Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Lab 2: Constraint Satisfaction Problems (CSPs) - Assignment

Open In Colab
# Setup - Run this cell first
from constraint import Problem, AllDifferentConstraint, ExactSumConstraint
import itertools

print("✓ Setup complete! python-constraint library loaded.")

Exercise 1: Sudoku with Binary Variables (40 points)

In the lecture, we implemented Sudoku where each cell variable had a domain of {1, 2, ..., 9}. Now we’ll use a binary encoding where:

  • Variables: Xr,c,vX_{r,c,v} for each row rr, column cc, and value vv (where r,c,v{0,1,...,8}r, c, v \in \{0, 1, ..., 8\})

  • Domain: {0,1}\{0, 1\} for each variable

  • Meaning: Xr,c,v=1X_{r,c,v} = 1 if cell (r,c)(r, c) contains value v+1v+1, otherwise Xr,c,v=0X_{r,c,v} = 0

Constraints for Binary Sudoku

  1. Each cell has exactly one value:

    v=08Xr,c,v=1,r,c{0,...,8}\sum_{v=0}^{8} X_{r,c,v} = 1, \quad \forall r, c \in \{0, ..., 8\}
  2. Each value appears exactly once in each row:

    c=08Xr,c,v=1,r,v{0,...,8}\sum_{c=0}^{8} X_{r,c,v} = 1, \quad \forall r, v \in \{0, ..., 8\}
  3. Each value appears exactly once in each column:

    r=08Xr,c,v=1,c,v{0,...,8}\sum_{r=0}^{8} X_{r,c,v} = 1, \quad \forall c, v \in \{0, ..., 8\}
  4. Each value appears exactly once in each 3×3 subgrid:

    i=02j=02X3m+i,3n+j,v=1,m,n{0,1,2},v{0,...,8}\sum_{i=0}^{2} \sum_{j=0}^{2} X_{3m+i, 3n+j, v} = 1, \quad \forall m, n \in \{0, 1, 2\}, \forall v \in \{0, ..., 8\}
def solve_sudoku_binary(puzzle):
    """
    Solve Sudoku using binary decision variables.
    
    Variables: X[r,c,v] = 1 if cell (r,c) contains value v+1
    
    Args:
        puzzle: 9x9 list where 0 represents empty cells
    
    Returns:
        Solved 9x9 grid or None if no solution
    """
    problem = Problem()
    
    # TODO: Create binary variables X[r,c,v] for r,c,v in 0..8
    # Variable name format: (row, col, value)
    # Hint: For each cell (r,c) and each possible value v (0-8):
    #   - If puzzle[r][c] == 0: empty cell, add binary variable with domain [0, 1]
    #   - If puzzle[r][c] == v + 1: given cell with this value, fix to [1]
    #   - Otherwise: given cell with different value, fix to [0]
    # YOUR CODE HERE
    pass
    
    # TODO: Constraint 1 - Each cell has exactly one value
    # For each cell (r,c), sum of X[r,c,v] over all v must equal 1
    # Use: problem.addConstraint(ExactSumConstraint(1), cell_vars)
    # YOUR CODE HERE
    pass
    
    # TODO: Constraint 2 - Each value appears exactly once in each row
    # For each row r and value v, sum of X[r,c,v] over all c must equal 1
    # YOUR CODE HERE
    pass
    
    # TODO: Constraint 3 - Each value appears exactly once in each column
    # For each column c and value v, sum of X[r,c,v] over all r must equal 1
    # YOUR CODE HERE
    pass
    
    # TODO: Constraint 4 - Each value appears exactly once in each 3x3 subgrid
    # For each box (box_r, box_c) and value v, sum of X[r,c,v] in the box must equal 1
    # Box cells: (box_r * 3 + i, box_c * 3 + j) for i,j in 0..2
    # YOUR CODE HERE
    pass
    
    # Solve
    solution = problem.getSolution()
    
    if solution:
        # Convert binary solution back to 9x9 grid
        result = [[0] * 9 for _ in range(9)]
        for r in range(9):
            for c in range(9):
                for v in range(9):
                    if solution[(r, c, v)] == 1:
                        result[r][c] = v + 1
                        break
        return result
    return None

print("✓ solve_sudoku_binary function defined!")
# Test the binary Sudoku solver
puzzle = [
    [5, 3, 0, 0, 7, 0, 0, 0, 0],
    [6, 0, 0, 1, 9, 5, 0, 0, 0],
    [0, 9, 8, 0, 0, 0, 0, 6, 0],
    [8, 0, 0, 0, 6, 0, 0, 0, 3],
    [4, 0, 0, 8, 0, 3, 0, 0, 1],
    [7, 0, 0, 0, 2, 0, 0, 0, 6],
    [0, 6, 0, 0, 0, 0, 2, 8, 0],
    [0, 0, 0, 4, 1, 9, 0, 0, 5],
    [0, 0, 0, 0, 8, 0, 0, 7, 9]
]

