Home > matpower6.0 > opf_hessfcn.m

opf_hessfcn

PURPOSE ^

OPF_HESSFCN Evaluates Hessian of Lagrangian for AC OPF.

SYNOPSIS ^

function Lxx = opf_hessfcn(x, lambda, cost_mult, om, Ybus, Yf, Yt, mpopt, il)

DESCRIPTION ^

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.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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