Openjij/greedy_trap.py

import time
from datetime import datetime
import numpy as np
import networkx as nx
import matplotlib.pyplot as plt
import seaborn as sns
from itertools import permutations
import openjij as oj
import os
import random

# ============ 輸出文件夾設置 (含時間戳) ============
TIMESTAMP = datetime.now().strftime("%Y%m%d_%H%M%S")
OUTPUT_DIR = f"./plots_output/plots_output_{TIMESTAMP}"
if not os.path.exists(OUTPUT_DIR):
    os.makedirs(OUTPUT_DIR)
    print(f"✓ 已建立輸出文件夾: {OUTPUT_DIR}")
# ==================================================

# ==================== 可調整參數區 ====================
# 問題規模設定
CITIES = 15         # 初始核心問題規模 (測試建議先用 12-15)
CITIES_BOUND = 10     # 問題規模的上下範圍 (單一數字控制：以 CITIES 為中心，例如 15±5，即 10~20)
RANDOM = False  
RANDOM_SEED = 42   
COORD_RANGE = (0.0, 10.0)  

# 【核心創新】：環境變異數控制與複雜度壓力測試參數
N_LIST = list(range(CITIES - CITIES_BOUND, CITIES + CITIES_BOUND + 1, 4))  # 自動產生測試的問題規模範圍
COORD_STD = 15.0  # 空間分佈變異數 (越大代表城市分佈越不均勻、越崎嶇)
RISK_STD = 10.0    # 擾動風險變異數 (越大代表某些路段特別危險，H-infinity 衝突極大)

# 演算法執行設定
NUM_RUNS = 1      # 正式比較的執行次數
NUM_READS = 20     # 🌟 每次退火採樣讀取的解數量 (主測試)
NUM_READS_COMPLEXITY = 5  # 複雜度測試的讀取數量 (較少以加快速度)
NUM_READS_TSP = 1  # TSP 測試的讀取數量
BETA = 10     # 🌟 稍微調高 (原為 10.0)，讓低溫結冰得更紮實
SWEEPS_MAIN_TEST = 10000      # 主要比較測試的退火步數
SWEEPS_COMPLEXITY_TEST = 200 # 複雜度擴展測試的退火步數
SWEEPS_HEATMAP_TEST = 200    # 熱力圖測試的退火步數
HEATMAP_RUNS = 5             # 熱力圖測試每種情況的平均次數

# QUBO 參數
PENALTY = 1000.0       # 🌟 從 3000 大幅降回 500 (解開高爾夫球場效應)
BIG_PENALTY = 2000.0  # 🌟 起點約束

# TSP 求解參數
EXACT_LIMIT = 8  # 小於等於這個數量的城市會使用暴力搜尋求解最佳路徑 (確保結果正確性)，超過則使用 nearest neighbor heuristic

# 魯棒優化參數 (H-infinity Robust QUBO)
USE_ROBUST = True  # 強制開啟，這是本論文的核心
GAMMA = 0.5        
SIGMA = 1.0        
ALPHA = 1.0       

# 故障訊號參數 (Fault Signal)
ENABLE_FAULT_SIGNAL = True   # 是否加入故障訊號測試
FAULT_LAMBDA = 50000.0          # 故障懲罰權重 (lambda)，設高一點演算法才會怕
FAULT_PROBABILITY = 0.0     # 兩城市間發生故障/禁飛的機率 (0.1 代表 10% 的路徑斷線)

# 🌟 蜘蛛網參數 (Spider Web Terrain)
SPIDER_WEB_BACKGROUND = 0.3  # 背景值：越小越安全，越大表示背景越危險
SPIDER_WEB_SILK = 0.08         # 蜘蛛絲值：最安全的路徑
SPIDER_WEB_CUT = 0.98          # 切斷點值：最危險的路徑
SPIDER_WEB_SPOKES = 5          # 輻射線數量
SPIDER_WEB_RINGS = 3           # 同心圓數量
SPIDER_WEB_CUT_BLOCKS = 15     # 切斷點數量

# 執行區塊控制
ENABLE_MULTI_UAV = False            # 🌟 是否啟用多機 mTSP 模式 (True=雙機, False=單機 TSP)
RUN_MAIN_TEST = False            # 是否執行主要演算法(SA vs SQA)比較測試
RUN_TSP_TEST = False             # 是否執行獨立單機 TSP (Traveling Salesman Problem) 測試
RUN_5NODE_PERFECT_TRAP = True   # 🌟 [新增] 是否執行 5 節點完美陷阱對比測試 (更複雜的陷阱)
COMPARE_H_INFINITY = False        # 是否比較有無 H-infinity 避障算法的效果 (將輸出 2x3 或 3x2 圖表)
RUN_COMPLEXITY_TEST = False     # 是否執行複雜度擴展(N_LIST)測試
SHOW_TERRAIN_PLOTS = False       # 是否顯示能量地形圖與驗證報告
ENABLE_DISTURBANCE_MATRIX = True  # 是否生成擾動矩陣 (F 矩陣)
ENABLE_QUANTUM_MINEFIELD = True  # 是否產生極端能量障礙與欺騙陷阱 (若為 False，則產生一般隨機地圖)
PLOT_ENV_DIS_FAULT = True        # 是否繪製擾動矩陣與故障矩陣對比圖
# =====================================================

def test_5node_perfect_trap():
    """
    🕷️ 5 節點完美陷阱測試
    演示更複雜的陷阱：SQA 的正解 vs Greedy 的致命觸發器
    
    問題配置：
    - SQA 正解：0→4→3→2→1→0 (安全路徑，無擾動)
    - Greedy 誘餌：0→1→3→... (通過短邊 1-3 上當)
    """
    print("\n" + "="*70)
    print("🕷️ 5 節點完美陷阱測試：SQA vs Greedy")
    print("="*70)
    
    N = 5
    
    # 預設所有路徑都是死亡高牆 (F=100)
    F = np.full((N, N), 100.0)
    # 預設距離為 10
    D = np.full((N, N), 10.0)

    # --- 1. SQA 的正解安全路徑 (0 -> 4 -> 3 -> 2 -> 1 -> 0) ---
    safe_path = [(0, 4), (4, 3), (3, 2), (2, 1), (1, 0)]
    for u, v in safe_path:
        F[u, v] = F[v, u] = 0.0

    # 設定正解的距離
    D[0, 4] = D[4, 0] = 10.0
    D[4, 3] = D[3, 4] = 5.0
    D[3, 2] = D[2, 3] = 10.0
    D[2, 1] = D[1, 2] = 10.0
    D[1, 0] = D[0, 1] = 1.0  # 第一步的共同起點

    # --- 2. 佈置致命捷徑誘餌 (Greedy 專用) ---
    # 在節點 1 和 3 之間開一條 0 擾動、距離超短的捷徑
    F[1, 3] = F[3, 1] = 0.0
    D[1, 3] = D[3, 1] = 1.0  # 致命誘惑！

    np.fill_diagonal(F, 0)
    np.fill_diagonal(D, 0)
    
    # 建立 Fault_Mat：在安全路徑 (4→3) 上放置故障阻礙
    # 逼迫 SQA 展現鯁棒性，找到迴避故障的替代路徑
    Fault_Mat = np.zeros((N, N))
    
    print(f"\n【問題設置】5 節點、1 台 UAV")
    print(f"\n【參數設置】")
    print(f"  α (ALPHA) = {ALPHA}")
    print(f"  σ (SIGMA) = {SIGMA}")
    print(f"  γ (GAMMA) = {GAMMA}")
    disturbance_weight = ALPHA * (SIGMA / (GAMMA**2))
    print(f"  擾動權重係數 = α × (σ/γ²) = {ALPHA} × ({SIGMA}/{GAMMA**2}) = {disturbance_weight:.2f}")
    
    print(f"\n距離矩陣 D:")
    print(D)
    print(f"\n擾動矩陣 F:")
    print(F)
    print(f"\n故障矩陣 Fault_Mat:")
    print(Fault_Mat)
    
    # 構建單機 TSP QUBO
    Q = build_robust_qubo_tsp(
        N, D, F, Fault_Mat=Fault_Mat,
        alpha=ALPHA, gamma=GAMMA, sigma=SIGMA,
        ENABLE_FAULT_SIGNAL=True, lambda_val=FAULT_LAMBDA,
        A=PENALTY  # 使用全局 PENALTY
    )
    
    print(f"\n【QUBO 模型】{len(Q)} 項")
    
    # --- 方案 1: Greedy 演算法 ---
    print("\n[測試 1/3] Greedy 近鄰演算法...")
    # Debug: 顯示Greedy使用的參數
    print(f"  [Debug] 使用參數: α={ALPHA}, σ={SIGMA}, γ={GAMMA}")
    greedy_path, greedy_cost, greedy_time = solve_greedy_tsp(
        N, D, F, Fault_Mat=Fault_Mat,
        alpha=ALPHA, gamma=GAMMA, sigma=SIGMA,
        ENABLE_FAULT_SIGNAL=True, lambda_val=FAULT_LAMBDA
    )
    print(f"  ✓ Greedy 路徑: {greedy_path}")
    print(f"  ✓ 成本: {greedy_cost:.2f} | 耗時: {greedy_time:.6f}s")
    
    # Debug: 驗證是否真的使用ALPHA計算
    disturbance_weight = ALPHA * (SIGMA / (GAMMA**2))
    trap_edge_cost = D[1, 3] + disturbance_weight * (F[1, 3]**2)
    print(f"  [Debug] 陷阱邊 (1,3) 成本 = {D[1,3]:.2f} + {disturbance_weight:.2f} × {F[1,3]:.2f}² = {trap_edge_cost:.2f}")
    
    # 判斷 Greedy 是否掉入陷阱
    if 3 in greedy_path and 1 in greedy_path:
        # 檢查是否經過短邊 (1,3)
        path_str = str(greedy_path)
        if ('1' in path_str and '3' in path_str) or ('3' in path_str and '1' in path_str):
            print(f"  ⚠️ Greedy 已掉入陷阱（被短邊 1→3 吸引）")
        else:
            print(f"  ✓ Greedy 找到安全路徑")
    else:
        print(f"  ✓ Greedy 找到安全路徑（未掉入陷阱）")
    
    # --- 方案 2: SA 求解 ---
    print("\n[測試 2/3] Simulated Annealing (SA) 求解...")
    Q_norm = {k: v / PENALTY for k, v in Q.items()}  # 正規化
    sa_sampler = oj.SASampler()
    
    sa_results = []
    for run_idx in range(3):
        response_sa = sa_sampler.sample_qubo(
            Q_norm,
            num_reads=NUM_READS_TSP,
            num_sweeps=SWEEPS_MAIN_TEST,
            beta_min=0.1,
            beta_max=BETA
        )
        best_sample = response_sa.first.sample
        sa_path = decode_solution(best_sample, N)
        
        if sa_path is not None:
            sa_cost = calculate_cost(
                sa_path, D, F, Fault_Mat=Fault_Mat,
                alpha=ALPHA, gamma=GAMMA, sigma=SIGMA,
                ENABLE_FAULT_SIGNAL=True, lambda_val=FAULT_LAMBDA
            )
            sa_results.append({'path': sa_path, 'cost': sa_cost, 'run': run_idx+1})
            # Debug: 如果是第一次運行，顯示α值
            if run_idx == 0:
                print(f"  [Debug] SA 運行時 α={ALPHA}")
            print(f"  Run {run_idx+1}: 路徑={sa_path}, 成本={sa_cost:.2f}")
    
    # --- 方案 3: SQA 求解 ---
    print("\n[測試 3/3] Simulated Quantum Annealing (SQA) 求解...")
    sqa_sampler = oj.SQASampler()
    
    # 建立 SQA 的退火排程
    TOTAL_SWEEPS = SWEEPS_MAIN_TEST
    smooth_sched = get_smooth_sqa_schedule(beta=BETA, total_sweeps=TOTAL_SWEEPS, num_steps=50)
    
    sqa_results = []
    for run_idx in range(3):
        response_sqa = sqa_sampler.sample_qubo(
            Q,
            schedule=smooth_sched,
            num_reads=NUM_READS_TSP
        )
        best_sample = response_sqa.first.sample
        sqa_path = decode_solution(best_sample, N)
        
        if sqa_path is not None:
            sqa_cost = calculate_cost(
                sqa_path, D, F, Fault_Mat=Fault_Mat,
                alpha=ALPHA, gamma=GAMMA, sigma=SIGMA,
                ENABLE_FAULT_SIGNAL=True, lambda_val=FAULT_LAMBDA
            )
            sqa_results.append({'path': sqa_path, 'cost': sqa_cost, 'run': run_idx+1})
            # Debug: 如果是第一次運行，顯示α值
            if run_idx == 0:
                print(f"  [Debug] SQA 運行時 α={ALPHA}")
            print(f"  Run {run_idx+1}: 路徑={sqa_path}, 成本={sqa_cost:.2f}")
        else:
            print(f"  Run {run_idx+1}: [警告] decode_solution 返回 None")
    
    # --- 結果彙總 ---
    print("\n" + "="*70)
    print("📊 5 節點陷阱測試結果彙總")
    print("="*70)
    
    # 最佳解
    best_greedy = {'path': greedy_path, 'cost': greedy_cost}
    best_sa = min(sa_results, key=lambda x: x['cost']) if sa_results else None
    best_sqa = min(sqa_results, key=lambda x: x['cost']) if sqa_results else None
    
    print(f"\nGreedy:   路徑={best_greedy['path']}, 成本={best_greedy['cost']:.2f}")
    if best_sa:
        print(f"SA:       路徑={best_sa['path']}, 成本={best_sa['cost']:.2f}")
    if best_sqa:
        print(f"SQA:      路徑={best_sqa['path']}, 成本={best_sqa['cost']:.2f}")
    
    # 計算改進
    if best_sa:
        sa_improvement = (greedy_cost - best_sa['cost']) / greedy_cost * 100
        print(f"\nSA 相對 Greedy 的改進: {sa_improvement:.1f}%")
    
    if best_sqa:
        sqa_improvement = (greedy_cost - best_sqa['cost']) / greedy_cost * 100
        print(f"SQA 相對 Greedy 的改進: {sqa_improvement:.1f}%")
    
    print("="*70 + "\n")
    
    # 🌟 新增：繪製算法成本對比柱狀圖
    if PLOT_ENV_DIS_FAULT:
        algorithms = ['Greedy', 'SA', 'SQA']
        costs = [greedy_cost]
        colors = ['#ff6b6b', '#4ecdc4', '#45b7d1']
        
        if best_sa:
            costs.append(best_sa['cost'])
        else:
            costs.append(0)
        
        if best_sqa:
            costs.append(best_sqa['cost'])
        else:
            costs.append(0)
        
        fig, ax = plt.subplots(figsize=(10, 6))
        bars = ax.bar(algorithms, costs, color=colors, edgecolor='black', linewidth=2, alpha=0.8)
        
        # 在柱子上方標記成本值
        for bar, cost in zip(bars, costs):
            height = bar.get_height()
            if cost > 0:
                if cost > 100000:
                    label = f'{cost:.0f}'
                else:
                    label = f'{cost:.2f}'
                ax.text(bar.get_x() + bar.get_width()/2., height,
                       label,
                       ha='center', va='bottom', fontsize=14, fontweight='bold')
        
        ax.set_ylabel('Cost', fontsize=14, fontweight='bold')
        ax.set_title('5-Node Perfect Trap: Algorithm Cost Comparison', fontsize=16, fontweight='bold')
        ax.grid(axis='y', alpha=0.3, linestyle='--')
        
        plt.tight_layout()
        plt.savefig(os.path.join(OUTPUT_DIR, "trap_test_cost_comparison.png"), dpi=300)
        print(f"🎨 算法成本對比圖已儲存: {os.path.join(OUTPUT_DIR, 'trap_test_cost_comparison.png')}")
    
    # 🌟 新增：繪製擾動與故障矩陣對比圖
    if PLOT_ENV_DIS_FAULT:
        # 使用測試過程中實際使用的 Fault_Mat（不要重新創建）
        plot_disturbance_and_fault_matrices(D, F, Fault_Mat)
    
    return {
        'greedy': best_greedy,
        'sa': best_sa,
        'sqa': best_sqa,
    }


