Test_maxloadlim

TEST_MAXLOADLIM

0001 %TEST_MAXLOADLIM
0002 
0003 %   MATPOWER
0004 %   Copyright (c) 2015-2016, Power Systems Engineering Research Center (PSERC)
0005 %   by Camille Hamon
0006 %
0007 %   This file is part of MATPOWER.
0008 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0009 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0010 
0011 classdef Test_maxloadlim < matlab.unittest.TestCase
0012     
0013     properties
0014         systems = {'case2','case9'};
0015         load_dir = struct('case2',[0 1],...
0016             'case9',[0 0 0 0 1 0 0 0 0;
0017             0 0 0 0 0 0 1 0 0;
0018             0 0 0 0 0 0 0 0 1;
0019             0 0 0 0 1 0 1 0 1;
0020             0 0 0 0 1 0 1 0 0;
0021             0 0 0 0 1 0 0 0 1;
0022             0 0 0 0 0 0 1 0 1]',...
0023             'case39',[eye(39);ones(1,39)]');
0024         gen_dir = struct('case9',[0 1 0;
0025             0 0 1;
0026             0 1 1]',...
0027             'case39',[zeros(9,1) eye(9);0 ones(1,9)]');
0028         % For the case 2 the theoretical result is P = E^2/(2*X)
0029         max_load_lims = struct('case2',5);
0030         threshold_MW = 1;
0031     end
0032     properties(TestParameter)
0033         idx_dir_ieee9 = num2cell(1:7);
0034         idx_dir_ieee39 = num2cell(1:40);
0035         gen_var_ieee9 = num2cell(1:3);
0036         gen_var_ieee39 = num2cell(1:10);
0037     end
0038     methods(TestClassSetup)
0039         function setuptests(testcase)
0040             % This is to turn off the warnings in the CPF when we get close
0041             % to the nose
0042             warning off MATLAB:singularMatrix
0043             warning off MATLAB:nearlySingularMatrix
0044         end
0045     end
0046     
0047     methods(TestClassTeardown)
0048         function teardownOnce(testcase)
0049             warning on MATLAB:singularMatrix
0050             warning on MATLAB:nearlySingularMatrix
0051         end
0052     end
0053     
0054     methods(Test)     
0055         function testAgainstCPF_case39(testCase,idx_dir_ieee39)
0056             define_constants;
0057             mpc = loadcase('case39');
0058             idx_nonzero_loads = mpc.bus(:,PD) > 0;
0059             dir_all = testCase.load_dir.('case39');
0060             dirCPF = dir_all(:,idx_dir_ieee39);
0061             dirCPF(~idx_nonzero_loads)=0;
0062             dirCPF2 = dirCPF(idx_nonzero_loads);
0063             if sum(dirCPF2) == 0
0064                 % case of load increase in non-zero load
0065                 testCase.verifyEqual(1,1);
0066             else
0067                 switch idx_dir_ieee39
0068                     case 9
0069                         msg = ['This case fails because the load at bus 9 is capacitive' ...
0070                             ' which means that some generators hit their lower reactive limits' ...
0071                             ' , a case that is not handled by my CPF'];
0072                         fprintf(1,[msg '\n']);
0073                     case 31
0074                         msg = ['This case fails because the load at bus 31 is small'...
0075                             ' and is connected to the slack bus that absorbs any load increase'...
0076                             ' so that the CPF procedure reaches the maximum number of iterations'...
0077                             ' before reaching the MLL'];
0078                         fprintf(1,[msg '\n']);
0079                 end
0080                 results_cpf = ch_runCPF('case39','',0,dirCPF2);
0081                 max_loads_cpf = results_cpf.bus(:,PD)*mpc.baseMVA;
0082                 results_mll = maxloadlim(mpc,dirCPF,'verbose',0);
0083                 max_loads_mll = results_mll.bus(:,PD);
0084                 testCase.verifyEqual(max_loads_mll,max_loads_cpf,'RelTol',0.05);
0085                 % NOTE: we used a 5% relative tolerance here, rather than
0086                 % an absolute tolerance, because the base case of IEEE39 is
0087                 % not feasible for the Q lims, which seem to bring about
0088                 % small discrepancies.
0089             end
0090         end
0091         
0092         function testAgainstCPF_case9(testCase,idx_dir_ieee9)
0093             define_constants;
0094             % Loading the case
0095             mpc = loadcase('case9');
0096             dir_all = testCase.load_dir.('case9');
0097             dir = dir_all(:,idx_dir_ieee9);
0098             results_mll = maxloadlim(mpc,dir,'verbose',0);
0099             max_loads_mll = results_mll.bus(:,PD);
0100             % Remember to set chooseStartPoint to 0 in ch_runCPF
0101             idx_nonzero_loads = mpc.bus(:,PD) > 0;
0102             dirCPF = dir(idx_nonzero_loads);
0103             results_cpf = ch_runCPF('case9static','loads568',0,dirCPF);
0104             max_loads_cpf = results_cpf.bus(:,PD)*mpc.baseMVA;           
0105             testCase.verifyEqual(max_loads_mll,max_loads_cpf,'AbsTol',testCase.threshold_MW);
0106         end
0107         
0108         function testAgainstMatpowerCPF_case9(testCase,idx_dir_ieee9)
0109             define_constants;
0110             % Loading the case
0111             mpc = loadcase('case9');
0112             dir_all = testCase.load_dir.('case9');
0113             dir = dir_all(:,idx_dir_ieee9);
0114             % Preparing the target case for Matpower CPF
0115             mpc_target = mpc;
0116             nonzero_loads = mpc_target.bus(:,PD) ~= 0;
0117             Q_P = mpc_target.bus(nonzero_loads,QD)./mpc_target.bus(nonzero_loads,PD);
0118             mpc_target.bus(:,PD) = mpc_target.bus(:,PD)+2*dir*mpc_target.baseMVA;
0119             mpc_target.bus(nonzero_loads,QD) = Q_P.*mpc_target.bus(nonzero_loads,PD);
0120             % Run the CPF with matpower
0121             [results,~] = runcpf(mpc,mpc_target,mpoption('out.all',0));
0122             % Extract the maximum loads
0123             max_loads_cpf = results.bus(:,PD);
0124             % Solve the maximum loadability limit without considering
0125             % reactive power limits
0126             results_mll = maxloadlim(mpc,dir,'use_qlim',0);
0127             % Extract the maximum loads
0128             max_loads_mll = results_mll.bus(:,PD);
0129             % We compare with a precision of 0.5MW
0130             testCase.verifyEqual(max_loads_mll,max_loads_cpf,'AbsTol',testCase.threshold_MW);
0131         end
0132         
0133         function testAgainstMatpowerCPF_case39(testCase,idx_dir_ieee39)
0134             define_constants;
0135             % Loading the case
0136             mpc = loadcase('case39');
0137             idx_nonzero_loads = mpc.bus(:,PD) > 0;
0138             dir_all = testCase.load_dir.('case39');
0139             dir = dir_all(:,idx_dir_ieee39);
0140             dir(~idx_nonzero_loads)=0;
0141             if sum(dir) == 0 || idx_dir_ieee39 == 31
0142                 % The code does not currently support load increase at
0143                 % nonzero loads.
0144                 % The MATPOWER CPF takes long time for increase at bus 31
0145                 % which is the slack bus.
0146                 testCase.verifyEqual(1,1);
0147             else
0148                 % Preparing the target case for Matpower CPF
0149                 mpc_target = mpc;
0150                 nonzero_loads = mpc_target.bus(:,PD) ~= 0;
0151                 Q_P = mpc_target.bus(nonzero_loads,QD)./mpc_target.bus(nonzero_loads,PD);
0152                 mpc_target.bus(:,PD) = mpc_target.bus(:,PD)+2*dir*mpc_target.baseMVA;
0153                 mpc_target.bus(nonzero_loads,QD) = Q_P.*mpc_target.bus(nonzero_loads,PD);
0154                 % Run the CPF with matpower
0155                 [results,~] = runcpf(mpc,mpc_target,mpoption('out.all',0));
0156                 % Extract the maximum loads
0157                 max_loads_cpf = results.bus(:,PD);
0158                 % Solve the maximum loadability limit without considering
0159                 % reactive power limits
0160                 results_mll = maxloadlim(mpc,dir,'use_qlim',0);
0161                 % Extract the maximum loads
0162                 max_loads_mll = results_mll.bus(:,PD);
0163                 % We compare with a precision of 0.5MW
0164                 testCase.verifyEqual(max_loads_mll,max_loads_cpf,'AbsTol',testCase.threshold_MW);
0165             end
0166         end
0167         
0168         function testAgainstTheoretical_case2(testCase)
0169             define_constants;
0170             % Loading the case
0171             mpc = loadcase('case2');
0172             dir = testCase.load_dir.('case2');
0173             res_maxloadlim = maxloadlim(mpc,dir);
0174             % Get nonzero loads
0175             idx_nonzero_loads = res_maxloadlim.bus(:,PD) > 0;
0176             max_loads = res_maxloadlim.bus(idx_nonzero_loads,PD)/res_maxloadlim.baseMVA;
0177             max_loads_theo = testCase.max_load_lims.('case2');
0178             testCase.verifyEqual(max_loads,max_loads_theo,'AbsTol',testCase.threshold_MW);
0179         end
0180         
0181         function testVarGen_case9(testCase,idx_dir_ieee9,gen_var_ieee9)
0182             define_constants;
0183             % Loading the case
0184             mpc = loadcase('case9');
0185             idx_nonzero_loads = mpc.bus(:,PD) > 0;
0186             dir_all = testCase.load_dir.('case9');
0187             dir_load = dir_all(:,idx_dir_ieee9);
0188             dir_load(~idx_nonzero_loads)=0;
0189             dir_var_gen_all = testCase.gen_dir.('case9');
0190             dir_var_gen = dir_var_gen_all(:,gen_var_ieee9);
0191             idx_var_gen = find(dir_var_gen);
0192             dir_var_gen = dir_var_gen(idx_var_gen);
0193             % Normalizing with respect to both loads and gens
0194             gen_load_dir = [dir_load;dir_var_gen];
0195             dir_load = dir_load/norm(gen_load_dir);
0196             dir_var_gen = dir_var_gen/norm(dir_var_gen);
0197             % Find MLL in the direction of load and gen increase
0198             results_with_gens = maxloadlim(mpc,dir_load,'verbose',0,'idx_var_gen',idx_var_gen,'dir_var_gen',dir_var_gen);
0199             % Set gens to their values in previous results and re-run in
0200             % load space only
0201             mpc2 = mpc;
0202             mpc2.gen(:,PG) = results_with_gens.gen(:,PG);
0203             dir_load2 = dir_load/norm(dir_load);
0204             % Find MLL in the direction of load and gen increase
0205             results_without_gens = maxloadlim(mpc2,dir_load2,'verbose',0);
0206             % Compare the maximum loads
0207             mll_with_gen = results_with_gens.bus(:,PD);
0208             mll_without_gen = results_without_gens.bus(:,PD);
0209             testCase.verifyEqual(mll_with_gen,mll_without_gen,'AbsTol',testCase.threshold_MW);
0210         end
0211         
0212         function testVarGen_case39(testCase,idx_dir_ieee39,gen_var_ieee39)
0213             define_constants;
0214             % Loading the case
0215             mpc = loadcase('case39');
0216             idx_nonzero_loads = mpc.bus(:,PD) > 0;
0217             dir_all = testCase.load_dir.('case39');
0218             dir_load = dir_all(:,idx_dir_ieee39);
0219             dir_load(~idx_nonzero_loads)=0;
0220             dir_var_gen_all = testCase.gen_dir.('case39');
0221             % Switching rows 1 and 2 (gens 30 and 31) because row 1
0222             % is zero and corresponds to the slack bus which is gen 31
0223             dir_var_gen_all([1 2],:) = dir_var_gen_all([2 1],:);
0224             dir_var_gen = dir_var_gen_all(:,gen_var_ieee39);
0225             idx_var_gen = find(dir_var_gen);
0226             dir_var_gen = dir_var_gen(idx_var_gen);
0227             % Normalizing with respect to both loads and gens
0228             gen_load_dir = [dir_load;dir_var_gen];
0229             dir_load = dir_load/norm(gen_load_dir);
0230             dir_var_gen = dir_var_gen/norm(dir_var_gen);
0231             if sum(dir_load) == 0
0232                 % case of load increase in non-zero load
0233                 testCase.verifyEqual(1,1);
0234             else
0235                 % Find MLL in the direction of load and gen increase
0236                 results_with_gens = maxloadlim(mpc,dir_load,'verbose',0,'idx_var_gen',idx_var_gen,'dir_var_gen',dir_var_gen);
0237                 % Set gens to their values in previous results and re-run in
0238                 % load space only
0239                 mpc2 = mpc;
0240                 mpc2.gen(:,PG) = results_with_gens.gen(:,PG);
0241                 dir_load2 = dir_load/norm(dir_load);
0242                 % Find MLL in the direction of load and gen increase
0243                 results_without_gens = maxloadlim(mpc2,dir_load2,'verbose',0);
0244                 % Compare the maximum loads
0245                 mll_with_gen = results_with_gens.bus(:,PD);
0246                 mll_without_gen = results_without_gens.bus(:,PD);
0247                 testCase.verifyEqual(mll_with_gen,mll_without_gen,'RelTol',0.05);
0248             end
0249         end
0250     end
0251 end