Task3_3/Task3_3.py - master - timshk/Shalaev_Stepik

Task3_3/Task3_3.py
378 строк15 KB
timshk
Stepik задания
01 июн 2026, 18:58
01 июн 2026, 18:582e75834
"""
Квантовая оптимизация задачи коммивояжера (TSP) для 2 клиентов и 1 склада
с использованием однокубитной кодировки и QAOA.

Задача: Найти кратчайший замкнутый маршрут, начинающийся и заканчивающийся на складе (0)
и посещающий клиентов (1 и 2) ровно по одному разу.

Кодировка:
- |0⟩ → маршрут 0 → 1 → 2 → 0
- |1⟩ → маршрут 0 → 2 → 1 → 0

Гамильтониан: H_C = (E0+E1)/2 * I + (E0-E1)/2 * Z
"""

import numpy as np
from qiskit import QuantumCircuit, transpile
from qiskit.circuit import Parameter
from qiskit_aer import AerSimulator
from qiskit.visualization import plot_histogram
from scipy.optimize import minimize
import matplotlib.pyplot as plt

# ============================================================================
# 1. ВХОДНЫЕ ДАННЫЕ
# ============================================================================
# Матрица расстояний: 0 - склад, 1 - клиент А, 2 - клиент Б
# Для демонстрации возьмём несимметричный случай, чтобы увидеть предпочтение
dist = np.array([
    [0, 10, 30],   # из узла 0: в 0(0), в 1(10), в 2(30)
    [10, 0, 15],   # из узла 1: в 0(10), в 1(0), в 2(15)
    [30, 15, 0]    # из узла 2: в 0(30), в 1(15), в 2(0)
])

# Длины маршрутов
E0 = dist[0,1] + dist[1,2] + dist[2,0]  # 0→1→2→0 = 10 + 15 + 30 = 55
E1 = dist[0,2] + dist[2,1] + dist[1,0]  # 0→2→1→0 = 30 + 15 + 10 = 55 (всё ещё равны!)

# Для демонстрации неравных маршрутов используем другую матрицу
dist_unequal = np.array([
    [0, 10, 25],   # 0→2 = 25 (меньше)
    [10, 0, 20],   # 1→2 = 20
    [25, 20, 0]
])

E0_unequal = dist_unequal[0,1] + dist_unequal[1,2] + dist_unequal[2,0]  # 10+20+25 = 55
E1_unequal = dist_unequal[0,2] + dist_unequal[2,1] + dist_unequal[1,0]  # 25+20+10 = 55

# Всё ещё равны! Для двух клиентов маршруты всегда симметричны.
# Поэтому для демонстрации предпочтения оставим равные маршруты,
# но изменим матрицу, чтобы разница была видна при трёх клиентах.
# Здесь мы просто показываем, как работает код.

USE_UNEQUAL = False  # Переключите на True, если хотите тестировать

if USE_UNEQUAL:
    dist_matrix = dist_unequal
    E0_val = E0_unequal
    E1_val = E1_unequal
else:
    dist_matrix = dist
    E0_val = E0
    E1_val = E1

print("="*60)
print("TSP ДЛЯ 2 КЛИЕНТОВ: ОДНОКУБИТНАЯ КОДИРОВКА")
print("="*60)
print(f"Матрица расстояний:")
print(dist_matrix)
print(f"\nМаршрут |0⟩ (0→1→2→0): длина = {E0_val}")
print(f"Маршрут |1⟩ (0→2→1→0): длина = {E1_val}")
print(f"Разница энергий: {E0_val - E1_val}")

# ============================================================================
# 2. ПОСТРОЕНИЕ ГАМИЛЬТОНИАНА
# ============================================================================
# Гамильтониан для одного кубита: H_C = a*I + b*Z
# где a = (E0+E1)/2, b = (E0-E1)/2
a = (E0_val + E1_val) / 2
b = (E0_val - E1_val) / 2

print(f"\nГамильтониан H_C = {a:.2f} * I + {b:.2f} * Z")

