Linear Programming: Optimisation at the Edges
Subscribe to my newsletter to hear about my latest posts. No spam, I promise.
In continuous optimisation problems, it is impossible to try all combinations. But what if a problem could be solved simply by evaluating the corners of the feasible region? This is exactly what linear programming aims to do.
In the previous parts of this series, we have applied no constraints to the number of tables and chairs produced, allowing it to vary within a positive range. Let’s now assume that the company takes its current machine and employee costs as constant. Production is now constrained by the number of hours of carpentry and varnishing.
The company has a total of
- Chairs:
hour of carpentry, hours of varnishing - Tables:
hours of carpentry, hours of varnishing
Formalising constraints we get:
- Carpentry constraint:
- Varnishing constraint:
Let’s specify the prices of each:
- Price per table:
- Price per chair:
With respect to these constraints, we are aiming to maximise revenue, defined by the following function (based on the unit prices):
Assuming that any table or chair produced will be sold (this is a very popular brand).
How can we optimise this without trying all possible combinations?
As both the objective function and the constraints are linear (i.e., their graph is a straight line), we can leverage linear programming to evaluate only a few of the solutions.
The constraints define the feasible region of the problem, shown in the chart below:
Code to generate the plot
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import ConvexHull
# Define the inequalities as functions
def ineq1(C):
return (200 - C) / 2 # From C + 2T ≤ 200
def ineq2(C):
return (300 - 2 * C) / 2 # From 2C + 2T ≤ 300
# Objective function for total revenue (for contour plot)
def objective(C, T):
return 250 * C + 400 * T # Revenue function
# Generate a grid of C and T values
= np.linspace(0, 200, 400)
C = np.linspace(0, 200, 400)
T = np.meshgrid(C, T)
C_grid, T_grid
# Calculate the corresponding T values for each inequality
= ineq1(C)
T1 = ineq2(C)
T2
# Create a mask for the feasible region
= T_grid <= ineq1(C_grid)
ineq1_mask = T_grid <= ineq2(C_grid)
ineq2_mask = (C_grid >= 0) & (T_grid >= 0)
non_negative_mask = ineq1_mask & ineq2_mask & non_negative_mask
feasible_region
# Calculate the objective function within the feasible region
= objective(C_grid, T_grid)
Z = np.ma.array(Z, mask=~feasible_region)
Z_masked
# Plot the feasible region with contour plot
=(10, 8))
plt.figure(figsize0, np.minimum(T1, T2)), color='gray', alpha=0.5, label='Feasible Region')
plt.fill_between(C, np.maximum(
# Plot the lines defining the inequalities
=r'$C + 2T \leq 200$ - Carpentry Hours')
plt.plot(C, T1, label=r'$2C + 2T \leq 300$ - Varnishing Hours')
plt.plot(C, T2, label
# Identify and plot the vertices of the feasible region
= np.array([[0, 100], [100, 50], [150, 0], [0, 0]]) # Vertices in (C, T)
vertices
# Plot the feasible region boundary
= ConvexHull(vertices)
hull for simplex in hull.simplices:
0], vertices[simplex, 1], 'k-')
plt.plot(vertices[simplex,
# Plot the vertices
0], vertices[:, 1], 'ro', markersize=8, label='Vertices')
plt.plot(vertices[:,
# Labels and legend
'Chairs (C)', fontsize=15)
plt.xlabel('Tables (T)', fontsize=15)
plt.ylabel(0, 201, step=10), fontsize=12)
plt.xticks(np.arange(0, 201, step=10), fontsize=12)
plt.yticks(np.arange(=12)
plt.legend(fontsize'Feasible Region', fontsize=15)
plt.title(0, 200)
plt.xlim(0, 200)
plt.ylim( plt.show()
As the constraints are linear, the shape they form is a complex polytope, a pointy shape in a multidimensional space. Most importantly, this shape has corners. Because the objective function is also linear, it will necessarily intersect with the feasible region at one of its corners.
Below is the feasible region with a contour plot of different values of the objective function. As you can see, the intersection of the objective function and the feasible region is at a corner of the feasible region.
Code to generate the plot
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import ConvexHull
# Define the inequalities as functions
def ineq1(C):
return (200 - C) / 2 # From C + 2T ≤ 200
def ineq2(C):
return (300 - 2 * C) / 2 # From 2C + 2T ≤ 300
# Objective function for total revenue (for contour plot)
def objective(C, T):
return 250 * C + 400 * T # Revenue function
# Generate a grid of C and T values
= np.linspace(0, 200, 400)
C = np.linspace(0, 200, 400)
T = np.meshgrid(C, T)
C_grid, T_grid
# Calculate the corresponding T values for each inequality
= ineq1(C)
T1 = ineq2(C)
T2
# Create a mask for the feasible region
= T_grid <= ineq1(C_grid)
ineq1_mask = T_grid <= ineq2(C_grid)
ineq2_mask = (C_grid >= 0) & (T_grid >= 0)
non_negative_mask = ineq1_mask & ineq2_mask & non_negative_mask
feasible_region
# Calculate the objective function within the feasible region
= objective(C_grid, T_grid)
Z = np.ma.array(Z, mask=~feasible_region)
Z_masked
# Plot the feasible region with contour plot
=(10, 8))
plt.figure(figsize0, np.minimum(T1, T2)), color='gray', alpha=0.5, label='Feasible Region')
plt.fill_between(C, np.maximum(
# Plot the lines defining the inequalities
=r'$C + 2T \leq 200$ - Carpentry Hours')
plt.plot(C, T1, label=r'$2C + 2T \leq 300$ - Varnishing Hours')
plt.plot(C, T2, label
# Identify and plot the vertices of the feasible region
= np.array([[0, 100], [100, 50], [150, 0], [0, 0]]) # Vertices in (C, T)
vertices
# Plot the feasible region boundary
= ConvexHull(vertices)
hull for simplex in hull.simplices:
0], vertices[simplex, 1], 'k-')
plt.plot(vertices[simplex,
# Plot the vertices
0], vertices[:, 1], 'ro', markersize=8, label='Vertices')
plt.plot(vertices[:,
# Add contour plot for the objective function
= plt.contour(C, T, Z, levels=10, colors='blue', alpha=0.7,label="Objective function")
contours =True, fontsize=12, fmt='%1.0f')
plt.clabel(contours, inline
# Labels and legend
'Chairs (C)', fontsize=15)
plt.xlabel('Tables (T)', fontsize=15)
plt.ylabel(0, 201, step=10), fontsize=12)
plt.xticks(np.arange(0, 201, step=10), fontsize=12)
plt.yticks(np.arange(=12)
plt.legend(fontsize'Feasible Region with Revenue Contours', fontsize=15)
plt.title(0, 200)
plt.xlim(0, 200)
plt.ylim( plt.show()
From this observation, there is no need to evaluate all combinations; one simply has to evaluate the corners of the feasible region.
To identify the corners of the feasible region, we consider the constraints:
- Carpentry hours:
- Varnishing hours:
- Non-negativity:
,
Solving these constraints, we find the vertices (corners) of the feasible region:
Origin point:
At
:
- From the carpentry constraint:
- From the varnishing constraint:
- The more restrictive is
, so we have the point
- At
:
- From the carpentry constraint:
- From the varnishing constraint:
- The more restrictive is
, so we have the point
- Intersection point of the two constraints:
- From the carpentry constraint:
- Substitute into the varnishing constraint:
Simplify: Substitute back into : So the intersection point is
Therefore, the corners of the feasible region are at:
We then evaluate the objective function at each corner:
At
:At
:At
:At
:
The maximum revenue occurs at
This type of shortcut is only possible when both the objective function and the constraints are linear.
This may still seem abstract when talking about chairs and tables. But many other problems can also be solved using the same framework:
- Scheduling matinenance staff shifts in a hotel
- Combination of raw materials in industrial food production
- Machine time allocation in a factory
Final Thoughts
The efficiency of linear programming reminds us that the more information we have about a problem, the less computation we need to solve it. This information about the problem (meta-information) can be leveraged to devise smart algorithms taking shortcuts.
But what happens when you cannot make any assumptions? What if you can only observe the output of the objective function? If this sounds exciting, head over to this series’ next post on black-box optimisation problems.