Home > matpower7.1 > lib > t > t_pf_radial.m

t_pf_radial

PURPOSE ^

T_PF_RADIAL Tests for distribution power flow solvers.

SYNOPSIS ^

function t_pf_radial(quiet)

DESCRIPTION ^

T_PF_RADIAL  Tests for distribution power flow solvers.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function t_pf_radial(quiet)
0002 %T_PF_RADIAL  Tests for distribution power flow solvers.
0003 
0004 if nargin < 1
0005     quiet = 0;
0006 end
0007 
0008 t_begin(329, quiet);
0009 
0010 if quiet
0011     verbose = 0;
0012 else
0013     verbose = 1;
0014 end
0015 if have_feature('octave')
0016     if have_feature('octave', 'vnum') >= 4
0017         file_in_path_warn_id = 'Octave:data-file-in-path';
0018     else
0019         file_in_path_warn_id = 'Octave:load-file-in-path';
0020     end
0021     s1 = warning('query', file_in_path_warn_id);
0022     warning('off', file_in_path_warn_id);
0023 end
0024 mpopt = mpoption('out.all', 0, 'verbose', verbose);
0025 
0026 %% define named indices into bus, gen, branch matrices
0027 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0028     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0029 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0030     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0031     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0032 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0033     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0034     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0035 
0036 %% Test Distribution Power Flow
0037 mpopt = mpoption(mpopt, 'pf.alg', 'NR', 'pf.radial.max_it', 100);
0038 casefile = {
0039     'case4_dist'
0040     'case18'
0041     'case22'
0042     'case33bw'
0043     'case69'
0044     'case85'
0045     'case141'
0046     };
0047 method = {
0048     'PQSUM' 'Power Summation'
0049     'ISUM'  'Current Summation'
0050     'YSUM'  'Admittance Summation'
0051     };
0052 
0053 % Original test cases (no PV buses)
0054 iter_expected = {
0055            'none'    'NR' 'PQSUM' 'ISUM' 'YSUM'
0056      'case4_dist'      3      12     12     11
0057          'case18'      4       9     12     10
0058          'case22'      3       3      6      6
0059        'case33bw'      3       5      8      8
0060          'case69'      4       5      8      8
0061          'case85'      4       5      9      9
0062         'case141'      3       4      7      7
0063         };
0064 iter_nopv = zeros(size(casefile,1),size(method,1)+1);
0065 for i = 1:length(casefile)
0066     % which row in iter_expected is for current casefile
0067     row = find(strcmp(casefile{i}, iter_expected(:,1)) == 1);
0068     mpc = loadcase(casefile{i});
0069     % solve it with NR
0070     r = runpf(mpc, mpopt);
0071     iter_nopv(i,1) = r.iterations;
0072     % which columng in iter_expected is for current method
0073     col = find(strcmp('NR', iter_expected(1,:)) == 1);
0074     t = ['Newton''s method, ' casefile{i} ' : '];
0075     t_is(r.iterations, iter_expected{row,col}, 6, [t 'iterations']);
0076     for j = 1:size(method,1)
0077         mpopt1 = mpoption(mpopt,'pf.alg',method{j,1});
0078         r1 = runpf(mpc,mpopt1);
0079         iter_nopv(i,j+1) = r1.iterations;
0080         % which columng in iter_expected is for current method
0081         col = find(strcmp(method{j,1}, iter_expected(1,:)) == 1);
0082         t = [method{j,2} ', ' casefile{i} ' : '];
0083         t_ok(r1.success, [t 'success']);
0084         t_is(r1.iterations, iter_expected{row,col}, 6, [t 'iterations']);
0085         t_is(r1.bus, r.bus, 6, [t 'bus']);
0086         t_is(r1.gen, r.gen, 6, [t 'gen']);
0087         t_is(r1.branch, r.branch, 6, [t 'branch']);
0088     end
0089 end
0090 
0091 % Test cases with added PV buses
0092 clear iter_expected
0093 iter_expected{1} = {
0094 % pf.radial.vcorr = 0
0095            'none'    'NR' 'PQSUM' 'ISUM' 'YSUM'
0096      'case4_dist'      3      12     12     11
0097          'case18'      4      20     16     13
0098          'case22'      4      11     25     25
0099        'case33bw'      4       8     16     16
0100          'case69'      4      12     27     27
0101          'case85'      5       9     17     17
0102         'case141'      5       9     18     18
0103 };
0104 iter_expected{2} = {
0105 % pf.radial.vcorr = 1
0106            'none'    'NR' 'PQSUM' 'ISUM' 'YSUM'
0107      'case4_dist'      3       7      6      7
0108          'case18'      4      10     12     12
0109          'case22'      4      12     14     14
0110        'case33bw'      4       8     11     11
0111          'case69'      4      13     17     17
0112          'case85'      5       9     13     13
0113         'case141'      5       9     12     12
0114 };
0115 iter_pv = zeros(size(casefile,1),size(method,1)+1,2);
0116 for i = 1:length(casefile)
0117     % which row in iter_expected is for current casefile
0118     row = find(strcmp(casefile{i}, iter_expected{1}(:,1)) == 1);
0119     mpc = loadcase(casefile{i});
0120     pv = mpc.bus(:,BUS_TYPE) == PV;
0121     if all(~pv)
0122         %%% ADD PV BUSES %%%
0123         N = 5; % How many PV buses should I add?
0124         % Solve the case without PV buses with Newton method
0125         r = runpf(mpc, mpopt);
0126         % Find the last N buses with lowest voltage
0127         [Vm, B] = sort(r.bus(:,VM));
0128         B = B(1:N);
0129         Vm = Vm(1:N);
0130         % add PV generators at buses in B
0131         % set VG 0.05 pu bigger then voltage at buses in B
0132         mpc.gen = repmat(mpc.gen,N+1,1);
0133         mpc.gen(2:end,GEN_BUS) = mpc.bus(B, BUS_I);
0134         mpc.gen(2:end,VG) = Vm + 0.05;
0135         mpc.bus(B,BUS_TYPE) = 2;
0136         if isfield(mpc,'gencost')
0137             mpc.gencost = repmat(mpc.gencost,N+1,1);
0138         end
0139         %%% END OF ADD PV BUSES %%%
0140     end
0141     r = runpf(mpc, mpopt);
0142     iter_pv(i,1,1) = r.iterations;
0143     iter_pv(i,1,2) = r.iterations;
0144     % which columng in iter_expected is for current method
0145     col = find(strcmp('NR', iter_expected{1}(1,:)) == 1);
0146     t = ['Newton''s method, ' casefile{i} ' : '];
0147     t_is(r.iterations, iter_expected{1}{row,col}, 6, [t 'iterations']);
0148     for j = 1:size(method,1)
0149         for vcorr = 0:1
0150             mpopt1 = mpoption(mpopt,'pf.alg',method{j,1},'pf.radial.vcorr',vcorr);
0151             r1 = runpf(mpc,mpopt1);
0152             iter_pv(i,j+1,vcorr+1) = r1.iterations;
0153             % which columng in iter_expected is for current method
0154             col = find(strcmp(method{j,1}, iter_expected{vcorr+1}(1,:)) == 1);
0155             t = [method{j,2} ', vcorr = ' num2str(vcorr) ', ' casefile{i} ' : '];
0156             t_ok(r1.success, [t 'success']);
0157             t_is(r1.iterations, iter_expected{vcorr+1}{row,col}, 6, [t 'iterations']);
0158             t_is(r1.bus, r.bus, 6, [t 'bus']);
0159             t_is(r1.gen, r.gen, 6, [t 'gen']);
0160             t_is(r1.branch, r.branch, 6, [t 'branch']);
0161         end
0162     end
0163 end
0164 
0165 t_end;
0166 
0167 if verbose > 0
0168     fprintf('\nITERATIONS: Original test cases (no PV buses), pf.radial.vcorr = 0\n');
0169     print_iterations(casefile,method,iter_nopv)
0170     fprintf('\nITERATIONS: Test cases with added PV buses, pf.radial.vcorr = 0\n');
0171     print_iterations(casefile,method,iter_pv(:,:,1))
0172     fprintf('\nITERATIONS: Test cases with added PV buses, pf.radial.vcorr = 1\n');
0173     print_iterations(casefile,method,iter_pv(:,:,2))
0174 end
0175 
0176 if have_feature('octave')
0177     warning(s1.state, file_in_path_warn_id);
0178 end
0179 
0180 function print_iterations(casefile,method,iterations)
0181 fprintf('%22s','NR');
0182 for j = 1:size(method,1)
0183     fprintf('%7s',method{j,1});
0184 end
0185 fprintf('\n');
0186 for i = 1:length(casefile)
0187     fprintf('%15s',casefile{i});
0188     for j = 1:size(method,1)+1
0189         fprintf('%7i',iterations(i,j));
0190     end
0191     fprintf('\n');
0192 end

Generated on Fri 09-Oct-2020 11:21:31 by m2html © 2005