The introduction of second-order cone relaxation constraints significantly increases the solution time

[ty72216]

Hello professor, in the process of modeling network optimization and reconstruction with the goal of maximizing power supply recovery and optimizing power flow as the main focus, the second-order cone relaxation constraint is not considered. This problem is modeled as a MIP, which can be solved within a few seconds. However, after introducing the second-order cone relaxation constraint, the solution time significantly increased, and even 6 hours were not enough to complete the solution. There may be some unreasonable aspects in the model construction.
%% Second-order cone relaxation
con=[con,V2(:,1).*I>=Pij.^2+Qij.^2];
May I ask you how to improve the model to speed up the solving process and obtain the global optimal solution as much as possible?
Here is my code:

%%
clc
clear
V2_up=1.05^2;
load xSWij
tic
%% Parameters
NN=33;%Number of nodes
NL=37;%Number of branches (including tie lines)
NL1=32;
M=10000;
UB = 10; % Voltage base kV
SB = 10; %Power base MVA
ZB = UB^2/SB; % Impedance base ohm
for readdata=1:1
% Line impedances
BranchData=[1,2,0.00575259116172393,0.00293244885684409;2,3,0.0307595167324284,0.0156667639990117;3,4,0.0228356655660625,0.0116299673811859;4,5,0.0237777927519847,0.0121103898534774;5,6,0.0510994811437299,0.0441115179103993;6,7,0.0116798814042811,0.0386084968641515;7,8,0.0443860450374230,0.0146684835371073;8,9,0.0642643047350938,0.0461704713630771;9,10,0.0651378001392601,0.0461704713630771;10,11,0.0122663711756499,0.00405551437648650;11,12,0.0233597628085623,0.00772419507398506;12,13,0.0915922323797259,0.0720633708437217;13,14,0.0337917936354629,0.0444796338307266;14,15,0.0368739845615927,0.0328184701851062;15,16,0.0465635442949519,0.0340039282336176;16,17,0.0804239697121708,0.107377542183589;17,18,0.0456713311321249,0.0358133115708193;2,19,0.0102323747345198,0.00976443076800212;19,20,0.0938508419247845,0.0845668336290739;20,21,0.0255497405718650,0.0298485858109407;21,22,0.0442300637152505,0.0584805173089354;3,23,0.0281515090257032,0.0192356166503198;23,24,0.0560284909243828,0.0442425422210243;24,25,0.0559037058666447,0.0437434019900721;6,26,0.0126656833604117,0.00645138748505699;26,27,0.0177319567045764,0.00902819892734764;27,28,0.0660736880722955,0.0582559042050069;28,29,0.0501760717164684,0.0437122057256376;29,30,0.0316642084010292,0.0161284687126425;30,31,0.0607952801299761,0.0600840053008693;31,32,0.0193728802138317,0.0225798561976995;32,33,0.0212758523443369,0.0330805188063561;8,21,0.124800000000000,0.124800000000000;9,15,0.124800000000000,0.124800000000000;12,22,0.124800000000000,0.124800000000000;18,33,0.0312000000000000,0.0312000000000000;25,29,0.0312000000000000,0.0312000000000000;];
% Predicted loads for each node
NodeData=[1,0,0;2,0.0100000000000000,0.00600000000000000;3,0.00900000000000000,0.00400000000000000;4,0.0120000000000000,0.00800000000000000;5,0.00600000000000000,0.00300000000000000;6,0.00600000000000000,0.00200000000000000;7,0.0200000000000000,0.0100000000000000;8,0.0200000000000000,0.0100000000000000;9,0.00600000000000000,0.00200000000000000;10,0.00600000000000000,0.00200000000000000;11,0.00450000000000000,0.00300000000000000;12,0.00600000000000000,0.00350000000000000;13,0.00600000000000000,0.00350000000000000;14,0.0120000000000000,0.00800000000000000;15,0.00600000000000000,0.00100000000000000;16,0.00600000000000000,0.00200000000000000;17,0.00600000000000000,0.00200000000000000;18,0.00900000000000000,0.00400000000000000;19,0.00900000000000000,0.00400000000000000;20,0.00900000000000000,0.00400000000000000;21,0.00900000000000000,0.00400000000000000;22,0.00900000000000000,0.00400000000000000;23,0.00900000000000000,0.00500000000000000;24,0.0420000000000000,0.0200000000000000;25,0.0420000000000000,0.0200000000000000;26,0.00600000000000000,0.00250000000000000;27,0.00600000000000000,0.00250000000000000;28,0.00600000000000000,0.00200000000000000;29,0.0120000000000000,0.00700000000000000;30,0.0200000000000000,0.0600000000000000;31,0.0150000000000000,0.00700000000000000;32,0.0210000000000000,0.0100000000000000;33,0.00600000000000000,0.00400000000000000;];
end
ntime=1; % Single time slice
pl=zeros(33,ntime);ql=zeros(33,ntime);
for i=1:33
pl(i,:)=NodeData(i,2);
ql(i,:)=NodeData(i,3);
end
AM=zeros(NN);
for i=1:NL
AM(BranchData(i,1),BranchData(i,2))=1;
AM(BranchData(i,2),BranchData(i,1))=1;
end % Adjacency matrix, if ij are connected, AM(i,j)=AM(j,i)=1;
% However, it is clear that the direction of all branches is from nodes with smaller numbers to larger ones
% Thus, for node i, looking at row i, those before AM(i,i) are flowing towards node i, and those after are flowing out from node i
r=BranchData(:,3); %Resistance
x=BranchData(:,4); %Reactance
DGn=3; % number of DG
% Nodes where DGs can be connected
DGcan=[7,16,23];
w=ones(NN,1);
a=[1 10 11 17 27];
b=[4 6 7 18 19 21 23 26 32];
a=a+1;
b=b+1;
for i=1:NN
if ismember(i,a)
w(i)=100;
else if ismember(i,b)
w(i)=10;
end
end
end

