Running a 4-index MILP model in OPTModel Procedure

[thossain]

Hi! I am trying to run a 4 index MILP in OPTModel procedure. I am using SAS Enterprise Guise 7.1. I have a number a{i,f,p} which is the supplied amount. I have a dataset from where I can read this data. This is how I have done it. I am just entering two constraints here which have a {i,f,p} in them. But I have many others. Same goes for many other of my parameters which are not listed as they will make the code large here. The problem is, if I do it this way, it was able to create the a{i,f,p} and populate it with the required data from the input dataset. But, it becomes an error while writing the constraint. As I have used the {i,f} as a set, in the constraint it can not recognize them separately anymore. I am trying to create them as a Node but having trouble with this issue. Is there any easier way to do this? Can I somehow populate my number a{i,f,p} directly without creating the {i,f} set? 

Any help would be highly appreciated. Thank you!

proc optmodel;
	set <str> Fields;
	set <str> Depots init {};
	set <str> Feedstock;
	set <str> Price;
	set <str, str> FF; /* this is the set of fields and feedstock to create the number a{i,f,p} with three indices*/

	num aifp {FF,Price} init 0;

read data ORout.Fields into Fields = [fips] x y;
read data ORout.Fields into Depots = [fips] x y;
read data ORout.Price into Price = [Pr] gp = number ;
read data ORout.vfcost into Feedstock = [Feedstock] vf = variable cf = fixed 
	fs=fst ds=dst Qf=que pf=prep mf = dmloss;
read data ORout.supply into FF = [fips Feedstock]
	{p in Price} < aifp[fips, Feedstock, p] = col ("Pr"||p) >;


var Assign{Fields, Depots, Feedstock, Price} binary;
var Build {Depots} binary;
var Amount{Fields, Depots, Feedstock, Price} integer >= 0;


/********Some Constraints************/

/****Amount supplied from a field can not be greater than the available supply of that field****/
con Supply_constraints{i in Fields, j in Depots, f in Feedstock, p in Price}:
		Amount[i,j,f,p] <= aifp[i,f,p]* Assign[i,j,f,p];

/***Considering the dry matter loss and deducting that mount from field supply******/
/***Supplied amount and recived amount will not be exactly same due to dry matter loss***/
con Dry_Matter {j in Depots, f in Feedstock, p in Price}:
		sum{i in Fields}(1-mf[f])*aifp[i,f,p] = sum{i in Fields} Amount [i,j,f,p];

quit;

[RobPratt]

If I understand correctly, FF is a subset of Fields cross Feedstock.  Here is a "dense" approach that treats a[i,f,p] as 0 for <i,f> pairs that do not appear in FF:

   num aifp {Fields,Feedstock,Price} init 0;

   read data ORout.supply into [fips Feedstock]
      {p in Price} < aifp[fips, Feedstock, p] = col ("Pr"||p) >;

   var Assign{Fields, Depots, Feedstock, Price} binary;

   var Amount{Fields, Depots, Feedstock, Price} integer >= 0;

   con Supply_constraints{i in Fields, j in Depots, f in Feedstock, p in Price}:
      Amount[i,j,f,p] <= aifp[i,f,p]* Assign[i,j,f,p];

   con Dry_Matter {j in Depots, f in Feedstock, p in Price}:
      sum{i in Fields}(1-mf[f])*aifp[i,f,p] = sum{i in Fields} Amount [i,j,f,p];

Alternatively, here is a more efficient "sparse" approach, as in the Sparse Modeling documentation example:

   num aifp {FF,Price} init 0;

   read data ORout.supply into FF = [fips Feedstock]
      {p in Price} < aifp[fips, Feedstock, p] = col ("Pr"||p) >;

   var Assign{FF, Depots, Price} binary;

   var Amount{FF, Depots, Price} integer >= 0;

   con Supply_constraints{<i,f> in FF, j in Depots, p in Price}:
      Amount[i,f,j,p] <= aifp[i,f,p]* Assign[i,f,j,p];

   con Dry_Matter {j in Depots, f in Feedstock, p in Price}:
      sum{<i,(f)> in FF}(1-mf[f])*aifp[i,f,p] = sum{<i,(f)> in FF} Amount [i,f,j,p];

