## Sunday, March 4, 2018

### Matlab: Problem Based Optimization

MATLAB has a new way to express optimization problems. Until now they used a matrix-based API (this is now called: "Solver Based Optimization"). In 2017b they added "Problem Based Optimization" which is kind of a algebraic modeling language (like AMPL or GAMS) inside Matlab. Let's revisit a simple problem related to warehousing [3].

#### Mathematical Model

We use the following indices:
• $$p$$: products
• $$f$$: factories
• $$w$$: warehouses
• $$s$$: stores
We introduce the following variables:
• $$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$$.
The data associated with the model is:
• $$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
The optimization model looks like:

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

This implementation uses the old matrix based API. The code is somewhat removed from the original problem as we need to use quite some non-trivial code to deal with the details in populating the solver matrices. This is still a very simple model: for more complex models this amount of code and its complexity will only increase. See [1] for more details.

  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 60 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' distfw(f,w) 'distance factory-warehouse' distsw(s,w) 'distance store-warehouse' ; 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); distfw(f,w) = abs(xloc(f)-xloc(w))+abs(yloc(f)-yloc(w)); distsw(s,w) = abs(xloc(s)-xloc(w))+abs(yloc(s)-yloc(w)); 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; 

There is a notable difference in number of different language constructs we use. In Matlab we use:

1. Different versions of sum
2. A loop
3. Diag, transpose ('), element wise division (./)
4. Squeeze
In GAMS we basically only use sum.

#### References

1. Factory, Warehouse, Sales Allocation Model: Solver-Based, http://www.mathworks.com/help/optim/ug/factory-warehouse-sales-allocation-model.html
2. Factory, Warehouse, Sales Allocation Model: Model-Based, http://www.mathworks.com/help/optim/examples/factory-warehouse-sales-allocation-model.html
3. Matlab vs GAMS: Integer Programming, http://yetanothermathprogrammingconsultant.blogspot.com/2016/10/matlab-vs-gams-integer-programming.html