%% Variables
% Control variables under fault scenarios
DGp=sdpvar(DGn,1,'full'); % DG active power. Dimension: DGn nodes; Power values are in per unit.
DGq=sdpvar(DGn,1,'full');
Ah=binvar(DGn,ntime,'full');

Pij=sdpvar(NL,1,'full'); % Branch active power, from nodes with smaller numbers to larger ones
Qij=sdpvar(NL,1,'full'); % Branch reactive power, from nodes with smaller numbers to larger ones
V1=sdpvar(NN,1,'full'); % Substation node voltage
V2=sdpvar(NL,2,'full'); % Load node voltage squared terms
I=sdpvar(NL,1,'full'); % Branch load current squared terms

fl_PF = sdpvar(NL,1,'full'); % Virtual fault flow
fn_PF = sdpvar(NN,1,'full');
% Switch state variables
sxyi = binvar(NL,1,'full');
sxyj = binvar(NL,1,'full');

% Recovery of node load after failure
pll=sdpvar(NN,1,'full');
qll=sdpvar(NN,1,'full');

puti=sdpvar(1,1,'full'); % Active power absorbed by the substation node from the upper-level grid
quti=sdpvar(1,1,'full'); % Reactive power absorbed by the substation node from the upper-level grid

x1 = sdpvar(NN,1,'full'); % Optimal active load shedding
y1 = sdpvar(NN,1,'full'); % Optimal reactive load shedding
q = binvar(NN,1,'full'); % Divide into island node variables
lxy=binvar(NL,1,'full'); % Branch disconnection state variable
q1=sdpvar(DGn,1,'full'); %
Fij=sdpvar(37,ntime,'full'); % Virtual Commodity flow

