22 November 2022 (1,512 words)

Gekko

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

GEKKO is a Python package for machine learning, optimization of mixed-integer, and differential algebraic equations.

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

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

Figure 2. Import dependencies
from gekko import GEKKO
import pandas as pd
import os.path
import json
import numpy as np

Data file

The data for Model 9 is shown in Figure 3. We're using the json format that we used in some previous models. We specify the apopt solver, which is built-in to Gekko and suitable for solving linear and mixed integer problems. Other built-in solvers include bpopt (best for systems biology applications) and ipopt (a non-linear solver generally best for problems with large numbers of degrees of freedom or when starting without a good initial guess).

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

Declarations

As shown in Figure 5, we create the model and assign various data to the Model. We maintain a code style that is similar to our previous models in this series.

Note that the model declaration includes a remote parameter. Gekko models can be solved either locally or on a remote server. For our model, specifying the remote server leads to successfully solving the model, but attempting to retrieve the dual values produces an error because our code expects to find the values in the local temporary folder. We haven't found a way to retrieve the dual values when using the remote server option.

Figure 5. Declarations
Model = GEKKO(name=Data['Name'], remote=False)

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

Coefficients = Data['Coefficients']
Model.Products = 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, starts with defining the variables. We explicitly define the objective function variable, unlike most other optimization libraries. We then define the production variables using a dictionary. While this approach is not required by Gekko, it is convenient.

The constraint and objective function definitions are very similar to previous models. Note that we specify the optimization direction using Model.Maximize on the objective function variable.

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
Model.TotalMargin = Model.Var()
Model.Production  = dict(map(lambda p: (p, Model.Var(lb=Model.VarLBounds, ub=Model.VarUBounds)), Model.Products))

Model.PeopleHours = Model.Equation(sum(Model.People[p] * Model.Production[p] for p in Model.Products) ≤ Model.Hours)
Model.MaterialUsage = Model.Equation(sum(Model.Materials[p] * Model.Production[p] for p in Model.Products) ≤ Model.kg)
Model.SalesRelationship = Model.Equation(sum(Model.Sales[p] * Model.Production[p] for p in Model.Products) ≤ Model.SalesLimit)

Model.Equation(Model.TotalMargin == sum(Model.Margin[p] * Model.Production[p] for p in Model.Products))
Model.Maximize(Model.TotalMargin)

Solve model

As shown in Figure 7, we convert our string name for the solver into an engine number, as expected by Gekko. We also define some options, including for the level of output that the solver produces during the solve process. We want to display dual prices in the results, so we need to specify a diagnostic level sufficient to enable us to do that. The diagnostic level can be varied from 0 to 5.

Figure 7. Solve model
Model.options.MAX_TIME = Model.TimeLimit
Model.options.DIAGLEVEL = 2      # Enable extraction of dual prices

if Model.Engine == 'apopt':
    EngineNum = 1
elif Model.Engine == 'bpopt':
    EngineNum = 2
elif Model.Engine == 'ipopt':  # ipopt will be used if other solvers are not available
    EngineNum = 3

try:
    Success = True
    Model.solve(solver=EngineNum, linear=1, disp=True, debug=True)
except:
    Success = False

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
Optimal = False

if Success:
    Optimal = True
    WriteSolution = True
    StatusText = 'Optimal'
else:
    StatusText = 'Unsuccessful'

Write output

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

Figure 9. Write output
print(Data['Name'],'\n')
print('Status:', StatusText)
print('Solver:', Model.Engine, '\n')

if WriteSolution:
    print(f"Total margin = ${Model.TotalMargin.value[0]:,.2f}\n")
    pd.options.display.float_format = "{:,.4f}".format
    ProductResults = pd.DataFrame()
    for p in Model.Products:
        ProductResults.loc[p, 'Production'] = Model.Production[p].value[0]
    display(ProductResults)
    
    ConstraintStatus = pd.DataFrame(columns=['Slack', 'Dual'])
    Duals = np.loadtxt(Model.path+'/apm_lam.txt')                       # Read dual prices from temporary folder
    ResultFilename = os.path.join(Model.path, 'results.json')
    with open(ResultFilename, 'r') as f:
        Results = json.load(f)                                          # Read slack values from temporary folder    
    for c in range(3):
        ConstraintStatus.loc[c] = [Results['slk_'+str(c+1)][0], Duals[c]]
    display(ConstraintStatus)    
else:
    print('No solution loaded\n')

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. Also note that the dual prices having the opposite sign compared with previous models.

Figure 10. Solution
Boutique pottery shop - Model 9 

Status: Optimal
Solver: apopt 

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.3845

Finally, unlike most other libraries, we need to clean up after the solve is complete, as shown in Figure 11. The clean up process removes the solver's temporary folder.

Figure 11. Clean up
Model.cleanup()

Evaluation of this model

This model is very similar to 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.

Despite the similarities, we tend to prefer Pyomo for solving linear programs, simply because Pyomo is focussed on optimization modelling, while Gekko is focussed more on other applications.

Next steps

In subsequent articles we'll repeat the model implementation using our other selected libraries. Next on the list is CVXPY, which is a modelling language for convex optimization problems.

Conclusion

In this article we built the Production mix model using the Gekko library. Compared with the Pyomo models, the code is quite similar. Gekko is a capable modelling library that is easy to use. However, it is focussed more on machine learning and related applications, rather than optimization. Since Gekko offers no significant advantage compared with Pyomo, we tend to prefer Pyomo over Gekko for this type of model.

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

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