# Sensitivity analysis of a linear program

This tutorial explains how to use the lp_sensitivity_report function to create sensitivity reports like those that are produced by the Excel Solver. This is most often used in introductory classes to linear programming.

In brief, sensitivity analysis of a linear program is about asking two questions:

1. Given an optimal solution, how much can I change the objective coefficients before a different solution becomes optimal?
2. Given an optimal solution, how much can I change the right-hand side of a linear constraint before a different solution becomes optimal?

JuMP provides a function, lp_sensitivity_report, to help us compute these values, but this tutorial extends that to create more informative tables in the form of a DataFrame.

## Setup

This tutorial uses the following packages:

using JuMP
import HiGHS
import DataFrames

as well as this small linear program:

model = Model(HiGHS.Optimizer)
@variable(model, x >= 0)
@variable(model, 0 <= y <= 3)
@variable(model, z <= 1)
@objective(model, Min, 12x + 20y - z)
@constraint(model, c1, 6x + 8y >= 100)
@constraint(model, c2, 7x + 12y >= 120)
@constraint(model, c3, x + y <= 20)
optimize!(model)
solution_summary(model; verbose = true)
* Solver : HiGHS

* Status
Termination status : OPTIMAL
Primal status      : FEASIBLE_POINT
Dual status        : FEASIBLE_POINT
Result count       : 1
Has duals          : true
Message from the solver:
"kHighsModelStatusOptimal"

* Candidate solution
Objective value      : 2.04000e+02
Objective bound      : 0.00000e+00
Relative gap         : Inf
Dual objective value : 2.04000e+02
Primal solution :
x : 1.50000e+01
y : 1.25000e+00
z : 1.00000e+00
Dual solution :
c1 : 2.50000e-01
c2 : 1.50000e+00
c3 : 0.00000e+00

* Work counters
Solve time (sec)   : 2.74897e-04
Simplex iterations : 2
Barrier iterations : 0
Node count         : -1


Can you identify:

• The objective coefficient of each variable?
• The right-hand side of each constraint?
• The optimal primal and dual solutions?

## Sensitivity reports

Now let's call lp_sensitivity_report:

report = lp_sensitivity_report(model)
SensitivityReport(Dict{ConstraintRef, Tuple{Float64, Float64}}(c1 : 6 x + 8 y ≥ 100.0 => (-4.0, 2.857142857142857), y ≤ 3.0 => (-1.75, Inf), c2 : 7 x + 12 y ≥ 120.0 => (-3.3333333333333335, 4.666666666666667), y ≥ 0.0 => (-Inf, 1.25), z ≤ 1.0 => (-Inf, Inf), x ≥ 0.0 => (-Inf, 15.0), c3 : x + y ≤ 20.0 => (-3.75, Inf)), Dict{VariableRef, Tuple{Float64, Float64}}(y => (-4.0, 0.5714285714285714), x => (-0.3333333333333333, 3.0), z => (-Inf, 1.0)))

It returns a SensitivityReport object, which maps:

• Every variable reference to a tuple (d_lo, d_hi)::Tuple{Float64,Float64}, explaining how much the objective coefficient of the corresponding variable can change by, such that the original basis remains optimal.
• Every constraint reference to a tuple (d_lo, d_hi)::Tuple{Float64,Float64}, explaining how much the right-hand side of the corresponding constraint can change by, such that the basis remains optimal.

Both tuples are relative, rather than absolute. So, given an objective coefficient of 1.0 and a tuple (-0.5, 0.5), the objective coefficient can range between 1.0 - 0.5 an 1.0 + 0.5.

For example:

report[x]
(-0.3333333333333333, 3.0)

indicates that the objective coefficient on x, that is, 12, can decrease by -0.333 or increase by 3.0 and the primal solution (15, 1.25) will remain optimal. In addition:

report[c1]
(-4.0, 2.857142857142857)

