MAXLOADLIM computes the maximum loadability limit in one direction. It
uses dispatchable loads in MATPOWER
RESULTS = MAXLOADLIM(MPC,DIR_MLL) returns the results from the
optimization problem looking for the maximum loadability limit in
the direction of load increase DIR_MLL. DIR_MLL defines the directions
of load increases for all buses. For buses with zero loads, the
direction of load increases must be zero. RESULTS contains all fields
returned from the runopf MATPOWER function. It also contains the
following additional fields:
* dir_mll: the direction of load increase used as input.
* stab_marg: the stability margin to the maximum loadability point from
the base case defined in the input MPC.
* bif: information about the bifurcation at the MLL point.
RESULTS = MAXLOADLIM(MPC,DIR_MLL,NAME,VALUE) uses the options defined
by the pair NAME,VALUE. The currently supported options are
* 'verbose': 1 or 0 (Default). If set to 1, a summary of the results
at the maximum loadability limit is printed.
* 'use_qlim': 1 (Default) or 0. Enforces or not the reactive power
limits of the generators.
* 'Vlims_bus_nb': [] (Default) or array of integers. By default, the
bus voltage limits are not enforced. This option allows for defining
a set of buses at which the voltage limits are enforced.
See also PREPARE_MAXLOADLIM, POSTPROC_MAXLOADLIM, PRINT_MAXLOADLIM,
RUNOPF.0001 function results = maxloadlim(mpc,dir_mll,varargin) 0002 % MAXLOADLIM computes the maximum loadability limit in one direction. It 0003 % uses dispatchable loads in MATPOWER 0004 % RESULTS = MAXLOADLIM(MPC,DIR_MLL) returns the results from the 0005 % optimization problem looking for the maximum loadability limit in 0006 % the direction of load increase DIR_MLL. DIR_MLL defines the directions 0007 % of load increases for all buses. For buses with zero loads, the 0008 % direction of load increases must be zero. RESULTS contains all fields 0009 % returned from the runopf MATPOWER function. It also contains the 0010 % following additional fields: 0011 % * dir_mll: the direction of load increase used as input. 0012 % * stab_marg: the stability margin to the maximum loadability point from 0013 % the base case defined in the input MPC. 0014 % * bif: information about the bifurcation at the MLL point. 0015 % 0016 % RESULTS = MAXLOADLIM(MPC,DIR_MLL,NAME,VALUE) uses the options defined 0017 % by the pair NAME,VALUE. The currently supported options are 0018 % * 'verbose': 1 or 0 (Default). If set to 1, a summary of the results 0019 % at the maximum loadability limit is printed. 0020 % * 'use_qlim': 1 (Default) or 0. Enforces or not the reactive power 0021 % limits of the generators. 0022 % * 'Vlims_bus_nb': [] (Default) or array of integers. By default, the 0023 % bus voltage limits are not enforced. This option allows for defining 0024 % a set of buses at which the voltage limits are enforced. 0025 % 0026 % See also PREPARE_MAXLOADLIM, POSTPROC_MAXLOADLIM, PRINT_MAXLOADLIM, 0027 % RUNOPF. 0028 0029 % MATPOWER 0030 % Copyright (c) 2015-2016, Power Systems Engineering Research Center (PSERC) 0031 % by Camille Hamon 0032 % 0033 % This file is part of MATPOWER/mx-maxloadlim. 0034 % Covered by the 3-clause BSD License (see LICENSE file for details). 0035 % See https://github.com/MATPOWER/mx-maxloadlim/ for more info. 0036 0037 define_constants; 0038 0039 %% Checking the options, if any 0040 input_checker = inputParser; 0041 0042 default_verbose = 0; 0043 verbose_levels = [0;1]; 0044 check_verbose = @(x)(isnumeric(x) && isscalar(x) && any(x == verbose_levels)); 0045 addParameter(input_checker,'verbose',default_verbose,check_verbose); 0046 0047 % Direction of change for generators 0048 default_dir_var_gen = []; 0049 check_dir_var_gen = @(x)(isempty(x) || (isnumeric(x) && iscolumn(x))); 0050 addParameter(input_checker,'dir_var_gen',default_dir_var_gen,check_dir_var_gen); 0051 0052 % Generator numbers of the variable generators; 0053 default_idx_var_gen = []; 0054 check_idx_var_gen = @(x)(isempty(x) || (isnumeric(x) && iscolumn(x))); 0055 addParameter(input_checker,'idx_var_gen',default_idx_var_gen,check_idx_var_gen); 0056 0057 input_checker.KeepUnmatched = true; 0058 parse(input_checker,varargin{:}); 0059 0060 options = input_checker.Results; 0061 0062 %% Iterate the process when considering variable generators 0063 cur_stab_marg = 0; 0064 idx_var_gen = options.idx_var_gen; 0065 dir_var_gen = options.dir_var_gen; 0066 nb_var_gen = length(idx_var_gen); 0067 repeat = 1; 0068 iter = 0; 0069 iter_max = nb_var_gen; 0070 settings = varargin; 0071 0072 while iter <= iter_max && repeat 0073 if options.verbose 0074 fprintf(1,'Beginning of iteration %d\n',iter); 0075 end 0076 iter = iter + 1; 0077 0078 %% Prepare the matpower case for the maximum loadability limit problem 0079 % We remove reactive power limits of the slack bus since, by 0080 % assumption, the slack bus is a strong grid. 0081 [ref,~] = bustypes(mpc.bus,mpc.gen); 0082 gen_ref = ismember(mpc.gen(:,GEN_BUS),ref); 0083 mpc.gen(gen_ref,QMAX) = 9999; 0084 mpc.gen(gen_ref,QMIN) = -9999; 0085 mpc_vl = prepare_maxloadlim(mpc,dir_mll,settings{:}); 0086 0087 %% Run opf 0088 % Turning off the printing and initializing from the base case 0089 mpopt = mpoption('verbose',options.verbose,'opf.init_from_mpc',1); 0090 mpopt = mpoption(mpopt,'out.all',0); 0091 % Decreasing the threshold for the relative complementarity constraints 0092 mpopt = mpoption(mpopt,'mips.comptol',1e-8); 0093 % Change solver 0094 mpopt = mpoption(mpopt,'opf.ac.solver','MIPS'); 0095 % Execute opf 0096 results = runopf(mpc_vl,mpopt); 0097 0098 %% Post-processing 0099 results = postproc_maxloadlim(results,dir_mll); 0100 % update the stability margin with that of the current iteration 0101 results.stab_marg = results.stab_marg+cur_stab_marg; 0102 cur_stab_marg = results.stab_marg; 0103 0104 %% Check if it stopped because of a variable generator reached its PMAX 0105 % We check the Lagrangian multiplier of Pg<=PMAX 0106 if isempty(idx_var_gen) 0107 gen_hit_pmax = 0; 0108 else 0109 gen_hit_pmax = (abs(results.gen(idx_var_gen,PMAX)-results.gen(idx_var_gen,PG)) < 5e-4) | ... 0110 (results.var.mu.u.Pg(idx_var_gen) > 1e-4); %1e-4 for numerical error 0111 end 0112 if sum(gen_hit_pmax) == 0 0113 % The OPF did not stop because PG = PMAX so we are done. 0114 repeat = 0; 0115 else 0116 % We remove it from the direction of change and the set of variable 0117 % generators 0118 idx_var_gen(gen_hit_pmax) = []; 0119 dir_var_gen(gen_hit_pmax) = []; 0120 % We update the settings that are passed down to the function 0121 % preparing the constraints 0122 var_gen_set = find(strcmp(settings,'idx_var_gen')); 0123 settings{var_gen_set+1} = idx_var_gen; 0124 dir_gen_set = find(strcmp(settings,'dir_var_gen')); 0125 settings{dir_gen_set+1} = dir_var_gen; 0126 % We update mpc with the results from the current iteration 0127 mpc.bus(:,[PD QD]) = results.bus(:,[PD QD]); 0128 mpc.gen(:,PG) = results.gen(:,PG); 0129 mpc.bus(:,[VM VA]) = results.bus(:,[VM VA]); 0130 % % Remove OPF infor 0131 % mpc.gen(:,MU_PMAX:MU_QMIN) = []; 0132 end 0133 end 0134 0135 %% Printing 0136 if options.verbose 0137 print_maxloadlim(mpc,results); 0138 end