Home > matpower7.0 > lib > t > t_jacobian.m

t_jacobian

PURPOSE ^

T_JACOBIAN Numerical tests of partial derivative code.

SYNOPSIS ^

function t_jacobian(quiet)

DESCRIPTION ^

T_JACOBIAN  Numerical tests of partial derivative code.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function t_jacobian(quiet)
0002 %T_JACOBIAN  Numerical tests of partial 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(64, 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 mpc = loadcase(casefile);
0031 [baseMVA, bus, gen, branch, success, et] = runpf(mpc, mpopt);
0032 
0033 %% switch to internal bus numbering and build admittance matrices
0034 [i2e, bus, gen, branch] = ext2int(bus, gen, branch);
0035 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0036 Sbus = makeSbus(baseMVA, bus, gen);  %% net injected power in p.u.
0037 Ybus_full   = full(Ybus);
0038 Yf_full     = full(Yf);
0039 Yt_full     = full(Yt);
0040 Vm = bus(:, VM);
0041 Va = bus(:, VA) * pi/180;
0042 V = Vm .* exp(1j * Va);
0043 Vr = real(V);
0044 Vi = imag(V);
0045 f = branch(:, F_BUS);       %% list of "from" buses
0046 t = branch(:, T_BUS);       %% list of "to" buses
0047 nl = length(f);
0048 nb = length(V);
0049 VV = V * ones(1, nb);
0050 SS = Sbus * ones(1, nb);
0051 pert = 1e-8;
0052 
0053 
0054 %%-----  run tests for polar, then cartesian, coordinate voltages  -----
0055 for vcart = 0:1
0056     %%-----  create perturbed voltages  -----
0057     if vcart        %% cartesian coordinate voltages (V1=Vr, V2=Vi)
0058         coord = 'cartesian';
0059         vv = {'r', 'i'};
0060         V1p = (Vr*ones(1,nb) + pert*eye(nb,nb)) + 1j * Vi * ones(1,nb);
0061         V2p = Vr*ones(1,nb) + 1j * (Vi*ones(1,nb) + pert*eye(nb,nb));
0062     else            %% polar coordinate voltages (V1=Va, V2=Vm)
0063         coord = 'polar';
0064         vv = {'a', 'm'};
0065         V1p = (Vm*ones(1,nb)) .* (exp(1j * (Va*ones(1,nb) + pert*eye(nb,nb))));
0066         V2p = (Vm*ones(1,nb) + pert*eye(nb,nb)) .* (exp(1j * Va) * ones(1,nb));
0067     end
0068 
0069     %%-----  check dImis_dV code  -----
0070     %% full matrices
0071     [dImis_dV1_full, dImis_dV2_full] = dImis_dV(Sbus, Ybus_full, V, vcart);
0072 
0073     %% sparse matrices
0074     [dImis_dV1, dImis_dV2] = dImis_dV(Sbus, Ybus, V, vcart);
0075     dImis_dV1_sp = full(dImis_dV1);
0076     dImis_dV2_sp = full(dImis_dV2);
0077 
0078     %% compute numerically to compare
0079     num_dImis_dV1 = full( ( Ybus * V1p - conj(SS./V1p) - Ybus*VV + conj(SS./VV)) / pert );
0080     num_dImis_dV2 = full( ( Ybus * V2p - conj(SS./V2p) - Ybus*VV + conj(SS./VV)) / pert );
0081 
0082     t_is(dImis_dV1_sp, num_dImis_dV1, 5, sprintf('%s - dImis_dV%s (sparse)', coord, vv{1}));
0083     t_is(dImis_dV2_sp, num_dImis_dV2, 5, sprintf('%s - dImis_dV%s (sparse)', coord, vv{2}));
0084     t_is(dImis_dV1_full, num_dImis_dV1, 5, sprintf('%s - dImis_dV%s (full)', coord, vv{1}));
0085     t_is(dImis_dV2_full, num_dImis_dV2, 5, sprintf('%s - dImis_dV%s (full)', coord, vv{2}));
0086 
0087 
0088     %%-----  check dSbus_dV code  -----
0089     %% full matrices
0090     [dSbus_dV1_full, dSbus_dV2_full] = dSbus_dV(Ybus_full, V, vcart);
0091 
0092     %% sparse matrices
0093     [dSbus_dV1, dSbus_dV2] = dSbus_dV(Ybus, V, vcart);
0094     dSbus_dV1_sp = full(dSbus_dV1);
0095     dSbus_dV2_sp = full(dSbus_dV2);
0096 
0097     %% compute numerically to compare
0098     num_dSbus_dV1 = full( (V1p .* conj(Ybus * V1p) - VV .* conj(Ybus * VV)) / pert );
0099     num_dSbus_dV2 = full( (V2p .* conj(Ybus * V2p) - VV .* conj(Ybus * VV)) / pert );
0100 
0101     t_is(dSbus_dV1_sp, num_dSbus_dV1, 5, sprintf('%s - dSbus_dV%s (sparse)', coord, vv{1}));
0102     t_is(dSbus_dV2_sp, num_dSbus_dV2, 5, sprintf('%s - dSbus_dV%s (sparse)', coord, vv{2}));
0103     t_is(dSbus_dV1_full, num_dSbus_dV1, 5, sprintf('%s - dSbus_dV%s (full)', coord, vv{1}));
0104     t_is(dSbus_dV2_full, num_dSbus_dV2, 5, sprintf('%s - dSbus_dV%s (full)', coord, vv{2}));
0105 
0106     %%-----  check dIbr_dV code  -----
0107     %% full matrices
0108     [dIf_dV1_full, dIf_dV2_full, dIt_dV1_full, dIt_dV2_full, If, It] = dIbr_dV(branch, Yf_full, Yt_full, V, vcart);
0109 
0110     %% sparse matrices
0111     [dIf_dV1, dIf_dV2, dIt_dV1, dIt_dV2, If, It] = dIbr_dV(branch, Yf, Yt, V, vcart);
0112     dIf_dV1_sp = full(dIf_dV1);
0113     dIf_dV2_sp = full(dIf_dV2);
0114     dIt_dV1_sp = full(dIt_dV1);
0115     dIt_dV2_sp = full(dIt_dV2);
0116 
0117     %% compute numerically to compare
0118     num_dIf_dV1 = full( (Yf * V1p - Yf * VV) / pert );
0119     num_dIf_dV2 = full( (Yf * V2p - Yf * VV) / pert );
0120     num_dIt_dV1 = full( (Yt * V1p - Yt * VV) / pert );
0121     num_dIt_dV2 = full( (Yt * V2p - Yt * VV) / pert );
0122 
0123     t_is(dIf_dV1_sp, num_dIf_dV1, 5, sprintf('%s - dIf_dV%s (sparse)', coord, vv{1}));
0124     t_is(dIf_dV2_sp, num_dIf_dV2, 5, sprintf('%s - dIf_dV%s (sparse)', coord, vv{2}));
0125     t_is(dIt_dV1_sp, num_dIt_dV1, 5, sprintf('%s - dIt_dV%s (sparse)', coord, vv{1}));
0126     t_is(dIt_dV2_sp, num_dIt_dV2, 5, sprintf('%s - dIt_dV%s (sparse)', coord, vv{2}));
0127     t_is(dIf_dV1_full, num_dIf_dV1, 5, sprintf('%s - dIf_dV%s (full)', coord, vv{1}));
0128     t_is(dIf_dV2_full, num_dIf_dV2, 5, sprintf('%s - dIf_dV%s (full)', coord, vv{2}));
0129     t_is(dIt_dV1_full, num_dIt_dV1, 5, sprintf('%s - dIt_dV%s (full)', coord, vv{1}));
0130     t_is(dIt_dV2_full, num_dIt_dV2, 5, sprintf('%s - dIt_dV%s (full)', coord, vv{2}));
0131 
0132     %%-----  check dSbr_dV code  -----
0133     %% full matrices
0134     [dSf_dV1_full, dSf_dV2_full, dSt_dV1_full, dSt_dV2_full, Sf, St] = dSbr_dV(branch, Yf_full, Yt_full, V, vcart);
0135 
0136     %% sparse matrices
0137     [dSf_dV1, dSf_dV2, dSt_dV1, dSt_dV2, Sf, St] = dSbr_dV(branch, Yf, Yt, V, vcart);
0138     dSf_dV1_sp = full(dSf_dV1);
0139     dSf_dV2_sp = full(dSf_dV2);
0140     dSt_dV1_sp = full(dSt_dV1);
0141     dSt_dV2_sp = full(dSt_dV2);
0142 
0143     %% compute numerically to compare
0144     V1pf = V1p(f,:);
0145     V2pf = V2p(f,:);
0146     V1pt = V1p(t,:);
0147     V2pt = V2p(t,:);
0148     Sf2 = (V(f)*ones(1,nb)) .* conj(Yf * VV);
0149     St2 = (V(t)*ones(1,nb)) .* conj(Yt * VV);
0150     S1pf = V1pf .* conj(Yf * V1p);
0151     S2pf = V2pf .* conj(Yf * V2p);
0152     S1pt = V1pt .* conj(Yt * V1p);
0153     S2pt = V2pt .* conj(Yt * V2p);
0154 
0155     num_dSf_dV1 = full( (S1pf - Sf2) / pert );
0156     num_dSf_dV2 = full( (S2pf - Sf2) / pert );
0157     num_dSt_dV1 = full( (S1pt - St2) / pert );
0158     num_dSt_dV2 = full( (S2pt - St2) / pert );
0159 
0160     t_is(dSf_dV1_sp, num_dSf_dV1, 5, sprintf('%s - dSf_dV%s (sparse)', coord, vv{1}));
0161     t_is(dSf_dV2_sp, num_dSf_dV2, 5, sprintf('%s - dSf_dV%s (sparse)', coord, vv{2}));
0162     t_is(dSt_dV1_sp, num_dSt_dV1, 5, sprintf('%s - dSt_dV%s (sparse)', coord, vv{1}));
0163     t_is(dSt_dV2_sp, num_dSt_dV2, 5, sprintf('%s - dSt_dV%s (sparse)', coord, vv{2}));
0164     t_is(dSf_dV1_full, num_dSf_dV1, 5, sprintf('%s - dSf_dV%s (full)', coord, vv{1}));
0165     t_is(dSf_dV2_full, num_dSf_dV2, 5, sprintf('%s - dSf_dV%s (full)', coord, vv{2}));
0166     t_is(dSt_dV1_full, num_dSt_dV1, 5, sprintf('%s - dSt_dV%s (full)', coord, vv{1}));
0167     t_is(dSt_dV2_full, num_dSt_dV2, 5, sprintf('%s - dSt_dV%s (full)', coord, vv{2}));
0168 
0169     %%-----  check dAbr_dV code  -----
0170     %% full matrices
0171     [dAf_dV1_full, dAf_dV2_full, dAt_dV1_full, dAt_dV2_full] = ...
0172                             dAbr_dV(dSf_dV1_full, dSf_dV2_full, dSt_dV1_full, dSt_dV2_full, Sf, St);
0173     %% sparse matrices
0174     [dAf_dV1, dAf_dV2, dAt_dV1, dAt_dV2] = ...
0175                             dAbr_dV(dSf_dV1, dSf_dV2, dSt_dV1, dSt_dV2, Sf, St);
0176     dAf_dV1_sp = full(dAf_dV1);
0177     dAf_dV2_sp = full(dAf_dV2);
0178     dAt_dV1_sp = full(dAt_dV1);
0179     dAt_dV2_sp = full(dAt_dV2);
0180 
0181     %% compute numerically to compare
0182     num_dAf_dV1 = full( (abs(S1pf).^2 - abs(Sf2).^2) / pert );
0183     num_dAf_dV2 = full( (abs(S2pf).^2 - abs(Sf2).^2) / pert );
0184     num_dAt_dV1 = full( (abs(S1pt).^2 - abs(St2).^2) / pert );
0185     num_dAt_dV2 = full( (abs(S2pt).^2 - abs(St2).^2) / pert );
0186 
0187     t_is(dAf_dV1_sp, num_dAf_dV1, 4, sprintf('%s - dAf_dV%s (sparse)', coord, vv{1}));
0188     t_is(dAf_dV2_sp, num_dAf_dV2, 4, sprintf('%s - dAf_dV%s (sparse)', coord, vv{2}));
0189     t_is(dAt_dV1_sp, num_dAt_dV1, 4, sprintf('%s - dAt_dV%s (sparse)', coord, vv{1}));
0190     t_is(dAt_dV2_sp, num_dAt_dV2, 4, sprintf('%s - dAt_dV%s (sparse)', coord, vv{2}));
0191     t_is(dAf_dV1_full, num_dAf_dV1, 4, sprintf('%s - dAf_dV%s (full)', coord, vv{1}));
0192     t_is(dAf_dV2_full, num_dAf_dV2, 4, sprintf('%s - dAf_dV%s (full)', coord, vv{2}));
0193     t_is(dAt_dV1_full, num_dAt_dV1, 4, sprintf('%s - dAt_dV%s (full)', coord, vv{1}));
0194     t_is(dAt_dV2_full, num_dAt_dV2, 4, sprintf('%s - dAt_dV%s (full)', coord, vv{2}));
0195 end
0196 
0197 t_end;
Generated on Mon 24-Jun-2019 15:58:45 by m2html © 2005