Linear programming: Theory and applications

Linear optimization main concepts and implementation in Python

Photo by Patrick Fore on Unsplash

Numerical optimization is a fundamental tool in quantitative decision-making processes. To optimize a given goal, one must select values for interrelated decision variables under some prescribed set of circumstances. The format of the mathematical equations that describe the objective and the constraints of the problem is used to distinguish optimization problems between two major categories: Linear and Nonlinear Programming.

Management sciences and operations research make extensive use of linear models, whereas nonlinear programming problems tend to arise naturally in the physical sciences and engineering (Nocedal & Wright, 2006). In linear problems, as the name suggests, the objective(s) and constraints are described by linear functions only, which will be the focus of the current article.

Throughout this article, some of the main theoretical aspects of linear programming will be covered, besides applications in classical problems using Python. To do this, we will use the libraries scipy and pyomo (Bynum et al., 2021).

You may find here a short summary of the remainder of this text:

• Problem statement
• Feasible space and solution techniques
• The product mix problem
• The transportation problem
• Further reading
• Conclusions
• References

You can find the complete code in this GIT repository. Now enjoy the ride!

Problem statement

When formulating an optimization problem, one must define an objective that is a function of a vector decision variables x and might be subject to some equality and inequality constraints, which are functions of x as well. The objective can be defined either in a minimization or maximization sense although the former is the most usual. Notice that by simply multiplying the coefficients cᵢ of a maximization objective, it can be re-formulated in a minimization sense.

In linear problems, the decision variable space must be limited somehow. Otherwise, decision variables would converge to positive or negative infinity values according to the objective function. Constraints might be formulated by either equality or inequality relationships. They are usually stated as rows in the problem matrix. Let us distinguish the matrix of equality constraints A_eq from the matrix of inequality constraints A_ub.

Lower and upper boundaries for each component of x might be explicit in the formulation, which reduces the search space. As a convention, due to solution techniques, usually lower bounds for decision variables are by default equal to zero. This leads to a general problem formulation such as the following.

Linear problem in generic formulation. (Image by the author).

Inequality constraints might be re-formulated as equality constraints by adding non-negative slack variables. If this occurs and all decision variables are stated non-negative, we say the linear program is in the standard form.

Inequality constraint re-formulated as equality by including a nonnegative slack variable. (Image by the author).

And can then be stated as the following.

Linear programming problem in standard form. (Image by the author).

Notice that inequalities work as equalities if their corresponding slack variables are equal to zero. We will dive deeper into that in the next section when discussing the feasible space and solution techniques.

Feasible space and solution techniques

Let us introduce a simple example to explore the concepts of this section.

Example problem. (Image by the author).

Which can be represented in a two-dimensional space as the following.

Notice that in each vertice of the feasible region (colored), a subset of constraints is active. Moreover, observe that the optimal solution in this example lies in a vertice. Recall that, when inequalities are active, their corresponding slack variables for the problem in the standard form are equal to zero. In linear problems, we will be interested in finding which constraints are active in a vertice that contains the optimal solution.

A generalization for higher dimensional problems is that, if an optimal solution exists, the set of optimal points can be a single vertice (unique solution), an edge or face, or even the entire feasible set (multiple solutions). The problem has no solution if the feasible set is empty (the infeasible case) or if the objective function is unbounded on the feasible region (the unbounded case).

In a formal definition, suppose the matrix A of the problem in its standard form has full row rank. In practice, preprocessing usually is applied to remove some redundancies. Reformulation by adding slack, surplus, and artificial variables is also used to force A to satisfy this condition (Nocedal & Wright, 2006). If we can identify a feasible point x with at most m nonzero components and the corresponding matrix B (m × m) of the columns from A matching nonzero indices in x is nonsingular, x is called a basic feasible point or a basic feasible solution. The variables set to zero are denoted nonbasic variables, while the remaining are denoted basic variables.

Back to the feasible region… The feasible set defined by the linear constraints is a polytope (in the two-variable example provided it is a polygon). Algebraically, all its vertices correspond to basic feasible points previously described. We can, therefore, make a connection between the algebraic and geometric viewpoints of the problem and gain insight into understanding how the simplex method works.