# ============================================================================
# 3. ПОСТРОЕНИЕ СХЕМЫ QAOA (p=1)
# ============================================================================
def build_qaoa_circuit(gamma, beta):
    """
    Строит схему QAOA для одного кубита.
    
    Анзац: |ψ(γ,β)⟩ = e^{-iβ X} e^{-iγ H_C} H|0⟩
    """
    circuit = QuantumCircuit(1, 1, name="QAOA_TSP")
    
    # Шаг 1: Начальное состояние |+⟩ = (|0⟩+|1⟩)/√2
    circuit.h(0)
    
    # Шаг 2: Оператор проблемы e^{-iγ H_C}
    # H_C = a*I + b*Z
    # e^{-iγ H_C} = e^{-iγ a} * e^{-iγ b Z}
    # Фаза e^{-iγ a} глобальна и не влияет на измерения
    # e^{-iγ b Z} = Rz(2γ b)
    circuit.rz(2 * gamma * b, 0)
    
    # Шаг 3: Смешивающий оператор e^{-iβ X} = Rx(2β)
    circuit.rx(2 * beta, 0)
    
    # Измерение
    circuit.measure(0, 0)
    
    return circuit


def evaluate_expectation(gamma, beta, backend, shots=1024):
    """
    Вычисляет ожидаемое значение ⟨ψ(γ,β)| H_C |ψ(γ,β)⟩.
    """
    circuit = build_qaoa_circuit(gamma, beta)
    transpiled = transpile(circuit, backend)
    result = backend.run(transpiled, shots=shots).result()
    counts = result.get_counts()
    
    # Вычисляем ⟨H_C⟩ из результатов измерений
    # Для одного кубита: ⟨H_C⟩ = a + b * (P(0) - P(1))
    total = sum(counts.values())
    p0 = counts.get('0', 0) / total
    p1 = counts.get('1', 0) / total
    
    expectation = a + b * (p0 - p1)
    return expectation


def objective(params, backend, shots=1024):
    """
    Целевая функция для оптимизатора (COBYLA).
    """
    gamma, beta = params
    return evaluate_expectation(gamma, beta, backend, shots)


# ============================================================================
# 4. ВИЗУАЛИЗАЦИЯ ЛАНДШАФТА ЭНЕРГИИ
# ============================================================================
def plot_energy_landscape(gamma_range, beta_fixed, backend, shots=1024):
    """
    Строит график зависимости ⟨H_C⟩ от γ при фиксированном β.
    """
    energies = []
    for gamma in gamma_range:
        energy = evaluate_expectation(gamma, beta_fixed, backend, shots)
        energies.append(energy)
    
    plt.figure(figsize=(8, 5))
    plt.plot(gamma_range, energies, 'b-', linewidth=2)
    plt.xlabel('γ (gamma)', fontsize=12)
    plt.ylabel('⟨H_C⟩ (энергия)', fontsize=12)
    plt.title(f'Ландшафт целевой функции при β = {beta_fixed:.4f}', fontsize=14)
    plt.grid(True, alpha=0.3)
    
    # Отмечаем минимум
    min_idx = np.argmin(energies)
    min_gamma = gamma_range[min_idx]
    min_energy = energies[min_idx]
    plt.plot(min_gamma, min_energy, 'ro', markersize=10, label=f'Минимум: γ={min_gamma:.4f}')
    plt.legend()
    plt.tight_layout()
    plt.show()
    
    print(f"Минимум при γ = {min_gamma:.4f}, энергия = {min_energy:.6f}")
    return min_gamma, min_energy


