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

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