Note the difference in order of indices for both Assign and Amount, with i,f,j,p instead of i,j,f,p.

[thossain]

Thank you so much for your reply! I did try the first method and it worked! However, I am having issues with the feasibility of the problem. It is coming as infeasible. I tried to reduce the complexity by making aifp two dimensional a[I,f] or a[I,p]. For both tries, it is showing me in the log that,
NOTE: The row activity (0) of constraint Assign_fields['37013',CS] violates its lower bound (1). Basically for the very first row of the data, Assign_fields violates. I don't understand why is this happening. Even when I used aifp[I,f,p], it showed me the same error and was infeasible.

I am attaching the datalines and constraints here. How can I get over the infeasibility? Any suggestions? 

Thanks in advance!

data test;
	input fips $6. CS SW;
	datalines;
37013  33000 37000
37031  31000 29000
37095  44000 48000
37097  29000 21000
37137  13000 13000
37141  19000 20000
37171  16000 17000
37177  24000 27000
37187  18000 83000
37197  17000 10000
45027  48000 50000
45031  18000 18000
45061  29000 26000
45075  53000 57000
45085  35000 36000
;
proc optmodel;

/***Declaring the sets****/
	set <str> Fields;
	set <str> Depots init {};
	set <str> Feedstock;
 
/*Distance parameters*/;
    num x{Fields union Depots};
	num y{Fields union Depots};
	num dij{i in Fields, j in Depots}=sqrt((x[i] - x[j])^2 + (y[i] 
		- y[j])^2);

/*Field cost parameters*/
	num supply{Fields,Feedstock} init 0;
	num gp = 70; /*Farmgate price including grower payment and harvest and collection cost*/
	num fs{Feedstock} init 0; /*field site storage*/
  

/*Transportation cost parameters*/
	num vf {Feedstock} init 0;
	num cf {Feedstock} init 0;
	num tijf {i in Fields, j in Depots,f in Feedstock}= cf[f] + vf[f]*dij[i,j];
 	num C = 150000;
 	num N = 504000;
 	num pf{Feedstock} init 0;
 	num ds{Feedstock} init 0;
 	num Qf{Feedstock} init 0;

read data ORout.Fields into Fields = [fips] x y;
read data ORout.Depots into Depots = [fips] x y;
read data ORout.vfcost into Feedstock = [Feed] vf = variable cf = fixed 
	fs=fst ds=dst Qf=que pf=prep;
read data test into [fips]{f in Feedstock}< supply[fips,f]=col(f)>;

var Assign{Fields, Depots, Feedstock} binary;
var Build {Depots} binary;


/*****Objective function********/
min Cost = sum{i in Fields, j in Depots, f in Feedstock} (gp+ fs[f] + ds[f] + tijf[i,j,f] + Qf[f] + 
           pf[f])* Assign[i,j,f]*supply [i,f]+ sum{j in Depots} N*Build[j];

/****Each field has to be assigned to a depot*/
con Assign_fields{i in Fields, f in Feedstock}: 
		sum{j in Depots} Assign[i,j,f]=1;
con Build_depot{i in Fields, j in Depots,f in Feedstock}: 
		Assign[i,j,f] <=Build[j];

/****Total amount supplied to a depot has to be less than its capacity*******/
con Depot_Capacity{j in Depots}: 
		sum {i in Fields, f in Feedstock}Supply[i,f]*Assign[i,j,f]<=C*Build[j];
	
solve obj Cost with MILP/Primalin;
quit;

[RobPratt]

Can you please provide the other three data sets as well?

[thossain]

I am attaching the two files for depot and field location. I tried at the beginning with many other constraints and another integer variable X[i,j,f,p] which is supplied feedstock for price f from field i to depot j. But then it was coming as infeasible so made the model much simpler now. But still it is not working. The constant values are:

num C = 150000
num N = 504000

I am also attaching how the model was at the beginning, just to make it clear. Although, the field and depot location were different for this one. The vfcost dataset is also with the code. Thank you so much! I am having a really hard time solving these! 

data ORout.vfcost;
input feed $4. fixed variable fst dst que prep dmloss;
datalines;
 CS  3.42 0.114 3.97 0.88 1.21 19.4 0.12
 SW  3.42 0.114 3.02 0.88 1.34 18.77 0.08
 ;

