Production mix - Model 11, SciPy

10 March 2023

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

SciPy is an open source Python package for a wide range of scientific computing applications, including optimization, integration, interpolation, eigenvalue problems, algebraic equations, differential equations, statistics and many other classes of problems.

Our objective is to compare a linear programming model built using SciPy with the same model built using Pyomo.

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 model

The model is available on GitHub.

Formulation for Model 11

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

Model 11 Python code

Import dependencies

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

from scipy.optimize import linprog
import pandas as pd
import numpy as np
import os.path
import json

Data file

The data for Model 11 is shown in Figure 2. We're using the json format that we used in some previous models.

As our engine we specify the HiGHS solver, which is a fast, open source linear solver that is packaged with SciPy.

{
    "Name": "Boutique pottery shop - Model 11",
    "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": "highs",
    "TimeLimit": 60
}

Get data

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

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

Declarations

As shown in Figure 4, we assign our data to suitable data structures. Unlike most other models in this series, we don't assign the data to a Model object. Instead, we'll define the objective function and constraints using lists, as expected by SciPy's linprog function. Even so, the general structure of this code is similar to previous models.

Name = Data['Name']
Hours = Data['Hours']
kg = Data['kg']
SalesLimit = Data['SalesLimit']
VarInitial = Data['VarInitial']   # Not used
VarLBounds = Data['VarLBounds']
VarUBounds = Data['VarUBounds']
Engine = Data['Engine']
TimeLimit = Data['TimeLimit']

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

Define the model

Like our CVXPY model, SciPy defines a linear program using of matrix notation. Note that the variables are implicit, rather than being explicitly declared. We will be able to access the variable values later, using the built-in name x.

As shown in Figure 5, we start by defining our parameters as empty numpy arrays (i.e., arrays of zeros, where the array length is the number of products). Then we populate the arrays using the data we loaded from the json file.

Note that SciPy always minimizes, while we want to maximize, so we need to change the sign of our objective function coefficients.

Finally, we create the objective function, constraints' left-hand side, and constraints' right-hand side lists. Note that all of our constraints are in the form ≤, which we'll specify using the A_ub and b_ub notation in the next section.

Margin = np.zeros(NumProducts)
People = np.zeros(NumProducts)
Materials = np.zeros(NumProducts)
Sales = np.zeros(NumProducts)
for p in Products:
    i = int(p)-1
    Margin[i]    = -Coefficients[p]['Margin']  # Need to negate, as SciPy always minimizes, but we want to maximize
    People[i]    =  Coefficients[p]['People']
    Materials[i] =  Coefficients[p]['Materials']
    Sales[i]     =  Coefficients[p]['Sales']

ObJCoeff = Margin
Constraints = [People, Materials, Sales]
rhs = [Hours, kg, SalesLimit];

Solve model

As shown in Figure 6, we specify and solve the model in a single step. We also include options like the solver engine and time limit.

Model = linprog(c = ObJCoeff, A_ub = Constraints, b_ub = rhs, bounds = [(VarLBounds, VarUBounds)], method = Engine, options = {'time_limit': TimeLimit})

SciPy's linprog function expects that a linear program is expressed in the form shown in Figure 7.

That is:

• Objective function coefficients (transposed) multiplied by the x variables.
• Inequality constraints expressed as upper bounds.
• Equality constraints.
• Lower and upper bounds on the variables.

Therefore, we specify the c, A_ub, b_ub, and bounds parameters in our call to the linprog function. We have no equality constraints in this model, so we omit the A_eq and b_eq parameters.

Process results

The code for processing the solver result, as shown in Figure 8, is similar to the code for previous models. linprog returns a Boolean that indicates whether the model solved successfully.

WriteSolution = False
Optimal = False
Condition = Model.success

if Condition:
    Optimal = True
    WriteSolution = True

Write output

The code for writing the output, as shown in Figure 9, is similar to our previous models. The main difference is the syntax for extracting the variable values, along with the slack and dual values, from the solver.

print(Name, '\n')
print('Status: ', Model.success, '\n')

if WriteSolution:
    print(f"Total margin = ${-Model.fun:,.2f}\n")  # Need to negate, as we're maximizing
    pd.options.display.float_format = "{:,.4f}".format
    ProductResults = pd.DataFrame()
    for p in Products:
        ProductResults.loc[p, 'Production'] = Model.x[int(p)-1]
    display(ProductResults)
    ConstraintStatus = pd.DataFrame(columns=['Slack', 'Dual'])
    for c in range(len(Constraints)):
        ConstraintStatus.loc[c] = [Model.slack[c], Model['ineqlin']['marginals'][c]]
    display(ConstraintStatus)
else:
    print('No solution loaded\n')
    print('Model:')
    print(Model)

When we find an optimal solution, the output is as shown in Figure 10. This output is similar to previous models, except for the model's name.

Note that, since SciPy always minimizes, we need to change the sign of the objective function to represent our maximization objective. Similarly, the signs of the slack and dual values are inverted (though we have printed them as reported).

    Boutique pottery shop - Model 11 

    Status:  True 

    Total margin = $3,076.92

    Production
    1       6.4103
    2      12.8205 

    Slack      Dual
    0   41.6667   -0.0000
    1    0.0000   -6.1538
    2    0.0000  -15.3846
  

Evaluation of this model

This SciPy model has more in common with our CVXPY model than with our Pyomo models. That is, both CVXPY and SciPy use a matrix notation that reflects the standard mathematical definition of a linear program. Although some mathematical purists prefer matrix notation, from a modelling perspective it can be less intuitive. Even so, despite the syntax of SciPy being markedly different to Pyomo, the general structure of the Python program is familiar.

SciPy is a capable modelling library, with a scope that is much broader than optimization. Consequently, because SciPy is not focussed solely on optimization modelling, it lacks the specific modelling support that makes Pyomo easier to use. Therefore, for a linear programming model, there is no reason to use SciPy compared with Pyomo, though SciPy may be suitable for other types of modelling.

Next steps

In the next article we'll conclude the Production mix series with a summary of the libraries we've used and our experience of working with the libraries.

Conclusion

In this article we built the Production mix model using the SciPy library. Compared with the Pyomo models, the code is broadly similar though the model specification is markedly different to Pyomo.

For some types of modelling, SciPy would be a capable and appropriate choice. Here we're solving a straightforward linear programming model, for which SciPy offers no significant advantage compared with Pyomo. Therefore, we tend to prefer Pyomo over SciPy for this type of model.

In the next article, we'll summarize and compare the optimization libraries that we've used in this series of articles.

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