DEV Community: Thana B.

Implementation of Optimization of Red Blood Cell Shape Using Spherical Harmonics Parameterization

Thana B. — Sun, 06 Oct 2024 13:58:06 +0000

Introduction

Red blood cells (RBCs) exhibit a distinctive biconcave disc shape, essential for their role in gas exchange and traversing microvasculature. Modeling and optimizing this shape can provide valuable insights into the cell's mechanical properties and behaviors under physiological and pathological conditions.

This article presents a comprehensive mathematical model for optimizing the RBC shape using spherical harmonics parameterization under axial symmetry. We then provide a detailed Python implementation, explicitly mapping each mathematical expression to corresponding code blocks. The implementation utilizes the CasADi library for symbolic computations and numerical optimization.

Mathematical Formulation for RBC Shape Optimization with Fixed Posture and Symmetry

1. Decision Variables

Spherical Harmonics Coefficients (Axially Symmetric):

The shape of the RBC is represented using only the coefficients for axisymmetric spherical harmonics:

C_{L,0} \quad \text{for } L \in {0, 1, \ldots, L_{\text{max}}}

This reduces the number of decision variables significantly since only the coefficients corresponding to $K = 0$ are retained due to axial symmetry.

2. Objective Function

The objective is to minimize the total energy of the RBC membrane:

\text{Minimize} \quad E_{\text{total}}(C_{L,0}) = E_{\text{bending}}(C_{L,0}) + E_{\text{shear}}(C_{L,0}) + E_{\text{coupling}}(C_{L,0}) + E_{\text{viscous}}(C_{L,0}) + E_{\text{thermal}}(C_{L,0}) + E_{\text{hydro}}(C_{L,0})

Where each energy term is calculated using the reduced set of coefficients. The energy components are:

Bending Energy

E_{\text{bending}} = \sum_{i=1}^{N} \frac{k_b}{2} \left( 2H_i - C_0 \right)^2 \cdot A_i

$k_b$ : Bending modulus.
$H_i$ : Mean curvature at point $i$ .
$C_0$ : Spontaneous curvature (assumed zero in this model).
$A_i$ : Area element at point $i$ .

Shear Energy

E_{\text{shear}} = \sum_{i=1}^{N} \frac{K_s}{2} \left( \lambda_{1,i}^2 + \lambda_{2,i}^2 - 2 \right) \cdot A_i

$K_s$ : Shear modulus.
$λ1,i,λ2,i\lambda_{1,i}, \lambda_{2,i}$ : Principal stretches at point $i$ .

Coupling Energy

E_{\text{coupling}} = \sum_{i=1}^{N} \gamma (H_i - C_0) \left( \lambda_{1,i}^2 + \lambda_{2,i}^2 \right) \cdot A_i

$γ\gamma$ : Coupling constant.

Viscous Energy

E_{\text{viscous}} = \sum_{i=1}^{N} \eta \left( \frac{\partial H_i}{\partial t} \right)^2 \cdot A_i

$η\eta$ : Viscosity coefficient.
$∂Hi∂t\frac{\partial H_i}{\partial t}$ : Temporal change in curvature.

Thermal Fluctuation Energy

E_{\text{thermal}} = \sum_{i=1}^{N} F_{\text{thermal},i} \cdot A_i

$Fthermal,iF_{\text{thermal},i}$ : Stochastic thermal fluctuation force at point $i$ .

Hydrodynamic Forces

E_{\text{hydro}} = \sum_{i=1}^{N} P_{\text{flow},i} \cdot \mathbf{n}_i \cdot A_i

$Pflow,iP_{\text{flow},i}$ : Pressure or shear force from fluid at point $i$ .
$ni\mathbf{n}_i$ : Normal vector at point $i$ .

3. Constraints

3.1 Volume Constraint

Ensure the RBC volume matches physiological values:

V(C_{L,0}) = \sum_{i=1}^{N} V_i = V_{\text{target}} \quad \text{(e.g., } V_{\text{target}} \approx 90 \, \mu\text{m}^3 \text{)}

3.2 Surface Area Constraint

Ensure the surface area remains constant:

A(C_{L,0}) = \sum_{i=1}^{N} A_i = A_{\text{target}} \quad \text{(e.g., } A_{\text{target}} \approx 145 \, \mu\text{m}^2 \text{)}

3.3 Area Difference Constraint

Account for the area difference between the inner and outer leaflets:

\Delta A(C_{L,0}) = \sum_{i=1}^{N} 2H_i \cdot A_i = \Delta A_{\text{target}} \quad \text{(e.g., } \Delta A_{\text{target}} \approx 0.6 - 1.0 \text{)}

3.4 Reduced Volume Constraint

Maintain a reduced volume less than 1 for a non-spherical shape:

\frac{V(C_{L,0})}{\left( \frac{4}{3} \pi \left( \frac{A(C_{L,0})}{4 \pi} \right)^{3/2} \right)} \quad \text{and } v \in [0.6, 0.7]

3.5 Symmetry Constraints

Impose symmetry on the coefficients:

C_{L,K} = 0 \quad \text{for all } K \neq 0

Python Implementation Using CasADi

We implement the mathematical model using Python and the CasADi library, which allows for symbolic differentiation and optimization. Below, each code block is mapped to the corresponding mathematical expression.

1. Import Necessary Libraries

import casadi as ca  # CasADi for symbolic computations
import numpy as np   # NumPy for numerical operations
import matplotlib.pyplot as plt  # For plotting
from mpl_toolkits.mplot3d import Axes3D  # For 3D plotting

2. Define Constants and Parameters

Set the physical constants, target values, and discretization parameters as per physiological data.

# Physical constants
k_b = 2e-19    # Bending modulus (J)
K_s = 5e-2     # Shear modulus (N/m), may be simplified
gamma = 1e-9   # Coupling constant
C_0 = 0        # Spontaneous curvature (assumed zero)

# Target values for constraints
V_target = 90e-18      # Target volume (m³)
A_target = 145e-12     # Target surface area (m²)
Delta_A_target = 0.85  # Target area difference
v_target = 0.7         # Target reduced volume

# Spherical harmonics parameters
L_max = 6    # Maximum degree for spherical harmonics

# Discretization parameters
N_theta = 200    # Number of discretization points for theta
theta_values = np.linspace(0, np.pi, N_theta)  # θ grid points
dtheta = theta_values[1] - theta_values[0]

3. Set Up the Optimization Problem

Decision Variables:

Mathematical Expression:

C_{L,0} \quad \text{for } L = 0 \text{ to } L_{\text{max}}

Code Implementation:

# Create an optimization problem
opt = ca.Opti()

# Define Theta as a parameter
theta_mx = opt.parameter(N_theta)
opt.set_value(theta_mx, theta_values)
cos_theta = ca.cos(theta_mx)
sin_theta = ca.sin(theta_mx)

# Define Decision Variables
C_L0_var = opt.variable(L_max + 1)  # Spherical harmonics coefficients

4. Compute Legendre Polynomials

Compute $PL(cos⁡θ)P_L(\cos\theta)$ for $L = 0$ to $LmaxL_{\text{max}}$ .

Mathematical Expression:

P_L(\cos\theta)

Code Implementation:

# Function to compute Legendre polynomials iteratively
def legendre_P_iterative(L_max, x):
    P = []
    P.append(ca.MX.ones(x.shape))  # P0 = 1
    if L_max >= 1:
        P.append(x)  # P1 = x
        for L in range(2, L_max + 1):
            P_L = ((2 * L - 1) * x * P[L - 1] - (L - 1) * P[L - 2]) / L
            P.append(P_L)
    return P

# Compute Legendre Polynomials
P_L_list = legendre_P_iterative(L_max, cos_theta)

5. Shape Representation $r(θ)r(\theta)$

Mathematical Expression:

r(\theta) = \sum_{L=0}^{L_{\text{max}}} C_{L,0} \, P_L(\cos\theta)

Code Implementation:

# Compute r(θ) symbolically
r_theta_var = C_L0_var[0] * P_L_list[0]
for L in range(1, L_max + 1):
    r_theta_var += C_L0_var[L] * P_L_list[L]

# Compute mean curvature H(θ)
numerator_var = (
    dr_dtheta_var**2
    - r_theta_var * d2r_dtheta2_var
    + r_theta_var * dr_dtheta_var * ca.tan(theta_mx)
)
denominator_var = (dr_dtheta_var**2 + r_theta_var**2) ** (1.5)
H_i_var = numerator_var / denominator_var

7. Define Energy Components

Bending Energy

Mathematical Expression:

E_{\text{bending}} = \sum_{i=1}^{N} \frac{k_b}{2} \left( 2H_i - C_0 \right)^2 A_i

Code Implementation:

# Bending Energy
E_bending_var = (k_b / 2.0) * ca.sum1((2.0 * H_i_var - C_0) ** 2 * dA_var)

Shear Energy

Assuming negligible shear deformation:

Mathematical Expression:

E_{\text{shear}} = 0

Code Implementation:

# Shear Energy (simplified)
E_shear_var = 0

Coupling Energy

Mathematical Expression:

E_{\text{coupling}} = \sum_{i=1}^{N} \gamma (H_i - C_0) \cdot A_i

Code Implementation:

# Coupling Energy
E_coupling_var = gamma * ca.sum1((H_i_var - C_0) * dA_var)

Total Energy

Mathematical Expression:

E_{\text{total}} = E_{\text{bending}} + E_{\text{shear}} + E_{\text{coupling}}

Code Implementation:

# Total Energy
E_total_var = E_bending_var + E_shear_var + E_coupling_var

8. Define Constraints

Volume Constraint

Mathematical Expression:

\sum_{i=1}^{N} V_i = V_{\text{target}}

Where:

V_i = 2\pi r_i^2 \sin\theta_i \Delta\theta

Code Implementation:

# Volume Constraint
V_var = 2 * np.pi * ca.sum1((r_theta_var**2) * sin_theta * dtheta)
opt.subject_to(V_var == V_target)

Surface Area Constraint

Mathematical Expression:

\sum_{i=1}^{N} A_i = A_{\text{target}}

Code Implementation:

# Surface Area Constraint
A_var = 2.0 * ca.pi * ca.sum1(r_theta_var * sin_theta * ca.sqrt(E_var) * dtheta)
opt.subject_to(A_var == A_target)

Reduced Volume Constraint

Mathematical Expression:

\frac{V}{\left( \frac{4}{3} \pi R_0^3 \right)} = v_{\text{target}}, \quad R_0 = \sqrt{\frac{A}{4\pi}}

Code Implementation:

# Reduced Volume Constraint
v_var = V_var / ((4.0 / 3.0) * ca.pi * (A_var / (4.0 * ca.pi)) ** 1.5)
opt.subject_to(v_var == v_target)

Area Difference Constraint

Mathematical Expression:

\Delta A = \sum_{i=1}^{N} 2H_i A_i = \Delta A_{\text{target}}

Code Implementation:

# Area Difference Constraint
Delta_A_var = ca.sum1(2.0 * H_i_var * dA_var)
opt.subject_to(Delta_A_var == Delta_A_target)

Symmetry Constraints

Only coefficients $C_{L,0}$ are considered, ensuring axial symmetry.

Code Implementation:

Implicitly enforced by only defining $C_{L,0}$ variables (no $C_{L,K}$ for $\neq 0$ ).

Positivity Constraint for $r(θ)r(\theta)$

Ensure physical validity of the shape.

Mathematical Expression:

r(\theta_i) > 0 \quad \forall i

Code Implementation:

# Ensure positive radius
opt.subject_to(r_theta_var > 0)

Biconcavity Enforcement

To ensure a biconcave shape, enforce $drdθ<0\frac{dr}{d\theta} < 0$ around $θ=π2\theta = \frac{\pi}{2}$ .

Mathematical Expression:

\frac{dr}{d\theta} \bigg|_{\theta = \frac{\pi}{2}} < 0

Code Implementation:

# Enforce biconcavity around θ = π/2
delta_theta = np.pi / 20        # Range around θ = π/2
theta_mid = np.pi / 2
theta_lower = theta_mid - delta_theta
theta_upper = theta_mid + delta_theta

# Indices within specified θ range
theta_mid_indices = np.where(
    (theta_values >= theta_lower) & (theta_values <= theta_upper)
)[0]

epsilon = 1e-6  # Small positive value

# Enforce dr/dθ <= -epsilon in the specified θ range
for idx in theta_mid_indices:
    opt.subject_to(dr_dtheta_var[idx] <= -epsilon)

9. Set Objective and Solve

Code Implementation:

# Set the objective function
opt.minimize(E_total_var)

# Set initial guesses for coefficients
C_L0_initial = np.zeros(L_max + 1)
C_L0_initial[0] = bound_scale  # Approximate radius
opt.set_initial(C_L0_var, C_L0_initial)

# Set solver options
p_opts = {"expand": True}
s_opts = {"max_iter": 10000, "tol": 1e-10, "acceptable_tol": 1e-8}

# Specify solver
opt.solver("ipopt", p_opts, s_opts)

# Solve the optimization problem
try:
    sol = opt.solve()
except RuntimeError as e:
    print("Solver failed:", e)
    sol = opt.debug   # Retrieve solution even if solver fails

# Retrieve the optimal coefficients
C_L0_opt = sol.value(C_L0_var)

10. Post-processing and Visualization

Compute Optimized Shape

Mathematical Expression:

r_{\text{opt}}(\theta) = \sum_{L=0}^{L_{\text{max}}} C_{L,0}^{\text{opt}} \, P_L(\cos\theta)

Code Implementation:

# Compute the optimal shape r(θ) using the optimal coefficients
r_theta_opt = C_L0_opt[0] * P_L_list[0]
for L in range(1, L_max + 1):
    r_theta_opt += C_L0_opt[L] * P_L_list[L]
r_theta_opt = sol.value(r_theta_opt)

Generate 3D Coordinates

Convert $ropt(θ)r_{\text{opt}}(\theta)$ to Cartesian coordinates for visualization.

Mathematical Expressions:

\begin{align*} x &= r(\theta) \sin\theta \cos\phi \ y &= r(\theta) \sin\theta \sin\phi \ z &= r(\theta) \cos\theta \end{align*}

Code Implementation:

# Additional discretization for φ
N_phi = 200
phi_values = np.linspace(0, 2 * np.pi, N_phi)
theta_grid, phi_grid = np.meshgrid(theta_values, phi_values)

# Replicate r(θ) across φ
r_grid = np.tile(r_theta_opt, (N_phi, 1))

# Compute the 3D coordinates (x, y, z)
x_coords = r_grid * np.sin(theta_grid) * np.cos(phi_grid)
y_coords = r_grid * np.sin(theta_grid) * np.sin(phi_grid)
z_coords = r_grid * np.cos(theta_grid)

Plot the Optimized Shape

Code Implementation:

# Plot the 3D shape
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(111, projection="3d")
ax.plot_surface(
    x_coords,
    y_coords,
    z_coords,
    rstride=4,
    cstride=4,
    color='red',
    edgecolor="none",
    alpha=0.8
)
ax.set_xlabel("x (m)")
ax.set_ylabel("y (m)")
ax.set_zlabel("z (m)")
ax.set_title("Optimized Red Blood Cell Shape (3D)")

# Ensure equal scaling on all axes
ax.set_box_aspect([np.ptp(x_coords), np.ptp(y_coords), np.ptp(z_coords)])

plt.tight_layout()
plt.show()

11. Validate Results

Ensure that the optimized shape satisfies the constraints.

Code Implementation:

# Compute final volume and area
V_final = sol.value(V_var)
A_final = sol.value(A_var)
v_final = sol.value(v_var)
Delta_A_final = sol.value(Delta_A_var)

print(f"Final Volume: {V_final:.3e} m³")
print(f"Target Volume: {V_target:.3e} m³")
print(f"Final Surface Area: {A_final:.3e} m²")
print(f"Target Surface Area: {A_target:.3e} m²")
print(f"Final Reduced Volume: {v_final:.3f}")
print(f"Target Reduced Volume: {v_target:.3f}")
print(f"Final Delta A: {Delta_A_final:.3e}")
print(f"Target Delta A: {Delta_A_target:.3e}")

Discussion

Mapping Mathematical Expressions to Code

In the above implementation, we meticulously map each mathematical expression from the model to corresponding code blocks. This explicit correspondence ensures that the code accurately reflects the theoretical formulations and that each component is correctly implemented.

Decision Variables: The spherical harmonics coefficients $C_{L,0}$ are represented by C_L0_var.
Shape Representation: The radial function $r(θ)r(\theta)$ is constructed using the Legendre polynomials and the coefficients.
Geometric Quantities: Derivatives, mean curvature, and area elements are computed symbolically, matching the mathematical definitions.
Energy Components: The bending and coupling energies are calculated as per their mathematical expressions, and summed to form the total energy.
Constraints: Volume, surface area, reduced volume, area difference, and symmetry constraints are enforced using opt.subject_to.
Optimization: The total energy is minimized using IPOPT via CasADi.
Visualization: The optimized shape is visualized in 3D, providing a graphical representation of the results.
Validation: Final values are compared to target constraints to verify the optimization's success.

Simplifications and Assumptions

Shear Energy: Set to zero, assuming negligible shear deformation for simplicity.
Viscous, Thermal, and Hydrodynamic Energies: Omitted in this static optimization model.
Spontaneous Curvature ( $C_0$ ): Assumed zero, but can be adjusted if needed.

Advantages of Using CasADi

Symbolic Differentiation: Enables precise computation of derivatives required for curvature and energy calculations.
Optimization Framework: Provides tools to define variables, objectives, and constraints in an intuitive manner.
Solver Integration: Seamlessly integrates with IPOPT and other solvers for efficient numerical optimization.

Conclusion

This article presented a detailed mathematical model for optimizing the shape of a red blood cell using spherical harmonics parameterization under axial symmetry. By mapping each mathematical expression to corresponding code blocks, we provided a transparent and precise implementation.

The model minimizes the membrane's total energy while satisfying physiological constraints, resulting in an optimized biconcave shape consistent with actual RBCs. The approach demonstrates the power of combining mathematical modeling with computational tools like CasADi to explore and analyze complex biological systems.

Applying Monte Carlos Model Modification to Handle Infeasibility Internally in Stochastic Programming

Thana B. — Sun, 06 Oct 2024 10:07:27 +0000

Introduction

When dealing with stochastic programming and Monte Carlo simulation-based scenario generation, you may encounter infeasible scenarios that complicate the optimization process. One effective strategy to address this issue is to modify the optimization model to handle infeasibility internally. This approach eliminates the need to adjust scenario probabilities on the fly and allows the optimization to proceed smoothly.

In this guide, I'll teach you how to apply this strategy, focusing on two main methods:

Introducing Penalty Terms or Slack Variables: Adjusting the model to allow for constraint violations at a cost.
Using Chance Constraints: Replacing hard constraints with probabilistic ones that can tolerate a small acceptable level of constraint violation.

We'll go through each method in detail, explaining the concepts, providing mathematical formulations, and illustrating with examples in the context of the Stochastic Knapsack Problem.

1. Introducing Penalty Terms or Slack Variables

Concept Overview

Goal: Allow the model to handle infeasible scenarios internally by permitting constraint violations but penalizing them in the objective function.

How It Works:

Slack Variables: Add variables to represent the extent of constraint violations (slack or surplus) and include associated penalties in the objective function.
Penalty Terms: Directly incorporate penalties for constraint violations into the objective function.

This approach ensures that the optimization process considers the cost of violating constraints and seeks solutions that minimize both the original objective and the penalties.

Applying Penalty Terms to the Stochastic Knapsack Problem

Original Stochastic Knapsack Model

Let's recall the original formulation:

Decision Variables:
- $xi∈0,1x_i \in {0, 1}$ : Whether to include item $i$ in the knapsack.
Parameters:
- $v_i^s$ : Value of item $i$ in scenario $s$ .
- $w_i^s$ : Weight of item $i$ in scenario $s$ .
- $W$ : Knapsack capacity.
- $p_s$ : Probability of scenario $s$ .
Objective Function:

\text{Maximize } Z = \sum_{s=1}^{n} p_s \left( \sum_{i=1}^{n_{\text{items}}} v_i^s x_i \right)

Constraints:

\sum_{i=1}^{n_{\text{items}}} w_i^s x_i \leq W, \quad \forall s = 1, 2, \dots, n

Modified Model with Penalty Terms

We introduce penalties for exceeding the knapsack capacity in any scenario.

Additional Variables:
- $δs≥0\delta_s \geq 0$ : Amount by which the capacity constraint is violated in scenario $s$ .
Modified Constraints:

\sum_{i=1}^{n_{\text{items}}} w_i^s x_i \leq W + \delta_s, \quad \forall s = 1, 2, \dots, n

Penalizing Constraint Violations:
- We add a penalty term to the objective function that includes the sum of $δs\delta_s$ multiplied by a penalty coefficient.
Penalty Coefficient:
- Let $M$ be a large positive number representing the penalty per unit of constraint violation.
Modified Objective Function:

\text{Maximize } Z = \sum_{s=1}^{n} p_s \left( \sum_{i=1}^{n_{\text{items}}} v_i^s x_i - M \delta_s \right)

Model Interpretation

The model seeks to maximize the expected total value of selected items while minimizing the expected penalty from constraint violations.
The penalty $M$ should be set high enough to discourage constraint violations unless absolutely necessary.

Implementation Steps

Set Penalty Coefficient $M$ :
- Choose an appropriate value for $M$ .
- $M$ should be large enough to make constraint violations undesirable but not so large that it causes numerical instability.
Modify Constraints:
- Add slack variables $δs\delta_s$ to the capacity constraints for each scenario.
Add Penalty Terms to Objective Function:
- Subtract the weighted penalty $\delta_s$ in the objective function.
Adjust Variable Domains:
- Ensure that $δs≥0\delta_s \geq 0$ for all scenarios.
Solve the Modified Optimization Problem:
- Use an appropriate optimization solver capable of handling the modified model.

Example

Let's work through a simplified example with numerical values.

Scenario Data:

Items: 2 items.
Scenarios: 3 scenarios.
Knapsack Capacity: $W = 10$ units.

Scenario $s$	$p_s$	$v_1^s$	$w_1^s$	$v_2^s$	$w_2^s$
1	0.333	7	6	5	5
2	0.333	8	7	6	6
3	0.333	6	8	4	7

Original Constraints:

Scenario 3 is infeasible because the minimum possible total weight ( $w_1^3 + w_2^3 = 8 + 7 = 15$ ) exceeds $W = 10$ .

Modified Model with Penalty Terms:

Introduce $δs≥0\delta_s \geq 0$ for $s = 1, 2, 3$ .
Modified Constraints:

x_1 + 5 x_2 \leq 10 + \delta_1 \