print("Original Puzzle:")
for row in puzzle:
    print(row)

print("\nSolving with binary variables...")
solution = solve_sudoku_binary(puzzle)

if solution:
    print("\n✓ Sudoku solved successfully!")
    print("\nSolution:")
    for row in solution:
        print(row)
else:
    print("No solution found.")

Exercise 1 Questions


Exercise 2: Traveling Salesman Problem (TSP) (30 points)

The Traveling Salesman Problem (TSP) asks: given a list of cities and distances between them, what is the shortest route that visits each city exactly once and returns to the starting city?

Binary CSP Formulation

  • Variables: Xi,jX_{i,j} for each pair of cities (i,j)(i, j) where iji \neq j

  • Domain: {0,1}\{0, 1\}

  • Meaning: Xi,j=1X_{i,j} = 1 if the tour goes directly from city ii to city jj

Constraints

  1. Leave each city exactly once:

    jiXi,j=1,i\sum_{j \neq i} X_{i,j} = 1, \quad \forall i
  2. Enter each city exactly once:

    ijXi,j=1,j\sum_{i \neq j} X_{i,j} = 1, \quad \forall j
  3. Subtour elimination: Prevent disconnected cycles (we’ll use a position-based approach)

  4. Total distance constraint:

    ijidi,jXi,jC\sum_{i} \sum_{j \neq i} d_{i,j} \cdot X_{i,j} \leq C
def solve_tsp(distances, max_cost=None):
    """
    Solve the Traveling Salesman Problem using CSP with binary variables.
    
    Uses a position-based formulation to avoid subtours:
    - Variables: pos[i] = position of city i in the tour (0 to n-1)
    - Variables: X[i,j] = 1 if edge (i,j) is in the tour
    
    Args:
        distances: dict mapping (i,j) to distance between cities i and j
        max_cost: optional maximum total tour cost
    
    Returns:
        (tour, total_distance) or None if no solution
    """
    # Extract cities from distance matrix
    cities = set()
    for (i, j) in distances:
        cities.add(i)
        cities.add(j)
    cities = sorted(list(cities))
    n = len(cities)
    
    problem = Problem()
    
    # TODO: Position variables - assign position in tour to each city
    # Fix city 0 at position 0 (to break symmetry)
    # Other cities get positions from range(1, n)
    # Hint: problem.addVariable(('pos', city), [0]) for first city
    #       problem.addVariable(('pos', city), range(1, n)) for others
    # YOUR CODE HERE
    pass
    
    # TODO: All positions must be different
    # Use: problem.addConstraint(AllDifferentConstraint(), position_vars)
    # YOUR CODE HERE
    pass
    
    # TODO: Edge variables X[i,j] = 1 if we travel from i to j
    # For each pair of distinct cities (i,j), add binary variable ('edge', i, j)
    # YOUR CODE HERE
    pass
    
    # TODO: Constraint - Leave each city exactly once
    # For each city i, sum of X[i,j] over all j != i must equal 1
    # YOUR CODE HERE
    pass
    
    # TODO: Constraint - Enter each city exactly once
    # For each city j, sum of X[i,j] over all i != j must equal 1
    # YOUR CODE HERE
    pass
    
    # TODO: Link edge variables to position variables
    # If X[i,j] = 1, then pos[j] = (pos[i] + 1) mod n
    # Define: edge_position_constraint(edge_val, pos_i, pos_j, n=n)
    #   - If edge_val == 1: return (pos_i + 1) % n == pos_j
    #   - Otherwise: return True (no constraint)
    # Apply this constraint for each edge variable
    # YOUR CODE HERE
    pass
    
    # Optional: Maximum cost constraint (already implemented)
    if max_cost is not None:
        def cost_constraint(*edge_values):
            total = 0
            edge_list = [(i, j) for i in cities for j in cities if i != j]
            for idx, (i, j) in enumerate(edge_list):
                if edge_values[idx] == 1:
                    total += distances.get((i, j), distances.get((j, i), float('inf')))
            return total <= max_cost
        
        edge_vars = [('edge', i, j) for i in cities for j in cities if i != j]
        problem.addConstraint(cost_constraint, edge_vars)
    
    # Solve
    solution = problem.getSolution()
    
    if solution:
        # Reconstruct tour from positions
        tour = [None] * n
        for city in cities:
            pos = solution[('pos', city)]
            tour[pos] = city
        
        # Calculate total distance
        total_dist = 0
        for idx in range(n):
            i = tour[idx]
            j = tour[(idx + 1) % n]
            total_dist += distances.get((i, j), distances.get((j, i), 0))
        
        return tour, total_dist
    
    return None

