LP 3: Larger-Scale LP Models

Packages Used: No new packages are used in this notebook.

It is convenient when creating larger-scale LP models to use arrays for specifying individual scalar decision variables and coefficients instead of individual scalar values. Set notation can be used to indicate iteration over a set of array indices as long as the order in which the iteration occurs does not matter; when order is important, sequences should be used.

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

Everything is the same except that n-element arrays will be used instead of individual scalars for the sale price, purchase price, and purchase quantity, where n is the number of items purchased and $ N = \{ 1, \dots, n \}$$ = \{ 1, 2 \} $ is the set of array indices:

$ \begin{eqnarray*} \quad \mbox{Maximize} \quad \sum_{i \in N} \left( p_i - c_i \right) x_i \\ \quad \mbox{subject to} \quad \sum_{i \in N} c_i x_i &\le& 120 &\quad\mbox{:Budget constraint}\\ \sum_{i \in N} x_i &\le& 30 &\quad\mbox{:Weight constraint}\\ x_i &\ge& 0 \mbox{,} \end{eqnarray*} $

which can be then be formulated in JuMP as follows:

using JuMP, GLPK
p = [6, 15]
c = [2, 6]
N = 1:length(p)
model = Model(GLPK.Optimizer)
@variable(model, 0 <= x[i = N] )                             # UB = ∞ for each xᵢ
@objective(model, Max, sum((p[i] - c[i])*x[i] for i ∈ N) )   # \in<TAB> ⟹ ∈
@constraint(model, sum(c[i]*x[i] for i ∈ N) <= 120 )         # Budget constraint
@constraint(model, sum(x[i] for i ∈ N) <= 30 )               # Weight constraint
print(model)
optimize!(model)
Πᵒ = objective_value(model)
xᵒ = [value(x[i]) for i ∈ N]
Πᵒ, xᵒ

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

(195.0, [14.999999999999998, 15.0])

If both items have different maximum order quantity restrictions, modeled as UBs of 40 and 10:

using JuMP, GLPK
u = [40, 10]                                                 
p = [6, 15]
c = [2, 6]
N = 1:length(p)
model = Model(GLPK.Optimizer)
@variable(model, 0 <= x[i = N] <= u[i] )                     # Different UB for each xᵢ
@objective(model, Max, sum((p[i] - c[i])*x[i] for i ∈ N) )
@constraint(model, sum(c[i]*x[i] for i ∈ N) <= 120 )
@constraint(model, sum(x[i] for i ∈ N) <= 30 )
print(model)
optimize!(model)
Πᵒ = objective_value(model)
xᵒ = [value(x[i]) for i ∈ N]
Πᵒ, xᵒ

$$ \begin{aligned} \max\quad & 4 x_{1} + 9 x_{2}\\ \text{Subject to} \quad & 2 x_{1} + 6 x_{2} \leq 120.0\\ & x_{1} + x_{2} \leq 30.0\\ & x_{1} \geq 0.0\\ & x_{2} \geq 0.0\\ & x_{1} \leq 40.0\\ & x_{2} \leq 10.0\\ \end{aligned} $$

(170.0, [20.0, 10.0])

Ex 2: Multi-Item Arbitrage Opportunity (cont. from Ex 1)

Continuing with Ex 1, if there were three other items avaiable in addition to coffee and peppers:

using DataFrames
df = DataFrame(Item = ["Coffee", "Peppers", "Ginger", "Tea", "Candy"], 
    Max_Order = [40, 10, 2, 5, 15], Sale_Price = [6, 15, 8, 25, 4], 
    Purchase_Price = [2, 6, 1, 5, 2])

5 rows × 4 columns

	Item	Max_Order	Sale_Price	Purchase_Price
	String	Int64	Int64	Int64
1	Coffee	40	6	2
2	Peppers	10	15	6
3	Ginger	2	8	1
4	Tea	5	25	5
5	Candy	15	4	2

using JuMP, GLPK
u,p,c = df.Max_Order, df.Sale_Price, df.Purchase_Price
N = 1:length(p)
model = Model(GLPK.Optimizer)
@variable(model, 0 <= x[i = N] <= u[i] )
@objective(model, Max, sum((p[i] - c[i])*x[i] for i ∈ N) )
@constraint(model, sum(c[i]*x[i] for i ∈ N) <= 120 )
@constraint(model, sum(x[i] for i ∈ N) <= 30 )
print(model)
optimize!(model)
df.Purchase_Quantity = value.(x)                       # Add output to DataFrame
println(" \nTotal Profit: ", objective_value(model))   # Format output
df

