Home > matpower4.0 > grad_copf.m

grad_copf

PURPOSE ^

------------------------------ deprecated ------------------------------

SYNOPSIS ^

function [df, dg, d2f] = grad_copf(x, om, Ybus, Yf, Yt, Afeq, bfeq, Af, bf, mpopt, il)

DESCRIPTION ^

------------------------------  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.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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);

Generated on Mon 26-Jan-2015 14:56:45 by m2html © 2005