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

Generated on Fri 20-Mar-2015 18:23:34 by m2html © 2005