% Definition of virtual voltage upper limit without considering current
V1_sup=sdpvar(NN,1,'full');
V2_sup=sdpvar(NL,2,'full');
Pij_sup=sdpvar(NL,1,'full');
Qij_sup=sdpvar(NL,1,'full');
DGp_sup=sdpvar(DGn,1,'full');
DGq_sup=sdpvar(DGn,1,'full');
puti_sup=sdpvar(1,1,'full');
quti_sup=sdpvar(1,1,'full');

%% Define virtual fault flow constraints
con = [];
for i = 1:NL
if xSWij(i,3)==1
con = [con,(sxyi(i,1)-1)*M+fn_PF(xSWij(i,1),1)<=fl_PF(i,1)];
con = [con,(1-sxyi(i,1))*M+fn_PF(xSWij(i,1),1)>=fl_PF(i,1)];
end
if xSWij(i,4)==1
con = [con,(sxyj(i,1)-1)*M+fn_PF(xSWij(i,2),1)<=fl_PF(i,1)];
con = [con,(1-sxyj(i,1))*M+fn_PF(xSWij(i,2),1)>=fl_PF(i,1)];
end
end
%Substation nodes are not affected by faults
con = [con,fn_PF(1,:)==1];
con=[con,V1(1,:)==1]; %Select the voltage of node 1 as the standard voltage
%% DistFlow
%Branch voltage constraint
for k=2:NN
con=[con,V2_up>=V1(k,:)];
con=[con,0.95^2<=V1(k,:)];
end

for i = 1:NN
for j = 1:NN
for l = 1:NL
if BranchData(l,1)==i&&BranchData(l,2)==j
con = [con,V2(l,1)>=0];
con = [con,V2(l,1)<=lxy(l,1)*M];
con = [con,V2(l,1)>=(lxy(l,1)-1)*M+V1(i,1)];
con = [con,V2(l,1)<=V1(i,1)];
con = [con,V2(l,2)>=0];
con = [con,V2(l,2)<=lxy(l,1)*M];
con = [con,V2(l,2)>=(lxy(l,1)-1)M+V1(j,1)];
con = [con,V2(l,2)<=V1(j,1)];
con = [con,V2(l,2)==V2(l,1)-2(r(l,1)*Pij(l,1)+x(l,1)*Qij(l,1))+(r(l,1)^2+x(l,1)^2)*I(l,1)];
end
end
end
end

%% Second-order cone relaxation
con=[con,V2(:,1).*I>=Pij.^2+Qij.^2];
%%

for i=1:NN
Pnode=0; %Inflow-outflow=0
Qnode=0;
Pnode=Pnode-pll(i);
Qnode=Qnode-qll(i);
for j=1:NN
if AM(i,j)==1
if j<i
for l=1:NL
if (BranchData(l,1)==j)&&(BranchData(l,2)==i)
ll=l;
end
end
Pnode=Pnode+Pij(ll)-r(ll)*I(ll);
Qnode=Qnode+Qij(ll)-x(ll)*I(ll);
end
if j>i
for l=1:NL
if (BranchData(l,1)==i)&&(BranchData(l,2)==j)
ll=l;
end
end
Pnode=Pnode-Pij(ll)-r(ll)*I(ll);
Qnode=Qnode-Qij(ll)-x(ll)*I(ll);
end
end
end

if i==1
    Pnode=Pnode+puti;
    Qnode=Qnode+quti;
end
for j=1:DGn
    if i==DGcan(j) 
        Pnode=Pnode+DGp(j);
        Qnode=Qnode+DGq(j);

    end
end
con=[con,Pnode==0];
con=[con,Qnode==0];

end

% Branch power and current transmission constraints

con=[con,(Pij(2:NL,1)<=0.5)&(Pij(2:NL,1)>=-0.5)];
con=[con,(Qij(2:NL,1)<=0.4)&(Qij(2:NL,1)>=-0.4)];
con=[con,(I<=0.46)&(I>=0)];

con=[con,(Pij(2:NL,1)<=0.5)&(Pij(1,1)>=0)];
on=[con,(Qij(2:NL,1)<=0.4)&(Qij(1,1)>=0)];

