Transpiler Tutorial for TSP with QUBO solver

# install jijmodeling-transpiler
!pip install jijmodeling-transpiler
!pip install openjij

import jijmodeling as jm
import openjij as oj
import jijmodeling_transpiler as jmt
import numpy as np
import matplotlib.pyplot as plt

d = jm.Placeholder("d", ndim=2)
n = d.len_at(0, latex="n")

i, j, t = jm.Element("i", n), jm.Element("j", n), jm.Element("t", n)

x = jm.BinaryVar("x", shape=(n, n))

problem = jm.Problem("TSP")
# Objective function
problem += jm.sum([i, j], d[i, j]*jm.sum(t, x[i, t]*x[j, (t+1) % n]))

# Constraints
jmC = jm.Constraint
problem += jmC("one-time", x[i, :].sum() == 1, forall=i)
problem += jmC("one-city", x[:, t].sum() == 1, forall=t)

problem

$\begin{array}{cccc}\text{Problem:} & \text{TSP} & & \\& & \min \quad \displaystyle \sum_{i = 0}^{n - 1} \sum_{j = 0}^{n - 1} d_{i, j} \cdot \sum_{t = 0}^{n - 1} x_{i, t} \cdot x_{j, \left(t + 1\right) \bmod n} & \\\text{{s.t.}} & & & \\ & \text{one-city} & \displaystyle \sum_{\ast_{0} = 0}^{n - 1} x_{\ast_{0}, t} = 1 & \forall t \in \left\{0,\ldots,n - 1\right\} \\ & \text{one-time} & \displaystyle \sum_{\ast_{1} = 0}^{n - 1} x_{i, \ast_{1}} = 1 & \forall i \in \left\{0,\ldots,n - 1\right\} \\\text{{where}} & & & \\& x & 2\text{-dim binary variable}\\\end{array}$

Transpile to QUBO

# Prepare data to create instance
def random_2d_tsp(n: int):
    x = np.random.uniform(0, 1, n)
    y = np.random.uniform(0, 1, n)
    XX, YY = np.meshgrid(x, y)
    distance = np.sqrt((XX - XX.T)**2 + (YY - YY.T)**2)
    return distance, (x, y)

num_city = 10
distance, (pos_x, pos_y) = random_2d_tsp(n=num_city)

plt.scatter(pos_x, pos_y)
plt.xlim(0, 1)
plt.ylim(0, 1)
plt.show()

# Substitute data to the mathematical model
compiled_instnace = jmt.core.compile_model(problem, {"d": distance})

# Transpile to QUBO
qubo_builder = jmt.core.pubo.transpile_to_pubo(compiled_instnace)

Why PUBO?

The transpile_to_pubo function is designed to transform a given mathematical model into a Polynomial Unconstrained Binary Optimization (PUBO) formulation. The term "PUBO" stands for Polynomial Unconstrained Binary Optimization, a more generalized form of Quadratic Unconstrained Binary Optimization (QUBO).

multipliers = {"one-time": 1.0, "one-city": 1.0}
qubo, offset = qubo_builder.get_qubo_dict(multipliers=multipliers)

Optimize by QUBO SA Sampler

sampler = oj.SASampler()
response = sampler.sample_qubo(qubo)

Decode the response from openjij

decoded = jmt.core.pubo.decode_from_openjij(response, qubo_builder, compiled_instnace)

feasibles = decoded.feasible()
objectives = np.array(feasibles.evaluation.objective)
lowest_index = np.argmin(objectives)
print(f"Lowest solution index: {lowest_index}, Lowest objective value: {objectives[lowest_index]}")

Lowest solution index: 0, Lowest objective value: 4.261925958221673

# check solution
nonzero_indices, nonzero_values, shape = feasibles.record.solution["x"][lowest_index]
print("indices: ", nonzero_indices)
print("values: ", nonzero_values)

indices:  ([0, 1, 2, 3, 4, 5, 6, 7, 8, 9], [4, 0, 1, 7, 5, 3, 6, 2, 8, 9])
values:  [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]

Draw tour

def convert_to_tour(nonzero_indices, num_city):
    tour = np.zeros(num_city, dtype=int)
    i_value, t_value = nonzero_indices
    tour[t_value] = i_value
    tour = np.append(tour, [tour[0]])
    return tour

tour = convert_to_tour(nonzero_indices, num_city)
tour

array([1, 2, 7, 5, 0, 4, 6, 3, 8, 9, 1])

plt.scatter(pos_x, pos_y)
plt.plot(pos_x[tour], pos_y[tour], "-")
plt.show()

Optimize the weights of penalties.

num_iter = 10
num_reads = 1
alpha = 1.2
results = []
A = np.max(distance)
Ahistory = [A]
best_tour = None
best_objective = np.inf
for i in range(num_iter):
    # re-put the parameters of TSP and make QUBO
    qubo, _ = qubo_builder.get_qubo_dict(multipliers={"one-city": A, "one-time": A})
    response = sampler.sample_qubo(qubo, num_reads=num_reads)
    decoded = jmt.core.pubo.decode_from_openjij(response, qubo_builder, compiled_instnace)
    feasibles = decoded.feasible()
    if len(feasibles.evaluation.objective) == 0:
        # If no feasible solution is obtained, increase the penalty parameter and try again.
        one_time_violaiton = np.mean(decoded.evaluation.constraint_violations["one-time"])
        one_city_violaiton = np.mean(decoded.evaluation.constraint_violations["one-city"])
        A += one_city_violaiton + one_time_violaiton
    else:
        argmin_e = np.argmin(feasibles.evaluation.objective)
        nonzero_indices, _, _ = feasibles.record.solution["x"][argmin_e]
        tour = convert_to_tour(nonzero_indices, num_city)

        if feasibles.evaluation.objective[argmin_e] < best_objective:
            best_objective = feasibles.evaluation.objective[argmin_e]
            best_tour = tour

        dij_dik = []
        for t, _i in enumerate(tour):
            _j = tour[(t-1)%num_city]
            _dij = distance[_i, _j]
            _k = tour[(t+1)%num_city]
            _dik = distance[_i, _k]
            dij_dik.append(_dij + _dik)
        # Update rule using the structure of TSP
        A = alpha * np.max(dij_dik)/2
    Ahistory.append(A)
    results.append(decoded)

objective = [np.mean(r.evaluation.objective) for r in results]
fig = plt.figure()

ax1 = fig.add_subplot(111)

ln1 = ax1.plot(Ahistory[:-1], "--", label="A")

ax2 = ax1.twinx()
ln2 = ax2.plot(objective, c="r", label="Objective")


h1, l1 = ax1.get_legend_handles_labels()
h2, l2 = ax2.get_legend_handles_labels()
ax1.legend(h1+h2, l1+l2, loc='upper right', fontsize=15)

ax1.set_xlabel('Step', fontsize=15)
ax1.set_ylabel(r'$A$', fontsize=15)
ax2.set_ylabel('Objective', fontsize=15)

Text(0, 0.5, 'Objective')

# plot the best solution
plt.scatter(pos_x, pos_y)
plt.plot(pos_x[best_tour], pos_y[best_tour], "-")
plt.show()