x_1 + 6 x_2 \leq 10 + \delta_2 \

x_1 + 7 x_2 \leq 10 + \delta_3 \

Modified Objective Function:

\text{Maximize } Z = \sum_{s=1}^{3} p_s \left( v_1^s x_1 + v_2^s x_2 - M \delta_s \right)

Let's set $M = 100$ .

Decision Variables:

$xi∈0,1x_i \in {0, 1}$ for $i = 1, 2$ .
$δs≥0\delta_s \geq 0$ for $s = 1, 2, 3$ .

Solving the Model:

Use an optimization solver to find the values of $x1,x2,δ1,δ2,δ3x_1, x_2, \delta_1, \delta_2, \delta_3$ that maximize $Z$ while satisfying the constraints.

Interpretation of Results:

The solver will likely select items that maximize total value while minimizing or avoiding penalties.
In Scenario 3, some penalty may be incurred if constraint violation is unavoidable.

Considerations When Using Penalty Terms

Choosing Penalty Coefficient $M$ :
- Too Small: The model may accept constraint violations too readily.
- Too Large: Numerical issues may arise, or the model may become too conservative.
Penalties Reflect Real Costs:
- If possible, set $M$ to reflect the actual cost or impact of violating the constraint.
Interpretability:
- Solutions with constraint violations should be acceptable within the context of the problem.

2. Using Chance Constraints

Concept Overview

Goal: Allow constraints to be violated in a controlled manner by requiring them to be satisfied with a specified probability level.

How It Works:

Chance Constraints: Constraints are reformulated to require that they are met in all but a small fraction $ϵ\epsilon$ of scenarios.
Probability Level: Decision-makers set $ϵ\epsilon$ , the maximum acceptable probability of constraint violation.

Applying Chance Constraints to the Stochastic Knapsack Problem

Model Formulation

Decision Variables:

$xi∈0,1x_i \in {0, 1}$ : Decision variables for selecting items.

Parameters:

As before: $v_i^s, w_i^s, W, p_s$ .

Chance Constraint:

Require that the capacity constraint is satisfied with at least probability $\epsilon$ :

P\left( \sum_{i=1}^{n_{\text{items}}} w_i^s x_i \leq W \right) \geq 1 - \epsilon

This means that the weight constraint can be violated in at most $ϵ×100%\epsilon \times 100\%$ of the scenarios.

Reformulating the Chance Constraint

Since we have a discrete set of scenarios, the chance constraint can be expressed as:

\sum_{s \in S_{\text{feasible}}} p_s \geq 1 - \epsilon

Where $SfeasibleS_{\text{feasible}}$ is the set of scenarios where the constraint is satisfied.

Integer Programming Formulation

Introduce binary variables $y_s$ to indicate whether the constraint is satisfied in scenario $s$ :

$y_s = 1$ if the constraint is satisfied in scenario $s$ .
$y_s = 0$ if the constraint is violated in scenario $s$ .

Modified Constraints:

Capacity Constraint with Indicator Variables:

$\sum_{i=1}^{n_{\text{items}}} w_i^s x_i \leq W + M (1 - y_s), \quad \forall s = 1, 2, \dots, n$
- $M$ is a large constant that effectively deactivates the constraint when $y_s = 0$ .
Chance Constraint:
$\sum_{s=1}^{n} p_s y_s \geq 1 - \epsilon$

Objective Function:
Maximize expected value as before:
$\text{Maximize } Z = \sum_{s=1}^{n} p_s \left( \sum_{i=1}^{n_{\text{items}}} v_i^s x_i \right)$

Variable Domains:
$xi∈0,1x_i \in {0, 1}$ $ys∈0,1y_s \in {0, 1}$

Explanation

For each scenario $s$ , if $y_s = 1$ , the capacity constraint must be satisfied.
If $y_s = 0$ , the capacity constraint can be violated.
The chance constraint ensures that the probability of scenarios where $y_s = 1$ (i.e., the constraint is satisfied) is at least $\epsilon$ .

Example

Using the same scenario data as before:

Scenario $s$	$p_s$	$v_1^s$	$w_1^s$	$v_2^s$	$w_2^s$
1	0.333	7	6	5	5
2	0.333	8	7	6	6
3	0.333	6	8	4	7

Suppose we can accept a 10% chance of violating the capacity constraint ( $ϵ=0.10\epsilon = 0.10$ ).

Variables:

x1,x2∈0,1x_1, x_2 \in {0, 1}

y1,y2,y3∈0,1y_1, y_2, y_3 \in {0, 1}

Constraints:

Capacity Constraints:

$x_1 + 5 x_2 \leq 10 + M (1 - y_1)$ $x_1 + 6 x_2 \leq 10 + M (1 - y_2)$ $x_1 + 7 x_2 \leq 10 + M (1 - y_3)$
Chance Constraint:

$y_1 + 0.333 y_2 + 0.333 y_3 \geq 0.90$
Variable Domains:

$xi∈0,1x_i \in {0, 1}$ $ys∈0,1y_s \in {0, 1}$

Objective Function:

Same as before, without penalties:

\text{Maximize } Z = \sum_{s=1}^{3} 0.333 \left( v_1^s x_1 + v_2^s x_2 \right)

M Value:

Choose a sufficiently large $M$ to deactivate the constraint when $y_s = 0$ .

Solving the Model:

Use an optimization solver that can handle mixed-integer programming problems.

Interpretation of Results:

The model will select $x_1$ and $x_2$ to maximize expected value while ensuring that the capacity constraint is satisfied in at least 90% of the total probability.
The capacity constraint can be violated in scenarios totaling up to 10% probability.

Considerations When Using Chance Constraints

Setting $ϵ\epsilon$ :
- Decide on the acceptable risk level for constraint violations based on practical considerations.
Computational Complexity:
- The model becomes a mixed-integer program, which may increase computational requirements.
Feasibility:
- Ensure that the model remains feasible under the specified chance constraint.

Comparing the Two Methods

Penalty Terms vs. Chance Constraints

Penalty Terms:
- Flexibility: Violations are penalized but allowed if they lead to a better overall objective value.
- Penalty Calibration: Requires careful selection of penalty coefficients.
- Objective Balance: Directly balances maximizing the objective and minimizing constraint violations.
Chance Constraints:
- Risk Control: Provides explicit control over the acceptable probability of constraint violations.
- Binary Decisions: Uses binary variables to enforce or relax constraints in specific scenarios.
- Complexity: May be more complex due to the introduction of additional binary variables.

Choosing Between Methods

Penalty Terms:
- Suitable when the cost of constraint violations can be quantified or estimated.
- Appropriate if occasional violations are acceptable and their impact is known.
Chance Constraints:
- Ideal when there is a need to limit the probability of constraint violations explicitly.
- Useful in risk management contexts where constraints must be satisfied with high confidence.

Implementation Tips

General Advice

Understanding the Context:
- Consider the practical implications of constraint violations in your specific problem.
- Understand whether penalties or controlled violation probabilities align better with decision-maker preferences.
Solver Selection:
- Use optimization solvers capable of handling mixed-integer programming (e.g., CPLEX, Gurobi).
- Ensure the solver can handle the size and complexity of your modified model.
Parameter Tuning:
- Penalty Coefficient $M$ :
- Test different values to observe their impact on the solution.
- Avoid excessively large values that may cause numerical issues.
- Chance Constraint Level $ϵ\epsilon$ :
- Experiment with different risk levels to find one that balances feasibility and optimality.

Testing and Validation

Small-Scale Testing:
- Start with a simplified version of your problem to test the modified model.
Scenario Analysis:
- Analyze how the solutions change with different penalty coefficients or chance constraint levels.
Sensitivity Analysis:
- Assess the sensitivity of the optimal solution to changes in parameters.
Feasibility Checks:
- Verify that the solutions satisfy constraints in the required proportion of scenarios.

Documentation

Model Assumptions:
- Clearly document assumptions made in modifying the model.
Parameter Choices:
- Record how penalty coefficients and chance constraint levels were selected.
Solution Interpretation:
- Provide interpretation of the results, especially regarding constraint violations.

Conclusion

By modifying your optimization model to handle infeasibility internally, you avoid the need to adjust scenario probabilities during optimization. Both introducing penalty terms or slack variables and using chance constraints are effective methods to achieve this.

Penalty Terms or Slack Variables enable the model to penalize constraint violations, balancing the objective value with the cost of violations.
Chance Constraints allow for a controlled level of constraint violations by specifying the acceptable probability of such events.

Implementing these methods involves adjusting your mathematical model, carefully selecting parameters, and using appropriate optimization solvers. By following the steps outlined above, you can apply these techniques to your stochastic programming problems and effectively manage infeasibility within your optimization model.

Monte Carlo Simulation-Based Scenario Generation in Stochastic Programming: Addressing Uncertainty in the Knapsack Problem

Thana B. — Sun, 06 Oct 2024 09:55:21 +0000

Introduction

In the realm of optimization and decision-making, the Knapsack Problem stands as a quintessential example, highlighting the challenges of selecting the most valuable combination of items under a weight constraint. Traditionally, this problem assumes known and fixed weights and values for the items. However, real-world scenarios often involve uncertainties due to factors like market fluctuations, measurement errors, or logistical variances.

To effectively handle such uncertainties, Stochastic Programming integrates randomness directly into optimization models. A practical technique within this framework is Monte Carlo Simulation-Based Scenario Generation, which creates numerous possible scenarios by randomly sampling from historical data. Each scenario represents a potential state of the world, capturing the variability in item weights and values.

This comprehensive article delves into the Monte Carlo simulation technique in the context of stochastic programming, emphasizing its value in solving real-world problems. We will explore foundational concepts, explain the methodology step-by-step, and provide detailed examples using the Knapsack Problem. Additionally, we address the challenge of infeasible scenarios that may arise during scenario generation and present strategies to handle them effectively within the optimization model.

The Value of Monte Carlo Simulation in Real-World Problems

Embracing Uncertainty in Decision Making

In practical applications, decision-makers often grapple with uncertainty in key parameters influencing outcomes:

Supply Chain Management: Demand and supply levels fluctuate due to seasonal changes or market trends.
Finance: Investment returns vary with market volatility.
Engineering: Material properties may deviate from specifications due to manufacturing tolerances.

Ignoring these uncertainties can lead to decisions that are suboptimal or even detrimental when actual conditions deviate from expectations.

Limitations of Deterministic Models

Deterministic models assume precise knowledge of all parameters, which is rarely the case in reality. These models may:

Overlook Risks: By not accounting for variability, they ignore potential adverse outcomes.
Yield Fragile Solutions: Solutions that perform well under assumed conditions but poorly when parameters vary.
Fail to Capture Realistic Scenarios: Missing the complexity and randomness inherent in real-world situations.

Advantages of Stochastic Programming with Monte Carlo Simulation

Stochastic Programming incorporates uncertainty by considering multiple possible scenarios, each representing a different set of parameter values. Monte Carlo Simulation enhances this approach by:

Reflecting Realistic Variability: Generating scenarios based on historical data captures the true nature of uncertainties.
Providing Probabilistic Insights: Each scenario has an associated probability, allowing for expected value calculations and risk assessments.
Enabling Robust Decision-Making: Solutions are evaluated across numerous scenarios, ensuring they perform well under varying conditions.

Understanding Monte Carlo Simulation in This Context

Concept of Monte Carlo Simulation

Monte Carlo Simulation is a computational method that uses random sampling to model and analyze complex systems influenced by uncertainty. Key features include:

Random Sampling from Historical Data: Drawing values from historical records to reflect actual observed variations.
Scenario Generation: Creating numerous possible combinations of parameters to simulate different future states.
Probability Assignment: Assuming each randomly generated scenario has an equal chance of occurring when sampling is unbiased.

Application to Stochastic Programming

In stochastic programming:

Data Collection: Historical data for uncertain parameters (e.g., item weights and values) are gathered.
Random Sampling: Parameters are randomly selected from historical data for each item in each scenario.
Scenario Creation: Combining these random selections to form a complete scenario.
Repetition: Generating a large number of scenarios by repeating the process.
Equal Probability: Assigning an equal probability to each scenario, typically $1n\frac{1}{n}$ where $n$ is the number of scenarios.

The Knapsack Problem

1. Traditional Deterministic Model

Problem Definition

Objective: Maximize the total value of items placed in a knapsack without exceeding its weight capacity.
Assumptions: The weights and values of items are known and fixed.

Mathematical Formulation

Let:

$n$ = number of items
$v_i$ = value of item $i$
$w_i$ = weight of item $i$
$W$ = maximum weight capacity of the knapsack
$xi∈0,1x_i \in {0,1}$ = decision variable (1 if item $i$ is included, 0 otherwise)

Objective Function:

Maximize total value:

\text{Maximize } Z = \sum_{i=1}^{n} v_i x_i

Constraint:

Total weight cannot exceed capacity:

\sum_{i=1}^{n} w_i x_i \leq W

Limitations

Ignores uncertainty in $w_i$ and $v_i$ .
May result in solutions that are infeasible or suboptimal when parameters vary in reality.

2. Stochastic Programming Using Monte Carlo Simulation-Based Scenario Generation

Step-by-Step Guide

Step 1: Collect Historical Data

Gather historical data for each item's weight and value. For example:

Item 1:
- Weights: [12, 13, 11, 12]
- Values: [40, 42, 39, 41]
Item 2:
- Weights: [8, 9, 7, 8]
- Values: [30, 32, 28, 31]

Step 2: Random Sampling from Historical Data

For each scenario:

Item-wise Random Selection:
- Randomly select a weight and value for Item 1 from its historical data.
- Randomly select a weight and value for Item 2 from its historical data.
- Continue this process for all items.
Combine Selections:
- The selected weights and values for all items form one scenario.

Step 3: Repeat to Generate Multiple Scenarios

Number of Scenarios: Decide on a suitable number $n$ (e.g., 10,000).
Process:
- Repeat the random sampling process $n$ times.
- Each iteration generates a new scenario with randomly selected weights and values for all items.

Step 4: Assign Equal Probabilities

Since each scenario is generated through unbiased random sampling, assign a probability $ps=1np_s = \frac{1}{n}$ to each scenario $s$ .

Step 5: Formulate the Stochastic Optimization Model

Decision Variables:

$xi∈0,1x_i \in {0, 1}$ : Whether to include item $i$ in the knapsack.

Objective Function:

Maximize the expected total value across all scenarios:

\text{Maximize } Z = \sum_{s=1}^{n} p_s \left( \sum_{i=1}^{n} v_i^s x_i \right) = \frac{1}{n} \sum_{s=1}^{n} \left( \sum_{i=1}^{n} v_i^s x_i \right)

Constraints:

Weight constraint must be satisfied in each scenario: $∑i=1nwisxi≤W,∀s=1,2,…,n\sum_{i=1}^{n} w_i^s x_i \leq W, \quad \forall s = 1, 2, \dots, n$

Step 6: Solve the Model

Use optimization software capable of handling large-scale stochastic programs (e.g., CPLEX, Gurobi).
Techniques like Sample Average Approximation (SAA) can be employed to approximate the expected value.

Handling Infeasible Scenarios in Monte Carlo Simulation

Understanding Infeasibility

Infeasible Scenario:

A scenario in which the constraints of the optimization problem cannot be satisfied by any feasible solution, including the option of taking no action (e.g., selecting no items).

Causes of Infeasibility:

Unrealistic Parameter Values: Random sampling might generate parameter combinations that are not realistic or physically possible.
Excessive Uncertainty: High variability in the input data can lead to extreme values making the problem infeasible.
Model Misalignment: The model may not correctly capture the relationships and dependencies between variables, leading to contradictions.

Implications of Infeasible Scenarios

Optimization Failure: The presence of infeasible scenarios can cause the optimization problem to have no solution.
Computational Issues: Solvers may not handle infeasibility gracefully, leading to errors or excessive computation times.
Decision-Making Challenges: Without feasible solutions, decision-makers cannot derive actionable insights from the model.

Strategies to Handle Infeasible Scenarios

1. Modifying the Optimization Model to Handle Infeasibility Internally

Instead of adjusting probabilities or filtering scenarios, modify the model to handle infeasible scenarios internally. This can be done by:

Introducing Penalty Terms or Slack Variables: Allowing constraint violations at a cost.
Using Chance Constraints: Replacing hard constraints with probabilistic ones that tolerate a small acceptable level of constraint violation.

Applying Model Modification to the Knapsack Problem

Introducing Penalty Terms or Slack Variables

Concept Overview

Goal: Allow the model to handle infeasible scenarios internally by permitting constraint violations but penalizing them in the objective function.

How It Works:

Slack Variables: Add variables to represent the extent of constraint violations (slack or surplus) and include associated penalties in the objective function.
Penalty Terms: Directly incorporate penalties for constraint violations into the objective function.

This approach ensures that the optimization process considers the cost of violating constraints and seeks solutions that balance maximizing the objective with minimizing penalties.

Applying Penalty Terms to the Stochastic Knapsack Problem

Original Stochastic Knapsack Model

Decision Variables:
- $xi∈0,1x_i \in {0, 1}$ : Whether to include item $i$ in the knapsack.
Parameters:
- $v_i^s$ : Value of item $i$ in scenario $s$ .
- $w_i^s$ : Weight of item $i$ in scenario $s$ .
- $W$ : Knapsack capacity.
- $p_s$ : Probability of scenario $s$ .
Objective Function:
$\text{Maximize } Z = \sum_{s=1}^{n} p_s \left( \sum_{i=1}^{n} v_i^s x_i \right)$
Constraints:
$\sum_{i=1}^{n} w_i^s x_i \leq W, \quad \forall s = 1, 2, \dots, n$

Modified Model with Penalty Terms

We introduce penalties for exceeding the knapsack capacity in any scenario.

Additional Variables:
- $δs≥0\delta_s \geq 0$ : Amount by which the capacity constraint is violated in scenario $s$ .
Modified Constraints:
$\sum_{i=1}^{n} w_i^s x_i \leq W + \delta_s, \quad \forall s = 1, 2, \dots, n$
Penalizing Constraint Violations:
- Add a penalty term to the objective function, including the sum of $δs\delta_s$ multiplied by a penalty coefficient $M$ .
Modified Objective Function:
$\text{Maximize } Z = \sum_{s=1}^{n} p_s \left( \sum_{i=1}^{n} v_i^s x_i - M \delta_s \right)$

Model Interpretation

The model seeks to maximize the expected total value of selected items while minimizing the expected penalty from constraint violations.
The penalty $M$ should be set high enough to discourage constraint violations unless absolutely necessary.

Implementation Steps

Set Penalty Coefficient $M$ :
- Choose an appropriate value for $M$ that reflects the cost of violating the constraint.
- $M$ should be large enough to make constraint violations undesirable but not so large as to cause numerical instability.
  - Choosing Penalty Coefficient $M$ :
    - Too Small: The model may accept constraint violations too readily.
    - Too Large: Numerical issues may arise, or the model may become too conservative.
Modify Constraints:
- Add slack variables $δs\delta_s$ to the capacity constraints for each scenario.
Add Penalty Terms to Objective Function:
- Incorporate $\delta_s$ in the objective function.
Adjust Variable Domains:
- Ensure that $δs≥0\delta_s \geq 0$ for all scenarios.
Solve the Modified Optimization Problem:
- Use an optimization solver capable of handling the modified model.

Example

Let's work through a simplified example with numerical values.

Scenario Data:

Items: 2 items.
Scenarios: 3 scenarios.
Knapsack Capacity: $W = 10$ units.

Scenario $s$	$p_s$	$v_1^s$	$w_1^s$	$v_2^s$	$w_2^s$
1	0.333	7	6	5	5
2	0.333	8	7	6	6
3	0.333	6	8	4	7

Original Constraints:

Scenario 3 is infeasible because the total weight of any combination exceeds $W = 10$ .

Modified Model with Penalty Terms:

Introduce $δs≥0\delta_s \geq 0$ for $s = 1, 2, 3$ .
Modified Constraints:

$\begin{cases} 6 x_1 + 5 x_2 \leq 10 + \delta_1 \ 7 x_1 + 6 x_2 \leq 10 + \delta_2 \ 8 x_1 + 7 x_2 \leq 10 + \delta_3 \ \end{cases}$
Modified Objective Function:

$\text{Maximize } Z = \sum_{s=1}^{3} p_s \left( v_1^s x_1 + v_2^s x_2 - M \delta_s \right)$
Let's set $M = 100$ .

Decision Variables:

$xi∈0,1x_i \in {0, 1}$ for $i = 1, 2$ .
$δs≥0\delta_s \geq 0$ for $s = 1, 2, 3$ .

Solving the Model:

Use an optimization solver to find the values of $x1,x2,δ1,δ2,δ3x_1, x_2, \delta_1, \delta_2, \delta_3$ that maximize $Z$ while satisfying the constraints.

Interpretation of Results:

The solver will seek to maximize the expected value of selected items while minimizing penalties from any constraint violations.
In Scenario 3, some penalty may be incurred if constraint violation is unavoidable.

Considerations When Using Penalty Terms

Choosing Penalty Coefficient $M$ :
- Too Small: The model may accept constraint violations too readily.
- Too Large: Numerical issues may arise, or the model may become overly conservative.
Penalties Reflect Real Costs:
- If possible, set $M$ to reflect the actual cost or impact of violating the constraint.
Interpretability:
- Solutions with constraint violations should be acceptable within the context of the problem.

Conclusion

Monte Carlo Simulation-Based Scenario Generation offers a practical and effective means of handling uncertainty in optimization problems like the Knapsack Problem. However, infeasible scenarios can arise, complicating the optimization process. By modifying the optimization model to handle infeasibility internally—through penalty terms or chance constraints—we can:

Maintain All Scenarios: Include both feasible and infeasible scenarios in the analysis.
Avoid On-the-Fly Probability Adjustments: Eliminate the need to adjust scenario probabilities during optimization.
Balance Objective Maximization and Feasibility: Find solutions that optimize performance while considering the cost of constraint violations.

Key Takeaways:

Data-Driven Decision-Making: Using historical data and Monte Carlo simulation captures real-world variability.
Model Flexibility: Modifying the model to handle infeasibility allows for robust and adaptable optimization.
Improved Decision Quality: Considering penalties or acceptable levels of constraint violations leads to more practical and implementable solutions.

Practical Implications:

Applicability: The approach is valuable in various fields, including supply chain management, finance, and engineering, where uncertainty and infeasibility are common challenges.
Scalability: While model modifications increase complexity, they enable the handling of larger and more realistic problems.
Strategic Advantage: Organizations can develop strategies that are resilient to uncertainties and feasible under a wider range of conditions.

This comprehensive article has combined foundational concepts with practical strategies to address infeasibility in stochastic programming using Monte Carlo simulation. By understanding and applying these methods, decision-makers can enhance the robustness and reliability of their optimization models in the face of uncertainty.

