Theory of LP and the Simplex method¶

Colab notebook

First, we will briefly discuss the primal-dual theory with a few examples. Consider the Factory problem that was discussed earlier:

\[\begin{split} \begin{align} &\text{max}&&\\ &\qquad z=40𝑥_1 + 50𝑥_2\\ &\text{subject to}&&\\ &\qquad 1x_1 + 2𝑥_2 \le 40\\ &\qquad 4𝑥_1 + 3𝑥_2 \le 120\\ &\qquad x_1, x_2 \ge 0.\\ \end{align} \end{split}\]

!pip install gurobipy
from gurobipy import *
m = Model()
x1 = m.addVar(lb=0, vtype = GRB.CONTINUOUS, name='x1') 
x2 = m.addVar(lb=0, vtype = GRB.CONTINUOUS, name='x2')
m.setObjective(40*x1+50*x2, GRB.MAXIMIZE)
m.addConstr(1*x1+2*x2<=40, name='c1')
m.addConstr(4*x1+3*x2<=120, name= 'c2')
m.optimize()
print('*'*100)
for var in m.getVars(): # descision variable
    print(var.varName, '=', var.x, (var.obj,var.SAObjLow, var.SAObjUp, var.RC))
print('*'*100)
print('optimal total revenue:', m.objVal)
print('*'*100)
for con in m.getConstrs(): # constraints
    print(con.ConstrName, ': slack =', con.slack,', shadow price=',
          con.pi,',', (con.RHS, con.SARHSLow, con.SARHSUp))
print('*'*100)
print('*'*100)

Collecting gurobipy
  Downloading gurobipy-9.5.0-cp37-cp37m-manylinux2014_x86_64.whl (11.5 MB)
     |████████████████████████████████| 11.5 MB 4.8 MB/s 
?25hInstalling collected packages: gurobipy
Successfully installed gurobipy-9.5.0
Restricted license - for non-production use only - expires 2023-10-25
Gurobi Optimizer version 9.5.0 build v9.5.0rc5 (linux64)
Thread count: 1 physical cores, 2 logical processors, using up to 2 threads
Optimize a model with 2 rows, 2 columns and 4 nonzeros
Model fingerprint: 0x3a526911
Coefficient statistics:
  Matrix range     [1e+00, 4e+00]
  Objective range  [4e+01, 5e+01]
  Bounds range     [0e+00, 0e+00]
  RHS range        [4e+01, 1e+02]
Presolve time: 0.01s
Presolved: 2 rows, 2 columns, 4 nonzeros

Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    9.0000000e+31   3.250000e+30   9.000000e+01      0s
       2    1.3600000e+03   0.000000e+00   0.000000e+00      0s

Solved in 2 iterations and 0.02 seconds (0.00 work units)
Optimal objective  1.360000000e+03
****************************************************************************************************
x1 = 24.0 (40.0, 25.0, 66.66666666666667, 0.0)
x2 = 8.0 (50.0, 30.0, 80.0, 0.0)
****************************************************************************************************
optimal total revenue: 1360.0
****************************************************************************************************
c1 : slack = 0.0 , shadow price= 16.0 , (40.0, 30.0, 80.0)
c2 : slack = 0.0 , shadow price= 6.0 , (120.0, 60.0, 160.0)
****************************************************************************************************
****************************************************************************************************

Here, the optimal solution is

\[ (x^*,z)=((24,8),1360). \]

In addition, the Python solution shows that

\[\begin{split} \begin{align} \qquad c_1 &= 40~(25, 66.667),\\ \qquad c_2 &= 50~(30, 80),\\ \qquad b_1 &= 40~(30, 80),~\text{with shadow price } 16, \\ \qquad b_2 &= 120~(60, 160),~\text{with shadow price } 6,\\ \qquad s_1 &= 0, \\ \qquad s_2 &= 0. \\ \end{align} \end{split}\]

That is, the cost coefficients, \(\mathbf c = (c_1, c_2)\), have values 40 and 50, and the ranges in which they are allowed to change without affecting the optimal \(x^*\) are \((25, 66.667)\) and \((25, 66.667)\). Similarly, the RHS constraints, \(\mathbf b = (b_1, b_2)\), have values of 40 and 120 and can change to values \((30, 80)\) and \((60, 160)\) without affecting the optimal solution mix. Also, the shadow prices of the RHS constraints are 16 and 6, and there is no slack.

Is it possible to infer the optimal \(z\) or at least bound its value without solving the LP model?

The above objective is \(40𝑥_1 + 50𝑥_2\) can be bounded using the first constraint and the fact that both decision variables are non-negative

\[ 40𝑥_1 + 50𝑥_2<=40\times(x_1 + 2𝑥_2)=40x_1 + 80𝑥_2\le 1600. \]

Can we do better?

Let’s use the second constraint,

\[ 40𝑥_1 + 50𝑥_2<=50/3\times(4𝑥_1 + 3𝑥_2)=66.67x_1 + 50𝑥_2=2000 \]

but, here, the value 2000 is higher than the best upper bound we have so far, which is 1600.

Systematically, we can write that

\[ 40𝑥_1 + 50𝑥_2\le d_1x_1 +d_2x_2 \le h, \]

and let \(h\) be the upper bound on the maximum of the objective. The trick is that we will use the constraint equations to infer \(d_1, d_2\) and \(h\). That is, we multiply the first constraint by \(v_1\ge0\), the second by \(v_2\ge0\), and then add the two:

