0001 function t_hessian(quiet)
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 if nargin < 1
0034 quiet = 0;
0035 end
0036
0037 t_begin(44, quiet);
0038
0039 casefile = 'case30';
0040
0041
0042 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0043 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0044 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0045 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0046 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0047
0048
0049 mpopt = mpoption('verbose', 0, 'out.all', 0);
0050 [baseMVA, bus, gen, branch, success, et] = runpf(casefile, mpopt);
0051
0052
0053 [i2e, bus, gen, branch] = ext2int(bus, gen, branch);
0054 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0055 Vm = bus(:, VM);
0056 Va = bus(:, VA) * pi/180;
0057 V = Vm .* exp(1j * Va);
0058 f = branch(:, F_BUS);
0059 t = branch(:, T_BUS);
0060 nl = length(f);
0061 nb = length(V);
0062 Cf = sparse(1:nl, f, ones(nl, 1), nl, nb);
0063 Ct = sparse(1:nl, t, ones(nl, 1), nl, nb);
0064 pert = 1e-8;
0065
0066
0067 t = ' - d2Sbus_dV2 (complex power injections)';
0068 lam = 10 * rand(nb, 1);
0069 num_Haa = zeros(nb, nb);
0070 num_Hav = zeros(nb, nb);
0071 num_Hva = zeros(nb, nb);
0072 num_Hvv = zeros(nb, nb);
0073 [dSbus_dVm, dSbus_dVa] = dSbus_dV(Ybus, V);
0074 [Haa, Hav, Hva, Hvv] = d2Sbus_dV2(Ybus, V, lam);
0075 for i = 1:nb
0076 Vap = V;
0077 Vap(i) = Vm(i) * exp(1j * (Va(i) + pert));
0078 [dSbus_dVm_ap, dSbus_dVa_ap] = dSbus_dV(Ybus, Vap);
0079 num_Haa(:, i) = (dSbus_dVa_ap - dSbus_dVa).' * lam / pert;
0080 num_Hva(:, i) = (dSbus_dVm_ap - dSbus_dVm).' * lam / pert;
0081
0082 Vmp = V;
0083 Vmp(i) = (Vm(i) + pert) * exp(1j * Va(i));
0084 [dSbus_dVm_mp, dSbus_dVa_mp] = dSbus_dV(Ybus, Vmp);
0085 num_Hav(:, i) = (dSbus_dVa_mp - dSbus_dVa).' * lam / pert;
0086 num_Hvv(:, i) = (dSbus_dVm_mp - dSbus_dVm).' * lam / pert;
0087 end
0088
0089 t_is(full(Haa), num_Haa, 4, ['Haa' t]);
0090 t_is(full(Hav), num_Hav, 4, ['Hav' t]);
0091 t_is(full(Hva), num_Hva, 4, ['Hva' t]);
0092 t_is(full(Hvv), num_Hvv, 4, ['Hvv' t]);
0093
0094
0095 t = ' - d2Sbr_dV2 (complex power flows)';
0096 lam = 10 * rand(nl, 1);
0097
0098 num_Gfaa = zeros(nb, nb);
0099 num_Gfav = zeros(nb, nb);
0100 num_Gfva = zeros(nb, nb);
0101 num_Gfvv = zeros(nb, nb);
0102 num_Gtaa = zeros(nb, nb);
0103 num_Gtav = zeros(nb, nb);
0104 num_Gtva = zeros(nb, nb);
0105 num_Gtvv = zeros(nb, nb);
0106 [dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St] = dSbr_dV(branch, Yf, Yt, V);
0107 [Gfaa, Gfav, Gfva, Gfvv] = d2Sbr_dV2(Cf, Yf, V, lam);
0108 [Gtaa, Gtav, Gtva, Gtvv] = d2Sbr_dV2(Ct, Yt, V, lam);
0109 for i = 1:nb
0110 Vap = V;
0111 Vap(i) = Vm(i) * exp(1j * (Va(i) + pert));
0112 [dSf_dVa_ap, dSf_dVm_ap, dSt_dVa_ap, dSt_dVm_ap, Sf_ap, St_ap] = ...
0113 dSbr_dV(branch, Yf, Yt, Vap);
0114 num_Gfaa(:, i) = (dSf_dVa_ap - dSf_dVa).' * lam / pert;
0115 num_Gfva(:, i) = (dSf_dVm_ap - dSf_dVm).' * lam / pert;
0116 num_Gtaa(:, i) = (dSt_dVa_ap - dSt_dVa).' * lam / pert;
0117 num_Gtva(:, i) = (dSt_dVm_ap - dSt_dVm).' * lam / pert;
0118
0119 Vmp = V;
0120 Vmp(i) = (Vm(i) + pert) * exp(1j * Va(i));
0121 [dSf_dVa_mp, dSf_dVm_mp, dSt_dVa_mp, dSt_dVm_mp, Sf_mp, St_mp] = ...
0122 dSbr_dV(branch, Yf, Yt, Vmp);
0123 num_Gfav(:, i) = (dSf_dVa_mp - dSf_dVa).' * lam / pert;
0124 num_Gfvv(:, i) = (dSf_dVm_mp - dSf_dVm).' * lam / pert;
0125 num_Gtav(:, i) = (dSt_dVa_mp - dSt_dVa).' * lam / pert;
0126 num_Gtvv(:, i) = (dSt_dVm_mp - dSt_dVm).' * lam / pert;
0127 end
0128
0129 t_is(full(Gfaa), num_Gfaa, 4, ['Gfaa' t]);
0130 t_is(full(Gfav), num_Gfav, 4, ['Gfav' t]);
0131 t_is(full(Gfva), num_Gfva, 4, ['Gfva' t]);
0132 t_is(full(Gfvv), num_Gfvv, 4, ['Gfvv' t]);
0133
0134 t_is(full(Gtaa), num_Gtaa, 4, ['Gtaa' t]);
0135 t_is(full(Gtav), num_Gtav, 4, ['Gtav' t]);
0136 t_is(full(Gtva), num_Gtva, 4, ['Gtva' t]);
0137 t_is(full(Gtvv), num_Gtvv, 4, ['Gtvv' t]);
0138
0139
0140 t = ' - d2Ibr_dV2 (complex currents)';
0141 lam = 10 * rand(nl, 1);
0142
0143 num_Gfaa = zeros(nb, nb);
0144 num_Gfav = zeros(nb, nb);
0145 num_Gfva = zeros(nb, nb);
0146 num_Gfvv = zeros(nb, nb);
0147 num_Gtaa = zeros(nb, nb);
0148 num_Gtav = zeros(nb, nb);
0149 num_Gtva = zeros(nb, nb);
0150 num_Gtvv = zeros(nb, nb);
0151 [dIf_dVa, dIf_dVm, dIt_dVa, dIt_dVm, If, It] = dIbr_dV(branch, Yf, Yt, V);
0152 [Gfaa, Gfav, Gfva, Gfvv] = d2Ibr_dV2(Yf, V, lam);
0153 [Gtaa, Gtav, Gtva, Gtvv] = d2Ibr_dV2(Yt, V, lam);
0154 for i = 1:nb
0155 Vap = V;
0156 Vap(i) = Vm(i) * exp(1j * (Va(i) + pert));
0157 [dIf_dVa_ap, dIf_dVm_ap, dIt_dVa_ap, dIt_dVm_ap, If_ap, It_ap] = ...
0158 dIbr_dV(branch, Yf, Yt, Vap);
0159 num_Gfaa(:, i) = (dIf_dVa_ap - dIf_dVa).' * lam / pert;
0160 num_Gfva(:, i) = (dIf_dVm_ap - dIf_dVm).' * lam / pert;
0161 num_Gtaa(:, i) = (dIt_dVa_ap - dIt_dVa).' * lam / pert;
0162 num_Gtva(:, i) = (dIt_dVm_ap - dIt_dVm).' * lam / pert;
0163
0164 Vmp = V;
0165 Vmp(i) = (Vm(i) + pert) * exp(1j * Va(i));
0166 [dIf_dVa_mp, dIf_dVm_mp, dIt_dVa_mp, dIt_dVm_mp, If_mp, It_mp] = ...
0167 dIbr_dV(branch, Yf, Yt, Vmp);
0168 num_Gfav(:, i) = (dIf_dVa_mp - dIf_dVa).' * lam / pert;
0169 num_Gfvv(:, i) = (dIf_dVm_mp - dIf_dVm).' * lam / pert;
0170 num_Gtav(:, i) = (dIt_dVa_mp - dIt_dVa).' * lam / pert;
0171 num_Gtvv(:, i) = (dIt_dVm_mp - dIt_dVm).' * lam / pert;
0172 end
0173
0174 t_is(full(Gfaa), num_Gfaa, 4, ['Gfaa' t]);
0175 t_is(full(Gfav), num_Gfav, 4, ['Gfav' t]);
0176 t_is(full(Gfva), num_Gfva, 4, ['Gfva' t]);
0177 t_is(full(Gfvv), num_Gfvv, 4, ['Gfvv' t]);
0178
0179 t_is(full(Gtaa), num_Gtaa, 4, ['Gtaa' t]);
0180 t_is(full(Gtav), num_Gtav, 4, ['Gtav' t]);
0181 t_is(full(Gtva), num_Gtva, 4, ['Gtva' t]);
0182 t_is(full(Gtvv), num_Gtvv, 4, ['Gtvv' t]);
0183
0184
0185 t = ' - d2ASbr_dV2 (squared apparent power flows)';
0186 lam = 10 * rand(nl, 1);
0187
0188 num_Gfaa = zeros(nb, nb);
0189 num_Gfav = zeros(nb, nb);
0190 num_Gfva = zeros(nb, nb);
0191 num_Gfvv = zeros(nb, nb);
0192 num_Gtaa = zeros(nb, nb);
0193 num_Gtav = zeros(nb, nb);
0194 num_Gtva = zeros(nb, nb);
0195 num_Gtvv = zeros(nb, nb);
0196 [dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St] = dSbr_dV(branch, Yf, Yt, V);
0197 [dAf_dVa, dAf_dVm, dAt_dVa, dAt_dVm] = ...
0198 dAbr_dV(dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St);
0199 [Gfaa, Gfav, Gfva, Gfvv] = d2ASbr_dV2(dSf_dVa, dSf_dVm, Sf, Cf, Yf, V, lam);
0200 [Gtaa, Gtav, Gtva, Gtvv] = d2ASbr_dV2(dSt_dVa, dSt_dVm, St, Ct, Yt, V, lam);
0201 for i = 1:nb
0202 Vap = V;
0203 Vap(i) = Vm(i) * exp(1j * (Va(i) + pert));
0204 [dSf_dVa_ap, dSf_dVm_ap, dSt_dVa_ap, dSt_dVm_ap, Sf_ap, St_ap] = ...
0205 dSbr_dV(branch, Yf, Yt, Vap);
0206 [dAf_dVa_ap, dAf_dVm_ap, dAt_dVa_ap, dAt_dVm_ap] = ...
0207 dAbr_dV(dSf_dVa_ap, dSf_dVm_ap, dSt_dVa_ap, dSt_dVm_ap, Sf_ap, St_ap);
0208 num_Gfaa(:, i) = (dAf_dVa_ap - dAf_dVa).' * lam / pert;
0209 num_Gfva(:, i) = (dAf_dVm_ap - dAf_dVm).' * lam / pert;
0210 num_Gtaa(:, i) = (dAt_dVa_ap - dAt_dVa).' * lam / pert;
0211 num_Gtva(:, i) = (dAt_dVm_ap - dAt_dVm).' * lam / pert;
0212
0213 Vmp = V;
0214 Vmp(i) = (Vm(i) + pert) * exp(1j * Va(i));
0215 [dSf_dVa_mp, dSf_dVm_mp, dSt_dVa_mp, dSt_dVm_mp, Sf_mp, St_mp] = ...
0216 dSbr_dV(branch, Yf, Yt, Vmp);
0217 [dAf_dVa_mp, dAf_dVm_mp, dAt_dVa_mp, dAt_dVm_mp] = ...
0218 dAbr_dV(dSf_dVa_mp, dSf_dVm_mp, dSt_dVa_mp, dSt_dVm_mp, Sf_mp, St_mp);
0219 num_Gfav(:, i) = (dAf_dVa_mp - dAf_dVa).' * lam / pert;
0220 num_Gfvv(:, i) = (dAf_dVm_mp - dAf_dVm).' * lam / pert;
0221 num_Gtav(:, i) = (dAt_dVa_mp - dAt_dVa).' * lam / pert;
0222 num_Gtvv(:, i) = (dAt_dVm_mp - dAt_dVm).' * lam / pert;
0223 end
0224
0225 t_is(full(Gfaa), num_Gfaa, 2, ['Gfaa' t]);
0226 t_is(full(Gfav), num_Gfav, 2, ['Gfav' t]);
0227 t_is(full(Gfva), num_Gfva, 2, ['Gfva' t]);
0228 t_is(full(Gfvv), num_Gfvv, 2, ['Gfvv' t]);
0229
0230 t_is(full(Gtaa), num_Gtaa, 2, ['Gtaa' t]);
0231 t_is(full(Gtav), num_Gtav, 2, ['Gtav' t]);
0232 t_is(full(Gtva), num_Gtva, 2, ['Gtva' t]);
0233 t_is(full(Gtvv), num_Gtvv, 2, ['Gtvv' t]);
0234
0235
0236 t = ' - d2ASbr_dV2 (squared real power flows)';
0237 lam = 10 * rand(nl, 1);
0238
0239 num_Gfaa = zeros(nb, nb);
0240 num_Gfav = zeros(nb, nb);
0241 num_Gfva = zeros(nb, nb);
0242 num_Gfvv = zeros(nb, nb);
0243 num_Gtaa = zeros(nb, nb);
0244 num_Gtav = zeros(nb, nb);
0245 num_Gtva = zeros(nb, nb);
0246 num_Gtvv = zeros(nb, nb);
0247 [dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St] = dSbr_dV(branch, Yf, Yt, V);
0248 [dAf_dVa, dAf_dVm, dAt_dVa, dAt_dVm] = ...
0249 dAbr_dV(real(dSf_dVa), real(dSf_dVm), real(dSt_dVa), real(dSt_dVm), real(Sf), real(St));
0250 [Gfaa, Gfav, Gfva, Gfvv] = d2ASbr_dV2(real(dSf_dVa), real(dSf_dVm), real(Sf), Cf, Yf, V, lam);
0251 [Gtaa, Gtav, Gtva, Gtvv] = d2ASbr_dV2(real(dSt_dVa), real(dSt_dVm), real(St), Ct, Yt, V, lam);
0252 for i = 1:nb
0253 Vap = V;
0254 Vap(i) = Vm(i) * exp(1j * (Va(i) + pert));
0255 [dSf_dVa_ap, dSf_dVm_ap, dSt_dVa_ap, dSt_dVm_ap, Sf_ap, St_ap] = ...
0256 dSbr_dV(branch, Yf, Yt, Vap);
0257 [dAf_dVa_ap, dAf_dVm_ap, dAt_dVa_ap, dAt_dVm_ap] = ...
0258 dAbr_dV(real(dSf_dVa_ap), real(dSf_dVm_ap), real(dSt_dVa_ap), real(dSt_dVm_ap), real(Sf_ap), real(St_ap));
0259 num_Gfaa(:, i) = (dAf_dVa_ap - dAf_dVa).' * lam / pert;
0260 num_Gfva(:, i) = (dAf_dVm_ap - dAf_dVm).' * lam / pert;
0261 num_Gtaa(:, i) = (dAt_dVa_ap - dAt_dVa).' * lam / pert;
0262 num_Gtva(:, i) = (dAt_dVm_ap - dAt_dVm).' * lam / pert;
0263
0264 Vmp = V;
0265 Vmp(i) = (Vm(i) + pert) * exp(1j * Va(i));
0266 [dSf_dVa_mp, dSf_dVm_mp, dSt_dVa_mp, dSt_dVm_mp, Sf_mp, St_mp] = ...
0267 dSbr_dV(branch, Yf, Yt, Vmp);
0268 [dAf_dVa_mp, dAf_dVm_mp, dAt_dVa_mp, dAt_dVm_mp] = ...
0269 dAbr_dV(real(dSf_dVa_mp), real(dSf_dVm_mp), real(dSt_dVa_mp), real(dSt_dVm_mp), real(Sf_mp), real(St_mp));
0270 num_Gfav(:, i) = (dAf_dVa_mp - dAf_dVa).' * lam / pert;
0271 num_Gfvv(:, i) = (dAf_dVm_mp - dAf_dVm).' * lam / pert;
0272 num_Gtav(:, i) = (dAt_dVa_mp - dAt_dVa).' * lam / pert;
0273 num_Gtvv(:, i) = (dAt_dVm_mp - dAt_dVm).' * lam / pert;
0274 end
0275
0276 t_is(full(Gfaa), num_Gfaa, 2, ['Gfaa' t]);
0277 t_is(full(Gfav), num_Gfav, 2, ['Gfav' t]);
0278 t_is(full(Gfva), num_Gfva, 2, ['Gfva' t]);
0279 t_is(full(Gfvv), num_Gfvv, 2, ['Gfvv' t]);
0280
0281 t_is(full(Gtaa), num_Gtaa, 2, ['Gtaa' t]);
0282 t_is(full(Gtav), num_Gtav, 2, ['Gtav' t]);
0283 t_is(full(Gtva), num_Gtva, 2, ['Gtva' t]);
0284 t_is(full(Gtvv), num_Gtvv, 2, ['Gtvv' t]);
0285
0286
0287 t = ' - d2AIbr_dV2 (squared current magnitudes)';
0288 lam = 10 * rand(nl, 1);
0289
0290 num_Gfaa = zeros(nb, nb);
0291 num_Gfav = zeros(nb, nb);
0292 num_Gfva = zeros(nb, nb);
0293 num_Gfvv = zeros(nb, nb);
0294 num_Gtaa = zeros(nb, nb);
0295 num_Gtav = zeros(nb, nb);
0296 num_Gtva = zeros(nb, nb);
0297 num_Gtvv = zeros(nb, nb);
0298 [dIf_dVa, dIf_dVm, dIt_dVa, dIt_dVm, If, It] = dIbr_dV(branch, Yf, Yt, V);
0299 [dAf_dVa, dAf_dVm, dAt_dVa, dAt_dVm] = ...
0300 dAbr_dV(dIf_dVa, dIf_dVm, dIt_dVa, dIt_dVm, If, It);
0301 [Gfaa, Gfav, Gfva, Gfvv] = d2AIbr_dV2(dIf_dVa, dIf_dVm, If, Yf, V, lam);
0302 [Gtaa, Gtav, Gtva, Gtvv] = d2AIbr_dV2(dIt_dVa, dIt_dVm, It, Yt, V, lam);
0303 for i = 1:nb
0304 Vap = V;
0305 Vap(i) = Vm(i) * exp(1j * (Va(i) + pert));
0306 [dIf_dVa_ap, dIf_dVm_ap, dIt_dVa_ap, dIt_dVm_ap, If_ap, It_ap] = ...
0307 dIbr_dV(branch, Yf, Yt, Vap);
0308 [dAf_dVa_ap, dAf_dVm_ap, dAt_dVa_ap, dAt_dVm_ap] = ...
0309 dAbr_dV(dIf_dVa_ap, dIf_dVm_ap, dIt_dVa_ap, dIt_dVm_ap, If_ap, It_ap);
0310 num_Gfaa(:, i) = (dAf_dVa_ap - dAf_dVa).' * lam / pert;
0311 num_Gfva(:, i) = (dAf_dVm_ap - dAf_dVm).' * lam / pert;
0312 num_Gtaa(:, i) = (dAt_dVa_ap - dAt_dVa).' * lam / pert;
0313 num_Gtva(:, i) = (dAt_dVm_ap - dAt_dVm).' * lam / pert;
0314
0315 Vmp = V;
0316 Vmp(i) = (Vm(i) + pert) * exp(1j * Va(i));
0317 [dIf_dVa_mp, dIf_dVm_mp, dIt_dVa_mp, dIt_dVm_mp, If_mp, It_mp] = ...
0318 dIbr_dV(branch, Yf, Yt, Vmp);
0319 [dAf_dVa_mp, dAf_dVm_mp, dAt_dVa_mp, dAt_dVm_mp] = ...
0320 dAbr_dV(dIf_dVa_mp, dIf_dVm_mp, dIt_dVa_mp, dIt_dVm_mp, If_mp, It_mp);
0321 num_Gfav(:, i) = (dAf_dVa_mp - dAf_dVa).' * lam / pert;
0322 num_Gfvv(:, i) = (dAf_dVm_mp - dAf_dVm).' * lam / pert;
0323 num_Gtav(:, i) = (dAt_dVa_mp - dAt_dVa).' * lam / pert;
0324 num_Gtvv(:, i) = (dAt_dVm_mp - dAt_dVm).' * lam / pert;
0325 end
0326
0327 t_is(full(Gfaa), num_Gfaa, 3, ['Gfaa' t]);
0328 t_is(full(Gfav), num_Gfav, 3, ['Gfav' t]);
0329 t_is(full(Gfva), num_Gfva, 3, ['Gfva' t]);
0330 t_is(full(Gfvv), num_Gfvv, 2, ['Gfvv' t]);
0331
0332 t_is(full(Gtaa), num_Gtaa, 3, ['Gtaa' t]);
0333 t_is(full(Gtav), num_Gtav, 3, ['Gtav' t]);
0334 t_is(full(Gtva), num_Gtva, 3, ['Gtva' t]);
0335 t_is(full(Gtvv), num_Gtvv, 2, ['Gtvv' t]);
0336
0337 t_end;