Book distribution#
Click –> Live Code to activate live coding on this page!
Problem#
We’ll consider planning the shipment of books from distribution centers to stores where they are needed. There are three distribution centers at different cities: Amsterdam, Den Haag and Rotterdam. They have \(250\), \(130\) and \(235\) \(\text{books}\). There are four stores in Haarlem, Utrecht, Delft and Breda. They ordered \(70\), \(230\), \(240\) and \(70\) \(\text{books}\). All stores can receive books from all distribution centers, as shown below:
There are the following costs in \(\text{€ of transportation per book}\):
From \ To |
Haarlem |
Utrecht |
Delft |
Breda |
|---|---|---|---|---|
Amsterdam |
\( 7 \) |
\( 11 \) |
\( 16 \) |
\( 26 \) |
Den Haag |
\( 7 \) |
\( 13 \) |
\( 5 \) |
\( 20 \) |
Rotterdam |
\(16 \) |
\( 28 \) |
\( 7 \) |
\( 12 \) |
The goal is to minimize the shipping costs while meeting demand.
Model#
We need to define our model in the form of a linear constrained optimization model (2.1) and to apply scipy.optimize.linprog in matrix notation (2.2).
We’ll define the model as follows:
Design variables: how many books go from which distribution center to which store
Objective function: minimum shippings costs
Inequality constraint functions: don’t run out of stock in warehouses
Equality constraint functions: meet demand of stores
Bounds: book are only transported from distribution centers to stores.
Find best solution manually#
Try and adjust the book transports yourself. Can you find the optimal solution?
Design variables#
Let’s start with our design variables. In this case that could be the amount of books \(x\) transported from each distribution center \(i = \left[ 1 \left(\text{Amsterdam}\right),2 \left(\text{Den Haag}\right),3 \left( \text{Rotterdam}\right) \right] \) to each store \(j = \left[ 1 \left( \text{Haarlem}\right),2 \left(\text{Utrecht}\right),3 \left(\text{Delft}\right),4 \left( \text{Breda} \right)\right] \):
Alternatively, the design variables can be expressed as follows:
Objective function#
Now we can define the costs as the sum of the amount of books transported with the cost per transport as in (2.1)in the form of \(\mathop {\min }\limits_x f\left(x_{ij}\right) \):
This can be converted to matrix notation. In case of the design variables defined as in (2.2) in the form \(\mathop {\min }\limits_x {c^T}x\) with \(c\):
Inequality constraints#
Let’s continue with the inequality constraints, which should deal with that the stock of the distribution centers should always be bigger than \(0\). Or, stated differently, the sum of the amount of books transported out of each distribution center should be small or equal than the stock of each distribution center. These can be defined in the form \({{g}}\left(x_{ij}\right) \le 0\) as:
Or in matrix notation in the form \({A_{ub}}x \le {b_{ub}}\) as:
Equality constraints#
Finally, let’s define the inequality constraints, which should deal with that the each store receives the correct amount of books. Or, stated differently, the sum amount of paper books transported to each store should be equal to the demand of that store. These can be defined in the from \({h}\left(x_{ij}\right) = 0\) as:
Or in matrix notation in the form \({A_{eq}}x = {b_{eq}}\) as:
Bounds#
The number of books cannot be negative (assuming that book are only transported from distribution centers to stores). The maximum number of books transported could be the number of book available in a distribution center, but this is already defined by the constraint functions. Therefore, the bounds can be defined as:
Method#
Now let’s solve this problem using an optimization method.
Import libraries#
As this problem is higher-dimensional, we cannot easily plot the solution space. Therefore, we won’t import matplotlib.
import scipy as sp
import numpy as np
import scipy as sp
import numpy as np
Define variables#
As before, the (length of) design variable itself doesn’t have to be specified. So there’s actually nothing be be done here.
Define objective function#
The objective function was defined in (2.6), which gives:
c = np.array([7, 11, 16, 26, 7, 13, 5, 20, 16, 28, 7, 12])
Define constraint functions#
The constraint functions were defined in (2.8) and (2.10), which can be coded as follows:
A_ub = np.array([[1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1]])
b_ub = np.array([250, 130, 235])
A_eq = np.array([[1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0],
[0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0],
[0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0],
[0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1]])
b_eq = np.array([70, 230, 240, 70])
Bounds#
The bounds were defined in (2.11) and result in:
bounds = np.array([[0, None],
[0, None],
[0, None],
[0, None],
[0, None],
[0, None],
[0, None],
[0, None],
[0, None],
[0, None],
[0, None],
[0, None]])
Solve the problem#
result = sp.optimize.linprog(c, A_ub, b_ub, A_eq, b_eq, bounds)
print(result)
message: Optimization terminated successfully. (HiGHS Status 7: Optimal)
success: True
status: 0
fun: 5380.0
x: [ 2.000e+01 2.300e+02 0.000e+00 0.000e+00 5.000e+01
0.000e+00 8.000e+01 0.000e+00 0.000e+00 0.000e+00
1.600e+02 7.000e+01]
nit: 7
lower: residual: [ 2.000e+01 2.300e+02 0.000e+00 0.000e+00
5.000e+01 0.000e+00 8.000e+01 0.000e+00
0.000e+00 0.000e+00 1.600e+02 7.000e+01]
marginals: [ 0.000e+00 0.000e+00 1.100e+01 1.600e+01
0.000e+00 2.000e+00 0.000e+00 1.000e+01
7.000e+00 1.500e+01 0.000e+00 0.000e+00]
upper: residual: [ inf inf inf inf
inf inf inf inf
inf inf inf inf]
marginals: [ 0.000e+00 0.000e+00 0.000e+00 0.000e+00
0.000e+00 0.000e+00 0.000e+00 0.000e+00
0.000e+00 0.000e+00 0.000e+00 0.000e+00]
eqlin: residual: [ 0.000e+00 0.000e+00 0.000e+00 0.000e+00]
marginals: [ 9.000e+00 1.300e+01 7.000e+00 1.200e+01]
ineqlin: residual: [ 0.000e+00 0.000e+00 5.000e+00]
marginals: [-2.000e+00 -2.000e+00 -0.000e+00]
mip_node_count: 0
mip_dual_bound: 0.0
mip_gap: 0.0
Exercise#
Click –> Live Code to activate live coding on this page.