Production mix - Model 4, Pyomo json file

16 August 2022

In this article we continue the Python Production mix series, using the Pyomo library. Specifically, we build Model 4, which changes Model 3 to:

• Import the data from an external json file.
• Read the data into the Model object, rather than into separate objects.

These changes reflect features that we may need to include in an operational model.

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 4

For this model, we're using the same general formulation that we used for Model 2.

Model 4 Python code

Import dependencies

The first task is to import the libraries that are needed for our program. As shown in Figure 1, in addition to the Pyomo and pandas libraries, we import the os and json libraries – which we'll use to import the data from the json file.

import pyomo.environ as pyo
import pandas as pd
import os.path
import json

Get data

The main difference between Model 4 and Model 3 is the data format. Specifically, as shown in Figure 2, we have the data in a json file. Json files are commonly used for storing and exchanging data, so it is useful to explore how to handle json data in a Python model. For more information about the json file format, see www.json.org.

{
    "Name": "Boutique pottery shop - Model 4",
    "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}
    },
    "VarInitial": 0,
    "VarLBounds": 0,
    "VarUBounds": 100,
    "Engine": "cbc",
    "TimeLimit": 60
}

Note that to edit a json file in Jupyter Lab, you'll need to right-click on the file and select Open with > Editor.

To load the json file, we use the os and json libraries, as shown in Figure 3. This code loads all the json file data into a single object, which we'll parse in the next section.

DataFilename = os.path.join('.', 'productiondata4.json')
with open(DataFilename, 'r') as f:
    Data = json.load(f)

Declarations

The json file contains the same data as the Python file we used for Model 3. The difference is that we need to parse the json data into Python data structures before we can use it.

As shown in Figure 4, we continue to use a Pyomo concrete model. Note that we include the model's name in the declaration, so we can use Model.name in the output, even though we don't explicitly define Model.name.

Model = pyo.ConcreteModel(name = Data['Name'])

Model.Hours = pyo.Param(within = pyo.NonNegativeReals, initialize = Data['Hours'])
Model.kg = pyo.Param(within = pyo.NonNegativeReals, initialize = Data['kg'])
Model.SalesLimit = pyo.Param(within = pyo.NonNegativeReals, initialize = Data['SalesLimit'])
Model.VarInitial = pyo.Param(within = pyo.NonNegativeReals, initialize = Data['VarInitial'])
Model.VarLBounds = pyo.Param(within = pyo.NonNegativeReals, initialize = Data['VarLBounds'])
Model.VarUBounds = pyo.Param(within = pyo.NonNegativeReals, initialize = Data['VarUBounds'])
Model.Engine = pyo.Param(within = pyo.Any, initialize = Data['Engine'])
Model.TimeLimit = pyo.Param(within = pyo.NonNegativeReals, initialize = Data['TimeLimit'])

Coefficients = Data['Coefficients']
Model.Products = pyo.Set(initialize = list(Coefficients.keys()))  # Pyomo Set rather than Python set

Model.People = pyo.Param(Model.Products, within = pyo.NonNegativeReals, mutable = True)
Model.Materials = pyo.Param(Model.Products, within = pyo.NonNegativeReals, mutable = True)
Model.Sales = pyo.Param(Model.Products, within = pyo.Reals, mutable = True)
Model.Margin = pyo.Param(Model.Products, within = pyo.Reals, mutable = True)

for p in Model.Products:    
    Model.People[p] = Coefficients[p]['People']
    Model.Materials[p] = Coefficients[p]['Materials']
    Model.Sales[p] = Coefficients[p]['Sales']
    Model.Margin[p] = Coefficients[p]['Margin']

