VRP through constraint programming

Question

I am with a problem that consists of determining whether a VRP instance is feasible, this question is a continuation of an older thread. In order to do it, I was thinking of using Constraint Programming (CP). Since this is my first time using CP, I am facing some problems, specifically, in what concerns modeling the subcycle constraints in CP. Before I show the CP model, let's define the VRP notations to substantiate the model. Let

• $D(V, A)$ be a complete digraph;
• $V$ be the set of nodes;
• $0 \in V$ be the depot node;
• $A$ be the set of arcs;
• $t_a$ be the arc $a$ metric traversing time;
• $T$ be a vehicle's time limit; and
• $M$ be the set of available vehicles.

Note that we are dealing with the Distance Constrained VRP (DCVRP) variant. Since we just want to find a feasible solution, if such exists, there's no need for optimizing anything, hence model this problem as a decision problem suffices. Let

• $x_a \in \mathbb{B}$ be equals to $1$ if arc $a \in A$ is used in the solution, and $0$ otherwise;
• $\mu_i \in \mathbb{R}_{+}$ be equals to the amount of incurred time at node $i \in V\backslash\{0\}$;

Below is presented a CP model: $$ \sum_{a \in \delta^{+}(i)} x_a = \sum_{a \in \delta^{-}(i)} x_a \quad \forall i \in V \quad (1)\\ \sum_{a \in \delta^{+}(i)} x_a = 1 \quad \forall i \in V \backslash \{0\} \quad (2)\\ \sum_{a \in \delta^{+}(0)} x_a \leqslant |M| \quad (3)\\ u_i \geqslant t_{0i} x_{0i} \quad \forall i \in V \backslash \{0\} \quad (4)\\ u_j \geqslant u_i + t_{ij} x_{ij} - T (1 - x_{ij}) \quad \forall i, j \in V\backslash \{0\} \quad (5)\\ u_i \leqslant T - t_{i0} x_{i0} \quad \forall i \in V\backslash \{0\} \quad (6) $$

The constraints set (1) represents the degree constraints. The set of constraints (2) states that every customer must be visited exactly once. The set of constraints (3) forces that every solution use at most $|M|$ routes. And, the constraints sets (4-6) are the MTZ constraints to avoid subcycles disconnected from the depot, and at the same time guarantee that a route will consume at most $T$ of time.

I know I can model the constraints sets (1-3) easily, however, I can not state the same for the sets of constraints (4-6) since I do not know if it is possible to create non-discrete variables in CP models, as I said, this is my first time using CP. Therefore, I would like to know if it is possible to represent the constraints (4-6) in CP or it would require a workaround, and if yes what. Also, I would like to know if anyone has ever solved this problem, determine whether a VRP instance is feasible or not, through CP.

Many thanks for your time and regards.

UPDATE 1: I added the constraints set (2). And stated that the times $t_a$ $\forall a \in A$ are metric.

Answer 1

The model proposed by the OP is not a CP model, but a MIP model. Although in theory it would be possible to implement this model directly in CP, in practice it would perform terrible; you would be much better off implementing this model in a MIP solver, and strengthening it with valid inequalities (and use proper subtour separation routines instead of weak MTZ constraints).

There are at least 2 common CP formulations: one formulation uses 'next' variables, the other formulation uses scheduling constructs (e.g. interval variables). Here I'll only state the former formulation.

To simplify the CP formulation, let's reformulate the vertex set $V=\{0,n+1\}\cup V'$. Vertex $0$ is the starting depot, vertex $n+1$ is the ending depot (starting and ending depots can have the same physical location, but we just represent them as 2 different vertices). Vertices $V'=\{1,2,\dots,n\}$ are the customers that need to be visited.

For the CP model, we create $|M|$ copies of the destination depot. Let $D$ denote the set of destination depot copies. The model uses 2 types of variables:

• Variable $next_i$ denotes the next node visited immediately after node $i$
• Variable $S_i$ represents the arrival time at node $i$

Variables and domain definitions: \begin{align} & next_i\in V'\cup D &\forall i\in V'\\ & S_i \in [t_{0,i},T-t_{i,n+1}] &\forall i\in V'\\ & S_i \in [0,T] &\forall i\in D \end{align} Note that we define a $next_i$ variable for every customer, and a arrival time variable for every customer and destination depot. The notation $[x,y]$ defines the continuous time interval from $x$ to $y$ (inclusive).

Constraints:

\begin{align} &\textit{allDifferent}(next_i|i\in V') &\\ & S_{next_i} = S_i+t_{i,next_i} & \forall i\in V' \end{align} The allDifferent constraint is a global constraint that forces every decision variable to assume pairwise different values. In other words, this constraint ensures that for $i,j\in V'$ with $i\neq j$, it holds that $next_i\neq next_j$. Practically the allDifferent constraint ensures that each customer is visited at most once. In CP a feasible solution requires that each variable is assigned exactly one value from its domain, so this will ensure that every customer is visited exactly once.

The second constraint updates the arrival times. Note that in CP, we can index a variable with another variable! Again, this constraint ensures that vehicles return to the depot in time (before $T$), and that there are no subtours.

Note that the above formulation is very compact and actually easier that an equivalent ILP formulation. Unfortunately I don't think this CP formulation will work well to solve the problem that the OP intends to solve, namely verifying whether a single feasible solution exists. If a feasible solution exists, the CP solver might do a reasonable job at finding it. However, if no solution exists, the solver will likely do a rather poor job at proving that no such solution exists. If you have additional constraints, such as tight time windows on the customers, this would actually help the CP solver as it would narrow down the domains of the $S_i$ variables. Moreover, in your implementation you should probably set the 'inference level' of the allDifferent constraint to 'extended' (assuming the CP solver supports different inference levels for the allDifferent constraint)

Remarks:

1. The above CP model can be easily extended to support time-windows (i.e. VRPTW), by shrinking the domains of the $S_i$ variables to the corresponding time-windows. To allow a vehicle to arrive early at a customer and wait till the beginning of its time-window, the equality sign ($=$) in the 2nd constraint must be replaced by an inequality ($\geq$).

2. The above CP formulation assumes that the triangle-inequality holds. If those do not hold, the initial domains of the $S_i$ variables for $i\in V'$ must be adjusted to prevent cutting-off feasible solutions. E.g. the left-hand side of the domain should be changed from $t_{0,i}$ to the shortest path from $0$ to $i$. Similar for the RHS.

Answer 2

Can't speak on CP to model this but we have solved this by employing Graph Theory to model all the given constraints of yours and more on real data.

Your constraints resemble the Hamiltonian Paths Problem and in our implementation we didn't have to go as far as employing CP, MIP or any solver because we understand the instance. The way I see it, it is a problem of managing multiple moving agents in a network. Since you asked, we used a basis function and developed a Graph Theory based solver in-house for this.

The response time we are getting is exceptional, there is visible optimization and we are able to derive even the combinations to better utilize the fleet while following all the constraints one would like to apply. Hope it helped and gave you an alternative.