Home > matpower5.0 > 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 %   $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

Generated on Mon 26-Jan-2015 15:21:31 by m2html © 2005