Openjij/single_uav.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

# ============ 輸出文件夾設置 (含時間戳) ============
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 = 50         # 初始核心問題規模 (測試建議先用 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 = 10      # 正式比較的執行次數
BETA = 1     # 🌟 稍微調高 (原為 10.0)，讓低溫結冰得更紮實
SWEEPS_MAIN_TEST = 1000      # 主要比較測試的退火步數
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 = 10.0       

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

# 執行區塊控制
RUN_MAIN_TEST = False            # 是否執行主要演算法(SA vs SQA)比較測試
RUN_TSP_TEST = True             # 🌟 是否執行單機 TSP (Traveling Salesman Problem) 測試
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 generate_quantum_minefield(N, random_seed=42):
    """
    生成具有「極大能量障礙」、「欺騙陷阱」與「峽谷路障」的極端測試矩陣。
    目標：逼迫 SA 掉入陷阱或卡死在路障前，而 SQA 能穿透勢壘找到繞道。
    """
    np.random.seed(random_seed)
    
    # 1. 築起能量高牆 (The Wall)
    D = np.random.uniform(30.0, 40.0, (N, N))
    F = np.random.uniform(20.0, 30.0, (N, N))
    np.fill_diagonal(D, 0)
    np.fill_diagonal(F, 0)
    
    # 2. 鑿出「黃金隧道 (Global Minimum)」
    golden_path = list(range(N))
    for i in range(N):
        u = golden_path[i]
        v = golden_path[(i + 1) % N]
        D[u, v] = 2.0   
        F[u, v] = 0.0   
        
    # 3. 佈置「致命陷阱 (Deceptive Traps)」
    for i in range(N):
        trap_node = (i + int(N/2)) % N  
        if trap_node != golden_path[(i + 1) % N] and trap_node != i:
            # 致命誘惑：眼前只要 0.1，但進去後就是 30.0 以上的高牆
            D[i, trap_node] = 0.1  
            F[i, trap_node] = 0.1  

    # ==========================================
    # 4. [新增] 在峽谷中砸下路障，並建立「穿隧繞道」
    # ==========================================
    # 挑選黃金隧道中的幾個節點作為路障發生地 (例如 1/3 和 2/3 處)
    obstacle_indices = [N // 3, 2 * N // 3] 
    
    for idx in obstacle_indices:
        u = golden_path[idx]
        v = golden_path[(idx + 1) % N]
        
        # [砸下巨石] 切斷原本順暢的黃金路線
        D[u, v] = 50.0
        F[u, v] = 50.0  # 模擬 I_fault 突發故障
        
        # [建立繞道] 徵用一個原本是陷阱的點作為安全繞道 (Bypass)
        bypass_node = (u + int(N/2) + 1) % N 
        
        # 第一段：高能勢壘 (阻擋 SA)
        # SA 在 u 點時看到：往陷阱只要 0.1，往繞道要 12.0。SA 必定會選陷阱然後死掉。
        D[u, bypass_node] = 12.0 
        F[u, bypass_node] = 0.0  
        
        # 第二段：柳暗花明 (專屬 SQA 的出路)
        # SQA 會看破全局，穿透 12.0 的障礙，走這一步接回黃金大道
        D[bypass_node, v] = 1.0  
        F[bypass_node, v] = 0.0  
        
        # 確保 Bypass 節點不要保留原本的陷阱死路設定
        # 關閉 Bypass 到其他隨機點的低成本誘惑
        for j in range(N):
            if j != v and j != bypass_node:
                D[bypass_node, j] = np.random.uniform(30.0, 40.0)
                F[bypass_node, j] = np.random.uniform(20.0, 30.0)

    print(f"⚠️ [極限測試] 已成功生成 N={N} 的量子雷區地貌（含動態路障）！")
    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 generate_disturbance_matrix(num_cities, seed=None):
    if seed is None: seed = RANDOM_SEED
    np.random.seed(seed)
    F = np.random.rand(num_cities, num_cities)
    F = (F + F.T) / 2
    np.fill_diagonal(F, 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 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

    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, 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)

    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)

    for i in range(1, N):
        vars_city = [idx_mtsp(k, i, p, N) for k in range(2) 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)

    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

