Solver Max - Production mix - Model 6, Pyomo abstract

21 September 2022

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

Declare the model as a Pyomo pyo.AbstractModel, rather than as a pyo.ConcreteModel.
Read the data from a dat file rather than a json file.

These changes show that, contrary to how abstract and concrete models are portrayed in most blogs, there is actually little difference between abstract and concrete Pyomo models.

Articles in this series

Articles in the Python Production mix series:

Download the models

The model is available on GitHub.

Formulation for Model 6

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

Model 6 Python code

Import dependencies

The first task is to import the libraries that are needed for our program. As shown in Figure 1, we aren't using a json file, so we don't need the os and json libraries that we imported in Model 5.

Figure 1. Import dependencies

import pyomo.environ as pyo
import pandas as pd

Data file

The data for Model 6 is shown in Figure 2. Unlike Model 5, where we read json format data and assign it to Pyomo Param data structures, here we define the data directly in param data structures. Note the different capitalization of the two structures – uppercase for parameters defined in the model, and lowercase for parameters defined in a dat file. The same convention applies to sets.

Figure 2. External data file, productiondata6.dat

param: Name  := 'Boutique pottery shop - Model 6';
param: Hours := 250;
param: kg    := 500;
param: SalesLimit := 0;
param: VarInitial := 0;
param: VarLBounds := 0;
param: VarUBounds := 100;
param: Engine := "cbc";
param: TimeLimit := 60;

param:  Products: People  Materials  Sales  Margin  :=
    Discs     12.50   18.00      -2      80.00
    Orbs      10.00   30.00       1     200.00;

Note that, unlike concrete models, we don't load the data at this stage. This is the essential difference between the two types of Pyomo model. That is:

pyo.ConcreteModel. Data is loaded before the constraints and objective function are defined.
pyo.AbstractModel. The constraints are objective function are defined first, with the data loaded only when we solve the model.

Declarations

We declare an pyo.AbstractModel, rather than a pyo.ConcreteModel. Because we have not loaded the data, the declarations are empty sets and parameters, as shown in Figure 3. Though we can declare the domains, such as pyo.NonNegativeReals.

Figure 3. Declarations

Model = pyo.AbstractModel()

Model.Products = pyo.Set()  # Pyomo Set rather than Python set

Model.Name = pyo.Param(within = pyo.Any)
Model.Hours = pyo.Param(within = pyo.NonNegativeReals)
Model.kg = pyo.Param(within = pyo.NonNegativeReals)
Model.SalesLimit = pyo.Param(within = pyo.NonNegativeReals)
Model.VarInitial = pyo.Param(within = pyo.NonNegativeReals)
Model.VarLBounds = pyo.Param(within = pyo.NonNegativeReals)
Model.VarUBounds = pyo.Param(within = pyo.NonNegativeReals)
Model.Engine = pyo.Param(within = pyo.Any)
Model.TimeLimit = pyo.Param(within = pyo.NonNegativeReals)
Model.People = pyo.Param(Model.Products, within = pyo.NonNegativeReals) 
Model.Materials = pyo.Param(Model.Products, within = pyo.NonNegativeReals)
Model.Sales = pyo.Param(Model.Products, within = pyo.Reals)
Model.Margin = pyo.Param(Model.Products, within = pyo.Reals)

Define the model

The model definition, as shown in Figure 4, is the same as for Model 5.

In a more complex model, the concrete and abstract model definitions may not be the same. Specifically, in a concrete mode, we may use the data to make decisions within the constraint and objective function rules. That is not possible in a pure abstract model – though we could read some data, to inform decision making, but there is typically little advantage in adopting a hybrid approach.

Figure 4. Define the model

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 5, we create an instance of an abstract model by loading the data just before we solve the model. Having created an instance, we need to refer to our Instance object rather than our Model object. In abstract models, a common error is to refer to the model definition rather than its instance.

Figure 5. Solve model

Instance = Model.create_instance("productiondata6.dat")
Solver = pyo.SolverFactory(pyo.value(Instance.Engine))
if pyo.value(Instance.Engine) == 'cbc':
    Solver.options['seconds'] = pyo.value(Instance.TimeLimit)
elif pyo.value(Instance.Engine) == 'glpk':
    Solver.options['tmlim'] = pyo.value(Instance.TimeLimit)

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

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

Process results

The code for processing the solver result, as shown in Figure 6, is the same as for Model 5 except that we refer to Instance rather than Model.

Figure 6. Process results

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
        Instance.solutions.load_from(Results)
        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 7, is almost the same as for Model 5, except that, again, we refer to Instance rather than Model.

Figure 7. Write output

print(pyo.value(Instance.Name), '\n')
print('Status:', Results.solver.termination_condition)
print('Solver:', pyo.value(Instance.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 = ${Instance.TotalMargin():,.2f}\n')
    pd.options.display.float_format = "{:,.4f}".format
    ProductResults = pd.DataFrame()
    for p in Instance.Products:
        ProductResults.loc[p, 'Production'] = round(pyo.value(Instance.Production[p]), 4)
    print(ProductResults, '\n')

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

When we find an optimal solution, the output is shown in Figure 8. This output is the same as for Model 5.

Figure 8. Solution

    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

Evaluation of this model

Model 6 is our final Pyomo model in this series of articles. This model translates the concrete Model 5 into an abstract model.

Most blog articles about Pyomo models suggest that there is a substantial difference between concrete and abstract models. A model like our Model 1 is often used as a typical concrete model, while a model like our Model 6 is used as a typical abstract model. When presented that way, there is indeed a substantial difference between the two types of model.

However, a fairer comparison is between our Model 5 and Model 6. When presented this way, there is actually little difference between the two model types. As the Pyomo documentation says:

Python programmers will probably prefer to write concrete models, while users of some other algebraic modeling languages may tend to prefer to write abstract models. The choice is largely a matter of taste; some applications may be a little more straightforward using one or the other.
Abstract versus Concrete models

Next steps

So far in this series, we've explored the Pyomo library quite extensively. Next, we'll continue our Production mix series by implementing the model using other Python modelling libraries. Specifically, we'll use the PuLP, OR Tools, Gekko, CVXPY, and SciPy libraries to implement the same model. By using different libraries, we'll get some insight into their similarities and differences, along with their relative strengths and weaknesses.

Conclusion

This article completes the Pyomo part of this series. Compared with models 1 to 5, which use Pyomo's concrete model type, here we use Pyomo's abstract model type. It turns out that there isn't much difference between the two model types, with the choice between them being largely a matter of taste.

In the next article, we'll start using different libraries to build the same Production mix model.

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 6, Pyomo abstract