Home > matpower5.0 > extras > sdp_pf > sdpopf_solver.m

sdpopf_solver

PURPOSE ^

SDPOPF_SOLVER A semidefinite programming relaxtion of the OPF problem

SYNOPSIS ^

function [results,success,raw] = sdpopf_solver(om, mpopt)

DESCRIPTION ^

SDPOPF_SOLVER A semidefinite programming relaxtion of the OPF problem

   [RESULTS,SUCCESS,RAW] = SDPOPF_SOLVER(OM, MPOPT)

   Inputs are an OPF model object and a MATPOWER options vector.

   Outputs are a RESULTS struct, SUCCESS flag and RAW output struct.

   RESULTS is a MATPOWER case struct (mpc) with the usual baseMVA, bus
   branch, gen, gencost fields, along with the following additional
   fields:
       .f           final objective function value
       .mineigratio Minimum eigenvalue ratio measure of rank condition
                    satisfaction
       .zero_eval   Zero eigenvalue measure of optimal voltage profile
                    consistency

   SUCCESS     1 if solver converged successfully, 0 otherwise

   RAW         raw output
       .xr     final value of optimization variables
           .A  Cell array of matrix variables for the dual problem
           .W  Cell array of matrix variables for the primal problem
       .pimul  constraint multipliers on 
           .lambdaeq_sdp   Active power equalities
           .psiU_sdp       Generator active power upper inequalities
           .psiL_sdp       Generator active power lower inequalities
           .gammaeq_sdp    Reactive power equalities
           .gammaU_sdp     Reactive power upper inequalities
           .gammaL_sdp     Reactive power lower inequalities
           .muU_sdp        Square of voltage magnitude upper inequalities
           .muL_sdp        Square of voltage magnitude lower inequalities
       .info   solver specific termination code

   See also OPF.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated on Mon 26-Jan-2015 15:21:31 by m2html © 2005