Good day to you. From a previous model with cplex (MIP 1), I use M0 as a parameter.
Total_P = 8
M0 = np.array([[0,1,0,1,0,0,0,2]]).reshape(Total_P,1)
Then, I tried to develop a different model (MIP 2) with adjusted objectives and constraints from MIP2, however, the M0 now becomes a decision variable and use for several constraints (#constraint 38). So, I created a code as follows:
Total_P = 8
M0 = np.empty((Total_P,1), dtype = object)
for l in range(0,Total_P):
M0[l][0]= mdl.integer_var(name='M0'+ str(',') + str(l+1))
for l in range(0,Total_P):
if M0[l][0] is not None:
mdl.add_constraint((M0[1][0] + M0[5][0]) == 1)
mdl.add_constraint((M0[3][0] + M0[6][0]) == 1)
mdl.add_constraint((M0[0][0] + M0[2][0] + M0[4][0] + M0[7][0]) == 2)
After running the cplex with the same parameters input:
The result from MIP 1: obj_lambda=213
The result from MIP 2: obj_lamba= 205.0 and the M0 is [[0,1,1,0,0,0,1,1]]
The correct result of M0 from the MIP 2 should be [[0,1,0,1,0,0,0,2]] as in the M0 (parameter in MIP 1) and obj_lambda=213.
Here is the full code of MIP 1:
import cplex
from docplex.mp.model import Model
import numpy as np
mdl = Model(name='Scheduling MIP 1')
inf = cplex.infinity
bigM= 10000
Total_P = 8 # number of places
Total_T = 6 # number of transitions
r1 = 100 #processing time for PM1 (second)
r2 = 200 #processing time for PM2 (second)
v = 3 #robot moving time (second)
w = 5 #loading or unloading time (second)
M0 = np.array([[0,1,1,0,0,0,1,1]]).reshape(Total_P,1)
M0_temp = M0.reshape(1,Total_P)
M0_final = M0_temp.tolist()*Total_T
h = np.array([v+w,w+r1,v+w,w+r2,v+w,0,0,0]).reshape(Total_P,1)
#Parameter
AT =np.array([[1,-1,0,0,0,0],
[0,1,-1,0,0,0],
[0,0,1,-1,0,0],
[0,0,0,1,-1,0],
[0,0,0,0,1,-1],
[0,-1,1,0,0,0],
[0,0,0,-1,1,0],
[-1,1,-1,1,-1,1]])
AT_temp = AT.transpose()
places = np.array(['p1','p2','p3','p4','p5','p6','p7','p8'])
P_conflict = np.empty((1,Total_P), dtype = object)
P_zero = np.empty((1,Total_P), dtype = object)
#Define the place without conflict place
CP = np.count_nonzero(AT, axis=1, keepdims=True) #calculate the nonzero elements for each row
P_conflict = []
P_zero = []
for a in range(0, len(CP)):
if CP[a].item(0) > 2:
P_conflict.append(places[a])
else:
P_zero.append(places[a])
print('p_zero', np.shape(P_zero))
obj_lambda = [ ]
obj_lambda = mdl.continuous_var(lb = 0, ub=inf, name='obj_lambda')
x = np.array(mdl.continuous_var_list(Total_T, 0, inf, name='x')).reshape(Total_T,)
ind_x = np.where(AT[0] == 1)
print('ind_x', ind_x)
def get_index_value(input):
data = []
for l in range(len(P_zero)):
ind_x = np.where(AT[l] == input)
get_value = x[ind_x]
data.append(get_value)
return data
x_in = get_index_value(1)
x_out = get_index_value(-1)
#constraint 16
for l in range(0,len(P_zero)): #only for P non conflict
#if x_in is not None and x_out is not None and obj_lambda is not None:
mdl.add_constraint(x_out[l][0]-x_in[l][0] >= h[l][0] - M0[l][0]*obj_lambda)
#Decision Var
Z = np.empty((Total_T, Total_T), dtype = object)
for k in range(Total_T):
for i in range(Total_T):
Z[k][i] = mdl.binary_var(name='Z' + str(k+1) + str(',') + str(i+1))
storage_ZAT = []
for k in range(Total_T):
ZAT = np.matmul(Z[k].reshape(1,Total_T),AT_temp)
storage_ZAT.append(ZAT) #storage_ZAT = np.append(storage_ZAT, ZAT, axis=0)
ZAT_final = np.asarray(storage_ZAT).reshape(Total_T,Total_P)
M = np.empty((Total_T, Total_P), dtype = object)
for k in range(0,Total_T):
for l in range (0,Total_P):
if k == Total_T-1:
M[Total_T-1][l] = M0_final[0][l]
else:
M[k][l] = mdl.continuous_var(name='M' + str(k + 1) + str(',') + str(l + 1))
M_prev = np.empty((Total_T, Total_P), dtype = object)
if M is not None:
for k in range(0,Total_T):
for l in range (0,Total_P):
if k is not 0:
M_prev[k][l] = M[k-1][l]
else:
M_prev[0][l] = M0_final[0][l]
#Constraint 17
for k in range(Total_T):
for l in range(Total_P):
mdl.add_constraint(M[k][l] == M_prev[k][l] + ZAT_final[k][l])
#Constraint 18
mdl.add_constraints(mdl.sum(Z[k][i] for k in range(Total_T)) == 1 for i in range(Total_T))
# Constraint 19
mdl.add_constraints(mdl.sum(Z[k][i] for i in range(Total_T)) == 1 for k in range(Total_T))
# # Parameters
VW_temp = [[v + w]]
VW_temp = VW_temp*Total_T
VW = np.array(VW_temp) #eshape(Total_T,)
#Define S
S = np.array(mdl.continuous_var_list(Total_T, 0, inf, name='S')).reshape(Total_T,1)
#Constraint 20
for k in range(Total_T-1):
mdl.add_constraint(S[k][0] - S[k+1][0] <= -VW[k][0])
# # Constraint 21
mdl.add_constraint(S[Total_T-1][0] - (S[0][0] + obj_lambda) <=-VW[Total_T-1][0])
x_temp = x.reshape(Total_T,1)
print('x_temp',x_temp)
# Constraint 22
for k in range(Total_T):
for i in range(Total_T):
mdl.add_constraint(S[k][0] - x_temp[i][0] <= (1-Z[k][i])*bigM) #why x_temp? because it is the reshape of x
# Constraint 23
for k in range(Total_T):
for i in range(Total_T):
mdl.add_constraint(S[k][0] - x_temp[i][0] >= (Z[k][i]-1)*bigM)
mdl.minimize(obj_lambda)
mdl.print_information()
solver = mdl.solve() #(log_output=True)
if solver is not None:
mdl.print_solution()
miu = 1 / obj_lambda.solution_value
print('obj_lamba=', obj_lambda.solution_value)
print('miu =', miu)
else:
print("Solver is error")
Here is the code of MIP 2:
import cplex
from docplex.mp.model import Model
import numpy as np
mdl = Model(name='Scheduling MILP 2')
inf = cplex.infinity
bigM= 10000
Total_P = 8 # number of places
Total_T = 6 # number of transitions
r1 = 100 #processing time for PM1 (second)
r2 = 200 #processing time for PM2 (second)
v = 3 #robot moving time (second)
w = 5 #loading or unloading time (second)
h = np.array([v+w,w+r1,v+w,w+r2,v+w,0,0,0]).reshape(Total_P,1)
#Parameter
AT =np.array([[1,-1,0,0,0,0],
[0,1,-1,0,0,0],
[0,0,1,-1,0,0],
[0,0,0,1,-1,0],
[0,0,0,0,1,-1],
[0,-1,1,0,0,0],
[0,0,0,-1,1,0],
[-1,1,-1,1,-1,1]])
AT_temp = AT.transpose()
places = np.array(['p1','p2','p3','p4','p5','p6','p7','p8'])
P_conflict = np.empty((1,Total_P), dtype = object)
P_zero = np.empty((1,Total_P), dtype = object)
#Define the place without conflict place
CP = np.count_nonzero(AT, axis=1, keepdims=True) #calculate the nonzero elements for each row
P_conflict = []
P_zero = []
for a in range(0, len(CP)):
if CP[a].item(0) > 2:
P_conflict.append(places[a])
else:
P_zero.append(places[a])
print('p_zero', P_zero)
y = [ ] #miu
y = mdl.continuous_var(lb = 0, ub=inf, name='miu') #error: docplex.mp.utils.DOcplexException: Variable obj_lambda cannot be used as denominator of 1
x = np.array(mdl.continuous_var_list(Total_T, 0, inf, name='x')).reshape(Total_T,)
ind_x = np.where(AT[0] == 1)
def get_index_value(input):
data = []
for l in range(len(P_zero)):
ind_x = np.where(AT[l] == input)
get_value = x[ind_x]
data.append(get_value)
return data
x_in_hat = get_index_value(1)
x_out_hat = get_index_value(-1)
M0 = np.empty((Total_P,1), dtype = object)
for l in range(0,Total_P):
M0[l][0]= mdl.integer_var(name='M0'+ str(',') + str(l+1))
M0_temp= M0.reshape(1,Total_P)
M0_final = M0_temp.tolist()*Total_T
for k in range(Total_T):
for l in range(Total_P):
if M0_final[k][l] is not None:
mdl.add_constraint((M0_final[0][1] + M0_final[0][5]) == 1)
mdl.add_constraint((M0_final[0][3] + M0_final[0][6]) == 1)
mdl.add_constraint((M0_final[0][0] + M0_final[0][2] + M0_final[0][4] + M0_final[0][7]) == 2)
#constraint 30
for l in range(0,len(P_zero)): #only for P non conflict
mdl.add_constraint(x_out_hat[l][0]-x_in_hat[l][0] >= h[l][0]*y - M0[l][0])
#Decision Var
Z = np.empty((Total_T, Total_T), dtype = object)
for k in range(Total_T):
for i in range(Total_T):
# if k == 0 and i == 0:
# Z[0][0]=1
# else:
Z[k][i] = mdl.binary_var(name='Z' + str(k+1) + str(',') + str(i+1))
#
storage_ZAT = []
for k in range(Total_T):
ZAT = np.matmul(Z[k].reshape(1,Total_T),AT_temp)
storage_ZAT.append(ZAT) #storage_ZAT = np.append(storage_ZAT, ZAT, axis=0)
ZAT_final = np.asarray(storage_ZAT).reshape(Total_T,Total_P)
M = np.empty((Total_T, Total_P), dtype = object)
for k in range(0,Total_T):
for l in range (0,Total_P):
if k == Total_T-1:
M[Total_T-1][l] = M0_final[0][l]
else:
M[k][l] = mdl.continuous_var(name='M' + str(k + 1) + str(',') + str(l + 1))
M_prev = np.empty((Total_T, Total_P), dtype = object)
if M is not None:
for k in range(0,Total_T):
for l in range (0,Total_P):
if k is not 0:
M_prev[k][l] = M[k-1][l]
else:
M_prev[0][l] = M0_final[0][l]
#
#Constraint 31
for k in range(Total_T):
for l in range(Total_P):
mdl.add_constraint(M[k][l] == M_prev[k][l] + ZAT_final[k][l])
#
#Constraint 32
mdl.add_constraints(mdl.sum(Z[k][i] for k in range(Total_T)) == 1 for i in range(Total_T))
# Constraint 33
mdl.add_constraints(mdl.sum(Z[k][i] for i in range(Total_T)) == 1 for k in range(Total_T))
# #
# # # # Parameters
VW_temp = [[v + w]]
VW_temp = VW_temp*Total_T
VW = np.array(VW_temp) #eshape(Total_T,)
#
S_hat = np.array(mdl.continuous_var_list(Total_T, 0, inf, name='S_hat')).reshape(Total_T,1)
#Constraint 34
for k in range(0,Total_T-1):
if k < (Total_T-1):
mdl.add_constraint(S_hat[k][0] - S_hat[k+1][0] <= -VW[k][0]*y)
if k == (Total_T - 1):
mdl.add_constraint(S_hat[k][0] - (S_hat[0][0] + 1) <= -VW[Total_T - 1][0] * y)# Constraint 35
#
#
x_temp_hat = x.reshape(Total_T,1)
#
# Constraint 36
for k in range(Total_T):
for i in range(Total_T):
mdl.add_constraint((S_hat[k][0] - x_temp_hat[i][0]) <= (1-Z[k][i])*bigM) #why x_temp? because it is the reshape of x
# Constraint 37
for k in range(Total_T):
for i in range(Total_T):
mdl.add_constraint((S_hat[k][0] - x_temp_hat[i][0]) >= (Z[k][i]-1)*bigM)
mdl.maximize(y)
mdl.print_information()
solver = mdl.solve()
if solver is not None:
mdl.print_solution()
obj_lambda = 1/y.solution_value
print('miu =', y.solution_value)
print('obj_lamba=', obj_lambda)
print('M0 final=', M0_temp)
else:
print("Solver is error")
In this case, could anyone tell me: did I define the M0 as a decision variable and put them into the constraints correctly, please?
Thank you.