print("✓ solve_tsp function defined!")
# Test TSP with a small example: 4 cities
# 
#     A ---10--- B
#     |  \      |
#    15   25   35
#     |      \  |
#     D ---20--- C
#

distances = {
    ('A', 'B'): 10, ('B', 'A'): 10,
    ('A', 'C'): 25, ('C', 'A'): 25,
    ('A', 'D'): 15, ('D', 'A'): 15,
    ('B', 'C'): 35, ('C', 'B'): 35,
    ('B', 'D'): 30, ('D', 'B'): 30,
    ('C', 'D'): 20, ('D', 'C'): 20,
}

print("Cities: A, B, C, D")
print("\nDistance Matrix:")
print("     A    B    C    D")
print(f"A    -   10   25   15")
print(f"B   10    -   35   30")
print(f"C   25   35    -   20")
print(f"D   15   30   20    -")

print("\nSolving TSP...")
result = solve_tsp(distances)

if result:
    tour, total_dist = result
    print(f"\n✓ Optimal tour found!")
    print(f"Tour: {' → '.join(tour)} → {tour[0]}")
    print(f"Total distance: {total_dist}")
else:
    print("No solution found.")
# Test with a larger example: 5 cities with coordinates
import math

# City coordinates
city_coords = {
    'Paris': (2.35, 48.85),
    'London': (-0.12, 51.51),
    'Berlin': (13.40, 52.52),
    'Rome': (12.50, 41.90),
    'Madrid': (-3.70, 40.42)
}

# Calculate Euclidean distances (simplified, not actual km)
def euclidean_dist(c1, c2):
    return round(math.sqrt((c1[0]-c2[0])**2 + (c1[1]-c2[1])**2), 2)

distances_europe = {}
cities_list = list(city_coords.keys())
for i, city1 in enumerate(cities_list):
    for city2 in cities_list:
        if city1 != city2:
            distances_europe[(city1, city2)] = euclidean_dist(
                city_coords[city1], city_coords[city2]
            )

print("European Cities TSP")
print("=" * 40)
print("\nCity Coordinates:")
for city, coord in city_coords.items():
    print(f"  {city}: {coord}")

print("\nSolving TSP...")
result = solve_tsp(distances_europe)

if result:
    tour, total_dist = result
    print(f"\n✓ Tour found!")
    print(f"Tour: {' → '.join(tour)} → {tour[0]}")
    print(f"Total distance: {total_dist:.2f} units")
else:
    print("No solution found.")

Exercise 2 Questions


Exercise 3: Graph Coloring Problem (30 points)

The Graph Coloring Problem asks: given a graph, assign colors to vertices such that no two adjacent vertices share the same color, using at most kk colors.

Binary CSP Formulation

  • Variables: Xv,cX_{v,c} for each vertex vv and color cc

  • Domain: {0,1}\{0, 1\}

  • Meaning: Xv,c=1X_{v,c} = 1 if vertex vv is assigned color cc

Constraints

  1. Each vertex gets exactly one color:

    c=1kXv,c=1,vV\sum_{c=1}^{k} X_{v,c} = 1, \quad \forall v \in V
  2. Adjacent vertices have different colors:

    Xu,c+Xv,c1,(u,v)E,c{1,...,k}X_{u,c} + X_{v,c} \leq 1, \quad \forall (u,v) \in E, \forall c \in \{1, ..., k\}

This means if vertex uu has color cc, then vertex vv cannot have color cc (and vice versa).

def solve_graph_coloring(vertices, edges, num_colors):
    """
    Solve the Graph Coloring Problem using CSP with binary variables.
    
    Variables: X[v,c] = 1 if vertex v has color c
    
    Args:
        vertices: list of vertex names
        edges: list of (u, v) tuples representing edges
        num_colors: number of colors available
    
    Returns:
        dict mapping vertex to color, or None if no solution
    """
    problem = Problem()
    colors = list(range(1, num_colors + 1))
    
    # TODO: Create binary variables X[v,c] for each vertex and color
    # For each vertex v and color c, add variable (v, c) with domain [0, 1]
    # YOUR CODE HERE
    pass
    
    # TODO: Constraint 1 - Each vertex gets exactly one color
    # For each vertex v, sum of X[v,c] over all colors must equal 1
    # Use: problem.addConstraint(ExactSumConstraint(1), vertex_vars)
    # YOUR CODE HERE
    pass
    
    # TODO: Constraint 2 - Adjacent vertices cannot have same color
    # For each edge (u,v) and each color c: X[u,c] + X[v,c] <= 1
    # Define: adjacent_constraint(xu_c, xv_c): return xu_c + xv_c <= 1
    # Apply to variables [(u, c), (v, c)] for each edge and color
    # YOUR CODE HERE
    pass
    
    # Solve
    solution = problem.getSolution()
    
    if solution:
        # Extract coloring from binary solution
        coloring = {}
        for v in vertices:
            for c in colors:
                if solution[(v, c)] == 1:
                    coloring[v] = c
                    break
        return coloring
    
    return None