con = [con,DGp<=0.1,DGp>=0];
con = [con,DGq<=0.01,DGq>=0];

%% Linear relaxation of optimal cutting load at nodes

for i = 1:NN
con = [con,x1(i,1)>=0];
con = [con,x1(i,1)<=1];
end
for i = 1:NN
con = [con,y1(i,1)>=0];
con = [con,y1(i,1)<=1];
end

con = [con,q == fn_PF];
con = [con,pll>=pl.(x1-(1-q)M),pll<=x1.pl];
con = [con,pll<=qM,pll>=-qM];
con = [con,qll>=ql.(y1-(1-q)M),qll<=y1.ql];
con = [con,qll<=qM,qll>=-qM];

con=[con,-Msxyi<=Pij,Msxyi>=Pij];
con=[con,-Msxyj<=Pij,Msxyj>=Pij];
con=[con,-Msxyi<=Qij,Msxyi>=Qij];
con=[con,-Msxyj<=Qij,Msxyj>=Qij];
con=[con,-Msxyi<=I,Msxyi>=I];
con=[con,-Msxyj<=I,Msxyj>=I];

con=[con,(fl_PF(:,:)>=0)&(fl_PF(:,:)<=1)];
con=[con,(puti<=0.5)&(puti>=0)];
con=[con,(quti<=0.5)&(quti>=0)];

con = [con,sum(pll)==sum(DGp)+puti-sum(I.*r)];
% con = [con,sum(qll)==sum(DGq)+quti-sum(I.*x)];

%% Radiative constraint

for k=2:33
    if k ~=DGcan
        node_out=find(BranchData(:,1)==k);
        node_in=find(BranchData(:,2)==k);
        con=[con,sum(Fij(node_in))-sum(Fij(node_out))== x1(k)];
    end

end
for k = 2:33
   node_out=find(BranchData(:,1)==k);
   node_in=find(BranchData(:,2)==k);
    for j=1:DGn
    if k==DGcan(j)  
       con=[con,sum(Fij(node_in))-sum(Fij(node_out))<=M*Ah(j,1)+x1(k)];
       con=[con,sum(Fij(node_in))-sum(Fij(node_out))>=-M*Ah(j,1)+x1(k)];
    end
   end
end

con=[con,lxy(:)<=sxyi(:)];
con=[con,lxy(:)<=sxyj(:)];
con=[con,lxy(:)>=sxyi(:)+sxyj(:)-1];
con=[con,sum(q1)==sum(x1)];
con=[con,Fij<=sum(x1)];
con=[con,-Fij<=sum(x1)];
con = [con,sum(lxy)==sum(q(2:33))-sum(Ah)];

for i = 1:NL
con=[con,-Mlxy(i)<=Fij(i)];
con=[con,Fij(i)<=Mlxy(i)];
end
% define fault
con = [con,fl_PF(1,1)==0];
%% Voltage tight constraint
con=[con,V1_sup(1,:)==1];
for k=2:NN
con=[con,V2_up>=V1_sup(k,:)];
con=[con,0.95^2<=V1_sup(k,:)];
end
for i = 1:NN
for j = 1:NN
for l = 1:NL
if BranchData(l,1)==i&&BranchData(l,2)==j
con = [con,V2_sup(l,1)>=0];
con = [con,V2_sup(l,1)<=lxy(l,1)*M];
con = [con,V2_sup(l,1)>=(lxy(l,1)-1)*M+V1_sup(i,1)];
con = [con,V2_sup(l,1)<=V1_sup(i,1)];
con = [con,V2_sup(l,2)>=0];
con = [con,V2_sup(l,2)<=lxy(l,1)*M];
con = [con,V2_sup(l,2)>=(lxy(l,1)-1)M+V1_sup(j,1)];
con = [con,V2_sup(l,2)<=V1_sup(j,1)];
con = [con,V2_sup(l,2)==V2_sup(l,1)-2(r(l,1)*Pij(l,1)+x(l,1)*Qij(l,1))];
end
end
end
end

