1 Introduction↩︎

Finite domain propagation is a powerful approach to solving complex combinatorial satisfaction and optimization problems. The core of a finite domain propagation solver are the propagators: given a constraint \(c\) and a domain \(D\) representing the set of possible values of the variables in \(c\), they remove the values of the variables which cannot take part in any solution subject to \(c\). The propagation solver interleaves the fix-point computation of the propagators with search to find a solution. In this paper, we study how we should propagate constraints of the form \(x - y \leq d\).

Difference (logic) constraints (also called temporal constraints), that is \(x - y \leq d\), where \(x\) and \(y\) are variables and \(d\) is an integer constant, are one of the simplest forms of constraints. They are well-studied, with efficient algorithms known for their satisfaction and implication, because of their connection to shortest path algorithms.

But traditionally, finite domain propagation solvers do not make use of these algorithms. Difference constraints are represented as individual propagators just like any other constraints. Unlike most constraints, propagation on difference constraints is complete. That is, if \(C\) is a set of difference constraints, then a propagation solver starting from domain \(D\) will fail if \(C \wedge D\) is unsatisfiable.

But the order of the propagation of the constraints can make a significant difference. In a sense, the propagation engine implements a version of Fords labelling and scanning algorithm (see [1], page 558) that also checks for negative cycles. The potential bad behaviour of propagation on constraints of difference is well recognized.

Half-reified constraints [2] of difference \(b \Rightarrow x - y \leq d\), that is if \(b\) then the constraint \(x - y \leq d\) holds, are very useful in expressing scheduling and other problems. Finite domain propagation engines typically implement them as weak bounds propagators that only propagate to either propagate the constraint \(x - y \leq d\) when \(b\) is known true or set \(\lnot b\) if the bounds of \(x\) and \(y\) force it to be false. Finite domain propagation engines are not complete for checking implication of half-reified difference constraints by difference constraints.

SAT Modulo Theories (SMT) solvers [3], [4] treat difference constraints quite differently than finite domain propagators. Effectively, they treat fully reified difference constraints \(b \Leftrightarrow x - y \leq d\), equivalent to \(\left(b \Rightarrow x - y \leq d \right) \wedge \left(\neg b \Rightarrow y - x \leq -d-1\right)\), where \(b\) acts as the “name” of the constraint in the SMT solver. They use specialized shortest path algorithms [5] to determine the unsatisfiability of the difference constraints, which are consequences of the current Boolean evaluation. For example, if \(b\) is true then \(x - y \leq d\) and if \(b\) is false, then \(y - x \leq -d- 1\). They also determine all the difference constraints that are entailed or disentailed by the difference constraints which are consequences of the current Boolean evaluation.

In this paper, we investigate how we can build a global propagator for difference constraints for use in a finite domain propagation engine, that treats all difference constraints together, that can determine all consequences of half-reified constraints. As we shall see, the questions that arise in a finite domain propagation context are different from those in SMT, and there are different trade-offs for the implementation. To use this propagator in a lazy clause generation solver [6], we need to be able to explain its propagators, just like in an SMT solver.

An earlier much shorter version of this paper was originally published in PPDP2008 [7]. The contributions of this paper are:

• improved algorithms for bounds propagation using global difference logic,

• the development of hybrid propagation approaches that combine the global propagator with individual difference constraint propagators,

• the first implementation of a global difference logic propagator including explanations in a general purpose finite domain solver,

• the first evaluation of different explanation options, and

• experiments illustrating that treating difference logic constraint globally can be substantially more efficient than treating them as individual constraints.

The remainder of the paper is organized as follows. 2 introduces our notation for finite domain solving, and illustrate the usual propagators for difference constraints. 3 introduces notation for graphs and shortest paths, and gives the fundamental theorem relating difference constraint satisfaction and implication to shortest paths. We then explain the state-of-the-art incremental algorithms for difference constraints. Then, 4 defines the capabilities and default implementation of a global difference constraint propagator, before discussing how to improve it by handling bounds specially, and how to explain the resulting propagations. 5 shows the integration of the global difference logic propagator in a modern lazy clause generation solver including preprocessing and simplification. Finally, 6 gives experimental results, 7 discusses related work, and 8 concludes the paper.

2 Finite Domain Solvers↩︎

We now introduce our terminology for finite domain propagation based solving. Let \({\cal V}\) be a set of (integer) variables. For theory purposes we will treat Boolean variables as 0-1 integers, but for clarity of exposition use \(\bot = 0\) and \(\top = 1\) to represent values taken by Boolean variables. In an abuse of notation \(b\) and \(\neg b\) for Boolean variable \(b\) are treated as shorthand for the expressions \(b = \top\) and \(b = \bot\), respectively. A valuation, \(\theta\), is a mapping of variables to values, denoted \(\{x_1 \mapsto d_1, \ldots, x_n \mapsto d_n\}\). We define \(\operatorname{vars}(\theta) = \{x_1, \ldots, x_n\}\). We can apply a valuation to a variable \(\theta(x_i)\) to return the value \(d_i\), and extend application of valuations \(\theta\) to arbitrary expressions involving \(\operatorname{vars}(\theta)\) in the obvious way.

A primitive constraint, \(c\), is a set of valuations over a set of variables \(\operatorname{vars}(c)\). A valuation \(\theta\) is a solution of \(c\) if \(\{ x \mapsto \theta(x) ~|~ x \in vars(c) \} \in c\). A constraint \(C\) is a conjunction of primitive constraints, which we often treat as a set. A valuation \(\theta\) is a solution of constraint \(C\) if it is a solution for each \(c \in C\). We write \(C_1 \models C_2\) if every solution of \(C_1\) is a solution of \(C_2\). We extend this notation to valuations, writing \(\theta \models C\) if \(\bigwedge_{i=1}^n x_i = d_i \models C\), where \(\theta = \{x_1 \mapsto d_1, \ldots, x_n \mapsto d_n\}\).

An atomic constraint is a unary constraint (we can restrict to the forms \(x = d, x \neq d, x \geq d, x \leq d\)), or \(\mathit{false}\). A domain \(D\) is a conjunction of atomic constraints over \(vars(D)\). \(D\) is a false domain if it has no solutions. We use notation \(D(x) = \{ \theta(x) ~|~ \theta \text{~is a solution of~} D \}\). We use range notation \(\left[l\;..\;u\right] = \{ d ~|~ l \leq d \leq u\}\). A domain \(D\) is a range domain if for each variable \(x \in vars(D)\) then \(\exists l,u: D(x) = \left[l\;..\;u\right]\). We will only be interested in range domains for the purpose of this paper. We can map a valuation \(\theta\) to a domain \(D_0 = \wedge_{x \in vars(D)} x = \theta(x)\).