# ---------------------------------------------------------
# 🌟 [新增] 貪婪演算法解決 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
    邏輯：每次挑選目前累積成本較低的 UAV，讓它飛往剩餘未拜訪城市中「當下綜合成本最低」的那個。
    Makespan 只計算距離（D），不包括風險和故障成本。
    """
    N = D.shape[0]
    unvisited = set(range(1, N))
    
    p1, p2 = [], []
    curr1, curr2 = 0, 0
    cost1, cost2 = 0.0, 0.0  # 綜合成本（用於決策）
    dist_cost1, dist_cost2 = 0.0, 0.0  # 純距離成本（用於 Makespan）
    
    # 計算兩點之間的綜合成本 (與 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

    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

# ---------------------------------------------------------
# 🌟 [新增] 產生平滑的 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):
    slots = []
    for k in range(2):
        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)
    return slots[0], slots[1]

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):
    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
        res_sa = sa_sampler.sample_qubo(Q, num_reads=10, num_sweeps=SWEEPS_MAIN_TEST)
        mk, risk, p1, p2 = get_makespan_and_risk(res_sa.first.sample, 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=10)
        mk, risk, p1, p2 = get_makespan_and_risk(res_std.first.sample, 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
        max_possible_cost = max_d + max_f
        dyn_penalty = max_possible_cost * 12.0
        dyn_big_penalty = dyn_penalty * 5.0
        
        Q, _ = build_robust_qubo(D, F, Fault_Mat=Fault_Mat, penalty=dyn_penalty, big_penalty=dyn_big_penalty)
        
        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=5)
            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(F, Fault_Mat):
    """將擾動矩陣 (F) 與 故障矩陣 (Fault_Mat) 並排繪製熱力圖"""
    if Fault_Mat is None:
        return
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
    fig.suptitle("Environment Disturbance vs. Fault Map", fontsize=16, fontweight='bold')
    
    # 畫擾動矩陣 F (連續值)
    sns.heatmap(F, ax=axes[0], cmap="YlOrRd", annot=False)
    axes[0].set_title("Disturbance Matrix (F)", fontsize=14)
    axes[0].set_xlabel("City index"); axes[0].set_ylabel("City index")
    
    # 畫故障矩陣 Fault_Mat (0或1)
    sns.heatmap(Fault_Mat, ax=axes[1], cmap="Reds", cbar=False, linewidths=0.5, linecolor='lightgray')
    axes[1].set_title("Fault Matrix (0: Normal, 1: Fault)", fontsize=14)
    axes[1].set_xlabel("City index"); axes[1].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"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("="*50)
    if ENABLE_QUANTUM_MINEFIELD:
        print("🚀 生成極限量子雷區測試...")
        D, F = generate_quantum_minefield(CITIES, random_seed=RANDOM_SEED)
        
        # 根據 ENABLE_DISTURBANCE_MATRIX 參數決定是否使用擾動矩陣
        if not ENABLE_DISTURBANCE_MATRIX:
            F = np.zeros((CITIES, CITIES))
            print("📍 擾動矩陣已禁用，使用零矩陣")
        
        # 為極限雷區建立環狀顯示用實體座標 (符合原有 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("🚀 生成一般隨機任務環境...")
        D, city_coords = generate_distance_matrix(CITIES, random=RANDOM)
        if ENABLE_DISTURBANCE_MATRIX:
            F = generate_disturbance_matrix(CITIES)
        else:
            F = np.zeros((CITIES, CITIES))
            print("📍 擾動矩陣已禁用，使用零矩陣")
        # 您可選擇是否在此呼叫原有的干擾產生器
        # D, F = inject_deceptive_trap(D, F, alpha=ALPHA, gamma=GAMMA, sigma=SIGMA)
        # D, F = inject_hidden_shortcut(D, F)
    
    # 🌟 新增：生成故障矩陣
    if ENABLE_FAULT_SIGNAL:
        print("🚨 故障訊號 (Fault Signal) 已啟用！正在生成隨機禁飛路段...")
        Fault_Mat = generate_fault_matrix(CITIES, prob=FAULT_PROBABILITY, seed=RANDOM_SEED)
    else:
        Fault_Mat = None
        
    # 原本的其他環境干擾暫時維持 (也可以選擇註解掉，因為 generate_quantum_minefield 已經很極端)
    # 為了保持之前的介面一致，這裡先保留動態懲罰值計算
    
    # ======== 🌟 新增：動態自適應懲罰值 ========
    # 找出矩陣中最遙遠的距離與最大的擾動
    max_d = np.max(D)
    max_f = np.max(F) if F is not None else 0
    max_possible_cost = max_d + max_f
    
    # 根據最大可能成本，動態設定懲罰值 (通常設為 10 倍 ~ 15 倍最穩定)
    PENALTY = max_possible_cost * 12.0
    BIG_PENALTY = PENALTY * 5.0

    Q = None
    
    if RUN_MAIN_TEST:
        Q, _ = build_robust_qubo(D, F, Fault_Mat=Fault_Mat, penalty=PENALTY, big_penalty=BIG_PENALTY)
        
    if RUN_MAIN_TEST:
        print(f"🔧 [自動調校] 偵測到單趟最大成本: {max_possible_cost:.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)
            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)
            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)
            
            # 正式評估
            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)} 項")
        
        # --- 方案 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=1, num_sweeps=SWEEPS_MAIN_TEST)
            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")
        
    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])
        
        # 繪製圖表：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()