userfcn_direction_mll_formulation

 USERFCN_DIRECTION_MLL_FORMULATION adds one variable and as many
 constraints as dispatchable loads to enforce the load increase direction

0001 function om = userfcn_direction_mll_formulation(om,args)
0002 % USERFCN_DIRECTION_MLL_FORMULATION adds one variable and as many
0003 % constraints as dispatchable loads to enforce the load increase direction
0004 
0005 %   MATPOWER
0006 %   Copyright (c) 2015-2016, Power Systems Engineering Research Center (PSERC)
0007 %   by Camille Hamon
0008 %
0009 %   This file is part of MATPOWER.
0010 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0011 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0012 
0013 define_constants;
0014 mpc = get_mpc(om);
0015 dir_mll = mpc.dir_mll;
0016 
0017 % identify dispatchable loads
0018 idx_vl = isload(mpc.gen);
0019 n_vl = sum(idx_vl);
0020 n_g = size(mpc.gen,1)-n_vl;
0021 if length(dir_mll) ~= n_vl
0022     error_msg = ['The number of dispatchable loads is not equal to the '...
0023         'length of the direction vector'];
0024     error(error_msg);
0025 end
0026 % Add the amount of load increase alpha with 0 <= alpha <= inf
0027 om = add_vars(om,'alpha',1,0,0,inf);
0028 
0029 %% Load increase
0030 % Add the constraint for enforcing the direction of load increase
0031 Pl0 = -mpc.gen(idx_vl,PG)/mpc.baseMVA;
0032 % finding the internal indices of the dispatchable loads
0033 int_idx_disp_loads = find(idx_vl); 
0034 idx_A_dirmll_i = [1:n_vl 1:n_vl]';
0035 idx_A_dirmll_j = [int_idx_disp_loads' (n_g+n_vl+1)*ones(1,n_vl)]';
0036 vals_A_dirmll = [ones(n_vl,1);dir_mll];
0037 A_dirmll = sparse(idx_A_dirmll_i,idx_A_dirmll_j,vals_A_dirmll,n_vl,n_g+n_vl+1);
0038 om = add_constraints(om,'dir_mll',A_dirmll,-Pl0,-Pl0,{'Pg','alpha'});
0039 % Add cost of alpha to -1 to maximize loads in the given direction
0040 om = add_costs(om,'alpha_cost',struct('Cw',-1),{'alpha'});
0041 
0042 %% Generator changes
0043 % Add the constraint for enforcing the direction of generation change
0044 idx_var_gen = find(mpc.dir_var_gen_all);
0045 if ~isempty(idx_var_gen)
0046     % Add the amount of generator change beta with 0 <= beta <= inf
0047     om = add_vars(om,'beta',1,0,0,inf);
0048     Pg0 = mpc.gen(idx_var_gen,PG)/mpc.baseMVA;
0049     nb_var_gen = length(idx_var_gen);
0050     idx_A_var_gen_i = [1:nb_var_gen 1:nb_var_gen]'; % Constraint number
0051     idx_A_var_gen_j = [idx_var_gen;(n_g+n_vl+1)*ones(nb_var_gen,1)]; % Generator number and beta column (column in constraint matrix)
0052     vals_A_var_gen = [ones(nb_var_gen,1);-nonzeros(mpc.dir_var_gen_all)];
0053     A_var_gen = sparse(idx_A_var_gen_i,idx_A_var_gen_j,vals_A_var_gen,nb_var_gen,n_g+n_vl+1);
0054     om = add_constraints(om,'dir_var_gen',A_var_gen,...
0055         Pg0,Pg0,{'Pg','beta'});
0056     % Add cost of beta to -1 to maximize the generator change in a given
0057     % direction
0058     om = add_costs(om,'beta_cost',struct('Cw',-1e3),{'beta'});
0059     % Add constraint for beta <= alpha
0060     A_beta_alpha = [-1 1];
0061     om = add_constraints(om,'beta_alpha',A_beta_alpha,-Inf,0,{'alpha','beta'});
0062 end