12 December 2022 (1,449 words)

CVXPY

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

CVXPY is a Python package for representing and solving convex optimization problems.

Our objective is to compare a model built using CVXPY 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 10

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 10 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 CVXPY library, which we've previously installed, along with some other libraries that we'll use.

Figure 2. Import dependencies
import cvxpy as cp
import pandas as pd
import numpy as np
import os.path
import json

Data file

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

As our engine we specify the cvxopt solver, which is a general purpose convex solver suitable for linear problems. Other solvers can also be used. Our code allows for the use of cbc, glop, glpk, or cvxopt – all of which we have previously installed. Each of these solvers returns the same solution, except that cbc fails to return dual values. When using cbc with Pyomo, we needed to specifically request dual values before running the model (see Production mix - Model 5), but we have not been able to find a similar process with CVXPY.

Figure 3. External data file, productiondata9.json
{
    "Name": "Boutique pottery shop - Model 10",
    "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": "cvxopt",
    "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('.', 'productiondata10.json')
with open(DataFilename, 'r') as f:
    Data = json.load(f)

Declarations

As shown in Figure 5, we assign various data to suitable data structures. Unlike most other models in this series, we don't assign the data to a Model object. This is because when declaring the CVXPY model we must include the objective function and constraints – which we'll do later. Even so, the general structure of this code is very similar to previous models.

Figure 5. Declarations
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)

Margin = np.zeros(NumProducts)
People = np.zeros(NumProducts)
Materials = np.zeros(NumProducts)
Sales = np.zeros(NumProducts)
for p in Products:
    i = int(p)
    Margin[i]    = Coefficients[p]['Margin']
    People[i]    = Coefficients[p]['People']
    Materials[i] = Coefficients[p]['Materials']
    Sales[i]     = Coefficients[p]['Sales']

Define the model

The model definition, as shown in Figure 6, starts with defining the variables and the objective function.

The constraints are defined as a list. First we create an empty list, then append each constraint to the list. The variable bounds are explicitly specified as constraints.

Note that CVXPY uses the @ symbol to represent matrix multiplication. Unlike most other optimization libraries, we don't need to specify the domain over which we're summing.

We've adopted a simple style in this example, though we could use the def structure that we've used in some previous models. When developing more complex models, it can be useful to adopt a def structure, to make the model more flexible and easier to understand.

Figure 6. Define the model
Production = cp.Variable(NumProducts)   # Variables

objective = cp.Maximize(cp.sum(Margin @ Production))   # Objectve function
constraints = []   # Constraints
constraints += [cp.sum(People @ Production) <= Hours]
constraints += [cp.sum(Materials @ Production) <= kg]
constraints += [cp.sum(Sales @ Production) <= SalesLimit]

constraints += [Production >= VarLBounds]   # Bounds on variables
constraints += [Production <= VarUBounds]

Solve model

As shown in Figure 7, we convert our string name for the solver into an engine object, as expected by CVXPY. We also define options that specify how verbose the solver output is and a maximum time limit.

Note that we can obtain a list of installed solvers that CVXPY can use via the code: print(cp.installed_solvers()).

Figure 7. Solve model
Model = cp.Problem(objective, constraints)

if Engine == 'cbc':
    EngineObj = cp.CBC
elif Engine == 'glop':
    EngineObj = cp.GLOP
elif Engine == 'glpk':
    EngineObj = cp.GLPK
elif Engine == 'cvxopt':
    EngineObj = cp.CVXOPT

Result = Model.solve(solver=EngineObj, verbose=True, max_seconds=TimeLimit)

Process results

The code for processing the solver result, as shown in Figure 8, is similar to the code for previous models.

Figure 8. Process results
WriteSolution = False
OptimalWriteSolution = False
Optimal = False
Condition = Model.status

if Condition == 'optimal':
    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 slack and dual values from the solver output files.

Figure 9. Write output
print(Name, '\n')
print('Status:', Model.status)
print('Solver:', Engine, '\n')

if WriteSolution:
    print(f"Total margin = ${objective.value:,.2f}\n")
    pd.options.display.float_format = "{:,.4f}".format
    ProductResults = pd.DataFrame()
    for p in Products:
        ProductResults.loc[p, 'Production'] = Production[int(p)].value
    display(ProductResults)
    
    ConstraintStatus = pd.DataFrame(columns=['Slack', 'Dual'])
    for c in range(3):
        ConstraintStatus.loc[c] = [constraints[c].expr.value, constraints[c].dual_value]
    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's name.

Figure 10. Solution
Boutique pottery shop - Model 10 

Status: optimal
Solver: cvxopt 

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

This CVXPY model is similar to our Pyomo models. Although the syntax of the two libraries is markedly different, the general structure of the model definitions and solution process is familiar.

CVXPY is a capable modelling library. It is focussed on general convex optimization, offering access to some solvers that are not supported by Pyomo – like CVXOPT (a package for solving convex optimization models) and OSQP (a package for solving convex quadratic models) – which may be important in some situations. For our linear programming model, there is no reason to use CPXPY compared with Pyomo.

Next steps

In the next article we'll repeat the model implementation using our last selected library SciPy, which is a comprehensive collection of tools and algorithms for a wide range of scientific computing tasks.

Conclusion

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

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

In the next article, we'll build the Production mix model using SciPy.

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