Solver Max - Production mix - Model 2, Pyomo separate data

14 July 2022

In this article we continue the Python Production mix series, using the Pyomo library. Specifically, we build Model 2, which improves Model 1 by extracting the hard coded coefficients and putting them into Python data structures.

Separately the data from the model is an important step in developing a practical design that can be used for operational models.

Articles in this series

Articles in the Python Production mix series:

Download the models

The model is available on GitHub.

Formulation for Model 2

For this model, by separating the data from the model, we're using the general (rather than specific) formulation, as shown in Figure 1.

Figure 1. General formulation

\begin{alignat}{1} &\text{Objective} \\ & \quad \text{Maximize} \ fTotalMargin &= \quad & \sum_{p=1}^n \text{dMargin}_{p} \times vProduction_{p} \tag{1} \\ \\ &\text{Subject to} \\ & \quad \sum_{p=1}^n \text{dPeople}_{p} \times vProduction_{p} &\le \quad & \text{dHours} \tag{2}\\ & \quad \sum_{p=1}^n \text{dMaterials}_{p} \times vProduction_{p} \quad &\le \quad & \text{dkg} \tag{3}\\ & \quad \sum_{p=1}^n \text{dSales}_{p} \times vProduction_{n} &\le \quad & 0 \tag{4}\\ \\ &\text{Variables} \\ & \quad vProduction_{p} &= \quad & \text{Units of each product} \quad &\forall \ p \in P \tag{5}\\ \\ &\text{Data} && \\ & \quad \text{dMargin}_{p} &= \quad &\text{Profit margin on each unit} \quad &\forall \ p \in P \tag{6}\\ & \quad \text{dPeople}_{p} &= \quad &\text{Staff hours to make each unit} \quad &\forall \ p \in P \tag{7}\\ & \quad \text{dMaterials}_{p} \quad &= \quad &\text{kg of materials to make each unit} \quad &\forall \ p \in P \tag{8}\\ & \quad \text{dSales}_{p} &= \quad &\text{Relative sales volume for each product} \quad &\forall \ p \in P \tag{9}\\ \\ &\text{Sets} \\ & \quad P &\ &\text{Product} \tag{10} \\ \\ &\text{Dimensions} \\ & \quad n &\ &\text{Number of products} \tag{11} \\ \\ &\text{Bounds} \\ & \quad \text{Non-negative} &\ &\quad \tag{12} \\ \end{alignat}

Model 2 Python code

Import dependencies

The first task is to import the libraries that are needed for our program. As shown in Figure 2, in addition to the Pyomo library, we import the pandas library. We'll use pandas when writing the model's output.

Figure 2. Import dependencies

import pyomo.environ as pyo
import pandas as pd

Declarations

As shown in Figure 3, we continue to use a Pyomo concrete model.

Figure 3. Declarations

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

Hours = 250
kg = 500
SalesLimit = 0

Coefficients = {
    'Discs' : {'People': 12.50, 'Materials': 18.00, 'Sales': -2.00, 'Margin':  80.00},  
    'Orbs'  : {'People': 10.00, 'Materials': 30.00, 'Sales':  1.00, 'Margin': 200.00},
}

Products = Coefficients.keys()

VarInitial = 0
VarBounds = (0, 100)

Unlike Model 1, here we separate the data from the model. Therefore, we need to declare and populate some data structures.

We start with simple constants, called Hours, kg, and SalesLimit, which will be the right-hand side values of our constraints.

Next, we declare and populate a Python Dictionary called Coefficients to contain the coefficients for each variable. Each row of the dictionary represents a product, enabling us to easily add more products. Within each row, we have the People, Materials, and Sales coefficients that we'll use in the constraints, along with the Margin coefficient that we'll use in the objective function. If we added a new constraint, then we would add the associated variable coefficients as a new column in this dictionary.

Finally, we create a dictionary keys object called Products, which contains the dictionary's top-level keys – i.e. Discs and Orbs. We'll use these keys to select the coefficients for each product.

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. In Model 1, we had a variable for each product. Here we define a set of variables, indexed by Products. Since Products contains two keys, two variables are created. We also initialize each of the variables, just to show how it is done, though it doesn't matter in this model.
Constraints. Rather than specifying individual product variables, for the left-hand side of each constraint we sum over all Products. For example, when the index p has the value Discs, the term Coefficients[p]['People'] has value 12.50 (see Figure 3). The right-hand side of each constraint is the value we specified in the declarations section above.
Objective function. The objective function is defined in a similar way.

Figure 4. Define the model

Model.Production = pyo.Var(Products, domain = pyo.NonNegativeReals, initialize = VarInitial, bounds = VarBounds)

