Skip to content

Modeling Complex Optimization Problems

In this tutorial we will demonstrate how to model a more complex NP-hard optimization problem, the Traveling Salesperson Problem (TSP). Together we will walk through all the necessary modeling steps and explain how to handle and implement the variables, constraints, environments, and objective functions for an advanced Model. Before you begin, we recommend reviewing the Model Introduction to get familiar with the core concepts.

Example: Modeling the Traveling Salesperson Problem (TSP)

The Traveling Salesperson Problem (TSP) is a classical optimization problem that involves finding the shortest possible route visiting each city exactly once and returning to the original starting point. Although easy to understand, the TSP becomes computationally challenging: as the number of cities grows, the number of possible routes increases exponentially, making exact solutions infeasible for larger scenarios. However, it serves a wide range of practical applications in logistics planning, transportation networks, microchip circuit design, manufacturing, and DNA sequencing.

For this use case we will define a network of 7 cities. The network of cities can be visualized as a graph, with cities as nodes and distances as weighted edges between them.

import networkx as nx
import matplotlib.pyplot as plt
from itertools import combinations

# Graph initialization
G = nx.Graph()

# Define cities as nodes and add them to the Graph
cities = ["A", "B", "C", "D", "E", "F", "G"]
G.add_nodes_from(cities)

# Distances as edges
# (Missing edges implicitly have a high cost)
distances = {
    ("A", "B"): 29, ("A", "C"): 20, ("A", "D"): 21, ("A", "E"): 16, ("A", "F"): 31, ("A", "G"): 100,
    ("B", "C"): 15, ("B", "D"): 29, ("B", "E"): 28, ("B", "F"): 40, ("B", "G"): 72, ("C", "D"): 15,
    ("C", "E"): 14, ("C", "F"): 25, ("C", "G"): 81, ("D", "E"): 4,  ("D", "F"): 12, ("D", "G"): 92,
    ("E", "F"): 16, ("E", "G"): 94, ("F", "G"): 95
}

# Add distances to the Graph as edges
# Add non-existing connections as edges with a default distance of 10000
for city1, city2 in combinations(cities, 2):
    distance = distances.get((city1, city2)) or distances.get((city2, city1)) or 10000
    G.add_edge(city1, city2, weight=distance)

# Graph visualization
pos = nx.circular_layout(G)
nx.draw(G, pos, with_labels=True, node_color='lightblue', node_size=700, font_weight='bold')
nx.draw_networkx_edge_labels(G, pos, edge_labels=nx.get_edge_attributes(G, 'weight'))
plt.title("Traveling Salesperson Problem")
plt.axis('off')
plt.show()

Creating the Model

Let's build our optimization model using the Model class. To begin, initialize the Model instance and optionally assign it a clear name that is descriptive of our task.

from luna_quantum import Model, Sense

# Create the model
model = Model("TSP Example")

# Set the Model to minimize its objective function
model.set_sense(Sense.Min)

Defining the Variables and an Environment

In the Traveling Salesperson Problem (TSP) example, we formulate the task as a binary optimization problem, focusing on determining the optimal visiting order of cities. Since each city must appear exactly once in the route, the encoding assigns each city to a unique position within the tour. Thus, for a problem involving 7 cities, there are precisely 7 positions in the tour.

To represent this mathematically, we define binary decision variables x[u][j] ∈ {0, 1} for each possible pairing of the city u and position j. A decision variable x[u][j] = 1 indicates that city u is visited at position j. Consequently, for this instance with 7 cities and 7 positions, we require a total of 49 binary decision variables.

from luna_quantum import Variable, Vtype

# Declare decision variables within the model's environment
variables = {}
with model.environment:
    for node in G.nodes():
        for path_order in range(len(G.nodes())):
            var_name = f"{node}_{path_order}"
            variables[var_name] = Variable(name=var_name, vtype=Vtype.Binary)

Creating an Objective Function

In the TSP, the objective function represents the total distance traveled, which we aim to minimize. To construct this objective mathematically, we sum over all consecutive position pairs (j, j+1) within the tour. For each of these pairs, we then sum over all possible city pairs (u, v) to consider every potential transition from city u at position j to city v at position j+1. By multiplying each possible city-to-city transition by its corresponding distance and the associated decision variables, we obtain the total route distance. The resulting expression effectively captures the total distance contributed by each city transition along the entire route.

num_nodes = len(G.nodes())
objective_terms = []

for u, v, weight in G.edges(data='weight'):
    for pos in range(num_nodes):
        u_at_pos = variables[f"{u}_{pos}"]
        v_next_pos = variables[f"{v}_{(pos + 1) % num_nodes}"]
        objective_terms.append(u_at_pos * v_next_pos * weight)

model.objective = sum(objective_terms)

Adding Constraints

Constraints restrict the feasible solution space defined by the decision variables. In the TSP, the following logical constraints must hold:

While the objective function accurately represents the total traveled distance we aim to minimize, it cannot operate effectively in isolation. Without constraints, the optimization model might produce unreasonable or invalid solutions, such as assigning multiple cities to the same position or skipping cities entirely. To ensure a valid solution, the objective function must be accompanied by carefully defined constraints:

