Sunday, October 2, 2016

Matlab vs GAMS: Integer Programming

In a previous post I translated a small GAMS model (model number one from the GAMS model library) into Matlab and compared the codes side-by-side. In this post we take slightly more complicated MIP model written in Matlab and convert it into GAMS.   

The problem discussed here (1) involves a network of factories, warehouses and sales outlets. We need to find the least expensive flow of products from factories to warehouses to stores. One particularity is that each store gets its products from one warehouse.

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:

\[\boxed{\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 vs GAMS Implementation

Here we compare the Matlab implementation from (1) to a GAMS transcription of the model.

Matlab code GAMS code

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

 
 

sets
  k
'locations' /factory1*factory20,
           warehouse1*warehouse20,store1*store40/

  f(k)
'factories' /factory1*factory20/
  w(k)
'warehouses' /warehouse1*warehouse20/
  s(k)
'stores' /store1*store40/

;

xyloc = randperm(N2,F+W+S); % unique locations of facilities

[xloc,yloc] = ind2sub([N N],xyloc);
Pick unique locations on a grid
Not 100% the same: just random locations. No guarantee of uniqueness.

parameters xloc(k),yloc(k);
xloc(k) = uniformint(1,20);
yloc(k) = uniformint(1,20);

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;

 
 

set p 'products' /product1*product20/;
parameters

   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'
;
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 = 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

 
The indexing is cleaner thanks to \(f,s,w\) being subsets of \(k\).

parameters
   distfw(f,w)
'distance factory-warehouse'
   distsw(s,w)
'distance store-warehouse'

;
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));

obj1 = zeros(P,F,W); % Allocate arrays

obj2 = zeros(S,W);

% Generate the entries of obj1 and obj2.

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

% Combine the entries into one vector.

obj = [obj1(:);obj2(:)]; % obj is the objective function vector

 
  positive variable x(p,f,w) 'flow factory to warehouse';
binary variable y(s,w) 'assign warehouse to store'
;
variable z 'objective variable'
;

equation obj 'objective'
;
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));

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

Matrix Aineq holding the coefficients for the inequality constraints is large and sparse, so use a sparse matrix instead of a dense one.   
 

equation prodcap(f,p) 'production capacity';
prodcap(f,p)..
   
sum
(w,x(p,f,w)) =l= pcap(f,p);

% 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

 
 

equation whcap(w) 'warehouse capacity';
whcap(w)..
 
sum
((p,s),(d(s,p)/turn(p))*y(s,w)) =l= wcap(w);

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

Matrix Aeq is for the coefficients of the equality constraints. Also stored as a sparse matrix
 

equation demand(p,w);
demand(p,w)..
  
sum(f,x(p,f,w)) =e= sum
(s, d(s,p)*y(s,w));

% 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

 
 

equation assign(s) 'assign wh to store';
assign(s)..
  
sum
(w,y(s,w)) =e= 1;

intcon = P*F*W+1:length(obj);

lb = zeros(length(obj),1);

ub = Inf(length(obj),1);

ub(P*F*W+1:end) = 1;

[solution,fval,exitflag,output] = intlinprog(obj,intcon,Aineq,bineq,Aeq,beq,lb,ub);

 
 

model m /all/;
option
optcr=0;
solve
m minimizing z using mip;

The conclusion is that when we compare an equation based modeling system compared to a matrix based system, the equation based formulation is really more compact and readable when the model equations become more complex.

References