def generate_quantum_minefield(N, random_seed=42):
    """
    🕷️ 蜘蛛網地貌生成器：直接繪製版 (Direct Drawing)
    """
    # ========== 蜘蛛網參數設定（使用全局參數） ==========
    num_spokes = SPIDER_WEB_SPOKES
    num_rings = SPIDER_WEB_RINGS
    num_anchor_threads = "all"
    background_val = SPIDER_WEB_BACKGROUND
    silk_val = SPIDER_WEB_SILK
    cut_val = SPIDER_WEB_CUT
    num_cut_blocks = SPIDER_WEB_CUT_BLOCKS
    cut_block_width = 2

    rng = np.random.default_rng(random_seed)
    M = np.full((N, N), background_val, dtype=float)

    cx = N / 2 + rng.uniform(-0.8, 0.8)
    cy = N / 2 + rng.uniform(-0.8, 0.8)
    base_radius = 0.34 * N

    # ========== 繪圖輔助函數 ==========
    def set_pixel(arr, x, y, value, mode='min'):
        xi, yi = int(round(x)), int(round(y))
        if 0 <= xi < arr.shape[1] and 0 <= yi < arr.shape[0]:
            if mode == 'min':
                arr[yi, xi] = min(arr[yi, xi], value)
            else:
                arr[yi, xi] = max(arr[yi, xi], value)

    def draw_line_4connected(arr, x0, y0, x1, y1, value, mode='min'):
        x = int(round(x0))
        y = int(round(y0))
        x1 = int(round(x1))
        y1 = int(round(y1))
        set_pixel(arr, x, y, value, mode=mode)
        dx = x1 - x
        dy = y1 - y
        sx = 1 if dx > 0 else -1 if dx < 0 else 0
        sy = 1 if dy > 0 else -1 if dy < 0 else 0
        dx = abs(dx)
        dy = abs(dy)
        err = 0
        if dx >= dy:
            while x != x1 or y != y1:
                moved = False
                if x != x1:
                    x += sx
                    set_pixel(arr, x, y, value, mode=mode)
                    moved = True
                err += dy
                while err >= dx and y != y1 and dx != 0:
                    y += sy
                    set_pixel(arr, x, y, value, mode=mode)
                    err -= dx
                    moved = True
                if not moved and y != y1:
                    y += sy
                    set_pixel(arr, x, y, value, mode=mode)
        else:
            while x != x1 or y != y1:
                moved = False
                if y != y1:
                    y += sy
                    set_pixel(arr, x, y, value, mode=mode)
                    moved = True
                err += dx
                while err >= dy and x != x1 and dy != 0:
                    x += sx
                    set_pixel(arr, x, y, value, mode=mode)
                    err -= dy
                    moved = True
                if not moved and x != x1:
                    x += sx
                    set_pixel(arr, x, y, value, mode=mode)

    def draw_polyline_4connected(arr, pts, value, mode='min'):
        for (x0, y0), (x1, y1) in zip(pts[:-1], pts[1:]):
            draw_line_4connected(arr, x0, y0, x1, y1, value, mode=mode)

    def boundary_intersection_from_center(cx, cy, angle, W, H):
        dx = np.cos(angle)
        dy = np.sin(angle)
        candidates = []
        if abs(dx) > 1e-9:
            t = (0 - cx) / dx
            y = cy + t * dy
            if t > 0 and 0 <= y <= H - 1:
                candidates.append((t, 0, y))
            t = ((W - 1) - cx) / dx
            y = cy + t * dy
            if t > 0 and 0 <= y <= H - 1:
                candidates.append((t, W - 1, y))
        if abs(dy) > 1e-9:
            t = (0 - cy) / dy
            x = cx + t * dx
            if t > 0 and 0 <= x <= W - 1:
                candidates.append((t, x, 0))
            t = ((H - 1) - cy) / dy
            x = cx + t * dx
            if t > 0 and 0 <= x <= W - 1:
                candidates.append((t, x, H - 1))
        candidates.sort(key=lambda z: z[0])
        _, x, y = candidates[0]
        return x, y

    base_angles = np.linspace(0, 2 * np.pi, num_spokes, endpoint=False)
    angles = np.sort(base_angles + rng.normal(0, 0.08, num_spokes))

    ring_fracs = np.linspace(0.60, 1.0, num_rings)  # 🌟 改為 0.60，讓中心區域更小
    ring_fracs += rng.normal(0, 0.02, num_rings)
    ring_fracs = np.sort(np.clip(ring_fracs, 0.55, 1.0))

    nodes = []
    for rf in ring_fracs:
        ring_nodes = []
        for ang in angles:
            local_r = base_radius * rf * (1 + rng.normal(0, 0.03))
            wobble = 0.35 * np.sin(2 * ang + 3 * rf) + rng.normal(0, 0.12)
            x = cx + (local_r + wobble) * np.cos(ang)
            y = cy + (local_r + wobble) * np.sin(ang)
            ring_nodes.append((x, y))
        nodes.append(ring_nodes)

    for j in range(num_spokes):
        pts = [(cx, cy)] + [nodes[i][j] for i in range(num_rings)]
        draw_polyline_4connected(M, pts, silk_val, mode='min')

    for i in range(num_rings):
        ring_pts = nodes[i] + [nodes[i][0]]
        draw_polyline_4connected(M, ring_pts, silk_val, mode='min')

    outer_nodes = nodes[-1]
    anchor_indices = np.arange(num_spokes)
    for idx in anchor_indices:
        x0, y0 = outer_nodes[idx]
        angle = np.arctan2(y0 - cy, x0 - cx)
        xb, yb = boundary_intersection_from_center(cx, cy, angle, N, N)
        draw_line_4connected(M, x0, y0, xb, yb, silk_val, mode='min')

    silk_pixels = np.argwhere(M <= silk_val + 1e-9)
    chosen = []
    min_separation = 6.0
    for idx in rng.permutation(len(silk_pixels)):
        y, x = silk_pixels[idx]
        if all((x - px) ** 2 + (y - py) ** 2 >= min_separation ** 2 for py, px in chosen):
            chosen.append((y, x))
        if len(chosen) >= num_cut_blocks:
            break

    for y, x in chosen:
        x0 = max(0, x)
        x1 = min(N, x + cut_block_width)
        y0 = max(0, y)
        y1 = min(N, y + cut_block_width)
        M[y0:y1, x0:x1] = np.maximum(M[y0:y1, x0:x1], cut_val)

    # ========== 轉換為 QUBO 相容格式 ==========
    # 🌟 改進：背景的 F 值直接受 background_val 驅動
    # background_val 越大 -> F 值越大（越危險）
    # background_val 越小 -> F 值越小（越安全）
    
    F = np.zeros((N, N), dtype=float)
    for i in range(N):
        for j in range(N):
            if M[i, j] <= silk_val + 1e-9:
                F[i, j] = 0.0                  # 蜘蛛絲：安全
            elif M[i, j] >= cut_val - 1e-9:
                F[i, j] = 10.0                 # 切斷點：最危險
            else:
                # 中間區域（背景）: 0.08 < M < 0.98
                # F 值正比於 background_val，讓參數改動有實際效果
                # background_val 小 (如 0.1) -> F 最大約 8.0
                # background_val 大 (如 0.5) -> F 最大約 8.0
                # 但要讓改變更明顯，用 background_val 直接決定最大值
                f_max_for_background = 10.0 * background_val / cut_val  # 這樣 background_val 變化會直接影響 F
                F[i, j] = f_max_for_background * (M[i, j] - silk_val) / (background_val - silk_val)

    # 城市坐標
    city_coords = []
    all_node_positions = [(cx, cy)] + [pos for ring in nodes for pos in ring]
    for i in range(N):
        if i < len(all_node_positions):
            city_coords.append(all_node_positions[i])
        else:
            x = rng.uniform(0, N)
            y = rng.uniform(0, N)
            city_coords.append((x, y))

    city_coords = np.array(city_coords)
    D = np.zeros((N, N), dtype=float)
    for i in range(N):
        for j in range(i + 1, N):
            dist = np.linalg.norm(city_coords[i] - city_coords[j])
            D[i, j] = dist
            D[j, i] = dist
    
    return D, F

def verify_minefield_stats(D, F):
    """
    數值驗證：計算「黃金捷徑」與「一般隨機路徑」的成本差異，
    證明地圖中確實存在極端的高牆與深谷。
    """
    N = len(D)
    
    # 1. 計算黃金捷徑 (0 -> 1 -> 2 -> ... -> N-1 -> 0) 的總成本
    golden_path = list(range(N))
    golden_cost = 0
    for i in range(N):
        u = golden_path[i]
        v = golden_path[(i + 1) % N]
        golden_cost += (D[u, v] + F[u, v])
        
    # 2. 計算一條「不小心踩到高牆」的隨機路徑成本
    np.random.seed(99)
    random_path = np.random.permutation(N).tolist()
    random_cost = 0
    for i in range(N):
        u = random_path[i]
        v = random_path[(i + 1) % N]
        random_cost += (D[u, v] + F[u, v])

    print("\n🔍 [雷區地貌驗證報告]")
    print(f"✅ 黃金隧道總成本: {golden_cost:.2f} (演算法的終極目標)")
    print(f"❌ 隨機踩雷總成本: {random_cost:.2f} (一般 SA 容易落入的下場)")
    print(f"⚖️ 難度倍率: 隨機路徑的成本是黃金捷徑的 {random_cost/golden_cost:.1f} 倍！\n")

def plot_energy_landscape_heatmap(D, F):
    """
    視覺化驗證：繪製 D+F 總成本矩陣的熱力圖。
    """
    Total_Cost_Matrix = D + F
    
    plt.figure(figsize=(8, 6))
    # 使用 'hot' 顏色地圖，顏色越亮(黃/白)代表成本越高(高牆)，越暗(黑/紅)代表成本越低(深谷)
    plt.imshow(Total_Cost_Matrix, cmap='hot', interpolation='nearest')
    plt.colorbar(label='Total Energy Cost (D + F)')
    plt.title('Quantum Minefield Energy Landscape')
    plt.xlabel('To City')
    plt.ylabel('From City')
    
    # 標示出黃金路線 (次對角線)
    plt.tight_layout()
    plt.savefig(os.path.join(OUTPUT_DIR, "energy_landscape_heatmap.png"), dpi=300)
    print(f"🎨 能量地形圖已儲存: {os.path.join(OUTPUT_DIR, 'energy_landscape_heatmap.png')}")
    # plt.show() 將在最後統一呼叫

def plot_energy_landscape_3d(D, F):
    """
    視覺化驗證二：繪製 D+F 總成本矩陣的 3D 能量地形圖。
    """
    Total_Cost_Matrix = D + F
    N = len(D)
    
    X, Y = np.meshgrid(range(N), range(N))
    
    fig = plt.figure(figsize=(10, 8))
    ax = fig.add_subplot(111, projection='3d')
    
    # 繪製 3D 表面
    surf = ax.plot_surface(X, Y, Total_Cost_Matrix, cmap='hot', edgecolor='none', alpha=0.9)
    fig.colorbar(surf, label='Total Energy Cost (D + F)', shrink=0.5, aspect=5)
    
    ax.set_title('3D Quantum Minefield Energy Landscape')
    ax.set_xlabel('To City')
    ax.set_ylabel('From City')
    ax.set_zlabel('Cost')
    
    plt.tight_layout()
    plt.savefig(os.path.join(OUTPUT_DIR, "energy_landscape_3d.png"), dpi=300)
    print(f"🎨 3D 能量地形圖已儲存: {os.path.join(OUTPUT_DIR, 'energy_landscape_3d.png')}")
    # plt.show() 將在最後統一呼叫

# ----------------- 基礎輔助函數 -----------------
def generate_distance_matrix(num_cities, random=True, seed=None, coord_range=None):
    if seed is None: seed = RANDOM_SEED
    if coord_range is None: coord_range = COORD_RANGE
    if not random: np.random.seed(seed)
    else: np.random.seed(None)
    low, high = coord_range
    coords = np.random.uniform(low, high, size=(num_cities, 2))
    D = np.zeros((num_cities, num_cities), dtype=float)
    for i in range(num_cities):
        for j in range(i + 1, num_cities):
            dist = np.linalg.norm(coords[i] - coords[j])
            D[i, j] = dist
            D[j, i] = dist
    return np.round(D, 2), coords

def calculate_distance_matrix(city_coords):
    """
    根據給定的城市座標計算距離矩陣 D
    
    Args:
        city_coords: shape (N, 2) 的座標陣列
    
    Returns:
        D: shape (N, N) 的距離矩陣
    """
    num_cities = len(city_coords)
    D = np.zeros((num_cities, num_cities), dtype=float)
    for i in range(num_cities):
        for j in range(i + 1, num_cities):
            dist = np.linalg.norm(city_coords[i] - city_coords[j])
            D[i, j] = dist
            D[j, i] = dist
    return np.round(D, 2)

