Home > matpower5.1 > extras > sdp_pf > testGlobalOpt.m

testGlobalOpt

PURPOSE ^

TESTGLOABALOPT A sufficient condition for OPF global optimality

SYNOPSIS ^

function [globalopt,comp,Apsd] = testGlobalOpt(mpc,mpopt)

DESCRIPTION ^

TESTGLOABALOPT A sufficient condition for OPF global optimality

   [GLOBALOPT,COMP,APSD] = TESTGLOBALOPT(MPC,MPOPT)

   Evaluates a sufficient condition for global optimality of a solution to
   an OPF problem based on the KKT conditions of a semidefinite relaxation
   of the OPF problem. The guarantee of global optimality requires
   satisfaction of the complementarity condition (trace(A*W) = 0) and the
   positive semidefinite A matrix condition (A >= 0). Note that failure
   to satisfy these conditions does not necessarily indicate that the
   solution is not globally optimal, but satisfaction gurantees global
   optimality. See [1] for further details.

   Inputs:
       MPC : A solved MATPOWER case (e.g., mpc = runopf(mpc,mpopt)).
       MPOPT : The options struct used to solve the MATPOWER case. If not
           specified, it is assumed to be the default mpoption.

   Outputs:
       GLOBALOPT : Flag indicating whether the complementarity condition
           and positive semidefinite A matrix constraints are satisfied to
           the tolerances specified in mpopt.
       COMP : Value of the complementarity condition trace(A*W). A value
           equal to zero within the tolerance sqrt(ZEROEVAL_TOL) satisfies
           the complementarity condition.
       APSD : Flag indicating if the A matrix is positive semidefinite to
           within the tolerance of ZEROEVAL_TOL specified in mpopt.
           Positive semidefiniteness is determined using a Cholesky
           factorization.

   Note that this function does not currently handle DC Lines or
   angle-difference constraints. Piecewise-linear cost functions should
   work, but were not extensively tested.

   Note that tight solution tolerances are typically required for good
   results.

   Positive branch resistances are not required. However, enforcing a 
   small branch resistance may result in better results. With a small 
   minimum branch resistance, global optimality is confirmed for all test
   cases up to and including the 118-bus system.


   Example: Both the complementarity and positive semidefinite conditions
   are satisfied for the IEEE 14-bus system.

           define_constants;
           mpc = loadcase('case14');
           mpc.branch(:,BR_R) = max(mpc.branch(:,BR_R), 1e-4);
           mpopt = mpoption('opf.ac.solver', 'MIPS', ...
               'opf.violation', 1e-10, 'mips.gradtol', 1e-10, ...
               'mips.comptol', 1e-10, 'mips.costtol', 1e-10);
           results = runopf(mpc, mpopt);
           [globalopt, comp, Apsd] = testGlobalOpt(results, mpopt)

 [1] D.K. Molzahn, B.C. Lesieutre, and C.L. DeMarco, "A Sufficient
     Condition for Global Optimality of Solutions to the Optimal Power
     Flow Problem," To appear in IEEE Transactions on Power Systems 
     (Letters).

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [globalopt,comp,Apsd] = testGlobalOpt(mpc,mpopt)
0002 %TESTGLOABALOPT A sufficient condition for OPF global optimality
0003 %
0004 %   [GLOBALOPT,COMP,APSD] = TESTGLOBALOPT(MPC,MPOPT)
0005 %
0006 %   Evaluates a sufficient condition for global optimality of a solution to
0007 %   an OPF problem based on the KKT conditions of a semidefinite relaxation
0008 %   of the OPF problem. The guarantee of global optimality requires
0009 %   satisfaction of the complementarity condition (trace(A*W) = 0) and the
0010 %   positive semidefinite A matrix condition (A >= 0). Note that failure
0011 %   to satisfy these conditions does not necessarily indicate that the
0012 %   solution is not globally optimal, but satisfaction gurantees global
0013 %   optimality. See [1] for further details.
0014 %
0015 %   Inputs:
0016 %       MPC : A solved MATPOWER case (e.g., mpc = runopf(mpc,mpopt)).
0017 %       MPOPT : The options struct used to solve the MATPOWER case. If not
0018 %           specified, it is assumed to be the default mpoption.
0019 %
0020 %   Outputs:
0021 %       GLOBALOPT : Flag indicating whether the complementarity condition
0022 %           and positive semidefinite A matrix constraints are satisfied to
0023 %           the tolerances specified in mpopt.
0024 %       COMP : Value of the complementarity condition trace(A*W). A value
0025 %           equal to zero within the tolerance sqrt(ZEROEVAL_TOL) satisfies
0026 %           the complementarity condition.
0027 %       APSD : Flag indicating if the A matrix is positive semidefinite to
0028 %           within the tolerance of ZEROEVAL_TOL specified in mpopt.
0029 %           Positive semidefiniteness is determined using a Cholesky
0030 %           factorization.
0031 %
0032 %   Note that this function does not currently handle DC Lines or
0033 %   angle-difference constraints. Piecewise-linear cost functions should
0034 %   work, but were not extensively tested.
0035 %
0036 %   Note that tight solution tolerances are typically required for good
0037 %   results.
0038 %
0039 %   Positive branch resistances are not required. However, enforcing a
0040 %   small branch resistance may result in better results. With a small
0041 %   minimum branch resistance, global optimality is confirmed for all test
0042 %   cases up to and including the 118-bus system.
0043 %
0044 %
0045 %   Example: Both the complementarity and positive semidefinite conditions
0046 %   are satisfied for the IEEE 14-bus system.
0047 %
0048 %           define_constants;
0049 %           mpc = loadcase('case14');
0050 %           mpc.branch(:,BR_R) = max(mpc.branch(:,BR_R), 1e-4);
0051 %           mpopt = mpoption('opf.ac.solver', 'MIPS', ...
0052 %               'opf.violation', 1e-10, 'mips.gradtol', 1e-10, ...
0053 %               'mips.comptol', 1e-10, 'mips.costtol', 1e-10);
0054 %           results = runopf(mpc, mpopt);
0055 %           [globalopt, comp, Apsd] = testGlobalOpt(results, mpopt)
0056 %
0057 % [1] D.K. Molzahn, B.C. Lesieutre, and C.L. DeMarco, "A Sufficient
0058 %     Condition for Global Optimality of Solutions to the Optimal Power
0059 %     Flow Problem," To appear in IEEE Transactions on Power Systems
0060 %     (Letters).
0061 
0062 %   MATPOWER
0063 %   Copyright (c) 2013-2015 by Power System Engineering Research Center (PSERC)
0064 %   by Daniel Molzahn, PSERC U of Wisc, Madison
0065 %   and Ray Zimmerman, PSERC Cornell
0066 %
0067 %   $Id: testGlobalOpt.m 2644 2015-03-11 19:34:22Z ray $
0068 %
0069 %   This file is part of MATPOWER.
0070 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0071 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0072 
0073 %% define named indices into data matrices
0074 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0075     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0076 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0077     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0078     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0079 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
0080 
0081 %% Load in options, give warnings
0082 
0083 if nargin < 2
0084     mpopt = mpoption;
0085 end
0086 
0087 if nargin == 0
0088     error('testGlobalOpt: testGlobalOpt requires a solved optimal power flow problem.');
0089 end
0090 
0091 mpc = loadcase(mpc);
0092 mpc = ext2int(mpc);
0093 
0094 ignore_angle_lim    = mpopt.opf.ignore_angle_lim;
0095 binding_lagrange    = mpopt.sdp_pf.bind_lagrange;   %% Tolerance for considering a Lagrange multiplier to indicae a binding constraint
0096 zeroeval_tol        = mpopt.sdp_pf.zeroeval_tol;    %% Tolerance for considering an eigenvalue in LLeval equal to zero
0097 comp_tol            = sqrt(zeroeval_tol);       %% Tolerance for considering trace(A*W) equal to zero
0098 
0099 if size(mpc.bus,2) < LAM_P
0100     error('testGlobalOpt: testGlobalOpt requires a solved optimal power flow problem.');
0101 end
0102 
0103 if toggle_dcline(mpc, 'status')
0104     error('testGlobalOpt: DC lines are not implemented in testGlobalOpt');
0105 end
0106 
0107 nbus = size(mpc.bus,1);
0108 Sbase = mpc.baseMVA;
0109 
0110 if ~ignore_angle_lim && (any(mpc.branch(:,ANGMIN) ~= -360) || any(mpc.branch(:,ANGMAX) ~= 360))
0111     warning('testGlobalOpt: Angle difference constraints are not implemented in testGlobalOpt. Ignoring angle difference constraints.');
0112 end
0113 
0114 
0115 %% Create matrices
0116 
0117 [Yk,Yk_,Mk,Ylineft,Ylinetf,Y_lineft,Y_linetf] = makesdpmat(mpc);
0118 
0119 %% Form the A matrix
0120 
0121 lambda = mpc.bus(:,LAM_P)*Sbase;
0122 gamma = mpc.bus(:,LAM_Q)*Sbase;
0123 
0124 % Convert voltage Lagrange multiplier to be in terms of voltage squared
0125 % dLdV^2 * dV^2dV = dLdV --> dLdV^2 = dLdV / (2*V)
0126 mu_V = mpc.bus(:,MU_VMAX) - mpc.bus(:,MU_VMIN);
0127 mu = mu_V ./ (2*mpc.bus(:,VM));
0128 
0129 A = sparse(2*nbus,2*nbus);
0130 
0131 % Add bus terms
0132 for k=1:nbus
0133     A = A + lambda(k)*Yk(k) + gamma(k)*Yk_(k) + mu(k)*Mk(k);
0134 end
0135 
0136 % Add branch terms
0137 Pft = mpc.branch(:,PF) / Sbase;
0138 Ptf = mpc.branch(:,PT) / Sbase;
0139 Qft = mpc.branch(:,QF) / Sbase;
0140 Qtf = mpc.branch(:,QT) / Sbase;
0141 Sft = sqrt(Pft.^2 + Qft.^2);
0142 Stf = sqrt(Ptf.^2 + Qtf.^2);
0143 
0144 dLdSft = mpc.branch(:,MU_SF) * Sbase;
0145 dLdStf = mpc.branch(:,MU_ST) * Sbase;
0146 
0147 dLdSft_over_Sft = dLdSft ./ (Sft);
0148 dLdStf_over_Stf = dLdStf ./ (Stf);
0149 dLdPft = Pft .* dLdSft_over_Sft;
0150 dLdQft = Qft .* dLdSft_over_Sft;
0151 dLdPtf = Ptf .* dLdStf_over_Stf;
0152 dLdQtf = Qtf .* dLdStf_over_Stf;
0153 
0154 binding_branches_ft = find(abs(dLdSft) > binding_lagrange);
0155 binding_branches_tf = find(abs(dLdStf) > binding_lagrange);
0156 
0157 if ~isempty(binding_branches_ft)
0158     for ft_idx=1:length(binding_branches_ft)
0159         k = binding_branches_ft(ft_idx);
0160         if strcmp(upper(mpopt.opf.flow_lim), 'S')       % apparent power limits
0161             
0162             A = A + dLdPft(k)*Ylineft(k) ...
0163                   + dLdQft(k)*Y_lineft(k);
0164         elseif strcmp(upper(mpopt.opf.flow_lim), 'P')   % active power limits
0165             % Active power Lagrange multipliers are stored in mpc.branch(:,MU_SF)
0166             % No need for conversion
0167             
0168             A = A + dLdSft(k) * Ylineft(k);
0169 
0170         else
0171             error('testGlobalOpt: Only ''S'' and ''P'' options are currently implemented for MPOPT.opf.flow_lim\n');
0172         end
0173     end
0174 end
0175 
0176 if ~isempty(binding_branches_tf)
0177     for tf_idx=1:length(binding_branches_tf)
0178         k = binding_branches_tf(tf_idx);
0179         if strcmp(upper(mpopt.opf.flow_lim), 'S')       % apparent power limits
0180 
0181             A = A + dLdPtf(k)*Ylinetf(k) ...
0182                   + dLdQtf(k)*Y_linetf(k);
0183 
0184         elseif strcmp(upper(mpopt.opf.flow_lim), 'P')   % active power limits
0185             % Active power Lagrange multipliers are stored in mpc.branch(:,MU_SF)
0186             % No need for conversion
0187 
0188             A = A + dLdStf(k) * Ylinetf(k);
0189 
0190         else
0191             error('testGlobalOpt: Only ''S'' and ''P'' options are currently implemented for MPOPT.opf.flow_lim\n');
0192         end
0193 
0194     end
0195 end
0196 
0197 %% Calculate trace(A*W)
0198 
0199 V = mpc.bus(:,VM).*(cos(mpc.bus(:,VA)*pi/180)+1i*sin(mpc.bus(:,VA)*pi/180));
0200 x = [real(V); imag(V)];
0201 
0202 comp = 0;
0203 
0204 [Arow,Acol,Aval] = find(A);
0205 
0206 for i=1:length(Aval)
0207     comp = comp + Aval(i)*x(Arow(i))*x(Acol(i));
0208 end
0209 
0210 
0211 %% Determine if A is psd using a Cholesky decomposition
0212 
0213 Abar = A + zeroeval_tol*speye(size(A));
0214 per = amd(Abar);
0215 [junk, p] = chol(Abar(per,per));
0216 
0217 Apsd = p == 0;
0218 
0219 %% Calculate globalopt flag
0220 
0221 if Apsd && abs(comp) < comp_tol
0222     globalopt = 1;
0223 else
0224     globalopt = 0;
0225 end

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