from itertools import product from collections import OrderedDict from ortools.linear_solver import pywraplp from ortools.sat.python import cp_model import numpy as np from utils import get_m_value, parse_data # n_opt, n_mach, operation_data, machine_data = parse_data("gfjsp_10_10_1.txt") n_opt, n_mach, operation_data, machine_data = parse_data("gfjsp_10_5_1.txt") operation_data["0"]["h"] = 0 operation_data["1"]["h"] = 0 operation_data["2"]["h"] = 0 operation_data["3"]["h"] = 20 infinity = 1.0e6 # declae the CP-SAT model model = cp_model.CpModel() solver = cp_model.CpSolver() # define the parameters # minimum lag between the starting time of operation i and the ending time of operation j para_lmin = np.full((n_opt, n_opt), dtype=object, fill_value=-np.inf) # maximum lag between the starting time of operation i and the ending time of operation j para_lmax = np.full((n_opt, n_opt), dtype=object, fill_value=np.inf) # processing time of operation i in machine m # para_p = np.full((n_opt, n_mach), dtype=object, fill_value=np.inf) para_p = np.full((n_mach, n_opt), dtype=object, fill_value=np.inf) # the shape of h in the original file is (n_machine) while the shape of para_h in the # paper is (n_opt, n_mach) # maximum holding time of operation i in machine m para_h = np.empty((n_opt, n_mach), dtype=object) # para_h = np.empty(n_mach, dtype=object) # mapping of operation i to machine m # 20 for the furnaces, 0 for Cutting, Pressing, and Forging # 0: Cutting; 1: Pressing; 2: Forging; 3: Furnace holding_time_dict = { "0": 0, "1": 0, "2": 0, "3": 20, } # weight of operation i in machine m # para_w = np.empty((n_opt, n_mach), dtype=object) para_w = np.full((n_mach, n_opt), dtype=object, fill_value=np.inf) # input/output delay time between two consecutive operations in mahcine m para_delta = np.empty((n_mach), dtype=object) # setup time of machine m when processing operation i before j # para_a = np.full((n_opt, n_opt, n_mach), dtype=object, fill_value=np.inf) para_a = np.full((n_mach, n_opt, n_opt), dtype=object, fill_value=np.inf) # capacity of machine para_mach_capacity = np.empty((n_mach), dtype=object) for m in range(n_mach): # capacity of machine is a set of constant numbers para_mach_capacity[m] = machine_data[str(m)]["c"] # input/output delay time between two consecutive operations in mahcine m # delta(m): loading and unloading time of machine m (=1 for all machines) para_delta[m] = 1 # set up time of machine m when processing operation i before j # a(i,j,m): setup time of machine m when processing operation i before j (aijm = -inf if there # is no setups) for idx_setup, setup_data in enumerate(machine_data[str(m)]["setup_data"][0]): para_a[m, int(setup_data[0]), int(setup_data[1])] = setup_data[2] # maximum holding time of operation i in machine m para_h[:, m] = holding_time_dict[str(machine_data[str(m)]["t"])] # lag time for i in range(n_opt): for idx_lag, lag_data in enumerate(operation_data[str(i)]["lag"]): # minimum lag between the starting time of operation i and the ending time of operation j para_lmin[i, int(lag_data[0])] = lag_data[1] # maximum lag between the starting time of operation i and the ending time of operation j para_lmax[i, int(lag_data[0])] = lag_data[2] for idx_pw, pw_data in enumerate(operation_data[str(i)]["pw"]): # operation_data[str(1)]["pw"] # # the shape of para_p in the original file is the transpose of the shape of para_p # para_p[i, int(pw_data[0])] = pw_data[1] # # the shape of para_w in the original file is the transpose of the shape of para_w # para_w[i, int(pw_data[0])] = pw_data[2] # the shape of para_p in the original file is the transpose of the shape of para_p para_p[int(pw_data[0]), i] = pw_data[1] # the shape of para_w in the original file is the transpose of the shape of para_w para_w[int(pw_data[0]), i] = pw_data[2] # reformat the shape of para_p and para_w para_p = para_p.T para_w = para_w.T # # reshape the shape of para_a # para_a = np.einsum("mij->ijm", para_a) # the big M value big_m = get_m_value(para_p=para_p, para_h=para_h, para_lmin=para_lmin, para_a=para_a) # big_m = 1.0e4 para_p_horizon = np.copy(para_p) para_p_horizon[para_p_horizon == np.inf] = 0 para_h_horizon = np.copy(para_h) para_h_horizon[para_h_horizon == np.inf] = 0 para_lmax_horizon = np.copy(para_lmax) para_lmax_horizon[para_lmax_horizon == np.inf] = 0 horizon = ( np.sum(para_p_horizon, axis=1) + np.sum(para_h_horizon, axis=1) + np.sum(para_lmax_horizon, axis=1) ) horizon = int(np.sum(horizon)) + 1 print(f"horizon: {horizon}") # replace infinities with infinity from the solver para_lmin[para_lmin == -np.inf] = -infinity para_lmin[para_lmin == np.inf] = infinity para_lmax[para_lmax == np.inf] = infinity para_lmax[para_lmax == -np.inf] = -infinity para_p[para_p == np.inf] = infinity para_p[para_p == -np.inf] = -infinity para_w[para_w == np.inf] = infinity para_w[para_w == -np.inf] = -infinity para_a[para_a == np.inf] = infinity para_a[para_a == -np.inf] = -infinity # ===================================================== # convert all the numpy arrays with integers data type para_lmin = para_lmin.astype(int) para_lmax = para_lmax.astype(int) para_p = para_p.astype(int) para_w = para_w.astype(int) para_a = para_a.astype(int) para_delta = para_delta.astype(int) para_mach_capacity = para_mach_capacity.astype(int) # take only part of the data for testing n_opt_selected = 8 n_opt = n_opt_selected para_lmin = para_lmin[:n_opt_selected, :n_opt_selected] para_lmax = para_lmax[:n_opt_selected, :n_opt_selected] para_p = para_p[:n_opt_selected, :] para_h = para_h[:n_opt_selected, :] para_w = para_w[:n_opt_selected, :] para_a = para_a[:n_opt_selected, :n_opt_selected, :] # create variables # scaler = 1.0e4 # horizon *= scaler # if operation i is processed by machine m var_y = np.empty((n_opt, n_mach), dtype=object) for i, m in product(range(n_opt), range(n_mach)): var_y[i, m] = model.NewBoolVar(f"y_{i}_{m}") # if operation i is processed before operation j var_x = np.empty((n_opt, n_opt), dtype=object) for i, j in product(range(n_opt), range(n_opt)): var_x[i, j] = model.NewBoolVar(f"x_{i}_{j}") # if operation i is processed before operation j on machine m var_z = np.empty((n_opt, n_opt, n_mach), dtype=object) for i, j, m in product(range(n_opt), range(n_opt), range(n_mach)): var_z[i, j, m] = model.NewBoolVar(f"z_{i}_{j}_{m}") # starting time of operation i var_s = np.empty((n_opt), dtype=object) for i in range(n_opt): var_s[i] = model.NewIntVar(0, horizon, f"s_{i}") # completion time of operation i var_c = np.empty((n_opt), dtype=object) for i in range(n_opt): var_c[i] = model.NewIntVar(0, horizon, f"c_{i}") # make span var_c_max = model.NewIntVar(0, horizon, "C_max") # define the variables in table of section 3 for i in range(n_opt): # eq. (3) # !!!!!!!!! expr = [para_p[i, m] * var_y[i, m] for i in range(n_mach)] expr = [para_p[i, m] * var_y[i, m] for m in range(n_mach)] model.Add(var_c[i] >= var_s[i] + sum(expr)) # eq. (4) # !!!!!!!!!! expr = [(para_p[i, m] + para_h[i, m]) * var_y[i, m] for i in range(n_mach)] expr = [(para_p[i, m] + para_h[i, m]) * var_y[i, m] for m in range(n_mach)] model.Add(var_c[i] <= var_s[i] + sum(expr)) # eq. (5) # !!!!!!!!!!!!!! # for m in range(n_mach): # # sum of y_im = 1 # model.Add(sum([var_y[i, m] for i in range(n_opt)]) == 1) for i in range(n_opt): # sum of y_im = 1 model.Add(sum([var_y[i, m] for m in range(n_mach)]) == 1) for i, j in product(range(n_opt), range(n_opt)): # eq. (6) # minimum lag between the starting time of operation i and the ending time of operation j model.Add(var_s[j] >= var_c[i] + para_lmin[i, j]) # eq. (7) # maximum lag between the starting time of operation i and the ending time of operation j model.Add(var_s[j] <= var_c[i] + para_lmax[i, j]) # TODO: remove this # setting all the variables to 0 for now para_a[:, :, :] = 0 # eq. (22) and (23) # this new version reimplements the constraints in eq. (22) and (23) in the paper # https://developers.google.com/optimization/cp/channeling for i, j, m in product(np.arange(n_opt), np.arange(n_opt), np.arange(n_mach)): # https://github.com/d-krupke/cpsat-primer # model.Add(x + z == 10).OnlyEnforceIf([b2, b3.Not()]) # eq. (22) # left part of the implication bool_b1 = model.NewBoolVar(f"bool_b1_{i}_{j}_{m}") model.Add(var_y[i, m] + var_y[j, m] == 2).OnlyEnforceIf(bool_b1) model.Add(var_y[i, m] + var_y[j, m] != 2).OnlyEnforceIf(bool_b1.Not()) # model.AddBoolAnd(var_y[i, m], var_y[j, m]).OnlyEnforceIf(bool_b1) # model.AddAtMostOne([var_y[i, m], var_y[j, m]]).OnlyEnforceIf(bool_b1.Not()) # right part of the implication bool_b2 = model.NewBoolVar(f"bool_b2_{i}_{j}_{m}") bool_a1 = model.NewBoolVar(f"bool_a1_{i}_{j}_{m}") bool_a2 = model.NewBoolVar(f"bool_a2_{i}_{j}_{m}") model.Add(var_s[j] >= var_c[i] + para_a[m, i, j]).OnlyEnforceIf(bool_a1) model.Add(var_s[j] < var_c[i] + para_a[m, i, j]).OnlyEnforceIf(bool_a1.Not()) model.Add(var_s[i] >= var_c[j] + para_a[m, j, i]).OnlyEnforceIf(bool_a2) model.Add(var_s[i] < var_c[j] + para_a[m, j, i]).OnlyEnforceIf(bool_a2.Not()) # model.AddBoolOr(bool_a1, bool_a2).OnlyEnforceIf(bool_b2) # model.AddBoolAnd(bool_a1.Not(), bool_a2.Not()).OnlyEnforceIf(bool_b2.Not()) model.Add(bool_a1 + bool_a2 == 1).OnlyEnforceIf(bool_b2) model.Add(bool_a1 + bool_a2 != 1).OnlyEnforceIf(bool_b2.Not()) # the implication model.AddImplication(bool_b1, bool_b2) # # eq. (23) # # left part of the implication is the same as eq. (22), so we can reuse bool_b1 # # right part of the implication # bool_b3 = model.NewBoolVar(f"bool_b3_{i}_{j}_{m}") # max_eq23_first = model.NewIntVar(0, horizon, f"max_eq23_{i}_{j}_{m}_first") # model.AddMaxEquality(max_eq23_first, [0, var_c[i] - var_s[i] - var_c[j] + var_s[j]]) # max_eq23_second = model.NewIntVar(0, horizon, f"max_eq23_{i}_{j}_{m}_second") # model.AddMaxEquality( # max_eq23_second, [0, var_c[j] - var_s[j] - var_c[i] + var_s[i]] # ) # # define the intermediate variable for var_s[j] >= var_s[i] + max_eq23_first + para_delta[m] # bool_a3 = model.NewBoolVar(f"bool_a3_{i}_{j}_{m}") # model.Add(var_s[j] >= var_s[i] + max_eq23_first + para_delta[m]).OnlyEnforceIf( # bool_a3 # ) # model.Add(var_s[j] < var_s[i] + max_eq23_first + para_delta[m]).OnlyEnforceIf( # bool_a3.Not() # ) # # define the intermediate variable for var_s[i] >= var_s[j] + max_eq23_second + para_delta[m] # bool_a4 = model.NewBoolVar(f"bool_a4_{i}_{j}_{m}") # model.Add(var_s[i] >= var_s[j] + max_eq23_second + para_delta[m]).OnlyEnforceIf( # bool_a4 # ) # model.Add(var_s[i] < var_s[j] + max_eq23_second + para_delta[m]).OnlyEnforceIf( # bool_a4.Not() # ) # model.AddBoolOr(bool_a3, bool_a4).OnlyEnforceIf(bool_b3) # model.AddBoolAnd(bool_a3.Not(), bool_a4.Not()).OnlyEnforceIf(bool_b3.Not()) # # the implication # model.AddImplication(bool_b1, bool_b3) # https://developers.google.com/optimization/cp/channeling # TODO: num_t can be a huge number, leading to memory error # sum_time = sum(para_p) + sum(para_h) + sum(para_a) # this is wrong sum_time = horizon # num_t = int(np.sum(sum_time) / 1.0e4) num_t = int(sum_time / 1.0e0) var_u = np.empty((n_opt, n_mach, num_t), dtype=object) for i, m, t in product(range(n_opt), range(n_mach), range(num_t)): var_u[i, m, t] = model.NewIntVar(0, np.max(para_w), f"u_{i}_{m}_{t}") # eq. (25) for idx_m, m in enumerate(np.arange(n_mach)): for idx_t, t in enumerate(np.arange(num_t)): constr_25 = 0 for idx_i, i in enumerate(np.arange(n_opt)): # constr_25 += var_y[idx_i, idx_m] * var_u[idx_i, idx_m, idx_t] # xy = model.NewIntVar(0, 8*5, 'xy') # model.AddMultiplicationEquality(xy, [x, y]) yu = model.NewIntVar(0, np.max(para_w), f"yu_{idx_i}_{idx_m}_{idx_t}") model.AddMultiplicationEquality( yu, [var_y[idx_i, idx_m], var_u[idx_i, idx_m, idx_t]] ) constr_25 += yu model.Add(constr_25 <= para_mach_capacity[idx_m]) # eq. (24) for idx_m, m in enumerate(np.arange(n_mach)): for idx_i, i in enumerate(np.arange(n_opt)): for idx_t, t in enumerate(np.arange(num_t)): bool_list = [] bool_x7 = model.NewBoolVar(f"bool_x7_{i}_{idx_t}") model.Add(var_s[i] <= t).OnlyEnforceIf(bool_x7) bool_list.append(bool_x7) bool_x8 = model.NewBoolVar(f"bool_x8_{i}_{idx_t}") model.Add(var_c[i] >= t).OnlyEnforceIf(bool_x8) bool_list.append(bool_x8) bool_x9_and = model.NewBoolVar(f"bool_x9_{i}_{idx_t}") # model.AddBoolAnd(bool_x9_and == 1).OnlyEnforceIf([bool_x7, bool_x8]) # when bool_x7 and bool_x8 are both true, var_u is is w_im model.Add(var_u[idx_i, idx_m, idx_t] == para_w[idx_i, idx_m]).OnlyEnforceIf( bool_x9_and ) model.Add(var_u[idx_i, idx_m, idx_t] == 0).OnlyEnforceIf(bool_x9_and.Not()) # solve the model status = solver.Solve(model) # cp_model.OPTIMAL: 4 # cp_model.FEASIBLE: 2 # cp_model.INFEASIBLE: 3 if status == cp_model.OPTIMAL or status == cp_model.FEASIBLE: print(f"the optimal makespan is {solver.ObjectiveValue()}") else: print(f"status is {status}") print(f"{solver.SufficientAssumptionsForInfeasibility()}") solver.ObjectiveValue() # [VarX[_].solution_value() for _ in range(X.shape[0])] if status == cp_model.INFEASIBLE: # print infeasible boolean variables index print( "SufficientAssumptionsForInfeasibility = " f"{solver.SufficientAssumptionsForInfeasibility()}" ) # print infeasible boolean variables infeasibles = solver.SufficientAssumptionsForInfeasibility() for i in infeasibles: print("Infeasible constraint: %d" % model.GetBoolVarFromProtoIndex(i))