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).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