Frequency Assignment Problem - QAOA¶

Standard QAOA with annealing schedule (without optimizer)¶

In [1]:
import warnings

import matplotlib.patches as patches
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.lines import Line2D
from qiskit import QuantumCircuit
from qiskit.primitives import StatevectorSampler
from qiskit.quantum_info import SparsePauliOp, Statevector

warnings.filterwarnings('ignore', category=DeprecationWarning)

np.random.seed(42)

N = 8                     # number of transmitters
positions = np.random.rand(N, 2) * 80   # 80×80 meters
freq_count = 3
num_qubits = N * freq_count   # qubits (one-hot encoding)
lambda_penalty = 30.0          # penalty for invalid assignment
interf_range = 30.0            # maximum interference range

P_LAYERS = 4  # Number of QAOA layers
shots = 8192

print(f"Number of transmitters: {N}")
print(f"Number of frequencies: {freq_count}")
print(f"Total qubits: {num_qubits}")
print(f"QAOA layers: {P_LAYERS}")
print("\nTransmitter positions:")
for i, pos in enumerate(positions):
    print(f"  Transmitter {i}: x={pos[0]:.2f}m, y={pos[1]:.2f}m")

Number of transmitters: 8
Number of frequencies: 3
Total qubits: 24
QAOA layers: 4

Transmitter positions:
  Transmitter 0: x=29.96m, y=76.06m
  Transmitter 1: x=58.56m, y=47.89m
  Transmitter 2: x=12.48m, y=12.48m
  Transmitter 3: x=4.65m, y=69.29m
  Transmitter 4: x=48.09m, y=56.65m
  Transmitter 5: x=1.65m, y=77.59m
  Transmitter 6: x=66.60m, y=16.99m
  Transmitter 7: x=14.55m, y=14.67m
In [2]:
def draw_transmitter(ax, x, y, color, label, size=8):
    """Draw transmitter as triangle with antenna"""
    ax.plot([x, x], [y, y + size*0.8], color='black', linewidth=3, zorder=5)
    triangle = patches.RegularPolygon((x, y), 3, radius=size*0.6, 
                                     orientation=0, 
                                     facecolor=color, 
                                     edgecolor='black', 
                                     linewidth=2,
                                     zorder=4)
    ax.add_patch(triangle)
    ax.text(x, y-size*0.3, label, fontsize=12, color='blue', 
            ha='center', va='center', weight='bold', zorder=6)

fig, ax = plt.subplots(figsize=(10, 10))
for i, (x, y) in enumerate(positions):
    draw_transmitter(ax, x, y, '#D3D3D3', str(i))
    circle = plt.Circle((x, y), 50, color='#B4E7F5', alpha=0.3, zorder=1)
    ax.add_patch(circle)

ax.set_title('Transmitter Positions (before frequency assignment)', fontsize=14, weight='bold')
ax.set_xlabel('Distance [m]', fontsize=12)
ax.set_ylabel('Distance [m]', fontsize=12)
ax.grid(True, alpha=0.3)
ax.set_aspect('equal')
ax.legend(['Transmitter', 'Interference zone (r=50m)'], loc='upper right')
plt.tight_layout()
plt.show()

distance_matrix = np.zeros((N, N))
for i in range(N):
    for j in range(N):
        distance_matrix[i, j] = np.linalg.norm(positions[i] - positions[j])

print("Distances between transmitters [m]:")
print("     ", end="")
for j in range(N):
    print(f"T{j:2d}   ", end="")
print()
for i in range(N):
    print(f"T{i:2d}: ", end="")
    for j in range(N):
        if i == j:
            print("  -   ", end="")
        else:
            print(f"{distance_matrix[i,j]:5.1f} ", end="")
    print()
No description has been provided for this image 
Distances between transmitters [m]:
     T 0   T 1   T 2   T 3   T 4   T 5   T 6   T 7   
T 0:   -    40.1  65.9  26.2  26.6  28.4  69.5  63.3 
T 1:  40.1   -    58.1  58.0  13.6  64.2  31.9  55.1 
T 2:  65.9  58.1   -    57.4  56.7  66.0  54.3   3.0 
T 3:  26.2  58.0  57.4   -    45.2   8.8  81.1  55.5 
T 4:  26.6  13.6  56.7  45.2   -    50.9  43.8  53.7 
T 5:  28.4  64.2  66.0   8.8  50.9   -    88.8  64.2 
T 6:  69.5  31.9  54.3  81.1  43.8  88.8   -    52.1 
T 7:  63.3  55.1   3.0  55.5  53.7  64.2  52.1   -

Pauli Hamiltonian Construction¶

One-hot Encoding and Binary Variables¶

Each transmitter $i$ is assigned $K$ qubits (where $K$ = number of frequencies). The qubit at position $q_{i,k} = i \cdot K + k$ represents:

  • $|1\rangle$: transmitter $i$ uses frequency $k$
  • $|0\rangle$: transmitter $i$ does not use frequency $k$

Mapping Binary Variables to Pauli Operators¶

To convert a classical optimization problem to a quantum Hamiltonian, we use the relationship between binary variables $x \in \{0, 1\}$ and the Pauli $Z$ operator:

$$x = \frac{1 - Z}{2}$$

where $Z|0\rangle = |0\rangle$ and $Z|1\rangle = -|1\rangle$. Therefore:

  • State $|0\rangle$: $Z = +1 \Rightarrow x = 0$
  • State $|1\rangle$: $Z = -1 \Rightarrow x = 1$

One-hot Constraint¶

For each transmitter, exactly one qubit must be in state $|1\rangle$. This constraint is formulated as a penalty term:

$$H_{\text{constraint}} = \lambda \sum_{i=0}^{N-1} \left( \sum_{k=0}^{K-1} x_{i,k} - 1 \right)^2$$

By substituting $x = \frac{1-Z}{2}$ and expanding the square, we obtain three types of terms:

Quadratic terms $Z_k Z_{k'}$ (cross-interactions between qubits of the same transmitter): $$\text{coefficient} = +2\lambda$$

Linear terms $Z_k$ (individual qubits): $$\text{coefficient} = -2\lambda$$

Constant term (identity, energy offset): $$\text{coefficient} = +4\lambda$$

For $K=3$ frequencies, each transmitter has: 3 cross terms + 3 linear terms + 1 offset = 7 Pauli terms.

Interference Terms¶

If two transmitters $i$ and $j$ are within interference range and use the same frequency $f$, we add a penalty:

$$H_{\text{interf}} = \sum_{(i,j): d_{ij} < r} \sum_f w_{ij} \cdot x_{i,f} \cdot x_{j,f}$$

where $w_{ij}$ is a weight dependent on distance (here using $20 \cdot \frac{r_{\text{max}}}{d_{ij}}$ dB). After substitution, this produces $Z_i Z_j$ terms with coefficient $\frac{w_{ij}}{4}$.

In [3]:
pauli_terms = []

print(f"Penalty coefficient λ = {lambda_penalty}")
print(f"Interference range = {interf_range}m")

# ONE-HOT CONSTRAINT: (Σ_k Z_k - 1)^2
print("Exactly one bit |1⟩ per transmitter")
print("Method: (Σ Z - 1)^2 with CORRECT signs")

constraint_count = 0
for i in range(N):
    base = i * freq_count

    # Quadratic terms Z_i * Z_j: +2*lambda
    for k1 in range(freq_count):
        for k2 in range(k1 + 1, freq_count):
            label = ['I'] * num_qubits
            label[base + k1] = 'Z'
            label[base + k2] = 'Z'
            pauli_string = ''.join(label)
            pauli_terms.append((pauli_string, lambda_penalty * 2.0))
            constraint_count += 1

    # Linear terms Z_i: -2*lambda
    for k in range(freq_count):
        label = ['I'] * num_qubits
        label[base + k] = 'Z'
        pauli_string = ''.join(label)
        pauli_terms.append((pauli_string, lambda_penalty * (-2.0)))
        constraint_count += 1

    # Constant term: +4*lambda
    label = ['I'] * num_qubits
    pauli_string = ''.join(label)
    pauli_terms.append((pauli_string, lambda_penalty * 4.0))
    constraint_count += 1

    if i == 0:
        print(
            f"  Transmitter {i}: 3 cross(+) + 3 linear(-) + 1 offset = 7 terms")

print(f"Total constraint terms: {constraint_count}")

# INTERFERENCE TERMS
interference_count = 0
interference_details = []

for i in range(N):
    for j in range(i + 1, N):
        dist = distance_matrix[i, j]
        if dist >= interf_range:
            continue

        loss_db = 20 * (interf_range / dist)

        for f in range(freq_count):
            p1 = i * freq_count + f
            p2 = j * freq_count + f
            label = ['I'] * num_qubits
            label[p1] = 'Z'
            label[p2] = 'Z'
            pauli_string = ''.join(label)
            coef = loss_db / 4.0
            pauli_terms.append((pauli_string, coef))
            interference_count += 1

            interference_details.append({
                'tx1': i,
                'tx2': j,
                'freq': f,
                'dist': dist,
                'loss': loss_db,
                'coef': coef
            })

print(f"Total interference terms: {interference_count}")
print("\nExamples:")
for idx, detail in enumerate(interference_details[:5]):
    print(
        f"  Transmitter {detail['tx1']} <-> {detail['tx2']}, "
        f"f={detail['freq']}: dist={detail['dist']:.1f}m, loss={detail['loss']:.1f}dB"
    )

Penalty coefficient λ = 30.0
Interference range = 30.0m
Exactly one bit |1⟩ per transmitter
Method: (Σ Z - 1)^2 with CORRECT signs
  Transmitter 0: 3 cross(+) + 3 linear(-) + 1 offset = 7 terms
Total constraint terms: 56
Total interference terms: 18

Examples:
  Transmitter 0 <-> 3, f=0: dist=26.2m, loss=22.9dB
  Transmitter 0 <-> 3, f=1: dist=26.2m, loss=22.9dB
  Transmitter 0 <-> 3, f=2: dist=26.2m, loss=22.9dB
  Transmitter 0 <-> 4, f=0: dist=26.6m, loss=22.6dB
  Transmitter 0 <-> 4, f=1: dist=26.6m, loss=22.6dB
In [4]:
print(f"Total: {constraint_count} terms")
print("\nFirst transmitter (T0) - all 7 terms:")
constraint_idx = 0
for i in range(1):  # Only first transmitter for detailed view
    base = i * freq_count
    print(f"\nTransmitter {i} constraint terms:")

    # Cross terms
    for k1 in range(freq_count):
        for k2 in range(k1 + 1, freq_count):
            pauli_word, coef = pauli_terms[constraint_idx]
            highlighted = [f"q{base+k1}:Z", f"q{base+k2}:Z"]
            print(
                f"  [{constraint_idx:2d}] {pauli_word} | coef={coef:6.1f} | Cross: {highlighted}"
            )
            constraint_idx += 1

    # Linear terms
    for k in range(freq_count):
        pauli_word, coef = pauli_terms[constraint_idx]
        highlighted = [f"q{base+k}:Z"]
        print(
            f"  [{constraint_idx:2d}] {pauli_word} | coef={coef:6.1f} | Linear: {highlighted}"
        )
        constraint_idx += 1

    # Constant term
    pauli_word, coef = pauli_terms[constraint_idx]
    print(
        f"  [{constraint_idx:2d}] {pauli_word} | coef={coef:6.1f} | Offset (Identity)"
    )
    constraint_idx += 1

# Display ALL interference terms
print(f"Total: {interference_count} terms")
print("\nAll interference terms (full listing):")
interf_start_idx = constraint_count
for idx, detail in enumerate(interference_details):
    pauli_idx = interf_start_idx + idx
    pauli_word, coef = pauli_terms[pauli_idx]
    q1 = detail['tx1'] * freq_count + detail['freq']
    q2 = detail['tx2'] * freq_count + detail['freq']
    print(
        f"  [{pauli_idx:2d}] T{detail['tx1']} <-> T{detail['tx2']} f{detail['freq']}: "
        f"q{q1}:Z x q{q2}:Z | "
        f"dist={detail['dist']:4.1f}m, loss={detail['loss']:5.1f}dB, coef={coef:5.2f}"
    )
    if idx < len(interference_details) - 1:
        # Show the actual Pauli word every 5 terms for verification
        if idx % 5 == 4:
            print(f"       Pauli word: {pauli_word}")

# Matrix visualization of Pauli operators
print("Legend: Z=Pauli-Z, I=Identity\n")

# Create a visual matrix for first few terms
print("Constraint terms (T0 example):")
print("     ", end="")
for q in range(num_qubits):
    print(f"q{q:2d}", end=" ")
print()

for idx in range(min(7, len(pauli_terms))):
    pauli_word, coef = pauli_terms[idx]
    print(f"[{idx:2d}] ", end="")
    for q, char in enumerate(pauli_word):
        if char == 'Z':
            print(f" Z ", end=" ")
        else:
            print(f" · ", end=" ")
    print(f"| coef={coef:6.1f}")

print("\nInterference terms (first 10):")
print("     ", end="")
for q in range(num_qubits):
    print(f"q{q:2d}", end=" ")
print()

for idx in range(constraint_count, min(constraint_count + 10,
                                       len(pauli_terms))):
    pauli_word, coef = pauli_terms[idx]
    print(f"[{idx:2d}] ", end="")
    for q, char in enumerate(pauli_word):
        if char == 'Z':
            print(f" Z ", end=" ")
        else:
            print(f" · ", end=" ")
    print(f"| coef={coef:6.1f}")

# Summary statistics
z_counts = []
for pauli_word, _ in pauli_terms:
    z_count = pauli_word.count('Z')
    z_counts.append(z_count)

print(f"Total Pauli terms: {len(pauli_terms)}")
print(f"Average Z operators per term: {np.mean(z_counts):.2f}")
print(f"Terms with 1 Z: {z_counts.count(1)} (linear)")
print(f"Terms with 2 Z: {z_counts.count(2)} (quadratic)")
print(f"Terms with 0 Z: {z_counts.count(0)} (constant)")

# Export all Pauli words to file
with open('pauli_words_complete.txt', 'w') as f:
    f.write("CONSTRAINT TERMS:\n")
    for idx in range(constraint_count):
        pauli_word, coef = pauli_terms[idx]
        f.write(f"[{idx:3d}] {pauli_word} | coef={coef:8.2f}\n")
        print(f"[{idx:3d}] {pauli_word} | coef={coef:8.2f}\n")

    f.write("\n\nINTERFERENCE TERMS:\n")
    for idx in range(constraint_count, len(pauli_terms)):
        pauli_word, coef = pauli_terms[idx]
        detail_idx = idx - constraint_count
        if detail_idx < len(interference_details):
            detail = interference_details[detail_idx]
            f.write(f"[{idx:3d}] {pauli_word} | coef={coef:8.2f} | "
                    f"T{detail['tx1']}<->T{detail['tx2']} f{detail['freq']}\n")
            print(f"[{idx:3d}] {pauli_word} | coef={coef:8.2f} | "
                  f"T{detail['tx1']}<->T{detail['tx2']} f{detail['freq']}\n")
        else:
            f.write(f"[{idx:3d}] {pauli_word} | coef={coef:8.2f}\n")

Total: 56 terms

First transmitter (T0) - all 7 terms:

Transmitter 0 constraint terms:
  [ 0] ZZIIIIIIIIIIIIIIIIIIIIII | coef=  60.0 | Cross: ['q0:Z', 'q1:Z']
  [ 1] ZIZIIIIIIIIIIIIIIIIIIIII | coef=  60.0 | Cross: ['q0:Z', 'q2:Z']
  [ 2] IZZIIIIIIIIIIIIIIIIIIIII | coef=  60.0 | Cross: ['q1:Z', 'q2:Z']
  [ 3] ZIIIIIIIIIIIIIIIIIIIIIII | coef= -60.0 | Linear: ['q0:Z']
  [ 4] IZIIIIIIIIIIIIIIIIIIIIII | coef= -60.0 | Linear: ['q1:Z']
  [ 5] IIZIIIIIIIIIIIIIIIIIIIII | coef= -60.0 | Linear: ['q2:Z']
  [ 6] IIIIIIIIIIIIIIIIIIIIIIII | coef= 120.0 | Offset (Identity)
Total: 18 terms

All interference terms (full listing):
  [56] T0 <-> T3 f0: q0:Z x q9:Z | dist=26.2m, loss= 22.9dB, coef= 5.72
  [57] T0 <-> T3 f1: q1:Z x q10:Z | dist=26.2m, loss= 22.9dB, coef= 5.72
  [58] T0 <-> T3 f2: q2:Z x q11:Z | dist=26.2m, loss= 22.9dB, coef= 5.72
  [59] T0 <-> T4 f0: q0:Z x q12:Z | dist=26.6m, loss= 22.6dB, coef= 5.65
  [60] T0 <-> T4 f1: q1:Z x q13:Z | dist=26.6m, loss= 22.6dB, coef= 5.65
       Pauli word: IZIIIIIIIIIIIZIIIIIIIIII
  [61] T0 <-> T4 f2: q2:Z x q14:Z | dist=26.6m, loss= 22.6dB, coef= 5.65
  [62] T0 <-> T5 f0: q0:Z x q15:Z | dist=28.4m, loss= 21.2dB, coef= 5.29
  [63] T0 <-> T5 f1: q1:Z x q16:Z | dist=28.4m, loss= 21.2dB, coef= 5.29
  [64] T0 <-> T5 f2: q2:Z x q17:Z | dist=28.4m, loss= 21.2dB, coef= 5.29
  [65] T1 <-> T4 f0: q3:Z x q12:Z | dist=13.6m, loss= 44.0dB, coef=10.99
       Pauli word: IIIZIIIIIIIIZIIIIIIIIIII
  [66] T1 <-> T4 f1: q4:Z x q13:Z | dist=13.6m, loss= 44.0dB, coef=10.99
  [67] T1 <-> T4 f2: q5:Z x q14:Z | dist=13.6m, loss= 44.0dB, coef=10.99
  [68] T2 <-> T7 f0: q6:Z x q21:Z | dist= 3.0m, loss=199.2dB, coef=49.81
  [69] T2 <-> T7 f1: q7:Z x q22:Z | dist= 3.0m, loss=199.2dB, coef=49.81
  [70] T2 <-> T7 f2: q8:Z x q23:Z | dist= 3.0m, loss=199.2dB, coef=49.81
       Pauli word: IIIIIIIIZIIIIIIIIIIIIIIZ
  [71] T3 <-> T5 f0: q9:Z x q15:Z | dist= 8.8m, loss= 68.0dB, coef=17.00
  [72] T3 <-> T5 f1: q10:Z x q16:Z | dist= 8.8m, loss= 68.0dB, coef=17.00
  [73] T3 <-> T5 f2: q11:Z x q17:Z | dist= 8.8m, loss= 68.0dB, coef=17.00
Legend: Z=Pauli-Z, I=Identity

Constraint terms (T0 example):
     q 0 q 1 q 2 q 3 q 4 q 5 q 6 q 7 q 8 q 9 q10 q11 q12 q13 q14 q15 q16 q17 q18 q19 q20 q21 q22 q23 
[ 0]  Z   Z   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·  | coef=  60.0
[ 1]  Z   ·   Z   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·  | coef=  60.0
[ 2]  ·   Z   Z   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·  | coef=  60.0
[ 3]  Z   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·  | coef= -60.0
[ 4]  ·   Z   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·  | coef= -60.0
[ 5]  ·   ·   Z   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·  | coef= -60.0
[ 6]  ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·  | coef= 120.0

Interference terms (first 10):
     q 0 q 1 q 2 q 3 q 4 q 5 q 6 q 7 q 8 q 9 q10 q11 q12 q13 q14 q15 q16 q17 q18 q19 q20 q21 q22 q23 
[56]  Z   ·   ·   ·   ·   ·   ·   ·   ·   Z   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·  | coef=   5.7
[57]  ·   Z   ·   ·   ·   ·   ·   ·   ·   ·   Z   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·  | coef=   5.7
[58]  ·   ·   Z   ·   ·   ·   ·   ·   ·   ·   ·   Z   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·  | coef=   5.7
[59]  Z   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   Z   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·  | coef=   5.6
[60]  ·   Z   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   Z   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·  | coef=   5.6
[61]  ·   ·   Z   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   Z   ·   ·   ·   ·   ·   ·   ·   ·   ·  | coef=   5.6
[62]  Z   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   Z   ·   ·   ·   ·   ·   ·   ·   ·  | coef=   5.3
[63]  ·   Z   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   Z   ·   ·   ·   ·   ·   ·   ·  | coef=   5.3
[64]  ·   ·   Z   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   Z   ·   ·   ·   ·   ·   ·  | coef=   5.3
[65]  ·   ·   ·   Z   ·   ·   ·   ·   ·   ·   ·   ·   Z   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·  | coef=  11.0
Total Pauli terms: 74
Average Z operators per term: 1.46
Terms with 1 Z: 24 (linear)
Terms with 2 Z: 42 (quadratic)
Terms with 0 Z: 8 (constant)
[  0] ZZIIIIIIIIIIIIIIIIIIIIII | coef=   60.00

[  1] ZIZIIIIIIIIIIIIIIIIIIIII | coef=   60.00

[  2] IZZIIIIIIIIIIIIIIIIIIIII | coef=   60.00

[  3] ZIIIIIIIIIIIIIIIIIIIIIII | coef=  -60.00

[  4] IZIIIIIIIIIIIIIIIIIIIIII | coef=  -60.00

[  5] IIZIIIIIIIIIIIIIIIIIIIII | coef=  -60.00

[  6] IIIIIIIIIIIIIIIIIIIIIIII | coef=  120.00

[  7] IIIZZIIIIIIIIIIIIIIIIIII | coef=   60.00

[  8] IIIZIZIIIIIIIIIIIIIIIIII | coef=   60.00

[  9] IIIIZZIIIIIIIIIIIIIIIIII | coef=   60.00

[ 10] IIIZIIIIIIIIIIIIIIIIIIII | coef=  -60.00

[ 11] IIIIZIIIIIIIIIIIIIIIIIII | coef=  -60.00

[ 12] IIIIIZIIIIIIIIIIIIIIIIII | coef=  -60.00

[ 13] IIIIIIIIIIIIIIIIIIIIIIII | coef=  120.00

[ 14] IIIIIIZZIIIIIIIIIIIIIIII | coef=   60.00

[ 15] IIIIIIZIZIIIIIIIIIIIIIII | coef=   60.00

[ 16] IIIIIIIZZIIIIIIIIIIIIIII | coef=   60.00

[ 17] IIIIIIZIIIIIIIIIIIIIIIII | coef=  -60.00

[ 18] IIIIIIIZIIIIIIIIIIIIIIII | coef=  -60.00

[ 19] IIIIIIIIZIIIIIIIIIIIIIII | coef=  -60.00

[ 20] IIIIIIIIIIIIIIIIIIIIIIII | coef=  120.00

[ 21] IIIIIIIIIZZIIIIIIIIIIIII | coef=   60.00

[ 22] IIIIIIIIIZIZIIIIIIIIIIII | coef=   60.00

[ 23] IIIIIIIIIIZZIIIIIIIIIIII | coef=   60.00

[ 24] IIIIIIIIIZIIIIIIIIIIIIII | coef=  -60.00

[ 25] IIIIIIIIIIZIIIIIIIIIIIII | coef=  -60.00

[ 26] IIIIIIIIIIIZIIIIIIIIIIII | coef=  -60.00

[ 27] IIIIIIIIIIIIIIIIIIIIIIII | coef=  120.00

[ 28] IIIIIIIIIIIIZZIIIIIIIIII | coef=   60.00

[ 29] IIIIIIIIIIIIZIZIIIIIIIII | coef=   60.00

[ 30] IIIIIIIIIIIIIZZIIIIIIIII | coef=   60.00

[ 31] IIIIIIIIIIIIZIIIIIIIIIII | coef=  -60.00

[ 32] IIIIIIIIIIIIIZIIIIIIIIII | coef=  -60.00

[ 33] IIIIIIIIIIIIIIZIIIIIIIII | coef=  -60.00

[ 34] IIIIIIIIIIIIIIIIIIIIIIII | coef=  120.00

[ 35] IIIIIIIIIIIIIIIZZIIIIIII | coef=   60.00

[ 36] IIIIIIIIIIIIIIIZIZIIIIII | coef=   60.00

[ 37] IIIIIIIIIIIIIIIIZZIIIIII | coef=   60.00

[ 38] IIIIIIIIIIIIIIIZIIIIIIII | coef=  -60.00

[ 39] IIIIIIIIIIIIIIIIZIIIIIII | coef=  -60.00

[ 40] IIIIIIIIIIIIIIIIIZIIIIII | coef=  -60.00

[ 41] IIIIIIIIIIIIIIIIIIIIIIII | coef=  120.00

[ 42] IIIIIIIIIIIIIIIIIIZZIIII | coef=   60.00

[ 43] IIIIIIIIIIIIIIIIIIZIZIII | coef=   60.00

[ 44] IIIIIIIIIIIIIIIIIIIZZIII | coef=   60.00

[ 45] IIIIIIIIIIIIIIIIIIZIIIII | coef=  -60.00

[ 46] IIIIIIIIIIIIIIIIIIIZIIII | coef=  -60.00

[ 47] IIIIIIIIIIIIIIIIIIIIZIII | coef=  -60.00

[ 48] IIIIIIIIIIIIIIIIIIIIIIII | coef=  120.00

[ 49] IIIIIIIIIIIIIIIIIIIIIZZI | coef=   60.00

[ 50] IIIIIIIIIIIIIIIIIIIIIZIZ | coef=   60.00

[ 51] IIIIIIIIIIIIIIIIIIIIIIZZ | coef=   60.00

[ 52] IIIIIIIIIIIIIIIIIIIIIZII | coef=  -60.00

[ 53] IIIIIIIIIIIIIIIIIIIIIIZI | coef=  -60.00

[ 54] IIIIIIIIIIIIIIIIIIIIIIIZ | coef=  -60.00

[ 55] IIIIIIIIIIIIIIIIIIIIIIII | coef=  120.00

[ 56] ZIIIIIIIIZIIIIIIIIIIIIII | coef=    5.72 | T0<->T3 f0

[ 57] IZIIIIIIIIZIIIIIIIIIIIII | coef=    5.72 | T0<->T3 f1

[ 58] IIZIIIIIIIIZIIIIIIIIIIII | coef=    5.72 | T0<->T3 f2

[ 59] ZIIIIIIIIIIIZIIIIIIIIIII | coef=    5.65 | T0<->T4 f0

[ 60] IZIIIIIIIIIIIZIIIIIIIIII | coef=    5.65 | T0<->T4 f1

[ 61] IIZIIIIIIIIIIIZIIIIIIIII | coef=    5.65 | T0<->T4 f2

[ 62] ZIIIIIIIIIIIIIIZIIIIIIII | coef=    5.29 | T0<->T5 f0

[ 63] IZIIIIIIIIIIIIIIZIIIIIII | coef=    5.29 | T0<->T5 f1

[ 64] IIZIIIIIIIIIIIIIIZIIIIII | coef=    5.29 | T0<->T5 f2

[ 65] IIIZIIIIIIIIZIIIIIIIIIII | coef=   10.99 | T1<->T4 f0

[ 66] IIIIZIIIIIIIIZIIIIIIIIII | coef=   10.99 | T1<->T4 f1

[ 67] IIIIIZIIIIIIIIZIIIIIIIII | coef=   10.99 | T1<->T4 f2

[ 68] IIIIIIZIIIIIIIIIIIIIIZII | coef=   49.81 | T2<->T7 f0

[ 69] IIIIIIIZIIIIIIIIIIIIIIZI | coef=   49.81 | T2<->T7 f1

[ 70] IIIIIIIIZIIIIIIIIIIIIIIZ | coef=   49.81 | T2<->T7 f2

[ 71] IIIIIIIIIZIIIIIZIIIIIIII | coef=   17.00 | T3<->T5 f0

[ 72] IIIIIIIIIIZIIIIIZIIIIIII | coef=   17.00 | T3<->T5 f1

[ 73] IIIIIIIIIIIZIIIIIZIIIIII | coef=   17.00 | T3<->T5 f2

Pauli Terms to Ising Model Conversion¶

Ising Hamiltonian Structure¶

A Pauli Hamiltonian containing only $Z$ operators can be directly mapped to the Ising model:

$$H_{\text{Ising}} = \sum_i h_i \sigma_i^z + \sum_{i<j} J_{ij} \sigma_i^z \sigma_j^z + C$$

where:

  • $h_i$ are single-qubit (local) fields
  • $J_{ij}$ are two-qubit (interaction) couplings
  • $C$ is the constant energy offset

The pauli_to_ising Algorithm¶

The function classifies each Pauli term by the number of $Z$ operators in the string:

# of $Z$ Term Type Example Storage
0 Identity IIIIII offset $\leftarrow$ offset + coef
1 Linear IIZIII $h[(i,)] \leftarrow h[(i,)]$ + coef
2 Quadratic IZIIZI $J[(i,j)] \leftarrow J[(i,j)]$ + coef

When the same position appears multiple times (e.g., multiple constraints generate the same $Z_i$ term), coefficients are accumulated.

Significance for QAOA Circuit¶

The decomposition into $h$ and $J$ enables efficient implementation of the cost Hamiltonian:

  • Linear terms $h_i$ → RZ(2γ·h_i) rotation on qubit $i$
  • Quadratic terms $J_{ij}$ → RZZ(2γ·J_ij) two-qubit gate on qubits $i, j$
In [5]:
def pauli_to_ising(pauli_terms, num_qubits):
    """
    Convert Pauli terms to Ising model format (h, J, offset).
    """
    h = {}  # Single-qubit terms
    J = {}  # Two-qubit terms
    offset = 0.0

    for pauli_string, coef in pauli_terms:
        z_positions = [i for i, char in enumerate(pauli_string) if char == 'Z']

        if len(z_positions) == 0:
            offset += coef
        elif len(z_positions) == 1:
            i = z_positions[0]
            key = (i,)
            h[key] = h.get(key, 0.0) + coef
        elif len(z_positions) == 2:
            i, j = z_positions
            key = (min(i, j), max(i, j))
            J[key] = J.get(key, 0.0) + coef

    return h, J, offset


# Convert Pauli terms to Ising format
h, J, offset = pauli_to_ising(pauli_terms, num_qubits)

print("Ising model conversion:")
print(f"  h terms (single-qubit): {len(h)}")
print(f"  J terms (two-qubit): {len(J)}")
print(f"  Offset: {offset:.2f}")

Ising model conversion:
  h terms (single-qubit): 24
  J terms (two-qubit): 42
  Offset: 960.00

QAOA Circuit¶

The qaoa_circuit function implements standard QAOA with the following structure:

  1. Initialization: Hadamard gates create a uniform superposition $|+\rangle^{\otimes n}$
  2. $P$ layers alternating:
    • Cost Hamiltonian $U_C(\gamma) = e^{-i\gamma H_C}$: RZ and RZZ gates from the Ising model
    • Mixer Hamiltonian $U_M(\beta) = e^{-i\beta H_M}$: RX rotations (X-mixer)

Parameters $\gamma$ and $\beta$ are normalized by the maximum coefficient $w_{\text{max}}$ for numerical stability.

Custom QAOA theory¶

is here at this site

Link: qiskit.cz QAOA commuting gates QUBO

In [6]:
def qaoa_circuit(gammas,
                 betas,
                 h,
                 J,
                 num_qubits,
                 mixer_type='x',
                 use_barriers=True):
    """
    Standard QAOA circuit with cost and mixer Hamiltonians.

    Parameters:
    - gammas: Parameters for cost Hamiltonian (list of length P)
    - betas: Parameters for mixer Hamiltonian (list of length P)
    - h: Single-qubit terms {(i,): coefficient}
    - J: Two-qubit terms {(i,j): coefficient}
    - num_qubits: Number of qubits
    - mixer_type: 'x' for standard X mixer, 'xy' for XY mixer
    """
    qc = QuantumCircuit(num_qubits)

    # Normalization
    h_vals = list(h.values()) if h else [0]
    J_vals = list(J.values()) if J else [0]
    wmax = max(np.max(np.abs(h_vals)), np.max(
        np.abs(J_vals))) if (h_vals or J_vals) else 1.0
    if wmax == 0:
        wmax = 1.0

    P = len(gammas)

    # Initial state: uniform superposition
    for i in range(num_qubits):
        qc.h(i)

    if use_barriers:
        qc.barrier()

    # P layers of QAOA
    for layer in range(P):
        gamma = gammas[layer]
        beta = betas[layer]

        # ---------- COST HAMILTONIAN ----------
        for (i,), v in h.items():
            if abs(v) > 1e-15:
                qc.rz(2 * gamma * v / wmax, i)

        for (i, j), vij in J.items():
            if abs(vij) > 1e-15:
                qc.rzz(2 * gamma * vij / wmax, i, j)

        if use_barriers:
            qc.barrier()

        # ---------- MIXER HAMILTONIAN ----------
        if mixer_type == 'x':
            for i in range(num_qubits):
                qc.rx(-2 * beta, i)
        elif mixer_type == 'xy':
            for i in range(num_qubits - 1):
                qc.rxx(-2 * beta, i, i + 1)
                qc.ryy(-2 * beta, i, i + 1)

        if use_barriers:
            qc.barrier()

    return qc


def energy_ising(bitstring, h, J, offset):
    """Calculate Ising energy for a given bitstring."""
    x = [int(b) for b in bitstring]
    s = [1 - 2 * xi for xi in x]  # Convert to spins

    energy = offset

    for (i,), v in h.items():
        energy += v * s[i]

    for (i, j), vij in J.items():
        energy += vij * s[i] * s[j]

    return energy


def samples_dict(counts, num_qubits):
    """Convert Qiskit counts to standardized samples dictionary."""
    samples = {}
    for bitstring, count in counts.items():
        if len(bitstring) < num_qubits:
            bitstring = bitstring.zfill(num_qubits)
        samples[bitstring] = count
    return samples

Annealing Schedule¶

Instead of optimizing parameters $\gamma, \beta$, we use a fixed annealing schedule:

  • $\beta$: decreases from 1 to 0 (mixer dominates at the beginning)
  • $\gamma$: increases from 0 to 1 (cost Hamiltonian dominates at the end)

This strategy mimics the adiabatic transition and often provides reasonable results without classical optimization.

In [7]:
# Fixed annealing schedule
betas = np.linspace(1, 0, P_LAYERS)   # Start high, decrease
gammas = np.linspace(0, 1, P_LAYERS)  # Start low, increase

print(f"Annealing schedule (P={P_LAYERS}):")
print(f"  Gammas: {gammas}")
print(f"  Betas:  {betas}")

# Build QAOA circuit
qc = qaoa_circuit(gammas, betas, h, J, num_qubits, mixer_type='x', use_barriers=True)
qc.measure_all()

print(f"  Depth: {qc.depth()}")
print(f"  Gates: {qc.size()}")

# Run circuit
sampler = StatevectorSampler()
job = sampler.run([qc], shots=shots)
result = job.result()
counts = result[0].data.meas.get_counts()

# Convert to samples
samples = samples_dict(counts, num_qubits)
print(f"Total unique bitstrings: {len(samples)}")

Annealing schedule (P=4):
  Gammas: [0.         0.33333333 0.66666667 1.        ]
  Betas:  [1.         0.66666667 0.33333333 0.        ]
  Depth: 38
  Gates: 408
Total unique bitstrings: 8131

Results Analysis and Decoding¶

One-hot Encoding Validation¶

Measured bitstrings are divided into blocks of $K$ qubits. A valid solution has exactly one 1 bit in each block.

Energy Calculation¶

For a bitstring $\mathbf{x}$, the Ising energy is computed as: $$E(\mathbf{x}) = \sum_i h_i s_i + \sum_{i<j} J_{ij} s_i s_j + C$$

where $s_i = 1 - 2x_i$ (mapping $\{0,1\} \to \{+1,-1\}$).

In [8]:
def is_valid_assignment(bitstring, n_transmitters, n_frequencies):
    """Check if bitstring represents valid one-hot encoding."""
    for i in range(n_transmitters):
        chunk = bitstring[i * n_frequencies:(i + 1) * n_frequencies]
        if chunk.count('1') != 1:
            return False
    return True


def decode_assignment(bitstring, n_transmitters, n_frequencies):
    """Decode bitstring to frequency assignment."""
    assignment = []
    for i in range(n_transmitters):
        chunk = bitstring[i * n_frequencies:(i + 1) * n_frequencies]
        freq = chunk.find('1')
        if freq == -1:
            freq = 0
        assignment.append(freq)
    return assignment


def count_collisions(assignment, distance_matrix, interf_range):
    """Count interference collisions."""
    collisions = 0
    collision_pairs = []
    n = len(assignment)
    for i in range(n):
        for j in range(i + 1, n):
            if assignment[i] == assignment[j] and distance_matrix[
                    i, j] < interf_range:
                collisions += 1
                collision_pairs.append((i, j))
    return collisions, collision_pairs


# Analyze all samples
print("Analyzing sampled solutions...")

valid_solutions = {}
invalid_count = 0

for bitstring, count in samples.items():
    if not is_valid_assignment(bitstring, N, freq_count):
        invalid_count += count
        continue

    energy = energy_ising(bitstring, h, J, offset)
    assignment = decode_assignment(bitstring, N, freq_count)
    collisions, collision_pairs = count_collisions(assignment, distance_matrix,
                                                   interf_range)

    valid_solutions[bitstring] = {
        'count': count,
        'energy': energy,
        'assignment': assignment,
        'collisions': collisions,
        'collision_pairs': collision_pairs
    }

valid_count = sum(s['count'] for s in valid_solutions.values())
print("\nResults:")
print(f"  Valid solutions: {len(valid_solutions)}")
print(f"  Valid shots: {valid_count}/{shots} ({valid_count/shots*100:.1f}%)")
print(
    f"  Invalid shots: {invalid_count}/{shots} ({invalid_count/shots*100:.1f}%)"
)

Analyzing sampled solutions...

Results:
  Valid solutions: 179
  Valid shots: 188/8192 (2.3%)
  Invalid shots: 8004/8192 (97.7%)
In [9]:
# Find best solution (minimum energy)
if len(valid_solutions) > 0:
    best_bitstring = min(valid_solutions.keys(),
                         key=lambda k: valid_solutions[k]['energy'])
    best_solution = valid_solutions[best_bitstring]

    print(f"\nBitstring: {best_bitstring}")
    print(
        f"Frequency: {best_solution['count']}/{shots} ({best_solution['count']/shots*100:.1f}%)"
    )
    print(f"Energy: {best_solution['energy']:.2f}")
    print(f"Collisions: {best_solution['collisions']}")

    assignment = best_solution['assignment']
    colors = ['#FFB3BA', '#BAFFC9', '#BAE1FF']
    freq_names = ['f0 (red)', 'f1 (green)', 'f2 (blue)']

    print("\nFrequency Assignment:")
    for i, f in enumerate(assignment):
        chunk = best_bitstring[i * freq_count:(i + 1) * freq_count]
        print(f"  Transmitter {i}: {chunk} → {freq_names[f]}")

    if best_solution['collision_pairs']:
        print("\nCollision details:")
        for i, j in best_solution['collision_pairs']:
            dist = distance_matrix[i, j]
            print(
                f"  T{i} <-> T{j}: {freq_names[assignment[i]]}, dist={dist:.1f}m"
            )
else:
    print("No valid solutions found!")
    assignment = [0] * N

Bitstring: 010100100001001100010001
Frequency: 1/8192 (0.0%)
Energy: -94.46
Collisions: 0

Frequency Assignment:
  Transmitter 0: 010 → f1 (green)
  Transmitter 1: 100 → f0 (red)
  Transmitter 2: 100 → f0 (red)
  Transmitter 3: 001 → f2 (blue)
  Transmitter 4: 001 → f2 (blue)
  Transmitter 5: 100 → f0 (red)
  Transmitter 6: 010 → f1 (green)
  Transmitter 7: 001 → f2 (blue)
In [10]:
# Show top solutions
print("\nTop 10 solutions by energy:")
print("-" * 90)

sorted_solutions = sorted(valid_solutions.items(),
                          key=lambda x: x[1]['energy'])[:10]

for i, (bitstring, sol) in enumerate(sorted_solutions, 1):
    print(f"{i:2d}. Energy: {sol['energy']:8.2f} | "
          f"Collisions: {sol['collisions']} | "
          f"Freq: {sol['count']:4d}/{shots} | "
          f"Assignment: {sol['assignment']}")

# Find zero-collision solutions
zero_collision = [
    (bs, sol) for bs, sol in valid_solutions.items() if sol['collisions'] == 0
]
if zero_collision:
    print(f"\nFound {len(zero_collision)} solution(s) with ZERO collisions!")
    best_zero = min(zero_collision, key=lambda x: x[1]['energy'])
    print(
        f"  Best zero-collision: Energy={best_zero[1]['energy']:.2f}, Assignment={best_zero[1]['assignment']}"
    )
else:
    print("\nNo zero-collision solutions found in samples.")

Top 10 solutions by energy:
------------------------------------------------------------------------------------------
 1. Energy:   -94.46 | Collisions: 0 | Freq:    1/8192 | Assignment: [1, 0, 0, 2, 2, 0, 1, 2]
 2. Energy:   -94.46 | Collisions: 0 | Freq:    1/8192 | Assignment: [0, 1, 2, 1, 2, 2, 2, 0]
 3. Energy:   -94.46 | Collisions: 0 | Freq:    1/8192 | Assignment: [2, 2, 0, 0, 1, 1, 0, 2]
 4. Energy:   -94.46 | Collisions: 0 | Freq:    1/8192 | Assignment: [2, 1, 2, 0, 0, 1, 0, 1]
 5. Energy:   -94.46 | Collisions: 0 | Freq:    1/8192 | Assignment: [1, 2, 1, 2, 0, 0, 0, 0]
 6. Energy:   -94.46 | Collisions: 0 | Freq:    1/8192 | Assignment: [1, 0, 1, 2, 2, 0, 0, 2]
 7. Energy:   -94.46 | Collisions: 0 | Freq:    1/8192 | Assignment: [1, 2, 2, 2, 0, 0, 1, 0]
 8. Energy:   -94.46 | Collisions: 0 | Freq:    1/8192 | Assignment: [2, 1, 1, 1, 0, 0, 1, 2]
 9. Energy:   -94.46 | Collisions: 0 | Freq:    1/8192 | Assignment: [0, 2, 2, 1, 1, 2, 2, 0]
10. Energy:   -94.46 | Collisions: 0 | Freq:    1/8192 | Assignment: [0, 1, 1, 2, 2, 1, 0, 0]

Found 13 solution(s) with ZERO collisions!
  Best zero-collision: Energy=-94.46, Assignment=[1, 0, 0, 2, 2, 0, 1, 2]
In [11]:
# Probability distribution from statevector
qc_statevector = qaoa_circuit(gammas, betas, h, J, num_qubits, mixer_type='x', use_barriers=False)
final_state = Statevector(qc_statevector)
probs = final_state.probabilities()

plt.figure(figsize=(14, 5))

plt.subplot(1, 2, 1)
top_50 = np.argsort(probs)[-50:][::-1]
plt.bar(range(len(top_50)), probs[top_50], color='#A8DADC', edgecolor='#457B9D', linewidth=1.5)
plt.xlabel('State rank', fontsize=12)
plt.ylabel('Probability', fontsize=12)
plt.title('Top 50 most probable states', fontsize=14, weight='bold')
plt.grid(True, alpha=0.3, axis='y')

plt.subplot(1, 2, 2)
plt.hist(probs[probs > 1e-6], bins=50, color='#F4A261', edgecolor='#E76F51', linewidth=1.5)
plt.xlabel('Probability', fontsize=12)
plt.ylabel('Number of states', fontsize=12)
plt.title('Probability distribution', fontsize=14, weight='bold')
plt.yscale('log')
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Top states from statevector
print("\nTop 10 most probable states (statevector):")
top_indices = np.argsort(probs)[-10:][::-1]
for rank, idx in enumerate(top_indices, 1):
    binary = format(idx, f'0{num_qubits}b')
    prob = probs[idx]
    is_valid = is_valid_assignment(binary, N, freq_count)
    if is_valid:
        energy = energy_ising(binary, h, J, offset)
        print(f"  {rank:2d}. |{binary}⟩  P={prob:.4f} ({prob*100:.2f}%)  E={energy:.1f}  ✓")
    else:
        print(f"  {rank:2d}. |{binary}⟩  P={prob:.4f} ({prob*100:.2f}%)  INVALID")
No description has been provided for this image 
Top 10 most probable states (statevector):
   1. |010000001001100001100010⟩  P=0.0000 (0.00%)  INVALID
   2. |100000001001010001010100⟩  P=0.0000 (0.00%)  INVALID
   3. |001000100100010010010001⟩  P=0.0000 (0.00%)  INVALID
   4. |001000100100010100010001⟩  P=0.0000 (0.00%)  INVALID
   5. |001000001001010010010100⟩  P=0.0000 (0.00%)  INVALID
   6. |100000001001100001100010⟩  P=0.0000 (0.00%)  INVALID
   7. |001000100100001100001010⟩  P=0.0000 (0.00%)  INVALID
   8. |010000001001100100100010⟩  P=0.0000 (0.00%)  INVALID
   9. |001000100100001010001010⟩  P=0.0000 (0.00%)  INVALID
  10. |100000001001010010010100⟩  P=0.0000 (0.00%)  INVALID
In [12]:
# Final visualization
colors = ['#FFB3BA', '#BAFFC9', '#BAE1FF']
freq_names = ['f0 (red)', 'f1 (green)', 'f2 (blue)']
freq_colors = [colors[f] for f in assignment]

fig, ax = plt.subplots(figsize=(12, 10))

# Interference lines
for i in range(N):
    for j in range(i + 1, N):
        dist = distance_matrix[i, j]
        if dist < interf_range and assignment[i] == assignment[j]:
            ax.plot([positions[i, 0], positions[j, 0]],
                    [positions[i, 1], positions[j, 1]],
                    color='#FF6B6B',
                    linewidth=3,
                    alpha=0.7,
                    zorder=2)
        elif dist < interf_range:
            ax.plot([positions[i, 0], positions[j, 0]],
                    [positions[i, 1], positions[j, 1]],
                    color='#95E1D3',
                    linewidth=1.5,
                    linestyle='--',
                    alpha=0.5,
                    zorder=1)

# Transmitters
for i, (x, y) in enumerate(positions):
    draw_transmitter(ax, x, y, freq_colors[i], str(i), size=6)
    circle = plt.Circle((x, y),
                        interf_range,
                        color=freq_colors[i],
                        alpha=0.3,
                        zorder=0)
    ax.add_patch(circle)

legend_elements = [
    Line2D([0], [0],
           marker='^',
           color='w',
           label=f'Frequency {i} ({freq_names[i].split()[0]})',
           markerfacecolor=colors[i],
           markersize=15,
           markeredgecolor='black',
           markeredgewidth=2) for i in range(freq_count)
]
legend_elements.extend([
    Line2D([0], [0],
           color='#FF6B6B',
           linewidth=3,
           label='Interference (collision)'),
    Line2D([0], [0],
           color='#95E1D3',
           linewidth=2,
           linestyle='--',
           label='In range (OK)')
])

ax.legend(handles=legend_elements, loc='upper left', fontsize=11)
ax.set_title('QAOA Result: Frequency Assignment', fontsize=16, weight='bold')
ax.set_xlabel('Distance [m]', fontsize=12)
ax.set_ylabel('Distance [m]', fontsize=12)
ax.grid(True, alpha=0.3)
ax.set_aspect('equal')
plt.tight_layout()
plt.show()
No description has been provided for this image