Production mix - Model 1, Pyomo concrete

1 July 2022

In this article we build Model 1 of the Python Production mix series, using the Pyomo library. The Production mix model relates to our hypothetical boutique pottery business, which is described in more detail in the article Python optimization Rosetta Stone.

Our objective for this article is to build and explain the workings of a simple Pyomo example.

Although Model 1 is like many examples of a Pyomo model that you may find on the web and in textbooks, it certainly does not represent best (or even good) practice. In subsequent articles, we'll incrementally improve the model, leading to a structure that is more suitable for operational models.

Even so, Model 1 is an easy-to-understand place to start our exploration of building optimization models in Python. It is worth understanding this model before moving on to more sophisticated versions.

Articles in this series

Articles in the Python Production mix series:

• Python optimization Rosetta Stone
• Production mix - Model 1, Pyomo concrete
• Production mix - Model 2, Pyomo separate data
• Production mix - Model 3, Pyomo external data
• Production mix - Model 4, Pyomo json file
• Production mix - Model 5, Pyomo using def
• Production mix - Model 6, Pyomo abstract
• Production mix - Model 7, PuLP
• Production mix - Model 8, OR-Tools
• Production mix - Model 9, Gekko
• Production mix - Model 10, CVXPY
• Production mix - Model 11, SciPy
• Production mix - Conclusions

Download the models

The model is available on GitHub.

Formulation for Model 1

For this first model, we're using a very direct coding approach, implementing the specific formulation (described in Python optimization Rosetta Stone), with hard coded data, as shown in Figure 1.

Model 1 Python code

Import dependencies

The first task is to import the libraries that are needed for our program. In this case, as shown in Figure 2, the only dependency is the Pyomo library.

import pyomo.environ as pyo

For some models we'll use several libraries, such as NumPy, pandas, json, os, matplotlib, and others. But this is a simple model, so we don't need anything else.

Note that we import the library using the alias pyo. In our code, we'll use the pyo prefix for all Pyomo objects. This style of explicitly referring to library objects can seem verbose, especially as it isn't necessary for small, simple models. But for larger and more complex models, using explicit references produces clearer code and avoids potential conflicts between libraries.

For example, if a model needs to calculate the sine of a variable that represents an angle, then we could try using the sin function from the math library. But that won't work in a Pyomo model because it is the solver, rather than Python, that needs to do the calculation. Instead, we need to use the sine function from the Pyomo library, like pyo.sin. Including the prefix makes it clear exactly which library we're intending to use.

Declarations

Next, we declare the model. Pyomo has two types of model:

• Concrete. The model is defined using a known set of data values.
• Abstract. The model is defined using only symbols, without knowing the data values. The values are supplied when the model is solved.

In this case, as shown in Figure 3, we declare our model as a concrete model. We'll explore an abstract implementation of the Production mix model in a later article.

Model = pyo.ConcreteModel(name = 'Boutique pottery shop - Model 1')

The declaration creates a Pyomo object that we'll use to access various features of the library. Note that the object's name, "Model", is arbitrary – though beware that the object name is case sensitive.

It is a good practice to give each model a descriptive name. Therefore, when we declare the model, we give it the name 'Boutique pottery shop - Model 1'. The model's name can be shown in the results, making it easy to identify where the results come from.

As we'll be hard coding the data values in this model's objective function and constraints, there are no further data structures to declare.

Define the model

With the preliminaries done, we can now define the model. We do this, as shown in Figure 4, by adding variables, constraints, and an objective function to the Model object:

• Variables. We have variables for the quantity of each product, Discs and Orbs. The variables represent physical production quantities, so we specify the variables to be in the domain pyo.NonNegativeReals, meaning that they can take any value greater than or equal to zero. We'll solve this model as a linear program, so part units (like 3.76 Discs) are acceptable. If we wanted to limit the production variables to integer quantities, then we could specify the domain as pyo.NonNegativeIntegers. A variety of other Predefined Virtual Sets are also available.
• Constraints. We express the People, Materials, and Sales constraints as expressions expr. The expressions closely follow the specific formulation shown in Figure 1, with hard coded coefficients for each variable term and the right-hand side of the constraints.
• Objective function. The objective function is also defined using an expression, though it does not have a right-hand side. We specify the sense of the objective to indicate whether we want to pyo.minimize or pyo.maximize the objective function.