Model.PeopleHours = pyo.Constraint(expr = sum(Coefficients[p]['People'] * Model.Production[p] for p in Products) <= Hours)
Model.MaterialUsage = pyo.Constraint(expr = sum(Coefficients[p]['Materials'] * Model.Production[p] for p in Products) <= kg)
Model.SalesRelationship = pyo.Constraint(expr = sum(Coefficients[p]['Sales'] * Model.Production[p] for p in Products) <= SalesLimit)

Model.TotalMargin = pyo.Objective(expr = sum(Coefficients[p]['Margin'] * Model.Production[p] for p in Products), sense = pyo.maximize)

Check the model definition

It is important to check that our model definition works correctly. Python has a "pretty print" module, pprint, for outputting data structures in a formatted way. We can call pprint on our Pyomo Model object by simply running the line Model.pprint(). The output is shown in Figure 5.

As a quick check, the important parts in the pprint output are:

We have two variables, Production[Discs] and Production[Orbs].
The objective function is: 80.0*Production[Discs] + 200.0*Production[Orbs].
The Material constraint is: 18.0*Production[Discs] + 30.0*Production[Orbs], with a right-hand side of 500.0. The other constraints are similar.

Each of these outputs is as expected.

Figure 5. Output of Model.pprint()

1 Set Declarations
    Production_index : Size=1, Index=None, Ordered=False
    Key  : Dimen : Domain : Size : Members
    None :     1 :    Any :    2 : {'Discs', 'Orbs'}

    1 Var Declarations
    Production : Size=2, Index=Production_index
    Key   : Lower : Value : Upper : Fixed : Stale : Domain
    Discs :     0 :   0.0 :  None : False : False : NonNegativeReals
    Orbs :     0 :   0.0 :  None : False : False : NonNegativeReals

    1 Objective Declarations
    TotalMargin : Size=1, Index=None, Active=True
    Key  : Active : Sense    : Expression
    None :   True : maximize : 80.0*Production[Discs] + 200.0*Production[Orbs]

    3 Constraint Declarations
    Material : Size=1, Index=None, Active=True
    Key  : Lower : Body                                           : Upper : Active
    None :  -Inf : 18.0*Production[Discs] + 30.0*Production[Orbs] : 500.0 :   True
    People : Size=1, Index=None, Active=True
    Key  : Lower : Body                                           : Upper : Active
    None :  -Inf : 12.5*Production[Discs] + 10.0*Production[Orbs] : 250.0 :   True
    Sales : Size=1, Index=None, Active=True
    Key  : Lower : Body                                      : Upper : Active
    None :  -Inf : -2.0*Production[Discs] + Production[Orbs] :   0.0 :   True

    6 Declarations: Production_index Production People Material Sales TotalMargin

Solve model

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

This code is the same as Model 1.

Figure 6. Solve model

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

Write output

As for Model 1, we output the model's name, the solve status, and the objective function value – as shown in Figure 7.

But instead of writing the variables individually, here we use a pandas DataFrame to collate the variable values. Note that we access the variable values using pyo.value(Model.Production[p]), where p in an index of the Products keys.

Having collated the results in a DataFrame, we simply display the DataFrame, which writes a nicely formatted table of results.

Figure 7. Write solution

print(Model.name, '\n')
print('Status: ', Results.solver.status, '\n')
print(f'Total margin = ${Model.TotalMargin():,.2f}\n')
ProductResults = pd.DataFrame()
for p in Products:
    ProductResults.loc[p, 'Production'] = round(pyo.value(Model.Production[p]), 2)
display(ProductResults)

The output is shown in Figure 8. Apart from formatting, the solution is the same as for Model 1.

Figure 8. Solution

    Boutique pottery shop - Model 2

    Status:  ok 

    Total margin = $3,076.92

    Production
    Discs        6.41
    Orbs        12.82

Evaluation of this model

Model 2 is an improvement compared with Model 1. Pyomo enables us to create and use data structures that are a step towards a more general model design.

Tools like pprint() help us validate that our model works correctly. Such validation will become more important as our model design becomes more sophisticated.

There is still more that we can do to improve the model, which we'll continuing exploring in the next article of this series.

Next steps

Although we've separated the data from the model, the data is still contained in the same file. In the next article, we'll put the data in an external file. We'll also improve the way we handle the solution process.

Conclusion

This article continues our exploration of optimization modelling in Python. We've improved Model 1, with a clear separation of the data from the model.

In the next article, we'll continue 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.

Python
Pyomo

Production mix - Model 2, Pyomo separate data