def generate_disturbance_matrix(num_cities, seed=None):
    """
    將蜘蛛網拓樸結構純粹應用在擾動矩陣 F 上。
    
    邏輯：
      - 所有「沒有蜘蛛絲連接」的路段都極度危險 (F=10.0)
      - 有「蜘蛛絲」連接的路段都是安全的 (F=0.0)
      - 蜘蛛網由輻射線 (Spokes) 和同心圓 (Rings) 組成
    """
    if seed is None: 
        seed = RANDOM_SEED
    np.random.seed(seed)
    
    N = num_cities
    
    # 預設：所有「沒有蜘蛛絲」的路徑都是極度危險的 (F=10.0)
    F = np.full((N, N), 10.0) 
    np.fill_diagonal(F, 0)

    if N <= 3:
        return np.zeros((N, N))

    # 動態計算蜘蛛網的形狀 (例如 N=31 -> 大約 5條輻射線, 6圈)
    spokes = max(3, int(np.sqrt(N - 1)))
    rings = (N - 1) // spokes

    valid_edges = set()

    # 1. 中心 Depot (0) 連到第一圈的城市
    for s in range(1, spokes + 1):
        if s < N:
            valid_edges.add((0, s))

    # 2. 輻射線 (向外延伸的蜘蛛絲)
    for r in range(rings - 1):
        for s in range(1, spokes + 1):
            curr_node = 1 + r * spokes + (s - 1)
            next_node = curr_node + spokes 
            if next_node < N:
                valid_edges.add((curr_node, next_node))

    # 3. 同心圓 (橫向連接的蜘蛛絲)
    for r in range(rings):
        for s in range(1, spokes + 1):
            curr_node = 1 + r * spokes + (s - 1)
            next_s = s + 1 if s <= spokes else 1
            next_node = 1 + r * spokes + (next_s - 1)
            
            if curr_node < N and next_node < N and s < spokes:
                 valid_edges.add((curr_node, next_node))
            elif curr_node < N and s == spokes: 
                 first_node_in_ring = 1 + r * spokes
                 if first_node_in_ring < N:
                     valid_edges.add((curr_node, first_node_in_ring))

    # 4. 處理無法整除的「剩餘邊緣城市」
    for i in range(1 + rings * spokes, N):
        random_outer = np.random.randint(1 + (rings-1)*spokes, 1 + rings*spokes)
        valid_edges.add((random_outer, i))

    # 🌟 5. 將有蜘蛛絲連接的邊，設為安全路徑 (F = 0)
    for u, v in valid_edges:
        F[u, v] = 0.0
        F[v, u] = 0.0

    return F

def generate_fault_matrix(num_cities, prob=0.1, seed=None):
    """產生 0 或 1 的故障矩陣 I_fault(i,j)"""
    if seed is not None: 
        np.random.seed(seed)
    # 產生 0 或 1 的矩陣
    Fault_Mat = np.random.choice([0, 1], size=(num_cities, num_cities), p=[1-prob, prob])
    # 確保矩陣對稱 (i 到 j 故障，j 到 i 也故障)
    Fault_Mat = np.maximum(Fault_Mat, Fault_Mat.T)
    np.fill_diagonal(Fault_Mat, 0.0)
    return Fault_Mat

def generate_spider_web_system(N, num_obstacles=8):
    """
    這會產生一個真正的蜘蛛網拓樸：
    F = 0 代表蜘蛛絲路徑
    F = 10 代表網格外的危險區
    Obstacle_Mat = 1 代表蜘蛛絲斷掉的地方
    """
    # 初始化：全域都是危險區 (F=10)
    F = np.full((N, N), 10.0)
    np.fill_diagonal(F, 0)
    Obstacle_Mat = np.zeros((N, N))

    # 計算結構：中心(1) + 圈數(rings) * 輻射線(spokes)
    spokes = max(3, int(np.sqrt(N - 1)))
    rings = (N - 1) // spokes
    
    web_edges = []

    # A. 建立蜘蛛網連線邏輯
    # 1. 中心連向第一圈
    for s in range(1, spokes + 1):
        if s < N: web_edges.append((0, s))
    # 2. 輻射線 (向外延伸)
    for r in range(rings - 1):
        for s in range(1, spokes + 1):
            curr = 1 + r * spokes + (s - 1)
            nxt = curr + spokes
            if nxt < N: web_edges.append((curr, nxt))
    # 3. 同心圓 (環狀連接)
    for r in range(rings):
        for s in range(1, spokes + 1):
            curr = 1 + r * spokes + (s - 1)
            next_s = (s % spokes) + 1
            nxt = 1 + r * spokes + (next_s - 1)
            if curr < N and nxt < N: web_edges.append((curr, nxt))

    # B. 套用到 F 矩陣：有連線的地方 F=0
    for u, v in web_edges:
        F[u, v] = F[v, u] = 0.0

    # C. 在這些「蜘蛛絲」上隨機抽樣放障礙物
    if len(web_edges) > num_obstacles:
        obs_links = random.sample(web_edges, num_obstacles)
        for u, v in obs_links:
            Obstacle_Mat[u, v] = Obstacle_Mat[v, u] = 1.0

    return F, Obstacle_Mat

def idx_mtsp(k, i, p, N):
    return k * (N * N) + i * N + p

def build_robust_qubo(D, F, Fault_Mat=None, gamma=GAMMA, sigma=SIGMA, alpha=ALPHA, penalty=PENALTY, big_penalty=BIG_PENALTY, lambda_val=FAULT_LAMBDA):
    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

    # 🌟 根據ENABLE_MULTI_UAV決定機數
    num_uav = 2 if ENABLE_MULTI_UAV else 1
    
    for k in range(num_uav):
        for p in range(N):
            q = (p + 1) % N
            for i in range(N):
                for j in range(N):
                    dij, fij = D[i, j], F[i, j]
                    if dij == 0: continue
                    # 這裡就是 H-infinity 的核心風險項
                    risk_term = (sigma / (gamma**2)) * (fij**2)
                    total_weight = dij + (alpha * risk_term)
                    
                    # 🌟 新增：加入故障訊號懲罰
                    if ENABLE_FAULT_SIGNAL and Fault_Mat is not None:
                        total_weight += lambda_val * Fault_Mat[i, j]
                        
                    u, v = idx_mtsp(k, i, p, N), idx_mtsp(k, j, q, N)
                    addQ(u, v, total_weight)

    # 1. 每個時間點 p 只能去一個城市
    for k in range(num_uav):
        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)

    # 2. 每個城市 i (除了起點) 都必須被拜訪剛好一次
    for i in range(1, N):
        if ENABLE_MULTI_UAV:
            vars_city = [idx_mtsp(k, i, p, N) for k in range(num_uav) for p in range(N)]
        else:
            vars_city = [idx_mtsp(0, i, p, N) for p in range(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)

    # 3. 起點約束 (第一步與最後一步必須是 0)
    for k in range(num_uav):
        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

# ---------------------------------------------------------
# 🌟 [新增] 貪婪演算法解決 Robust mTSP
# ---------------------------------------------------------
def solve_greedy_mtsp(D, F, Fault_Mat=None, gamma=GAMMA, sigma=SIGMA, alpha=ALPHA, lambda_val=FAULT_LAMBDA):
    """
    貪婪演算法 (Nearest Neighbor) 解決 Robust mTSP
    邏輯：
      - ENABLE_MULTI_UAV=True: 兩台 UAV 輪流選最便宜的城市
      - ENABLE_MULTI_UAV=False: 單台 UAV 飛遍所有城市（回傳 p1, []）
    """
    N = D.shape[0]
    unvisited = set(range(1, N))
    
    # 計算兩點之間的綜合成本 (與 QUBO 的權重公式一致)
    def edge_cost(u, v):
        """綜合成本：用於演算法決策"""
        if u == v: return float('inf')
        w = D[u, v] + alpha * (sigma / (gamma**2)) * (F[u, v]**2)
        if ENABLE_FAULT_SIGNAL and Fault_Mat is not None:
            w += lambda_val * Fault_Mat[u, v]
        return w

    if ENABLE_MULTI_UAV:
        # 🌟 多機版本：兩台 UAV 動態分配
        p1, p2 = [], []
        curr1, curr2 = 0, 0
        cost1, cost2 = 0.0, 0.0  # 綜合成本（用於決策）
        dist_cost1, dist_cost2 = 0.0, 0.0  # 純距離成本（用於 Makespan）
        
        while unvisited:
            # 決定換哪台無人機飛 (選目前綜合成本較低的，以平衡 Makespan)
            if cost1 <= cost2:
                # 找對 UAV 1 來說最便宜的下一個城市
                next_node = min(unvisited, key=lambda v: edge_cost(curr1, v))
                cost1 += edge_cost(curr1, next_node)
                dist_cost1 += D[curr1, next_node]  # 只加距離
                curr1 = next_node
                p1.append(next_node)
                unvisited.remove(next_node)
            else:
                # 找對 UAV 2 來說最便宜的下一個城市
                next_node = min(unvisited, key=lambda v: edge_cost(curr2, v))
                cost2 += edge_cost(curr2, next_node)
                dist_cost2 += D[curr2, next_node]  # 只加距離
                curr2 = next_node
                p2.append(next_node)
                unvisited.remove(next_node)
                
        # 最後兩台都要飛回起點 (Depot 0) - 只計算距離
        dist_cost1 += D[curr1, 0]
        dist_cost2 += D[curr2, 0]
        
        # Makespan = 兩台的最大飛行距離
        makespan = max(dist_cost1, dist_cost2)
        return p1, p2, makespan
    
    else:
        # 🌟 單機版本：一台 UAV 飛遍所有城市
        p1 = []
        curr1 = 0
        dist_cost1 = 0.0  # 純距離成本
        
        while unvisited:
            # UAV 1 一個人包辦所有城市
            next_node = min(unvisited, key=lambda v: edge_cost(curr1, v))
            dist_cost1 += D[curr1, next_node]  # 只加距離
            curr1 = next_node
            p1.append(next_node)
            unvisited.remove(next_node)
                
        # 最後飛回起點 (Depot 0)
        dist_cost1 += D[curr1, 0]
        
        # Makespan => 單機 TSP 的總距離
        makespan = dist_cost1
        
        # 🌟 回傳 p1, [], makespan
        # 空陣列 [] 會被當作 p2，後續所有程式碼都會自動把它的成本計為 0
        return p1, [], makespan

# ---------------------------------------------------------
# 🌟 [新增] 產生平滑的 SQA Schedule
# ---------------------------------------------------------
def get_smooth_sqa_schedule(beta, total_sweeps, num_steps=100):
    """產生標準 SQA 的平滑退火排程 (s 從 0.0 遞增到 1.0)"""
    schedule = []
    sweeps_per_step = max(1, int(total_sweeps / num_steps))
    for s in np.linspace(0.0, 1.0, num_steps):
        schedule.append([float(s), float(beta), int(sweeps_per_step)])
    return schedule

def decode_slots(sample, N):
    """
    根據 ENABLE_MULTI_UAV 設定解碼 UAV 的 slots。
    
    - 多機模式 (ENABLE_MULTI_UAV=True): 解碼兩台 UAV，回傳 (slots[0], slots[1])
    - 單機模式 (ENABLE_MULTI_UAV=False): 只解碼一台 UAV，回傳 (slots[0], [])
    """
    num_uav = 2 if ENABLE_MULTI_UAV else 1
    slots = []
    for k in range(num_uav):
        row = []
        for p in range(N):
            chosen = 0
            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)
    
    if num_uav == 1:
        return slots[0], []  # 單機：u2 真的是空的
    return slots[0], slots[1]  # 多機：正確回傳兩台的 slots

def decode_single_uav(sample, N):
    """
    單機模式：直接從 slot 順序讀出路徑
    
    QUBO 變數 x[i, p] 的意思是「城市 i 在第 p 步被拜訪」，
    順序本身就是解，不需要 repair 再重排。
    
    Args:
        sample: QUBO 採樣結果
        N: 城市數
    
    Returns:
        path: 訪問順序 [city1, city2, ..., cityN-1] (不含depot 0)
    """
    path = []
    for p in range(N):
        for i in range(N):
            if sample.get(idx_mtsp(0, i, p, N), 0) == 1:
                if i != 0:          # 跳過 depot
                    path.append(i)
                break
    
    # 若有城市沒出現（非法解），補上缺少的城市
    missing = set(range(1, N)) - set(path)
    path += sorted(missing)
    return path

def repair_routes_from_slots(u1_slots, u2_slots, N, D):
    count1, count2 = [0]*N, [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):
        if count1[city] > count2[city]: 
            assign1.append(city)
        elif count2[city] > count1[city]: 
            assign2.append(city)
        else:
            # 平票：計算該城市到兩個 UAV 分配城市群的平均距離
            if assign1:
                avg_dist1 = np.mean([D[city, c] for c in assign1])
            else:
                avg_dist1 = float('inf')
            
            if assign2:
                avg_dist2 = np.mean([D[city, c] for c in assign2])
            else:
                avg_dist2 = float('inf')
            
            # 選擇距離較近的 UAV
            if avg_dist1 <= avg_dist2: 
                assign1.append(city)
            else: 
                assign2.append(city)
    return sorted(assign1), sorted(assign2)

def uav_cost(path, D):
    if not path or 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 uav_disturbance_energy(path, F, gamma=GAMMA, sigma=SIGMA, alpha=ALPHA):
    if not path or len(path) == 0: return 0.0
    risk = alpha * (sigma / (gamma**2)) * (F[0, path[0]]**2)
    for i in range(len(path)-1): risk += alpha * (sigma / (gamma**2)) * (F[path[i], path[i+1]]**2)
    risk += alpha * (sigma / (gamma**2)) * (F[path[-1], 0]**2)
    return float(risk)

def best_order_for_cities(cities, D, exact_limit=EXACT_LIMIT):
    cities = list(cities)
    if len(cities) <= 1: return cities, uav_cost(cities, D)
    if len(cities) <= exact_limit:
        best_perm, best_cost = None, float('inf')
        for perm in permutations(cities):
            c = uav_cost(list(perm), D)
            if c < best_cost: best_cost, best_perm = c, list(perm)
        return best_perm, best_cost
    else: # 簡單的 nearest neighbor heuristic 當作 fallback
        rem = set(cities)
        curr = min(rem, key=lambda x: D[0, x])
        rem.remove(curr)
        route = [curr]
        while rem:
            nxt = min(rem, key=lambda x: D[curr, x])
            route.append(nxt)
            rem.remove(nxt)
            curr = nxt
        return route, uav_cost(route, D)

