Home > matpower7.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-2019, 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/mx-sdp_pf.
0068 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0069 %   See https://github.com/MATPOWER/mx-sdp_pf/ 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

Generated on Fri 09-Oct-2020 11:21:31 by m2html © 2005