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

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
