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
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
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
0075 if isempty(cost_mult)
0076 cost_mult = 1;
0077 end
0078
0079
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
0089 nb = size(bus, 1);
0090 nl = size(branch, 1);
0091 ng = size(gen, 1);
0092 nxyz = length(x);
0093
0094
0095 if nargin < 8
0096 il = (1:nl);
0097 end
0098 nl2 = length(il);
0099
0100
0101 Pg = x(vv.i1.Pg:vv.iN.Pg);
0102 Qg = x(vv.i1.Qg:vv.iN.Qg);
0103
0104
0105 gen(:, PG) = Pg * baseMVA;
0106 gen(:, QG) = Qg * baseMVA;
0107
0108
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
0122 d2f_dPg2 = sparse(ng, 1);
0123 d2f_dQg2 = sparse(ng, 1);
0124 ipolp = find(pcost(:, MODEL) == POLYNOMIAL);
0125 d2f_dPg2(ipolp) = baseMVA^2 * polycost(pcost(ipolp, :), Pg(ipolp)*baseMVA, 2);
0126 if ~isempty(qcost)
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
0134 if ~isempty(N)
0135 nw = size(N, 1);
0136 r = N * x - rh;
0137 iLT = find(r < -kk);
0138 iEQ = find(r == 0 & kk == 0);
0139 iGT = find(r > kk);
0140 iND = [iLT; iEQ; iGT];
0141 iL = find(dd == 1);
0142 iQ = find(dd == 2);
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;
0149 M = sparse(iND, iND, mm(iND), nw, nw);
0150 diagrr = sparse(1:nw, 1:nw, rr, nw, nw);
0151
0152
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
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
0172 nmu = length(lambda.ineqnonlin) / 2;
0173 if nmu
0174 muF = lambda.ineqnonlin(1:nmu);
0175 muT = lambda.ineqnonlin((1:nmu)+nmu);
0176 else
0177 muF = zeros(0,1);
0178 muT = zeros(0,1);
0179 end
0180 if upper(mpopt.opf.flow_lim(1)) == 'I'
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);
0186 t = branch(il, T_BUS);
0187 Cf = sparse(1:nl2, f, ones(nl2, 1), nl2, nb);
0188 Ct = sparse(1:nl2, t, ones(nl2, 1), nl2, nb);
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'
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
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
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
0216 [fp, dfp] = opf_costfcn(xp, om);
0217 [fm, dfm] = opf_costfcn(xm, om);
0218
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;