main2.py - master - LadyChat2015/stepik

main2.py
180 строк7 KB
LadyChat2015
create: main2.py
05 апр 2026, 17:13
05 апр 2026, 17:13ce9e6a2
import numpy as np
from qiskit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp, Statevector
import matplotlib.pyplot as plt

# =====================================================
# Простая задача: 2 груза, 2 машины
# Груз 0: вес 10, Груз 1: вес 15
# Машина 0: вместимость 25, Машина 1: вместимость 20
# Оптимальное решение: оба груза в машину 0
# Переменные: x00, x01, x10, x11
# =====================================================

# Правильная QUBO матрица
# Цель: минимизировать количество машин + штрафы за нарушения
PENALTY = 100

Q = np.array([
    # x00   x01   x10   x11
    [ -1,    0,    0,    0],  # x00: стимул использовать машину 0
    [  0,   -1,    0,    0],  # x01: стимул использовать машину 1
    [  0,    0,   -1,    0],  # x10: стимул использовать машину 0
    [  0,    0,    0,   -1],  # x11: стимул использовать машину 1
])

# Добавляем штраф: каждый груз должен быть назначен ровно в одну машину
# Для груза 0: (x00 + x01 - 1)^2 = x00 + x01 - 2*x00*x01 + 1
Q[0,0] += PENALTY * (-2 + 1)  # -2 + 1 = -1
Q[1,1] += PENALTY * (-2 + 1)  # -2 + 1 = -1
Q[0,1] += PENALTY * 2
Q[1,0] += PENALTY * 2

# Для груза 1: (x10 + x11 - 1)^2
Q[2,2] += PENALTY * (-2 + 1)  # -2 + 1 = -1
Q[3,3] += PENALTY * (-2 + 1)  # -2 + 1 = -1
Q[2,3] += PENALTY * 2
Q[3,2] += PENALTY * 2

print("QUBO матрица:")
print(Q)

def qubo_to_ising(Q):
    n = Q.shape[0]
    terms = {"I"*n: 0}
    for i in range(n):
        for j in range(i, n):
            coeff = Q[i,j] + (Q[j,i] if i!=j else 0)
            if i == j:
                terms["I"*n] += coeff * 0.5
                z = ["I"]*n
                z[i] = "Z"
                terms["".join(z)] = terms.get("".join(z),0) - coeff * 0.5
            else:
                terms["I"*n] += coeff * 0.25
                zi = ["I"]*n
                zi[i] = "Z"
                terms["".join(zi)] = terms.get("".join(zi),0) - coeff * 0.25
                zj = ["I"]*n
                zj[j] = "Z"
                terms["".join(zj)] = terms.get("".join(zj),0) - coeff * 0.25
                zij = ["I"]*n
                zij[i] = "Z"
                zij[j] = "Z"
                terms["".join(zij)] = terms.get("".join(zij),0) + coeff * 0.25
    paulis, coeffs = [], []
    for p,c in terms.items():
        if abs(c) > 1e-6:
            paulis.append(p)
            coeffs.append(c)
    return SparsePauliOp(paulis, coeffs)

ham = qubo_to_ising(Q)
print(f"\nГамильтониан: {len(ham.paulis)} слагаемых")

def evaluate_energy(x, ham):
    energy = 0
    n = ham.num_qubits
    for pauli, coeff in zip(ham.paulis, ham.coeffs):
        sign = 1
        for i,ch in enumerate(pauli):
            if ch == 'Z' and x[i] == 1:
                sign = -sign
        energy += coeff * sign
    return energy

# Точное решение перебором
print("\n=== Точное решение (полный перебор) ===")
best_energy = float('inf')
best_x = None
solutions = []
for i in range(16):
    x = [int(b) for b in format(i, '04b')]
    e = evaluate_energy(x, ham)
    solutions.append((x, e))
    if e < best_energy:
        best_energy = e
        best_x = x
    print(f"{x} -> E={e:.4f}")

print(f"\nЛучшая битовая строка: {best_x}")
print(f"Энергия: {best_energy:.4f}")

# Интерпретация
x00, x01, x10, x11 = best_x
print(f"\nГруз 0 → машина {'0' if x00 else '1' if x01 else 'НЕТ'}")
print(f"Груз 1 → машина {'0' if x10 else '1' if x11 else 'НЕТ'}")

# QAOA с сеточным поиском
print("\n=== QAOA (сеточный поиск) ===")
best_gamma = 0
best_beta = 0
best_energy_qaoa = float('inf')

gamma_vals = np.linspace(0, np.pi, 20)
beta_vals = np.linspace(0, np.pi, 20)
energies_grid = np.zeros((len(gamma_vals), len(beta_vals)))

for i, gamma in enumerate(gamma_vals):
    for j, beta in enumerate(beta_vals):
        qc = QuantumCircuit(4)
        qc.h(range(4))
        for pauli, coeff in zip(ham.paulis, ham.coeffs):
            if pauli == "I"*4: continue
            z_pos = [idx for idx,ch in enumerate(pauli) if ch == 'Z']
            if len(z_pos) == 1:
                qc.rz(2*gamma*coeff, z_pos[0])
            elif len(z_pos) == 2:
                qc.cx(z_pos[0], z_pos[1])
                qc.rz(2*gamma*coeff, z_pos[1])
                qc.cx(z_pos[0], z_pos[1])
        for qubit in range(4):
            qc.rx(2*beta, qubit)
        sv = Statevector(qc)
        energy = np.real(sv.expectation_value(ham))
        energies_grid[i, j] = energy
        if energy < best_energy_qaoa:
            best_energy_qaoa = energy
            best_gamma = gamma
            best_beta = beta

print(f"Лучшие параметры: γ={best_gamma:.4f}, β={best_beta:.4f}")
print(f"Энергия QAOA: {best_energy_qaoa:.4f}")

# Финальное состояние
qc_final = QuantumCircuit(4)
qc_final.h(range(4))
for pauli, coeff in zip(ham.paulis, ham.coeffs):
    if pauli == "I"*4: continue
    z_pos = [idx for idx,ch in enumerate(pauli) if ch == 'Z']
    if len(z_pos) == 1:
        qc_final.rz(2*best_gamma*coeff, z_pos[0])
    elif len(z_pos) == 2:
        qc_final.cx(z_pos[0], z_pos[1])
        qc_final.rz(2*best_gamma*coeff, z_pos[1])
        qc_final.cx(z_pos[0], z_pos[1])
for qubit in range(4):
    qc_final.rx(2*best_beta, qubit)

sv = Statevector(qc_final)
probs = sv.probabilities_dict()
most_likely = max(probs, key=probs.get)
bits = most_likely[::-1]
print(f"\nНаиболее вероятное состояние QAOA: {bits}")
print(f"Вероятность: {probs[most_likely]:.4f}")

bits_list = [int(b) for b in bits]
x00, x01, x10, x11 = bits_list
print(f"\nГруз 0 → машина {'0' if x00 else '1' if x01 else 'НЕТ'}")
print(f"Груз 1 → машина {'0' if x10 else '1' if x11 else 'НЕТ'}")

# Визуализация
plt.figure(figsize=(10, 6))
plt.contourf(gamma_vals, beta_vals, energies_grid.T, levels=50, cmap='viridis')
plt.colorbar(label='Energy')
plt.plot(best_gamma, best_beta, 'r*', markersize=15, label='Best found')
plt.xlabel('γ')
plt.ylabel('β')
plt.title('QAOA Energy Landscape')
plt.legend()
plt.show()