print("✓ solve_graph_coloring function defined!")
# Test 1: Simple triangle graph (requires 3 colors)
#
#      A
#     / \
#    B---C
#

print("Test 1: Triangle Graph")
print("=" * 40)

vertices_triangle = ['A', 'B', 'C']
edges_triangle = [('A', 'B'), ('B', 'C'), ('A', 'C')]

print("Graph: A-B-C-A (triangle)")
print("\nTrying with 2 colors...")
result = solve_graph_coloring(vertices_triangle, edges_triangle, 2)
if result:
    print(f"✓ Coloring found: {result}")
else:
    print("✗ No solution with 2 colors (expected!)")

print("\nTrying with 3 colors...")
result = solve_graph_coloring(vertices_triangle, edges_triangle, 3)
if result:
    print(f"✓ Coloring found: {result}")
# Test 2: Map coloring - Simplified map of some European countries
#
#        Norway
#          |
#       Sweden -- Finland
#          |
#       Denmark
#
# Adjacencies: Norway-Sweden, Sweden-Finland, Sweden-Denmark

print("\nTest 2: Nordic Countries Map Coloring")
print("=" * 40)

countries = ['Norway', 'Sweden', 'Finland', 'Denmark']
borders = [
    ('Norway', 'Sweden'),
    ('Sweden', 'Finland'),
    ('Sweden', 'Denmark')
]

print("Countries:", countries)
print("Borders:", borders)

print("\nTrying with 2 colors...")
result = solve_graph_coloring(countries, borders, 2)
if result:
    print(f"✓ Coloring found: {result}")
else:
    print("✗ No solution with 2 colors")

print("\nTrying with 3 colors...")
result = solve_graph_coloring(countries, borders, 3)
if result:
    print(f"✓ Coloring found: {result}")
# Test 3: Petersen Graph (famous graph requiring exactly 3 colors)
# The Petersen graph has 10 vertices and 15 edges

print("\nTest 3: Petersen Graph (Classic Graph Theory Example)")
print("=" * 40)

# Outer pentagon: 0-1-2-3-4-0
# Inner star: 5-7-9-6-8-5
# Spokes: 0-5, 1-6, 2-7, 3-8, 4-9

petersen_vertices = list(range(10))
petersen_edges = [
    # Outer pentagon
    (0, 1), (1, 2), (2, 3), (3, 4), (4, 0),
    # Inner star (pentagram)
    (5, 7), (7, 9), (9, 6), (6, 8), (8, 5),
    # Spokes connecting outer to inner
    (0, 5), (1, 6), (2, 7), (3, 8), (4, 9)
]

print(f"Vertices: {petersen_vertices}")
print(f"Number of edges: {len(petersen_edges)}")

print("\nTrying with 2 colors...")
result = solve_graph_coloring(petersen_vertices, petersen_edges, 2)
if result:
    print(f"✓ Coloring found!")
else:
    print("✗ No solution with 2 colors")

print("\nTrying with 3 colors...")
result = solve_graph_coloring(petersen_vertices, petersen_edges, 3)
if result:
    print(f"✓ Coloring found!")
    # Display by color groups
    color_groups = {1: [], 2: [], 3: []}
    for v, c in result.items():
        color_groups[c].append(v)
    for c, verts in color_groups.items():
        print(f"  Color {c}: vertices {verts}")

Exercise 3 Questions

# TODO: Test Australia map coloring
# Uncomment and complete the code below

australia_states = ['WA', 'NT', 'SA', 'QLD', 'NSW', 'VIC', 'TAS']
australia_borders = [
    ('WA', 'NT'), ('WA', 'SA'),
    ('NT', 'SA'), ('NT', 'QLD'),
    ('SA', 'QLD'), ('SA', 'NSW'), ('SA', 'VIC'),
    ('QLD', 'NSW'),
    ('NSW', 'VIC')
]

# Your code here:
# Test with 3 colors
# result_3 = solve_graph_coloring(...)

# Test with 4 colors  
# result_4 = solve_graph_coloring(...)