Production mix - Model 5, Pyomo using def

5 September 2022

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

• Define the constraints and objective function using def function blocks.
• Output the slack values and dual prices (also known as shadow prices) for each constraint.

These changes give us more control over how the model is defined and provide more information about the solution.

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 5

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

Model 5 Python code

Import dependencies

The first task is to import the libraries that are needed for our program. The dependencies, as shown in Figure 1, are the same as for Model 4.

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

Get data

The data for Model 5, as shown in Figure 2, is identical to the data for Model 4, except for the model's name.

{
  "Name": "Boutique pottery shop - Model 5",
  "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('.', 'productiondata5.json')
with open(DataFilename, 'r') as f:
    Data = json.load(f)

Declarations

The declarations, as shown in Figure 4, are identical to Model 4.

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']

Define the model

The model definition is shown in Figure 5. Our definition of the Production variables is the same as for Model 4, but the definitions of the constraints and the objective function are very different.

Specifically, we use a def function to define each of the constraints and the objective function. In this model, each definition consists of three lines of code:

• First line: Declare a function. We need to give the function a name and pass at least the Model object to the function. In more complex models we would also pass indices and possibly other objects.
• Second line: We return a rule for the constraint or objective function. The rules are the same as the expressions we defined in Model 4. In more complex models, this line would expand to multiple lines, with logic to decide what to return. For example, a common task is to skip specific instances of a constraint.
• Third line: Here we call the function, indicating which rule to use and whether we're defining a pyo.Constraint or pyo.Objective. This line doesn't need to be immediately after the function definition, though it is common practice to do so. For complex models, which have long functions, it may be clearer to define the functions and then have their calls collated in a subsequent block of code.

The advantage of using def functions is that they give us much greater control over the definitions. For all but simple models, the standard way to define constraints and the objective function in Pyomo is to use def functions.

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

def rule_hours(Model):
    return sum(Model.People[p] * Model.Production[p] for p in Model.Products) <= Model.Hours
Model.PeopleHours = pyo.Constraint(rule = rule_hours)

def rule_usage(Model):
    return sum(Model.Materials[p] * Model.Production[p] for p in Model.Products) <= Model.kg
Model.MaterialUsage = pyo.Constraint(rule = rule_usage)

def rule_sales(Model):
    return sum(Model.Sales[p] * Model.Production[p] for p in Model.Products) <= Model.SalesLimit
Model.SalesRelationship = pyo.Constraint(rule = rule_sales)

def rule_Obj(Model):
    return sum(Model.Margin[p] * Model.Production[p] for p in Model.Products)
Model.TotalMargin = pyo.Objective(rule = rule_Obj, sense = pyo.maximize)

Solve model

As shown in Figure 6, the code for solving the model is almost the same as for Model 4. The exception is that we've added a line that tells the solver to collate Model.dual prices for the constraints, which we'll print in the output.

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)

Model.dual = pyo.Suffix(direction = pyo.Suffix.IMPORT)

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

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 4, except that:

• We define the default data frame format to have 4 decimal places (i.e., using pd.options.display.float_format), so we have additional precision for calculations that we do below.
• We've added a section to populate a ConstraintStatus data frame that contains the slack values and dual prices for each constraint.

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')
    pd.options.display.float_format = "{:,.4f}".format
    ProductResults = pd.DataFrame()
    for p in Model.Products:
        ProductResults.loc[p, 'Production'] = pyo.value(Model.Production[p])
        print(ProductResults, '\n')

    ConstraintStatus = pd.DataFrame(columns=['lSlack', 'uSlack', 'Dual'])
    for c in Model.component_objects(pyo.Constraint, active = True):
        ConstraintStatus.loc[c.name] = [c.lslack(), c.uslack(), Model.dual[c]]
        print(ConstraintStatus)
else:
    print('No solution loaded\n')
    print('Model:')
    Model.pprint()

When we find an optimal solution, the output is shown in Figure 9.

    Boutique pottery shop - Model 5 

    Status: optimal
    Solver: cbc 

    Total margin = $3,076.92

    Production
    Discs      6.4103
    Orbs      12.8205 

    lSlack  uSlack    Dual
    PeopleHours           inf 41.6667 -0.0000
    MaterialUsage         inf -0.0000  6.1538
    SalesRelationship     inf -0.0000 15.3846
  

The table of slack values and dual prices can be interpreted as follows:

• lSlack is the difference between the constraint's lower bound and its value in the solution. None of our constarints have lower bounds, so the differences are infinite.
• uSlack is the difference between the constraint's upper bound and its value in the solution. For the PeopleHours constraint, uSlack = 41.6667. We can verify this value by substituting the solution's variable values into the constraint, so the left-hand side is: 12.50*6.4103 + 10.00*12.8205 = 208.3338, which is 41.6662 less than the right-hand side of 250 (give-or-take a small rounding difference). If we do a similar calculation for the MaterialsUsage and SalesRelationship constraints, the result is zero – indicating that those constraints are binding.
• The dual prices indicate the marginal change in the objective function for a unit change in a constraint's right-hand side, all else being equal. For the non-binding PeopleHours constraint, the dual price is zero because a change in the right-hand side value has no impact on the solution's objective function value. The other two constraints are binding, so their dual prices are non-zero. For example, if we change the available materials to 501 kg, and re-solve the model, then the objective function increases by 6.16 to 3,083.08 (again, allowing for a small rounding difference). Similarly, changing the available materials to 499 kg reduces the objective function value by 6.15 to 3,070.77.

Evaluation of this model

Model 5 is our final concrete model in this series of articles.

This model contains the essential elements of a Pyomo optimization model, including importing data, creating data structures, defining the model, solving the model, processing the result, and producing outputs that display the solution and help us understand the solution.

Every model is different, so other models may need to alter or extend the design of this model. Nonetheless, this model provides a good template from which to develop other Pyomo models.

Next steps

We could further extend this model, particularly by adding more error checking and data validation. For an operational model, such features may be important.

But we'll leave such details for another time. Instead, the next article will focus on an alternative implementation, using a Pyomo "abstract" model. Subsequent articles will then implement the Production mix model using other Python modelling libraries.

Conclusion

This article continues our exploration of optimization modelling in Python. Compared with Model 4, we define the constraints and objective function using functions, which gives us much greater control over the model definition. We also output additional information about the solution, specifically the slack values and dual prices.

In the next article, we'll look at an alternative way of defining the model, via a Pyomo abstract model.

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