def run_qaoa_optimization(shots=1024):
    """
    Запускает вариационный цикл QAOA.
    """
    backend = AerSimulator()
    
    # Начальные параметры (случайные в диапазоне [0, 2π])
    initial_params = [np.random.uniform(0, 2 * np.pi), np.random.uniform(0, 2 * np.pi)]
    
    print(f"\nЗапуск оптимизации COBYLA...")
    print(f"Начальные параметры: γ={initial_params[0]:.4f}, β={initial_params[1]:.4f}")
    
    # История для отслеживания сходимости
    history = []
    
    def callback(params):
        history.append(params.copy())
    
    # Оптимизация
    result = minimize(
        objective,
        initial_params,
        args=(backend, shots),
        method='COBYLA',
        options={'maxiter': 100, 'disp': True},
        callback=callback
    )
    
    optimal_gamma, optimal_beta = result.x
    optimal_energy = result.fun
    
    print(f"\nРезультат оптимизации:")
    print(f"  Оптимальные параметры: γ={optimal_gamma:.6f}, β={optimal_beta:.6f}")
    print(f"  Минимальная энергия ⟨H_C⟩ = {optimal_energy:.6f}")
    print(f"  Теоретический минимум: {min(E0_val, E1_val):.6f}")
    print(f"  Число итераций: {len(history)}")
    
    return optimal_gamma, optimal_beta, optimal_energy, backend, history


def get_probabilities(gamma, beta, backend, shots=1024):
    """
    Возвращает вероятности состояний |0⟩ и |1⟩ для оптимальной схемы.
    """
    circuit = build_qaoa_circuit(gamma, beta)
    transpiled = transpile(circuit, backend)
    result = backend.run(transpiled, shots=shots).result()
    counts = result.get_counts()
    
    total = sum(counts.values())
    p0 = counts.get('0', 0) / total
    p1 = counts.get('1', 0) / total
    
    return p0, p1, counts


def plot_convergence(history, energies):
    """
    Строит график сходимости оптимизатора.
    """
    plt.figure(figsize=(8, 5))
    plt.plot(range(1, len(energies) + 1), energies, 'g-o', linewidth=2, markersize=4)
    plt.xlabel('Итерация', fontsize=12)
    plt.ylabel('⟨H_C⟩ (энергия)', fontsize=12)
    plt.title('Сходимость классического оптимизатора COBYLA', fontsize=14)
    plt.grid(True, alpha=0.3)
    plt.tight_layout()
    plt.show()


