Rolling horizon inventory problem
AnsweredHello,
I am trying to model an inventory optimization problem in a rolling horizon manner (A.K.A model predictive control)
This is my model with description:
Note: each item of stock has 14 days shelf life and after that, it is disposed incurring a wastage.
Here is my attempted code in Gurobipy:
And this the error I get:
I was wondering if someone could take a look at my code and tell me what is wrong?any help would be highly appreciated. Thank you for your time.
Soroush
-
-
Gurobi Staff
Hi Soroush,
You are trying to access \(\texttt{x[t+l]}\) and you loop over \(t \in \{0,\dots 49\}, l\in \{0,\dots,5\}\). So at some point you are trying to access \(\texttt{x[45+5]}\) which is not defined, because you have only 50 \(\texttt{x}\) variables with indices \(0,\dots,49\).
Best regards,
Jaromił1 -
-
Dear Jaromił,
Thank you so much for your fast response! I appreciate it. I will probably run into other errors and get back to you.
BR,
Soroush
0 -
-
Dear Jaromił,
I have figured out the indexes problem somehow, however, the solution I get seems to be wrong.
I believe there might be something wrong with my code of the objective function:m.setObjective(gp.quicksum((x[t+l]-r)**2 + b[t+l] for l in range(1,P)) , GRB.MINIMIZE)
If we were to solve the problem by hand for the "first term" in objective function we would have:(x[0+1]-3)^2since x[1] = x[0] + u[0] - d[0]= 0 + u[0] - 1therefore in order to minimize the objective function:(u[0] -1 - 3)^2we have to choose u[0] = 4, however, the optimal solution by Gurobi is u[0] = 1.0Here is my modified code block:# Prediction horizon and time constants
I do appreciate any help.
P = 6 #prediction horizon(days)
K = P-1 # counter(days)
T = 55 # Total time horizon(days)
# Parameters
r = 3 #desired inventory level(constant)
# Control variable limits
u_min=0
u_max = 7
# state and input variables initialization (variables)
x0 = 0
d =[1, 0, 3, 4, 5, 4, 0, 0, 4, 3, 2, 5, 4, 0, 1, 1, 6, 1, 1, 0, 4, 1,
0, 6, 3, 6, 2, 1, 6, 0, 6, 0, 3, 3, 2, 4, 0, 2, 6, 4, 1, 6, 6, 5,
4, 0, 6, 0, 4, 0, 1, 0, 3, 2, 0, 1, 5, 4, 5, 5, 3, 4, 6, 3, 3, 1,
0, 3, 6, 1] #for debuging purpose we fixed the demand vector
# Create a new model
m = Model("MPCmodel")
# Create variables
x = m.addVars(T, vtype=GRB.INTEGER, name="x")
b = m.addVars(T, vtype=GRB.INTEGER, name="b")
#w = m.addVars(T+1, vtype=GRB.INTEGER, name="w")
u = m.addVars(T, lb=u_min , ub=u_max, vtype=GRB.INTEGER, name="u")
x[0] = x0
for t in range(50):
# Constraints
m.addConstrs((x[k+t] + u[k+t] - d[k+t] == x[k+t+1] for k in range(K)),"c1")
m.addConstrs((d[k+t] - x[k+t] == b[k+t] for k in range(K)),"c2")
m.addConstrs((b[k+t] >=0 for k in range(K)),"c3")
#m.addConstrs((max(u[k+t-14] - gp.quicksum(d[k+t-n] for n in range(1,15)), 0) == w[k+t] for k in range(K) for t in range(14,T+1)),"c3")
# Objective function
m.setObjective(gp.quicksum((x[t+l]-r)**2 + b[t+l] for l in range(1,P)) , GRB.MINIMIZE)
m.write('MPCmodel.lp')
m.write('model.mps')
m.optimize()
BR,Soroush0 -
-
-
Hi Soroush,
Your objective is indented to be in the \(\texttt{for}\)-loop. You should un-indent it and add an addition \(\texttt{for}\)-loop to the quicksum.
# Objective function
m.setObjective(gp.quicksum((x[t+l]-r)**2 + b[t+l] for l in range(1,P) for t in range(50)) , GRB.MINIMIZE)You also should have a look at the LP file MPCmodel.lp you are writing to check whether your model looks like what you would actually expect.
If you think that there is a better optimal solution available, you should try providing this solution as a MIP start. This way, Gurobi will tell you whether your MIP start is feasible and if yes, then it will tell you the objective value of your feasible solution point.
Moreover, in your code there are 2 \(\max\) statements missing which are present in the formulation you posted. You can model \(\max\) terms via Gurobi's general constraint feature. Note that you have to introduce auxiliary variables for each expression in the \(\max\) term, i.e., in your case you have to model
\[\begin{align*}
z_{k+t} &= d_{k+t} - x_{k+t}\\
b_{k+t} &= \max\{z_{k+t},0\}\\
y_{k+t} &= u_{k+t-14} - \sum d_{k+t-14}\\
w_{k+t} &= \max\{y_{k+t},0\}
\end{align*}\]
where \(z,y\) are free continuous variables.Best regards,
Jaromił0 -
-
Dear Jaromil,
Thank you so much for you help, it means a lot to me, considering you probably have to answer a lot of questions daily. I do appreciate it; however, my problem persists, maybe giving a more detailed explanation would help.
As I posted my model, this model tries to minimize the deviation of inventory level from a desired setpoint “r”, minimize the number of lost orders (or backlogs) and minimize the wasted items during a planning horizon.
The algorithm to solve this problem is as follows:
- At time t solve the optimization problem J with regards to the respective constraints(uploaded model on top of the page);
- feed u(t) to the system
- iterate with t = t+1
I have coded up the problem using your recommendations and yet I have three main questions:
- Considering the auxiliary variables, I get this error:
TypeError: must be real number, not gurobipy.LinExpr
- As the algorithm says, is there any way to get the first u(t) for each t and feed it to the next iteration
- The objective function is a summation of two “positive” terms, why there is a couple of minus x[i] in the MPCmodel.lp file, and why there is a denominator of 2 at the end?
I was wondering if maybe this explanation can clarify my problem.
Once again, I really appreciate your time, consideration, and help.
My final code with your recommendations:# Prediction horizon and time constants
P = 6 #prediction horizon(days)
K = P-1 # counter(days)
T = 55 # Total time horizon(days)
# Parameters
r = 3 #desired inventory level(constant)
# Control variable limits
u_min=0
u_max = 5
# state and input variables initialization (variables)
x0 = 0
d =[ 1, 0, 3, 4, 5, 4, 0, 0, 4, 3, 2, 5, 4, 0, 1, 1, 6, 1, 1, 0, 4, 1,
0, 6, 3, 6, 2, 1, 6, 0, 6, 0, 3, 3, 2, 4, 0, 2, 6, 4, 1, 6, 6, 5,
4, 0, 6, 0, 4, 0, 1, 0, 3, 2, 0, 1, 5, 4, 5, 5, 3, 4, 6, 3, 3, 1,
0, 3, 6, 1
]
# Create a new model
m = Model("MPCmodel")
# Create variables
x = m.addVars(T, vtype=GRB.INTEGER, name="x")
z = m.addVar(T,lb=-GRB.INFINITY,ub=GRB.INFINITY,vtype=GRB.CONTINUOUS,name="z_aux_var")
y = m.addvar(T,lb=-GRB.INFINITY,ub=GRB.INFINITY,vtype=GRB.CONTINUOUS,name="y_aux_var")
m.addVar((z == d[k+t]-x[k+t] for k in range(K)), name="aux_var1")
m.addVar(((y == u[k+t-14] - gp.quicksum(d[k+t-n] for n in range(1,15)), 0) for k in range(K) for t in range(14,T)), name="aux_var2")
b = m.addVars(T, vtype=GRB.INTEGER, name="b")
w = m.addVars(T, vtype=GRB.INTEGER, name="w")
u = m.addVars(T, lb=u_min , ub=u_max, vtype=GRB.INTEGER, name="u")
x[0] = x0
for t in range(50):
# Constraints
m.addConstrs((x[k+t] + u[k+t] - d[k+t] == x[k+t+1] for k in range(K)),"c1")
m.addConstrs((b[k+t] == max_(0,z) for k in range(K)),"c2")
m.addConstrs((w[k+t] == max_(0,y)for k in range(K)),"c3")
# Objective function
m.setObjective(gp.quicksum((x[t+l]-r)**2 + b[t+l] for l in range(1,P) for t in range(50)) , GRB.MINIMIZE)
m.write('MPCmodel.lp')
m.write('model.mps')
m.optimize()0 -
-
-
Hi Soroush,
- Considering the auxiliary variables, I get this error:
TypeError: must be real number, not gurobipy.LinExpr
You have to respect the method's arguments. For the auxiliary variables \(z,y\), you should use the addVars method, because you want to add multiple variables.
The terms
m.addVar((z == d[k+t]-x[k+t] for k in range(K)), name="aux_var1")
m.addVar(((y == u[k+t-14] - gp.quicksum(d[k+t-n] for n in range(1,15)), 0) for k in range(K) for t in range(14,T)), name="aux_var2")are not variables but constraints. Thus, you should use the addConstrs method. Moreover, these auxiliary constraints have to be added in the \(\texttt{t for}\)-loop. The additional condition
for t in range(14,T)
cannot be used in the \(\texttt{t for}\)-loop because you would have two \(\texttt{t}\) variables. But you can use an \(\texttt{if}\)-condition.
import gurobipy as gp
from gurobipy import GRB
import numpy as np
# Prediction horizon and time constants
P = 6 #prediction horizon(days)
K = P-1 # counter(days)
T = 55 # Total time horizon(days)
# Parameters
r = 3 #desired inventory level(constant)
# Control variable limits
u_min=0
u_max = 5
# state and input variables initialization (variables)
x0 = 0
d =[1, 0, 3, 4, 5, 4, 0, 0, 4, 3, 2, 5, 4, 0, 1, 1, 6, 1, 1, 0, 4, 1,
0, 6, 3, 6, 2, 1, 6, 0, 6, 0, 3, 3, 2, 4, 0, 2, 6, 4, 1, 6, 6, 5,
4, 0, 6, 0, 4, 0, 1, 0, 3, 2, 0, 1, 5, 4, 5, 5, 3, 4, 6, 3, 3, 1,
0, 3, 6, 1]
# Create a new model
m = gp.Model("MPCmodel")
# Create variables
x = m.addVars(T, vtype=GRB.INTEGER, name="x")
z = m.addVars(T,lb=-GRB.INFINITY,ub=GRB.INFINITY,vtype=GRB.CONTINUOUS,name="z_aux_var")
y = m.addVars(T,lb=-GRB.INFINITY,ub=GRB.INFINITY,vtype=GRB.CONTINUOUS,name="y_aux_var")
b = m.addVars(T, vtype=GRB.INTEGER, name="b")
w = m.addVars(T, vtype=GRB.INTEGER, name="w")
u = m.addVars(T, lb=u_min , ub=u_max, vtype=GRB.INTEGER, name="u")
x[0] = x0
for t in range(50):
# Constraints
m.addConstrs((x[k+t] + u[k+t] - d[k+t] == x[k+t+1] for k in range(K)),"c1")
m.addConstrs((z[k+t] == d[k+t]-x[k+t] for k in range(K)), name="aux_constr1")
m.addConstrs((b[k+t] == gp.max_(0,z[k+t]) for k in range(K)),"c2")
if t >= 14:
m.addConstrs((y[k+t] == u[k+t-14] - gp.quicksum(d[k+t-n] for n in range(1,15)) for k in range(K)), name="aux_constr2")
m.addConstrs((w[k+t] == gp.max_(0,y[k+t]) for k in range(K)),"c3")
# Objective function
m.setObjective(gp.quicksum((x[t+l]-r)**2 + b[t+l] for l in range(1,P) for t in range(50)) , GRB.MINIMIZE)
m.write('MPCmodel.lp')
m.optimize()As the algorithm says, is there any way to get the first u(t) for each t and feed it to the next iteration
You can extract solution information via a variable's X attribute. An example would be
m.optimize()
for t in T:
print("%s: %f"%(u[t].VarName, u[t].X))So you can extract solution values of certain variables, feed them to your model and adjust it accordingly, and resolve your model.
The objective function is a summation of two “positive” terms, why there is a couple of minus x[i] in the MPCmodel.lp file, and why there is a denominator of 2 at the end?
Your objective is the sum of many terms of the form \((x_{t+l} - r)^2 + b_{t+l}\). Gurobi separates expressions into a linear and a quadratic part so the term \((x_{t+l} - r)^2 + b_{t+l}\) becomes \(b_{t+l} - 2\cdot r \cdot x_{t+l} + r^2 + x_{t+l}^2\). Summing over all \(t,l\) results in the objective function you see in the LP file.
The division by \(2\) appears, because Gurobi uses the standard form for quadratic models, which uses a symmetric matrix Q divided by 2. For additional insights about (nonconvex) quadratic models in Gurobi, I would recommend to have a look at our Non-Convex webinar.
Best regards,
Jaromił0 -
-
Dear Jaromił,
I am deeply thankful for this help and your perfect clarification. It solved my problem.
Thank you so much for your detailed answer, it is probably the best experience of mine on any kind of forum or something like that.
Again thank you and Gurobi for your amazing support.BR,
Soroush
0 -
Please sign in to leave a comment.
Comments
7 comments