from base_optimizer.optimizer_common import * def list_range(start, end=None): return list(range(start)) if end is None else list(range(start, end)) @timer_wrapper def gurobi_optimizer(pcb_data, component_data, feeder_data, reduction=True, partition=True, initial=False, hinter=True): # data preparation: convert data to index component_list, nozzle_list = defaultdict(int), defaultdict(int) component_feeder = defaultdict(int) cpidx_2_part, nzidx_2_nozzle, cpidx_2_nzidx = {}, {}, {} arg_slot_rng = None if len(feeder_data) == 0 else [feeder_data.iloc[0].slot, feeder_data.iloc[-1].slot] for idx, data in component_data.iterrows(): part, nozzle = data.part, data.nz cpidx_2_part[idx] = part nz_key = [key for key, val in nzidx_2_nozzle.items() if val == nozzle] nz_idx = len(nzidx_2_nozzle) if len(nz_key) == 0 else nz_key[0] nzidx_2_nozzle[nz_idx] = nozzle component_list[part] = 0 component_feeder[part] = data.fdn cpidx_2_nzidx[idx] = nz_idx for _, data in pcb_data.iterrows(): idx = component_data[component_data.part == data.part].index.tolist()[0] nozzle = component_data.loc[idx].nz nozzle_list[nozzle] += 1 component_list[data.part] += 1 part_feederbase = defaultdict(int) if feeder_data is not None: for _, data in feeder_data.iterrows(): idx = -1 for idx, part_ in cpidx_2_part.items(): if data.part == part_: break assert idx != -1 part_feederbase[idx] = data.slot # part index - slot ratio = 1 if reduction else 2 I, J = len(cpidx_2_part.keys()), len(nzidx_2_nozzle.keys()) # === determine the hyper-parameter of L === # first phase: calculate the number of heads for each type of nozzle nozzle_heads = defaultdict(int) for nozzle in nozzle_list.keys(): nozzle_heads[nozzle] = 1 while sum(nozzle_heads.values()) != max_head_index: max_cycle_nozzle = None for nozzle, head_num in nozzle_heads.items(): if max_cycle_nozzle is None or nozzle_list[nozzle] / head_num > nozzle_list[max_cycle_nozzle] / \ nozzle_heads[max_cycle_nozzle]: max_cycle_nozzle = nozzle assert max_cycle_nozzle is not None nozzle_heads[max_cycle_nozzle] += 1 nozzle_comp_points = defaultdict(list) for part, points in component_list.items(): idx = component_data[component_data.part == part].index.tolist()[0] nozzle = component_data.loc[idx].nz nozzle_comp_points[nozzle].append([part, points]) level = 1 if len(component_list) == 1 or len(component_list) % max_head_index == 0 else 2 part_assignment, cycle_assignment = [], [] def aux_func(info): return max(map(lambda points: max([p[1] for p in points]), info)) pre_objbst, pre_changetime = None, None def terminate_condition(mdl, where): if where == GRB.Callback.MIP: objbst, objbnd = mdl.cbGet(GRB.Callback.MIP_OBJBST), mdl.cbGet(GRB.Callback.MIP_OBJBND) changetime = mdl.cbGet(GRB.Callback.RUNTIME) nonlocal pre_objbst, pre_changetime # condition: value change if pre_objbst and abs(pre_objbst - objbst) < 1e-3: if pre_changetime and changetime - pre_changetime > 100 * (1 - objbnd / objbst): mdl.terminate() else: pre_changetime = changetime pre_objbst = objbst def recursive_assign(assign_points, nozzle_compo_points, cur_level, total_level) -> int: def func(points): return map(lambda points: max([p[1] for p in points]), points) if cur_level > total_level and sum(func(nozzle_compo_points.values())) == 0: return 0 elif assign_points <= 0 and cur_level == 1: return -1 # backtrack elif assign_points <= 0 or cur_level > total_level: return 1 # fail nozzle_compo_points_cpy = copy.deepcopy(nozzle_compo_points) prev_assign = 0 for part in part_assignment[cur_level - 1]: if part != -1: prev_assign += 1 head_idx = 0 for nozzle, head in nozzle_heads.items(): while head: min_idx = -1 for idx, (part, points) in enumerate(nozzle_compo_points_cpy[nozzle]): if points >= assign_points and ( min_idx == -1 or points < nozzle_compo_points_cpy[nozzle][min_idx][1]): min_idx = idx part_assignment[cur_level - 1][head_idx] = -1 if min_idx == -1 else \ nozzle_compo_points_cpy[nozzle][min_idx][0] if min_idx != -1: nozzle_compo_points_cpy[nozzle][min_idx][1] -= assign_points head -= 1 head_idx += 1 cycle_assignment[cur_level - 1] = assign_points for part in part_assignment[cur_level - 1]: if part != -1: prev_assign -= 1 if prev_assign == 0: res = 1 else: points = min(len(pcb_data) // max_head_index + 1, aux_func(nozzle_compo_points_cpy.values())) res = recursive_assign(points, nozzle_compo_points_cpy, cur_level + 1, total_level) if res == 0: return 0 elif res == 1: # All cycles have been completed, but there are still points left to be allocated return recursive_assign(assign_points - 1, nozzle_compo_points, cur_level, total_level) # second phase: (greedy) recursive search to assign points for each cycle set and obtain an initial solution while True: part_assignment = [[-1 for _ in range(max_head_index)] for _ in range(level)] cycle_assignment = [-1 for _ in range(level)] points = min(len(pcb_data) // max_head_index + 1, max(component_list.values())) if recursive_assign(points, nozzle_comp_points, 1, level) == 0: break level += 1 L = len(cycle_assignment) if partition else len(pcb_data) S = ratio * sum(component_feeder.values()) * 2 if len(feeder_data) == 0 else arg_slot_rng[-1] - arg_slot_rng[0] + 1 # the available feeder num M = len(pcb_data) # a sufficiently large number (number of placement points) HC = [[0 for _ in range(J)] for _ in range(I)] for i in range(I): for j in range(J): HC[i][j] = 1 if cpidx_2_nzidx[i] == j else 0 mdl = Model('SMT') mdl.setParam('Seed', 0) mdl.setParam('OutputFlag', hinter) # set whether output the debug information mdl.setParam('TimeLimit', 3600 * 3) mdl.setParam('PoolSearchMode', 2) mdl.setParam('PoolSolutions', 3e2) mdl.setParam('PoolGap', 1e-4) # mdl.setParam('MIPFocus', 2) # mdl.setParam("Heuristics", 0.5) # Use only if other methods, including exploring the tree with the default settings, do not yield a viable solution # mdl.setParam("ZeroObjNodes", 100) # === Decision Variables === x = mdl.addVars(list_range(I), list_range(S), list_range(max_head_index), list_range(L), vtype=GRB.BINARY, name='x') y = mdl.addVars(list_range(I), list_range(max_head_index), list_range(L), vtype=GRB.BINARY, name='y') v = mdl.addVars(list_range(S), list_range(max_head_index), list_range(L), vtype=GRB.BINARY, name='v') c = mdl.addVars(list_range(I), list_range(max_head_index), list_range(L), vtype=GRB.INTEGER, name='c') mdl.addConstrs( c[i, h, l] <= component_list[cpidx_2_part[i]] for i in range(I) for h in range(max_head_index) for l in range(L)) # todo: the condition for upper limits of feeders exceed 1 f = {} for i in range(I): if i not in part_feederbase.keys(): for s in range(S): f[s, i] = mdl.addVar(vtype=GRB.BINARY, name='f_' + str(s) + '_' + str(i)) else: for s in range(S): f[s, i] = 1 if part_feederbase[i] == s + arg_slot_rng[0] else 0 p = mdl.addVars(list_range(-(max_head_index - 1) * ratio, S), list_range(L), vtype=GRB.BINARY, name='p') z = mdl.addVars(list_range(J), list_range(max_head_index), list_range(L), vtype=GRB.BINARY) d = mdl.addVars(list_range(L), list_range(max_head_index), vtype=GRB.INTEGER, name='d') d_plus = mdl.addVars(list_range(J), list_range(max_head_index), list_range(L), vtype=GRB.INTEGER, name='d_plus') d_minus = mdl.addVars(list_range(J), list_range(max_head_index), list_range(L), vtype=GRB.INTEGER, name='d_minus') max_cycle = math.ceil(len(pcb_data) / max_head_index) PU = mdl.addVars(list_range(-(max_head_index - 1) * ratio, S), list_range(L), vtype=GRB.INTEGER, name='PU') WL = mdl.addVars(list_range(L), vtype=GRB.INTEGER, ub=max_cycle, name='WL') NC = mdl.addVars(list_range(max_head_index), vtype=GRB.INTEGER, name='NC') part_2_cpidx = defaultdict(int) for idx, part in cpidx_2_part.items(): part_2_cpidx[part] = idx if initial: # initial some variables to speed up the search process # ensure the priority of the workload assignment cycle_index = sorted(range(len(cycle_assignment)), key=lambda k: cycle_assignment[k], reverse=True) part_list = [] for cycle in cycle_index: cycle_part = part_assignment[cycle] for part in cycle_part: if part != -1 and part not in part_list: part_list.append(part) slot = 0 for part in part_list: if feeder_data is not None: while slot in feeder_data.keys(): slot += 1 # skip assigned feeder slot if part_2_cpidx[part] in part_feederbase.keys(): continue part_feederbase[part_2_cpidx[part]] = slot # f[slot, part_2_cpidx[part]].Start = 1 slot += 1 for idx, cycle in enumerate(cycle_index): WL[idx].Start = cycle_assignment[cycle] for h in range(max_head_index): part = part_assignment[cycle][h] if part == -1: continue i = part_2_cpidx[part] y[i, h, idx].Start = 1 v[part_feederbase[i], h, idx].Start = 1 # === Objective === mdl.setObjective(Fit_cy * quicksum(WL[l] for l in range(L)) + 2 * Fit_nz * quicksum( NC[h] for h in range(max_head_index)) + Fit_pu * quicksum( PU[s, l] for s in range(-(max_head_index - 1) * ratio, S) for l in range(L))) # === Constraint === if not partition: mdl.addConstrs(WL[l] <= 1 for l in range(L)) # work completion mdl.addConstrs(c[i, h, l] == WL[l] * y[i, h, l] for i in range(I) for h in range(max_head_index) for l in range(L)) # mdl.addConstrs( # c[i, h, l] <= max_cycle * y[i, h, l] for i in range(I) for h in range(max_head_index) for l in range(L)) # mdl.addConstrs(c[i, h, l] <= WL[l] for i in range(I) for h in range(max_head_index) for l in range(L)) # mdl.addConstrs( # c[i, h, l] >= WL[l] - max_cycle * (1 - y[i, h, l]) for i in range(I) for h in range(max_head_index) for l in # range(L)) mdl.addConstrs( quicksum(c[i, h, l] for h in range(max_head_index) for l in range(L)) == component_list[cpidx_2_part[i]] for i in range(I)) # variable constraint mdl.addConstrs(quicksum(y[i, h, l] for i in range(I)) <= 1 for h in range(max_head_index) for l in range(L)) # simultaneous pick for s in range(-(max_head_index - 1) * ratio, S): rng = list(range(max(0, -math.floor(s / ratio)), min(max_head_index, math.ceil((S - s) / ratio)))) for l in range(L): mdl.addConstr(quicksum(v[s + h * ratio, h, l] for h in rng) <= max_head_index * p[s, l]) mdl.addConstr(quicksum(v[s + h * ratio, h, l] for h in rng) >= p[s, l]) mdl.addConstrs(PU[s, l] == p[s, l] * WL[l] for s in range(-(max_head_index - 1) * ratio, S) for l in range(L)) # mdl.addConstrs(PU[s, l] <= max_cycle * p[s, l] for s in range(-(max_head_index - 1) * ratio, S) for l in range(L)) # mdl.addConstrs(PU[s, l] <= WL[l] for s in range(-(max_head_index - 1) * ratio, S) for l in range(L)) # mdl.addConstrs( # PU[s, l] >= WL[l] - max_cycle * (1 - p[s, l]) for s in range(-(max_head_index - 1) * ratio, S) for l in # range(L)) # nozzle change mdl.addConstrs( z[j, h, l] - z[j, h, l + 1] == d_plus[j, h, l] - d_minus[j, h, l] for l in range(L - 1) for j in range(J) for h in range(max_head_index)) mdl.addConstrs(z[j, h, 0] - z[j, h, L - 1] == d_plus[j, h, L - 1] - d_minus[j, h, L - 1] for j in range(J) for h in range(max_head_index)) mdl.addConstrs( 2 * d[l, h] == quicksum(d_plus[j, h, l] for j in range(J)) + quicksum(d_minus[j, h, l] for j in range(J)) for l in range(L) for h in range(max_head_index)) mdl.addConstrs(NC[h] == quicksum(d[l, h] for l in range(L)) for h in range(max_head_index)) mdl.addConstrs(quicksum(y[i, h, l] for i in range(I) for h in range(max_head_index)) * M >= WL[l] for l in range(L)) # nozzle-component compatibility mdl.addConstrs( y[i, h, l] <= quicksum(HC[i][j] * z[j, h, l] for j in range(J)) for i in range(I) for h in range(max_head_index) for l in range(L)) # available number of feeder mdl.addConstrs(quicksum(f[s, i] for s in range(S)) <= component_feeder[cpidx_2_part[i]] for i in range(I)) # available number of nozzle mdl.addConstrs(quicksum(z[j, h, l] for h in range(max_head_index)) <= max_head_index for j in range(J) for l in range(L)) # upper limit for occupation for feeder slot mdl.addConstrs(quicksum(f[s, i] for i in range(I)) <= 1 for s in range(S)) mdl.addConstrs( quicksum(v[s, h, l] for s in range(S)) >= quicksum(y[i, h, l] for i in range(I)) for h in range(max_head_index) for l in range(L)) # others mdl.addConstrs(quicksum(z[j, h, l] for j in range(J)) <= 1 for h in range(max_head_index) for l in range(L)) mdl.addConstrs( quicksum(x[i, s, h, l] for h in range(max_head_index) for l in range(L)) >= f[s, i] for i in range(I) for s in range(S)) mdl.addConstrs( quicksum(x[i, s, h, l] for h in range(max_head_index) for l in range(L)) <= M * f[s, i] for i in range(I) for s in range(S)) # mdl.addConstrs( # f[s, i] >= x[i, s, h, l] for s in range(S) for i in range(I) for h in range(max_head_index) for l in range(L)) # # mdl.addConstrs( # quicksum(x[i, s, h, l] for h in range(max_head_index) for l in range(L)) >= f[s, i] for s in # range(S) for i in range(I)) # the constraints to speed up the search process mdl.addConstrs( quicksum(x[i, s, h, l] for i in range(I) for s in range(S)) <= 1 for h in range(max_head_index) for l in range(L)) if reduction: # mdl.addConstrs(WL[l] >= WL[l + 1] for l in range(L - 1)) mdl.addConstr(quicksum(WL[l] for l in range(L)) <= sum(cycle_assignment)) mdl.addConstr(quicksum(WL[l] for l in range(L)) >= math.ceil(len(pcb_data) / max_head_index)) mdl.addConstrs(quicksum(z[j, h, l] for j in range(J) for h in range(max_head_index)) >= quicksum( z[j, h, l + 1] for j in range(J) for h in range(max_head_index)) for l in range(L - 1)) mdl.addConstrs(y[i, h, l] <= WL[l] for i in range(I) for h in range(max_head_index) for l in range(L)) mdl.addConstrs(v[s, h, l] <= WL[l] for s in range(S) for h in range(max_head_index) for l in range(L)) mdl.addConstrs( x[i, s, h, l] >= y[i, h, l] + v[s, h, l] - 1 for i in range(I) for s in range(S) for h in range(max_head_index) for l in range(L)) mdl.addConstrs( x[i, s, h, l] <= y[i, h, l] for i in range(I) for s in range(S) for h in range(max_head_index) for l in range(L)) mdl.addConstrs( x[i, s, h, l] <= v[s, h, l] for i in range(I) for s in range(S) for h in range(max_head_index) for l in range(L)) # === search process === mdl.update() # mdl.write('mdl.lp') if hinter: print('num of constrs: ', str(len(mdl.getConstrs())), ', num of vars: ', str(len(mdl.getVars()))) mdl.optimize(terminate_condition) # mdl.optimize() # === result generation === nozzle_assign, component_assign = [], [] feeder_assign, cycle_assign = [], [] if mdl.Status == GRB.OPTIMAL or mdl.Status == GRB.INTERRUPTED or mdl.Status == GRB.TIME_LIMIT: # === selection from solution pool === component_pos, component_avg_pos = defaultdict(list), defaultdict(list) for _, data in pcb_data.iterrows(): component_index = component_data[component_data.part == data.part].index.tolist()[0] component_pos[component_index].append([data.x, data.y]) for i in component_pos.keys(): component_pos[i] = sorted(component_pos[i], key=lambda pos: (pos[0], pos[1])) component_avg_pos[i] = [sum(map(lambda pos: pos[0], component_pos[i])) / len(component_pos[i]), sum(map(lambda pos: pos[1], component_pos[i])) / len(component_pos[i])] min_dist, solution_number = None, -1 for sol_counter in range(mdl.SolCount): nozzle_assign, component_assign = [], [] feeder_assign, cycle_assign = [], [] mdl.Params.SolutionNumber = sol_counter pos_counter = defaultdict(int) dist = 0 cycle_placement, cycle_points = defaultdict(list), defaultdict(list) for l in range(L): if abs(WL[l].Xn) <= 1e-4: continue cycle_placement[l], cycle_points[l] = [-1] * max_head_index, [None] * max_head_index for h in range(max_head_index): for l in range(L): if abs(WL[l].Xn) <= 1e-4: continue pos_list = [] for i in range(I): if abs(y[i, h, l].Xn) <= 1e-4: continue for _ in range(round(WL[l].Xn)): pos_list.append(component_pos[i][pos_counter[i]]) pos_counter[i] += 1 cycle_placement[l][h] = i cycle_points[l][h] = [sum(map(lambda pos: pos[0], pos_list)) / len(pos_list), sum(map(lambda pos: pos[1], pos_list)) / len(pos_list)] for l in range(L): if abs(WL[l].Xn) <= 1e-4: continue if min_dist is None or dist < min_dist: min_dist = dist solution_number = sol_counter mdl.Params.SolutionNumber = solution_number # === 更新吸嘴、元件、周期数优化结果 === for l in range(L): nozzle_assign.append([-1 for _ in range(max_head_index)]) component_assign.append([-1 for _ in range(max_head_index)]) feeder_assign.append([-1 for _ in range(max_head_index)]) cycle_assign.append(round(WL[l].Xn)) if abs(WL[l].Xn) <= 1e-4: continue for h in range(max_head_index): for i in range(I): if abs(y[i, h, l].Xn - 1) < 1e-4: component_assign[-1][h] = i for j in range(J): if HC[i][j]: nozzle_assign[-1][h] = j for s in range(S): if abs(v[s, h, l].Xn - 1) < 1e-4 and component_assign[l][h] != -1: feeder_assign[l][h] = s # === 更新供料器分配结果 == component_head = defaultdict(int) for i in range(I): cycle_num = 0 for l, component_cycle in enumerate(component_assign): for head, component in enumerate(component_cycle): if component == i: component_head[i] += cycle_assign[l] * head cycle_num += cycle_assign[l] component_head[i] /= cycle_num # 不同元件的加权拾取贴装头 average_pos = 0 for _, data in pcb_data.iterrows(): average_pos += (data.x - component_head[part_2_cpidx[data.part]] * head_interval) average_pos /= len(pcb_data) # 实际贴装位置的加权平均 average_slot = 0 for l in range(L): if abs(WL[l].Xn) <= 1e-4: continue min_slot, max_slot = None, None for head in range(max_head_index): if abs(WL[l].Xn) <= 1e-4 or feeder_assign[l][head] == -1: continue slot = feeder_assign[l][head] - head * ratio if min_slot is None or slot < min_slot: min_slot = slot if max_slot is None or slot > max_slot: max_slot = slot average_slot += (max_slot - min_slot) * cycle_assign[l] average_slot /= sum(cycle_assign) start_slot = round((average_pos + stopper_pos[0] - slotf1_pos[0]) / slot_interval + average_slot / 2) + 1 for l in range(L): if abs(WL[l].Xn) <= 1e-4: continue for h in range(max_head_index): for s in range(S): if abs(v[s, h, l].Xn - 1) < 1e-4 and component_assign[l][h] != -1: feeder_assign[l][h] = start_slot + s * (2 if ratio == 1 else 1) if hinter: print('total cost = {}'.format(mdl.objval)) print('cycle = {}, nozzle change = {}, pick up = {}'.format(quicksum(WL[l].Xn for l in range(L)), quicksum( NC[h].Xn for h in range(max_head_index)), quicksum( PU[s, l].Xn for s in range(-(max_head_index - 1) * ratio, S) for l in range(L)))) print('workload: ') for l in range(L): print(WL[l].Xn, end=', ') print('') print('result') print('nozzle assignment: ', nozzle_assign) print('component assignment: ', component_assign) print('feeder assignment: ', feeder_assign) print('cycle assignment: ', cycle_assign) return component_assign, feeder_assign, cycle_assign