$$ \begin{aligned} \max\quad & 4 x_{1} + 9 x_{2} + 7 x_{3} + 20 x_{4} + 2 x_{5}\\ \text{Subject to} \quad & 2 x_{1} + 6 x_{2} + x_{3} + 5 x_{4} + 2 x_{5} \leq 120.0\\ & x_{1} + x_{2} + x_{3} + x_{4} + x_{5} \leq 30.0\\ & x_{1} \geq 0.0\\ & x_{2} \geq 0.0\\ & x_{3} \geq 0.0\\ & x_{4} \geq 0.0\\ & x_{5} \geq 0.0\\ & x_{1} \leq 40.0\\ & x_{2} \leq 10.0\\ & x_{3} \leq 2.0\\ & x_{4} \leq 5.0\\ & x_{5} \leq 15.0\\ \end{aligned} $$

 
Total Profit: 256.0

5 rows × 5 columns

	Item	Max_Order	Sale_Price	Purchase_Price	Purchase_Quantity
	String	Int64	Int64	Int64	Float64
1	Coffee	40	6	2	13.0
2	Peppers	10	15	6	10.0
3	Ginger	2	8	1	2.0
4	Tea	5	25	5	5.0
5	Candy	15	4	2	0.0

Ex 3: Production Mix

Three different chemical products are produced using four different raw materials as input. Up to 100, 40, and 80 liters of each product can be sold for \$10, 20, and 30 per L, respectively. Producing 10 liters of the first product requires 2, 3, and 5 L of inputs 1-3, respectively; 10 liters of the second requires 7 and 3 L of inputs three and four; 5 liters of the third product requires 2 L of inputs one and two and 1 L of input four. Up to 40, 20, 30, and 80 L of inputs 1-4, respectively, are available at the cost of \\$4, 3, 5, and 1 per L. Determine how many liters of each product should be produced.

Solution: First step is to make the fixed proportions commensurate by determining the quantity of each input needed to produce 1 L of each product:

A = [2 3 5 0; 0 0 7 3; 2 2 0 1] ./ [10, 10, 5]

3×4 Matrix{Float64}:
 0.2  0.3  0.5  0.0
 0.0  0.0  0.7  0.3
 0.4  0.4  0.0  0.2

Next, the following LP is formulated:

$ \begin{eqnarray*} \quad \mbox{Maximize} \quad \sum_{i \in N} p_i x_i &-& \sum_{j \in M} c_j \sum_{i \in N} a_{ij} x_i \\ \quad \mbox{subject to} \quad \sum_{i \in N} a_{ij} x_i &\le& b_j \mbox{,} \quad j \in M \\ 0 \le x_i &\le& u_i \mbox{,} \end{eqnarray*} $

where the decision variables $x_i, p_i, u_i, i \in N$, represent the production (L), sale price (\$/L), and upper bound (L) of each product in the set $N$ of products; $c_j, bj, j \in M$, represent the cost (\\$/L) and upper bound (L) of each input in the set $M$ of inputs; and $a{ij}$ represents the liters of input j required to produce one liter of product i.

p = [10,20,30]      # Product sale price ($/L)
c = [4,3,5,1]       # Input cost ($/L)
u = [100,40,80]     # Product upper bound (L)
b = [40,20,30,80]   # Input upper bound (L)

using JuMP, GLPK, DataFrames
N = 1:length(p)
M = 1:length(c)
model = Model(GLPK.Optimizer)
@variable(model, 0. <= x[i = N] <= u[i])
@objective(model, Max, sum(p[i]*x[i] for i ∈ N) - sum(c[j]*A[i,j]*x[i] for i ∈ N, j ∈ M) )
@constraint(model, con[j = M], sum(A[i,j]*x[i] for i ∈ N) <= b[j] )   # Label constriants
print(model)
optimize!(model)
Πᵒ = objective_value(model)
xᵒ = [value(x[i]) for i ∈ N]
Πᵒ, xᵒ

$$ \begin{aligned} \max\quad & 5.799999999999999 x_{1} + 16.2 x_{2} + 27 x_{3}\\ \text{Subject to} \quad & 0.2 x_{1} + 0.4 x_{3} \leq 40.0\\ & 0.3 x_{1} + 0.4 x_{3} \leq 20.0\\ & 0.5 x_{1} + 0.7 x_{2} \leq 30.0\\ & 0.3 x_{2} + 0.2 x_{3} \leq 80.0\\ & x_{1} \geq 0.0\\ & x_{2} \geq 0.0\\ & x_{3} \geq 0.0\\ & x_{1} \leq 100.0\\ & x_{2} \leq 40.0\\ & x_{3} \leq 80.0\\ \end{aligned} $$

