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.

from pyomo.environ import *
from pyomo.opt import *
opt = solvers.SolverFactory("glpk")

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:

model.dual = Suffix(direction=Suffix.IMPORT)

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

results = opt.solve(model)
model.dual[model.c1],model.dual[model.c2],model.dual[model.c3]

(0.0, 1.5, 1.0)

model.x.get_values()

{1: 2.0, 2: 6.0}

model.z.expr()

36.0

Now suppose we can increase the production capacity at either one of the three plants by one unit. Which one should we choose in order to increase the profit the most?

Answer: The shadow price associated to the second constraint (plant 2) is the highest, so we should choose this one. Increasing the b_2 from 12 to 13 will increase the profit from 36 to 37.5. (Check this directly by changing the code above.)

---

Next, let us solve the explicit dual, which we set up in the usual way.

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:

dual.y.get_values()

{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!

dual.z.expr()

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.

model.dual[model.c1], model.dual[model.c2], model.dual[model.c3]

(0.0, 1.5, 1.0)

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

dual.dual[dual.c1], dual.dual[dual.c2]

(2.0, 6.0)

The dual of the dual is the primal problem!