0001 function t_jacobian(quiet)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033 if nargin < 1
0034 quiet = 0;
0035 end
0036
0037 t_begin(28, quiet);
0038
0039 casefile = 'case30';
0040
0041
0042 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0043 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0044 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0045 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0046 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0047
0048
0049 mpopt = mpoption('verbose', 0, 'out.all', 0);
0050 mpc = loadcase(casefile);
0051 [baseMVA, bus, gen, branch, success, et] = runpf(mpc, mpopt);
0052
0053
0054 [i2e, bus, gen, branch] = ext2int(bus, gen, branch);
0055 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0056 Ybus_full = full(Ybus);
0057 Yf_full = full(Yf);
0058 Yt_full = full(Yt);
0059 Vm = bus(:, VM);
0060 Va = bus(:, VA) * pi/180;
0061 V = Vm .* exp(1j * Va);
0062 f = branch(:, F_BUS);
0063 t = branch(:, T_BUS);
0064 nl = length(f);
0065 nb = length(V);
0066 pert = 1e-8;
0067
0068
0069
0070 [dSbus_dVm_full, dSbus_dVa_full] = dSbus_dV(Ybus_full, V);
0071
0072
0073 [dSbus_dVm, dSbus_dVa] = dSbus_dV(Ybus, V);
0074 dSbus_dVm_sp = full(dSbus_dVm);
0075 dSbus_dVa_sp = full(dSbus_dVa);
0076
0077
0078 Vmp = (Vm*ones(1,nb) + pert*eye(nb,nb)) .* (exp(1j * Va) * ones(1,nb));
0079 Vap = (Vm*ones(1,nb)) .* (exp(1j * (Va*ones(1,nb) + pert*eye(nb,nb))));
0080 num_dSbus_dVm = full( (Vmp .* conj(Ybus * Vmp) - V*ones(1,nb) .* conj(Ybus * V*ones(1,nb))) / pert );
0081 num_dSbus_dVa = full( (Vap .* conj(Ybus * Vap) - V*ones(1,nb) .* conj(Ybus * V*ones(1,nb))) / pert );
0082
0083 t_is(dSbus_dVm_sp, num_dSbus_dVm, 5, 'dSbus_dVm (sparse)');
0084 t_is(dSbus_dVa_sp, num_dSbus_dVa, 5, 'dSbus_dVa (sparse)');
0085 t_is(dSbus_dVm_full, num_dSbus_dVm, 5, 'dSbus_dVm (full)');
0086 t_is(dSbus_dVa_full, num_dSbus_dVa, 5, 'dSbus_dVa (full)');
0087
0088
0089
0090 [dSf_dVa_full, dSf_dVm_full, dSt_dVa_full, dSt_dVm_full, Sf, St] = dSbr_dV(branch, Yf_full, Yt_full, V);
0091
0092
0093 [dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St] = dSbr_dV(branch, Yf, Yt, V);
0094 dSf_dVa_sp = full(dSf_dVa);
0095 dSf_dVm_sp = full(dSf_dVm);
0096 dSt_dVa_sp = full(dSt_dVa);
0097 dSt_dVm_sp = full(dSt_dVm);
0098
0099
0100 Vmpf = Vmp(f,:);
0101 Vapf = Vap(f,:);
0102 Vmpt = Vmp(t,:);
0103 Vapt = Vap(t,:);
0104 Sf2 = (V(f)*ones(1,nb)) .* conj(Yf * V*ones(1,nb));
0105 St2 = (V(t)*ones(1,nb)) .* conj(Yt * V*ones(1,nb));
0106 Smpf = Vmpf .* conj(Yf * Vmp);
0107 Sapf = Vapf .* conj(Yf * Vap);
0108 Smpt = Vmpt .* conj(Yt * Vmp);
0109 Sapt = Vapt .* conj(Yt * Vap);
0110
0111 num_dSf_dVm = full( (Smpf - Sf2) / pert );
0112 num_dSf_dVa = full( (Sapf - Sf2) / pert );
0113 num_dSt_dVm = full( (Smpt - St2) / pert );
0114 num_dSt_dVa = full( (Sapt - St2) / pert );
0115
0116 t_is(dSf_dVm_sp, num_dSf_dVm, 5, 'dSf_dVm (sparse)');
0117 t_is(dSf_dVa_sp, num_dSf_dVa, 5, 'dSf_dVa (sparse)');
0118 t_is(dSt_dVm_sp, num_dSt_dVm, 5, 'dSt_dVm (sparse)');
0119 t_is(dSt_dVa_sp, num_dSt_dVa, 5, 'dSt_dVa (sparse)');
0120 t_is(dSf_dVm_full, num_dSf_dVm, 5, 'dSf_dVm (full)');
0121 t_is(dSf_dVa_full, num_dSf_dVa, 5, 'dSf_dVa (full)');
0122 t_is(dSt_dVm_full, num_dSt_dVm, 5, 'dSt_dVm (full)');
0123 t_is(dSt_dVa_full, num_dSt_dVa, 5, 'dSt_dVa (full)');
0124
0125
0126
0127 [dAf_dVa_full, dAf_dVm_full, dAt_dVa_full, dAt_dVm_full] = ...
0128 dAbr_dV(dSf_dVa_full, dSf_dVm_full, dSt_dVa_full, dSt_dVm_full, Sf, St);
0129
0130 [dAf_dVa, dAf_dVm, dAt_dVa, dAt_dVm] = ...
0131 dAbr_dV(dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St);
0132 dAf_dVa_sp = full(dAf_dVa);
0133 dAf_dVm_sp = full(dAf_dVm);
0134 dAt_dVa_sp = full(dAt_dVa);
0135 dAt_dVm_sp = full(dAt_dVm);
0136
0137
0138 num_dAf_dVm = full( (abs(Smpf).^2 - abs(Sf2).^2) / pert );
0139 num_dAf_dVa = full( (abs(Sapf).^2 - abs(Sf2).^2) / pert );
0140 num_dAt_dVm = full( (abs(Smpt).^2 - abs(St2).^2) / pert );
0141 num_dAt_dVa = full( (abs(Sapt).^2 - abs(St2).^2) / pert );
0142
0143 t_is(dAf_dVm_sp, num_dAf_dVm, 4, 'dAf_dVm (sparse)');
0144 t_is(dAf_dVa_sp, num_dAf_dVa, 4, 'dAf_dVa (sparse)');
0145 t_is(dAt_dVm_sp, num_dAt_dVm, 4, 'dAt_dVm (sparse)');
0146 t_is(dAt_dVa_sp, num_dAt_dVa, 4, 'dAt_dVa (sparse)');
0147 t_is(dAf_dVm_full, num_dAf_dVm, 4, 'dAf_dVm (full)');
0148 t_is(dAf_dVa_full, num_dAf_dVa, 4, 'dAf_dVa (full)');
0149 t_is(dAt_dVm_full, num_dAt_dVm, 4, 'dAt_dVm (full)');
0150 t_is(dAt_dVa_full, num_dAt_dVa, 4, 'dAt_dVa (full)');
0151
0152
0153
0154 [dIf_dVa_full, dIf_dVm_full, dIt_dVa_full, dIt_dVm_full, If, It] = dIbr_dV(branch, Yf_full, Yt_full, V);
0155
0156
0157 [dIf_dVa, dIf_dVm, dIt_dVa, dIt_dVm, If, It] = dIbr_dV(branch, Yf, Yt, V);
0158 dIf_dVa_sp = full(dIf_dVa);
0159 dIf_dVm_sp = full(dIf_dVm);
0160 dIt_dVa_sp = full(dIt_dVa);
0161 dIt_dVm_sp = full(dIt_dVm);
0162
0163
0164 num_dIf_dVm = full( (Yf * Vmp - Yf * V*ones(1,nb)) / pert );
0165 num_dIf_dVa = full( (Yf * Vap - Yf * V*ones(1,nb)) / pert );
0166 num_dIt_dVm = full( (Yt * Vmp - Yt * V*ones(1,nb)) / pert );
0167 num_dIt_dVa = full( (Yt * Vap - Yt * V*ones(1,nb)) / pert );
0168
0169 t_is(dIf_dVm_sp, num_dIf_dVm, 5, 'dIf_dVm (sparse)');
0170 t_is(dIf_dVa_sp, num_dIf_dVa, 5, 'dIf_dVa (sparse)');
0171 t_is(dIt_dVm_sp, num_dIt_dVm, 5, 'dIt_dVm (sparse)');
0172 t_is(dIt_dVa_sp, num_dIt_dVa, 5, 'dIt_dVa (sparse)');
0173 t_is(dIf_dVm_full, num_dIf_dVm, 5, 'dIf_dVm (full)');
0174 t_is(dIf_dVa_full, num_dIf_dVa, 5, 'dIf_dVa (full)');
0175 t_is(dIt_dVm_full, num_dIt_dVm, 5, 'dIt_dVm (full)');
0176 t_is(dIt_dVa_full, num_dIt_dVa, 5, 'dIt_dVa (full)');
0177
0178 t_end;