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 % $Id: testGlobalOpt.m 2274 2014-01-17 15:56:29Z ray $ 0064 % by Daniel Molzahn, PSERC U of Wisc, Madison 0065 % and Ray Zimmerman, PSERC Cornell 0066 % Copyright (c) 2013-2014 by Power System Engineering Research Center (PSERC) 0067 % 0068 % This file is part of MATPOWER. 0069 % See http://www.pserc.cornell.edu/matpower/ for more info. 0070 % 0071 % MATPOWER is free software: you can redistribute it and/or modify 0072 % it under the terms of the GNU General Public License as published 0073 % by the Free Software Foundation, either version 3 of the License, 0074 % or (at your option) any later version. 0075 % 0076 % MATPOWER is distributed in the hope that it will be useful, 0077 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0078 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0079 % GNU General Public License for more details. 0080 % 0081 % You should have received a copy of the GNU General Public License 0082 % along with MATPOWER. If not, see <http://www.gnu.org/licenses/>. 0083 % 0084 % Additional permission under GNU GPL version 3 section 7 0085 % 0086 % If you modify MATPOWER, or any covered work, to interface with 0087 % other modules (such as MATLAB code and MEX-files) available in a 0088 % MATLAB(R) or comparable environment containing parts covered 0089 % under other licensing terms, the licensors of MATPOWER grant 0090 % you additional permission to convey the resulting work. 0091 0092 %% define named indices into data matrices 0093 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ... 0094 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus; 0095 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ... 0096 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ... 0097 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch; 0098 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost; 0099 0100 %% Load in options, give warnings 0101 0102 if nargin < 2 0103 mpopt = mpoption; 0104 end 0105 0106 if nargin == 0 0107 error('testGlobalOpt: testGlobalOpt requires a solved optimal power flow problem.'); 0108 end 0109 0110 mpc = loadcase(mpc); 0111 mpc = ext2int(mpc); 0112 0113 ignore_angle_lim = mpopt.opf.ignore_angle_lim; 0114 binding_lagrange = mpopt.sdp_pf.bind_lagrange; %% Tolerance for considering a Lagrange multiplier to indicae a binding constraint 0115 zeroeval_tol = mpopt.sdp_pf.zeroeval_tol; %% Tolerance for considering an eigenvalue in LLeval equal to zero 0116 comp_tol = sqrt(zeroeval_tol); %% Tolerance for considering trace(A*W) equal to zero 0117 0118 if size(mpc.bus,2) < LAM_P 0119 error('testGlobalOpt: testGlobalOpt requires a solved optimal power flow problem.'); 0120 end 0121 0122 if toggle_dcline(mpc, 'status') 0123 error('testGlobalOpt: DC lines are not implemented in testGlobalOpt'); 0124 end 0125 0126 nbus = size(mpc.bus,1); 0127 Sbase = mpc.baseMVA; 0128 0129 if ~ignore_angle_lim && (any(mpc.branch(:,ANGMIN) ~= -360) || any(mpc.branch(:,ANGMAX) ~= 360)) 0130 warning('testGlobalOpt: Angle difference constraints are not implemented in testGlobalOpt. Ignoring angle difference constraints.'); 0131 end 0132 0133 0134 %% Create matrices 0135 0136 [Yk,Yk_,Mk,Ylineft,Ylinetf,Y_lineft,Y_linetf] = makesdpmat(mpc); 0137 0138 %% Form the A matrix 0139 0140 lambda = mpc.bus(:,LAM_P)*Sbase; 0141 gamma = mpc.bus(:,LAM_Q)*Sbase; 0142 0143 % Convert voltage Lagrange multiplier to be in terms of voltage squared 0144 % dLdV^2 * dV^2dV = dLdV --> dLdV^2 = dLdV / (2*V) 0145 mu_V = mpc.bus(:,MU_VMAX) - mpc.bus(:,MU_VMIN); 0146 mu = mu_V ./ (2*mpc.bus(:,VM)); 0147 0148 A = sparse(2*nbus,2*nbus); 0149 0150 % Add bus terms 0151 for k=1:nbus 0152 A = A + lambda(k)*Yk(k) + gamma(k)*Yk_(k) + mu(k)*Mk(k); 0153 end 0154 0155 % Add branch terms 0156 Pft = mpc.branch(:,PF) / Sbase; 0157 Ptf = mpc.branch(:,PT) / Sbase; 0158 Qft = mpc.branch(:,QF) / Sbase; 0159 Qtf = mpc.branch(:,QT) / Sbase; 0160 Sft = sqrt(Pft.^2 + Qft.^2); 0161 Stf = sqrt(Ptf.^2 + Qtf.^2); 0162 0163 dLdSft = mpc.branch(:,MU_SF) * Sbase; 0164 dLdStf = mpc.branch(:,MU_ST) * Sbase; 0165 0166 dLdSft_over_Sft = dLdSft ./ (Sft); 0167 dLdStf_over_Stf = dLdStf ./ (Stf); 0168 dLdPft = Pft .* dLdSft_over_Sft; 0169 dLdQft = Qft .* dLdSft_over_Sft; 0170 dLdPtf = Ptf .* dLdStf_over_Stf; 0171 dLdQtf = Qtf .* dLdStf_over_Stf; 0172 0173 binding_branches_ft = find(abs(dLdSft) > binding_lagrange); 0174 binding_branches_tf = find(abs(dLdStf) > binding_lagrange); 0175 0176 if ~isempty(binding_branches_ft) 0177 for ft_idx=1:length(binding_branches_ft) 0178 k = binding_branches_ft(ft_idx); 0179 if strcmp(upper(mpopt.opf.flow_lim), 'S') % apparent power limits 0180 0181 A = A + dLdPft(k)*Ylineft(k) ... 0182 + dLdQft(k)*Y_lineft(k); 0183 elseif strcmp(upper(mpopt.opf.flow_lim), 'P') % active power limits 0184 % Active power Lagrange multipliers are stored in mpc.branch(:,MU_SF) 0185 % No need for conversion 0186 0187 A = A + dLdSft(k) * Ylineft(k); 0188 0189 else 0190 error('testGlobalOpt: Only ''S'' and ''P'' options are currently implemented for MPOPT.opf.flow_lim\n'); 0191 end 0192 end 0193 end 0194 0195 if ~isempty(binding_branches_tf) 0196 for tf_idx=1:length(binding_branches_tf) 0197 k = binding_branches_tf(tf_idx); 0198 if strcmp(upper(mpopt.opf.flow_lim), 'S') % apparent power limits 0199 0200 A = A + dLdPtf(k)*Ylinetf(k) ... 0201 + dLdQtf(k)*Y_linetf(k); 0202 0203 elseif strcmp(upper(mpopt.opf.flow_lim), 'P') % active power limits 0204 % Active power Lagrange multipliers are stored in mpc.branch(:,MU_SF) 0205 % No need for conversion 0206 0207 A = A + dLdStf(k) * Ylinetf(k); 0208 0209 else 0210 error('testGlobalOpt: Only ''S'' and ''P'' options are currently implemented for MPOPT.opf.flow_lim\n'); 0211 end 0212 0213 end 0214 end 0215 0216 %% Calculate trace(A*W) 0217 0218 V = mpc.bus(:,VM).*(cos(mpc.bus(:,VA)*pi/180)+1i*sin(mpc.bus(:,VA)*pi/180)); 0219 x = [real(V); imag(V)]; 0220 0221 comp = 0; 0222 0223 [Arow,Acol,Aval] = find(A); 0224 0225 for i=1:length(Aval) 0226 comp = comp + Aval(i)*x(Arow(i))*x(Acol(i)); 0227 end 0228 0229 0230 %% Determine if A is psd using a Cholesky decomposition 0231 0232 Abar = A + zeroeval_tol*speye(size(A)); 0233 per = amd(Abar); 0234 [junk, p] = chol(Abar(per,per)); 0235 0236 Apsd = p == 0; 0237 0238 %% Calculate globalopt flag 0239 0240 if Apsd && abs(comp) < comp_tol 0241 globalopt = 1; 0242 else 0243 globalopt = 0; 0244 end