def get_makespan_and_risk(sample, N, D, F):
    """
    根據 ENABLE_MULTI_UAV 決定解碼方式：
    - 單機模式：直接讀出路徑，無需 repair 和 best_order
    - 多機模式：decode slots -> repair -> best_order
    """
    if not ENABLE_MULTI_UAV:
        # 單機模式：直接從時間步順序解碼
        p1 = decode_single_uav(sample, N)
        c1 = uav_cost(p1, D)
        dist_energy = uav_disturbance_energy(p1, F)
        return c1, dist_energy, p1, []
    else:
        # 多機模式：需要進行 repair 和最佳排序
        u1, u2 = decode_slots(sample, N)
        a1, a2 = repair_routes_from_slots(u1, u2, N, D)
        p1, c1 = best_order_for_cities(a1, D)
        p2, c2 = best_order_for_cities(a2, D)
        makespan = max(c1, c2)
        dist_energy = uav_disturbance_energy(p1, F) + uav_disturbance_energy(p2, F)
        return makespan, dist_energy, p1, p2

# ====================================================================
# 🌟 [新增] 單機 TSP (Traveling Salesman Problem) 版本函數群
# ====================================================================
# 這些函數提供單機TSP的完整解決方案，用於與mTSP版本對比
# ====================================================================

def solve_greedy_tsp(CITIES, D, F, Fault_Mat, alpha, gamma, sigma, ENABLE_FAULT_SIGNAL, lambda_val):
    """
    貪婪演算法 (Nearest Neighbor) 解決單機 TSP
    邏輯：從起點 0 開始，每次貪心地選擇「綜合成本最低」的未拜訪城市，最後回到起點。
    特點：這種短視的貪心策略容易被「誘騙陷阱」欺騙，是 Greedy 的典型弱點。
    """
    start_time = time.time()
    unvisited = set(range(1, CITIES))  # 1 到 CITIES-1
    path = [0]                         # 起點為 0
    total_cost = 0.0

    # 綜合評估函數（與 QUBO 的權重公式一致）
    def edge_cost(u, v):
        if u == v: 
            return float('inf')
        w = D[u, v] + alpha * (sigma / (gamma**2)) * (F[u, v]**2)
        if ENABLE_FAULT_SIGNAL and Fault_Mat is not None:
            w += lambda_val * Fault_Mat[u, v]
        return w

    current_node = 0
    while unvisited:
        # 找下一個綜合成本最小的城市（局部最佳化）
        next_node = min(unvisited, key=lambda city: edge_cost(current_node, city))
        total_cost += edge_cost(current_node, next_node)
        path.append(next_node)
        unvisited.remove(next_node)
        current_node = next_node
    
    # 最後一哩路：被迫從最後一個點飛回起點（這就是 Greedy 的致命傷！）
    total_cost += edge_cost(current_node, 0)
    path.append(0)
    
    exec_time = time.time() - start_time
    return path, total_cost, exec_time


def build_robust_qubo_tsp(CITIES, D, F, Fault_Mat, alpha, gamma, sigma, ENABLE_FAULT_SIGNAL, lambda_val, A):
    """
    為單機 TSP 構建 Robust QUBO 模型
    
    變數索引化：x[v,i] 代表「城市 v 在第 i 個時間步被拜訪」
    約束：
      1. 每個城市恰好被拜訪一次
      2. 每個時間步恰好有一個城市被拜訪
    目標函數：最小化總路徑成本（包括距離與 H-infinity 風險）
    """
    Q = {}
    N_visit = CITIES - 1  # 需要拜訪的城市數量（不含起點 0）
    
    # 變數 index 映射：將 (city, step) 轉為 1D index
    def get_idx(city, step):
        # city 是 1~CITIES-1， step 是 0~N_visit-1
        return (city - 1) * N_visit + step

    # ===== 約束條件 A：每個城市只能在一個時間點被拜訪 (去且只去一次) =====
    for v in range(1, CITIES):
        for i in range(N_visit):
            idx = get_idx(v, i)
            Q[(idx, idx)] = Q.get((idx, idx), 0) - A
            for j in range(i + 1, N_visit):
                idx2 = get_idx(v, j)
                Q[(idx, idx2)] = Q.get((idx, idx2), 0) + 2 * A

    # ===== 約束條件 B：每個時間點只能排一個城市 =====
    for i in range(N_visit):
        for v in range(1, CITIES):
            idx = get_idx(v, i)
            Q[(idx, idx)] = Q.get((idx, idx), 0) - A
            for u in range(v + 1, CITIES):
                idx2 = get_idx(u, i)
                Q[(idx, idx2)] = Q.get((idx, idx2), 0) + 2 * A

    # ===== 綜合評估函數 =====
    def edge_cost(u, v):
        w = D[u, v] + alpha * (sigma / (gamma**2)) * (F[u, v]**2)
        if ENABLE_FAULT_SIGNAL and Fault_Mat is not None:
            w += lambda_val * Fault_Mat[u, v]
        return w

    # ===== 目標函數：計算路徑總成本 =====
    # (a) 從起點 0 出發的第一步
    for v in range(1, CITIES):
        idx = get_idx(v, 0)
        Q[(idx, idx)] = Q.get((idx, idx), 0) + edge_cost(0, v)
        
    # (b) 中間的連續路徑 (step i 到 step i+1)
    for i in range(N_visit - 1):
        for u in range(1, CITIES):
            for v in range(1, CITIES):
                if u != v:
                    idx_u = get_idx(u, i)
                    idx_v = get_idx(v, i + 1)
                    if idx_u < idx_v:
                        Q[(idx_u, idx_v)] = Q.get((idx_u, idx_v), 0) + edge_cost(u, v)
                    else:
                        Q[(idx_v, idx_u)] = Q.get((idx_v, idx_u), 0) + edge_cost(u, v)
                        
    # (c) 最後一步：從最後一個城市回到起點 0
    for u in range(1, CITIES):
        idx = get_idx(u, N_visit - 1)
        Q[(idx, idx)] = Q.get((idx, idx), 0) + edge_cost(u, 0)

    return Q


