constraint-programming

Constraint Programming

You are an expert in constraint programming for supply chain optimization. Your goal is to help solve complex combinatorial problems using constraint propagation, global constraints, and intelligent search strategies that excel where traditional MIP struggles.

Initial Assessment

Before applying constraint programming, understand:

1. Problem Characteristics

   • Type? (scheduling, allocation, sequencing, assignment)
   • Constraints complex or logical? (if-then, all-different, etc.)
   • Variables: discrete domains?
   • Many feasibility constraints vs. optimization?

2. Why Constraint Programming

   • MIP formulation too weak?
   • Scheduling with disjunctive resources?
   • Complex logical constraints?
   • Need to find feasible solution quickly?

3. Problem Size

   • Number of variables?
   • Domain sizes?
   • Number of constraints?
   • Time limit for solving?

4. Technical Environment

   • CP solver available? (OR-Tools, Gecode, MiniZinc)
   • Need integrated with other systems?
   • Optimization or just feasibility?

---

Constraint Programming Fundamentals

Core Concepts

Variables: Decision variables with discrete domains

# Example: Variable x can take values 1, 2, or 3
x ∈ {1, 2, 3}

Constraints: Relations that must hold between variables

# Example: x + y ≤ 10
# Example: AllDifferent([x, y, z])

Propagation: Reduce variable domains based on constraints

# If x + y ≤ 10 and x = 8, then propagate y ≤ 2

Search: Systematic exploration with backtracking

# Try x = 1, check consistency, recurse
# If inconsistent, backtrack and try x = 2

---

Job Shop Scheduling with CP

Implementation with OR-Tools CP-SAT

from ortools.sat.python import cp_model
import matplotlib.pyplot as plt
import numpy as np
from typing import List, Tuple, Dict

