LP 2: Modeling and Sensitivity Analysis

Packages Used: Functions from the following packages are used in this notebook for the first time:

• GLPK
• JuMP

1. Modeling

Clearly, as shown in LP 1, the Nelder-Mead and gradient-based methods are not effective procedures for solving problems with a linear objective function. Specialized procedures like the simplex method (1947) can be used to solve these problems effectively.

Ex 1: Coffee and Peppers Arbitrage Opportunity (cont. from Ex 5 in LP-1)

In LP-1, first, through displaying the feasible region and objective function, the optimal solution was determined to be to purchase 15 pounds of both coffee and peppers. Both constrained nonlinear optimization and linear programming were used to solve this problem. The Nelder-Mead method used for constrained nonlinear optimization had difficulty finding the known optimal solution for some starting points, while the (revised) simplex method used for linear programming easily found the known optimal. As the size of problems gets larger (e.g., more items to consider for purchasing), this difference gets worse and is the reason that specialized techniques like the simplex method have been developed for constrained linear optimization problems.

Question: You are planning an international trip and have the opportunity to purchase some estate-grown coffee and high-quality dried peppers while there that you can sell upon your return to a friend who supplies several local restaurants. You can sell the coffee and peppers to your friend for \$6 and \\$15 per pound, respectively, and can purchase them at your destination for \$2 and \\$6 per pound, respectively. You have a budget of \$120 to spend on these purchases, and you have estimated that you will be able to pack up to 30 pounds, in total, of both items in your knapsack since there is a maximum weight limit of 50 pounds on the one allowed checked bag and your clothes and the knapsack together weigh 20 pounds. How many pounds of coffee and peppers should you purchase?

Solution Approach:

Objective will be to try to maximize the total profit from a mix of both items. Unit profits for coffee and peppers are \$4 and \\$9 per pound, respectively, and so total profit will be the linear function

$ \begin{eqnarray*} \quad \mbox{Total profit:} \quad \mathit{\Pi} &=& \pi_1 x_1 + \pi_2 x_2 \\ &=& \left( p_1 - c_1 \right) x_1 + \left( p_2 - c_2 \right) x_2 \\ &=& ( 6 - 2) x_1 + ( 15 - 6 ) x_2 \\ &=& 4x_1 + 9x_2 \end{eqnarray*}$,

where the decision variables $x_1$ and $x_2$ are pounds of coffee and peppers purchased, $p_1$ and $p_2$ the sale price, and $c_1$ and $c_2$ the purchase price, respectiveley. The budget and weight constaints are the following:

$ \begin{eqnarray*} \quad\mbox{Budget constraint:} \quad &c_1 x_1 + c_2 x_2 &\le& 120 = 2x_1 + 6x_2 \le 120 \\ \mbox{Weight constraint:} \quad &x_1 + x_2 &\le& 30 \end{eqnarray*}$

Also, the decision variable are restricted from being negative with the lower-bound constraints $x_1 \ge 0$ and $x_2 \ge 0$.

Create Model using JuMP and Solve it Using GLPK

The package JuMP provides a convenient algebraic modeling language for creating mathematical programming models, and the package GLPK is an open-source solver that implements the simplex method:

# Construct model
using JuMP, GLPK
model = Model(GLPK.Optimizer)
@variable(model, 0 <= x₁ )
@variable(model, 0 <= x₂ )
@objective(model, Max, 4x₁ + 9x₂ )
@constraint(model, 2x₁ + 6x₂ <= 120 )
@constraint(model, x₁ + x₂ <= 30 )
print(model)

$$ \begin{aligned} \max\quad & 4 x₁ + 9 x₂\\ \text{Subject to} \quad & 2 x₁ + 6 x₂ \leq 120.0\\ & x₁ + x₂ \leq 30.0\\ & x₁ \geq 0.0\\ & x₂ \geq 0.0\\ \end{aligned} $$

# Solve model and display results
optimize!(model)
Πᵒ = objective_value(model)
x₁ᵒ = value(x₁)
x₂ᵒ = value(x₂)
Πᵒ, x₁ᵒ, x₂ᵒ

(195.0, 14.999999999999998, 15.0)

Display Model and Solution

using Plots
plot(Shape([0, 30, 15, 0],[0, 0, 15, 20]) , opacity=.25, label="Feasible region",
    xlims=(0,65), ylims=(0,35), xlabel="Coffee (lb)", ylabel="Peppers (lb)", 
    aspect_ratio=:equal)
xrng = 0:.1:65
fline(x₁,a₁,a₂,b) = abs(a₂) > sqrt(eps()) ? (b - a₁*x₁)/a₂ : x₁   # Line function
plot!(xrng, [fline(i,2,6,120) for i in xrng], linewidth=2, label="Budget constraint")
plot!(xrng, [fline(i,1,1,30) for i in xrng], linewidth=2, label="Weight constraint")
plot!(xrng, [fline(i,4,9,Πᵒ) for i in xrng], linewidth=2, label="Objective function")
scatter!([x₁ᵒ], [x₂ᵒ], label="Optimal solution")