Variants of the simplex method are probably the most widely used solution techniques of linear programming. They explore the fundamental theorem of linear programming (Luenberger & Ye, 2008):

1. If there is a feasible solution, there is a basic feasible solution.
2. If there is an optimal feasible solution, there is an optimal basic feasible solution.

Thus, if an optimal solution exists, there is a vertice that contains an optimal solution. Aiming to find such a vertice, steps in the simplex method explore adjacent vertices, for which the set of basic indices differs in exactly one component until no further improvement is achieved. The reader interested in understanding better how to obtain a starting feasible point, prove the optimality of a solution, and how to iterate between two consecutive vertices might refer to Nocedal & Wright (2006) or Winston & Goldberg (2004).

Interior point methods are also widely used, especially for large linear programs. Interior-point methods share common features that distinguish them from the simplex method. Each interior-point iteration is expensive to compute and can make significant progress toward the solution, while the simplex method usually requires a larger number of inexpensive iterations (Nocedal & Wright, 2006). Some optimization solvers might also use variants of the simplex method to refine solutions obtained by interior point methods.

Both the simplex and interior point methods are represented in the figure below.

Graphical representation of the simplex method (left) and an interior point method (right). (Image by the author).

Now done with some important theoretical aspects, in the next section, we will implement our first problem using Python. To solve it, we will first use the Python package scipy which has wrappers for the open-source solver HiGHS. Furthermore, we will implement the same problem using pyomo (Bynum et al., 2021) and solve it with the CBC solver. Both solvers will make use of variants of the simplex method to obtain their solutions.

The product mix problem

In the product mix problem, we must allocate limited I resources to produce J products maximizing expected returns. Each product j has a known unitary profit margin cⱼ and requires a known amount of each resource fᵢⱼ. Consider the availability of each resource bᵢ. Our decision variables represent how much of each product should be produced xⱼ. We can formulate the problem as the following.

Product mix problem. (Image by the author).

From this structure, we can easily convert the problem to a matrix form, in which the matrix A_ub has elements fᵢⱼ in row i and column j.

Consider an example in which we have three reactants {A, B, C} used to produce three products {D, E, F}. The availability of each reactant is respectively 8000, 3000, and 2000. The profit margin of each product is respectively 2.15, 1.34, 1.72. The proportions are given in the table below.

+---------------------+--------+-------+------+
| Reactant \ Product  |   D    |   E   |   F  |
+---------------------+--------+-------+------+
| A                   |  7/10  |  1/3  |  1/2 |
| B                   |  1/5   |  2/3  |  1/6 |
| C                   |  1/10  |  0    |  1/3 |
+---------------------+--------+-------+------+

Let us start implementing this problem by importing useful libraries.

import numpy as np
from scipy.optimize import linprog

And let us instantiate numpy arrays to represent our parameters.

# Pseudo costs
margins = np.array([2.15, 1.34, 1.72])
c = - margins

# Proportions
A = np.array([
    [7/10, 1/3, 1/2],
    [1/5, 2/3, 1/6],
    [1/10, 0.0, 1/3]
])

# Availability
b = np.array([8000.0, 3000.0, 2500.0])

We are almost done now. All we need to do next is to call the scipy function linprog to solve the problem.

sol = linprog(c, A_ub=A, b_ub=b)
print(sol)

And it should return a solution with properties fun (the objective function) and x (the vector of decision variables).

This problem is more intuitive to formulate using matrix notation than others might be. For more complex problems, algebraic modeling languages (AML) can be very helpful. An amazing open-source alternative available in Python is pyomo (Bynum et al., 2021). So let us apply it to the product mix problem.

import pyomo.environ as pyo

There are two approaches for modeling a problem in pyomo: Abstract and Concrete models. In the first approach, the algebraic expressions of the problem are defined before some data values are supplied, whereas, in the second, the model instance is created immediately as its elements are defined. You can find more about these approaches in the library documentation or in the book by Bynum et al. (2021). Throughout this article, we will adopt the Concrete model formulation.

model = pyo.ConcreteModel()

In this problem there are two sets:

• Resources I
• Products J

When using pyomo we can create these sets as modeling components, which has some benefits of set operations implemented. It goes as the following.

