Mathematical Model
We use the following indices:
- \(p\): products
- \(f\): factories
- \(w\): warehouses
- \(s\): stores
- \(x_{p,f,w} \ge 0\): shipments of product \(p\) from factory \(f\) to warehouse \(w\),
- \(y_{s,w} \in \{0,1\}\): links each store \(s\) to a single warehouse \(w\).
- \(pcost_{f,p}\): unit production cost
- \(tcost_{p}\): unit transportation cost
- \(dist_{f,w}, dist_{s,w}\): distances
- \(pcap_{f,p}\): factory production capacities
- \(wcap_{w}\): warehouse capacity
- \(d_{s,p}\): demand for product \(p\) at store \(s\)
- \(turn_{p}\): product turnover rate
MIP Model |
---|
\[\large{\begin{align} \min&\sum_{p,f,w} (pcost_{f,p}+tcost_p\cdot dist_{f,w})\cdot x_{p,f,w}+\sum_{s,w,p} d_{s,p}\cdot tcost_p \cdot dist_{s,w} \cdot y_{s,w}\\ &\sum_w x_{p,f,w} \le pcap_{f,p} \>\>\forall f,p \>\> \text{(production capacity)}\\ &\sum_f x_{p,f,w} = \sum_s d_{s,p}\cdot y_{s,w} \>\>\forall p,w \>\>\text{(demand)}\\ &\sum_{p,s} \frac{d_{s,p}}{turn_{p}} y_{s,w} \le wcap_{w}\>\>\forall w \>\>\text{(warehouse capacity)}\\ &\sum_w y_{s,w} = 1\>\>\forall s \>\>\text{(one warehouse for a store)}\\ &x_{p,f,w} \ge 0 \\ &y_{s,w} \in \{0,1\} \end{align}}\] |
Matlab: Solver Based Implementation
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 | %----------------------------------------------------- % Generate a random problem: Facility Locations %----------------------------------------------------- rng(1) % for reproducibility N = 20; % N from 10 to 30 seems to work. Choose large values with caution. N2 = N*N; f = 0.05; % density of factories w = 0.05; % density of warehouses s = 0.1; % density of sales outlets F = floor(f*N2); % number of factories W = floor(w*N2); % number of warehouses S = floor(s*N2); % number of sales outlets xyloc = randperm(N2,F+W+S); % unique locations of facilities [xloc,yloc] = ind2sub([N N],xyloc); %----------------------------------------------------- % Generate Random Capacities, Costs, and Demands %----------------------------------------------------- P = 20; % 20 products % Production costs between 20 and 100 pcost = 80*rand(F,P) + 20; % Production capacity between 500 and 1500 for each product/factory pcap = 1000*rand(F,P) + 500; % Warehouse capacity between P*400 and P*800 for each product/warehouse wcap = P*400*rand(W,1) + P*400; % Product turnover rate between 1 and 3 for each product turn = 2*rand(1,P) + 1; % Product transport cost per distance between 5 and 10 for each product tcost = 5*rand(1,P) + 5; % Product demand by sales outlet between 200 and 500 for each % product/outlet d = 300*rand(S,P) + 200; %--------------------------------------------------------- % Generate Objective and Constraint Matrices and Vectors %--------------------------------------------------------- obj1 = zeros(P,F,W); % Allocate arrays obj2 = zeros(S,W); distfw = zeros(F,W); % Allocate matrix for factory-warehouse distances for ii = 1:F for jj = 1:W distfw(ii,jj) = abs(xloc(ii) - xloc(F + jj)) + abs(yloc(ii) ... - yloc(F + jj)); end end distsw = zeros(S,W); % Allocate matrix for sales outlet-warehouse distances for ii = 1:S for jj = 1:W distsw(ii,jj) = abs(xloc(F + W + ii) - xloc(F + jj)) ... + abs(yloc(F + W + ii) - yloc(F + jj)); end end for ii = 1:P for jj = 1:F for kk = 1:W obj1(ii,jj,kk) = pcost(jj,ii) + tcost(ii)*distfw(jj,kk); end end end for ii = 1:S for jj = 1:W obj2(ii,jj) = distsw(ii,jj)*sum(d(ii,:).*tcost); end end obj = [obj1(:);obj2(:)]; % obj is the objective function vector matwid = length(obj); Aineq = spalloc(P*F + W,matwid,P*F*W + S*W); % Allocate sparse Aeq bineq = zeros(P*F + W,1); % Allocate bineq as full % Zero matrices of convenient sizes: clearer1 = zeros(size(obj1)); clearer12 = clearer1(:); clearer2 = zeros(size(obj2)); clearer22 = clearer2(:); % First the production capacity constraints counter = 1; for ii = 1:F for jj = 1:P xtemp = clearer1; xtemp(jj,ii,:) = 1; % Sum over warehouses for each product and factory xtemp = sparse([xtemp(:);clearer22]); % Convert to sparse Aineq(counter,:) = xtemp'; % Fill in the row bineq(counter) = pcap(ii,jj); counter = counter + 1; end end % Now the warehouse capacity constraints vj = zeros(S,1); % The multipliers for jj = 1:S vj(jj) = sum(d(jj,:)./turn); % A sum of P elements end for ii = 1:W xtemp = clearer2; xtemp(:,ii) = vj; xtemp = sparse([clearer12;xtemp(:)]); % Convert to sparse Aineq(counter,:) = xtemp'; % Fill in the row bineq(counter) = wcap(ii); counter = counter + 1; end Aeq = spalloc(P*W + S,matwid,P*W*(F+S) + S*W); % Allocate as sparse beq = zeros(P*W + S,1); % Allocate vectors as full counter = 1; % Demand is satisfied: for ii = 1:P for jj = 1:W xtemp = clearer1; xtemp(ii,:,jj) = 1; xtemp2 = clearer2; xtemp2(:,jj) = -d(:,ii); xtemp = sparse([xtemp(:);xtemp2(:)]'); % Change to sparse row Aeq(counter,:) = xtemp; % Fill in row counter = counter + 1; end end % Only one warehouse for each sales outlet: for ii = 1:S xtemp = clearer2; xtemp(ii,:) = 1; xtemp = sparse([clearer12;xtemp(:)]'); % Change to sparse row Aeq(counter,:) = xtemp; % Fill in row beq(counter) = 1; counter = counter + 1; end %----------------------------------------------------- % Bound Constraints and Integer Variables %----------------------------------------------------- intcon = P*F*W+1:length(obj); lb = zeros(length(obj),1); ub = Inf(length(obj),1); ub(P*F*W+1:end) = 1; opts = optimoptions('intlinprog','Display','off','PlotFcn',@optimplotmilp); %----------------------------------------------------- % Solve the Problem %----------------------------------------------------- [solution,fval,exitflag,output] = intlinprog(obj,intcon,... Aineq,bineq,Aeq,beq,lb,ub,opts); |
Matlab: Problem Based Implementation
This implementation is using the new Matlab Problem Based Optimization language features [2]. In my opinion this is a big improvement. When reading this code we at least recognize the original problem. In addition we see the code is more compact.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 | %----------------------------------------------------- % Generate a random problem: Facility Locations %----------------------------------------------------- rng(1) % for reproducibility N = 20; % N from 10 to 30 seems to work. Choose large values with caution. N2 = N*N; f = 0.05; % density of factories w = 0.05; % density of warehouses s = 0.1; % density of sales outlets F = floor(f*N2); % number of factories W = floor(w*N2); % number of warehouses S = floor(s*N2); % number of sales outlets xyloc = randperm(N2,F+W+S); % unique locations of facilities [xloc,yloc] = ind2sub([N N],xyloc); %----------------------------------------------------- % Generate Random Capacities, Costs, and Demands %----------------------------------------------------- P = 20; % 20 products % Production costs between 20 and 100 pcost = 80*rand(F,P) + 20; % Production capacity between 500 and 1500 for each product/factory pcap = 1000*rand(F,P) + 500; % Warehouse capacity between P*400 and P*800 for each product/warehouse wcap = P*400*rand(W,1) + P*400; % Product turnover rate between 1 and 3 for each product turn = 2*rand(1,P) + 1; % Product transport cost per distance between 5 and 10 for each product tcost = 5*rand(1,P) + 5; % Product demand by sales outlet between 200 and 500 for each % product/outlet d = 300*rand(S,P) + 200; %----------------------------------------------------- % Generate Variables and Constraints %----------------------------------------------------- distfw = zeros(F,W); % Allocate matrix for factory-warehouse distances for ii = 1:F for jj = 1:W distfw(ii,jj) = abs(xloc(ii) - xloc(F + jj)) + abs(yloc(ii) ... - yloc(F + jj)); end end distsw = zeros(S,W); % Allocate matrix for sales outlet-warehouse distances for ii = 1:S for jj = 1:W distsw(ii,jj) = abs(xloc(F + W + ii) - xloc(F + jj)) ... + abs(yloc(F + W + ii) - yloc(F + jj)); end end x = optimvar('x',P,F,W,'LowerBound',0); y = optimvar('y',S,W,'Type','integer','LowerBound',0,'UpperBound',1); capconstr = sum(x,3) <= pcap'; demconstr = squeeze(sum(x,2)) == d'*y; warecap = sum(diag(1./turn)*(d'*y),1) <= wcap'; salesware = sum(y,2) == ones(S,1); %----------------------------------------------------- % Create Problem and Objective %----------------------------------------------------- factoryprob = optimproblem; objfun1 = sum(sum(sum(x,3).*(pcost'),2),1); objfun2 = 0; for p = 1:P objfun2 = objfun2 + tcost(p)*sum(sum(squeeze(x(p,:,:)).*distfw)); end r = sum(distsw.*y,2); % r is a length s vector v = d*(tcost(:)); objfun3 = sum(v.*r); factoryprob.Objective = objfun1 + objfun2 + objfun3; factoryprob.Constraints.capconstr = capconstr; factoryprob.Constraints.demconstr = demconstr; factoryprob.Constraints.warecap = warecap; factoryprob.Constraints.salesware = salesware; %----------------------------------------------------- % Solve the Problem %----------------------------------------------------- opts = optimoptions('intlinprog','Display','off','PlotFcn',@optimplotmilp); [sol,fval,exitflag,output] = solve(factoryprob,opts); |
GAMS Implementation
For completeness, here is the GAMS representation [3]:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 | sets k 'locations' /factory1*factory20, warehouse1*warehouse20, store1*store40/ f(k) 'factories' /factory1*factory20/ w(k) 'warehouses' /warehouse1*warehouse20/ s(k) 'stores' /store1*store40/ p 'products' /product1*product20/ ; parameters xloc(k) 'locations: x-coordinates' yloc(k) 'locations: y-coordinates' pcost(f,p) 'unit production cost' pcap(f,p) 'production capacity' wcap(w) 'warehouse capacity' turn(p) 'product turnover rate' tcost(p) 'unit transportation cost' d(s,p) 'demand' ; xloc(k) = uniformint(1,20); yloc(k) = uniformint(1,20); pcost(f,p) = uniform(20,100); pcap(f,p) = uniform(500,1500); wcap(w) = uniform(400*card(p),800*card(p)); turn(p) = uniform(1,3); tcost(p) = uniform(5,10); d(s,p) = uniform(200,500); positive variable x(p,f,w) 'flow factory to warehouse'; binary variable y(s,w) 'assign warehouse to store'; variable z 'objective variable'; equations obj 'objective' prodcap(f,p) 'production capacity' whcap(w) 'warehouse capacity' demand(p,w) assign(s) 'assign wh to store' ; obj.. z =e= sum((p,f,w), (pcost(f,p)+tcost(p)*distfw(f,w))*x(p,f,w)) + sum((s,w,p), d(s,p)*tcost(p)*distsw(s,w)*y(s,w)); prodcap(f,p).. sum(w,x(p,f,w)) =l= pcap(f,p); whcap(w).. sum((p,s),(d(s,p)/turn(p))*y(s,w)) =l= wcap(w); demand(p,w).. sum(f,x(p,f,w)) =e= sum(s, d(s,p)*y(s,w)); assign(s).. sum(w,y(s,w)) =e= 1; model m /all/; option optcr=0; solve m minimizing z using mip; |
It is instructive to see the differences in approaches. Here are some highlights:
Expression | Matlab (Model Based) | GAMS |
---|---|---|
\[\sum_{p,f,w} pcost_{f,p} \cdot x_{p,f,w}\] | sum(sum(sum(x,3).*(pcost'),2),1) |
sum((p,f,w), pcost(f,p)*x(p,f,w)) |
\[\sum_{p,f,w} tcost_p \cdot distfw_{f,w} \cdot x_{p,f,w}\] | objfun2 = 0; for p = 1:P objfun2 = objfun2 + tcost(p)*sum(sum(squeeze(x(p,:,:)).*distfw)); end |
sum((p,f,w), tcost(p)*distfw(f,w)*x(p,f,w)) |
\[\sum_w x_{p,f,w} \le pcap_{f,p}\] | sum(x,3) <= pcap' |
sum(w,x(p,f,w)) =l= pcap(f,p) |
\[\sum_f x_{p,f,w} = \sum_s d_{s,p}\cdot y_{s,w}\] | squeeze(sum(x,2)) == d'*y |
sum(f,x(p,f,w)) =e= sum(s, d(s,p)*y(s,w)) |
\[\sum_{p,s} \frac{d_{s,p}}{turn_{p}} y_{s,w} \le wcap_{w}\] | sum(diag(1./turn)*(d'*y),1) <= wcap' |
sum((p,s),(d(s,p)/turn(p))*y(s,w)) =l= wcap(w) |
\[\sum_w y_{s,w} = 1\] | sum(y,2) == ones(S,1) |
sum(w,y(s,w)) =e= 1 |
There is a notable difference in number of different language constructs we use. In Matlab we use:
- Different versions of sum
- A loop
- Diag, transpose ('), element wise division (./)
- Squeeze
In GAMS we basically only use sum.
References
- Factory, Warehouse, Sales Allocation Model: Solver-Based, http://www.mathworks.com/help/optim/ug/factory-warehouse-sales-allocation-model.html
- Factory, Warehouse, Sales Allocation Model: Model-Based, http://www.mathworks.com/help/optim/examples/factory-warehouse-sales-allocation-model.html
- Matlab vs GAMS: Integer Programming, http://yetanothermathprogrammingconsultant.blogspot.com/2016/10/matlab-vs-gams-integer-programming.html