Observations:

1. Solution procedure:

For this problem, the linear programming approach is preferred to penalty-constrained optimization; this is due to its linear objective function.

The penalty-constrained optimization (see the result, below, from Ex 5 in LP 1) did not reach the optimal solution of [15, 15] when starting from [0, 0], and the number of function evaluations increased as the starting point was placed closer to the optimal solution. Clearly, the general-purpose Nelder-Mead nonlinear optimization algorithm has difficulties when an objective function is a flat plane, as is the case for a linear objective function.

The linear programming solution procedure (the simplex method) only pivots from one corner point of the feasible region to another, eventually reaching the optimal solution, which, because the objective function is linear, is always at one of the corner points. For example, the simplex method needs to make only two pivots to reach [15, 15] starting at [0, 0], while penalty-constrained optimization required 99 iterations starting at [1, 1] (see Ex 5 in LP 1). (There are specialized interior point linear programming methods that do search inside the feasible region; these methods are more efficient than the simplex method but not as useful when used inside of mixed-integer programming problem solvers.)

2. Continuous variables:

_The quantities $x_i$ are assumed to be continuous, real-valued decision variables._

The size of individual coffee beans and peppers is small enough relative to the size of the knapsack that representing the quantities of both as continuous variables seems reasonable. If only a few units of each item are being considered, then just rounding a continuous variable may lead to a nonoptimal result; instead, integer-valued variables could be used to represent the quantities of each item, which would then require solving a mixed-integer linear program (MILP), which is usually much more difficult to solve.

For example, if the coffee and peppers only could be purchased in 4-pound bags, then rounding the continuous result to get the number of 4-bags would give a result that exceeds the 30-pound weight limit:

round(x₁ᵒ/4), round(x₂ᵒ/4),
4*round(x₁ᵒ/4) + 4*round(x₂ᵒ/4),
4*round(x₁ᵒ/4) + 4*round(x₂ᵒ/4) <= 30

(4.0, 4.0, 32.0, false)

In order to correctly account for the 4-pound bags, the model can include integer restrictions on the $x_i$, which makes it a MILP, and where the quantities $x_i$ now represent the numbers of bags instead of pounds purchased.

3. Adding constraints:

Adding a constraint to model never improves the solution.

Adding a constraint either reduces or does not change a solution. The solution is reduced if the additional constraint is binding; that is, if the constraint cuts off the previous solution from the feasible region, the solution is unchanged if the additional constraint is nonbinding.

Extension: Additional Constraints

Examples of additional constraints that can be added to the model include the following [Brown and Dell, 2017]:

1. Bounded proportion

For each pound of coffee, there must be at least five pounds of peppers purchased.

@constraint(model, con1, x₁ <= x₂/5 )   # Add constraint to model and then solve
print(model)
optimize!(model)
Πᵒ = objective_value(model)
x₁ᵒ = value(x₁)
x₂ᵒ = value(x₂)
Πᵒ, x₁ᵒ, x₂ᵒ, x₂ᵒ/x₁ᵒ

$$ \begin{aligned} \max\quad & 4 x₁ + 9 x₂\\ \text{Subject to} \quad & 2 x₁ + 6 x₂ \leq 120.0\\ & x₁ + x₂ \leq 30.0\\ & x₁ - 0.2 x₂ \leq 0.0\\ & x₁ \geq 0.0\\ & x₂ \geq 0.0\\ \end{aligned} $$

(183.75, 3.7499999999999982, 18.75, 5.000000000000003)

2. Fixed proportion

There must be exactly three pounds of coffee purchased for every two pounds of peppers purchased.

delete(model, con1)                     # Remove previous constraint
@constraint(model, con2, x₁/3 == x₂/2 )
print(model)
optimize!(model)
Πᵒ = objective_value(model)
x₁ᵒ = value(x₁)
x₂ᵒ = value(x₂)
Πᵒ, x₁ᵒ, x₂ᵒ, x₂ᵒ/x₁ᵒ

$$ \begin{aligned} \max\quad & 4 x₁ + 9 x₂\\ \text{Subject to} \quad & 0.3333333333333333 x₁ - 0.5 x₂ = 0.0\\ & 2 x₁ + 6 x₂ \leq 120.0\\ & x₁ + x₂ \leq 30.0\\ & x₁ \geq 0.0\\ & x₂ \geq 0.0\\ \end{aligned} $$

(180.0, 18.0, 12.0, 0.6666666666666666)

3. Rates

Only three hours are available for shopping, and it takes an hour to purchase six pounds of coffee and an hour to purchase four pounds of peppers.