Unlike Model 3, we assign each item of data to a Pyomo pyo.Param object. We do this so that all the model's data is part of the Model object. While this isn't strictly necessary, it does enable us to create more consistent model definitions (as we'll do in the next section).

To parse the single data values, we use the field names specified in the file. For example, to get the model's name, having read the json file into Data, we use Data['Name']. We access the other single data values in a similar way.

We read the coefficients structure in one step. Then we create a products index by getting the keys of the json Coefficients structure (i.e., "Discs" and "Orbs"). We use the keys to define a Pyomo pyo.Set. Note that a Pyomo Set (uppercase S) differs from a Python set (lowercase s).

Before we populate the objects that are indexed by product – i.e., Model.People, Model.Materials, Model.Sales, and Model.Margin – we declare them as "mutable", meaning that they can be changed. We need to do this as they are initially empty, with values assigned in the block that follows.

In contrast, the other values are immutable, by default. For example, if we add the line Model.Hours = 1000 at the end of Figure 3, then we would get a runtime error telling us that Model.Hours is immutable, so its value cannot be changed. To make the value of Model.Hours changeable, we would need to add , mutable = True to its definition.

Finally, most of the parameters are defined as NonNegativeReals, meaning that they can take any value greater than or equal to zero. The exceptions are Model.Engine, which is pyo.Any because it contains text, and the values Model.Sales and Model.Margin, which we allow to be positive or negative.

Define the model

The model definition, as shown in Figure 5, is similar to Model 3. The differences are that we use the Model object throughout, rather than using a mixure of Python data structures and the Model object.

For example, in Model 3 we used Coefficients[p]['People'] * Model.Production[p], while in this model we use Model.People[p] * Model.Production[p]. In a simple model, like this one, the distinction between these two approaches is not significant. For more complex models, the code is simpler, clearer, and easier to modify, when we use the Model object consistently throughout the model definition.

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

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

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

Solve model

As shown in Figure 6, the code for solving the model is almost the same as for Model 3 – except that we use the Model object rather than individual Python data structures. Note that we need to use, for example, pyo.value(Model.Engine) to get the value, rather than just using Model.Engine.

Solver = pyo.SolverFactory(pyo.value(Model.Engine))

if pyo.value(Model.Engine) == 'cbc':
    Solver.options['seconds'] = pyo.value(Model.TimeLimit)
elif pyo.value(Model.Engine) == 'glpk':
    Solver.options['tmlim'] = pyo.value(Model.TimeLimit)

Results = Solver.solve(Model, load_solutions = False, tee = False)

Process results

The code for processing the solver result, as shown in Figure 7, is the same as for Model 3.

WriteSolution = False
Optimal = False
LimitStop = False
Condition = Results.solver.termination_condition

if Condition == pyo.TerminationCondition.optimal:
    Optimal = True
if Condition == pyo.TerminationCondition.maxTimeLimit or Condition == pyo.TerminationCondition.maxIterations:
    LimitStop = True
if Optimal or LimitStop:
    try:
        WriteSolution = True
        Model.solutions.load_from(Results)  # Defer loading results until now, in case there is no solution to load
        SolverData = Results.Problem._list
        SolutionLB = SolverData[0].lower_bound
        SolutionUB = SolverData[0].upper_bound
    except:
        WriteSolution = False

Write output

The code for writing the output, as shown in Figure 8, is almost the same as for Model 3 – again, except that we access the data values using, for example, pyo.value(Model.Engine).

print(Model.name, '\n')
print('Status:', Results.solver.termination_condition)
print('Solver:', pyo.value(Model.Engine), '\n')

if LimitStop:  # Indicate how close we are to a solution
    print('Objective bounds')
    print('----------------')
if SolutionLB is None:
    print('Lower:      None')
else:
    print(f'Lower: {SolutionLB:9,.2f}')
if SolutionUB is None:
    print('Upper:      None\n')
else:
    print(f'Upper: {SolutionUB:9,.2f}\n')
if WriteSolution:
    print(f'Total margin = ${Model.TotalMargin():,.2f}\n')
    ProductResults = pd.DataFrame()
    for p in Model.Products:
        ProductResults.loc[p, 'Production'] = round(pyo.value(Model.Production[p]), 2)
    print(ProductResults)
else:
    print('No solution loaded\n')
    print('Model:')
    Model.pprint()

When we find an optimal solution, the output is shown in Figure 9. This output is the same as for Model 3, apart from the model's name.

    Boutique pottery shop - Model 4 

    Status: optimal
    Solver: cbc 

    Total margin = $3,076.92

    Production
    Discs        6.41
    Orbs        12.82
  

Evaluation of this model

Model 4 explores an alternative external data file format, along with having a more consistent model definition.

Being able to handle different data formats is an important part of building operational models. Similarly, having a consistent model definition makes the building and maintenance of operational models easier and less error prone.

Next steps

We have one more concrete model, which we'll explore in the next article of this series.

That is, we'll look at using a more flexible way of defining constraints and the objective, using def function blocks. Using functions will give us more control over how a constraint or the objective function is defined, including making decisions about the parameters to use on a case-by-case basis, or whether we should skip defining a specific instance based on potentially complex criteria.

In addition, we'll output the slack values and dual prices (also known as shadow prices) for each constraint, to provide more information about the solution.

Conclusion

This article continues our exploration of optimization modelling in Python. Compared with Model 3, we loaded data from an external json file rather than using a Python data file, and we used Pyomo objects throughout the model, for greater consistency. 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.