How to assign a new variable equal to min of two other variable (Gurobi in Python)
Ongoing
Hi,
I have been trying to write a variable say c = min(a, b) and then this variable enters the objective function. The problem is I can't use c = min_(a,b) because then c becomes a gurobipy.GenExpr and I can't do summation in the objective function (unsupported operand type(s) for +: 'GenExpr' and 'GenExpr').
How to fix this? I have been trying this issue for a long time and still can't figure it out. Thank you so much in advance for your help!
Since the whole code is very long, here I'm only sharing a part of it first to give you context (x, y, z are all variables, other are parameters).
-
Hi Ye,
1. MIN, MAX Constraints:
you are looking for a min_ constraint, as described here.
There are also max_ constraints, described here, which generally work similarly.
Assuming that you have defined your variables someplace earlier in your code, its part could look as follows:
if P[2,i]>(1-alpha)*P[0,i]+P[4,i]:
y1_1[0,i] = 0
m.addConstr(y2_1[0,i] = min_(s2_1[0,i],q))
m.addConstr(x1_1[0,i] = max_(0,s1_1[0,i]-z_1_x1[0,i]))
m.addConstr(x2_1[0,i] = max_(0, s2_1[0,i]-q-z_2_x1[0,i]))On a side note:
y1_1[0,i] = 0
will put an integer 0 in your dictionary y1_1, where you store variables otherwise. Is that the expected behavior?
2. Sum in objective function:
If your R1 is a Gurobi's tupledict, you may generate a sum of all variables contained therein by:
R1.sum()
Hope this helps.
Best regards
Jonasz0 -
Hi Ye,
In addition to what Jonasz already stated, the general constraints methods accept Var objects only. This means that you cannot directly provide a sum or difference of two variable to the min_ / max_ function and have to introduce additional auxiliary variables.
[...]
auxVar = m.addVar(lb=-GRB.INFINITY,name="auxVar_%d"%(i))
m.addVar(auxVar == s1_1[0,i]-z_1_x1[0,i], name="auxConstr_%d"%(i))
m.addConstr(x1_1[0,i] = max_(0,auxVar)Moreover, the \(\texttt{if}\)-clauses will work only if your \(P\) object holds pre-defined parameter values. It will not work if \(P\) are optimization variables. In the case that \(P\) are indeed optimization variables, you have to introduce additional binary variables and formulate your model with the help of conditional statements. Please have a look at How do I model conditional statements in Gurobi? for more details.
Best regards,
Jaromił0 -
Thank you so much, Jonasz and Jaromił! I really appreciated it!
- Jonasz, yes, I tried writing into constraints before, but I didn't know the trick of writing the objective function as R1.sum(). Also thanks for pointing out y1_1[0,i] = 0, I also wrote it into a constraint to avoid further issues.
- Jaromił, thanks so much for pointing this out. I was struggling with this as well! Thanks for the suggestion!
Best,
Ye
0 -
I do have another question regarding writing the constraint:
y1_1[i] = min( s1_1[i], max(q-s2_1[i], 0) )I wrote it into two constraints by adding two additional variables y1_10 and y_11, but it seems when I added two more constrains, then the model became unsolvable. I wrote as follows:m_ML2.addConstr(y1_10 == q - s2_1[i])
m_ML2.addConstr(y1_11 == gp.max_(y1_10,0))
m_ML2.addConstr(y1_1[i] == gp.min_(s1_1[i], y1_11) )The error is: Model is infeasible or unbounded. When I only used one constraint it worked, but it's not the constraint I wanted.m_ML2.addConstr(y1_1[i] == q - s2_1[i]) #only write this one worksI would really appreciate it if someone could answer this question. It is really confusing why I cannot specify several constraints related to one variable.Thanks,Ye0 -
Did you explicitly set the auxiliary variables to be free by setting their lower bound to \(-\infty\)?
auxVar = m.addVar(lb=-GRB.INFINITY,name="auxVar")You can check the Knowledge Base article How do I resolve the error "model is infeasible or unbounded"?
If you did this already, could you please provide a minimal reproducible example showing the issue?
Best regards,
Jaromił0 -
Thanks, Jaromił for answering this question! I tried what you suggested, the model is still not solvable when T=3 (it works when T=2). Here is the minimal reproducible example showing the issue:
import numpy as np
import gurobipy as gp
from gurobipy import GRB, abs_
T = 3
q = 98107
alpha = 0.154
beta = 0.02
P0 = np.random.normal(30, 4, T)
P = np.array((P0, P0-10, P0, P0-20, P0-25))
w1 = np.random.normal(28937.89109,7875.353526,T+1)#from data, mean, sd, size
w2 = np.random.normal(62906.31683,9551.438123,T+1) #from data
try:
m_ML2 = gp.Model("Model")
# define decision variables
y1_10 = m_ML2.addVar(vtype = GRB.CONTINUOUS, name = 'y10', lb=-GRB.INFINITY)
x1_1 = m_ML2.addVars(T, name = 'x1_1', lb = 0)
x2_1 = m_ML2.addVars(T, name = 'x2_1', lb = 0)
y1_1 = m_ML2.addVars(T, name = 'y1_1', lb = 0)
y2_1 = m_ML2.addVars(T, name = 'y2_1', lb = 0)
z1_1 = m_ML2.addVars(T, name = 'z1_1', lb = 0)
z2_1 = m_ML2.addVars(T, name = 'z2_1', lb = 0)
s1_1 = m_ML2.addVars(T, name = 's1_1', lb = 0)
s2_1 = m_ML2.addVars(T, name = 's2_1', lb = 0)
R1 = m_ML2.addVars(T, name = 'R1')
#period 0
m_ML2.addConstr(z1_1[0] == 5000)
m_ML2.addConstr(z2_1[0] == 10000)
m_ML2.addConstr(y2_1[0] == 0)
m_ML2.addConstr(x2_1[0] == w2[0])
m_ML2.addConstr(y1_1[0] == 0 )
m_ML2.addConstr(x1_1[0] == w1[0])
#reward
m_ML2.addConstr(R1[0] == P[0,0]*(y2_1[0]+(1-alpha)*y1_1[0])+P[1,0]*x2_1[0]+P[2,0]*x1_1[0]-P[4,0]*(q-y1_1[0]-y2_1[0])-P[3,0]*(w1[0]+w2[0]))
for i in range(1,T):
m_ML2.addConstr(z1_1[i] == 5000)
m_ML2.addConstr(z2_1[i] == 10000)
m_ML2.addConstr(s1_1[i] == w1[i])
m_ML2.addConstr(s2_1[i] == w2[i]+z1_1[i-1]+z2_1[i-1])
m_ML2.addConstr(y2_1[i] == gp.min_(q,s2_1[i]) )
m_ML2.addConstr(x2_1[i] == 0)
# m_ML2.addConstr(y1_1[i] == q - y2_1[i]) #only using this one works instead of the following two constrsm_ML2.addConstr(y1_10 == q - y2_1[i])
m_ML2.addConstr(y1_1[i] == gp.min_(s1_1[i], y1_10) )
m_ML2.addConstr(x1_1[i] == 0 )
#reward
m_ML2.addConstr(R1[i] == P[0,i]*(y2_1[i]+(1-alpha)*y1_1[i])+P[1,i]*x2_1[i]+P[2,i]*x1_1[i]-P[3,i]*(z1_1[i]+z2_1[i])-P[4,i]*(q-y1_1[i]-y2_1[i])-P[3,i]*(w1[i]+w2[i]))
m_ML2.setObjective(R1.sum(), GRB.MAXIMIZE)
# Optimize model
m_ML2.optimize()
for v in m_ML2.getVars():
print('%s %g' % (v.varName, v.x))print('Obj: %g' % m_ML2.objVal)
except gp.GurobiError as e:
print('Error code ' + str(e.errno) + ': ' + str(e))except AttributeError:
print('Encountered an attribute error')0 -
If you disable DualReductions then compute an IIS
m_ML2.setParam("DualReductions",0)
[...]
m_ML2.optimize()
m_ML2.computeIIS()
m_ML2.write("iis.ilp")and have a look at the IIS file \(\texttt{iis.ilp}\) you see that
R0: z1_1[0] = 5000
R1: z2_1[0] = 10000
R10: - z1_1[0] - z2_1[0] + s2_1[1] = 59970.24755004315
R12: y10 + y2_1[1] = 98107
R18: - z1_1[1] - z2_1[1] + s2_1[2] = 78968.96449813887
R20: y10 + y2_1[2] = 98107
Bounds
y10 free
y2_1[1] free
y2_1[2] free
z1_1[0] free
z2_1[0] free
s2_1[1] free
s2_1[2] free
General Constraints
GC0: y2_1[1] = MIN ( s2_1[1] , 98107 )
GC2: y2_1[2] = MIN ( s2_1[2] , 98107 )this set of constraints and variables cannot be fulfilled. It is because
\[\begin{align}
s_{2,1,1} &= 59970.24755004315 + 10000 + 5000 = 74970.24755004315\\
s_{2,1,2} &= 78968.96449813887 + 10000 + 5000 = 93968.96449813887\\
y_{2,1,1} &= \min(74970.24755004315, 98107) = 74970.24755004315\\
y_{2,1,2} &= \min(93968.96449813887, 98107) = 93968.96449813887\\
\end{align}\]
And then you also finally have\[\begin{align}
y_{10} &= 98107 - 74970.24755004315\\
y_{10} &= 98107 - 93968.96449813887
\end{align}\]
which obviously cannot hold.Now it is up to you to double check your model and understand why this happens. My guess is that you probably need more than one \(y_{10}\). You probably need \(y_{10i}\) with one \(i\) for every \(\min\) term.
Best regards,
Jaromił0 -
Thank you so much, Jaromił!!! I redefine y_10 as a vector then it works! I thought the new value of y_10 will overwrite the previous value so that's why I only used one variable y_10 before.
The method you adopted to find out what went wrong is very helpful and it can help me debug in the future! I wonder how you found iis.ilp file. Since I am using Gurobi in Spyder (using Anaconda), when I wrote the following,
m_ML2.setParam("DualReductions",0)
[...]
m_ML2.optimize()
m_ML2.computeIIS()
m_ML2.write("iis.ilp")I only have this in the output console instead of what you got (some specific numerical values):
Computing Irreducible Inconsistent Subsystem (IIS)...
Constraints | Bounds | Runtime
Min Max Guess | Min Max Guess |
--------------------------------------------------------------------------
0 27 - 0 27 - 0s
9 9 - 1 1 - 0sIIS computed: 9 constraints, 1 bounds
IIS runtime: 0.07 seconds
Encountered an attribute errorThanks a lot!
Best,
Ye
0 -
Hi Ye,
m_ML2.write("iis.ilp")will generate the \(\texttt{iis.ilp}\) file in the working directory of your Spyder/Anaconda installation. It is easiest to just search your machine for it. It might be possible that you are missing writing permissions. Then your Spyder/Anaconda environment will not generate any files. In this case, it is best to execute your IDE as admin.
Please note that the provided code only works when your model is guaranteed to be infeasible, because an IIS can only be computed for infeasible models. A correct way of implementing an IIS computation would be
m_ML2.optimize()
if (m_ML2.Status == GRB.INFEASIBLE):
m_ML2.computeIIS()
m_ML2.write("iis.ilp")Best regards,
Jaromił0 -
Jaromił, thank you so much for your answers, it worked! You really helped me move forward with my research project!!!
One more question regarding your previous response on How do I model conditional statements in Gurobi? I encountered a very similar situation.
The problem looks like this:
m.addConstr(x == c1*x1+c2*x2+c3*x3)
m.addConstr( x>=y+eps-M*(1-b) )
m.addConstr( x<=y+ M*b )
The only difference is that x is a linear expression of several optimization vars: x1, x2, x3. I have the error: Unsupported type (<class 'numpy.ndarray'>) for LinExpr addition argument.
I would very much appreciate your help! Thank you in advance!
Best,
Ye
0 -
Hi Ye,
Could you share how you define \(x\)? The error indicates that a numpy array of parameters is used instead of a variable tupledict. Could you please try to construct a minimal working example reproducing the issue?
Best regards,
Jaromił0 -
Thank you so much, Jaromił! I really appreciate your response! When I was trying to construct a minimal working example, I found out that I used a vector instead of a scaler for the coefficient for my x. Then the big M approach you introduced in the post (How do I model conditional statements in Gurobi?) worked!
One more question about resetting the model. I need to run this model more than one time. Since each period has different prices (and other exogenous parameters) so the model has different parameters in each run. When I run the model a lot of times, then the speed becomes very slow, it took me multiple days to run one time. I wonder if this is something I can improve upon.
I used resetParams() to reset parameters for each run:
m_ML2 = gp.Model()
m_ML2.reset(0) # since I am running model multiple times via looping, I reset the model for each run
m_ML2.resetParams()The issue looks like this (time here in the last column can take up to thousands of seconds):
Attached is the minimal working example. When Tmax is a small number (around d<30), it works. But if Tmax= 101 and n>2, then this situation occurs.
import random
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
import gurobipy as gp
import numpy as np
from gurobipy import GRB, abs_
#set the numer of run times here
Tmax = 100
nmax = 1rng = np.random.default_rng(123)
q = 98107
alpha = 0.154
beta = 0.02
Twarmup = 101
eps = 0.00002
M = 1000000
W0_factor = np.array([[-0.050083, -0.13667, 0.0020002], #c
[0.62053, 0.36688, 0.8887], # ar1
[0.52539, 0.39767, -0.70058], # ma1
[4.6059, 3.7998, 0.05894]]) # var
#Pi_O Pi_M Pi_F
W0_OTC = np.array([[14.617, -34.846, -8.940, 1.037, -0.089], #c
[0.165, -0.274, -0.183, 0.806, -0.183], # ar1
[0.358, 0.342, 0.454, -0.253, 0.454], # ma1
[0.11, 0.006, -0.027, 0.026, -0.027], # sma
[0.876, 2.207, 2.284, 0.005, 0.023], #Pi_O
[0.673, 0.471, 0.328, -0.007, 0.003], #Pi_M
[0.237, 1.101, 0.191, 0.001, 0.002], #Pi_F
[4.839, 47.704, 16.58, 0.020227, 0.0016579]]) #var
#p_c p_1 p_2 c_h c_p
# Cholesky decomposition matrix
cov = [[5, 0.5, 0.4, 0.3, 0.2],
[0.5, 4, 0.4, 0.2, 0.1],
[0.4, 0.5, 3, 0.3, 0.2],
[0.3, 0.2, 0.15, 2, 0.1],
[0.3, 0.2, 0.15, 0.1, 1]]
L = np.linalg.cholesky(cov)V_ML = np.zeros((Tmax-1))
for T in range(1, Tmax):
print(T)
V2 = np.zeros((nmax))
for n in range(nmax):
Pi = np.zeros((3, Tmax))
P = np.zeros((5, Tmax))
pF = np.zeros((1, Tmax))
tempPi = np.zeros((3, Tmax))
temp = np.zeros((5,Tmax+1))
w1 = np.random.normal(28937.89109,7875.353526,T+1)#from data, mean, sd, size
w2 = np.random.normal(62906.31683,9551.438123,T+1) #from data# Warm up factor price
u_O = np.random.normal(0,np.sqrt(W0_factor[3,0]),Twarmup+53)
u_M = np.random.normal(0,np.sqrt(W0_factor[3,1]),Twarmup+53)
u_F = np.random.normal(0,np.sqrt(W0_factor[3,2]),Twarmup+53)
Pi_O = np.zeros(Twarmup+53)
Pi_M = np.zeros(Twarmup+53)
Pi_F = np.zeros(Twarmup+53)
for t in range(2,len(u_O)):
Pi_O[t] = W0_factor[0,0] + (1 + W0_factor[1,0])*Pi_O[t-1] - W0_factor[1,0]*Pi_O[t-2] + u_O[t] + W0_factor[2,0]*u_O[t-1]
for t in range(2,len(u_M)):
Pi_M[t] = W0_factor[0,1] + (1 + W0_factor[1,1])*Pi_M[t-1] - W0_factor[1,1]*Pi_M[t-2] + u_M[t] + W0_factor[2,1]*u_M[t-1]
for t in range(2,len(u_F)):
Pi_F[t] = W0_factor[0,2] + (1 + W0_factor[1,2])*Pi_F[t-1] - W0_factor[1,2]*Pi_F[t-2] + u_F[t] + W0_factor[2,2]*u_F[t-1]
Pi_O = Pi_O + 25
Pi_M = Pi_M + 55
Pi_F = Pi_F + 2.5
Pi_O[Pi_O<0]=0
Pi_M[Pi_M<0]=0
Pi_F[Pi_F<0]=0
#warm up OTC prices
w_C = np.random.normal(0,np.sqrt(W0_OTC[7,0]),Twarmup+53)
w_1 = np.random.normal(0,np.sqrt(W0_OTC[7,1]),Twarmup+53)
w_2 = np.random.normal(0,np.sqrt(W0_OTC[7,2]),Twarmup+53)
w_H = np.random.normal(0,np.sqrt(W0_OTC[7,3]),Twarmup+53)
w_P = np.random.normal(0,np.sqrt(W0_OTC[7,4]),Twarmup+53)
P_C = np.zeros(Twarmup+53)
P_1 = np.zeros(Twarmup+53)
P_2 = np.zeros(Twarmup+53)
C_H = np.zeros(Twarmup+53)
C_P = np.zeros(Twarmup+53)
for t in range(53,len(w_C)): ## may be rewritten as function later
P_C[t] = W0_OTC[0,0] + P_C[t-1]*W0_OTC[1,0] + Pi_O[t]*W0_OTC[4,0] + Pi_M[t]*W0_OTC[5,0] + Pi_F[t]*W0_OTC[6,0] + w_C[t] + W0_OTC[3,0]*w_C[t-52] + W0_OTC[2,0]*w_C[t-1] + W0_OTC[2,0]*W0_OTC[3,0]*w_C[t-53]
for t in range(53,len(w_1)):
P_1[t] = W0_OTC[0,1] + P_1[t-1]*W0_OTC[1,1] + Pi_O[t]*W0_OTC[4,1] + Pi_M[t]*W0_OTC[5,1] + Pi_F[t]*W0_OTC[6,1] + w_1[t] + W0_OTC[3,1]*w_1[t-52] + W0_OTC[2,1]*w_1[t-1] + W0_OTC[2,1]*W0_OTC[3,1]*w_1[t-53]
for t in range(53,len(w_2)):
P_2[t] = W0_OTC[0,2] + P_2[t-1]*W0_OTC[1,2] + Pi_O[t]*W0_OTC[4,2] + Pi_M[t]*W0_OTC[5,2] + Pi_F[t]*W0_OTC[6,2] + w_2[t] + W0_OTC[3,2]*w_2[t-52] + W0_OTC[2,2]*w_2[t-1] + W0_OTC[2,2]*W0_OTC[3,2]*w_2[t-53]
for t in range(53,len(w_H)):
C_H[t] = W0_OTC[0,3] + C_H[t-1]*W0_OTC[1,3] + Pi_O[t]*W0_OTC[4,3] + Pi_M[t]*W0_OTC[5,3] + Pi_F[t]*W0_OTC[6,3] + w_H[t] + W0_OTC[3,3]*w_H[t-52] + W0_OTC[2,3]*w_H[t-1] + W0_OTC[2,3]*W0_OTC[3,3]*w_H[t-53]
for t in range(53,len(w_P)):
C_P[t] = W0_OTC[0,4] + C_P[t-1]*W0_OTC[1,4] + Pi_O[t]*W0_OTC[4,4] + Pi_M[t]*W0_OTC[5,4] + Pi_F[t]*W0_OTC[6,4] + w_P[t] + W0_OTC[3,4]*w_P[t-52] + W0_OTC[2,4]*w_P[t-1] + W0_OTC[2,4]*W0_OTC[3,4]*w_P[t-53]
P_C[P_C<0]=0
P_1[P_1<0]=0
P_2[P_2<0]=0
C_H[C_H<0]=0
C_P[C_P<0]=0Pi[:,0] = np.stack((Pi_O[-1], Pi_M[-1], Pi_F[-1]))
#OTC prices fist period
P[:,0] = np.stack((P_C[-1], P_2[-1], P_1[-1], C_H[-1], C_P[-1]))
pF[0,0] = C_H[-1]+rng.normal(0,C_H[-1]/100) #mean and standard deviation, forward market price for the first period of the spot market price
# generate 2nd to 100 period prices
W_factor = np.array([[-0.050083, 1+0.62053, -0.62053],
[-0.13667, 1+0.36688, -0.36688],
[0.0020002, 1+0.8887, -0.8887]])
tempPi[0,0] = W_factor[0,:]@np.array([1,Pi_O[-1],0]) #inner product
tempPi[1,0] = W_factor[1,:]@np.array([1, Pi_M[-1], 0])
tempPi[2,0] = W_factor[2,:]@np.array([1, Pi_F[-1], 0])
tempPi[tempPi<0]=0
#generate forecast for the OTC
W_OTC = np.array([[14.617, 0.16547, 0.87606, 0.67273, 0.23732], #P_C
[-8.9395, -0.18333, 2.2838, 0.32753, 0.19068], #P_2
[-34.846, -0.27437, 2.2071, 0.47132, 1.1006], #P_1
[1.0374, 0.80573, 0.004734, -0.0065202, 0.0012593], #C_H
[-0.089304, -0.18324, 0.022836, 0.0032739, 0.0019056]]) #C_P
temp[0,0] = W_OTC[0,:]@np.hstack((1, P_C[-1], tempPi[:,0])) #P_C
temp[1,0] = W_OTC[1,:]@np.hstack((1, P_2[-1], tempPi[:,0])) #P_2
temp[2,0] = W_OTC[2,:]@np.hstack((1, P_1[-1], tempPi[:,0])) #P_1
temp[3,0] = W_OTC[3,:]@np.hstack((1, C_H[-1], tempPi[:,0])) #C_H
temp[4,0] = W_OTC[4,:]@np.hstack((1, C_P[-1], tempPi[:,0])) #C_P
temp[temp<0]=0
for i in range(1,T):
#Factor prices DGP
inno_Pi_O = np.random.normal(0,np.sqrt(4.6059),T+1)
inno_Pi_M = np.random.normal(0,np.sqrt(3.7998),T+1)
inno_Pi_F = np.random.normal(0,np.sqrt(0.05894),T+1)
MA_factor = np.array([0.52539, 0.39767, -0.70058])
if i == 1:
Pi[0,i] = W_factor[0,:]@np.hstack([1, Pi[0,i-1], 0]) + inno_Pi_O[i] + MA_factor[0]*inno_Pi_O[i-1]
Pi[1,i] = W_factor[1,:]@np.hstack([1, Pi[1,i-1], 0]) + inno_Pi_M[i] + MA_factor[1]*inno_Pi_M[i-1]
Pi[2,i] = W_factor[2,:]@np.hstack([1, Pi[2,i-1], 0]) + inno_Pi_F[i] + MA_factor[2]*inno_Pi_F[i-1]
else:
Pi[0,i] = W_factor[0,:]@np.hstack([1, Pi[0,i-1], Pi[0,i-2]])+ inno_Pi_O[i] + MA_factor[0]*inno_Pi_O[i-1]
Pi[1,i] = W_factor[1,:]@np.hstack([1, Pi[1,i-1], Pi[1,i-2]]) + inno_Pi_M[i] + MA_factor[1]*inno_Pi_M[i-1]
Pi[2,i] = W_factor[2,:]@np.hstack([1, Pi[2,i-1], Pi[2,i-2]]) + inno_Pi_F[i] + MA_factor[2]*inno_Pi_F[i-1]
Pi[Pi<0]=0
# generate forecast/expectation for the next period's factor price, t=2
# use the estimated factor price parameters
tempPi[0,i] = W_factor[0,:]@np.array([1, Pi[0,i], Pi[0,i-1]])#forecast for the ith period
tempPi[1,i] = W_factor[1,:]@np.array([1, Pi[1,i], Pi[1,i-1]])
tempPi[2,i] = W_factor[2,:]@[1, Pi[2,i], Pi[2,i-1]]
tempPi[tempPi<0]=0
#OTC prices DGP
#generate forecast for the OTC
inno_P_C = np.random.normal(0,np.sqrt(4.839),T+52)
inno_P_2 = np.random.normal(0,np.sqrt(16.58),T+52)
inno_P_1 = np.random.normal(0,np.sqrt(47.704),T+52)
inno_C_H = np.random.normal(0,np.sqrt(0.020227),T+52)
inno_C_P = np.random.normal(0,np.sqrt(0.0016579),T+52)
MA_OTC = np.array([0.35802, 0.45427, 0.34237, -0.25309, 0.45424])
SMA_OTC = np.array([0.11019, -0.027154, 0.0063966, 0.024557, -0.027257])
P[0,i] = W_OTC[0,:]@np.hstack([1, P[0,i-1], tempPi[:,i]]) + inno_P_C[i+52] + MA_OTC[0]*inno_P_C[i+51] + SMA_OTC[0]*inno_P_C[i] + MA_OTC[0]*SMA_OTC[0]*inno_P_C[i-1] #P_C
P[1,i] = W_OTC[1,:]@np.hstack([1, P[1,i-1], tempPi[:,i]]) + inno_P_2[i+52] + MA_OTC[1]*inno_P_2[i+51] + SMA_OTC[1]*inno_P_2[i] + MA_OTC[1]*SMA_OTC[1]*inno_P_2[i-1] #P_2
P[2,i] = W_OTC[2,:]@np.hstack([1, P[2,i-1], tempPi[:,i]]) + inno_P_1[i+52] + MA_OTC[2]*inno_P_1[i+51] + SMA_OTC[2]*inno_P_1[i] + MA_OTC[2]*SMA_OTC[2]*inno_P_1[i-1] #P_1
P[3,i] = W_OTC[3,:]@np.hstack([1, P[3,i-1], tempPi[:,i]]) + inno_C_H[i+52] + MA_OTC[3]*inno_C_H[i+51] + SMA_OTC[3]*inno_C_H[i] + MA_OTC[3]*SMA_OTC[3]*inno_C_H[i-1] #C_H
P[4,i] = W_OTC[4,:]@np.hstack([1, P[4,i-1], tempPi[:,i]]) + inno_C_P[i+52] + MA_OTC[4]*inno_C_P[i+51] + SMA_OTC[4]*inno_C_P[i] + MA_OTC[4]*SMA_OTC[4]*inno_C_P[i-1] #C_P
P[P<0]=0
#forward market
pF[0,i] = P[3,i]+rng.normal(0,P[3,i]/100) #mean and standard deviation
pF[pF<0]=0
#generate forecast for the OTC
temp[0,i] = W_OTC[0,:]@np.hstack([1, P[0,i], tempPi[:,i]]) #P_C
temp[1,i] = W_OTC[1,:]@np.hstack([1, P[1,i], tempPi[:,i]]) #P_2
temp[2,i] = W_OTC[2,:]@np.hstack([1, P[2,i], tempPi[:,i]]) #P_1
temp[3,i] = W_OTC[3,:]@np.hstack([1, P[3,i], tempPi[:,i]]) #C_H
temp[4,i] = W_OTC[4,:]@np.hstack([1, P[4,i], tempPi[:,i]]) #C_P
temp[temp<0]=0
#Forecast for last period
#T = 100
tempPi[0,T] = W_factor[0,:]@np.hstack([1, Pi[0,T], Pi[0,T-1]]) #forecast for the ith period
tempPi[1,T] = W_factor[1,:]@np.hstack([1, Pi[1,T], Pi[1,T-1]])
tempPi[2,T] = W_factor[2,:]@np.stack([1, Pi[2,T], Pi[2,T-1]])
tempPi[tempPi<0]=0
inno_P_C = np.random.normal(0,np.sqrt(4.839),T+52+1)
inno_P_2 = np.random.normal(0,np.sqrt(16.58),T+52+1)
inno_P_1 = np.random.normal(0,np.sqrt(47.704),T+52+1)
inno_C_H = np.random.normal(0,np.sqrt(0.020227),T+52+1)
inno_C_P = np.random.normal(0,np.sqrt(0.0016579),T+52+1)
MA_OTC = np.array([0.35802, 0.45427, 0.34237, -0.25309, 0.45424])
SMA_OTC = np.array([0.11019, -0.027154, 0.0063966, 0.024557, -0.027257])
P[0,T] = W_OTC[0,:]@np.hstack([1, P[0,T-1], tempPi[:,T]]) + inno_P_C[T+52] + MA_OTC[0]*inno_P_C[T+51] + SMA_OTC[0]*inno_P_C[T] + MA_OTC[0]*SMA_OTC[0]*inno_P_C[T-1] #P_C
P[1,T] = W_OTC[1,:]@np.hstack([1, P[1,T-1], tempPi[:,T]]) + inno_P_2[T+52] + MA_OTC[1]*inno_P_2[T+51] + SMA_OTC[1]*inno_P_2[T] + MA_OTC[1]*SMA_OTC[1]*inno_P_2[T-1] #P_2
P[2,T] = W_OTC[2,:]@np.hstack([1, P[2,T-1], tempPi[:,T]]) + inno_P_1[T+52] + MA_OTC[2]*inno_P_1[T+51] + SMA_OTC[2]*inno_P_1[T] + MA_OTC[2]*SMA_OTC[2]*inno_P_1[T-1] #P_1
P[3,T] = W_OTC[3,:]@np.hstack([1, P[3,T-1], tempPi[:,T]]) + inno_C_H[T+52] + MA_OTC[3]*inno_C_H[T+51] + SMA_OTC[3]*inno_C_H[T] + MA_OTC[3]*SMA_OTC[3]*inno_C_H[T-1] #C_H
P[4,T] = W_OTC[4,:]@np.hstack([1, P[4,T-1], tempPi[:,T]]) + inno_C_P[T+52] + MA_OTC[4]*inno_C_P[T+51] + SMA_OTC[4]*inno_C_P[T] + MA_OTC[4]*SMA_OTC[4]*inno_C_P[T-1] #C_P
P[P<0]=0
#forward market
pF[0,T] = P[3,T]+rng.normal(0,P[3,T]/100) #mean and standard deviation
pF[pF<0]=0
# simulate demand: constant
avg_input = 1#106453
d = np.ones((1, Tmax))*35*avg_input
#calculate feature prices:
Feature = np.dot(L, P)
X_feature = np.vstack((np.ones((1,Tmax)), P, Feature))
X_feature = X_feature.T
#feature for procurement decisions
X_feature_buy = np.vstack((np.ones((1,Tmax)), P, pF,Feature))
X_feature_buy = X_feature_buy.Ttry:
m_ML2 = gp.Model()
m_ML2.reset(0) # since I am running model multiple times via looping, I reset the model for each run
m_ML2.resetParams()#define feature coefficients
alpha0 = m_ML2.addVars(T, name = 'alpha0', lb=-GRB.INFINITY)
alpha1 = m_ML2.addVars(T, name = 'alpha1', lb=-GRB.INFINITY)
alpha2 = m_ML2.addVars(T, name = 'alpha2', lb=-GRB.INFINITY)
alpha3 = m_ML2.addVars(T, name = 'alpha3', lb=-GRB.INFINITY)
alpha4 = m_ML2.addVars(T, name = 'alpha4', lb=-GRB.INFINITY)
alpha5 = m_ML2.addVars(T, name = 'alpha5', lb=-GRB.INFINITY)
alpha6 = m_ML2.addVars(T, name = 'alpha6', lb=-GRB.INFINITY)
alpha7 = m_ML2.addVars(T, name = 'alpha7', lb=-GRB.INFINITY)
alpha8 = m_ML2.addVars(T, name = 'alpha8', lb=-GRB.INFINITY)
alpha9 = m_ML2.addVars(T, name = 'alpha9', lb=-GRB.INFINITY)
alpha10 = m_ML2.addVars(T, name = 'alpha10', lb=-GRB.INFINITY)
alpha0_abs = m_ML2.addVars(T, name = 'alpha0_abs', lb=0)
alpha1_abs = m_ML2.addVars(T, name = 'alpha1_abs', lb=0)
alpha2_abs = m_ML2.addVars(T, name = 'alpha2_abs', lb=0)
alpha3_abs = m_ML2.addVars(T, name = 'alpha3_abs', lb=0)
alpha4_abs = m_ML2.addVars(T, name = 'alpha4_abs', lb=0)
alpha5_abs = m_ML2.addVars(T, name = 'alpha5_abs', lb=0)
alpha6_abs = m_ML2.addVars(T, name = 'alpha6_abs', lb=0)
alpha7_abs = m_ML2.addVars(T, name = 'alpha7_abs', lb=0)
alpha8_abs = m_ML2.addVars(T, name = 'alpha8_abs', lb=0)
alpha9_abs = m_ML2.addVars(T, name = 'alpha9_abs', lb=0)
alpha10_abs = m_ML2.addVars(T, name = 'alpha10_abs', lb=0)
beta0 = m_ML2.addVars(T, name = 'beta0', lb=-GRB.INFINITY)
beta1 = m_ML2.addVars(T, name = 'beta1', lb=-GRB.INFINITY)
beta2 = m_ML2.addVars(T, name = 'beta2', lb=-GRB.INFINITY)
beta3 = m_ML2.addVars(T, name = 'beta3', lb=-GRB.INFINITY)
beta4 = m_ML2.addVars(T, name = 'beta4', lb=-GRB.INFINITY)
beta5 = m_ML2.addVars(T, name = 'beta5', lb=-GRB.INFINITY)
beta6 = m_ML2.addVars(T, name = 'beta6', lb=-GRB.INFINITY)
beta7 = m_ML2.addVars(T, name = 'beta7', lb=-GRB.INFINITY)
beta8 = m_ML2.addVars(T, name = 'beta8', lb=-GRB.INFINITY)
beta9 = m_ML2.addVars(T, name = 'beta9', lb=-GRB.INFINITY)
beta10 = m_ML2.addVars(T, name = 'beta10', lb=-GRB.INFINITY)
beta0_abs = m_ML2.addVars(T, name = 'beta0_abs', lb=0)
beta1_abs = m_ML2.addVars(T, name = 'beta1_abs', lb=0)
beta2_abs = m_ML2.addVars(T, name = 'beta2_abs', lb=0)
beta3_abs = m_ML2.addVars(T, name = 'beta3_abs', lb=0)
beta4_abs = m_ML2.addVars(T, name = 'beta4_abs', lb=0)
beta5_abs = m_ML2.addVars(T, name = 'beta5_abs', lb=0)
beta6_abs = m_ML2.addVars(T, name = 'beta6_abs', lb=0)
beta7_abs = m_ML2.addVars(T, name = 'beta7_abs', lb=0)
beta8_abs = m_ML2.addVars(T, name = 'beta8_abs', lb=0)
beta9_abs = m_ML2.addVars(T, name = 'beta9_abs', lb=0)
beta10_abs = m_ML2.addVars(T, name = 'beta10_abs', lb=0)
gamma0 = m_ML2.addVars(T, name = 'gamma0', lb=-GRB.INFINITY)
gamma1 = m_ML2.addVars(T, name = 'gamma1', lb=-GRB.INFINITY)
gamma2 = m_ML2.addVars(T, name = 'gamma2', lb=-GRB.INFINITY)
gamma3 = m_ML2.addVars(T, name = 'gamma3', lb=-GRB.INFINITY)
gamma4 = m_ML2.addVars(T, name = 'gamma4', lb=-GRB.INFINITY)
gamma5 = m_ML2.addVars(T, name = 'gamma5', lb=-GRB.INFINITY)
gamma6 = m_ML2.addVars(T, name = 'gamma6', lb=-GRB.INFINITY)
gamma7 = m_ML2.addVars(T, name = 'gamma7', lb=-GRB.INFINITY)
gamma8 = m_ML2.addVars(T, name = 'gamma8', lb=-GRB.INFINITY)
gamma9 = m_ML2.addVars(T, name = 'gamma9', lb=-GRB.INFINITY)
gamma10 = m_ML2.addVars(T, name = 'gamma10', lb=-GRB.INFINITY)
gamma0_abs = m_ML2.addVars(T, name = 'gamma0_abs', lb=0)
gamma1_abs = m_ML2.addVars(T, name = 'gamma1_abs', lb=0)
gamma2_abs = m_ML2.addVars(T, name = 'gamma2_abs', lb=0)
gamma3_abs = m_ML2.addVars(T, name = 'gamma3_abs', lb=0)
gamma4_abs = m_ML2.addVars(T, name = 'gamma4_abs', lb=0)
gamma5_abs = m_ML2.addVars(T, name = 'gamma5_abs', lb=0)
gamma6_abs = m_ML2.addVars(T, name = 'gamma6_abs', lb=0)
gamma7_abs = m_ML2.addVars(T, name = 'gamma7_abs', lb=0)
gamma8_abs = m_ML2.addVars(T, name = 'gamma8_abs', lb=0)
gamma9_abs = m_ML2.addVars(T, name = 'gamma9_abs', lb=0)
gamma10_abs = m_ML2.addVars(T, name = 'gamma10_abs', lb=0)
#define hedging feature coefficients
zeta0 = m_ML2.addVars(T, name = 'zeta0', lb=-GRB.INFINITY)
zeta1 = m_ML2.addVars(T, name = 'zeta1', lb=-GRB.INFINITY)
zeta2 = m_ML2.addVars(T, name = 'zeta2', lb=-GRB.INFINITY)
zeta3 = m_ML2.addVars(T, name = 'zeta3', lb=-GRB.INFINITY)
zeta4 = m_ML2.addVars(T, name = 'zeta4', lb=-GRB.INFINITY)
zeta5 = m_ML2.addVars(T, name = 'zeta5', lb=-GRB.INFINITY)
zeta6 = m_ML2.addVars(T, name = 'zeta6', lb=-GRB.INFINITY)
zeta7 = m_ML2.addVars(T, name = 'zeta7', lb=-GRB.INFINITY)
zeta8 = m_ML2.addVars(T, name = 'zeta8', lb=-GRB.INFINITY)
zeta9 = m_ML2.addVars(T, name = 'zeta9', lb=-GRB.INFINITY)
zeta10 = m_ML2.addVars(T, name = 'zeta10', lb=-GRB.INFINITY)
zeta11 = m_ML2.addVars(T, name = 'zeta11', lb=-GRB.INFINITY)
zeta0_abs = m_ML2.addVars(T, name = 'zeta0_abs', lb=0)
zeta1_abs = m_ML2.addVars(T, name = 'zeta1_abs', lb=0)
zeta2_abs = m_ML2.addVars(T, name = 'zeta2_abs', lb=0)
zeta3_abs = m_ML2.addVars(T, name = 'zeta3_abs', lb=0)
zeta4_abs = m_ML2.addVars(T, name = 'zeta4_abs', lb=0)
zeta5_abs = m_ML2.addVars(T, name = 'zeta5_abs', lb=0)
zeta6_abs = m_ML2.addVars(T, name = 'zeta6_abs', lb=0)
zeta7_abs = m_ML2.addVars(T, name = 'zeta7_abs', lb=0)
zeta8_abs = m_ML2.addVars(T, name = 'zeta8_abs', lb=0)
zeta9_abs = m_ML2.addVars(T, name = 'zeta9_abs', lb=0)
zeta10_abs = m_ML2.addVars(T, name = 'zeta10_abs', lb=0)
zeta11_abs = m_ML2.addVars(T, name = 'zeta11_abs', lb=0)
#define procurement decisions
q_spot = m_ML2.addVars(T, vtype = GRB.BINARY, name = 'q_spot')
q_fwd = m_ML2.addVars(T, vtype = GRB.BINARY, name = 'q_fwd')
Zeta_threshold = m_ML2.addVars(T, name = 'Zeta_threshold')
# define auxiliary decision variables
x1_10 = m_ML2.addVars(T, name = 'x1_10', lb=-GRB.INFINITY)
x2_10 = m_ML2.addVars(T, name = 'x2_10', lb=-GRB.INFINITY)
y1_10 = m_ML2.addVars(T, name = 'y1_10', lb=-GRB.INFINITY)
y1_11 = m_ML2.addVars(T, name = 'y1_11', lb=-GRB.INFINITY)
y1_12 = m_ML2.addVars(T, name = 'y1_11', lb=-GRB.INFINITY)
y1_13 = m_ML2.addVars(T, name = 'y1_11', lb=-GRB.INFINITY)
y2_10 = m_ML2.addVars(T, name = 'y2_10', lb=-GRB.INFINITY)
y2_11 = m_ML2.addVars(T, name = 'y2_11', lb=-GRB.INFINITY)
y2_12 = m_ML2.addVars(T, name = 'y2_12', lb=-GRB.INFINITY)
x1_1 = m_ML2.addVars(T, name = 'x1_1')
x2_1 = m_ML2.addVars(T, name = 'x2_1')
y1_1 = m_ML2.addVars(T, name = 'y1_1')
y2_1 = m_ML2.addVars(T, name = 'y2_1')
z1_1 = m_ML2.addVars(T, name = 'z1_1')
z2_1 = m_ML2.addVars(T, name = 'z2_1')
s1_1 = m_ML2.addVars(T, name = 's1_1')
s2_1 = m_ML2.addVars(T, name = 's2_1')
#Define threholds
z_2_x1 = m_ML2.addVars(T, name = 'z_2_x1')
z_1_x1 = m_ML2.addVars(T, name = 'z_1_x1')
z_1_y1 = m_ML2.addVars(T, name = 'z_1_y1')
R_ML = m_ML2.addVars(T, name = 'R_ML')
#period 0
m_ML2.addConstr(s1_1[0] == w1[0])
m_ML2.addConstr(s2_1[0] == w2[0])
#constraint for coefficents
m_ML2.addConstr(alpha0_abs[0] == abs_(alpha0[0]))
m_ML2.addConstr(alpha1_abs[0] == abs_(alpha1[0]))
m_ML2.addConstr(alpha2_abs[0] == abs_(alpha2[0]))
m_ML2.addConstr(alpha3_abs[0] == abs_(alpha3[0]))
m_ML2.addConstr(alpha4_abs[0] == abs_(alpha4[0]))
m_ML2.addConstr(alpha5_abs[0] == abs_(alpha5[0]))
m_ML2.addConstr(alpha6_abs[0] == abs_(alpha6[0]))
m_ML2.addConstr(alpha7_abs[0] == abs_(alpha7[0]))
m_ML2.addConstr(alpha8_abs[0] == abs_(alpha8[0]))
m_ML2.addConstr(alpha9_abs[0] == abs_(alpha9[0]))
m_ML2.addConstr(alpha10_abs[0] == abs_(alpha10[0]))
m_ML2.addConstr(beta0_abs[0] == abs_(beta0[0]))
m_ML2.addConstr(beta1_abs[0] == abs_(beta1[0]))
m_ML2.addConstr(beta2_abs[0] == abs_(beta2[0]))
m_ML2.addConstr(beta3_abs[0] == abs_(beta3[0]))
m_ML2.addConstr(beta4_abs[0] == abs_(beta4[0]))
m_ML2.addConstr(beta5_abs[0] == abs_(beta5[0]))
m_ML2.addConstr(beta6_abs[0] == abs_(beta6[0]))
m_ML2.addConstr(beta7_abs[0] == abs_(beta7[0]))
m_ML2.addConstr(beta8_abs[0] == abs_(beta8[0]))
m_ML2.addConstr(beta9_abs[0] == abs_(beta9[0]))
m_ML2.addConstr(beta10_abs[0] == abs_(beta10[0]))
m_ML2.addConstr(gamma0_abs[0] == abs_(gamma0[0]))
m_ML2.addConstr(gamma1_abs[0] == abs_(gamma1[0]))
m_ML2.addConstr(gamma2_abs[0] == abs_(gamma2[0]))
m_ML2.addConstr(gamma3_abs[0] == abs_(gamma3[0]))
m_ML2.addConstr(gamma4_abs[0] == abs_(gamma4[0]))
m_ML2.addConstr(gamma5_abs[0] == abs_(gamma5[0]))
m_ML2.addConstr(gamma6_abs[0] == abs_(gamma6[0]))
m_ML2.addConstr(gamma7_abs[0] == abs_(gamma7[0]))
m_ML2.addConstr(gamma8_abs[0] == abs_(gamma8[0]))
m_ML2.addConstr(gamma9_abs[0] == abs_(gamma9[0]))
m_ML2.addConstr(gamma10_abs[0] == abs_(gamma10[0]))
m_ML2.addConstr(zeta0_abs[0] >= zeta0[0])
m_ML2.addConstr(zeta1_abs[0] >= zeta1[0])
m_ML2.addConstr(zeta2_abs[0] >= zeta2[0])
m_ML2.addConstr(zeta3_abs[0] >= zeta3[0])
m_ML2.addConstr(zeta4_abs[0] >= zeta4[0])
m_ML2.addConstr(zeta5_abs[0] >= zeta5[0])
m_ML2.addConstr(zeta6_abs[0] >= zeta6[0])
m_ML2.addConstr(zeta7_abs[0] >= zeta7[0])
m_ML2.addConstr(zeta8_abs[0] >= zeta8[0])
m_ML2.addConstr(zeta9_abs[0] >= zeta9[0])
m_ML2.addConstr(zeta10_abs[0] >= zeta10[0])
m_ML2.addConstr(zeta11_abs[0] >= zeta11[0])
m_ML2.addConstr(zeta0_abs[0] >= -zeta0[0])
m_ML2.addConstr(zeta1_abs[0] >= -zeta1[0])
m_ML2.addConstr(zeta2_abs[0] >= -zeta2[0])
m_ML2.addConstr(zeta3_abs[0] >= -zeta3[0])
m_ML2.addConstr(zeta4_abs[0] >= -zeta4[0])
m_ML2.addConstr(zeta5_abs[0] >= -zeta5[0])
m_ML2.addConstr(zeta6_abs[0] >= -zeta6[0])
m_ML2.addConstr(zeta7_abs[0] >= -zeta7[0])
m_ML2.addConstr(zeta8_abs[0] >= -zeta8[0])
m_ML2.addConstr(zeta9_abs[0] >= -zeta9[0])
m_ML2.addConstr(zeta10_abs[0] >= -zeta10[0])
m_ML2.addConstr(zeta11_abs[0] >= -zeta11[0])
#purchasing decision first
m_ML2.addConstr(q_spot[0] == 1)
m_ML2.addConstr(Zeta_threshold[0] == zeta0[0]*X_feature_buy[0,0]+zeta1[0]*X_feature_buy[0,1]+zeta2[0]*X_feature_buy[0,2]+zeta3[0]*X_feature_buy[0,3]\
+zeta4[0]*X_feature_buy[0,4]+zeta5[0]*X_feature_buy[0,5]+zeta6[0]*X_feature_buy[0,6]+zeta7[0]*X_feature_buy[0,7]\
+zeta8[0]*X_feature_buy[0,8]+zeta9[0]*X_feature_buy[0,9]+zeta10[0]*X_feature_buy[0,10]+zeta11[0]*X_feature_buy[0,11])
m_ML2.addConstr(Zeta_threshold[0] >= pF[0,0] + eps - M*(1-q_fwd[0])) # if the forward price pF is greater than the threshold, don't buy
m_ML2.addConstr(Zeta_threshold[0] <= pF[0,0] + M*q_fwd[0]) # pF[0] here is to buy corn&soy for the next eta period
#Threshold Upperbar: Holding vs. Selling
#\bar{z}_2
if temp[0,0]+temp[4,0]-P[1,0]-P[3,0]<=0:
m_ML2.addConstr(z_2_x1[0] == 0)
else:
m_ML2.addConstr(z_2_x1[0] == alpha0[0]*X_feature[0,0]+alpha1[0]*X_feature[0,1]+alpha2[0]*X_feature[0,2]
+alpha3[0]*X_feature[0,3]+alpha4[0]*X_feature[0,4]+alpha5[0]*X_feature[0,5]+alpha6[0]*X_feature[0,6]
+alpha7[0]*X_feature[0,7]+alpha8[0]*X_feature[0,8]+alpha9[0]*X_feature[0,9]+alpha10[0]*X_feature[0,10])
#\bar{z}_1
if P[2,0] <= min(temp[1,0],temp[0,0]+temp[4,0])-P[3,0]:
m_ML2.addConstr(z_1_x1[0] == beta0[0]*X_feature[0,0]+beta1[0]*X_feature[0,1]+beta2[0]*X_feature[0,2]+beta3[0]*X_feature[0,3]
+beta4[0]*X_feature[0,4]+beta5[0]*X_feature[0,5]+beta6[0]*X_feature[0,6]+beta7[0]*X_feature[0,7]
+beta8[0]*X_feature[0,8]+beta9[0]*X_feature[0,9]+beta10[0]*X_feature[0,10])
elif P[2,0]>=max(temp[1,0],temp[0,0]+temp[4,0])-P[3,0]:
m_ML2.addConstr(z_1_x1[0] == 0)
else:
m_ML2.addConstr(z_1_x1[0] == beta0[0]*X_feature[0,0]+beta1[0]*X_feature[0,1]+beta2[0]*X_feature[0,2]+beta3[0]*X_feature[0,3]
+beta4[0]*X_feature[0,4]+beta5[0]*X_feature[0,5]+beta6[0]*X_feature[0,6]+beta7[0]*X_feature[0,7]
+beta8[0]*X_feature[0,8]+beta9[0]*X_feature[0,9]+beta10[0]*X_feature[0,10])
#Underbar: Holding vs. Fulfilling
#\underbar{z}_1
if (1-alpha)* P[0,0]+P[4,0]<=min(temp[1,0],temp[0,0]+temp[4,0])-P[3,0]:
m_ML2.addConstr(z_1_y1[0] == gamma0[0]*X_feature[0,0]+gamma1[0]*X_feature[0,1]+gamma2[0]*X_feature[0,2]+gamma3[0]*X_feature[0,3]
+gamma4[0]*X_feature[0,4]+gamma5[0]*X_feature[0,5]+gamma6[0]*X_feature[0,6]+gamma7[0]*X_feature[0,7]
+gamma8[0]*X_feature[0,8]+gamma9[0]*X_feature[0,9]+gamma10[0]*X_feature[0,10])
elif (1-alpha)* P[0,0]+P[4,0]>=max(temp[1,0],temp[0,0]+temp[4,0])-P[3,0]:
m_ML2.addConstr(z_1_y1[0] == 0)
else:
m_ML2.addConstr(z_1_y1[0] == gamma0[0]*X_feature[0,0]+gamma1[0]*X_feature[0,1]+gamma2[0]*X_feature[0,2]+gamma3[0]*X_feature[0,3]
+gamma4[0]*X_feature[0,4]+gamma5[0]*X_feature[0,5]+gamma6[0]*X_feature[0,6]+gamma7[0]*X_feature[0,7]
+gamma8[0]*X_feature[0,8]+gamma9[0]*X_feature[0,9]+gamma10[0]*X_feature[0,10])
#decision variables 5 cases
if P[0,0]+P[4,0]<P[1,0]:
if P[2,0]>(1-alpha)*P[0,0]+P[4,0]:
m_ML2.addConstr(y1_1[0] == 0.)
m_ML2.addConstr(y2_1[0] == 0.)
m_ML2.addConstr(x1_10[0] == s1_1[0]-z_1_x1[0]) #auxiliary var to x1_1
m_ML2.addConstr(x1_1[0] == gp.max_(0., x1_10[0])) #the general constraints methods accept Var objects only, not sum/minus b/w two vars.
m_ML2.addConstr(x2_1[0] == s2_1[0])
else:
m_ML2.addConstr(y1_10[0] == s1_1[0]-z_1_y1[0])
m_ML2.addConstr(y1_11[0] == gp.max_(y1_10[0],0.))
m_ML2.addConstr(y1_1[0] == gp.min_(y1_10[0],q))
m_ML2.addConstr(x1_10[0] == s1_1[0]-z_1_x1[0]-q)
m_ML2.addConstr(x1_1[0] == gp.max_(0., x1_10[0]))
m_ML2.addConstr(x2_1[0] == s2_1[0])
else:
if P[2,0]>(1-alpha)*P[0,0]+P[4,0]:#case 2a
m_ML2.addConstr(y1_1[0] == 0.)
m_ML2.addConstr(y2_1[0] == gp.min_(s2_1[0], q))
m_ML2.addConstr(x1_10[0] == s1_1[0]-z_1_x1[0])
m_ML2.addConstr(x1_1[0] == gp.max_(0., x1_10[0]))
m_ML2.addConstr(x2_10[0] == s2_1[0]-q-z_2_x1[0])
m_ML2.addConstr(x2_1[0] == gp.max_(0., x2_10[0]))
elif P[2,0]>P[1,0]-alpha* P[0,0]: #case 2b
m_ML2.addConstr(y1_10[0] == s1_1[0]-z_1_y1[0])
m_ML2.addConstr(y1_11[0] == gp.max_(y1_10[0],0.))
m_ML2.addConstr(y1_12[0] == q-s2_1[0])
m_ML2.addConstr(y1_13[0] == gp.max_(y1_12[0],0.))
m_ML2.addConstr(y1_1[0] == gp.min_(y1_11[0], y1_13[0]))
m_ML2.addConstr(y2_1[0] == gp.min_(s2_1[0], q))
m_ML2.addConstr(x1_10[0] == s1_1[0]-z_1_x1[0]-y1_1[0])
m_ML2.addConstr(x1_1[0] == gp.max_(0., x1_10[0]))
m_ML2.addConstr(x2_10[0] == s2_1[0]-q-z_2_x1[0])
m_ML2.addConstr(x2_1[0] == gp.max_(0., x2_10[0]))
else: #case 2c
m_ML2.addConstr(y1_10[0] == s1_1[0]-z_1_y1[0])
m_ML2.addConstr(y1_11[0] == gp.max_(y1_10[0], 0.))
m_ML2.addConstr(y1_1[0] == gp.min_(q, y1_11[0]))
m_ML2.addConstr(y2_10[0] == q -s1_1[0]+z_1_y1[0])
m_ML2.addConstr(y2_11[0] == gp.min_(q, y2_10[0]))
m_ML2.addConstr(y2_12[0] == gp.max_(0., y2_11[0]))
m_ML2.addConstr(y2_1[0] == gp.min_(s2_1[0], y2_12[0]))
m_ML2.addConstr(z1_1[0] == s1_1[0]-x1_1[0]-y1_1[0])
m_ML2.addConstr(z2_1[0] == s2_1[0]-x2_1[0]-y2_1[0])
#reward
m_ML2.addConstr(R_ML[0] == -P[3,0]*q_spot[0]-pF[0,0]*q_fwd[0]+ P[0,0]*(y2_1[0]+(1-alpha)*y1_1[0])+P[1,0]*x2_1[0]+P[2,0]*x1_1[0]-P[3,0]*(z1_1[0]+z2_1[0])-P[4,0]*(q-y1_1[0]-y2_1[0])-P[3,0]*(w1[0]+w2[0])
- alpha0_abs[0] - alpha1_abs[0] - alpha2_abs[0]- alpha3_abs[0]- alpha4_abs[0]- alpha5_abs[0]- alpha6_abs[0]
- alpha7_abs[0]- alpha8_abs[0]- alpha9_abs[0]- alpha10_abs[0]- beta0_abs[0] - beta1_abs[0] - beta2_abs[0]
- beta3_abs[0] - beta4_abs[0] - beta5_abs[0] - beta6_abs[0] - beta7_abs[0] - beta8_abs[0] - beta9_abs[0]
- beta10_abs[0]- gamma0_abs[0] - gamma1_abs[0] - gamma2_abs[0] - gamma3_abs[0] - gamma4_abs[0] - gamma5_abs[0]
- gamma6_abs[0] - gamma7_abs[0] - gamma8_abs[0] - gamma9_abs[0] - gamma10_abs[0] - zeta0_abs[0] - zeta1_abs[0] - zeta2_abs[0] - zeta3_abs[0] - zeta4_abs[0]
- zeta5_abs[0] - zeta6_abs[0] - zeta7_abs[0] - zeta8_abs[0] - zeta9_abs[0] - zeta10_abs[0] - zeta11_abs[0])
#period i
for i in range(1,T):
m_ML2.addConstr(s1_1[i] == w1[i])
m_ML2.addConstr(s2_1[i] == w2[i]+z1_1[i-1]+z2_1[i-1])
m_ML2.addConstr(alpha0_abs[i] == abs_(alpha0[i]))
m_ML2.addConstr(alpha1_abs[i] == abs_(alpha1[i]))
m_ML2.addConstr(alpha2_abs[i] == abs_(alpha2[i]))
m_ML2.addConstr(alpha3_abs[i] == abs_(alpha3[i]))
m_ML2.addConstr(alpha4_abs[i] == abs_(alpha4[i]))
m_ML2.addConstr(alpha5_abs[i] == abs_(alpha5[i]))
m_ML2.addConstr(alpha6_abs[i] == abs_(alpha6[i]))
m_ML2.addConstr(alpha7_abs[i] == abs_(alpha7[i]))
m_ML2.addConstr(alpha8_abs[i] == abs_(alpha8[i]))
m_ML2.addConstr(alpha9_abs[i] == abs_(alpha9[i]))
m_ML2.addConstr(alpha10_abs[i] == abs_(alpha10[i]))
m_ML2.addConstr(beta0_abs[i] == abs_(beta0[i]))
m_ML2.addConstr(beta1_abs[i] == abs_(beta1[i]))
m_ML2.addConstr(beta2_abs[i] == abs_(beta2[i]))
m_ML2.addConstr(beta3_abs[i] == abs_(beta3[i]))
m_ML2.addConstr(beta4_abs[i] == abs_(beta4[i]))
m_ML2.addConstr(beta5_abs[i] == abs_(beta5[i]))
m_ML2.addConstr(beta6_abs[i] == abs_(beta6[i]))
m_ML2.addConstr(beta7_abs[i] == abs_(beta7[i]))
m_ML2.addConstr(beta8_abs[i] == abs_(beta8[i]))
m_ML2.addConstr(beta9_abs[i] == abs_(beta9[i]))
m_ML2.addConstr(beta10_abs[i] == abs_(beta10[i]))
m_ML2.addConstr(gamma0_abs[i] == abs_(gamma0[i]))
m_ML2.addConstr(gamma1_abs[i] == abs_(gamma1[i]))
m_ML2.addConstr(gamma2_abs[i] == abs_(gamma2[i]))
m_ML2.addConstr(gamma3_abs[i] == abs_(gamma3[i]))
m_ML2.addConstr(gamma4_abs[i] == abs_(gamma4[i]))
m_ML2.addConstr(gamma5_abs[i] == abs_(gamma5[i]))
m_ML2.addConstr(gamma6_abs[i] == abs_(gamma6[i]))
m_ML2.addConstr(gamma7_abs[i] == abs_(gamma7[i]))
m_ML2.addConstr(gamma8_abs[i] == abs_(gamma8[i]))
m_ML2.addConstr(gamma9_abs[i] == abs_(gamma9[i]))
m_ML2.addConstr(gamma10_abs[i] == abs_(gamma10[i]))
m_ML2.addConstr(zeta0_abs[i] >= zeta0[i])
m_ML2.addConstr(zeta1_abs[i] >= zeta1[i])
m_ML2.addConstr(zeta2_abs[i] >= zeta2[i])
m_ML2.addConstr(zeta3_abs[i] >= zeta3[i])
m_ML2.addConstr(zeta4_abs[i] >= zeta4[i])
m_ML2.addConstr(zeta5_abs[i] >= zeta5[i])
m_ML2.addConstr(zeta6_abs[i] >= zeta6[i])
m_ML2.addConstr(zeta7_abs[i] >= zeta7[i])
m_ML2.addConstr(zeta8_abs[i] >= zeta8[i])
m_ML2.addConstr(zeta9_abs[i] >= zeta9[i])
m_ML2.addConstr(zeta10_abs[i] >= zeta10[i])
m_ML2.addConstr(zeta11_abs[i] >= zeta11[i])
m_ML2.addConstr(zeta0_abs[i] >= -zeta0[i])
m_ML2.addConstr(zeta1_abs[i] >= -zeta1[i])
m_ML2.addConstr(zeta2_abs[i] >= -zeta2[i])
m_ML2.addConstr(zeta3_abs[i] >= -zeta3[i])
m_ML2.addConstr(zeta4_abs[i] >= -zeta4[i])
m_ML2.addConstr(zeta5_abs[i] >= -zeta5[i])
m_ML2.addConstr(zeta6_abs[i] >= -zeta6[i])
m_ML2.addConstr(zeta7_abs[i] >= -zeta7[i])
m_ML2.addConstr(zeta8_abs[i] >= -zeta8[i])
m_ML2.addConstr(zeta9_abs[i] >= -zeta9[i])
m_ML2.addConstr(zeta10_abs[i] >= -zeta10[i])
m_ML2.addConstr(zeta11_abs[i] >= -zeta11[i])
#procurement decision
#here, we write the optimal decision embedded in the constraint, we can also use if condition to make optimal decision
m_ML2.addConstr(Zeta_threshold[i] == zeta0[i]*X_feature_buy[i,0]+zeta1[i]*X_feature_buy[i,1]+zeta2[i]*X_feature_buy[i,2]+zeta3[i]*X_feature_buy[i,3]\
+zeta4[i]*X_feature_buy[i,4]+zeta5[i]*X_feature_buy[i,5]+zeta6[i]*X_feature_buy[i,6]+zeta7[i]*X_feature_buy[i,7]\
+zeta8[i]*X_feature_buy[i,8]+zeta9[i]*X_feature_buy[i,9]+zeta10[i]*X_feature_buy[i,10]+zeta11[i]*X_feature_buy[i,11])
m_ML2.addConstr(q_spot[i] + q_fwd[i-1] == 1)
#Here is the if statement
m_ML2.addConstr(Zeta_threshold[i] >= pF[0,i] + eps - M*(1-q_fwd[i])) # if the forward price pF is greater than the threshold, don't buy
m_ML2.addConstr(Zeta_threshold[i] <= pF[0,i] + M*q_fwd[i]) # pF[i] here is to buy corn&soy for the next eta period
#Threshold Upperbar: Holding vs. Selling
#\bar{z}_2
if temp[0,i]+temp[4,i]-P[1,i]-P[3,i]<=0:
m_ML2.addConstr(z_2_x1[i] == 0)
else:
m_ML2.addConstr(z_2_x1[i] == alpha0[i]*X_feature[i,0]+alpha1[i]*X_feature[i,1]+alpha2[i]*X_feature[i,2]+alpha3[i]*X_feature[i,3]
+alpha4[i]*X_feature[i,4]+alpha5[i]*X_feature[i,5]+alpha6[i]*X_feature[i,6]+alpha7[i]*X_feature[i,7]
+alpha8[i]*X_feature[i,8]+alpha9[i]*X_feature[i,9]+alpha10[i]*X_feature[i,10])
#\bar{z}_1
if P[2,i] <= min(temp[1,i],temp[0,i]+temp[4,i])-P[3,i]:
m_ML2.addConstr(z_1_x1[i] == beta0[i]*X_feature[i,0]+beta1[i]*X_feature[i,1]+beta2[i]*X_feature[i,2]+beta3[i]*X_feature[i,3]
+beta4[i]*X_feature[i,4]+beta5[i]*X_feature[i,5]+beta6[i]*X_feature[i,6]+beta7[i]*X_feature[i,7]
+beta8[i]*X_feature[i,8]+beta9[i]*X_feature[i,9]+beta10[i]*X_feature[i,10])
elif P[2,i]>=max(temp[1,i],temp[0,i]+temp[4,i])-P[3,i]:
m_ML2.addConstr(z_1_x1[i] == 0)
else:
m_ML2.addConstr(z_1_x1[i] == beta0[i]*X_feature[i,0]+beta1[i]*X_feature[i,1]+beta2[i]*X_feature[i,2]+beta3[i]*X_feature[i,3]
+beta4[i]*X_feature[i,4]+beta5[i]*X_feature[i,5]+beta6[i]*X_feature[i,6]+beta7[i]*X_feature[i,7]
+beta8[i]*X_feature[i,8]+beta9[i]*X_feature[i,9]+beta10[i]*X_feature[i,10])
#Underbar: Holding vs. Fulfilling
#\underbar{z}_1
if (1-alpha)* P[0,i]+P[4,i]<=min(temp[1,i],temp[0,i]+temp[4,i])-P[3,i]:
m_ML2.addConstr(z_1_y1[i] == gamma0[i]*X_feature[i,0]+gamma1[i]*X_feature[i,1]+gamma2[i]*X_feature[i,2]+gamma3[i]*X_feature[i,3]
+gamma4[i]*X_feature[i,4]+gamma5[i]*X_feature[i,5]+gamma6[i]*X_feature[i,6]+gamma7[i]*X_feature[i,7]
+gamma8[i]*X_feature[i,8]+gamma9[i]*X_feature[i,9]+gamma10[i]*X_feature[i,10])
elif (1-alpha)* P[0,i]+P[4,i]>=max(temp[1,i],temp[0,i]+temp[4,i])-P[3,i]:
m_ML2.addConstr(z_1_y1[i] == 0)
else:
m_ML2.addConstr(z_1_y1[i] == gamma0[i]*X_feature[i,0]+gamma1[i]*X_feature[i,1]+gamma2[i]*X_feature[i,2]+gamma3[i]*X_feature[i,3]
+gamma4[i]*X_feature[i,4]+gamma5[i]*X_feature[i,5]+gamma6[i]*X_feature[i,6]+gamma7[i]*X_feature[i,7]
+gamma8[i]*X_feature[i,8]+gamma9[i]*X_feature[i,9]+gamma10[i]*X_feature[i,10])
#decision variables 5 cases
if P[0,i]+P[4,i]<P[1,i]:
if P[2,i]>(1-alpha)*P[0,i]+P[4,i]:
m_ML2.addConstr(y1_1[i] == 0.)
m_ML2.addConstr(y2_1[i] == 0.)
m_ML2.addConstr(x1_10[i] == s1_1[i]-z_1_x1[i]) #auxiliary var to x1_1
m_ML2.addConstr(x1_1[i] == gp.max_(0., x1_10[i])) #the general constraints methods accept Var objects only, not sum/minus b/w two vars.
m_ML2.addConstr(x2_1[i] == s2_1[i])
else:
m_ML2.addConstr(y1_10[i] == s1_1[i]-z_1_y1[i])
m_ML2.addConstr(y1_11[i] == gp.max_(y1_10[i],0.))
m_ML2.addConstr(y1_1[i] == gp.min_(y1_10[i],q))
m_ML2.addConstr(x1_10[i] == s1_1[i]-z_1_x1[i]-q)
m_ML2.addConstr(x1_1[i] == gp.max_(0., x1_10[i]))
m_ML2.addConstr(x2_1[i] == s2_1[i])
else:
if P[2,i]>(1-alpha)*P[0,i]+P[4,i]:#case 2a
m_ML2.addConstr(y1_1[i] == 0.)
m_ML2.addConstr(y2_1[i] == gp.min_(s2_1[i], q))
m_ML2.addConstr(x1_10[i] == s1_1[i]-z_1_x1[i])
m_ML2.addConstr(x1_1[i] == gp.max_(0., x1_10[i]))
m_ML2.addConstr(x2_10[i] == s2_1[i]-q-z_2_x1[i])
m_ML2.addConstr(x2_1[i] == gp.max_(0., x2_10[i]))
elif P[2,i]>P[1,i]-alpha* P[0,i]: #case 2b
m_ML2.addConstr(y1_10[i] == s1_1[i]-z_1_y1[i])
m_ML2.addConstr(y1_11[i] == gp.max_(y1_10[i],0.))
m_ML2.addConstr(y1_12[i] == q-s2_1[i])
m_ML2.addConstr(y1_13[i] == gp.max_(y1_12[i],0.))
m_ML2.addConstr(y1_1[i] == gp.min_(y1_11[i], y1_13[i]))
m_ML2.addConstr(y2_1[i] == gp.min_(s2_1[i], q))
m_ML2.addConstr(x1_10[i] == s1_1[i]-z_1_x1[i]-y1_1[i])
m_ML2.addConstr(x1_1[i] == gp.max_(0., x1_10[i]))
m_ML2.addConstr(x2_10[i] == s2_1[i]-q-z_2_x1[i])
m_ML2.addConstr(x2_1[i] == gp.max_(0., x2_10[i]))
else: #case 2c
m_ML2.addConstr(y1_10[i] == s1_1[i]-z_1_y1[i])
m_ML2.addConstr(y1_11[i] == gp.max_(y1_10[i], 0.))
m_ML2.addConstr(y1_1[i] == gp.min_(q, y1_11[i]))
m_ML2.addConstr(y2_10[i] == q -s1_1[i]+z_1_y1[i])
m_ML2.addConstr(y2_11[i] == gp.min_(q, y2_10[i]))
m_ML2.addConstr(y2_12[i] == gp.max_(0., y2_11[i]))
m_ML2.addConstr(y2_1[i] == gp.min_(s2_1[i], y2_12[i]))
m_ML2.addConstr(z1_1[i] == s1_1[i]-x1_1[i]-y1_1[i])
m_ML2.addConstr(z2_1[i] == s2_1[i]-x2_1[i]-y2_1[i])
#reward
m_ML2.addConstr(R_ML[i] == -P[3,i]*q_spot[i]-pF[0,i]*q_fwd[i] + P[0,i]*(y2_1[i]+(1-alpha)*y1_1[i])+P[1,i]*x2_1[i]+P[2,i]*x1_1[i]-P[3,i]*(z1_1[i]+z2_1[i])-P[4,i]*(q-y1_1[i]-y2_1[i])-P[3,i]*(w1[i]+w2[i])
- alpha0_abs[i] - alpha1_abs[i] - alpha2_abs[i]- alpha3_abs[i]- alpha4_abs[i]- alpha5_abs[i]- alpha6_abs[i]
- alpha7_abs[i]- alpha8_abs[i]- alpha9_abs[i]- alpha10_abs[i]- beta0_abs[i] - beta1_abs[i] - beta2_abs[i]
- beta3_abs[i] - beta4_abs[i] - beta5_abs[i] - beta6_abs[i] - beta7_abs[i] - beta8_abs[i] - beta9_abs[i]
- beta10_abs[i]- gamma0_abs[i] - gamma1_abs[i] - gamma2_abs[i] - gamma3_abs[i] - gamma4_abs[i] - gamma5_abs[i]
- gamma6_abs[i] - gamma7_abs[i] - gamma8_abs[i] - gamma9_abs[i] - gamma10_abs[i] - zeta0_abs[i] - zeta1_abs[i] - zeta2_abs[i] - zeta3_abs[i] - zeta4_abs[i]
- zeta5_abs[i] - zeta6_abs[i] - zeta7_abs[i] - zeta8_abs[i] - zeta9_abs[i] - zeta10_abs[i] - zeta11_abs[i])
m_ML2.setObjective(R_ML.sum(), GRB.MAXIMIZE)
# Optimize model
m_ML2.optimize()
V2[n] = m_ML2.objVal
for v in m_ML2.getVars():
print('%s %g' % (v.varName, v.x))
print('Obj: %g' % m_ML2.objVal)
except gp.GurobiError as e:
print('Error code ' + str(e.errno) + ': ' + str(e))
except AttributeError:
print('Encountered an attribute error')V_ML[T-1] = np.mean(V2)
# Plot
import numpy as np
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
V_MLSort = np.sort(V_ML).T
x = np.array([np.linspace(2,Tmax, Tmax-1)]).T
plt.plot(x,V_MLSort, label = "ML")
plt.title("Performance of Each Policy")
plt.ylabel("Accumulated Profit")
plt.xlabel("Time Horizon (T)")
ax.legend()
plt.show()
0 -
You are solving many models with randomly generated data. It is possible that for some data, the model is harder than it is on average. It may have to do with numerical issues arising from random coefficients or a particular data set is simply really hard to solve.
You are currently solving up to the default gap of 0.01% which is very small given that you objective function values are > 1e+07. Do you really need such a tight MIPGap? If not, you should try setting the desired MIPGap to, e.g., 0.1%
m_ML2.setParam("MIPGap",0.001)This may already solve your issue however it is not guaranteed that it does. There may be random data sets which make your model hard to solve up to 0.1% gap as well. You could try to improve the numerics of your model a bit. This may help, cf. Guidelines for Numerical Issues.
What you also could try is to write every model as an LP/MPS file via the write method and have a closer look at a model which takes a long time to solve. Once you identified one such case, you can have a closer look at this particular models numerics and maybe find parameters to improve its performance. You can then try whether the parameters you found also work in general for your whole loop.
Best regards,
Jaromił0 -
Thank you so much, Jaromił! I tried all of your suggestions, and they are very helpful.
I still have one question:
Adjusting the MIPGap is very useful; I changed the gap to 0.2%. But still, there are some instances where the gap is very hard to narrow down (It took 10 hours for one instance to converge to a 0.7% gap.) I wonder if I can set for such instance to a larger gap, e.g., set a 0.7% gap for this kind of harder-to-solve instance (the one that takes more than 10 hours to solve), but keep a 0.2% gap for the rest of the instances?
Really appreciated!
Best,
Ye
0 -
Hi Ye,
Changing the MIPGap parameter "on the fly" is not supported. Still, there is a solution to your situation.
You can implement a callback, where you query the current best solution (\(z_P\)) and best bound (\(z_D\)). You check the current runtime in the callback and if the runtime is too large, you can compute the MIPGap value manually with the use of the best solution and best bound values. You can then manually terminate the optimization process from inside the callback function, once some termination criterion of your choice has been met.
A simplistic code snippet might look something like
def mycallback(model, where):
if where == GRB.Callback.MIPNODE:
runtime = model.cbGet(GRB.Callback.RUNTIME)
bestobj = model.cbGet(GRB.Callback.MIPSOL_OBJBST)
bestbnd = model.cbGet(GRB.Callback.MIPSOL_OBJBND)
# compute mipgap as described in the documentation
# add some termination criteria based on runtime + mipgap
# and terminate model if one of the criteria is metAn alternative would be to query the log message in a MESSAGE callback and extract the MIPGap value computed by Gurobi. Then based on this value together with some runtime criterion, terminate the optimization from inside the MESSAGE callback.
Best regards,
Jaromił0 -
Hi Jaromił,
Thanks a lot for your solution! However, when I implemented the code, sometimes I have the following message: Encountered an attribute error. I don't know if I made a mistake somewhere:
I put the following code before I define my model:
from gurobipy import *
def mycallback(model, where):
if where == GRB.Callback.MIPNODE:
runtime = model.cbGet(GRB.Callback.RUNTIME)
bestobj = model.cbGet(GRB.Callback.MIPSOL_OBJBST)
bestbnd = model.cbGet(GRB.Callback.MIPSOL_OBJBND)
gap = (bestbnd-bestobj)/bestobj# compute mipgap as described in the documentation
if runtime >= 10 and gap <= 0.6:
model.terminate()
[...my model]
m_ML2.optimize(mycallback)Then I have the following message which I didn't see when I didn't use the callback:
Thank you so much!
Best,
Ye
0 -
Hi Ye,
Your model is infeasible or unbounded. This means that you cannot access the solution attributes. Thus you get an attribute error. I guess it is similar to the issue we discussed in this thread before.
You should also check whether the bestbnd and bestobj you get are finite values to avoid numerical errors. To check whether a feasible point is available, you can check whether MIPNODE_SOLCNT callback attribute is > 0.
Best regards,
Jaromił0 -
Hi Jaromił,
Thank you so much! Yes, you are right, I used the iis.ilp file to find out that I forgot to relax the bound to be negative infinite. One minor issue is that iis.ilp file cannot overwrite, each time I need to rename the file if I want to produce the most updated ilp file when I run the code.
Also, it seems that the Callback code I am writing doesn't work. Under the condition
if runtime >= 200 and gap <= 0.6:
model.terminate()#My model still keeps running even if those two conditions have been satisfied (see picture below).
My code looks like this:
from gurobipy import *
def mycallback(model, where):
if where == GRB.Callback.MIPNODE:
count = model.cbGet(GRB.Callback.MIPNODE_SOLCNT)
runtime = model.cbGet(GRB.Callback.RUNTIME)
bestobj = model.cbGet(GRB.Callback.MIPSOL_OBJBST)
bestbnd = model.cbGet(GRB.Callback.MIPSOL_OBJBND)
gap = (bestbnd-bestobj)/bestobj# compute mipgap as described in the documentation
print("Count",count)
if runtime >= 200 and gap <= 0.6:
model.terminate()# add some termination criteria based on runtime + mipgap
# and terminate model if one of the criteria is met
[then I define my model...]Thank you so much!
Best,
Ye
0 -
Hi Ye,
You are in the MIPNODE callback, but you are trying to access to MIPSOL_OBJBST and MIPSOL_OBJBND values. This won't work. You have to access MIPNODE_OBJBST and MIPNODE_OBJBND.
I just noticed that I made the mistake originally, so this one is on me :)
Just to be sure, you should take the absolute values of the nominator and denominator when computing the MIPGap value, see MIPGap. There are also many special cases when computing the MIPGap such as "what to do when the denominator is 0", but I don't think that these apply in your case.
Best regards,
Jaromił0 -
Hi Jaromił,
Thanks for responding! I tried the following but still didn't work
from gurobipy import *
def mycallback(model, where):
if where == GRB.Callback.MIPNODE:
count = model.cbGet(GRB.Callback.MIPNODE_SOLCNT)
runtime = model.cbGet(GRB.Callback.RUNTIME)
if where == GRB.Callback.MIPSOL_OBJBST:
bestobj = model.cbGet(GRB.Callback.MIPSOL_OBJBST)
if where == GRB.Callback.MIPSOL_OBJBND:
bestbnd = model.cbGet(GRB.Callback.MIPSOL_OBJBND)
gap = gp.abs_(bestbnd-bestobj)/gp.abs_(bestobj)# compute mipgap as described in the documentation
if runtime >= 200 and gap <= 0.6:
model.terminate()# add some termination criteria based on runtime + mipgap
[my model...]Could you help me write an example? In each run, the model termination criteria are the runtime is >200s and the mipgap is less than 0.6%.
I check all the documents related to MIPGap , MIPNODE_SOLCNT callback attribute and other related ones, but still didn't find a suitable example that works for me.
Thank you so much!
Best,
Ye
0
Please sign in to leave a comment.
Comments
20 comments