Production mix - Model 3, Pyomo external data

2 August 2022

In this article we continue the Python Production mix series, using the Pyomo library. Specifically, we build Model 3, which improves Model 2 by:

• Extracting the data from the model, and putting it in an external file.
• Implementing better handling of the solve process.
• Expanding the output to include additional information.

Each of these improvements are important steps in developing a practical design that can be used for operational models.

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 3

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

Model 3 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 library, we import the pandas library. We'll use pandas when writing the model's output.

import pyomo.environ as pyo
import pandas as pd

Get data

In this model we've put the data in an external Python file, as shown in Figure 2. The file contains the data we used for Model 2, plus the model's name and options for the solver.

Name  = 'Boutique pottery shop - Model 3'
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)

Engine = 'cbc'
TimeLimit = 60   # seconds

To read the data, we simply import it as we would any other Python code, as shown in Figure 3.

from productiondata3 import *

Declarations

As shown in Figure 4, we continue to use a Pyomo concrete model. The model's name comes from the external file. Since we put all the data in an external file, we need nothing further in this section.

Model = pyo.ConcreteModel(name = Name)

Define the model

The model definition, as shown in Figure 5, is identical to Model 2.

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)

Solve model

In the external data file, we specify the solver to use and a time limit option for the solver. Since different solvers may have different options, we need to translate the generic time limit data value into the chosen solver's specific option. We do this in Figure 6, where we first specify the solver to use (via the Engine object from the data file), then specify the option seconds or tmlim, representing a limit on the solve time, for solvers CBC and GLPK respectively. Note that some solvers specify a time limit in milliseconds, rather than seconds, so we may also need to translate the units to match what the solver is expecting (e.g. TimeLimit * 1000).

Note that this simply model solves in a fraction of a second, so a 60 second time limit does nothing. But for models that take longer to solve, setting a time limit can be useful.

When we call Solver.solve, we defer loading the solution by specifying load_solutions = False. We do to avoid the program crashing if there is no solution to load.

We've also added the option tee, which displays the solver's output as it runs. This can be useful for seeing how the solver is progressing towards a solution, or for diagnosing problems with the solve process. Try changing the option's value to True to see what it does.

Solver = pyo.SolverFactory(Engine)

if Engine == 'cbc':
    Solver.options['seconds'] = TimeLimit
elif Engine == 'glpk':
    Solver.options['tmlim'] = TimeLimit

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

Process results

Before writing any output, we need to process the result of the solve to see what happened, as shown in Figure 7.

In this case, we want to decide if there is a solution to write, and the condition in which the solver stopped – so we create three Boolean variables to help us make decisions about what to output.

We get the solver's condition from its termination_condition property, from which we can decide whether the result is optimal or whether the solver stopped on either a time limit or the maximum number of iterations. There are other possible stopping conditions, such as the solver encountering an error – in those conditions, we usually won't have a solution to write.

Note that various solvers are not consistent in how they report hitting a time limit or maximum number of iterations. Here we make no distinction, though sometimes it might matter which condition occurs.

If the solver stopped on a time or iterations limit, it may not have found a solution. Therefore, we try to load a solution. If there is no solution to load, then this attempt will fail, so we run the except clause, where we indicate that there is no solution to write. If the test passes, then we have a solution to write. We also get the lower and upper bounds on the objective function – as these can be helpful for indicating how close a feasible solution is to being optimal.

If the solver stopped with an optimal solution, then it should return a solution for us to output – though we allow for the possibility that this case also fails.

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

Given information about the outcome of the solve process, we can now decide what output to write, as shown in Figure 8:

• Irrespective of the solver's condition, we output the model's name, the solve status, and the solver we used.
• If the solver stopped on a time or iterations limit, then we can write some potentially useful information about the bounds on the objective function value.
• If there is a solution to write, then we write the same output that we wrote for Model 2.
• If there is no solution available to write, then we say so. We also output the model, as it might be helpful in diagnosing whatever issues prevented the solver from finding a solution (setting the tee option to True, mentioned above, may also be useful).

print(Model.name, '\n')
print('Status:', Results.solver.termination_condition)
print('Solver:', 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 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 2, apart from the model's name.

    Boutique pottery shop - Model 3 

    Status: optimal
    Solver: cbc 

    Total margin = $3,076.92

    Production
    Discs        6.41
    Orbs        12.82
  

We can make this model infeasible by, for example, setting the lower bound on the variables to 30. We do that by changing the data file to say VarBounds = (30, 100) (and saving the file). If we run the model again, then the output is as shown in Figure 10, where we state that the model is infeasible, that no solution has been loaded, and we output the model (including the bounds of 30 and 100 on the variables).

    Boutique pottery shop - Model 3 

    Status: infeasible
    Solver: cbc 

    No solution loaded

    Model:
    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 :    30 :     0 :   100 : False : False : NonNegativeReals
    Orbs :    30 :     0 :   100 : 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
    MaterialUsage : Size=1, Index=None, Active=True
    Key  : Lower : Body                                           : Upper : Active
    None :  -Inf : 18.0*Production[Discs] + 30.0*Production[Orbs] : 500.0 :   True
    PeopleHours : Size=1, Index=None, Active=True
    Key  : Lower : Body                                           : Upper : Active
    None :  -Inf : 12.5*Production[Discs] + 10.0*Production[Orbs] : 250.0 :   True
    SalesRelationship : 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 PeopleHours MaterialUsage SalesRelationship TotalMargin
  

Remember to change the bounds back to 0 and 100 in the data file.

Evaluation of this model

Model 3 is a significant improvement compared with Model 2.

Having an external data file allows us to easily change the data independently of the model. That is, we can either overwrite the data file with new data, or we can have several data files from which we choose one to solve.

Having attempted to solve the model, this version gives us more control over the output. There may or may not be a solution available to write out, so we need to handle either outcome. We can also choose what information to write, depending on the solver's condition.

Next steps

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

In the next article, we'll put the data in an external json file, which we may receive if the data is created by some automated process. We'll also read the data into the Model object, rather than separate objects, and we'll make the model definition more consistent.

Conclusion

This article continues our exploration of optimization modelling in Python. We've improved Model 2, but there's still more to do. 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.