Please help in below facility location code

[#1]

Dear SAS community members,

I need some help in below code. i wrote below code to solve a facility location problem but the result differs from the one i get from Excel solver. Could you please help explaining what is the wrong in the below code:

proc optmodel;

set DC = {'BO','NA','PR','SP','WO'};
set City = {'BO','BR','CO','HA','MN','NA','NH','NL','PO','PR','SP','WO'};

number transoprtationcost {DC,City} = [
0 37 42 82 34
93 65 106 59 68
69 33 105 101 72
98 103 73 27 66
55 20 92 93 60
37 0 72 79 41
128 137 94 63 98
95 113 57 57 71
62 48 104 127 85
42 72 0 68 38
82 79 68 0 47
34 41 38 47 0
];
number fixedcost {DC}=[10000 10000 10000 10000 10000];
number demand {City}=[425 12 43 125 110 86 129 28 66 320 220 182];
number supply {DC}=[2000 2000 2000 2000 2000];
number M= 1000000000;

var X{DC,City} integer >=0;
var Y{DC} binary;

min totalcost = sum {d in DC, c in City}
X[d,c]*transoprtationcost[d,c]+sum{d in DC}Y[d]*fixedcost[d];

con supplycons {d in DC}:sum{c in City} X[d,c]<=supply[d];
con demandcons {c in City}:sum{d in DC} X[d,c]>=demand[c];
con conserv {d in DC, c in City}:X[d,c]-M*Y[d]<=0;
con numberofDC:sum {d in DC} Y[d]=1;

solve;

print X Y totalcost;

Quit;

Regards,

Ahmed Taher;

[#2 · Accepted solution]

Dear All,

after a lot of tries finally i got the right soultion for this problem which will be as follows:

proc optmodel;

set DC = {'BO','NA','PR','SP','WO'};
set City = {'BO','BR','CO','HA','MN','NA','NH','NL','PO','PR','SP','WO'};

number transoprtationcost {DC,City} = [
0 93 69 98 55 37 128 95 62 42 82 34
37 65 33 103 20 0 137 113 48 72 79 41
42 106 105 73 92 72 94 57 104 0 68 38
82 59 101 27 93 79 63 57 127 68 0 47
34 68 72 66 60 41 98 71 85 38 47 0
];
number fixedcost {DC}=[10000 10000 10000 10000 10000];
number demand {City}=[425 12 43 125 110 86 129 28 66 320 220 182];
number supply {DC}=[1746 1746 1746 1746 1746];
number M= 100000;

var X{DC,City} integer >=0;
var Y{DC} binary;

minimize totalcost = sum {d in DC}sum{c in City}
X[d,c]*transoprtationcost[d,c]+sum{n in DC}fixedcost[n]*Y[n];

con supplycons {d in DC}:sum{c in City} X[d,c]<=supply[d];
con demandcons {c in City}:sum{d in DC} X[d,c]>=demand[c];
con conserv {d in DC, c in City}:X[d,c]-M*Y[d]<=0;
/*con numberofDC:sum {d in DC} Y[d]=1;*/

solve;

print X totalcost;

Quit;

[#3]

I have a couple of suggestions for improvement.

1. An arbitrarily large value like 100000 for M can cause numerical difficulties for optimization solvers.  In this case, the supplycons constraint implies that X[d,c] <= supply[d], so you can instead use the following tighter constraint:

con conserv {d in DC, c in City}:X[d,c]-supply[d]*Y[d]<=0;

2. You can in fact avoid the conserv constraint altogether if you tighten the supplycons constraint as follows:

con supplycons {d in DC}:sum{c in City} X[d,c]<=supply[d]*Y[d];

The resulting model (without conserv and with the new supplycons) is tighter, smaller, and more numerically stable.  For your small instance, even the original model solves quickly.  But these best practices can make a big difference for larger instances.

[#4]

Dear Rob,

First of all happy new year. 

I have a question regarding the above code I was asking how to add a constraints or a logical link to this code that force the source to ship only quantities that are divided by 5. 

Regards,

Thanks

[#5]

Do you want each X[d,c] to be a multiple of 5, or sum{c in City} X[d,c] to be a multiple of 5 for each d, or something else?

In any case, you can introduce a family of integer variables Z[] and constraints that force the expression to equal 5*Z[].

[#6]

Thanks a lot for your quick reply, however I didn't get the point. Would you please give me more explanation? For your question to be more clear in my question I want each X[d,c] to be a multiple of 5

Regards,
Ahmed Taher

[#7]

var Z{DC,City} integer >= 0;
con Five {d in DC, c in City}:
   X[d,c] = 5*Z[d,c];

[#8]

Thanks a lot Man. You are the hero of these communities.

[#9]

Glad to help.  By the way, here is an alternative approach that declares X as an implicit variable instead:

/*var X{DC,City} integer >=0;*/
var Z{DC,City} integer >= 0;
/*con Five {d in DC, c in City}:*/
/*   X[d,c] = 5*Z[d,c];*/
impvar X {d in DC, c in City} = 5*Z[d,c];

[#10]

Dear Boss @RobPratt RobPratt 🙂 

I need your help again in the below code. my issue is that i know the logical thinking of how to solve the issue so I know the how to solve it mathematically but since I'm not that godd in SAS programming I am asking for you kind help in programming. in the below code

I want to add constraint that says that I need DC 1 and DC 2 are open by default , and I need the optimization to choose the best DC from 'Potential DC 3', 'Potential DC 4', 'Potential DC 5'  in addition to DC 1 and 2 which are forced to be open as mentioned before. 

also I want to add a level of service constraint that says 75% of demand should be within 50 miles of a DC and the max average distance from a DC is <= 60 mile.

Your help is always appreciated. 

Proc optmodel;

set Plants = {'Plant 1','Contractor'};
set CentralWH = {'DC 1',	'DC 2',	'Potential DC 3',	'Potential DC 4',	'Potential DC 5'};
set RegionalWH = {'SO 1','SO 2','SO 3','SO 4','SO 5','SO 6','SO 7','SO 8'};

number outbounddist {Plants,CentralWH} = [
77	48	41	57	94
54	33	21	94	44
];

number inbounddist {CentralWH,RegionalWH} = [
62	55	114	106	51	76	122	80
63	92	90	50	88	47	136	47
76	85	56	65	129	40	89	100
42	128	74	70	45	56	81	131
77	135	118	59	124	93	41	94
];
number fixedcostWH {CentralWH} = [100 100 80 70 60];
number varcostWH {CentralWH} = [10 10 25 30 50];
number varcostplant {Plants} = [10 50];
number Demand {RegionalWH} = [121	96	166	175	210	130	180	160];
number Supply {Plants} = [1250 200];

var X {Plants,CentralWH} integer >=0;
var Z {CentralWH,RegionalWH} integer >=0;
var Y {CentralWH} binary;

min totalcost = sum {p in Plants,c in CentralWH} X[p,c]*1.5*outbounddist[p,c]+
                sum {p in Plants,c in CentralWH} X[p,c]*varcostplant[p]+
				sum {c in CentralWH,r in RegionalWH} Z[c,r]*1*inbounddist[c,r]+
				sum {c in CentralWH,r in RegionalWH} Z[c,r]*varcostWH[c]+
				sum {c in CentralWH} Y[c]*fixedcostWH[c];

 con supplycons {p in Plants}: sum {c in CentralWH} X[p,c]<=Supply[p];
 /*con supplycons1{p in Plants}: sum {c in CentralWH} X[p,'DC 1']>=1;*/
 /*con supplycons2{p in Plants}: sum {c in CentralWH} X[p,'DC 2']>=1;*/
 con demandcons {r in RegionalWH}: sum {c in CentralWH} Z[c,r]>=Demand[r];
 con nonstock   {c in CentralWH}: sum {p in Plants}X[p,c]-sum {r in RegionalWH}Z[c,r]=0;
 con conserva   {p in Plants,c in CentralWH}:X[p,c]-Y[c]*Supply[p]<=0;
 con numberofDC : sum {c in CentralWH} Y[c]>=1;

 solve;

 print X Y Z totalcost;

 quit;

[#11]

Sorry man for a lot of question.

Also i want to ask what if i want to add a constraint says that in case of demand shortage there will be a penality of X$ will be charged to the plants for each item shortage.

[#12]

To force DC 1 and DC 2 to be open, you can declare a constraint:

con FixOpen {c in {'DC 1','DC 2'}}:
   Y[c] = 1;

Or use the FIX statement:

for {c in {'DC 1','DC 2'}}
   fix Y[c] = 1;

To force exactly one potential DC to be open:

con ExactlyOnePotentialOpen:
   sum {c in {'Potential DC 3','Potential DC 4','Potential DC 5'}} Y[c] = 1;

Alternatively, you can scan for 'Potential' as the first word in the name:

con ExactlyOnePotentialOpen:
   sum {c in CentralWH: scan(c,1) = 'Potential'} Y[c] = 1;

For the other constraints you want, can you please write out the mathematical expressions?  The benefit of an optimization modeling language is that the corresponding constraint declarations will look very similar to the mathematics.