A propagator \(p(c)\) for constraint \(c\) is an inference algorithm, it maps a domain \(D\) to a conjunction of literals \(p(c)(D)\), where \(D \wedge c \models p(c)(D)\). We assume each propagator is checking. That is, if \(\forall x \in vars(c). |D(x)| = 1\), then \(p(c)(D) \wedge D \models \mathit{false}\) if and only if \(\theta_D\) is not a solution of \(c\). A domain propagator \(dom(c)\), for constraint \(c\), is the strongest possible propagator for \(c\): after propagation, the domain \(D' = D \wedge dom(c)(D)\) only contains values in \(d \in D'(x), \forall x \in {\cal V}\), where there is a solution \(\theta\) of \(D \wedge c\) and \(\theta(x) = d\). A propagation solver \(prop(\mathcal{P},D)\) applied to a set of propagators \(\mathcal{P}\) and a domain \(D\) repeatedly applies the propagators \(p \in \mathcal{P}\) adding new atomic constraints from propagators to \(D\) until it reaches a domain \(D'\) where \(D' \models D' \wedge p(D')\) for all \(p \in \mathcal{P}\), and returns \(D'\). In an abuse of notation we write \(prop(C,D)\) to refer to \(prop(\{ dom(c) \mid c \in C\}, D)\).

A constraint satisfaction problem (CSP) is a constraint \(C\), often broken into a domain constraint and the remainder \(C \leftrightarrow D \wedge C'\).

In lazy clause generation (LCG) solvers [6] propagators are also required to give explanations for each new consequence \(l \in p(c)(D)\) (i.e. where \(D \not\models l\)). That is, an explanation clause \(e \equiv l_1 \wedge \dots \wedge l_n \Rightarrow l\), such that \(\forall 1 \leq i \leq n, D \models l_i\) and \(c \models e\). LCG solvers, like SAT solvers, create an implication graph, where every new consequence is attached to a reason. On failure, this is used to create a nogood by repeatedly replacing literals in the explanation of failure until only one literal that became true after the last decision remains. This nogood is guaranteed to generate new propagation information. See [6] for more details.

3 Solving Difference Constraints↩︎

Difference constraints are highly connected to shortest paths, a statement we shall formalize shortly. In this section we give our graph notation, and then explain how the difference constraints algorithms work.

3.1 Graphs, Paths and Potential Functions↩︎

A weighted directed graph \(G = (V,E)\) is made up of vertices \(V\) and a set \(E\) of weighted directed edges \((u, v, d)\) from vertex \(u \in V\) to vertex \(v \in V\) with weight \(d\). We also use the notation \(u \stackrel{d}{\to} v\) to denote the edge \((u,v,d)\). A path \(P\) from \(v_0\) to \(v_k\) in graph \(G\), denoted \(v_0 \rightsquigarrow v_k\), is a sequence of edges \(e_1, \ldots, e_k\) where \(e_i = (v_{i-1}, v_i, d_i) \in E\). A simple path \(P\) is a path where \(v_i \neq v_j, 0 \leq i < j \leq k\). A (simple) cycle \(P\) is a path \(P\) where \(v_0 = v_k\) and \(v_i \neq v_j, 0 \leq i < j \leq k \wedge (i \neq 0 \vee j \neq k)\). The path weight of a path \(P\), denoted \(w(P)\) is \(\Sigma_{i=1}^k d_i\).

Let \(G\) be a graph without negative weight cycles, that is, without a cycle \(P\) where \(w(P) < 0\). Then, we can define a shortest path from \(v_0\) to \(v_k\), which we denote by \(SP(v_0,v_k)\), as the (simple) path \(P\) from \(v_0\) to \(v_k\) such that \(w(P)\) is minimized. Let \(wSP(x,y) = w(SP(x,y))\) or \(+\infty\) if no path exists from \(x\) to \(y\). Given a graph \(G\) and vertex \(x\) define the functions \(\delta_x^\leftarrow, \delta_x^\rightarrow: V\rightarrow \mathbb{R}\) as \[\delta_x^\leftarrow (y) = wSP(x,y)\qquad \text{and}\qquad \delta_x^\rightarrow (y) = wSP(y,x)\enspace.\]

Let \(G\) be a graph without negative weight cycles. Then \(\pi\) is a valid potential function for \(G\) if \(\pi(u) + d - \pi(v)\ge 0\) for every edge \((u,v,d)\) in \(G\).

There are many algorithms (see e.g. [8]) for detecting negative weight cycles in a weighted directed graph. They either detect a cycle or determine a valid potential function for the graph. The standard approach is the Bellman-Ford algorithm, which is \({\cal O}(|V||E|)\)

Given a valid potential function \(\pi\) for a graph \(G=(V,E)\), we can define the reduced cost graph \(rc(G)\) as \((V, \{ (x,y,\pi(x) + d - \pi(y) ~|~ (x,y,d) \in E \})\). All weights in the reduced cost graph are non-negative, and we can recover the original path length \(w(P)\) for path \(P\) from \(x\) to \(y\) from paths in the reduced cost graph. This follows from \(w(P) = w + \pi(y) - \pi(x)\), where \(w\) is the weight of the corresponding path in the reduced cost graph. Since edges in the reduced cost graph are non-negative, we can use Dijkstra’s algorithm to calculate shortest paths in the reduced cost graph in time \({\cal O}(|V|\log |V| + |E|)\) instead of \({\cal O}(|V||E|)\) assuming the use of a Fibonacci heap.

3.2 Difference Constraints↩︎

Difference constraints are a well-studied class of constraints (e.g. [9], [10]). Difference constraints are of the form \(x - y \leq d\). Note that we can encode bounds constraints of the form \(x \geq l\) and \(x \leq u\) by selecting a dummy variable \(v_0\) to represent the fixed value 0, and encoding them as \(v_0 - x \leq -l\) and \(x - v_0 \leq u\) respectively. We call these encoded bounds constraints. We can map difference constraints to a weighted directed graph. For the remainder of this paper, let \(v_0\) denote the dummy variable.

The following well-known result characterizes how the constraint graph can be used for satisfiability and implication checking of difference constraints.

Ramalingam et al. [11] define efficient algorithms for satisfiability of difference constraints after incrementally adding or deleting a constraint. Cotton and Maler [5] define efficient incremental algorithms for difference constraints. The satisfiability algorithm for incremental addition is identical to that of [11], while they also give an algorithm for checking implication of difference constraints on addition. For our purposes, we are only interested in incremental addition algorithms so we will use the formulation of Cotton and Maler.

The incremental satisfaction algorithm for addition of [5], [11], shown in [alg:incsat], relies on maintaining a potential function \(\pi\) for the constraint graph \(G_C\). In a sense, it is an incremental Bellman-Ford algorithm. When a new constraint \(u - v \leq d\) is added, the edge \(u \stackrel{d}{\to} v\) is added to \(G_C\) and a new potential function \(\pi'\) is calculated or unsatisfiability is detected. The algorithm is \({\cal O}(n \log n + m)\) for \(m\) difference constraints on \(n\) variables (using a Fibonacci heap to implement argmin).

The implication algorithm of [5], shown in [alg:incimp], simply checks, for each difference constraint \(x - y \leq d\) of interest, whether the new edge has created a path from \(x\) to \(y\) of length \(\leq d\). It makes use of the potential function previously calculated to compute shortest paths on the reduced cost graph using Dijkstra, rather than using a more expensive algorithm that handles negative weight edges. The algorithm is \({\cal O}(n \log n + m + p)\) for \(m\) difference constraints on \(n\) variables and \(p\) constraints to check for implication.

4 Difference Constraint Propagation↩︎

In this section, we explain how we can create a global propagator for all difference constraints and reified difference constraints.

4.1 A Global Difference Constraint Propagator↩︎

We can use IncSat and IncImp to create a global propagator as follows. The global propagator needs to support the following operations:

• Add a new difference constraints \(x - y \leq d\)

• Add a new bound \(x \geq l\) or \(x \leq u\)

• Add a new half-reified difference constraint \(b \Rightarrow u - v \leq d'\)

• Mark a state and backtrack to a previous marked state

Adding a new difference constraint \(x - y \leq d\) causes the propagator to run IncSat, and fails if this determines unsatisfiability. If it does not, then it runs IncImp to determine if any difference constraints \(x - y \leq d' \in C'\) become entailed. For each half-reified difference logic constraint \(b \Rightarrow y - x \leq -d'-1\), if we find that \(x - y \leq d'\) is entailed, then we set \(b=0\). This is the same as an SMT difference constraint solver modulo the use of half reification.

This is not enough, since we are also interested in any new bounds that result from the constraint addition. A corollary of 1 gives the key insight.

To extract new bounds, the propagator needs to determine for each variable \(z\) the weight of a shortest path from \(v_0\) to \(z\) via the new edge \((x,y,d)\), i.e. \(\delta^\leftarrow_x(v_0) + d +\delta^\to_y(z)\). The negation of this is a lower bound on \(z\). It also determines the weight of a shortest path from \(z\) to \(v_0\) via the new edge \((x,y,d)\), \(\delta^\leftarrow_x(z) + d +\delta^\to_y(v_0)\), which is an upper bound on \(z\). The propagator updates the bounds of variables which have changed regarding these weights.

Adding a new bound \(x \leq u\) or \(x \geq l\) is simply adding a new difference constraint \(x - v_0 \leq u\) or \(v_0 - x \leq -l\). This is then treated as described in the previous paragraphs.

Adding a new half-reified constraint \(b \Rightarrow u - v \leq d'\) is a new feature not usually present in SMT solvers, since in this context all reified constraints are usually known from the beginning of solving. We need to check if the constraint is disentailed with the current set of difference constraints. To do so, we calculate \(wSP(u,v)\) and \(wSP(v,u)\) by using Dijkstra’s algorithm on the reduced cost graph. Note that half-reified bounds constraints \(b \Rightarrow u \leq d'\) or \(b \Rightarrow-v \leq d'\) can be handled by examining the domain \(D\) directly.

For backtracking, we need to restore the state of the difference constraint propagator to that at an earlier time. The only state of the solver is the set of difference constraints posted, and the set of posted half-reified difference constraints, as well as the valid potential function. The potential function \(\pi\) remains valid on backtracking since it simply ensures that \(\pi(x) + d - \pi(y) \geq 0\) for all edges \((x,y,d) \in G\), and removing edges does not invalidate this. Hence, the potential function need not be trailed, as noted by Wang et al. [12]. Backtracking must simply remove the representation of the constraints added since the marked state.

Under the assumptions that the domain \(D\) is a range domain and no Boolean variable appears twice in the set of difference constraints, bounds constraints, and (added) half-reified difference constraints \(C\), the global propagator defined in this subsection is a domain propagator for (the conjunction) \(C\).

4.2 Handling Bounds Constraints Better↩︎

While in the SMT context bounds are simply another form of difference constraint, for a propagation engine bounds updates are much more frequent than difference constraint additions. Hence, it is worth treating them separately. We will not use the dummy variable \(v_0\) to encode bounds constraints but treat them directly. The basis of this treatment is the following theorem.

The above theorem implies we can check the satisfiability of bounds constraints using propagation on the difference constraints, and we can split the implication into two checks: just using the difference constraints and just using the computed bounds.

At first, let us consider bounds updating. We define an algorithm that simultaneously considers all bounds changes since the last time the global propagator was run, as opposed to adding them one by one, which is required by the base approach.

Since the calculation of the new lower and upper bounds is symmetric, we describe only the algorithm IncLB (see [alg:BoundCheck]) for the lower bounds.

For the new lower bounds of a variable \(x\), we want to know if the negative of the weight of the shortest path to \(x\) from \(v_0\), \(- \delta^\to_{v_0}(x)\) is greater than the current lower bound of \(x\), \(\min D(x)\). Because we do not represent a dummy variable \(v_0\) in the constraint graph, we compute a valid potential function value \(\pi(v_0)\) for \(v_0\) regarding \(\pi\) in the first step. For this calculation, we consider only variables whose current lower bound is greater than their lower bound resulting from the last run of IncLB.

At the second step, we run Dijkstra on the reduced cost graph with a starting priority queue of the variables which have changed lower bound since the last run of the propagator. The initial value for variable \(x\) is the reduced cost of the “imaginary” edge between \(v_0\) and \(x\). At last, we create the bounds constraints for calculated new bounds.

Our Dijkstra’s algorithm does not explore all variables in the graph, it visits only variables \(x\) for which \(-\delta^\to_{v_0}(x) > \min D(x)\) holds, that is where a new lower bound has been found.

Note that IncLB (and IncUB) requires \({\cal O}(n \log n + m)\) time and \({\cal O}(n + m)\) space for \(m\) difference constraints (not bounds constraints) on \(n\) variables. Incremental propagation of difference constraints using a FIFO queue of propagators¹ can require \({\cal O}(nm)\) time.

None

We can thus check satisfiability of addition of bounds constraints by running IncLB on the new lower bounds and IncUB on the new upper bounds and seeing if in the result any variable has an empty domain. Afterwards, we can check new implications of difference constraints caused by the addition of bounds constraints by simply checking for each \((x - y \leq d) \in C'\) if it is implied by bounds, i.e. \(\max D(x) - \min D(y) \leq d\).

We check satisfiability on addition of a new difference constraints \(u - v \leq d\) by running IncSat on the graph \(G_C\) containing only the difference constraints (not bounds constraints), and then perform bounds propagation by determining the possibly new lower bound for \(v\) given by \(-d + \min D(u)\) and doing bounds propagation with IncLB and the possibly new upper bound for \(u\) of \(d + \max D(v)\) and doing bounds propagation with IncUB. We can then check implication by first checking if the bounds imply some difference constraint in \(C'\), and then running IncImp on the graph \(G_{C\cup \{u -v \leq d\}}\).

Adding a new half-reified difference constraint \(b \Rightarrow u - v \leq d'\) requires us to check if it is disentailed by the domain \(D\), and if not we calculate \(wSP(u,v)\) and \(wSP(v,u)\) using the reduced cost graph of \(G_C\).

4.3 Explanations↩︎

In an LCG solver, each propagation needs to be justified by an explanation that gives the reason for the propagation in clausal form. Difference logic relies on three different kinds of propagations.

For the simplicity of explaining, we extend the notion of an edge \(s \stackrel{d}{\to} t\) in the constraint graph \(G_C\) with a Boolean variable \(b\) and depict it with \(s\mathrel{\substack{d\\\longrightarrow\\(b)}} t\). In the case the edge originates from a different constraint \(s - t \leq d\) then \(b = \top\). In the other case, \(b\) is the same Boolean variable as from the corresponding half-reified constraint represented by the edge. In addition, we denote \(\mathit{pred}(t)=s\) the predecessor of \(t\) on the shortest path to \(t\) and \(\mathit{pred}^\rightarrow(t)=b\) the Boolean variable of the edge \(\mathit{pred}(t)\mathrel{\substack{d\\\longrightarrow\\(b)}} t\), both can be determined in line 16 of [alg:BoundCheck], when a new lower bound of \(t\) is detected.

4.3.1 Explaining Bound Propagations↩︎

The simplest forms of propagation are bound propagations in IncLB and IncUB. While the algorithms explore bound propagations jointly for the whole constraint graph \(G_C = (V, E)\), each individual bound can be explained locally by the arc that caused the bound change. Note this explanation is more reusable than more complex explanations. For lower bound changes in IncLB, if \(\gamma(t)\) is updated in line 16 of [alg:BoundCheck] along edge \(s\mathrel{\substack{d\\\longrightarrow\\(b)}} t\), the explanation for the lower bound update of \(t\) is that the edge is active together with the bound of the predecessor: \[\mathit{pred}(t)\geq-\delta_{v_0}^\rightarrow(\mathit{pred}(t))\land \mathit{pred}(t)^\rightarrow \Rightarrow t\geq-\delta_{v_0}^\rightarrow(t))\]

The upper bound updates are explain analogously following edges in opposite direction.

4.3.2 Explaining Cycles of Negative Length↩︎

The other two types of propagations both propagate Boolean variables used in implications to \(\bot\). The simpler case is that a cycle of negative length is found at the addition of a new edge in IncSat. This either triggers a conflict if it is the consequence of adding a new edge \(u\mathrel{\substack{d\\\longrightarrow\\(b)}} v\) with \(b\mapsto\top\) or it is applied to half-reified edges after execution of IncImp to propagate and explain setting \(b\mapsto\bot\). Similar to the previous case, the predecessor \(\mathit{pred}(t)=s\) and the Boolean implying the preceding edge \(\mathit{pred}^\rightarrow(t)=b\) is stored in line 10 of [alg:incsat]. \(\gamma(u)<0\) in line 11 entails that there is a path of negative length from \(u\) back to \(u\), so the chain of predecessors in this loop is unrolled for the explanation: \[\mathit{pred}^\rightarrow(u)\land\mathit{pred}^\rightarrow(\mathit{pred}(u))\land\ldots \Rightarrow\neg b\]

The edge under consideration from \(u\) to \(v\) closing the cycle will never be added to the predecessor relation, and no other cycles of negative length exist when applying IncSat, so the explanation chain is guaranteed to end at \(v\) which forms the root of the predecessor tree. Interestingly, Cotton and Maler [5] who define [alg:incsat] do not really describe how to explain negative cycles, though they use this in their implementation.

4.3.3 Explaining Boolean Propagations Based on Bounds↩︎

Finally, there is a third type of propagation that sets Boolean variables \(b\) of disentailed constraints \(b \Rightarrow x - y \leq d'\) where the constraint is disentailed by the bounds on \(x\) and \(y\), so \(\min D(x) - \max D(y) > d'\). The direct explanation (Simple) for this propagation is the following: \[\label{eq:expl95simple} x \geq\min D(x)\land y\leq\max D(y) \Rightarrow\neg b\tag{1}\]

However, this explanation is not necessarily as strong as it could be. The propagation occurs as soon as \(\min D(x)-\max D(y)>d'\), while the explanation does not consider how much above \(d'\) the difference actually is. Therefore, the explanation can be lifted (Lifted): \[\mathit{expl}(\neg b)=\begin{cases} x\geq\max D(y)+d'+1\land y\leq\max D(y) \Rightarrow\neg b & \text{if }\min D(x) \text{ was updated}\\ x\geq\min D(x)\land y\leq\min D(y)-d'-1 \Rightarrow\neg b & \text{otherwise} \end{cases}\] Now the most recently moved bound (or the lower bound if both were updated) is included with the least restrictive value that causes propagation of the Boolean variable.

Finally, a third option under consideration is to produce the lifted explanation lazily (Lazy), i.e., only when actually needed, and based on literals already in use by the solver (so lifting might occur on both bounds ensuring that the combination of both bounds is still a valid explanation).

4.4 Possible Variations and Optimizations↩︎

2 also provides us with other ways of building difference constraint propagators by mixing the standard approach of a propagator per difference constraint and reified difference constraint with the global approach.

Standard propagators will perform bounds propagation and implication by bounds whereas the global propagator performs a consistency check of the system of difference constraints. Hence we can combine the global propagator which only considers addition of difference constraints \(x - y \leq d\) and half-reified difference constraints \(b \Rightarrow u - v \leq d'\) with the usual propagators for these constraints. By prioritizing the global propagator before the standard propagators we can avoid the worst case behavior of 1 determine implication purely by difference constraints, and use the standard incremental queueing for bounds propagators/implication by bounds.

We can further improve upon the algorithms above by taking into account fixedness of variables, and reasoning better on which difference constraints in \(C'\) can be implied by a constraint addition.

Any fixed variable \(x\) where \(D(x) = \{d\}\) acts similarly to the dummy variable \(v_0\) since any path \(x \rightsquigarrow y\) implies a path \(v_0 \stackrel{-d}{\rightarrow} x \rightsquigarrow y\), and similarly \(y \rightsquigarrow x\) implies the path \(y \rightsquigarrow x \stackrel{d}{\to} v_0\). Once a variable is fixed and we have updated the bounds caused by the fixing, then the variable plays no further role. We can ignore fixed variables \(x\) in IncSat, IncImp, IncLB and IncUB without compromising correctness using a result analogous to 1. Hence edges to and from fixed variables need not be considered in these algorithms.

After every addition of a constraint we need to check the constraints in \(C'\) for implication. We can do better than checking every one if we keep track of what changes have been made by the latest constraints addition.

Cotton and Maler [5] show how we can restrict the shortest path calculations in IncImp to those that may actually decrease a shortest path. There is a shorter path from \(x\) to \(y\) via new edge \((u,v,d)\) if and only if there is a shorter path from \(x\) to \(v\) or \(u\) to \(y\) via this edge. If we calculate \(\delta_v^\leftarrow\) (as opposed to \(\delta_u^\leftarrow\)) and \(\delta_u^\rightarrow\) we can modify Dijkstra’s algorithm to restrict attention to those nodes that give a shorter path using the edge \((u,v,d)\) (see [5] for details). We can calculate the weight of a shortest path from \(x\) to \(y\) via \((u,v,d)\) as \(\delta_v^\leftarrow(x) - d + \delta_u^\rightarrow(y)\).

Note that \(wSP(x,y) = w + \pi(y) - \pi(x)\) where \(w\) is the weight of a shortest path from \(x\) to \(y\) in the reduced cost graph. Clearly \(w \geq 0\) and hence \(wSP(x,y) \geq \pi(y) - \pi(x)\). Hence if \(\pi(y) - \pi(x) > d\) we know that \(x - y \leq d\) is not implied by the constraints. Note also that \(\pi\) only changes by reducing the values at particular nodes (See IncSAT). Hence if \(\pi(y) - \pi(x) > d\) then at some future time \(\pi'(y) - \pi'(x) > d\) unless the potential function value at \(y\) changed.

Given both the above observations we can improve the checking of implication of constraints in \(C'\) by only checking constraints \(x - y \leq d \in C'\) if \(\delta_u^\rightarrow(y)\) made use of the edge \((u,v,d)\) in its calculation, and if \(\pi(y) - \pi(x) > d\) in the past, then \(\pi(y)\) has changed.

5 Difference Logic in Huub↩︎

Huub [13], [14] is a modern CP solver developed in Rust that is based on the IPASIR-UP [15] interface to modern SAT solvers. This section provides specifics about the interaction of a difference logic component in a modern CP solver, about pre-processing and simplification, and proposes several control options that are evaluated in the experiment section.

5.1 Difference Logic Representation↩︎

In Huub, difference logic is represented by a graph \(G=(V,E,E')\), where \(V\subseteq\mathcal{V}\) is a subset of the variables, and \(E\) is the set of active edges \(u\stackrel{d}{\to} v\) representing active difference constraints. Additionally, \(E'\) is a set of implied edges \(u \mathrel{\substack{d\\\longrightarrow\\(b)}} v\) representing the implied (half-reified) difference constraint \(b\Rightarrow u-v\leq d\) for an unknown Boolean variable \(b\). These edges will either be dropped if \(b\mapsto\bot\), or moved to the active edges \(E\) if \(b\mapsto\top\). Note that this representation is enough to cover all kinds of difference constraints, including (half-)reified equality and inequality constraints based on the following mapping:

• \(x-y\leq d\): Directly added to \(E\).

• \(b\Rightarrow x-y\leq d\): Directly added to \(E'\).

• \(b\Leftrightarrow x-y\leq d\): Transformed to \(b \Rightarrow x-y\leq d\) and \(\neg b\Rightarrow y-x\leq-d-1\).

• \(x-y=d\): One variable is a transformation of the other (called a view in the solver), no need to represent as a constraint.

• \(b\Rightarrow x-y=d\): Transformed to \(b\Rightarrow x-y\leq d\) and \(b\Rightarrow y-x\leq-d\).

• \(x-y\neq d\): Transformed to \(b\Rightarrow x-y\leq d-1\) and \(\neg b\Rightarrow y-x\leq-d-1\) for a new Boolean variable \(b\).

• \(b\Rightarrow x-y\neq d\): Transformed to \(b\Rightarrow c\lor e\), \(\neg c\lor\neg e\), \(c\Rightarrow x-y\leq d-1\) and \(e\Rightarrow y-x\leq-d-1\) for new Boolean variables \(c\) and \(e\).

• \(b\Leftrightarrow x-y=d\): Transformed to \(b\Rightarrow x-y=d\) and \(\neg b\Rightarrow x-y\neq d\), then further transformation as defined above.

Therefore, all versions of potentially (half-)reified difference, equality, and inequality constraints on two variables can be represented in difference logic. However, for equality and inequality, additional Boolean variables need to be introduced, and multiple difference constraints are generated. This might hurt performance, but it might also allow better propagation within the difference logic graph. Therefore, a control parameter for the difference logic component is introduced:

Level 0:

Difference logic is deactivated.

Level 1:

Include global and (half-)reified difference constraints.

Level 2:

Additionally include half-reified equality constraints.

Level 3:

Additionally include inequality and reified equality constraints.

Depending on the selected level, the included constraints are captured when the model is built. Internally, each node \(v\) is associated with a list of incoming active edges \(E_v^+\), outgoing active edges \(E_v^-\), incoming implied edges \({E'}_v^+\), and outgoing implied edges \({E'}_v^-\). Additionally, each Boolean variable \(b\) that is used in implied edges is associated with its implied edges \(E_b\). During search, these edges dynamically move from implied to active, or are dropped from the implied edges. However, these changes need to be reverted if the solver backtracks. For this purpose, two types of trailed lists are used:

Trailed Appending:

Active edges are stored in lists that trail the current active length \(\ell\). A new element is stored at index \(\ell\) (the array is 0-indexed), and backtracking resets the active length \(\ell\). Therefore, recently added elements are not part of the active array any more after backtracking. Note that they are not actually removed, but just overwritten on future extensions. This design allows both setting implied edges to active, and adding difference constraints from external sources.

Trailed Closing:

Implied edges are stored in a type of trailed todo-list. Elements with lower index than the trailed value \(c\) are marked as closed, while other elements are still considered open. If an element \(i\) with \(i\geq c\) is closed, its position is swapped with \(c\), and \(c\) is increased by 1. On backtracking, \(c\) is decreased again, moving the element back into the active section. Note that the order of active elements is not preserved, but closing and backtracking are done in constant time. Additionally, since closed elements are kept at the beginning of the list, trailed appending could be combined with this type for dynamic addition of implied edges. Also consider that each implied edge \(u \mathrel{\substack{d\\\longrightarrow\\(b)}} v\) is contained in three lists: outgoing from \(u\), incoming to \(v\), and implied by \(b\). Therefore, each edge is associated with its current index in all three lists, updating them on every swap. Then, from each list, closing the edge can be triggered in constant time for all three lists containing the edge.

While in theory a priority queue that allows dynamic priority updates or even a Fibonacci heap is required to achieve the given complexity results, in practice a simple priority queue where duplicate entries become stale and are skipped has shown to be more efficient and is used in the implementation.

5.2 Initial Simplification↩︎

A modern CP solver typically solves a problem in two major stages. First, a simplification stage tries to reduce the problem as much as possible without going into search, by propagating constraints at the root level, identifying variables that are just projections of other variables, and simplifying constraints when possible. Only when no further simplification is possible, search is started. Since simplification is only done once, it is very important for the performance of a solver to simplify as much as possible. Difference logic allows simplification that goes beyond what propagators for the individual constraints could do, and may be able to significantly reduce the complexity of the difference logic graph, or even prove infeasibility without starting search. Already before simplification starts for the difference logic component, Huub identifies variables that are views of other variables, i.e., \(y\) is a linear transformation of another variable \(x\) defined by \(x=k\cdot y+c\).

The first time difference logic is called, the whole graph needs to be checked for consistency, and the potential function needs to be defined. While starting with an empty graph and iteratively adding all active edges would be possible, it is more efficient to deal with the whole graph at once, applying the steps in [alg:init].

The first step is to run the Bellman Ford algorithm from an imaginary initial vertex \(v_0\) that is connected via an edge \(v_0\stackrel{0}{\to} v\) to each other vertex in \(V\). It computes the shortest path from \(v_0\) to each vertex \(v\in V\) in \({\cal O}(|V|\cdot|E|)\) while being able to deal with negative edge weights. First, it allows us to detect cycles of negative length (which would lead to a shortest path of \(-\infty\)), leading to a root node proof of infeasibility. Second, setting \(\pi(v)\) to the length of the shortest path to \(v\) is a valid potential function, which is returned if no failure is detected [16].

Next, an initial propagation of bounds is performed using IncLB and IncUB setting the sets of nodes to evaluate \(\mathcal{V}_l\) and \(\mathcal{V}_u\) to the set of all nodes in the difference logic graph \(V\). This entails removing edges from \(E'\) that are implied by bounds, or setting \(b\mapsto\bot\) for edges in \(E'\) where the inverse is implied. Initial propagation of Boolean variables evaluates all Boolean variables \(b\) that imply edges in \(E'\) which already have a fixed value (either from the propagation of bounds or externally). Implied edges are dropped if \(b\mapsto\bot\), or moved to \(E\) by calling IncSat if \(b\mapsto\top\).

In a next step, nodes are removed if they already have a fixed value, or if there are no edges connecting them to the remainder of the graph (neither in \(E\) nor \(E'\)). In the first case, implied edges might still exist, e.g., \(u\) is fixed, but \(u \mathrel{\substack{d\\\longrightarrow\\(b)}} v\in E'\) for an unknown Boolean variable \(b\). Now the implied edge reduces to an implied bound constraint for \(v\), \(b\Rightarrow v\geq \theta(u)-d\), which is passed back to the solver to be handled directly (analogously \(b\Rightarrow u\leq \theta(v)+d\) if \(v\) is fixed).

The last major step in the initialization is to compute all pairwise shortest paths in the active graph by repeated execution of Dijkstra’s algorithm, which together with the initialization of the potential function by Bellman Ford constitutes Johnson’s algorithm [16]. It runs in time complexity \(O(|V|^2\log|V|+|V||E|)\), which can be very significant for large graphs, but can lead to substantial simplifications of the difference logic graph which benefit all subsequent evaluations during search. The simplifications performed with the resulting matrix of shortest paths between all vertices \(\mathcal{D}_{uv}\) are:

Removal of redundant global edges:

\(u\stackrel{d}{\to} v\) is removed if \(d>\mathcal{D}_{uv}\), additionally, at most one edge from \(u\) to \(v\) is kept if multiple edges \(u\stackrel{\mathcal{D}_{uv}}{\to} v\) are defined.

Removal of redundant implied edges:

\(u\mathrel{\substack{d\\\longrightarrow\\(b)}} v\) is removed if \(d\geq\mathcal{D}_{uv}\).

Fixing Boolean variables to \(\bot\):

Addition of \(u\mathrel{\substack{d\\\longrightarrow\\(b)}} v\) would lead to a negative cycle if \(-d>\mathcal{D}_{uv}\), leading to fixing \(b\mapsto\bot\).

Node unification:

\(\mathcal{D}_{uu}=0\) entails that there is a cycle of length \(0\). The only valid solution to such a cycle is to set \(u-v=d\) for each edge \(u\stackrel{d}{\to}v\) in the cycle. Therefore, variable \(v\) can be replaced with a view \(v\mapsto u-d\), collapsing all variables in the cycle to a single variable.

Finally, nodes that got fixed or isolated during the previous step are removed again, and Boolean variables not connected to any implied edges any more (because all of them got removed) are also removed from the infrastructure in the graph. Afterwards, the first iterative simplification is triggered.

5.3 Iterative Simplification↩︎

Since simplification of other propagators can allow further simplification in difference logic, these iterative simplification steps are executed until no propagator can simplify any further. The iterative simplification performs the following steps:

• Merge unified variables

• Propagation

• Removal of nodes and Boolean variables

• Check for subsumption

Already before initialization, variables that are just linear transformations of other variables are replaced by views. However, further views can be identified during simplification both by the difference logic component and other propagators. Additionally, in difference logic, multiple variables that are views of the same underlying variable can be represented by the same node in the difference logic graph. Assume \(v=u+c\), even if \(u\) and \(v\) are originally not connected by any difference constraints, they can be represented by the same node in the difference logic graph. Simply replace every edge \(v\stackrel{d}{\to}w\) by \(u\stackrel{d-c}{\to} w\), analogously for edges in the opposite direction and implied edges. Only views with a multiplicative factor like \(v=2\cdot u\) are represented as separate nodes, while all additive transformations are merged. This check is performed first in each iterative simplification.

At the end of propagation in each iterative simplification, again fixed and isolated nodes are removed, and irrelevant Boolean variables are removed. Finally, if all nodes have been removed, subsumption is returned. Otherwise, new iterative simplification might be triggered. If the difference logic component is not subsumed, is finally transformed into its propagator form for the search stage of the solver. In this stage, nodes are no longer removed, as all changes now need to be trailed to be able to revert them. Only the addition of edges to \(E\) and the removal of edges from \(E'\) are possible from this moment on as described in the trailing infrastructure.

5.4 Propagation↩︎

During search, propagation of the difference logic component is triggered every time a bound of any node changes, or any of the Boolean variables that imply edges in \(E'\) is fixed. These changes are communicated by advisors [17]: The difference logic component registers for updates of all node variables and relevant Boolean variables, and gets notified of each of their changes. Therefore, no exhaustive search of bound changes over the nodes is required, but \(\mathcal{V}_l\) and \(\mathcal{V}_u\) for IncLB and IncUB as well as the set of fixed Boolean variables \(\mathcal{B}\) are directly built based on advisor calls.

The propagation of bound changes first executes IncLB and IncUB, followed by checking all implied edges connected to nodes with bound updates if they are implied by the bounds (remove them from \(E'\)) or falsified by the bounds (remove them from \(E'\) and propagate the corresponding Boolean \(b\mapsto\bot\)).

The propagation of Boolean changes entails dropping an edge from \(E'\) if the Boolean is set to \(\bot\), or moving it from \(E'\) to \(E\) if set to \(\top\), which leads to checking for conflicts and potentially updating \(\pi\) via IncSat, and optionally checking the impact on further implied edges via IncImp. Note that this check via IncImp can be turned on or off via a parameter. While skipping this step can lead to conflicts that are only discovered later in the search (e.g., IncImp would propagate \(b\mapsto\bot\) for \(x\mathrel{\substack{d\\\longrightarrow\\(b)}} y\), instead \(b\) is later set to true and IncSat runs into a conflict), running IncImp is the most costly part of propagation, and the benefit of skipping it can outweigh the cost even though IncImp implements the improvement by Cotton and Maler [5] calculating \(\delta_v^\leftarrow\) and \(\delta_u^\rightarrow\) as described in 4.4.

Additionally, propagation of bound changes and fixed Boolean variables is independent of each other even though they both operate on the same graph and one kind of propagation might trigger the other as a consequence. This entails that the propagations can be triggered separately, and therefore with different priorities PrioBounds and PrioBools.

6 Experiments↩︎

We consider three implementations in Huub in the experiments: sp: only separate propagators for each constraint (the standard approach); gp: a global propagator using IncSat and IncImp on difference and reified difference constraints, and IncLB and IncUB to handle bounds. gp+s: gp plus the additional simplification stage using Johnson’s algorithm followed by the removal of redundant edges and nodes. We additionally evaluate several different parameter settings for gp+s. We do not consider the base implementation of 4.1 since it is an order of magnitude slower than gp when there are bounds constraints.

All the experiments were carried out on a cluster equipped with Intel Xeon Silver 4314 (2.4 GHz, 24 MB Cache, no hyperthreading) running Ubuntu 22.04 set up for maximum repeatability [18]. Each experiment uses a single thread and at most \(12.8\) GB of RAM. Huub [13] is evaluated in a development build² and compared to Chuffed 0.13.2 [19] and Google OR-Tools CP-SAT 9.15.6755 [20] as the two of the best performing solvers in recent MiniZinc challenges [21]. We specifically compare to other solvers that are applicable to general problems specified in MiniZinc, while other solvers implementing a global difference propagator (see 7) are not general purpose finite domain solvers, but more specialized to specific domains. Full tables of the results on all instances are available in the online supplementary material [22].

6.1 Worst case behavior for separate propagators↩︎

8 shows a problem where the use of single propagators could be terrible compared to the global propagator gp just for bounds propagation. In the first experiment, we validate this experimentally, comparing the separate propagators sp and the global propagator gp on this problem. sp uses the standard implementation of linear inequality propagators in Huub. For different number of variables \(n\), we applied consecutive lower bound updates of \(y_0\) from \((k-1)n\) to \(kn\), \(k \in 1..10\), and measured the time needed to reach these 10 fixpoints. The results are show in 4. The gp times are all below \(70\) ms, while the sp times show the expected cubic growth. One can see from the results that very bad propagation behavior can occur in practical solvers when using single propagators, while a global propagator only needs to run once to reach a fixpoint.

5 shows the comparison in logarithmic scale, including not just the propagation time, but also the initialization and simplification time. For gp, initialization and simplification actually takes the majority of the time, but still stays below \(0.25\) s for all experiments. For sp, initialization time is similar to gp in absolute values, but barely notable compared to the huge propagation time.

Additionally, we compare gp+s, which includes the full simplification using Johnson’s algorithm. While this algorithm also has cubic runtime like sp, runtime is still orders of magnitude lower compared to sp. First, the cubic cost in this case occurs only once during simplification, while for sp it repeatedly occurs in propagation. Second, for gp+s the cubic runtime occurs within one algorithm in the propagator, while sp shows a cubic invokation of different propagators which is orders of magnitude more costly. Therefore, even for worst-case examples, simplification shows acceptable runtime cost. Note that a similar behavior to gp+s occurs in a modified version of the worst-case example where all non-zero arcs are negative due to the initial computation of the potential function using Bellman-Ford.

6.2 Benchmark Problems↩︎

In the following sections, we evaluate three sets of benchmark instances to evaluate the effect of several configuration choices for difference logic, and to compare the performance with and without difference logic.

To test difference logic on a wide range of different problems, we evaluate the repository of MiniZinc challenge benchmarks [23] that contains all previous problems used in the yearly MiniZinc challenge since 2008. The problems are filtered to remove those that do not contain any difference constraints as captured by the propagator in this paper for any instance. Problems occurring in different years are only kept in their most recent appearance. The resulting benchmark set contains 626 instances and is run with the same timeout as available in the MiniZinc challenge, which is 20 minutes.

Further we use two specific problems that contain difference constraints as a major component of their model for evaluation. The second benchmark is a set of complex resource-constrained project scheduling benchmarks using producer / consumer constraints ProdCons [24] and precedences. The model is taken directly from this reference using the dAfter decomposition. The instances range from 10 to 1000 activities using 5 resources. The small instance set for analysing different configurations using 10 instances per size contains 120 instances in total, while the full data set contains 1080 instances.

The third benchmark in the evaluation is RCPSPmax [25], which is a variant of resource-constraint project scheduling with minimal and maximal time lags. Already the decision variant of the problem is NP-hard. Again, there is a small data set using 10 instances per size for 120 instances in total, and a full data set with 1130 instances from the CD, UBO, and SM datasets in literature ranging from 10 to 1000 activities using 5 resources (per dataset and size at most 100 instances).

6.3 Evaluation of Different Configurations↩︎

In this section we evaluate a range of difference logic configurations on the whole MiniZinc challenge benchmark and the small instance sets for ProdCons and RCPSPmax. The following choices are evaluated in all possible 162 combinations:

• Level of captured constraints: 1 (global and (half-)reified difference constraints), 2 (+ half-reified equalities), 3 (+ inequalities and reified equalities)

• Priorities: Lowest (0), medium (2), highest (4) individually for bound propagation and Boolean propagation

• Implied check: Use of IncImp for early detection of conflicts (on or off)

• Bound propagation level: Lazy + lifted (Lazy), Eager + lifted (Lifted), Eager not lifted (Simple)

6.3.1 MiniZinc Challenge Benchmarks↩︎

Since there are 162 configurations to evaluate, we use the tool mzn-bench³ to calculate a score of each configuration in comparison to the baseline configuration (level 0, difference logic is off, i.e. the standard solver sp) using the scoring system of the MiniZinc challenge [23]. For each instance, this score awards a point to the solver providing a better solution status (optimal / unsatisfiable if the other solver did not reach a proof, better objective if no solver reached a proof), or splits the points proportionally to the runtime if both solvers provide a proof of optimality or unsatisfiability).

For each of the following boxplots, we collect all results sharing a particular configuration setting (with varying other configurations) and show minimum, maximum, quartiles (Q1, Q3), and median in each graph. A dotted line at 310 points indicates the amount of points needed to match the performance of the baseline configuration.

6 shows the impact of the level of constraints captured by difference logic. Clearly, level 3 is outperformed by the other levels, with even Q3 below the baseline. Level 1 outperforms the other levels, indicating that the extra effort to include equality and inequality constraints by introducing extra variables generates an overhead that does not pay off by reaching better propagation in the difference logic graph.

7 shows the results for bound propagation priorities, while 8 shows Boolean propagation priorities. The results are very interesting as they show that for bound propagation the lowest setting is beneficial, outperforming the baseline in most cases, while medium and highest mostly reduce the number of points. In contrast, for Boolean propagation the lowest setting is clearly inferior, and the highest setting, even considering some negative outliers, can provide the best results. These results show that separating the priorities of the different propagation stages is valuable for the performance of difference logic.

9 shows the difference between activation and deactivation of IncImp. It is clear that activation mostly reduces the scores, while deactivation shows much better results. This indicates that the overhead of the IncImp check is in general not worth the effort even though earlier propagation of Boolean variables can be achieved.

Finally, 10 shows the comparison of the different types of explanations for Boolean variables propagated by bound changes. While Lifted and Lazy show similar performance, it is clear that Simple, not relying on any lifting, is inferior in most configurations. However, it has a very wide span including both the highest and lowest overall scores.

Next, to validate the individual findings and check for any strong interactions between different settings, we look at the top configurations overall. The top 5 configurations, all above 320 points, all agree on IncImp off and lowest priority for bound propagation. Four of them use highest Boolean propagation priority, one medium. Two configurations use level 1, three use level 2. The explanation types are mixed with 2 Lazy, 2 Lifted, and 1 Simple. These match well with the individual evaluations above.

6.3.2 Scheduling Benchmarks↩︎

The configuration evaluation was conducted on a very diverse set of benchmarks, ensuring robustness of the results. However, since the importance of difference logic can vary for different problem domains, we also evaluate all configurations on the small data sets for ProdCons and RCPSPmax to see the effects in more detail on problem domains where difference logic plays a larger role.

Note that for ProdCons all configurations outperform the baseline, while for RCPSPmax 157 out of 162 configurations outperform the baseline, showing a clear advantage of using difference logic on these problems. This section summarizes findings on the benchmarks, for the corresponding graphs seen 9 and 10.

For both problems there are no additional constraints to capture at levels 2 or 3, therefore the results only show noise.

More differences occur for priorities, especially for ProdCons. For bound propagation, the lowest setting actually has the highest median of points, but also the lowest Q1, Q3, and maximum, showing a less clear overall picture. For Boolean propagation, the picture from the general benchmark is actually reversed, with the lowest setting outperforming the medium and highest setting. Priority settings seem to have less impact for RCPSPmax, with a light trend of lower scores for lowest priorities in both categories. The best overall results are achieved with highest bound and medium Boolean priority. These diverse results show that the best priority setting is more problem-dependent.

IncImp off is again clearly ahead compared to the activated setting. While the result is already strong for ProdCons, for RCPSPmax it is actually so strong that all results using IncImp score strictly below all results deactivating it, making this finding the most strongly supported finding among the configuration experiments.

Finally, the choice for explanations is again more problem-specific. ProdCons shows highest points for Simple, closely followed by Lifted, both clearly ahead of Lazy, while RCPSPmax shows Lazy ahead of Simple and Lifted last.

Overall, we choose the configuration with level 1, lowest bound priority, highest Boolean priority, no use of IncImp, and Lifted explanations as our final configuration, since this configuration excels in the diverse MiniZinc challenge benchmark (\(320.95\) points, maximum achieved \(321.95\)), combines all preferred configuration options from the individual investigations on this benchmark, but also performs well on the specific problems (\(67.17\) points with maximum achieved \(67.87\) on ProdCons, \(62.12\) points with maximum achieved \(65.44\) on RCPSPmax). For RCPSPmax the top scores require higher priorities for bound propagation and Lazy explanations, but these options did not show consistent enough performance to be selected overall.

6.4 Comparison of Top Configurations on Full Instance Set↩︎

This section shows the comparison of Huub with and without difference logic to other state-of-the-art solvers on the MiniZinc benchmark and the full benchmarks on ProdCons and RCPSPmax. Specifically, Huub (sp) denotes running without difference logic, Huub (gp+s) denotes running the best configuration of difference logic according to the previous section including simplification, and Huub (gp) denotes the same configuration, but with initial simplification deactivated. For ProdCons and RCPSPmax, all versions of Huub are run with the flag --reason-eager 1 to avoid triggering a rare bug in the interface to CaDiCaL on unsatisfiability proofs.

1 shows the comparison on the MiniZinc challenge benchmark. The different columns show the number of optimally solved instances, their average runtime, the number of instances proven unsatisfiable, their average runtime, the number of instances where a satisfying solution was found, but timeout occurred (#to), the number of satisfaction instances solved, their average runtime, the number of unknown instances (no solution found within runtime, Unkn.), the number of instances aborted due to errors, and the points according to the MiniZinc challenge in direct comparison to Huub (sp).

The evaluation shows that CP-SAT, the winner of the MiniZinc challenge, is still the strongest solver on this general benchmark. As it is a portfolio solver combining different strategies, it is well suited for this benchmark as different strategies are expected to excel at different groups of instances. While it has the highest number of optimal solutions and solved satisfaction instances, as well as the highest number of points, it is outperformed by all three versions of Huub on the number of unsatisfiable instances, and Huub with difference logic has the lowest number of unknown instances on the benchmark. Further, difference logic reduces the gap between Huub and CP-SAT in points by around one third, showing the valuable contribution of this component on the general benchmark.

Without initial simplification, Huub with difference logic is still better than the original in points, but only just, showing the importance of these initial steps in the overall performance. Overall, difference logic does not increase the number of optimal or proven unsatisfiable instances on this benchmark, but the benefit seems to come from finding better solutions faster on more complex instances. Chuffed is inferior to the other competitors, showing fewer optimal solutions, more unknown instances, and more errors.

2 shows the comparison on the full ProdCons benchmark, where Huub already outperforms the competitors without difference logic. Using difference logic, the performance is greatly boosted. More instances are solved to optimality in shorter time, more instances are proven unsatisfiable in shorter time, and only 2 unknown instances remain. These results show the great benefits of difference logic on such scheduling benchmarks. The comparison to the version without simplification also shows the importance of this component. Without simplification, there is still a large improvement in the number of optimal solutions, however, there is also an increase in the number of unknown solutions. The majority of the improvements of using difference logic require working on the simplified graph.

Finally, we compare the different options on RCPSPmax in 3. Again, Huub can outperform the competitors already without difference logic, but further improve with difference logic. Again, difference logic provides more optimal and proven unsatisfiable solutions in less runtime and fewer unknown instances. The benefit of simplification is mostly seen in the number of unknown instances. Remarkably, difference logic with simplification is so efficient in detecting unsatisfiability on this benchmark that the average time is below \(0.01\) second, so most of these instances are identified during simplification.

Overall, the results show that difference logic makes Huub an even stronger competitor on a very general set of instances, and helps the solver to excel on scheduling benchmarks where it shows excellent performance in comparison to other state-of-the-art solvers. The simplification stage proves to be an important part of the practical success of difference logic.

7 Related Work↩︎

Early work on difference logic constraints (also called temporal constraints) was done by Dechter et al. [10], which introduced the simple temporal constraints network for difference logic constraints and the formal framework for reasoning about them using shortest path algorithms. This network is the same directed constraint graph \(G_C\) as defined in 1 on page , in which they also encode bounds constraints on variables with a dummy node. This work formed a basis for later work for finite domain solvers, e.g., IBM ILOG CP Optimizer.

As mentioned in the introduction, SMT solvers [3], [4] like, e.g., Z3 and cvc5, have been treating difference logic constraints on a global level and reasoning about them using specialized shortest path algorithms [5] whether the current system is satisfiable or the entailment or disentailment of reified constraints. Our work adapts these SMT methods for finite domain solvers with and without clause learning.

In the context of CP, the previous published version of this work [7] in 2008 introduced the global difference logic propagator to the CP community and implemented it in the G12 finite domain solver engine G12?. Before that difference logic constraints were modelled as separated linear constraints.

Schutt and Stuckey [24] implemented the global difference logic propagator [7] in the Lazy Clause Generation [26] solver Chuffed [19], [27] for solving problems with producer/consumer constraints, but do not provide details how any propagation is explained, except they refer to the difference logic theory propagator [5] and mentioned that bounds literals \(\left\llbracket \pm x\leq d \right\rrbracket\) are treated separately as described in [7]. They only considered reified constraints \(b \leftrightarrow x - y\leq d\), which needed to be statically added to the model before solving, whereas our solver allows dynamic additions of half-reified and reified constraints during the search.

More recently, Hebrard [28] introduces the new CP/SAT hybrid solver Tempo currently specialized for disjunctive scheduling, in which they implement a global difference logic propagator. Tempo lazily creates bound literals in a branch only when necessary and they only exist in the subtree rooted at this branch, whereas in our solver they currently exist from the start, but the Huub architecture allows for dynamic trailed additions, which we can explore in the future. They also separate lower and upper bound propagation and difference logic propagation as we do. However, their bound propagators uses a version of Dijkstra with a simple FIFO queue, which directly works on negative edges and has the worst case runtime complexity of \({\cal O}(nm)\) opposed to ours \({\cal O}(n\log n + m)\).

Their difference logic propagator supports dynamic additions of half-reified difference logic constraints, in which two constraints can also be expressed in a XOR-relationship by using the same Boolean variable, but one in the negative form. They use this relationship to model the non-overlap of two tasks sharing the same disjunctive resource. For example, the relationship between two tasks A and B models that either A must finish before the start of B or B must finish before the start of A. Our propagator can support the same relationship, but we also consider equality and inequality constraints of two variables. Implication of such half-reified constraints are only checked regarding the variable’s bounds and their explanation is the same as our explanation Lifted. They note that checking implication via the difference logic graph is too costly, but do not provide any empirical results as we do.

They empirically evaluated Tempo on various disjunctive scheduling problems with the objective of minimizing the makespan, in which they use the global difference logic propagator and the edge-finding rule from the global disjunctive propagator (also called unary or no-overlap propagator) to identify new difference logic constraints (“edges”) between tasks, which are then shared with the difference logic propagator. In our case, we do not have any propagator that infers new difference logic constraints and shares them with the global difference logic propagator at the moment. However, our propagator is implemented in a more general solver and evaluated on a wider set of problems involving different objective functions.

Aries [29] is a dedicated clause learning solver for solving only disjunctive scheduling problems minimizing the makespan that uses a global difference propagator. The non-overlapping constraints between two tasks sharing the same disjunctive resource are modelled by two reified difference constraints using the same binary variable. The makespan is modelled with difference logic constraints with any other task, which is the same in our case if the model declares inequalities between makespan and the end of each task. It uses an adaption of the Bellman-Ford algorithm described in [30] to propagate the bounds, which has a worst-case runtime complexity of \({\cal O}(mn)\) compared to \({\cal O}(n\log n + m)\). Their bounds’ propagation is equivalent to our Simple explanation.

IBM ILOG CP Optimizer [31] is a specialized solver for planning and scheduling that uses a global difference propagator for propagating the time bounds of the start and end of (optional) interval variables representing the start and end time of an activity. Since the solver is a closed-source software only limited information about the implementation and internal workings are known. In [31], they mentioned that their difference logic graph is based on the graph in [10] extended by the presence statuses of interval variables. Bound constraints are explicitly handled in [10] compared to our separation. The improved Bellman-Ford algorithm [32] is used for the initial propagation of the graph, whereas an extended version of the algorithm described in [30] is used for subsequent propagation. Both have a worst-case runtime complexity \({\cal O}(nm)\), whereas our algorithm for the subsequent propagation has only \({\cal O}(n\log n + m)\). Compared to our implementation, they also consider optional interval variables leading to different logic constraints, in that the variables can be optional. They perform time-bound adjustment on those optional variables if their presence status is implied by the other, even in the case if there presence status is unknown. For that reason, they maintain a global “implication” graph about the presence statuses of all optional interval variables. Furthermore, they identify negative cycles of such difference constraints, for which the presence status of the involved optional variables are implied for each other, and set the presence statuses to “absent” in such a case. Our propagator would be able to handle difference logic constraints with optional variables \(x-y\leq d\), but can only propagate them when they are present. In our case, such a constraint would be modelled by a half-reified constraint \(c \Rightarrow x-y \leq d\) and the constraint \(c = a\wedge b\) where the Boolean variables \(a\) and \(b\) represent the presence status of the optional variables \(x\) and \(y\).

Google OR-Tools CP-SAT solver [20] is an LCG solver that can also make use of a simplex algorithm for solving LP relaxations of linear constraints in the model. While difference logic constraints are internally represented as two-variables per inquality (TVPI) constraints (i.e. \(ax + by \leq d\) where \(a\) and \(b\) are constant numbers) and no global difference logic propagator exists, it applies a transitive closure algorithm on unit TVPIs (i.e. TVPI’s where \(a\) and \(b\) is \(1\) or \(-1\)) to tighten the bound on the constraint (i.e. \(d\)), but only at the root level. For that purpose, the algorithm builds a directly acyclic graph by adding two edges \(x\to y\) and \(y\to x\) for each constraint and then computes the transitive closure using a specific variable order. However, it comes with two caveats: First, if there is a cycle amongst the set of constraints then the algorithm is not executed. Second, the algorithm has a work limit for computing the transitive closure, once the limit is reached the algorithm stops. We note that the cycle detection does not consider the length of it and, therefore, does not identify infeasible system. Beside this algorithm, difference logic constraints and their (half-)reified versions are treated as separated propagators, but propagated in an efficient order to minimize the workload using their specific variable ordering. Regarding the internal LP solver, CP-SAT does not automatically add difference logic constraints to the LP relaxation, but decides it depending on the linear constraints in the model. We also ran CP-SAT with settings that always add the difference logic constraints to the LP relaxation, but the results did not improve.

To sum up, the works [24], [28], [29], [31] implement a version of the global difference logic propagator in a finite domain constraint solver with or without clause learning, which reasons about the consistency of the system, implication of (half-)reified constraints, and the variables’ bounds. All of them are specific implementation of a group of problems and only tested on them: disjunctive scheduling [28], [29], scheduling problems involving producer/consumer constraints [24], and planning and scheduling problems [31], whereas our propagator is implemented in a general purpose solver and tested on a wide range of problems that include difference logic constraints. In the case of [24], [31], implementation details are not known.

8 Conclusion↩︎

Difference constraints appear widely in constraint programming models, and the default representation as individual propagators has some known bad behavior. In this paper we explore how a global treatment of difference constraints can improve propagation. We explore a number of different design choices. Overall we find that applying the global propagator at lowest priority for bounds, highest priority for booleans, not running the IncImp algorithm, and applying initial simplification leads to the best results overall.

We note that while SMT [3], [4] solvers treat difference constraints globally, the naive importing of their methodology is impractical for finite domain propagation solving because of the relative importance of bounds propagation for FD. For problems with bounds constraints the approach outlined in Section 4.1 is an order of magnitude slower than gp.

The global difference solver opens the possibility of improving the behaviour of CP solvers by instrumenting globals to take advantage of known difference constraints and/or learn new difference logic consequences of propagation, as in the disjunctive propagator for Tempo. This remains as important future work.

9 Configuration Experiments on ProdCons↩︎

This section shows the graphs for the configuration experiments on the small ProdCons instance set. All configurations outperform the baseline which requires more than 60 points.

10 Configuration Experiments on RCPSPmax↩︎

This section shows the graphs for the configuration experiments on the small RCPSPmax instance set. Almost all configurations outperform the baseline which requires more than 55 points.

11 Submission Information↩︎

Author Contributions↩︎

This paper is an extension of an original conference paper [7] written by Thibaut Feydy, Andreas Schutt and Peter J. Stuckey; Thibaut Feydy has meanwhile left academia. The extended paper details the full new implementation in a new learning solver Huub [14]. Lucas Kletzander developed the difference logic implementation in Huub with help from Jip J. Dekker, the lead developer of Huub. Lucas performed experiments and wrote Section 5. The remaining authors were all involved in discussions on algorithmic design and edited the paper substantially.

Compliance with Ethical Standards↩︎

As a commentary paper there are no conflicts of interest, or other complex ethical standards to apply. Generative AI was only used to generate scripts for data collection and generation of figures and tables. All scripts and outputs were checked by the authors for correctness.

Competing Interests↩︎

Not applicable, although Peter J. Stuckey sits on the Editorial board for Constraints.

Data Availability Declaration↩︎

The solvers used in the evaluation are available as open source. The detailed result tables are available as online supplementary material.

References↩︎

---

1. A FIFO queue is the usual case, and for a LIFO queue it is easy to create even worse behavior.↩︎

2. https://github.com/kletzi/huub/tree/feat/difference_logic_benchmark↩︎

3. https://github.com/MiniZinc/mzn-bench↩︎