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

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