Model.Discs = pyo.Var(domain = pyo.NonNegativeReals)
Model.Orbs = pyo.Var(domain = pyo.NonNegativeReals)

Model.PeopleHours = pyo.Constraint(expr = 12.50 * Model.Discs + 10.00 * Model.Orbs <= 250)
Model.MaterialUsage = pyo.Constraint(expr = 18.00 * Model.Discs + 30.00 * Model.Orbs <= 500)
Model.SalesRelationship = pyo.Constraint(expr = -2.00 * Model.Discs + 1.00 * Model.Orbs <= 0)

Model.TotalMargin = pyo.Objective(expr = 80.00 * Model.Discs + 200.00 * Model.Orbs, sense = pyo.maximize)

Solve model

Now we can solve the model, as shown in Figure 5.

Pyomo does not have any built-in solvers – that's why we installed solvers for Pyomo when we set up our Python modelling environment, as described in the article Setting up a Python modelling environment. Therefore, we need to specify a solver – which we do in the first line. Our model is a linear program, so the CBC solver is appropriate. We could also use the GLPK solver, simply by changing 'cbc' to 'glpk'.

The second line solves the model, assigning the outcome to a Results object.

Solver = pyo.SolverFactory('cbc')
Results = Solver.solve(Model)

Write output

All going well, the model will solve successfully. So now we can write the results, as shown in Figure 6.

Note that we should check that the solver was successful in finding at least a feasible solution. We should also capture and handle any errors that may occur. However, for the sake of simplicity we'll forgo those checks in this model – we'll add that level of sophistication in later models.

Our model is simple, so there isn't much output required. In general, we want to know:

• Model name. This is especially useful if the results are copied elsewhere, so that we can easily identify which model the results came from.
• Solver status. Different modelling libraries and different solvers may return a variety of status codes. Here we simply write the solver's status value. In more complex models, it is helpful to interpret the status code, rather than simply writing it to the results.
• Objective function value. We report the objective function value, after solving, using our name for the objective function Model.TotalMargin(). It is also useful to include the $ units of the objective function, to help with interpretation.
• Variable values. Similarly, we write the values of the variables.

print(Model.name, '\n')
print('Status: ', Results.solver.status, '\n')
print(f'Total margin = ${Model.TotalMargin():,.2f}')
print(f'Production of discs = {Model.Discs():6.2f}')
print(f'Production of orbs  = {Model.Orbs():6.2f}')

The output is shown in Figure 7. The "ok" status means that an optimal solution was found by the CBC solver. We note that production of Orbs is exactly twice the production of Discs, indicating that the Sales constraint is binding. This solution is the same as that found by the original Excel model.

'Boutique pottery shop - Model 1'

Status:  ok

Total margin = $3,076.92
Production of discs =   6.41
Production of orbs  =  12.82

Evaluation of this model

This is a very simply linear programming model. The entire model – preliminary setup, defining the model, solving the model, and writing the results – needs only 15 lines of Python code.

By using a concrete model, with hard coded coefficients, the Pyomo library makes the modelling process very easy. In particular, the constraint and objective function expressions are an almost literal implementation of the specific formulation.

Many examples of Pyomo concrete models in blogs and textbooks look much like Model 1. Often, such models are presented as if this is how you should build models in Python. However, in practice, we would never build an operational model like Model 1.

Hard coding coefficients directly in the constraint and objective function expressions is a poor practice that makes changing data unnecessarily difficult, especially if the data comes from a database or other process. For a larger model, hard coding coefficients is – at best – cumbersome. In addition, changing hard coded values manually is risky, with high potential for inadvertently changing (or breaking) the model definition rather than just changing the data.

The model is also very rigid, with each product's variable individually named. For larger models, variables are usually defined as indexed sets rather than being named individually.

Consequently, while this model design is useful for introducing some Python modelling concepts, it lacks the flexibility and robustness needed for an operational model.

Next steps

To start addressing the limitations of Model 1, in the next article we'll separate the data from the model definition. That is, we'll replace the hard coded coefficients with Python data structures designed to contain the values.

Conclusion

This article starts our exploration of optimization modelling in Python. Model 1 represents the Production mix situation as a concrete Pyomo model, with hard coded coefficients. While this type of model is common, it is not good practice. In the next article, we'll start improving the design.

If you would like to know more about this model, or you want help with your own models, then please contact us.