## Linear optimization main concepts and implementation in Python

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:

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

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

**as well. The objective can be defined either in a**

*x**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.

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*.

And can then be stated as the following.

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.

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

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

**to satisfy this condition (Nocedal & Wright, 2006). If we can identify a feasible point**

*A***with at most**

*x**m*nonzero components and the corresponding matrix

**(**

*B**m*×

*m*) of the columns from

**matching nonzero indices in**

*A***is nonsingular,**

*x***is called a**

*x**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):

- If there is a
*feasible solution*, there is a*basic feasible solution*. - 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.

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.

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.

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:

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.

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.

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:

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.

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:

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.

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.