prepare_maxloadlim

 PREPARE_MAXSLL prepares the mpc case for computing the maximum
 loadability limit by adding the relevant constraints and variables.
   MPC_VL = PREPARE_MAXLL(MPC,DIR_MLL) returns the matpower case MPC_VL
   prepared from MPC by transforming all loads to dispatchable loads,
   adding a field for the direction DIR_MLL of load increase and adapting 
   limits for the OPF formulation.
   MPC_VL = PREPARE_MAXLL(MPC,DIR_MLL,Name,Value) does the same with additional
   options specified in pairs Name,Value. The two supported options are as 
   follows:
     * '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.

0001 function mpc_vl = prepare_maxloadlim(mpc,dir_mll,varargin)
0002 % PREPARE_MAXSLL prepares the mpc case for computing the maximum
0003 % loadability limit by adding the relevant constraints and variables.
0004 %   MPC_VL = PREPARE_MAXLL(MPC,DIR_MLL) returns the matpower case MPC_VL
0005 %   prepared from MPC by transforming all loads to dispatchable loads,
0006 %   adding a field for the direction DIR_MLL of load increase and adapting
0007 %   limits for the OPF formulation.
0008 %   MPC_VL = PREPARE_MAXLL(MPC,DIR_MLL,Name,Value) does the same with additional
0009 %   options specified in pairs Name,Value. The two supported options are as
0010 %   follows:
0011 %     * 'use_qlim': 1 (Default) or 0. Enforces or not the reactive power
0012 %     limits of the generators.
0013 %     * 'Vlims_bus_nb': [] (Default) or array of integers. By default, the
0014 %     bus voltage limits are not enforced. This option allows for defining
0015 %     a set of buses at which the voltage limits are enforced.
0016 
0017 %   MATPOWER
0018 %   Copyright (c) 2015-2016, Power Systems Engineering Research Center (PSERC)
0019 %   by Camille Hamon
0020 %
0021 %   This file is part of MATPOWER/mx-maxloadlim.
0022 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0023 %   See https://github.com/MATPOWER/mx-maxloadlim/ for more info.
0024 
0025 define_constants;
0026 n_gen = size(mpc.gen,1);
0027 mpc0 = mpc;
0028 %% Checking the options, if any
0029 input_checker = inputParser;
0030 
0031 % Q-lims
0032 default_qlim = 1;
0033 check_qlim = @(x)(isnumeric(x) && isscalar(x));
0034 addParameter(input_checker,'use_qlim',default_qlim,check_qlim);
0035 
0036 % Enfore V-lims
0037 default_vlim = [];
0038 check_vlim = @(x)(isnumeric(x) && all(floor(x) == ceil(x))); % expects array of integer values (bus numbers)
0039 addParameter(input_checker,'Vlims_bus_nb',default_vlim,check_vlim);
0040 
0041 % Direction of change for generators
0042 default_dir_var_gen = [];
0043 check_dir_var_gen = @(x)(isempty(x) || (isnumeric(x) && iscolumn(x)));
0044 addParameter(input_checker,'dir_var_gen',default_dir_var_gen,check_dir_var_gen);
0045 
0046 % Generator numbers of the variable generators;
0047 default_idx_var_gen = [];
0048 check_idx_var_gen = @(x)(isempty(x) || (isnumeric(x) && iscolumn(x)));
0049 addParameter(input_checker,'idx_var_gen',default_idx_var_gen,check_idx_var_gen);
0050 
0051 % Parse
0052 input_checker.KeepUnmatched = true;
0053 parse(input_checker,varargin{:});
0054 options = input_checker.Results;
0055 
0056 %% CHECKS
0057 % Check whether the number of directions of load increases is equal to the
0058 % number of buses
0059 if size(mpc.bus,1) ~= length(dir_mll)
0060     error_msg = ['The number of directions of load increases ' ...
0061         'is not equal to the number of buses'];
0062     error(error_msg);
0063 end
0064 
0065 % Check whether load increases have been defined for zero loads
0066 idx_zero_loads = mpc.bus(:,PD) == 0;
0067 if sum(dir_mll(idx_zero_loads))>0
0068     error('Directions of load increases cannot be defined for zero loads.');
0069 end
0070 
0071 % Check whether the number of variable generators is equal to the number of
0072 % elements in the direction of change of generators
0073 if size(options.idx_var_gen,1) ~= size(options.dir_var_gen,1)
0074     error('The number of variable generators does not match the direction vector.');
0075 end
0076 
0077 % Make sure that the slack bus is not included among the variable
0078 % generators
0079 [ref_init, ~] = bustypes(mpc.bus, mpc.gen);
0080 idx_gen_slack = find(mpc.gen(1:n_gen,GEN_BUS) == ref_init);
0081 if sum(ismember(options.idx_var_gen,idx_gen_slack)) ~= 0
0082     error('The direction vector cannot include changes at the slack bus');
0083 end
0084 
0085 %% Preparation of the case mpc_vl
0086 % Initialise with a power flow with q-lims considered
0087 mpopt = mpoption('pf.enforce_q_lims', options.use_qlim,'verbose',0,'out.all',0);
0088 mpc = runpf(mpc,mpopt);
0089 % Reset the bus types after the power flow
0090 mpc.bus(:,BUS_TYPE) = mpc0.bus(:,BUS_TYPE);
0091 % Convert all loads to dispatchable
0092 mpc_vl = load2disp(mpc);
0093 
0094 % Extract the part of dir_mll corresponding to nonzero loads
0095 dir_mll = dir_mll(mpc.bus(:, PD) > 0);
0096 
0097 % Add a field to mpc_vl for the load increase
0098 mpc_vl.dir_mll = dir_mll;
0099 
0100 % Create a vector with the direction for the variable generators with as
0101 % many elements as the number of gens + number of sheddable loads
0102 % (necessary for conversion to internal indexing).
0103 dir_var_gen_all = sparse(options.idx_var_gen,1,options.dir_var_gen,...
0104     size(mpc_vl.gen,1),1);
0105 % Add a field for the generators
0106 mpc_vl.dir_var_gen = options.dir_var_gen;
0107 mpc_vl.idx_var_gen = options.idx_var_gen;
0108 mpc_vl.dir_var_gen_all = dir_var_gen_all;
0109 
0110 % Adjust the Pmin of dispatchable loads to make them negative enough so
0111 % that the max load lim can be found
0112 idx_vl = isload(mpc_vl.gen);
0113 mpc_vl.gen(idx_vl,PMIN) = 300*mpc_vl.gen(idx_vl,PMIN);
0114 % Adjust Qmin so that Qmin/Pmin is the power factor of the load, if
0115 % inductive, and change Qmax if the load is capacitive
0116 idx_vl_inductive = idx_vl & mpc_vl.gen(:,QMAX) == 0;
0117 idx_vl_capacitive = idx_vl & mpc_vl.gen(:,QMIN) == 0;
0118 tanphi_vl_ind = mpc_vl.gen(idx_vl_inductive,QG)./mpc_vl.gen(idx_vl_inductive,PG);
0119 tanphi_vl_cap = mpc_vl.gen(idx_vl_capacitive,QG)./mpc_vl.gen(idx_vl_capacitive,PG);
0120 mpc_vl.gen(idx_vl_inductive,QMIN) = mpc_vl.gen(idx_vl_inductive,PMIN).*tanphi_vl_ind;
0121 mpc_vl.gen(idx_vl_capacitive,QMAX) = mpc_vl.gen(idx_vl_capacitive,PMIN).*tanphi_vl_cap;
0122 % Make the cost zero
0123 mpc_vl.gencost(:,COST:end) = 0;
0124 % Make the non variable generators not dispatchable
0125 % Note, we look only for the real PV buses, i.e. we do not consider the
0126 % dispatchable loads in this search. Hence the search over 1:n_gen
0127 [ref, pv, pq] = bustypes(mpc_vl.bus, mpc_vl.gen);
0128 idx_gen_pv = find(ismember(mpc_vl.gen(1:n_gen,GEN_BUS),pv));
0129 % Note we filter out non variable generators with nonzero direction of
0130 % increase from the pv buses
0131 idx_non_var_pv = setdiff(idx_gen_pv,options.idx_var_gen(options.dir_var_gen~=0));
0132 mpc_vl.gen(idx_non_var_pv,PMIN) = mpc_vl.gen(idx_non_var_pv,PG);
0133 mpc_vl.gen(idx_non_var_pv,PMAX) = mpc_vl.gen(idx_non_var_pv,PG);
0134 % Raise the flow limits so that they are not binding
0135 mpc_vl.branch(:,RATE_A) = 9999;%1e5;
0136 % Raise the slack bus limits so that they are not binding
0137 idx_gen_slack = find(mpc_vl.gen(1:n_gen,GEN_BUS) == ref);
0138 mpc_vl.gen(idx_gen_slack,[QMAX,PMAX]) = 9999;
0139 % Decrease the slack bus limits so that they are not binding
0140 mpc_vl.gen(idx_gen_slack,[QMIN,PMIN]) = -9999;
0141 % Change the voltage constraints of the PQ buses so that they are not
0142 % binding
0143 mpc_vl.bus(pq,VMIN) = 0.01;
0144 mpc_vl.bus(pq,VMAX) = 10;
0145 % Lock the voltages of the slack bus
0146 mpc_vl.bus(ref,VMAX) = mpc_vl.gen(idx_gen_slack(1),VG);
0147 mpc_vl.bus(ref,VMIN) = mpc_vl.gen(idx_gen_slack(1),VG);
0148 % Put Vmax = Vset and low Vmin for all pv buses
0149 for bb = 1:length(pv)
0150     idx_gen_at_bb = find(ismember(mpc_vl.gen(1:n_gen,GEN_BUS),pv(bb)));
0151     mpc_vl.bus(pv(bb),VMAX) = mpc_vl.gen(idx_gen_at_bb(1),VG);
0152     mpc_vl.bus(pv(bb),VMIN) = 0.01;
0153 %     mpc_vl.bus(pv(bb),VMAX) = mpc_vl.gen(idx_gen_at_bb(1),VG);
0154 %     mpc_vl.bus(pv(bb),VMIN) = mpc_vl.gen(idx_gen_at_bb(1),VG);
0155 end
0156 % If we do not consider Qlim, increase Qmax and decrease Qmin
0157 % of all generators to arbitrarily large values
0158 if ~options.use_qlim
0159     mpc_vl.gen(idx_gen_pv,QMAX) = 9999;
0160     mpc_vl.gen(idx_gen_pv,QMIN) = -9999;
0161     for bb = 1:length(pv)
0162         idx_gen_at_bb = find(ismember(mpc_vl.gen(1:n_gen,GEN_BUS),pv(bb)));
0163         mpc_vl.bus(pv(bb),VMAX) = mpc_vl.gen(idx_gen_at_bb(1),VG);
0164         mpc_vl.bus(pv(bb),VMIN) = mpc_vl.gen(idx_gen_at_bb(1),VG);
0165     end
0166 end
0167 if ~isempty(options.Vlims_bus_nb)
0168     idx_gen_vlim = find(ismember(mpc_vl.gen(:,GEN_BUS),options.Vlims_bus_nb));
0169     mpc_vl.bus(options.Vlims_bus_nb,VMAX) = mpc_vl.gen(idx_gen_vlim,VG);
0170     mpc_vl.bus(options.Vlims_bus_nb,VMIN) = mpc_vl.gen(idx_gen_vlim,VG);
0171 end
0172 
0173 % Convert external to internal indexing
0174 mpc_vl = add_userfcn(mpc_vl, 'ext2int', @userfcn_direction_mll_ext2int);
0175 % Build the constraint for enforcing the direction of load increase
0176 mpc_vl = add_userfcn(mpc_vl, 'formulation', @userfcn_direction_mll_formulation);