(1998.0, [0.0, 40.0, 50.0])

Note: In the above model, @constraint(model, con[j = M], is used to label each constraint as con1 to con4 so that they can be referenced; for example (see below), each constraint is referenced to determine its "shadow price."

Sensitivity Analysis

In order to provide more information for sensitivity analysis of the solution, the following additional data can be reported:

yᵒ = zeros(length(M))
for j ∈ M
    for i ∈ N
       yᵒ[j] += A[i,j]*xᵒ[i]
    end
end
yᵒ

4-element Vector{Float64}:
 20.0
 20.0
 28.0
 22.0

or, equivalently:

yᵒ = [sum([sum(A[i,j]*xᵒ[i]) for i ∈ N]) for j ∈ M]         # Input requirements (L)

4-element Vector{Float64}:
 20.0
 20.0
 28.0
 22.0

println(" \nTotal Profit: ", objective_value(model), "\n")
df = DataFrame(Product = N, Price = p, Quantity = xᵒ, UB = u)
df2 = DataFrame(Input = M, Cost = c, Quantity = yᵒ, UB = b)
show(df, show_row_number=false), println("\n"), show(df2, show_row_number=false);

 
Total Profit: 1998.0

3×4 DataFrame
 Product  Price  Quantity  UB    
 Int64    Int64  Float64   Int64 
─────────────────────────────────
       1     10       0.0    100
       2     20      40.0     40
       3     30      50.0     80

4×4 DataFrame
 Input  Cost   Quantity  UB    
 Int64  Int64  Float64   Int64 
───────────────────────────────
     1      4      20.0     40
     2      3      20.0     20
     3      5      28.0     30
     4      1      22.0     80

Since the upper bound of 20 for input 2 is binding, increasing the bound should increase profits; for example, increasing the bound to 21 increases profit by \$67.50, from \\$1998 to \$2065.50. This is also reflected in the shadow price of each constraint:

df2.Shadow_Price = [shadow_price(con[j]) for j ∈ M]
show(df2, show_row_number=false)

4×5 DataFrame
 Input  Cost   Quantity  UB     Shadow_Price 
 Int64  Int64  Float64   Int64  Float64      
─────────────────────────────────────────────
     1      4      20.0     40           0.0
     2      3      20.0     20          67.5
     3      5      28.0     30           0.0
     4      1      22.0     80           0.0

Reading Matrix Data from a File

The following can be used to read the matrix A from the CSV file LP-3-Data.csv:

0.2,0.3,0.5,0
0,0,0.7,0.3
0.4,0.4,0,0.2

using DataFrames, CSV
A = Matrix(DataFrame(CSV.File("LP-3-Data.csv", header=false)))

3×4 Matrix{Float64}:
 0.2  0.3  0.5  0.0
 0.0  0.0  0.7  0.3
 0.4  0.4  0.0  0.2

Ex 4: Employee Staffing

Question: Wally World adventure park is deciding how many full-time employees to hire for the upcoming season. Each full-time employee should work for five consecutive days each week, followed by two days off. They have estimated that their staffing requirements for Sunday through Saturday each week are 150, 160, 110, 175, 130, 150, and 190 full-time employees, respectively. Additional part-time staff will be hired for peak periods during the season. Determine the minimum number of full-time employees needed.

Solution approach: The decision variable $x_i$ represents the number of employees that start their five consecutive days of work on day $i$, for $i = 1$ corresponding to Sunday. Since consecutive days will be defined, sequences of indices $i$ for $x_i$ will be used instead of set notation. The most difficult part of formulation a model is getting the sequences of $x_i$'s to wrap around day 7 (Saturday) to continue at day 1 (Sunday). The basic model will be

$ \begin{eqnarray} \quad \mbox{Minimize} \quad \sum_{i = 1}^{7} x_i \\ \quad \mbox{subject to} \quad \sum_{\mbox{j-start ?}}^{\mbox{j-end ?}} x_{\mbox{[x-index ?]}} &\ge& b_i \mbox{,} &\quad i = 1, \dots, 7 \\ x_i &\ge& 0 \mbox{.} \end{eqnarray} $

Since employees should work five consectitive days, the only employees not able to work day $i$ would be those that start on days $i + 1$ and $i + 2$, the summation starting limit ("j-start") should be $i + 3$, or 3 days in the future, and the ending limit ("j-end") $i + 3 + (5 - 1)$ = $i + 7$ days in the future. Since $i > 5 \implies $$j > 7$, the index for $x$ ("x-index") should be $j - 7 \left\lfloor\frac{j - 1}{7} \right\rfloor$ so that, for example, $j = 8 \implies $$ 8 - 7 \left\lfloor\frac{7}{7} \right\rfloor = $$1$. It is a good idea to test the indices, where k is the x-index:

for i = 1:7
    for j = i+3:i+7
        k = Int(j - 7*floor((j - 1)/7))
        @show i,j,k
    end
end

(i, j, k) = (1, 4, 4)
(i, j, k) = (1, 5, 5)
(i, j, k) = (1, 6, 6)
(i, j, k) = (1, 7, 7)
(i, j, k) = (1, 8, 1)
(i, j, k) = (2, 5, 5)
(i, j, k) = (2, 6, 6)
(i, j, k) = (2, 7, 7)
(i, j, k) = (2, 8, 1)
(i, j, k) = (2, 9, 2)
(i, j, k) = (3, 6, 6)
(i, j, k) = (3, 7, 7)
(i, j, k) = (3, 8, 1)
(i, j, k) = (3, 9, 2)
(i, j, k) = (3, 10, 3)
(i, j, k) = (4, 7, 7)
(i, j, k) = (4, 8, 1)
(i, j, k) = (4, 9, 2)
(i, j, k) = (4, 10, 3)
(i, j, k) = (4, 11, 4)
(i, j, k) = (5, 8, 1)
(i, j, k) = (5, 9, 2)
(i, j, k) = (5, 10, 3)
(i, j, k) = (5, 11, 4)
(i, j, k) = (5, 12, 5)
(i, j, k) = (6, 9, 2)
(i, j, k) = (6, 10, 3)
(i, j, k) = (6, 11, 4)
(i, j, k) = (6, 12, 5)
(i, j, k) = (6, 13, 6)
(i, j, k) = (7, 10, 3)
(i, j, k) = (7, 11, 4)
(i, j, k) = (7, 12, 5)
(i, j, k) = (7, 13, 6)
(i, j, k) = (7, 14, 7)

The complete model formulation is then:

$ \begin{eqnarray} \quad \mbox{Minimize} \quad \sum_{i = 1}^{7} x_i \\ \quad \mbox{subject to} \quad \sum_{j = i+3}^{i+7} x_{j\,-\,7 \left\lfloor\frac{j - 1}{7} \right\rfloor} &\ge& b_i \mbox{,} &\quad i = 1, \dots, 7 \\ x_i &\ge& 0 \mbox{,} \end{eqnarray} $

which can be then be formulated in JuMP as follows:

using JuMP, GLPK
b = [150, 160, 110, 175, 130, 150, 190]
model = Model(GLPK.Optimizer)
@variable(model, 0 <= x[1:7] )
@objective(model, Min, sum(x[i] for i = 1:7) )
for i = 1:7      # Can use FOR-loop since constraints are not being labeled
    @constraint(model, sum(x[Int(j - 7*floor((j - 1)/7))] for j = i+3:i+7) >= b[i] )
end
print(model)
optimize!(model)
Nᵒ = objective_value(model)
xᵒ = value.(x)
Nᵒ, xᵒ

$$ \begin{aligned} \min\quad & x_{1} + x_{2} + x_{3} + x_{4} + x_{5} + x_{6} + x_{7}\\ \text{Subject to} \quad & x_{1} + x_{4} + x_{5} + x_{6} + x_{7} \geq 150.0\\ & x_{1} + x_{2} + x_{5} + x_{6} + x_{7} \geq 160.0\\ & x_{1} + x_{2} + x_{3} + x_{6} + x_{7} \geq 110.0\\ & x_{1} + x_{2} + x_{3} + x_{4} + x_{7} \geq 175.0\\ & x_{1} + x_{2} + x_{3} + x_{4} + x_{5} \geq 130.0\\ & x_{2} + x_{3} + x_{4} + x_{5} + x_{6} \geq 150.0\\ & x_{3} + x_{4} + x_{5} + x_{6} + x_{7} \geq 190.0\\ & x_{1} \geq 0.0\\ & x_{2} \geq 0.0\\ & x_{3} \geq 0.0\\ & x_{4} \geq 0.0\\ & x_{5} \geq 0.0\\ & x_{6} \geq 0.0\\ & x_{7} \geq 0.0\\ \end{aligned} $$

(225.0, [0.0, 35.0, 0.0, 65.0, 50.0, 0.0, 75.0])

In general, the solution will not be integer-valued. As a rule of thumb, rounding the (continuous) decision variables in an LP solution when their values are greater than 20 will not result in too much error, but for values less than 20, the error can be significant, and a MILP model should be used with decision variables restricted to being an integer. In this example, taking the ceiling of any fractional values will ensure that the solution remains feasible, but at the cost of increasing the number of employees required; instead, the MILP model will make the correct internal decisions with respect to rounding that will result in the least possible increase in the number of employees.

Ex 5: LP to Minimize the Sum of Absolute Deviations

This is the same problem we solved using univariate nonlinear optimization in NL Opt 1.

using Plots
p = [1, 2, 4, 6, 9]
n = length(p)
LB = minimum(p) - 1
UB = maximum(p) + 1
dist1(x,p) = abs.(x .- p)
sumdist1(x) = sum(dist1(x,p))
plot(sumdist1, LB, UB, legend=false)
scatter!(p, [sumdist1(x) for x in p], xticks=p)

using Statistics
median(p)

4.0

using Optim
xᵒ = optimize(sumdist1, LB, UB).minimizer
TDᵒ = sumdist1(xᵒ)
TDᵒ, xᵒ

(12.000000027215956, 3.999999972784044)

Instead of as a MP, this unconstrained nonlinear optimization problem is represented mathematically as follows:

$ \quad\begin{eqnarray} x^* &=& \mathrm{arg}\underset{x}{\operatorname{min}} \left\{ \sum_{i=1}^{n} \lvert p_i - x \rvert : L\!B \leq x \leq U\!B \right\} \\ T\!D^* &=& \sum_{i=1}^{n} \lvert p_i - x^* \rvert \end{eqnarray}$

It cannot be directly formulated as an LP because of the absolute values used to determine the distances. This can be formulated as an LP using a change of variables.

Absolute Value: Use a change of variables to represent absolute values in an MP.

The variables $r_i \ge 0$ and $s_i \ge 0$ are introduced to represent the cases where the unrestricted $x$ is to the right or left of $p_i$, respectively:

$ \quad\left.\begin{matrix} x \ge p_i: r_i = x - p_i \mbox{ and } s_i = 0 \\ x < p_i: s_i = p_i - x \mbox{ and } r_i = 0 \end{matrix} \right\} \implies r_i + s_i = \lvert x - p_i \rvert,\, x - r_i + s_i = p_i, \mbox{ and } r_i, s_i \ge 0 $,

