增加了HSMO整线优化方法，读取数据增加了供料器部分

base_optimizer/optimizer_common.py

@@ -22,6 +22,7 @@ import pandas as pd
 import matplotlib.pyplot as plt
 import matplotlib
 import traceback
+import openpyxl
 
 matplotlib.use('TkAgg')
 
@@ -67,7 +68,16 @@ t_fix_camera_check = 0.12  # 固定相机检测时间
 T_pp, T_tr, T_nc, T_pl = 2, 5, 25, 0
 
 # 时间参数 （数据拟合获得）
-Fit_cy, Fit_nz, Fit_pu, Fit_pl, Fit_mv = 0.326, 0.870, 0.159, 0.041, 0.001
+# Fit_cy, Fit_nz, Fit_pu, Fit_pl, Fit_mv = 0.326, 0.870, 0.159, 0.041, 0.030
+Fit_cy, Fit_nz, Fit_pu, Fit_pl, Fit_mv = 0.326, 0.870, 0.159, 0.041, 0.035
+
+
+class Point:
+    def __init__(self, _x, _y, _r=0, _h=None):
+        self.x = _x
+        self.y = _y
+        self.r = _r
+        self.h = _h
 
 
 class OptResult:
@@ -108,9 +118,12 @@ class OptInfo:
         print(f'-Pick time: {self.pickup_time: .3f}, Pick distance: {self.pickup_distance: .3f}')
         print(f'-Place time: {self.place_time: .3f}, Place distance: {self.place_distance: .3f}')
         print(
-            f'-Round time: {self.total_time - self.place_time - self.place_time: .3f}, Round distance: '
+            f'-Round time: {self.total_time - self.operation_time - self.pickup_time - self.place_time: .3f}, Round distance: '
             f'{self.total_distance - self.pickup_distance - self.place_distance: .3f}')
 