1. Each city u must be visited exactly once, regardless of the position j. This means that for each city u, the sum of decision variables over all positions j must equal one:

for node in G.nodes():
    city_expression = sum(variables[f"{node}_{path_order}"] for path_order in range(len(G.nodes())))
    model.add_constraint(city_expression == 1, f"city_constraint_{node}")

2. Each position j must be assigned exactly one city u. That is, for each position j, the sum of decision variables over all cities must equal one:

for path_order in range(len(G.nodes())):
    position_expression = sum(variables[f"{node}_{path_order}"] for node in G.nodes())
    model.constraints += position_expression == 1, f"path_order_constraint_{path_order}"

Great! We have just managed to create our first TSP optimization problem. We have defined the decision variables of the Model inside its environment and have built an objective function together with the needed constraints.

Access and Analyze the TSP Model

Now that we have built an optimization model, we can inspect different aspects of the Model.

The summary, returns information of the objective function, the constraints, bounds and types of the variables. Check out the Model's summary with: 

print(model)

Output: 

Model: TSP Example
Minimize
  29 * A_0 * B_1 + 20 * A_0 * C_1 + 21 * A_0 * D_1 + 16 * A_0 * E_1 
  + 31 * A_0 * F_1 + 100 * A_0 * G_1 + 29 * A_1 * B_2 + 20 * A_1 * C_2 
  + 21 * A_1 * D_2 + 16 * A_1 * E_2 + 31 * A_1 * F_2 + 100 * A_1 * G_2 
  + 29 * A_2 * B_3 + 20 * A_2 * C_3 + 21 * A_2 * D_3 + 16 * A_2 * E_3 
  ...
  + 16 * E_1 * F_2 + 94 * E_1 * G_2 + 16 * E_2 * F_3 + 94 * E_2 * G_3 
  + 16 * E_3 * F_4 + 94 * E_3 * G_4 + 16 * E_4 * F_5 + 94 * E_4 * G_5 
  + 16 * E_5 * F_6 + 94 * E_5 * G_6 + 16 * E_6 * F_0 + 94 * E_6 * G_0 
  + 95 * F_0 * G_1 + 95 * F_1 * G_2 + 95 * F_2 * G_3 + 95 * F_3 * G_4 
  + 95 * F_4 * G_5 + 95 * F_5 * G_6 + 95 * F_6 * G_0
Subject To
  city_constraint_A: A_0 + A_1 + A_2 + A_3 + A_4 + A_5 + A_6 == 1
  city_constraint_B: B_0 + B_1 + B_2 + B_3 + B_4 + B_5 + B_6 == 1
  city_constraint_C: C_0 + C_1 + C_2 + C_3 + C_4 + C_5 + C_6 == 1
  city_constraint_D: D_0 + D_1 + D_2 + D_3 + D_4 + D_5 + D_6 == 1
  city_constraint_E: E_0 + E_1 + E_2 + E_3 + E_4 + E_5 + E_6 == 1
  city_constraint_F: F_0 + F_1 + F_2 + F_3 + F_4 + F_5 + F_6 == 1
  city_constraint_G: G_0 + G_1 + G_2 + G_3 + G_4 + G_5 + G_6 == 1
  path_order_constraint_0: A_0 + B_0 + C_0 + D_0 + E_0 + F_0 + G_0 == 1
  path_order_constraint_1: A_1 + B_1 + C_1 + D_1 + E_1 + F_1 + G_1 == 1
  path_order_constraint_2: A_2 + B_2 + C_2 + D_2 + E_2 + F_2 + G_2 == 1
  path_order_constraint_3: A_3 + B_3 + C_3 + D_3 + E_3 + F_3 + G_3 == 1
  path_order_constraint_4: A_4 + B_4 + C_4 + D_4 + E_4 + F_4 + G_4 == 1
  path_order_constraint_5: A_5 + B_5 + C_5 + D_5 + E_5 + F_5 + G_5 == 1
  path_order_constraint_6: A_6 + B_6 + C_6 + D_6 + E_6 + F_6 + G_6 == 1
Bounds
  0 <= A_0 <= 1
  0 <= A_1 <= 1
  0 <= A_2 <= 1
  0 <= A_3 <= 1
  0 <= A_4 <= 1
  0 <= A_5 <= 1
  0 <= A_6 <= 1
  0 <= B_0 <= 1
  0 <= B_1 <= 1
  ...
  0 <= F_6 <= 1
  0 <= G_0 <= 1
  0 <= G_1 <= 1
  0 <= G_2 <= 1
  0 <= G_3 <= 1
  0 <= G_4 <= 1
  0 <= G_5 <= 1
  0 <= G_6 <= 1