Goal Programming: Balancing Multiple Objectives with Flexibility

Thana B. — Sun, 06 Oct 2024 09:54:34 +0000

In the field of operations research, decision-makers often face complex problems involving multiple, sometimes conflicting, objectives. Traditional optimization techniques, like linear programming, focus on optimizing a single objective function subject to a set of rigid constraints. However, real-world scenarios rarely adhere to such simplicity. Constraints can be flexible, and goals may need to be balanced rather than strictly optimized. Goal programming emerges as a powerful tool in these situations, enabling decision-makers to accommodate multiple objectives and find feasible solutions even when traditional models indicate infeasibility.

The Importance of Goal Programming

In practice, not all constraints are absolute; many represent soft constraints or goals that we aim to achieve if possible. Traditional linear programming treats constraints as rigid requirements that must be strictly satisfied for a solution to be feasible. If these constraints are conflicting or overly restrictive, the model may become infeasible, offering no solution.

Goal programming addresses this issue by introducing flexibility into the model without altering the original problem instance. It allows decision-makers to prioritize and balance multiple objectives, accommodating the inherent "softness" of real-world constraints. By permitting controlled deviations from certain constraints, goal programming transforms an unsolvable problem into a solvable one, enabling users to fine-tune their models effectively.

Benefits of Goal Programming

Flexibility: Accommodates soft constraints, reflecting real-world conditions where not all requirements are absolute.
Feasibility: Transforms infeasible models into solvable ones by allowing deviations from less critical constraints.
Prioritization: Enables ranking objectives according to their importance, ensuring critical goals are addressed first.
Tunability: Allows users to adjust goal priorities without changing the original problem instance, providing a powerful tool for model refinement.

How to Formulate Goal Constraints

Formulating goal constraints involves converting soft constraints and objectives into a form that allows for deviations, measured by deviational variables. This process introduces flexibility into the model, enabling the accommodation of conflicting objectives and the balancing of multiple goals.

Steps to Formulate Goal Constraints

Identify the Goals: List all desired outcomes and soft constraints.
Write the Requirement as a Regular Constraint: Formulate each goal as a standard constraint (equality or inequality).
Reformulate as a Goal Constraint: Introduce deviational variables to measure underachievement and overachievement.

The general method for formulating goal constraints is summarized in the table below:

Desired Situation	Formulation of Goal Constraint	Contribution to Objective Function
$LHS≤RHS\text{LHS} \leq \text{RHS}$	$LHS+d−−d+=RHS\text{LHS} + d^- - d^+ = \text{RHS}$	Minimize $d^+$
$LHS=RHS\text{LHS} = \text{RHS}$	$LHS+d−−d+=RHS\text{LHS} + d^- - d^+ = \text{RHS}$	Minimize $d^- + d^+$
$LHS≥RHS\text{LHS} \geq \text{RHS}$	$LHS+d−−d+=RHS\text{LHS} + d^- - d^+ = \text{RHS}$	Minimize $d^-$

Establish the Objective Function: Create an objective function that minimizes the weighted sum of the deviational variables, reflecting the priorities of each goal.

Explanation of Deviational Variables

Negative Deviational Variable ( $d^-$ ): Measures the amount by which a goal is underachieved.
Positive Deviational Variable ( $d^+$ ): Measures the amount by which a goal is overachieved.

Note: For any particular goal, only one of the deviational variables ( $d^-$ or $d^+$ ) will be positive in the optimal solution, as a goal cannot be both underachieved and overachieved simultaneously.

Example: Allocating Diamonds to Stores

Problem Statement

The owner of a chain of jewelry stores has to decide how to distribute parts of a new shipment of diamonds to five stores in a region. The first three stores are located in shopping malls. The following conditions have to be observed:

Absolutely Necessary (Hard Constraints):

(a) Allocate between 1,000 and 1,200 carats in total to the five stores.
(b) Store 5 must receive at least 300 carats of diamonds.

Desired Properties of the Allocation (Soft Constraints and Objectives):

(c) The stores in the malls (Stores 1, 2, and 3) should receive at least 80% of all the diamonds, if possible.
(d) The allocations to the stores in the malls should be equal to each other, if possible.
(e) The probabilities of theft in the stores have been estimated to be 0.1%, 0.1%, 9%, 2%, and 3% for Stores 1 to 5, respectively. The owner would like to minimize the expected loss due to theft.

Goal Priorities:

Priority 1: requirement (e) is 25 times as important as requirement (d)
Priority 2: requirement (d) is 2 times as important as requirement (c)
Priority 3: requirement (c) as a baseline of priority. Ensure mall stores receive at least 80% of the diamonds.

Traditional Linear Programming Model

In the traditional model, all constraints, including the soft constraints and objectives, are treated as hard constraints.

Decision Variables

$x_i$ : Number of carats allocated to Store $i$ (for $i = 1, 2, 3, 4, 5$ ).

Objective Function

Minimize Expected Loss (e): $Z=0.001x1+0.001x2+0.09x3+0.02x4+0.03x5\text{Minimize } Z = 0.001 x_1 + 0.001 x_2 + 0.09 x_3 + 0.02 x_4 + 0.03 x_5$

Constraints

Hard Constraints (Absolutely Necessary):
- Total Allocation Bounds: $\leq x_1 + x_2 + x_3 + x_4 + x_5 \leq 1200$
- Minimum Allocation to Store 5: $x5≥300x_5 \geq 300$
Supposed to be Soft Constraints:
- Mall Stores Allocation (c): $x1+x2+x3≥0.8(x1+x2+x3+x4+x5)x_1 + x_2 + x_3 \geq 0.8(x_1 + x_2 + x_3 + x_4 + x_5)$
- Equal Allocation Among Mall Stores (d): $x_1 = x_2 = x_3$

Shortcomings of the Traditional Model

Infeasibility: Due to the rigidity of hard constraints and conflicting objectives, the model may become infeasible, offering no solution.
Lack of Flexibility: Cannot prioritize goals or allow for deviations, leaving the decision-maker without actionable solutions.

Goal Programming Approach

Goal programming introduces flexibility by converting objectives and soft constraints into goal constraints with deviational variables, allowing for prioritization and balancing of goals.

Step 1: Convert Soft Constraints and Objectives into Goal Constraints

(c) Mall Stores Allocation Goal

Original Constraint: $x1+x2+x3≥0.8(x1+x2+x3+x4+x5)x_1 + x_2 + x_3 \geq 0.8(x_1 + x_2 + x_3 + x_4 + x_5)$
Rewriting as Goal Constraint: $x_1 + x_2 + x_3 + d_1^- - d_1^+ = 0.8(x_1 + x_2 + x_3 + x_4 + x_5)$
Simplify: $0.2(x_1 + x_2 + x_3) - 0.8 x_4 - 0.8 x_5 + d_1^- - d_1^+ = 0$
Contribution to Objective Function: Minimize $d_1^-$ (underachievement).

(d) Equal Allocation Among Mall Stores

For Store 1 and Store 2:

Original Constraint: $x_1 = x_2$
Rewriting as Goal Constraint: $x_1 - x_2 + d_2^- - d_2^+ = 0$
Contribution to Objective Function: Minimize $d_2^- + d_2^+$ (total deviation).

For Store 1 and Store 3:

Original Constraint: $x_1 = x_3$
Rewriting as Goal Constraint: $x_1 - x_3 + d_3^- - d_3^+ = 0$
Contribution to Objective Function: Minimize $d_3^- + d_3^+$ .

(e) Minimize Expected Loss Due to Theft

Original Objective Function: $Z=0.001x1+0.001x2+0.09x3+0.02x4+0.03x5\text{Minimize } Z = 0.001 x_1 + 0.001 x_2 + 0.09 x_3 + 0.02 x_4 + 0.03 x_5$
Setting a Target (Aspiration Level): Aim for minimal expected loss (e.g., $Z = 0$ , recognizing it may not be attainable).
Rewriting as Goal Constraint: $0.001 x_1 + 0.001 x_2 + 0.09 x_3 + 0.02 x_4 + 0.03 x_5 + d_5^- - d_5^+ = 0$
Contribution to Objective Function: Minimize $d_5^+$ (overachievement).

Step 2: Assign Weights Based on Goal Priorities

Priority 1 (Minimize Expected Loss): Weight $W_5 = 50$ (25 times more important than Goal 2).
Priority 2 (Equal Allocation): Weight $W_2 = W_3 = 2$ (2 times more important than Goal 3).
Priority 3 (Mall Stores Allocation): Weight $W_1 = 1$ .

Step 3: Formulate the Goal Programming Mathematical Model

Decision Variables

$xi≥0x_i \geq 0$ : Allocation to Store $i$ .
$dj−≥0,dj+≥0d_j^- \geq 0, \quad d_j^+ \geq 0$ : Deviational variables for Goal $j$ .

Objective Function

Minimize the weighted sum of deviational variables:

\text{Minimize } Z = 50 d_5^+ + 2 (d_2^- + d_2^+ + d_3^- + d_3^+) + d_1^-

Constraints

Hard Constraints (Absolutely Necessary):

Total Allocation Bounds: $\leq x_1 + x_2 + x_3 + x_4 + x_5 \leq 1200$
Minimum Allocation to Store 5: $x5≥300x_5 \geq 300$

Goal Constraints:

Mall Stores Allocation (c): $0.2 x_1 + 0.2 x_2 + 0.2 x_3 - 0.8 x_4 - 0.8 x_5 + d_1^- - d_1^+ = 0$
Equal Allocation Between Stores 1 and 2 (d): $x_1 - x_2 + d_2^- - d_2^+ = 0$
Equal Allocation Between Stores 1 and 3 (d): $x_1 - x_3 + d_3^- - d_3^+ = 0$
Expected Loss Minimization (e): $0.001 x_1 + 0.001 x_2 + 0.09 x_3 + 0.02 x_4 + 0.03 x_5 + d_5^- - d_5^+ = 0$

Non-Negativity Constraints:

x_i \geq 0, \qquad d_j^- \geq 0, \qquad d_j^+ \geq 0

Step 4: Solve the Goal Programming Model

Using an appropriate optimization solver, we find the values of $x_i$ and $dj±d_j^\pm$ that minimize $Z$ .

Interpretation of Results

Assuming the optimal solution is:

Allocations:
- $x_1 = 466.69$ carats
- $x_2 = 233.31$ carats
- $x_3 = 0$ carats
- $x_4 = 0$ carats
- $x_5 = 300$ carats
Total Allocation:
$x_1 + x_2 + x_3 + x_4 + x_5 = 1000 \text{ carats (meeting the minimum requirement)}$
Deviational Variables:
- $valued_1^- = \text{Positive value}$ : Indicating underachievement of the mall stores allocation goal.
- $d_1^+ = 0$ .
- $valuesd_2^-, d_2^+, d_3^-, d_3^+ = \text{Positive values}$ : Indicating deviations from equal allocation among mall stores.
- $d_5^+ = 9.7$ : Overachievement of expected loss.
- $d_5^- = 0$ .

Advantages of Goal Programming in This Example

Feasibility: Provides a solution where the traditional model may be infeasible due to conflicting constraints.
Prioritization: Allows the decision-maker to focus on the most critical objectives.
Flexibility: Accommodates deviations from less critical goals to achieve the most important ones.
Tunability: Users can adjust weights to explore different trade-offs without altering the problem instance.

Conclusion

Goal programming is a valuable method for addressing complex decision-making problems involving multiple objectives and soft constraints. By incorporating deviational variables and prioritizing goals, it allows decision-makers to find feasible solutions without altering the original problem instance.

Key Takeaways:

Translates Soft Constraints into Goal Constraints: Allows for controlled deviations, reflecting real-world flexibility.
Balances Conflicting Objectives: Prioritizes goals according to their importance.
Maintains Feasibility: Provides solutions even when traditional models fail due to infeasibility.
Adjustable: Enables fine-tuning of models through weight adjustments without changing the problem itself.

Goal programming empowers decision-makers to find practical solutions in complex environments, aligning closely with their objectives and the realities of the constraints they face.

Reference:

(Springer Texts In Business And Economics) H. A. Eiselt, Carl-Louis Sandblom - Operations Research_ A Model-Based Approach-Springer (2022)

A Practical Review: Solving Vehicle Routing Problems with OR-Tools and SCIP

Thana B. — Sun, 06 Oct 2024 09:31:04 +0000

Introduction

Vehicle Routing Problems (VRP) are a class of combinatorial optimization problems that are crucial in the fields of logistics and transportation. The objective is to determine the optimal set of routes for a fleet of vehicles to traverse in order to deliver goods to a given set of customers. The challenge lies in balancing multiple constraints such as vehicle capacities, delivery time windows, and minimizing the total travel time or distance.

In this blog, we will explore how to solve VRP using two powerful optimization tools: Google OR-Tools and SCIP. Google OR-Tools is an open-source software suite for optimization, providing a range of solvers for different types of problems. SCIP, on the other hand, is a framework for Constraint Integer Programming and Mixed Integer Programming, known for its flexibility and performance.

We will implement the VRP with time windows using both OR-Tools and SCIP, compare their performance, and share insights from our experience with these tools.

Problem Statement

The specific VRP we are tackling involves a fleet of vehicles that must deliver goods to a set of customers within specified time windows. The problem can be summarized as follows:

Time Matrix: A matrix representing the travel time between each pair of locations.
Time Windows: Each location (including the depot) has a time window within which the delivery must be made.
Demands: Each customer has a demand that must be satisfied by the vehicles.
Vehicle Capacities: Each vehicle has a maximum capacity that cannot be exceeded.
Number of Vehicles: The fleet consists of a fixed number of vehicles.
Depot: The starting and ending point for all vehicles.

Note that waiting time is allowed and vehicles are not necessary to depart from its depot at 0th minute.

The objective is to minimize the total travel time while satisfying all the constraints. The problem is modeled and solved using Google OR-Tools, and the solution is compared with various models implemented in SCIP.

Tools Overview

Google OR-Tools

Brief Introduction to OR-Tools

Google OR-Tools is an open-source software suite for optimization, developed by Google. It provides a range of solvers for different types of optimization problems, including linear programming, mixed-integer programming, constraint programming, and routing problems. OR-Tools is designed to be easy to use and integrate into various applications, making it a popular choice for both academic research and industrial applications.

Key Features and Advantages

Ease of Use: OR-Tools provides a high-level API that simplifies the process of defining and solving optimization problems. This makes it accessible to users with varying levels of expertise in optimization.
Versatility: It supports a wide range of problem types, including VRP, scheduling, and resource allocation.
Performance: OR-Tools is optimized for performance, leveraging advanced algorithms and heuristics to find solutions quickly.
Cross-Platform: It is available on multiple platforms, including Windows, macOS, and Linux.
Integration: OR-Tools can be easily integrated with other Google services and tools, such as Google Cloud Platform.

Use Cases and Applications

Logistics and Transportation: Solving VRP, optimizing delivery routes, and managing fleet operations.
Scheduling: Creating efficient schedules for employees, machines, or other resources.
Resource Allocation: Optimizing the allocation of resources in various industries, such as manufacturing and healthcare.
Research and Education: Used in academic research and teaching to demonstrate optimization techniques and solve complex problems.

SCIP

Brief Introduction to SCIP

SCIP (Solving Constraint Integer Programs) is a framework for Constraint Integer Programming and Mixed Integer Programming. Developed by the Zuse Institute Berlin, SCIP is known for its flexibility and high performance. It combines techniques from constraint programming, mixed-integer programming, and SAT solving, making it a powerful tool for solving a wide range of optimization problems.

Key Features and Advantages

Flexibility: SCIP allows users to define custom constraints and objective functions, providing a high degree of control over the optimization process.
Performance: SCIP is one of the fastest non-commercial solvers available, often outperforming commercial solvers in benchmark tests.
Extensibility: Users can extend SCIP by adding custom plugins for specific problem types or optimization techniques.
Open Source: SCIP is open-source and freely available for academic and non-commercial use, making it accessible to a wide audience.
Comprehensive Documentation: SCIP provides extensive documentation and examples, helping users get started quickly and effectively.

Use Cases and Applications

Combinatorial Optimization: Solving complex combinatorial problems, such as VRP, scheduling, and resource allocation.
Operations Research: Used in various industries to optimize processes, reduce costs, and improve efficiency.
Academic Research: Widely used in academic research to develop and test new optimization algorithms and techniques.
Industrial Applications: Applied in industries such as manufacturing, logistics, and telecommunications to solve real-world optimization problems.

By leveraging the strengths of both OR-Tools and SCIP, we can effectively tackle the VRP with time windows and gain valuable insights into the performance and usability of these powerful optimization tools.

Implementation Details

Google OR-Tools Implementation

I modified code from
Vehicle Routing Problem with Time Windows and add capacity constraints to slightly increase complexity and located my code in the file named or-tools.py in the same directory with this file. For mathematician or operations research analyst, this may sound weird. When you write code with OR-tools, you don't formulate mathematical model to code. You just must read the documents and understand how it defines time, distance as cumulative variables and how to set constraints for nodes in a route. Please take a look at the here while reading this explanation.

"""Capacited Vehicles Routing Problem (CVRP) with Time Windows."""

from ortools.constraint_solver import routing_enums_pb2
from ortools.constraint_solver import pywrapcp


def create_data_model():
    """Stores the data for the problem."""
    data = {}
    data["time_matrix"] = [
        [0, 6, 9, 8, 7, 3, 6, 2, 3, 2, 6, 6, 4, 4, 5, 9, 7],
        [6, 0, 8, 3, 2, 6, 8, 4, 8, 8, 13, 7, 5, 8, 12, 10, 14],
        [9, 8, 0, 11, 10, 6, 3, 9, 5, 8, 4, 15, 14, 13, 9, 18, 9],
        [8, 3, 11, 0, 1, 7, 10, 6, 10, 10, 14, 6, 7, 9, 14, 6, 16],
        [7, 2, 10, 1, 0, 6, 9, 4, 8, 9, 13, 4, 6, 8, 12, 8, 14],
        [3, 6, 6, 7, 6, 0, 2, 3, 2, 2, 7, 9, 7, 7, 6, 12, 8],
        [6, 8, 3, 10, 9, 2, 0, 6, 2, 5, 4, 12, 10, 10, 6, 15, 5],
        [2, 4, 9, 6, 4, 3, 6, 0, 4, 4, 8, 5, 4, 3, 7, 8, 10],
        [3, 8, 5, 10, 8, 2, 2, 4, 0, 3, 4, 9, 8, 7, 3, 13, 6],
        [2, 8, 8, 10, 9, 2, 5, 4, 3, 0, 4, 6, 5, 4, 3, 9, 5],
        [6, 13, 4, 14, 13, 7, 4, 8, 4, 4, 0, 10, 9, 8, 4, 13, 4],
        [6, 7, 15, 6, 4, 9, 12, 5, 9, 6, 10, 0, 1, 3, 7, 3, 10],
        [4, 5, 14, 7, 6, 7, 10, 4, 8, 5, 9, 1, 0, 2, 6, 4, 8],
        [4, 8, 13, 9, 8, 7, 10, 3, 7, 4, 8, 3, 2, 0, 4, 5, 6],
        [5, 12, 9, 14, 12, 6, 6, 7, 3, 3, 4, 7, 6, 4, 0, 9, 2],
        [9, 10, 18, 6, 8, 12, 15, 8, 13, 9, 13, 3, 4, 5, 9, 0, 9],
        [7, 14, 9, 16, 14, 8, 5, 10, 6, 5, 4, 10, 8, 6, 2, 9, 0],
    ]
    data["time_windows"] = [
        (0, 5),  # depot
        (7, 12),  # 1
        (10, 15),  # 2
        (16, 18),  # 3
        (10, 13),  # 4
        (0, 5),  # 5
        (5, 10),  # 6
        (0, 4),  # 7
        (5, 10),  # 8
        (0, 3),  # 9
        (10, 16),  # 10
        (10, 15),  # 11
        (0, 5),  # 12
        (5, 10),  # 13
        (7, 8),  # 14
        (10, 15),  # 15
        (11, 15),  # 16
    ]
    data["demands"] = [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
    data["vehicle_capacities"] = [15, 15, 15, 15]
    data["num_vehicles"] = 4
    data["depot"] = 0
    return data


def print_solution(data, manager, routing, solution):
    """Prints solution on console."""
    print(f"Objective: {solution.ObjectiveValue()}")
    time_dimension = routing.GetDimensionOrDie("Time")
    total_time = 0
    for vehicle_id in range(data["num_vehicles"]):
        index = routing.Start(vehicle_id)
        plan_output = f"Route for vehicle {vehicle_id}:\n"
        while not routing.IsEnd(index):
            time_var = time_dimension.CumulVar(index)
            plan_output += (
                f"{manager.IndexToNode(index)}"
                f" Time({solution.Min(time_var)},{solution.Max(time_var)})"
                " -> "
            )
            index = solution.Value(routing.NextVar(index))
        time_var = time_dimension.CumulVar(index)
        plan_output += (
            f"{manager.IndexToNode(index)}"
            f" Time({solution.Min(time_var)},{solution.Max(time_var)})\n"
        )
        plan_output += f"Time of the route: {solution.Min(time_var)}min\n"
        print(plan_output)
        total_time += solution.Min(time_var)
    print(f"Total time of all routes: {total_time}min")


def main():
    """Solve the VRP with time windows."""
    # Instantiate the data problem.
    data = create_data_model()

    # Create the routing index manager.
    manager = pywrapcp.RoutingIndexManager(
        len(data["time_matrix"]), data["num_vehicles"], data["depot"]
    )

    # Create Routing Model.
    routing = pywrapcp.RoutingModel(manager)

    # Create and register a transit callback.
    def time_callback(from_index, to_index):
        """Returns the travel time between the two nodes."""
        # Convert from routing variable Index to time matrix NodeIndex.
        from_node = manager.IndexToNode(from_index)
        to_node = manager.IndexToNode(to_index)
        return data["time_matrix"][from_node][to_node]

    transit_callback_index = routing.RegisterTransitCallback(time_callback)

    # Define cost of each arc.
    routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)

    # Add Capacity constraint.
    def demand_callback(from_index):
        """Returns the demand of the node."""
        # Convert from routing variable Index to demands NodeIndex.
        from_node = manager.IndexToNode(from_index)
        return data["demands"][from_node]

    demand_callback_index = routing.RegisterUnaryTransitCallback(demand_callback)
    routing.AddDimensionWithVehicleCapacity(
        demand_callback_index,
        0,  # null capacity slack
        data["vehicle_capacities"],  # vehicle maximum capacities
        True,  # start cumul to zero
        "Capacity",
    )

    # Add Time Windows constraint.
    time = "Time"
    routing.AddDimension(
        transit_callback_index,
        30,  # allow waiting time
        30,  # maximum time per vehicle
        False,  # Don't force start cumul to zero.
        time,
    )
    time_dimension = routing.GetDimensionOrDie(time)
    # Add time window constraints for each location except depot.
    for location_idx, time_window in enumerate(data["time_windows"]):
        if location_idx == data["depot"]:
            continue
        index = manager.NodeToIndex(location_idx)
        time_dimension.CumulVar(index).SetRange(time_window[0], time_window[1])
    # Add time window constraints for each vehicle start node.
    depot_idx = data["depot"]
    for vehicle_id in range(data["num_vehicles"]):
        index = routing.Start(vehicle_id)
        time_dimension.CumulVar(index).SetRange(
            data["time_windows"][depot_idx][0], data["time_windows"][depot_idx][1]
        )

    # Instantiate route start and end times to produce feasible times.
    for i in range(data["num_vehicles"]):
        routing.AddVariableMinimizedByFinalizer(
            time_dimension.CumulVar(routing.Start(i))
        )
        routing.AddVariableMinimizedByFinalizer(time_dimension.CumulVar(routing.End(i)))

    # Setting first solution heuristic.
    search_parameters = pywrapcp.DefaultRoutingSearchParameters()
    search_parameters.first_solution_strategy = (
        routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC
    )

    # Solve the problem.
    solution = routing.SolveWithParameters(search_parameters)

    # Print solution on console.
    if solution:
        print_solution(data, manager, routing, solution)


