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
0046
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
0056 if isempty(cost_mult)
0057 cost_mult = 1;
0058 end
0059
0060
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
0070 nb = size(bus, 1);
0071 nl = size(branch, 1);
0072 ng = size(gen, 1);
0073 nxyz = length(x);
0074
0075
0076 if nargin < 8
0077 il = (1:nl);
0078 end
0079 nl2 = length(il);
0080
0081
0082 Pg = x(vv.i1.Pg:vv.iN.Pg);
0083 Qg = x(vv.i1.Qg:vv.iN.Qg);
0084
0085
0086 gen(:, PG) = Pg * baseMVA;
0087 gen(:, QG) = Qg * baseMVA;
0088
0089
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
0103 d2f_dPg2 = sparse(ng, 1);
0104 d2f_dQg2 = sparse(ng, 1);
0105 ipolp = find(pcost(:, MODEL) == POLYNOMIAL);
0106 d2f_dPg2(ipolp) = baseMVA^2 * polycost(pcost(ipolp, :), Pg(ipolp)*baseMVA, 2);
0107 if ~isempty(qcost)
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
0115 if ~isempty(N)
0116 nw = size(N, 1);
0117 r = N * x - rh;
0118 iLT = find(r < -kk);
0119 iEQ = find(r == 0 & kk == 0);
0120 iGT = find(r > kk);
0121 iND = [iLT; iEQ; iGT];
0122 iL = find(dd == 1);
0123 iQ = find(dd == 2);
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;
0130 M = sparse(iND, iND, mm(iND), nw, nw);
0131 diagrr = sparse(1:nw, 1:nw, rr, nw, nw);
0132
0133
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
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
0153 nmu = length(lambda.ineqnonlin) / 2;
0154 if nmu
0155 muF = lambda.ineqnonlin(1:nmu);
0156 muT = lambda.ineqnonlin((1:nmu)+nmu);
0157 else
0158 muF = zeros(0,1);
0159 muT = zeros(0,1);
0160 end
0161 if upper(mpopt.opf.flow_lim(1)) == 'I'
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);
0167 t = branch(il, T_BUS);
0168 Cf = sparse(1:nl2, f, ones(nl2, 1), nl2, nb);
0169 Ct = sparse(1:nl2, t, ones(nl2, 1), nl2, nb);
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'
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
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
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
0197 [fp, dfp] = opf_costfcn(xp, om);
0198 [fm, dfm] = opf_costfcn(xm, om);
0199
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;