def decode_solution(sample, CITIES):
    """
    將 QUBO/SA/SQA 的量子採樣結果解碼成實際的城市訪問順序
    
    返回值：
      - list: [0, city1, city2, ..., cityN, 0] 的完整路徑
      - None: 如果結果違反 TSP 約束（重複、漏掉城市等）
    """
    N_visit = CITIES - 1
    path = [-1] * N_visit
    
    # 讀取量子/模擬退火的計算結果
    for key, val in sample.items():
        if val == 1:
            city = (key // N_visit) + 1
            step = key % N_visit
            if step < N_visit:
                path[step] = city
    
    # 防呆：如果演算法違反 TSP 規則（沒走完或重複走），標記為無效解
    if -1 in path or len(set(path)) != N_visit:
        return None 
        
    return [0] + path + [0]  # 加上起點與終點 0


def calculate_cost(path, D, F, Fault_Mat, alpha, gamma, sigma, ENABLE_FAULT_SIGNAL, lambda_val):
    """
    計算給定路徑的總成本
    
    成本公式：Σ (distance + H-infinity risk + fault penalty)
    
    對於無效路徑（None），返回無限大。
    """
    if path is None:
        return float('inf')  # 不合法路徑給予無限大成本
    
    total_cost = 0.0
    for i in range(len(path) - 1):
        u, v = path[i], path[i+1]
        w = D[u, v] + alpha * (sigma / (gamma**2)) * (F[u, v]**2)
        if ENABLE_FAULT_SIGNAL and Fault_Mat is not None:
            w += lambda_val * Fault_Mat[u, v]
        total_cost += w
    return total_cost


# ============================
# 複雜度與環境控制功能 (新增)
# ============================
def run_comparative_evaluations(Q, N, D, F, s_star):
    print(f"\n⚔️ [階段三] 正式對決 (SA vs SQA, {NUM_RUNS} runs)")
    
    results = {
        'sa': {'mk': [], 'risk': [], 'energy': [], 'faults': [], 'dist': [], 'best_mk': float('inf'), 'best_p1': None, 'best_p2': None},
        'sqa': {'mk': [], 'risk': [], 'energy': [], 'faults': [], 'dist': [], 'best_mk': float('inf'), 'best_p1': None, 'best_p2': None}
    }
    
    sqa_sampler = oj.SQASampler()
    sa_sampler = oj.SASampler()
    
    # 🌟 核心修改：建立 Sweeps 的退火排程
    TOTAL_SWEEPS = SWEEPS_MAIN_TEST
    
    smooth_sched = get_smooth_sqa_schedule(beta=BETA, total_sweeps=TOTAL_SWEEPS, num_steps=100)
    
    for r in range(NUM_RUNS):
        # SA 改為明確設定初始和結束溫度，讓 openjij 自動管理排程
        res_sa = sa_sampler.sample_qubo(
            Q,
            num_reads=NUM_READS,
            num_sweeps=SWEEPS_MAIN_TEST,
            beta_min=0.1,      # 高溫起點
            beta_max=10.0      # 低溫終點 (正規化後的合理值)
        )
        # 🌟 改成隨機取一條，讓結果有差異
        all_samples_sa = list(res_sa.samples())
        chosen_sa = random.choice(all_samples_sa)
        mk, risk, p1, p2 = get_makespan_and_risk(chosen_sa, N, D, F)
        fault_count = count_faults_hit(p1, Fault_Mat) + count_faults_hit(p2, Fault_Mat) if ENABLE_FAULT_SIGNAL else 0
        
        # 🌟 新增：計算raw distance
        raw_dist = 0
        for path in [p1, p2]:
            if not path: continue
            full = [0] + path + [0]
            for i in range(len(full)-1):
                raw_dist += D[full[i], full[i+1]]
        
        results['sa']['mk'].append(mk)
        results['sa']['risk'].append(risk)
        results['sa']['energy'].append(res_sa.first.energy)
        results['sa']['faults'].append(fault_count)
        results['sa']['dist'].append(raw_dist)
        if mk < results['sa']['best_mk']:
            results['sa']['best_mk'] = mk
            results['sa']['best_p1'] = p1[:]  # 🌟 使用 [:] 深度複製，確保快照
            results['sa']['best_p2'] = p2[:]  # 🌟 使用 [:] 深度複製，確保快照
        print(f"      [SA 執行 {r+1}/{NUM_RUNS}] 路徑 1: [0, {', '.join(map(str, p1))}, 0], 路徑 2: [0, {', '.join(map(str, p2))}, 0] | 成本: {mk:.2f} | 故障數: {fault_count}")
        
        # Standard SQA
        t0 = time.time()
        res_std = sqa_sampler.sample_qubo(Q, schedule=smooth_sched, num_reads=NUM_READS)
        # 註：smooth_sched 的 beta 已在 get_smooth_sqa_schedule 中根據 BETA 參數調整
        # 🌟 改成隨機取一條，讓結果有差異
        all_samples_sqa = list(res_std.samples())
        chosen_sqa = random.choice(all_samples_sqa)
        mk, risk, p1, p2 = get_makespan_and_risk(chosen_sqa, N, D, F)
        fault_count = count_faults_hit(p1, Fault_Mat) + count_faults_hit(p2, Fault_Mat) if ENABLE_FAULT_SIGNAL else 0
        
        # 🌟 新增：計算raw distance
        raw_dist = 0
        for path in [p1, p2]:
            if not path: continue
            full = [0] + path + [0]
            for i in range(len(full)-1):
                raw_dist += D[full[i], full[i+1]]
        
        results['sqa']['mk'].append(mk)
        results['sqa']['risk'].append(risk)
        results['sqa']['energy'].append(res_std.first.energy)
        results['sqa']['faults'].append(fault_count)
        results['sqa']['dist'].append(raw_dist)
        if mk < results['sqa']['best_mk']:
            results['sqa']['best_mk'] = mk
            results['sqa']['best_p1'] = p1[:]  # 🌟 使用 [:] 深度複製，確保快照
            results['sqa']['best_p2'] = p2[:]  # 🌟 使用 [:] 深度複製，確保快照
        print(f"      [Std-SQA 執行 {r+1}/{NUM_RUNS}] 耗時: {time.time()-t0:.2f}s | 路徑 1: [0, {', '.join(map(str, p1))}, 0], 路徑 2: [0, {', '.join(map(str, p2))}, 0] | 成本: {mk:.2f} | 故障數: {fault_count}")

    return results

def create_basic_distribution_charts(results, greedy_mk=None, greedy_risk=None, greedy_energy=None, greedy_faults=None, greedy_dist=None):
    # 🌟 修改：在 COMPARE_H_INFINITY 模式下忽略 greedy 数据（因为会导致显示过于复杂）
    if COMPARE_H_INFINITY:
        greedy_mk = greedy_risk = greedy_energy = greedy_faults = greedy_dist = None
        
    if COMPARE_H_INFINITY:
        # 根据是否有greedy数据调整列数 (3列: Risk, Faults, Distance)
        num_cols = 3
        fig, axes = plt.subplots(2, num_cols, figsize=(15, 10))
        title_suffix = " (+ Greedy)" if greedy_mk is not None else ""
        fig.suptitle(r"Performance Distribution: Without vs With $H_\infty$ Obstacle Avoidance" + title_suffix, fontsize=16, fontweight='bold', y=1.05)
        
        row_labels = ['Without Obstacle Avoidance', r'With $H_\infty$ Obstacle Avoidance']
        res_keys = ['baseline', 'robust']
        
        for row, r_key in enumerate(res_keys):
            res = results[r_key]
            
            if greedy_mk is not None:
                labels = ['Greedy', 'SA', 'SQA']
                colors = ['#2ca02c', 'gray', 'steelblue']
                risk_data = [greedy_risk if greedy_risk is not None else [0.0], res['sa']['risk'], res['sqa']['risk']]
                fault_data = [greedy_faults if greedy_faults is not None else [0], res['sa']['faults'], res['sqa']['faults']]
                dist_data = [greedy_dist if greedy_dist is not None else [0.0], res['sa']['dist'], res['sqa']['dist']]
            else:
                labels = ['SA', 'SQA']
                colors = ['gray', 'steelblue']
                risk_data = [res['sa']['risk'], res['sqa']['risk']]
                fault_data = [res['sa']['faults'], res['sqa']['faults']]
                dist_data = [res['sa']['dist'], res['sqa']['dist']]
            
            parts1 = axes[row, 0].violinplot(risk_data, showmeans=True)
            axes[row, 0].set_title(f"[{row_labels[row]}] Risk Distribution")
            axes[row, 0].set_ylabel("Risk Penalty")

            parts2 = axes[row, 1].violinplot(fault_data, showmeans=True)
            axes[row, 1].set_title(f"[{row_labels[row]}] Faults Hit Distribution")
            axes[row, 1].set_ylabel("Number of Faults")

            parts3 = axes[row, 2].violinplot(dist_data, showmeans=True)
            axes[row, 2].set_title(f"[{row_labels[row]}] Raw Distance")
            axes[row, 2].set_ylabel("Distance Units")
            
            parts_list = [parts1, parts2, parts3]
            
            for col, parts in enumerate(parts_list):
                ax = axes[row, col]
                ax.set_xticks(np.arange(1, len(labels) + 1))
                ax.set_xticklabels(labels, fontsize=11)
                ax.grid(alpha=0.3, axis='y')
                for pc, color in zip(parts['bodies'], colors):
                    pc.set_facecolor(color)
                    pc.set_alpha(0.7)

    else:
        # 单行模式：3个subplot (Risk, Faults, Raw Distance) - Makespan和Energy已单独分離
        fig, axes = plt.subplots(1, 3, figsize=(15, 5))
        title_suffix = " (+ Greedy)" if greedy_mk is not None else ""
        fig.suptitle("Performance Distribution Comparison" + title_suffix, fontsize=16, fontweight='bold', y=1.05)
        
        if greedy_mk is not None:
            labels = ['Greedy', 'SA', 'SQA']
            colors = ['#2ca02c', 'gray', 'steelblue']
            risk_data = [greedy_risk if greedy_risk is not None else [0.0], results['sa']['risk'], results['sqa']['risk']]
            fault_data = [greedy_faults if greedy_faults is not None else [0], results['sa']['faults'], results['sqa']['faults']]
            dist_data = [greedy_dist if greedy_dist is not None else [0.0], results['sa']['dist'], results['sqa']['dist']]
        else:
            labels = ['SA', 'SQA']
            colors = ['gray', 'steelblue']
            risk_data = [results['sa']['risk'], results['sqa']['risk']]
            fault_data = [results['sa']['faults'], results['sqa']['faults']]
            dist_data = [results['sa']['dist'], results['sqa']['dist']]
        
        # 1. Risk 分布
        parts1 = axes[0].violinplot(risk_data, showmeans=True)
        axes[0].set_title(r"$H_\infty$ Disturbance Risk Distribution")
        axes[0].set_ylabel("Risk Penalty")
        
        # 2. Faults 分布
        parts2 = axes[1].violinplot(fault_data, showmeans=True)
        axes[1].set_title("Faults Hit Distribution")
        axes[1].set_ylabel("Number of Faults")
        
        # 3. Raw Distance 分布
        parts3 = axes[2].violinplot(dist_data, showmeans=True)
        axes[2].set_title("Raw Total Distance (D)")
        axes[2].set_ylabel("Distance Units")
        
        # 設定提琴圖外觀
        parts_list = [parts1, parts2, parts3]
        for i, ax in enumerate(axes):
            ax.set_xticks(np.arange(1, len(labels) + 1))
            ax.set_xticklabels(labels, fontsize=11)
            ax.grid(alpha=0.3, axis='y')
            parts = parts_list[i]
            for pc, color in zip(parts['bodies'], colors):
                pc.set_facecolor(color)
                pc.set_alpha(0.7)

    plt.tight_layout(rect=[0, 0.03, 1, 0.95])
    plt.savefig(os.path.join(OUTPUT_DIR, "distribution_comparison.png"), dpi=300)
    print(f"\n🎨 圖表已繪製: {os.path.join(OUTPUT_DIR, 'distribution_comparison.png')}")
    # plt.show() (Moved to the end)

def create_makespan_energy_charts(results, greedy_mk=None, greedy_energy=None):
    """繪製第二分布圖：Makespan 與 Energy 分布對比 (2個子圖)
    注：Energy 分布不包含 Greedy 數據
    """
    if COMPARE_H_INFINITY:
        greedy_mk = greedy_energy = None
        
    if COMPARE_H_INFINITY:
        # 2x2 模式用於 H-infinity 對比
        fig, axes = plt.subplots(2, 2, figsize=(14, 10))
        fig.suptitle(r"Makespan & Energy Distribution: Without vs With $H_\infty$ Obstacle Avoidance", fontsize=16, fontweight='bold', y=1.05)
        
        row_labels = ['Without Obstacle Avoidance', r'With $H_\infty$ Obstacle Avoidance']
        res_keys = ['baseline', 'robust']
        
        for row, r_key in enumerate(res_keys):
            res = results[r_key]
            
            if greedy_mk is not None:
                mk_labels = ['Greedy', 'SA', 'SQA']
                mk_colors = ['#2ca02c', 'gray', 'steelblue']
                mk_data = [greedy_mk, res['sa']['mk'], res['sqa']['mk']]
                
                # Energy 不包含 Greedy
                energy_labels = ['SA', 'SQA']
                energy_colors = ['gray', 'steelblue']
                energy_data = [res['sa']['energy'], res['sqa']['energy']]
            else:
                mk_labels = ['SA', 'SQA']
                mk_colors = ['gray', 'steelblue']
                mk_data = [res['sa']['mk'], res['sqa']['mk']]
                
                energy_labels = ['SA', 'SQA']
                energy_colors = ['gray', 'steelblue']
                energy_data = [res['sa']['energy'], res['sqa']['energy']]
            
            # Makespan 圖
            parts1 = axes[row, 0].violinplot(mk_data, showmeans=True)
            axes[row, 0].set_title(f"[{row_labels[row]}] Makespan Distribution")
            axes[row, 0].set_ylabel("Distance Cost")
            axes[row, 0].set_xticks(np.arange(1, len(mk_labels) + 1))
            axes[row, 0].set_xticklabels(mk_labels, fontsize=11)
            axes[row, 0].grid(alpha=0.3, axis='y')
            for pc, color in zip(parts1['bodies'], mk_colors):
                pc.set_facecolor(color)
                pc.set_alpha(0.7)

            # Energy 圖 (不包含 Greedy)
            parts2 = axes[row, 1].violinplot(energy_data, showmeans=True)
            axes[row, 1].set_title(f"[{row_labels[row]}] Energy Distribution")
            axes[row, 1].set_ylabel("QUBO Energy")
            axes[row, 1].set_xticks(np.arange(1, len(energy_labels) + 1))
            axes[row, 1].set_xticklabels(energy_labels, fontsize=11)
            axes[row, 1].grid(alpha=0.3, axis='y')
            for pc, color in zip(parts2['bodies'], energy_colors):
                pc.set_facecolor(color)
                pc.set_alpha(0.7)

    else:
        # 單行模式：2個子圖 (Makespan, Energy)
        fig, axes = plt.subplots(1, 2, figsize=(14, 5))
        title_suffix = " (Makespan incl. Greedy, Energy excl. Greedy)" if greedy_mk is not None else ""
        fig.suptitle("Makespan & Energy Distribution Comparison" + title_suffix, fontsize=16, fontweight='bold', y=1.05)
        
        if greedy_mk is not None:
            # Makespan 包含 Greedy
            mk_labels = ['Greedy', 'SA', 'SQA']
            mk_colors = ['#2ca02c', 'gray', 'steelblue']
            mk_data = [greedy_mk, results['sa']['mk'], results['sqa']['mk']]
            
            # Energy 不包含 Greedy
            energy_labels = ['SA', 'SQA']
            energy_colors = ['gray', 'steelblue']
            energy_data = [results['sa']['energy'], results['sqa']['energy']]
        else:
            mk_labels = ['SA', 'SQA']
            mk_colors = ['gray', 'steelblue']
            mk_data = [results['sa']['mk'], results['sqa']['mk']]
            
            energy_labels = ['SA', 'SQA']
            energy_colors = ['gray', 'steelblue']
            energy_data = [results['sa']['energy'], results['sqa']['energy']]
        
        # 1. Makespan 分布 (包含 Greedy)
        parts1 = axes[0].violinplot(mk_data, showmeans=True)
        axes[0].set_title("Makespan (Distance Cost) Distribution")
        axes[0].set_ylabel("Distance Cost")
        axes[0].set_xticks(np.arange(1, len(mk_labels) + 1))
        axes[0].set_xticklabels(mk_labels, fontsize=11)
        axes[0].grid(alpha=0.3, axis='y')
        for pc, color in zip(parts1['bodies'], mk_colors):
            pc.set_facecolor(color)
            pc.set_alpha(0.7)
        
        # 2. Energy 分布 (不包含 Greedy)
        parts2 = axes[1].violinplot(energy_data, showmeans=True)
        axes[1].set_title("QUBO System Raw Energy")
        axes[1].set_ylabel("System Energy")
        axes[1].set_xticks(np.arange(1, len(energy_labels) + 1))
        axes[1].set_xticklabels(energy_labels, fontsize=11)
        axes[1].grid(alpha=0.3, axis='y')
        for pc, color in zip(parts2['bodies'], energy_colors):
            pc.set_facecolor(color)
            pc.set_alpha(0.7)

    plt.tight_layout(rect=[0, 0.03, 1, 0.95])
    plt.savefig(os.path.join(OUTPUT_DIR, "makespan_energy_comparison.png"), dpi=300)
    print(f"🎨 Makespan & Energy 分布圖已繪製: {os.path.join(OUTPUT_DIR, 'makespan_energy_comparison.png')}")
    # plt.show() (Moved to the end)


def plot_real_route_comparison(N, coords, p1_greedy, p2_greedy, p1_sa, p2_sa, p1_sqa, p2_sqa):
    """
    視覺化 Greedy 與 SA 與 SQA 的實體飛行路線對比 (使用真實座標)
    """
    fig, axes = plt.subplots(1, 3, figsize=(20, 7))
    fig.suptitle("UAV Actual Physical Trajectory: Greedy vs SA vs SQA", fontsize=16, fontweight='bold')
    
    G = nx.Graph()
    G.add_nodes_from(range(N))
    
    # 🌟 關鍵修改：使用真實的 (x, y) 座標作為節點位置 🌟
    pos = {i: (coords[i][0], coords[i][1]) for i in range(N)}
    
    titles = ["Greedy Algorithm", "Classical SA", "SQA"]
    routes_list = [(p1_greedy, p2_greedy), (p1_sa, p2_sa), (p1_sqa, p2_sqa)]
    
    for ax, title, (p1, p2) in zip(axes, titles, routes_list):
        ax.set_title(title, fontsize=14)
        
        # 畫節點 (依照真實地理位置散佈)
        nx.draw_networkx_nodes(G, pos, nodelist=[0], node_color='red', node_shape='s', node_size=300, ax=ax, label='Depot')
        nx.draw_networkx_nodes(G, pos, nodelist=range(1, N), node_color='skyblue', node_size=150, ax=ax)
        nx.draw_networkx_labels(G, pos, font_size=8, ax=ax)
        
        # 畫路線
        def add_edges(path, color, style):
            if not path: return
            edges = [(0, path[0])] + [(path[i], path[i+1]) for i in range(len(path)-1)] + [(path[-1], 0)]
            nx.draw_networkx_edges(G, pos, edgelist=edges, edge_color=color, style=style, width=2.0, arrows=True, ax=ax)

        add_edges(p1, 'blue', 'solid')
        add_edges(p2, 'darkorange', 'dashed')
        
    plt.tight_layout()
    plt.savefig(os.path.join(OUTPUT_DIR, "route_trajectory_comparison.png"), dpi=300)
    print(f"🎨 飛行軌跡對比圖已繪製: {os.path.join(OUTPUT_DIR, 'route_trajectory_comparison.png')}")

# ============================
# 複雜度與環境控制功能 (新增)
# =============================
def inject_hidden_shortcut(D, F):
    """
    注入隱藏的黃金捷徑 (Hidden Shortcut)：
    創造一條極度狹窄且成本極低的路線 (例如 4 -> 5 -> 6 -> 7)。
    如果沒照順序走，成本會極高；一旦走對，總成本會大幅下降。
    """
    N = D.shape[0]
    if N < 8:
        return D, F  # 城市太少無法做長捷徑

    # 定義黃金路線的節點 (假設是 4, 5, 6, 7)
    golden_nodes = [4, 5, 6, 7]

    # 1. 築起高牆：先把這幾個城市之間的所有連線，都變成極高成本 (阻止 SA 亂湊)
    for i in golden_nodes:
        for j in golden_nodes:
            if i != j:
                D[i, j] = 8.0
                F[i, j] = 0.8  # 高風險

    # 2. 挖出深谷：只開通 4->5, 5->6, 6->7 這條唯一且完美的捷徑
    for i in range(len(golden_nodes) - 1):
        n1 = golden_nodes[i]
        n2 = golden_nodes[i + 1]
        
        # 距離極短，風險為 0
        D[n1, n2] = 0.05
        D[n2, n1] = 0.05
        F[n1, n2] = 0.0
        F[n2, n1] = 0.0

    print("\n✨ [隱藏捷徑已佈署] 演算法將挑戰尋找極狹窄的黃金路線 (4->5->6->7)！")
    return D, F

def inject_deceptive_trap(D, F, alpha=10.0, gamma=0.5, sigma=1.0):
    """
    注入欺騙性陷阱 (Deceptive Trap)：
    改造矩陣，創造一條「看似完美的捷徑」，測試演算法是否會陷入局部陷阱。
    """
    N = D.shape[0]
    if N < 4:
        return D, F  # 城市太少無法做陷阱

    # 定義陷阱節點
    trap_start = 1
    trap_end = 2
    safe_detour = 3

    # ==========================================
    # 陷阱 1：致命捷徑 (距離極度誘人，但風險爆表)
    # ==========================================
    D[trap_start, trap_end] = 0.1
    D[trap_end, trap_start] = 0.1
    F[trap_start, trap_end] = 0.99
    F[trap_end, trap_start] = 0.99

    # ==========================================
    # 陷阱 2：安全繞路 (距離較遠，但完全無風險)
    # ==========================================
    # 路線: trap_start -> safe_detour -> trap_end
    D[trap_start, safe_detour] = 4.0
    D[safe_detour, trap_start] = 4.0
    D[safe_detour, trap_end] = 4.0
    D[trap_end, safe_detour] = 4.0

    F[trap_start, safe_detour] = 0.01
    F[safe_detour, trap_start] = 0.01
    F[safe_detour, trap_end] = 0.01
    F[trap_end, safe_detour] = 0.01

    
    # 計算並印出真實成本，讓您確認陷阱是否成立
    # 真實成本 = 距離 + alpha * (sigma / gamma^2) * F^2
    risk_multiplier = alpha * (sigma / (gamma**2))
    
    shortcut_cost = 0.1 + risk_multiplier * (0.99**2)
    detour_cost = (4.0 + risk_multiplier * (0.01**2)) + (4.0 + risk_multiplier * (0.01**2))
    
    print(f"   ➤ 致命捷徑 (1->2) 表觀距離: 0.1 | 隱藏真實成本: {shortcut_cost:.2f}")
    print(f"   ➤ 安全繞路 (1->3->2) 表觀距離: 8.0 | 隱藏真實成本: {detour_cost:.2f}")
    
    return D, F

def generate_controlled_matrices(num_cities, coord_std, risk_std, seed=42):
    """根據指定的標準差(變異數)生成 D 與 F 矩陣，精準控制問題複雜度"""
    np.random.seed(seed)
    
    # 1. 生成距離矩陣 D (使用常態分佈控制空間聚集度)
    coords = np.random.normal(loc=0.0, scale=coord_std, size=(num_cities, 2))
    D = np.zeros((num_cities, num_cities), dtype=float)
    for i in range(num_cities):
        for j in range(i + 1, num_cities):
            dist = np.linalg.norm(coords[i] - coords[j])
            D[i, j] = D[j, i] = dist
            
    # 2. 生成擾動矩陣 F (使用常態分佈控制風險極端值)
    F = np.abs(np.random.normal(loc=1.0, scale=risk_std, size=(num_cities, num_cities)))
    F = (F + F.T) / 2  # 確保對稱
    np.fill_diagonal(F, 0.0)
    
    return np.round(D, 2), np.round(F, 2)

def decode_and_eval(sample, N, D):
    # 快速解碼與貪婪/暴力求成本
    u1, u2 = [], []
    for p in range(N):
        for i in range(N):
            if sample.get(idx_mtsp(0, i, p, N), 0) == 1: u1.append(i)
            if sample.get(idx_mtsp(1, i, p, N), 0) == 1: u2.append(i)
    
    # 簡單修復 (去除重複，補齊缺失)
    c1, c2 = set(u1), set(u2)
    a1, a2 = [], []
    for city in range(1, N):
        if city in c1 and city not in c2: 
            a1.append(city)
        elif city in c2 and city not in c1: 
            a2.append(city)
        else:
            # 平票：計算該城市到兩個 UAV 分配城市群的平均距離
            if a1:
                avg_dist1 = np.mean([D[city, c] for c in a1])
            else:
                avg_dist1 = float('inf')
            
            if a2:
                avg_dist2 = np.mean([D[city, c] for c in a2])
            else:
                avg_dist2 = float('inf')
            
            # 選擇距離較近的 UAV
            if avg_dist1 <= avg_dist2: 
                a1.append(city)
            else: 
                a2.append(city)
            
    # 評估成本
    def eval_cost(cities):
        if not cities: return 0
        if len(cities) <= EXACT_LIMIT:
            bc = float('inf')
            for p in permutations(cities):
                c = D[0, p[0]] + sum(D[p[i], p[i+1]] for i in range(len(p)-1)) + D[p[-1], 0]
                bc = min(bc, c)
            return bc
        # 貪婪
        rem, path, curr = set(cities), [], 0
        while rem:
            nxt = min(rem, key=lambda x: D[curr, x])
            path.append(nxt)
            rem.remove(nxt)
            curr = nxt
        return D[0, path[0]] + sum(D[path[i], path[i+1]] for i in range(len(path)-1)) + D[path[-1], 0]

    return max(eval_cost(a1), eval_cost(a2))

def run_complexity_scaling():
    print(f"\n🚀 開始複雜度壓力測試...")
    print(f"⚙️ 空間變異數(COORD_STD)={COORD_STD}, 風險變異數(RISK_STD)={RISK_STD}")
    
    std_means, std_errs = [], []
    
    sampler = oj.SQASampler()
    
    for N in N_LIST:
        print(f"\n📊 測試規模 N={N} ...")
        D, F = generate_controlled_matrices(N, COORD_STD, RISK_STD)
        
        # 🌟 新增：生成故障矩陣
        Fault_Mat = generate_fault_matrix(N, prob=FAULT_PROBABILITY) if ENABLE_FAULT_SIGNAL else None
        
        # 動態計算懲罰值（基於路徑上界）
        max_d = np.max(D)
        max_f = np.max(F) if F is not None else 0
        coeff = ALPHA * SIGMA / (GAMMA ** 2)
        max_edge_weight = max_d + coeff * (max_f ** 2)  # 單條邊最大成本
        objective_upper_bound = CITIES * max_edge_weight  # 整條路徑上界
        
        # PENALTY 必須大於整條路徑的總成本上界
        dyn_penalty = objective_upper_bound * 2.0
        dyn_big_penalty = dyn_penalty * 3.0
        
        print(f'max_edge_weight = {max_edge_weight:.1f}')
        print(f'objective_upper_bound = {objective_upper_bound:.1f}')
        print(f'PENALTY = {dyn_penalty:.1f}')
        print(f'BIG_PENALTY = {dyn_big_penalty:.1f}')
        
        Q, _ = build_robust_qubo(D, F, Fault_Mat=Fault_Mat, penalty=dyn_penalty, big_penalty=dyn_big_penalty)
        # 🌟 QUBO 正規化
        Q = {k: v / dyn_penalty for k, v in Q.items()}
        
        std_results = []
        
        # 總步數統一為 SWEEPS_COMPLEXITY_TEST 步
        std_sched = get_smooth_sqa_schedule(BETA, total_sweeps=SWEEPS_COMPLEXITY_TEST, num_steps=20)
                
        # 這裡設定每次 N 跑的次數，與原本 NUM_RUNS 獨立，建議 10 次
        runs_for_scale = 10 
        for r in range(runs_for_scale):
            # Standard SQA
            res_std = sampler.sample_qubo(Q, schedule=std_sched, num_reads=NUM_READS_COMPLEXITY)
            std_results.append(decode_and_eval(res_std.first.sample, N, D))
            
        mean_std, err_std = np.mean(std_results), np.std(std_results)
        
        std_means.append(mean_std)
        std_errs.append(err_std)
        
        print(f"   SQA 平均: {mean_std:.2f} ± {err_std:.2f}")
        
    return std_means, std_errs

def plot_crossover(std_means, std_errs):
    plt.figure(figsize=(10, 6))
    plt.title(f"Complexity Scaling: Algorithm Performance vs Problem Size (N)\n(Coord Std={COORD_STD}, Risk Std={RISK_STD})", fontsize=14, fontweight='bold')
    
    # Draw line and error bars
    plt.errorbar(N_LIST, std_means, yerr=std_errs, fmt='-o', color='steelblue', 
                 linewidth=2.5, capsize=5, markersize=8, label='SQA')
    
    plt.xlabel("Problem Scale (Number of Cities $N$)", fontsize=12)
    plt.ylabel("Optimized Makespan (Distance Cost)", fontsize=12)
    plt.xticks(N_LIST)
    plt.grid(alpha=0.4, linestyle='--')
    plt.legend(fontsize=11)
    
    plt.tight_layout()
    plt.savefig(os.path.join(OUTPUT_DIR, "complexity_crossover.png"), dpi=300)
    print(f"🎨 Complexity scaling analysis saved to: {os.path.join(OUTPUT_DIR, 'complexity_crossover.png')}")
    # plt.show() (Moved to the end)

def count_faults_hit(path, Fault_Mat):
    """計算路徑中踩到多少個故障點 (故障路段數)"""
    if Fault_Mat is None or path is None:
        return 0
    fault_count = 0
    for i in range(len(path) - 1):
        city_from = path[i]
        city_to = path[i + 1]
        if Fault_Mat[city_from, city_to] > 0.5:  # 故障點
            fault_count += 1
    return fault_count

def plot_disturbance_and_fault_matrices(D, F, Fault_Mat):
    """將距離矩陣 (D)、擾動矩陣 (F) 與 故障矩陣 (Fault_Mat) 並排繪製熱力圖"""
    if Fault_Mat is None:
        return
    fig, axes = plt.subplots(1, 3, figsize=(18, 5))
    fig.suptitle("TSP Problem Setup: Distance, Disturbance, and Fault Map", fontsize=16, fontweight='bold')
    
    # 畫距離矩陣 D (連續值)
    sns.heatmap(D, ax=axes[0], cmap="Blues", annot=False)
    axes[0].set_title("Distance Matrix (D)", fontsize=14)
    axes[0].set_xlabel("City index"); axes[0].set_ylabel("City index")
    
    # 畫擾動矩陣 F (連續值)
    sns.heatmap(F, ax=axes[1], cmap="YlOrRd", annot=False)
    axes[1].set_title("Disturbance Matrix (F)", fontsize=14)
    axes[1].set_xlabel("City index"); axes[1].set_ylabel("City index")
    
    # 畫故障矩陣 Fault_Mat (0或1)
    sns.heatmap(Fault_Mat, ax=axes[2], cmap="Reds", cbar=False, linewidths=0.5, linecolor='lightgray')
    axes[2].set_title("Fault Matrix (0: Normal, 1: Fault)", fontsize=14)
    axes[2].set_xlabel("City index"); axes[2].set_ylabel("City index")
    
    plt.tight_layout()
    plt.savefig(os.path.join(OUTPUT_DIR, "disturbance_fault_matrices.png"), dpi=300)
    print(f"🎨 距離、擾動和故障矩陣圖已儲存: {os.path.join(OUTPUT_DIR, 'disturbance_fault_matrices.png')}")
    # plt.show() (Moved to the end)

def plot_comprehensive_best_comparison(greedy_mks, sa_mks, sqa_mks, greedy_p1, greedy_p2, sa_p1, sa_p2, sqa_p1, sqa_p2, D, F, Fault_Mat):
    """繪製 2x2 Subplot: 綜合比較最佳解的 Makespan、故障數、純距離、純擾動 (加入 Greedy)"""
    if Fault_Mat is None:
        return
        
    def calc_raw_metrics(p1, p2):
        dist_total, disturb_total, fault_hits = 0, 0, 0
        for path in [p1, p2]:
            if not path: continue
            full = [0] + path + [0]
            for i in range(len(full)-1):
                u, v = full[i], full[i+1]
                dist_total += D[u, v]
                disturb_total += F[u, v]
                fault_hits += Fault_Mat[u, v]
        return dist_total, disturb_total, fault_hits
        
    # 計算三者的真實物理指標
    greedy_dist, greedy_disturb, greedy_faults = calc_raw_metrics(greedy_p1, greedy_p2)
    sa_dist, sa_disturb, sa_faults = calc_raw_metrics(sa_p1, sa_p2)
    sqa_dist, sqa_disturb, sqa_faults = calc_raw_metrics(sqa_p1, sqa_p2)

    fig, axes = plt.subplots(2, 2, figsize=(15, 11))
    fig.suptitle("Comprehensive Performance: Greedy vs SA vs SQA", fontsize=18, fontweight='bold')
    
    labels = ['Greedy', 'Classical SA', 'SQA']
    colors = ['#2ca02c', '#4C72B0', '#DD8452']  # 綠色(Greedy), 藍色(SA), 橘色(SQA)
    
    def plot_bar(ax, vals, title, ylabel, is_int=False):
        bars = ax.bar(labels, vals, color=colors, edgecolor='black', linewidth=1.2)
        ax.set_title(title, fontsize=14, fontweight='bold')
        ax.set_ylabel(ylabel, fontsize=12)
        if is_int:
            ax.yaxis.set_major_locator(plt.MaxNLocator(integer=True))
        for bar, v in zip(bars, vals):
            yval = bar.get_height()
            text_str = f'{int(v)}' if is_int else f'{v:.1f}'
            ax.text(bar.get_x() + bar.get_width()/2.0, yval + (yval*0.01), text_str, 
                    ha='center', va='bottom', fontweight='bold', fontsize=12)

    plot_bar(axes[0, 0], [greedy_mks, sa_mks, sqa_mks], "1. Best Makespan (Total Objective Cost)", "Cost Score")
    plot_bar(axes[0, 1], [greedy_faults, sa_faults, sqa_faults], "2. Faults Hit (Safety Penalty)", "Number of Hits", is_int=True)
    axes[0, 1].set_title("2. Faults Hit (Safety Penalty)", color='darkred', fontweight='bold')
    plot_bar(axes[1, 0], [greedy_dist, sa_dist, sqa_dist], "3. Total Raw Distance (D)", "Distance Units")
    plot_bar(axes[1, 1], [greedy_disturb, sa_disturb, sqa_disturb], "4. Total Environmental Disturbance (F)", "Disturbance Level")

    plt.tight_layout(rect=[0, 0.03, 1, 0.95])
    plt.savefig(os.path.join(OUTPUT_DIR, "comprehensive_best_comparison.png"), dpi=300)
    print(f"🎨 三柱綜合最佳對比圖已儲存: {os.path.join(OUTPUT_DIR, 'comprehensive_best_comparison.png')}")
    # plt.show() (Moved to the end)

# ====================================================================
# 🌟 [新增] 單機 TSP 結果繪圖函數
# ====================================================================
def plot_tsp_results_comparison(tsp_results, city_coords=None):
    """
    繪製單機 TSP 的 Greedy vs SA 對比圖表
    
    包含：
    1. 成本對比柱狀圖
    2. 路徑可視化（如果提供城市座標）
    """
    if tsp_results is None:
        return
    
    # ===== 圖表 1: 成本對比 =====
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
    fig.suptitle("Single Machine TSP: Greedy vs Simulated Annealing", fontsize=16, fontweight='bold')
    
    # 成本對比
    algorithms = ['Greedy\n(Myopic)', 'SA\n(Metaheuristic)']
    costs = [
        tsp_results['comparison']['greedy_cost'],
        tsp_results['comparison']['sa_cost']
    ]
    colors = ['#2ca02c', '#4C72B0']
    bars = ax1.bar(algorithms, costs, color=colors, edgecolor='black', linewidth=1.5, width=0.5)
    
    # 添加數值標籤
    for bar, cost in zip(bars, costs):
        height = bar.get_height()
        ax1.text(bar.get_x() + bar.get_width()/2., height,
                f'{cost:.2f}',
                ha='center', va='bottom', fontweight='bold', fontsize=12)
    
    # 添加改進百分比
    improvement = tsp_results['comparison']['improvement']
    ax1.text(0.5, max(costs) * 0.5, 
            f"改進: {improvement:.1f}%",
            ha='center', fontsize=14, fontweight='bold',
            bbox=dict(boxstyle='round', facecolor='yellow', alpha=0.7))
    
    ax1.set_ylabel("Total Path Cost", fontsize=12)
    ax1.set_title("Cost Comparison", fontsize=13, fontweight='bold')
    ax1.grid(axis='y', alpha=0.3)
    
    # ===== 圖表 2: 路徑長度對比 =====
    greedy_path = tsp_results['greedy']['path']
    sa_path = tsp_results['sa']['path']
    
    ax2.text(0.5, 0.9, "Single Machine TSP Path Summary", 
            ha='center', fontsize=13, fontweight='bold', transform=ax2.transAxes)
    
    # Greedy 路徑信息
    if greedy_path is not None:
        ax2.text(0.05, 0.75, f"Greedy Path ({len(greedy_path)-2} cities):", 
                fontsize=11, fontweight='bold', transform=ax2.transAxes)
        greedy_str = f"[{greedy_path[0]}, {', '.join(map(str, greedy_path[1:-1]))}, {greedy_path[-1]}]"
        ax2.text(0.05, 0.65, greedy_str, fontsize=9, family='monospace', 
                transform=ax2.transAxes, wrap=True)
    else:
        ax2.text(0.05, 0.75, "Greedy Path: Invalid/None", 
                fontsize=11, fontweight='bold', color='red', transform=ax2.transAxes)
    
    ax2.text(0.05, 0.55, f"Cost: {tsp_results['comparison']['greedy_cost']:.2f}", 
            fontsize=10, color='#2ca02c', fontweight='bold', transform=ax2.transAxes)
    
    # SA 路徑信息
    if sa_path is not None:
        ax2.text(0.05, 0.40, f"SA Path ({len(sa_path)-2} cities):", 
                fontsize=11, fontweight='bold', transform=ax2.transAxes)
        sa_str = f"[{sa_path[0]}, {', '.join(map(str, sa_path[1:-1]))}, {sa_path[-1]}]"
        ax2.text(0.05, 0.30, sa_str, fontsize=9, family='monospace', 
                transform=ax2.transAxes, wrap=True)
    else:
        ax2.text(0.05, 0.40, "SA Path: Invalid/None", 
                fontsize=11, fontweight='bold', color='red', transform=ax2.transAxes)
    
    ax2.text(0.05, 0.20, f"Cost: {tsp_results['comparison']['sa_cost']:.2f}", 
            fontsize=10, color='#4C72B0', fontweight='bold', transform=ax2.transAxes)
    
    ax2.axis('off')
    
    plt.tight_layout()
    plt.savefig(os.path.join(OUTPUT_DIR, "tsp_comparison.png"), dpi=300, bbox_inches='tight')
    print(f"🎨 單機 TSP 對比圖已儲存: {os.path.join(OUTPUT_DIR, 'tsp_comparison.png')}")
    
    # ===== 圖表 3: 如果有城市座標，繪製路徑圖 =====
    if city_coords is not None and len(city_coords) > 0 and greedy_path is not None and sa_path is not None:
        fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
        fig.suptitle("Single Machine TSP: Path Visualization", fontsize=16, fontweight='bold')
        
        # 繪製 Greedy 路徑
        ax1.scatter(city_coords[:, 0], city_coords[:, 1], s=100, c='red', zorder=3, label='Cities')
        greedy_path_coords = city_coords[greedy_path]
        ax1.plot(greedy_path_coords[:, 0], greedy_path_coords[:, 1], 'g-', linewidth=2, label='Greedy Path')
        ax1.plot(greedy_path_coords[0, 0], greedy_path_coords[0, 1], 'go', markersize=12, label='Start/End')
        
        # 添加城市標籤
        for i, (x, y) in enumerate(city_coords):
            ax1.text(x, y+1, str(i), ha='center', fontsize=9, fontweight='bold')
        
        ax1.set_title(f"Greedy Path (Cost: {tsp_results['comparison']['greedy_cost']:.2f})", 
                     fontsize=12, fontweight='bold')
        ax1.set_xlabel("X Coordinate")
        ax1.set_ylabel("Y Coordinate")
        ax1.legend(loc='best')
        ax1.grid(alpha=0.3)
        ax1.set_aspect('equal')
        
        # 繪製 SA 路徑
        ax2.scatter(city_coords[:, 0], city_coords[:, 1], s=100, c='red', zorder=3, label='Cities')
        sa_path_coords = city_coords[sa_path]
        ax2.plot(sa_path_coords[:, 0], sa_path_coords[:, 1], 'b-', linewidth=2, label='SA Path')
        ax2.plot(sa_path_coords[0, 0], sa_path_coords[0, 1], 'bo', markersize=12, label='Start/End')
        
        # 添加城市標籤
        for i, (x, y) in enumerate(city_coords):
            ax2.text(x, y+1, str(i), ha='center', fontsize=9, fontweight='bold')
        
        ax2.set_title(f"SA Path (Cost: {tsp_results['comparison']['sa_cost']:.2f})", 
                     fontsize=12, fontweight='bold')
        ax2.set_xlabel("X Coordinate")
        ax2.set_ylabel("Y Coordinate")
        ax2.legend(loc='best')
        ax2.grid(alpha=0.3)
        ax2.set_aspect('equal')
        
        plt.tight_layout()
        plt.savefig(os.path.join(OUTPUT_DIR, "tsp_path_visualization.png"), dpi=300, bbox_inches='tight')
        print(f"🎨 單機 TSP 路徑可視化已儲存: {os.path.join(OUTPUT_DIR, 'tsp_path_visualization.png')}")
    elif sa_path is None or greedy_path is None:
        print("⚠️ 警告: SA 路徑無效(違反 TSP 約束)，略過路徑可視化圖表")

# ============================
# 參數記錄功能
# ============================
def save_parameters_to_file(output_dir):
    """將所有參數記錄到 .txt 文件"""
    params_file = os.path.join(output_dir, "parameters_log.txt")
    
    with open(params_file, 'w', encoding='utf-8') as f:
        f.write("=" * 80 + "\n")
        f.write("實驗參數設定記錄\n")
        f.write("=" * 80 + "\n\n")
        f.write(f"生成時間: {TIMESTAMP}\n\n")
        
        # 問題規模設定
        f.write("[問題規模設定]\n")
        f.write(f"CITIES = {CITIES}  # 初始核心問題規模\n")
        f.write(f"CITIES_BOUND = {CITIES_BOUND}  # 問題規模的上下範圍\n")
        f.write(f"RANDOM = {RANDOM}\n")
        f.write(f"RANDOM_SEED = {RANDOM_SEED}\n")
        f.write(f"COORD_RANGE = {COORD_RANGE}\n")
        f.write(f"N_LIST = {N_LIST}\n\n")
        
        # 環境變異數控制
        f.write("[環境變異數控制與複雜度測試參數]\n")
        f.write(f"COORD_STD = {COORD_STD}  # 空間分佈變異數\n")
        f.write(f"RISK_STD = {RISK_STD}  # 擾動風險變異數\n\n")
        
        # 演算法執行設定
        f.write("[演算法執行設定]\n")
        f.write(f"NUM_RUNS = {NUM_RUNS}  # 正式比較的執行次數\n")
        f.write(f"BETA = {BETA}  # 退火溫度參數\n")
        f.write(f"SWEEPS_MAIN_TEST = {SWEEPS_MAIN_TEST}  # 主要測試的退火步數\n")
        f.write(f"SWEEPS_COMPLEXITY_TEST = {SWEEPS_COMPLEXITY_TEST}  # 複雜度測試的退火步數\n")
        f.write(f"HEATMAP_RUNS = {HEATMAP_RUNS}  # 熱力圖測試每種情況的平均次數\n\n")
        
        # QUBO參數
        f.write("[QUBO參數]\n")
        f.write(f"PENALTY = {PENALTY}  # 約束違反懲罰權重\n")
        f.write(f"BIG_PENALTY = {BIG_PENALTY}  # 起點約束權重\n")
        f.write(f"EXACT_LIMIT = {EXACT_LIMIT}  # TSP精確求解上限\n\n")
        
        # 魯棒優化參數
        f.write("[魯棒優化參數 (H-infinity)]\n")
        f.write(f"USE_ROBUST = {USE_ROBUST}  # 是否啟用魯棒優化\n")
        f.write(f"GAMMA = {GAMMA}  # H-infinity 阻尼因子\n")
        f.write(f"SIGMA = {SIGMA}  # H-infinity 擾動標度\n")
        f.write(f"ALPHA = {ALPHA}  # H-infinity 權重係數\n\n")
        
        # 故障訊號參數
        f.write("[故障訊號參數]\n")
        f.write(f"ENABLE_FAULT_SIGNAL = {ENABLE_FAULT_SIGNAL}  # 是否加入故障訊號\n")
        f.write(f"FAULT_LAMBDA = {FAULT_LAMBDA}  # 故障懲罰權重\n")
        f.write(f"FAULT_PROBABILITY = {FAULT_PROBABILITY}  # 故障發生機率\n\n")
        
        # 執行區塊控制
        f.write("[執行區塊控制]\n")
        f.write(f"ENABLE_MULTI_UAV = {ENABLE_MULTI_UAV}  # 多機 mTSP 模式 (True=2台UAV, False=單機TSP)\n")
        f.write(f"RUN_MAIN_TEST = {RUN_MAIN_TEST}  # 執行主要演算法比較\n")
        f.write(f"RUN_TSP_TEST = {RUN_TSP_TEST}  # 執行獨立單機 TSP 測試\n")
        f.write(f"COMPARE_H_INFINITY = {COMPARE_H_INFINITY}  # 比較H-infinity效果\n")
        f.write(f"RUN_COMPLEXITY_TEST = {RUN_COMPLEXITY_TEST}  # 執行複雜度測試\n")
        f.write(f"SHOW_TERRAIN_PLOTS = {SHOW_TERRAIN_PLOTS}  # 顯示能量地形圖\n")
        f.write(f"ENABLE_DISTURBANCE_MATRIX = {ENABLE_DISTURBANCE_MATRIX}  # 生成擾動矩陣\n")
        f.write(f"ENABLE_QUANTUM_MINEFIELD = {ENABLE_QUANTUM_MINEFIELD}  # 產生量子雷區\n")
        f.write(f"PLOT_ENV_DIS_FAULT = {PLOT_ENV_DIS_FAULT}  # 繪製環境矩陣\n\n")
        
        f.write("=" * 80 + "\n")
        f.write("📝 參數記錄已保存\n")
        f.write("=" * 80 + "\n")
        
    print(f"📝 參數記錄文件已生成: {params_file}")
    return params_file

# ============================
# 主程式執行入口
# ===========================
if __name__ == "__main__":
    # 立即生成參數記錄文件
    save_parameters_to_file(OUTPUT_DIR)
    
    # 🌟 [新增] 顯示當前運行模式
    print("="*60)
    print("🎯 當前運行模式:")
    if ENABLE_MULTI_UAV:
        print("   ✓ 多機 mTSP 模式 (2台 UAV)")
    else:
        print("   ✓ 單機 TSP 模式 (1台 UAV)")
    print("="*60)
    
    print("="*50)
    if ENABLE_QUANTUM_MINEFIELD:
        print("<EFBFBD> 生成極限量子雷區測試...")
        D, F = generate_quantum_minefield(CITIES, random_seed=RANDOM_SEED)
        
        # 根據 ENABLE_DISTURBANCE_MATRIX 參數決定是否使用擾動矩陣
        if not ENABLE_DISTURBANCE_MATRIX:
            F = np.zeros((CITIES, CITIES))
            print("📍 擾動矩陣已禁用，使用零矩陣")
        
        # 為極限雷區建立障礙物矩陣
        Obstacle_Mat = np.zeros((CITIES, CITIES))
        
        # 為極限雷區建立環狀顯示用實體座標 (符合原有 layout 的精神)
        city_coords = []
        for i in range(CITIES):
            angle = 2 * np.pi * i / CITIES
            city_coords.append((50 + 40 * np.cos(angle), 50 + 40 * np.sin(angle)))
        
        # 進行數值驗證與畫圖
        if SHOW_TERRAIN_PLOTS:
            verify_minefield_stats(D, F)
            plot_energy_landscape_heatmap(D, F)
            plot_energy_landscape_3d(D, F)
    else:
        print("🚀 生成一般隨機任務環境...")
        # 1. 物理維度：隨機生成真實世界的城市座標，並計算物理距離 D
        city_coords = np.random.uniform(COORD_RANGE[0], COORD_RANGE[1], (CITIES, 2))
        D = calculate_distance_matrix(city_coords)
        
        # 2. 風險維度 + 障礙物維度：將蜘蛛網拓樸結構應用到擾動矩陣 F，並隨機放置障礙物
        if ENABLE_DISTURBANCE_MATRIX:
            # 🕷️ 同時生成蜘蛛網擾動矩陣 F 和障礙物矩陣 Obstacle_Mat
            F, Obstacle_Mat = generate_spider_web_with_obstacles(CITIES, num_obstacles=10)
        else:
            F = np.zeros((CITIES, CITIES))
            Obstacle_Mat = np.zeros((CITIES, CITIES))
            print("📍 擾動矩陣已禁用，使用零矩陣")
    
    # 3. 故障維度：使用障礙物矩陣 (與下游代碼相容性，Obstacle_Mat 映射為 Fault_Mat)
    Fault_Mat = Obstacle_Mat  # 向後兼容：障礙物矩陣作為故障矩陣
    # 可選：如果要啟用故障訊號，取消下面註釋
    # if ENABLE_FAULT_SIGNAL:
    #     print("🚨 故障訊號 (Fault Signal) 已啟用！正在生成隨機禁飛路段...")
    #     Fault_Mat = generate_fault_matrix(CITIES, prob=FAULT_PROBABILITY, seed=RANDOM_SEED)
    
    # ======== 🌟 改進：基於路徑成本上界的動態懲罰值 ========
    # 找出矩陣中最遙遠的距離與最大的擾動
    max_d = np.max(D)
    max_f = np.max(F) if F is not None else 0
    coeff = ALPHA * SIGMA / (GAMMA ** 2)
    max_edge_weight = max_d + coeff * (max_f ** 2)   # 單條邊最大成本（含風險項）
    objective_upper_bound = CITIES * max_edge_weight  # 整條路徑成本上界

    PENALTY = objective_upper_bound * 2.0             # 約束懲罰 > 整條路徑上界
    BIG_PENALTY = PENALTY * 3.0                       # 起點約束更強
    
    print(f"\n動態懲罰計算:")
    print(f"  max_edge_weight = {max_edge_weight:.1f}")
    print(f"  objective_upper_bound = {objective_upper_bound:.1f}")
    print(f"  PENALTY = {PENALTY:.1f}")
    print(f"  BIG_PENALTY = {BIG_PENALTY:.1f}\n")

    Q = None
    
    if RUN_MAIN_TEST:
        Q, _ = build_robust_qubo(D, F, Fault_Mat=Fault_Mat, penalty=PENALTY, big_penalty=BIG_PENALTY)
        # 🌟 QUBO 正規化：讓能量尺度回到 [0, 1] 附近
        Q = {k: v / PENALTY for k, v in Q.items()}
        
    if RUN_MAIN_TEST:
        print(f"🔧 [自動調校] 偵測到整條路徑成本上界: {objective_upper_bound:.2f}")
        print(f"🔧 [自動調校] 動態 PENALTY 設為: {PENALTY:.2f}, BIG_PENALTY 設為: {BIG_PENALTY:.2f}")
        
        if COMPARE_H_INFINITY:
            print("==================================================")
            print("🚀 [開始 H-infinity 避障算法對比測試]")
            print("==================================================")
            
            # --- 基準測試 (無避障策略，不考慮障礙物風險) ---
            print("\n--- 基準測試 (無避障策略 - Algorithm Ignores Obstacles) ---")
            Q_baseline, _ = build_robust_qubo(D, F, Fault_Mat=Fault_Mat, alpha=0.0, penalty=PENALTY, big_penalty=BIG_PENALTY)
            # 🌟 QUBO 正規化
            Q_baseline = {k: v / PENALTY for k, v in Q_baseline.items()}
            results_base = run_comparative_evaluations(Q_baseline, CITIES, D, F, 0.5)

            # --- 穩健測試 (有避障策略，考慮障礙物風險) ---
            print("\n--- 穩健測試 (有避障策略 - Algorithm Avoids Obstacles) ---")
            Q_robust, _ = build_robust_qubo(D, F, Fault_Mat=Fault_Mat, alpha=ALPHA, penalty=PENALTY, big_penalty=BIG_PENALTY)
            # 🌟 QUBO 正規化
            Q_robust = {k: v / PENALTY for k, v in Q_robust.items()}
            results_rob = run_comparative_evaluations(Q_robust, CITIES, D, F, 0.5)
            
            # Pack results for plotting
            results = {
                'baseline': results_base,
                'robust': results_rob
            }
            results_for_academic = results_rob

        else:
            print("==================================================")
            print(f"🚀 Robust mTSP Optimization (mTSP N={CITIES})")
            print("==================================================")
            
            print("\n[初始化] 距離矩陣 D:")
            print(D)
            print("\n[初始化] 擾動矩陣 F:")
            print(F)
            
            Q, _ = build_robust_qubo(D, F, Fault_Mat=Fault_Mat, penalty=PENALTY, big_penalty=BIG_PENALTY)
            # 🌟 QUBO 正規化：讓能量尺度回到 [0, 1] 附近
            Q = {k: v / PENALTY for k, v in Q.items()}
            
            # 正式評估
            results = run_comparative_evaluations(Q, CITIES, D, F, 0.5)
            results_for_academic = results
        
    # ====================================================================
    # 🌟 [新增] 單機 TSP 測試區塊
    # ====================================================================
    tsp_results = None
    if RUN_TSP_TEST:
        print("\n" + "="*60)
        print("🚀 [單機 TSP 測試] Traveling Salesman Problem (Single UAV)")
        print("="*60)
        
        # --- 方案 1: Greedy 近鄰演算法 ---
        print("\n[TSP 測試 1/3] 執行 Greedy 近鄰演算法...")
        tsp_greedy_path, tsp_greedy_cost, tsp_greedy_time = solve_greedy_tsp(
            CITIES, D, F, Fault_Mat=Fault_Mat if ENABLE_FAULT_SIGNAL else None,
            alpha=ALPHA, gamma=GAMMA, sigma=SIGMA, 
            ENABLE_FAULT_SIGNAL=ENABLE_FAULT_SIGNAL, lambda_val=FAULT_LAMBDA
        )
        print(f"  ✓ Greedy TSP 成本: {tsp_greedy_cost:.2f} | 耗時: {tsp_greedy_time:.4f}s")
        print(f"  路徑: {tsp_greedy_path}")
        
        # --- 方案 2: 構建 QUBO 模型用於 SA 和 SQA ---
        print("\n[TSP 測試 2/3] 構建單機 TSP 的 QUBO 模型...")
        Q_tsp = build_robust_qubo_tsp(
            CITIES, D, F, Fault_Mat=Fault_Mat if ENABLE_FAULT_SIGNAL else None,
            alpha=ALPHA, gamma=GAMMA, sigma=SIGMA, 
            ENABLE_FAULT_SIGNAL=ENABLE_FAULT_SIGNAL, lambda_val=FAULT_LAMBDA,
            A=PENALTY
        )
        print(f"  ✓ QUBO 矩陣規模: {len(Q_tsp)} 項")
        # 🌟 QUBO 正規化
        Q_tsp = {k: v / PENALTY for k, v in Q_tsp.items()}
        
        # --- 方案 3: 使用 SA 求解 TSP ---
        print("\n[TSP 測試 3/3] 使用 Simulated Annealing 求解 TSP...")
        sa_sampler = oj.SASampler()
        
        sa_tsp_results = []
        for run_idx in range(min(3, NUM_RUNS)):  # 只跑 3 次示意
            response_sa = sa_sampler.sample_qubo(
                Q_tsp,
                num_reads=NUM_READS_TSP,
                num_sweeps=SWEEPS_MAIN_TEST,
                beta_min=0.1,      # 高溫起點
                beta_max=10.0      # 低溫終點 (正規化後的合理值)
            )
            best_sample = response_sa.first.sample
            tsp_path = decode_solution(best_sample, CITIES)
            
            # 跳過無效解
            if tsp_path is None:
                print(f"  ⚠️ SA TSP Run {run_idx+1}: 無效解(違反 TSP 約束)")
                continue
            
            tsp_cost = calculate_cost(tsp_path, D, F, Fault_Mat=Fault_Mat if ENABLE_FAULT_SIGNAL else None,
                                     alpha=ALPHA, gamma=GAMMA, sigma=SIGMA, 
                                     ENABLE_FAULT_SIGNAL=ENABLE_FAULT_SIGNAL, lambda_val=FAULT_LAMBDA)
            sa_tsp_results.append({
                'path': tsp_path,
                'cost': tsp_cost,
                'energy': response_sa.first.energy
            })
            print(f"  ✓ SA TSP Run {run_idx+1}: 成本={tsp_cost:.2f}, 能量={response_sa.first.energy:.4f}")
        
        # 找最佳 SA 結果
        if len(sa_tsp_results) > 0:
            best_sa_tsp = min(sa_tsp_results, key=lambda x: x['cost'])
            print(f"\n  🏆 SA TSP 最佳結果: {best_sa_tsp['cost']:.2f}")
            print(f"  路徑: {best_sa_tsp['path']}")
        else:
            print("\n  ❌ SA TSP: 未找到有效解，使用 Greedy 作為備選")
            best_sa_tsp = {
                'path': tsp_greedy_path,
                'cost': tsp_greedy_cost,
                'energy': float('inf')
            }
        
        # --- 彙總結果 ---
        tsp_results = {
            'greedy': {
                'path': tsp_greedy_path,
                'cost': tsp_greedy_cost,
                'time': tsp_greedy_time
            },
            'sa': best_sa_tsp,
            'comparison': {
                'greedy_cost': tsp_greedy_cost,
                'sa_cost': best_sa_tsp['cost'],
                'improvement': ((tsp_greedy_cost - best_sa_tsp['cost']) / tsp_greedy_cost * 100) if tsp_greedy_cost > 0 else 0
            }
        }
        
        # 打印對比摘要
        print("\n" + "="*60)
        print("📊 單機 TSP 求解摘要:")
        print("="*60)
        print(f"Greedy 隨機演算法:  成本 = {tsp_results['comparison']['greedy_cost']:.2f}")
        print(f"Simulated Annealing: 成本 = {tsp_results['comparison']['sa_cost']:.2f}")
        print(f"優化改進: {tsp_results['comparison']['improvement']:.1f}%")
        print("="*60 + "\n")
        
    if RUN_5NODE_PERFECT_TRAP:
        print("\n🌟 執行 5 節點完美陷阱對比測試...")
        trap_results = test_5node_perfect_trap()
        print(f"✓ 測試完成，結果已記錄")
    
    s_m, s_e = None, None
    if RUN_COMPLEXITY_TEST:
        # Stage 4: Complexity stress testing and crossover analysis
        s_m, s_e = run_complexity_scaling()
        
    print("\n🏁 所有運算完成！開始繪製圖表...")
    
    # 🌟 新增：繪製 TSP 結果圖表
    if RUN_TSP_TEST and tsp_results is not None:
        print("\n📊 繪製單機 TSP 結果對比圖...")
        plot_tsp_results_comparison(tsp_results, city_coords=np.array(city_coords))
    
    if RUN_MAIN_TEST:
        # 🌟 新增：執行 Greedy 演算法 (提前計算以便傳給繪圖函數)
        t0_greedy = time.time()
        greedy_p1, greedy_p2, greedy_makespan = solve_greedy_mtsp(
            D, F, Fault_Mat=Fault_Mat if ENABLE_FAULT_SIGNAL else None
        )
        elapsed_time_greedy = time.time() - t0_greedy
        
        # 計算 Greedy 的故障數和raw distance
        greedy_fault_count = 0
        greedy_raw_dist = 0
        if ENABLE_FAULT_SIGNAL:
            greedy_fault_count = count_faults_hit(greedy_p1, Fault_Mat) + count_faults_hit(greedy_p2, Fault_Mat)
        
        # 🌟 新增：計算greedy的raw distance
        for path in [greedy_p1, greedy_p2]:
            if not path: continue
            full = [0] + path + [0]
            for i in range(len(full)-1):
                greedy_raw_dist += D[full[i], full[i+1]]
        
        # 輸出 Greedy 的結果
        print(f"      [Greedy 演算法] 耗時: {elapsed_time_greedy:.4f}s | 路徑 1: [0, {', '.join(map(str, greedy_p1))}, 0], 路徑 2: [0, {', '.join(map(str, greedy_p2))}, 0] | 成本: {greedy_makespan:.2f} | 故障數: {greedy_fault_count}")
        
        # 繪製圖表一：Risk、Faults、Distance 分布對比
        create_basic_distribution_charts(results, greedy_mk=[greedy_makespan], greedy_risk=[0.0], greedy_faults=[greedy_fault_count], greedy_dist=[greedy_raw_dist])
        
        # 繪製圖表：Makespan、Energy 分布對比
        create_makespan_energy_charts(results, greedy_mk=[greedy_makespan], greedy_energy=[greedy_makespan])
        
        # 繪製圖表：路徑對比圖
        if COMPARE_H_INFINITY:
            target_res = results['robust']
        else:
            target_res = results

        # 🌟 新增：1. 畫出擾動與故障矩陣對比熱力圖
        if ENABLE_FAULT_SIGNAL and PLOT_ENV_DIS_FAULT:
            plot_disturbance_and_fault_matrices(F, Fault_Mat)

        best_sa_1 = target_res['sa']['best_p1']
        best_sa_2 = target_res['sa']['best_p2']
        best_qa_1 = target_res['sqa']['best_p1']
        best_qa_2 = target_res['sqa']['best_p2']
        qa_makespan = target_res['sqa']['best_mk']
        
        sa_makespan = target_res['sa']['best_mk']
        
        # 🌟 修改：畫出「三柱」綜合指標對比圖 (Greedy vs SA vs SQA)
        if ENABLE_FAULT_SIGNAL:
            plot_comprehensive_best_comparison(
                greedy_mks=greedy_makespan,
                sa_mks=sa_makespan, 
                sqa_mks=qa_makespan, 
                greedy_p1=greedy_p1,
                greedy_p2=greedy_p2,
                sa_p1=best_sa_1, 
                sa_p2=best_sa_2, 
                sqa_p1=best_qa_1, 
                sqa_p2=best_qa_2, 
                D=D, 
                F=F, 
                Fault_Mat=Fault_Mat
            )
            
        plot_real_route_comparison(CITIES, city_coords, greedy_p1, greedy_p2, best_sa_1, best_sa_2, best_qa_1, best_qa_2)
        
    if RUN_COMPLEXITY_TEST and s_m is not None:
        # Plot Chart Three: Complexity Stress Testing and Crossover Analysis
        plot_crossover(s_m, s_e)
    
    print("\n" + "="*60)
    print("🏁 所有實驗與圖表繪製完成！")
    print("="*60)
    print(f"\n📁 所有圖表已儲存在: {os.path.abspath(OUTPUT_DIR)}")
    print("\n📊 已生成的圖表文件:")
    if os.path.exists(OUTPUT_DIR):
        files = sorted(os.listdir(OUTPUT_DIR))
        for i, fname in enumerate(files, 1):
            if fname.endswith('.png'):
                print(f"   {i}. {fname}")
    print("\n📈 圖表說明:")
    if RUN_TSP_TEST:
        print("   • tsp_comparison.png - 單機 TSP: Greedy vs SA 成本對比")
        print("   • tsp_path_visualization.png - 單機 TSP: 實際路徑可視化")
    if RUN_MAIN_TEST:
        print("   • distribution_comparison.png - mTSP: Risk/Faults/Distance 分布")
        print("   • makespan_energy_comparison.png - mTSP: Makespan & Energy 分布")
        print("   • comprehensive_best_comparison.png - mTSP: Greedy/SA/SQA 綜合對比")
        print("   • real_route_comparison.png - mTSP: 最優路徑可視化")
    if ENABLE_FAULT_SIGNAL and PLOT_ENV_DIS_FAULT:
        print("   • disturbance_fault_matrices.png - 擾動與故障矩陣熱力圖")
    print("\n💡 提示: 所有圖表窗口將同時顯示，請在各窗口完成檢視後關閉。")
    print("="*60 + "\n")
    
    # 同時顯示所有圖表窗口
    plt.show()