+        print(f'-Round & place time per cycle: {(self.total_time - self.pickup_time - self.operation_time) * 1000.0 / (self.cycle_counter + 1e-10): .3f}, ', end='')
+        print(f'-Round & place distance per cycle: {(self.total_distance - self.pickup_distance) / (self.cycle_counter + 1e-10): .3f}')
+
         minutes, seconds = int(self.total_time // 60), int(self.total_time) % 60
         millisecond = int((self.total_time - minutes * 60 - seconds) * 60)
 
@@ -856,3 +869,10 @@ def random_division(num, div):
 
 def list_range(start, end=None):
     return list(range(start)) if end is None else list(range(start, end))
+
+
+def kth_indices_partition(num, kth):
+    if len(num) > kth:
+        return np.argpartition(num, kth)
+    else:
+        return np.array(range(len(num)))

base_optimizer/optimizer_interface.py

@@ -15,26 +15,22 @@ def base_optimizer(machine_index, pcb_data, component_data, feeder_data, params,
 
     if params.machine_optimizer == 'cell-division':  # 基于元胞分裂的遗传算法
         component_result, cycle_result, feeder_slot_result = optimizer_celldivision(pcb_data, component_data)
-        placement_result, head_sequence = greedy_placement_route_generation(component_data, pcb_data, component_result,
-                                                                            cycle_result, feeder_slot_result)
+        placement_result, head_sequence = place_allocate_sequence_route_generation(component_data, pcb_data,
+                                                                                   component_result, cycle_result,
+                                                                                   feeder_slot_result)
     elif params.machine_optimizer == 'feeder-priority':  # 基于基座扫描的供料器优先算法
         component_result, cycle_result, feeder_slot_result = feeder_priority_assignment(component_data, pcb_data,
                                                                                         feeder_data)
-        placement_result, head_sequence = greedy_placement_route_generation(component_data, pcb_data, component_result,
-                                                                            cycle_result, feeder_slot_result)
-        # placement_result, head_sequence = beam_search_route_generation(component_data, pcb_data, component_result,
-        #                                                                cycle_result, feeder_slot_result)
-        # placement_result, head_sequence = scan_based_placement_route_generation(component_data, pcb_data,
-        #                                                                         component_result, cycle_result,
-        #                                                                         feeder_slot_result)
-
+        placement_result, head_sequence = scan_based_placement_route_generation(component_data, pcb_data,
+                                                                                   component_result, cycle_result,
+                                                                                   feeder_slot_result)
     elif params.machine_optimizer == 'hybrid-genetic':  # 基于拾取组的混合遗传算法
         component_result, cycle_result, feeder_slot_result, placement_result, head_sequence = optimizer_hybrid_genetic(
             pcb_data, component_data, hinter=hinter)
 
     elif params.machine_optimizer == 'aggregation':  # 基于batch-level的整数规划 + 启发式算法
         component_result, cycle_result, feeder_slot_result, placement_result, head_sequence = optimizer_aggregation(
-            component_data, pcb_data)
+            component_data, pcb_data, hinter=hinter)
     elif params.machine_optimizer == 'genetic-scanning':
         component_result, cycle_result, feeder_slot_result, placement_result, head_sequence = optimizer_genetic_scanning(
             component_data, pcb_data, hinter=hinter)
@@ -43,15 +39,16 @@ def base_optimizer(machine_index, pcb_data, component_data, feeder_data, params,
             component_data, pcb_data, hinter=hinter)
     elif params.machine_optimizer == "two-phase":
         component_result, feeder_slot_result, cycle_result = gurobi_optimizer(pcb_data, component_data, feeder_data,
-                                                                              initial=True, partition=True,
+                                                                              initial=True, partition=False,
                                                                               reduction=True, hinter=hinter)
 
-        placement_result, head_sequence = greedy_placement_route_generation(component_data, pcb_data, component_result,
-                                                                            cycle_result, feeder_slot_result)
+        placement_result, head_sequence = place_allocate_sequence_route_generation(component_data, pcb_data,
+                                                                                   component_result, cycle_result,
+                                                                                   feeder_slot_result)
     else:
         raise 'machine optimizer method ' + params.method + ' is not existed'
 
-    print('----- Placement machine ' + str(machine_index) + ' ----- ')
+    # print('----- Placement machine ' + str(machine_index) + ' ----- ')
     opt_res = OptResult(component_result, cycle_result, feeder_slot_result, placement_result, head_sequence)
     # 估算贴装用时
     info = placement_info_evaluation(component_data, pcb_data, opt_res, hinter=False)
@@ -67,4 +64,6 @@ def base_optimizer(machine_index, pcb_data, component_data, feeder_data, params,
             f'result/{params.filename[:-4]}-{params.line_optimizer}-M0{machine_index} {params.save_suffix}',
             component_data, pcb_data, opt_res)
 
+    # output_optimize_result(f'{params.filename[:-4]}', component_data, pcb_data, opt_res)
+
     return info

base_optimizer/result_analysis.py

@@ -355,15 +355,12 @@ def output_optimize_result(file_path, component_data, pcb_data, optimizer_result
     if 'desc' not in output_data.columns:
         column_index = int(np.where(output_data.columns.values.reshape(-1) == 'part')[0][0])
         output_data.insert(loc=column_index + 1, column='desc', value='')
-    file_dir = file_path[:file_path.rfind('/') + 1]
-    if not os.path.exists(file_dir):
-        os.makedirs(file_dir)
 
-    output_data.to_excel(file_path + '.xlsx', sheet_name='tb1', float_format='%.3f', na_rep='')
+    output_data.to_csv('result/' + file_path + '.txt', sep='\t', float_format='%.3f', header=False, index=False)
 
 
 def optimization_assign_result(component_data, pcb_data, optimizer_result, nozzle_hinter=False, component_hinter=False, 
-                               feeder_hinter=False):
+                               feeder_hinter=False, placement_hinter=False):
     if nozzle_hinter:
         columns = ['H{}'.format(i + 1) for i in range(max_head_index)] + ['cycle']
 
@@ -428,6 +425,29 @@ def optimization_assign_result(component_data, pcb_data, optimizer_result, nozzl
         print(feedr_assign)
         print('')
 
+    if placement_hinter:
+        columns = ['H{}'.format(i + 1) for i in range(max_head_index)] + ['cycle']
+
+        placement_assign = pd.DataFrame(columns=columns)
+        for cycle, _ in enumerate(optimizer_result.placement_assign):
+            placement_assign.loc[cycle, 'cycle'] = 1
+            for head in range(max_head_index):
+                point = optimizer_result.placement_assign[cycle][head]
+                if point != -1:
+                    placement_assign.loc[cycle, 'H{}'.format(head + 1)] = 'P{}'.format(point)
+                else:
+                    placement_assign.loc[cycle, 'H{}'.format(head + 1)] = ''
+
+        headseq_assign = pd.DataFrame(columns=columns)
+        for cycle, headseq in enumerate(optimizer_result.head_sequence):
+            headseq_assign.loc[cycle, 'cycle'] = 1
+            for head in range(len(headseq)):
+                headseq_assign.loc[cycle, 'H{}'.format(head + 1)] = 'H{}'.format(headseq[head])
+
+        print(placement_assign)
+        print(headseq_assign)
+        print('')
+
 
 def placement_info_evaluation(component_data, pcb_data, optimizer_result, hinter=False):
     # === 优化结果参数 ===
@@ -555,28 +575,43 @@ def placement_info_evaluation(component_data, pcb_data, optimizer_result, hinter
 
             # 贴装路径
             if optimizer_result.placement_assign and optimizer_result.head_sequence:
+                head_angle = [0 for _ in range(max_head_index)]
                 for head in optimizer_result.head_sequence[cycle]:
                     index = optimizer_result.placement_assign[cycle][head]
                     if index == -1:
                         continue
                     mount_pos.append([pcb_data.iloc[index]['x'] - head * head_interval + stopper_pos[0],
                                       pcb_data.iloc[index]['y'] + stopper_pos[1]])
-                    mount_angle.append(pcb_data.iloc[index]['r'])
+                    head_angle[head] = pcb_data.iloc[index]['r']
 
                 # 单独计算贴装路径
                 for cntPoints in range(len(mount_pos) - 1):
                     info.place_distance += max(abs(mount_pos[cntPoints][0] - mount_pos[cntPoints + 1][0]),
                                                abs(mount_pos[cntPoints][1] - mount_pos[cntPoints + 1][1]))
 
+                if mount_pos[0][0] < mount_pos[-1][0]:
+                    mount_pos = reversed(mount_pos)
+
                 # 考虑R轴预旋转，补偿同轴角度转动带来的额外贴装用时
-                info.operation_time += head_rotary_time(mount_angle[0])  # 补偿角度转动带来的额外贴装用时
                 info.operation_time += t_nozzle_put * nozzle_put_counter + t_nozzle_pick * nozzle_pick_counter
                 for idx, pos in enumerate(mount_pos):
                     info.operation_time += t_place
-                    move_time = max(axis_moving_time(cur_pos[0] - pos[0], 0), axis_moving_time(cur_pos[1] - pos[1], 1))
+
                     if idx == 0:
+                        move_time = max(axis_moving_time(cur_pos[0] - pos[0], 0),
+                                        axis_moving_time(cur_pos[1] - pos[1], 1))
                         info.round_time += move_time
                     else:
+
+                        cur_head = optimizer_result.head_sequence[cycle][idx]
+                        side_head = cur_head - 1 if cur_head % 2 else cur_head + 1
+                        if optimizer_result.head_sequence[cycle][idx - 1] != side_head:
+                            move_time = max(axis_moving_time(cur_pos[0] - pos[0], 0),
+                                            axis_moving_time(cur_pos[1] - pos[1], 1))
+                        else:
+                            move_time = max(axis_moving_time(cur_pos[0] - pos[0], 0),
+                                            axis_moving_time(cur_pos[1] - pos[1], 1),
+                                            head_rotary_time(head_angle[cur_head] - head_angle[side_head]))
                         info.place_time += move_time
 
                     info.total_distance += max(abs(cur_pos[0] - pos[0]), abs(cur_pos[1] - pos[1]))
@@ -594,5 +629,3 @@ def placement_info_evaluation(component_data, pcb_data, optimizer_result, hinter
     return info
 
 
-
-

base_optimizer/smopt_aggregation.py

@@ -1,12 +1,12 @@
 from base_optimizer.optimizer_common import *
-
+from base_optimizer.smtopt_route import *
 
 def list_range(start, end=None):
     return list(range(start)) if end is None else list(range(start, end))
 
 
 @timer_wrapper
-def optimizer_aggregation(component_data, pcb_data):
+def optimizer_aggregation(component_data, pcb_data, hinter=True):
     # === phase 0: data preparation ===
     M = 1000                                # a sufficient large number
     a, b = 1, 6                             # coefficient
@@ -43,7 +43,7 @@ def optimizer_aggregation(component_data, pcb_data):
 
     # === phase 1: mathematical model solver ===
     mdl = Model('SMT')
-    mdl.setParam('OutputFlag', 0)
+    mdl.setParam('OutputFlag', hinter)
 
     # === Decision Variables ===
     # the number of components of type i that are placed by nozzle type j on placement head k
@@ -104,8 +104,7 @@ def optimizer_aggregation(component_data, pcb_data):
     mdl.setParam("TimeLimit", 100)
 
     mdl.optimize()
-
-    if mdl.Status == GRB.OPTIMAL:
+    if mdl.Status == GRB.OPTIMAL or mdl.Status == GRB.TIME_LIMIT:
         print('total cost = {}'.format(mdl.objval))
 
         # convert cp model solution to standard output
@@ -160,25 +159,9 @@ def optimizer_aggregation(component_data, pcb_data):
         feeder_slot_result = feeder_assignment(component_data, pcb_data, component_result, cycle_result)
 
         # === phase 2: heuristic method ===
-        mount_point_pos = defaultdict(list)
-        for pcb_idx, data in pcb_data.iterrows():
-            part = data['part']
-            part_index = component_data[component_data['part'] == part].index.tolist()[0]
-            mount_point_pos[part_index].append([data['x'], data['y'], pcb_idx])
-
-        for index_ in mount_point_pos.keys():
-            mount_point_pos[index_].sort(key=lambda x: (x[1], x[0]))
-
-        for cycle_idx, _ in enumerate(cycle_result):
-            for _ in range(cycle_result[cycle_idx]):
-                placement_result.append([-1 for _ in range(max_head_index)])
-                for head in range(max_head_index):
-                    if component_result[cycle_idx][head] == -1:
-                        continue
-                    index_ = component_result[cycle_idx][head]
-                    placement_result[-1][head] = mount_point_pos[index_][-1][2]
-                    mount_point_pos[index_].pop()
-                head_sequence.append(dynamic_programming_cycle_path(pcb_data, placement_result[-1], feeder_slot_result[cycle_idx]))
+        placement_result, head_sequence = greedy_level_placing_route_generation(component_data, pcb_data,
+                                                                                component_result, cycle_result,
+                                                                                feeder_slot_result)
 
     else:
         warnings.warn('No solution found!', UserWarning)

base_optimizer/smopt_celldivision.py

@@ -33,7 +33,7 @@ def convert_cell_2_result(pcb_data, component_data, component_cell, population):
 
     wl = [0 for _ in range(max_head_index)]     # workload
 
-    e1, e2, e3 = 1, 0.5, 1. / 6
+    e1, e2, e3 = 1, 2, 1. / 6
 
     component_result, cycle_result, feeder_slot_result = [], [], []
     for index in population:
@@ -101,12 +101,16 @@ def optimizer_celldivision(pcb_data, component_data, hinter=True):
     golden_section = 0.618
 
     # 获取元件元胞
-    point_num = len(pcb_data)
-    component_cell = pd.DataFrame({'index': np.arange(len(component_data)), 'points': np.zeros(len(component_data), dtype=int)})
-    for point_cnt in range(point_num):
-        part = pcb_data.loc[point_cnt, 'part']
-        index = np.where(component_data['part'].values == part)
-        component_cell.loc[index[0], 'points'] += 1
+    feeder_num = sum(component_data['fdn'])
+    component_cell = pd.DataFrame({'index': np.arange(feeder_num), 'points': np.zeros(feeder_num, dtype=int)})
+    cell_index = 0
+    for part_index, data in component_data.iterrows():
+        total_points, division_points = data.points, math.ceil(data.points / data.fdn)
+        for _ in range(data.fdn):
+            component_cell.loc[cell_index, 'index'] = part_index
+            component_cell.loc[cell_index, 'points'] = min(division_points, total_points)
+            total_points -= division_points
+            cell_index += 1
     component_cell = component_cell[~component_cell['points'].isin([0])]
 
     # component_cell.sort_values(by = "points" , inplace = True, ascending = False)

base_optimizer/smopt_hybridgenetic.py

@@ -1,5 +1,5 @@
 from base_optimizer.optimizer_common import *
-
+from base_optimizer.smtopt_route import *
 
 def dynamic_programming_cycle_path(cycle_placement, cycle_points):
     head_sequence = []
@@ -231,8 +231,8 @@ def cal_individual_val(component_nozzle, component_point_pos, designated_nozzle,
 
 def convert_individual_2_result(component_data, component_point_pos, designated_nozzle, pickup_group,
                                 pickup_group_cycle, pair_group, feeder_lane, individual):
+
     component_result, cycle_result, feeder_slot_result = [], [], []
-    placement_result, head_sequence_result = [], []
 
     # === 记录不同元件对应的槽位 ===
     feeder_part_arrange = defaultdict(list)
@@ -273,26 +273,7 @@ def convert_individual_2_result(component_data, component_point_pos, designated_
 
                 pickup_cycle_result[idx][head] -= cycle
 
-    component_point_index = defaultdict(int)
-    for cycle_set in range(len(cycle_result)):
-        for cycle in range(cycle_result[cycle_set]):
-            placement_result.append([-1 for _ in range(max_head_index)])
-            mount_point = [[0, 0] for _ in range(max_head_index)]
-            for head in range(max_head_index):
-                part_index = component_result[cycle_set][head]
-                if part_index == -1:
-                    continue
-
-                part = component_data.iloc[part_index]['part']
-                point_info = component_point_pos[part][component_point_index[part]]
-                
-                placement_result[-1][head] = point_info[2]
-                mount_point[head] = point_info[0:2]
-
-                component_point_index[part] += 1
-            head_sequence_result.append(dynamic_programming_cycle_path(placement_result[-1], mount_point))
-
-    return component_result, cycle_result, feeder_slot_result, placement_result, head_sequence_result
+    return component_result, cycle_result, feeder_slot_result
 
 
 @timer_wrapper
@@ -510,11 +491,6 @@ def optimizer_hybrid_genetic(pcb_data, component_data, hinter=True):
             pop_val.append(val)         # val is related to assembly time
 
         for _ in range(n_generations):
-            # idx = np.argmin(pop_val)
-            # if len(best_pop_val) == 0 or pop_val[idx] < best_pop_val[-1]:
-            #     best_individual = copy.deepcopy(population[idx])
-            # best_pop_val.append(pop_val[idx])
-
             # min-max convert
             max_val = 1.5 * max(pop_val)
             convert_pop_val = list(map(lambda v: max_val - v, pop_val))
@@ -565,6 +541,15 @@ def optimizer_hybrid_genetic(pcb_data, component_data, hinter=True):
             pbar.update(1)
 
     best_individual = population[np.argmin(pop_val)]
+    component_result, cycle_result, feeder_slot_result = convert_individual_2_result(component_data,
+                                                                                     component_point_pos,
+                                                                                     designated_nozzle, pickup_group,
+                                                                                     pickup_group_cycle, pair_group,
+                                                                                     feeder_lane, best_individual)
+    placement_result, head_sequence_result = place_cluster_greedy_route_generation(component_data, pcb_data,
+                                                                                   component_result, cycle_result,
+                                                                                   feeder_slot_result)
+
+
+    return component_result, cycle_result, feeder_slot_result, placement_result, head_sequence_result
 
-    return convert_individual_2_result(component_data, component_point_pos, designated_nozzle, pickup_group,
-                                       pickup_group_cycle, pair_group, feeder_lane, best_individual)

base_optimizer/smopt_mathmodel.py

@@ -1,4 +1,5 @@
 from base_optimizer.optimizer_common import *
+from base_optimizer.smtopt_route import *
 
 
 def head_task_model(component_data, pcb_data, hinter=True):
@@ -6,11 +7,10 @@ def head_task_model(component_data, pcb_data, hinter=True):
     mdl = Model('pick_route')
     mdl.setParam('Seed', 0)
     mdl.setParam('OutputFlag', hinter)          # set whether output the debug information
-    mdl.setParam('TimeLimit', 600)
+    mdl.setParam('TimeLimit', 1000)
 
     H = max_head_index
     I = len(component_data)
-    S = len(component_data)
     K = len(pcb_data)
 
     nozzle_type, component_type = [], []
@@ -29,12 +29,16 @@ def head_task_model(component_data, pcb_data, hinter=True):
     M = 10000
     CompOfNozzle = [[0 for _ in range(J)] for _ in range(I)]    # Compatibility
 
-    component_point = [0 for _ in range(I)]
+    component_point, component_fdn = [0 for _ in range(I)], [0 for _ in range(I)]
+
     for _, data in pcb_data.iterrows():
         idx = component_data[component_data.part == data.part].index.tolist()[0]
         nozzle = component_data.iloc[idx].nz
         CompOfNozzle[idx][nozzle_type.index(nozzle)] = 1
         component_point[idx] += 1
+        component_fdn[idx] = component_data.iloc[idx].fdn
+
+    S = sum(component_fdn)
 
     # objective related
     g = mdl.addVars(list_range(K), vtype=GRB.BINARY)
@@ -94,7 +98,7 @@ def head_task_model(component_data, pcb_data, hinter=True):
         range(-(H - 1) * r, S) for k in range(K))
 
     # feeder related
-    mdl.addConstrs(quicksum(f[s, i] for s in range(S)) <= 1 for i in range(I))
+    mdl.addConstrs(quicksum(f[s, i] for s in range(S)) <= component_fdn[i] for i in range(I))
     mdl.addConstrs(quicksum(f[s, i] for i in range(I)) <= 1 for s in range(S))
     mdl.addConstrs(
         quicksum(x[i, k, h] * y[s, k, h] for h in range(H) for k in range(K)) >= f[s, i] for i in range(I) for s in range(S))
@@ -110,7 +114,7 @@ def head_task_model(component_data, pcb_data, hinter=True):
     # objective
     mdl.setObjective(Fit_cy * quicksum(g[k] for k in range(K)) + Fit_nz * quicksum(
         d[k, h] for h in range(H) for k in range(K)) + Fit_pu * quicksum(
-        e[s, k] for s in range(-(H - 1) * r, S) for k in range(K)) + Fit_mv * head_interval * quicksum(u[k] for k in range(K)),
+        e[s, k] for s in range(-(H - 1) * r, S) for k in range(K)) + Fit_mv * quicksum(u[k] for k in range(K)),
                      GRB.MINIMIZE)
 
     mdl.optimize()
@@ -119,19 +123,21 @@ def head_task_model(component_data, pcb_data, hinter=True):
     for k in range(K):
         if abs(g[k].x) < 1e-6:
             continue
+        component_assign, feeder_slot_assign = [-1 for _ in range(H)], [-1 for _ in range(H)]
 
-        component_result.append([-1 for _ in range(H)])
-        feeder_slot_result.append([-1 for _ in range(H)])
-        cycle_result.append(1)
         for h in range(H):
             for i in range(I):
                 if abs(x[i, k, h].x) > 1e-6:
-                    component_result[-1][h] = i
+                    component_assign[h] = i
 
             for s in range(S):
                 if abs(y[s, k, h].x) > 1e-6:
-                    feeder_slot_result[-1][h] = slot_start + s * interval_ratio - 1
+                    feeder_slot_assign[h] = slot_start + s * interval_ratio - 1
 
+        if sum(component_assign) != -H:
+            component_result.append(component_assign)
+            feeder_slot_result.append(feeder_slot_assign)
+            cycle_result.append(1)
 
     if hinter:
         print(component_result)
@@ -143,7 +149,7 @@ def place_route_model(component_data, pcb_data, component_result, feeder_slot_re
     mdl = Model('place_route')
     mdl.setParam('Seed', 0)
     mdl.setParam('OutputFlag', hinter)  # set whether output the debug information
-    mdl.setParam('TimeLimit', 10)
+    mdl.setParam('TimeLimit', 1000)
 
     component_type = []
     for _, data in component_data.iterrows():
@@ -227,7 +233,7 @@ def place_route_model(component_data, pcb_data, component_result, feeder_slot_re
             else:
                 # there are components on the head
                 mdl.addConstrs(quicksum(w[p, q, k, a] for a in A_from(h) for q in range(P)) + quicksum(
-                    w[q, p, k, a] for a in A_to(h) for q in range(P)) + y[p, k, h] + z[p, k, h] <= 4 *
+                    w[q, p, k, a] for a in A_to(h) for q in range(P)) + y[p, k, h] + z[p, k, h] <= 2 *
                                CompOfPoint[component_result[k][h]][p] for p in range(P))
 
     # each head corresponds to a maximum of one point in each cycle
@@ -362,7 +368,9 @@ def place_route_model(component_data, pcb_data, component_result, feeder_slot_re
 def optimizer_mathmodel(component_data, pcb_data, hinter=True):
 
     component_result, cycle_result, feeder_slot_result = head_task_model(component_data, pcb_data, hinter)
-    placement_result, head_sequence = place_route_model(component_data, pcb_data, component_result, feeder_slot_result)
-    # placement_result, head_sequence = greedy_placement_route_generation(component_data, pcb_data, component_result,
-    #                                                                     cycle_result)
+    # placement_result, head_sequence = place_route_model(component_data, pcb_data, component_result, feeder_slot_result)
+    placement_result, head_sequence = place_allocate_sequence_route_generation(component_data, pcb_data,
+                                                                               component_result, cycle_result,
+                                                                               feeder_slot_result)
+
     return component_result, cycle_result, feeder_slot_result, placement_result, head_sequence

base_optimizer/smopt_twophase.py

@@ -1,4 +1,6 @@
 from base_optimizer.optimizer_common import *
+from base_optimizer.smtopt_route import *
+from base_optimizer.result_analysis import *
 
 
 def list_range(start, end=None):
@@ -80,12 +82,13 @@ def gurobi_optimizer(pcb_data, component_data, feeder_data, reduction=True, part
             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
+#            condition: value change
+            if abs(objbst - 1e+100) > 1:    # 避免未找到可行解提前退出
+                if pre_objbst and abs(pre_objbst - objbst) < 1e-3:
+                    if pre_changetime and changetime - pre_changetime > 90 * (1 - objbnd / objbst):
+                        mdl.terminate()
+                else:
+                    pre_changetime = changetime
 
             pre_objbst = objbst
 
@@ -147,7 +150,7 @@ def gurobi_optimizer(pcb_data, component_data, feeder_data, reduction=True, part
         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
+    S = ratio * sum(component_feeder.values()) 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):
@@ -157,12 +160,12 @@ def gurobi_optimizer(pcb_data, component_data, feeder_data, reduction=True, part
     mdl = Model('SMT')
     mdl.setParam('Seed', 0)
     mdl.setParam('OutputFlag', hinter)  # set whether output the debug information
-    mdl.setParam('TimeLimit', 3600 * 3)
+    mdl.setParam('TimeLimit', 3600)
     mdl.setParam('PoolSearchMode', 2)
-    mdl.setParam('PoolSolutions', 3e2)
+    mdl.setParam('PoolSolutions', 100)
     mdl.setParam('PoolGap', 1e-4)
     # mdl.setParam('MIPFocus', 2)
-    # mdl.setParam("Heuristics", 0.5)
+    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)
@@ -178,7 +181,6 @@ def gurobi_optimizer(pcb_data, component_data, feeder_data, reduction=True, part
         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():
@@ -227,7 +229,7 @@ def gurobi_optimizer(pcb_data, component_data, feeder_data, reduction=True, part
                 continue
 
             part_feederbase[part_2_cpidx[part]] = slot
-            # f[slot, part_2_cpidx[part]].Start = 1
+            f[slot, part_2_cpidx[part]].Start = 1
             slot += 1
 
         for idx, cycle in enumerate(cycle_index):
@@ -334,13 +336,12 @@ def gurobi_optimizer(pcb_data, component_data, feeder_data, reduction=True, part
         in range(L))
 
     if reduction:
-        # mdl.addConstrs(WL[l] >= WL[l + 1] for l in range(L - 1))
+        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))
 
@@ -365,140 +366,142 @@ def gurobi_optimizer(pcb_data, component_data, feeder_data, reduction=True, part
     # mdl.optimize()
 
     # === result generation ===
-    nozzle_assign, component_assign = [], []
-    feeder_assign, cycle_assign = [], []
+    opt_res_list = defaultdict(OptResult)
     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)
+        component_pos = defaultdict(list[Point])
         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])
+            component_pos[component_index].append(Point(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])]
+        for part in component_pos.keys():
+            component_pos[part] = sorted(component_pos[part], key=lambda pos: (pos.x, pos.y))
 
-        min_dist, solution_number = None, -1
+        min_dist, solution_number = None, 0
         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)
+            opt_res = OptResult()
+            # == 转换标准的贴装头分配的解 ===
             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
+                opt_res.cycle_assign.append(round(WL[l].Xn))
+                opt_res.component_assign.append([-1] * max_head_index)
+                opt_res.feeder_slot_assign.append([-1] * 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 h in range(max_head_index):
                     for i in range(I):
                         if abs(y[i, h, l].Xn) <= 1e-4:
                             continue
+                        opt_res.component_assign[-1][h] = i
 
-                        for _ in range(round(WL[l].Xn)):
-                            pos_list.append(component_pos[i][pos_counter[i]])
-                            pos_counter[i] += 1
+                    for s in range(S):
+                        if abs(v[s, h, l].Xn - 1) < 1e-4 and opt_res.component_assign[-1][h] != -1:
+                            opt_res.feeder_slot_assign[-1][h] = s
 
-                        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
+            # 根据贴装头位置，转换供料器槽位
+            cp_avg_head, cp_sum_cycle = defaultdict(float), defaultdict(int)
+            for cycle, component_assign in enumerate(opt_res.component_assign):
+                for head, part in enumerate(component_assign):
+                    if part == -1:
+                        continue
+                    cp_avg_head[part] += opt_res.cycle_assign[cycle] * head
+                    cp_sum_cycle[part] += opt_res.cycle_assign[cycle]
 
-            if min_dist is None or dist < min_dist:
-                min_dist = dist
-                solution_number = sol_counter
+            for part, head in cp_avg_head.items():
+                cp_avg_head[part] = head / cp_sum_cycle[part]
 
-        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)])
+            avg_position = sum([data.x - cp_avg_head[part_2_cpidx[data.part]] * head_interval for _, data in
+                                pcb_data.iterrows()]) / len(pcb_data)
+            avg_slot = 0
+            D_PU, D_PL, D_BW, D_FW = 0, 0, 0, 0
+            for cycle, slots in enumerate(opt_res.feeder_slot_assign):
+                min_slot, max_slot = max_slot_index, 0
+                for head, slot in enumerate(slots):
+                    if slot == -1:
+                        continue
+                    min_slot = min(min_slot, slot - head * ratio)
+                    max_slot = max(max_slot, slot - head * ratio)
+                avg_slot += (max_slot - min_slot) * opt_res.cycle_assign[cycle]
+                D_PU += (max_slot - min_slot) * slot_interval * opt_res.cycle_assign[cycle]     # 拾取路径
 
-            cycle_assign.append(round(WL[l].Xn))
-            if abs(WL[l].Xn) <= 1e-4:
-                continue
+            avg_slot /= sum(opt_res.cycle_assign)
+            start_slot = round((avg_position + stopper_pos[0] - slotf1_pos[0]) / slot_interval + avg_slot / 2) + 1
 
-            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 cycle in range(len(opt_res.feeder_slot_assign)):
+                for head in range(max_head_index):
+                    if (slot := opt_res.feeder_slot_assign[cycle][head]) == -1:
+                        continue
+                    opt_res.feeder_slot_assign[cycle][head] = start_slot + slot * (2 if ratio == 1 else 1)
 
-                        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
+            component_pos_counter = defaultdict(int)
+            cycle_place_pos = defaultdict(list[Point])
             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 cycle in range(len(opt_res.cycle_assign)):
+                    if (part := opt_res.component_assign[cycle][head]) == -1:
+                        continue
 
+                    avg_place_pos = Point(0, 0, _h=head)
+                    for counter in range(round(opt_res.cycle_assign[cycle])):
+                        avg_place_pos.x = (1 - 1.0 / (counter + 1)) * avg_place_pos.x + (
+                                component_pos[part][component_pos_counter[part]].x - head * head_interval) / (
+                                                  counter + 1)
+                        avg_place_pos.y = (1 - 1.0 / (counter + 1)) * avg_place_pos.y + component_pos[part][
+                            component_pos_counter[part]].y / (counter + 1)
+                        component_pos_counter[part] += 1
+                    avg_place_pos.x += stopper_pos[0]
+                    avg_place_pos.y += stopper_pos[1]
+                    cycle_place_pos[cycle].append(avg_place_pos)
+
+            for cycle in range(len(opt_res.cycle_assign)):
+                min_slot, max_slot = max_slot_index, 0
+                for head in range(max_head_index):
+                    if (slot := opt_res.feeder_slot_assign[cycle][head]) == -1:
+                        continue
+                    min_slot = min(min_slot, slot - head * interval_ratio)
+                    max_slot = max(max_slot, slot - head * interval_ratio)
+                # cycle_place_pos[cycle] = sorted(cycle_place_pos[cycle], key=lambda pt: pt.x)
+
+                pick_pos = Point(slotf1_pos[0] + (min_slot + max_slot) / 2 * slot_interval, slotf1_pos[1])
+                _, seq = dynamic_programming_cycle_route(cycle_place_pos[cycle], pick_pos)
+                head_position = [Point(0, 0) for _ in range(max_head_index)]
+                for point in cycle_place_pos[cycle]:
+                    head_position[point.h] = point
+
+                for idx in range(len(seq) - 1):
+                    h1, h2 = seq[idx], seq[idx + 1]
+                    D_PL += max(abs(head_position[h1].x - head_position[h2].x),
+                                abs(head_position[h1].y - head_position[h2].y)) * opt_res.cycle_assign[cycle]
+
+
+            # opt_res.placement_assign, opt_res.head_sequence = scan_based_placement_route_generation(component_data,
+            #                                                                                         pcb_data,
+            #                                                                                         opt_res.component_assign,
+            #                                                                                         opt_res.cycle_assign,
+            #                                                                                         opt_res.feeder_slot_assign,
+            #                                                                                         hinter=False)
+            # info = placement_info_evaluation(component_data, pcb_data, opt_res)
+            # print(f'{info.place_distance + info.pickup_distance: .3f}\t{D_PL + D_PU: .3f}')
+            opt_res_list[sol_counter] = opt_res
+
+    solution_number = 0
+    mdl.Params.SolutionNumber = 0
+    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):
-            if abs(WL[l].Xn) <= 1e-4:
-                continue
+            print(WL[l].Xn, end=', ')
 
-            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)
+        print('')
+        print('result')
+        print('component assignment: ', opt_res_list[solution_number].component_assign)
+        print('feeder assignment: ', opt_res_list[solution_number].feeder_slot_assign)
+        print('cycle assignment: ', opt_res_list[solution_number].cycle_assign)
 
-        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
+    return opt_res_list[solution_number].component_assign, opt_res_list[solution_number].feeder_slot_assign, \
+           opt_res_list[solution_number].cycle_assign