From 97dda65da902ae2bbf54f819a178d2900a97bcd4 Mon Sep 17 00:00:00 2001 From: ken910606 Date: Wed, 10 Dec 2025 21:16:17 +0800 Subject: [PATCH] first upload --- m-tsp-sqa.py | 460 ++++++++++++++++++++++++++++++++++++++ openjijtest.py | 587 +++++++++++++++++++++++++++++++++++++++++++++++++ sqa_tsp.py | 248 +++++++++++++++++++++ 3 files changed, 1295 insertions(+) create mode 100644 m-tsp-sqa.py create mode 100644 openjijtest.py create mode 100644 sqa_tsp.py diff --git a/m-tsp-sqa.py b/m-tsp-sqa.py new file mode 100644 index 0000000..8a186c4 --- /dev/null +++ b/m-tsp-sqa.py @@ -0,0 +1,460 @@ +import numpy as np +import matplotlib.pyplot as plt +from collections import Counter +import math +import openjij as oj +from itertools import permutations + +# ============================ +# 預設距離矩陣：10-city TSP +# ============================ +D = np.array([ + [0, 2, 8, 5, 7, 6, 3, 9, 4, 2], + [2, 0, 5, 7, 3, 8, 4, 6, 9, 1], + [8, 5, 0, 2, 6, 7, 3, 4, 1, 5], + [5, 7, 2, 0, 4, 3, 8, 6, 2, 7], + [7, 3, 6, 4, 0, 5, 9, 2, 7, 3], + [6, 8, 7, 3, 5, 0, 4, 7, 9, 6], + [3, 4, 3, 8, 9, 4, 0, 5, 6, 2], + [9, 6, 4, 6, 2, 7, 5, 0, 8, 3], + [4, 9, 1, 2, 7, 9, 6, 8, 0, 5], + [2, 1, 5, 7, 3, 6, 2, 3, 5, 0] +]) + +N = D.shape[0] +penalty = 20 + +def idx(i, t, N): + return i * N + t + +# ============================ +# KDE 函數 +# ============================ +def kde_1d(samples, num_points=200): + if len(samples) == 0: + return np.array([]), np.array([]) + + xs = np.linspace(min(samples), max(samples), num_points) + n = len(samples) + if n < 2: + return xs, np.zeros_like(xs) + std = np.std(samples) + if std == 0: + std = 1.0 + h = 1.06 * std * (n ** (-1/5)) + if h == 0: + h = 1.0 + ys = [] + inv_sqrt_2pi = 1.0 / math.sqrt(2.0 * math.pi) + for x in xs: + s = 0.0 + for xi in samples: + z = (x - xi) / h + s += math.exp(-0.5 * z * z) * inv_sqrt_2pi + ys.append(s / (n * h)) + return xs, np.array(ys) + +# ============================ +# Route diversity matrix 函數 +# ============================ +def route_diversity_matrix(routes, N): + freq = np.zeros((N, N), dtype=float) + valid_routes = 0 + for route in routes: + if len(route) == N and all(0 <= city < N for city in route): + for t, city in enumerate(route): + if 0 <= city < N: # 確保city有效 + freq[city, t] += 1 + valid_routes += 1 + if valid_routes > 0: + freq /= valid_routes + return freq + +# ============================ +# 建立 QUBO：位置編碼 + 固定起點 0 +# ============================ +def build_qubo(D, penalty, fix_start=True): + N = D.shape[0] + Q = {} + + # 距離項 + for t in range(N): + t2 = (t + 1) % N + for i in range(N): + for j in range(N): + if i != j: + u = idx(i, t, N) + v = idx(j, t2, N) + Q[(u, v)] = Q.get((u, v), 0) + D[i, j] + + # 約束A: 每個 time slot 有且只有一個城市 + for t in range(N): + for i in range(N): + u = idx(i, t, N) + Q[(u, u)] = Q.get((u, u), 0) - penalty + for j in range(i + 1, N): + v = idx(j, t, N) + Q[(u, v)] = Q.get((u, v), 0) + 2 * penalty + + # 約束B: 每個城市被訪問一次 + for i in range(N): + for t in range(N): + u = idx(i, t, N) + Q[(u, u)] = Q.get((u, u), 0) - penalty + for t2 in range(t + 1, N): + v = idx(i, t2, N) + Q[(u, v)] = Q.get((u, v), 0) + 2 * penalty + + # 固定起點：t=0 必須是 city 0 + if fix_start: + big = 9999.0 + # 禁止 i!=0 在 t=0 + for i in range(1, N): + u = idx(i, 0, N) + Q[(u, u)] = Q.get((u, u), 0) + big + # 獎勵 city 0 在 t=0 + u0 = idx(0, 0, N) + Q[(u0, u0)] = Q.get((u0, u0), 0) - big + + return Q, N + +# ============================ +# Route decode & cost +# ============================ +def decode_route(sample, N): + route = [] + for t in range(N): + city_for_t = None + for i in range(N): + if sample.get(idx(i, t, N), 0) == 1: + city_for_t = i + break + if city_for_t is None: + city_for_t = -1 + route.append(city_for_t) + return route + +def compute_cost(route, D): + if -1 in route: + return 99999 # 無效路徑懲罰 + N = len(route) + total = 0 + for k in range(N): + a = route[k] + b = route[(k + 1) % N] + total += D[a, b] + return total + +# ============================ +# 兩台UAV TSP分配策略 +# ============================ +def split_route_for_two_uavs(route): + """ + 將單一TSP路徑分配給兩台UAV + 策略：交替分配城市，確保負載平衡 + """ + if len(route) == 0: + return [], [] + + # 移除無效城市 + valid_route = [city for city in route if city >= 0 and city < N] + + if len(valid_route) == 0: + return [], [] + + # 確保起點是0 + if valid_route[0] != 0: + if 0 in valid_route: + valid_route.remove(0) + valid_route = [0] + valid_route + else: + valid_route = [0] + valid_route + + # 策略1: 奇偶分配 + uav1_path = [] + uav2_path = [] + + for i, city in enumerate(valid_route): + if i % 2 == 0: + uav1_path.append(city) + else: + uav2_path.append(city) + + return uav1_path, uav2_path + +def compute_two_uav_cost(uav1_path, uav2_path, D): + """計算兩台UAV的總成本""" + def path_cost(path): + if len(path) <= 1: + return 0 + cost = 0 + for i in range(len(path)): + current = path[i] + next_city = path[(i + 1) % len(path)] + cost += D[current, next_city] + return cost + + cost1 = path_cost(uav1_path) + cost2 = path_cost(uav2_path) + return cost1 + cost2, cost1, cost2 + +# ============================ +# 經典 TSP 解法：暴力枚舉 (限制小規模) +# ============================ +def classic_tsp_two_uav(D, max_permutations=1000): + """ + 經典TSP求解並分配給兩台UAV + 由於10城市的排列數太大(9!=362880)，限制搜索數量 + """ + cities = list(range(1, len(D))) # 排除起點0 + best_total_cost = float('inf') + best_route = None + best_uav1 = None + best_uav2 = None + + count = 0 + for perm in permutations(cities): + if count >= max_permutations: + break + count += 1 + + route = [0] + list(perm) # 固定0為起點 + + # 分配給兩台UAV + uav1_path, uav2_path = split_route_for_two_uavs(route) + total_cost, cost1, cost2 = compute_two_uav_cost(uav1_path, uav2_path, D) + + if total_cost < best_total_cost: + best_total_cost = total_cost + best_route = route + best_uav1 = uav1_path + best_uav2 = uav2_path + + return best_uav1, best_uav2, best_total_cost, best_route + +# ============================ +# 跑 SA or SQA 多次 +# ============================ +def run_algorithm(name, sampler, Q, D, N, num_runs=20, num_reads=20): + all_costs = [] + all_routes = [] + all_uav1_paths = [] + all_uav2_paths = [] + all_uav_costs = [] + + print(f"\n=== Running {name} ===") + for r in range(num_runs): + result = sampler.sample_qubo(Q, num_reads=num_reads) + best = result.first.sample + route = decode_route(best, N) + cost = compute_cost(route, D) + + # 分配給兩台UAV + uav1_path, uav2_path = split_route_for_two_uavs(route) + total_uav_cost, cost1, cost2 = compute_two_uav_cost(uav1_path, uav2_path, D) + + all_costs.append(cost) + all_routes.append(tuple(route)) + all_uav1_paths.append(tuple(uav1_path)) + all_uav2_paths.append(tuple(uav2_path)) + all_uav_costs.append(total_uav_cost) + + print(f"{name} Run {r+1:02d}: Route={route}") + print(f" TSP Cost={cost}, UAV1={uav1_path}(cost={cost1}), UAV2={uav2_path}(cost={cost2}), Total={total_uav_cost}") + + return all_costs, all_routes, all_uav1_paths, all_uav2_paths, all_uav_costs + +# ============================ +# 創建綜合結果圖表 +# ============================ +def create_comprehensive_charts(sa_results, sqa_results, classic_results): + """創建綜合分析圖表""" + fig = plt.figure(figsize=(18, 12)) + + # 解包結果 + sa_costs, sa_routes, sa_uav1, sa_uav2, sa_uav_costs = sa_results + sqa_costs, sqa_routes, sqa_uav1, sqa_uav2, sqa_uav_costs = sqa_results + classic_uav1, classic_uav2, classic_cost, classic_route = classic_results + + # 1. TSP成本分布比較 (左上) + ax1 = plt.subplot(2, 3, 1) + bins = range(min(sa_costs + sqa_costs), max(sa_costs + sqa_costs) + 2) + ax1.hist(sa_costs, bins=bins, alpha=0.6, label="SA", density=True, color='lightblue') + ax1.hist(sqa_costs, bins=bins, alpha=0.6, label="SQA", density=True, color='lightcoral') + + # 添加KDE + if len(sa_costs) > 1: + xs_sa, ys_sa = kde_1d(sa_costs) + ax1.plot(xs_sa, ys_sa, label="SA KDE", color='blue', linewidth=2) + if len(sqa_costs) > 1: + xs_sqa, ys_sqa = kde_1d(sqa_costs) + ax1.plot(xs_sqa, ys_sqa, label="SQA KDE", color='red', linewidth=2) + + ax1.axvline(x=classic_cost, color='green', linestyle='--', + label=f'Classic: {classic_cost}', linewidth=2) + ax1.set_title("TSP Cost Distribution") + ax1.set_xlabel("Tour Cost") + ax1.set_ylabel("Density") + ax1.legend() + ax1.grid(True, alpha=0.3) + + # 2. UAV總成本分布比較 (右上) + ax2 = plt.subplot(2, 3, 2) + uav_bins = range(min(sa_uav_costs + sqa_uav_costs), max(sa_uav_costs + sqa_uav_costs) + 2) + ax2.hist(sa_uav_costs, bins=uav_bins, alpha=0.6, label="SA UAV", density=True, color='lightblue') + ax2.hist(sqa_uav_costs, bins=uav_bins, alpha=0.6, label="SQA UAV", density=True, color='lightcoral') + + classic_uav_total, _, _ = compute_two_uav_cost(classic_uav1, classic_uav2, D) + ax2.axvline(x=classic_uav_total, color='green', linestyle='--', + label=f'Classic UAV: {classic_uav_total}', linewidth=2) + ax2.set_title("Two-UAV Total Cost Distribution") + ax2.set_xlabel("Total UAV Cost") + ax2.set_ylabel("Density") + ax2.legend() + ax2.grid(True, alpha=0.3) + + # 3. 箱型圖比較 (中上) + ax3 = plt.subplot(2, 3, 3) + box_data = [sa_costs, sqa_costs, sa_uav_costs, sqa_uav_costs] + bp = ax3.boxplot(box_data, labels=["SA\nTSP", "SQA\nTSP", "SA\nUAV", "SQA\nUAV"], + patch_artist=True, showmeans=True) + + colors = ['lightblue', 'lightcoral', 'skyblue', 'salmon'] + for patch, color in zip(bp['boxes'], colors): + patch.set_facecolor(color) + + ax3.set_title("Cost Comparison Box Plot") + ax3.set_ylabel("Cost") + ax3.grid(True, alpha=0.3) + + # 4. Route diversity heatmap - SA (左下) + ax4 = plt.subplot(2, 3, 4) + sa_mat = route_diversity_matrix(sa_routes, N) + im1 = ax4.imshow(sa_mat, aspect='auto', origin='lower', cmap='viridis') + ax4.set_title("SA Route Diversity") + ax4.set_xlabel("Position") + ax4.set_ylabel("City") + plt.colorbar(im1, ax=ax4, shrink=0.8) + + # 5. Route diversity heatmap - SQA (右下) + ax5 = plt.subplot(2, 3, 5) + sqa_mat = route_diversity_matrix(sqa_routes, N) + im2 = ax5.imshow(sqa_mat, aspect='auto', origin='lower', cmap='viridis') + ax5.set_title("SQA Route Diversity") + ax5.set_xlabel("Position") + ax5.set_ylabel("City") + plt.colorbar(im2, ax=ax5, shrink=0.8) + + # 6. 最佳路徑可視化 (中下) + ax6 = plt.subplot(2, 3, 6) + + # 找出最佳量子解 + best_sa_idx = sa_costs.index(min(sa_costs)) + best_sqa_idx = sqa_costs.index(min(sqa_costs)) + best_sa_route = list(sa_routes[best_sa_idx]) + best_sqa_route = list(sqa_routes[best_sqa_idx]) + + # 簡單的圓形佈局 + angles = np.linspace(0, 2*np.pi, N, endpoint=False) + x_pos = np.cos(angles) + y_pos = np.sin(angles) + + # 繪製城市 + ax6.scatter(x_pos, y_pos, s=200, c='red', zorder=5) + for i, (x, y) in enumerate(zip(x_pos, y_pos)): + ax6.annotate(str(i), (x, y), ha='center', va='center', + fontsize=10, fontweight='bold', color='white', zorder=6) + + # 選擇更好的路徑來顯示 + if min(sa_costs) <= min(sqa_costs): + display_route = best_sa_route + route_name = "SA" + route_cost = min(sa_costs) + else: + display_route = best_sqa_route + route_name = "SQA" + route_cost = min(sqa_costs) + + # 繪製路徑 + if all(0 <= city < N for city in display_route): + route_x = [x_pos[city] for city in display_route] + [x_pos[display_route[0]]] + route_y = [y_pos[city] for city in display_route] + [y_pos[display_route[0]]] + ax6.plot(route_x, route_y, 'b-', linewidth=2, alpha=0.8, zorder=3) + + ax6.set_title(f"Best {route_name} Route\nCost: {route_cost}") + ax6.set_xlim(-1.3, 1.3) + ax6.set_ylim(-1.3, 1.3) + ax6.set_aspect('equal') + ax6.axis('off') + + plt.tight_layout() + plt.savefig("tsp_two_uav_comprehensive_analysis.png", dpi=300, bbox_inches='tight') + print(f"\n📊 綜合分析圖表已保存: tsp_two_uav_comprehensive_analysis.png") + plt.show() + +# ============================ +# Main +# ============================ +if __name__ == "__main__": + print("🚀 開始10城市兩台UAV TSP分析...") + print(f"Distance Matrix ({N}x{N}):") + print(D) + + Q, N = build_qubo(D, penalty, fix_start=True) + + # Samplers + sqa_sampler = oj.SQASampler() + sa_sampler = oj.SASampler() + + # 參數設定 + NUM_RUNS = 15 + NUM_READS = 20 + + # 運行量子退火算法 + print("🔥 Running Quantum Annealing Algorithms...") + sa_results = run_algorithm("SA", sa_sampler, Q, D, N, NUM_RUNS, NUM_READS) + sqa_results = run_algorithm("SQA", sqa_sampler, Q, D, N, NUM_RUNS, NUM_READS) + + # 運行經典TSP解法 + print("\n🎯 Running Classic TSP (limited search)...") + classic_results = classic_tsp_two_uav(D, max_permutations=5000) + classic_uav1, classic_uav2, classic_total_cost, classic_route = classic_results + + print(f"Classic TSP Results:") + print(f" Best Route: {classic_route}") + print(f" UAV1 Path: {classic_uav1}") + print(f" UAV2 Path: {classic_uav2}") + print(f" Total Cost: {classic_total_cost}") + + # 統計分析 + sa_costs, sa_routes, sa_uav1, sa_uav2, sa_uav_costs = sa_results + sqa_costs, sqa_routes, sqa_uav1, sqa_uav2, sqa_uav_costs = sqa_results + + print(f"\n📊 Algorithm Performance Summary:") + print(f"{'Algorithm':<15} {'TSP Cost':<20} {'UAV Total Cost':<20}") + print(f"{'-'*55}") + print(f"{'SA':<15} min={min(sa_costs):<3} mean={np.mean(sa_costs):<6.1f} " + f"min={min(sa_uav_costs):<3} mean={np.mean(sa_uav_costs):<6.1f}") + print(f"{'SQA':<15} min={min(sqa_costs):<3} mean={np.mean(sqa_costs):<6.1f} " + f"min={min(sqa_uav_costs):<3} mean={np.mean(sqa_uav_costs):<6.1f}") + print(f"{'Classic':<15} {classic_total_cost:<23} {classic_total_cost:<20}") + + # 創建綜合圖表 + print("\n🎨 Creating comprehensive analysis charts...") + create_comprehensive_charts(sa_results, sqa_results, classic_results) + + # 最終比較 + best_sa_uav = min(sa_uav_costs) + best_sqa_uav = min(sqa_uav_costs) + + print(f"\n🏆 Final Two-UAV TSP Results:") + print(f"SA best UAV cost: {best_sa_uav}") + print(f"SQA best UAV cost: {best_sqa_uav}") + print(f"Classic UAV cost: {classic_total_cost}") + + if best_sa_uav <= best_sqa_uav and best_sa_uav <= classic_total_cost: + print("🥇 Winner: SA") + elif best_sqa_uav <= best_sa_uav and best_sqa_uav <= classic_total_cost: + print("🥇 Winner: SQA") + else: + print("🥇 Winner: Classic Algorithm") \ No newline at end of file diff --git a/openjijtest.py b/openjijtest.py new file mode 100644 index 0000000..cbd3aba --- /dev/null +++ b/openjijtest.py @@ -0,0 +1,587 @@ +import numpy as np +import matplotlib.pyplot as plt +from itertools import permutations +import math +import openjij as oj +import random + +# ============================ +# 1. KDE 函數 +# ============================ +def kde_1d(samples, num_points=200): + if len(samples) == 0: + return np.array([]), np.array([]) + + xs = np.linspace(min(samples), max(samples), num_points) + n = len(samples) + if n < 2: + return xs, np.zeros_like(xs) + + std = np.std(samples) + if std == 0: + std = 1.0 + + h = 1.06 * std * (n ** (-1/5)) + ys = [] + inv_sqrt_2pi = 1.0 / math.sqrt(2 * math.pi) + + for x in xs: + s = 0.0 + for xi in samples: + z = (x - xi) / h + s += math.exp(-0.5 * z * z) * inv_sqrt_2pi + ys.append(s / (n * h)) + + return xs, np.array(ys) + + +# ============================ +# 2. route diversity matrix（支援可變長度路徑） +# ============================ +def route_diversity_matrix(routes, N): + """ + routes: list of routes，每條 route 是城市序列，例如 [0, 2, 5, 1, 9] + 我們只用前 N 個位置來統計（超過 N 的部分忽略） + """ + freq = np.zeros((N, N), dtype=float) + valid_routes = 0 + + for route in routes: + if len(route) == 0: + continue + valid_routes += 1 + for t in range(min(len(route), N)): + city = route[t] + if 0 <= city < N: + freq[city, t] += 1 + + if valid_routes > 0: + freq /= valid_routes + + return freq + + +# ============================ +# 3. MTSP 變數編碼：x[k,i,p] → index +# ============================ +def idx_mtsp(k, i, p, N): + """ x[k,i,p] → 1D index """ + return k * (N * N) + i * N + p + + +# ============================ +# 4. 建立 2-UAV M-TSP QUBO（城市數不固定） +# ============================ +def build_mtsp_qubo_variable(D, penalty=20.0, big_penalty=9999.0): + """ + - 城市: 0..N-1，其中 0 為 depot + - 兩台 UAV: k = 0,1 + - 每台 UAV 有 N 個 slot: p = 0..N-1 + - 變數: x_{k,i,p} + QUBO包含： + 距離成本 + 位置約束 + 城市約束 + 起點約束 + """ + D = np.asarray(D) + N = D.shape[0] + Q = {} + + def addQ(u, v, w): + if w == 0: + return + if u > v: + u, v = v, u + Q[(u, v)] = Q.get((u, v), 0.0) + w + + # ---- 距離成本 ---- + for k in range(2): + for p in range(N): + q = (p + 1) % N + for i in range(N): + for j in range(N): + dij = D[i, j] + if dij == 0: + continue + u = idx_mtsp(k, i, p, N) + v = idx_mtsp(k, j, q, N) + addQ(u, v, dij) + + # ---- 位置約束：每 slot 只能一個城市 ---- + for k in range(2): + for p in range(N): + vars_pos = [idx_mtsp(k, i, p, N) for i in range(N)] + + for u in vars_pos: + addQ(u, u, -penalty) + + for a in range(N): + for b in range(a+1, N): + addQ(vars_pos[a], vars_pos[b], 2*penalty) + + # ---- 城市約束：1..N-1 每城恰好一次 ---- + for i in range(1, N): + vars_city = [] + for k in range(2): + for p in range(N): + vars_city.append(idx_mtsp(k, i, p, N)) + + for u in vars_city: + addQ(u, u, -penalty) + + L = len(vars_city) + for a in range(L): + for b in range(a+1, L): + addQ(vars_city[a], vars_city[b], 2*penalty) + + # ---- 起點約束：每台 UAV 的 p=0 必須是 city 0 ---- + for k in range(2): + for i in range(1, N): + addQ(idx_mtsp(k, i, 0, N), idx_mtsp(k, i, 0, N), big_penalty) + addQ(idx_mtsp(k, 0, 0, N), idx_mtsp(k, 0, 0, N), -big_penalty) + + return Q, N + + +# ============================ +# 5. 解碼 slot 序列（不修剪） +# ============================ +def decode_slots(sample, N): + """ + 回傳： + u1_slots, u2_slots + 各自長度 N，裡面是城市 index（可能重複，也可能很多 0） + """ + slots = [] + for k in range(2): + row = [] + for p in range(N): + chosen = 0 # 預設當作待在 depot + for i in range(N): + if sample.get(idx_mtsp(k, i, p, N), 0) == 1: + chosen = i + break + row.append(chosen) + slots.append(row) + return slots[0], slots[1] + + +# ============================ +# 6. 根據 slots 做 repair → 產生合法 uav1 / uav2 path +# ============================ +def repair_routes_from_slots(u1_slots, u2_slots, N): + """ + 依據 slots 傾向資訊，把城市分配給兩台 UAV， + 並用暴力 TSP 在各自子集合裡找最短路徑，確保： + - 每個城市 1..N-1 被恰好一台 UAV 拜訪 + - 兩台 UAV 路徑皆為 0 -> ... -> 0（成本計算時處理） + """ + + # 統計每個城市在兩台 UAV slots 中出現次數 + count1 = [0]*N + count2 = [0]*N + for c in u1_slots: + if 0 <= c < N: + count1[c] += 1 + for c in u2_slots: + if 0 <= c < N: + count2[c] += 1 + + # 分配城市：看誰出現次數多，平手則讓城市較少的一邊拿 + assign1 = [] + assign2 = [] + for city in range(1, N): # city 0 = depot 不分配 + if count1[city] > count2[city]: + assign1.append(city) + elif count2[city] > count1[city]: + assign2.append(city) + else: + # 平手 → 給目前負載較少的 UAV + if len(assign1) <= len(assign2): + assign1.append(city) + else: + assign2.append(city) + + return sorted(assign1), sorted(assign2) + + +# ============================ +# 7. 單 UAV 成本 & 子集合最佳排列 +# ============================ +def uav_cost(path, D): + if not path: + return 0.0 + path = [c for c in path if 0 <= c < len(D)] + if len(path) == 0: + return 0.0 + + cost = D[0, path[0]] + for i in range(len(path)-1): + cost += D[path[i], path[i+1]] + cost += D[path[-1], 0] + return float(cost) + + +def best_order_for_cities(cities, D): + """ + 對給定 city set，暴力搜尋最短 0 -> perm -> 0 的路徑 + cities: list of city index (1..N-1) + """ + if len(cities) <= 1: + return list(cities), uav_cost(list(cities), D) + + best_perm = None + best_cost = float('inf') + for perm in permutations(cities): + c = uav_cost(list(perm), D) + if c < best_cost: + best_cost = c + best_perm = list(perm) + return best_perm, best_cost + + +# ============================ +# 8. 利用 QUBO 取樣 + repair 得到合法解 +# ============================ +def sample_and_repair(sampler, Q, N, D): + """ + 一次取樣 + repair，即可保證得到合法 2-UAV M-TSP 解。 + 不需要反覆重試。 + """ + result = sampler.sample_qubo(Q, num_reads=30) + sample = result.first.sample + + # 1) 先解出 slots + u1_slots, u2_slots = decode_slots(sample, N) + + # 2) 依 slots 傾向分配城市給兩台 UAV + assign1, assign2 = repair_routes_from_slots(u1_slots, u2_slots, N) + + # 3) 各自子集合內找最短路徑 + u1_path, c1 = best_order_for_cities(assign1, D) + u2_path, c2 = best_order_for_cities(assign2, D) + + total = c1 + c2 + + return u1_path, u2_path, c1, c2, total + + + +# ============================ +# 9. Classic 2-UAV MTSP 暴力搜尋（全域最佳） +# ============================ +def classic_mtsp_two_uav_bruteforce(D): + N = D.shape[0] + cities = list(range(1, N)) + best_total = float('inf') + best_u1, best_u2 = None, None + best_route = None + + def cost_of(path): + return uav_cost(path, D) + + for perm in permutations(cities): + for cut in range(len(perm)+1): + u1 = list(perm[:cut]) + u2 = list(perm[cut:]) + total = cost_of(u1) + cost_of(u2) + if total < best_total: + best_total = total + best_u1 = u1 + best_u2 = u2 + best_route = [0] + list(perm) + + return best_u1, best_u2, best_total, best_route + + +# ============================ +# 10. 共用：計算兩 UAV 總成本 +# ============================ +def compute_two_uav_cost(uav1_path, uav2_path, D): + c1 = uav_cost(uav1_path, D) + c2 = uav_cost(uav2_path, D) + return c1 + c2, c1, c2 + + +# ============================ +# 11. SA & SQA：每 run 使用 sample_and_repair +# ============================ +def run_mtsp_sa_sqa(Q, N, D, NUM_RUNS=10): + sa_costs = [] + sa_routes = [] + sa_u1_list = [] + sa_u2_list = [] + sa_uav_costs = [] # store tuple (c1, c2) + + sqa_costs = [] + sqa_routes = [] + sqa_u1_list = [] + sqa_u2_list = [] + sqa_uav_costs = [] # store tuple (c1, c2) + + sa = oj.SASampler() + sqa = oj.SQASampler() + + print("\n🔥 Running SA (with repair, always feasible)...") + for r in range(NUM_RUNS): + u1, u2, c1, c2, total = sample_and_repair(sa, Q, N, D) + + full_route = [0] + u1 + u2 + + print(f"SA Run {r+1:02d}: UAV1={u1}, UAV2={u2}, " + f"UAV1_cost={c1:.2f}, UAV2_cost={c2:.2f}, Total={total:.2f}") + + sa_costs.append(total) + sa_routes.append(full_route) + sa_u1_list.append(u1) + sa_u2_list.append(u2) + sa_uav_costs.append((c1, c2)) # <-- store both costs + + print("\n🔥 Running SQA (with repair, always feasible)...") + for r in range(NUM_RUNS): + u1, u2, c1, c2, total = sample_and_repair(sqa, Q, N, D) + + full_route = [0] + u1 + u2 + + print(f"SQA Run {r+1:02d}: UAV1={u1}, UAV2={u2}, " + f"UAV1_cost={c1:.2f}, UAV2_cost={c2:.2f}, Total={total:.2f}") + + sqa_costs.append(total) + sqa_routes.append(full_route) + sqa_u1_list.append(u1) + sqa_u2_list.append(u2) + sqa_uav_costs.append((c1, c2)) # <-- store both costs + + sa_results = (sa_costs, sa_routes, sa_u1_list, sa_u2_list, sa_uav_costs) + sqa_results = (sqa_costs, sqa_routes, sqa_u1_list, sqa_u2_list, sqa_uav_costs) + + return sa_results, sqa_results + +# ============================ +# 12. 綜合分析圖表（你指定的版本） +# ============================ +def create_comprehensive_charts(sa_results, sqa_results, classic_results, D, N): + """創建綜合分析圖表""" + fig = plt.figure(figsize=(18, 12)) + + # 解包結果 + sa_costs, sa_routes, sa_uav1, sa_uav2, sa_uav_costs = sa_results + sqa_costs, sqa_routes, sqa_uav1, sqa_uav2, sqa_uav_costs = sqa_results + classic_uav1, classic_uav2, classic_cost, classic_route = classic_results + + # --------------------------------------------------------- + # 1. MTSP 成本分布比較 (左上) - 系統總成本 + # --------------------------------------------------------- + ax1 = plt.subplot(2, 3, 1) + bins = range(int(min(sa_costs + sqa_costs)) - 1, + int(max(sa_costs + sqa_costs)) + 2) + ax1.hist(sa_costs, bins=bins, alpha=0.6, label="SA Total", density=True, color='lightblue') + ax1.hist(sqa_costs, bins=bins, alpha=0.6, label="SQA Total", density=True, color='lightcoral') + + # 添加KDE + if len(sa_costs) > 1: + xs_sa, ys_sa = kde_1d(sa_costs) + ax1.plot(xs_sa, ys_sa, label="SA KDE", color='blue', linewidth=2) + if len(sqa_costs) > 1: + xs_sqa, ys_sqa = kde_1d(sqa_costs) + ax1.plot(xs_sqa, ys_sqa, label="SQA KDE", color='red', linewidth=2) + + ax1.axvline(x=classic_cost, color='green', linestyle='--', + label=f'Classic: {classic_cost}', linewidth=2) + ax1.set_title("Total System Cost Distribution") + ax1.set_xlabel("Total Cost (UAV1 + UAV2)") + ax1.set_ylabel("Density") + ax1.legend() + ax1.grid(True, alpha=0.3) + + # --------------------------------------------------------- + # 2. 單一 UAV 成本分布比較 (右上) - 修正為「個別 UAV」視角 + # --------------------------------------------------------- + ax2 = plt.subplot(2, 3, 2) + + # 展平數據：將 [(c1, c2), (c1, c2)...] 變成 [c1, c2, c1, c2...] + # 這樣可以看到「單台無人機」通常飛多遠 + sa_single_costs = [c for costs in sa_uav_costs for c in costs] + sqa_single_costs = [c for costs in sqa_uav_costs for c in costs] + + # 計算 classic 的單機成本 + c1_classic = uav_cost(classic_uav1, D) + c2_classic = uav_cost(classic_uav2, D) + + # 設定 bins + all_single = sa_single_costs + sqa_single_costs + uav_bins = range(int(min(all_single)) - 1, int(max(all_single)) + 2) + + ax2.hist(sa_single_costs, bins=uav_bins, alpha=0.6, label="SA Single UAV", density=True, color='skyblue') + ax2.hist(sqa_single_costs, bins=uav_bins, alpha=0.6, label="SQA Single UAV", density=True, color='salmon') + + ax2.axvline(x=c1_classic, color='darkgreen', linestyle=':', label='Classic UAV1') + ax2.axvline(x=c2_classic, color='limegreen', linestyle=':', label='Classic UAV2') + + ax2.set_title("Individual UAV Cost Distribution") + ax2.set_xlabel("Cost per UAV") + ax2.set_ylabel("Density") + ax2.legend() + ax2.grid(True, alpha=0.3) + + # --------------------------------------------------------- + # 3. 箱型圖比較 (中上) - 修正重點 + # --------------------------------------------------------- + ax3 = plt.subplot(2, 3, 3) + + # 數據 1 & 2: 系統總成本 (Total Cost) + # 數據 3 & 4: 單機成本 (Individual UAV Cost) - 使用上面展平的數據 + # 這樣 Box Plot 會有明顯的高低落差（總成本高，單機成本低） + + box_data = [sa_costs, sqa_costs, sa_single_costs, sqa_single_costs] + + bp = ax3.boxplot( + box_data, + labels=["SA\nTotal", "SQA\nTotal", "SA\nSingle", "SQA\nSingle"], + patch_artist=True, + showmeans=True + ) + + colors = ['lightblue', 'lightcoral', 'skyblue', 'salmon'] + for patch, color in zip(bp['boxes'], colors): + patch.set_facecolor(color) + + ax3.set_title("Cost Comparison: Total vs Single UAV") + ax3.set_ylabel("Cost") + ax3.grid(True, alpha=0.3) + + # --------------------------------------------------------- + # 4. Route diversity heatmap - SA (左下) + # --------------------------------------------------------- + ax4 = plt.subplot(2, 3, 4) + sa_mat = route_diversity_matrix(sa_routes, N) + im1 = ax4.imshow(sa_mat, aspect='auto', origin='lower', cmap='viridis') + ax4.set_title("SA Route Diversity") + ax4.set_xlabel("Position") + ax4.set_ylabel("City") + plt.colorbar(im1, ax=ax4, shrink=0.8) + + # --------------------------------------------------------- + # 5. Route diversity heatmap - SQA (右下) + # --------------------------------------------------------- + ax5 = plt.subplot(2, 3, 5) + sqa_mat = route_diversity_matrix(sqa_routes, N) + im2 = ax5.imshow(sqa_mat, aspect='auto', origin='lower', cmap='viridis') + ax5.set_title("SQA Route Diversity") + ax5.set_xlabel("Position") + ax5.set_ylabel("City") + plt.colorbar(im2, ax=ax5, shrink=0.8) + + # --------------------------------------------------------- + # 6. 最佳路徑可視化 (中下) + # --------------------------------------------------------- + ax6 = plt.subplot(2, 3, 6) + + # 找出 SA 或 SQA 的最佳解 + best_sa_idx = sa_costs.index(min(sa_costs)) + best_sqa_idx = sqa_costs.index(min(sqa_costs)) + + best_sa_u1 = sa_uav1[best_sa_idx] + best_sa_u2 = sa_uav2[best_sa_idx] + + best_sqa_u1 = sqa_uav1[best_sqa_idx] + best_sqa_u2 = sqa_uav2[best_sqa_idx] + + # 用成本較好的那組 + if min(sa_costs) <= min(sqa_costs): + u1 = best_sa_u1 + u2 = best_sa_u2 + route_name = "SA" + route_cost = min(sa_costs) + else: + u1 = best_sqa_u1 + u2 = best_sqa_u2 + route_name = "SQA" + route_cost = min(sqa_costs) + + # --- 繪製城市 --- + angles = np.linspace(0, 2*np.pi, N, endpoint=False) + x_pos = np.cos(angles) + y_pos = np.sin(angles) + + ax6.scatter(x_pos, y_pos, s=200, c='red', zorder=5) + for i, (x, y) in enumerate(zip(x_pos, y_pos)): + ax6.annotate(str(i), (x, y), ha='center', va='center', + fontsize=10, fontweight='bold', color='white', zorder=6) + + # --- UAV1: 藍色 --- + if len(u1) > 0: + pts = [0] + u1 + [0] # 加入 depot + xs = [x_pos[c] for c in pts] + ys = [y_pos[c] for c in pts] + ax6.plot(xs, ys, '-o', color='blue', linewidth=2, alpha=0.8, + label=f'UAV1 ({len(u1)})') + + # --- UAV2: 橘色 --- + if len(u2) > 0: + pts = [0] + u2 + [0] # 加入 depot + xs = [x_pos[c] for c in pts] + ys = [y_pos[c] for c in pts] + ax6.plot(xs, ys, '-o', color='orange', linewidth=2, alpha=0.8, + label=f'UAV2 ({len(u2)})') + + ax6.set_title(f"Best {route_name} MTSP Solution\nTotal Cost: {route_cost:.2f}") + ax6.set_xlim(-1.3, 1.3) + ax6.set_ylim(-1.3, 1.3) + ax6.set_aspect('equal') + ax6.axis('off') + ax6.legend(loc='upper right') + + plt.tight_layout() + plt.savefig("mtsp_two_uav_comprehensive_analysis.png", dpi=300, bbox_inches='tight') + print(f"\n📊 綜合分析圖表已保存: mtsp_two_uav_comprehensive_analysis.png") + plt.show() + +# ============================ +# 13. main：整合執行 +# ============================ +if __name__ == "__main__": + # 距離矩陣：10 城市 + D = np.array([ + [0, 2, 8, 5, 7, 6, 3, 9, 4, 2], + [2, 0, 5, 7, 3, 8, 4, 6, 9, 1], + [8, 5, 0, 2, 6, 7, 3, 4, 1, 5], + [5, 7, 2, 0, 4, 3, 8, 6, 2, 7], + [7, 3, 6, 4, 0, 5, 9, 2, 7, 3], + [6, 8, 7, 3, 5, 0, 4, 7, 9, 6], + [3, 4, 3, 8, 9, 4, 0, 5, 6, 2], + [9, 6, 4, 6, 2, 7, 5, 0, 8, 3], + [4, 9, 1, 2, 7, 9, 6, 8, 0, 5], + [2, 1, 5, 7, 3, 6, 2, 3, 5, 0] + ]) + N = D.shape[0] + + print("🚀 2-UAV M-TSP (Quantum + Repair vs Classic) 開始運行...\n") + print(f"Distance Matrix ({N}x{N}):") + print(D) + + # 建立 QUBO + Q, _ = build_mtsp_qubo_variable(D, penalty=20.0) + + # SA / SQA，每 run 一定是合法解 + NUM_RUNS = 10 + sa_results, sqa_results = run_mtsp_sa_sqa(Q, N, D, NUM_RUNS=NUM_RUNS) + + # Classic 暴力全域最佳 + print("\n🎯 Running Classic 2-UAV MTSP (Bruteforce ALL solutions)...") + classic_uav1, classic_uav2, classic_total_cost, classic_route = classic_mtsp_two_uav_bruteforce(D) + print(f"Classic UAV1={classic_uav1}, UAV2={classic_uav2}, Cost={classic_total_cost}") + print(f"Classic merged route (0 + perm): {classic_route}") + + # Summary + sa_costs, _, _, _, _ = sa_results + sqa_costs, _, _, _, _ = sqa_results + print("\n📊 Result Summary:") + print(f"SA best cost = {min(sa_costs):.2f}, mean = {np.mean(sa_costs):.2f}") + print(f"SQA best cost = {min(sqa_costs):.2f}, mean = {np.mean(sqa_costs):.2f}") + print(f"Classic global optimum = {classic_total_cost:.2f}") + + # 畫圖 + classic_results = (classic_uav1, classic_uav2, classic_total_cost, classic_route) + print("\n🎨 Creating comprehensive analysis charts...") + create_comprehensive_charts(sa_results, sqa_results, classic_results, D, N) + + print("\n🏁 Done.") diff --git a/sqa_tsp.py b/sqa_tsp.py new file mode 100644 index 0000000..3819b61 --- /dev/null +++ b/sqa_tsp.py @@ -0,0 +1,248 @@ +import numpy as np +import matplotlib.pyplot as plt +from collections import Counter +import math +import openjij as oj + +# ============================ +# 啟發距離矩陣：10-city TSP +# ============================ +D = np.array([ + [0, 2, 8, 5, 7, 6, 3, 9, 4, 2], + [2, 0, 5, 7, 3, 8, 4, 6, 9, 1], + [8, 5, 0, 2, 6, 7, 3, 4, 1, 5], + [5, 7, 2, 0, 4, 3, 8, 6, 2, 7], + [7, 3, 6, 4, 0, 5, 9, 2, 7, 3], + [6, 8, 7, 3, 5, 0, 4, 7, 9, 6], + [3, 4, 3, 8, 9, 4, 0, 5, 6, 2], + [9, 6, 4, 6, 2, 7, 5, 0, 8, 3], + [4, 9, 1, 2, 7, 9, 6, 8, 0, 5], + [2, 1, 5, 7, 3, 6, 2, 3, 5, 0] +]) + +N = D.shape[0] +penalty = 20 # 罰則 > max(D) 即可 + +def idx(i, t, N): + return i * N + t + +# ============================ +# 建立 QUBO：位置編碼 + 固定起點 0 +# ============================ +def build_qubo(D, penalty, fix_start=True): + N = D.shape[0] + Q = {} + + # 距離項 + for t in range(N): + t2 = (t + 1) % N + for i in range(N): + for j in range(N): + if i != j: + u = idx(i, t, N) + v = idx(j, t2, N) + Q[(u, v)] = Q.get((u, v), 0) + D[i, j] + + # 約束A: 每個 time slot 有且只有一個城市 + for t in range(N): + for i in range(N): + u = idx(i, t, N) + Q[(u, u)] = Q.get((u, u), 0) - penalty + for j in range(i + 1, N): + v = idx(j, t, N) + Q[(u, v)] = Q.get((u, v), 0) + 2 * penalty + + # 約束B: 每個城市被訪問一次 + for i in range(N): + for t in range(N): + u = idx(i, t, N) + Q[(u, u)] = Q.get((u, u), 0) - penalty + for t2 in range(t + 1, N): + v = idx(i, t2, N) + Q[(u, v)] = Q.get((u, v), 0) + 2 * penalty + + # 固定起點：t=0 必須是 city 0 + if fix_start: + big = 9999.0 + # 禁止 i!=0 在 t=0 + for i in range(1, N): + u = idx(i, 0, N) + Q[(u, u)] = Q.get((u, u), 0) + big + # 獎勵 city 0 在 t=0 + u0 = idx(0, 0, N) + Q[(u0, u0)] = Q.get((u0, u0), 0) - big + + return Q, N + +# ============================ +# Route decode & cost +# ============================ +def decode_route(sample, N): + route = [] + for t in range(N): + city_for_t = None + for i in range(N): + if sample.get(idx(i, t, N), 0) == 1: + city_for_t = i + break + # 若有問題就塞 -1，方便debug + if city_for_t is None: + city_for_t = -1 + route.append(city_for_t) + return route + +def compute_cost(route, D): + N = len(route) + total = 0 + for k in range(N): + a = route[k] + b = route[(k + 1) % N] + total += D[a, b] + return total + +# ============================ +# 簡單 KDE 實作（Gaussian kernel） +# ============================ +def kde_1d(samples, num_points=200): + xs = np.linspace(min(samples), max(samples), num_points) + n = len(samples) + if n < 2: + return xs, np.zeros_like(xs) + std = np.std(samples) + if std == 0: + std = 1.0 + # Silverman's rule + h = 1.06 * std * (n ** (-1/5)) + if h == 0: + h = 1.0 + ys = [] + inv_sqrt_2pi = 1.0 / math.sqrt(2.0 * math.pi) + for x in xs: + s = 0.0 + for xi in samples: + z = (x - xi) / h + s += math.exp(-0.5 * z * z) * inv_sqrt_2pi + ys.append(s / (n * h)) + return xs, np.array(ys) + +# ============================ +# Route diversity heatmap +# ============================ +def route_diversity_matrix(routes, N): + # freq[i, t] = city i 在位置 t 出現次數 + freq = np.zeros((N, N), dtype=float) + for route in routes: + for t, city in enumerate(route): + if 0 <= city < N: + freq[city, t] += 1 + # normalize to probabilities + if len(routes) > 0: + freq /= len(routes) + return freq + +# ============================ +# 跑 SA or SQA 多次 +# ============================ +def run_algorithm(name, sampler, Q, D, N, num_runs=20, num_reads=20): + all_costs = [] + all_routes = [] + + print(f"\n=== Running {name} ===") + for r in range(num_runs): + result = sampler.sample_qubo(Q, num_reads=num_reads) + best = result.first.sample + route = decode_route(best, N) + cost = compute_cost(route, D) + all_costs.append(cost) + all_routes.append(tuple(route)) + print(f"{name} Run {r+1:02d}: Route={route}, Cost={cost}") + + return all_costs, all_routes + +# ============================ +# Main +# ============================ +if __name__ == "__main__": + Q, N = build_qubo(D, penalty, fix_start=True) + + # Samplers + sqa_sampler = oj.SQASampler() + sa_sampler = oj.SASampler() + + # 可自行調整次數與 num_reads + NUM_RUNS = 20 + NUM_READS = 20 + + sqa_costs, sqa_routes = run_algorithm("SQA", sqa_sampler, Q, D, N, NUM_RUNS, NUM_READS) + sa_costs, sa_routes = run_algorithm("SA", sa_sampler, Q, D, N, NUM_RUNS, NUM_READS) + + # ---------------------------- + # 基本統計 + # ---------------------------- + print("\n=== Summary ===") + print(f"SQA: min={min(sqa_costs)}, max={max(sqa_costs)}, mean={np.mean(sqa_costs):.2f}, std={np.std(sqa_costs):.2f}") + print(f"SA : min={min(sa_costs)}, max={max(sa_costs)}, mean={np.mean(sa_costs):.2f}, std={np.std(sa_costs):.2f}") + + # ============================ + # 圖 1: Cost histogram + KDE + # ============================ + fig1, ax1 = plt.subplots(figsize=(8,5)) + bins = range(min(sa_costs + sqa_costs), max(sa_costs + sqa_costs) + 2) + + ax1.hist(sa_costs, bins=bins, alpha=0.4, label="SA", density=True) + ax1.hist(sqa_costs, bins=bins, alpha=0.4, label="SQA", density=True) + + xs_sa, ys_sa = kde_1d(sa_costs) + xs_sqa, ys_sqa = kde_1d(sqa_costs) + ax1.plot(xs_sa, ys_sa, label="SA KDE") + ax1.plot(xs_sqa, ys_sqa, label="SQA KDE") + + ax1.set_title("Tour Cost Distribution (SA vs SQA)") + ax1.set_xlabel("Tour Cost") + ax1.set_ylabel("Density") + ax1.legend() + fig1.tight_layout() + fig1.savefig("tsp_cost_hist_kde_sa_vs_sqa.png") + + # ============================ + # 圖 2: Violin plot + Box plot + # ============================ + fig2, axs2 = plt.subplots(1, 2, figsize=(10,5)) + + # Violin + axs2[0].violinplot([sa_costs, sqa_costs], positions=[1,2], showmeans=True) + axs2[0].set_xticks([1,2]) + axs2[0].set_xticklabels(["SA", "SQA"]) + axs2[0].set_title("Violin Plot of Tour Costs") + + # Box plot + axs2[1].boxplot([sa_costs, sqa_costs], labels=["SA", "SQA"], showmeans=True) + axs2[1].set_title("Box Plot of Tour Costs") + + fig2.tight_layout() + fig2.savefig("tsp_cost_violin_box_sa_vs_sqa.png") + + # ============================ + # 圖 3: Route diversity heatmap + # ============================ + sqa_mat = route_diversity_matrix(sqa_routes, N) + sa_mat = route_diversity_matrix(sa_routes, N) + + fig3, axs3 = plt.subplots(1, 2, figsize=(12,5)) + + im0 = axs3[0].imshow(sa_mat, aspect='auto', origin='lower') + axs3[0].set_title("Route Diversity (SA)") + axs3[0].set_xlabel("Position t") + axs3[0].set_ylabel("City index") + fig3.colorbar(im0, ax=axs3[0]) + + im1 = axs3[1].imshow(sqa_mat, aspect='auto', origin='lower') + axs3[1].set_title("Route Diversity (SQA)") + axs3[1].set_xlabel("Position t") + axs3[1].set_ylabel("City index") + fig3.colorbar(im1, ax=axs3[1]) + + fig3.tight_layout() + fig3.savefig("tsp_route_diversity_heatmap_sa_vs_sqa.png") + + plt.show() \ No newline at end of file