We revisit the elementary "WYNDOR GLASS CO." example and show

  • how to solve the dual problem explicitly,
  • how to access dual variables already provided by the solver.

We start as in the previous examples and set up the primal problem.

In [1]:
from pyomo.environ import *
from pyomo.opt import *
opt = solvers.SolverFactory("glpk")
In [2]:
model = ConcreteModel()

model.x = Var([1,2], within=NonNegativeReals)

model.c1 = Constraint(expr = model.x[1] <= 4)
model.c2 = Constraint(expr = 2*model.x[2] <= 12)
model.c3 = Constraint(expr = 3*model.x[1] + 2*model.x[2] <= 18)
model.z = Objective(expr = 3*model.x[1] + 5*model.x[2], sense=maximize)

We need to tell Pyomo that is should import dual variable values provided by the solver:

In [3]:
model.dual = Suffix(direction=Suffix.IMPORT)

Now we solve the problem and access the results as usual.

In [4]:
results = opt.solve(model)
In [5]:
model.x.get_values()
Out[5]:
{1: 2.0, 2: 6.0}
In [6]:
model.z.expr()
Out[6]:
36.0

Let us first solve the explicit dual, which we set up in the usual way.

In [7]:
dual = ConcreteModel()

dual.y = Var([1,2,3], within=NonNegativeReals)

dual.c1 = Constraint(expr = dual.y[1] + 3*dual.y[3] >= 3)
dual.c2 = Constraint(expr = 2*dual.y[2] + 2*dual.y[3] >= 5)
dual.z = Objective(expr = 4*dual.y[1] + 12*dual.y[2] + 18*dual.y[3], sense=minimize)

dual.dual = Suffix(direction=Suffix.IMPORT)
results = opt.solve(dual)

These are the values of the dual variables:

In [8]:
dual.y.get_values()
Out[8]:
{1: 0.0, 2: 1.5, 3: 1.0}

And this is the dual cost function. We see clearly that dual cost equals primal profit!

In [9]:
dual.z.expr()
Out[9]:
36.0

We can also ask the solver for the primal problem for dual variables ("shadow prices"). Note that there is one dual variable per constraint, so the dual variables are indexed by their corresponding constraint.

In [10]:
model.dual[model.c1], model.dual[model.c2], model.dual[model.c3]
Out[10]:
(0.0, 1.5, 1.0)

And we can ask for the dual variables of the dual problem:

In [11]:
dual.dual[dual.c1], dual.dual[dual.c2]
Out[11]:
(2.0, 6.0)

The dual of the dual is the primal problem!