Now, we will finally delve into the code of how to use gurobi with machine learning models.
The code that will be explained below is found in the following link: https://github.com/joseortegalabra/gurobi-ml-optimization-prices/blob/main/5_optimization/2_multiple_ml_models.ipynb
This is the improvement mentioned above, where multiple Machine learning models are added but it can still be run with the free gurobi license as they are all simple models (linear regressions)
Note: in the repo there are many more codes for data generation, data exploration, training of different models and development of optimization models (some simpler with a single ML model and others more complex with more complex ML models than require a license). By reading the repo readme you can delve into what other examples were developed, due to time and being outside the objective of the post, it was not mentioned or explained here.
Below are the steps to create the optimizer combining classical operations research and machine learning models (these steps are the ones followed in the codes that are in the repository, but everyone is free to order their code in the way they want. makes more sense to it).
1. Define sets
In this case, a single set is defined, the regions in which the product is sold.
It is important to remember that everything is defined based on dataframes (to use the gurobi pandas package), so sets are defined as indexes of dataframes
list_regions = ['Great_Lakes',
'Midsouth',
'Northeast',
'Northern_New_England',
'Plains',
'SouthCentral',
'Southeast',
'West']
regions = list_regionsindex_regions = pd.Index(regions)
2. Create the optimizer
Call an instance of the gurobi Model class with which you can create the optimizer which will then be solved and obtain its optimal solution
m = gp.Model(name = "Avocado_Price_Allocation")
3. Load optimizer parameters
In this example, the parameters that are loaded are the following: minimum and maximum price of the product independent of the region, the total supply of products considering all the regions (as these parameters are transversal to the sets, they are defined as a numerical python variable). On the other hand, the minimum and maximum demand in each region is also known; product supply; inventory costs and shipping costs to each region (as they are different parameters for each set, they are defined as a Series of pandas)
# a_min, a_max: min and max price of product A
a_min = 0
a_max = 2# b_min(r), b_max(r): min and max historical products send to each region (value get from historical data)
b_min = data_min_delivery
b_max = data_max_delivery
# B: supply product
B = 25
# c_waste: cost of waste product
c_waste = 0.1
# c_transport(r): cost transport for each region
c_transport = pd.Series(
{
"Great_Lakes": 0.3,
"Midsouth": 0.1,
"Northeast": 0.4,
"Northern_New_England": 0.5,
"SouthCentral": 0.3,
"Southeast": 0.2,
"West": 0.2,
"Plains": 0.2,
}, name='transport_cost')
c_transport = c_transport.loc[regions]
4. Load features from machine learning models that are not decision variables
The features of machine learning models may be decision variables of the optimizer or they may not be and may be fixed values.
Thus, its operation is more similar to a parameter (because they are fixed values) that are not explicitly defined in the optimizer formulation as they are hidden as fixed features of the machine learning model.
instance_ml_model = pd.DataFrame(
data={
"peak": peak_or_not
},
index=regions
)
This is how it is observed for each region, the fixed values of the features that are not decision variables. In the next steps this dataframe will be completed with the features of the ML model that are decision variables of the optimizer.
5. Define decision variables
Define all the optimizer decision variables. In this example, price, supply, demand prediction, inventory quantity and therefore also the quantity sold are defined as decision variables. All these decision variables are defined for each region.
Additionally, note that when defining the decision variables, the lower and upper limits that each of these variables can take were also defined; These limits may also have been defined in the form of constraints, the effect on the optimization results and speed of calculation should not be affected, this is detailed later.
# p(r): price. feature of machine learning model
p = gppd.add_vars(m, index_regions, name = "price", lb = a_min, ub = a_max) # bounds prices# x(r): supply
x = gppd.add_vars(m, index_regions, name = "x", lb = b_min, ub= b_max) # bounds supply - using historical data
# s(r): solds given a certain price
s = gppd.add_vars(m, index_regions, name = "s")
# u(r): inventary. units not sold. waste.
u = gppd.add_vars(m, index_regions, name = "w")
# d(r): demand. output of machine learning model
d = gppd.add_vars(m, index_regions, lb = -gp.GRB.INFINITY, name = "demand") # BY DEFULT LOWER BOUND IS ZERO
6. Define “classical constraints” of an optimizer
Next, the “classical constraints” of an optimizer are defined, that is, those that are not related to machine learning models.
6.1 Add the Supply Constraint
Make sure that the total number of avocados supplied is equal to B.
m.addConstr(x.sum() == B, name = 'supply')
Obs: “addConstr” writes the restrictions implicitly, that is, as if the mathematical formula were being written (defining inequalities or equalities in pandas datafames) without using gurobi functions
6.2 Add Constraints That Define Sales Quantity
The sales quantity is the minimum of the allocated quantity and the predicted demand.
gppd.add_constrs(m, s, gp.GRB.LESS_EQUAL, x, name = 'solds <= supply')
gppd.add_constrs(m, s, gp.GRB.LESS_EQUAL, d, name = 'solds <= demand')
Obs 1: “add_constrs” writes the constraints separating the left side, inequality, and the right side where the inequality must be defined using gurobi codes.
Obs 2: you can see that constraints can be saved with names. In addition, they could be saved in a python variable to be able to access said restriction easily in case you want to perform some analysis of this
6.3 Add the Wastage Constraints
Define the predicted unsold number of avocados in each region, given by the supplied quantity that is not sold. For each region r
gppd.add_constrs(m, u, gp.GRB.EQUAL, x - s, name = 'waste')
6.4 Model update — add the constraint to gurobi model
In Gurobi it is necessary to execute a model.update() to load the constraints described previously. For a data scientist this step is very similar to the model.compile() that is used when writing neural networks.
Compiling the optimization model is very important. When carrying out this step, Gurobi internally is responsible for transforming the written restrictions into an optimal structure for solving the problem. This can be observed when calling the method that solves the optimization model, which performs a first presolve step where what is described above is carried out. For example, to define the lower and upper limits of a variable, this can be done when defining the decision variable or it can be done by creating a constraint for that purpose; In the pre-solving stage, the modeling is transformed into the ideal structure to be determined by gurobi.
m.update()
7. Add constraints that are machine learning models
The restrictions associated with machine learning models are written below. Important, before writing the constraints, you must first create a dataframe that contains both the fixed features and the features that are decision variables of the optimizer. This dataframe represents the instance of the data that will be used to perform inference with the machine learning models. Since there are decision variables in the process, the optimizer will move these variables, changing the model prediction while seeking to maximize/minimize the objective function.
After having the instance defined, the machine learning model is loaded and the constraint is added. In the example code you can see that the restriction is saved in a python variable to be able to access it and print the status.
# load model
model_ml = dict_models[region]## add model to predict the demand for each region with the SAME MODEL
pred_constr = add_predictor_constr(gp_model = m,
predictor = model_ml,
input_vars = instance,
output_vars = d[region], # filter decision variable for the element of the set region,
name = f'model_predict_{region}'
)
pred_constr.print_stats()
An instance of a machine learning model can be seen as follows:
Here you can see that the region (the set defined in this optimization problem) is stored in the dataframe index. The features that are not decision variables take a fixed value while the features that are decision variables take a “gurobi var” value.
Obs 1: In this instance I have a single row since I want to use a model to predict the demand for a specific region (West region). Additionally, the option is available to use the same model (expanding the dataframe rows with the instance for the model) to predict one more set, for example the weather. Thus, to predict the demand for the product in the region r = west and for all times t using the same machine learning model, it is enough to add more rows to the dataframe instance, in this case, add more rows to represent the different times t.
Obs 2: The previous comment indicates that a dataframe with multiple rows can be used to represent a set of data that a machine learning model uses to generate an inference. On the other hand, if there are two or more sets, all that is needed is to create more input dataframes where each of them has a fixed set and only one set moves along the rows.
8. Define Objective Function
The goal is to maximize the net revenue, which is the product of price and quantity, minus costs over all regions.
m.setObjective((p * s).sum() - c_waste * u.sum() - (c_transport * x).sum(),
gp.GRB.MAXIMIZE)
9. Solve optimization problem
The objective is quadratic since we take the product of price and the predicted sales, both of which are variables. Maximizing a quadratic term is said to be non-convex, and we specify this by setting the value of the Gurobi NonConvex parameter to be 2
m.Params.NonConvex = 2
m.optimize()
Additionally, you can review whether a feasible solution was obtained or not. The gurobi code that solves the optimization problem always runs without problems (so as not to stop a python pipeline if an infeasible model is missed). Then, it is necessary to review whether a feasible solution or an infeasibility was obtained. Thus it is necessary to print the status of the solution (or read the prints obtained when executing the previous line of code).
#### know the status of the model - 2 a optimal solution was founded
# docu: https://www.gurobi.com/documentation/current/refman/optimization_status_codes.html#sec:StatusCodes
m.Status
A status equal to 2 means that a feasible solution was obtained
10. Save solution to a dataframe
In order to access the optimal values of a decision variable, it must be done through: “gppd.X”
decision_var.gppd.X
For this example:
# create dataframe with index
solution = pd.DataFrame(index=index_regions)# save optimal values
solution["Price"] = p.gppd.X
solution["Historical_Max_Price"] = data_historical_max_price # this is informative value get from historical data
solution["Allocated"] = x.gppd.X
solution["Sold"] = s.gppd.X
solution["Wasted"] = u.gppd.X
solution["Pred_demand"] = d.gppd.X
# round values
solution = solution.round(3)
# get value objetive function
opt_revenue = m.ObjVal
This is how the following dataframe is obtained with the optimal values