where each $r_i + s_i$ is a term in the objective function, each $x - r_i + s_i = p_i$ a constraint, and $r_i, s_i \ge 0$ are lower-bounds on the new variables.

using JuMP, GLPK
p = [1, 2, 4, 6, 9]
N = 1:length(p)
model = Model(GLPK.Optimizer)
@variable(model, x )
@variable(model, 0 <= r[i=N] )
@variable(model, 0 <= s[i=N] )
@objective(model, Min, sum(r[i] + s[i] for i ∈ N) )
for i ∈ N
    @constraint(model, x - r[i] + s[i] == p[i] )
end
print(model)
optimize!(model)
TDᵒ = objective_value(model)
xᵒ = value.(x)
TDᵒ, xᵒ

$$ \begin{aligned} \min\quad & r_{1} + s_{1} + r_{2} + s_{2} + r_{3} + s_{3} + r_{4} + s_{4} + r_{5} + s_{5}\\ \text{Subject to} \quad & x - r_{1} + s_{1} = 1.0\\ & x - r_{2} + s_{2} = 2.0\\ & x - r_{3} + s_{3} = 4.0\\ & x - r_{4} + s_{4} = 6.0\\ & x - r_{5} + s_{5} = 9.0\\ & r_{1} \geq 0.0\\ & r_{2} \geq 0.0\\ & r_{3} \geq 0.0\\ & r_{4} \geq 0.0\\ & r_{5} \geq 0.0\\ & s_{1} \geq 0.0\\ & s_{2} \geq 0.0\\ & s_{3} \geq 0.0\\ & s_{4} \geq 0.0\\ & s_{5} \geq 0.0\\ \end{aligned} $$

(12.0, 4.0)