0001 function Lxx = opf_hessfcn(x, lambda, cost_mult, om, Ybus, Yf, Yt, mpopt, il)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
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
0054 if isempty(cost_mult)
0055 cost_mult = 1;
0056 end
0057
0058
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
0068 nb = size(bus, 1);
0069 nl = size(branch, 1);
0070 ng = size(gen, 1);
0071 nxyz = length(x);
0072
0073
0074 if nargin < 8
0075 il = (1:nl);
0076 end
0077 nl2 = length(il);
0078
0079
0080 Pg = x(vv.i1.Pg:vv.iN.Pg);
0081 Qg = x(vv.i1.Qg:vv.iN.Qg);
0082
0083
0084 gen(:, PG) = Pg * baseMVA;
0085 gen(:, QG) = Qg * baseMVA;
0086
0087
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
0101 d2f_dPg2 = sparse(ng, 1);
0102 d2f_dQg2 = sparse(ng, 1);
0103 ipolp = find(pcost(:, MODEL) == POLYNOMIAL);
0104 d2f_dPg2(ipolp) = baseMVA^2 * polycost(pcost(ipolp, :), Pg(ipolp)*baseMVA, 2);
0105 if ~isempty(qcost)
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
0113 if ~isempty(N)
0114 nw = size(N, 1);
0115 r = N * x - rh;
0116 iLT = find(r < -kk);
0117 iEQ = find(r == 0 & kk == 0);
0118 iGT = find(r > kk);
0119 iND = [iLT; iEQ; iGT];
0120 iL = find(dd == 1);
0121 iQ = find(dd == 2);
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;
0128 M = sparse(iND, iND, mm(iND), nw, nw);
0129 diagrr = sparse(1:nw, 1:nw, rr, nw, nw);
0130
0131
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
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
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
0156 nmu = length(lambda.ineqnonlin) / 2;
0157 if nmu
0158 muF = lambda.ineqnonlin(1:nmu);
0159 muT = lambda.ineqnonlin((1:nmu)+nmu);
0160 else
0161 muF = zeros(0,1);
0162 muT = zeros(0,1);
0163 end
0164 if upper(mpopt.opf.flow_lim(1)) == 'I'
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);
0170 t = branch(il, T_BUS);
0171 Cf = sparse(1:nl2, f, ones(nl2, 1), nl2, nb);
0172 Ct = sparse(1:nl2, t, ones(nl2, 1), nl2, nb);
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'
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
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
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
0200 [fp, dfp] = opf_costfcn(xp, om);
0201 [fm, dfm] = opf_costfcn(xm, om);
0202
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;