for i=1:NN
Pnode=0;
Qnode=0;
Pnode=Pnode-pll(i);
Qnode=Qnode-qll(i);
for j=1:NN
if AM(i,j)==1
if j<i
for l=1:NL
if (BranchData(l,1)==j)&&(BranchData(l,2)==i)
ll=l;
end
end
Pnode=Pnode+Pij_sup(ll);
Qnode=Qnode+Qij_sup(ll);
end
if j>i
for l=1:NL
if (BranchData(l,1)==i)&&(BranchData(l,2)==j)
ll=l;
end
end
Pnode=Pnode-Pij_sup(ll);
Qnode=Qnode-Qij_sup(ll);
end
end
end

if i==1
    Pnode=Pnode+puti_sup;
    Qnode=Qnode+quti_sup;
end
for j=1:DGn
    if i==DGcan(j) 
        Pnode=Pnode+DGp_sup(j);
        Qnode=Qnode+DGq_sup(j);

    end
end
con=[con,Pnode==0];
con=[con,Qnode==0];

end
%Branch power constraint
%Substation connection branch
con=[con,Pij_sup(1,1)<=0.3715,Pij_sup(1,1)>=0];
con=[con,Qij_sup(1,1)<=0.23,Qij_sup(1,1)>=0];

con=[con,Pij_sup(2:NL,1)<=0.3715,Pij_sup(2:NL,1)>=-0.3715];
con=[con,Qij_sup(2:NL,1)<=0.23,Qij_sup(2:NL,1)>=-0.23];

%DG output constraint
con = [con,DGp_sup>=0,DGp_sup<=0.1];
con = [con,DGq_sup>=0,DGq_sup<=0.01];
%Substation output constraints
con=[con,(puti_sup<=0.3715)&(puti_sup>=0)];
con=[con,(quti_sup<=0.3715)&(quti_sup>=0)];
%Power balance constraint
con = [con,sum(pll)==sum(DGp_sup)+puti_sup];
con = [con,sum(qll)==sum(DGq_sup)+quti_sup];

con=[con,-Msxyi<=Pij_sup,Msxyi>=Pij_sup];
con=[con,-Msxyj<=Pij_sup,Msxyj>=Pij_sup];
con=[con,-Msxyi<=Qij_sup,Msxyi>=Qij_sup];
con=[con,-Msxyj<=Qij_sup,Msxyj>=Qij_sup];

con=[con,V2<=V2_sup];
%% obj
obj=sum(-pll.*w)+sum(q)*0.0001; %Node load restored to maximum power supply
%% Solving models
ops=sdpsettings('verbose', 1, 'solver', 'cplex');
sol=optimize(con,obj,ops);
ops.cplex.timelimit=10800;
if sol.problem == 0
disp('succcessful solved');
else
disp('error');
yalmiperror(sol.problem)
end
toc
%%
Pij=value(Pij);
x1=value(x1);
q=value(q);
Ah=value(Ah);
pll=value(pll);
V2=value(V2);
sxyi=value(sxyi);
sxyj=value(sxyj);
Fij=value(Fij);
lxy=value(lxy);
DGp=value(DGp);
I=value(I);
V2_sup=value(V2_sup);

[johanlofberg]

MISOCP is much less mature than MILP and can thus take much longer to solve

[ty72216]

Are there any ways to accelerate the speed of MISOCP solving?

[johanlofberg]

Nothing magic. Try other solver (mosek/gurobi) and all the algorithmic options available, or try various high-level strategies such as iteratively adding violated SOCP cuts

[ty72216]

Thank you for your reply!

[ty72216]

Installing the latest version of mosek but without the SCOP module, the solver(mosek) not applicable . Is it because the yalmip version is incompatible?

[johanlofberg]

You must use an old version of YALMIP, as I don't recognize the print-out from current version.

Also, if you try to use mosek, but haven't installed the SOCP module, it will fail, of course.