if __name__ == "__main__":
    main()

The code consists of many components, I will explain each component step-by-step.

Routing Index Manager

manager = pywrapcp.RoutingIndexManager(
    len(data["time_matrix"]), data["num_vehicles"], data["depot"]
)

The RoutingIndexManager is responsible for managing the indexing of nodes (locations) and vehicles. It creates a mapping between the problem nodes and internal indices used by the solver. This mapping is crucial for retrieving parameters/values from the data dictionary, which is indexed by node numbers (indices of the numpy array) rather than the internal indices used by Google OR-Tools.

Creating the Routing Model

routing = pywrapcp.RoutingModel(manager)

The RoutingModel is the core of the routing solver. It uses the RoutingIndexManager to manage the nodes and vehicles. When you initialize this object, think of it as declaring the decision variable $x_{ijk}$ , which represents the decision of whether a vehicle $k$ moves from node $i$ to node $j$ .

Registering the Time Callback

def time_callback(from_index, to_index):
    """Returns the travel time between the two nodes."""
    from_node = manager.IndexToNode(from_index)
    to_node = manager.IndexToNode(to_index)
    return data["time_matrix"][from_node][to_node]

transit_callback_index = routing.RegisterTransitCallback(time_callback)

The time_callback function returns the travel time between two nodes. The RegisterTransitCallback method registers this callback with the routing model, allowing it to be used for calculating the travel time between nodes. The transit_callback_index is an identifier for this callback.

Setting the Objective Function

routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)

This line sets the objective function as the sum of all travel times for the routes taken by the vehicles. Mathematically, it can be expressed as:

\sum_{i \in N} \sum_{j \in N} \sum_{k \in V} c_{ij} \cdot x_{ijk}

where $c_{ij}$ is the travel time from node $i$ to node $j$ , and $x_{ijk}$ is the decision variable indicating whether vehicle $k$ moves from node $i$ to node $j$ .

Registering the Demand Callback

def demand_callback(from_index):
    """Returns the demand of the node."""
    from_node = manager.IndexToNode(from_index)
    return data["demands"][from_node]

demand_callback_index = routing.RegisterUnaryTransitCallback(demand_callback)

The demand_callback function returns the demand of a node. The RegisterUnaryTransitCallback method registers this callback with the routing model, allowing it to be used for calculating the cumulative demand at each node. The demand_callback_index is an identifier for this callback.

Adding the Capacity Constraint

routing.AddDimensionWithVehicleCapacity(
    demand_callback_index,
    0,  # null capacity slack
    data["vehicle_capacities"],  # vehicle maximum capacities
    True,  # start cumul to zero
    "Capacity",
)

This line adds a capacity constraint to the routing model. The AddDimensionWithVehicleCapacity method sets up a rule to cumulate the number of delivered units at each node. The cumulative variable for each vehicle cannot exceed its capacity limit, specified by data["vehicle_capacities"].

Adding the Time Windows Constraint

time = "Time"
routing.AddDimension(
    transit_callback_index,
    30,  # allow waiting time
    30,  # maximum time per vehicle
    False,  # Don't force start cumul to zero.
    time,
)
time_dimension = routing.GetDimensionOrDie(time)

This section adds a time windows constraint to the routing model. The AddDimension method sets up a rule to track the cumulative travel time for each vehicle. The parameters are as follows:

transit_callback_index: The callback function for calculating travel time.
30: The maximum allowed waiting time.
30: The maximum travel time per vehicle.
False: Do not force the start cumulative value to zero.
time: The name of the dimension.

The GetDimensionOrDie method retrieves the time dimension, which is used to add time window constraints for each location.

Adding Time Window Constraints for Each Location

for location_idx, time_window in enumerate(data["time_windows"]):
    if location_idx == data["depot"]:
        continue
    index = manager.NodeToIndex(location_idx)
    time_dimension.CumulVar(index).SetRange(time_window[0], time_window[1])

This loop adds time window constraints for each location except the depot. The SetRange method sets the allowable time window for each location, ensuring that deliveries are made within the specified time windows.

Adding Time Window Constraints for Each Vehicle Start Node

depot_idx = data["depot"]
for vehicle_id in range(data["num_vehicles"]):
    index = routing.Start(vehicle_id)
    time_dimension.CumulVar(index).SetRange(
        data["time_windows"][depot_idx][0], data["time_windows"][depot_idx][1]
    )

This loop adds time window constraints for the start node of each vehicle. The SetRange method ensures that each vehicle starts its route within the specified time window for the depot.

Instantiating Route Start and End Times

for i in range(data["num_vehicles"]):
    routing.AddVariableMinimizedByFinalizer(
        time_dimension.CumulVar(routing.Start(i))
    )
    routing.AddVariableMinimizedByFinalizer(time_dimension.CumulVar(routing.End(i)))

These lines instantiate the start and end times for each vehicle's route, ensuring that the solver produces feasible times for the routes.

Optimizing

# Setting first solution heuristic.
search_parameters = pywrapcp.DefaultRoutingSearchParameters()
search_parameters.first_solution_strategy = (
    routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC
)

# Solve the problem.
solution = routing.SolveWithParameters(search_parameters)

# Print solution on console.
if solution:
    print_solution(data, manager, routing, solution)

This section sets the search parameters for the solver, specifying the first solution heuristic as PATH_CHEAPEST_ARC. The SolveWithParameters method solves the problem using the specified search parameters. If a solution is found, the print_solution function is called to print the solution on the console.

SCIP Implementations

Node-Based Model

Here is the mathematical formulation, I implemented in this code.

from pyscipopt import Model, quicksum
import networkx as nx

def create_data_model():
    """Stores the data for the problem."""
    data = {}
    # Time matrix representing travel times between locations (consistent with time windows)
    data["time_matrix"] = [
        [0, 6, 9, 8, 7, 3, 6, 2, 3, 2, 6, 6, 4, 4, 5, 9, 7],
        [6, 0, 8, 3, 2, 6, 8, 4, 8, 8, 13, 7, 5, 8, 12, 10, 14],
        [9, 8, 0, 11, 10, 6, 3, 9, 5, 8, 4, 15, 14, 13, 9, 18, 9],
        [8, 3, 11, 0, 1, 7, 10, 6, 10, 10, 14, 6, 7, 9, 14, 6, 16],
        [7, 2, 10, 1, 0, 6, 9, 4, 8, 9, 13, 4, 6, 8, 12, 8, 14],
        [3, 6, 6, 7, 6, 0, 2, 3, 2, 2, 7, 9, 7, 7, 6, 12, 8],
        [6, 8, 3, 10, 9, 2, 0, 6, 2, 5, 4, 12, 10, 10, 6, 15, 5],
        [2, 4, 9, 6, 4, 3, 6, 0, 4, 4, 8, 5, 4, 3, 7, 8, 10],
        [3, 8, 5, 10, 8, 2, 2, 4, 0, 3, 4, 9, 8, 7, 3, 13, 6],
        [2, 8, 8, 10, 9, 2, 5, 4, 3, 0, 4, 6, 5, 4, 3, 9, 5],
        [6, 13, 4, 14, 13, 7, 4, 8, 4, 4, 0, 10, 9, 8, 4, 13, 4],
        [6, 7, 15, 6, 4, 9, 12, 5, 9, 6, 10, 0, 1, 3, 7, 3, 10],
        [4, 5, 14, 7, 6, 7, 10, 4, 8, 5, 9, 1, 0, 2, 6, 4, 8],
        [4, 8, 13, 9, 8, 7, 10, 3, 7, 4, 8, 3, 2, 0, 4, 5, 6],
        [5, 12, 9, 14, 12, 6, 6, 7, 3, 3, 4, 7, 6, 4, 0, 9, 2],
        [9, 10, 18, 6, 8, 12, 15, 8, 13, 9, 13, 3, 4, 5, 9, 0, 9],
        [7, 14, 9, 16, 14, 8, 5, 10, 6, 5, 4, 10, 8, 6, 2, 9, 0],
    ]
    data["time_windows"] = [
        (0, 5),     # 0 - depot
        (7, 12),    # 1
        (10, 15),   # 2
        (16, 18),   # 3
        (10, 13),   # 4
        (0, 5),     # 5
        (5, 10),    # 6
        (0, 4),     # 7
        (5, 10),    # 8
        (0, 3),     # 9
        (10, 16),   # 10
        (10, 15),   # 11
        (0, 5),     # 12
        (5, 10),    # 13
        (7, 8),     # 14
        (10, 15),   # 15
        (11, 15),   # 16
    ]
    data["demands"] = [0, 1, 1, 2, 4, 2, 4, 8, 8, 1, 2, 1, 2, 4, 4, 8, 8]
    data["vehicle_capacities"] = [15, 15, 15, 15]  # Capacities for each vehicle
    data["num_vehicles"] = 4
    data["depot"] = 0
    return data

def main():
    """Solve the VRPTW problem using PySCIPOpt with variable depot departure times."""
    # Retrieve data
    data = create_data_model()
    time_matrix = data["time_matrix"]
    time_windows = data["time_windows"]
    demands = data["demands"]
    vehicle_capacities = data["vehicle_capacities"]
    num_vehicles = data["num_vehicles"]
    depot = data["depot"]

    n = len(time_matrix)  # Number of nodes
    N = range(n)  # Set of nodes
    V = [i for i in N if i != depot]  # Customers excluding depot
    K = range(num_vehicles)  # Set of vehicles

    # Calculate M for Big-M constraints (maximum possible time)
    max_time = max(tw[1] for tw in time_windows) + max(
        time_matrix[i][j] for i in N for j in N
    )

    # Initialize the model
    model = Model("VRPTW with Variable Depot Departure Times")

    # Variables
    x = {}  # Binary variables: x[i,j,k] == 1 if vehicle k travels from i to j
    t = {}  # Arrival time at node i by any vehicle
    t_start = {}  # Departure time from depot for vehicle k
    t_end = {}    # Return time to depot for vehicle k
    w = {}  # Waiting time at node i by vehicle k
    y = {}  # Binary variables: y[i,k] == 1 if customer i is visited by vehicle k

    # Decision variables
    for k in K:
        capacity_k = vehicle_capacities[k]
        # Departure time from depot
        t_start[k] = model.addVar(lb=time_windows[depot][0], ub=time_windows[depot][1], vtype="I", name=f"t_start_{k}")
        # Return time to depot
        t_end[k] = model.addVar(lb=time_windows[depot][0], ub=None, vtype="I", name=f"t_end_{k}")
        for i in N:
            if i != depot:
                y[i, k] = model.addVar(vtype="B", name=f"y_{i}_{k}")
            for j in N:
                if i != j:
                    x[i, j, k] = model.addVar(vtype="B", name=f"x_{i}_{j}_{k}")

    for i in N:
        if i != depot:
            t[i] = model.addVar(lb=time_windows[i][0], ub=time_windows[i][1], vtype="I", name=f"t_{i}")

    for i in V:
        for k in K:
            w[i, k] = model.addVar(lb=0, ub=max_time, vtype="I", name=f"w_{i}_{k}")

    # Objective Function: Minimize total operation time of all vehicles
    model.setObjective(
        quicksum(t_end[k] - t_start[k] for k in K),
        "minimize",
    )

    # Constraints

    # 1. Each customer is visited exactly once
    for i in V:
        model.addCons(
            quicksum(y[i, k] for k in K) == 1,
            name=f"VisitOnce_{i}"
        )

    # 2. For each vehicle, ensure it departs from and returns to the depot correctly
    for k in K:
        model.addCons(
            quicksum(x[depot, j, k] for j in N if j != depot) == 1,
            name=f"DepartDepot_{k}"
        )
        model.addCons(
            quicksum(x[i, depot, k] for i in N if i != depot) == 1,
            name=f"ReturnDepot_{k}"
        )

    # 3. Link x and y variables
    for k in K:
        for i in V:
            model.addCons(
                quicksum(x[i, j, k] for j in N if j != i) == y[i, k],
                name=f"LinkXY_{i}_{k}"
            )
            model.addCons(
                quicksum(x[j, i, k] for j in N if j != i) == y[i, k],
                name=f"LinkXY_Inflow_{i}_{k}"
            )

    # 4. Flow conservation constraints for nodes (excluding depot)
    for k in K:
        for i in V:
            inflow = quicksum(x[j, i, k] for j in N if j != i)
            outflow = quicksum(x[i, j, k] for j in N if j != i)
            model.addCons(
                inflow - outflow == 0,
                name=f"FlowConservation_{i}_{k}"
            )

    # 5. Capacity constraints for each vehicle
    for k in K:
        capacity_k = vehicle_capacities[k]
        model.addCons(
            quicksum(demands[i] * y[i, k] for i in V) <= capacity_k,
            name=f"VehicleCapacity_{k}"
        )

    # 6. Time window constraints and time propagation constraints
    for k in K:
        # Departure time from depot is within time window
        model.addCons(
            t_start[k] >= time_windows[depot][0],
            name=f"DepotTimeWindowStart_{k}"
        )
        model.addCons(
            t_start[k] <= time_windows[depot][1],
            name=f"DepotTimeWindowEnd_{k}"
        )

        # Time propagation constraints
        for i in N:
            for j in N:
                if i != j:
                    if i == depot and j != depot:
                        # From depot to first customer
                        model.addCons(
                            t[j] >= t_start[k] + time_matrix[i][j] - max_time * (1 - x[i, j, k]),
                            name=f"TimePropagationDepot_{i}_{j}_{k}"
                        )
                    elif i != depot and j != depot:
                        # Between customers
                        model.addCons(
                            t[j] >= t[i] + w[i, k] + time_matrix[i][j] - max_time * (1 - x[i, j, k]),
                            name=f"TimePropagation_{i}_{j}_{k}"
                        )
                    elif i != depot and j == depot:
                        # Return to depot
                        model.addCons(
                            t_end[k] >= t[i] + w[i, k] + time_matrix[i][j] - max_time * (1 - x[i, j, k]),
                            name=f"TimePropagationReturnDepot_{i}_{j}_{k}"
                        )

        # Waiting time constraints
        for i in V:
            # Waiting time is non-negative
            model.addCons(
                w[i, k] >= 0,
                name=f"WaitingTimeNonNegative_{i}_{k}"
            )
            # If customer i is not visited by vehicle k, waiting time is zero
            model.addCons(
                w[i, k] <= max_time * y[i, k],
                name=f"WaitingTimeZeroIfNotVisited_{i}_{k}"
            )

    # Time windows for customer arrivals
    for i in V:
        model.addCons(
            t[i] >= time_windows[i][0],
            name=f"TimeWindowStart_{i}"
        )
        model.addCons(
            t[i] <= time_windows[i][1],
            name=f"TimeWindowEnd_{i}"
        )

    # Waiting time zero if not visited (already added)

    # Solve the model
    model.optimize()

    # Retrieve and print the solution
    status = model.getStatus()
    if status == "optimal" or status == "bestsollimit":
        print_solution(data, x, t, t_start, t_end, w, K, N, depot, model)
    else:
        print("No optimal solution found.")
        print(f"Solver status: {status}")

def print_solution(data, x, t, t_start, t_end, w, K, N, depot, model):
    """Prints solution in the specified format."""
    solution = model.getBestSol()
    total_time = 0

    print(f"Objective: {model.getSolObjVal(solution):.1f}")
    for k in K:
        vehicle_id = k
        # Build the route for vehicle k
        edges = [(i, j) for i in N for j in N if i != j and model.getSolVal(solution, x.get((i, j, k), 0)) > 0.5]
        if not edges:
            continue  # Skip if vehicle k is not used
        G = nx.DiGraph()
        G.add_edges_from(edges)

        if depot not in G.nodes:
            continue  # Skip if vehicle k is not used

        plan_output = f"Route for vehicle {vehicle_id}:\n"

        index = depot
        time = model.getSolVal(solution, t_start[k])
        time_str = f"Time({time:.1f},{time:.1f})"
        plan_output += f"{index} {time_str} -> "
        total_route_time = time

        while True:
            successors = list(G.successors(index))
            if not successors or successors[0] == depot:
                break
            next_node = successors[0]
            arrival_time = model.getSolVal(solution, t[next_node])
            arrival_time_str = f"Time({arrival_time:.1f},{arrival_time:.1f})"
            plan_output += f"{next_node} {arrival_time_str} -> "
            total_route_time = arrival_time
            index = next_node

        # Return to depot
        return_time = model.getSolVal(solution, t_end[k])
        arrival_time_str = f"Time({return_time:.1f},{return_time:.1f})"
        plan_output += f"{depot} {arrival_time_str}\n"
        route_duration = return_time - model.getSolVal(solution, t_start[k])
        plan_output += f"Time of the route: {route_duration:.1f}min\n"
        print(plan_output)
        total_time += route_duration

    print(f"Total time of all routes: {total_time:.1f}min")

if __name__ == "__main__":
    main()

Sets and Indices

$N$ : Set of all nodes (including depot)
$V$ : Set of customer nodes (excluding depot)
$K$ : Set of vehicles
$\in N$ : Indices for nodes
$\in K$ : Index for vehicles

Parameters

$c_{ij}$ : Travel time from node $i$ to node $j$
$a_i, b_i]$ : Time window for node $i$
$d_i$ : Demand at node $i$
$Q_k$ : Capacity of vehicle $k$
$M$ : A large constant (Big-M)

Decision Variables

$x_{ijk}$ : Binary variable, 1 if vehicle $k$ travels from node $i$ to node $j$ , 0 otherwise
$t_i$ : Arrival time at node $i$
$t_{start_k}$ : Departure time from depot for vehicle $k$
$t_{end_k}$ : Return time to depot for vehicle $k$
$w_{ik}$ : Waiting time at node $i$ by vehicle $k$
$y_{ik}$ : Binary variable, 1 if customer $i$ is visited by vehicle $k$ , 0 otherwise

Objective Function

Minimize the total operation time of all vehicles:

\text{Minimize} \quad \sum_{k \in K} (t_{end_k} - t_{start_k})

Constraints

1. Each customer is visited exactly once

\sum_{k \in K} y_{ik} = 1 \quad \forall i \in V

2. Each vehicle departs from and returns to the depot exactly once

\sum_{j \in N, j \neq i} x_{0jk} = 1 \quad \forall k \in K

\sum_{i \in N, i \neq j} x_{i0k} = 1 \quad \forall k \in K

3. Link $x$ and $y$ variables

\sum_{j \in N, j \neq i} x_{ijk} = y_{ik} \quad \forall i \in V, \forall k \in K

\sum_{j \in N, j \neq i} x_{jik} = y_{ik} \quad \forall i \in V, \forall k \in K

4. Flow conservation constraints for nodes (excluding depot)

\sum_{j \in N, j \neq i} x_{jik} - \sum_{j \in N, j \neq i} x_{ijk} = 0 \quad \forall i \in V, \forall k \in K

5. Capacity constraints for each vehicle

\sum_{i \in V} d_i y_{ik} \leq Q_k \quad \forall k \in K

6. Time window constraints and time propagation constraints

Departure time from depot is within time window

a_0 \leq t_{\text{start}_k} \leq b_0 \quad \forall k \in K

Time propagation constraints

From depot to first customer:

$t_j \geq t_{start_k} + c_{0j} - M(1 - x_{0jk}) \quad \forall j \in V, \forall k \in K$
Between customers:

$t_j \geq t_i + w_{ik} + c_{ij} - M(1 - x_{ijk}) \quad \forall i, j \in V, i \neq j, \forall k \in K$
Return to depot:

$t_{end_k} \geq t_i + w_{ik} + c_{i0} - M(1 - x_{i0k}) \quad \forall i \in V, \forall k \in K$

Flow-based Model

Here is the mathematical formulation, I implemented in this code.

from pyscipopt import Model, quicksum
import networkx as nx

