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
Posted by at

No comments:

Post a Comment

Subscribe to: Post Comments (Atom)