\[ v_1(1x_1 + 2𝑥_2)+v_2(4𝑥_1 + 3𝑥_2)\le 40v_1+120v_2 \]

or

\[ (v_1+4v_2)x_1+(2v_2+3v_2)x_2\le 40v_1+120v_2. \]

In the above notation:

\[\begin{split} \begin{align} &d_1=v_1+4v_2, \\ &d_2=2v_2+3v_2,\qquad \text{and}\\ &h=40v_1+120v_2. \end{align} \end{split}\]

How do we choose the best coefficients \(v_1\), and \(v_2\)? We must ensure that \(d_1\ge 40\) and \(d_2\ge 50\), and we want \(h\) to be as small as possible under these constraints. This is again an LP model which is called the dual to the primal set

\[\begin{split} \begin{align} &\text{min}\\ &\qquad h=40v_1 + 120v_2\\ &\text{s.t.}\\ &\qquad 1v_1 + 4v_2 \ge 40\\ &\qquad 2v_1 + 3v_2 \ge 50\\ &\qquad v_1, v_2 \ge 0. \end{align} \end{split}\]

In general, the dual to primal LP is another LP model that is derived from it in the following way:

Each variable in the primal becomes a constraint in the dual
Each constraint in the primal becomes a variable in the dual
The objective direction is inversed: maximum in the primal becomes minimum in the dual, and vice versa.

Hence, for the max primal

\[\begin{split} \begin{align} &\max \\ &\qquad\mathbf c\cdot\mathbf x,\\ &\text{s.t.}\\ &\qquad\mathbf A\mathbf x\le \mathbf b,\\ &\qquad\mathbf x \ge 0, \end{align} \end{split}\]

the corresponding dual, is

\[\begin{split} \begin{align} &\min \\ &\qquad\mathbf b\cdot\mathbf v,\\ &\text{s.t.}\\ &\qquad\mathbf A^T\mathbf v\ge \mathbf c,\\ &\qquad\mathbf v\ge 0. \end{align} \end{split}\]

The interpretation is that we solve for \(\mathbf v\), the shadow prices of the primal, by constraining the shadow prices with the cost coefficients, \(\mathbf c\).

Solving for \(v\) using Python, we find that the optimal is

\[ (v^*,h)=((16,6),1360). \]

m = Model()    
x1 = m.addVar(lb=0, vtype = GRB.CONTINUOUS, name='v1') 
x2 = m.addVar(lb=0, vtype = GRB.CONTINUOUS, name='v2')
m.setObjective(40*x1+120*x2, GRB.MINIMIZE)
m.addConstr(1*x1+4*x2>=40, name='c1')
m.addConstr(2*x1+3*x2>=50, name= 'c2')
m.optimize()
print('*'*100)
for var in m.getVars(): # descision variable
    print(var.varName, '=', var.x, (var.obj,var.SAObjLow, var.SAObjUp, var.RC))
print('*'*100)
print('optimal total revenue:', m.objVal)
print('*'*100)
for con in m.getConstrs(): # constraints
    print(con.ConstrName, ': slack =', con.slack,', shadow price=',
          con.pi,',', (con.RHS, con.SARHSLow, con.SARHSUp))

Gurobi Optimizer version 9.5.0 build v9.5.0rc5 (linux64)
Thread count: 1 physical cores, 2 logical processors, using up to 2 threads
Optimize a model with 2 rows, 2 columns and 4 nonzeros
Model fingerprint: 0x8c4006b8
Coefficient statistics:
  Matrix range     [1e+00, 4e+00]
  Objective range  [4e+01, 1e+02]
  Bounds range     [0e+00, 0e+00]
  RHS range        [4e+01, 5e+01]
Presolve time: 0.01s
Presolved: 2 rows, 2 columns, 4 nonzeros

Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    0.0000000e+00   4.500000e+01   0.000000e+00      0s
       2    1.3600000e+03   0.000000e+00   0.000000e+00      0s

Solved in 2 iterations and 0.02 seconds (0.00 work units)
Optimal objective  1.360000000e+03
****************************************************************************************************
v1 = 16.0 (40.0, 30.0, 80.0, 0.0)
v2 = 6.0 (120.0, 60.0, 160.0, 0.0)
****************************************************************************************************
optimal total revenue: 1360.0
****************************************************************************************************
c1 : slack = 0.0 , shadow price= 24.0 , (40.0, 25.0, 66.66666666666666)
c2 : slack = 0.0 , shadow price= 8.0 , (50.0, 30.0, 80.0)

In addition, as shown above

\[\begin{split} \begin{align} \qquad b_1 &= 40~(30,80),\\ \qquad b_2 &= 120~(60,160),\\ \qquad c_1 &= 40~(25, 66.667),~\text{with shadow price } 24, \\ \qquad c_2 &= 50~(30, 80),~\text{with shadow price } 8,\\ \qquad s_1 &= 0, \\ \qquad s_2 &= 0. \\ \end{align} \end{split}\]

The dual’s decision variables, \(\mathbf v\), are the primal’s shadow prices and the dual’s \(\mathbf b\) and \(\mathbf c\) correspond with their primal values. Lastly, the dual’s shadow prices are the primal’s decision variables.

The primal-dual correspondence gives us more flexibility in solving the LP model. In cases where the dual is simpler, we can solve it instead of the primal.

Few other properties emerge from the primal-dual relationship:

OR/MS in Python

Theory of LP and the Simplex method¶

Weak duality¶

Complementary slackness¶

KKT conditions for optimality¶

Improving search¶

A few remarks:¶

The Simplex Algorithm¶