0001 function t_jacobian(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(64, 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 mpc = loadcase(casefile);
0031 [baseMVA, bus, gen, branch, success, et] = runpf(mpc, mpopt);
0032
0033
0034 [i2e, bus, gen, branch] = ext2int(bus, gen, branch);
0035 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0036 Sbus = makeSbus(baseMVA, bus, gen);
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);
0046 t = branch(:, T_BUS);
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
0055 for vcart = 0:1
0056
0057 if vcart
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
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
0070
0071 [dImis_dV1_full, dImis_dV2_full] = dImis_dV(Sbus, Ybus_full, V, vcart);
0072
0073
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
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
0089
0090 [dSbus_dV1_full, dSbus_dV2_full] = dSbus_dV(Ybus_full, V, vcart);
0091
0092
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
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
0107
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
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
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
0133
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
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
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
0170
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
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
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;