0001 function t_hessian(quiet)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013 if nargin < 1
0014 quiet = 0;
0015 end
0016
0017 t_begin(96, quiet);
0018
0019 casefile = 'case30';
0020
0021
0022 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0023 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0024 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0025 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0026 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0027
0028
0029 mpopt = mpoption('verbose', 0, 'out.all', 0);
0030 [baseMVA, bus, gen, branch, success, et] = runpf(casefile, mpopt);
0031
0032
0033 [i2e, bus, gen, branch] = ext2int(bus, gen, branch);
0034 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0035 Sbus = makeSbus(baseMVA, bus, gen);
0036 Vm = bus(:, VM);
0037 Va = bus(:, VA) * pi/180;
0038 V = Vm .* exp(1j * Va);
0039 Vr = real(V);
0040 Vi = imag(V);
0041 f = branch(:, F_BUS);
0042 t = branch(:, T_BUS);
0043 nl = length(f);
0044 nb = length(V);
0045 Cf = sparse(1:nl, f, ones(nl, 1), nl, nb);
0046 Ct = sparse(1:nl, t, ones(nl, 1), nl, nb);
0047 pert = 1e-8;
0048
0049
0050 for vcart = 0:1
0051
0052 if vcart
0053 coord = 'cartesian';
0054 vv = {'rr', 'ri', 'ir', 'ii'};
0055 V1p = (Vr*ones(1,nb) + pert*eye(nb,nb)) + 1j * Vi * ones(1,nb);
0056 V2p = Vr*ones(1,nb) + 1j * (Vi*ones(1,nb) + pert*eye(nb,nb));
0057 else
0058 coord = 'polar';
0059 vv = {'aa', 'av', 'va', 'vv'};
0060 V1p = (Vm*ones(1,nb)) .* (exp(1j * (Va*ones(1,nb) + pert*eye(nb,nb))));
0061 V2p = (Vm*ones(1,nb) + pert*eye(nb,nb)) .* (exp(1j * Va) * ones(1,nb));
0062 end
0063
0064
0065
0066 t = ' - d2Imis_dV2 (complex current injections)';
0067 lam = 10 * rand(nb, 1);
0068 num_H11 = zeros(nb, nb);
0069 num_H12 = zeros(nb, nb);
0070 num_H21 = zeros(nb, nb);
0071 num_H22 = zeros(nb, nb);
0072 [dImis_dV1, dImis_dV2] = dImis_dV(Sbus, Ybus, V, vcart);
0073 [H11, H12, H21, H22] = d2Imis_dV2(Sbus, Ybus, V, lam, vcart);
0074 for i = 1:nb
0075 V1p = V;
0076 V2p = V;
0077 if vcart
0078 V1p(i) = (Vr(i) + pert) + 1j * Vi(i);
0079 V2p(i) = Vr(i) + 1j * (Vi(i) + pert);
0080 else
0081 V1p(i) = Vm(i) * exp(1j * (Va(i) + pert));
0082 V2p(i) = (Vm(i) + pert) * exp(1j * Va(i));
0083 end
0084 [dImis_dV1_1p, dImis_dV2_1p] = dImis_dV(Sbus, Ybus, V1p, vcart);
0085 num_H11(:, i) = (dImis_dV1_1p - dImis_dV1).' * lam / pert;
0086 num_H21(:, i) = (dImis_dV2_1p - dImis_dV2).' * lam / pert;
0087
0088 [dImis_dV1_2p, dImis_dV2_2p] = dImis_dV(Sbus, Ybus, V2p, vcart);
0089 num_H12(:, i) = (dImis_dV1_2p - dImis_dV1).' * lam / pert;
0090 num_H22(:, i) = (dImis_dV2_2p - dImis_dV2).' * lam / pert;
0091 end
0092
0093 t_is(full(H11), num_H11, 4, sprintf('%s - H%s%s', coord, vv{1}, t));
0094 t_is(full(H12), num_H12, 4, sprintf('%s - H%s%s', coord, vv{2}, t));
0095 t_is(full(H21), num_H21, 4, sprintf('%s - H%s%s', coord, vv{3}, t));
0096 t_is(full(H22), num_H22, 4, sprintf('%s - H%s%s', coord, vv{4}, t));
0097
0098
0099 t = ' - d2Sbus_dV2 (complex power injections)';
0100 lam = 10 * rand(nb, 1);
0101 num_H11 = zeros(nb, nb);
0102 num_H12 = zeros(nb, nb);
0103 num_H21 = zeros(nb, nb);
0104 num_H22 = zeros(nb, nb);
0105 [dSbus_dV1, dSbus_dV2] = dSbus_dV(Ybus, V, vcart);
0106 [H11, H12, H21, H22] = d2Sbus_dV2(Ybus, V, lam, vcart);
0107 for i = 1:nb
0108 V1p = V;
0109 V2p = V;
0110 if vcart
0111 V1p(i) = (Vr(i) + pert) + 1j * Vi(i);
0112 V2p(i) = Vr(i) + 1j * (Vi(i) + pert);
0113 else
0114 V1p(i) = Vm(i) * exp(1j * (Va(i) + pert));
0115 V2p(i) = (Vm(i) + pert) * exp(1j * Va(i));
0116 end
0117 [dSbus_dV1_1p, dSbus_dV2_1p] = dSbus_dV(Ybus, V1p, vcart);
0118 num_H11(:, i) = (dSbus_dV1_1p - dSbus_dV1).' * lam / pert;
0119 num_H21(:, i) = (dSbus_dV2_1p - dSbus_dV2).' * lam / pert;
0120
0121 [dSbus_dV1_2p, dSbus_dV2_2p] = dSbus_dV(Ybus, V2p, vcart);
0122 num_H12(:, i) = (dSbus_dV1_2p - dSbus_dV1).' * lam / pert;
0123 num_H22(:, i) = (dSbus_dV2_2p - dSbus_dV2).' * lam / pert;
0124 end
0125
0126 t_is(full(H11), num_H11, 4, sprintf('%s - H%s%s', coord, vv{1}, t));
0127 t_is(full(H12), num_H12, 4, sprintf('%s - H%s%s', coord, vv{2}, t));
0128 t_is(full(H21), num_H21, 4, sprintf('%s - H%s%s', coord, vv{3}, t));
0129 t_is(full(H22), num_H22, 4, sprintf('%s - H%s%s', coord, vv{4}, t));
0130
0131
0132 t = ' - d2Ibr_dV2 (complex currents)';
0133 lam = 10 * rand(nl, 1);
0134
0135 num_Gf11 = zeros(nb, nb);
0136 num_Gf12 = zeros(nb, nb);
0137 num_Gf21 = zeros(nb, nb);
0138 num_Gf22 = zeros(nb, nb);
0139 num_Gt11 = zeros(nb, nb);
0140 num_Gt12 = zeros(nb, nb);
0141 num_Gt21 = zeros(nb, nb);
0142 num_Gt22 = zeros(nb, nb);
0143 [dIf_dV1, dIf_dV2, dIt_dV1, dIt_dV2, If, It] = dIbr_dV(branch, Yf, Yt, V, vcart);
0144 [Gf11, Gf12, Gf21, Gf22] = d2Ibr_dV2(Yf, V, lam, vcart);
0145 [Gt11, Gt12, Gt21, Gt22] = d2Ibr_dV2(Yt, V, lam, vcart);
0146 for i = 1:nb
0147 V1p = V;
0148 V2p = V;
0149 if vcart
0150 V1p(i) = (Vr(i) + pert) + 1j * Vi(i);
0151 V2p(i) = Vr(i) + 1j * (Vi(i) + pert);
0152 else
0153 V1p(i) = Vm(i) * exp(1j * (Va(i) + pert));
0154 V2p(i) = (Vm(i) + pert) * exp(1j * Va(i));
0155 end
0156 [dIf_dV1_1p, dIf_dV2_1p, dIt_dV1_1p, dIt_dV2_1p, If_1p, It_1p] = ...
0157 dIbr_dV(branch, Yf, Yt, V1p, vcart);
0158 num_Gf11(:, i) = (dIf_dV1_1p - dIf_dV1).' * lam / pert;
0159 num_Gf21(:, i) = (dIf_dV2_1p - dIf_dV2).' * lam / pert;
0160 num_Gt11(:, i) = (dIt_dV1_1p - dIt_dV1).' * lam / pert;
0161 num_Gt21(:, i) = (dIt_dV2_1p - dIt_dV2).' * lam / pert;
0162
0163 [dIf_dV1_2p, dIf_dV2_2p, dIt_dV1_2p, dIt_dV2_2p, If_2p, It_2p] = ...
0164 dIbr_dV(branch, Yf, Yt, V2p, vcart);
0165 num_Gf12(:, i) = (dIf_dV1_2p - dIf_dV1).' * lam / pert;
0166 num_Gf22(:, i) = (dIf_dV2_2p - dIf_dV2).' * lam / pert;
0167 num_Gt12(:, i) = (dIt_dV1_2p - dIt_dV1).' * lam / pert;
0168 num_Gt22(:, i) = (dIt_dV2_2p - dIt_dV2).' * lam / pert;
0169 end
0170
0171 t_is(full(Gf11), num_Gf11, 4, sprintf('%s - Gf%s%s', coord, vv{1}, t));
0172 t_is(full(Gf12), num_Gf12, 4, sprintf('%s - Gf%s%s', coord, vv{2}, t));
0173 t_is(full(Gf21), num_Gf21, 4, sprintf('%s - Gf%s%s', coord, vv{3}, t));
0174 t_is(full(Gf22), num_Gf22, 4, sprintf('%s - Gf%s%s', coord, vv{4}, t));
0175
0176 t_is(full(Gt11), num_Gt11, 4, sprintf('%s - Gt%s%s', coord, vv{1}, t));
0177 t_is(full(Gt12), num_Gt12, 4, sprintf('%s - Gt%s%s', coord, vv{2}, t));
0178 t_is(full(Gt21), num_Gt21, 4, sprintf('%s - Gt%s%s', coord, vv{3}, t));
0179 t_is(full(Gt22), num_Gt22, 4, sprintf('%s - Gt%s%s', coord, vv{4}, t));
0180
0181
0182 t = ' - d2Sbr_dV2 (complex power flows)';
0183 lam = 10 * rand(nl, 1);
0184
0185 num_Gf11 = zeros(nb, nb);
0186 num_Gf12 = zeros(nb, nb);
0187 num_Gf21 = zeros(nb, nb);
0188 num_Gf22 = zeros(nb, nb);
0189 num_Gt11 = zeros(nb, nb);
0190 num_Gt12 = zeros(nb, nb);
0191 num_Gt21 = zeros(nb, nb);
0192 num_Gt22 = zeros(nb, nb);
0193 [dSf_dV1, dSf_dV2, dSt_dV1, dSt_dV2, Sf, St] = dSbr_dV(branch, Yf, Yt, V, vcart);
0194 [Gf11, Gf12, Gf21, Gf22] = d2Sbr_dV2(Cf, Yf, V, lam, vcart);
0195 [Gt11, Gt12, Gt21, Gt22] = d2Sbr_dV2(Ct, Yt, V, lam, vcart);
0196 for i = 1:nb
0197 V1p = V;
0198 V2p = V;
0199 if vcart
0200 V1p(i) = (Vr(i) + pert) + 1j * Vi(i);
0201 V2p(i) = Vr(i) + 1j * (Vi(i) + pert);
0202 else
0203 V1p(i) = Vm(i) * exp(1j * (Va(i) + pert));
0204 V2p(i) = (Vm(i) + pert) * exp(1j * Va(i));
0205 end
0206 [dSf_dV1_1p, dSf_dV2_1p, dSt_dV1_1p, dSt_dV2_1p, Sf_1p, St_1p] = ...
0207 dSbr_dV(branch, Yf, Yt, V1p, vcart);
0208 num_Gf11(:, i) = (dSf_dV1_1p - dSf_dV1).' * lam / pert;
0209 num_Gf21(:, i) = (dSf_dV2_1p - dSf_dV2).' * lam / pert;
0210 num_Gt11(:, i) = (dSt_dV1_1p - dSt_dV1).' * lam / pert;
0211 num_Gt21(:, i) = (dSt_dV2_1p - dSt_dV2).' * lam / pert;
0212
0213 [dSf_dV1_2p, dSf_dV2_2p, dSt_dV1_2p, dSt_dV2_2p, Sf_2p, St_2p] = ...
0214 dSbr_dV(branch, Yf, Yt, V2p, vcart);
0215 num_Gf12(:, i) = (dSf_dV1_2p - dSf_dV1).' * lam / pert;
0216 num_Gf22(:, i) = (dSf_dV2_2p - dSf_dV2).' * lam / pert;
0217 num_Gt12(:, i) = (dSt_dV1_2p - dSt_dV1).' * lam / pert;
0218 num_Gt22(:, i) = (dSt_dV2_2p - dSt_dV2).' * lam / pert;
0219 end
0220
0221 t_is(full(Gf11), num_Gf11, 4, sprintf('%s - Gf%s%s', coord, vv{1}, t));
0222 t_is(full(Gf12), num_Gf12, 4, sprintf('%s - Gf%s%s', coord, vv{2}, t));
0223 t_is(full(Gf21), num_Gf21, 4, sprintf('%s - Gf%s%s', coord, vv{3}, t));
0224 t_is(full(Gf22), num_Gf22, 4, sprintf('%s - Gf%s%s', coord, vv{4}, t));
0225
0226 t_is(full(Gt11), num_Gt11, 4, sprintf('%s - Gt%s%s', coord, vv{1}, t));
0227 t_is(full(Gt12), num_Gt12, 4, sprintf('%s - Gt%s%s', coord, vv{2}, t));
0228 t_is(full(Gt21), num_Gt21, 4, sprintf('%s - Gt%s%s', coord, vv{3}, t));
0229 t_is(full(Gt22), num_Gt22, 4, sprintf('%s - Gt%s%s', coord, vv{4}, t));
0230
0231
0232 t = ' - d2Abr_dV2 (squared current magnitudes)';
0233 lam = 10 * rand(nl, 1);
0234
0235 num_Gf11 = zeros(nb, nb);
0236 num_Gf12 = zeros(nb, nb);
0237 num_Gf21 = zeros(nb, nb);
0238 num_Gf22 = zeros(nb, nb);
0239 num_Gt11 = zeros(nb, nb);
0240 num_Gt12 = zeros(nb, nb);
0241 num_Gt21 = zeros(nb, nb);
0242 num_Gt22 = zeros(nb, nb);
0243 d2If_dV2 = @(V, mu)d2Ibr_dV2(Yf, V, mu, vcart);
0244 d2It_dV2 = @(V, mu)d2Ibr_dV2(Yt, V, mu, vcart);
0245 [dIf_dV1, dIf_dV2, dIt_dV1, dIt_dV2, If, It] = dIbr_dV(branch, Yf, Yt, V, vcart);
0246 [dAf_dV1, dAf_dV2, dAt_dV1, dAt_dV2] = ...
0247 dAbr_dV(dIf_dV1, dIf_dV2, dIt_dV1, dIt_dV2, If, It);
0248 [Gf11, Gf12, Gf21, Gf22] = d2Abr_dV2(d2If_dV2, dIf_dV1, dIf_dV2, If, V, lam);
0249 [Gt11, Gt12, Gt21, Gt22] = d2Abr_dV2(d2It_dV2, dIt_dV1, dIt_dV2, It, V, lam);
0250 for i = 1:nb
0251 V1p = V;
0252 V2p = V;
0253 if vcart
0254 V1p(i) = (Vr(i) + pert) + 1j * Vi(i);
0255 V2p(i) = Vr(i) + 1j * (Vi(i) + pert);
0256 else
0257 V1p(i) = Vm(i) * exp(1j * (Va(i) + pert));
0258 V2p(i) = (Vm(i) + pert) * exp(1j * Va(i));
0259 end
0260 [dIf_dV1_1p, dIf_dV2_1p, dIt_dV1_1p, dIt_dV2_1p, If_1p, It_1p] = ...
0261 dIbr_dV(branch, Yf, Yt, V1p, vcart);
0262 [dAf_dV1_1p, dAf_dV2_1p, dAt_dV1_1p, dAt_dV2_1p] = ...
0263 dAbr_dV(dIf_dV1_1p, dIf_dV2_1p, dIt_dV1_1p, dIt_dV2_1p, If_1p, It_1p);
0264 num_Gf11(:, i) = (dAf_dV1_1p - dAf_dV1).' * lam / pert;
0265 num_Gf21(:, i) = (dAf_dV2_1p - dAf_dV2).' * lam / pert;
0266 num_Gt11(:, i) = (dAt_dV1_1p - dAt_dV1).' * lam / pert;
0267 num_Gt21(:, i) = (dAt_dV2_1p - dAt_dV2).' * lam / pert;
0268
0269 [dIf_dV1_2p, dIf_dV2_2p, dIt_dV1_2p, dIt_dV2_2p, If_2p, It_2p] = ...
0270 dIbr_dV(branch, Yf, Yt, V2p, vcart);
0271 [dAf_dV1_2p, dAf_dV2_2p, dAt_dV1_2p, dAt_dV2_2p] = ...
0272 dAbr_dV(dIf_dV1_2p, dIf_dV2_2p, dIt_dV1_2p, dIt_dV2_2p, If_2p, It_2p);
0273 num_Gf12(:, i) = (dAf_dV1_2p - dAf_dV1).' * lam / pert;
0274 num_Gf22(:, i) = (dAf_dV2_2p - dAf_dV2).' * lam / pert;
0275 num_Gt12(:, i) = (dAt_dV1_2p - dAt_dV1).' * lam / pert;
0276 num_Gt22(:, i) = (dAt_dV2_2p - dAt_dV2).' * lam / pert;
0277 end
0278
0279 t_is(full(Gf11), num_Gf11, 2.5, sprintf('%s - Gf%s%s', coord, vv{1}, t));
0280 t_is(full(Gf12), num_Gf12, 2.5, sprintf('%s - Gf%s%s', coord, vv{2}, t));
0281 t_is(full(Gf21), num_Gf21, 2.5, sprintf('%s - Gf%s%s', coord, vv{3}, t));
0282 t_is(full(Gf22), num_Gf22, 2.5, sprintf('%s - Gf%s%s', coord, vv{4}, t));
0283
0284 t_is(full(Gt11), num_Gt11, 2.5, sprintf('%s - Gt%s%s', coord, vv{1}, t));
0285 t_is(full(Gt12), num_Gt12, 2.5, sprintf('%s - Gt%s%s', coord, vv{2}, t));
0286 t_is(full(Gt21), num_Gt21, 2.5, sprintf('%s - Gt%s%s', coord, vv{3}, t));
0287 t_is(full(Gt22), num_Gt22, 2.5, sprintf('%s - Gt%s%s', coord, vv{4}, t));
0288
0289
0290 t = ' - d2Abr_dV2 (squared apparent power flows)';
0291 lam = 10 * rand(nl, 1);
0292
0293 num_Gf11 = zeros(nb, nb);
0294 num_Gf12 = zeros(nb, nb);
0295 num_Gf21 = zeros(nb, nb);
0296 num_Gf22 = zeros(nb, nb);
0297 num_Gt11 = zeros(nb, nb);
0298 num_Gt12 = zeros(nb, nb);
0299 num_Gt21 = zeros(nb, nb);
0300 num_Gt22 = zeros(nb, nb);
0301 d2Sf_dV2 = @(V, mu)d2Sbr_dV2(Cf, Yf, V, mu, vcart);
0302 d2St_dV2 = @(V, mu)d2Sbr_dV2(Ct, Yt, V, mu, vcart);
0303 [dSf_dV1, dSf_dV2, dSt_dV1, dSt_dV2, Sf, St] = dSbr_dV(branch, Yf, Yt, V, vcart);
0304 [dAf_dV1, dAf_dV2, dAt_dV1, dAt_dV2] = ...
0305 dAbr_dV(dSf_dV1, dSf_dV2, dSt_dV1, dSt_dV2, Sf, St);
0306 [Gf11, Gf12, Gf21, Gf22] = d2Abr_dV2(d2Sf_dV2, dSf_dV1, dSf_dV2, Sf, V, lam);
0307 [Gt11, Gt12, Gt21, Gt22] = d2Abr_dV2(d2St_dV2, dSt_dV1, dSt_dV2, St, V, lam);
0308 for i = 1:nb
0309 V1p = V;
0310 V2p = V;
0311 if vcart
0312 V1p(i) = (Vr(i) + pert) + 1j * Vi(i);
0313 V2p(i) = Vr(i) + 1j * (Vi(i) + pert);
0314 else
0315 V1p(i) = Vm(i) * exp(1j * (Va(i) + pert));
0316 V2p(i) = (Vm(i) + pert) * exp(1j * Va(i));
0317 end
0318 [dSf_dV1_1p, dSf_dV2_1p, dSt_dV1_1p, dSt_dV2_1p, Sf_1p, St_1p] = ...
0319 dSbr_dV(branch, Yf, Yt, V1p, vcart);
0320 [dAf_dV1_1p, dAf_dV2_1p, dAt_dV1_1p, dAt_dV2_1p] = ...
0321 dAbr_dV(dSf_dV1_1p, dSf_dV2_1p, dSt_dV1_1p, dSt_dV2_1p, Sf_1p, St_1p);
0322 num_Gf11(:, i) = (dAf_dV1_1p - dAf_dV1).' * lam / pert;
0323 num_Gf21(:, i) = (dAf_dV2_1p - dAf_dV2).' * lam / pert;
0324 num_Gt11(:, i) = (dAt_dV1_1p - dAt_dV1).' * lam / pert;
0325 num_Gt21(:, i) = (dAt_dV2_1p - dAt_dV2).' * lam / pert;
0326
0327 [dSf_dV1_2p, dSf_dV2_2p, dSt_dV1_2p, dSt_dV2_2p, Sf_2p, St_2p] = ...
0328 dSbr_dV(branch, Yf, Yt, V2p, vcart);
0329 [dAf_dV1_2p, dAf_dV2_2p, dAt_dV1_2p, dAt_dV2_2p] = ...
0330 dAbr_dV(dSf_dV1_2p, dSf_dV2_2p, dSt_dV1_2p, dSt_dV2_2p, Sf_2p, St_2p);
0331 num_Gf12(:, i) = (dAf_dV1_2p - dAf_dV1).' * lam / pert;
0332 num_Gf22(:, i) = (dAf_dV2_2p - dAf_dV2).' * lam / pert;
0333 num_Gt12(:, i) = (dAt_dV1_2p - dAt_dV1).' * lam / pert;
0334 num_Gt22(:, i) = (dAt_dV2_2p - dAt_dV2).' * lam / pert;
0335 end
0336
0337 t_is(full(Gf11), num_Gf11, 2.5, sprintf('%s - Gf%s%s', coord, vv{1}, t));
0338 t_is(full(Gf12), num_Gf12, 2.5, sprintf('%s - Gf%s%s', coord, vv{2}, t));
0339 t_is(full(Gf21), num_Gf21, 2.5, sprintf('%s - Gf%s%s', coord, vv{3}, t));
0340 t_is(full(Gf22), num_Gf22, 2.5, sprintf('%s - Gf%s%s', coord, vv{4}, t));
0341
0342 t_is(full(Gt11), num_Gt11, 2.5, sprintf('%s - Gt%s%s', coord, vv{1}, t));
0343 t_is(full(Gt12), num_Gt12, 2.5, sprintf('%s - Gt%s%s', coord, vv{2}, t));
0344 t_is(full(Gt21), num_Gt21, 2.5, sprintf('%s - Gt%s%s', coord, vv{3}, t));
0345 t_is(full(Gt22), num_Gt22, 2.5, sprintf('%s - Gt%s%s', coord, vv{4}, t));
0346
0347
0348 t = ' - d2Abr_dV2 (squared real power flows)';
0349 lam = 10 * rand(nl, 1);
0350
0351 num_Gf11 = zeros(nb, nb);
0352 num_Gf12 = zeros(nb, nb);
0353 num_Gf21 = zeros(nb, nb);
0354 num_Gf22 = zeros(nb, nb);
0355 num_Gt11 = zeros(nb, nb);
0356 num_Gt12 = zeros(nb, nb);
0357 num_Gt21 = zeros(nb, nb);
0358 num_Gt22 = zeros(nb, nb);
0359 d2Sf_dV2 = @(V, mu)d2Sbr_dV2(Cf, Yf, V, mu, vcart);
0360 d2St_dV2 = @(V, mu)d2Sbr_dV2(Ct, Yt, V, mu, vcart);
0361 [dSf_dV1, dSf_dV2, dSt_dV1, dSt_dV2, Sf, St] = dSbr_dV(branch, Yf, Yt, V, vcart);
0362 [dAf_dV1, dAf_dV2, dAt_dV1, dAt_dV2] = ...
0363 dAbr_dV(real(dSf_dV1), real(dSf_dV2), real(dSt_dV1), real(dSt_dV2), real(Sf), real(St));
0364 [Gf11, Gf12, Gf21, Gf22] = d2Abr_dV2(d2Sf_dV2, real(dSf_dV1), real(dSf_dV2), real(Sf), V, lam);
0365 [Gt11, Gt12, Gt21, Gt22] = d2Abr_dV2(d2St_dV2, real(dSt_dV1), real(dSt_dV2), real(St), V, lam);
0366 for i = 1:nb
0367 V1p = V;
0368 V2p = V;
0369 if vcart
0370 V1p(i) = (Vr(i) + pert) + 1j * Vi(i);
0371 V2p(i) = Vr(i) + 1j * (Vi(i) + pert);
0372 else
0373 V1p(i) = Vm(i) * exp(1j * (Va(i) + pert));
0374 V2p(i) = (Vm(i) + pert) * exp(1j * Va(i));
0375 end
0376 [dSf_dV1_1p, dSf_dV2_1p, dSt_dV1_1p, dSt_dV2_1p, Sf_1p, St_1p] = ...
0377 dSbr_dV(branch, Yf, Yt, V1p, vcart);
0378 [dAf_dV1_1p, dAf_dV2_1p, dAt_dV1_1p, dAt_dV2_1p] = ...
0379 dAbr_dV(real(dSf_dV1_1p), real(dSf_dV2_1p), real(dSt_dV1_1p), real(dSt_dV2_1p), real(Sf_1p), real(St_1p));
0380 num_Gf11(:, i) = (dAf_dV1_1p - dAf_dV1).' * lam / pert;
0381 num_Gf21(:, i) = (dAf_dV2_1p - dAf_dV2).' * lam / pert;
0382 num_Gt11(:, i) = (dAt_dV1_1p - dAt_dV1).' * lam / pert;
0383 num_Gt21(:, i) = (dAt_dV2_1p - dAt_dV2).' * lam / pert;
0384
0385 [dSf_dV1_2p, dSf_dV2_2p, dSt_dV1_2p, dSt_dV2_2p, Sf_2p, St_2p] = ...
0386 dSbr_dV(branch, Yf, Yt, V2p, vcart);
0387 [dAf_dV1_2p, dAf_dV2_2p, dAt_dV1_2p, dAt_dV2_2p] = ...
0388 dAbr_dV(real(dSf_dV1_2p), real(dSf_dV2_2p), real(dSt_dV1_2p), real(dSt_dV2_2p), real(Sf_2p), real(St_2p));
0389 num_Gf12(:, i) = (dAf_dV1_2p - dAf_dV1).' * lam / pert;
0390 num_Gf22(:, i) = (dAf_dV2_2p - dAf_dV2).' * lam / pert;
0391 num_Gt12(:, i) = (dAt_dV1_2p - dAt_dV1).' * lam / pert;
0392 num_Gt22(:, i) = (dAt_dV2_2p - dAt_dV2).' * lam / pert;
0393 end
0394
0395 t_is(full(Gf11), num_Gf11, 2.5, sprintf('%s - Gf%s%s', coord, vv{1}, t));
0396 t_is(full(Gf12), num_Gf12, 2.5, sprintf('%s - Gf%s%s', coord, vv{2}, t));
0397 t_is(full(Gf21), num_Gf21, 2.5, sprintf('%s - Gf%s%s', coord, vv{3}, t));
0398 t_is(full(Gf22), num_Gf22, 2.5, sprintf('%s - Gf%s%s', coord, vv{4}, t));
0399
0400 t_is(full(Gt11), num_Gt11, 2.5, sprintf('%s - Gt%s%s', coord, vv{1}, t));
0401 t_is(full(Gt12), num_Gt12, 2.5, sprintf('%s - Gt%s%s', coord, vv{2}, t));
0402 t_is(full(Gt21), num_Gt21, 2.5, sprintf('%s - Gt%s%s', coord, vv{3}, t));
0403 t_is(full(Gt22), num_Gt22, 2.5, sprintf('%s - Gt%s%s', coord, vv{4}, t));
0404 end
0405
0406 t_end;