def create_data_model():
    """Stores the data for the problem."""
    data = {}
    # Time matrix representing travel times between locations (consistent with time windows)
    data["time_matrix"] = [
        [0, 6, 9, 8, 7, 3, 6, 2, 3, 2, 6, 6, 4, 4, 5, 9, 7],
        [6, 0, 8, 3, 2, 6, 8, 4, 8, 8, 13, 7, 5, 8, 12, 10, 14],
        [9, 8, 0, 11, 10, 6, 3, 9, 5, 8, 4, 15, 14, 13, 9, 18, 9],
        [8, 3, 11, 0, 1, 7, 10, 6, 10, 10, 14, 6, 7, 9, 14, 6, 16],
        [7, 2, 10, 1, 0, 6, 9, 4, 8, 9, 13, 4, 6, 8, 12, 8, 14],
        [3, 6, 6, 7, 6, 0, 2, 3, 2, 2, 7, 9, 7, 7, 6, 12, 8],
        [6, 8, 3, 10, 9, 2, 0, 6, 2, 5, 4, 12, 10, 10, 6, 15, 5],
        [2, 4, 9, 6, 4, 3, 6, 0, 4, 4, 8, 5, 4, 3, 7, 8, 10],
        [3, 8, 5, 10, 8, 2, 2, 4, 0, 3, 4, 9, 8, 7, 3, 13, 6],
        [2, 8, 8, 10, 9, 2, 5, 4, 3, 0, 4, 6, 5, 4, 3, 9, 5],
        [6, 13, 4, 14, 13, 7, 4, 8, 4, 4, 0, 10, 9, 8, 4, 13, 4],
        [6, 7, 15, 6, 4, 9, 12, 5, 9, 6, 10, 0, 1, 3, 7, 3, 10],
        [4, 5, 14, 7, 6, 7, 10, 4, 8, 5, 9, 1, 0, 2, 6, 4, 8],
        [4, 8, 13, 9, 8, 7, 10, 3, 7, 4, 8, 3, 2, 0, 4, 5, 6],
        [5, 12, 9, 14, 12, 6, 6, 7, 3, 3, 4, 7, 6, 4, 0, 9, 2],
        [9, 10, 18, 6, 8, 12, 15, 8, 13, 9, 13, 3, 4, 5, 9, 0, 9],
        [7, 14, 9, 16, 14, 8, 5, 10, 6, 5, 4, 10, 8, 6, 2, 9, 0],
    ]
    data["time_windows"] = [
        (0, 5),     # 0 - depot
        (7, 12),    # 1
        (10, 15),   # 2
        (16, 18),   # 3
        (10, 13),   # 4
        (0, 5),     # 5
        (5, 10),    # 6
        (0, 4),     # 7
        (5, 10),    # 8
        (0, 3),     # 9
        (10, 16),   # 10
        (10, 15),   # 11
        (0, 5),     # 12
        (5, 10),    # 13
        (7, 8),     # 14
        (10, 15),   # 15
        (11, 15),   # 16
    ]
    data["demands"] = [0, 1, 1, 2, 4, 2, 4, 8, 8, 1, 2, 1, 2, 4, 4, 8, 8]
    data["vehicle_capacities"] = [15, 15, 15, 15]  # Capacities for each vehicle
    data["num_vehicles"] = 4
    data["depot"] = 0
    return data

def main():
    """Solve the VRPTW problem using PySCIPOpt with a flow-based model."""
    # Retrieve data
    data = create_data_model()
    time_matrix = data["time_matrix"]
    time_windows = data["time_windows"]
    demands = data["demands"]
    vehicle_capacities = data["vehicle_capacities"]
    num_vehicles = data["num_vehicles"]
    depot = data["depot"]

    n = len(time_matrix)  # Number of nodes
    N = range(n)  # Set of nodes
    V = [i for i in N if i != depot]  # Customers excluding depot
    K = range(num_vehicles)  # Set of vehicles

    # Calculate M for Big-M constraints (maximum possible time)
    max_time = max(tw[1] for tw in time_windows) + max(
        time_matrix[i][j] for i in N for j in N
    )

    # Initialize the model
    model = Model("VRPTW with Flow-Based Model")

    # Variables
    x = {}  # Binary variables: x[i,j,k] == 1 if vehicle k travels from i to j
    f = {}  # Continuous variables: flow of demand on arc (i,j) by vehicle k
    t = {}  # Arrival time at node i by any vehicle
    t_start = {}  # Departure time from depot for vehicle k
    t_end = {}    # Return time to depot for vehicle k
    w = {}  # Waiting time at node i by vehicle k
    y = {}  # Binary variables: y[i,k] == 1 if customer i is visited by vehicle k

    # Decision variables
    for k in K:
        capacity_k = vehicle_capacities[k]
        # Departure time from depot
        t_start[k] = model.addVar(lb=time_windows[depot][0], ub=time_windows[depot][1], vtype="I", name=f"t_start_{k}")
        # Return time to depot
        t_end[k] = model.addVar(lb=time_windows[depot][0], ub=None, vtype="I", name=f"t_end_{k}")
        for i in N:
            if i != depot:
                y[i, k] = model.addVar(vtype="B", name=f"y_{i}_{k}")
            for j in N:
                if i != j:
                    x[i, j, k] = model.addVar(vtype="B", name=f"x_{i}_{j}_{k}")
                    # Flow variable
                    f[i, j, k] = model.addVar(lb=0, vtype="I", name=f"f_{i}_{j}_{k}")

    for i in N:
        if i != depot:
            t[i] = model.addVar(lb=time_windows[i][0], ub=time_windows[i][1], vtype="I", name=f"t_{i}")

    for i in V:
        for k in K:
            w[i, k] = model.addVar(lb=0, ub=max_time, vtype="I", name=f"w_{i}_{k}")

    # Objective Function: Minimize total operation time of all vehicles
    model.setObjective(
        quicksum(t_end[k] - t_start[k] for k in K),
        "minimize",
    )

    # Constraints

    # 1. Each customer is visited exactly once
    for i in V:
        model.addCons(
            quicksum(y[i, k] for k in K) == 1,
            name=f"VisitOnce_{i}"
        )

    # 2. For each vehicle, ensure it departs from and returns to the depot correctly
    for k in K:
        model.addCons(
            quicksum(x[depot, j, k] for j in N if j != depot) == 1,
            name=f"DepartDepot_{k}"
        )
        model.addCons(
            quicksum(x[i, depot, k] for i in N if i != depot) == 1,
            name=f"ReturnDepot_{k}"
        )

    # 3. Link x and y variables
    for k in K:
        for i in V:
            model.addCons(
                quicksum(x[i, j, k] for j in N if j != i) == y[i, k],
                name=f"LinkXY_{i}_{k}"
            )
            model.addCons(
                quicksum(x[j, i, k] for j in N if j != i) == y[i, k],
                name=f"LinkXY_Inflow_{i}_{k}"
            )

    # 4. Flow conservation constraints (using flow variables)
    for k in K:
        capacity_k = vehicle_capacities[k]
        # Capacity constraints on arcs
        for i in N:
            for j in N:
                if i != j:
                    model.addCons(
                        f[i, j, k] <= capacity_k * x[i, j, k],
                        name=f"CapacityArc_{i}_{j}_{k}"
                    )
        # Flow conservation at depot
        inflow_f = quicksum(f[j, depot, k] for j in N if j != depot)
        outflow_f = quicksum(f[depot, j, k] for j in N if j != depot)
        total_demand_k = quicksum(demands[i] * y[i, k] for i in V)
        model.addCons(
            inflow_f - outflow_f == -total_demand_k,
            name=f"FlowConservation_Depot_{k}"
        )
        # Flow conservation at customers
        for i in V:
            inflow_f = quicksum(f[j, i, k] for j in N if j != i)
            outflow_f = quicksum(f[i, j, k] for j in N if j != i)
            model.addCons(
                inflow_f - outflow_f == demands[i] * y[i, k],
                name=f"FlowConservation_{i}_{k}"
            )

    # 5. Time window constraints and time propagation constraints
    for k in K:
        # Departure time from depot is within time window
        model.addCons(
            t_start[k] >= time_windows[depot][0],
            name=f"DepotTimeWindowStart_{k}"
        )
        model.addCons(
            t_start[k] <= time_windows[depot][1],
            name=f"DepotTimeWindowEnd_{k}"
        )

        # Time propagation constraints
        for i in N:
            for j in N:
                if i != j:
                    if i == depot and j != depot:
                        # From depot to first customer
                        model.addCons(
                            t[j] >= t_start[k] + time_matrix[i][j] - max_time * (1 - x[i, j, k]),
                            name=f"TimePropagationDepot_{i}_{j}_{k}"
                        )
                    elif i != depot and j != depot:
                        # Between customers
                        model.addCons(
                            t[j] >= t[i] + w[i, k] + time_matrix[i][j] - max_time * (1 - x[i, j, k]),
                            name=f"TimePropagation_{i}_{j}_{k}"
                        )
                    elif i != depot and j == depot:
                        # Return to depot
                        model.addCons(
                            t_end[k] >= t[i] + w[i, k] + time_matrix[i][j] - max_time * (1 - x[i, j, k]),
                            name=f"TimePropagationReturnDepot_{i}_{j}_{k}"
                        )

        # Waiting time constraints
        for i in V:
            # Waiting time is non-negative
            model.addCons(w[i, k] >= 0,
                name=f"WaitingTimeNonNegative_{i}_{k}"
            )
            # If customer i is not visited by vehicle k, waiting time is zero
            model.addCons(
                w[i, k] <= max_time * y[i, k],
                name=f"WaitingTimeZeroIfNotVisited_{i}_{k}"
            )

    # Time windows for customer arrivals
    for i in V:
        model.addCons(
            t[i] >= time_windows[i][0],
            name=f"TimeWindowStart_{i}"
        )
        model.addCons(
            t[i] <= time_windows[i][1],
            name=f"TimeWindowEnd_{i}"
        )

    # Solve the model
    model.optimize()

    # Retrieve and print the solution
    status = model.getStatus()
    if status == "optimal" or status == "bestsollimit":
        print_solution(data, x, t, t_start, t_end, w, K, N, depot, model)
    else:
        print("No optimal solution found.")
        print(f"Solver status: {status}")

def print_solution(data, x, t, t_start, t_end, w, K, N, depot, model):
    """Prints solution in the specified format."""
    solution = model.getBestSol()
    total_time = 0

    print(f"Objective: {model.getSolObjVal(solution):.1f}")
    for k in K:
        vehicle_id = k
        # Build the route for vehicle k
        edges = [(i, j) for i in N for j in N if i != j and model.getSolVal(solution, x.get((i, j, k), 0)) > 0.5]
        if not edges:
            continue  # Skip if vehicle k is not used
        G = nx.DiGraph()
        G.add_edges_from(edges)

        if depot not in G.nodes:
            continue  # Skip if vehicle k is not used

        plan_output = f"Route for vehicle {vehicle_id}:\n"

        index = depot
        time = model.getSolVal(solution, t_start[k])
        time_str = f"Time({time:.1f},{time:.1f})"
        plan_output += f"{index} {time_str} -> "
        total_route_time = time

        while True:
            successors = list(G.successors(index))
            if not successors or successors[0] == depot:
                break
            next_node = successors[0]
            arrival_time = model.getSolVal(solution, t[next_node])
            arrival_time_str = f"Time({arrival_time:.1f},{arrival_time:.1f})"
            plan_output += f"{next_node} {arrival_time_str} -> "
            total_route_time = arrival_time
            index = next_node

        # Return to depot
        return_time = model.getSolVal(solution, t_end[k])
        arrival_time_str = f"Time({return_time:.1f},{return_time:.1f})"
        plan_output += f"{depot} {arrival_time_str}\n"
        route_duration = return_time - model.getSolVal(solution, t_start[k])
        plan_output += f"Time of the route: {route_duration:.1f}min\n"
        print(plan_output)
        total_time += route_duration

    print(f"Total time of all routes: {total_time:.1f}min")

if __name__ == "__main__":
    main()

Sets and Indices

$N$ : Set of all nodes (including depot)
$V$ : Set of customer nodes (excluding depot)
$K$ : Set of vehicles
$\in N$ : Indices for nodes
$\in K$ : Index for vehicles

Parameters

$c_{ij}$ : Travel time from node $i$ to node $j$
$a_i, b_i]$ : Time window for node $i$
$d_i$ : Demand at node $i$
$Q_k$ : Capacity of vehicle $k$
$M$ : A large constant (Big-M)

Decision Variables

$x_{ijk}$ : Binary variable, 1 if vehicle $k$ travels from node $i$ to node $j$ , 0 otherwise
$f_{ijk}$ : Continuous variable, flow of demand on arc $(i, j)$ by vehicle $k$
$t_i$ : Arrival time at node $i$
$t_{start_k}$ : Departure time from depot for vehicle $k$
$t_{end_k}$ : Return time to depot for vehicle $k$
$w_{ik}$ : Waiting time at node $i$ by vehicle $k$
$y_{ik}$ : Binary variable, 1 if customer $i$ is visited by vehicle $k$ , 0 otherwise

Objective Function

Minimize the total operation time of all vehicles:

\text{Minimize} \quad \sum_{k \in K} (t_{end_k} - t_{start_k})

Constraints

1. Each customer is visited exactly once

\sum_{k \in K} y_{ik} = 1 \quad \forall i \in V

2. Each vehicle departs from and returns to the depot exactly once

\sum_{j \in N, j \neq i} x_{0jk} = 1 \quad \forall k \in K

\sum_{i \in N, i \neq j} x_{i0k} = 1 \quad \forall k \in K

3. Link $x$ and $y$ variables

\sum_{j \in N, j \neq i} x_{ijk} = y_{ik} \quad \forall i \in V, \forall k \in K

\sum_{j \in N, j \neq i} x_{jik} = y_{ik} \quad \forall i \in V, \forall k \in K

4. Flow conservation constraints (using flow variables)

Capacity constraints on arcs:

$f_{ijk} \leq Q_k \cdot x_{ijk} \quad \forall i, j \in N, i \neq j, \forall k \in K$
Flow conservation at depot:

$\sum_{j \in N, j \neq 0} f_{j0k} - \sum_{j \in N, j \neq 0} f_{0jk} = -\sum_{i \in V} d_i \cdot y_{ik} \quad \forall k \in K$
Flow conservation at customers:

$\sum_{j \in N, j \neq i} f_{jik} - \sum_{j \in N, j \neq i} f_{ijk} = d_i \cdot y_{ik} \quad \forall i \in V, \forall k \in K$

5. Time window constraints and time propagation constraints

Departure time from depot is within time window:

$a_0 \leq t_{start_k} \leq b_0 \quad \forall k \in K$
Time propagation constraints:
- From depot to first customer: $tj≥tstartk+c0j−M(1−x0jk)∀j∈V,∀k∈Kt_j \geq t_{start_k} + c_{0j} - M(1 - x_{0jk}) \quad \forall j \in V, \forall k \in K$
- Between customers: $tj≥ti+wik+cij−M(1−xijk)∀i,j∈V,i≠j,∀k∈Kt_j \geq t_i + w_{ik} + c_{ij} - M(1 - x_{ijk}) \quad \forall i, j \in V, i \neq j, \forall k \in K$
- Return to depot: $tendk≥ti+wik+ci0−M(1−xi0k)∀i∈V,∀k∈Kt_{end_k} \geq t_i + w_{ik} + c_{i0} - M(1 - x_{i0k}) \quad \forall i \in V, \forall k \in K$
Waiting time constraints:
- Waiting time is non-negative: $wik≥0∀i∈V,∀k∈Kw_{ik} \geq 0 \quad \forall i \in V, \forall k \in K$
- If customer $i$ is not visited by vehicle $k$ , waiting time is zero: $wik≤M⋅yik∀i∈V,∀k∈Kw_{ik} \leq M \cdot y_{ik} \quad \forall i \in V, \forall k \in K$

6. Time windows for customer arrivals

a_i \leq t_i \leq b_i \quad \forall i \in V

MTZ Model

Here is the mathematical formulation, I implemented in this code.

from pyscipopt import Model, quicksum
import networkx as nx

def create_data_model():
    """Stores the data for the problem."""
    data = {}
    # Time matrix representing travel times between locations (consistent with time windows)
    data["time_matrix"] = [
        [0, 6, 9, 8, 7, 3, 6, 2, 3, 2, 6, 6, 4, 4, 5, 9, 7],
        [6, 0, 8, 3, 2, 6, 8, 4, 8, 8, 13, 7, 5, 8, 12, 10, 14],
        [9, 8, 0, 11, 10, 6, 3, 9, 5, 8, 4, 15, 14, 13, 9, 18, 9],
        [8, 3, 11, 0, 1, 7, 10, 6, 10, 10, 14, 6, 7, 9, 14, 6, 16],
        [7, 2, 10, 1, 0, 6, 9, 4, 8, 9, 13, 4, 6, 8, 12, 8, 14],
        [3, 6, 6, 7, 6, 0, 2, 3, 2, 2, 7, 9, 7, 7, 6, 12, 8],
        [6, 8, 3, 10, 9, 2, 0, 6, 2, 5, 4, 12, 10, 10, 6, 15, 5],
        [2, 4, 9, 6, 4, 3, 6, 0, 4, 4, 8, 5, 4, 3, 7, 8, 10],
        [3, 8, 5, 10, 8, 2, 2, 4, 0, 3, 4, 9, 8, 7, 3, 13, 6],
        [2, 8, 8, 10, 9, 2, 5, 4, 3, 0, 4, 6, 5, 4, 3, 9, 5],
        [6, 13, 4, 14, 13, 7, 4, 8, 4, 4, 0, 10, 9, 8, 4, 13, 4],
        [6, 7, 15, 6, 4, 9, 12, 5, 9, 6, 10, 0, 1, 3, 7, 3, 10],
        [4, 5, 14, 7, 6, 7, 10, 4, 8, 5, 9, 1, 0, 2, 6, 4, 8],
        [4, 8, 13, 9, 8, 7, 10, 3, 7, 4, 8, 3, 2, 0, 4, 5, 6],
        [5, 12, 9, 14, 12, 6, 6, 7, 3, 3, 4, 7, 6, 4, 0, 9, 2],
        [9, 10, 18, 6, 8, 12, 15, 8, 13, 9, 13, 3, 4, 5, 9, 0, 9],
        [7, 14, 9, 16, 14, 8, 5, 10, 6, 5, 4, 10, 8, 6, 2, 9, 0],
    ]
    data["time_windows"] = [
        (0, 5),     # 0 - depot
        (7, 12),    # 1
        (10, 15),   # 2
        (16, 18),   # 3
        (10, 13),   # 4
        (0, 5),     # 5
        (5, 10),    # 6
        (0, 4),     # 7
        (5, 10),    # 8
        (0, 3),     # 9
        (10, 16),   # 10
        (10, 15),   # 11
        (0, 5),     # 12
        (5, 10),    # 13
        (7, 8),     # 14
        (10, 15),   # 15
        (11, 15),   # 16
    ]
    data["demands"] = [0, 1, 1, 2, 4, 2, 4, 8, 8, 1, 2, 1, 2, 4, 4, 8, 8]
    data["vehicle_capacities"] = [15, 15, 15, 15]  # Capacities for each vehicle
    data["num_vehicles"] = 4
    data["depot"] = 0
    return data

def main():
    """Solve the VRPTW problem using PySCIPOpt with MTZ constraints."""
    # Retrieve data
    data = create_data_model()
    time_matrix = data["time_matrix"]
    time_windows = data["time_windows"]
    demands = data["demands"]
    vehicle_capacities = data["vehicle_capacities"]
    num_vehicles = data["num_vehicles"]
    depot = data["depot"]

    n = len(time_matrix)  # Number of nodes
    N = range(n)  # Set of nodes
    V = [i for i in N if i != depot]  # Customers excluding depot
    K = range(num_vehicles)  # Set of vehicles

    # Calculate M for Big-M constraints (maximum possible time)
    max_time = max(tw[1] for tw in time_windows) + max(
        time_matrix[i][j] for i in N for j in N
    )

    # Initialize the model
    model = Model("VRPTW with MTZ Constraints")

    # Variables
    x = {}  # Binary variables: x[i,j,k] == 1 if vehicle k travels from i to j
    t = {}  # Arrival time at node i by vehicle k
    t_start = {}  # Departure time from depot for vehicle k
    t_end = {}    # Return time to depot for vehicle k
    y = {}  # Binary variables: y[i,k] == 1 if vehicle k visits node i
    u = {}  # MTZ variables for subtour elimination

    # Decision variables
    for k in K:
        capacity_k = vehicle_capacities[k]
        # Departure time from depot
        t_start[k] = model.addVar(lb=time_windows[depot][0], ub=time_windows[depot][1], vtype="I", name=f"t_start_{k}")
        # Return time to depot
        t_end[k] = model.addVar(lb=time_windows[depot][0], ub=None, vtype="I", name=f"t_end_{k}")
        u[depot, k] = model.addVar(lb=0, ub=n, vtype="I", name=f"u_{depot}_{k}")  # u at depot is 0
        model.addCons(u[depot, k] == 0, name=f"MTZ_depot_{k}")
        for i in N:
            if i != depot:
                y[i, k] = model.addVar(vtype="B", name=f"y_{i}_{k}")
                u[i, k] = model.addVar(lb=1, ub=n, vtype="I", name=f"u_{i}_{k}")
                t[i, k] = model.addVar(lb=time_windows[i][0], ub=time_windows[i][1], vtype="I", name=f"t_{i}_{k}")
            for j in N:
                if i != j:
                    x[i, j, k] = model.addVar(vtype="B", name=f"x_{i}_{j}_{k}")

    # Objective Function: Minimize total operation time of all vehicles
    model.setObjective(
        quicksum(t_end[k] - t_start[k] for k in K),
        "minimize",
    )

    # Constraints

    # 1. Each customer is visited exactly once
    for i in V:
        model.addCons(
            quicksum(y[i, k] for k in K) == 1,
            name=f"VisitOnce_{i}"
        )

    # 2. For each vehicle, ensure it departs from and returns to the depot correctly
    for k in K:
        model.addCons(
            quicksum(x[depot, j, k] for j in N if j != depot) == 1,
            name=f"DepartDepot_{k}"
        )
        model.addCons(
            quicksum(x[i, depot, k] for i in N if i != depot) == 1,
            name=f"ReturnDepot_{k}"
        )

    # 3. Flow conservation constraints for nodes (excluding depot)
    for k in K:
        for i in V:
            model.addCons(
                quicksum(x[i, j, k] for j in N if j != i) - quicksum(x[j, i, k] for j in N if j != i) == 0,
                name=f"FlowConservation_{i}_{k}"
            )

    # 4. Link x and y variables
    for k in K:
        for i in V:
            model.addCons(
                quicksum(x[i, j, k] for j in N if j != i) == y[i, k],
                name=f"LinkXY_{i}_{k}"
            )

    # 5. Capacity constraints for each vehicle
    for k in K:
        capacity_k = vehicle_capacities[k]
        model.addCons(
            quicksum(demands[i] * y[i, k] for i in V) <= capacity_k,
            name=f"VehicleCapacity_{k}"
        )

    # 6. Time window constraints and time propagation constraints
    for k in K:
        # Departure time from depot is within time window
        model.addCons(
            t_start[k] >= time_windows[depot][0],
            name=f"DepotTimeWindowStart_{k}"
        )
        model.addCons(
            t_start[k] <= time_windows[depot][1],
            name=f"DepotTimeWindowEnd_{k}"
        )

        # Time propagation constraints
        for i in N:
            if i != depot:
                # Arrival time at node i for vehicle k
                model.addCons(
                    t[i, k] >= time_windows[i][0],
                    name=f"TimeWindowStart_{i}_{k}"
                )
                model.addCons(
                    t[i, k] <= time_windows[i][1],
                    name=f"TimeWindowEnd_{i}_{k}"
                )
            for j in N:
                if i != j:
                    travel_time = time_matrix[i][j]
                    if i == depot:
                        # From depot to first customer
                        model.addCons(
                            t[j, k] >= t_start[k] + travel_time - max_time * (1 - x[i, j, k]),
                            name=f"TimePropagationDepot_{i}_{j}_{k}"
                        )
                    elif j == depot:
                        # From last customer to depot
                        model.addCons(
                            t_end[k] >= t[i, k] + travel_time - max_time * (1 - x[i, j, k]),
                            name=f"TimePropagationReturnDepot_{i}_{j}_{k}"
                        )
                    else:
                        # Between customers
                        model.addCons(
                            t[j, k] >= t[i, k] + travel_time - max_time * (1 - x[i, j, k]),
                            name=f"TimePropagation_{i}_{j}_{k}"
                        )

    # 7. Subtour elimination (MTZ) constraints
    for k in K:
        for i in N:
            if i != depot:
                model.addCons(
                    u[i, k] >= 1,
                    name=f"MTZ_LB_{i}_{k}"
                )
                model.addCons(
                    u[i, k] <= n,
                    name=f"MTZ_UB_{i}_{k}"
                )
        for i in N:
            for j in N:
                if i != j and i != depot and j != depot:
                    model.addCons(
                        u[i, k] - u[j, k] + n * x[i, j, k] <= n - 1,
                        name=f"MTZ_{i}_{j}_{k}"
                    )

    # 8. Ensure that if a customer is not visited by vehicle k, its time variables are inactive
    for k in K:
        for i in V:
            model.addCons(
                y[i, k] <= quicksum(x[j, i, k] for j in N if j != i),
                name=f"VisitIndicator_{i}_{k}"
            )

    # Solve the model
    model.optimize()

    # Retrieve and print the solution
    status = model.getStatus()
    if status == "optimal" or status == "bestsollimit":
        print_solution(data, x, t_start, t_end, t, y, K, N, depot, model)
    else:
        print("No optimal solution found.")
        print(f"Solver status: {status}")

