4 October 2022 (1,325 words)

PuLP

In this article we continue the Python Production mix series. Specifically, we build Model 7 using the PuLP library.

PuLP, like Pyomo, is a COIN-OR project. The library is freely available, the code is open source, and it is widely used.

Our objective is to compare a model built using PuLP with the same model built using Pyomo.

Articles in this series

Articles in the Python Production mix series:

Download the model

The Python code and data for this model are available in the following files:

The "Jupyter notebook" file contains a formatted combination of Python code and markdown code – this file should be opened and run in Jupyter Lab. We describe setting up our Python environment, including Jupyter Lab and various Python libraries, in the article Set up a Python modelling environment.

The "Python code" file is a plain text file containing only the Python code of this model. Download this file if you have a non-Jupyter environment for running Python programs.

The "Data" file is a plain text file containing the model's data, in json format.

The model files are also available on GitHub.

Formulation for Model 7

For this model, we're using the same general formulation that we used for previous models in this series, as shown in Figure 1.

Figure 1. General formulation
\begin{alignat}{1} &\text{Objective} \\ & \quad \text{Maximize} \ fTotalMargin &= \quad & \sum_{p=1}^n dMargin_{p} \times vProduction_{p} \tag{1} \\ &\text{Subject to} \\ & \quad \sum_{p=1}^n dPeople_{p} \times vProduction_{p} &\le & dHours \tag{2}\\ & \quad \sum_{p=1}^n dMaterials_{p} \times vProduction_{p} \quad &\le & dkg \tag{3}\\ & \quad \sum_{p=1}^n dSales_{p} \times vProduction_{n} &\le & 0 \tag{4}\\ &\text{Variables} \\ & \quad vProduction_{p} &= & \text{Units of each product} \tag{5} &\begin{array}{l} \forall \ p \in \{1 \ldots n\} \end{array} \\ &\text{Data} && \\ & \quad dMargin_{p} &= &\text{Profit margin on each unit} \tag{6} &\begin{array}{l} \forall \ p \in \{1 \ldots n\} \\ \end{array} \\ & \quad dPeople_{p} &= &\text{Staff hours to make each unit} \tag{7} &\begin{array}{l} \forall \ p \in \{1 \ldots n\} \\ \end{array} \\ & \quad dMaterials_{p} \quad &= &\text{kg of materials to make each unit} \tag{8} &\begin{array}{l} \forall \ p \in \{1 \ldots n\} \\ \end{array} \\ & \quad dSales_{p} &= &\text{Relative sales volume for each product} \quad \tag{9} &\begin{array}{l} \forall \ p \in \{1 \ldots n\} \\ \end{array} \\ &\text{Indices} \\ & \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 7 Python code

Import dependencies

The first task is to import the libraries that are needed for our program. As shown in Figure 2, we import the pulp library, which we've previously installed, along with some other libraries that we'll use.

Figure 2. Import dependencies
import pulp as pu
import pandas as pd
import os.path
import json

Data file

The data for Model 7 is shown in Figure 3. We're using the json format that we used in some previous models. The only difference is that PuLP does not allow spaces in the model name, so we use underscores instead.

Figure 3. External data file, productiondata7.json
{
    "Name": "Boutique_pottery_shop_Model_7",
    "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": 10
}

Get data

We import the data from the json file using the code shown in Figure 4. This code is the same as the code we used for previous json files, apart from the filename.

Figure 4. Get data
DataFilename = os.path.join('.', 'productiondata7.json')
with open(DataFilename, 'r') as f:
    Data = json.load(f)

Declarations

As shown in Figure 5, we declare the model as a linear program, and specify that it is a maximization problem. The data is assigned to a Model object, in a way that is similar to previous models. The syntax is simpler than we used in, for example, our Pyomo Model 5, though Pyomo offers a richer object model.

Figure 5. Declarations
Model = pu.LpProblem(Data['Name'], pu.LpMaximize)

Model.Hours = Data['Hours']
Model.kg = Data['kg']
Model.SalesLimit = Data['SalesLimit']
Model.VarInitial = Data['VarInitial']
Model.VarLBounds = Data['VarLBounds']
Model.VarUBounds = Data['VarUBounds']
Model.Engine = Data['Engine']
Model.TimeLimit = Data['TimeLimit']

Coefficients = Data['Coefficients']
Model.Products = list(Coefficients.keys())

Model.People = {}
Model.Materials = {}
Model.Sales = {}
Model.Margin = {}

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, as shown in Figure 6, is similar to how we defined the Pyomo Model 5 and Model 6. That is, we use def functions to define each constraint and the objective.