means that the right-hand side of the c1 constraint (100), can decrease by 4 units, or increase by 2.85 units, aand the primal solution (15, 1.25) will remain optimal.

## Variable sensitivity

By themselves, the tuples aren't informative. Let's put them in context by collating a range of other information about a variable:

function variable_report(xi)
return (
name = name(xi),
lower_bound = has_lower_bound(xi) ? lower_bound(xi) : -Inf,
value = value(xi),
upper_bound = has_upper_bound(xi) ? upper_bound(xi) : Inf,
reduced_cost = reduced_cost(xi),
obj_coefficient = coefficient(objective_function(model), xi),
allowed_decrease = report[xi],
allowed_increase = report[xi],
)
end
variable_report (generic function with 1 method)

Calling our function on x:

x_report = variable_report(x)
(name = "x", lower_bound = 0.0, value = 15.0, upper_bound = Inf, reduced_cost = 0.0, obj_coefficient = 12.0, allowed_decrease = -0.3333333333333333, allowed_increase = 3.0)

That's a bit hard to read, so let's call this on every variable in the model and put things into a DataFrame:

variable_df =
DataFrames.DataFrame(variable_report(xi) for xi in all_variables(model))

3 rows × 8 columns

namelower_boundvalueupper_boundreduced_costobj_coefficientallowed_decreaseallowed_increase
StringFloat64Float64Float64Float64Float64Float64Float64
1x0.015.0Inf0.012.0-0.3333333.0
2y0.01.253.00.020.0-4.00.571429
3z-Inf1.01.0-1.0-1.0-Inf1.0

Great! That looks just like the reports in Excel.

## Constraint sensitivity

We can do something similar with constraints:

function constraint_report(ci)
return (
name = name(ci),
value = value(ci),
rhs = normalized_rhs(ci),
slack = normalized_rhs(ci) - value(ci),
allowed_decrease = report[ci],
allowed_increase = report[ci],
)
end

c1_report = constraint_report(c1)
(name = "c1", value = 100.0, rhs = 100.0, slack = 0.0, shadow_price = -0.25, allowed_decrease = -4.0, allowed_increase = 2.857142857142857)

That's a bit hard to read, so let's call this on every variable in the model and put things into a DataFrame:

constraint_df = DataFrames.DataFrame(
constraint_report(ci) for (F, S) in list_of_constraint_types(model) for
ci in all_constraints(model, F, S) if F == AffExpr
)

3 rows × 7 columns

StringFloat64Float64Float64Float64Float64Float64
1c1100.0100.00.0-0.25-4.02.85714
2c2120.0120.00.0-1.5-3.333334.66667
3c316.2520.03.750.0-3.75Inf

## Analysis questions

Now we can use these dataframes to ask questions of the solution.

For example, we can find basic variables by looking for variables with a reduced cost of 0:

basic = filter(row -> iszero(row.reduced_cost), variable_df)

2 rows × 8 columns

namelower_boundvalueupper_boundreduced_costobj_coefficientallowed_decreaseallowed_increase
StringFloat64Float64Float64Float64Float64Float64Float64
1x0.015.0Inf0.012.0-0.3333333.0
2y0.01.253.00.020.0-4.00.571429

and non-basic variables by looking for non-zero reduced costs:

non_basic = filter(row -> !iszero(row.reduced_cost), variable_df)

1 rows × 8 columns

namelower_boundvalueupper_boundreduced_costobj_coefficientallowed_decreaseallowed_increase
StringFloat64Float64Float64Float64Float64Float64Float64
1z-Inf1.01.0-1.0-1.0-Inf1.0

we can also find constraints that are binding by looking for zero slacks:

binding = filter(row -> iszero(row.slack), constraint_df)

2 rows × 7 columns

StringFloat64Float64Float64Float64Float64Float64
1c1100.0100.00.0-0.25-4.02.85714
2c2120.0120.00.0-1.5-3.333334.66667

binding2 = filter(row -> !iszero(row.shadow_price), constraint_df)