t_hessian

T_HESSIAN  Numerical tests of 2nd derivative code.

0001 function t_hessian(quiet)
0002 %T_HESSIAN  Numerical tests of 2nd derivative code.
0003 
0004 %   MATPOWER
0005 %   Copyright (c) 2004-2016, Power Systems Engineering Research Center (PSERC)
0006 %   by Ray Zimmerman, PSERC Cornell
0007 %   and Baljinnyam Sereeter, Delft University of Technology
0008 %
0009 %   This file is part of MATPOWER.
0010 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0011 %   See https://matpower.org for more info.
0012 
0013 if nargin < 1
0014     quiet = 0;
0015 end
0016 
0017 t_begin(96, quiet);
0018 
0019 casefile = 'case30';
0020 
0021 %% define named indices into bus, gen, branch matrices
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 %% run powerflow to get solved case
0029 mpopt = mpoption('verbose', 0, 'out.all', 0);
0030 [baseMVA, bus, gen, branch, success, et] = runpf(casefile, mpopt);
0031 
0032 %% switch to internal bus numbering and build admittance matrices
0033 [i2e, bus, gen, branch] = ext2int(bus, gen, branch);
0034 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0035 Sbus = makeSbus(baseMVA, bus, gen);  %% net injected power in p.u.
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);       %% list of "from" buses
0042 t = branch(:, T_BUS);       %% list of "to" buses
0043 nl = length(f);
0044 nb = length(V);
0045 Cf = sparse(1:nl, f, ones(nl, 1), nl, nb);      %% connection matrix for line & from buses
0046 Ct = sparse(1:nl, t, ones(nl, 1), nl, nb);      %% connection matrix for line & to buses
0047 pert = 1e-8;
0048 
0049 %%-----  run tests for polar, then cartesian, coordinate voltages  -----
0050 for vcart = 0:1
0051     %%-----  create perturbed voltages  -----
0052     if vcart        %% cartesian coordinate voltages (V1=Vr, V2=Vi)
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            %% polar coordinate voltages (V1=Va, V2=Vm)
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     %%-----  check d2Imis_dV2 code  -----
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);       %% perturb Vr
0079             V2p(i) = Vr(i) + 1j * (Vi(i) + pert);       %% perturb Vi
0080         else
0081             V1p(i) = Vm(i) * exp(1j * (Va(i) + pert));  %% perturb Va
0082             V2p(i) = (Vm(i) + pert) * exp(1j * Va(i));  %% perturb Vm
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     %%-----  check d2Sbus_dV2 code  -----
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);       %% perturb Vr
0112             V2p(i) = Vr(i) + 1j * (Vi(i) + pert);       %% perturb Vi
0113         else
0114             V1p(i) = Vm(i) * exp(1j * (Va(i) + pert));  %% perturb Va
0115             V2p(i) = (Vm(i) + pert) * exp(1j * Va(i));  %% perturb Vm
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     %%-----  check d2Ibr_dV2 code  -----
0132     t = ' - d2Ibr_dV2 (complex currents)';
0133     lam = 10 * rand(nl, 1);
0134     % lam = [1; zeros(nl-1, 1)];
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);       %% perturb Vr
0151             V2p(i) = Vr(i) + 1j * (Vi(i) + pert);       %% perturb Vi
0152         else
0153             V1p(i) = Vm(i) * exp(1j * (Va(i) + pert));  %% perturb Va
0154             V2p(i) = (Vm(i) + pert) * exp(1j * Va(i));  %% perturb Vm
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     %%-----  check d2Sbr_dV2 code  -----
0182     t = ' - d2Sbr_dV2 (complex power flows)';
0183     lam = 10 * rand(nl, 1);
0184     % lam = [1; zeros(nl-1, 1)];
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);       %% perturb Vr
0201             V2p(i) = Vr(i) + 1j * (Vi(i) + pert);       %% perturb Vi
0202         else
0203             V1p(i) = Vm(i) * exp(1j * (Va(i) + pert));  %% perturb Va
0204             V2p(i) = (Vm(i) + pert) * exp(1j * Va(i));  %% perturb Vm
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     %%-----  check d2Abr_dV2 code  -----
0232     t = ' - d2Abr_dV2 (squared current magnitudes)';
0233     lam = 10 * rand(nl, 1);
0234     % lam = [1; zeros(nl-1, 1)];
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);       %% perturb Vr
0255             V2p(i) = Vr(i) + 1j * (Vi(i) + pert);       %% perturb Vi
0256         else
0257             V1p(i) = Vm(i) * exp(1j * (Va(i) + pert));  %% perturb Va
0258             V2p(i) = (Vm(i) + pert) * exp(1j * Va(i));  %% perturb Vm
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     %%-----  check d2Abr_dV2 code  -----
0290     t = ' - d2Abr_dV2 (squared apparent power flows)';
0291     lam = 10 * rand(nl, 1);
0292     % lam = [1; zeros(nl-1, 1)];
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);       %% perturb Vr
0313             V2p(i) = Vr(i) + 1j * (Vi(i) + pert);       %% perturb Vi
0314         else
0315             V1p(i) = Vm(i) * exp(1j * (Va(i) + pert));  %% perturb Va
0316             V2p(i) = (Vm(i) + pert) * exp(1j * Va(i));  %% perturb Vm
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     %%-----  check d2Abr_dV2 code  -----
0348     t = ' - d2Abr_dV2 (squared real power flows)';
0349     lam = 10 * rand(nl, 1);
0350     % lam = [1; zeros(nl-1, 1)];
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);       %% perturb Vr
0371             V2p(i) = Vr(i) + 1j * (Vi(i) + pert);       %% perturb Vi
0372         else
0373             V1p(i) = Vm(i) * exp(1j * (Va(i) + pert));  %% perturb Va
0374             V2p(i) = (Vm(i) + pert) * exp(1j * Va(i));  %% perturb Vm
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;