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

Generated on Fri 16-Dec-2016 12:45:37 by m2html © 2005