Objectives
-
Formulate and solve dual problems using Julia
-
Extract basis information from solutions using Julia
Section 6.3 Solving Dual Problems with Julia
In Section 6.2, one main takeaway was the fact that solving the primal problem also gives us a solution to the dual problem. Here, we will see how to explicitly formulate and solve dual problems using Julia.
Subsection 6.3.1 Finding Dual Values and Certificates using Julia
m = Model(HiGHS.Optimizer)
@variable(m, x₁ ≥ 0)
@variable(m, x₂ ≥ 0)
@variable(m, x₃ ≥ 0)
@objective(m, Max, 3x₁ + 4x₂ + 7x₃)
@constraint(m, c1, x₁ + 2x₂ + 3x₃ ≤ 12)
@constraint(m, c2, 4x₁ + 3x₂ + 2x₃ ≤ 18)
print(m)
And as we will see below, including the terms
c1 and c2 in the constraints is important. The result of this is the LOP as expected.
Max 3 x₁ + 4 x₂ + 7 x₃ Subject to c1 : x₁ + 2 x₂ + 3 x₃ ≤ 12 c2 : 4 x₁ + 3 x₂ + 2 x₃ ≤ 18 x₁ ≥ 0 x₂ ≥ 0 x₃ ≥ 0
And then solving this with the following:
set_silent(m)
optimize!(m)
is_solved_and_feasible(m)
The results is
true and finding the optimal values of the variables and objective with:
value(x₁), value(x₂), value(x₃)
results in \((3.0, 0.0, 3.0)\text{,}\) and for the objective use
objective_value(m) which results in \(30.0\text{.}\)
Since we have also used constants
c1 and c2 in the constraints, we can now get the values of these as well with value(c1) and value(c2). These return 12.0 and 18.0 respectively. These are the values of the left-hand sides of the constraints at optimality also since they equal the right-hand sides, the slack variables for each of these are 0.0.
As we saw in Section 6.2, the solution to the dual is in the objective row of the optimal tableau. However, it doesn’t appear that we have the tableau. Let’s first solve the dual problem from (6.2.3) explicitly with the following:
md = Model(HiGHS.Optimizer)
@variable(md, y₁ ≥ 0)
@variable(md, y₂ ≥ 0)
@objective(md, Min, 12y₁ + 18y₂)
@constraint(md, c1, y₁ + 4y₂ ≥ 3)
@constraint(md, c2, 2y₁ + 3y₂ ≥ 4)
@constraint(md, c3, 3y₁ + 2y₂ ≥ 7)
print(md)
and the result is
Min 12 y₁ + 18 y₂ Subject to c1 : y₁ + 4 y₂ ≥ 3 c2 : 2 y₁ + 3 y₂ ≥ 4 c3 : 3 y₁ + 2 y₂ ≥ 7 y₁ ≥ 0 y₂ ≥ 0
And then solving this with the following:
set_silent(md)
optimize!(md)
is_solved_and_feasible(md)
The results is
true and finding the optimal values of the variables and objective with:
value(y₁), value(y₂)
results in \((2.1999999999999997, 0.2)\text{,}\) and for the objective use
objective_value(md) which results in \(30.0\text{,}\) which is the same value as that of the primal problem. (Yes, the strong duality theorem comes through.) Again, we can find the values of the constraints with value(c1), value(c2), and value(c3), which return 3.0, 5.0, and 7.0 respectively. Since the first and third values are the same as the right-hand sides, the slack variables for these constraints are 0.0, while the second constraint has a slack variable value of 1.0.
It seems unnecessary to have to solved the dual problem in light of the work in Section 6.2 and its true that we don’t need to. If instead, we enter the following code:
dual(c1),dual(c2)
The result is
(-2.2, -0.19999999999999973). The dual function returns the value of the slack variable coefficients in the objective function written in dictionary form. That is, returning to the final tableau in (6.2.2), the objective row can be written as
\begin{equation*}
\begin{aligned} 10 z \amp = 300 - 10 x_2 - 22x_4 - 2x_5 \amp\amp \\ \amp\amp\amp\text{or} \\ z \amp = 30 - x_2 - \frac{22}{10} x_4 - \frac{2}{10} x_5 \end{aligned}
\end{equation*}
The values of the dual variables using Julia results in the coefficients of the slack variables that correspond to the two constraints and in fact the result of the
dual command is the negative of the values of the dual variables found directly.
You may notice that
dual gave the coefficients of \(x_4\) and \(x_5\text{,}\) but not \(x_2\text{.}\) We need a little more work for that.
Subsection 6.3.2 Dual Values of Regular Variables in Julia
The dual values of the regular variables (those that are not slack or surplus variables) can also be found using Julia. These correspond to the constraints in the dual problem. To find these, we can use the following code:
reduced_cost(x₁), reduced_cost(x₂), reduced_cost(x₃)
which results in
(-0.0, -1.0, -0.0). This is again the negative of the dual variable values for the three constraints in the dual problem. These provide the certificate of optimality for the primal problem as we saw in Section 6.2 in terms of the objective row of the final tableau (to the left of the vertical line).
Subsection 6.3.3 Finding a basis using Julia
The primal and dual example that we saw above showed that we can solve the primal problem and get the values of the dual problem. However, another useful part of the solution that we have seen over the first few chapters of this book is the basis of the solution.
1
The basis will also play a significant role in determining uniqueness of solutions.
Although the basis is important in the simplex method, often packages like
HiGHS may not use the simplex tableau as we have in this book. We provide a possible alternative in Chapter 8. However, with some work, we can get the basis out of the solution using JuMP and HiGHS.
The following code will show which of the original variables are basic:
[
xi => get_attribute(xi, MOI.VariableBasisStatus()) for
xi in all_variables(m)
]
and for the Primal problem above (this is the
m model), this results in
3-element Vector{Pair{VariableRef, MathOptInterface.BasisStatusCode}}:
x₁ => MathOptInterface.BASIC
x₂ => MathOptInterface.NONBASIC_AT_LOWER
x₃ => MathOptInterface.BASIC
which shows that \(x_1\) and \(x_3\) are basic and \(x_2\) is non-basic. The
NONBASIC_AT_LOWER indicates that the variable is basic and is at its lower bound of zero.
Similarly, we can get the basis status of the constraints, which correspond to the dual variables, with the following code:
[ name(c) => MOI.get(m, MOI.ConstraintBasisStatus(), c)
for c in all_constraints(m, include_variable_in_set_constraints = false)]
which returns:
2-element Vector{Pair{String, MathOptInterface.BasisStatusCode}}:
"x₄" => MathOptInterface.NONBASIC
"x₅" => MathOptInterface.NONBASIC
Indicating that both constraints are non-basic.