model.I = pyo.Set(initialize=["A", "B", "C"])
model.J = pyo.Set(initialize=["D", "E", "F"])

As we are using the ConcreteModel approach, we must provide values for each component immediately when instantiating it. For sets, this is done via the initialize keyword argument.

Now let us create the fixed parameters of our model. In pyomo, we can use the Param class to store parameters. Parameters can be defined over sets. In this case, we must pass the given sets as the first arguments in the definition. As for other pyomo components, multiple sets can be passed in this statement if the element defined (variable, parameter, expression, or constraint) is indexed by more than one set. This definition is in the Python *args style.

# Availability
model.b = pyo.Param(model.I, initialize=dict(zip(model.I, b)))

# Margins (now as a maximization objective)
model.c = pyo.Param(model.J, initialize=dict(zip(model.J, margins)))

# Proportions
proportions = {}
for k, i in enumerate(model.I):
    for l, j in enumerate(model.J):
        proportions[i, j] = A[k, l]

model.f = pyo.Param(model.I, model.J, initialize=proportions)

Now, let us instantiate the decision variables and constraints. This is done using the Var and Constraint components respectively. Notice that I specified the domain of x as nonnegative reals. This can be done in pyomo via the within keyword argument.

# Decision variables
model.x = pyo.Var(model.J, within=pyo.NonNegativeReals)


# Availability constraint
def availability_rule_cstr(model, i):
    return sum(model.f[i, j] * model.x[j] for j in model.J) <= model.b[i]

model.cstr_available = pyo.Constraint(model.I, rule=availability_rule_cstr)

Observe that the constraint was instantiated passing the rule keyword argument. The rule must be a callable that receives the model as the first argument followed by indices of its corresponding domain.

We are almost done! Let us define the objective function.

# Objective function
def obj_func(model):
    return sum(model.x[j] * model.c[j] for j in model.J)

model.obj = pyo.Objective(rule=obj_func, sense=pyo.maximize)

At last, let us instantiate a solver and use it to solve our model. In this example, I used the open-source solver CBC. You can download CBC binaries from this link. You can also find an installation tutorial here. As the CBC executable is included in the PATH variable of my system, I can instantiate the solver without specifying the path to an executable file. If yours is not, parse the keyword argument “executable” with the path to your executable file.

cbc = pyo.SolverFactory("cbc")
sol = cbc.solve(model, tee=False)

And we are done! Now you can access any modeling component and check its values by using the print or display methods. As expected, the solution using pyomo and CBC returned exactly the same values as the previous one using scipy.

In the following section, let us solve another problem slightly more complex: the transportation problem.

The transportation problem

In the transportation problem, we must match suppliers I of limited capacity bᵢ to a set of customers J to meet their known demands dⱼ. Each pair (i, j) has a known supply cost cᵢⱼ. Our goal is to minimize the total cost by defining how much of each demand is supplied by each supplier. We, therefore, consider decision variables xᵢⱼ to represent those quantities. If we consider conditions such that the total supply capacity meets total demand, we can formulate the problem as the following.

Transportation problem. (Image by the author).

Let us consider an example with four customers {A, B, C, D} and three suppliers {1, 2, 3}. The demand of each customer is respectively 5, 15, 13, and 17. The capacity of each supplier is respectively 14, 26, and 11. Now consider the following costs matrix.

+---+------+------+------+------+
|   |   A  |    B |    C |    D |
+---+------+------+------+------+
| 1 |   10 |    5 |   20 |   12 |
| 2 |   12 |    7 |   12 |   19 |
| 3 |    6 |   12 |   16 |   17 |
+---+------+------+------+------+

Let us instantiate these parameters as Python objects.

costs = pd.DataFrame({
    "A": [10, 12, 6],
    "B": [5, 7, 12],
    "C": [20, 12, 16],
    "D": [12, 19, 17],
}, index=[1, 2, 3])


availability = {1: 14, 2: 26, 3: 11}
demands = {"A": 5, "B": 15, "C": 13, "D": 17}

And let us start implementing our pyomo model.

model = pyo.ConcreteModel()

Once again we have two sets:

• Suppliers I
• Customers J