class CPJobShopScheduler:
    """
    Job Shop Scheduling using Constraint Programming

    Problem: Schedule jobs on machines minimizing makespan
    Each job has sequence of operations on specific machines
    """

    def __init__(self,
                 jobs: List[List[Tuple[int, int]]],
                 time_limit_seconds: int = 60):
        """
        Initialize CP Job Shop Scheduler

        jobs: list of jobs, each job is list of (machine, duration)
              Example: [[(0, 3), (1, 2), (2, 2)],  # Job 0
                       [(0, 2), (2, 1), (1, 4)]]   # Job 1
        """

        self.jobs = jobs
        self.n_jobs = len(jobs)
        self.n_machines = max(max(op[0] for op in job) for job in jobs) + 1
        self.time_limit = time_limit_seconds

        # CP model
        self.model = cp_model.CpModel()

        # Decision variables
        self.job_starts = []  # [job][operation] -> start time variable
        self.job_ends = []    # [job][operation] -> end time variable
        self.job_intervals = []  # [job][operation] -> interval variable

        self.makespan = None
        self.solution = None

    def build_model(self):
        """Build CP model for job shop scheduling"""

        # Upper bound on time horizon
        horizon = sum(duration for job in self.jobs
                     for machine, duration in job)

        print(f"Building CP model...")
        print(f"Jobs: {self.n_jobs}, Machines: {self.n_machines}")
        print(f"Time Horizon: {horizon}")

        # Create variables for each operation
        for job_id, job in enumerate(self.jobs):
            job_start_vars = []
            job_end_vars = []
            job_interval_vars = []

            for op_id, (machine, duration) in enumerate(job):
                # Suffix for variable names
                suffix = f'_j{job_id}_o{op_id}_m{machine}'

                # Start time variable
                start_var = self.model.NewIntVar(0, horizon, f'start{suffix}')

                # End time variable
                end_var = self.model.NewIntVar(0, horizon, f'end{suffix}')

                # Interval variable (start, duration, end)
                interval_var = self.model.NewIntervalVar(
                    start_var, duration, end_var, f'interval{suffix}'
                )

                job_start_vars.append(start_var)
                job_end_vars.append(end_var)
                job_interval_vars.append(interval_var)

            self.job_starts.append(job_start_vars)
            self.job_ends.append(job_end_vars)
            self.job_intervals.append(job_interval_vars)

        # Precedence constraints: operations within job must be sequential
        for job_id in range(self.n_jobs):
            for op_id in range(len(self.jobs[job_id]) - 1):
                self.model.Add(
                    self.job_ends[job_id][op_id] <=
                    self.job_starts[job_id][op_id + 1]
                )

        # Disjunctive constraints: operations on same machine cannot overlap
        machine_to_intervals = [[] for _ in range(self.n_machines)]

        for job_id, job in enumerate(self.jobs):
            for op_id, (machine, duration) in enumerate(job):
                machine_to_intervals[machine].append(
                    self.job_intervals[job_id][op_id]
                )

        # NoOverlap constraint for each machine
        for machine in range(self.n_machines):
            if machine_to_intervals[machine]:
                self.model.AddNoOverlap(machine_to_intervals[machine])

        # Objective: minimize makespan
        self.makespan = self.model.NewIntVar(0, horizon, 'makespan')

        # Makespan is max end time of all jobs
        for job_id in range(self.n_jobs):
            last_op_idx = len(self.jobs[job_id]) - 1
            self.model.Add(
                self.makespan >= self.job_ends[job_id][last_op_idx]
            )

        self.model.Minimize(self.makespan)

        print("CP model built successfully!")

    def solve(self) -> Dict:
        """
        Solve the CP model

        Returns: solution dictionary
        """

        print(f"\nSolving with time limit: {self.time_limit}s...")

        # Create solver
        solver = cp_model.CpSolver()
        solver.parameters.max_time_in_seconds = self.time_limit

        # Optional: set number of workers for parallel solving
        solver.parameters.num_search_workers = 8

        # Solve
        status = solver.Solve(self.model)

        # Extract solution
        if status == cp_model.OPTIMAL or status == cp_model.FEASIBLE:
            print(f"\nSolution found!")
            print(f"Status: {'OPTIMAL' if status == cp_model.OPTIMAL else 'FEASIBLE'}")
            print(f"Makespan: {solver.Value(self.makespan)}")
            print(f"Solve time: {solver.WallTime():.2f}s")

            # Extract schedule
            schedule = []
            for job_id in range(self.n_jobs):
                for op_id, (machine, duration) in enumerate(self.jobs[job_id]):
                    start = solver.Value(self.job_starts[job_id][op_id])
                    end = solver.Value(self.job_ends[job_id][op_id])

                    schedule.append({
                        'job': job_id,
                        'operation': op_id,
                        'machine': machine,
                        'start': start,
                        'end': end,
                        'duration': duration
                    })

            self.solution = {
                'status': 'OPTIMAL' if status == cp_model.OPTIMAL else 'FEASIBLE',
                'makespan': solver.Value(self.makespan),
                'schedule': schedule,
                'solve_time': solver.WallTime(),
                'lower_bound': solver.BestObjectiveBound(),
                'gap': (solver.Value(self.makespan) - solver.BestObjectiveBound()) /
                       solver.Value(self.makespan) * 100
            }

            return self.solution

        else:
            print("No solution found!")
            return {
                'status': 'INFEASIBLE' if status == cp_model.INFEASIBLE else 'UNKNOWN',
                'makespan': None,
                'schedule': [],
                'solve_time': solver.WallTime()
            }

    def plot_gantt(self):
        """Visualize schedule as Gantt chart"""

        if not self.solution or not self.solution['schedule']:
            print("No solution to visualize!")
            return

        schedule = self.solution['schedule']

        fig, ax = plt.subplots(figsize=(14, 8))

        colors = plt.cm.Set3(np.linspace(0, 1, self.n_jobs))

        for task in schedule:
            ax.barh(
                task['machine'],
                task['duration'],
                left=task['start'],
                height=0.6,
                color=colors[task['job']],
                edgecolor='black',
                linewidth=1.5
            )

            # Add job label
            ax.text(
                task['start'] + task['duration'] / 2,
                task['machine'],
                f"J{task['job']}\nO{task['operation']}",
                ha='center',
                va='center',
                fontsize=9,
                fontweight='bold'
            )

        ax.set_xlabel('Time', fontsize=12)
        ax.set_ylabel('Machine', fontsize=12)
        ax.set_title(
            f"Job Shop Schedule - Makespan: {self.solution['makespan']}",
            fontsize=14
        )
        ax.set_yticks(range(self.n_machines))
        ax.set_yticklabels([f'M{i}' for i in range(self.n_machines)])
        ax.grid(True, axis='x', alpha=0.3)

        # Legend
        legend_elements = [
            plt.Rectangle((0,0),1,1, fc=colors[i],
                         edgecolor='black', label=f'Job {i}')
            for i in range(self.n_jobs)
        ]
        ax.legend(handles=legend_elements, loc='upper right')

        plt.tight_layout()
        plt.show()


# Example usage
if __name__ == "__main__":
    # Classic 6x6 job shop problem (Fisher and Thompson, 1963)
    jobs = [
        [(0, 1), (1, 3), (2, 6), (3, 7), (4, 3), (5, 6)],
        [(1, 8), (0, 5), (2, 10), (4, 10), (5, 10), (3, 4)],
        [(0, 5), (1, 4), (2, 8), (4, 9), (3, 1), (5, 7)],
        [(1, 5), (0, 5), (2, 5), (3, 3), (4, 8), (5, 9)],
        [(2, 9), (1, 3), (3, 5), (4, 4), (5, 3), (0, 1)],
        [(1, 3), (2, 3), (4, 9), (3, 10), (5, 4), (0, 1)]
    ]

    # Create scheduler
    scheduler = CPJobShopScheduler(jobs, time_limit_seconds=30)

    # Build and solve
    scheduler.build_model()
    result = scheduler.solve()

    # Print results
    if result['status'] in ['OPTIMAL', 'FEASIBLE']:
        print(f"\nMakespan: {result['makespan']}")
        print(f"Lower Bound: {result['lower_bound']}")
        print(f"Gap: {result['gap']:.2f}%")

        # Visualize
        scheduler.plot_gantt()