top of page

Linear Solving In Python

Updated: Jan 16, 2024


Introduction

Earlier this year, I took a class in finances in excel class on Coursera and if you know me I love learning how to solve real world problems programmatically and mathematically.


One of the concepts I got really good at was working through linear and non-linear optimization problems. They are used to solve problems to help manage what products to sell based on time of production, cost of materials, and return on investment. (Sounds fun right?)




Half way through the class, I learned how to use Excel Solvers and was super interested in how I would be able to write my own.

This tutorial is an introduction to using python to start solving these problems programmatically and how we can start including this code in larger projects.




The Basics


Linear equations is defined as:

A linear equation is an algebraic equation where each term has an exponent of 1 and when this equation is graphed, it always results in a straight line.


When you see the equations for linear equations, they look like ( y = 8x - 9 ). These are the classic ( y = mx + b ) type of equation to understand the slope and the intercept for the linear line on the graph.


Now if we are trying to solve an optimization problem for linear equations there are usually several inputs and constraints to the problem.


The vocabulary that is useful to understand is:

  • Constraints These are generally equations to show the area of interest in a problem.

  • Objective This is what we are trying to solve. This can be seen in trying to find the maximum or minimum of a given value.


So usually these linear equation problems can be described as based on these given constraints what is the (max/min) value we can find. In future tutorials we will go over how we can express these problems mathematically, but in this intro we are just going to focus on the use of Python and how to get started with the problems.





Getting Started


For our first problem to solve, we are going to start by expressing 3 different constraints and then find the maximum x and y value that can be achieved given the constraints. First we want to list the constrain and then the objective

Constraints:

  • 2y+-1x<=8

  • y+2x<=14

  • -1y+2x<=10

Objective:

  • max (x, y)




Graphing the Equations


The first thing we want to do is make a graph of the equation. This helps to visualize the problem we are facing. In Python we want to do this with matplotlib.

We first import the package as well as numpy. We create the graph's space and then we can start plotting the lines in the graph.

The lines y1, y2, y3 need to be expressed in the format y=mx+b where y is an array of position for every x position in the graph. We then plot each of our linear lines and define them in the key. Last we give more context to the graph providing labels for the x and y axis as well as defining the view and key.


import matplotlib.pyplot as plt
import numpy as np

graph_x = np.linspace(0, 100, 2000)

y1 = (8 - graph_x * -1 ) / 2
y2 = (14 - graph_x * 2)
y3 = (10 - graph_x * 2) / -1

plt.plot(graph_x, y1, label=r'$-X+2Y\leq8$')
plt.plot(graph_x, y2, label=r'$2X+Y\leq14$')
plt.plot(graph_x, y3, label=r'$2X-Y\leq10$')

plt.xlim((0, 12))
plt.ylim((0, 12))
plt.xlabel('x values')
plt.ylabel('y values')

plt.legend(bbox_to_anchor=(.65, 0.95), loc=2, borderaxespad=0.1)

Running the application we can now see the graph displayed below.




Creating the Linear Equation Solver


The linear solver is the bit of code this tutorial is all about. This is the bit of code that based on the constraints of the graph, can figure out the values we want to find based on a given parameter.

For this tutorial we are going to be using pyomo. This package can do linear and nonlinear equations and the best part is its free. There are other solvers out there but most of them need a subscription / membership and if you are anything like me, open source is always better.


Installation

First we need to add the package to our python library.


pip install pyomo

The next step if you have to download the installer you are looking to use as part of the pyomo package. The one that I chose to use was GLPK 4.65. You want to download the package off of source forge and then unzip it and place it on your computer. Next you need to add the path to the PATH environment variable and restart your machines to make the new path work. If everything worked correctly we can now write the code to implement the solver.


import pyomo.environ as pyo
from pyomo.environ import *
from pyomo.opt import SolverFactory

model = pyo.ConcreteModel()
model.x = pyo.Var(bounds=(0,10))
model.y = pyo.Var(bounds=(0,10))

y = model.y
x = model.x

model.C1 = pyo.Constraint(expr= -x+2*y<=8)
model.C2 = pyo.Constraint(expr= 2*x+y<=14)
model.C3 = pyo.Constraint(expr= 2*x-y<=10)

model.obj = pyo.Objective(expr= x+y, sense=maximize)

opt = SolverFactory('glpk')
opt.solve(model)

x_value = pyo.value(x)
y_value = pyo.value(y)


Displaying The Results




The last thing I want to display is the location of the plot for the optimal solution in the graph. We can do this by first creating a shaded region for the area that covers the area within our constraint. To add this in the code above for our graph, add the following lines before plotting the the lines on the graph.


d = np.linspace(-2,16,300)
x,y = np.meshgrid(d,d)
plt.imshow( ((y>=0) & (x>=0) & (y<=14-2*x) & (2*y<=8-x*-1) & (-1*y<=10-2*x)).astype(int) , 
                extent=(x.min(),x.max(),y.min(),y.max()),origin="lower", cmap="Greys", alpha = 0.3);

What we do in this code is create a mesh grid in numpy and then within that space only show grid points that are inside the constraints we have listed in our linear lines. These are statements and need to be written as such unlike the array we create in the graphing sample above.


The last bit of code we want to add is the plot and a label for our Solver's solution. we want to take the references of the output and use them in the graph by adding these lines.


plt.plot(x_value, y_value, 'o')
plt.text(x_value, y_value, f'({x_value}, {y_value})')
plt.show()

This is the output of the graph with the desired results:





Conclusion



This is only an introduction to how to use solvers and to be able to express these problems. In the next tutorial we will show how these problems are expressed in different ways and how the solvers can be used in real world scenarios.

bottom of page