VRP through constraint programming

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.

• I don't believe your model is valid for multiple vehicles. I can see at least a couple of problems with the MTZ constraints. The easiest to see is this: if the optimal solution would send vehicle 1 through nodes $i$ and $j$ in that order (not necessarily consecutive), and vehicle 2 through nodes $j$ and $i$ in that order, (4) would render that solution infeasible ($u_i < u_j < u_i$).
– prubin
Aug 17, 2021 at 18:13
• To answer your question about continuous variables, I think most CP solvers can accommodate a mix of integer and continuous variables. If nothing else, you can approximate the arc transit times with integers (for instance, by using rational approximations and then multiplying everything by the least common denominator).
– prubin
Aug 17, 2021 at 18:16
• You're trying to write a CP model as you would write a MILP model. There's little chance that a CP solver would work better than a MILP solver in this case. VRP problems are rather modeled in CP with next/pred/start/load variable AFAIK. Unfortunately, I don't know any reference in English about this Aug 17, 2021 at 20:03
• Yes, you understood my second comment. It might not be necessary, depending on the CP solver. I also support the comment by @fontanf that CP models are written entirely differently from MIP models. In particular, the constraints (and variable types) available for a CP model depend on the solver you use (whereas all MIP solvers understand linear constraints and continuous/integer/binary variables).
– prubin
Aug 17, 2021 at 21:10
• To expand on my previous remark, if I were using CP for this I might start with a variable $w_i$ for each non-depot node $i$, with domain $M$, telling me which vehicle will visit node $i$. That is not something you would do with a MIP model.
– prubin
Aug 17, 2021 at 21:13

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.

• Very thank you, this is exactly what I am looking for. However, concerning the constraints $S_{next_i} \geqslant S_i + d_{i,next_i}$, can you indicate me some CP solver that allows represent such inequality? As I said, this is my first time using CP. Aug 18, 2021 at 0:13
• @MatheusDiógenesAndrade As far as I know, any CP solver can do this (IBM CP Solver/ Choco/Gecode/...). Indexing a variable by another variable is called an element constraint in CP. See here for an example. Aug 18, 2021 at 0:35
• @MatheusDiógenesAndrade I'd recommend using IBM's CP modeling for python (DOCplex.CP) as their Python CP interface is very simple, well designed and intuitive. Here are example models. In particular, look at the facility model where the element constraint is used. Aug 18, 2021 at 0:35
• very thank you. Aug 18, 2021 at 0:37
• @JorisKinable you mentioned an interval variable based approach for modeling the problem. Since I deal with the VRPTW out of pure intuition this is the approach that I took, it seems to struggle a lot on even finding feasible solutions. Have you done any comparison with the model that you proposed? I suspect that a model without interval variables would require way less variables with more compact domains. Mar 18 at 19:41

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.

• Welcome to OR.SE. I have noticed that all of your answers have included promotional links, which in general is only allowed in restricted situations. Please see here for further details on disclosure of affiliation; otherwise your posts will come across as spam. Aug 24, 2021 at 11:50
• I am new to this platform and didn't know about this. Thanks for the insight. Aug 24, 2021 at 15:27
• My intention though was to educate the right stuff and provide actual resources for the burning issues which now won't happen because of such objections. Aug 24, 2021 at 15:36