opf_hessfcn

OPF_HESSFCN  Evaluates Hessian of Lagrangian for AC OPF.
   LXX = OPF_HESSFCN(X, LAMBDA, COST_MULT, OM, YBUS, YF, YT, MPOPT, IL)

   Hessian evaluation function for AC optimal power flow, suitable
   for use with MIPS or FMINCON's interior-point algorithm.

   Inputs:
     X : optimization vector
     LAMBDA (struct)
       .eqnonlin : Lagrange multipliers on power balance equations
       .ineqnonlin : Kuhn-Tucker multipliers on constrained branch flows
     COST_MULT : (optional) Scale factor to be applied to the cost
          (default = 1).
     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
     MPOPT : MATPOWER options struct
     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:
     LXX : Hessian of the Lagrangian.

   Examples:
       Lxx = opf_hessfcn(x, lambda, cost_mult, om, Ybus, Yf, Yt, mpopt);
       Lxx = opf_hessfcn(x, lambda, cost_mult, om, Ybus, Yf, Yt, mpopt, il);

   See also OPF_COSTFCN, OPF_CONSFCN.

0001 function Lxx = opf_hessfcn(x, lambda, cost_mult, om, Ybus, Yf, Yt, mpopt, il)
0002 %OPF_HESSFCN  Evaluates Hessian of Lagrangian for AC OPF.
0003 %   LXX = OPF_HESSFCN(X, LAMBDA, COST_MULT, OM, YBUS, YF, YT, MPOPT, IL)
0004 %
0005 %   Hessian evaluation function for AC optimal power flow, suitable
0006 %   for use with MIPS or FMINCON's interior-point algorithm.
0007 %
0008 %   Inputs:
0009 %     X : optimization vector
0010 %     LAMBDA (struct)
0011 %       .eqnonlin : Lagrange multipliers on power balance equations
0012 %       .ineqnonlin : Kuhn-Tucker multipliers on constrained branch flows
0013 %     COST_MULT : (optional) Scale factor to be applied to the cost
0014 %          (default = 1).
0015 %     OM : OPF model object
0016 %     YBUS : bus admittance matrix
0017 %     YF : admittance matrix for "from" end of constrained branches
0018 %     YT : admittance matrix for "to" end of constrained branches
0019 %     MPOPT : MATPOWER options struct
0020 %     IL : (optional) vector of branch indices corresponding to
0021 %          branches with flow limits (all others are assumed to be
0022 %          unconstrained). The default is [1:nl] (all branches).
0023 %          YF and YT contain only the rows corresponding to IL.
0024 %
0025 %   Outputs:
0026 %     LXX : Hessian of the Lagrangian.
0027 %
0028 %   Examples:
0029 %       Lxx = opf_hessfcn(x, lambda, cost_mult, om, Ybus, Yf, Yt, mpopt);
0030 %       Lxx = opf_hessfcn(x, lambda, cost_mult, om, Ybus, Yf, Yt, mpopt, il);
0031 %
0032 %   See also OPF_COSTFCN, OPF_CONSFCN.
0033 
0034 %   MATPOWER
0035 %   $Id: opf_hessfcn.m 2328 2014-05-23 18:43:00Z ray $
0036 %   by Ray Zimmerman, PSERC Cornell
0037 %   and Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Autonoma de Manizales
0038 %   Copyright (c) 1996-2010 by Power System Engineering Research Center (PSERC)
0039 %
0040 %   This file is part of MATPOWER.
0041 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0042 %
0043 %   MATPOWER is free software: you can redistribute it and/or modify
0044 %   it under the terms of the GNU General Public License as published
0045 %   by the Free Software Foundation, either version 3 of the License,
0046 %   or (at your option) any later version.
0047 %
0048 %   MATPOWER is distributed in the hope that it will be useful,
0049 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0050 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0051 %   GNU General Public License for more details.
0052 %
0053 %   You should have received a copy of the GNU General Public License
0054 %   along with MATPOWER. If not, see <http://www.gnu.org/licenses/>.
0055 %
0056 %   Additional permission under GNU GPL version 3 section 7
0057 %
0058 %   If you modify MATPOWER, or any covered work, to interface with
0059 %   other modules (such as MATLAB code and MEX-files) available in a
0060 %   MATLAB(R) or comparable environment containing parts covered
0061 %   under other licensing terms, the licensors of MATPOWER grant
0062 %   you additional permission to convey the resulting work.
0063 
0064 %%----- initialize -----
0065 %% define named indices into data matrices
0066 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0067     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0068     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0069 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0070     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0071     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0072 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
0073 
0074 %% default args
0075 if isempty(cost_mult)
0076     cost_mult = 1;
0077 end
0078 
0079 %% unpack data
0080 mpc = get_mpc(om);
0081 [baseMVA, bus, gen, branch, gencost] = ...
0082     deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch, mpc.gencost);
0083 cp = get_cost_params(om);
0084 [N, Cw, H, dd, rh, kk, mm] = deal(cp.N, cp.Cw, cp.H, cp.dd, ...
0085                                     cp.rh, cp.kk, cp.mm);
0086 vv = get_idx(om);
0087 
0088 %% unpack needed parameters
0089 nb = size(bus, 1);          %% number of buses
0090 nl = size(branch, 1);       %% number of branches
0091 ng = size(gen, 1);          %% number of dispatchable injections
0092 nxyz = length(x);           %% total number of control vars of all types
0093 
0094 %% set default constrained lines
0095 if nargin < 8
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 %% reconstruct V
0109 Va = zeros(nb, 1);
0110 Va = x(vv.i1.Va:vv.iN.Va);
0111 Vm = x(vv.i1.Vm:vv.iN.Vm);
0112 V = Vm .* exp(1j * Va);
0113 nxtra = nxyz - 2*nb;
0114 pcost = gencost(1:ng, :);
0115 if size(gencost, 1) > ng
0116     qcost = gencost(ng+1:2*ng, :);
0117 else
0118     qcost = [];
0119 end
0120 
0121 %% ----- evaluate d2f -----
0122 d2f_dPg2 = sparse(ng, 1);               %% w.r.t. p.u. Pg
0123 d2f_dQg2 = sparse(ng, 1);               %% w.r.t. p.u. Qg
0124 ipolp = find(pcost(:, MODEL) == POLYNOMIAL);
0125 d2f_dPg2(ipolp) = baseMVA^2 * polycost(pcost(ipolp, :), Pg(ipolp)*baseMVA, 2);
0126 if ~isempty(qcost)          %% Qg is not free
0127     ipolq = find(qcost(:, MODEL) == POLYNOMIAL);
0128     d2f_dQg2(ipolq) = baseMVA^2 * polycost(qcost(ipolq, :), Qg(ipolq)*baseMVA, 2);
0129 end
0130 i = [vv.i1.Pg:vv.iN.Pg vv.i1.Qg:vv.iN.Qg]';
0131 d2f = sparse(i, i, [d2f_dPg2; d2f_dQg2], nxyz, nxyz);
0132 
0133 %% generalized cost
0134 if ~isempty(N)
0135     nw = size(N, 1);
0136     r = N * x - rh;                 %% Nx - rhat
0137     iLT = find(r < -kk);            %% below dead zone
0138     iEQ = find(r == 0 & kk == 0);   %% dead zone doesn't exist
0139     iGT = find(r > kk);             %% above dead zone
0140     iND = [iLT; iEQ; iGT];          %% rows that are Not in the Dead region
0141     iL = find(dd == 1);             %% rows using linear function
0142     iQ = find(dd == 2);             %% rows using quadratic function
0143     LL = sparse(iL, iL, 1, nw, nw);
0144     QQ = sparse(iQ, iQ, 1, nw, nw);
0145     kbar = sparse(iND, iND, [   ones(length(iLT), 1);
0146                                 zeros(length(iEQ), 1);
0147                                 -ones(length(iGT), 1)], nw, nw) * kk;
0148     rr = r + kbar;                  %% apply non-dead zone shift
0149     M = sparse(iND, iND, mm(iND), nw, nw);  %% dead zone or scale
0150     diagrr = sparse(1:nw, 1:nw, rr, nw, nw);
0151     
0152     %% linear rows multiplied by rr(i), quadratic rows by rr(i)^2
0153     w = M * (LL + QQ * diagrr) * rr;
0154     HwC = H * w + Cw;
0155     AA = N' * M * (LL + 2 * QQ * diagrr);
0156     d2f = d2f + AA * H * AA' + 2 * N' * M * QQ * sparse(1:nw, 1:nw, HwC, nw, nw) * N;
0157 end
0158 d2f = d2f * cost_mult;
0159 
0160 %%----- evaluate Hessian of power balance constraints -----
0161 nlam = length(lambda.eqnonlin) / 2;
0162 lamP = lambda.eqnonlin(1:nlam);
0163 lamQ = lambda.eqnonlin((1:nlam)+nlam);
0164 [Gpaa, Gpav, Gpva, Gpvv] = d2Sbus_dV2(Ybus, V, lamP);
0165 [Gqaa, Gqav, Gqva, Gqvv] = d2Sbus_dV2(Ybus, V, lamQ);
0166 d2G = [
0167     real([Gpaa Gpav; Gpva Gpvv]) + imag([Gqaa Gqav; Gqva Gqvv]) sparse(2*nb, nxtra);
0168     sparse(nxtra, 2*nb + nxtra)
0169 ];
0170 
0171 %%----- evaluate Hessian of flow constraints -----
0172 nmu = length(lambda.ineqnonlin) / 2;
0173 if nmu
0174     muF = lambda.ineqnonlin(1:nmu);
0175     muT = lambda.ineqnonlin((1:nmu)+nmu);
0176 else    %% keep dimensions of empty matrices/vectors compatible
0177     muF = zeros(0,1);   %% (required to avoid problems when using Knitro
0178     muT = zeros(0,1);   %%  on cases with all lines unconstrained)
0179 end
0180 if upper(mpopt.opf.flow_lim(1)) == 'I'      %% current
0181     [dIf_dVa, dIf_dVm, dIt_dVa, dIt_dVm, If, It] = dIbr_dV(branch(il,:), Yf, Yt, V);
0182     [Hfaa, Hfav, Hfva, Hfvv] = d2AIbr_dV2(dIf_dVa, dIf_dVm, If, Yf, V, muF);
0183     [Htaa, Htav, Htva, Htvv] = d2AIbr_dV2(dIt_dVa, dIt_dVm, It, Yt, V, muT);
0184 else
0185   f = branch(il, F_BUS);    %% list of "from" buses
0186   t = branch(il, T_BUS);    %% list of "to" buses
0187   Cf = sparse(1:nl2, f, ones(nl2, 1), nl2, nb);     %% connection matrix for line & from buses
0188   Ct = sparse(1:nl2, t, ones(nl2, 1), nl2, nb);     %% connection matrix for line & to buses
0189   [dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St] = dSbr_dV(branch(il,:), Yf, Yt, V);
0190   if upper(mpopt.opf.flow_lim(1)) == 'P'    %% real power
0191     [Hfaa, Hfav, Hfva, Hfvv] = d2ASbr_dV2(real(dSf_dVa), real(dSf_dVm), real(Sf), Cf, Yf, V, muF);
0192     [Htaa, Htav, Htva, Htvv] = d2ASbr_dV2(real(dSt_dVa), real(dSt_dVm), real(St), Ct, Yt, V, muT);
0193   else                                      %% apparent power
0194     [Hfaa, Hfav, Hfva, Hfvv] = d2ASbr_dV2(dSf_dVa, dSf_dVm, Sf, Cf, Yf, V, muF);
0195     [Htaa, Htav, Htva, Htvv] = d2ASbr_dV2(dSt_dVa, dSt_dVm, St, Ct, Yt, V, muT);
0196   end
0197 end
0198 d2H = [
0199     [Hfaa Hfav; Hfva Hfvv] + [Htaa Htav; Htva Htvv] sparse(2*nb, nxtra);
0200     sparse(nxtra, 2*nb + nxtra)
0201 ];
0202 
0203 %%-----  do numerical check using (central) finite differences  -----
0204 if 0
0205     nx = length(x);
0206     step = 1e-5;
0207     num_d2f = sparse(nx, nx);
0208     num_d2G = sparse(nx, nx);
0209     num_d2H = sparse(nx, nx);
0210     for i = 1:nx
0211         xp = x;
0212         xm = x;
0213         xp(i) = x(i) + step/2;
0214         xm(i) = x(i) - step/2;
0215         % evaluate cost & gradients
0216         [fp, dfp] = opf_costfcn(xp, om);
0217         [fm, dfm] = opf_costfcn(xm, om);
0218         % evaluate constraints & gradients
0219         [Hp, Gp, dHp, dGp] = opf_consfcn(xp, om, Ybus, Yf, Yt, mpopt, il);
0220         [Hm, Gm, dHm, dGm] = opf_consfcn(xm, om, Ybus, Yf, Yt, mpopt, il);
0221         num_d2f(:, i) = cost_mult * (dfp - dfm) / step;
0222         num_d2G(:, i) = (dGp - dGm) * lambda.eqnonlin   / step;
0223         num_d2H(:, i) = (dHp - dHm) * lambda.ineqnonlin / step;
0224     end
0225     d2f_err = full(max(max(abs(d2f - num_d2f))));
0226     d2G_err = full(max(max(abs(d2G - num_d2G))));
0227     d2H_err = full(max(max(abs(d2H - num_d2H))));
0228     if d2f_err > 1e-6
0229         fprintf('Max difference in d2f: %g\n', d2f_err);
0230     end
0231     if d2G_err > 1e-5
0232         fprintf('Max difference in d2G: %g\n', d2G_err);
0233     end
0234     if d2H_err > 1e-6
0235         fprintf('Max difference in d2H: %g\n', d2H_err);
0236     end
0237 end
0238 
0239 Lxx = d2f + d2G + d2H;