Decision Optimization

import cplex
from docplex.mp.model import Model
import numpy as np

mdl = Model(name='Scheduling')
inf = cplex.infinity
bigM= 1000000

Total_P = 8 
Total_T = 6

r1 = 100 
r2 = 200

v = 1 
w = 1 #l

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) 
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])

y = [ ] #miu
obj_lambda = [ ]
y = mdl.continuous_var(lb = 0, ub=inf, name='miu')
#miu = 1/obj_lambda


x = np.array(mdl.integer_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(Total_P):
    M0[l][0]= mdl.binary_var(name='M0' + str(l+1) + str(',') + str(l+1))

# Constraint 38 (initial marking constraint

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] == 1)

M0_temp = M0.reshape(1,Total_P)
M0_final = M0_temp.tolist()*Total_T

#constraint 30
for l in range(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)    
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.integer_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.integer_var_list(

Total_T, 0, inf, name='S_hat')).reshape(Total_T,1)

#Constraint 34

for k in range(Total_T-1):
     mdl.add_constraint(S_hat[k][0] - S_hat[k+1][0] <= -VW[k][0]*y)

# Constraint 35
mdl.add_constraint(S_hat[

Total_T-1][0] - (S_hat[0][0] + 1) <= -VW[Total_T-1][0]*y)

mdl.maximize(y)
mdl.print_information()
solver = mdl.solve() 
if solver is not None:
    mdl.print_solution()
else:
    print("Solver is error")

obj_lambda = 1/y.solution_value

if mdl.solve():
  print('time', mdl.solve_details.time)