------------------------------ deprecated ------------------------------
OPF solvers based on CONSTR & LPCONSTR to be removed in future version.
--------------------------------------------------------------------------
GRAD_COPF Evaluates gradient of objective function and constraints for OPF.
[DF, DG, D2F] = GRAD_COPF(X, OM, YBUS, YF, YT, AFEQ, BFEQ, AF, BF, MPOPT)
[DF, DG, D2F] = GRAD_COPF(X, OM, YBUS, YF, YT, AFEQ, BFEQ, AF, BF, MPOPT, IL)
Gradient (and Hessian) evaluation routine for AC optimal power flow costs
and constraints, suitable for use with CONSTR.
Inputs:
X : optimization vector
OM : OPF model object
YBUS : bus admittance matrix
YF : admittance matrix for "from" end of constrained branches
YT : admittance matrix for "to" end of constrained branches
AFEQ : sparse A matrix for linear equality constraints
BFEQ : right hand side vector for linear equality constraints
AF : sparse A matrix for linear inequality constraints
BF : right hand side vector for linear inequality constraints
MPOPF : MATPOWER options vector
IL : (optional) vector of branch indices corresponding to
branches with flow limits (all others are assumed to be
unconstrained). The default is [1:nl] (all branches).
YF and YT contain only the rows corresponding to IL.
Outputs:
DF : gradient of objective function (column vector)
DG : constraint gradients, column j is gradient of G(j)
see FUN_COPF for order of elements in G
D2F : (optional) Hessian of objective function (sparse matrix)
See also FUN_COPF.0001 function [df, dg, d2f] = grad_copf(x, om, Ybus, Yf, Yt, Afeq, bfeq, Af, bf, mpopt, il) 0002 %------------------------------ deprecated ------------------------------ 0003 % OPF solvers based on CONSTR & LPCONSTR to be removed in future version. 0004 %-------------------------------------------------------------------------- 0005 %GRAD_COPF Evaluates gradient of objective function and constraints for OPF. 0006 % [DF, DG, D2F] = GRAD_COPF(X, OM, YBUS, YF, YT, AFEQ, BFEQ, AF, BF, MPOPT) 0007 % [DF, DG, D2F] = GRAD_COPF(X, OM, YBUS, YF, YT, AFEQ, BFEQ, AF, BF, MPOPT, IL) 0008 % 0009 % Gradient (and Hessian) evaluation routine for AC optimal power flow costs 0010 % and constraints, suitable for use with CONSTR. 0011 % 0012 % Inputs: 0013 % X : optimization vector 0014 % OM : OPF model object 0015 % YBUS : bus admittance matrix 0016 % YF : admittance matrix for "from" end of constrained branches 0017 % YT : admittance matrix for "to" end of constrained branches 0018 % AFEQ : sparse A matrix for linear equality constraints 0019 % BFEQ : right hand side vector for linear equality constraints 0020 % AF : sparse A matrix for linear inequality constraints 0021 % BF : right hand side vector for linear inequality constraints 0022 % MPOPF : MATPOWER options vector 0023 % IL : (optional) vector of branch indices corresponding to 0024 % branches with flow limits (all others are assumed to be 0025 % unconstrained). The default is [1:nl] (all branches). 0026 % YF and YT contain only the rows corresponding to IL. 0027 % 0028 % Outputs: 0029 % DF : gradient of objective function (column vector) 0030 % DG : constraint gradients, column j is gradient of G(j) 0031 % see FUN_COPF for order of elements in G 0032 % D2F : (optional) Hessian of objective function (sparse matrix) 0033 % 0034 % See also FUN_COPF. 0035 0036 % MATPOWER 0037 % $Id: grad_copf.m,v 1.19 2010/04/27 18:55:02 ray Exp $ 0038 % by Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Autonoma de Manizales 0039 % and Ray Zimmerman, PSERC Cornell 0040 % Copyright (c) 1996-2010 by Power System Engineering Research Center (PSERC) 0041 % 0042 % This file is part of MATPOWER. 0043 % See http://www.pserc.cornell.edu/matpower/ for more info. 0044 % 0045 % MATPOWER is free software: you can redistribute it and/or modify 0046 % it under the terms of the GNU General Public License as published 0047 % by the Free Software Foundation, either version 3 of the License, 0048 % or (at your option) any later version. 0049 % 0050 % MATPOWER is distributed in the hope that it will be useful, 0051 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0052 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0053 % GNU General Public License for more details. 0054 % 0055 % You should have received a copy of the GNU General Public License 0056 % along with MATPOWER. If not, see <http://www.gnu.org/licenses/>. 0057 % 0058 % Additional permission under GNU GPL version 3 section 7 0059 % 0060 % If you modify MATPOWER, or any covered work, to interface with 0061 % other modules (such as MATLAB code and MEX-files) available in a 0062 % MATLAB(R) or comparable environment containing parts covered 0063 % under other licensing terms, the licensors of MATPOWER grant 0064 % you additional permission to convey the resulting work. 0065 0066 %%----- initialize ----- 0067 %% define named indices into data matrices 0068 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ... 0069 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ... 0070 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen; 0071 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost; 0072 0073 %% default args 0074 if nargin < 11 0075 il = []; 0076 end 0077 0078 %% unpack data 0079 mpc = get_mpc(om); 0080 [baseMVA, bus, gen, branch, gencost] = ... 0081 deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch, mpc.gencost); 0082 cp = get_cost_params(om); 0083 [N, Cw, H, dd, rh, kk, mm] = deal(cp.N, cp.Cw, cp.H, cp.dd, ... 0084 cp.rh, cp.kk, cp.mm); 0085 vv = get_idx(om); 0086 0087 %% problem dimensions 0088 nb = size(bus, 1); %% number of buses 0089 nl = size(branch, 1); %% number of branches 0090 ng = size(gen, 1); %% number of dispatchable injections 0091 ny = getN(om, 'var', 'y'); %% number of piece-wise linear costs 0092 nxyz = length(x); %% total number of control vars of all types 0093 0094 %% set default constrained lines 0095 if isempty(il) 0096 il = (1:nl); %% all lines have limits by default 0097 end 0098 nl2 = length(il); %% number of constrained lines 0099 0100 %% grab Pg & Qg 0101 Pg = x(vv.i1.Pg:vv.iN.Pg); %% active generation in p.u. 0102 Qg = x(vv.i1.Qg:vv.iN.Qg); %% reactive generation in p.u. 0103 0104 %% put Pg & Qg back in gen 0105 gen(:, PG) = Pg * baseMVA; %% active generation in MW 0106 gen(:, QG) = Qg * baseMVA; %% reactive generation in MVAr 0107 0108 %%----- evaluate objective function ----- 0109 %% polynomial cost of P and Q 0110 % use totcost only on polynomial cost; in the minimization problem 0111 % formulation, pwl cost is the sum of the y variables. 0112 ipol = find(gencost(:, MODEL) == POLYNOMIAL); %% poly MW and MVAr costs 0113 xx = [ gen(:, PG); gen(:, QG)]; 0114 if ~isempty(ipol) 0115 f = sum( totcost(gencost(ipol, :), xx(ipol)) ); %% cost of poly P or Q 0116 else 0117 f = 0; 0118 end 0119 0120 %% piecewise linear cost of P and Q 0121 if ny > 0 0122 ccost = full(sparse(ones(1,ny), vv.i1.y:vv.iN.y, ones(1,ny), 1, nxyz)); 0123 f = f + ccost * x; 0124 else 0125 ccost = zeros(1, nxyz); 0126 end 0127 0128 %% generalized cost term 0129 if ~isempty(N) 0130 nw = size(N, 1); 0131 r = N * x - rh; %% Nx - rhat 0132 iLT = find(r < -kk); %% below dead zone 0133 iEQ = find(r == 0 & kk == 0); %% dead zone doesn't exist 0134 iGT = find(r > kk); %% above dead zone 0135 iND = [iLT; iEQ; iGT]; %% rows that are Not in the Dead region 0136 iL = find(dd == 1); %% rows using linear function 0137 iQ = find(dd == 2); %% rows using quadratic function 0138 LL = sparse(iL, iL, 1, nw, nw); 0139 QQ = sparse(iQ, iQ, 1, nw, nw); 0140 kbar = sparse(iND, iND, [ ones(length(iLT), 1); 0141 zeros(length(iEQ), 1); 0142 -ones(length(iGT), 1)], nw, nw) * kk; 0143 rr = r + kbar; %% apply non-dead zone shift 0144 M = sparse(iND, iND, mm(iND), nw, nw); %% dead zone or scale 0145 diagrr = sparse(1:nw, 1:nw, rr, nw, nw); 0146 0147 %% linear rows multiplied by rr(i), quadratic rows by rr(i)^2 0148 w = M * (LL + QQ * diagrr) * rr; 0149 0150 f = f + (w' * H * w) / 2 + Cw' * w; 0151 end 0152 0153 %%----- evaluate cost gradient ----- 0154 %% index ranges 0155 iPg = vv.i1.Pg:vv.iN.Pg; 0156 iQg = vv.i1.Qg:vv.iN.Qg; 0157 0158 %% polynomial cost of P and Q 0159 df_dPgQg = zeros(2*ng, 1); %% w.r.t p.u. Pg and Qg 0160 df_dPgQg(ipol) = baseMVA * polycost(gencost(ipol, :), xx(ipol), 1); 0161 df = zeros(nxyz, 1); 0162 df(iPg) = df_dPgQg(1:ng); 0163 df(iQg) = df_dPgQg((1:ng) + ng); 0164 0165 %% piecewise linear cost of P and Q 0166 df = df + ccost'; % As in MINOS, the linear cost row is additive wrt 0167 % any nonlinear cost. 0168 0169 %% generalized cost term 0170 if ~isempty(N) 0171 HwC = H * w + Cw; 0172 AA = N' * M * (LL + 2 * QQ * diagrr); 0173 df = df + AA * HwC; 0174 0175 %% numerical check 0176 if 0 %% 1 to check, 0 to skip check 0177 ddff = zeros(size(df)); 0178 step = 1e-7; 0179 tol = 1e-3; 0180 for k = 1:length(x) 0181 xx = x; 0182 xx(k) = xx(k) + step; 0183 ddff(k) = (fun_copf(xx, om, Ybus, Yf, Yt, Afeq, bfeq, Af, bf, mpopt, il) - f) / step; 0184 end 0185 if max(abs(ddff - df)) > tol 0186 idx = find(abs(ddff - df) == max(abs(ddff - df))); 0187 fprintf('\nMismatch in gradient\n'); 0188 fprintf('idx df(num) df diff\n'); 0189 fprintf('%4d%16g%16g%16g\n', [ 1:length(df); ddff'; df'; abs(ddff - df)' ]); 0190 fprintf('MAX\n'); 0191 fprintf('%4d%16g%16g%16g\n', [ idx'; ddff(idx)'; df(idx)'; abs(ddff(idx) - df(idx))' ]); 0192 fprintf('\n'); 0193 end 0194 end %% numeric check 0195 end 0196 0197 %% ---- evaluate cost Hessian ----- 0198 if nargout > 2 0199 pcost = gencost(1:ng, :); 0200 if size(gencost, 1) > ng 0201 qcost = gencost(ng+1:2*ng, :); 0202 else 0203 qcost = []; 0204 end 0205 0206 %% polynomial generator costs 0207 d2f_dPg2 = sparse(ng, 1); %% w.r.t. p.u. Pg 0208 d2f_dQg2 = sparse(ng, 1); %% w.r.t. p.u. Qg 0209 ipolp = find(pcost(:, MODEL) == POLYNOMIAL); 0210 d2f_dPg2(ipolp) = baseMVA^2 * polycost(pcost(ipolp, :), Pg(ipolp)*baseMVA, 2); 0211 if ~isempty(qcost) %% Qg is not free 0212 ipolq = find(qcost(:, MODEL) == POLYNOMIAL); 0213 d2f_dQg2(ipolq) = baseMVA^2 * polycost(qcost(ipolq, :), Qg(ipolq)*baseMVA, 2); 0214 end 0215 i = (pgbas:qgend)'; 0216 d2f = sparse(i, i, [d2f_dPg2; d2f_dQg2], nxyz, nxyz); 0217 0218 %% generalized cost 0219 if ~isempty(N) 0220 d2f = d2f + AA * H * AA' + 2 * N' * M * QQ * sparse(1:nw, 1:nw, HwC, nw, nw) * N; 0221 end 0222 end 0223 0224 %%----- evaluate partials of constraints ----- 0225 %% reconstruct V 0226 Va = x(vv.i1.Va:vv.iN.Va); 0227 Vm = x(vv.i1.Vm:vv.iN.Vm); 0228 V = Vm .* exp(1j * Va); 0229 0230 %% compute partials of injected bus powers 0231 [dSbus_dVm, dSbus_dVa] = dSbus_dV(Ybus, V); %% w.r.t. V 0232 neg_Cg = sparse(gen(:, GEN_BUS), 1:ng, -1, nb, ng); %% Pbus w.r.t. Pg 0233 %% Qbus w.r.t. Qg 0234 0235 %% compute partials of Flows w.r.t. V 0236 if mpopt(24) == 2 %% current 0237 [dFf_dVa, dFf_dVm, dFt_dVa, dFt_dVm, Ff, Ft] = dIbr_dV(branch(il,:), Yf, Yt, V); 0238 else %% power 0239 [dFf_dVa, dFf_dVm, dFt_dVa, dFt_dVm, Ff, Ft] = dSbr_dV(branch(il,:), Yf, Yt, V); 0240 end 0241 if mpopt(24) == 1 %% real part of flow (active power) 0242 dFf_dVa = real(dFf_dVa); 0243 dFf_dVm = real(dFf_dVm); 0244 dFt_dVa = real(dFt_dVa); 0245 dFt_dVm = real(dFt_dVm); 0246 Ff = real(Ff); 0247 Ft = real(Ft); 0248 end 0249 0250 %% squared magnitude of flow (of complex power or current, or real power) 0251 [df_dVa, df_dVm, dt_dVa, dt_dVm] = ... 0252 dAbr_dV(dFf_dVa, dFf_dVm, dFt_dVa, dFt_dVm, Ff, Ft); 0253 0254 %% index ranges 0255 iVa = vv.i1.Va:vv.iN.Va; 0256 iVm = vv.i1.Vm:vv.iN.Vm; 0257 iPg = vv.i1.Pg:vv.iN.Pg; 0258 iQg = vv.i1.Qg:vv.iN.Qg; 0259 nleq = size(Afeq, 1); 0260 nliq = size(Af, 1); 0261 0262 %% construct Jacobian of constraints 0263 dg = sparse(nxyz, 2*nb+2*nl2+nleq+nliq); 0264 %% equality constraints 0265 dg([iVa iVm], 1:2*nb) = [ 0266 real(dSbus_dVa), real(dSbus_dVm); %% P mismatch w.r.t Va, Vm 0267 imag(dSbus_dVa), imag(dSbus_dVm); %% Q mismatch w.r.t Va, Vm 0268 ]'; 0269 dg(iPg, 1:nb) = neg_Cg'; %% P mismatch w.r.t Pg 0270 dg(iQg, (1:nb) + nb) = neg_Cg'; %% Q mismatch w.r.t Qg 0271 dg(:, (1:nleq)+2*nb) = Afeq'; %% linear equalities 0272 %% inequality constraints 0273 dg([iVa iVm], (1:2*nl2)+2*nb+nleq) = [ 0274 df_dVa, df_dVm; %% "from" flow limit 0275 dt_dVa, dt_dVm; %% "to" flow limit 0276 ]'; 0277 dg(:, (1:nliq)+2*nb+2*nl2+nleq) = Af'; %% linear inequalities 0278 0279 %% use non-sparse matrices 0280 dg = full(dg);