And the corresponding syntax should be quite similar to the previous example.

model.I = pyo.Set(initialize=[1, 2, 3])
model.J = pyo.Set(initialize=["A", "B", "C", "D"])

Let us instantiate our parameters as pyomo components.

model.b = pyo.Param(model.I, initialize=availability)

model.d = pyo.Param(model.J, initialize=demands)

c = {(i, j): costs.loc[i, j] for i in costs.index for j in costs.columns}
model.c = pyo.Param(model.I, model.J, initialize=c)

Our decision variables…

model.x = pyo.Var(model.I, model.J, within=pyo.NonNegativeReals)

And this time we should have both equality and inequality constraints as each demand must be met (equality) while respecting suppliers' capacities (inequality).

def availability_rule(model, i):
    return sum(model.x[i, j] for j in model.J) <= model.b[i]

model.availability_constr = pyo.Constraint(model.I, rule=availability_rule)


def demand_rule(model, j):
    return sum(model.x[i, j] for i in model.I) == model.d[j]

model.demand_constr = pyo.Constraint(model.J, rule=demand_rule)

Let us define the objective function and once again use CBC to solve the problem. Notice that this time, I defined the objective using the expr keyword argument, which receives directly a pyomo expression.

model.obj = pyo.Objective(expr=sum(model.x[i, j] * model.c[i, j] for (i, j) in model.x), sense=pyo.minimize)
sol = cbc.solve(model, tee=False)

And in this example, I will export my results to a pandas DataFrame for better visualization.

results = pd.DataFrame(index=model.I, columns=model.J)
for i, j in model.x:
    results.loc[i, j] = model.x[i, j].value

Which should return something like this.

+---+---------+----------+----------+----------+
|   |    A    |      B   |      C   |      D   |
+---+---------+----------+----------+----------+
| 1 |     0.0 |      2.0 |      0.0 |     12.0 |
| 2 |     0.0 |     13.0 |     13.0 |      0.0 |
| 3 |     5.0 |      0.0 |      0.0 |      5.0 |
+---+---------+----------+----------+----------+

And in this scenario, all customers would be served with minimum total cost. But what if we had less availability of suppliers than total demand? The problem would then be infeasible and we would need to define which constraints should be violated to return some optimization result. In this problem, a simple suggestion would be to introduce artificial variables to meet the demands with unitary costs greater than those of attending from any of the original suppliers. Remember that the best strategy should be unique for each scenario discussed in real-world problems though.

Further reading

For a deeper understanding of the theoretical aspects of Linear Programming, I strongly advise reading the related chapters in the books by Luenberger & Ye (2008) and Nocedal & Wright (2006). The concept of duality can be especially useful due to sensitivity analysis, an economic interpretation of the problem, and solution techniques (just to mention a few applications).

Operations Research makes extensive use of numerical optimization. The book by Winston & Goldberg (2004) can be very helpful for those interested in more applications of optimization besides some related sciences.

Some real-world problems include discrete variables, to formulate disjunctive problem aspects and integrality. I wrote a few medium articles on the subject, which might be helpful. See here an introduction to branch & bound, the classical knapsack problem, and the job-shop scheduling problem.

Some problems also include nonlinear functions in either the objective or constraints. If that’s the case, you might want to check this other article:

Nonlinear programming: Theory and applications

Conclusions

Throughout this article, some of the most relevant theoretical aspects of linear optimization have been explained in detail and illustrated with two practical implementation examples: the product mix and the transportation problems. The Python library scipy was used to solve the product-mix problem using matrix notation and both problems were modeled using the Python AML pyomo and solved using open-source solvers. The code used is entirely available for readers in this GIT repository.

References

Bynum, M. L. et al., 2021. Pyomo-optimization modeling in python. Springer.

Luenberger, D. G. & Ye, Y., 2008. Linear and Nonlinear Programming. 3rd ed. Stanford: Springer.

Nocedal, J. & Wright, S. J., 2006. Numerical Optimization. 2nd ed. New York: Springer New York.

Winston, W. L. & Goldberg, J. B., 2004. Operations research: applications and algorithms. 4th ed. Belmont, CA: Thomson Brooks/Cole Belmont.