def print_solution(data, x, t_start, t_end, t, y, K, N, depot, model):
    """Prints solution in the specified format."""
    solution = model.getBestSol()
    total_time = 0

    print(f"Objective: {model.getSolObjVal(solution):.1f}")
    for k in K:
        vehicle_id = k
        # Build the route for vehicle k
        edges = [(i, j) for i in N for j in N if i != j and model.getSolVal(solution, x.get((i, j, k), 0)) > 0.5]
        if not edges:
            continue  # Skip if vehicle k is not used
        G = nx.DiGraph()
        G.add_edges_from(edges)

        if depot not in G.nodes:
            continue  # Skip if vehicle k is not used

        plan_output = f"Route for vehicle {vehicle_id}:\n"

        route = [depot]
        time_log = [f"Time({model.getSolVal(solution, t_start[k]):.1f})"]
        next_node = depot

        while True:
            successors = list(G.successors(next_node))
            if not successors:
                break
            next_node = successors[0]
            if next_node == depot:
                break
            route.append(next_node)
            arrival_time = model.getSolVal(solution, t[next_node, k])
            time_log.append(f"Time({arrival_time:.1f})")

        # Append depot at the end
        route.append(depot)
        time_log.append(f"Time({model.getSolVal(solution, t_end[k]):.1f})")

        # Print the route and times
        for idx, node in enumerate(route):
            plan_output += f"{node} {time_log[idx]} -> "
        plan_output = plan_output.rstrip(" -> ") + "\n"
        route_duration = model.getSolVal(solution, t_end[k]) - model.getSolVal(solution, t_start[k])
        plan_output += f"Time of the route: {route_duration:.1f}min\n"
        print(plan_output)
        total_time += route_duration

    print(f"Total time of all routes: {total_time:.1f}min")

if __name__ == "__main__":
    main()

Sets and Indices

$N$ : Set of all nodes (including depot)
$V$ : Set of customer nodes (excluding depot)
$K$ : Set of vehicles
$\in N$ : Indices for nodes
$\in K$ : Index for vehicles

Parameters

$c_{ij}$ : Travel time from node $i$ to node $j$
$a_i, b_i]$ : Time window for node $i$
$d_i$ : Demand at node $i$
$Q_k$ : Capacity of vehicle $k$
$M$ : A large constant (Big-M)

Decision Variables

$x_{ijk}$ : Binary variable, 1 if vehicle $k$ travels from node $i$ to node $j$ , 0 otherwise
$t_{ik}$ : Arrival time at node $i$ by vehicle $k$
$t_{start_k}$ : Departure time from depot for vehicle $k$
$t_{end_k}$ : Return time to depot for vehicle $k$
$y_{ik}$ : Binary variable, 1 if vehicle $k$ visits node $i$ , 0 otherwise
$u_{ik}$ : MTZ variable for subtour elimination

Objective Function

Minimize the total operation time of all vehicles:

\text{Minimize} \quad \sum_{k \in K} (t_{end_k} - t_{start_k})

Constraints

1. Each customer is visited exactly once

\sum_{k \in K} y_{ik} = 1 \quad \forall i \in V

2. Each vehicle departs from and returns to the depot exactly once

\sum_{j \in N, j \neq i} x_{0jk} = 1 \quad \forall k \in K

\sum_{i \in N, i \neq j} x_{i0k} = 1 \quad \forall k \in K

3. Flow conservation constraints for nodes (excluding depot)

\sum_{j \in N, j \neq i} x_{ijk} - \sum_{j \in N, j \neq i} x_{jik} = 0 \quad \forall i \in V, \forall k \in K

4. Link $x$ and $y$ variables

\sum_{j \in N, j \neq i} x_{ijk} = y_{ik} \quad \forall i \in V, \forall k \in K

5. Capacity constraints for each vehicle

\sum_{i \in V} d_i y_{ik} \leq Q_k \quad \forall k \in K

6. Time window constraints and time propagation constraints

Departure time from depot is within time window:

$a_0 \leq t_{start_k} \leq b_0 \quad \forall k \in K$
Time propagation constraints:
- From depot to first customer: $tjk≥tstartk+c0j−M(1−x0jk)∀j∈V,∀k∈Kt_{jk} \geq t_{start_k} + c_{0j} - M(1 - x_{0jk}) \quad \forall j \in V, \forall k \in K$
- Between customers: $tjk≥tik+cij−M(1−xijk)∀i,j∈V,i≠j,∀k∈Kt_{jk} \geq t_{ik} + c_{ij} - M(1 - x_{ijk}) \quad \forall i, j \in V, i \neq j, \forall k \in K$
- Return to depot: $tendk≥tik+ci0−M(1−xi0k)∀i∈V,∀k∈Kt_{end_k} \geq t_{ik} + c_{i0} - M(1 - x_{i0k}) \quad \forall i \in V, \forall k \in K$
Time window constraints for customer arrivals:
$a_i \leq t_{ik} \leq b_i \quad \forall i \in V, \forall k \in K$

7. Subtour elimination (MTZ) constraints

MTZ constraints for subtour elimination: $uik≥1∀i∈V,∀k∈Ku_{ik} \geq 1 \quad \forall i \in V, \forall k \in K$ $uik≤n∀i∈V,∀k∈Ku_{ik} \leq n \quad \forall i \in V, \forall k \in K$ $uik−ujk+nxijk≤n−1∀i,j∈V,i≠j,∀k∈Ku_{ik} - u_{jk} + n x_{ijk} \leq n - 1 \quad \forall i, j \in V, i \neq j, \forall k \in K$

8. Ensure that if a customer is not visited by vehicle $k$ , its time variables are inactive

y_{ik} \leq \sum_{j \in N, j \neq i} x_{jik} \quad \forall i \in V, \forall k \in K

Lazy Constraint Generation for Subtour Elimination

Here is the mathematical formulation, I implemented in this code.

from pyscipopt import Model, quicksum, Conshdlr, SCIP_RESULT
import networkx as nx

def create_data_model():
    """Stores the data for the problem."""
    data = {}
    # Time matrix representing travel times between locations (consistent with time windows)
    data["time_matrix"] = [
        [0, 6, 9, 8, 7, 3, 6, 2, 3, 2, 6, 6, 4, 4, 5, 9, 7],
        [6, 0, 8, 3, 2, 6, 8, 4, 8, 8, 13, 7, 5, 8, 12, 10, 14],
        [9, 8, 0, 11, 10, 6, 3, 9, 5, 8, 4, 15, 14, 13, 9, 18, 9],
        [8, 3, 11, 0, 1, 7, 10, 6, 10, 10, 14, 6, 7, 9, 14, 6, 16],
        [7, 2, 10, 1, 0, 6, 9, 4, 8, 9, 13, 4, 6, 8, 12, 8, 14],
        [3, 6, 6, 7, 6, 0, 2, 3, 2, 2, 7, 9, 7, 7, 6, 12, 8],
        [6, 8, 3, 10, 9, 2, 0, 6, 2, 5, 4, 12, 10, 10, 6, 15, 5],
        [2, 4, 9, 6, 4, 3, 6, 0, 4, 4, 8, 5, 4, 3, 7, 8, 10],
        [3, 8, 5, 10, 8, 2, 2, 4, 0, 3, 4, 9, 8, 7, 3, 13, 6],
        [2, 8, 8, 10, 9, 2, 5, 4, 3, 0, 4, 6, 5, 4, 3, 9, 5],
        [6, 13, 4, 14, 13, 7, 4, 8, 4, 4, 0, 10, 9, 8, 4, 13, 4],
        [6, 7, 15, 6, 4, 9, 12, 5, 9, 6, 10, 0, 1, 3, 7, 3, 10],
        [4, 5, 14, 7, 6, 7, 10, 4, 8, 5, 9, 1, 0, 2, 6, 4, 8],
        [4, 8, 13, 9, 8, 7, 10, 3, 7, 4, 8, 3, 2, 0, 4, 5, 6],
        [5, 12, 9, 14, 12, 6, 6, 7, 3, 3, 4, 7, 6, 4, 0, 9, 2],
        [9, 10, 18, 6, 8, 12, 15, 8, 13, 9, 13, 3, 4, 5, 9, 0, 9],
        [7, 14, 9, 16, 14, 8, 5, 10, 6, 5, 4, 10, 8, 6, 2, 9, 0],
    ]
    data["time_windows"] = [
        (0, 5),     # 0 - depot
        (7, 12),    # 1
        (10, 15),   # 2
        (16, 18),   # 3
        (10, 13),   # 4
        (0, 5),     # 5
        (5, 10),    # 6
        (0, 4),     # 7
        (5, 10),    # 8
        (0, 3),     # 9
        (10, 16),   # 10
        (10, 15),   # 11
        (0, 5),     # 12
        (5, 10),    # 13
        (7, 8),     # 14
        (10, 15),   # 15
        (11, 15),   # 16
    ]
    data["demands"] = [0, 1, 1, 2, 4, 2, 4, 8, 8, 1, 2, 1, 2, 4, 4, 8, 8]
    data["vehicle_capacities"] = [15, 15, 15, 15]  # Capacities for each vehicle
    data["num_vehicles"] = 4
    data["depot"] = 0
    return data

def main():
    """Solve the VRPTW problem using PySCIPOpt with subtour elimination as lazy constraints."""
    # Retrieve data
    data = create_data_model()
    time_matrix = data["time_matrix"]
    time_windows = data["time_windows"]
    demands = data["demands"]
    vehicle_capacities = data["vehicle_capacities"]
    num_vehicles = data["num_vehicles"]
    depot = data["depot"]

    n = len(time_matrix)  # Number of nodes
    N = range(n)  # Set of nodes
    V = [i for i in N if i != depot]  # Customers excluding depot
    K = range(num_vehicles)  # Set of vehicles

    # Calculate M for Big-M constraints (maximum possible time)
    max_time = max(tw[1] for tw in time_windows) + max(
        time_matrix[i][j] for i in N for j in N
    )

    # Initialize the model
    model = Model("VRPTW with Lazy Constraints")

    # Variables
    x = {}  # Binary variables: x[i,j,k] == 1 if vehicle k travels from i to j
    t = {}  # Arrival time at node i by vehicle k
    t_start = {}  # Departure time from depot for vehicle k
    t_end = {}    # Return time to depot for vehicle k
    y = {}  # Binary variables: y[i,k] == 1 if vehicle k visits node i

    # Decision variables
    for k in K:
        capacity_k = vehicle_capacities[k]
        # Departure time from depot
        t_start[k] = model.addVar(lb=time_windows[depot][0], ub=time_windows[depot][1], vtype="I", name=f"t_start_{k}")
        # Return time to depot
        t_end[k] = model.addVar(lb=time_windows[depot][0], ub=None, vtype="I", name=f"t_end_{k}")
        for i in N:
            if i != depot:
                y[i, k] = model.addVar(vtype="B", name=f"y_{i}_{k}")
                t[i, k] = model.addVar(lb=time_windows[i][0], ub=time_windows[i][1], vtype="I", name=f"t_{i}_{k}")
            for j in N:
                if i != j:
                    x[i, j, k] = model.addVar(vtype="B", name=f"x_{i}_{j}_{k}")

    # Objective Function: Minimize total operation time of all vehicles
    model.setObjective(
        quicksum(t_end[k] - t_start[k] for k in K),
        "minimize",
    )

    # Constraints

    # 1. Each customer is visited exactly once
    for i in V:
        model.addCons(
            quicksum(y[i, k] for k in K) == 1,
            name=f"VisitOnce_{i}"
        )

    # 2. For each vehicle, ensure it departs from and returns to the depot correctly
    for k in K:
        model.addCons(
            quicksum(x[depot, j, k] for j in N if j != depot) == 1,
            name=f"DepartDepot_{k}"
        )
        model.addCons(
            quicksum(x[i, depot, k] for i in N if i != depot) == 1,
            name=f"ReturnDepot_{k}"
        )

    # 3. Flow conservation constraints for nodes (excluding depot)
    for k in K:
        for i in V:
            model.addCons(
                quicksum(x[i, j, k] for j in N if j != i) - quicksum(x[j, i, k] for j in N if j != i) == 0,
                name=f"FlowConservation_{i}_{k}"
            )

    # 4. Link x and y variables
    for k in K:
        for i in V:
            model.addCons(
                quicksum(x[i, j, k] for j in N if j != i) == y[i, k],
                name=f"LinkXY_{i}_{k}"
            )

    # 5. Capacity constraints for each vehicle
    for k in K:
        capacity_k = vehicle_capacities[k]
        model.addCons(
            quicksum(demands[i] * y[i, k] for i in V) <= capacity_k,
            name=f"VehicleCapacity_{k}"
        )

    # 6. Time window constraints and time propagation constraints
    for k in K:
        # Departure time from depot is within time window
        model.addCons(
            t_start[k] >= time_windows[depot][0],
            name=f"DepotTimeWindowStart_{k}"
        )
        model.addCons(
            t_start[k] <= time_windows[depot][1],
            name=f"DepotTimeWindowEnd_{k}"
        )

        # Time propagation constraints
        for i in N:
            if i != depot:
                # Arrival time at node i for vehicle k
                model.addCons(
                    t[i, k] >= time_windows[i][0],
                    name=f"TimeWindowStart_{i}_{k}"
                )
                model.addCons(
                    t[i, k] <= time_windows[i][1],
                    name=f"TimeWindowEnd_{i}_{k}"
                )
            for j in N:
                if i != j:
                    travel_time = time_matrix[i][j]
                    if i == depot:
                        # From depot to first customer
                        model.addCons(
                            t[j, k] >= t_start[k] + travel_time - max_time * (1 - x[i, j, k]),
                            name=f"TimePropagationDepot_{i}_{j}_{k}"
                        )
                    elif j == depot:
                        # From last customer to depot
                        model.addCons(
                            t_end[k] >= t[i, k] + travel_time - max_time * (1 - x[i, j, k]),
                            name=f"TimePropagationReturnDepot_{i}_{j}_{k}"
                        )
                    else:
                        # Between customers
                        model.addCons(
                            t[j, k] >= t[i, k] + travel_time - max_time * (1 - x[i, j, k]),
                            name=f"TimePropagation_{i}_{j}_{k}"
                        )

    # Include the subtour elimination constraint handler
    class SubtourEliminationConshdlr(Conshdlr):
        def __init__(self, x, n, N, V, K, depot):
            super().__init__()
            self.x = x
            self.n = n
            self.N = N
            self.V = V
            self.K = K
            self.depot = depot

        def conscheck(self, constraints, solution, checkintegrality, checklprows, printreason, completely):
            # For integer solutions, we ensure no subtours exist per vehicle
            for k in self.K:
                selected_edges = [(i, j) for i in self.N for j in self.N if i != j and solution[self.x[i, j, k]] > 0.5]
                G = nx.DiGraph()
                G.add_edges_from(selected_edges)
                # Find strongly connected components (subtours)
                subtours = list(nx.strongly_connected_components(G))
                for subtour in subtours:
                    if self.depot in subtour:
                        continue  # Skip component containing depot
                    if len(subtour) < 2:
                        continue
                    # Subtour detected, infeasible
                    if printreason:
                        print(f"Subtour detected for vehicle {k}: {subtour}")
                    return {"result": SCIP_RESULT.INFEASIBLE}

            return {"result": SCIP_RESULT.FEASIBLE}

        def consenfolp(self, constraints, nusefulconss, solinfeasible):
            # Enforce LP solution
            violated = False
            for k in self.K:
                selected_edges = [(i, j) for i in self.N for j in self.N if i != j and self.model.getVal(self.x[i, j, k]) > 0.5]
                G = nx.DiGraph()
                G.add_edges_from(selected_edges)
                # Find strongly connected components (subtours)
                subtours = list(nx.strongly_connected_components(G))
                for subtour in subtours:
                    if self.depot in subtour:
                        continue  # Skip component containing depot
                    if len(subtour) < 2:
                        continue
                    # Subtour detected, add subtour elimination constraint
                    subtour_nodes = list(subtour)
                    expr = quicksum(self.x[i, j, k] for i in subtour_nodes for j in subtour_nodes if i != j)
                    self.model.addCons(expr <= len(subtour_nodes) - 1)
                    violated = True
            if violated:
                return {"result": SCIP_RESULT.CONSADDED}
            else:
                return {"result": SCIP_RESULT.FEASIBLE}

        def consenfeas(self, constraints, solution, objinfeasible):
            # Enforce constraints on feasible integer solutions
            violated = False
            for k in self.K:
                selected_edges = [(i, j) for i in self.N for j in self.N if i != j and self.model.getSolVal(solution, self.x[i, j, k]) > 0.5]
                G = nx.DiGraph()
                G.add_edges_from(selected_edges)
                # Find strongly connected components (subtours)
                subtours = list(nx.strongly_connected_components(G))
                for subtour in subtours:
                    if self.depot in subtour:
                        continue  # Skip component containing depot
                    if len(subtour) < 2:
                        continue
                    # Subtour detected, add subtour elimination constraint
                    subtour_nodes = list(subtour)
                    expr = quicksum(self.x[i, j, k] for i in subtour_nodes for j in subtour_nodes if i != j)
                    self.model.addCons(expr <= len(subtour_nodes) - 1)
                    violated = True
            if violated:
                return {"result": SCIP_RESULT.CONSADDED}
            else:
                return {"result": SCIP_RESULT.FEASIBLE}

    # Include the constraint handler
    conshdlr = SubtourEliminationConshdlr(x, n, N, V, K, depot)
    model.includeConshdlr(conshdlr, "SubtourEliminationConshdlr", "Eliminates subtours in VRPTW")

    # Solve the model
    model.optimize()

    # Retrieve and print the solution
    status = model.getStatus()
    if status == "optimal" or status == "bestsollimit":
        print_solution(data, x, t_start, t_end, t, y, K, N, depot, model)
    else:
        print("No optimal solution found.")
        print(f"Solver status: {status}")

def print_solution(data, x, t_start, t_end, t, y, K, N, depot, model):
    """Prints solution in the specified format."""
    solution = model.getBestSol()
    total_time = 0

    print(f"Objective: {model.getSolObjVal(solution):.1f}")
    for k in K:
        vehicle_id = k
        # Build the route for vehicle k
        edges = [(i, j) for i in N for j in N if i != j and model.getSolVal(solution, x.get((i, j, k), 0)) > 0.5]
        if not edges:
            continue  # Skip if vehicle k is not used
        G = nx.DiGraph()
        G.add_edges_from(edges)

        if depot not in G.nodes:
            continue  # Skip if vehicle k is not used

        plan_output = f"Route for vehicle {vehicle_id}:\n"

        route = [depot]
        time_log = [f"Time({model.getSolVal(solution, t_start[k]):.1f})"]
        next_node = depot

        while True:
            successors = list(G.successors(next_node))
            if not successors:
                break
            next_node = successors[0]
            if next_node == depot:
                break
            route.append(next_node)
            arrival_time = model.getSolVal(solution, t[next_node, k])
            time_log.append(f"Time({arrival_time:.1f})")

        # Append depot at the end
        route.append(depot)
        time_log.append(f"Time({model.getSolVal(solution, t_end[k]):.1f})")

        # Print the route and times
        for idx, node in enumerate(route):
            plan_output += f"{node} {time_log[idx]} -> "
        plan_output = plan_output.rstrip(" -> ") + "\n"
        route_duration = model.getSolVal(solution, t_end[k]) - model.getSolVal(solution, t_start[k])
        plan_output += f"Time of the route: {route_duration:.1f}min\n"
        print(plan_output)
        total_time += route_duration

    print(f"Total time of all routes: {total_time:.1f}min")

if __name__ == "__main__":
    main()

Sets and Indices

$N$ : Set of all nodes (including depot)
$V$ : Set of customer nodes (excluding depot)
$K$ : Set of vehicles
$\in N$ : Indices for nodes
$\in K$ : Index for vehicles

Parameters

$c_{ij}$ : Travel time from node $i$ to node $j$
$a_i, b_i]$ : Time window for node $i$
$d_i$ : Demand at node $i$
$Q_k$ : Capacity of vehicle $k$
$M$ : A large constant (Big-M)

Decision Variables

$x_{ijk}$ : Binary variable, 1 if vehicle $k$ travels from node $i$ to node $j$ , 0 otherwise
$t_{ik}$ : Arrival time at node $i$ by vehicle $k$
$t_{start_k}$ : Departure time from depot for vehicle $k$
$t_{end_k}$ : Return time to depot for vehicle $k$
$y_{ik}$ : Binary variable, 1 if vehicle $k$ visits node $i$ , 0 otherwise

Objective Function

Minimize the total operation time of all vehicles:

\text{Minimize} \quad \sum_{k \in K} (t_{end_k} - t_{start_k})

Constraints

1. Each customer is visited exactly once

\sum_{k \in K} y_{ik} = 1 \quad \forall i \in V

2. Each vehicle departs from and returns to the depot exactly once

\sum_{j \in N, j \neq i} x_{0jk} = 1 \quad \forall k \in K