data ORout.Price ;
input name $6. pr;
datalines; 
Pr60  60
Pr70  70
Pr80  80
Pr90  90
Pr100  100
;
data test;
	input fips $6. feed $4. Pr60 Pr70 Pr80 Pr90 Pr100;
	datalines;
37035  1  890 950 60 0 160
37035  2  0 330 330 330 330
37039  2  0 50 3 21 55
37041  1  1400 1600 1900 2200 2600
37041  2  11000 18000 19000 19000 19000
37043  2  0	31 58 660 970
37047  2  63000 63000 63000 16000 16000
37049  1  0	0 500 2300 2200
37049  2  31000 31000 31000 31000 15000
37051  2  23000 23000 23000 23000 23000
37055  2  2200 1900 1000 940 890
37057  1  5600 4500 4500 4500 4500
37057  2  0	0 18 240 240
37059  1  9100 6300 6400 6400 6400
37059  2  11 0 0 200 200
;

Proc optmodel;

set <str> Fields;
set <str> Depots init {};
set <str> Feedstock;
set <str> Price;

read data Fieldtest into Fields = [fips] x y;
read data Fieldtest into Depots = [fips] x y;
read data ORout.Price into Price = [name] gp = Pr ;
read data ORout.vfcost into Feedstock = [Feed] vf = variable cf = fixed 
	fs=fst ds=dst Qf=que pf=prep mf = dmloss;
read data test into [fips Feed]{p in Price}<aifp[fips,Feed,p]=col(p)>;

/***Defining the variables***/
var Assign{Fields, Depots, Feedstock, Price} binary;
var Build {Depots} binary;
var Amount{Fields, Depots, Feedstock, Price} integer >= 0;


/****Objective function*******/
min Cost = sum{i in Fields, j in Depots, f in Feedstock, p in Price} (gp[p]+ fs[f] + tijf[i,j,f] + Qf[f] + 
           ds[f] +  pf[f])* Amount [i,j,f,p] + sum{j in Depots} N*Build[j];


/********Constraints************/

/****Each field has to be assigned to a depot, field supply can't be stored and they have to be shipped****/
con Assign_fields{i in Fields, f in Feedstock, p in Price}: 
		sum{j in Depots} Assign[i,j,f,p]=1;

/***One field can not supply same type of biomass at different prices to the same depot location*****/
con Cost_constraints {i in Fields, j in Depots, f in Feedstock}:
		sum{p in Price}	Assign [i,j,f,p] <= 1;

/****Amount supplied from a field can not be greater than the available supply of that field****/
con Supply_constraints{i in Fields, j in Depots, f in Feedstock,p in Price}:
		Amount[i,j,f,p] <= aifp[i,f,p]* Assign[i,j,f,p];


/****Total amount supplied to a depot has to be less than its capacity*******/
con Depot_Capacity{j in Depots}:
		sum{i in Fields, f in feedstock, p in Price} Amount [i,j,f,p] <= C*Build [j];

/*** If a depot is built, it has to be supplied with a minimum amount of biomass****/
/***To increase efficiency***/
con Depot_utilization{j in Depots}:
		sum{i in Fields, f in feedstock, p in Price} Amount [i,j,f,p] >= U*C*Build [j];

	
solve obj Cost with MILP/Primalin;

print Amount;

quit;

[RobPratt]

I am able to replicate this behavior, and even the LP relaxation is infeasible.  To help diagnose the cause of infeasibility, you can run these statements:

      solve with lp relaxint / iis=true;
      expand / iis;

On my machine, this returns a set of 5 variables and 35 constraints that together are infeasible.

[thossain]

Thank you so much! I will try with this. Also, SAS is not letting me to put constraint on my integer decision variable which needs to be greater than zero. It is showing me the error that hard inequality are not allowed. And wants me to put an equal sign after the inequality. Is there any way that I can declare:
var Amount{Fields, Depots, Feedstock, Price} integer > 0;

Thank you again!

[RobPratt]

You are correct that OPTMODEL does not allow strict inequalities for variable bounds.  To enforce > 0 for an integer variable, you can specify >= 1.