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

Generated on Mon 24-Jun-2019 15:58:45 by m2html © 2005