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