delete(model, con2)
@constraint(model, con3, x₁/6 + x₂/4 <= 3 ) 
print(model)
optimize!(model)
Πᵒ = objective_value(model)
x₁ᵒ = value(x₁)
x₂ᵒ = value(x₂)
Πᵒ, x₁ᵒ, x₂ᵒ

$$ \begin{aligned} \max\quad & 4 x₁ + 9 x₂\\ \text{Subject to} \quad & 2 x₁ + 6 x₂ \leq 120.0\\ & x₁ + x₂ \leq 30.0\\ & 0.16666666666666666 x₁ + 0.25 x₂ \leq 3.0\\ & x₁ \geq 0.0\\ & x₂ \geq 0.0\\ \end{aligned} $$

(108.0, 0.0, 12.0)

4. Max percent of total

Coffee cannot exceed one-third of the total pounds of coffee and peppers purchased.

delete(model, con3)
@constraint(model, con4, x₁ <= (x₁ + x₂)/3 )
print(model)
optimize!(model)
Πᵒ = objective_value(model)
x₁ᵒ = value(x₁)
x₂ᵒ = value(x₂)
Πᵒ, x₁ᵒ, x₂ᵒ, x₁ᵒ/(x₁ᵒ + x₂ᵒ)

$$ \begin{aligned} \max\quad & 4 x₁ + 9 x₂\\ \text{Subject to} \quad & 2 x₁ + 6 x₂ \leq 120.0\\ & x₁ + x₂ \leq 30.0\\ & 0.6666666666666667 x₁ - 0.3333333333333333 x₂ \leq 0.0\\ & x₁ \geq 0.0\\ & x₂ \geq 0.0\\ \end{aligned} $$

(188.57142857142856, 8.57142857142857, 17.142857142857142, 0.3333333333333333)

5. Bound on the mix

The average carbon footprint of imported items cannot exceed twenty grams per pound. The footprint of coffee and peppers are eighty and five grams per pound, respectively.

delete(model, con4)
@constraint(model, con5, 80x₁ + 5x₂ <= 20(x₁ + x₂) )
print(model)
optimize!(model)
Πᵒ = objective_value(model)
x₁ᵒ = value(x₁)
x₂ᵒ = value(x₂)
Πᵒ, x₁ᵒ, x₂ᵒ, (80x₁ᵒ + 5x₂ᵒ)/(x₁ᵒ + x₂ᵒ)

$$ \begin{aligned} \max\quad & 4 x₁ + 9 x₂\\ \text{Subject to} \quad & 2 x₁ + 6 x₂ \leq 120.0\\ & x₁ + x₂ \leq 30.0\\ & 60 x₁ - 15 x₂ \leq 0.0\\ & x₁ \geq 0.0\\ & x₂ \geq 0.0\\ \end{aligned} $$

(184.6153846153846, 4.615384615384615, 18.461538461538463, 19.999999999999996)

2. Sensitivity Analysis

Sensitivity analysis for an LP involving looking at how changing any of the coefficients in the objective or a constraint would change the solution. It is useful because it gives a decision-maker guidance with respect to how much any of the data used to find a particular solution to a problem can change before there is a change in the optimal solution. This is useful because often, there is a lot of uncertainty regarding many of the coefficients used to formulate an LP (they may have been guesstimated), and sensitivity analysis for these coefficients can indicate whether additional effort is needed to provide a better estimate. Although there are analytical methods for sensitivity analysis specific to an LP, in the following, a more general numerical approach is used that can be used for other optimization models that usually do not have analytical-based methods available for sensitivity analysis.

Ex 2: Pepper Sales Price versus Purchase Quantity

Continuing with the Coffee and Peppers Arbitrage problem from Ex 1, instead of a firm \$15 per pound, assume that there has not been a firm indication of how much the peppers can be sold for. You can always get \\$5 per pound selling the peppers to a friend who runs a local market, but you have heard rumors that some restaurant owners will pay up \$45 per pound if there is a shortage of the peppers since they are not commonly available. Your guesstimate of the price you will receive when you get back from your trip is the following:

@show p₂ = sqrt(5*45);

p₂ = sqrt(5 * 45) = 15.0

Question: How much could the price received for peppers change without changing the plan to purchase 15 pounds of peppers?

Solution approach: Will first display the optimal quantity of peppers to purchase for prices between your LB and UB, and then determine the range of possible prices that do not change the current plan.

Note: It is important not to evaluate too many prices since each will require a full LP solution. First, look coarsely, then narrow the bounds to focus if more resolution is needed.

function price2quantity(p₂)
    model = Model(GLPK.Optimizer)
    @variable(model, 0 <= x₁ )
    @variable(model, 0 <= x₂ )
    @objective(model, Max, 4x₁ + (p₂ - 6)x₂ )
    @constraint(model, 2x₁ + 6x₂ <= 120 )
    @constraint(model, x₁ + x₂ <= 30 )
    optimize!(model)
    return value(x₂)