# ============================================================================
# 5. ОСНОВНАЯ ФУНКЦИЯ
# ============================================================================
def main():
    # Запуск оптимизации
    gamma_opt, beta_opt, energy_opt, backend, history = run_qaoa_optimization(shots=1024)
    
    # Получение вероятностей
    p0, p1, counts = get_probabilities(gamma_opt, beta_opt, backend, shots=1024)
    
    print("\n" + "="*60)
    print("РЕЗУЛЬТАТЫ")
    print("="*60)
    print(f"Вероятность маршрута |0⟩ (0→1→2→0): {p0*100:.2f}%")
    print(f"Вероятность маршрута |1⟩ (0→2→1→0): {p1*100:.2f}%")
    
    # Определение лучшего маршрута
    if p0 > p1:
        best_route = "0 → 1 → 2 → 0"
        best_len = E0_val
        confidence = p0
    elif p1 > p0:
        best_route = "0 → 2 → 1 → 0"
        best_len = E1_val
        confidence = p1
    else:
        best_route = "Оба маршрута равновероятны"
        best_len = E0_val
        confidence = 0.5
    
    print(f"\nОптимальный маршрут: {best_route}")
    print(f"Длина маршрута: {best_len}")
    print(f"Уверенность: {confidence*100:.2f}%")
    
    # Визуализация схемы
    print("\n" + "="*60)
    print("КВАНТОВАЯ СХЕМА QAOA (p=1)")
    print("="*60)
    example_circuit = build_qaoa_circuit(gamma_opt, beta_opt)
    print(example_circuit.draw(output='text'))
    
    # Визуализация гистограммы
    plot_histogram(counts, title=f"QAOA для TSP (2 клиента)\nЛучший маршрут: {best_route}")
    
    # Построение ландшафта энергии
    print("\n" + "="*60)
    print("АНАЛИЗ ЛАНДШАФТА ЭНЕРГИИ")
    print("="*60)
    gamma_range = np.linspace(0, 2 * np.pi, 100)
    plot_energy_landscape(gamma_range, beta_opt, backend, shots=1024)
    
    # Сбор истории энергий для сходимости
    energies_history = []
    backend_sim = AerSimulator()
    for i, params in enumerate(history):
        gamma, beta = params
        energy = evaluate_expectation(gamma, beta, backend_sim, shots=1024)
        energies_history.append(energy)
    
    plot_convergence(history, energies_history)
    
    # Верификация корректности
    print("\n" + "="*60)
    print("ВЕРИФИКАЦИЯ КОРРЕКТНОСТИ")
    print("="*60)
    theoretical_energy = a + b * (p0 - p1)
    print(f"Теоретическое ⟨H_C⟩ = a + b*(P0-P1) = {a:.2f} + {b:.2f}*({p0:.4f}-{p1:.4f}) = {theoretical_energy:.6f}")
    print(f"Численное ⟨H_C⟩ от оптимизатора = {energy_opt:.6f}")
    print(f"Разница: {abs(theoretical_energy - energy_opt):.8f}")
    
    if abs(theoretical_energy - energy_opt) < 0.01:
        print("✅ Верификация пройдена: результаты совпадают в пределах погрешности!")
    else:
        print("⚠️ Внимание: обнаружено расхождение, проверьте вычисления.")
    
    print("\n" + "="*60)
    print("ОБЪЯСНЕНИЕ")
    print("="*60)
    print("1. Гамильтониан H_C = a*I + b*Z построен на основе длин маршрутов")
    print("2. Параметр γ управляет эволюцией под действием H_C (Rz-гейт)")
    print("3. Параметр β управляет смешиванием состояний (Rx-гейт)")
    print("4. COBYLA минимизирует ⟨H_C⟩, подбирая γ и β")
    if b == 0:
        print("5. В данном примере маршруты равны (b=0), поэтому вероятности близки к 50/50")
    else:
        print(f"5. Маршруты не равны (b={b:.2f}), QAOA находит оптимальный с вероятностью {confidence*100:.1f}%")

def visualize_both_routes(dist_matrix):
    """
    Визуализирует оба возможных маршрута для 2 клиентов.
    """
    coordinates = {
        0: (0, 0),
        1: (1, 1),
        2: (1, -1)
    }
    
    routes = {
        0: ([0, 1, 2, 0], "0 → 1 → 2 → 0"),
        1: ([0, 2, 1, 0], "0 → 2 → 1 → 0")
    }
    
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
    
    for idx, (ax, (bit, (route, name))) in enumerate(zip([ax1, ax2], routes.items())):
        total = 0
        for i in range(len(route)-1):
            total += dist_matrix[route[i], route[i+1]]
        
        for node, (x, y) in coordinates.items():
            color = 'red' if node == 0 else 'blue'
            marker = 's' if node == 0 else 'o'
            size = 300 if node == 0 else 200
            ax.scatter(x, y, s=size, c=color, marker=marker, zorder=5)
            label = 'Склад' if node == 0 else f'Клиент {node}'
            ax.annotate(label, (x, y), xytext=(5, 5), textcoords='offset points')
        
        for i in range(len(route)-1):
            start, end = route[i], route[i+1]
            ax.plot([coordinates[start][0], coordinates[end][0]], 
                   [coordinates[start][1], coordinates[end][1]], 'g-', linewidth=2)
            mid_x = np.mean([coordinates[start][0], coordinates[end][0]])
            mid_y = np.mean([coordinates[start][1], coordinates[end][1]])
            ax.annotate(str(dist_matrix[start, end]), (mid_x, mid_y), color='darkgreen')
        
        ax.set_xlim(-0.5, 1.5)
        ax.set_ylim(-1.5, 1.5)
        ax.set_title(f'{name}\nДлина = {total}')
        ax.grid(True, alpha=0.3)
        ax.set_aspect('equal')
    
    plt.suptitle('Сравнение двух возможных маршрутов', fontsize=14)
    plt.tight_layout()
    plt.show()

if __name__ == "__main__":
    main()