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