Description of sdpopf

0001 function [results,success,raw] = sdpopf_solver(om, mpopt)
0002 %SDPOPF_SOLVER A semidefinite programming relaxtion of the OPF problem
0003 %
0004 %   [RESULTS,SUCCESS,RAW] = SDPOPF_SOLVER(OM, MPOPT)
0005 %
0006 %   Inputs are an OPF model object and a MATPOWER options vector.
0007 %
0008 %   Outputs are a RESULTS struct, SUCCESS flag and RAW output struct.
0009 %
0010 %   RESULTS is a MATPOWER case struct (mpc) with the usual baseMVA, bus
0011 %   branch, gen, gencost fields, along with the following additional
0012 %   fields:
0013 %       .f           final objective function value
0014 %       .mineigratio Minimum eigenvalue ratio measure of rank condition
0015 %                    satisfaction
0016 %       .zero_eval   Zero eigenvalue measure of optimal voltage profile
0017 %                    consistency
0018 %
0019 %   SUCCESS     1 if solver converged successfully, 0 otherwise
0020 %
0021 %   RAW         raw output
0022 %       .xr     final value of optimization variables
0023 %           .A  Cell array of matrix variables for the dual problem
0024 %           .W  Cell array of matrix variables for the primal problem
0025 %       .pimul  constraint multipliers on
0026 %           .lambdaeq_sdp   Active power equalities
0027 %           .psiU_sdp       Generator active power upper inequalities
0028 %           .psiL_sdp       Generator active power lower inequalities
0029 %           .gammaeq_sdp    Reactive power equalities
0030 %           .gammaU_sdp     Reactive power upper inequalities
0031 %           .gammaL_sdp     Reactive power lower inequalities
0032 %           .muU_sdp        Square of voltage magnitude upper inequalities
0033 %           .muL_sdp        Square of voltage magnitude lower inequalities
0034 %       .info   solver specific termination code
0035 %
0036 %   See also OPF.
0037 
0038 %   MATPOWER
0039 %   Copyright (c) 2013-2019, Power Systems Engineering Research Center (PSERC)
0040 %   by Daniel Molzahn, PSERC U of Wisc, Madison
0041 %   and Ray Zimmerman, PSERC Cornell
0042 %
0043 %   This file is part of MATPOWER/mx-sdp_pf.
0044 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0045 %   See https://github.com/MATPOWER/mx-sdp_pf/ for more info.
0046 
0047 % Unpack options
0048 ignore_angle_lim    = mpopt.opf.ignore_angle_lim;
0049 verbose             = mpopt.verbose;
0050 maxNumberOfCliques  = mpopt.sdp_pf.max_number_of_cliques;   %% Maximum number of maximal cliques
0051 eps_r               = mpopt.sdp_pf.eps_r;               %% Resistance of lines when zero line resistance specified
0052 recover_voltage     = mpopt.sdp_pf.recover_voltage;     %% Method for recovering the optimal voltage profile
0053 recover_injections  = mpopt.sdp_pf.recover_injections;  %% Method for recovering the power injections
0054 minPgendiff         = mpopt.sdp_pf.min_Pgen_diff;       %% Fix Pgen to midpoint of generator range if PMAX-PMIN < MINPGENDIFF (in MW)
0055 minQgendiff         = mpopt.sdp_pf.min_Qgen_diff;       %% Fix Qgen to midpoint of generator range if QMAX-QMIN < MINQGENDIFF (in MVAr)
0056 maxlinelimit        = mpopt.sdp_pf.max_line_limit;      %% Maximum line limit. Anything above this will be unconstrained. Line limits of zero are also unconstrained.
0057 maxgenlimit         = mpopt.sdp_pf.max_gen_limit;       %% Maximum reactive power generation limits. If Qmax or abs(Qmin) > maxgenlimit, reactive power injection is unlimited.
0058 ndisplay            = mpopt.sdp_pf.ndisplay;            %% Determine display frequency of diagonastic information
0059 cholopts.dense      = mpopt.sdp_pf.choldense;           %% Cholesky factorization options
0060 cholopts.aggressive = mpopt.sdp_pf.cholaggressive;      %% Cholesky factorization options
0061 binding_lagrange    = mpopt.sdp_pf.bind_lagrange;       %% Tolerance for considering a Lagrange multiplier to indicae a binding constraint
0062 zeroeval_tol        = mpopt.sdp_pf.zeroeval_tol;        %% Tolerance for considering an eigenvalue in LLeval equal to zero
0063 mineigratio_tol     = mpopt.sdp_pf.mineigratio_tol;     %% Tolerance for considering the rank condition satisfied
0064 
0065 if upper(mpopt.opf.flow_lim(1)) ~= 'S' && upper(mpopt.opf.flow_lim(1)) ~= 'P'
0066     error('sdpopf_solver: Only ''S'' and ''P'' options are currently implemented for MPOPT.opf.flow_lim when MPOPT.opf.ac.solver == ''SDPOPF''');
0067 end
0068 
0069 % set YALMIP options struct in SDP_PF (for details, see help sdpsettings)
0070 sdpopts = yalmip_options([], mpopt);
0071 
0072 %% define undocumented MATLAB function ismembc() if not available (e.g. Octave)
0073 if exist('ismembc')
0074     ismembc_ = @ismembc;
0075 else
0076     ismembc_ = @ismembc_octave;
0077 end
0078 
0079 if verbose > 0
0080     v = sdp_pf_ver('all');
0081     fprintf('SDPOPF Version %s, %s', v.Version, v.Date);
0082     if ~isempty(sdpopts.solver)
0083         fprintf('  --  Using %s.\n\n',upper(sdpopts.solver));
0084     else
0085         fprintf('\n\n');
0086     end
0087 end
0088 
0089 tic;
0090 
0091 %% define named indices into data matrices
0092 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0093     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0094 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0095     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0096     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0097 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0098     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0099     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0100 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
0101 
0102 
0103 %% Load mpc data
0104 
0105 mpc = get_mpc(om);
0106 
0107 if toggle_dcline(mpc, 'status')
0108     error('sdpopf_solver: DC lines are not implemented in SDP_PF');
0109 end
0110 
0111 nbus = size(mpc.bus,1);
0112 ngen = size(mpc.gen,1);
0113 nbranch = size(mpc.branch,1);
0114 bustype = mpc.bus(:,BUS_TYPE);
0115 
0116 if any(mpc.gencost(:,NCOST) > 3 & mpc.gencost(:,MODEL) == POLYNOMIAL)
0117     error('sdpopf_solver: SDPOPF is limited to quadratic cost functions.');
0118 end
0119 
0120 if size(mpc.gencost,1) > ngen && verbose >= 1 % reactive power costs are specified
0121     warning('sdpopf_solver: Reactive power costs are not implemented in SDPOPF. Ignoring reactive power costs.');
0122 end
0123 
0124 if ~ignore_angle_lim && (any(mpc.branch(:,ANGMIN) ~= -360) || any(mpc.branch(:,ANGMAX) ~= 360))
0125     warning('sdpopf_solver: Angle difference constraints are not implemented in SDPOPF. Ignoring angle difference constraints.');
0126 end
0127 
0128 % Enforce a minimum resistance for all branches
0129 mpc.branch(:,BR_R) = max(mpc.branch(:,BR_R),eps_r); 
0130 
0131 % Some of the larger system (e.g., case2746wp) have generators
0132 % corresponding to buses that have bustype == PQ. Change these
0133 % to PV buses.
0134 for i=1:ngen
0135     busidx = find(mpc.bus(:,BUS_I) == mpc.gen(i,GEN_BUS));
0136     if isempty(busidx) || ~(bustype(busidx) == PV || bustype(busidx) == REF)
0137         bustype(busidx) = PV;
0138     end
0139 end
0140 
0141 % Buses may be listed as PV buses without associated generators. Change
0142 % these buses to PQ.
0143 for i=1:nbus
0144     if bustype(i) == PV
0145         genidx = find(mpc.gen(:,GEN_BUS) == mpc.bus(i,BUS_I), 1);
0146         if isempty(genidx)
0147             bustype(i) = PQ;
0148         end
0149     end
0150 end
0151 
0152 % if any(mpc.gencost(:,MODEL) == PW_LINEAR)
0153 %     error('sdpopf_solver: Piecewise linear generator costs are not implemented in SDPOPF.');
0154 % end
0155 
0156 % The code does not handle the case where a bus has both piecewise-linear
0157 % generator cost functions and quadratic generator cost functions.
0158 for i=1:nbus
0159     genidx = find(mpc.gen(:,GEN_BUS) == mpc.bus(i,BUS_I));
0160     if ~all(mpc.gencost(genidx,MODEL) == PW_LINEAR) && ~all(mpc.gencost(genidx,MODEL) == POLYNOMIAL)
0161         error('sdpopf_solver: Bus %i has generators with both piecewise-linear and quadratic generators, which is not supported in this version of SDPOPF.',i);
0162     end
0163 end
0164 
0165 % Save initial bus indexing in a new column on the right side of mpc.bus
0166 mpc.bus(:,size(mpc.bus,2)+1) = mpc.bus(:,BUS_I);
0167 
0168 tsetup = toc;
0169 
0170 tic;
0171 %% Determine Maximal Cliques
0172 %% Step 1: Cholesky factorization to obtain chordal extension
0173 % Use a minimum degree permutation to obtain a sparse chordal extension.
0174 
0175 if maxNumberOfCliques ~= 1
0176     [Ainc] = makeIncidence(mpc.bus,mpc.branch);
0177     sparsity_pattern = abs(Ainc.'*Ainc) + speye(nbus,nbus);
0178     per = amd(sparsity_pattern,cholopts);
0179     [Rchol,p] = chol(sparsity_pattern(per,per),'lower');
0180 
0181     if p ~= 0
0182         error('sdpopf_solver: sparsity_pattern not positive definite!');
0183     end
0184 else
0185     per = 1:nbus;
0186 end
0187 
0188 % Rearrange mpc to the same order as the permutation per
0189 mpc.bus = mpc.bus(per,:);
0190 mpc.bus(:,BUS_I) = 1:nbus;
0191 bustype = bustype(per);
0192 
0193 for i=1:ngen
0194     mpc.gen(i,GEN_BUS) = find(per == mpc.gen(i,GEN_BUS));
0195 end
0196 
0197 [mpc.gen genidx] = sortrows(mpc.gen,1);
0198 mpc.gencost = mpc.gencost(genidx,:);
0199 
0200 for i=1:nbranch
0201     mpc.branch(i,F_BUS) = find(per == mpc.branch(i,F_BUS));
0202     mpc.branch(i,T_BUS) = find(per == mpc.branch(i,T_BUS));
0203 end
0204 % -------------------------------------------------------------------------
0205 
0206 Sbase = mpc.baseMVA;
0207 
0208 Pd = mpc.bus(:,PD) / Sbase;
0209 Qd = mpc.bus(:,QD) / Sbase;
0210 
0211 for i=1:nbus  
0212     genidx = find(mpc.gen(:,GEN_BUS) == i);
0213     if ~isempty(genidx)
0214         Qmax(i) = sum(mpc.gen(genidx,QMAX)) / Sbase;
0215         Qmin(i) = sum(mpc.gen(genidx,QMIN)) / Sbase;
0216     else
0217         Qmax(i) = 0;
0218         Qmin(i) = 0;
0219     end
0220 end
0221 
0222 Vmax = mpc.bus(:,VMAX);
0223 Vmin = mpc.bus(:,VMIN);
0224 
0225 Smax = mpc.branch(:,RATE_A) / Sbase;
0226 
0227 % For generators with quadratic cost functions (handle piecewise linear
0228 % costs elsewhere)
0229 c2 = zeros(ngen,1);
0230 c1 = zeros(ngen,1);
0231 c0 = zeros(ngen,1);
0232 for genidx=1:ngen
0233     if mpc.gencost(genidx,MODEL) == POLYNOMIAL
0234         if mpc.gencost(genidx,NCOST) == 3 % Quadratic cost function
0235             c2(genidx) = mpc.gencost(genidx,COST) * Sbase^2;
0236             c1(genidx) = mpc.gencost(genidx,COST+1) * Sbase^1;
0237             c0(genidx) = mpc.gencost(genidx,COST+2) * Sbase^0;
0238         elseif mpc.gencost(genidx,NCOST) == 2 % Linear cost function
0239             c1(genidx) = mpc.gencost(genidx,COST) * Sbase^1;
0240             c0(genidx) = mpc.gencost(genidx,COST+1) * Sbase^0;
0241         else
0242             error('sdpopf_solver: Only piecewise-linear and quadratic cost functions are implemented in SDPOPF');
0243         end
0244     end
0245 end
0246 
0247 if any(c2 < 0)
0248     error('sdpopf_solver: Quadratic term of generator cost function is negative. Must be convex (non-negative).');
0249 end
0250 
0251 maxlinelimit = maxlinelimit / Sbase;
0252 maxgenlimit = maxgenlimit / Sbase;
0253 minPgendiff = minPgendiff / Sbase;
0254 minQgendiff = minQgendiff / Sbase;
0255 
0256 % Create Ybus with new bus numbers
0257 [Y, Yf, Yt] = makeYbus(mpc);
0258 
0259 
0260 %% Step 2: Compute maximal cliques of chordal extension
0261 % Use adjacency matrix made from Rchol
0262 
0263 if maxNumberOfCliques ~= 1 && nbus > 3
0264     
0265     % Build adjacency matrix
0266     [f,t] = find(Rchol - diag(diag(Rchol)));
0267     Aadj = sparse(f,t,ones(length(f),1),nbus,nbus);
0268     Aadj = max(Aadj,Aadj');
0269 
0270     % Determine maximal cliques
0271     [MC,ischordal] = maxCardSearch(Aadj);
0272     
0273     if ~ischordal
0274         % This should never happen since the Cholesky decomposition should
0275         % always yield a chordal graph.
0276         error('sdpopf_solver: Chordal extension adjacency matrix is not chordal!');
0277     end
0278     
0279     for i=1:size(MC,2)
0280        maxclique{i} = find(MC(:,i));
0281     end
0282 
0283 else
0284     maxclique{1} = 1:nbus;
0285     nmaxclique = 1;
0286     E = [];
0287 end
0288 
0289 
0290 %% Step 3: Make clique tree and combine maximal cliques
0291 
0292 if maxNumberOfCliques ~= 1 && nbus > 3
0293     % Create a graph of maximal cliques, where cost between each maximal clique
0294     % is the number of shared buses in the corresponding maximal cliques.
0295     nmaxclique = length(maxclique);
0296     cliqueCost = sparse(nmaxclique,nmaxclique);
0297     for i=1:nmaxclique
0298         maxcliquei = maxclique{i};
0299         for k=i+1:nmaxclique
0300 
0301             cliqueCost(i,k) = sum(ismembc_(maxcliquei,maxclique{k}));
0302 
0303             % Slower alternative that doesn't use undocumented MATLAB function
0304             % cliqueCost(i,k) = length(intersect(maxcliquei,maxclique{k}));
0305 
0306         end
0307     end
0308     cliqueCost = max(cliqueCost,cliqueCost.');
0309 
0310     % Calculate the maximal spanning tree
0311     cliqueCost = -cliqueCost;
0312     [E] = prim(cliqueCost);
0313 
0314     % Combine maximal cliques
0315     if verbose >= 2
0316         [maxclique,E] = combineMaxCliques(maxclique,E,maxNumberOfCliques,ndisplay);
0317     else
0318         [maxclique,E] = combineMaxCliques(maxclique,E,maxNumberOfCliques,inf);
0319     end
0320     nmaxclique = length(maxclique);
0321 end
0322 tmaxclique = toc;
0323 
0324 %% Create SDP relaxation
0325 
0326 tic;
0327 
0328 constraints = [];
0329 
0330 % Control growth of Wref variables
0331 dd_blk_size = 50*nbus;
0332 dq_blk_size = 2*dd_blk_size;
0333 
0334 Wref_dd = zeros(dd_blk_size,2); % For terms like Vd1*Vd1 and Vd1*Vd2
0335 lWref_dd = dd_blk_size;
0336 matidx_dd = zeros(dd_blk_size,3);
0337 dd_ptr = 0;
0338 
0339 Wref_dq = zeros(dq_blk_size,2); % For terms like Vd1*Vd1 and Vd1*Vd2
0340 lWref_dq = dq_blk_size;
0341 matidx_dq = zeros(dq_blk_size,3);
0342 dq_ptr = 0;
0343 
0344 for i=1:nmaxclique
0345     nmaxcliquei = length(maxclique{i});
0346     for k=1:nmaxcliquei % row
0347         for m=k:nmaxcliquei % column
0348             
0349             % Check if this pair loop{i}(k) and loop{i}(m) has appeared in
0350             % any previous maxclique. If not, it isn't in Wref_dd.
0351             Wref_dd_found = 0;
0352             if i > 1
0353                 for r = 1:i-1
0354                    if any(maxclique{r}(:) == maxclique{i}(k)) && any(maxclique{r}(:) == maxclique{i}(m))
0355                        Wref_dd_found = 1;
0356                        break;
0357                    end
0358                 end
0359             end
0360 
0361             if ~Wref_dd_found
0362                 % If we didn't find this element already, add it to Wref_dd
0363                 % and matidx_dd
0364                 dd_ptr = dd_ptr + 1;
0365                 
0366                 if dd_ptr > lWref_dd
0367                     % Expand the variables by dd_blk_size
0368                     Wref_dd = [Wref_dd; zeros(dd_blk_size,2)];
0369                     matidx_dd = [matidx_dd; zeros(dd_ptr-1 + dd_blk_size,3)];
0370                     lWref_dd = length(Wref_dd(:,1));
0371                 end
0372                 
0373                 Wref_dd(dd_ptr,1:2) = [maxclique{i}(k) maxclique{i}(m)];
0374                 matidx_dd(dd_ptr,1:3) = [i k m];
0375             end
0376 
0377             Wref_dq_found = Wref_dd_found;
0378 
0379             if ~Wref_dq_found
0380                 % If we didn't find this element already, add it to Wref_dq
0381                 % and matidx_dq
0382                 dq_ptr = dq_ptr + 1;
0383                 
0384                 if dq_ptr > lWref_dq
0385                     % Expand the variables by dq_blk_size
0386                     Wref_dq = [Wref_dq; zeros(dq_blk_size,2)];
0387                     matidx_dq = [matidx_dq; zeros(dq_ptr-1 + dq_blk_size,3)];
0388                     lWref_dq = length(Wref_dq(:,1));
0389                 end
0390                 
0391                 Wref_dq(dq_ptr,1:2) = [maxclique{i}(k) maxclique{i}(m)];
0392                 matidx_dq(dq_ptr,1:3) = [i k m+nmaxcliquei];
0393             
0394                 if k ~= m % Already have the diagonal terms of the off-diagonal block
0395 
0396                     dq_ptr = dq_ptr + 1;
0397 
0398                     if dq_ptr > lWref_dq
0399                         % Expand the variables by dq_blk_size
0400                         Wref_dq = [Wref_dq; zeros(dq_blk_size,2)];
0401                         matidx_dq = [matidx_dq; zeros(dq_ptr-1 + dq_blk_size,3)];
0402                     end
0403 
0404                     Wref_dq(dq_ptr,1:2) = [maxclique{i}(m) maxclique{i}(k)];
0405                     matidx_dq(dq_ptr,1:3) = [i k+nmaxcliquei m];
0406                 end
0407             end
0408         end
0409     end
0410     
0411     if verbose >= 2 && mod(i,ndisplay) == 0
0412         fprintf('Loop identification: Loop %i of %i\n',i,nmaxclique);
0413     end
0414 end
0415 
0416 % Trim off excess zeros and empty matrix structures
0417 Wref_dd = Wref_dd(1:dd_ptr,:);
0418 matidx_dd = matidx_dd(1:dd_ptr,:);
0419 
0420 % Store index of Wref variables
0421 Wref_dd = [(1:length(Wref_dd)).' Wref_dd];
0422 Wref_dq = [(1:length(Wref_dq)).' Wref_dq];
0423 
0424 Wref_qq = Wref_dd;
0425 matidx_qq = zeros(size(Wref_qq,1),3);
0426 for i=1:size(Wref_qq,1) 
0427     nmaxcliquei = length(maxclique{matidx_dd(i,1)});
0428     matidx_qq(i,1:3) = [matidx_dd(i,1) matidx_dd(i,2)+nmaxcliquei matidx_dd(i,3)+nmaxcliquei];
0429 end
0430 
0431 
0432 %% Enforce linking constraints in dual
0433 
0434 for i=1:nmaxclique
0435     A{i} = [];
0436 end
0437 
0438 % Count the number of required linking constraints
0439 nbeta = 0;
0440 for i=1:size(E,1)
0441     overlap_idx = intersect(maxclique{E(i,1)},maxclique{E(i,2)});
0442     for k=1:length(overlap_idx)
0443         for m=k:length(overlap_idx)
0444             E11idx = find(maxclique{E(i,1)} == overlap_idx(k));
0445             E12idx = find(maxclique{E(i,1)} == overlap_idx(m));
0446             
0447             nbeta = nbeta + 3;
0448             
0449             if E11idx ~= E12idx
0450                 nbeta = nbeta + 1;
0451             end
0452         end
0453     end
0454     
0455     if verbose >= 2 && mod(i,ndisplay) == 0
0456         fprintf('Counting beta: %i of %i\n',i,size(E,1));
0457     end
0458 end
0459 
0460 % Make beta sdpvar
0461 if nbeta > 0
0462     beta = sdpvar(nbeta,1);
0463 end
0464 
0465 % Create the linking constraints defined by the maximal clique tree
0466 beta_idx = 0;
0467 for i=1:size(E,1)
0468     overlap_idx = intersect(maxclique{E(i,1)},maxclique{E(i,2)});
0469     nmaxclique1 = length(maxclique{E(i,1)});
0470     nmaxclique2 = length(maxclique{E(i,2)});
0471     for k=1:length(overlap_idx)
0472         for m=k:length(overlap_idx)
0473             E11idx = find(maxclique{E(i,1)} == overlap_idx(k));
0474             E12idx = find(maxclique{E(i,1)} == overlap_idx(m));
0475             E21idx = find(maxclique{E(i,2)} == overlap_idx(k));
0476             E22idx = find(maxclique{E(i,2)} == overlap_idx(m));
0477             
0478             beta_idx = beta_idx + 1;
0479             if ~isempty(A{E(i,1)})
0480                 A{E(i,1)} = A{E(i,1)} + 0.5*beta(beta_idx)*sparse([E11idx; E12idx], [E12idx; E11idx], [1; 1], 2*nmaxclique1, 2*nmaxclique1);
0481             else
0482                 A{E(i,1)} = 0.5*beta(beta_idx)*sparse([E11idx; E12idx], [E12idx; E11idx], [1; 1], 2*nmaxclique1, 2*nmaxclique1);
0483             end
0484             if ~isempty(A{E(i,2)})
0485                 A{E(i,2)} = A{E(i,2)} - 0.5*beta(beta_idx)*sparse([E21idx; E22idx], [E22idx; E21idx], [1; 1], 2*nmaxclique2, 2*nmaxclique2);
0486             else
0487                 A{E(i,2)} = -0.5*beta(beta_idx)*sparse([E21idx; E22idx], [E22idx; E21idx], [1; 1], 2*nmaxclique2, 2*nmaxclique2);
0488             end
0489             
0490             beta_idx = beta_idx + 1;
0491             A{E(i,1)} = A{E(i,1)} + 0.5*beta(beta_idx)*sparse([E11idx+nmaxclique1; E12idx+nmaxclique1],[E12idx+nmaxclique1; E11idx+nmaxclique1], [1;1], 2*nmaxclique1, 2*nmaxclique1);
0492             A{E(i,2)} = A{E(i,2)} - 0.5*beta(beta_idx)*sparse([E21idx+nmaxclique2; E22idx+nmaxclique2],[E22idx+nmaxclique2; E21idx+nmaxclique2], [1;1], 2*nmaxclique2, 2*nmaxclique2);
0493             
0494             beta_idx = beta_idx + 1;
0495             A{E(i,1)} = A{E(i,1)} + 0.5*beta(beta_idx)*sparse([E11idx; E12idx+nmaxclique1], [E12idx+nmaxclique1; E11idx], [1;1], 2*nmaxclique1, 2*nmaxclique1);
0496             A{E(i,2)} = A{E(i,2)} - 0.5*beta(beta_idx)*sparse([E21idx; E22idx+nmaxclique2], [E22idx+nmaxclique2; E21idx], [1;1], 2*nmaxclique2, 2*nmaxclique2);
0497             
0498             if E11idx ~= E12idx
0499                 beta_idx = beta_idx + 1;
0500                 A{E(i,1)} = A{E(i,1)} + 0.5*beta(beta_idx)*sparse([E11idx+nmaxclique1; E12idx], [E12idx; E11idx+nmaxclique1], [1;1], 2*nmaxclique1, 2*nmaxclique1);
0501                 A{E(i,2)} = A{E(i,2)} - 0.5*beta(beta_idx)*sparse([E21idx+nmaxclique2; E22idx], [E22idx; E21idx+nmaxclique2], [1;1], 2*nmaxclique2, 2*nmaxclique2);
0502             end
0503         end
0504     end
0505     
0506     if verbose >= 2 && mod(i,ndisplay) == 0
0507         fprintf('SDP linking: %i of %i\n',i,size(E,1));
0508     end
0509     
0510 end
0511 
0512 % For systems with only one maximal clique, A doesn't get defined above.
0513 % Explicitly define it here.
0514 if nmaxclique == 1
0515     A{1} = sparse(2*nbus,2*nbus);
0516 end
0517 
0518 %% Functions to build matrices
0519 
0520 [Yk,Yk_,Mk,Ylineft,Ylinetf,Y_lineft,Y_linetf] = makesdpmat(mpc);
0521 
0522 %% Create SDP variables for dual
0523 
0524 % lambda: active power equality constraints
0525 % gamma:  reactive power equality constraints
0526 % mu:     voltage magnitude constraints
0527 % psi:    generator upper and lower active power limits (created later)
0528 
0529 % Bus variables
0530 nlambda_eq = 0;
0531 ngamma_eq = 0;
0532 ngamma_ineq = 0;
0533 nmu_ineq = 0;
0534 
0535 for i=1:nbus
0536     
0537     nlambda_eq = nlambda_eq + 1;
0538     lambda_eq(i) = nlambda_eq;
0539     
0540     if bustype(i) == PQ || (Qmax(i) - Qmin(i) < minQgendiff)
0541         ngamma_eq = ngamma_eq + 1;
0542         gamma_eq(i) = ngamma_eq;
0543     else
0544         gamma_eq(i) = nan;
0545     end
0546         
0547     if (bustype(i) == PV || bustype(i) == REF) && (Qmax(i) - Qmin(i) >= minQgendiff) && (Qmax(i) <= maxgenlimit || Qmin(i) >= -maxgenlimit)
0548         ngamma_ineq = ngamma_ineq + 1;
0549         gamma_ineq(i) = ngamma_ineq;
0550     else
0551         gamma_ineq(i) = nan;
0552     end
0553     
0554     nmu_ineq = nmu_ineq + 1;
0555     mu_ineq(i) = nmu_ineq;
0556 
0557     psiU_sdp{i} = [];
0558     psiL_sdp{i} = [];
0559     
0560     pwl_sdp{i} = [];
0561 end
0562 
0563 lambdaeq_sdp = sdpvar(nlambda_eq,1);
0564 
0565 gammaeq_sdp = sdpvar(ngamma_eq,1);
0566 gammaU_sdp = sdpvar(ngamma_ineq,1);
0567 gammaL_sdp = sdpvar(ngamma_ineq,1);
0568 
0569 muU_sdp = sdpvar(nmu_ineq,1);
0570 muL_sdp = sdpvar(nmu_ineq,1);
0571 
0572 % Branch flow limit variables
0573 % Hlookup: [line_idx ft]
0574 if upper(mpopt.opf.flow_lim(1)) == 'S'
0575     Hlookup = zeros(2*sum(Smax ~= 0 & Smax < maxlinelimit),2);
0576 end
0577 
0578 nlconstraint = 0;
0579 for i=1:nbranch
0580     if Smax(i) ~= 0 && Smax(i) < maxlinelimit
0581         nlconstraint = nlconstraint + 1;
0582         Hlookup(nlconstraint,:) = [i 1];
0583         nlconstraint = nlconstraint + 1;
0584         Hlookup(nlconstraint,:) = [i 0];
0585     end
0586 end
0587 
0588 if upper(mpopt.opf.flow_lim(1)) == 'S'
0589     for i=1:nlconstraint
0590         Hsdp{i} = sdpvar(3,3);
0591     end
0592 elseif upper(mpopt.opf.flow_lim(1)) == 'P'
0593     Hsdp = sdpvar(2*nlconstraint,1);
0594 end
0595 
0596 
0597 
0598 %% Build bus constraints
0599 
0600 for k=1:nbus
0601         
0602     if bustype(k) == PQ
0603         % PQ bus has equality constraints on P and Q
0604         
0605         % Active power constraint
0606         % -----------------------------------------------------------------
0607         A = addToA(Yk(k), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, -lambdaeq_sdp(lambda_eq(k),1), maxclique);
0608 
0609         % No constraints on lambda_eq_sdp
0610         % -----------------------------------------------------------------
0611         
0612         
0613         % Reactive power constraint
0614         % -----------------------------------------------------------------
0615         A = addToA(Yk_(k), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, -gammaeq_sdp(gamma_eq(k),1), maxclique);
0616         
0617         % No constraints on gamma_eq_sdp
0618         % -----------------------------------------------------------------
0619         
0620         
0621         % Create cost function
0622         % -----------------------------------------------------------------
0623         if exist('cost','var')
0624             cost = cost ... 
0625                 + lambdaeq_sdp(lambda_eq(k))*(-Pd(k)) ...
0626                 + gammaeq_sdp(gamma_eq(k))*(-Qd(k));
0627         else
0628             cost = lambdaeq_sdp(lambda_eq(k))*(-Pd(k)) ...
0629                    + gammaeq_sdp(gamma_eq(k))*(-Qd(k));
0630         end
0631         % -----------------------------------------------------------------
0632         
0633     elseif bustype(k) == PV || bustype(k) == REF 
0634         % PV and slack buses have upper and lower constraints on both P and
0635         % Q, unless the constraints are tighter than minPgendiff and
0636         % minQgendiff, in which case the inequalities are set to equalities
0637         % at the midpoint of the upper and lower limits.
0638         
0639         genidx = find(mpc.gen(:,GEN_BUS) == k);
0640         
0641         if isempty(genidx)
0642             error('sdpopf_solver: Generator for bus %i not found!',k);
0643         end
0644         
0645         % Active power constraints
0646         % -----------------------------------------------------------------
0647         % Make Lagrange multipliers for generator limits
0648         psiU_sdp{k} = sdpvar(length(genidx),1);
0649         psiL_sdp{k} = sdpvar(length(genidx),1);
0650         
0651         clear lambdai Pmax Pmin
0652         for i=1:length(genidx)
0653             
0654             % Does this generator have pwl costs?
0655             if mpc.gencost(genidx(i),MODEL) == PW_LINEAR
0656                 pwl = 1;
0657                 % If so, define Lagrange multipliers for each segment
0658                 pwl_sdp{k}{i} = sdpvar(mpc.gencost(genidx(i),NCOST)-1,1);
0659                 
0660                 % Get breakpoints for this curve
0661                 x_pw = mpc.gencost(genidx(i),-1+COST+(1:2:2*length(pwl_sdp{k}{i})+1)) / Sbase; % piecewise powers
0662                 c_pw = mpc.gencost(genidx(i),-1+COST+(2:2:2*length(pwl_sdp{k}{i})+2)); % piecewise costs
0663             else
0664                 pwl = 0;
0665                 pwl_sdp{k}{i} = [];
0666             end
0667             
0668             Pmax(i) = mpc.gen(genidx(i),PMAX) / Sbase;
0669             Pmin(i) = mpc.gen(genidx(i),PMIN) / Sbase;
0670             
0671             if (Pmax(i) - Pmin(i) >= minPgendiff)
0672                 
0673                 % Constrain upper and lower Lagrange multipliers to be positive
0674                 constraints = [constraints;
0675                     psiL_sdp{k}(i) >= 0;
0676                     psiU_sdp{k}(i) >= 0];
0677                 
0678                 if pwl
0679                     
0680                     lambdasum = 0;
0681                     for pwiter = 1:length(pwl_sdp{k}{i})
0682                         
0683                         m_pw = (c_pw(pwiter+1) - c_pw(pwiter)) / (x_pw(pwiter+1) - x_pw(pwiter)); % slope for this piecewise cost segement
0684                         
0685                         if exist('cost','var')
0686                             cost = cost - pwl_sdp{k}{i}(pwiter) * (m_pw*x_pw(pwiter) - c_pw(pwiter));
0687                         else
0688                             cost = -pwl_sdp{k}{i}(pwiter) * (m_pw*x_pw(pwiter) - c_pw(pwiter));
0689                         end
0690                         
0691                         lambdasum = lambdasum + pwl_sdp{k}{i}(pwiter) * m_pw;
0692 
0693                     end
0694                     
0695                     constraints = [constraints; 
0696                                    pwl_sdp{k}{i}(:) >= 0];
0697                     
0698                     if ~exist('lambdai','var')
0699                         lambdai = lambdasum + psiU_sdp{k}(i) - psiL_sdp{k}(i);
0700                     else
0701                         % Equality constraint on prices (equal marginal
0702                         % prices requirement) for all generators at a bus
0703 %                         constraints = [constraints;
0704 %                             lambdasum + psiU_sdp{k}(i) - psiL_sdp{k}(i) == lambdai];
0705                         
0706                         % Instead of a strict equality, numerical
0707                         % performance works better with tight inequalities
0708                         % for equal marginal price requirement with
0709                         % multiple piecewise generators at the same bus.
0710                         constraints = [constraints;
0711                             lambdasum + psiU_sdp{k}(i) - psiL_sdp{k}(i) - binding_lagrange <= lambdai <= lambdasum + psiU_sdp{k}(i) - psiL_sdp{k}(i) + binding_lagrange];
0712                     end
0713                     
0714                     constraints = [constraints; 
0715                                     sum(pwl_sdp{k}{i}) == 1];
0716                     
0717                 else
0718                     if c2(genidx(i)) > 0
0719                         % Cost function has nonzero quadratic term
0720 
0721                         % Make a new R matrix
0722                         R{k}{i} = sdpvar(2,1);
0723 
0724                         constraints = [constraints; 
0725                             [1 R{k}{i}(1);
0726                              R{k}{i}(1) R{k}{i}(2)] >= 0]; % R psd
0727 
0728                         if ~exist('lambdai','var')
0729                             lambdai = c1(genidx(i)) + 2*sqrt(c2(genidx(i)))*R{k}{i}(1) + psiU_sdp{k}(i) - psiL_sdp{k}(i);
0730                         else
0731                             % Equality constraint on prices (equal marginal
0732                             % prices requirement) for all generators at a bus
0733                             constraints = [constraints;
0734                                 c1(genidx(i)) + 2*sqrt(c2(genidx(i)))*R{k}{i}(1) + psiU_sdp{k}(i) - psiL_sdp{k}(i) == lambdai];
0735                         end
0736 
0737                         if exist('cost','var')
0738                             cost = cost - R{k}{i}(2);
0739                         else
0740                             cost = -R{k}{i}(2);
0741                         end
0742                     else
0743                         % Cost function is linear for this generator
0744                         R{k}{i} = zeros(2,1);
0745 
0746                         if ~exist('lambdai','var')
0747                             lambdai = c1(genidx(i)) + psiU_sdp{k}(i) - psiL_sdp{k}(i);
0748                         else
0749                             % Equality constraint on prices (equal marginal
0750                             % prices requirement) for all generators at a bus
0751                             constraints = [constraints;
0752                                 c1(genidx(i)) + psiU_sdp{k}(i) - psiL_sdp{k}(i) == lambdai];
0753                         end
0754 
0755                     end
0756                 end
0757                 
0758                 if exist('cost','var')
0759                     cost = cost ...
0760                            + psiL_sdp{k}(i)*Pmin(i) - psiU_sdp{k}(i)*Pmax(i);
0761                 else
0762                     cost = psiL_sdp{k}(i)*Pmin(i) - psiU_sdp{k}(i)*Pmax(i);
0763                 end
0764 
0765             else % Active power generator is constrained to a small range
0766                 % Set the power injection to the midpoint of the small range
0767                 % allowed by the problem, and use a free variable to force an
0768                 % equality constraint.
0769 
0770                 Pavg = (Pmin(i)+Pmax(i))/2;
0771                 Pd(k) = Pd(k) - Pavg;
0772 
0773                 % Add on the (known) cost of this generation
0774                 if pwl
0775                     % Figure out what this fixed generation costs
0776                     pwidx = find(Pavg <= x_pw,1);
0777                     if isempty(pwidx)
0778                         pwidx = length(x_pw);
0779                     end
0780                     pwidx = pwidx - 1;
0781                     
0782                     m_pw = (c_pw(pwidx+1) - c_pw(pwidx)) / (x_pw(pwidx+1) - x_pw(pwidx)); % slope for this piecewise cost segement
0783                     fix_gen_cost = m_pw * (Pavg - x_pw(pwidx)) + c_pw(pwidx);
0784                     
0785                     if exist('cost','var')
0786                         cost = cost + fix_gen_cost;
0787                     else
0788                         cost = fix_gen_cost;
0789                     end
0790                 else
0791                     c0(genidx(i)) = c2(genidx(i))*Pavg^2 + c1(genidx(i))*Pavg + c0(genidx(i));
0792                 end
0793                 
0794                 % No constraints on lambdaeq_sdp (free var)
0795             end
0796             
0797             if ~pwl
0798                 if exist('cost','var')
0799                     cost = cost + c0(genidx(i));
0800                 else
0801                     cost = c0(genidx(i));
0802                 end
0803             end
0804         end
0805                
0806         % If there are only generator(s) with tight active power
0807         % constraints at this bus, explicitly set an equality constraint to
0808         % the (modified) load demand directly.
0809         if exist('lambdai','var')
0810             lambda_aggregate = lambdai;
0811         else
0812             lambda_aggregate = lambdaeq_sdp(lambda_eq(k));
0813         end
0814         
0815         lambdaeq_sdp(k) = -lambda_aggregate;
0816         
0817         cost = cost + lambda_aggregate*Pd(k);
0818         
0819         A = addToA(Yk(k), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, lambda_aggregate, maxclique);
0820         % -----------------------------------------------------------------
0821         
0822 
0823         % Reactive power constraints
0824         % -----------------------------------------------------------------
0825         % Sum reactive power injections from all generators at each bus
0826         
0827         clear gamma_aggregate
0828         if (Qmax(k) - Qmin(k) >= minQgendiff)
0829             
0830             if Qmin(k) >= -maxgenlimit % This generator has a lower Q limit
0831                 
0832                 gamma_aggregate = -gammaL_sdp(gamma_ineq(k));
0833                 
0834                 if exist('cost','var')
0835                     cost = cost ...
0836                            + gammaL_sdp(gamma_ineq(k))*(Qmin(k) - Qd(k));
0837                 else
0838                     cost = gammaL_sdp(gamma_ineq(k))*(Qmin(k) - Qd(k));
0839                 end
0840             
0841                 % Constrain upper and lower lagrange multipliers to be positive
0842                 constraints = [constraints;
0843                     gammaL_sdp(gamma_ineq(k)) >= 0];        
0844             end
0845             
0846             if Qmax(k) <= maxgenlimit % We have an upper Q limit
0847                 
0848                 if exist('gamma_aggregate','var')
0849                     gamma_aggregate = gamma_aggregate + gammaU_sdp(gamma_ineq(k));
0850                 else
0851                     gamma_aggregate = gammaU_sdp(gamma_ineq(k));
0852                 end
0853                 
0854                 if exist('cost','var')
0855                     cost = cost ...
0856                            - gammaU_sdp(gamma_ineq(k))*(Qmax(k) - Qd(k));
0857                 else
0858                     cost = -gammaU_sdp(gamma_ineq(k))*(Qmax(k) - Qd(k));
0859                 end
0860             
0861                 % Constrain upper and lower lagrange multipliers to be positive
0862                 constraints = [constraints;
0863                     gammaU_sdp(gamma_ineq(k)) >= 0];    
0864             end
0865             
0866         else
0867             gamma_aggregate = -gammaeq_sdp(gamma_eq(k));
0868             
0869             if exist('cost','var')
0870                 cost = cost ...
0871                        + gammaeq_sdp(gamma_eq(k))*((Qmin(k)+Qmax(k))/2 - Qd(k));
0872             else
0873                 cost = gammaeq_sdp(gamma_eq(k))*((Qmin(k)+Qmax(k))/2 - Qd(k));
0874             end
0875             
0876             % No constraints on gammaeq_sdp (free var)
0877         end
0878         
0879         % Add to appropriate A matrices
0880         if (Qmax(k) <= maxgenlimit || Qmin(k) >= -maxgenlimit)
0881             A = addToA(Yk_(k), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, gamma_aggregate, maxclique);
0882         end
0883         % -----------------------------------------------------------------
0884         
0885     else
0886         error('sdpopf_solver: Invalid bus type');
0887     end
0888 
0889     
0890     % Voltage magnitude constraints
0891     % -----------------------------------------------------------------
0892     mu_aggregate = -muL_sdp(mu_ineq(k)) + muU_sdp(mu_ineq(k));    
0893     A = addToA(Mk(k), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, mu_aggregate, maxclique);
0894     
0895     % Constrain upper and lower Lagrange multipliers to be positive
0896     constraints = [constraints;
0897         muL_sdp(mu_ineq(k)) >= 0;
0898         muU_sdp(mu_ineq(k)) >= 0];
0899     
0900     % Create sdp cost function
0901     cost = cost + muL_sdp(mu_ineq(k))*Vmin(k)^2 - muU_sdp(mu_ineq(k))*Vmax(k)^2;
0902     % -----------------------------------------------------------------
0903     
0904     if verbose >= 2 && mod(k,ndisplay) == 0
0905         fprintf('SDP creation: Bus %i of %i\n',k,nbus);
0906     end
0907     
0908 end % Loop through all buses
0909 
0910 
0911 
0912 %% Build line constraints
0913 
0914 nlconstraint = 0;
0915 
0916 for i=1:nbranch
0917     if Smax(i) ~= 0 && Smax(i) < maxlinelimit
0918         if upper(mpopt.opf.flow_lim(1)) == 'S'  % Constrain MVA at both line terminals
0919             
0920             nlconstraint = nlconstraint + 1;
0921             cost = cost ...
0922                    - (Smax(i)^2*Hsdp{nlconstraint}(1,1) + Hsdp{nlconstraint}(2,2) + Hsdp{nlconstraint}(3,3));
0923                
0924             constraints = [constraints; Hsdp{nlconstraint} >= 0];
0925             
0926             A = addToA(Ylineft(i), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, 2*Hsdp{nlconstraint}(1,2), maxclique);
0927             A = addToA(Y_lineft(i), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, 2*Hsdp{nlconstraint}(1,3), maxclique);
0928             
0929             nlconstraint = nlconstraint + 1;
0930             cost = cost ...
0931                    - (Smax(i)^2*Hsdp{nlconstraint}(1,1) + Hsdp{nlconstraint}(2,2) + Hsdp{nlconstraint}(3,3));
0932                
0933             constraints = [constraints; Hsdp{nlconstraint} >= 0];
0934             
0935             A = addToA(Ylinetf(i), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, 2*Hsdp{nlconstraint}(1,2), maxclique);
0936             A = addToA(Y_linetf(i), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, 2*Hsdp{nlconstraint}(1,3), maxclique);
0937 
0938         elseif upper(mpopt.opf.flow_lim(1)) == 'P'  % Constrain active power flow at both terminals
0939             
0940             nlconstraint = nlconstraint + 1;
0941             cost = cost ...
0942                    - Smax(i)*Hsdp(nlconstraint);
0943 
0944             A = addToA(Ylineft(i), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, Hsdp(nlconstraint), maxclique);
0945             
0946             nlconstraint = nlconstraint + 1;
0947             cost = cost ...
0948                    - Smax(i)*Hsdp(nlconstraint);
0949             
0950             A = addToA(Ylinetf(i), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, Hsdp(nlconstraint), maxclique);
0951             
0952         else
0953             error('sdpopf_solver: Invalid line constraint option');
0954         end
0955     end
0956     
0957     if verbose >= 2 && mod(i,ndisplay) == 0
0958         fprintf('SDP creation: Branch %i of %i\n',i,nbranch);
0959     end
0960     
0961 end
0962 
0963 if upper(mpopt.opf.flow_lim(1)) == 'P'
0964     constraints = [constraints; Hsdp >= 0];
0965 end
0966 
0967 
0968 %% Formulate dual psd constraints
0969 
0970 Aconstraints = zeros(nmaxclique,1);
0971 for i=1:nmaxclique
0972     % Can multiply A by any non-zero scalar. This may affect numerics
0973     % of the solver.
0974     constraints = [constraints; 1*A{i} >= 0]; 
0975     Aconstraints(i) = length(constraints);
0976 end
0977 
0978 tsdpform = toc;
0979 
0980 %% Solve the SDP
0981 
0982 if recover_voltage == 2 || recover_voltage == 3 || recover_voltage == 4
0983     sdpopts = sdpsettings(sdpopts,'saveduals',1);
0984 end
0985 
0986 % Preserve warning settings
0987 S = warning;
0988 
0989 % sdpopts = sdpsettings(sdpopts,'sedumi.eps',0);
0990 
0991 % Run sdp solver
0992 sdpinfo = solvesdp(constraints, -cost, sdpopts); % Negative cost to convert maximization to minimization problem
0993 
0994 if sdpinfo.problem == 2 || sdpinfo.problem == -2 || sdpinfo.problem == -3
0995     error(yalmiperror(sdpinfo.problem));
0996 end
0997 if ~have_feature('octave') || have_feature('octave', 'vnum') >= 4.001
0998     %% (avoid bug in Octave 4.0.x, where warning state is left corrupted)
0999     warning(S);
1000 end
1001 tsolve = sdpinfo.solvertime;
1002 
1003 if verbose >= 2
1004     fprintf('Solver exit message: %s\n',sdpinfo.info);
1005 end
1006 
1007 objval = double(cost);
1008 
1009 %% Extract the optimal voltage profile from the SDP solution
1010 
1011 % Find the eigenvector associated with a zero eigenvalue of each A matrix.
1012 % Enforce consistency between terms in different nullspace eigenvectors
1013 % that refer to the same voltage component. Each eigenvector can be
1014 % multiplied by a complex scalar to enforce consistency: create a matrix
1015 % with elements of the eigenvectors. Two more binding constraints are
1016 % needed: one with a voltage magnitude at its upper or lower limit and one
1017 % with the voltage angle equal to zero at the reference bus.
1018 
1019 tic;
1020 
1021 % Calculate eigenvalues from dual problem (A matrices)
1022 if recover_voltage == 1 || recover_voltage == 3 || recover_voltage == 4
1023     dual_mineigratio = inf;
1024     dual_eval = nan(length(A),1);
1025     for i=1:length(A)
1026         [evc, evl] = eig(double(A{i}));
1027 
1028         % We only need the eigenvector associated with a zero eigenvalue
1029         evl = diag(evl);
1030         dual_u{i} = evc(:,abs(evl) == min(abs(evl)));
1031 
1032         % Calculate mineigratio
1033         dual_eigA_ratio = abs(evl(3) / evl(1));
1034         dual_eval(i) = 1/dual_eigA_ratio;
1035         
1036         if dual_eigA_ratio < dual_mineigratio
1037             dual_mineigratio = dual_eigA_ratio;
1038         end
1039     end
1040 end
1041 
1042 % Calculate eigenvalues from primal problem (W matrices)
1043 if recover_voltage == 2 || recover_voltage == 3  || recover_voltage == 4
1044     primal_mineigratio = inf;
1045     primal_eval = nan(length(A),1);
1046     for i=1:length(A)
1047         [evc, evl] = eig(double(dual(constraints(Aconstraints(i)))));
1048 
1049         % We only need the eigenvector associated with a non-zero eigenvalue
1050         evl = diag(evl);
1051         primal_u{i} = evc(:,abs(evl) == max(abs(evl)));
1052 
1053         % Calculate mineigratio
1054         primal_eigA_ratio = abs(evl(end) / evl(end-2));
1055         primal_eval(i) = 1/primal_eigA_ratio;
1056         
1057         if primal_eigA_ratio < primal_mineigratio
1058             primal_mineigratio = primal_eigA_ratio;
1059         end
1060     end
1061 end
1062 
1063 if recover_voltage == 3 || recover_voltage == 4
1064     recover_voltage_loop = [1 2];
1065 else
1066     recover_voltage_loop = recover_voltage;
1067 end
1068     
1069 for r = 1:length(recover_voltage_loop)
1070     
1071     if recover_voltage_loop(r) == 1
1072         if verbose >= 2
1073             fprintf('Recovering dual solution.\n');
1074         end
1075         u = dual_u;
1076         eval = dual_eval;
1077         mineigratio = dual_mineigratio;
1078     elseif recover_voltage_loop(r) == 2
1079         if verbose >= 2
1080             fprintf('Recovering primal solution.\n');
1081         end
1082         u = primal_u;
1083         eval = primal_eval;
1084         mineigratio = primal_mineigratio;
1085     else
1086         error('sdpopf_solver: Invalid recover_voltage');
1087     end
1088  
1089     % Create a linear equation to match up terms from different matrices
1090     % Each matrix gets two unknown parameters: alpha and beta.
1091     % The unknown vector is arranged as [alpha1; alpha2; ... ; beta1; beta2 ; ...]
1092     % corresponding to the matrices representing maxclique{1}, maxclique{2}, etc.
1093     Aint = [];
1094     Aslack = [];
1095     slackbus_idx = find(mpc.bus(:,2) == 3);
1096     for i=1:nmaxclique
1097         maxcliquei = maxclique{i};
1098         for k=i+1:nmaxclique
1099             
1100             temp = maxcliquei(ismembc_(maxcliquei, maxclique{k}));
1101             
1102             % Slower alternative that does not use undocumented MATLAB
1103             % functions:
1104             % temp = intersect(maxclique{i},maxclique{k});
1105             
1106             Aint = [Aint; i*ones(length(temp),1) k*ones(length(temp),1) temp(:)];
1107         end
1108 
1109         if any(maxcliquei == slackbus_idx)
1110             Aslack = i;
1111         end
1112     end
1113     if nmaxclique > 1
1114         L = sparse(2*length(Aint(:,1))+1,2*nmaxclique);
1115     else
1116         L = sparse(1,2);
1117     end
1118     Lrow = 0; % index the equality
1119     for i=1:size(Aint,1)
1120         % first the real part of the voltage
1121         Lrow = Lrow + 1;
1122         L(Lrow,Aint(i,1)) = u{Aint(i,1)}(maxclique{Aint(i,1)}(:) == Aint(i,3));
1123         L(Lrow,Aint(i,1)+nmaxclique) = -u{Aint(i,1)}(find(maxclique{Aint(i,1)}(:) == Aint(i,3))+length(maxclique{Aint(i,1)}(:)));
1124         L(Lrow,Aint(i,2)) = -u{Aint(i,2)}(maxclique{Aint(i,2)}(:) == Aint(i,3));
1125         L(Lrow,Aint(i,2)+nmaxclique) = u{Aint(i,2)}(find(maxclique{Aint(i,2)}(:) == Aint(i,3))+length(maxclique{Aint(i,2)}(:)));
1126 
1127         % then the imaginary part of the voltage
1128         Lrow = Lrow + 1;
1129         L(Lrow,Aint(i,1)) = u{Aint(i,1)}(find(maxclique{Aint(i,1)}(:) == Aint(i,3))+length(maxclique{Aint(i,1)}(:)));
1130         L(Lrow,Aint(i,1)+nmaxclique) = u{Aint(i,1)}(maxclique{Aint(i,1)}(:) == Aint(i,3));
1131         L(Lrow,Aint(i,2)) = -u{Aint(i,2)}(find(maxclique{Aint(i,2)}(:) == Aint(i,3))+length(maxclique{Aint(i,2)}(:)));
1132         L(Lrow,Aint(i,2)+nmaxclique) = -u{Aint(i,2)}(maxclique{Aint(i,2)}(:) == Aint(i,3));
1133     end
1134 
1135     % Add in angle reference constraint, which is linear
1136     Lrow = Lrow + 1;
1137     L(Lrow,Aslack) = u{Aslack}(find(maxclique{Aslack}(:) == slackbus_idx)+length(maxclique{Aslack}(:)));
1138     L(Lrow,Aslack+nmaxclique) = u{Aslack}(maxclique{Aslack}(:) == slackbus_idx);
1139 
1140     % Get a vector in the nullspace
1141     if nmaxclique == 1 % Will probably get an eigenvalue exactly at zero. Only a 2x2 matrix, just do a normal eig
1142         [LLevec,LLeval] = eig(full(L).'*full(L));
1143         LLeval = LLeval(1);
1144         LLevec = LLevec(:,1);
1145     else
1146         % sometimes get an error when the solution is "too good" and we have an
1147         % eigenvalue very close to zero. eigs can't handle this situation
1148         % eigs is necessary for large numbers of loops, which shouldn't
1149         % have exactly a zero eigenvalue. For small numbers of loops,
1150         % we can just do eig. First try eigs, if it fails, then try eig.
1151         try
1152             warning off
1153             [LLevec,LLeval] = eigs(L.'*L,1,'SM');
1154             if ~have_feature('octave') || have_feature('octave', 'vnum') >= 4.001
1155                 %% (avoid bug in Octave 4.0.x, where warning state is left corrupted)
1156                 warning(S);
1157             end
1158         catch
1159             [LLevec,LLeval] = eig(full(L).'*full(L));
1160             LLeval = LLeval(1);
1161             LLevec = LLevec(:,1);
1162             if ~have_feature('octave') || have_feature('octave', 'vnum') >= 4.001
1163                 %% (avoid bug in Octave 4.0.x, where warning state is left corrupted)
1164                 warning(S);
1165             end
1166         end
1167     end
1168 
1169     % Form a voltage vector from LLevec. Piece together the u{i} to form
1170     % a big vector U. This vector has one degree of freedom in its length.
1171     % When we get values that are inconsistent (due to numerical errors if
1172     % the SDP OPF satisifes the rank condition), take an average of the
1173     % values, weighted by the corresponding eigenvalue.
1174     U = sparse(2*nbus,nmaxclique);
1175     Ueval = sparse(2*nbus,nmaxclique); 
1176     for i=1:nmaxclique
1177         newentries = [maxclique{i}(:); maxclique{i}(:)+nbus];
1178         newvals_phasor = (u{i}(1:length(maxclique{i}))+1i*u{i}(length(maxclique{i})+1:2*length(maxclique{i})))*(LLevec(i)+1i*LLevec(i+nmaxclique));
1179         newvals = [real(newvals_phasor(:)); imag(newvals_phasor(:))];
1180 
1181         Ueval(newentries,i) = 1/eval(i);
1182         U(newentries,i) = newvals;
1183 
1184     end
1185     Ueval_sum = sum(Ueval,2);
1186     for i=1:2*nbus
1187         Ueval(i,:) = Ueval(i,:) / Ueval_sum(i);
1188     end
1189     U = full(sum(U .* Ueval,2));
1190 
1191     % Find the binding voltage constraint with the largest Lagrange multiplier
1192     muL_val = abs(double(muL_sdp(mu_ineq)));
1193     muU_val = abs(double(muU_sdp(mu_ineq)));
1194 
1195     [maxmuU,binding_voltage_busU] = max(muU_val);
1196     [maxmuL,binding_voltage_busL] = max(muL_val);
1197 
1198     if (maxmuU < binding_lagrange) && (maxmuL < binding_lagrange)
1199         if verbose > 1
1200             fprintf('No binding voltage magnitude. Attempting solution recovery from line-flow constraints.\n');
1201         end
1202 
1203         % No voltage magnitude limits are binding.
1204         % Try to recover the solution from binding line flow limits.
1205         % This should cover the vast majority of cases
1206 
1207         % First try to find a binding line flow limit.
1208         largest_H_norm = 0;
1209         largest_H_norm_idx = 0;
1210         for i=1:length(Hsdp)
1211             if norm(double(Hsdp{i})) > largest_H_norm
1212                 largest_H_norm = norm(double(Hsdp{i}));
1213                 largest_H_norm_idx = i;
1214             end
1215         end
1216 
1217         if largest_H_norm < binding_lagrange
1218             if verbose > 1
1219                 error('sdpopf_solver: No binding voltage magnitude or line-flow limits, and no PQ buses. Error recovering solution.');
1220             end
1221         else
1222             % Try to recover solution from binding line-flow limit
1223             if Hlookup(largest_H_norm_idx,2) == 1
1224                 [Yline_row, Yline_col, Yline_val] = find(Ylineft(Hlookup(largest_H_norm_idx,1)));
1225                 [Y_line_row, Y_line_col, Y_line_val] = find(Y_lineft(Hlookup(largest_H_norm_idx,1)));
1226             else
1227                 [Yline_row, Yline_col, Yline_val] = find(Ylinetf(Hlookup(largest_H_norm_idx,1)));
1228                 [Y_line_row, Y_line_col, Y_line_val] = find(Y_linetf(Hlookup(largest_H_norm_idx,1)));
1229             end
1230             trYlineUUT = 0;
1231             for i=1:length(Yline_val)
1232                 trYlineUUT = trYlineUUT + Yline_val(i) * U(Yline_row(i))*U(Yline_col(i));
1233             end
1234             
1235             if upper(mpopt.opf.flow_lim(1)) == 'S'
1236                 trY_lineUUT = 0;
1237                 for i=1:length(Y_line_val)
1238                     trY_lineUUT = trY_lineUUT + Y_line_val(i) * U(Y_line_row(i))*U(Y_line_col(i));
1239                 end
1240                 chi = (Smax(Hlookup(largest_H_norm_idx,1))^2 / ( trYlineUUT^2 + trY_lineUUT^2 ) )^(0.25);
1241             elseif upper(mpopt.opf.flow_lim(1)) == 'P'
1242                 % We want to scale U so that trace(chi^2*U*U.') = Smax
1243                 chi = sqrt(Smax(Hlookup(largest_H_norm_idx,1)) / trYlineUUT);
1244             end
1245             
1246             U = chi*U;
1247         end
1248         
1249     else % recover from binding voltage magnitude
1250         if maxmuU > maxmuL
1251             binding_voltage_bus = binding_voltage_busU;
1252             binding_voltage_mag = Vmax(binding_voltage_bus);
1253         else
1254             binding_voltage_bus = binding_voltage_busL;
1255             binding_voltage_mag = Vmin(binding_voltage_bus);
1256         end
1257         chi = binding_voltage_mag^2 / (U(binding_voltage_bus)^2 + U(binding_voltage_bus+nbus)^2);
1258         U = sqrt(chi)*U;
1259     end
1260 
1261     Vopt = (U(1:nbus) + 1i*U(nbus+1:2*nbus));
1262 
1263     if real(Vopt(slackbus_idx)) < 0
1264         Vopt = -Vopt;
1265     end
1266     
1267     if recover_voltage == 3 || recover_voltage == 4
1268         Sopt = Vopt .* conj(Y*Vopt);
1269         PQerr = norm([real(Sopt(bustype(:) == PQ))*100+mpc.bus(bustype(:) == PQ,PD); imag(Sopt(bustype(:) == PQ))*100+mpc.bus(bustype(:) == PQ,QD)], 2);
1270         if recover_voltage_loop(r) == 1
1271             Vopt_dual = Vopt;
1272             dual_PQerr = PQerr;
1273         elseif recover_voltage_loop(r) == 2
1274             Vopt_primal = Vopt;
1275             primal_PQerr = PQerr;
1276         end
1277     end
1278 end
1279 
1280 if recover_voltage == 3
1281     if dual_PQerr > primal_PQerr
1282         fprintf('dual_PQerr: %g MW. primal_PQerr: %g MW. Using primal solution.\n',dual_PQerr,primal_PQerr)
1283         Vopt = Vopt_primal;
1284     else
1285         fprintf('dual_PQerr: %g MW. primal_PQerr: %g MW. Using dual solution.\n',dual_PQerr,primal_PQerr)
1286         Vopt = Vopt_dual;
1287     end
1288 elseif recover_voltage == 4
1289     if primal_mineigratio > dual_mineigratio
1290         if verbose >= 2
1291             fprintf('dual_mineigratio: %g. primal_mineigratio: %g. Using primal solution.\n',dual_mineigratio,primal_mineigratio)
1292         end
1293         Vopt = Vopt_primal;
1294     else
1295         if verbose >= 2
1296             fprintf('dual_mineigratio: %g. primal_mineigratio: %g. Using dual solution.\n',dual_mineigratio,primal_mineigratio)
1297         end
1298         Vopt = Vopt_dual;
1299     end
1300 end
1301 
1302 % Store voltages
1303 for i=1:nbus
1304     mpc.bus(i,VM) = abs(Vopt(i));
1305     mpc.bus(i,VA) = angle(Vopt(i))*180/pi;
1306     mpc.gen(mpc.gen(:,GEN_BUS) == i,VG) = abs(Vopt(i));
1307 end
1308 
1309 
1310 %% Calculate power injections
1311 
1312 if recover_injections == 2
1313     % Calculate directly from the SDP solution
1314     
1315     for i=1:length(A)
1316         W{i} = dual(constraints(Aconstraints(i)));
1317     end
1318 
1319     Pinj = zeros(nbus,1);
1320     Qinj = zeros(nbus,1);
1321     for i=1:nbus
1322 
1323         % Active power injection
1324         Pinj(i) = recoverFromW(Yk(i), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, W, maxclique);
1325 
1326         % Reactive power injection
1327         Qinj(i) = recoverFromW(Yk_(i), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, W, maxclique);
1328 
1329         % Voltage magnitude
1330 %         Vmag(i)  = sqrt(recoverFromW(Mk(i), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, W, maxclique));
1331 
1332     end
1333     
1334     Pflowft = zeros(nbranch,1);
1335     Pflowtf = zeros(nbranch,1);
1336     Qflowft = zeros(nbranch,1);
1337     Qflowtf = zeros(nbranch,1);
1338     for i=1:nbranch
1339         Pflowft(i) = recoverFromW(Ylineft(i), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, W, maxclique) * Sbase;
1340         Pflowtf(i) = recoverFromW(Ylinetf(i), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, W, maxclique) * Sbase;
1341         Qflowft(i) = recoverFromW(Y_lineft(i), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, W, maxclique) * Sbase;
1342         Qflowtf(i) = recoverFromW(Y_linetf(i), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, W, maxclique) * Sbase;
1343     end
1344     
1345     
1346 elseif recover_injections == 1 
1347     % Calculate power injections from voltage profile
1348     Sopt = Vopt .* conj(Y*Vopt);
1349     Pinj = real(Sopt);
1350     Qinj = imag(Sopt);
1351    
1352     Vf = Vopt(mpc.branch(:,F_BUS));
1353     Vt = Vopt(mpc.branch(:,T_BUS));
1354     
1355     Sft = Vf.*conj(Yf*Vopt);
1356     Stf = Vt.*conj(Yt*Vopt);
1357     
1358     Pflowft = real(Sft) * Sbase;
1359     Pflowtf = real(Stf) * Sbase;
1360     Qflowft = imag(Sft) * Sbase;
1361     Qflowtf = imag(Stf) * Sbase;
1362 end
1363 
1364 % Store the original bus ordering
1365 busorder = mpc.bus(:,end);
1366 mpc.bus = mpc.bus(:,1:end-1);
1367     
1368 % Convert from injections to generation
1369 Qtol = 5e-2;
1370 mpc.gen(:,PG) = 0;
1371 mpc.gen(:,QG) = 0;
1372 for i=1:nbus
1373     
1374     if bustype(i) == PV || bustype(i) == REF
1375         % Find all generators at this bus
1376         genidx = find(mpc.gen(:,GEN_BUS) == i);
1377 
1378         % Find generator outputs from the known injections by solving an
1379         % economic dispatch problem at each bus. (Identifying generators at
1380         % their limits from non-zero dual variables is not numerically
1381         % reliable.)
1382         
1383         % Initialize all generators at this bus
1384         gen_not_limited_pw = [];
1385         gen_not_limited_q = [];
1386         Pmax = mpc.gen(genidx,PMAX);
1387         Pmin = mpc.gen(genidx,PMIN);
1388         remove_from_list = [];
1389         for k=1:length(genidx)
1390             % Check if this generator is fixed at its midpoint
1391             if (Pmax(k) - Pmin(k) < minPgendiff * Sbase)
1392                 mpc.gen(genidx(k),PG) = (Pmin(k)+Pmax(k))/2;
1393                 remove_from_list = [remove_from_list; k];
1394             else
1395                 if mpc.gencost(genidx(k),MODEL) == PW_LINEAR
1396                     gen_not_limited_pw = [gen_not_limited_pw; genidx(k)];
1397                     mpc.gen(genidx(k),PG) = mpc.gen(genidx(k),PMIN);
1398                 else
1399                     gen_not_limited_q = [gen_not_limited_q; genidx(k)];
1400                     mpc.gen(genidx(k),PG) = 0;
1401                 end
1402             end
1403         end
1404         Pmax(remove_from_list) = [];
1405         Pmin(remove_from_list) = [];
1406 
1407         % Figure out the remaining power injections to distribute among the
1408         % generators at this bus.
1409         Pgen_remain = Pinj(i) + mpc.bus(i,PD) / Sbase - sum(mpc.gen(genidx,PG)) / Sbase;
1410 
1411         % Use the equal marginal cost criterion to distribute Pgen_remain
1412         gen_not_limited = [gen_not_limited_q; gen_not_limited_pw];
1413         if length(gen_not_limited) == 1
1414             % If only one generator not at its limit, assign the remaining
1415             % generation to this generator
1416             mpc.gen(gen_not_limited,PG) = mpc.gen(gen_not_limited,PG) + Pgen_remain * Sbase;
1417         elseif length(genidx) == 1
1418             % If there is only one generator at this bus, assign the
1419             % remaining generation to this generator.
1420             mpc.gen(genidx,PG) = mpc.gen(genidx,PG) + Pgen_remain * Sbase;
1421         else % Multiple generators that aren't at their limits
1422             
1423             % We know the LMP at this bus. For that LMP, determine the
1424             % generation output for each generator.
1425             
1426             % For piecewise linear cost functions:
1427             % Stack Pgen_remain into the generators in the order of
1428             % cheapest cost.
1429             if ~isempty(gen_not_limited_pw)
1430                 while Pgen_remain > 1e-5
1431 
1432                     % Find cheapest available segment
1433                     lc = inf;
1434                     for geniter = 1:length(gen_not_limited_pw)
1435                         nseg = mpc.gencost(gen_not_limited_pw(geniter),NCOST)-1;
1436 
1437                         % Get breakpoints for this curve
1438                         x_pw = mpc.gencost(gen_not_limited_pw(geniter),-1+COST+(1:2:2*nseg+1)) / Sbase; % piecewise powers
1439                         c_pw = mpc.gencost(gen_not_limited_pw(geniter),-1+COST+(2:2:2*nseg+2)); % piecewise costs
1440 
1441                         for pwidx = 1:nseg
1442                             m_pw = (c_pw(pwidx+1) - c_pw(pwidx)) / (x_pw(pwidx+1) - x_pw(pwidx)); % slope for this piecewise cost segement
1443                             if x_pw(pwidx) > mpc.gen(gen_not_limited_pw(geniter),PG)/Sbase
1444                                 x_min = x_pw(pwidx-1);
1445                                 x_max = x_pw(pwidx);
1446                                 break;
1447                             end
1448                         end
1449 
1450                         if m_pw < lc
1451                             lc = m_pw;
1452                             lc_x_min = x_min;
1453                             lc_x_max = x_max;
1454                             lc_geniter = geniter;
1455                         end
1456 
1457                     end
1458                     if mpc.gen(gen_not_limited_pw(lc_geniter), PG) + min(Pgen_remain, lc_x_max-lc_x_min)*Sbase < Pmax(lc_geniter)
1459                         % If we don't hit the upper gen limit, add up to the
1460                         % breakpoint for this segment to the generation.
1461                         mpc.gen(gen_not_limited_pw(lc_geniter), PG) = mpc.gen(gen_not_limited_pw(lc_geniter), PG) + min(Pgen_remain, lc_x_max-lc_x_min)*Sbase;
1462                         Pgen_remain = Pgen_remain - min(Pgen_remain, lc_x_max-lc_x_min);
1463                     else
1464                         % If we hit the upper generation limit, set the
1465                         % generator to this limit and remove it from the list
1466                         % of genenerators below their limit.
1467                         Pgen_remain = Pgen_remain - (Pmax(lc_geniter) - mpc.gen(gen_not_limited_pw(lc_geniter), PG));
1468                         mpc.gen(gen_not_limited_pw(lc_geniter), PG) = Pmax(lc_geniter);
1469                         gen_not_limited_pw(lc_geniter) = [];
1470                         Pmax(lc_geniter) = [];
1471                         Pmin(lc_geniter) = [];
1472                     end
1473                 end
1474             end
1475             
1476             
1477             % For generators with quadratic cost functions:
1478             % Solve the quadratic program that is the economic dispatch
1479             % problem at this bus.
1480             if ~isempty(gen_not_limited_q)
1481                 
1482                 Pg_q = sdpvar(length(gen_not_limited_q),1);
1483                 cost_q = Pg_q.' * diag(c2(gen_not_limited_q)) * Pg_q + c1(gen_not_limited_q).' * Pg_q;
1484                 
1485                 % Specifying no generator costs leads to numeric problems
1486                 % in the solver. If no generator costs are specified,
1487                 % assign uniform costs for all generators to assign active
1488                 % power generation at this bus.
1489                 if isnumeric(cost_q)
1490                     cost_q = sum(Pg_q);
1491                 end
1492                 
1493                 constraints_q = [Pmin/Sbase <= Pg_q <= Pmax/Sbase;
1494                                  sum(Pg_q) == Pgen_remain];
1495                 sol_q = solvesdp(constraints_q,cost_q, sdpsettings(sdpopts,'Verbose',0,'saveduals',0));
1496                 
1497                 % This problem should have a feasible solution if we have a
1498                 % rank 1 solution to the OPF problem. However, it is
1499                 % possible that no set of power injections satisfies the
1500                 % power generation constraints if we recover the closest
1501                 % rank 1 solution from a higher rank solution to the OPF
1502                 % problem.
1503                 
1504                 constraint_q_err = checkset(constraints_q);
1505                 if sol_q.problem == 0 || max(abs(constraint_q_err)) < 1e-3 % If error in assigning active power generation is less than 0.1 MW, ignore it
1506                     mpc.gen(gen_not_limited_q, PG) = double(Pg_q) * Sbase;
1507                 elseif sol_q.problem == 4 || sol_q.problem == 1
1508                     % Solutions with errors of "numeric problems" or "infeasible" are
1509                     % typically pretty close. Use the solution but give a
1510                     % warning.
1511                     mpc.gen(gen_not_limited_q, PG) = double(Pg_q) * Sbase;
1512                     warning('sdpopf_solver: Numeric inconsistency assigning active power injection to generators at bus %i.',busorder(i));
1513                 else
1514                     % If there is any other error, just assign active power
1515                     % generation equally among all generators without
1516                     % regard for active power limits.
1517                     mpc.gen(gen_not_limited_q, PG) = (Pgen_remain / length(gen_not_limited_q))*Sbase;
1518                     warning('sdpopf_solver: Error assigning active power injection to generators at bus %i. Assiging equally among all generators at this bus.',busorder(i));
1519                 end
1520                 
1521             end
1522             
1523         end
1524 
1525         % Now assign reactive power generation. <-- Problem: positive
1526         % reactive power lower limits causes problems. Initialize all
1527         % generators to lower limits and increase as necessary.
1528         mpc.gen(genidx,QG) = mpc.gen(genidx,QMIN);
1529         Qremaining = (Qinj(i) + Qd(i))*Sbase - sum(mpc.gen(genidx,QG));
1530         genidx_idx = 1;
1531         while abs(Qremaining) > Qtol && genidx_idx <= length(genidx)
1532             if mpc.gen(genidx(genidx_idx),QG) + Qremaining <= mpc.gen(genidx(genidx_idx),QMAX) % && mpc.gen(genidx(genidx_idx),QG) + Qremaining >= mpc.gen(genidx(genidx_idx),QMIN)
1533                 mpc.gen(genidx(genidx_idx),QG) = mpc.gen(genidx(genidx_idx),QG) + Qremaining;
1534                 Qremaining = 0;
1535             else
1536                 Qremaining = Qremaining - (mpc.gen(genidx(genidx_idx),QMAX) - mpc.gen(genidx(genidx_idx),QG));
1537                 mpc.gen(genidx(genidx_idx),QG) = mpc.gen(genidx(genidx_idx),QMAX);
1538             end
1539             genidx_idx = genidx_idx + 1;
1540         end
1541         if abs(Qremaining) > Qtol
1542             % If the total reactive generation specified is greater than
1543             % the reactive generation available, assign the remainder to
1544             % equally to all generators at the bus.
1545             mpc.gen(genidx,QG) = mpc.gen(genidx,QG) + Qremaining / length(genidx);
1546             warning('sdpopf_solver: Inconsistency in assigning reactive power to generators at bus index %i.',i);
1547         end
1548     end
1549 end
1550 
1551 % Line flows (PF, PT, QF, QT)
1552 mpc.branch(:,PF) = Pflowft;
1553 mpc.branch(:,PT) = Pflowtf;
1554 mpc.branch(:,QF) = Qflowft;
1555 mpc.branch(:,QT) = Qflowtf;
1556 
1557 
1558 % Store bus Lagrange multipliers into results
1559 mpc.bus(:,LAM_P) = -double(lambdaeq_sdp) / Sbase;
1560 for i=1:nbus
1561     if ~isnan(gamma_eq(i))
1562         mpc.bus(i,LAM_Q) = -double(gammaeq_sdp(gamma_eq(i))) / Sbase;
1563     elseif Qmin(i) >= -maxgenlimit && Qmax(i) <= maxgenlimit % bus has both upper and lower Q limits
1564         mpc.bus(i,LAM_Q) = (double(gammaU_sdp(gamma_ineq(i))) - double(gammaL_sdp(gamma_ineq(i)))) / Sbase;
1565     elseif Qmax(i) <= maxgenlimit % bus has only upper Q limit
1566         mpc.bus(i,LAM_Q) = double(gammaU_sdp(gamma_ineq(i))) / Sbase;
1567     elseif Qmin(i) >= -maxgenlimit % bus has only lower Q limit
1568         mpc.bus(i,LAM_Q) = double(gammaL_sdp(gamma_ineq(i))) / Sbase;
1569     else % bus has no reactive power limits, so cost for reactive power is zero
1570         mpc.bus(i,LAM_Q) = 0;
1571     end
1572 end
1573 
1574 % Convert Lagrange multipliers to be in terms of voltage magntiude rather
1575 % than voltage magnitude squared.
1576 
1577 mpc.bus(:,MU_VMAX) = 2*abs(Vopt) .* double(muU_sdp);
1578 mpc.bus(:,MU_VMIN) = 2*abs(Vopt) .* double(muL_sdp);
1579 
1580 % Store gen Lagrange multipliers into results
1581 for i=1:nbus
1582     genidx = find(mpc.gen(:,GEN_BUS) == i);
1583     for k=1:length(genidx)
1584         if ~isnan(double(psiU_sdp{i}(k))) && ~isnan(double(psiU_sdp{i}(k)))
1585             mpc.gen(genidx(k),MU_PMAX) = double(psiU_sdp{i}(k)) / Sbase;
1586             mpc.gen(genidx(k),MU_PMIN) = double(psiL_sdp{i}(k)) / Sbase;
1587         else % Handle generators with tight generation limits
1588 
1589             % Calculate generalized Lagrange multiplier R for this generator
1590             % To do so, first calculate the dual of R (primal matrix)
1591             Pavg = (mpc.gen(genidx(k),PMAX)+mpc.gen(genidx(k),PMIN))/(2*Sbase);
1592             c0k = mpc.gencost(genidx(k),COST+2);
1593             alphak = c2(genidx(k))*Pavg^2 + c1(genidx(k))*Pavg + c0k;
1594             
1595             %  dualR = [-c1(genidx(k))*Pavg+alphak-c0k -sqrt(c2(genidx(k)))*Pavg;
1596             %          -sqrt(c2(genidx(k)))*Pavg 1] >= 0;
1597              
1598             % At the optimal solution, trace(dualR*R) = 0 and det(R) = 0
1599             %  -c1(genidx(k))*Pavg+alphak-c0k + R22 + -2*sqrt(c2(genidx(k)))*Pavg*R12 = 0 % trace(dualR*R) = 0
1600             %  R22 - R12^2 = 0                                                            % det(R) = 0
1601             
1602             % Rearranging these equations gives
1603             % -c1(genidx(k))*Pavg+alphak-c0k + R12^2 - 2*sqrt(c2(genidx(k)))*Pavg*R12 = 0
1604             
1605             % Solve with the quadratic equation
1606             a = 1;
1607             b = -2*sqrt(c2(genidx(k)))*Pavg;
1608             c = -c1(genidx(k))*Pavg+alphak-c0k;
1609             R12 = (-b+sqrt(b^2-4*a*c)) / (2*a);
1610             if abs(imag(R12)) > 1e-3
1611                 warning('sdpopf_solver: R12 has a nonzero imaginary component. Lagrange multipliers for generator active power limits may be incorrect.');
1612             end
1613             R12 = real(R12);
1614             
1615             % With R known, calculate psi
1616             % c1(genidx(i)) + 2*sqrt(c2(genidx(i)))*R{k}{i}(1) + psi == lambdai
1617             psi_tight = -double(lambdaeq_sdp(lambda_eq(i))) / Sbase - c1(genidx(k)) / Sbase - 2*sqrt(c2(genidx(k))/Sbase^2)*R12;
1618             if psi_tight < 0
1619                 mpc.gen(genidx(k),MU_PMAX) = 0;
1620                 mpc.gen(genidx(k),MU_PMIN) = -psi_tight;
1621             else
1622                 mpc.gen(genidx(k),MU_PMAX) = psi_tight;
1623                 mpc.gen(genidx(k),MU_PMIN) = 0;
1624             end
1625         end
1626     end
1627     if ~isempty(genidx)
1628         if ~isnan(gamma_eq(i))
1629             gamma_tight = double(gammaeq_sdp(gamma_eq(i))) / Sbase;
1630             if gamma_tight > 0 
1631                 mpc.gen(genidx(k),MU_QMAX) = 0;
1632                 mpc.gen(genidx(k),MU_QMIN) = double(gammaeq_sdp(gamma_eq(i))) / Sbase;
1633             else
1634                 mpc.gen(genidx(k),MU_QMAX) = -double(gammaeq_sdp(gamma_eq(i))) / Sbase;
1635                 mpc.gen(genidx(k),MU_QMIN) = 0;
1636             end
1637         else
1638             if any(Qmin(i) < -maxgenlimit) % bus has no lower Q limit
1639                 mpc.gen(genidx,MU_QMIN) = 0;
1640             else
1641                 mpc.gen(genidx,MU_QMIN) = double(gammaL_sdp(gamma_ineq(i))) / Sbase;
1642             end
1643 
1644             if any(Qmax(i) > maxgenlimit) % bus has no upper Q limit
1645                 mpc.gen(genidx,MU_QMAX) = 0;
1646             else
1647                 mpc.gen(genidx,MU_QMAX) = double(gammaU_sdp(gamma_ineq(i))) / Sbase;
1648             end
1649         end
1650     end
1651 end
1652 
1653 % Line flow Lagrange multipliers are calculated in terms of squared
1654 % active and reactive power line flows.
1655 %
1656 % Calculate line-flow limit Lagrange multipliers for MVA limits as follows.
1657 % temp = Hsdp{i} ./ Sbase;
1658 % 2*sqrt((temp(1,2)^2 + temp(1,3)^2))
1659 if nlconstraint > 0
1660     for i=1:nbranch
1661         Hidx = Hlookup(:,1) == i & Hlookup(:,2) == 1;
1662         if any(Hidx)
1663             if upper(mpopt.opf.flow_lim(1)) == 'S'
1664                 Hft = double(Hsdp{Hidx}) ./ Sbase;
1665                 mpc.branch(i,MU_SF) = 2*sqrt((Hft(1,2)^2 + Hft(1,3)^2));
1666                 Htf = double(Hsdp{Hlookup(:,1) == i & Hlookup(:,2) == 0}) ./ Sbase;
1667                 mpc.branch(i,MU_ST) = 2*sqrt((Htf(1,2)^2 + Htf(1,3)^2));
1668             elseif upper(mpopt.opf.flow_lim(1)) == 'P'
1669                 mpc.branch(i,MU_SF) = double(Hsdp(Hidx)) / Sbase;
1670                 mpc.branch(i,MU_ST) = double(Hsdp(Hlookup(:,1) == i & Hlookup(:,2) == 0)) / Sbase;
1671             end
1672         else
1673             mpc.branch(i,MU_SF) = 0;
1674             mpc.branch(i,MU_ST) = 0;
1675         end
1676     end
1677 else
1678     mpc.branch(:,MU_SF) = 0;
1679     mpc.branch(:,MU_ST) = 0;
1680 end
1681 
1682 % Angle limits are not implemented
1683 mpc.branch(:,MU_ANGMIN) = nan;
1684 mpc.branch(:,MU_ANGMAX) = nan;
1685 
1686 % Objective value
1687 mpc.f = objval;
1688 
1689 %% Convert back to original bus ordering
1690 
1691 mpc.bus(:,BUS_I) = busorder;
1692 mpc.bus = sortrows(mpc.bus,1);
1693 
1694 for i=1:ngen
1695     mpc.gen(i,GEN_BUS) = busorder(mpc.gen(i,GEN_BUS));
1696 end
1697 
1698 [mpc.gen genidx] = sortrows(mpc.gen,1);
1699 mpc.gencost = mpc.gencost(genidx,:);
1700 
1701 for i=1:nbranch
1702     mpc.branch(i,F_BUS) = busorder(mpc.branch(i,F_BUS));
1703     mpc.branch(i,T_BUS) = busorder(mpc.branch(i,T_BUS));
1704 end
1705 
1706 
1707 %% Generate outputs
1708 
1709 results = mpc;
1710 success = mineigratio > mineigratio_tol && LLeval < zeroeval_tol;
1711 
1712 raw.xr.A = A; 
1713 if recover_injections == 2
1714     raw.xr.W = W;
1715 else
1716     raw.xr.W = [];
1717 end
1718 
1719 raw.pimul.lambdaeq_sdp = lambdaeq_sdp;
1720 raw.pimul.gammaeq_sdp = gammaeq_sdp;
1721 raw.pimul.gammaU_sdp = gammaU_sdp;
1722 raw.pimul.gammaL_sdp = gammaL_sdp;
1723 raw.pimul.muU_sdp = muU_sdp;
1724 raw.pimul.muL_sdp = muL_sdp;
1725 raw.pimul.psiU_sdp = psiU_sdp;
1726 raw.pimul.psiL_sdp = psiL_sdp;
1727 
1728 raw.info = sdpinfo.problem;
1729 
1730 results.zero_eval = LLeval;
1731 results.mineigratio = mineigratio;
1732 
1733 results.sdpinfo = sdpinfo;
1734 
1735 toutput = toc;
1736 results.et = tsetup + tmaxclique + tsdpform + tsolve + toutput;
sdpopf_solver

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