end
p = [5:1:45;]                           # Will require 41 full LP solutions
q = [price2quantity(pi) for pi in p]
plot(p, q, legend=false, xlabel="Purchase Price (\$)", ylabel="Purchase Quantity (lb)")

using Base
d = diff(q)

40-element Vector{Float64}:
  0.0
  0.0
  0.0
  0.0
  0.0
 15.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  5.0
  ⋮
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0

idx = findall(d .> 0)

2-element Vector{Int64}:
  6
 13

idx = [1; idx[:]; length(p)]

4-element Vector{Int64}:
  1
  6
 13
 41

idx = [1; findall( abs.(diff(q)) .> sqrt(eps()) )[:]; length(p)]

4-element Vector{Int64}:
  1
  6
 13
 41

p[idx]

4-element Vector{Int64}:
  5
 10
 17
 45

plot(p, q, legend=false, xlabel="Purchase Price (\$)", ylabel="Purchase Quantity (lb)",
    xticks=p[idx], yticks=q[idx])

Answer: The price received for peppers could decrease to \$10 or increase to \\$17 without changing the plan to purchase 15 pounds of peppers (note: although the purchase quantity does not change, the total profit does change).

Ex 3: Budget versus Total Profit

Question: Continuing with Ex 1, what is the minimum budget for purchasing coffee and peppers that will maximize total profit?

Note: There is limit on who much the budget can be increased because of the other weight constraint of 30 pounds.

function budget2profit(bᴮ)
    model = Model(GLPK.Optimizer)
    @variable(model, 0 <= x₁ )
    @variable(model, 0 <= x₂ )
    @objective(model, Max, 4x₁ + 9x₂ )
    @constraint(model, 2x₁ + 6x₂ <= bᴮ )
    @constraint(model, x₁ + x₂ <= 30 )
    optimize!(model)
    return objective_value(model)
end
bᴮ = [0:1:450;]
Π = [budget2profit(bᴮi) for bᴮi in bᴮ]
plot(bᴮ, Π, legend=false, xlabel="Budget (\$)", ylabel="Total Profit (\$)")

idx = argmax(Π)
idx, bᴮ[idx], Π[idx]

(181, 180, 270.0)

str = string(" \nMinimum Budget of \$", string(bᴮ[idx]), 
    " Maximizes Total Profit at \$", string(Π[idx]))
scatter!([bᴮ[idx]], [Π[idx]], title=str)

3. Another Example

Ex 4: Diet Problem

Ozark Farms uses at least 800 pounds of feed daily that is a mixture of corn and soybean meal with the following compositions:

Feedstuff	Protein (lb/lb)	Fiber (lb/lb)	Cost (\$/lb)
Corn	.09	.02	.30
Soybean meal	.60	.06	.90

Dietary requirements of the feed are that it contain at least 30\% protein and at most 5\% fiber. Ozark Farms would like to determine the minimum-cost daily feed mix. [Adapted from H.A. Tata, 2002, pp 18-19]

using JuMP, GLPK
model = Model(GLPK.Optimizer)
@variable(model, 0 <= x₁ )                            # Corn in daily mix (lb)
@variable(model, 0 <= x₂ )                            # Soybean meal in daily mix (lb)
@objective(model, Min, .3x₁ + .9x₂ )
@constraint(model, x₁ + x₂ >= 800 )                   # Minimum daily feed amount
@constraint(model, .09x₁ + .6x₂ >= .3(x₁ + x₂) )      # At least 30% protein
@constraint(model, .02x₁ + .06x₂ <= .05(x₁ + x₂) )    # At most 5% fiber
print(model)
set_optimizer_attribute(model, "msg_lev", GLPK.GLP_MSG_ALL)
optimize!(model)
Cᵒ = objective_value(model)
x₁ᵒ = value(x₁)
x₂ᵒ = value(x₂)
Cᵒ, x₁ᵒ, x₂ᵒ

$$ \begin{aligned} \min\quad & 0.3 x₁ + 0.9 x₂\\ \text{Subject to} \quad & x₁ + x₂ \geq 800.0\\ & -0.21 x₁ + 0.3 x₂ \geq 0.0\\ & -0.030000000000000002 x₁ + 0.009999999999999995 x₂ \leq 0.0\\ & x₁ \geq 0.0\\ & x₂ \geq 0.0\\ \end{aligned} $$

GLPK Simplex Optimizer 5.0
3 rows, 2 columns, 6 non-zeros
      0: obj =  0.000000000e+000 inf =  8.000e+002 (1)
      2: obj =  4.376470588e+002 inf =  0.000e+000 (0)
OPTIMAL LP SOLUTION FOUND

(437.6470588235294, 470.5882352941176, 329.4117647058823)