Binary
  A_0 A_1 A_2 A_3 A_4 A_5 A_6 B_0 B_1 B_2 B_3 B_4 B_5 B_6 C_0 C_1 C_2 C_3 C_4 
  C_5 C_6 D_0 D_1 D_2 D_3 D_4 D_5 D_6 E_0 E_1 E_2 E_3 E_4 E_5 E_6 F_0 F_1 F_2 
  F_3 F_4 F_5 F_6 G_0 G_1 G_2 G_3 G_4 G_5 G_6

We can further analyze the objective function by directly accessing it via model.objective. The expression object allows us to retrieve the number of distinct variables involved in it via expression.num_variables().

expression = model.objective
print(f"Number of variables included in the expression: {expression.num_variables()}")
print(f"Offset of expression: {expression.get_offset()}")
Output: 
Number of variables included in the expression: 49
Offset of expression: 0.0

From the output, we observe that the number of variables used in the objective is 49, which matches the total number of decision variables of the optimization model we created. This confirms that every variable contributes to the objective.

print(f"Total number of constraints = {model.num_constraints()}")
Output: 
Total number of constraints = 14

We can confirm what we already saw from the Model's summary. There are 14 constraints in total: 7 constraints ensuring each city is visited exactly once, and 7 constraints ensuring each position is occupied by exactly one city.

To explore constraints individually, you can iterate through model.constraints. Each constraint object includes a name attribute, which helps identify its role in an optimization model.

for idx, constraint in enumerate(model.constraints):
    print(f"{constraint.name} -> {constraint}")
Output: 
city_constraint_A -> A_0 + A_1 + A_2 + A_3 + A_4 + A_5 + A_6 == 1
city_constraint_B -> B_0 + B_1 + B_2 + B_3 + B_4 + B_5 + B_6 == 1
city_constraint_C -> C_0 + C_1 + C_2 + C_3 + C_4 + C_5 + C_6 == 1
city_constraint_D -> D_0 + D_1 + D_2 + D_3 + D_4 + D_5 + D_6 == 1
city_constraint_E -> E_0 + E_1 + E_2 + E_3 + E_4 + E_5 + E_6 == 1
city_constraint_F -> F_0 + F_1 + F_2 + F_3 + F_4 + F_5 + F_6 == 1
city_constraint_G -> G_0 + G_1 + G_2 + G_3 + G_4 + G_5 + G_6 == 1
path_order_constraint_0 -> A_0 + B_0 + C_0 + D_0 + E_0 + F_0 + G_0 == 1
path_order_constraint_1 -> A_1 + B_1 + C_1 + D_1 + E_1 + F_1 + G_1 == 1
path_order_constraint_2 -> A_2 + B_2 + C_2 + D_2 + E_2 + F_2 + G_2 == 1
path_order_constraint_3 -> A_3 + B_3 + C_3 + D_3 + E_3 + F_3 + G_3 == 1
path_order_constraint_4 -> A_4 + B_4 + C_4 + D_4 + E_4 + F_4 + G_4 == 1
path_order_constraint_5 -> A_5 + B_5 + C_5 + D_5 + E_5 + F_5 + G_5 == 1
path_order_constraint_6 -> A_6 + B_6 + C_6 + D_6 + E_6 + F_6 + G_6 == 1

You can also inspect how specific variables contribute to an expression by querying their associated coefficients. The methods expression.get_linear(x), expression.get_quadratic(x, y), and expression.get_higher_order((x, y, z, ...)) allow you to access the prefactors of linear, quadratic, and higher-order terms, respectively.

Since the encoding of the TSP objective contains only quadratic terms, we focus on those. For example, to understand how the variable A_0, representing city A at position 0, interacts with other variables, we use expression.get_quadratic(x, y) to evaluate all pairwise interactions with A_0. 

expression = model.objective
var_x = variables["A_0"]
for var_y in variables.values():
      quad_prefactor = expression.get_quadratic(var_x, var_y)
      if quad_prefactor != 0:
        print(f"var_x = {var_x}, var_y = {var_y}")
        print(f"Quadratic prefactor = {quad_prefactor}")

Output: 

var_x = A_0: Binary, var_y = B_1: Binary
Quadratic prefactor = 29.0
var_x = A_0: Binary, var_y = C_1: Binary
Quadratic prefactor = 20.0
var_x = A_0: Binary, var_y = D_1: Binary
Quadratic prefactor = 21.0
var_x = A_0: Binary, var_y = E_1: Binary
Quadratic prefactor = 16.0
var_x = A_0: Binary, var_y = F_1: Binary
Quadratic prefactor = 31.0
var_x = A_0: Binary, var_y = G_1: Binary
Quadratic prefactor = 100.0

The results confirm that A_0 only interacts with variables corresponding to other cities placed at the next tour position, 1. Further, the quadratic prefactors are equal to the connecting distance between the two investigated cities. This matches the expected structure of the TSP, where only transitions between consecutive tour positions contribute to the cost, and the cost is equal to the added distance which must be traveled from city to city.

Solving the Optimization Problem

Once your model is defined with an objective and constraints, you're ready to solve the optimization problem using your preferred algorithm and backend. If you want to start straight away using your existing models, use our Translators to transfer your models into luna_quantum Model's and start right away!