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

t_pf

PURPOSE ^

T_PF Tests for power flow solvers.

SYNOPSIS ^

function t_pf(quiet)

DESCRIPTION ^

T_PF  Tests for power flow solvers.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function t_pf(quiet)
0002 %T_PF  Tests for power flow solvers.
0003 
0004 %   MATPOWER
0005 %   Copyright (c) 2004-2019, Power Systems Engineering Research Center (PSERC)
0006 %   by Ray Zimmerman, PSERC Cornell
0007 %
0008 %   This file is part of MATPOWER.
0009 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0010 %   See https://matpower.org for more info.
0011 
0012 if nargin < 1
0013     quiet = 0;
0014 end
0015 
0016 AC_alg = {'NR', 'NR-SC', 'NR-IP', 'NR-IC', 'FDXB', 'FDBX', 'GS'};
0017 AC_name = {
0018     'Newton (default, power-polar)',
0019     'Newton (power-cartesian)',
0020     'Newton (current-polar)',
0021     'Newton (current-cartesian)',
0022     'Fast Decoupled (XB)',
0023     'Fast Decoupled (BX)',
0024     'Gauss-Seidel'
0025 };
0026 
0027 t_begin(length(AC_alg)*44 + 14, quiet);
0028 
0029 casefile = 't_case9_pf';
0030 if quiet
0031     verbose = 0;
0032 else
0033     verbose = 1;
0034 end
0035 if have_fcn('octave')
0036     if have_fcn('octave', 'vnum') >= 4
0037         file_in_path_warn_id = 'Octave:data-file-in-path';
0038     else
0039         file_in_path_warn_id = 'Octave:load-file-in-path';
0040     end
0041     s1 = warning('query', file_in_path_warn_id);
0042     warning('off', file_in_path_warn_id);
0043 end
0044 mpopt0 = mpoption('out.all', 0, 'pf.tol', 1e-9, 'verbose', 0);
0045 mpopt0 = mpoption(mpopt0, 'verbose', verbose);
0046 
0047 %% define named indices into bus, gen, branch matrices
0048 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0049     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0050 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0051     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0052     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0053 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0054     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0055     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0056 
0057 %% network with islands
0058 mpc0 = loadcase(casefile);
0059 mpc0.gen(1, PG) = 60;
0060 mpc0.gen(1, [PMIN PMAX QMIN QMAX PG QG]) = mpc0.gen(1, [PMIN PMAX QMIN QMAX PG QG]) / 2;
0061 mpc0.gen = [mpc0.gen(1, :); mpc0.gen];
0062 mpc1 = mpc0;
0063 mpc  = mpc0;
0064 nb = size(mpc.bus, 1);
0065 mpc1.bus(:, BUS_I)      = mpc1.bus(:, BUS_I) + nb;
0066 mpc1.branch(:, F_BUS)   = mpc1.branch(:, F_BUS) + nb;
0067 mpc1.branch(:, T_BUS)   = mpc1.branch(:, T_BUS) + nb;
0068 mpc1.gen(:, GEN_BUS)    = mpc1.gen(:, GEN_BUS) + nb;
0069 mpc.bus         = [mpc.bus; mpc1.bus];
0070 mpc.branch      = [mpc.branch; mpc1.branch];
0071 mpc.gen         = [mpc.gen; mpc1.gen];
0072 mpc1 = mpc;
0073 
0074 %%-----  AC power flow  -----
0075 %% get solved AC power flow case from MAT-file
0076 load soln9_pf;      %% defines bus_soln, gen_soln, branch_soln
0077 
0078 %% run AC PF
0079 for k = 1:length(AC_alg)
0080     t = sprintf('AC PF - %s : ', AC_name{k});
0081     mpopt = mpoption(mpopt0, 'pf.alg', AC_alg{k});
0082     [baseMVA, bus, gen, branch, success, et] = runpf(casefile, mpopt);
0083     t_ok(success, [t 'success']);
0084     t_is(bus, bus_soln, 6, [t 'bus']);
0085     t_is(gen, gen_soln, 6, [t 'gen']);
0086     t_is(branch, branch_soln, 6, [t 'branch']);
0087 
0088     r = runpf(casefile, mpopt);
0089     t_ok(r.success, [t 'success']);
0090     t_is(r.bus, bus_soln, 6, [t 'bus']);
0091     t_is(r.gen, gen_soln, 6, [t 'gen']);
0092     t_is(r.branch, branch_soln, 6, [t 'branch']);
0093 
0094     %% check Qg distribution, when Qmin = Qmax
0095     t = sprintf('%s - check Qg : ', AC_alg{k});
0096     mpopt = mpoption(mpopt, 'pf.alg', AC_alg{k}, 'verbose', 0);
0097     mpc = loadcase(casefile);
0098     mpc.gen(1, [QMIN QMAX]) = [20 20];
0099     r = runpf(mpc, mpopt);
0100     t_ok(r.success, [t 'success']);
0101     t_is(r.gen(1, QG), 24.07, 2, [t 'single gen, Qmin = Qmax']);
0102 
0103     mpc.gen = [mpc.gen(1, :); mpc.gen];
0104     mpc.gen(1, [QMIN QMAX]) = [10 10];
0105     mpc.gen(2, [QMIN QMAX]) = [0 50];
0106     r = runpf(mpc, mpopt);
0107     t_ok(r.success, [t 'success']);
0108     t_is(r.gen(1:2, QG), [10; 14.07], 2, [t '2 gens, Qmin = Qmax for one']);
0109 
0110     mpc.gen(1, [QMIN QMAX]) = [10 10];
0111     mpc.gen(2, [QMIN QMAX]) = [-50 -50];
0112     r = runpf(mpc, mpopt);
0113     t_ok(r.success, [t 'success']);
0114     t_is(r.gen(1:2, QG), [42.03; -17.97], 2, [t '2 gens, Qmin = Qmax for both']);
0115 
0116     mpc.gen(1, [QMIN QMAX]) = [0 50];
0117     mpc.gen(2, [QMIN QMAX]) = [0 100];
0118     r = runpf(mpc, mpopt);
0119     t_ok(r.success, [t 'success']);
0120     t_is(r.gen(1:2, QG), [8.02; 16.05], 2, [t '2 gens, proportional']);
0121 
0122     mpc.gen(1, [QMIN QMAX]) = [-50 0];
0123     mpc.gen(2, [QMIN QMAX]) = [50 150];
0124     r = runpf(mpc, mpopt);
0125     t_ok(r.success, [t 'success']);
0126     t_is(r.gen(1:2, QG), [-50+8.02; 50+16.05], 2, [t '2 gens, proportional']);
0127 
0128     mpc.gen(1, [QMIN QMAX]) = [-50 Inf];
0129     mpc.gen(2, [QMIN QMAX]) = [50 150];
0130     r = runpf(mpc, mpopt);
0131     t_ok(r.success, [t 'success']);
0132     t_is(r.gen(1:2, QG), [-31.61; 55.68], 2, [t '2 gens, one infinite range']);
0133 
0134     mpc.gen(1, [QMIN QMAX]) = [-50 Inf];
0135     mpc.gen(2, [QMIN QMAX]) = [50 Inf];
0136     r = runpf(mpc, mpopt);
0137     t_ok(r.success, [t 'success']);
0138     t_is(r.gen(1:2, QG), [-33.12; 57.18], 2, [t '2 gens, both infinite range']);
0139 
0140     mpc.gen(1, [QMIN QMAX]) = [-50 Inf];
0141     mpc.gen(2, [QMIN QMAX]) = [-Inf 150];
0142     r = runpf(mpc, mpopt);
0143     t_ok(r.success, [t 'success']);
0144     t_is(r.gen(1:2, QG), [76.07; -52], 2, [t '2 gens, both infinite range']);
0145 
0146     mpc.gen(1, [QMIN QMAX]) = [-Inf Inf];
0147     mpc.gen(2, [QMIN QMAX]) = [-Inf Inf];
0148     r = runpf(mpc, mpopt);
0149     t_ok(r.success, [t 'success']);
0150     t_is(r.gen(1:2, QG), [12.03; 12.03], 2, [t '2 gens, both infinite range']);
0151 
0152     t = sprintf('%s - reactive generation allocation : ', AC_alg{k});
0153     mpc = loadcase(casefile);
0154     %% generator data
0155     %    bus    Pg    Qg    Qmax    Qmin    Vg    mBase    status    Pmax    Pmin    Pc1    Pc2    Qc1min    Qc1max    Qc2min    Qc2max    ramp_agc    ramp_10    ramp_30    ramp_q    apf
0156     mpc.gen = [
0157         1    0    0    300    -300    1    100    1    250    10    0    0    0    0    0    0    0    0    0    0    0;
0158         2    54    0    0    -5    1    100    1    300    10    0    0    0    0    0    0    0    0    0    0    0;
0159         2    54    0    5    -5    1    100    1    300    10    0    0    0    0    0    0    0    0    0    0    0;
0160         2    55    0    25     10    1    100    1    300    10    0    0    0    0    0    0    0    0    0    0    0;
0161         30    25    1    300    -300    1    100    1    270    10    0    0    0    0    0    0    0    0    0    0    0;
0162         30    30    2    300    -300    1    100    1    270    10    0    0    0    0    0    0    0    0    0    0    0;
0163         30    30    -3    300    -300    1    100    1    270    10    0    0    0    0    0    0    0    0    0    0    0;
0164     ];
0165     mpc.bus(3, BUS_TYPE) = PQ;
0166     r = runpf(mpc, mpopt);
0167     t_ok(r.success, [t 'success']);
0168     t_is(r.gen(2:4, QG), [-5; -5; 10] + [1; 2; 3]*1.989129794, 7, [t 'PV bus']);
0169     t_is(r.gen(5:7, QG), [1; 2; -3], 8, [t 'PQ bus']);
0170 
0171     %% network with islands
0172     t = sprintf('%s - network w/islands : AC PF : ', AC_alg{k});
0173     mpc = mpc1;
0174     r = runpf(mpc, mpopt);
0175     t_ok(r.success, [t 'success']);
0176     t_is(r.bus( 1:9,  VA), bus_soln(:, VA), 7, [t 'voltage angles 1']);
0177     t_is(r.bus(10:18, VA), bus_soln(:, VA), 7, [t 'voltage angles 2']);
0178     Pg = [gen_soln(1, PG)-30; 30; gen_soln(2:3, PG)];
0179     t_is(r.gen(1:4, PG), Pg, 6, [t 'active power generation 1']);
0180     t_is(r.gen(5:8, PG), Pg, 6, [t 'active power generation 2']);
0181 
0182     t = sprintf('%s - all buses isolated : ', AC_alg{k});
0183     mpc.bus(:, BUS_TYPE) = NONE;
0184     try
0185         r = runpf(mpc, mpopt);
0186         t_is(r.success, 0, 12, [t 'success = 0']);
0187     catch
0188         t_ok(0, [t 'unexpected fatal error']);
0189     end
0190 
0191     %% case 14 with Q limits
0192     t = sprintf('%s - pf.enforce_q_lims == 0 : ', AC_alg{k});
0193     mpc = loadcase('case14');
0194     mpc.gen(1, QMIN) = -10;
0195     mpc.gen(:, QMAX) = [10; 30; 29; 15; 15];
0196     bt0 = mpc.bus(:, BUS_TYPE);
0197     bt = bt0;
0198     mpopt = mpoption(mpopt, 'pf.enforce_q_lims', 0);
0199     r = runpf(mpc, mpopt);
0200     t_ok(r.success, [t 'success']);
0201     t_is(r.bus(:, BUS_TYPE), bt, 12, [t 'bus type']);
0202     t_is(r.gen(:, QG), [-16.549300542; 43.557100134; 25.075348495; 12.730944405; 17.623451366], 6, [t 'Qg']);
0203 
0204     t = sprintf('%s - pf.enforce_q_lims == 1 : ', AC_alg{k});
0205     mpopt = mpoption(mpopt, 'pf.enforce_q_lims', 1);
0206     r = runpf(mpc, mpopt);
0207     bt = bt0;
0208     bt([1 2 3 8]) = [PQ PQ REF PQ];
0209     t_is(r.success, 0, 12, [t 'success = 0']);
0210     t_is(r.bus(:, BUS_TYPE), bt, 12, [t 'bus type']);
0211     t_is(r.gen(:, QG), [-10; 30; 31.608422873; 16.420423190; 15], 4, [t 'Qg']);
0212 
0213     t = sprintf('%s - pf.enforce_q_lims == 2 : ', AC_alg{k});
0214     mpopt = mpoption(mpopt, 'pf.enforce_q_lims', 2);
0215     r = runpf(mpc, mpopt);
0216     bt = bt0;
0217     bt([1 2 3 6 8]) = [REF PQ PQ PQ PQ];
0218     t_ok(r.success, [t 'success']);
0219     t_is(r.bus(:, BUS_TYPE), bt, 12, [t 'bus type']);
0220     t_is(r.gen(:, QG), [-6.30936644; 30; 29; 15; 15], 6, [t 'Qg']);
0221 end
0222 
0223 %%-----  DC power flow  -----
0224 mpopt = mpoption(mpopt, 'verbose', verbose);
0225 %% get solved AC power flow case from MAT-file
0226 load soln9_dcpf;        %% defines bus_soln, gen_soln, branch_soln
0227 
0228 %% run DC PF
0229 t = 'DC PF : ';
0230 [baseMVA, bus, gen, branch, success, et] = rundcpf(casefile, mpopt);
0231 t_ok(success, [t 'success']);
0232 t_is(bus, bus_soln, 6, [t 'bus']);
0233 t_is(gen, gen_soln, 6, [t 'gen']);
0234 t_is(branch, branch_soln, 6, [t 'branch']);
0235 r = rundcpf(casefile, mpopt);
0236 t_ok(r.success, [t 'success']);
0237 t_is(r.bus, bus_soln, 6, [t 'bus']);
0238 t_is(r.gen, gen_soln, 6, [t 'gen']);
0239 t_is(r.branch, branch_soln, 6, [t 'branch']);
0240 
0241 %% network with islands
0242 t = sprintf('DC PF - network w/islands : ');
0243 mpc  = mpc1;
0244 %mpopt = mpoption(mpopt, 'out.bus', 1, 'out.gen', 1, 'out.all', -1, 'verbose', 2);
0245 r = rundcpf(mpc, mpopt);
0246 t_ok(r.success, [t 'success']);
0247 t_is(r.bus( 1:9,  VA), bus_soln(:, VA), 8, [t 'voltage angles 1']);
0248 t_is(r.bus(10:18, VA), bus_soln(:, VA), 8, [t 'voltage angles 2']);
0249 Pg = [gen_soln(1, PG)-30; 30; gen_soln(2:3, PG)];
0250 t_is(r.gen(1:4, PG), Pg, 8, [t 'active power generation 1']);
0251 t_is(r.gen(5:8, PG), Pg, 8, [t 'active power generation 1']);
0252 
0253 %% island without slack bus (catch singluar matrix?)
0254 t = sprintf('DC PF - network w/islands w/o slack : ');
0255 k = find(mpc.bus(:, BUS_TYPE) == REF);
0256 mpc.bus(k(2), BUS_TYPE) = PV;
0257 warn_state = warning;
0258 warning('off', 'all');  %% turn of (near-)singular matrix warnings
0259 r = rundcpf(mpc, mpopt);
0260 warning(warn_state);
0261 t_is(r.success, 0, 12, [t 'success = 0']);
0262 
0263 t_end;
0264 
0265 if have_fcn('octave')
0266     warning(s1.state, file_in_path_warn_id);
0267 end
Generated on Mon 24-Jun-2019 15:58:45 by m2html © 2005