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-2016, Power Systems Engineering Research Center (PSERC) 0064 % by Daniel Molzahn, PSERC U of Wisc, Madison 0065 % and Ray Zimmerman, PSERC Cornell 0066 % 0067 % This file is part of MATPOWER. 0068 % Covered by the 3-clause BSD License (see LICENSE file for details). 0069 % See http://www.pserc.cornell.edu/matpower/ for more info. 0070 0071 %% define named indices into data matrices 0072 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ... 0073 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus; 0074 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ... 0075 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ... 0076 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch; 0077 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost; 0078 0079 %% Load in options, give warnings 0080 0081 if nargin < 2 0082 mpopt = mpoption; 0083 end 0084 0085 if nargin == 0 0086 error('testGlobalOpt: testGlobalOpt requires a solved optimal power flow problem.'); 0087 end 0088 0089 mpc = loadcase(mpc); 0090 mpc = ext2int(mpc); 0091 0092 ignore_angle_lim = mpopt.opf.ignore_angle_lim; 0093 binding_lagrange = mpopt.sdp_pf.bind_lagrange; %% Tolerance for considering a Lagrange multiplier to indicae a binding constraint 0094 zeroeval_tol = mpopt.sdp_pf.zeroeval_tol; %% Tolerance for considering an eigenvalue in LLeval equal to zero 0095 comp_tol = sqrt(zeroeval_tol); %% Tolerance for considering trace(A*W) equal to zero 0096 0097 if size(mpc.bus,2) < LAM_P 0098 error('testGlobalOpt: testGlobalOpt requires a solved optimal power flow problem.'); 0099 end 0100 0101 if toggle_dcline(mpc, 'status') 0102 error('testGlobalOpt: DC lines are not implemented in testGlobalOpt'); 0103 end 0104 0105 nbus = size(mpc.bus,1); 0106 Sbase = mpc.baseMVA; 0107 0108 if ~ignore_angle_lim && (any(mpc.branch(:,ANGMIN) ~= -360) || any(mpc.branch(:,ANGMAX) ~= 360)) 0109 warning('testGlobalOpt: Angle difference constraints are not implemented in testGlobalOpt. Ignoring angle difference constraints.'); 0110 end 0111 0112 0113 %% Create matrices 0114 0115 [Yk,Yk_,Mk,Ylineft,Ylinetf,Y_lineft,Y_linetf] = makesdpmat(mpc); 0116 0117 %% Form the A matrix 0118 0119 lambda = mpc.bus(:,LAM_P)*Sbase; 0120 gamma = mpc.bus(:,LAM_Q)*Sbase; 0121 0122 % Convert voltage Lagrange multiplier to be in terms of voltage squared 0123 % dLdV^2 * dV^2dV = dLdV --> dLdV^2 = dLdV / (2*V) 0124 mu_V = mpc.bus(:,MU_VMAX) - mpc.bus(:,MU_VMIN); 0125 mu = mu_V ./ (2*mpc.bus(:,VM)); 0126 0127 A = sparse(2*nbus,2*nbus); 0128 0129 % Add bus terms 0130 for k=1:nbus 0131 A = A + lambda(k)*Yk(k) + gamma(k)*Yk_(k) + mu(k)*Mk(k); 0132 end 0133 0134 % Add branch terms 0135 Pft = mpc.branch(:,PF) / Sbase; 0136 Ptf = mpc.branch(:,PT) / Sbase; 0137 Qft = mpc.branch(:,QF) / Sbase; 0138 Qtf = mpc.branch(:,QT) / Sbase; 0139 Sft = sqrt(Pft.^2 + Qft.^2); 0140 Stf = sqrt(Ptf.^2 + Qtf.^2); 0141 0142 dLdSft = mpc.branch(:,MU_SF) * Sbase; 0143 dLdStf = mpc.branch(:,MU_ST) * Sbase; 0144 0145 dLdSft_over_Sft = dLdSft ./ (Sft); 0146 dLdStf_over_Stf = dLdStf ./ (Stf); 0147 dLdPft = Pft .* dLdSft_over_Sft; 0148 dLdQft = Qft .* dLdSft_over_Sft; 0149 dLdPtf = Ptf .* dLdStf_over_Stf; 0150 dLdQtf = Qtf .* dLdStf_over_Stf; 0151 0152 binding_branches_ft = find(abs(dLdSft) > binding_lagrange); 0153 binding_branches_tf = find(abs(dLdStf) > binding_lagrange); 0154 0155 if ~isempty(binding_branches_ft) 0156 for ft_idx=1:length(binding_branches_ft) 0157 k = binding_branches_ft(ft_idx); 0158 if strcmp(upper(mpopt.opf.flow_lim), 'S') % apparent power limits 0159 0160 A = A + dLdPft(k)*Ylineft(k) ... 0161 + dLdQft(k)*Y_lineft(k); 0162 elseif strcmp(upper(mpopt.opf.flow_lim), 'P') % active power limits 0163 % Active power Lagrange multipliers are stored in mpc.branch(:,MU_SF) 0164 % No need for conversion 0165 0166 A = A + dLdSft(k) * Ylineft(k); 0167 0168 else 0169 error('testGlobalOpt: Only ''S'' and ''P'' options are currently implemented for MPOPT.opf.flow_lim\n'); 0170 end 0171 end 0172 end 0173 0174 if ~isempty(binding_branches_tf) 0175 for tf_idx=1:length(binding_branches_tf) 0176 k = binding_branches_tf(tf_idx); 0177 if strcmp(upper(mpopt.opf.flow_lim), 'S') % apparent power limits 0178 0179 A = A + dLdPtf(k)*Ylinetf(k) ... 0180 + dLdQtf(k)*Y_linetf(k); 0181 0182 elseif strcmp(upper(mpopt.opf.flow_lim), 'P') % active power limits 0183 % Active power Lagrange multipliers are stored in mpc.branch(:,MU_SF) 0184 % No need for conversion 0185 0186 A = A + dLdStf(k) * Ylinetf(k); 0187 0188 else 0189 error('testGlobalOpt: Only ''S'' and ''P'' options are currently implemented for MPOPT.opf.flow_lim\n'); 0190 end 0191 0192 end 0193 end 0194 0195 %% Calculate trace(A*W) 0196 0197 V = mpc.bus(:,VM).*(cos(mpc.bus(:,VA)*pi/180)+1i*sin(mpc.bus(:,VA)*pi/180)); 0198 x = [real(V); imag(V)]; 0199 0200 comp = 0; 0201 0202 [Arow,Acol,Aval] = find(A); 0203 0204 for i=1:length(Aval) 0205 comp = comp + Aval(i)*x(Arow(i))*x(Acol(i)); 0206 end 0207 0208 0209 %% Determine if A is psd using a Cholesky decomposition 0210 0211 Abar = A + zeroeval_tol*speye(size(A)); 0212 per = amd(Abar); 0213 [junk, p] = chol(Abar(per,per)); 0214 0215 Apsd = p == 0; 0216 0217 %% Calculate globalopt flag 0218 0219 if Apsd && abs(comp) < comp_tol 0220 globalopt = 1; 0221 else 0222 globalopt = 0; 0223 end