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 %   Copyright (c) 1996-2015 by Power System Engineering Research Center (PSERC)
0036 %   by Ray Zimmerman, PSERC Cornell
0037 %   and Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Autonoma de Manizales
0038 %
0039 %   $Id: opf_hessfcn.m 2644 2015-03-11 19:34:22Z ray $
0040 %
0041 %   This file is part of MATPOWER.
0042 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0043 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0044 
0045 %%----- initialize -----
0046 %% define named indices into data matrices
0047 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0048     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0049     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0050 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0051     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0052     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0053 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
0054 
0055 %% default args
0056 if isempty(cost_mult)
0057     cost_mult = 1;
0058 end
0059 
0060 %% unpack data
0061 mpc = get_mpc(om);
0062 [baseMVA, bus, gen, branch, gencost] = ...
0063     deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch, mpc.gencost);
0064 cp = get_cost_params(om);
0065 [N, Cw, H, dd, rh, kk, mm] = deal(cp.N, cp.Cw, cp.H, cp.dd, ...
0066                                     cp.rh, cp.kk, cp.mm);
0067 vv = get_idx(om);
0068 
0069 %% unpack needed parameters
0070 nb = size(bus, 1);          %% number of buses
0071 nl = size(branch, 1);       %% number of branches
0072 ng = size(gen, 1);          %% number of dispatchable injections
0073 nxyz = length(x);           %% total number of control vars of all types
0074 
0075 %% set default constrained lines
0076 if nargin < 8
0077     il = (1:nl);            %% all lines have limits by default
0078 end
0079 nl2 = length(il);           %% number of constrained lines
0080 
0081 %% grab Pg & Qg
0082 Pg = x(vv.i1.Pg:vv.iN.Pg);  %% active generation in p.u.
0083 Qg = x(vv.i1.Qg:vv.iN.Qg);  %% reactive generation in p.u.
0084 
0085 %% put Pg & Qg back in gen
0086 gen(:, PG) = Pg * baseMVA;  %% active generation in MW
0087 gen(:, QG) = Qg * baseMVA;  %% reactive generation in MVAr
0088  
0089 %% reconstruct V
0090 Va = zeros(nb, 1);
0091 Va = x(vv.i1.Va:vv.iN.Va);
0092 Vm = x(vv.i1.Vm:vv.iN.Vm);
0093 V = Vm .* exp(1j * Va);
0094 nxtra = nxyz - 2*nb;
0095 pcost = gencost(1:ng, :);
0096 if size(gencost, 1) > ng
0097     qcost = gencost(ng+1:2*ng, :);
0098 else
0099     qcost = [];
0100 end
0101 
0102 %% ----- evaluate d2f -----
0103 d2f_dPg2 = sparse(ng, 1);               %% w.r.t. p.u. Pg
0104 d2f_dQg2 = sparse(ng, 1);               %% w.r.t. p.u. Qg
0105 ipolp = find(pcost(:, MODEL) == POLYNOMIAL);
0106 d2f_dPg2(ipolp) = baseMVA^2 * polycost(pcost(ipolp, :), Pg(ipolp)*baseMVA, 2);
0107 if ~isempty(qcost)          %% Qg is not free
0108     ipolq = find(qcost(:, MODEL) == POLYNOMIAL);
0109     d2f_dQg2(ipolq) = baseMVA^2 * polycost(qcost(ipolq, :), Qg(ipolq)*baseMVA, 2);
0110 end
0111 i = [vv.i1.Pg:vv.iN.Pg vv.i1.Qg:vv.iN.Qg]';
0112 d2f = sparse(i, i, [d2f_dPg2; d2f_dQg2], nxyz, nxyz);
0113 
0114 %% generalized cost
0115 if ~isempty(N)
0116     nw = size(N, 1);
0117     r = N * x - rh;                 %% Nx - rhat
0118     iLT = find(r < -kk);            %% below dead zone
0119     iEQ = find(r == 0 & kk == 0);   %% dead zone doesn't exist
0120     iGT = find(r > kk);             %% above dead zone
0121     iND = [iLT; iEQ; iGT];          %% rows that are Not in the Dead region
0122     iL = find(dd == 1);             %% rows using linear function
0123     iQ = find(dd == 2);             %% rows using quadratic function
0124     LL = sparse(iL, iL, 1, nw, nw);
0125     QQ = sparse(iQ, iQ, 1, nw, nw);
0126     kbar = sparse(iND, iND, [   ones(length(iLT), 1);
0127                                 zeros(length(iEQ), 1);
0128                                 -ones(length(iGT), 1)], nw, nw) * kk;
0129     rr = r + kbar;                  %% apply non-dead zone shift
0130     M = sparse(iND, iND, mm(iND), nw, nw);  %% dead zone or scale
0131     diagrr = sparse(1:nw, 1:nw, rr, nw, nw);
0132     
0133     %% linear rows multiplied by rr(i), quadratic rows by rr(i)^2
0134     w = M * (LL + QQ * diagrr) * rr;
0135     HwC = H * w + Cw;
0136     AA = N' * M * (LL + 2 * QQ * diagrr);
0137     d2f = d2f + AA * H * AA' + 2 * N' * M * QQ * sparse(1:nw, 1:nw, HwC, nw, nw) * N;
0138 end
0139 d2f = d2f * cost_mult;
0140 
0141 %%----- evaluate Hessian of power balance constraints -----
0142 nlam = length(lambda.eqnonlin) / 2;
0143 lamP = lambda.eqnonlin(1:nlam);
0144 lamQ = lambda.eqnonlin((1:nlam)+nlam);
0145 [Gpaa, Gpav, Gpva, Gpvv] = d2Sbus_dV2(Ybus, V, lamP);
0146 [Gqaa, Gqav, Gqva, Gqvv] = d2Sbus_dV2(Ybus, V, lamQ);
0147 d2G = [
0148     real([Gpaa Gpav; Gpva Gpvv]) + imag([Gqaa Gqav; Gqva Gqvv]) sparse(2*nb, nxtra);
0149     sparse(nxtra, 2*nb + nxtra)
0150 ];
0151 
0152 %%----- evaluate Hessian of flow constraints -----
0153 nmu = length(lambda.ineqnonlin) / 2;
0154 if nmu
0155     muF = lambda.ineqnonlin(1:nmu);
0156     muT = lambda.ineqnonlin((1:nmu)+nmu);
0157 else    %% keep dimensions of empty matrices/vectors compatible
0158     muF = zeros(0,1);   %% (required to avoid problems when using Knitro
0159     muT = zeros(0,1);   %%  on cases with all lines unconstrained)
0160 end
0161 if upper(mpopt.opf.flow_lim(1)) == 'I'      %% current
0162     [dIf_dVa, dIf_dVm, dIt_dVa, dIt_dVm, If, It] = dIbr_dV(branch(il,:), Yf, Yt, V);
0163     [Hfaa, Hfav, Hfva, Hfvv] = d2AIbr_dV2(dIf_dVa, dIf_dVm, If, Yf, V, muF);
0164     [Htaa, Htav, Htva, Htvv] = d2AIbr_dV2(dIt_dVa, dIt_dVm, It, Yt, V, muT);
0165 else
0166   f = branch(il, F_BUS);    %% list of "from" buses
0167   t = branch(il, T_BUS);    %% list of "to" buses
0168   Cf = sparse(1:nl2, f, ones(nl2, 1), nl2, nb);     %% connection matrix for line & from buses
0169   Ct = sparse(1:nl2, t, ones(nl2, 1), nl2, nb);     %% connection matrix for line & to buses
0170   [dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St] = dSbr_dV(branch(il,:), Yf, Yt, V);
0171   if upper(mpopt.opf.flow_lim(1)) == 'P'    %% real power
0172     [Hfaa, Hfav, Hfva, Hfvv] = d2ASbr_dV2(real(dSf_dVa), real(dSf_dVm), real(Sf), Cf, Yf, V, muF);
0173     [Htaa, Htav, Htva, Htvv] = d2ASbr_dV2(real(dSt_dVa), real(dSt_dVm), real(St), Ct, Yt, V, muT);
0174   else                                      %% apparent power
0175     [Hfaa, Hfav, Hfva, Hfvv] = d2ASbr_dV2(dSf_dVa, dSf_dVm, Sf, Cf, Yf, V, muF);
0176     [Htaa, Htav, Htva, Htvv] = d2ASbr_dV2(dSt_dVa, dSt_dVm, St, Ct, Yt, V, muT);
0177   end
0178 end
0179 d2H = [
0180     [Hfaa Hfav; Hfva Hfvv] + [Htaa Htav; Htva Htvv] sparse(2*nb, nxtra);
0181     sparse(nxtra, 2*nb + nxtra)
0182 ];
0183 
0184 %%-----  do numerical check using (central) finite differences  -----
0185 if 0
0186     nx = length(x);
0187     step = 1e-5;
0188     num_d2f = sparse(nx, nx);
0189     num_d2G = sparse(nx, nx);
0190     num_d2H = sparse(nx, nx);
0191     for i = 1:nx
0192         xp = x;
0193         xm = x;
0194         xp(i) = x(i) + step/2;
0195         xm(i) = x(i) - step/2;
0196         % evaluate cost & gradients
0197         [fp, dfp] = opf_costfcn(xp, om);
0198         [fm, dfm] = opf_costfcn(xm, om);
0199         % evaluate constraints & gradients
0200         [Hp, Gp, dHp, dGp] = opf_consfcn(xp, om, Ybus, Yf, Yt, mpopt, il);
0201         [Hm, Gm, dHm, dGm] = opf_consfcn(xm, om, Ybus, Yf, Yt, mpopt, il);
0202         num_d2f(:, i) = cost_mult * (dfp - dfm) / step;
0203         num_d2G(:, i) = (dGp - dGm) * lambda.eqnonlin   / step;
0204         num_d2H(:, i) = (dHp - dHm) * lambda.ineqnonlin / step;
0205     end
0206     d2f_err = full(max(max(abs(d2f - num_d2f))));
0207     d2G_err = full(max(max(abs(d2G - num_d2G))));
0208     d2H_err = full(max(max(abs(d2H - num_d2H))));
0209     if d2f_err > 1e-6
0210         fprintf('Max difference in d2f: %g\n', d2f_err);
0211     end
0212     if d2G_err > 1e-5
0213         fprintf('Max difference in d2G: %g\n', d2G_err);
0214     end
0215     if d2H_err > 1e-6
0216         fprintf('Max difference in d2H: %g\n', d2H_err);
0217     end
0218 end
0219 
0220 Lxx = d2f + d2G + d2H;