Solving a simple transportation problem using CVXOPT

12 November 2019

Problem Description

Assume that we are the manager of a support chain. Our company has two factory A and B, and each of them has 300 and 250 products, respectively. Now we would like to deliver them two three retail stores in different cities 1, 2, 3. There is a cost for each delivery as listed below.

  Store 1 Store 2 Store 3
Factory A 5 6 4
Factory B 6 3 7

Problem: How to design our transpotation strategy so that we can minimize the cost.


1. Formulate the linear programing(LP) problem.

Let’s first try to write down our problem as a optimization problem.

Let $x = [x_1, x_2, … , x_6]$ be the number of clothes delivering from source to destination cities, respectively.

Maybe let’s make it more clear:

\[x_1: A\to1 \\ x_2: A\to2 \\ x_3: A\to3 \\ x_4: B\to1 \\ x_5: B\to2 \\ x_6: B\to3\]

Our optimization problem can be formulate as:

\[\begin{align} minimize \ \ \ &5x_1 + 6x_2 + 4x_3 + 6x_4 + 3x_5 + 7x_6\\ s.t. \ \ \ & x_1 + x_4 = 200 \\ & x_2 + x_5 = 300 \\ & x_3 + x_6 = 250 \\ & x_1 + x_2 + x_3 \le 300 \\ & x_4 + x_5 + x_6 \le 500 \\ & x_1, x_2, x_3, x_4, x_5, x_6 \ge 0 \\ \end{align}\]

Now we can see that this is actually a Linear Programming (LP) problem, we will try to use the LP solver in CVXOPT package to solve it.

2. Convert to standard form.

We first check the documentation of CVXOPT’s LP solver, and we found that we need to convert the problem to the standard form before calling the function. CVXOPT’s LP solver uses the following standard form:

\[\begin{align} minimize \ \ \ &c^Tx\\ s.t. \ \ \ & Gx + s \le h \\ & Ax = b \\ &s \succeq 0 \end{align}\]

Hence, we need to rewrite our equation sets to match this standard form.

\[\begin{align} minimize \ \ \ &5x_1 + 6x_2 + 4x_3 + 6x_4 + 3x_5 + 7x_6\\ s.t. \ \ \ & x_1 + x_2 + x_3 \le 300 \\ & x_4 + x_5 + x_6 \le 500 \\ & -x_1 \le 0 \\ & -x_2 \le 0 \\ & -x_3 \le 0 \\ & -x_4 \le 0 \\ & -x_5 \le 0 \\ & -x_6 \le 0 \\ & x_1 + x_4 = 200 \\ & x_2 + x_5 = 300 \\ & x_3 + x_6 = 250 \\ \end{align}\]

Therefore, we have:

\[c = [5,6,4,6,3,7]\] \[G = \begin{vmatrix} 1& 1& 1& 0& 0& 0 \\ 0& 0& 0& 1& 1& 1 \\ -1& 0& 0& 0& 0& 0 \\ 0& -1& 0& 0& 0& 0 \\ 0& 0& -1& 0& 0& 0 \\ 0& 0& 0& -1& 0& 0 \\ 0& 0& 0& 0& -1& 0 \\ 0& 0& 0& 0& 0& -1 \end{vmatrix}\] \[h = [300, 500, 0,0,0,0,0,0]^T\] \[A = \begin{vmatrix} 1& 0& 0& 1& 0& 0 \\ 0& 1& 0& 0& 1& 0 \\ 0& 0& 1& 0& 0& 1 \end{vmatrix}\] \[b = [200, 300, 250]^T\]

3. Use CVXOPT to solve the LP problem.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
from cvxopt import matrix, solvers
import math
c = matrix([5.,6.,4.,6.,3.,7.])
G = matrix([[1., 0.,-1., 0., 0., 0., 0., 0.],
           [1., 0., 0.,-1., 0., 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., 1., 0., 0., 0., 0., 0.,-1.],])
h = matrix([300., 500., 0.,0.,0.,0.,0.,0.])
A = matrix([[1., 0., 0.],
           [0., 1., 0.],
           [0., 0., 1.],
           [1., 0., 0.],
           [0., 1., 0.],
           [0., 0., 1.],])
b = matrix([200., 300., 250.])
sol = solvers.lp(c, G, h, A, b)
print(sol['x'])
for i in range(6):
    print('x{}={}'.format(i+1, round(sol['x'][i])))

Output:

     pcost       dcost       gap    pres   dres   k/t
 0:  3.8500e+03  1.7500e+03  2e+03  5e-03  2e-16  1e+00
 1:  3.1521e+03  2.7932e+03  4e+02  9e-04  6e-15  1e+01
 2:  3.0593e+03  3.0367e+03  2e+01  5e-05  4e-16  1e+00
 3:  3.0501e+03  3.0498e+03  3e-01  7e-07  3e-16  1e-02
 4:  3.0500e+03  3.0500e+03  3e-03  7e-09  3e-16  1e-04
Optimal solution found.
[ 5.00e+01]
[ 2.92e-05]
[ 2.50e+02]
[ 1.50e+02]
[ 3.00e+02]
[ 4.59e-04]

x1=50
x2=0
x3=250
x4=150
x5=300
x6=0

4. Convert the LP problem to dual form

Now we take one step further. Notice that in the CVXOPT documentation, this LP solver actually solves the pair of primal and dual linear programs and it also shows the duel problem on the page. Now we will try to do this convertion by ourselves.

Recall that the original minimization problem is a linear programing problem with inequality constraints, its standard form is:

\[\begin{align} minimize: &c^Tx\\ s.t.: & Gx \le h \\ & Ax = b \\ \end{align}\]

Based on Lagrangian, we have:

\[\begin{align} L(x, \lambda, \mu) &= c^Tx + \lambda^T(Gx-h) + \mu^T(Ax - b) \\ &=c^Tx + \lambda^TGx - \lambda^Th + \mu^TAx-\mu^Tb \\ &=-\lambda^Th-\mu^Tb + (c + G^T\lambda + A^T\mu)^Tx \end{align}\]

where $\lambda \ge 0$.

Now we construct the dual form:

\[\begin{align} g(\lambda, \mu) &= \inf_{x}L(x, \lambda, \mu) \\ &= \inf_{x}(-\lambda^Th-\mu^Tb + (c + G^T\lambda + A^T\mu)^Tx) \\ &= \begin{cases} -\lambda^Th-\mu^Tb, & \text{if}\ c + G^T\lambda + A^T\mu =0 \\ -\infty, & \text{otherwise} \end{cases} \end{align}\]

Since the problem is only meaningful when $c + G^T\lambda + A^T\mu =0$, we can formulate the dual problem as:

\[\begin{align} maximize: &-h^T\lambda-b^T\mu\\ s.t.: & G^T\lambda + A^T\mu + c = 0 \\ & \lambda \ge 0 \\ \end{align}\]