\sum_{i \in N, i \neq j} x_{i0k} = 1 \quad \forall k \in K

3. Flow conservation constraints for nodes (excluding depot)

\sum_{j \in N, j \neq i} x_{ijk} - \sum_{j \in N, j \neq i} x_{jik} = 0 \quad \forall i \in V, \forall k \in K

4. Link $x$ and $y$ variables

\sum_{j \in N, j \neq i} x_{ijk} = y_{ik} \quad \forall i \in V, \forall k \in K

5. Capacity constraints for each vehicle

\sum_{i \in V} d_i y_{ik} \leq Q_k \quad \forall k \in K

6. Time window constraints and time propagation constraints

Departure time from depot is within time window:

$a_0 \leq t_{start_k} \leq b_0 \quad \forall k \in K$
Time propagation constraints:
- From depot to first customer: $tjk≥tstartk+c0j−M(1−x0jk)∀j∈V,∀k∈Kt_{jk} \geq t_{start_k} + c_{0j} - M(1 - x_{0jk}) \quad \forall j \in V, \forall k \in K$
- Between customers: $tjk≥tik+cij−M(1−xijk)∀i,j∈V,i≠j,∀k∈Kt_{jk} \geq t_{ik} + c_{ij} - M(1 - x_{ijk}) \quad \forall i, j \in V, i \neq j, \forall k \in K$
- Return to depot: $tendk≥tik+ci0−M(1−xi0k)∀i∈V,∀k∈Kt_{end_k} \geq t_{ik} + c_{i0} - M(1 - x_{i0k}) \quad \forall i \in V, \forall k \in K$
Time window constraints for customer arrivals:
$a_i \leq t_{ik} \leq b_i \quad \forall i \in V, \forall k \in K$

7. Subtour elimination constraints (Lazy Constraint Generation)

Subtour elimination constraints are dynamically added during the optimization process. For any detected subtour $S$ : $∑i∈S∑j∈S,j≠ixijk≤∣S∣−1∀k∈K\sum_{i \in S} \sum_{j \in S, j \neq i} x_{ijk} \leq |S| - 1 \quad \forall k \in K$

Dynamic Addition of Constraints

In this model, subtour elimination constraints are not added upfront but are dynamically generated during the optimization process. This approach is known as Lazy Constraint Generation. During the optimization, the solver periodically checks the current solution for any subtours. If a subtour is detected, a new constraint is added to the model to eliminate that subtour. This process continues until no more subtours are found, ensuring that the final solution is feasible and optimal.

Performance Comparison

Below is a table comparing the performance of different models and tools used to solve the Capacitated Vehicle Routing Problem with Time Windows (CVRPTW). The comparison includes execution time, lines of code, and the total time of all routes.

Model/Tool	Execution Time in secs	Lines of Code	Total Time of All Routes in mins (Objective Function Value)
OR-Tools	0.02	170	82
Node-based Model (SCIP)	0.54	283	81
Flow-based Model (SCIP)	21	294	81
MTZ Model (SCIP)	11	288	81
Lazy Constraint Generation for Subtour Elimination (SCIP)	4	339	81

Description

OR-Tools: The fastest among all, with an execution time of 0.02 seconds and 170 lines of code. However, the total time of all routes is slightly higher at 82 minutes.
Node-based Model (SCIP): Provides a slightly better solution with a total route time of 81 minutes and executes in 0.54 seconds. The implementation requires 283 lines of code.
Flow-based Model (SCIP): Takes significantly longer to execute at 21 seconds but achieves the same total route time of 81 minutes. The implementation consists of 294 lines of code.
MTZ Model (SCIP): Executes in 11 seconds with a total route time of 81 minutes and requires 288 lines of code.
Lazy Constraint Generation for Subtour Elimination (SCIP): Balances execution time and code complexity, taking 4 seconds to execute with 339 lines of code, and achieves a total route time of 81 minutes.

This comparison highlights the trade-offs between execution time, code complexity, and solution quality. OR-Tools is the fastest but slightly less optimal in terms of total route time. SCIP models, particularly the Node-based model, provide better solutions but at the cost of increased execution time and code complexity.

Personal Experience and Insights

Ease of Use and Code Complexity

OR-Tools

Fewer Lines of Code: One of the standout features of OR-Tools is its simplicity and ease of use. The API is designed to be user-friendly, which allows for rapid development and fewer lines of code. This makes OR-Tools an excellent choice for proof-of-concept (PoC) projects or quick benchmarks.
Good for PoC or Benchmark: Due to its simplicity, OR-Tools is ideal for quickly setting up and testing optimization problems. It allows you to focus on the problem at hand without getting bogged down by the intricacies of the underlying algorithms.

SCIP

More Control and Customization: SCIP offers a high degree of control and customization, allowing you to fine-tune the optimization process to meet specific requirements. This is particularly useful for complex problems where standard solvers may not perform well.
More Complex: The trade-off for this level of control is increased complexity. SCIP requires a deeper understanding of optimization techniques and more lines of code to implement. This can be a barrier for beginners but is invaluable for advanced users who need to customize their solutions.

Reliability and Support

OR-Tools

Risk of Google Discontinuing Projects: One of the potential drawbacks of using OR-Tools is the risk associated with Google's history of discontinuing projects. While OR-Tools is currently well-supported and actively maintained, there is always a possibility that Google may decide to discontinue the project in the future. This uncertainty can be a concern for users who rely on OR-Tools for critical applications.
Bad Documentation: Another challenge with OR-Tools is the quality of its documentation. While the official documentation provides a good starting point, it can be lacking in detail and clarity, making it difficult for users to fully understand and utilize all the features and capabilities of the tool. This can be particularly frustrating for new users who are trying to learn how to use OR-Tools effectively.
Bad Choice for Deployment: Given the potential risks and challenges associated with OR-Tools, it may not be the best choice for deployment in production environments. The uncertainty around long-term support and the limitations of the documentation can make it difficult to ensure the reliability and maintainability of applications built with OR-Tools.

SCIP

Good Choice for Deployment: In contrast, SCIP is a more reliable choice for deployment in production environments. As an open-source project with a strong academic and industrial user base, SCIP is less likely to be discontinued or lose support. Additionally, the comprehensive documentation and active community make it easier for users to get help and find solutions to any issues they may encounter.
Better Documentation: SCIP provides extensive and detailed documentation, including examples and tutorials, which can help users understand and effectively utilize the tool. This makes it easier to implement and maintain optimization models, particularly for complex and specialized problems.
Long-Term Support: The strong academic and industrial backing of SCIP ensures that it will continue to be supported and developed for the foreseeable future. This makes it a more reliable choice for long-term projects and critical applications where stability and support are essential.

Conclusion

In summary, while OR-Tools offers simplicity and ease of use, making it a good choice for quick proofs of concept and benchmarking, it may not be the best option for deployment due to potential risks and documentation challenges. On the other hand, SCIP provides greater control and customization, along with better documentation and long-term support, making it a more reliable choice for deployment in production environments.

Operations Research: How to Allow Users to Tune Model when an Instance Causes Unsolvable result

Thana B. — Sun, 06 Oct 2024 03:02:10 +0000

Goal Programming: Balancing Multiple Objectives with Flexibility

The Importance of Goal Programming

Benefits of Goal Programming

Flexibility: Accommodates soft constraints, reflecting real-world conditions where not all requirements are absolute.
Feasibility: Transforms infeasible models into solvable ones by allowing deviations from less critical constraints.
Prioritization: Enables ranking objectives according to their importance, ensuring critical goals are addressed first.
Tunability: Allows users to adjust goal priorities without changing the original problem instance, providing a powerful tool for model refinement.

How to Formulate Goal Constraints

Steps to Formulate Goal Constraints

Identify the Goals: List all desired outcomes and soft constraints.
Write the Requirement as a Regular Constraint: Formulate each goal as a standard constraint (equality or inequality).
Reformulate as a Goal Constraint: Introduce deviational variables to measure underachievement and overachievement.

The general method for formulating goal constraints is summarized in the table below:

Desired Situation	Formulation of Goal Constraint	Contribution to Objective Function
$LHS≤RHS\text{LHS} \leq \text{RHS}$	$LHS+d−−d+=RHS\text{LHS} + d^- - d^+ = \text{RHS}$	Minimize $d^+$
$LHS=RHS\text{LHS} = \text{RHS}$	$LHS+d−−d+=RHS\text{LHS} + d^- - d^+ = \text{RHS}$	Minimize $d^- + d^+$
$LHS≥RHS\text{LHS} \geq \text{RHS}$	$LHS+d−−d+=RHS\text{LHS} + d^- - d^+ = \text{RHS}$	Minimize $d^-$

Establish the Objective Function: Create an objective function that minimizes the weighted sum of the deviational variables, reflecting the priorities of each goal.

Explanation of Deviational Variables

Negative Deviational Variable ( $d^-$ ): Measures the amount by which a goal is underachieved.
Positive Deviational Variable ( $d^+$ ): Measures the amount by which a goal is overachieved.

Example: Allocating Diamonds to Stores

Problem Statement

Absolutely Necessary (Hard Constraints):

(a) Allocate between 1,000 and 1,200 carats in total to the five stores.
(b) Store 5 must receive at least 300 carats of diamonds.

Desired Properties of the Allocation (Soft Constraints and Objectives):

(c) The stores in the malls (Stores 1, 2, and 3) should receive at least 80% of all the diamonds, if possible.
(d) The allocations to the stores in the malls should be equal to each other, if possible.
(e) The probabilities of theft in the stores have been estimated to be 0.1%, 0.1%, 9%, 2%, and 3% for Stores 1 to 5, respectively. The owner would like to minimize the expected loss due to theft.

Goal Priorities:

Priority 1: requirement (e) is 25 times as important as requirement (d)
Priority 2: requirement (d) is 2 times as important as requirement (c)
Priority 3: requirement (c) as a baseline of priority. Ensure mall stores receive at least 80% of the diamonds.

Traditional Linear Programming Model

In the traditional model, all constraints, including the soft constraints and objectives, are treated as hard constraints.

Decision Variables

$x_i$ : Number of carats allocated to Store $i$ (for $i = 1, 2, 3, 4, 5$ ).

Objective Function

Minimize Expected Loss (e): $Z=0.001x1+0.001x2+0.09x3+0.02x4+0.03x5\text{Minimize } Z = 0.001 x_1 + 0.001 x_2 + 0.09 x_3 + 0.02 x_4 + 0.03 x_5$

Constraints

Hard Constraints (Absolutely Necessary):

Total Allocation Bounds: $\leq x_1 + x_2 + x_3 + x_4 + x_5 \leq 1200$
Minimum Allocation to Store 5: $x5≥300x_5 \geq 300$

Soft Constraints (Treated as Hard Constraints):

Mall Stores Allocation (c):

$x_1 + x_2 + x_3 \geq 0.8(x_1 + x_2 + x_3 + x_4 + x_5)$
Equal Allocation Among Mall Stores (d):

$x_1 = x_2 = x_3$

Shortcomings of the Traditional Model

Infeasibility: Due to the rigidity of hard constraints and conflicting objectives, the model may become infeasible, offering no solution.
Lack of Flexibility: Cannot prioritize goals or allow for deviations, leaving the decision-maker without actionable solutions.

Goal Programming Approach

Goal programming introduces flexibility by converting objectives and soft constraints into goal constraints with deviational variables, allowing for prioritization and balancing of goals.

Step 1: Convert Soft Constraints and Objectives into Goal Constraints

(c) Mall Stores Allocation Goal

Original Constraint: $x1+x2+x3≥0.8(x1+x2+x3+x4+x5)x_1 + x_2 + x_3 \geq 0.8(x_1 + x_2 + x_3 + x_4 + x_5)$
Rewriting as Goal Constraint: $x_1 + x_2 + x_3 + d_1^- - d_1^+ = 0.8(x_1 + x_2 + x_3 + x_4 + x_5)$
Simplify: $0.2(x_1 + x_2 + x_3) - 0.8 x_4 - 0.8 x_5 + d_1^- - d_1^+ = 0$
Contribution to Objective Function: Minimize $d_1^-$ (underachievement).

(d) Equal Allocation Among Mall Stores

For Store 1 and Store 2:

Original Constraint: $x_1 = x_2$
Rewriting as Goal Constraint: $x_1 - x_2 + d_2^- - d_2^+ = 0$
Contribution to Objective Function: Minimize $d_2^- + d_2^+$ (total deviation).

For Store 1 and Store 3:

Original Constraint: $x_1 = x_3$
Rewriting as Goal Constraint: $x_1 - x_3 + d_3^- - d_3^+ = 0$
Contribution to Objective Function: Minimize $d_3^- + d_3^+$ .

(e) Minimize Expected Loss Due to Theft

Original Objective Function: $Z=0.001x1+0.001x2+0.09x3+0.02x4+0.03x5\text{Minimize } Z = 0.001 x_1 + 0.001 x_2 + 0.09 x_3 + 0.02 x_4 + 0.03 x_5$
Setting a Target (Aspiration Level): Aim for minimal expected loss (e.g., $Z = 0$ , recognizing it may not be attainable).
Rewriting as Goal Constraint: $0.001 x_1 + 0.001 x_2 + 0.09 x_3 + 0.02 x_4 + 0.03 x_5 + d_5^- - d_5^+ = 0$
Contribution to Objective Function: Minimize $d_5^+$ (overachievement).

Step 2: Assign Weights Based on Goal Priorities

Priority 1 (Minimize Expected Loss): Weight $W_5 = 50$ (25 times more important than Goal 2).
Priority 2 (Equal Allocation): Weight $W_2 = W_3 = 2$ (2 times more important than Goal 3).
Priority 3 (Mall Stores Allocation): Weight $W_1 = 1$ .

Step 3: Formulate the Goal Programming Mathematical Model

Decision Variables

$xi≥0x_i \geq 0$ : Allocation to Store $i$ .
$dj−≥0,dj+≥0d_j^- \geq 0, \quad d_j^+ \geq 0$ : Deviational variables for Goal $j$ .

Objective Function

Minimize the weighted sum of deviational variables:

\text{Minimize } Z = 50 d_5^+ + 2 (d_2^- + d_2^+ + d_3^- + d_3^+) + d_1^-

Constraints

Hard Constraints (Absolutely Necessary):

Total Allocation Bounds: $\leq x_1 + x_2 + x_3 + x_4 + x_5 \leq 1200$
Minimum Allocation to Store 5: $x5≥300x_5 \geq 300$

Goal Constraints:

Mall Stores Allocation (c): $0.2 x_1 + 0.2 x_2 + 0.2 x_3 - 0.8 x_4 - 0.8 x_5 + d_1^- - d_1^+ = 0$
Equal Allocation Between Stores 1 and 2 (d): $x_1 - x_2 + d_2^- - d_2^+ = 0$
Equal Allocation Between Stores 1 and 3 (d): $x_1 - x_3 + d_3^- - d_3^+ = 0$
Expected Loss Minimization (e): $0.001 x_1 + 0.001 x_2 + 0.09 x_3 + 0.02 x_4 + 0.03 x_5 + d_5^- - d_5^+ = 0$

Non-Negativity Constraints:

x_i \geq 0, \qquad d_j^- \geq 0, \qquad d_j^+ \geq 0

Step 4: Solve the Goal Programming Model

Using an appropriate optimization solver, we find the values of $x_i$ and $dj±d_j^\pm$ that minimize $Z$ .

Interpretation of Results

Assuming the optimal solution is:

Allocations:
- $x_1 = 466.69$ carats
- $x_2 = 233.31$ carats
- $x_3 = 0$ carats
- $x_4 = 0$ carats
- $x_5 = 300$ carats
Total Allocation:
$x_1 + x_2 + x_3 + x_4 + x_5 = 1000 \text{ carats (meeting the minimum requirement)}$
Deviational Variables:
- $valued_1^- = \text{Positive value}$ : Indicating underachievement of the mall stores allocation goal.
- $d_1^+ = 0$ .
- $valuesd_2^-, d_2^+, d_3^-, d_3^+ = \text{Positive values}$ : Indicating deviations from equal allocation among mall stores.
- $d_5^+ = 9.7$ : Overachievement of expected loss.
- $d_5^- = 0$ .

Advantages of Goal Programming in This Example

Feasibility: Provides a solution where the traditional model may be infeasible due to conflicting constraints.
Prioritization: Allows the decision-maker to focus on the most critical objectives.
Flexibility: Accommodates deviations from less critical goals to achieve the most important ones.
Tunability: Users can adjust weights to explore different trade-offs without altering the problem instance.

Conclusion

Key Takeaways:

Translates Soft Constraints into Goal Constraints: Allows for controlled deviations, reflecting real-world flexibility.
Balances Conflicting Objectives: Prioritizes goals according to their importance.
Maintains Feasibility: Provides solutions even when traditional models fail due to infeasibility.
Adjustable: Enables fine-tuning of models through weight adjustments without changing the problem itself.

Goal programming empowers decision-makers to find practical solutions in complex environments, aligning closely with their objectives and the realities of the constraints they face.

Reference:

(Springer Texts In Business And Economics) H. A. Eiselt, Carl-Louis Sandblom - Operations Research_ A Model-Based Approach-Springer (2022)

Enforcing Driver Breaks After Departure Time in Google OR-Tools VRP

Thana B. — Thu, 03 Oct 2024 00:18:32 +0000

In the world of logistics and transportation, planning efficient routes for vehicles is crucial. However, equally important is ensuring that drivers adhere to regulations, such as taking mandatory breaks after certain hours of driving. In this blog post, we'll explore how to enforce break constraints after a vehicle departs from the depot using Google OR-Tools in a Capacitated Vehicle Routing Problem with Time Windows and Breaks (CVRPTWB).

We'll walk through a Python example where we define an ActiveTime dimension that tracks the time since each vehicle departed from the depot. By leveraging this dimension, we'll set up mandatory breaks for drivers after specific hours of active driving.

Problem Overview

We are tackling a Capacitated Vehicle Routing Problem with Time Windows and Breaks (CVRPTWB). In this problem, we need to plan routes for vehicles that have:

Capacity constraints: Each vehicle has a maximum load capacity.
Time windows: Each customer must be served within a specific time frame.
Driver breaks: Drivers must take breaks after certain hours of driving since they departed from the depot.

Our goal is to minimize the total distance traveled while adhering to these constraints.

Setting Up the Environment

First, ensure you have Google OR-Tools installed:

pip install ortools

Now, let's dive into the code.

Code Walkthrough

Importing Necessary Libraries

from ortools.constraint_solver import pywrapcp
from ortools.constraint_solver import routing_enums_pb2
import random
import pandas as pd

Creating the Data Model

We define a function create_data_model() to generate our problem's data, including locations, demands, and time windows.

def create_data_model():
    """Stores the data for the problem."""
    data = {}
    ...
    return data

Key parameters:

Number of orders: data["num_orders"]
Number of vehicles: data["num_vehicles"]
Vehicle capacities: data["vehicle_capacities"]
Time windows: data["time_windows"]
Service time per demand unit: data["time_per_demand_unit"]

Note: We use random generation for locations, demands, and time windows to simulate real-world variability. Due to this randomness, sometimes the problem may not have a feasible solution. If that happens, you might need to run the code multiple times.

Initializing the Routing Model

We create the routing index manager and the routing model.

manager = pywrapcp.RoutingIndexManager(
    data["num_orders"] + 1, data["num_vehicles"], data["depot"]
)
routing = pywrapcp.RoutingModel(manager)
parameters = pywrapcp.DefaultRoutingSearchParameters()

parameters.first_solution_strategy = (
    routing_enums_pb2.FirstSolutionStrategy.PARALLEL_CHEAPEST_INSERTION
)

Defining Distance and Demand Callbacks

We define the Manhattan distance and demand callbacks to calculate distances between nodes and demands at each node.

def manhattan_distance(x1, y1, x2, y2):
    return abs(x1 - x2) + abs(y1 - y2)

def transit_callback(i, j):
    ...

def demand_callback(i):
    ...

We register these callbacks with the routing model.

transit_index = routing.RegisterTransitCallback(transit_callback)
routing.SetArcCostEvaluatorOfAllVehicles(transit_index)

demand_index = routing.RegisterUnaryTransitCallback(demand_callback)
routing.AddDimensionWithVehicleCapacity(
    demand_index, 0, data["vehicle_capacities"], True, "Capacity"
)

Setting Service Times and Time Windows

We calculate service times based on demands and add the time dimension with time windows.

service_times = [
    data["time_per_demand_unit"] * data["demands"][manager.IndexToNode(i)]
    for i in range(routing.Size())
]

def time_callback(i, j):
    ...

time_index = routing.RegisterTransitCallback(time_callback)
routing.AddDimension(time_index, data["horizon"], data["horizon"], False, "Time")
time_dimension = routing.GetDimensionOrDie("Time")

We then add time window constraints for each customer and the depot.

for order in range(1, manager.GetNumberOfNodes()):
    start, end = data["time_windows"][order]
    index = manager.NodeToIndex(order)
    time_dimension.CumulVar(index).SetRange(start, end)

for vehicle_id in range(data["num_vehicles"]):
    index = routing.Start(vehicle_id)
    time_dimension.CumulVar(index).SetRange(
        data["time_windows"][data["depot"]][0], data["time_windows"][data["depot"]][1]
    )

Key Part 1: Defining the ActiveTime Dimension

To enforce breaks after a certain amount of active driving time since departing from the depot, we introduce an ActiveTime dimension.

Defining the ActiveTime Callback

def active_time_callback(i, j):
    loc_i = data["locations"][manager.IndexToNode(i)]
    loc_j = data["locations"][manager.IndexToNode(j)]
    return (
        int(manhattan_distance(loc_i[0], loc_i[1], loc_j[0], loc_j[1]) / data["speed"])
        + service_times[manager.IndexToNode(i)]
    )

This callback calculates the transit time between nodes, including the service time at the current node.

Adding the ActiveTime Dimension

We register the active_time_callback and add the ActiveTime dimension. Crucially, we set fix_start_cumul_to_zero to True, ensuring that the active time for each vehicle starts from zero at the depot.

active_time_index = routing.RegisterTransitCallback(active_time_callback)
routing.AddDimensionWithVehicleTransits(
    [active_time_index]*data["num_vehicles"], 0, data["horizon"], True, "ActiveTime"
)
active_time_dimension = routing.GetDimensionOrDie("ActiveTime")

Explanation:

AddDimensionWithVehicleTransits allows us to add a dimension that can have different transit callbacks for each vehicle. In our case, we use the same callback for all vehicles.
The fourth argument True is for fix_start_cumul_to_zero, which fixes the cumulative value at the start of each vehicle's route to zero. This is essential for accurately tracking the active time since departing from the depot.

Minimizing Time Variables

We also add time minimization constraints to help the solver find better solutions.

