A performance question of SDP in YALMIP

Cheng Gong

unread,

Aug 9, 2015, 10:07:14 AM8/9/15

to YALMIP

Hello!

Today I tried to convert the model of OPF Slover(byJavad Lavaei，http://www.ieor.berkeley.edu/~lavaei/Software.html),which is a SDP relaxition of large scale optimal power flow modeled by CVX,to YALMIP.

It runs correctly.But when the size of problem become larger ,its run time explode.

Here is his code:

--------------------------------------------------------------------------------------------------------------------------

cvx_begin

if strcmp(solver,'sdpt3')

cvx_solver sdpt3

cvx_solver_settings('maxit',120,'gaptol',1e-15,'inftol',1e-15,'steptol',1e-15);

elseif strcmp(solver,'sedumi')

cvx_solver sedumi

elseif strcmp(solver,'mosek')

cvx_solver mosek

end

cvx_precision best

variable VV(nlx,nc) complex

variable V2(nb,nc)

variable Sg(ng,nc) complex

expression sf(nl,nc)

expression st(nl,nc)

expression sfP(nl,nc)

expression stP(nl,nc)

expression Sb(nb,nc)

expression W(nb*nc,nb*nc)

for rr = 1 : nc

in = (nb*(rr-1)+1) : (nb*rr);

ex = setdiff(1 : nl, cc{rr});

W(in,in) = W(in,in) + sparse(fx, tx, VV(:,rr), nb, nb);

W(in,in) = W(in,in) + sparse(tx, fx, conj(VV(:,rr)), nb, nb);

W(in,in) = W(in,in) + sparse(1:nb, 1:nb, V2(:,rr), nb, nb);

sf(ex,rr) = conj(diag(Yf(ex, :) * W(in,in) * incidentF(ex, :).'));

st(ex,rr) = conj(diag(Yt(ex, :) * W(in,in) * incidentT(ex, :).'));

Sb(:,rr) = Pd + Qd*1i + conj(diag(Ybus{rr} * W(in,in)));

sfP(ex,rr) = conj(diag(YfP(ex, :) * W(in,in) * incidentF(ex, :).'));

stP(ex,rr) = conj(diag(YtP(ex, :) * W(in,in) * incidentT(ex, :).'));

sf(cc{rr},rr) = zeros(length(cc{rr}),1);

st(cc{rr},rr) = zeros(length(cc{rr}),1);

sfP(cc{rr},rr) = zeros(length(cc{rr}),1);

stP(cc{rr},rr) = zeros(length(cc{rr}),1);

end

obj = sum(c2.*real(Sg(:,1)).^2 + c1.*real(Sg(:,1)) + c0);

panG = sum(sum(imag(Sg)));

panQL = sum(sum(ndxL_prob .* abs(sfP + stP)));

cost = obj + ep(1)*panG + ep(2)*panQL;

minimize( cost );

subject to

for rr = 1 : nc

for kk = 1 : nBag

in = nb*(rr-1) + busN(bags{kk});

W(in,in) == hermitian_semidefinite(size(bags{kk},2));

end

Sb == incidentG.' * Sg;

V2 <= (Vmax.^2) * ones(1,nc);

V2 >= (Vmin.^2) * ones(1,nc);

abs(sf) .* edges <= SlmMax * ones(1,nc);

abs(st) .* edges <= SlmMax * ones(1,nc);

real(Sg) <= Pmax * ones(1,nc);

real(Sg) >= Pmin * ones(1,nc);

imag(Sg) <= Qmax * ones(1,nc);

imag(Sg) >= Qmin * ones(1,nc);

for rr = 2 : nc

abs(real(Sg(:,rr) - Sg(:,1))) <= correction;

end

cvx_end

--------------------------------------------------------------------------------------------------------------------------

here is mine

--------------------------------------------------------------------------------------------------------------------------

VV = sdpvar(nlx,nc,'full','complex');

V2 = sdpvar(nb,nc);

Sg = sdpvar(ng,nc,'full','complex');

sf = sdpvar(nl,nc,'full','complex');

st = sdpvar(nl,nc,'full','complex');

sfP = sdpvar(nl,nc,'full','complex');

stP = sdpvar(nl,nc,'full','complex');

Sb = sdpvar(nb,nc,'full','complex');

W = sdpvar(nb*nc,nb*nc,'hermitian','complex');

for rr = 1 : nc

in = (nb*(rr-1)+1) : (nb*rr);

ex = setdiff(1 : nl, cc{rr});

W(in,in) = W(in,in) + sparse(fx, tx, VV(:,rr), nb, nb);

W(in,in) = W(in,in) + sparse(tx, fx, conj(VV(:,rr)), nb, nb);

W(in,in) = W(in,in) + sparse(1:nb, 1:nb, V2(:,rr), nb, nb);

sf(ex,rr) = conj(diag(Yf(ex, :) * W(in,in) * incidentF(ex, :).'));

st(ex,rr) = conj(diag(Yt(ex, :) * W(in,in) * incidentT(ex, :).'));

Sb(:,rr) = Pd + Qd*1i + conj(diag(Ybus{rr} * W(in,in)));

sfP(ex,rr) = conj(diag(YfP(ex, :) * W(in,in) * incidentF(ex, :).'));

stP(ex,rr) = conj(diag(YtP(ex, :) * W(in,in) * incidentT(ex, :).'));

sf(cc{rr},rr) = zeros(length(cc{rr}),1);

st(cc{rr},rr) = zeros(length(cc{rr}),1);

sfP(cc{rr},rr) = zeros(length(cc{rr}),1);

stP(cc{rr},rr) = zeros(length(cc{rr}),1);

end

obj = sum(c2.*real(Sg(:,1)).^2 + c1.*real(Sg(:,1)) + c0);

panG = sum(sum(imag(Sg)));

panQL = sum(sum(ndxL_prob .* abs(sfP + stP)));

cost = obj + ep(1)*panG + ep(2)*panQL;

Constraints = [];

for rr = 1 : nc

for kk = 1 : nBag

in = nb*(rr-1) + busN(bags{kk});

Constraints = [Constraints,W(in,in) >= 0];

end

Constraints = [Constraints,Sb == incidentG.' * Sg];

Constraints = [Constraints,V2 <= (Vmax.^2) * ones(1,nc)];

Constraints = [Constraints,V2 >= (Vmin.^2) * ones(1,nc)];

Constraints = [Constraints,abs(sf) .* edges <= SlmMax * ones(1,nc)];

Constraints = [Constraints,abs(st) .* edges <= SlmMax * ones(1,nc)];

Constraints = [Constraints,real(Sg) <= Pmax * ones(1,nc)];

Constraints = [Constraints,real(Sg) >= Pmin * ones(1,nc)];

Constraints = [Constraints,imag(Sg) <= Qmax * ones(1,nc)];

Constraints = [Constraints,imag(Sg) >= Qmin * ones(1,nc)];

for rr = 2 : nc

Constraints = [Constraints,abs(real(Sg(:,rr) - Sg(:,1))) <= correction];

end

% options = sdpsettings('verbose',1,'solver','sdpt3');

results = optimize(Constraints,cost);

double(cost)

--------------------------------------------------------------------------------------------------------------------------

only the case<ieee 118 have the same magnitude of run time

i guess i have something wrong.

Thanks,

Gong.

Cheng Gong

unread,

Aug 9, 2015, 11:19:39 PM8/9/15

to YALMIP

Sorry,

I made a mistake.It does not run correctly.Althogh the result is close to the right one.

--------------------------------------------------------------------------------------------------------------------------

stop: primal infeas has deteriorated too much, 6.1e-05 0, 0, 1

Stop: progress is bad**

--------------------------------------------------------------------------------------------------------------------------

sorry for my rude and my poor english,i'll try again

and hope for your help

Thanks

Johan Löfberg

unread,

Aug 10, 2015, 2:54:38 AM8/10/15

to YALMIP

As a start, it looks to me that you think W will be what CVX calls an expression, i.e., a variable which really isn't a variable but only a place-holder for describing something.

In YALMIP, there are only variables, and nothing else. Hence, W is defined as a variable with a huge amount of variables, and then later you replace some of those variables with other. There will be a huge amount left from the original ones though (everything outside the "in" blocks).

One alternative which works it is simple block diagonal

W = [];

for rr = 1 : nc

        in = (nb*(rr-1)+1) : (nb*rr);
        ex = setdiff(1 : nl, cc{rr});

        block = sparse(fx, tx, VV(:,rr), nb, nb) + sparse(tx, fx, conj(VV(:,rr)), nb, nb) + sparse(1:nb, 1:nb, V2(:,rr), nb, nb);
        W = blkdiag(W, block);

Alternatively, you would think that you simply can replace the current definition by zeros, but unfortunately that does not work as MATLAB does not like to index an non-built-in class into an initial double. However, the following ugly piece of code works

W = double2sdpvar(zeros(nb*nc,nb*nc));

  for rr = 1 : nc
        in = (nb*(rr-1)+1) : (nb*rr);
        ex = setdiff(1 : nl, cc{rr});    
             
        W(in,in) = W(in,in) + sparse(fx, tx, VV(:,rr), nb, nb);
        W(in,in) = W(in,in) + sparse(tx, fx, conj(VV(:,rr)), nb, nb); 
        W(in,in) = W(in,in) + sparse(1:nb, 1:nb, V2(:,rr), nb, nb);

There are probably some other ways to do it, but those are two alternatives. The basic problem is that you can only index sdpvars into sdpvars (or into empty matrices). The basic pattern to code is to muild using concatentation, or start with some basic variable (which you do not have) and then add to that.

Hence, quick fix is to replace all "expression" with double2sdpvar(zeros(n,m))

Cheng Gong

unread,

Aug 10, 2015, 9:34:07 AM8/10/15

to YALMIP

It works,

Thanks for your answer .

I tried your second alternative,but the performance of formulate these place-holders also has considerable difference as the scale of the cases increasing.

case30 case118 case300

CVX : 0.032778 0.359951 1.477761

YALMIP: 0.074443 0.429680 15.966221

But it's enough to me.

Thanks again!

Gong.

Johan Löfberg

unread,

Aug 10, 2015, 9:38:59 AM8/10/15

to YALMIP

Don't understand what those numbers are. When I run case300, both cvx and yalmip take around 5-10 seconds

Johan Löfberg

unread,

Aug 10, 2015, 9:53:17 AM8/10/15

to YALMIP

btw, don't forget yalmip('clear') in your code, to make sure yalmip doesn't clog up internally with all old variables

Johan Löfberg

unread,

Aug 10, 2015, 9:56:53 AM8/10/15

to YALMIP

and these lines are redundant now as you initialized them to 0

Cheng Gong

unread,

Aug 11, 2015, 2:27:12 AM8/11/15

to YALMIP

The numbers I posted are result of this piece:

--------------------------------------------------------------------------------------------------------------------------

tic;

for rr = 1 : nc

in = (nb*(rr-1)+1) : (nb*rr);

ex = setdiff(1 : nl, cc{rr});

W(in,in) = W(in,in) + sparse(fx, tx, VV(:,rr), nb, nb);

W(in,in) = W(in,in) + sparse(tx, fx, conj(VV(:,rr)), nb, nb);

W(in,in) = W(in,in) + sparse(1:nb, 1:nb, V2(:,rr), nb, nb);

sf(ex,rr) = conj(diag(Yf(ex, :) * W(in,in) * incidentF(ex, :).'));

st(ex,rr) = conj(diag(Yt(ex, :) * W(in,in) * incidentT(ex, :).'));

Sb(:,rr) = Pd + Qd*1i + conj(diag(Ybus{rr} * W(in,in)));

sfP(ex,rr) = conj(diag(YfP(ex, :) * W(in,in) * incidentF(ex, :).'));

stP(ex,rr) = conj(diag(YtP(ex, :) * W(in,in) * incidentT(ex, :).'));

sf(cc{rr},rr) = zeros(length(cc{rr}),1);

st(cc{rr},rr) = zeros(length(cc{rr}),1);

sfP(cc{rr},rr) = zeros(length(cc{rr}),1);

stP(cc{rr},rr) = zeros(length(cc{rr}),1);

end

toc;

--------------------------------------------------------------------------------------------------------------------------

I run it on my old computer,with win7 x86 and MATLAB 2014a,in case 300.

The function sdpvar.subsasgn takes 13.154s.

But the total run time is YALMIP faster.

I hasn't set solver SDPT3 the same as CVX.The algorithms of SDPT3 is different.

Johan Löfberg

unread,

Aug 11, 2015, 2:35:33 AM8/11/15

to YALMIP

I would avoid the unnecessary data shuffling you are doing and write it as

temp = sparse(fx, tx, VV(:,rr), nb, nb) + sparse(tx, fx, conj(VV(:,rr)), nb, nb) + sparse(1:nb, 1:nb, V2(:,rr), nb, nb);        
W(in,in) = W(in,in) + temp;

and as I said, the last 4 rows are redundant

Cheng Gong

unread,

Aug 11, 2015, 7:57:12 AM8/11/15

to YALMIP

I've write it as you say and it gets faster.

Thanks,

Your answer has taught me a lot.

Reply all

Reply to author

Forward