Like with our Pyomo models, it is not necessary to use def functions in this way. However, in more complex models, this approach gives us more control over the definitions – espeically when we need to make decisions about what terms in include in a constraint or the objective function.

Note that at the end of each constraint we provide a name. This name is used in the model output.

Figure 6. Define the model
Model.Production = pu.LpVariable.dicts("Products", Model.Products, lowBound=Model.VarLBounds, upBound=Model.VarUBounds, cat=pu.LpContinuous)
for p in Model.Products:
    Model.Production[p].setInitialValue(Model.VarInitial)

def constraint_hours():
    return pu.lpSum([Model.People[p] * Model.Production[p] for p in Model.Products]) <= Model.Hours, 'PeopleHours'
Model += constraint_hours()

def constraint_usage():
    return pu.lpSum([Model.Materials[p] * Model.Production[p] for p in Model.Products]) <= Model.kg, 'MaterialUsage'
Model += constraint_usage()

def constraint_sales():
    return pu.lpSum([Model.Sales[p] * Model.Production[p] for p in Model.Products]) <= Model.SalesLimit, 'SalesRelationship'
Model += constraint_sales()

def objective_margin():
    return pu.lpSum([Model.Margin[p] * Model.Production[p] for p in Model.Products])
Model.setObjective(objective_margin())

Solve model

As shown in Figure 7, we define options specific to a solver – in this case, a time limt (in seconds) for either CBC or GLPK. We then solve the model and record the solve status.

Figure 7. Solve model
if Model.Engine == 'cbc':
    Solver = pu.PULP_CBC_CMD(timeLimit = Model.TimeLimit)
elif Model.Engine == 'glpk':
    Solver = pu.GLPK_CMD(timeLimit = Model.TimeLimit)

Status = Model.solve(Solver)

Process results

The code for processing the solver result, as shown in Figure 8, is similar to the code for Model 6 except that we've simpiflied it to only write a solution if the solver status is optimal.

Note that the Status we recorded above is an integer. We could use that value to process the results. Alternatively, as we do below, we can ask PuLP to provide a text status, such as "Optimal".

Figure 8. Process results
WriteSolution = False
Optimal = False
Condition = pu.LpStatus[Model.status]

if Condition == 'Optimal':
    Optimal = True
    WriteSolution = True

Write output

The code for writing the output, as shown in Figure 9, is very similar to our previous models.

Figure 9. Write output
print(Model.name, '\n')
print('Status:', pu.LpStatus[Model.status])
print('Solver:', Model.Engine, '\n')

if WriteSolution:
    print(f"Total margin = ${Model.objective.value():,.2f}\n")
    pd.options.display.float_format = "{:,.4f}".format
    ProductResults = pd.DataFrame()
    for p in Model.Products:
        ProductResults.loc[p, 'Production'] = pu.value(Model.Production[p])
    display(ProductResults)

    ConstraintStatus = pd.DataFrame(columns=['Slack', 'Dual'])
    for name, c in list(Model.constraints.items()):
        ConstraintStatus.loc[name] = [c.slack, c.pi]
    display(ConstraintStatus)
else:
    print('No solution loaded\n')
    print('Model:')
    print(Model)

When we find an optimal solution, the output is shown in Figure 10. This output is similar to previous models, except for the model name and the slack values are defined in only one direction.

Figure 10. Solution
Boutique_pottery_shop_Model_7 

Status: Optimal
Solver: cbc 

Total margin = $3,076.92

       Production
Discs      6.4103
Orbs      12.8205 

                      Slack    Dual
PeopleHours         41.6667 -0.0000
MaterialUsage       -0.0000  6.1538
SalesRelationship   -0.0000 15.3846

Evaluation of this model

There is a close similarity between this PuLP model and our Pyomo models. Although the syntax of the two libraries is somewhat different, the general structure of the model definitions and solution process is familiar. This isn't surprising, as both PuLP and Pyomo are COIN-OR projects.

Despite the similaries, we tend to prefer Pyomo simply because PuLP offers no significant advantages, while Pyomo provides easier access to a wider range of solvers – enabling us to solve a greater variety of model types.

Next steps

In subsequent articles we'll repeat the model implementation using our other selected libraries. Next on the list is OR-Tools, which takes quite a different approach to this type of modelling.

Conclusion

In this article we built the Production mix model using the PuLP library. Compared with the Pyomo models, the code is quite similar. PuLP is a capable modelling library that is easy to use. However, since it offers no significant advantage compared with Pyomo, we tend to prefer Pyomo over PuLP.

In the next article, we'll build the Production mix model using OR-Tools.

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