for i in range(routing.Size()):
    routing.AddVariableMinimizedByFinalizer(time_dimension.CumulVar(i))
    routing.AddVariableMinimizedByFinalizer(active_time_dimension.CumulVar(i))
for j in range(data["num_vehicles"]):
    routing.AddVariableMinimizedByFinalizer(time_dimension.CumulVar(routing.Start(j)))
    routing.AddVariableMinimizedByFinalizer(time_dimension.CumulVar(routing.End(j)))
    routing.AddVariableMinimizedByFinalizer(active_time_dimension.CumulVar(routing.Start(j)))
    routing.AddVariableMinimizedByFinalizer(active_time_dimension.CumulVar(routing.End(j)))

Key Part 2: Adding Vehicle Breaks Based on ActiveTime

Now that we have our ActiveTime dimension, we can set up breaks that occur after certain hours of active driving.

Defining Break Data

We define the break intervals that we want to enforce. For example, drivers must take a 30-minute break between 6-8 hours and another between 14-16 hours of active driving.

break_data = [
    # [start_min, start_max, duration] (in seconds)
    [6 * 3600, 8 * 3600, 1 * 1800],  # Morning break
    [14 * 3600, 16 * 3600, 1 * 1800],  # Afternoon break
]

Creating and Enforcing Break Intervals

We create break intervals for each vehicle and enforce that they must be taken.

solver = routing.solver()
breaks = {}
for vehicle in range(data["num_vehicles"]):
    breaks[vehicle] = []
    for i, data_break in enumerate(break_data):
        break_interval = solver.FixedDurationIntervalVar(
            data_break[0],
            data_break[1],
            data_break[2],
            False,  # is_optional set to False
            f"Vehicle {vehicle}'s Break {i}",
        )
        breaks[vehicle].append(break_interval)

    # Apply break constraints
    solver.Add(breaks[vehicle][0].PerformedExpr() == 1)
    solver.Add(breaks[vehicle][1].PerformedExpr() == 1)

    index = routing.Start(vehicle)
    solver.Add(time_dimension.CumulVar(index) < breaks[vehicle][0].StartExpr())

Linking Breaks to ActiveTime Dimension

This is the crucial step: We set the breaks to be based on the ActiveTime dimension rather than the default time dimension.

# Set breaks based on ActiveTime dimension
active_time_dimension.SetBreakIntervalsOfVehicle(
    breaks[vehicle], vehicle, service_times
)

Explanation:

By calling SetBreakIntervalsOfVehicle with our active_time_dimension, we ensure that breaks are scheduled based on the cumulative active time since the vehicle departed from the depot.
The service_times parameter accounts for the time spent servicing customers.

Optional: Forcing Longer Working Hours for Testing

Note: In our example, to test the breaks effectively, we add constraints to ensure vehicles work longer hours. This is for experimentation purposes and can be removed in actual applications.

# Please delete this when you apply to your problem
# This is added for this experiment to ensure results with long hours of working so we can test our breaks
for vehicle in range(data["num_vehicles"]):
    index = routing.Start(vehicle)
    solver.Add(time_dimension.CumulVar(index) <= 1 * 3600)
    index = routing.End(vehicle)
    solver.Add(time_dimension.CumulVar(index) >= 23 * 3600)

Allowing Orders to Be Skipped with Penalties

To provide flexibility in routing, we add penalties for skipping orders. This helps the solver make decisions that balance total distance with the need to serve all customers.

kPenalty = 10000000  # Penalty for skipping an order
for order in range(1, manager.GetNumberOfNodes()):
    routing.AddDisjunction([manager.NodeToIndex(order)], kPenalty)

Solving the Problem

We set up the solver and solve the problem.

solution = routing.SolveWithParameters(parameters)

Extracting and Printing the Solution

Finally, we extract the solution and print out the routes, including the break schedules.

Printing the Break Plans

print("==== BREAK PLAN ====")
for vehicle_id in range(data["num_vehicles"]):
    vehicle_breaks = breaks[vehicle_id]
    for break_interval in vehicle_breaks:
        if solution.PerformedValue(break_interval):
            print(
                f"{break_interval.Name()} Start({break_interval.StartExpr()}) "
                f"Duration({break_interval.DurationExpr()}) "
                f"End({break_interval.EndExpr()})"
            )

Printing the Route Plans

print("\n\n==== ROUTE PLAN ====")
total_distance = 0
total_load = 0
total_time = 0
capacity_dimension = routing.GetDimensionOrDie("Capacity")
time_dimension = routing.GetDimensionOrDie("Time")

for vehicle_id in range(data["num_vehicles"]):
    index = routing.Start(vehicle_id)
    plan_output = f"Route for vehicle {vehicle_id}:\n"
    route_distance = 0
    route_load = 0
    while not routing.IsEnd(index):
        ...
    print(plan_output)
    total_distance += route_distance
    total_load += route_load
    total_time += solution.Max(time_var)

print(f"Total Distance of all routes: {total_distance} units")
print(f"Total Load of all routes: {total_load}")
print(f"Total Time of all routes: {round(total_time/3600, 3)} hours")

Full Code

For completeness, here's the full code with comments:

# Capacitated Vehicle Routing Problem with Time Windows and Breaks (CVRPTWB)

from ortools.constraint_solver import pywrapcp
from ortools.constraint_solver import routing_enums_pb2
import random
import pandas as pd

# Data Model
def create_data_model():
    """Stores the data for the problem."""
    data = {}
    # Problem parameters corresponding to the mathematical model
    data["num_orders"] = 20  # |C| = 20 customer nodes
    data["num_vehicles"] = 3  # |V| = 3 vehicles
    data["depot"] = 0  # Depot node
    data["vehicle_capacities"] = [40] * data["num_vehicles"]  # Q = 40 for each vehicle
    data["speed"] = 100  # Vehicle speed (used in time calculations)

    # Time parameters
    data["time_per_demand_unit"] = 100  # Service time per demand unit (s_i = 100 * d_i)
    data["horizon"] = 24 * 3600  # H = 24 hours in seconds

    num_locations = data["num_orders"] + 1  # Total nodes |N| = |C| + 1 (depot)

    # Random seed for reproducibility
    use_deterministic_seed = False
    if use_deterministic_seed:
        random.seed(0)
    else:
        random.seed(None)

    # Randomly generate locations (x_i, y_i for each node i)
    kXMax = 100000
    kYMax = 100000
    data["locations"] = [
        (random.randint(0, kXMax), random.randint(0, kYMax))
        for _ in range(num_locations)
    ]

    # Random demands (d_i for each customer i)
    data["demands"] = [0] + [random.randint(1, 10) for _ in range(data["num_orders"])]

    # Time windows ([a_i, b_i] for each node i)
    data["time_windows"] = []
    kTWDuration = 5 * 3600  # 5 hours in seconds
    for _ in range(num_locations):
        start = random.randint(0, data["horizon"] - kTWDuration)
        data["time_windows"].append((start, start + kTWDuration))

    return data


# Instantiate the data problem
data = create_data_model()

# Initialize routing and search parameters
manager = pywrapcp.RoutingIndexManager(
    data["num_orders"] + 1, data["num_vehicles"], data["depot"]
)
routing = pywrapcp.RoutingModel(manager)
parameters = pywrapcp.DefaultRoutingSearchParameters()

# Set first solution strategy (not in mathematical model, but affects solution quality)
parameters.first_solution_strategy = (
    routing_enums_pb2.FirstSolutionStrategy.PARALLEL_CHEAPEST_INSERTION
)


# Define the Manhattan distance function (c_ij in the model)
def manhattan_distance(x1, y1, x2, y2):
    return abs(x1 - x2) + abs(y1 - y2)


# Callback for distance calculation (implements c_ij)
def transit_callback(i, j):
    loc_i = data["locations"][manager.IndexToNode(i)]
    loc_j = data["locations"][manager.IndexToNode(j)]
    return manhattan_distance(loc_i[0], loc_i[1], loc_j[0], loc_j[1])


transit_index = routing.RegisterTransitCallback(transit_callback)
routing.SetArcCostEvaluatorOfAllVehicles(transit_index)


# Adding capacity dimension constraints (Constraint 3 in the model)
def demand_callback(i):
    return data["demands"][manager.IndexToNode(i)]


demand_index = routing.RegisterUnaryTransitCallback(demand_callback)
routing.AddDimensionWithVehicleCapacity(
    demand_index, 0, data["vehicle_capacities"], True, "Capacity"
)


# Set service times (s_i in the model)
service_times = [
    data["time_per_demand_unit"] * data["demands"][manager.IndexToNode(i)]
    for i in range(routing.Size())
]


def time_callback(i, j):
    loc_i = data["locations"][manager.IndexToNode(i)]
    loc_j = data["locations"][manager.IndexToNode(j)]
    return (
        int(manhattan_distance(loc_i[0], loc_i[1], loc_j[0], loc_j[1]) / data["speed"])
        + service_times[manager.IndexToNode(i)]
    )


time_index = routing.RegisterTransitCallback(time_callback)
routing.AddDimension(time_index, data["horizon"], data["horizon"], False, "Time")
time_dimension = routing.GetDimensionOrDie("Time")


# === Key Part 1: Defining ActiveTime Dimension === #
def active_time_callback(i, j):
    loc_i = data["locations"][manager.IndexToNode(i)]
    loc_j = data["locations"][manager.IndexToNode(j)]
    return (
        int(manhattan_distance(loc_i[0], loc_i[1], loc_j[0], loc_j[1]) / data["speed"])
        + service_times[manager.IndexToNode(i)]
    )

active_time_index = routing.RegisterTransitCallback(active_time_callback)
routing.AddDimensionWithVehicleTransits(
    [active_time_index]*data["num_vehicles"], 0, data["horizon"], True, "ActiveTime"
)
active_time_dimension = routing.GetDimensionOrDie("ActiveTime")
# ================================================ #


# Adding time windows (Constraint 4 in the model)
for order in range(1, manager.GetNumberOfNodes()):
    start, end = data["time_windows"][order]
    index = manager.NodeToIndex(order)
    time_dimension.CumulVar(index).SetRange(start, end)

for vehicle_id in range(data["num_vehicles"]):
    index = routing.Start(vehicle_id)
    active_time_dimension.CumulVar(index).SetValue(0)

# Add time window constraints for each vehicle start node.
depot_idx = data["depot"]
for vehicle_id in range(data["num_vehicles"]):
    index = routing.Start(vehicle_id)
    time_dimension.CumulVar(index).SetRange(
        data["time_windows"][depot_idx][0], data["time_windows"][depot_idx][1]
    )

# Minimizing time variables (part of the objective function)
for i in range(routing.Size()):
    routing.AddVariableMinimizedByFinalizer(time_dimension.CumulVar(i))
for j in range(data["num_vehicles"]):
    routing.AddVariableMinimizedByFinalizer(time_dimension.CumulVar(routing.Start(j)))
    routing.AddVariableMinimizedByFinalizer(time_dimension.CumulVar(routing.End(j)))

# Minimizing ActiveTime variables
for i in range(routing.Size()):
    routing.AddVariableMinimizedByFinalizer(active_time_dimension.CumulVar(i))
for j in range(data["num_vehicles"]):
    routing.AddVariableMinimizedByFinalizer(
        active_time_dimension.CumulVar(routing.Start(j))
    )
    routing.AddVariableMinimizedByFinalizer(
        active_time_dimension.CumulVar(routing.End(j))
    )


# === Key Part 2: Adding Vehicle Breaks === #
solver = routing.solver()


break_data = [
    # [start_min, start_max, duration] (in seconds)
    [6 * 3600, 8 * 3600, 1 * 1800],  # Morning break
    [14 * 3600, 16 * 3600, 1 * 1800],  # Afternoon break
]

# Optional: Forcing long working hours for testing
# Please delete this when you apply to your problem
# This is added for this experiment to ensure results with long hours of working so we can test our breaks
for vehicle in range(data["num_vehicles"]):
    index = routing.Start(vehicle)
    solver.Add(time_dimension.CumulVar(index) <= 1 * 3600)
    index = routing.End(vehicle)
    solver.Add(time_dimension.CumulVar(index) >= 23 * 3600)

breaks = {}
for vehicle in range(data["num_vehicles"]):
    # Declare a list of break intervals for each vehicle
    breaks[vehicle] = []
    for i, data_break in enumerate(break_data):
        # Declare each possible break interval for the current vehicle
        break_interval = solver.FixedDurationIntervalVar(
            data_break[0],
            data_break[1],
            data_break[2],
            False,  # is_optional set to False
            f"Vehicle {vehicle}'s Break {i}",
        )
        breaks[vehicle].append(break_interval)

    # Apply break constraints
    solver.Add(breaks[vehicle][0].PerformedExpr() == 1)
    solver.Add(breaks[vehicle][1].PerformedExpr() == 1)

    index = routing.Start(vehicle)
    solver.Add(time_dimension.CumulVar(index) < breaks[vehicle][0].StartExpr())

    # Set breaks based on ActiveTime dimension
    active_time_dimension.SetBreakIntervalsOfVehicle(
        breaks[vehicle], vehicle, service_times
    )
# ========================================== #


# Adding penalty costs to allow skipping orders
kPenalty = 10000000  # Penalty for skipping an order
for order in range(1, manager.GetNumberOfNodes()):
    routing.AddDisjunction([manager.NodeToIndex(order)], kPenalty)

# Solve the problem
solution = routing.SolveWithParameters(parameters)

print("==== BREAK PLAN ====")
for vehicle_id in range(data["num_vehicles"]):
    # Get the break intervals for the current vehicle
    vehicle_breaks = breaks[vehicle_id]
    # Check each break interval of the vehicle
    for break_interval in vehicle_breaks:
        # Find whether the break interval was performed or not from the solution
        if solution.PerformedValue(break_interval):
            print(
                f"{break_interval.Name()} Start({break_interval.StartExpr()}) "
                f"Duration({break_interval.DurationExpr()}) "
                f"End({break_interval.EndExpr()})"
            )


print("\n\n==== ROUTE PLAN ====")
total_distance = 0
total_load = 0
total_time = 0
capacity_dimension = routing.GetDimensionOrDie("Capacity")
time_dimension = routing.GetDimensionOrDie("Time")

for vehicle_id in range(data["num_vehicles"]):
    index = routing.Start(vehicle_id)
    plan_output = f"Route for vehicle {vehicle_id}:\n"
    route_distance = 0
    route_load = 0
    while not routing.IsEnd(index):
        node_index = manager.IndexToNode(index)
        load_var = capacity_dimension.CumulVar(index)
        time_var = time_dimension.CumulVar(index)
        route_load = solution.Value(load_var)
        time_min = solution.Min(time_var)
        time_max = solution.Max(time_var)

        plan_output += f" {node_index} Load({route_load}) Time({round(time_min/3600, 3)},{round(time_max/3600, 3)}) -> "

        previous_index = index
        index = solution.Value(routing.NextVar(index))
        route_distance += routing.GetArcCostForVehicle(
            previous_index, index, vehicle_id
        )

    node_index = manager.IndexToNode(index)
    time_var = time_dimension.CumulVar(index)
    time_min = solution.Min(time_var)
    time_max = solution.Max(time_var)
    load_var = capacity_dimension.CumulVar(index)
    route_load = solution.Value(load_var)

    plan_output += (
        f" {node_index} Time({round(time_min/3600, 3)},{round(time_max/3600, 3)})\n"
    )
    plan_output += f"Distance of the route: {route_distance} units\n"
    plan_output += f"Load of the route: {route_load}\n"
    print(plan_output)
    total_distance += route_distance
    total_load += route_load
    total_time += solution.Max(time_var)


print(f"Total Distance of all routes: {total_distance} units")
print(f"Total Load of all routes: {total_load}")
print(f"Total Time of all routes: {round(total_time/3600, 3)} hours")

Validate Result

I wrote additional code to confirm breaks as below:

import pandas as pd

def seconds_to_hhmm(seconds):
    # Convert seconds to HH:MM format
    hours = int(seconds // 3600)
    minutes = int((seconds % 3600) // 60)
    return f"{hours}:{minutes:02d}"


events = []
for vehicle_id in range(data["num_vehicles"]):
    index = routing.Start(vehicle_id)

    break_remark = ""
    is_start = True

    while not routing.IsEnd(index):
        next_index = solution.Value(routing.NextVar(index))
        node_index = manager.IndexToNode(index)
        next_node_index = manager.IndexToNode(next_index)

        time_var = time_dimension.CumulVar(index)
        arrival_time = solution.Min(time_var)  # in seconds

        if is_start:
            depot_arrival_time = arrival_time
            is_start = False

        time_var = time_dimension.CumulVar(next_index)
        actual_next_arrival_time = solution.Min(time_var)  # in seconds

        service_time = service_times[node_index]  # in seconds
        start_wait = arrival_time + service_time
        expected_next_arrival_time = arrival_time + time_callback(index, next_index)
        transit_duration = expected_next_arrival_time - service_time - arrival_time
        waiting_duration = actual_next_arrival_time - expected_next_arrival_time
        start_transit = start_wait + waiting_duration
        if routing.IsStart(index):
            events.append(
                (
                    vehicle_id,
                    node_index,
                    "start_route",
                    arrival_time,
                    seconds_to_hhmm(arrival_time),
                    0,
                    seconds_to_hhmm(arrival_time - depot_arrival_time),
                    "",
                )
            )
        else:
            events.append(
                (
                    vehicle_id,
                    node_index,
                    "arrival_time",
                    arrival_time,
                    seconds_to_hhmm(arrival_time),
                    0,
                    seconds_to_hhmm(arrival_time - depot_arrival_time),
                    "",
                )
            )

        events.append(
            (
                vehicle_id,
                node_index,
                "start_transit",
                start_transit,
                seconds_to_hhmm(start_transit),
                seconds_to_hhmm(transit_duration),
                seconds_to_hhmm(start_transit - depot_arrival_time),
                "",
            )
        )

        index = next_index
        if waiting_duration > 0:
            end_wait = start_wait + waiting_duration

            for break_id, break_info in enumerate(break_data):
                overlap_start = max(start_wait - depot_arrival_time, break_info[0])
                overlap_end = min(
                    end_wait - depot_arrival_time, break_info[1] + break_info[2]
                )
                overlap_duration = overlap_end - overlap_start

                if overlap_duration > break_info[2]:
                    this_remark = f"break {break_id} "
                    break_remark =  break_remark + this_remark

            events.append(
                (
                    vehicle_id,
                    next_node_index,
                    "start_wait (idle + break)",
                    start_wait,
                    seconds_to_hhmm(start_wait),
                    seconds_to_hhmm(waiting_duration),
                    seconds_to_hhmm(start_wait - depot_arrival_time),
                    break_remark,
                )
            )

        break_remark = ""
        if routing.IsEnd(next_index):
            events.append(
                (
                    vehicle_id,
                    next_node_index,
                    "end_route",
                    actual_next_arrival_time,
                    seconds_to_hhmm(actual_next_arrival_time),
                    0,
                    seconds_to_hhmm(actual_next_arrival_time - depot_arrival_time),
                    "",
                )
            )

df_result = pd.DataFrame(
    events,
    columns=[
        "vehicle_id",
        "node",
        "event",
        "timestamp",
        "time",
        "duration",
        "since_start",
        "remark",
    ],
).sort_values(["vehicle_id", "timestamp"])
df_result = df_result[
    [
        "vehicle_id",
        "node",
        "event",
        "timestamp",
        "time",
        "since_start",
        "duration",
        "remark",
    ]
].copy()

Here is an example output:

Conclusion

By introducing an ActiveTime dimension and linking driver breaks to this dimension, we can effectively enforce breaks after specific hours of driving since departing from the depot in Google OR-Tools. This approach ensures compliance with regulations and enhances driver safety.

Remember that due to random data generation, you might need to run the code multiple times to obtain a feasible solution. In real-world applications, you'd use actual data, which would provide consistent results.

DEV Community: Thana B.

Implementation of Optimization of Red Blood Cell Shape Using Spherical Harmonics Parameterization

Introduction

Mathematical Formulation for RBC Shape Optimization with Fixed Posture and Symmetry

1. Decision Variables

2. Objective Function

Bending Energy

Shear Energy

Coupling Energy

Viscous Energy

Thermal Fluctuation Energy

Hydrodynamic Forces

3. Constraints

3.1 Volume Constraint

3.2 Surface Area Constraint

3.3 Area Difference Constraint

3.4 Reduced Volume Constraint

3.5 Symmetry Constraints

Python Implementation Using CasADi

1. Import Necessary Libraries

2. Define Constants and Parameters

3. Set Up the Optimization Problem

Decision Variables:

4. Compute Legendre Polynomials

5. Shape Representation r(θ)r(\theta) r(θ)

6. Compute Geometric Quantities

First Derivative drdθ\frac{dr}{d\theta} dθdr​

Second Derivative d2rdθ2\frac{d^2r}{d\theta^2} dθ2d2r​

Metric Coefficient EE E

Area Element dAdA dA

Mean Curvature HH H

7. Define Energy Components

Bending Energy

Shear Energy

Coupling Energy

Total Energy

8. Define Constraints

Volume Constraint

Surface Area Constraint

Reduced Volume Constraint

Area Difference Constraint

Symmetry Constraints

Positivity Constraint for r(θ)r(\theta) r(θ)

Biconcavity Enforcement

9. Set Objective and Solve

10. Post-processing and Visualization

Compute Optimized Shape

Generate 3D Coordinates

Plot the Optimized Shape

11. Validate Results

Discussion

Mapping Mathematical Expressions to Code

Simplifications and Assumptions

Advantages of Using CasADi

Conclusion

Applying Monte Carlos Model Modification to Handle Infeasibility Internally in Stochastic Programming

1. Introducing Penalty Terms or Slack Variables

Concept Overview

Applying Penalty Terms to the Stochastic Knapsack Problem

Original Stochastic Knapsack Model

Modified Model with Penalty Terms

Model Interpretation

Implementation Steps

Example

Considerations When Using Penalty Terms

2. Using Chance Constraints

Concept Overview

Applying Chance Constraints to the Stochastic Knapsack Problem

Model Formulation

Reformulating the Chance Constraint

Integer Programming Formulation

Explanation

Example

Considerations When Using Chance Constraints

Comparing the Two Methods

Penalty Terms vs. Chance Constraints

Choosing Between Methods

Implementation Tips

General Advice

Testing and Validation

5. Shape Representation $r(θ)r(\theta)$

First Derivative $drdθ\frac{dr}{d\theta}$

Second Derivative $d2rdθ2\frac{d^2r}{d\theta^2}$

Metric Coefficient $E$

Area Element $d A$

Mean Curvature $H$

Positivity Constraint for $r(θ)r(\theta)$