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.
0022 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0023 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0024 
0025 define_constants;
0026 n_gen = size(mpc.gen,1);
0027 
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)(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)(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 
0090 % Convert all loads to dispatchable
0091 mpc_vl = load2disp(mpc);
0092 
0093 % Extract the part of dir_mll corresponding to nonzero loads
0094 dir_mll = dir_mll(mpc.bus(:, PD) > 0);
0095 
0096 % Add a field to mpc_vl for the load increase
0097 mpc_vl.dir_mll = dir_mll;
0098 
0099 % Create a vector with the direction for the variable generators with as
0100 % many elements as the number of gens + number of sheddable loads
0101 % (necessary for conversion to internal indexing).
0102 dir_var_gen_all = sparse(options.idx_var_gen,1,options.dir_var_gen,...
0103     size(mpc_vl.gen,1),1);
0104 % Add a field for the generators
0105 mpc_vl.dir_var_gen = options.dir_var_gen;
0106 mpc_vl.idx_var_gen = options.idx_var_gen;
0107 mpc_vl.dir_var_gen_all = dir_var_gen_all;
0108 
0109 % Adjust the Pmin of dispatchable loads to make them negative enough so
0110 % that the max load lim can be found
0111 idx_vl = isload(mpc_vl.gen);
0112 mpc_vl.gen(idx_vl,PMIN) = 300*mpc_vl.gen(idx_vl,PMIN);
0113 % Adjust Qmin so that Qmin/Pmin is the power factor of the load, if
0114 % inductive, and change Qmax if the load is capacitive
0115 idx_vl_inductive = idx_vl & mpc_vl.gen(:,QMAX) == 0;
0116 idx_vl_capacitive = idx_vl & mpc_vl.gen(:,QMIN) == 0;
0117 tanphi_vl_ind = mpc_vl.gen(idx_vl_inductive,QG)./mpc_vl.gen(idx_vl_inductive,PG);
0118 tanphi_vl_cap = mpc_vl.gen(idx_vl_capacitive,QG)./mpc_vl.gen(idx_vl_capacitive,PG);
0119 mpc_vl.gen(idx_vl_inductive,QMIN) = mpc_vl.gen(idx_vl_inductive,PMIN).*tanphi_vl_ind;
0120 mpc_vl.gen(idx_vl_capacitive,QMAX) = mpc_vl.gen(idx_vl_capacitive,PMIN).*tanphi_vl_cap;
0121 % Make the cost zero
0122 mpc_vl.gencost(:,COST:end) = 0;
0123 % Make the non variable generators not dispatchable
0124 % Note, we look only for the real PV buses, i.e. we do not consider the
0125 % dispatchable loads in this search. Hence the search over 1:n_gen
0126 [ref, pv, pq] = bustypes(mpc_vl.bus, mpc_vl.gen);
0127 idx_gen_pv = find(ismember(mpc_vl.gen(1:n_gen,GEN_BUS),pv));
0128 % Note we filter out non variable generators with nonzero direction of
0129 % increase from the pv buses
0130 idx_non_var_pv = setdiff(idx_gen_pv,options.idx_var_gen(options.dir_var_gen~=0));
0131 mpc_vl.gen(idx_non_var_pv,PMIN) = mpc_vl.gen(idx_non_var_pv,PG);
0132 mpc_vl.gen(idx_non_var_pv,PMAX) = mpc_vl.gen(idx_non_var_pv,PG);
0133 % Raise the flow limits so that they are not binding
0134 mpc_vl.branch(:,RATE_A) = 9999;%1e5;
0135 % Raise the slack bus limits so that they are not binding
0136 idx_gen_slack = find(mpc_vl.gen(1:n_gen,GEN_BUS) == ref);
0137 mpc_vl.gen(idx_gen_slack,[QMAX,PMAX]) = 9999;
0138 % Decrease the slack bus limits so that they are not binding
0139 mpc_vl.gen(idx_gen_slack,[QMIN,PMIN]) = -9999;
0140 % Change the voltage constraints of the PQ buses so that they are not
0141 % binding
0142 mpc_vl.bus(pq,VMIN) = 0.01;
0143 mpc_vl.bus(pq,VMAX) = 10;
0144 % Lock the voltages of the slack bus
0145 mpc_vl.bus(ref,VMAX) = mpc_vl.gen(idx_gen_slack(1),VG);
0146 mpc_vl.bus(ref,VMIN) = mpc_vl.gen(idx_gen_slack(1),VG);
0147 % Put Vmax = Vset and low Vmin for all pv buses
0148 for bb = 1:length(pv)
0149     idx_gen_at_bb = find(ismember(mpc_vl.gen(1:n_gen,GEN_BUS),pv(bb)));
0150     mpc_vl.bus(pv(bb),VMAX) = mpc_vl.gen(idx_gen_at_bb(1),VG);
0151     mpc_vl.bus(pv(bb),VMIN) = 0.01;
0152 end
0153 % If we do not consider Qlim, increase Qmax and decrease Qmin
0154 % of all generators to arbitrarily large values
0155 if ~options.use_qlim
0156     mpc_vl.gen(idx_gen_pv,QMAX) = 9999;
0157     mpc_vl.gen(idx_gen_pv,QMIN) = -9999;
0158     for bb = 1:length(pv)
0159         idx_gen_at_bb = find(ismember(mpc_vl.gen(1:n_gen,GEN_BUS),pv(bb)));
0160         mpc_vl.bus(pv(bb),VMAX) = mpc_vl.gen(idx_gen_at_bb(1),VG);
0161         mpc_vl.bus(pv(bb),VMIN) = mpc_vl.gen(idx_gen_at_bb(1),VG);
0162     end
0163 end
0164 if ~isempty(options.Vlims_bus_nb)
0165     idx_gen_vlim = find(ismember(mpc_vl.gen(:,GEN_BUS),options.Vlims_bus_nb));
0166     mpc_vl.bus(options.Vlims_bus_nb,VMAX) = mpc_vl.gen(idx_gen_vlim,VG);
0167     mpc_vl.bus(options.Vlims_bus_nb,VMIN) = mpc_vl.gen(idx_gen_vlim,VG);
0168 end
0169 
0170 % Convert external to internal indexing
0171 mpc_vl = add_userfcn(mpc_vl, 'ext2int', @userfcn_direction_mll_ext2int);
0172 % Build the constraint for enforcing the direction of load increase
0173 mpc_vl = add_userfcn(mpc_vl, 'formulation', @userfcn_direction_mll_formulation);