conditions:

1. There exists $k \in N$ such that $(U, U \setminus \{k\}) \in \text{eff}_{\tau}$, ($k \in U) \in \text{pre}_{\tau'}$, and each $\tau'' \in T \setminus \{x_i, x_j\}$ does not change $U$ or $(U, U \setminus \{l\}) \in \text{eff}_{\tau''}$ for some $l \in N$.

2. There exists $k \in N$ such that $(U, U \cup \{k\}) \in \text{eff}_{\tau}$, ($k \notin U) \in \text{pre}_{\tau'}$, and each $\tau'' \in T \setminus \{x_i, x_j\}$ does not change $U$ or $(U, U \cup \{l\}) \in \text{eff}_{\tau''}$ for some $l \in N$.

**Proof.** For the first case, $x_i$ removes $k$ from $U$, and no other transition adds $k$ to $U$. Since $x_j$ requires $k$ to be in $U$, once we apply $x_i$, we cannot apply $x_j$ later. The other case is proved similarly. ◀

Before starting beam search in a neighborhood, we remove a transition $\tau$ from the model if there exists a transition $\tau'$ in the suffix such that $\tau$ and $\tau'$ satisfy one of the conditions in Proposition 2. Detecting such a pair of transitions is done once at the beginning by checking the expression trees representing the preconditions and effects of transitions.

## 3.3 Bandit-Based Depth Selection

Selecting the depth of a neighborhood, $d$, in line 3 of Algorithm 3 is non-trivial. If $d$ is too small, it is unlikely that an improving solution exists. However, if $d$ is too large, each neighborhood search takes a long time. We want to select $d$ such that the total cost improvement is maximized within the time limit.

We formulate the depth selection as the budgeted multi-armed bandit problem with continuous random costs [35]. We have the set of depths $D \subseteq \{2, ..., n\}$. If we select a depth $d \in D$ and perform search in line 7 of Algorithm 3, we obtain a new solution $x$ by spending search time $t$. As $t$ is assumed to take a value in $[0, 1]$ in the budgeted multi-armed bandit problem, we divide the actual time by the time limit $T$. If the cost $\text{cost}_x(S^0)$ is smaller than the current best solution cost $\bar{f}$, the reward is $r = (\bar{f} - \text{cost}_x(S^0))/\bar{f}$. If no better solution is found, the reward is $r = 0$. We call this process a round, and we repeat rounds until reaching the time limit $T$. We do not know the reward $r$ and time $t$ before finishing a round, so we use random variables $r_{dk}$ and $t_{dk}$ representing the reward and time if depth $d$ is used at round $k$. Let $a$ be a strategy that selects a depth $a_k$ in the round $k$. The number of rounds performed by $a$ by the time limit, $K_{aT}$, is also a random variable. The objective is to find a strategy $a$ that maximizes the total expected reward $\mathbb{E}[\sum_{k=1}^{K_{aT}} r_{ak}]$.

We use Budgeted-UCB [35]. At each round, if some depths in $D$ have not been selected before, Budgeted-UCB selects one of them. Otherwise, let $m_{dk}$ be the number of rounds where the depth $d$ is selected up to round $k-1$, and let $\bar{r}_{dk}$ and $\bar{t}_{dk}$ be the average reward and search time for $d$ up to round $k-1$. Budgeted-UCB selects the depth $d$ that maximizes

$$ \frac{\bar{r}_{dk}}{\bar{t}_{dk}} + \frac{\epsilon_{dk}}{\bar{t}_{dk}} + \frac{\epsilon_{dk}}{\bar{t}_{dk}} \min\{\bar{r}_{dk} + \epsilon_{dk}, 1\} \quad (5) $$

where $\epsilon_{dk} = \sqrt{\frac{2 \log (k-1)}{m_{dk}}}$ and $\lambda$ is a positive lower bound of the search time of each round.

In practice, we initialize $D = \{2, 4, 8, ..., 2^l, n\}$, where $n$ is the length of the initial feasible solution and $l$ is the maximum integer such that $2^l < n$. If we get a solution whose length $n'$ is different from $n$ at round $k$, we replace $n$ with $n'$ in $D$ using $m_{n',k+1} = m_{n,k+1}$, $\bar{r}_{n',k+1} = \bar{r}_{n,k+1}$, and $\bar{t}_{n',k+1} = \bar{t}_{n,k+1}$ and ignore depths greater than $n'$ in $D$. If multiple depths have not been selected before or have the same value, we select the minimum depth among them. For $\lambda$, we use the time of the first round divided by 10 while there is no guarantee that it is a lower bound. Thus, the theoretical analysis of Budgeted-UCB studied in the original paper [35] does not necessarily apply to our setting. In addition, while Xia et al. [35] assumed that the pairs $\{(r_{dk}, t_{dk})\}_{k=1}^{\infty}$ are i.i.d., we do not have such a guarantee.
## 3.4 Start Selection

Once LNBS determines the depth $d$ to use, it selects a starting point $i$, which induces the prefix and the suffix, in line 4 of Algorithm 3. We select $i$ considering the cost change in a neighborhood. Concretely, the cost change by partial path $\langle x_i, ..., x_{i+d-1} \rangle$ is defined as

$$ \delta_{di} = \text{cost}_{\langle x_1, ..., x_{i+d-1} \rangle}(S^0) - \text{cost}_{\langle x_1, ..., x_{i-1} \rangle}(S^0). \qquad (6) $$

Since the path cost is non-decreasing, $\delta_{di} \ge 0$. We ignore $i$ with $\delta_{di} = 0^1$ and select one uniformly at random from the remaining options.

Our second approach is to select $i$ based on the probability biased by $\delta_{di}$. As we explain in the next subsection, for each $d$ and $i$, the beam width $b$ is maintained. Since smaller $b_{di}$ leads to a shorter search time, we discount the probability of selecting $i$ by $b_{di}$. Concretely, given the depth $d$, we can select the starting point $i$ with the probability

$$ p_{di} = \frac{\delta_{di}/b_{di}}{\sum_{j=1}^{n-d+1} \delta_{dj}/b_{dj}}. \qquad (7) $$

## 3.5 Beam Width Selection

Given the depth $d$ and the starting point $i$, LNBS selects a beam width $b$ in line 5 of Algorithm 3. Here, we use a similar strategy to CABS: for each $d$ and $i$, we initialize the beam width $b_{di}$ to be 1 and update it to $2b_{di}$ after each round with $d$ and $i$. If we find an improved solution in line 7, we reset $b_{d'i'} = 1$ only for $d'$ and $i'$ such that $i' > i$ or $i' + d' < i + d$; if $i' \le i$ and $i' + d' \ge i + d$, the prefix and the suffix for the neighborhood induced by $i'$ and $d'$ do not change, and we know that a better partial path was not found with beam widths smaller than $b_{d'i'}$.<fnref n="2" /> If the neighborhood is exhausted, i.e., complete = $\top$ in line 7, we ignore the combination of $d$ and $i$ in lines 3 and 4 until a new solution is found and $b_{di}$ is reset to 1. Since the number of neighborhoods is finite, LNBS eventually exhausts all the neighborhoods and finds the optimal solution, which guarantees completeness.

# 4 Experimental Evaluation

We compare LNBS with CABS, DD-LNS, CP, and mixed-integer programming (MIP).

## 4.1 Experimental Settings

As benchmarks, we use the same problems and instances used by previous work [26]: the traveling salesperson problem with time windows (TSPTW), the capacitated vehicle routing problem (CVRP), the multi-commodity pickup and delivery traveling salesperson problem (m-PDTSP), the single machine total weighted tardiness problem ($1 \parallel \sum w_i T_i$), the talent scheduling problem (Talent), the simple assembly line balancing problem to minimize the number of stations (SALBP-1), the bin packing problem (BPP), and the graph-clear problem (GCP). We use the same DP, CP, and MIP models as the previous work [26].<fnref n="3" />

<footnote n="1">¹ Theoretically, a better solution may be found with such $i$ if the suffix is not empty because a partial path may change the state from which the suffix is applied, which may change the cost of the suffix.</footnote>

<footnote n="2">² Theoretically, a better solution may be found with beam width smaller than $b_{d'i'}$ if the updated primal bound changes the search behavior.</footnote>

<footnote n="3">³ https://github.com/Kurorororo/didp-models</footnote>
LNBS, CABS, and DD-LNS are implemented in didp-rs v0.3.2⁴ using Rust 1.65.0. In LNBS and DD-LNS, CABS is run first to find a feasible solution, and then LNBS and DD-LNS are run to improve the solution. For DD-LNS, we use $W = 1000$ and $p = 0.1$ following the original paper [18]. As we described above, LNBS removes conflicting transitions from a suffix, selects the depth using Budgeted-UCB, and geometrically increases the beam width for each neighborhood. To select the starting point of a partial path, we consider two approaches, uniform and biased sampling. While biased sampling achieves better solution quality in CVRP and m-PDTSP, uniform sampling is better in TSPTW, $1 \parallel \sum w_i T_i$, SALBP-1, MOSP, and GCP (see Appendix A). In what follows, we only show results from uniform sampling. In Appendix A, we also evaluate the importance of removing conflicting transitions and Budgeted-UCB. We confirm that these mechanisms significantly improve the solution quality. A more comprehensive ablation study is left for future work.

We implemented the DP models using didppy, a Python interface for didp-rs. We use IBM ILOG CP Optimizer 22.1.0 for the CP models and Gurobi Optimizer 9.5.0 for the MIP models. CP Optimizer is known to use LNS [27]. The DP, CP, and MIP models are implemented in Python 3.10.2. All experiments are run on an Intel Xeon Gold 6148 processor with a single thread, an 8 GB memory limit, and a time limit of 1800 seconds. For LNBS and DD-LNS, we take the median of 5 runs.

Following the previous work [26], we use the primal gap and the primal integral to measure the performance [7]. If an algorithm finds a solution $x^t$ with the cost $f(x^t)$ at time $t$, and the optimal or best-known solution cost is $f^*$, the primal gap at $t$ for the algorithm is

$$p(t) = \begin{cases} 1 & \text{if } f(x^t) \cdot f^* < 0 \\ 0 & \text{if } f(x^t) = f^* = 0 \\ \frac{|f(x^t) - f^*|}{\max\{|f(x^t)|, |f^*|\}} & \text{otherwise.} \end{cases} \quad (8)$$

If the algorithm does not find a solution at time $t$, then $p(t) = 1$. Let $t_1, ..., t_{l-1}$ be time points where a better solution is found, $t_0 = 0$, and $t_l = T$ where $T$ is the time limit. Then, the primal integral, $P(T)$ is defined as

$$P(T) = \sum_{i=1}^{l} p(t_{i-1}) \cdot (t_i - t_{i-1}). \quad (9)$$

We use the primal gap at the time limit, $p(T)$, and the primal integral, $P(T)$, as the measures.

## 4.2 Experimental Results

We show the number of instances where the optimality of a solution is proved, the average primal gap at the time limit, and the average primal integral for each problem in Table 1. We omit MIP in Table 1 because it is outperformed by CP in the primal gap and the primal integral for all problem types (see Appendix A). In the number of optimally solved instances, MIP is the best in CVRP (with 26 problems solved). As reported by previous work [26], MIP is good at finding an optimal solution for small instances, while it fails to find a feasible solution for larger instances, which results in poor performance in the primal gap and the primal integral on average. In other problems, LNBS solves more instances optimally than MIP except for BPP, where LNBS solves 1139 and MIP solves 1157.

<footnote n="4">⁴ https://didp.ai</footnote>
<table><caption>Table 1 Summary of the experimental result. “#” is the number of optimally solved instances, “gap” is the average primal gap at the time limit, and “p.i.” is the average primal integral.</caption><thead><tr><th rowspan="2"></th><th colspan="3">CP</th><th colspan="3">CABS</th><th colspan="3">DD-LNS</th><th colspan="3">LNBS</th></tr><tr><th>#</th><th>gap</th><th>p.i.</th><th>#</th><th>gap</th><th>p.i.</th><th>#</th><th>gap</th><th>p.i.</th><th>#</th><th>gap</th><th>p.i.</th></tr></thead><tbody><tr><td>TSPTW (340)</td><td>47</td><td>0.0259</td><td>49.0</td><td><b>257</b></td><td>0.0033</td><td>9.0</td><td>109</td><td>0.0100</td><td>23.7</td><td>241</td><td><b>0.0016</b></td><td><b>5.7</b></td></tr><tr><td>CVRP (207)</td><td>0</td><td>0.3174</td><td>601.1</td><td><b>6</b></td><td>0.1772</td><td>339.5</td><td>0</td><td>0.2504</td><td>461.6</td><td><b>6</b></td><td><b>0.1640</b></td><td><b>316.8</b></td></tr><tr><td>m-PDTSP (1178)</td><td><b>1049</b></td><td>0.0122</td><td>25.5</td><td>1030</td><td>0.0023</td><td>5.1</td><td>459</td><td>0.0102</td><td>21.7</td><td>1029</td><td><b>0.0022</b></td><td><b>5.0</b></td></tr><tr><td>1||&Sigma;<i>w</i><sub>i</sub><i>T</i><sub>i</sub> (375)</td><td>150</td><td><b>0.0009</b></td><td><b>3.5</b></td><td><b>286</b></td><td>0.0346</td><td>74.4</td><td>100</td><td>0.0409</td><td>83.5</td><td>275</td><td>0.0051</td><td>13.0</td></tr><tr><td>Talent (1000)</td><td>0</td><td>0.0081</td><td>29.3</td><td><b>232</b></td><td>0.0173</td><td>38.2</td><td>0</td><td>0.0602</td><td>114.6</td><td><b>232</b></td><td><b>0.0041</b></td><td><b>11.1</b></td></tr><tr><td>SALBP-1 (2100)</td><td>1584</td><td>0.0046</td><td>28.4</td><td><b>1799</b></td><td><b>0.0003</b></td><td><b>2.2</b></td><td>1507</td><td>0.0067</td><td>13.5</td><td>1682</td><td>0.0022</td><td>7.3</td></tr><tr><td>BPP (1615)</td><td><b>1234</b></td><td><b>0.0014</b></td><td>7.7</td><td>1159</td><td>0.0017</td><td><b>6.1</b></td><td>775</td><td>0.0190</td><td>35.4</td><td>1139</td><td>0.0021</td><td>8.1</td></tr><tr><td>MOSP (570)</td><td>437</td><td>0.0044</td><td>13.0</td><td><b>526</b></td><td><b>0.0000</b></td><td><b>0.4</b></td><td>353</td><td>0.0203</td><td>37.5</td><td>523</td><td>0.0002</td><td>0.7</td></tr><tr><td>GCP (135)</td><td>1</td><td>0.0151</td><td>44.3</td><td><b>103</b></td><td><b>0.0000</b></td><td><b>0.6</b></td><td>3</td><td>0.0009</td><td>2.7</td><td>102</td><td>0.0001</td><td><b>0.6</b></td></tr><tr><td colspan="13">Larger Instances</td></tr><tr><td>m-PDTSP (240)</td><td>77</td><td>0.1491</td><td>285.9</td><td><b>101</b></td><td>0.0694</td><td>153.2</td><td>79</td><td>0.1345</td><td>265.7</td><td>98</td><td><b>0.0652</b></td><td><b>146.6</b></td></tr><tr><td>MOSP (760)</td><td>0</td><td>0.0676</td><td>150.6</td><td><b>150</b></td><td><b>0.0002</b></td><td><b>4.4</b></td><td>0</td><td>0.0402</td><td>72.7</td><td>148</td><td>0.0025</td><td>10.4</td></tr><tr><td>GCP (50)</td><td><b>0</b></td><td>0.5289</td><td>1268.3</td><td><b>0</b></td><td><b>0.0013</b></td><td><b>10.8</b></td><td><b>0</b></td><td>0.0764</td><td>137.8</td><td><b>0</b></td><td>0.0038</td><td>19.5</td></tr></tbody></table>

Compared to CABS, LNBS achieves the better primal gap and the primal integral in TSPTW, CVRP, m-PDTSP, $1||\Sigma w_i T_i$, and Talent. We show the distribution of the primal gap in TSPTW, CVRP, $1||\Sigma w_i T_i$, and Talent in Figure 2. The primal integral has a similar trend to the primal gap in these problems (see Appendix A). In contrast, CABS is better than LNBS in SALBP-1, BPP, MOSP, and GCP. These results are consistent with the observation that LNS is effective for routing and scheduling problems in CP [24, 27]. In the routing problems (TSPTW, CVRP, and m-PDTSP), CABS is already better than CP, and LNBS is even better than CABS. In $1||\Sigma w_i T_i$, while LNBS shows a significant improvement from CABS (0.0051 from 0.0346) outperforming MIP (0.0188), CP is still the best. However, in Talent, LNBS outperforms CP while CABS does not. Overall, LNBS is better than CP in seven problems while CABS is better than CP in six problems.

In TSPTW, the difference in the primal gap between LNBS and CABS (0.0016 and 0.0033) seems small, but it is because they achieve almost the same primal gap in many instances: they optimally solve all 135 instances in the Dumas set [10] and achieve almost the same average primal gap in the AFG set [2], which has 50 instances, and GendreauDumas Extended set [16], which has 130 instances. However, in the OhlmannThomas set [31], which has 25 instances, no instance is optimally solved, and LNBS shows a significant improvement in the primal gap (0.0184 from 0.0395).

In the number of optimally solved instances, CABS is equal to or better than LNBS in all problems. While CABS always searches the entire state space graph, LNBS searches multiple neighborhoods, and the entire state space graph is just one of them. Nevertheless, LNBS proves the optimality of more than 93% of instances that are optimally solved by CABS.

DD-LNS performs worse than CABS and LNBS in all the problems. Previous work has reported that DD-LNS is effective for TSPTW [18]. Note however that the DD-LNS and LNBS results in Table 1 are not the results reported by Gillard and Schaus. To validate that our implementation and experimental settings do not handicap DD-LNS, we compare the results of DD-LNS with those of the original paper in TSPTW.<fnref n="5" /> Our DD-LNS implementation

<footnote n="5">⁵ https://github.com/xgillard/ijcai_22_DDLNS/blob/main/results/tsptw/ddlns/results_w1000_</footnote>
<figure><img src="image_219.png" /><figcaption>Figure 2 Distribution of the primal gap at the time limit. Higher and left is better.</figcaption></figure>

finds a better solution than the original in all instances with the time limit of 600 seconds. This difference is likely due to the difference between the DP models used by us (from Kuroiwa and Beck [26]) and Gillard and Schaus. In TSPTW, a solution is a tour that starts from a depot, visits each customer *j* within the time window [*a*<sub>*j*</sub>, *b*<sub>*j*</sub>], and returns to the depot, and an optimal solution minimizes the total travel time. In both DP models, state variables are the set of unvisited customers *U*, the current customer *i*, and the current time *t*, and each transition corresponds to visiting a customer or the depot. Each DP model has a dual bound (called RLB by Gillard and Schaus), a lower bound on the optimal solution cost. While Gillard and Schaus use a dual bound based on a minimum spanning tree, Kuroiwa and Beck use a simpler one based on the minimum travel time between customers. In addition, Kuroiwa and Beck use information that was not considered by Gillard and Schaus. First, they use dominance between states based on the current time: a state *S* dominates another state *S'* if *S*[*U*] = *S'*[*U*], *S*[*i*] = *S'*[*i*], and *S*[*t*] ≤ *S'*[*t*]. Furthermore, since the time to visit customer *j* is underestimated by *t* + *c*<sub>*ij*</sub><sup>*</sup>, where *c*<sub>*ij*</sub><sup>*</sup> is the shortest travel time from *i* to *j*, they
define state constraints $\forall j \in U, t + c_{ij}^* \le b_j$. The dominance and the state constraints are useful to prune states, which potentially explains the performance gap. Another difference is whether considering the time window constraints at the depot or not. In the benchmark instances used above<fnref n="6" /> (and in Kuroiwa and Beck [26]), a time window $[a_0, b_0]$ is defined for the depot. Gillard and Schaus explicitly model a required return to the depot within $[a_0, b_0]$ while Kuroiwa and Beck do not. However, in the benchmark instances, $b_0 \ge b_j + c_{j0}$ holds for all customers $j$, where $c_{j0}$ is the travel time from $j$ to the depot. Thus, if all customers are visited within the time windows, the depot can be reached within the time window, and explicit modeling of the depot return window is unnecessary.

## 4.3 Larger Instances

LNBS performs better than CABS in m-PDTSP and worse in MOSP and GCP, but the difference in the average primal gap is small. To evaluate the difference more clearly, we use larger instances for these problems.

In m-PDTSP, a vehicle visits all nodes in a graph, picks up some commodities at some nodes, and delivers them to others. Each commodity has a weight and the total weight of commodities that a vehicle can carry is limited by the capacity. In the benchmark set for m-PDTSP, three types of instances are used: Class 1, Class 2, and Class 3 [21], and Class 1 instances are generated from instances of the sequential ordering problem (SOP) [3]. We generate larger Class 1 instances by using 30 SOP instances in TSPLIB<fnref n="7" /> that were not used by the previous work. The original instances have at most 47 nodes, and the new instances have 42 to 378 nodes. We use the same methods as the previous work [21] with the maximum weight $q \in \{1, 5\}$ and the capacity $Q \in \{5q, 10q, 20q, 100q\}$, resulting in 240 instance in total.

In MOSP, an instance is represented by a matrix, and the original set uses at most $125 \times 125$ matrices. We add instances using $150 \times 150$ to $1000 \times 1000$ matrices [8, 12].

In GCP, an instance is represented by a graph. The original instance set uses random and random planar graphs with 20, 30, and 40 nodes. We generate 50 instances using random graphs with 100 and 200 nodes following the method used by previous work [29].

As shown in Table 1, LNBS clearly outperforms CABS in m-PDTSP in the primal gap and the primal integral, but CABS is better in MOSP and GCP. We also show the distribution of the primal gap in m-PDTSP and MOSP in Figure 2. The primal integral has a similar tendency to the primal gap, and the result for GCP is qualitatively similar to that of MOSP.

## 4.4 Analysis of Problem Characteristics

In routing problems, a solution is a route visiting all nodes in a graph, and its cost is the length of the route. In the DP models for these problems, each transition corresponds to visiting one node, and the cost of a partial path increases when a transition is applied. We expect that different partial solutions tend to have different costs, and it is relatively easy to find a better partial path; because the path costs are diverse, unless the current partial path is optimal, better partial paths are included in a partial state space graph with high density. In such a case, beam search is likely to find a better partial path although it searches in a fraction of the partial state space graph restricted by the beam width.

In contrast, in SALBP-1 and BPP, the problem is to pack weighted items into capacitated bins while minimizing the number of bins. In the DP models, each transition packs one item into a bin, and the cost increases only when a new bin is opened. In the DP models for

<footnote n="6">⁶ https://lopez-ibanez.eu/tsptw-instances</footnote>

<footnote n="7">⁷ http://comopt.ifi.uni-heidelberg.de/software/TSPLIB95/sop/</footnote>
<figure><img src="image_220.png" /><figcaption>(a) Problem types where LNBS has lower mean gap: TSPTW, CVRP, m-PDTSP, $1 \parallel \sum w_i T_i$, and Talent.</figcaption></figure>

<figure><img src="image_221.png" /><figcaption>(b) Problem types where CABS has lower mean gap: SALBP-1, BPP, MOSP, and GCP.</figcaption></figure>

**Figure 3** Entropy of the cost distribution over partial paths vs. the solution length in each problem instance. ‘LNBS better’ means LNBS finds a better solution, ‘TIE’ means that LNBS and CABS achieve the same solution cost, and ‘CABS better’ means that CABS finds a better solution.

MOSP and GCP, the cost is computed by taking the maximum weights of edges in a path, and it does not increase unless a new edge has a higher weight than the current maximum. Therefore, we expect that many partial paths tend to have the same cost in these problems, making it difficult to improve a solution by searching only a partial state space graph.

Based on these observations, we hypothesize that LNBS tends to perform better than CABS when path costs in a partial state space graph are diverse. To test this hypothesis, we evaluate the diversity of costs in a partial state space graph using entropy in information theory. Given a solution for a DyPDL model $x = \langle x_1, ..., x_n \rangle$, let $Y_{di}(x)$ be the set of solution paths whose prefix is $\langle x_1, ..., x_{i-1} \rangle$ and the suffix is $\langle x_{i+d}, ..., x_n \rangle$. Let $C = \{ \text{cost}_y(S^0) | y \in Y_{di}(x) \}$ be the set of the path costs. Then, the entropy of the path costs is defined as follows:

$$H(Y_{di}(x)) = - \sum_{c \in C} \frac{|\{y \in Y_{di}(x) | \text{cost}_y(S^0) = c\}|}{|Y_{di}(x)|} \log_2 \frac{|\{y \in Y_{di}(x) | \text{cost}_y(S^0) = c\}|}{|Y_{di}(x)|} \quad (10)$$

As this value gets larger, the cost distribution becomes more diverse, and we expect that LNBS will perform better than CABS. However, even if the entropy is large, if the problem itself is easy, both CABS and LNBS will find optimal or near-optimal solutions. To consider such cases, we also evaluate the length of the initial solution found by CABS.

We evaluate entropy and the length of the initial solution found for each problem instance used in Section 4.2. We first run CABS until it finds a feasible solution and record the length of the solution. Then, we remove the first eight transitions from the solution and enumerate
23:16 Large Neighborhood Beam Search for Domain-Independent Dynamic Programming

all feasible prefixes of the solution, i.e., we use $d = 8$ and $i = 1$. In Figure 3, we show a scatter plot of entropy and the solution length divided into two plots to emphasize the differences between the problem types where LNBS has a lower primal gap on average (Figure 3a) and those where CABS is better (Figure 3b). With low entropy, CABS dominates. For higher entropy, the solution length begins to play a factor: for short solutions, CABS and LNBS perform equally but for longer solutions, LNBS tends to perform better. Indeed, for the problems where CABS performs better on average (Figure 3b), the entropy is quite low (less than 3.5). This result suggests that the entropy of the cost distribution over partial paths is related to the performance of LNBS. Since this analysis is based on path costs, it is not directly applicable to LNS with tree search, where a solution corresponds to a leaf node. However, if we consider the factors of neighborhood size and cost distribution over leaf nodes, we may be able to apply this analysis to LNS for CP and MIP.

5 Related Work

As we discussed, LNBS can be considered a generalization of DD-LNS [18], which combines DDs and LNS. DP is closely related to DDs [23], and DDs have been actively used for combinatorial optimization [9]. For example, DDs are used to obtain bounds on the optimal objective value [4, 33], and heuristics based on DDs have been proposed [6]. Moreover, ddo, a general-purpose DD solver for combinatorial optimization has been developed [5, 19]. In CP, DDs are used for constraint propagation [1, 22]. Recently, HADDOCK, a modeling language of a DD based on a state transition system, was proposed for CP [17].

In state space search, there exist several methods that improve a solution path by searching in a partial state space graph, but they were not framed as LNS. In classical planning, plan neighborhood graph search (PNGS) first constructs a partial state space graph by performing local search from each state in a solution path and then finds the shortest path in the graph [30]. In sliding tile puzzles, iterative tunneling search with A\* (ITSA\*) iteratively expands a partial state space graph, which includes states close to a given path, and finds the shortest path in that graph [13]. Unlike the above two algorithms, Joint and local path A\* (LPA\*) [32] try to find a better partial path between two states in a given path using A\* [20]. While Joint and LPA\* fix the length of a partial path to remove and deterministically select a neighborhood, LNBS dynamically adjusts them and uses beam search instead of A\*.

6 Conclusion

We proposed large neighborhood beam search (LNBS), a state space search algorithm based on large neighborhood search (LNS) and beam search for domain-independent dynamic programming (DIDP). Our configuration of LNBS exploits the multi-armed bandit problem and random sampling to select a neighborhood. We proved that LNBS is complete. LNBS finds better quality solutions on average than the state-of-the-art DIDP solver, complete anytime beam search (CABS), in five out of the nine benchmark problems. In particular, LNBS performs well in routing and scheduling problems, and our analysis suggests that this performance is related to the diversity of the cost distribution over partial paths. A deeper investigation of the characteristics of the problems that make LNS effective in state space search and tree search is an interesting direction for future work. Based on such analysis, developing better configurations for LNBS may also be possible.
References

1. H. R. Andersen, T. Hadzic, J. N. Hooker, and P. Tiedemann. A constraint store based on multivalued decision diagrams. In *Principles and Practice of Constraint Programming – CP 2007*, pages 118–132, 2007. doi:10.1007/978-3-540-74970-7_11.

2. Norbert Ascheuer. Hamiltonian path problems in the on-line optimization of flexible manufacturing systems. PhD thesis, Technische Universität Berlin, 1995.

3. Norbert Ascheuer, Michael Jünger, and Gerhard Reinelt. A branch & cut algorithm for the asymmetric traveling salesman problem with precedence constraint. *Computational Optimization and Applications*, 17:25–42, 2000. doi:10.1023/A:1008779125567.

4. David Bergman, Andre A. Cire, Willem-Jan van Hoeve, and J. N. Hooker. Optimization bounds from binary decision diagrams. *INFORMS Journal on Computing*, 26(2):253–268, 2014. doi:10.1287/ijoc.2013.0561.

5. David Bergman, Andre A. Cire, Willem Jan Van Hoeve, and J. N. Hooker. Discrete optimization with decision diagrams. *INFORMS Journal on Computing*, 28(1):47–66, December 2016. doi:10.1287/ijoc.2015.0648.

6. David Bergman, Andre A. Cire, Willem-Jan van Hoeve, and Tallys Yunes. Bdd-based heuristics for binary optimization. *Journal of Heuristics*, 20(2):211–234, 2014. doi:10.1007/s10732-014-9238-1.

7. Timo Berthold. Measuring the impact of primal heuristics. *Operations Research Letters*, 41:611–614, 2013. doi:10.1016/j.orl.2013.08.007.

8. Marco Antonio Moreira De Carvalho and Nei Yoshihiro Soma. A breadth-first search applied to the minimization of the open stacks. *Journal of the Operational Research Society*, 66:936–946, June 2015. doi:10.1057/jors.2014.60.

9. Margarita P. Castro, Andre A. Cire, and J. Christopher Beck. Decision diagrams for discrete optimization: A survey of recent advances. *INFORMS Journal on Computing*, 34(4):2271–2295, 2022. doi:10.1287/ijoc.2022.1170.

10. Yvan Dumas, Jacques Desrosiers, Eric Gelinas, and Marius M Solomon. An optimal algorithm for the traveling salesman problem with time windows. *Operations Research*, 43(2):367–371, 1995. doi:10.1287/opre.43.2.367.

11. Stefan Edelkamp, Shahid Jabbar, and Alberto Lluch Lafuente. Cost-algebraic heuristic search. In *Proceedings of the 20th National Conference on Artificial Intelligence (AAAI)*, pages 1362–1367, 2005.

12. Rafael de Magalhães Dias Frinhani, Marco Antonio Moreira de Carvalho, and Nei Yoshihiro Soma. A pagerank-based heuristic for the minimization of open stacks problem. *PLoS ONE*, 13(8):1–24, 2018. doi:10.1371/journal.pone.0203076.

13. David A Furcy. ITSA*: Iterative tunneling search with A*. In *Proceedings of AAAI Workshop on Heuristic Search, Memory-Based Heuristics and Their Applications*, pages 21–26, 2006.

14. Maria Garcia de la Banda and Peter J. Stuckey. Dynamic programming to minimize the maximum number of open stacks. *INFORMS Journal on Computing*, 19(4):607–617, 2007. doi:10.1287/ijoc.1060.0205.

15. Maria Garcia de la Banda, Peter J. Stuckey, and Geoffrey Chu. Solving talent scheduling with dynamic programming. *INFORMS Journal on Computing*, 23(1):120–137, 2011. doi:10.1287/ijoc.1090.0378.

16. Michel Gendreau, Alain Hertz, Gilbert Laporte, and Mihnea Stan. A generalized insertion heuristic for the traveling salesman problem with time windows. *Operations Research*, 46(3):330–346, 1998. doi:10.1287/opre.46.3.330.

17. Rebecca Gentzel, Laurent Michel, and W.-J. van Hoeve. HADDOCK: A language and architecture for decision diagram compilation. In Helmut Simonis, editor, *Principles and Practice of Constraint Programming – CP 2020*, pages 531–547, 2020. doi:10.1007/978-3-030-58475-7_31.
18 Xavier Gillard and Pierre Schaus. Large neighborhood search with decision diagrams. In *Proceedings of the 21st International Joint Conference on Artificial Intelligence, IJCAI-22*, pages 4754-4760, 2022. doi:10.24963/ijcai.2022/659.

19 Xavier Gillard, Pierre Schaus, and Vianney Coppé. Ddo, a generic and efficient framework for mdd-based optimization. In *Proceedings of the 29th International Joint Conference on Artificial Intelligence, IJCAI-20*, pages 5243-5245, 2020. doi:10.24963/ijcai.2020/757.

20 Peter E. Hart, Nils J. Nilsson, and Bertram Raphael. A formal basis for the heuristic determination of minimum cost paths. *IEEE Transactions on Systems Science and Cybernetics*, 4(2):100-107, 1968. doi:10.1109/TSSC.1968.300136.

21 Hipólito Hernández-Pérez and Juan José Salazar-González. The multi-commodity one-to-one pickup-and-delivery traveling salesman problem. *European Journal of Operational Research*, 196:987-995, August 2009. doi:10.1016/j.ejor.2008.05.009.

22 Samid Hoda, Willem-Jan van Hoeve, and J. N. Hooker. A systematic approach to MDD-based constraint programming. In *Principles and Practice of Constraint Programming – CP 2010*, pages 266-280, 2010.

23 John N. Hooker. Decision diagrams and dynamic programming. In Carla Gomes and Meinolf Sellmann, editors, *Integration of AI and OR Techniques in Constraint Programming for Combinatorial Optimization Problems*, pages 94-110, 2013.

24 Siddhartha Jain and Pascal Van Hentenryck. Large neighborhood search for dial-a-ride problems. In *Principles and Practice of Constraint Programming – CP 2011*, pages 400-413, 2011.

25 Ryo Kuroiwa and J. Christopher Beck. Domain-independent dynamic programming: Generic state space search for combinatorial optimization. In *Proceedings of the 33rd International Conference on Automated Planning and Scheduling (ICAPS)*, 2023. doi:10.1609/icaps.v33i1.27200.

26 Ryo Kuroiwa and J. Christopher Beck. Solving domain-independent dynamic programming problems with anytime heuristic search. In *Proceedings of the 33rd International Conference on Automated Planning and Scheduling (ICAPS)*, 2023. doi:10.1609/icaps.v33i1.27201.

27 Philippe Laborie, Jérôme Rogerie, Paul Shaw, and Petr Vilím. IBM ILOG CP optimizer for scheduling. *Constraints*, 23(2):210-250, 2018. doi:10.1007/s10601-018-9281-x.

28 Shu Lin, Na Meng, and Wenxin Li. Optimizing constraint solving via dynamic programming. In *Proceedings of the 28th International Joint Conference on Artificial Intelligence, IJCAI-19*, pages 1146-1154, July 2019. doi:10.24963/ijcai.2019/160.

29 Michael Morin, Margarita P. Castro, Kyle E.C. Booth, Tony T. Tran, Chang Liu, and J. Christopher Beck. Intruder alert! Optimization models for solving the mobile robot graph-clear problem. *Constraints*, 23(3):335-354, 2018. doi:10.1007/s10601-018-9288-3.

30 Hootan Nakhost and Martin Müller. Action elimination and plan neighborhood graph search: Two algorithms for plan improvement. In *Proceedings of the 20th International Conference on Automated Planning and Scheduling (ICAPS)*, pages 121-128, 2010. doi:10.1609/icaps.v20i1.13402.

31 Jeffrey W. Ohlmann and Barrett W. Thomas. A compressed-annealing heuristic for the traveling salesman problem with time windows. *INFORMS Journal on Computing*, 19(1):80-90, 2007. doi:10.1287/ijoc.1050.0145.

32 Daniel Ratner and Ira Pohl. Joint and LPA*: Combination of approximation and search. In *Proceedings of the fifth National Conference on Artificial Intelligence (AAAI)*., pages 173-177, 1986. doi:10.5555/2887770.2887798.

33 Isaac Rudich, Quentin Cappart, and Louis-Martin Rousseau. Peel-And-Bound: Generating Stronger Relaxed Bounds with Multivalued Decision Diagrams. In *Principles and Practice of Constraint Programming – CP 2022*, pages 35:1-35:20, 2022. doi:10.4230/LIPIcs.CP.2022.35.
34 Paul Shaw. Using constraint programming and local search methods to solve vehicle routing problems. In *Principles and Practice of Constraint Programming – CP98*, volume 1520, pages 417–431, 1998. doi:10.1007/3-540-49481-2_30.

35 Yingce Xia, Xu-Dong Zhang, Nenghai Yu, Geoffrey Holmes, and Yan Liu. Budgeted bandit problems with continuous random costs. In *Proceedings of the Seventh Asian Conference on Machine Learning*, pages 317–332, 2015.

36 Weixiong Zhang. Complete anytime beam search. In *Proceedings of the Fifteenth National/Tenth Conference on Artificial Intelligence/Innovative Applications of Artificial Intelligence (AAAI-98/IAAI-98)*, pages 425–430, 1998.

A Additional Experimental Results

■ Table 2 Comparison of CP, MIP, and two configurations of LNBS. “#” is the number of optimally solved instances, “gap” is the average primal gap at the time limit, and “p.i.” is the average primal integral.

<table><thead><tr><th rowspan="2"></th><th colspan="3">CP</th><th colspan="3">MIP</th><th colspan="3">LNBS</th><th colspan="3">LNBS/bias</th></tr><tr><th>#</th><th>gap</th><th>p.i.</th><th>#</th><th>gap</th><th>p.i.</th><th>#</th><th>gap</th><th>p.i.</th><th>#</th><th>gap</th><th>p.i.</th></tr></thead><tbody><tr><td>TSPTW (340)</td><td>46</td><td>0.0275</td><td>52.3</td><td>227</td><td>0.2268</td><td>484.1</td><td>241</td><td><strong>0.0016</strong></td><td><strong>5.7</strong></td><td><strong>242</strong></td><td>0.0019</td><td>6.2</td></tr><tr><td>CVRP (207)</td><td>0</td><td>0.3174</td><td>601.1</td><td><strong>26</strong></td><td>0.5845</td><td>1157.4</td><td>6</td><td>0.1640</td><td>316.8</td><td>6</td><td><strong>0.1633</strong></td><td><strong>314.8</strong></td></tr><tr><td>m-PDTSP (1178)</td><td><strong>1050</strong></td><td>0.0121</td><td>25.4</td><td>945</td><td>0.0858</td><td>180.0</td><td>1029</td><td>0.0022</td><td>5.0</td><td>1029</td><td><strong>0.0021</strong></td><td><strong>4.7</strong></td></tr><tr><td>1||&Sigma;<i>w</i><sub>i</sub><i>T</i><sub>i</sub> (375)</td><td>150</td><td><strong>0.0003</strong></td><td><strong>2.4</strong></td><td>109</td><td>0.0188</td><td>75.6</td><td><strong>275</strong></td><td>0.0051</td><td>13.0</td><td><strong>275</strong></td><td>0.0056</td><td>14.5</td></tr><tr><td>Talent (1000)</td><td>7</td><td>0.0072</td><td>27.6</td><td>0</td><td>0.0573</td><td>152.6</td><td><strong>232</strong></td><td><strong>0.0041</strong></td><td><strong>11.1</strong></td><td>231</td><td>0.0042</td><td><strong>11.1</strong></td></tr><tr><td>SALBP-1 (2100)</td><td>1584</td><td>0.0046</td><td>28.4</td><td>1357</td><td>0.3447</td><td>634.6</td><td><strong>1682</strong></td><td><strong>0.0022</strong></td><td><strong>7.3</strong></td><td>1675</td><td><strong>0.0022</strong></td><td>7.5</td></tr><tr><td>BPP (1615)</td><td><strong>1234</strong></td><td><strong>0.0014</strong></td><td><strong>7.7</strong></td><td>1157</td><td>0.0385</td><td>85.9</td><td>1139</td><td>0.0021</td><td>8.1</td><td>1129</td><td>0.0021</td><td>9.0</td></tr><tr><td>MOSP (570)</td><td>437</td><td>0.0044</td><td>13.0</td><td>224</td><td>0.0394</td><td>100.4</td><td><strong>523</strong></td><td><strong>0.0002</strong></td><td><strong>0.7</strong></td><td><strong>523</strong></td><td><strong>0.0002</strong></td><td><strong>0.7</strong></td></tr><tr><td>GCP (135)</td><td>1</td><td>0.0151</td><td>44.3</td><td>23</td><td>0.1102</td><td>311.9</td><td><strong>102</strong></td><td><strong>0.0001</strong></td><td><strong>0.6</strong></td><td><strong>102</strong></td><td><strong>0.0001</strong></td><td>0.7</td></tr><tr><td colspan="14">Larger Instances</td></tr><tr><td>m-PDTSP (240)</td><td>77</td><td>0.1481</td><td>284.1</td><td>47</td><td>0.5811</td><td>1096.8</td><td><strong>98</strong></td><td>0.0652</td><td><strong>146.6</strong></td><td>97</td><td><strong>0.0647</strong></td><td>147.0</td></tr><tr><td>MOSP (760)</td><td>0</td><td>0.0675</td><td>150.4</td><td>0</td><td>0.8806</td><td>1599.4</td><td><strong>148</strong></td><td><strong>0.0025</strong></td><td><strong>10.4</strong></td><td><strong>148</strong></td><td>0.0027</td><td>10.7</td></tr><tr><td>GCP (50)</td><td><strong>0</strong></td><td>0.5287</td><td>1268.1</td><td><strong>0</strong></td><td>0.5306</td><td>977.8</td><td><strong>0</strong></td><td><strong>0.0038</strong></td><td><strong>19.5</strong></td><td><strong>0</strong></td><td>0.0061</td><td>21.2</td></tr></tbody></table>

■ Table 3 Comparison of LNBS variants. “No removing conflicts” does not remove conflicting transitions in the suffix. “No Budgeted-UCB” selects the depth uniformly at random instead of using Budgeted-UCB. “#” is the number of optimally solved instances, “gap” is the average primal gap at the time limit, and “p.i.” is the average primal integral.

<table><thead><tr><th rowspan="2"></th><th colspan="3">LNBS</th><th colspan="3">No removing conflicts</th><th colspan="3">No Budgeted-UCB</th></tr><tr><th>#</th><th>gap.</th><th>p.i.</th><th>#</th><th>gap.</th><th>p.i.</th><th>#</th><th>gap.</th><th>p.i.</th></tr></thead><tbody><tr><td>TSPTW (340)</td><td>241</td><td><strong>0.0016</strong></td><td><strong>5.7</strong></td><td>234</td><td>0.0035</td><td>10.6</td><td><strong>256</strong></td><td>0.0034</td><td>9.4</td></tr><tr><td>CVRP (207)</td><td><strong>6</strong></td><td><strong>0.1640</strong></td><td><strong>316.8</strong></td><td>5</td><td>0.1682</td><td>322.5</td><td><strong>6</strong></td><td>0.1767</td><td>338.7</td></tr><tr><td>1||&Sigma;<i>w</i><sub>i</sub><i>T</i><sub>i</sub> (375)</td><td>275</td><td><strong>0.0051</strong></td><td><strong>13.0</strong></td><td>268</td><td>0.0224</td><td>56.0</td><td><strong>287</strong></td><td>0.0340</td><td>74.1</td></tr></tbody></table>

We show the result of MIP in Table 2. MIP solves more instances than CP in TSPTW, CVRP, and GCP, but CP is better in the primal gap and the primal integral.

In addition, we compare two configurations of LNBS in Table 2. One selects the starting point of a partial path uniformly at random (LNBS), and another selects it according to the probability distribution biased by partial path costs in Equation (7) (LNBS/bias). LNBS/bias solves one more instance in TSPTW and outperforms LNBS in CVRP and m-PDTSP in the primal gap.
<figure><img src="image_222.png" /><figcaption>Figure 4 Distribution of the primal gap at the time limit. Higher and left is better.</figcaption></figure>

We evaluate the importance of other components of LNBS using a subset of the problems: TSPTW, CVRP, and $1 \ll \sum w_i T_i$. In Table 3, we compare two variants of LNBS, where conflicting transitions in a suffix are not removed ('No removing conflicts'), and the depth is selected uniformly at random instead of Budgeted-UCB ('No Budgeted-UCB'). The two variants perform worse in terms of the primal gap and the primal integral. While 'No Budgeted-UCB' solves more instances to optimality, it is because the largest depth, which makes LNBS the same as CABS, is more likely to be selected by uniform sampling. Indeed, the primal gap and the primal integral of 'No Budgeted-UCB' are close to those of CABS.

In Figure 4, we present the distribution of the primal gap over instances in m-PDTSP, SALBP-1, BPP, MOSP, Graph-Clear, and large instances in GCP. LNBS is slightly better in m-PDTSP, and CABS is better in SALBP-1, BPP, and large instances in GCP. Figures 5 and 6 show the distribution of the primal integral. The tendency is similar to that of the primal gap.
<figure><img src="image_223.png" /><figcaption>(a) TSPTW.</figcaption></figure>

<figure><img src="image_224.png" /><figcaption>(b) CVRP.</figcaption></figure>

<figure><img src="image_225.png" /><figcaption>(c) m-PDTSP.</figcaption></figure>

<figure><img src="image_226.png" /><figcaption>(d) Large instances in m-PDTSP.</figcaption></figure>

<figure><img src="image_227.png" /><figcaption>(e) $1 \parallel \sum w_i T_i$.</figcaption></figure>

<caption>(f) Talent.</caption>

■ **Figure 5** Distribution of the primal integral at the time limit in the problems where LNBS is better. Higher and left is better.
<figure><img src="image_228.png" /><figcaption>Figure 6 Distribution of the primal integral at the time limit in the problems where CABS is better. Higher and left is better.</figcaption></figure>
# MDD Archive for Boosting the Pareto Constraint

**Steve Malalel**

Université Côte d'Azur, CNRS, I3S, Nice, France

**Arnaud Malapert**

Université Côte d'Azur, CNRS, I3S, Nice, France

**Marie Pelleau**

Université Côte d'Azur, CNRS, I3S, Nice, France

**Jean-Charles Régin**

Université Côte d'Azur, CNRS, I3S, Nice, France

## Abstract

Multi-objective problems are frequent in the real world. In general they involve several incomparable objectives and the goal is to find a set of Pareto optimal solutions, i.e. solutions that are incomparable two by two. In order to better deal with these problems in CP the global constraint PARETO was developed by Schaus and Hartert to handle the relations between the objective variables and the current set of Pareto optimal solutions, called the archive. This constraint handles three operations: adding a new solution to the archive, removing solutions from the archive that are dominated by a new solution, and reducing the bounds of the objective variables. The complexity of these operations depends on the size of the archive. In this paper, we propose to use a multi-valued Decision Diagram (MDD) to represent the archive of Pareto optimal solutions. MDDs are a compressed representation of solution sets, which allows us to obtain a compressed and therefore smaller archive. We introduce several algorithms to implement the above operations on compressed archives with a complexity depending on the size of the archive. We show experimentally on bin packing and multi-knapsack problems the validity of our approach.

**2012 ACM Subject Classification** Applied computing → Multi-criterion optimization and decision-making; Theory of computation → Constraint and logic programming; Mathematics of computing → Decision diagrams

**Keywords and phrases** Constraint Programming, Global Constraint, MDD, Multi-Objective Problem, Pareto Constraint

**Digital Object Identifier** 10.4230/LIPIcs.CP.2023.24

**Funding** This work has been supported by the French government, through the 3IA Côte d'Azur Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-19-P3IA-0002.

## 1 Introduction

Multi-objective combinatorial optimization (MOCO) problems are present in many industrial applications [13, 14]. They involve several incomparable objectives represented by objective variables.

For the sake of clarity and without loss of generality, we will consider that we have to solve a problem where all objective variables must be minimized.

A solution $S_1$ of objective variables is dominated by another solution $S_2$ if for each objective variable $obj_i$ the value of $obj_i$ in $S_2$ is better than or equal to the value of $obj_i$ in $S_1$. For instance the solution (4, 6, 3, 1) is dominated by (4, 3, 2, 1) but it is not dominated by (1, 1, 1, 3). The set of non dominated solutions defines the set of Pareto optimal solutions. In MOCO, the goal is to compute that set of Pareto optimal solutions. Usually the set of non

<figure><img src="image_229.png" /></figure>

© Steve Malalel, Arnaud Malapert, Marie Pelleau, and Jean-Charles Régin;

licensed under Creative Commons License CC-BY 4.0

29th International Conference on Principles and Practice of Constraint Programming (CP 2023).

Editor: Roland H. C. Yap; Article No. 24; pp. 24:1–24:15

<figure><img src="image_230.png" /></figure>

LIPICS

Leibniz International Proceedings in Informatics

Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl Publishing, Germany
dominated solutions are saved in an archive that is maintained during the search for solutions. Two operations are involved: *insert* that manages the addition of a new non dominated solution and *delete* that removes the solutions that are dominated by the new solution.

In addition, solvers dealing with MOCO problems have to deal with another question: Is a new solution dominated by an archive solution?

In Constraint Programming, we answer this question by avoiding the generation of dominated solutions. To do this, we add a constraint to the problem that ensures that no dominated solution can be computed. This added constraint is called Pareto constraint and was proposed by Schaus and Hartert [11, 4]. It implements the ideas of Gavanelli [3]. This constraint reduces the bounds of the objective variables such that dominated solutions cannot be produced. This result is obtained by preventing a new solution from being dominated by a solution from the archive. In order to understand this process, let us create a tuple composed of the current minimum of all objective variables. This tuple will dominate all solutions that can be constructed from the current objective variables. Thus, if there is a solution in the archive that dominates this tuple then clearly it will dominate all future solutions and so we can stop the search. Now, consider the objective variable $obj_i$. If we replace in our tuple the value of $obj_i$ by its maximum possible value and if we found $S$, a solution of the archive, dominating this tuple, then $obj_i$ must take a value less than or equal to that of $S$ otherwise the future solutions will be dominated by $S$. By applying this process for each variable, Schaus and Hartert establish the bound consistency of the constraint. We will denote by *filter* this process. The complexity of this operation is not detailed in their paper. A simple implementation will require to traverse *n* (the number of objectives) times the archive. It is not straightforward to reach a time complexity linear in the size of the archive.

It may be tempting to use multi-valued decision diagrams (MDDs) for this constraint because MDDs are a compressive data structure for representing solution sets. Perez [7] defined the MDD representing the Pareto constraint, i.e. the set of tuples allowed by the constraint that are the possible future non-dominated solutions. As mentioned by Perez, this approach failed mainly because at the beginning the MDD compresses very strongly the set of solutions since everything is almost possible, then it will decompress because we only delete tuples and the chances of recompression are low.

In this article we propose to use MDDs to represent the archive, that are the current non dominated solutions. Currently the archive is often represented as lists. We can also use quad-trees but this is only efficient if we have few objectives [6], which is not our case of study. With lists (or quad-trees for that matter), the complexity of the operations *insert* and *delete* is linear with the size of the archive (i.e. the number of elements multiplied by the number of objectives). Representing the archive by an MDD will allow parts of common solutions to be merged. As with an MDD, all the solutions are treated globally, we can therefore hope to save time thanks to these groupings.

The operation *delete* will therefore potentially save time. The operation *insert* may be slower because MDDs are a heavier data structure than lists. However, this operation takes much less time than *delete* or *filter*. For this last operation, we will benefit from the global view of the MDD. We propose to improve the algorithm of Schaus and Hartet in two ways: we define an algorithm to find the tightest solutions (i.e. the largest possible value of an objective) faster and we introduce a variant of this algorithm that processes all objective variables at once, and not successively.

Our algorithms are based on the following idea: Consider the objective variable $obj_i$. Let us remove from the MDD that represents the archive all the values of the objective variables different from $obj_i$ that are less than or equal to the minimum of their variable. Then we
perform a reduction of the MDD in order to obtain MDD<sub>D</sub>. If MDD<sub>D</sub> is not empty then it means that there exist paths from root to tt in the MDD, that is solutions in the archive which will dominate any new solution involving some values of obj<sub>i</sub>. More precisely, there exist solutions of the form (v<sub>1</sub>, v<sub>2</sub>, ..., v<sub>i</sub>, ..., v<sub>n</sub>) such that $\forall j = 1...n, j \neq i : v_j \leq \min(obj_j)$. These solutions dominate (or are equal to) any future solution with obj<sub>i</sub> ≥ v<sub>i</sub>. Hence, the maximum possible value for obj<sub>i</sub> is the smallest value of obj<sub>i</sub> in MDD<sub>D</sub>. The obtained algorithms have a linear time complexity in the size of the MDD, which improves the algorithm of Schaus and Hartet.

The paper is organized as follows. First, we recall some concepts and definitions of Constraint Programming, Multi-Objective Optimization Problems and Multi-valued Decision Diagrams. Then, we introduce the representation of the archive by an MDD. We present how the *insert* and *delete* operations are implemented. We detail two algorithms for the operation *filter* that establishing the bound consistency of the Pareto Constraint associated with an MDD. Next, we experiment with these methods on bin packing and multi-knapsack problems. At last, we conclude.

# 2 Preliminaries

## 2.1 Constraint Programming

A finite constraint network $\mathcal{N}$ is defined as a set of $n$ variables $X = \{x_1, ..., x_n\}$, a set of current domains $\mathcal{D} = \{D(x_1), ..., D(x_n)\}$ where $D(x_i)$ is the finite set of possible values for variable $x_i$, and a set $\mathcal{C}$ of constraints between variables. If $x$ is a variable, then $x^{\min} = \min(D(x))$ and $x^{\max} = \max(D(x))$. We introduce the particular notation $\mathcal{D}_0 = \{D_0(x_1), ..., D_0(x_n)\}$ to represent the set of initial domains of $\mathcal{N}$ on which constraint definitions were stated. An element of $D_0(x_1) \times \dots \times D_0(x_n)$ on the ordered set $\mathcal{D}$ is called a tuple and is denoted $\tau$. In a tuple $\tau$, the assignment of the $i^{th}$ variable is denoted $\tau_i$. A solution of $\mathcal{N}$ is a tuple $\tau$ that satisfies all the constraints in $\mathcal{C}$. A constraint $C$ on the ordered set of variables $X(C) = (x_{i_1}, ..., x_{i_r})$ is a subset $T(C)$ of the Cartesian product $D_0(x_{i_1}) \times \dots \times D_0(x_{i_r})$ that specifies the allowed combinations of values for the variables $x_{i_1}, ..., x_{i_r}$.

## 2.2 Multi-Objective Optimization

A multi-objective problem in combinatorial optimization is a problem where several objectives have to be improved while satisfying constraints. For the sake of clarity and without loss of generality, these objectives are represented by integer variables and have to be minimized. We denote by $O = (obj_1, ..., obj_m)$ the ordered set of the objective variables in $X$, the set of variables of the whole problem. Then, this problem can be modeled as follows:

$$ \begin{array}{ll} \text{Minimize} & O \\ \text{Subject to} & \mathcal{C} \end{array} \qquad (1) $$

However, the minimization of several objectives simultaneously may seem ambiguous. Indeed, there is no order of priority between the different objectives and improving one often means degrading at least one of the others. This generally introduces the need to make compromises during the solving process: the goal is not to find only one optimal solution but a set of solutions that are considered *Pareto optimal*.

In the rest of the paper, we will focus only on the objective variables. As a consequence, a tuple will always be understood only in relation to the set $O$, and the same applies to a solution. Thus, for a tuple $\tau$, the assignment of the $i^{th}$ objective variable is denoted by $\tau_i$.
24:4 MDD Archive for Boosting the Pareto Constraint

The following definitions are taken from [11]:

**Definition 1 (Pareto dominance).** Let $\tau$ and $\tau'$ be two solutions of a multi-objective problem represented by a constraint network $\mathcal{N}$.
We say that $\tau$ dominates $\tau'$, denoted $\tau \prec \tau'$, if and only if:

$$ \begin{array}{l} \forall i \in [1 \dots m] : \tau_i \le \tau'_i \\ \wedge \quad \exists i \in [1 \dots m] : \tau_i < \tau'_i \end{array} \qquad (2) $$

We say that $\tau$ weakly-dominates $\tau'$, denoted $\tau \preceq \tau'$, if and only if $\forall i \in [1 \dots m] : \tau_i \le \tau'_i$.

**Definition 2 (Pareto optimality).** Let $\mathcal{S}$ be the set of all the feasible solutions of a multi-objective problem represented by a constraint network $\mathcal{N}$. A solution $\tau$ is Pareto optimal if and only if there is no solution $\tau'$ in $\mathcal{S}$ that dominates $\tau$:

$$ \mathcal{H}\tau' \in \mathcal{S} : \tau' \prec \tau \qquad (3) $$

**Definition 3 (Pareto set).** Let $\mathcal{N}$ be the constraint network representing a multi-objective problem and $\mathcal{S}$ be the set of all the feasible solutions of $\mathcal{N}$. The Pareto set of $\mathcal{N}$ is the set of all the Pareto optimal solutions in $\mathcal{S}$:

$$ \{\tau \in \mathcal{S} | \mathcal{H}\tau' \in \mathcal{S} : \tau' \prec \tau\} \qquad (4) $$

The search for the exact Pareto set of a multi-objective problem can be impossible to achieve in a reasonable time. This leads to search for an approximation of the Pareto set: the *archive*.

**Definition 4 (Archive).** An archive $\mathcal{A}$ is a set of solutions such that there is no solution $\tau'$ in the archive that dominates another solution $\tau$ in the archive. This property is known as the *domination-free property*:

$$ \tau \in \mathcal{A}, \mathcal{H}\tau' \in \mathcal{A} : \tau' \prec \tau \qquad (5) $$

A basic way to maintain an archive $\mathcal{A}$ is to verify if a solution $\tau$ found during the search is dominated by a solution of the archive:

- If $\tau$ is dominated by at least one solution, do not add it in $\mathcal{A}$.

- If $\tau$ is not dominated by any solution, add it and remove from $\mathcal{A}$ all the solutions dominated by $\tau$.

## 2.3 Pareto Constraint

We reformulate the definition introduced in [11] in order to avoid the notion of “next discovered solution”.

**Definition 5 (Pareto Constraint).** Let $X$ be a set of objective variables and $\mathcal{A}$ be an archive defined on $O$. A Pareto constraint is a constraint $C$ associated with $\mathcal{A}$ defined by $\text{PARETO}(O, \mathcal{A}) = \{\tau \text{ s.t. } \tau \text{ is a tuple on } O \text{ and } \mathcal{H}\tau' \in \mathcal{A} \text{ with } \tau' \preceq \tau\}$.

When using this constraint during the search for solutions, newly found solutions must be inserted in the archive. One could notice that this constraint prevents finding a solution $\tau$ such that $\tau' \in \mathcal{A}$ and $\tau = \tau'$.
$$ \begin{array}{c} D(obj_1) = \{2, 4\} \\ D(obj_2) = \{1, 2, 5\} \\ D(obj_3) = \{1, 3, 4, 5\} \\ D(obj_4) = \{2, 5, 6\} \\ \wedge \\ \boxed{\phantom{A}} \\ \begin{array}{c} (3, 1, 3, 1) \\ (2, 1, 4, 2) \end{array} \end{array} \implies obj_3 < 4 $$

<caption>Figure 1 Application of the propagator of the PARETO constraint.</caption>

We adapt the definition of ideal point of multi-objective problems to our purpose:

**Definition 6 (Ideal tuple).** Let $C = \text{PARETO}(O, \mathcal{A})$ be a Pareto constraint with $O = (obj_1, \dots, obj_n)$.

- The ideal tuple of $C$ denoted by $\tau^*(O)$ is the tuple composed of the best objective values, that is $(obj_1^{\min}, \dots, obj_n^{\min})$.

- The ideal tuple for the value $a$ of the variable $obj_i$ denoted by $\tau^*(O, i, a)$ is the tuple composed of $obj_i = a$ and the best objective values for the other objective variables, that is $(obj_1^{\min}, \dots, obj_{i-1}^{\min}, a, obj_{i+1}^{\min}, \dots, obj_n^{\min})$.

**Proposition 7.** Let $C = \text{PARETO}(O, \mathcal{A})$ be a Pareto constraint; the following two properties are equivalent:

- $C$ is consistent;

- $\tau^*(O)$ is not weakly-dominated by any tuple of $\mathcal{A}$

*Proof.* By definition $\tau^*(O)$ weakly-dominates any tuple defined on $O$, thus if $\tau^*(O)$ is not weakly-dominated then it is a possible solution and the constraint is consistent. Otherwise, there is no solution and $C$ is not consistent. ◀

We reformulate the filtering algorithm associated with the Pareto constraint given in [11].

**Proposition 8.** Let $C = \text{PARETO}(O, \mathcal{A})$ be a Pareto constraint. The value $a$ of the objective variable $obj_i$ is not consistent with $C$ if and only if $\exists \tau \in \mathcal{A}$ such that $\tau \preceq \tau^*(O, i, a)$.

*Proof.* ⇒ If the value $a$ of $obj_i$ is not consistent with $C$ then every tuple $\tau$ of $C$ with $\tau_i = a$ is weakly-dominated by a tuple of $\mathcal{A}$. Therefore $\tau^*(O, i, a)$ is weakly-dominated by a tuple of $\mathcal{A}$.

⇔ The tuple $\tau^*(O, i, a)$ weakly-dominates all the possible tuples $\tau$ of $C$ with $\tau_i = a$. Thus if this tuple is weakly-dominated then there is no tuple with $\tau_i = a$ consistent with the constraint and the value $a$ of $obj_i$ is not consistent with $C$. ◀

From this proposition we can identify all values of all variables that are inconsistent with the constraint and so we can establish the arc consistency of the constraint which is equivalent in our case to the bound consistency. Figure 1 gives an example of domain reduction.

## 2.4 Multi-valued Decision Diagram

The decision diagrams considered in this paper are reduced ordered multi-valued decision diagrams (MDD) [5, 12, 1], which are a generalization of binary decision diagrams [2]. They use a fixed variable ordering for canonical representation and shared sub-graphs for compression obtained by means of a reduction operation. An MDD is a rooted directed acyclic graph (DAG) used to represent some multi-valued functions $f: \{0 \dots d-1\}^n \to \text{true, false}$.
<figure><img src="image_231.png" /><figcaption>Figure 2 An MDD representing the tuple set {((1, 1), (1, 3), (2, 1), (2, 3), (3, 1), (3, 2))}.</figcaption></figure>

Given the *n* input variables, the DAG contains *n* + 1 layers of nodes, such that each variable is represented at a specific layer of the graph. Each node on a given layer has at most *d* outgoing arcs to nodes in the next layer of the graph. Each arc is labeled by its corresponding integer. The arc (*u*, *a*, *v*) is from node *u* to node *v* and labeled by *a*. Sometimes it is convenient to say that *v* is a child of *u*. The set of outgoing arcs from node *u* is denoted by $\omega^+(u)$. All outgoing arcs of the layer *n* reach *tt*, the true terminal node (the false terminal node is typically omitted). There is an equivalence between $f(a_1, \dots, a_n) = \text{true}$ and the existence of a path from the root node to the *tt* whose arcs are labeled $a_1, \dots, a_n$. Figure 2 shows an example of MDD and the kind of compression it can offer.

The reduction of an MDD is one of the most important operations that may reduce the MDD size by an exponential factor. It consists in removing nodes that have no successor and merging equivalent nodes, i.e., nodes having the same set of neighbors associated with the same labels. This means that only nodes of the same layer can be merged. Other operations used in this paper are the addition and deletion of tuples of an MDD. They can be performed with in-place operations provided by Perez and Régin [8].

The advantage of using MDDs instead of the usual data structures is their compression capability which is useful for reducing memory consumption. Moreover, this compression may also improve the time performance of algorithms computing on a set.

## 3 Pareto Constraint Using MDD

We propose to use $MDD_A$, an MDD, to represent the solution archive. Consider $\tau$ a new solution. Without loss of generality, we assume that $\tau$ is not weakly-dominated by any tuple of the archive. As mentioned in the introduction we need to implement the operations *insert* and *delete*:

- *insert*: $\tau$ has to be added to $MDD_A$.

- *delete*: All the tuples dominated by $\tau$ must be removed from $MDD_A$.

In addition we need to design algorithms for implementing the *filter* operation of the *PARETO* constraint.

### 3.1 Insert and Delete Operations

Adding $\tau$ to $MDD_A$ can be done thanks to the in-place addition operation [8]. The deletion from $MDD_A$ of all the tuples that are dominated by $\tau$ can be done by creating $MDD_{dom}(\tau)$, the MDD of all the tuples weakly-dominated by $\tau$. $MDD_{dom}(\tau)$ is really simple: each layer contains only one node, and for each layer $i \in [1... (n-1)]$ there are arcs labeled with all the values that belong to the range [$\tau_i...max(D_0(obj_i))$] between the node of layer *i* and the node of layer *i* + 1. Figure 3 shows an example of such an MDD. Then, the operation $MDD_A$
<figure><img src="image_232.png" /><figcaption>Figure 3 An MDD representing all the tuples weakly-dominated by (2, 5, 1) with $D_0(obj_1) = [1...5]$, $D_0(obj_2) = [1...6]$ and $D_0(obj_3) = [1...3]$.</figcaption></figure>

MDD<sub>dom</sub>(τ) can be performed in-place thanks to the in-place difference operator [8]. It should be noted that this operation also deletes the tuple τ from MDD<sub>A</sub>, so it is better to perform first the *delete* operation and then the *insert* operation.

## 3.2 Filtering algorithm of the Pareto Constraint

Let $C = \text{PARETO}(O, \text{MDD}_A)$ be a Pareto constraint whose archive is represented by an MDD. We present two methods that eliminate all values of the variables satisfying Proposition 8 (i.e., that are not consistent with $C$). That establishes the bound consistency of $C$. The first one has to be executed for each objective variable, and the second one uses the concept of the first method to filter all the objective variables at the same time.

### 3.2.1 Unidirectional Marking

This method operates only on one objective variable $obj_i$ at a time. It must be repeated for each objective variable to be complete.

MDD<sub>A</sub> is traversed using a Depth-First Search (DFS) procedure from *root* following two rules:

- For each layer $j \neq i$ it is possible to go only through arcs whose values are less than or equal to $obj_j^{min}$.

- For the layer $i$, it is possible to go only through arcs whose values are less than or equal to $obj_i^{max}$.

Each time the node $tt$ is reached, the value of the arc of the layer $i$ belonging to the current path is memorized if it is less than the previous memorized value. Then, $obj_i^{max}$ takes the lower value between the current upper bound and the memorized value minus one. Figure 4 shows two examples based on the same archive with different states for the domains. Algorithm 3.1 is a possible implementation of this filtering for all the objective variables.

**Proposition 9.** The unidirectional marking method eliminates all the values of $obj_i$ that are not consistent (c.f. Proposition 8).

*Proof.* The unidirectional marking method finds paths corresponding to the tuples that weakly-dominant $\tau^*(O, i, obj_i^{max})$ and it eliminates all the values $a$ such that $\tau^*(O, i, a)$ is weakly-dominated, by setting the maximum of $obj_i$ to $minV - 1$ where $minV = \min\{a$ such that $\tau^*(O, i, a)$ is weakly-dominated\}$. ◀

The time complexity of this method for one objective variable is linear in the size of the MDD, because it traverses the MDD with a DFS. However, as it is repeated for each objective variable and the size of the MDD depends on the number of objective variables, the overall
<figure><img src="image_233.png" /><figcaption>■ **Figure 4** Application of the unidirectional marking in $MDD_A$ for the objective variable $obj_3$. $MDD_A$ contains the set of tuples $\{(3, 3, 3, 3), (1, 4, 3, 3), (1, 3, 4, 4), (9, 2, 3, 4)\}$. All the nodes and arcs reached with the unidirectional marking method are represented with plain lines while those not reached with dotted lines.</figcaption></figure>

time complexity is then quadratic in the number of objectives. This method takes advantage of the compression offered by MDDs. Nonetheless, we can notice that a large part of the DFS is shared between the filtering of each objective variable, that is to say there are many repetitions. We then propose an improvement of this method that will execute only two DFSs.

### 3.2.2 Bidirectional Marking

The first step is to identify in $MDD_A$, for all $j \in (1...n)$, all the beginning of tuple $\tau_{(1...j)} = (\tau_1, ..., \tau_j)$ such that $\tau_{(1...j)}$ weakly-dominates $(obj_1^{min}, ..., obj_j^{min})$. This can be done by using a DFS from *root* to find the corresponding paths by following one rule: for the layer *j* it is possible to go only through arcs whose values are less than or equal to $obj_j^{min}$. All the nodes reached with this method are considered as marked from *root*.

**Proposition 10.** If *tt* is reached by the first step of the bidirectional marking method, then the PARETO constraint *C* is not consistent.

*Proof.* If *tt* is reached by the first step of the bidirectional marking method it means that there exists a path corresponding to a tuple that weakly-dominates $\tau^*(O)$. Then *C* is not consistent according to Proposition 7.

<figure><img src="image_234.png" /></figure>
<caption>■ Algorithm 3.1 unidirectional marking algorithm.</caption>

```c
// min is passed by reference
recursiveDFS(MDD_A, D, i, l, u, path[], current, min)
    1  u.isVisited ← true
    2  path[l] ← u
    3  for each arc (u, a, v) do
    4      if (l = i and a ≤ obji^max) or a ≤ obji^min then
    5      if l = i then current ← a
    6      if v = tt then
    7          for each node ∈ path do node.reachTt ← true
    8          if current ≤ min then min ← current
    9      else
   10          if v.reachTt and l ≥ i and current ≤ min then min ← current
   11          if not v.isVisited then
   12      recursiveDFS(MDD_A, D, i, l + 1, v, path, current, min)
```

```c
oneWayMarkingFiltering(O, D, MDD_A)
    12  for each obji ∈ O do
    13      for each node ∈ MDD_A do
    14          node.isVisited ← false
    15          node.reachTt ← false
    16  path[] ← Ø for each index
    17  currentValue ← MAX_INTEGERS
    18  minValue ← MAX_INTEGERS
    19  recursiveDFS(MDD_A, D, i, 1, root, path, currentValue, minValue)
   20  obji^max ← min(obji^max, minValue - 1)
```

The second step is similar to the first one. It consists of identifying in MDD<sub>A</sub>, for all $j \in (1...n)$, all the end of tuple $\tau_{(j...n)} = (\tau_j, ..., \tau_n)$ such that $\tau_{(j...n)}$ weakly-dominates $(obj_i^{min}, ..., obj_n^{min})$. This time the DFS starts from $tt$, takes arcs only in reverse, and follows the same rule as for the first step. All the nodes reached during this step are considered as marked from $tt$.

The inconsistent edge can now be identified: these are the arcs $(u, a, v)$ on the layer $i$, such that $u$ is marked from root and $v$ is marked from $tt$. More formally we have:

**Proposition 11.** Let $\Lambda_i$ be the set of all the arcs $(u, a, v)$ on the layer $i$, such that $u$ is marked from root and $v$ is marked from $tt$. The value $a$ of $obj_i$ satisfies Proposition 8 if and only if $\exists(u, a, v) \in \Lambda_i$

**Proof.** For all the arcs $(u, l, v)$ in $\Lambda_i$ there is at least one path going from root to $u$ that weakly-dominates $(obj_1^{min}, ..., obj_{i-1}^{min})$, and there is also at least one path going from $v$ to $tt$ that weakly-dominates $(obj_{i+1}^{min}, ..., obj_n^{min})$. It means that there is at least one tuple $\tau$ such that $\tau_i = l$ and $\tau \le \tau^*(O, i, l)$. Conversely, if there exists one tuple $\tau$ such that $\tau_i = l$ and $\tau \le \tau^*(O, i, l)$, then there is at least one path going from root to $u$ that weakly-dominates $(obj_1^{min}, ..., obj_{i-1}^{min})$, and there is also at least one path going from $v$ to $tt$ that weakly-dominates $(obj_{i+1}^{min}, ..., obj_n^{min})$. Therefore $(u, l, v)$ belongs to $\Lambda_i$ ◀
<figure><img src="image_235.png" /></figure>

tt is reached from *root*: the constraint is not consistent.

■ **Figure 5** Application of the bidirectional marking in MDD<sub>A</sub>. It represents the set of tuples {((7, 1, 1, 5), (2, 6, 1, 5), (2, 1, 4, 3), (3, 5, 2, 3))}. All the nodes and arcs reached with the bidirectional marking method are represented with plain lines while those not reached with dotted lines.

We immediately have:

► **Proposition 12.** *The bidirectional marking method finds and eliminates for each objective variable $obj_i$ all the value a such that $\tau^*(O, i, a)$ is weakly-dominated, by setting the maximum of $obj_i$ to $\min V - 1$ where $\min V = \min(\{a such that \tau^*(O, i, a) is weakly-dominated\})$.*

Algorithm 3.2 is a possible implementation of this method.

Figure 5 shows two examples based on the same archive with different states for the domains. In Figure 5 (a), $obj_1$ and $obj_4$ are not filtered because for both cases there is no arc between a node marked from *root* and a node marked from *tt* at their corresponding layer. Concerning $obj_2$ there is the arc (*c*, 5, *f*) that satisfies this condition so $obj_2^{max} = \min(8, 5 - 1)$ and the domain becomes {3, 4}. For $obj_3$, the arc (*e*, 4, *h*) is the only one that satisfies this condition so $obj_3^{max} = \min(7, 4 - 1)$ and the domain becomes {2, 3}. In Figure 5 (b), *tt* is reached from *root* during the first step of the bidirectional marking. Therefore there is no possible assignment with current domains.

The MDD is traversed two times with DFS and one time with the search of minimal value for each objective variable, then the time complexity of this method is linear in the size of the MDD.
<caption>Algorithm 3.2 Bidirectional marking algorithm.</caption>

topDownMarking(MDD<sub>A</sub>, D, l, u, markedFromRoot)
1. markedFromRoot[l].add(u)
2. for each arc (u, a, v) do
3. &nbsp;&nbsp;&nbsp;&nbsp;if a ≤ obj<sub>l</sub><sup>min</sup> and v ∉ markedFromRoot[l + 1] then
4. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;topDownMarking(MDD<sub>A</sub>, D, l + 1, v, markedFromRoot)
5. bottomUpMarking(MDD<sub>A</sub>, D, l, v, markedFromTt)
6. &nbsp;&nbsp;&nbsp;&nbsp;markedFromTt[l + 1].add(v)
7. &nbsp;&nbsp;&nbsp;&nbsp;for each arc (u, a, v) do
8. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;if a ≤ obj<sub>l</sub><sup>min</sup> and u ∉ markedFromTt[l] then
9. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;bottomUpMarking(MDD<sub>A</sub>, D, l - 1, u, markedFromTt)
10. twoWaysMarkingFiltering(O, D, MDD<sub>A</sub>)
11. &nbsp;&nbsp;&nbsp;&nbsp;size ← O.size
12. &nbsp;&nbsp;&nbsp;&nbsp;for i ∈ 1... size + 1 do
13. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;markedFromRoot[i] ← ∅
14. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;markedFromTt[i] ← ∅
15. topDownMarking(MDD<sub>A</sub>, D, 1, root, markedFromRoot)
16. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;// If tt si reached, it means that there is no
17. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;// more possible assignment with current domains.
18. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;if tt ∈ markedFromRoot(size + 1) then
19. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;backtrack
20. &nbsp;&nbsp;&nbsp;&nbsp;else
21. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;bottomUpMarking(MDD<sub>A</sub>, D, size, tt, markedFromTt)
22. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;for i ∈ 1... size do
23. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;minValue ← MAX_INTEG<sub>R</sub>
24. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;for each u ∈ markedFromRoot[i] do
25. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;for each arc (u, a, v) do
26. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;minValue ← a
27. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;obj<sub>i</sub><sup>max</sup> ← min(obj<sub>i</sub><sup>max</sup>, minValue - 1)

4 Experiments

The methods presented in this paper have been implemented in Java 17 using Choco-solver version 4.10.10 [9]. All the experiments were run in sequential on a machine with an Intel(R) Xeon(R) W-2175 CPU @ 2.50GHz using Ubuntu 20.04.6 LTS version 5.4.0-146-generic.

In these experiments, three implementations of the PARETO constraint are compared:

List: the PARETO constraint of Choco, using a list for representing the archive. When a new solution is inserted, all the solutions in the list are compared with this solution to determine if they must be removed from the list. Concerning the filtering, for each objective variable *obj*<sub>*i*</sub> the inconsistent values are found by comparing all the solutions in the list with the *dominated point* *DP*<sub>*i*</sub>, defined in [11, 3].

M-U: the PARETO constraint associated with an MDD for representing the archive and using the Unidirectional marking algorithm as filtering algorithm.
24:12 MDD Archive for Boosting the Pareto Constraint

<table><caption>Table 1 Time (s) comparison between the use of lists and MDDs for the PARETO constraint with 10 objectives. The search ends when all solutions are found or if 30000 solutions are found, or if it exceeds 30 minutes (TO).</caption><thead><tr><th rowspan="2">n</th><th rowspan="2">Data</th><th rowspan="2"># Solutions found</th><th rowspan="2">#Solutions in A</th><th colspan="3">Total time</th><th colspan="3">Filtering time</th><th colspan="3">Deletion and insertion time</th></tr><tr><th>List</th><th>M-U</th><th>M-B</th><th>List</th><th>M-U</th><th>M-B</th><th>List</th><th>M-U</th><th>M-B</th></tr></thead><tbody><tr><td rowspan="5">12</td><td>b1</td><td>6770</td><td>1687</td><td>11</td><td>46</td><td><b>9</b></td><td><b>3</b></td><td>23</td><td>8</td><td>6</td><td><b>0.7</b></td><td>0.8</td></tr><tr><td>b2</td><td>7443</td><td>3004</td><td>36</td><td>61</td><td><b>11</b></td><td>7</td><td>26</td><td><b>6</b></td><td>25</td><td><b>0.8</b></td><td>0.9</td></tr><tr><td>b3</td><td>7351</td><td>2651</td><td>29</td><td>55</td><td><b>11</b></td><td><b>6</b></td><td>23</td><td>7</td><td>19</td><td><b>0.6</b></td><td>0.7</td></tr><tr><td>b4</td><td>7754</td><td>3187</td><td>40</td><td>61</td><td><b>13</b></td><td><b>7</b></td><td>32</td><td>8</td><td>29</td><td><b>0.8</b></td><td>0.9</td></tr><tr><td>b5</td><td>8502</td><td>3284</td><td>49</td><td>129</td><td><b>18</b></td><td><b>10</b></td><td>65</td><td>13</td><td>35</td><td><b>0.9</b></td><td><b>0.9</b></td></tr><tr><td rowspan="5">16</td><td>b6</td><td>30000</td><td>4894</td><td>420</td><td>781</td><td><b>103</b></td><td>87</td><td>360</td><td><b>70</b></td><td>309</td><td><b>4</b></td><td><b>4</b></td></tr><tr><td>b7</td><td>30000</td><td>5805</td><td>708</td><td>TO</td><td><b>239</b></td><td>170</td><td>TO</td><td><b>166</b></td><td>487</td><td>TO</td><td><b>5</b></td></tr><tr><td>b8</td><td>30000</td><td>5015</td><td>615</td><td>901</td><td><b>159</b></td><td><b>106</b></td><td>460</td><td>116</td><td>478</td><td><b>3</b></td><td><b>3</b></td></tr><tr><td>b9</td><td>30000</td><td>7763</td><td>1336</td><td>1171</td><td><b>178</b></td><td>236</td><td>467</td><td><b>128</b></td><td>1055</td><td>6</td><td><b>5</b></td></tr><tr><td>b10</td><td>30000</td><td>7144</td><td>1150</td><td>TO</td><td><b>312</b></td><td><b>215</b></td><td>TO</td><td>236</td><td>881</td><td>TO</td><td><b>5</b></td></tr><tr><td rowspan="5">20</td><td>b11</td><td>30000</td><td>2085</td><td>265</td><td>504</td><td><b>115</b></td><td>100</td><td>247</td><td><b>67</b></td><td>123</td><td><b>5</b></td><td><b>5</b></td></tr><tr><td>b12</td><td>30000</td><td>3812</td><td>386</td><td>938</td><td><b>142</b></td><td>131</td><td>508</td><td><b>86</b></td><td>207</td><td><b>6</b></td><td><b>6</b></td></tr><tr><td>b13</td><td>30000</td><td>1198</td><td><b>143</b></td><td>827</td><td>159</td><td><b>56</b></td><td>439</td><td>77</td><td>17</td><td><b>3</b></td><td>4</td></tr><tr><td>b14</td><td>30000</td><td>2647</td><td>246</td><td>1430</td><td><b>205</b></td><td><b>107</b></td><td>709</td><td>108</td><td>56</td><td><b>3</b></td><td><b>3</b></td></tr><tr><td>b15</td><td>30000</td><td>2707</td><td>713</td><td>TO</td><td><b>406</b></td><td>427</td><td>TO</td><td><b>216</b></td><td>108</td><td>TO</td><td><b>5</b></td></tr></tbody></table>

M-B: the PARETO constraint associated with an MDD for representing the archive and using the Bidirectional marking algorithm as filtering algorithm.

The considered problems are the bin packing and the multi-criteria knapsack problems.

## 4.1 Bin Packing Problem

The bin packing problem is a problem where *n* items have to be placed into bins. Each item has *m* types of weight, and each bin has a limit for each type of weight. For each type of weight the objective is to minimize the maximum weight among all the bins, so there are *m* objectives. These objectives encourage an equitable distribution of weights. The datasets used involve items with weights randomly chosen between 1 and 40, and a limit of 120 for each type of weight for each bin. These items have to be distributed between 8 bins. The problem is modeled using the matrix-based symmetry-breaking constraints proposed by Salem and Kieffer [10].

In order to compare the different methods, we focused on the time taken to find at most 30000 solutions. Moreover, we measured the total time taken by the filtering of the Pareto constraint throughout the search (Filtering time), and the total time taken to maintain the domination-free property of the archive throughout the search (Deletion and Insertion time, or D&I time). Table 1 shows the evolution of these times depending on the number of items *n* while Table 2 shows this evolution depending on the number of objectives *m*.

The first thing that comes out of the results of Table 1 is that M-B generally performs better than the list for this problem. These performances seem to depend on the size of the archive: the larger it is, the more the compression of MDDs shows its advantage. When we look at the time spent on the different operations, we can notice that this time saving is mainly done during the D&I part. For example with data b9, which has the largest archive with 7763 solutions, M-B is 7.5 times faster than the list and its D&I part is 200 times faster than the D&I part of the list. Concerning the filtering part M-B is sometimes slower than the list but when this is the case, it is not much slower. Moreover, when M-B is faster on the filtering it can be up to almost 2 times faster as with data b15.

The experiments in Table 2 show another interesting phenomenon about the filtering
<table><caption>Table 2 Time (s) comparison between the use of lists and MDDs for the PARETO constraint with 16 items. The search ends when all solutions are found or if 30000 solutions are found, or if it exceeds 30 minutes (TO).</caption><thead><tr><th rowspan="2">m</th><th rowspan="2">Data</th><th rowspan="2"># Solutions found</th><th rowspan="2">#Solutions in A</th><th colspan="3">Total time</th><th colspan="3">Filtering time</th><th colspan="3">Deletion and insertion time</th></tr><tr><th>List</th><th>M-U</th><th>M-B</th><th>List</th><th>M-U</th><th>M-B</th><th>List</th><th>M-U</th><th>M-B</th></tr></thead><tbody><tr><td rowspan="5">5</td><td>b16</td><td>4619</td><td>121</td><td><b>7</b></td><td>17</td><td>9</td><td><b>0.5</b></td><td>5</td><td>2</td><td><b>0.1</b></td><td>0.2</td><td>0.2</td></tr><tr><td>b17</td><td>5452</td><td>640</td><td><b>14</b></td><td>58</td><td>19</td><td><b>4</b></td><td>28</td><td>9</td><td>1</td><td><b>0.3</b></td><td><b>0.3</b></td></tr><tr><td>b18</td><td>5679</td><td>382</td><td><b>3</b></td><td>12</td><td>5</td><td><b>0.6</b></td><td>6</td><td>2</td><td><b>0.3</b></td><td><b>0.3</b></td><td><b>0.3</b></td></tr><tr><td>b19</td><td>8078</td><td>597</td><td><b>16</b></td><td>52</td><td>20</td><td><b>4</b></td><td>24</td><td>9</td><td>2</td><td><b>0.6</b></td><td>0.7</td></tr><tr><td>b20</td><td>4658</td><td>234</td><td><b>3</b></td><td>6</td><td>4</td><td><b>0.3</b></td><td>2</td><td>1</td><td><b>0.1</b></td><td>0.2</td><td>0.2</td></tr><tr><td rowspan="5">10</td><td>b6</td><td>30000</td><td>4894</td><td>420</td><td>781</td><td><b>103</b></td><td>87</td><td>360</td><td><b>70</b></td><td>309</td><td><b>4</b></td><td><b>4</b></td></tr><tr><td>b7</td><td>30000</td><td>5805</td><td>708</td><td>TO</td><td><b>239</b></td><td>170</td><td>TO</td><td><b>166</b></td><td>487</td><td>TO</td><td><b>5</b></td></tr><tr><td>b8</td><td>30000</td><td>5015</td><td>615</td><td>901</td><td><b>159</b></td><td><b>106</b></td><td>460</td><td>116</td><td>478</td><td><b>3</b></td><td><b>3</b></td></tr><tr><td>b9</td><td>30000</td><td>7763</td><td>1336</td><td>1171</td><td><b>178</b></td><td>236</td><td>467</td><td><b>128</b></td><td>1055</td><td>6</td><td><b>5</b></td></tr><tr><td>b10</td><td>30000</td><td>7144</td><td>1150</td><td>TO</td><td><b>312</b></td><td><b>215</b></td><td>TO</td><td>236</td><td>881</td><td>TO</td><td><b>5</b></td></tr><tr><td rowspan="5">15</td><td>b21</td><td>30000</td><td>2242</td><td>78</td><td>290</td><td><b>44</b></td><td><b>16</b></td><td>179</td><td>30</td><td>53</td><td><b>2</b></td><td><b>2</b></td></tr><tr><td>b22</td><td>30000</td><td>4651</td><td>583</td><td>1005</td><td><b>116</b></td><td>105</td><td>545</td><td><b>91</b></td><td>456</td><td><b>4</b></td><td><b>4</b></td></tr><tr><td>b23</td><td>30000</td><td>3442</td><td>271</td><td>363</td><td><b>68</b></td><td>70</td><td>180</td><td><b>39</b></td><td>178</td><td><b>3</b></td><td><b>4</b></td></tr><tr><td>b24</td><td>30000</td><td>3379</td><td>296</td><td>1180</td><td><b>97</b></td><td><b>54</b></td><td>614</td><td>63</td><td>223</td><td>5</td><td><b>4</b></td></tr><tr><td>b25</td><td>30000</td><td>8421</td><td>1359</td><td>1213</td><td><b>121</b></td><td>305</td><td>532</td><td><b>77</b></td><td>1014</td><td>10</td><td><b>9</b></td></tr><tr><td rowspan="5">20</td><td>b26</td><td>30000</td><td>7145</td><td>911</td><td>1098</td><td><b>77</b></td><td>185</td><td>568</td><td><b>47</b></td><td>706</td><td>6</td><td><b>5</b></td></tr><tr><td>b27</td><td>30000</td><td>6828</td><td>1059</td><td>TO</td><td><b>149</b></td><td>202</td><td>TO</td><td><b>108</b></td><td>828</td><td>TO</td><td><b>7</b></td></tr><tr><td>b28</td><td>30000</td><td>5327</td><td>684</td><td>889</td><td><b>87</b></td><td>140</td><td>529</td><td><b>61</b></td><td>521</td><td><b>5</b></td><td><b>5</b></td></tr><tr><td>b29</td><td>30000</td><td>5791</td><td>819</td><td>TO</td><td><b>247</b></td><td>243</td><td>TO</td><td><b>176</b></td><td>516</td><td>TO</td><td><b>6</b></td></tr><tr><td>b30</td><td>30000</td><td>7142</td><td>864</td><td>855</td><td><b>106</b></td><td>178</td><td>342</td><td><b>63</b></td><td>650</td><td><b>7</b></td><td><b>7</b></td></tr></tbody></table>

part: compared to the list, the more objectives there are, the faster the filtering with M-B. For example with data b26 where there are 20 objectives, the filtering with M-B is almost 4 times faster than the filtering with the list. However, when there are few objectives the list is globally better as shown by the results for $m = 5$.

The results with M-U are not as good as those with M-B, and are even worse than those with the list. Indeed, even if M-U has an advantage on the D&I part compared to the list, the filtering takes too much time which negates totally the advantage.

## 4.2 Multi-Criteria Knapsack Problem

This problem is a variant of the knapsack problem with $n$ items: the goal is not to maximize only one type of profit but $m$ types of profit. Then, each item has $m$ values and there are $m$ objectives to maximize. The data sets used represent items with weights and values randomly chosen between 1 and 40. For each data set, the limit of the knapsack is equal to $(\sum_{i=1}^{n} w_i)/2$ with $w_i$ the weight of the $i$-th item.

We took the same types of measures as for the bin packing problem (i.e. the total time, the filtering time and the D&I time) but we only ran the methods with the list and M-B. Table 3 lists the results obtained.

The results between the methods are similar for this problem, the solving times are more or less equivalent for each instance. However, we can observe a behavior that we have already seen with the bin packing problem: M-B is generally better when the size of the archive is large. The results show also that for this problem the list is generally more efficient for the filtering while M-B is better for the D&I.
<table><caption>Table 3 Time (s) comparison between the use of lists and MDDs for the PARETO constraint with 10 objectives. The search ends when all solutions are found or if 10000 solutions are found.</caption><thead><tr><th rowspan="2">n</th><th rowspan="2">Data</th><th rowspan="2"># Solutions found</th><th rowspan="2">#Solutions in A</th><th colspan="2">Total time</th><th colspan="2">Filtering time</th><th colspan="2">Deletion and insertion time</th></tr><tr><th>List</th><th>M-B</th><th>List</th><th>M-B</th><th>List</th><th>M-B</th></tr></thead><tbody><tr><td rowspan="3">20</td><td>k1</td><td>1809</td><td>640</td><td><strong>182</strong></td><td>186</td><td><strong>1</strong></td><td>4</td><td>0.3</td><td><strong>0.2</strong></td></tr><tr><td>k2</td><td>3667</td><td>1750</td><td><strong>309</strong></td><td>317</td><td><strong>6</strong></td><td>13</td><td>4</td><td><strong>0.9</strong></td></tr><tr><td>k3</td><td>1661</td><td>1661</td><td><strong>134</strong></td><td>135</td><td><strong>0.5</strong></td><td>2</td><td><strong>0.2</strong></td><td><strong>0.2</strong></td></tr><tr><td rowspan="3">21</td><td>k4</td><td>10000</td><td>6714</td><td>786</td><td><strong>672</strong></td><td><strong>52</strong></td><td>97</td><td>175</td><td><strong>6</strong></td></tr><tr><td>k5</td><td>7806</td><td>3265</td><td><strong>865</strong></td><td>868</td><td><strong>27</strong></td><td>59</td><td>15</td><td><strong>2</strong></td></tr><tr><td>k6</td><td>4537</td><td>815</td><td>472</td><td><strong>467</strong></td><td><strong>5</strong></td><td>11</td><td><strong>1</strong></td><td><strong>1</strong></td></tr><tr><td rowspan="3">22</td><td>k7</td><td>3705</td><td>2750</td><td><strong>1011</strong></td><td>1050</td><td><strong>20</strong></td><td>35</td><td>8</td><td><strong>1</strong></td></tr><tr><td>k8</td><td>1895</td><td>601</td><td>450</td><td><strong>446</strong></td><td><strong>2</strong></td><td>6</td><td><strong>0.3</strong></td><td><strong>0.3</strong></td></tr><tr><td>k9</td><td>10000</td><td>5192</td><td>1123</td><td><strong>1073</strong></td><td><strong>59</strong></td><td>92</td><td>73</td><td><strong>6</strong></td></tr><tr><td rowspan="3">23</td><td>k10</td><td>10000</td><td>4561</td><td>1116</td><td><strong>1101</strong></td><td><strong>54</strong></td><td>99</td><td>66</td><td><strong>6</strong></td></tr><tr><td>k11</td><td>4634</td><td>1593</td><td><strong>1360</strong></td><td>1420</td><td><strong>21</strong></td><td>46</td><td>3</td><td><strong>1</strong></td></tr><tr><td>k12</td><td>5390</td><td>2816</td><td>1225</td><td><strong>1193</strong></td><td><strong>22</strong></td><td>40</td><td>11</td><td><strong>2</strong></td></tr><tr><td rowspan="3">24</td><td>k13</td><td>1791</td><td>823</td><td><strong>848</strong></td><td>873</td><td><strong>4</strong></td><td>12</td><td>0.5</td><td><strong>0.3</strong></td></tr><tr><td>k14</td><td>10000</td><td>5781</td><td>1402</td><td><strong>1335</strong></td><td><strong>77</strong></td><td>112</td><td>89</td><td><strong>5</strong></td></tr><tr><td>k15</td><td>10000</td><td>3528</td><td><strong>1545</strong></td><td>1570</td><td><strong>54</strong></td><td>92</td><td>45</td><td><strong>5</strong></td></tr></tbody></table>

# 5 Conclusion

In this paper we presented methods to use an MDD as an archive for the PARETO Constraint. The insertion of a new solution into the MDD and the deletion from the MDD of the solutions dominated by this new solution are made by applying classical operators. We presented two methods for establishing the bound consistency of this constraint: the unidirectional marking method and the bidirectional marking method. The second method is linear in the size of the MDD.

We have shown that, depending on the problem, using an MDD as an archive with the bidirectional marking method can be very effective compared to the classical list representation of the archive. The use of MDDs is particularly well suited to maintaining the dominance-free relation of the archive. It is also interesting for the filtering algorithm. These performances are even more important as the number of objectives increases.

# References

1 David Bergman, André A. Ciré, Willem-Jan van Hoeve, and John N. Hooker. *Decision Diagrams for Optimization*. Artificial Intelligence: Foundations, Theory, and Algorithms. Springer, 2016. doi:10.1007/978-3-319-42849-9.

2 Randal E. Bryant. Graph-based algorithms for boolean function manipulation. *IEEE Trans. Computers*, 35(8):677–691, 1986. doi:10.1109/TC.1986.1676819.

3 Marco Gavanelli. An algorithm for multi-criteria optimization in csp s. In *ECAI’02: Proceedings of the 15th European Conference on Artificial Intelligence*, pages 136–140, January 2002.

4 Renaud Hartert and Pierre Schaus. A support-based algorithm for the bi-objective pareto constraint. *Proceedings of the National Conference on Artificial Intelligence*, 4:2674–2679, June 2014. doi:10.1609/aaai.v28i1.9119.

5 T.Y.K. Kam and Robert K. Brayton. Multi-valued decision diagrams. Technical Report UCB/ERL M90/125, EECS Department, University of California, Berkeley, 1990. URL: http://www2.eecs.berkeley.edu/Pubs/TechRpts/1990/1671.html.

6 Sanaz Mostaghim and Jürgen Teich. *Quad-trees: A Data Structure for Storing Pareto Sets in Multiobjective Evolutionary Algorithms with Elitism*, pages 81–104. Springer London, London, 2005. doi:10.1007/1-84628-137-7_5.
7 Guillaume Perez. *Decision diagrams: constraints and algorithms*. PhD thesis, Université Côte d'Azur, 2017.

8 Guillaume Perez and Jean-Charles Régin. Constructions and in-place operations for mdds based constraints. In Claude-Guy Quimper, editor, *Integration of AI and OR Techniques in Constraint Programming*, pages 279–293, Cham, 2016. Springer International Publishing.

9 Charles Prud'homme and Jean-Guillaume Fages. Choco-solver: A java library for constraint programming. *Journal of Open Source Software*, 7(78):4708, 2022. doi:10.21105/joss.04708.

10 Khadija Salem and Yann Kieffer. An experimental study on symmetry breaking constraints impact for the one dimensional bin-packing problem. In *Conference: 2020 Federated Conference on Computer Science and Information Systems*, pages 317–326, September 2020. doi:10.15439/2020F19.

11 Pierre Schaus and Renaud Hartert. Multi-objective large neighborhood search. In Christian Schulte, editor, *Principles and Practice of Constraint Programming*, pages 611–627, Berlin, Heidelberg, 2013. Springer Berlin Heidelberg.

12 A. Srinivasan, T. Ham, S. Malik, and R. K. Brayton. Algorithms for discrete function manipulation. In *1990 IEEE International Conference on Computer-Aided Design. Digest of Technical Papers*, pages 92–95, 1990. doi:10.1109/ICCAD.1990.129849.

13 Taha Vafaeeenezhad, Reza Tavakkoli-Moghaddam, and Naoufel Cheikhrouhou. Multi-objective mathematical modeling for sustainable supply chain management in the paper industry. *Computers & Industrial Engineering*, 135:1092–1102, 2019. doi:10.1016/j.cie.2019.05.027.

14 Yihan Wang, Zongguo Wen, Jianguo Yao, and Christian Doh Dinga. Multi-objective optimization of synergic energy conservation and co2 emission reduction in china's iron and steel industry under uncertainty. *Renewable and Sustainable Energy Reviews*, 134:110128, 2020. doi:10.1016/j.rser.2020.110128.
<figure><img src="image_236.png" /></figure>
# Learning a Generic Value-Selection Heuristic Inside a Constraint Programming Solver

## Tom Marty

<figure><img src="image_237.png" /></figure>

<figure><img src="image_238.png" /></figure>

<figure><img src="image_239.png" /></figure>

Polytechnique Montréal, Montreal, Canada

Ecole Polytechnique, Palaiseau, France

## Tristan François

<figure><img src="image_240.png" /></figure>

Ecole Polytechnique, Palaiseau, France

## Pierre Tessier

<figure><img src="image_241.png" /></figure>

Ecole Polytechnique, Palaiseau, France

## Louis Gautier

<figure><img src="image_242.png" /></figure>

Ecole Polytechnique, Palaiseau, France

## Louis-Martin Rousseau

<figure><img src="image_243.png" /></figure>

<figure><img src="image_244.png" /></figure>

<figure><img src="image_245.png" /></figure>

Polytechnique Montréal, Montreal, Canada

## Quentin Cappart

<figure><img src="image_246.png" /></figure>

<figure><img src="image_247.png" /></figure>

<figure><img src="image_248.png" /></figure>

Polytechnique Montréal, Montreal, Canada

## Abstract

Constraint programming is known for being an efficient approach to solving combinatorial problems. Important design choices in a solver are the *branching heuristics*, designed to lead the search to the best solutions in a minimum amount of time. However, developing these heuristics is a time-consuming process that requires problem-specific expertise. This observation has motivated many efforts to use machine learning to automatically learn efficient heuristics without expert intervention. Although several generic *variable-selection heuristics* are available in the literature, the options for *value-selection heuristics* are more scarce. We propose to tackle this issue by introducing a generic learning procedure that can be used to obtain a value-selection heuristic inside a constraint programming solver. This has been achieved thanks to the combination of a *deep Q-learning* algorithm, a tailored *reward signal*, and a *heterogeneous graph neural network*. Experiments on *graph coloring*, *maximum independent set*, and *maximum cut* problems show that this framework competes with the well-known impact-based and activity-based search heuristics and can find solutions close to optimality without requiring a large number of backtracks.

**2012 ACM Subject Classification** Computing methodologies → Artificial intelligence; Computing methodologies → Machine learning

**Keywords and phrases** Branching heuristic, Deep reinforcement learning

**Digital Object Identifier** 10.4230/LIPIcs.CP.2023.25

**Supplementary Material**

*Software (Source Code): https://github.com/corail-research/SeaPearl.jl*
archived at **swh:1:dir:fd35f1159af064307cb36323831a0254e9d4060f**

**Funding** The research presented in this paper has been partially funded thanks to a NSERC Discovery Grant held by Quentin Cappart.

**Acknowledgements** We thank CIRRELT for providing essential computing resources. We extend our gratitude to all the people who contributed to the related open-source project: Léo Boisvert, Tom Sander, Ziad El Assal, Malik Attalah, and Marco Novaes.

<figure><img src="image_249.png" /></figure>

© Tom Marty, Tristan François, Pierre Tessier, Louis Gautier, Louis-Martin Rousseau, and Quentin Cappart;
licensed under Creative Commons License CC-BY 4.0

29th International Conference on Principles and Practice of Constraint Programming (CP 2023).

Editor: Roland H. C. Yap; Article No. 25; pp. 25:1–25:19

<figure><img src="image_250.png" /></figure>

LIPICS
Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl Publishing, Germany
# 1 Introduction

Combinatorial optimization has countless industrial applications, such as scheduling, routing, or finance. Unfortunately, most of these problems are NP-hard and, thereby, challenging to solve efficiently. It is why finding good solutions has motivated intense research efforts for many years. Traditional methods for tackling them are somehow based on a *search procedure*: A clever enumeration of the solution space is performed to find a feasible and possibly optimal solution. Among these methods, *constraint programming* (CP) is an exact procedure. It constitutes a popular approach as it offers the possibility to find the optimal solution or good feasible approximations by stopping the search early. An additional asset is its declarative paradigm in modeling, which makes the technology easier for the end-user to grasp. Introducing solver-agnostic modeling languages, such as MiniZinc [35] has greatly facilitated this aspect. Aligned with this goal, the propagation engine inside a CP solver is mostly hidden from the end-user. However, ensuring a generic search procedure is trickier as non-trivial heuristics must be designed to make the solving process efficient for an arbitrary problem. That being said, generic *variable-selection* and *value-selection* heuristics have been successfully designed. Notable examples are *impact-based search* [37] or *activity-based search* [31], but they require computationally intensive initialization and yield poor performance at the beginning of the search. This makes these methods not always appropriate for general use. As a concrete example, the current version of MiniZinc<fnref n="1" /> does not propose generic value-selection heuristics, except *in (out) domain* or *impact-based search*. In practice, heuristics are often designed thanks to problem-specific expert knowledge, which is often out of reach for end-users that do not have a solid background in artificial intelligence.

In another context, *machine learning* (ML) has been recently considered for automating the design of branching heuristics, both in constraint programming [11], mixed-integer programming [16, 23], or SAT solving [41]. Specifically, *reinforcement learning* (RL) [45] or *imitation learning* [22] approaches, often combined with *deep learning* [27], have gained special attention. Although this idea seems appealing, this is not an easy task to achieve in practice as several technical considerations must be taken into account in order to ensure both the efficiency and the genericity of the approach. In constraint programming, we identified three questions to resolve when learning a generic branching heuristic inside a solver. They are as follows:

1. *How to train the machine learning model?* An intuitive way is to leverage an RL agent that would explore the tree search by making branching decisions and rewarding it based on the quality of the solution found on a terminal node. This would typically be done with a *depth-first search* traversal of the tree for getting a certificate of optimality. However, as pointed out by several authors [38, 42], the backtracking operations inside a solver raise difficulties when formalizing the task as a Markov decision process and may require redefining it. Besides, this training scheme intensifies the *credit assignment problem* [32], ubiquitous in reinforcement learning.

2. *How to evaluate the quality of a value selection?* A core component of an RL environment is the *reward function*, which gives a score to each decision performed. The end goal for the agent is to perform a sequence of decisions leading to the best-accumulated sum of rewards. In our case, an intuitive solution would be to reward the agent according to the quality of the solution found. However, this information is only available at terminal nodes, and only a zero reward is provided in branching nodes. This is related to the *sparse reward* problematic, which is known to complicate the training process.

<footnote n="1">¹ https://www.minizinc.org/doc-2.7.0/en</footnote>
3. *How to learn from a CP model?* This question relates to the type of architecture that can obtain a value-selection heuristic from a search node (i.e., a partially solved CP model). A promising direction has been proposed by Gasse et al. [16] for binary mixed-integer programs. They introduced a bipartite graph linking variables and constraints (i.e., the two types of nodes) when a variable is involved in a given constraint. The subsequent architecture is a heterogeneous graph neural network. However, this encoding is not directly applicable in constraint programming, as a CP model generally involves non-binary variables and combinatorial constraints. This has been partially addressed by Chalumeau et al. [12], who introduced a tripartite graph where variables, values, and constraints are specific types of nodes. However, this approach lacks genericity as the method requires retraining when the number of variables changes.

To our knowledge, answering such questions is still an open challenge in the research community. This paper proposes to progress in this direction. It introduces a generic learning procedure that can be used to obtain a value-selection heuristic from a constraint programming model given as input. The approach has been designed to be generic in that it can be used for any CP model given as input. In practice, a specific way to extract *features* from a constraint should be designed for any available constraint, but this has to be done only once per constraint type. In this proof of concept, we limit our experiments to three combinatorial optimization problems, namely *graph coloring*, *maximum independent set*, and *maximum cut*. Specifically, we propose three main contributions, each dedicated to addressing one of the aforementioned difficulties. They are as follows: (1) a learning procedure, based on restarts, for training a reinforcement learning agent directly inside a CP solver, (2) a reward function able to assign non-zero intermediate rewards based on the propagation that has been carried out on the node, and (3) a neural architecture based on a tripartite graph representation and a heterogeneous graph neural network. Experimental results show that combining these three ideas enables the search to find good solutions without requiring many backtracks and competes with the well-known impact-based and activity-based search heuristics.

The paper is structured as follows. The next section presents other approaches related to our contribution. Then, Section 3 introduces succinctly technical background on reinforcement learning and graph neural networks. The core contributions are then presented in Section 4. Finally, Section 5 provides experimental results and closes with a discussion of the results.

## 2 Related Work

Bengio et al. [5] identified three ways to leverage machine learning for combinatorial optimization. First, *end-to-end* learning aims to solve the problem only with a trained ML model. This has been, for instance, considered for the traveling salesman problem [4, 25]. However, such an approach does not guarantee the validity nor optimality of the solution obtained. Second, *learning to configure* is dedicated to providing insights to a solver before its execution. This can be, for instance, the decision to linearize the problem in the context of quadratic programs [7] or to learn when a decomposition is appropriate [26]. This approach is also referred to as parameter tuning [21]. We refer to the initial survey for extended information about these two families of approaches. Third, *learning within a search procedure* uses machine learning within the solver. Our contribution belongs to this last category of methods. Although the idea of combining learning and searching for solving combinatorial optimization problems was already discussed in the nineties [36], it has re-emerged recently with the rise of deep learning. Most combinatorial optimization solvers are based on *branch-and-bound*
25:4 Learning a Generic Value-Selection Heuristic

and *backtracking*. In this context, ML is often used with branching rules to follow. *Imitation learning* [22] has been for instance used to replicate the expensive *strong branching* strategy for mixed-integer programming solvers [16, 23]. One limitation of imitation learning is that the performances are bounded by the performance of the imitated strategy, which remains heuristic and perfectible [43]. This opens the door for RL approaches that have the guarantee to find the best branching strategy eventually [29]. A branching strategy can be split into two challenging decisions, *variable-selection* and *value-selection*. Reinforcement learning approaches have been considered for both of them.

Concerning the learning for selecting the next variable to branch on, Song et al. [42] proposed to combine a *double deep Q-network* algorithm [49] with a graph neural network for carrying out this task. The approach is trained to minimize the expected number of nodes to reach a leaf node using the *first-fail principle*. Although this is a good proxy for pruning a maximum of infeasible solutions for a constraint satisfaction problem, it does not extend naturally to optimization variants, for which one should consider a trade-off between the quality of the solution found and the number of nodes required to reach that solution. Similarly, van Driel et al. [48] leveraged a graph neural network to initialize a variable-selection heuristic for *Chuffed*, a hybrid CP-SAT solver. In an online setting, Doolard and Yorke-Smith [15] also proposed to learn variable ordering heuristics where training time is included in the total solving time. Bandit-based learning approaches were also considered by Xia and Yap to automatically select search heuristics [50].

For the value-selection heuristic, Chu and Stuckey [13] introduced a scoring function which gives a score indicating how good an assignation is, given the current domain. A training phase is the carried out in a supervised manner to learn this scoring function. Cappart et al. [11] proposed to train a model with reinforcement learning outside the CP solver and to integrate the agent, once trained, subsequently in the solver. This has been achieved by reaping the benefits of a dynamic programming formulation of a combinatorial problem. An important limitation of this work is that no information related to the CP solver, such as the propagation achieved on a node, can be used to drive the decision. Chalumeau et al. [12] mitigated this issue by carrying out the learning inside the solver. The model is trained to find the optimal solution and to prove it with the least number of explored search nodes. However, this goal is disconnected from finding the best solution as quickly as possible and is practically hard to achieve, even with a good heuristic. A more realistic goal is to find a good solution quickly without closing the search. This is how the contribution of this paper is positioned.

We want to point out that learning *how to branch* is not the only way to leverage ML inside a combinatorial optimization solver. Related works have also been proposed on learning tight optimization bounds [9] or for accelerating column generation approaches [34]. A recurrent design choice is an architecture based on graph neural networks. We refer to the following survey for more information about combinatorial optimization with graph neural networks [10].

3 Technical Background

This section introduces the required background on reinforcement learning and graph neural
network to grasp the technical aspects of the paper.
## 3.1 Reinforcement Learning

Let $\langle S, A, T, R \rangle$ be a 4-tuple representing a Markov decision process where $S$ is the set of states in the environment, $A$ is the set of actions that the agent can do, $T: S \times A \to S$ is a transition function leading the agent from one state to another, given the action taken, and $R: S \times A \to \mathbb{R}$ is a reward function of taking an action from a specific state. The sequence $[s_1, \dots, s_T]$ from the initial state ($s_1$) of an agent towards a terminal state ($s_T$) is referred to as an *episode*. The returned reward within a partial episode $[s_t, \dots, s_T]$ can be formalized as follows: $G_t = \sum_{i=t}^{T} R(s_i, a_i)$. We intentionally omitted the discounting factor as we do not want to discount the late rewards in our application. The agent is governed by a policy $\pi: S \to A$, which indicates the action that must be taken on a given state. The agent's goal is to find the policy that will lead it to maximize the accumulated reward until a terminal state is reached. The core idea of reinforcement learning is to determine this policy by letting the agent interact with the environment and increasing the probability of taking action if it leads to high subsequent rewards. There are a plethora of reinforcement learning algorithms dedicated to this task, such as *trust region policy optimization* [40] or *soft actor-critic* [18]. We refer to *SpinningUp* website for explanations of the main algorithms [1].

This section presents the core principles of *deep Q-learning* [33], which is the algorithm used in this paper. The idea is to compute an *action-value* function $Q^\pi(s_t, a_t) = G_t$. Intuitively, this function gives the accumulated reward that the agent will obtain when performing the action *a* at state *s* while subsequently following a policy $\pi$. The output of this function for a specific action is referred to as a *Q-value*. Provided that the action-value function can be computed exactly, the optimal policy $\pi^*$ turns to be simply the selection of the action having the highest Q-value on a specific state: $\pi^* = \underset{\pi}{\operatorname{argmax}} Q^\pi(s, a), \forall (s, a) \in (S, A)$. Although the exact computation of Q-values can theoretically be performed, a specific value must be computed for each pair of states and actions, which is not tractable for realistic situations. It is why a tremendous amount of work has been carried out to approximate accurately and efficiently Q-values. Among them, deep Q-learning aims to provide a neural estimator $\hat{Q}(s, a, \theta) \approx Q(s, a)$, where $\theta$ is a tensor of parameters that must be learned during a training phase. This algorithm is commonly enriched with other mechanisms dedicated to speed-up or stabilizing the training process, such as the *double deep Q-network* variant [49] or *prioritized experience replay* [39]. Concerning the neural architecture, we opted for a graph neural network, which is explained in the next section.

## 3.2 Graph Neural Network

Intuitively, the goal of a *graph neural network* (GNN) is to embed information contained in a graph (e.g., the structure of the graph, spatial properties, features of the nodes, etc.) into a $d$-dimensional tensor for each node $u \in V$ of the graph. To do so, information on a node is iteratively refined by aggregating information from neighboring nodes. Each iteration of aggregation is referred to as a *layer* of the GNN and involves parameters that must be learned. Let $h_u^k \in \mathbb{R}^{d \times 1}$ be the tensor representation of node $u$ at layer $k$ of the GNN, $h_u^{k+1} \in \mathbb{R}^{l \times 1}$ be the tensor representation of this node at the next layer ($l$ being the dimension of a node at the layer $k+1$), and $\theta_1 \in \mathbb{R}^{l \times d}$ and $\theta_2 \in \mathbb{R}^{l \times d}$ be two matrices of parameters, respectively. Each GNN layer carries out the following update:

$$h_u^{k+1} = g\left(\theta_1 h_u^k \star \left(\bigoplus_{v \in N(u)} \theta_2 h_v^k\right)\right) \quad \forall u \in V \quad (1)$$

Three operations are involved in this update: (1) $\bigoplus$ is an *aggregation* operator that is dedicated to aggregating the information of neighbors (e.g., *mean-pooling* or *sum-pooling*), (2) $\star$ is a *merging* which enables to combine the information of a node with the ones from the
neighbors (e.g., a concatenation), and (3) *g* is an element-wise non-linear activation function, such as the ones commonly used in fully-connected neural networks (e.g., ReLU [17]). Without loss of generality, the *bias* term is not included in the equation. A concrete implementation of a GNN defines these three functions adequately. The training is conducted in a fully-connected neural network through back-propagation and an optimizer based on stochastic gradient descent such as Adam [24].

4 Learning a Value-Selection Heuristic Inside a Solver

This section presents how a value-selection heuristic can be learned with reinforcement
learning in a CP solver from a model given as input. This is the core contribution of the
paper. Three mechanisms are introduced: (1) a *training procedure based on restarts*, (2) a
*reward function leveraging propagation of domains*, and (3) a *heterogeneous graph neural
network architecture*. They are described individually in the next subsections. They have been
implemented in the recently introduced *SeaPearl.jl* solver [12]. Inspired by the architecture
of MiniCP [30], the main specificity of SeaPearl is to natively integrate support for learning
inside the search procedure. This greatly facilitates the prototyping of new search algorithms
based on learning.

**4.1 Restart-Based Training**

Generally speaking, the performance of a reinforcement learning agent is tightly correlated with the definition of an *episode*. This corresponds to the agent’s interactions with the CP solver’s search procedure and is related to the goal desired for the agent. Two options are discussed in this section, (1) an *episode based on depth-first search*, introduced by Chalumeau et al. [12], and (2) an *episode based on restarts*, which is our first contribution.

Building branching heuristics for solving exact combinatorial optimization problems often concurrently targets two objectives: *finding quickly good solutions* and *proving the optimality of a solution*. The approach of Chalumeau et al. [12] relies heavily on the second objective and aims to minimize the number of visited search nodes before proving optimality (e.g., closing the search). To do so, they defined a training episode as a complete solving process carried out by the depth-first search of a solver and penalized through the reward function the generation of each node. This is illustrated in the left picture of Figure 1. However, this approach suffers from an important difficulty. An episode only terminates when the search is completed, which is often intractable for realistic problems as it requires exploring an exponentially large search tree. This is especially problematic during training, where the heuristic is still mediocre. In addition, using a depth-first search algorithm in a Markov Decision Process (MDP) framework required additional considerations not considered by Chalumeau et al. [12]. For example, using a backtracking algorithm in a regular temporal MDP renders their method prone to the *credit assignment problem* [32]. These considerations have been pointed out by Scavuzzo et al. [38] for mixed-integer programming.

Unlike this approach, we propose to train the model to find high-quality solutions quickly.
To do so, we followed the approach proposed by Cappart et al. [11]: an episode is defined as
a *single dive* in the search tree. No backtrack is allowed; the episode stops when a complete
solution is found or a failure is generated. Once the episode is terminated, a restart from
the root node is performed, and a new episode is generated, hence the name of *restart-based
episode*. This is illustrated in the right picture of Figure 1. One limitation of Cappart et
al. [11] is that episodes are executed outside the CP solver during the training and cannot
use the information updated during propagation for the branching. Inspired by Song et
<figure><img src="image_251.png" /><figcaption>Figure 1 The two training procedures (left: depth-first search [12], right: restart-based - ours).</figcaption></figure>

al. [42] for variable-selection heuristics, we addressed this limitation by executing each episode inside the solver during the training. Formally, this requires defining the dynamics of the environment as a Markov Decision Process (i.e., a tuple $\langle S, A, T, R \rangle$, see Section 3.1). It is defined as follows.

**Set of states** Let $\mathcal{P} = \langle X, D(X), C, O \rangle$ be the expression of a combinatorial optimization problem (COP), defined by its variables ($X$), the related domains ($D$), its constraints ($C$), and an objective function ($O$). Each state $s_t \in S$ is defined as the pair $s_t = (\mathcal{P}_t, x_t)$, where $\mathcal{P}_t$ is a partially solved COP (i.e., some variables may have been assigned), and $x_t \in X$ is a variable selected for branching, at step $t$ of the episode. The initial state $s_1 \in S$ corresponds to the situation after the execution of the fix-point at the root node. A terminal node is reached either if all the variables are assigned ($\forall x \in X: |D_t(x)| = 1$), or if a failure is detected ($\exists x \in X: |D_t(x)| = 0$). The variable selected for branching is obtained through a standard heuristic such as *first-fail*.

**Set of actions** Given a state $s_t = (\mathcal{P}_t, x_t)$, an action $a_t$ corresponds to the selection of a value $v \in D(x_t)$ for branching at step $t$. Finding the most promising value to branch on is the problem addressed in this paper.

**Transition function** Given a state $s_t = (\mathcal{P}_t, x_t)$ and an action $a_t = v$, the transition function executes three successive operations. First, it assigns the value $v$ to the variable $x$ (i.e., $D(x_{t+1}) = v$). Second, it executes the fix-point on $\mathcal{P}_t$ in order to prune the domains (i.e., $\mathcal{P}_{t+1} = \text{fixPoint}(\mathcal{P}_t)$). Third, it selects the next variable to branch on (i.e., $x_{t+1} = \text{nextVariable}(\mathcal{P}_{t+1})$). This results in a new state $s_{t+1} = (\mathcal{P}_{t+1}, x_{t+1})$. Integrating the propagation inside the transition is one important difference with Cappart et al. [11].

**Reward function** The function is defined separately in Section 4.2.

Concerning the training, we opted for a *double deep Q-learning* algorithm [49], known to perform well for discrete action spaces. However, other RL algorithms could also be used. We compared our restart-based training procedure using a simple terminal reward based on the solution's score with the backtracking-based approach of Chalumeau et al. [12] using their reward at each step (penalty of 1 for each explored node). We selected the *maximum independent set problem* for this comparison with instances with 50 nodes. Results are presented using performance profiles [14] in Figure 2. A detailed explanation of the experimental protocol is proposed in Section 5.

We evaluated both methods on two metrics matching the objective for which they were specifically trained. We look at the value of the solution obtained after a single dive (Figure 2a) in the tree search and the number of nodes visited to prove optimality using a depth-first search (Figure 2b). As expected, we observe that the agent trained with the *restart-based learning* strategy allows good results regarding the optimality gap for the first solution found
after a single dive. Remarkably, our method yields a comparable ability to prove optimality compared to Chalumeau et al. [12], whose primary aim was specifically to solve the problem in the minimum number of nodes. This last result has to be mitigated as both RL-based methods lie in the range of the random strategy (shaded blue area).

<figure><img src="image_252.png" /><figcaption>(a) Score of the first solution obtained.</figcaption></figure>

<figure><img src="image_253.png" /><figcaption>(b) Number of node visited until optimality.</figcaption></figure>

Figure 2 Comparison of both training methods on maximum independent set (50 nodes). As a non-learned baseline, we added the performances of an agent performing only random decisions. Training is carried out on randomly generated Barabási-Albert graphs [2]; we selected this type of distribution as the generated graphs are known to mimic human-made and natural organizations. The evaluation is performed on 20 other graphs following the same distribution.

Finally, as shown in Figure 2a, it is important to notice that the optimality gap returned by our method is still non-negligible at the first solution obtained. The complexity of a combinatorial problem lies mainly in closing this gap, which is why backtracking is required. Experiments with backtracking are proposed in Section 5.

## 4.2 Propagation-Based Reward

The definition of our reward must be aligned with our objective of *finding quickly good solutions* for the combinatorial problem. Based on our training procedure, an intuitive function is to reward the agent proportionally to the solution quality found at the end of an episode. In case of an infeasible solution found, a penalty can be given. The main drawback of this rewarding scheme is that this information is only available at terminal nodes, and no reward is provided in branching nodes. This is related to the *sparse reward* problem, which complicates the training process [47]. To address this challenge, one should find a way to give informative intermediate rewards along the solving process. To this end, we propose a new rewarding scheme based on the domain reduction of the objective variable (i.e., the variable that must be minimized or maximized). This reduction happens either thanks to the branching assignment or the application of the fix-point. There are two main components: (1) an *intermediate reward ($r^{\text{mid}}$)* collected at branching nodes, and (2) *terminal reward ($r^{\text{end}}$)* collected only at the end of an episode.

Assuming a minimization problem, the intermediate reward follows two principles: each domain reduction of the largest values of the domain is rewarded, and each domain reduction of the lowest values of the domain is penalized. It is important to note that following these principles does not guarantee the discovery of a good solution at the end of the branch. The rationale is to lead the agent to a situation where the minimum cost can be *eventually* obtained while removing costly solutions. It is formalized in Equations (2) to (4), where $r_t^{\text{mid}}$ is the reward obtained at step $t$, and is illustrated in Figure 3. As shown in Equation (5),
the terminal reward is set to -1 if the leaf node corresponds to an infeasible solution and 0 if it is feasible. Finally, the total reward ($r^{acc}$) accumulated during an episode of T steps is the sum of all intermediate rewards with the final term, as proposed in Equation (6).

$$r_t^{ub} = \#\{v \in D_t(x^{obj}) \mid v \notin D_{t+1}(x^{obj}) \wedge v > \max(D_t(x^{obj}))\} \quad (2)$$

$$r_t^{lb} = \#\{v \in D_t(x^{obj}) \mid v \notin D_{t+1}(x^{obj}) \wedge v < \min(D_t(x^{obj}))\} \quad (3)$$

$$r_t^{mid} = \frac{r_t^{ub} - r_t^{lb}}{|D_1(x^{obj})|} \quad (4)$$

$$r_t^{end} = -1 \text{ if unfeasible solution found (0 otherwise)} \quad (5)$$

$$r^{acc} = \left( \sum_{t=1}^{T-1} r_t^{mid} \right) + r_T^{end} \quad (6)$$

<figure><img src="image_254.png" /><figcaption>Figure 3 Intermediate reward when four values are pruned from the domain.</figcaption></figure>

An experimental analysis of this new reward scheme (propagation-based reward) is carried out for the graph coloring, maximum cut, and maximum independent set problems; we look at the quality of the solution found after a single dive in the search tree. As a baseline, we consider a reward (score reward) that only gives a value at terminal nodes ($r_T^{end}$) without an intermediate reward. Besides, we also consider the solutions returned by a random value-selection heuristic as a baseline. Figure 4 shows the evolution of the quality of the first solution returned (y-axis, averaged on 20 instances of the validation step) with the training time (number of episodes in the x-axis) using for training our restart-based search strategy defined in Section 4.1. Instances are Barabási-Albert randomly generated graphs with 50 nodes. Except for the rewarding scheme, the other parts of the architecture are unchanged. We observe that the propagation-based reward provides a more stable training (Figure 4a) and can converge to a better model or, at least, to an equally good model as the terminal score reward (Figures 4b and 4c).

It should be noted that depending on the problem, the reward signal may remain sparse inside episodes even with our definition; this explains the discrepancy across the three class problems. Indeed, constraint propagation might take several steps to reach the objective variable, meaning that for related intermediate decisions, no value will be pruned from the domain of the objective variable. The graph coloring problem is thus the problem for which taking these intermediate rewards is the most beneficial. Indeed, any previously unused color added will negatively impact the domain of the objective function, yielding an insightful negative reward. Conversely, branching on the maximum independent set problem does not consistently impact the objective function domain through the mechanism of constraint propagation, particularly at the beginning of the search. Our method yields no worse result than the usual reward signal in this setting. This worst-case scenario empirically validates the robustness of this reward.
25:10 Learning a Generic Value-Selection Heuristic

<figure><img src="image_255.png" /><figcaption>■ **Figure 4** Training curve for the two rewarding schemes, each validation step corresponds to performing a single dive in the search tree, the score obtained refers to the quality of the solution found on the leaf node.</figcaption></figure>

4.3 Heterogeneous Graph Neural Network Architecture

An important part of the framework is the neural network architecture that we designed to
perform a prediction of the next value to branch on. A high-level representation is proposed
in Figure 5. Four steps are carried out: (1) a *CP model encoder*, (2) a *graph neural network encoder*, (3) a *neural network decoder*, and (4) an *action-selection policy*. They are detailed
in the next subsections.

<figure><img src="image_256.png" /><figcaption>Figure 5 High-level overview of the neural architecture designed.</figcaption></figure>

Step 1: CP Model Encoder

The core idea is to learn for any CP model given as input, unlike Cappart et al. [11], who require a specific encoding for each combinatorial problem. This has been achieved for mixed-integer programs thanks to a bipartite graph representation [16] and by Chalumeau et al. [12] for CP models thanks to a tripartite graph. This last work does not leverage any feature related to the variables, values, or constraints. We built upon this last approach by adding such features. Specifically, let $\mathcal{P} = \langle X, D(X), C, O \rangle$ be the combinatorial problem we want to encode. The idea consists in building a simple undirected graph $G(V_1, V_2, V_3, f_1, f_2, f_3, E_1, E_2)$ encoding all the information of $\mathcal{P}_t$ from a state $s_t = (\mathcal{P}_t, x_t)$. In this representation, $V_1, V_2$, and $V_3$ are three sets of vertices, $f_1, f_2$, and $f_3$ are three sets of feature vectors, and $E_1$ with $E_2$ are two distinct sets of edges. This yields a graph with three types of nodes decorated with features. The first part of the encoding we propose is as follows: (1) each variable, constraint, and value corresponds to a specific type of node ($V_1 = X, V_2 = C$, and $V_3 = D$), (2) each time a variable $x \in V_1$ is involved in a constraint $c \in V_2$, an edge $(x, c) \in E_1$ is added
between both nodes, (3) each time a value $v \in V_3$ is in the domain of a variable $x \in V_1$, an edge $(v, x) \in E_2$ is added between both nodes. This gives a tripartite graph representation of a CP model generically. This is illustrated in Figure 6. The second part of the encoding is to add features to each node. Intuitively, the features will provide meaningful information and thus improve the quality of the model. The features we considered are proposed below. We note that we can easily extend this encoding by integrating new features.

1. Features attached to variables ($f_1$): the current domain size, the initial domain size, a binary indication if the variable is already assigned, and a binary indication if the variable corresponds to the objective.

2. Features attached to constraints ($f_2$): the constraint type (one-hot encoding), and a binary indication if the constraint propagation has reduced domains.

3. Features attached to values ($f_3$): its numerical value.

<figure><img src="image_257.png" /><figcaption>Figure 6 Representation computed by the CP encoder on a simple example.</figcaption></figure>

## Step 2: Graph Neural Network Encoder

Once the CP model has been encoded as a graph, the next step is to embed this representation as a latent vector of features for each node of the graph (see Section 3.2). We propose to carry out this operation with a graph neural network. Unlike the standard prediction scheme presented in Equation (1), our graph has three types of nodes. For this reason, we opted for a heterogeneous architecture. Concretely, a specific convolution is carried out for each node type. The architecture is detailed in Equations (7) to (9), where $\oplus$ is the sum-pooling or mean-pooling aggregation, operator $(\cdot, \|, \cdot)$ is a concatenation of vectors, $N_x(n)$ is the set of neighbouring nodes of $n$ from $V_1$ (variable), $N_c(n)$ is the set of neighbouring nodes of $n$ from $V_2$ (constraint), $N_v(n)$ is the set of neighbouring nodes of $n$ from $V_3$ (value), $\theta_{1,...,10}^k$ are weight matrices at layer $k$, and $g$ is the leakyReLU activation function [28]. Another difference with the canonical GNN equation is the integration of skip connections ($h_x^0, h_c^0$, and $h_v^0$) allowing to keep at each layer information from the input features. This technique is ubiquitous in deep convolutional networks such as in ResNet [20]. Finally, the initial embedding are initialized as follows: $h_x^0 = \theta_{11}f_1$, $h_c^0 = \theta_{12}f_2$, and $h_v^0 = \theta_{13}f_3$, where $\theta_{11,...,13}$
25:12 Learning a Generic Value-Selection Heuristic

are new weight matrices.

$$h_x^{k+1} = g\left(\theta_1^k h_x^0 \parallel \theta_2^k h_x^k \parallel \left( \bigoplus_{c \in N_c(x)} \theta_3^k h_c^k \right) \parallel \left( \bigoplus_{v \in N_v(x)} \theta_4^k h_v^k \right)\right) \quad \forall x \in V_1 \quad (7)$$

$$h_c^{k+1} = g\left(\theta_5^k h_c^0 \parallel \theta_6^k h_c^k \parallel \left( \bigoplus_{x \in N_x(c)} \theta_7^k h_x^k \right)\right) \quad \forall c \in V_2 \quad (8)$$

$$h_v^{k+1} = g\left(\theta_8^k h_v^0 \parallel \theta_9^k h_v^k \parallel \left( \bigoplus_{x \in N_x(v)} \theta_{10}^k h_x^k \right)\right) \quad \forall v \in V_3 \quad (9)$$

Step 3: Neural Network Decoder

At this step, a d-dimensional tensor is obtained for each graph node. Let $x \in V_1$ be the node representing the current variable selected for branching, and $V_x \subseteq V_3$ the subset of nodes representing the values available for $x$ (i.e., the values that are in the domain of the variable). The goal of the decoder is to predict a Q-value (see Section 3.1) for each $v \in V_x$. The computation is formalized in Equation (10), where $h_x^K$ and $h_v^K$ are the node embedding of variable $x$ and value $v$, respectively, after $K$ iterations of the GNN architecture. The functions $\varphi_x: \mathbb{R}^d \to \mathbb{R}^l, \varphi_v: \mathbb{R}^d \to \mathbb{R}^l, \varphi_q: \mathbb{R}^{2l} \to \mathbb{R}$ are fully-connected neural networks. Such a Q-value must be computed for each value $v \in V_x$. It is internally done thanks to matrix operations, allowing a more efficient computation.

$$\hat{Q}(h_x^K, h_v^K) = \varphi_q \left( \varphi_x(h_x^K) \parallel \varphi_v(h_v^K) \right) \quad \forall v \in V_x \quad (10)$$

Step 4: Action-Selection Policy

Once all the Q-values have been computed for the current variable, the policy is defined by an explorer that can decide to exploit the approximated Q-values by greedily choosing the best action as shown in Equation (11) or decide to select unpromising action associated with a lower Q-value (for example, by selecting a random action with probability $\epsilon$). This behavior derives from the trade-off between exploitation and exploration, which is necessary for early learning when the estimates of Q-values are poor, and when only a few states have been visited. Once trained, the Q-values should represent the branching choice leading to the best decision according to the reward of Equation (6).

$$\pi(v|x) = \operatorname*{argmax}_{v \in V_x} \hat{Q}(h_x^K, h_v^K) \quad (11)$$

Assembling all the pieces, this architecture gives a generic approach to obtaining a data-driven value-selection heuristic inside a CP solver. Concerning the search strategy used for evaluation (which is different from the restart-based one used for training), we propose to embed our predictions inside an iterative limited discrepancy search (ILDS) [19]. This strategy is commonly used when we are confident in the quality of the heuristic. The core idea is to restrict the number of branching choices deviating from the heuristic (i.e., a discrepancy). By doing so, the search will explore a subset of solutions expected to be good while giving a chance to reconsider the value-heuristic selection which is nevertheless prone to errors. This mechanism is enriched with a procedure that iteratively increases the number of discrepancies allowed once a level has been explored.

5 Experiments

The goal of this section is to evaluate the quality of the learned value-selection heuristic and the efficiency of the approach. Three combinatorial optimization problems are considered: graph coloring (COL), maximum independent set (MIS), and maximum cut (MAXCUT).
## 5.1 Experimental Protocol

Three configurations for the distribution of the problems generated are proposed for each problem: *small* (20 to 30 nodes), *medium* (40 to 50 nodes), and *large* (80 to 100 nodes) instances, except for MAXCUT which was already challenging for the *medium* size. Training is carried out on randomly generated Barabási-Albert graph [2] with a density factor varying between 4 and 15 according to the size of the instances. A specific model is trained for each configuration of each combinatorial problem. The training is done using randomly generated instances. Evaluation is then performed on 20 new graphs following the same distributions. The models are trained on an Nvidia Tesla V100 32Go GPU until convergence. It took up to 72 hours of training time for the most difficult cases (*graph coloring* with 80 nodes) and less than 1 hour for the simplest cases (*graph coloring* with 20 nodes). Each operation of the CP solver during training and evaluation is carried out on a CPU Intel Xeon Silver 4116 at 2.10GHz. The approach has been implemented in *Julia* and is integrated into the solver *SeaPearl*. The implementation is available on GitHub with BSD 3-Clause licence<fnref n="2" />.

We compared our approach (*Learned*, ILDS) with two other generic value selection heuristics: *impact-based search* (*Impact*) [37] and *activity-based search* (*Activity*) [31]. The standard *minDomain* heuristic is used for the variable selection. Comparisons with Chalumeau et al. [12] have been provided in Section 4.1. As it has been highlighted that this approach is not suited to find good solutions quickly, it is not included again in the next experiments. Each approach is evaluated with a fixed node budget depending on the parameters of the distribution used to generate the problems. For our approach, the performance obtained after the first dive in the tree search is also monitored (*Learned*, 1st dive). As *Impact* and *Activity* are online learning methods, they perform similarly to a random selection at the beginning of the search. For this reason, the performance obtained after the first dive in the tree search with such methods is omitted. Finally, we also included a comparison with a random selection using DFS with the same node budget (*Random*). Finally, the optimal cost (*OPT*) has been obtained with an exact approach without any restriction on the budget.

## 5.2 Quantitative Results

Table 1 summarizes the main results of our approach. As a general comment, our approach can find solutions of superior quality given a node budget or find the optimal solution by exploring fewer nodes than the baselines. Interestingly, our approach (*Learned*, ILDS) can learn a branching strategy giving high-quality solutions, even without backtracking (1st dive). For instance, a single dive for *maximum cut* with 50 nodes yields almost instantly a solution with an optimality gap of 0.16, whereas a depth-first search with a random selection (*Random*, DFS) required 19 seconds and roughly 53,000 nodes explored to find a solution with the same gap. Within this same budget, (*Learned*, ILDS) significantly improves the solution and achieves an optimality gap of 0.09. It is worth highlighting that (*Learned*, ILDS) took 130 seconds to explore 38,744 nodes and has, thereby, an exploration rate slower than the other methods. This significantly increased execution time is mainly because calling the graph neural network architecture (Section 4.3) at each tree search node is much more computationally expensive than calling a simple heuristic. This difficulty is further discussed in Section 5.3.

Concerning *Activity* and *Impact* heuristics, they yield no improvement on *graph coloring* compared to a random strategy. This can be explained by the fact that this class of problem has many possible combinations of variables and values for branching. This requires a

<footnote n="2">² https://github.com/corail-research/SeaPearl.jl</footnote>
25:14 Learning a Generic Value-Selection Heuristic

<table><caption>Table 1 Results for the three problems given a fixed node budget. The average result (rounded) on the 20 test instances is reported for each configuration. *Gap* indicates the optimality gap, *Node* gives the number of nodes explored before finding the best solution within the budget, and *Time* gives the time (seconds) before finding this solution.</caption><thead><tr><th rowspan="2">Size</th><th rowspan="2">OPT</th><th colspan="4">Learned</th><th colspan="4">Activity-Based</th><th colspan="4">Impact-Based</th><th colspan="4">Random</th><th rowspan="2">Budget</th></tr><tr><th>1<sup>st</sup> dive</th><th colspan="3">ILDS</th><th colspan="3">DFS</th><th colspan="3">DFS</th><th colspan="3">DFS</th><th colspan="3">DFS</th></tr><tr><th></th><th></th><th>Gap</th><th>Gap</th><th>Node</th><th>Time</th><th>Gap</th><th>Node</th><th>Time</th><th>Gap</th><th>Node</th><th>Time</th><th>Gap</th><th>Node</th><th>Time</th><th>Gap</th><th>Node</th><th>Time</th><th></th></tr></thead><tbody><tr><td rowspan="3">COL</td><td>20</td><td>5.05</td><td>0.06</td><td>0</td><td>27</td><td>&lt; 1</td><td>0</td><td>378</td><td>&lt; 1</td><td>0</td><td>374</td><td>&lt; 1</td><td>0</td><td>378</td><td>&lt; 1</td><td>10<sup>3</sup></td><td></td></tr><tr><td>40</td><td>7.90</td><td>0.08</td><td>0</td><td>104</td><td>&lt; 1</td><td>0</td><td>1,664</td><td>&lt; 1</td><td>0</td><td>1732</td><td>&lt; 1</td><td>0</td><td>1735</td><td>&lt; 1</td><td>10<sup>4</sup></td><td></td></tr><tr><td>80</td><td>8.75</td><td>0.06</td><td>0</td><td>120</td><td>1</td><td>0</td><td>7,051</td><td>2</td><td>0</td><td>7,057</td><td>2</td><td>0</td><td>7,211</td><td>2</td><td>10<sup>5</sup></td><td></td></tr><tr><td rowspan="3">MIS</td><td>30</td><td>9.90</td><td>0.08</td><td>0</td><td>88</td><td>&lt; 1</td><td>0</td><td>215</td><td>&lt; 1</td><td>0</td><td>297</td><td>&lt; 1</td><td>0</td><td>293</td><td>&lt; 1</td><td>10<sup>3</sup></td><td></td></tr><tr><td>50</td><td>15.00</td><td>0.09</td><td>0</td><td>539</td><td>1</td><td>0</td><td>5,807</td><td>1</td><td>0</td><td>7,474</td><td>1</td><td>0</td><td>8,942</td><td>1</td><td>10<sup>4</sup></td><td></td></tr><tr><td>100</td><td>21.70</td><td>0.20</td><td>0.02</td><td>28,392</td><td>253</td><td>0.09</td><td>35,536</td><td>7</td><td>0.10</td><td>38,154</td><td>8</td><td>0.10</td><td>41,774</td><td>9</td><td>10<sup>5</sup></td><td></td></tr><tr><td rowspan="2">MAXCUT</td><td>20</td><td>46.70</td><td>0.15</td><td>0.03</td><td>3,714</td><td>5</td><td>0.04</td><td>4,635</td><td>1</td><td>0.03</td><td>5,959</td><td>2</td><td>0.04</td><td>4877</td><td>1</td><td>10<sup>4</sup></td><td></td></tr><tr><td>50</td><td>222.00</td><td>0.16</td><td>0.09</td><td>38,744</td><td>130</td><td>0.17</td><td>44,664</td><td>14</td><td>0.17</td><td>47,970</td><td>17</td><td>0.17</td><td>53,110</td><td>19</td><td>10<sup>5</sup></td><td></td></tr></tbody></table>

significantly larger number of explored nodes to initialize these two heuristics efficiently. For the two other problems, characterized by a binary domain for the values to branch on, Activity and Impact provide significantly better results than the random strategy, which is the expected behavior. In all the tested situations, (Learned, ILDS) provides the best optimality gap within the node budget. Additional results are proposed in Figure 7 using performance profiles [14] for the two hardest situations (100 for maximum independent set, and 50 for maximum cut) given a node budget of 100 or 1000 nodes.

<figure><img src="image_258.png" /><figcaption>(a) Maximum cut with 50 nodes.</figcaption></figure>

<figure><img src="image_259.png" /><figcaption>(b) Maximum independent set with 100 nodes.</figcaption></figure>

<caption>Figure 7 Best solutions found within a restricted node budget on largest instances for the three problems considered. We set a small budget to evaluate the ability of each approach to find quickly a good solution, which is the objective aimed by this work. The performance profile ratio is computed using the optimal solution as a reference. Within the same maximal number of nodes visited (1000), we observe that (Learned, ILDS) dominate all the other methods. Besides, we still perform better than the baselines when restricting ten times the budget for ILDS-Learned.</caption>

5.3 Discussions and Opportunities of Further Research

The previous experiments showcased the promise of this framework to quickly find good solutions towards a generic value-selection heuristic inside a CP solver. There are nonetheless open challenges that must be considered for practical use. Four of them are discussed.
**Challenge 1: Scalability of the Representation**

Our approach faces a double penalty regarding its *scaling* capability: as the problem grows larger, the tripartite representation increases significantly in size, which results in a longer computation time required to make one branching decision. This impacts both training and evaluation. Additionally, the number of nodes (and, therefore, decisions to be made) in the search tree grows exponentially with the problem size, exacerbating the aforementioned phenomenon. Consequently, our approach is penalized twice due to the exponential behavior of combinatorial problems. As a concrete example, *graph coloring* instances with 80 nodes require 72 hours of training on a GPU, while only 1 hour is required for the smallest instances. An interesting research direction to mitigate this difficulty is to build a mechanism to compact the representation, for instance, thanks to network pruning tools [51] or with *transfer learning*. Another idea is to call the model only in a few nodes, in a similar fashion as Cappart et al. [9] did for *decision-diagram-based branch-and-bound* [6]. On a lower level of computation, standard constraint programming solvers perform sequential decisions and are therefore optimized for CPU architecture. Concerning the training, it is carried out on a GPU. In the current implementation, each branching decision requires loading the entire tripartite graph on the Video RAM, which is inefficient. We believe much work could be done to optimize this CPU/GPU architecture, for instance by delegating other operations on the GPUs, such as the propagation of few constraints [8, 46].

**Challenge 2: Tackling Highly Constrained Problems**

The experiments proposed in the paper considered combinatorial problems where the difficulty lay in finding the *best* solution. Still, it was easy to find a *feasible* solution, even of poor quality. We empirically observed that the learning performance largely depends on the abundance of feasible solutions in the search space. This is explained by the definition of the reward, which is based on the propagation occurring on the objective variable (see $r_t^{\text{mid}}$ in Section 4.2). However, when feasible solutions are not easily obtained, such as in highly constrained problems, the reward signal becomes less informative. Addressing such combinatorial problems remains an open challenge. We believe an extension of the reward signal can address this in order to handle other situations.

**Challenge 3: Learning a Combined Variable/Value Heuristic**

Although this work proposes to learn a value-selection heuristic, learning how to branch on variables has already been considered in the literature [42]. An interesting research direction is to adapt this architecture to learn a variable-selection and a value-selection heuristic in a unified way. A possible direction is to consider a model with a double-head decoder, the first for selecting the variable and the second for selecting the value. On the training aspect, two reinforcement learning agents could be trained, with an incentive to cooperate with the information sharing [44].

**Challenge 4: Proving the Optimality of a Solution**

The goal pursued in this paper is to find the best solution as quickly as possible. Another direction is to guide the search to speed-up the optimality proof. It is what has been proposed by Chalumeau et al. [12]. In practice, finding good solutions and proving optimality are complementary aspects inside a constraint programming solver and should be both considered. Possible directions to do so could be to redefine the reward function appropriately or to revise the definition of an episode, as proposed by Scavuzzo et al. with TreeMDPs [38].
# 6 Conclusion

The efficiency of constraint programming solvers is partially due to the branching heuristics used to guide the search. In practice, value-selection heuristics are often designed thanks to problem-specific expert knowledge, often out of reach for non-practitioners. In this paper, we proposed a method based on reinforcement learning for obtaining such a heuristic, thanks to historical data, characterized by problem instances following the same distribution of the one that must be solved. This has been achieved thanks to a restart-based training procedure, a non-sparse reward signal, and a heterogeneous graph neural network architecture. Experiments on three combinatorial optimization problems show that the framework can find better solutions close to optimality in fewer nodes visited than other generic baselines. Several limitations and challenges (e.g., tractability for larger or real-world instances, transfer learning, sparsity of the reward signal) have been identified, and addressing them is part of future work. We also plan to consider other combinatorial problems, such as the ones proposed in XCSP3 competitions [3].

# References

1. Joshua Achiam. Spinning up as a deep RL researcher, October 2018. URL: spinningup.openai.com/en/latest/spinningup/spinningup.html.

2. Réka Albert and Albert-László Barabási. Statistical mechanics of complex networks. *Reviews of modern physics*, 74(1):47, 2002.

3. Gilles Audemard, Christophe Lecoutre, and Emmanuel Lonca. Proceedings of the 2022 XCSP3 competition. *arXiv preprint arXiv:2209.00917*, 2022.

4. Irwan Bello, Hieu Pham, Quoc V. Le, Mohammad Norouzi, and Samy Bengio. Neural combinatorial optimization with reinforcement learning. *arXiv preprint arXiv:1611.09940*, 2016.

5. Yoshua Bengio, Andrea Lodi, and Antoine Prouvost. Machine learning for combinatorial optimization: a methodological tour d'horizon. *European Journal of Operational Research*, 290(2):405-421, 2021.

6. David Bergman, Andre A. Cire, Willem-Jan van Hoeve, and John N. Hooker. Discrete optimization with decision diagrams. *INFORMS Journal on Computing*, 28(1):47-66, 2016.

7. Pierre Bonami, Andrea Lodi, and Giulia Zarpellon. Learning a classification of mixed-integer quadratic programming problems. In *International Conference on the Integration of Constraint Programming, Artificial Intelligence, and Operations Research*, pages 595-604. Springer, 2018.

8. Federico Campeotto, Alessandro Dal Palu, Agostino Dovier, Ferdinando Fioretto, and Enrico Pontelli. Exploring the use of GPUs in constraint solving. In *Practical Aspects of Declarative Languages: 16th International Symposium, PADL 2014, San Diego, CA, USA, January 20-21, 2014. Proceedings 16*, pages 152-167. Springer, 2014.

9. Quentin Cappart, David Bergman, Louis-Martin Rousseau, Isabeau Prémont-Schwarz, and Augustin Parjadis. Improving variable orderings of approximate decision diagrams using reinforcement learning. *INFORMS Journal on Computing*, 2022.

10. Quentin Cappart, Didier Chételat, Elias B. Khalil, Andrea Lodi, Christopher Morris, and Petar Velickovic. Combinatorial optimization and reasoning with graph neural networks. *Journal of Machine Learning Research*, 24(130):1-61, 2023. URL: http://jmlr.org/papers/v24/21-0449.html.

11. Quentin Cappart, Thierry Moisan, Louis-Martin Rousseau, Isabeau Prémont-Schwarz, and Andre A. Cire. Combining reinforcement learning and constraint programming for combinatorial optimization. In *Proceedings of the AAAI Conference on Artificial Intelligence*, volume 35, pages 3677-3687, 2021.

12. Félix Chalumeau, Ilan Coulon, Quentin Cappart, and Louis-Martin Rousseau. Seapearl: A constraint programming solver guided by reinforcement learning. In *International Conference on Machine Learning*.
on *Integration of Constraint Programming, Artificial Intelligence, and Operations Research*, pages 392–409. Springer, 2021.

13. Geoffrey Chu and Peter J. Stuckey. Learning value heuristics for constraint programming. In *International Conference on Integration of Artificial Intelligence and Operations Research Techniques in Constraint Programming for Combinatorial Optimization Problems 2015*, pages 108–123. Springer, 2015.

14. Elizabeth D. Dolan and Jorge J. Moré. Benchmarking optimization software with performance profiles. *Mathematical programming*, 91(2):201–213, 2002.

15. Floris Doolaard and Neil Yorke-Smith. Online learning of variable ordering heuristics for constraint optimisation problems. *Annals of Mathematics and Artificial Intelligence*, pages 1–30, 2022.

16. Maxime Gasse, Didier Chételat, Nicola Ferroni, Laurent Charlin, and Andrea Lodi. Exact combinatorial optimization with graph convolutional neural networks. *Advances in Neural Information Processing Systems*, 32, 2019.

17. Xavier Glorot, Antoine Bordes, and Yoshua Bengio. Deep sparse rectifier neural networks. In *Proceedings of the fourteenth international conference on artificial intelligence and statistics*, pages 315–323. JMLR Workshop and Conference Proceedings, 2011.

18. Tuomas Haarnoja, Aurick Zhou, Pieter Abbeel, and Sergey Levine. Soft actor-critic: Off-policy maximum entropy deep reinforcement learning with a stochastic actor. In *International conference on machine learning*, pages 1861–1870. PMLR, 2018.

19. William D. Harvey and Matthew L. Ginsberg. Limited discrepancy search. In *Proceedings of the 14th international joint conference on Artificial intelligence*, volume 1, pages 607–613, 1995.

20. Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recognition. In *Proceedings of the IEEE conference on computer vision and pattern recognition*, pages 770–778, 2016.

21. Holger H. Hoos. Automated algorithm configuration and parameter tuning. In *Autonomous search*, pages 37–71. Springer, 2011.

22. Ahmed Hussein, Mohamed Medhat Gaber, Eyad Elyan, and Chrisina Jayne. Imitation learning: A survey of learning methods. *ACM Computing Surveys (CSUR)*, 50(2):1–35, 2017.

23. Elias Khalil, Pierre Le Bodic, Le Song, George Nemhauser, and Bistra Dilkina. Learning to branch in mixed integer programming. In *Proceedings of the AAAI Conference on Artificial Intelligence*, volume 30, 2016.

24. Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. *arXiv preprint arXiv:1412.6980*, 2014.

25. Wouter Kool, Herke van Hoof, and Max Welling. Attention, learn to solve routing problems! In *International Conference on Learning Representations*, 2019. URL: https://openreview.net/forum?id=ByxBFsRqYm.

26. Markus Kruber, Marco E Lübbecke, and Axel Parmentier. Learning when to use a decomposition. In *International conference on AI and OR techniques in constraint programming for combinatorial optimization problems*, pages 202–210. Springer, 2017.

27. Yann LeCun, Yoshua Bengio, and Geoffrey Hinton. Deep learning. *Nature*, 521(7553):436–444, 2015.

28. Andrew L. Maas, Awni Y. Hannun, Andrew Y. Ng, et al. Rectifier nonlinearities improve neural network acoustic models. In *Proc. icml*, volume 30, page 3. Atlanta, Georgia, USA, 2013.

29. Nina Mazyavkina, Sergey Sviridov, Sergei Ivanov, and Evgeny Burnaev. Reinforcement learning for combinatorial optimization: A survey. *Computers & Operations Research*, 134:105400, 2021.

30. Laurent Michel, Pierre Schaus, and Pascal van Hentenryck. MiniCP: a lightweight solver for constraint programming. *Mathematical Programming Computation*, 13:133–184, 2021.
25:18 Learning a Generic Value-Selection Heuristic

31 Laurent Michel and Pascal van Hentenryck. Activity-based search for black-box constraint programming solvers. In *International Conference on Integration of Artificial Intelligence and Operations Research Techniques in Constraint Programming*, pages 228–243. Springer, 2012.

32 Marvin Minsky. Steps toward artificial intelligence. *Proceedings of the IRE*, 49(1):8–30, 1961.
doi:10.1109/JRPROC.1961.287775.

33 Volodymyr Mnih, Koray Kavukcuoglu, David Silver, Alex Graves, Ioannis Antonoglou, Daan Wierstra, and Martin Riedmiller. Playing atari with deep reinforcement learning. *arXiv preprint arXiv:1312.5602*, 2013.

34 Mouad Morabit, Guy Desaulniers, and Andrea Lodi. Machine-learning-based column selection for column generation. *Transportation Science*, 55(4):815–831, 2021.

35 Nicholas Nethercote, Peter J. Stuckey, Ralph Becket, Sebastian Brand, Gregory J. Duck, and Guido Tack. Minizinc: Towards a standard CP modelling language. In *International Conference on Principles and Practice of Constraint Programming*, pages 529–543. Springer, 2007.

36 Jean-Yves Potvin, Danny Dubé, and Christian Robillard. A hybrid approach to vehicle routing using neural networks and genetic algorithms. *Applied Intelligence*, 6(3):241–252, 1996.

37 Philippe Refalo. Impact-based search strategies for constraint programming. In *International Conference on Principles and Practice of Constraint Programming*, pages 557–571. Springer, 2004.

38 Lara Scavuzzo, Feng Yang Chen, Didier Chételat, Maxime Gasse, Andrea Lodi, Neil Yorke-Smith, and Karen Aardal. Learning to branch with tree MDPs. *arXiv preprint arXiv:2205.11107*, 2022.

39 Tom Schaul, John Quan, Ioannis Antonoglou, and David Silver. Prioritized experience replay. *arXiv preprint arXiv:1511.05952*, 2015.

40 John Schulman, Sergey Levine, Pieter Abbeel, Michael Jordan, and Philipp Moritz. Trust region policy optimization. In *International conference on machine learning*, pages 1889–1897, 2015.

41 Daniel Selsam and Nikolaj Bjørner. Guiding high-performance SAT solvers with unsat-core predictions. In *International Conference on Theory and Applications of Satisfiability Testing*, pages 336–353. Springer, 2019.

42 Wen Song, Zhiguang Cao, Jie Zhang, Chi Xu, and Andrew Lim. Learning variable ordering heuristics for solving constraint satisfaction problems. *Engineering Applications of Artificial Intelligence*, 109:104603, 2022.

43 Haoran Sun, Wenbo Chen, Hui Li, and Le Song. Improving learning to branch via reinforcement learning. In *Learning Meets Combinatorial Algorithms at NeurIPS2020*, 2020.

44 Peter Sunehag, Guy Lever, Audrunas Gruslys, Wojciech Marian Czarnecki, Vinicius Zambaldi, Max Jaderberg, Marc Lanctot, Nicolas Sonnerat, Joel Z. Leibo, Karl Tuyls, and Thore Graepel. Value-decomposition networks for cooperative multi-agent learning, 2017. arXiv:1706.05296.

45 Richard S. Sutton and Andrew G. Barto. *Reinforcement learning: An introduction*. MIT press, 2018.

46 Fabio Tardivo and Agostino Dovier. Constraints propagation on GPU: A case study for alldifferent. In *Proceedings of the 37th Italian Conference on Computational Logic*, 2022.

47 Alexander Trott, Stephan Zheng, Caiming Xiong, and Richard Socher. Keeping your distance: Solving sparse reward tasks using self-balancing shaped rewards. *Advances in Neural Information Processing Systems*, 32, 2019.

48 Ronald van Driel, Emir Demirović, and Neil Yorke-Smith. Learning variable activity initialisation for lazy clause generation solvers. In *Integration of Constraint Programming, Artificial Intelligence, and Operations Research: 18th International Conference, CPAIOR 2021, Vienna, Austria, July 5–8, 2021, Proceedings 18*, pages 62–71. Springer, 2021.

49 Hado van Hasselt, Arthur Guez, and David Silver. Deep reinforcement learning with double Q-learning. In *Proceedings of the AAAI conference on artificial intelligence*, volume 30, 2016.
50 Wei Xia and Roland Yap. Learning robust search strategies using a bandit-based approach. In *Proceedings of the AAAI Conference on Artificial Intelligence*, volume 32, 2018.

51 Xin Yu, Thiago Serra, Srikumar Ramalingam, and Shandian Zhe. The combinatorial brain surgeon: Pruning weights that cancel one another in neural networks. In *International Conference on Machine Learning*, pages 25668–25683. PMLR, 2022.
<figure><img src="image_260.png" /></figure>
# Proof Logging for Smart Extensional Constraints

Matthew J. McIlree

<figure><img src="image_261.png" /></figure>

University of Glasgow, UK

Ciaran McCreesh

<figure><img src="image_262.png" /></figure>

University of Glasgow, UK

## Abstract

Proof logging provides an auditable way of guaranteeing that a solver has produced a correct answer using sound reasoning. This is standard practice for Boolean satisfiability solving, but for constraint programming, a challenge is that every propagator must be able to justify all inferences it performs. Here we demonstrate how to support proof logging for a wide range of previously uncertified global constraints. We do this by showing how to justify every inference that could be performed by the propagation algorithms for two families of generalised extensional constraint: “Smart Table” and “Regular Language Membership”.

**2012 ACM Subject Classification** Theory of computation → Logic and verification; Theory of computation → Discrete optimization

**Keywords and phrases** Constraint programming, proof logging, extensional constraints

**Digital Object Identifier** 10.4230/LIPIcs.CP.2023.26

**Supplementary Material** Software: https://doi.org/10.5281/zenodo.8132457

**Funding** Part of this work was done while the authors were participating in a programme at the Simons Institute for the Theory of Computing.

Ciaran McCreesh: Supported by a Royal Academy of Engineering research fellowship.

**Acknowledgements** The authors would like to thank Jakob Nordström for several helpful discussions regarding pseudo-Boolean encodings and unit propagation.

## 1 Introduction

A proof log for a problem-solving algorithm provides a verifiable certificate that the result is correct, and also an auditable record of the steps taken to obtain that result. In the field of Boolean satisfiability (SAT), proof-logging has become an expected capability of modern solvers, with a series of proof formats including DRAT [12, 13], LRAT [7], and FRAT [2] widely accepted for independent verification. A similar standard practice has not yet been adopted for Constraint Programming (CP) due to the difficulties of creating easily verifiable proofs for the more expressive formulations and reasoning present in this paradigm. However, recent work by Gocht et al. [10] has shown how the VeriPB proof system [3, 9] can be used to certify the reasoning carried out for several important expressive global constraints, offering a strong candidate for a complete, general CP proof logging method. In this setting, every propagator for a global constraint must be able to do two things: describe its semantics in a lower-level *pseudo-Boolean* format, and justify any reasoning it carries out using either *cutting planes* [6] or *reverse unit propagation* (RUP) [11] reasoning. Describing a constraint's semantics is a well-understood problem, but justifying propagation is not. Gocht et al. demonstrated proof logging for a range of global constraints and propagation strengths, including bounds-consistent integer linear inequalities and domain-consistent table constraints, but it remains an open question whether every global constraint propagation algorithm can be justified in this manner.

<figure><img src="image_263.png" /></figure>

© Matthew J. McIlree and Ciaran McCreesh; licensed under Creative Commons License CC-BY 4.0

29th International Conference on Principles and Practice of Constraint Programming (CP 2023).

Editor: Roland H. C. Yap; Article No. 26; pp. 26:1–26:17

<figure><img src="image_264.png" /></figure>

LIPICS

Leibniz International Proceedings in Informatics

Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl Publishing, Germany
This paper shows that the reasoning carried out by domain-consistent propagators for the SmartTable [16] and Regular [19] constraints can similarly be proof logged efficiently inside the VeriPB proof system. As well as being useful in their own right, these two constraints provide the necessary building blocks for implementing many others, since they allow for efficient strong propagation for extensional constraints that cannot be expressed efficiently as a conventional table. For example, SmartTable can be used to implement the Lex, AtMostOne and NotAllEqual constraints [16], whilst Regular can be used to implement Stretch, Geost and DiffN [15, 19]. Together, these two new constraints bring us a lot closer to fully auditable combinatorial solving, particularly in areas such as workforce scheduling where legal restrictions apply, and where algorithmic decisions can have a large effect upon people's livelihoods.

## 2 An Overview of How Proof Logging Works

The VeriPB proof system is built upon pseudo-Boolean (PB) problems. A PB model is a restricted constraint satisfaction problem, where all variables $x_i$ have domain $\{0, 1\}$ and all constraints are integer linear inequalities of the format $\sum_i c_i \ell_i \ge A$ where $c_i, A \in \mathbb{Z}$, and each $\ell_i$ is a pseudo-Boolean *literal*: either a variable $x$ or its negation $\bar{x} = 1 - x$. For clarity and by convention we will always refer to CP variables with upper-case symbols and PB variables with lower-case symbols. Additionally, to save space and make the meaning clearer we will make use of *reified* PB constraint shorthands:

- $(\ell_1 \land \dots \land \ell_m) \implies \sum_i c_i x_i \ge A$ means $\sum_{i=1}^n c_i x_i + J \bar{\ell}_1 + \dots + J \bar{\ell}_m \ge A$ where $J$ is chosen to be suitably large, e.g. $J = A - \sum_{i=1}^n \min(c_i, 0)$.

- $\sum_i c_i x_i \ge A \implies (\ell_1 \land \dots \land \ell_m)$ means $\{\sum_{i=1}^n -c_i x_i + K \bar{\ell}_t \ge -A + 1 : 1 \le t \le m\}$ where $K$ is chosen to be suitably large, e.g. $K = -A + 1 - \sum_{i=1}^n \min(c_i, 0)$.

- $(\ell_1 \land \dots \land \ell_m) \iff \sum_i c_i x_i \ge A$ means both of the above together.

PB constraints can also be viewed as a superset of conjunctive-normal form (CNF) for SAT problems, since a clause such as $x_1 \lor \bar{x}_2 \lor x_3$ is always equivalent to a PB constraint such as $x_1 + \bar{x}_2 + x_3 \ge 1$. This means that it is straightforward to generalise known SAT encodings of CP constraints to PB. The concept of *unit propagation* from SAT can also be generalised to a PB setting. We say a PB constraint propagates a literal $\ell$ if it cannot possibly be satisfied unless $\ell = 1$. For the “at least 1” constraints that can be viewed as clauses, this is the same as SAT unit propagation, but in general it is enforcing integer bounds consistency [5] on the PB inequalities. The (generalised) notion of unit propagation for PB constraints is therefore repeatedly setting any literals that propagate until either a contradiction or a fixed point is reached.

An auditable constraint solver following the proof-logging methodology of Gocht et al. [10] must be able to produce a pseudo-Boolean representation of the problem being solved, and then justify its backtracking and inference steps in a proof log by outputting further PB constraints derivable by one of the proof rules allowed by the VeriPB proof checker. The full proof system associated with VeriPB is described in detail by Bogaerts et al. [3], but for the purposes of this paper the only rule that is needed is the *reverse unit propagation* (RUP) rule, which says that a PB constraint *D* is derivable from a set of PB constraints *F* if applying the pseudo-Boolean generalisation of unit propagation (repeated integer bounds consistency) to *F* ∪ *D* results in a contradiction. Intuitively, it is a constraint that “obviously” follows once it is pointed out to the verifier.

This paper will omit discussion of how the variables and constraints of CP problem can be represented as PB constraints and variables in general, focussing specifically on how Regular and SmartTable can be represented. We will assume that we have access to a
collection of properly defined PB variables $x_{bneg}, x_{b0}, x_{b1}, x_{b2}...$ representing bits in the two's-complement encoding of a CP variable $X$ (with suitable constraints to enforce the domain), but also that we can freely use flags of the form $x_{=v}$, which are defined to equal 1 if and only if $X = v$. Explanation of how these encodings can be achieved is given by Gocht et al. [10].

Similarly, we will omit detailed discussion of how complete proofs of solutions, unsatisfiability, optimality, or enumeration are achieved. We simply note that each proof is essentially a description of the solver's backtracking search tree. Any time the solver backtracks, a RUP constraint that encodes the negated conjunction of the currently guessed assignments is emitted. A new propagation algorithm can fit into this existing framework by ensuring that the encoding of any specific filtering and failure inferences it makes are visible to the verifier on reverse unit propagation of the backtrack constraint. So for example, if a solver has guessed $A = 2, B = 1$; and a propagator is able to infer $C \neq 1$ then the proof logging methodology requires that the PB constraint

$$ \bar{a}_{=2} + \bar{b}_{=1} + \bar{c}_{=1} \geq 1 \qquad (1) $$

is somehow derived, which can be interpreted as

$$ \bar{a}_{=2} \lor \bar{b}_{=1} \lor \bar{c}_{=1} \quad \text{i.e.} \quad A = 2 \land B = 1 \implies C \neq 1. \qquad (2) $$

This is in many ways similar to how lazy clause generating solvers work [18], except that the new constraints introduced must be justified, rather than asserted.

# 3 Proof Logging for Smart Table Constraints

A smart table constraint generalises the idea of a table constraint to allow wildcards and comparisons, allowing for compact representations for a much larger set of relations, whilst still retaining an efficient domain-consistent propagator [16, 20, 21]. There are several ways to define SmartTable (also called HybridTable), depending on which restriction types are allowed within tuples [4]. In this work we restrict our focus to unary comparisons, binary comparisons, and set membership, and define SmartTable as follows.

## ► Definition 1. Smart Table Constraint

Let $X$ be a sequence of finite-domain variables $\langle X_1, \dots, X_n \rangle$. A smart entry constraint is a unary or binary constraint in one of three possible forms:

1. $\langle \text{var} \rangle \langle \text{op} \rangle a$

2. $\langle \text{var} \rangle \in S$ or $\langle \text{var} \rangle \notin S$

3. $\langle \text{var} \rangle \langle \text{op} \rangle \langle \text{var} \rangle$

where a is an integer constant, $S$ is a set of integer constants, $\langle \text{op} \rangle$ denotes an operator in the set $\{ \langle \text{<}, \text{\textless;}, \text{=}, \text{\textless;}, \text{\textgreater} \}$, and $\langle \text{var} \rangle$ denotes a variable. A smart tuple is a set of these smart entry constraints, and a smart table is a set of smart tuples. When a variable does not appear in a given tuple, it is implicitly unrestricted (equivalent to a wildcard entry).

For a given smart table $T$, where the scope of every smart entry constraint is a subset of $X$, a smart table constraint $\text{SmartTable}(X, T)$ requires that there is at least one smart tuple in $T$ where every smart entry constraint holds. It can hence be thought of as a special case of a disjunction of conjunctions of simple constraints.
## 3.1 PB Encodings for Smart Table Constraints

Recall that in order to support proof logging for a constraint, we must be able to describe its semantics in a lower-level PB form. Let $T$ be a smart table with smart tuples $\sigma_1, \dots, \sigma_k$. We can encode the **SmartTable** constraint by making use of some auxiliary PB variables. For each smart tuple $\sigma_i$ in the table we will have a selector PB variable $s_i$ that controls whether the tuple is active (i.e. not yet infeasible). Then for each smart entry $C_j \in \sigma$ we will have a PB variable $e_{ij}$ that is set to 1 if and only the constraint $C_j$ is satisfied, which we denote by $\mathbf{satisfied}(C_j)$.

Since smart entries are themselves simple binary and unary constraints, with all of them having PB encodings that already well understood, the process of encoding this satisfied property is straightforward. Inequalities on one or two variables can be encoded by imposing the inequality on the evaluation of the bit encodings. For example, for two CP variable $A$ and $B$ both with domain $\{1, \dots, 7\}$, and recalling that PB variables $a_{bi}$ and $b_{bi}$ denote the $i_{th}$ bits in the encoding of the variables $A$ and $B$, we would have

$$ \mathbf{satisfied}(A < B) := -a_{b0} - 2a_{b1} - 4a_{b2} + b_{b0} + 2b_{b1} + 4b_{b2} \geq 1. \quad (3) $$

Other inequality constraints can be handled similarly, adding or subtracting constants from the right-hand side as necessary. For $\in$ and $\notin$ it is simply a case of imposing “at least $n$” constraints on the flags representing the disallowed values, e.g.

$$ \mathbf{satisfied}(A \in \{1, 3, 5\}) := \bar{a}_{=2} + \bar{a}_{=4} + \bar{a}_{=6} + \bar{a}_{=7} \geq 4, \quad (4) $$

$$ \mathbf{satisfied}(B \notin \{1, 3, 5\}) := \bar{b}_{=1} + \bar{b}_{=3} + \bar{b}_{=5} \geq 3. \quad (5) $$

Finally, for $=$ and $\neq$ we can make direct use of the $x_{=v}$ flag if the right-hand side is a value. If the right-hand side is instead another variable we make use of further intermediate flags to express the relation in terms of strict or non-strict inequalities, so:

$$ \mathbf{satisfied}(A = B) := f_{A \geq B} + f_{A \leq B} \geq 2, \quad (6) $$

$$ \mathbf{satisfied}(A \neq B) := f_{A > B} + f_{A < B} \geq 1, \quad (7) $$

where the flags are enforced with two-way implication, e.g.

$$ f_{A<B} \iff \mathbf{satisfied}(A < B). \quad (8) $$

Obviously, the negation of these encodings can be expressed by encoding the negation of the constraint, e.g.

$$ \neg \mathbf{satisfied}(A < B) := \mathbf{satisfied}(A \geq B), \quad (9) $$

$$ \neg \mathbf{satisfied}(A \neq B) := \mathbf{satisfied}(A = B). \quad (10) $$

With the procedures for each of these smart entry constraint encodings in place, and recalling that we know how to reify arbitrary PB constraints on a literal, the PB model for smart table can then be created according to the following specification.

### ► Encoding 1. Smart Table PB Encoding

1. for all $\sigma_i \in T$

2. for all $C \in \sigma_i$

3. $e_C \iff \mathbf{satisfied}(C)$

4. $s_i \iff \sum_{C \in \sigma_i} e_C \geq |\sigma_i|$

5. $s_1 + s_2 + \dots + s_{|T|} \geq 1$
## 3.2 Justifying Smart Table Propagation

The goal of the propagator for the smart table constraint, as with many global constraint propagators, is to achieve domain consistency. As a running example, suppose we had a smart table constraint on three variables, A, B, C, all with domains {1, 2, 3}, defined by the following table T with two tuples:

$$ T := \{\sigma_1, \sigma_2\}; \qquad \sigma_1 = \{(A < B), (A \in \{1, 2\}), (C = 3)\}, \qquad (11) $$

$$ \sigma_2 = \{(A = B), (A \neq 1), (B \geq C)\}. \qquad (12) $$

By inspection, we can observe that the literal $B = 1$ does not have any support on this constraint, as it would contradict $A < B$ on the first tuple and $A = B, A \neq 1$ together on the second. Therefore, a propagation algorithm that achieves domain consistency should immediately remove the value 1 from the domain of B.

The propagation algorithm described by Mairy et al. [16] uses a “Simple Tabular Reduction” strategy to eliminate such unsupported values. Essentially it initialises a set sl (for “supportless”) with every variable value/pair that is possible given the current domains of variables, and then iterates through the tuples, removing any values that they in turn support. The literals that are left in sl are not supported by any tuple and hence the ones that should be eliminated. A key observation is that each smart tuple $\sigma$, as a conjunction of simple constraints, can be thought of as a small CSP $P_{\sigma}$ in its own right, and so the supported values are just all the solutions for this problem. It is shown that as long as the constraint graph of $P_{\sigma}$ is acyclic it can be effectively filtered as a collection of tree CSPs via a two pass filtering process. Due to the restricted structure of the smart tuples, the result of this is a collection of copies of the domains where the only values present are those that appear in some solution (global consistency for $P_{\sigma}$), and so this can be used to remove values from sl in polynomial time.

For proof logging, the main concern is then how to justify the elimination of these unsupported values. In particular, for a sequence of guesses $\mathcal{G}$ and newly derived unsupported assignment $X \neq v$, we would like to log the PB constraint

$$ \wedge \mathcal{G} \implies \bar{x}_{v} \geq 1. \qquad (13) $$

In some situations this might follow directly by reverse unit propagation, as is always the case for proof-logging the classical table constraint as shown by Gocht et al. [10]. For smart tables, however, we can show that more work is required in general. Going back to the previous example and following Encoding 1, we note that if we try to derive $\bar{b}_{=1} \geq 1$ by reverse unit propagation, we do not reach a contradiction after propagating its negation $b_{=1}$. This is despite $B = 1$ lacking support in the table. To see why this is the case, note that the only way unit propagation of a non-auxiliary literal such as $b_{=1}$ could falsify the whole model would be for it to eventually result in all the selectors $s_i$ being falsified. In turn, the only way that could happen is for there to be at least one $e_{iC}$ falsified for every $i \in \{1, \dots, |T|\}$. We can observe that this does not happen on propagation of $b_{=1}$. The smart entry constraints that are responsible for eliminating ($B = 1$)'s support are $A < B$ in the first tuple and both $A = B, A \neq 1$ in the second. From Encoding 1, the first of these is present in the PB model in the form of two constraints:

$$ e_{A<B} \implies b_{b0} + 2b_{b1} - a_{b0} - 2a_{b1} \geq 1; \qquad (14) $$

$$ \bar{e}_{A<B} \implies -b_{b0} - 2b_{b1} + a_{b0} + 2a_{b1} \geq 1. \qquad (15) $$
26:6 Proof Logging for Smart Extensional Constraints

Assuming the bit encodings are correctly defined, propagating $b_{=1}$ should propagate $b_{b0}$ and $\bar{b}_{b1}$, meaning (14) is reduced to

$$e_{A<B} \implies 1 - a_{b0} - 2a_{b1} \geq 1. \qquad (16)$$

Now it might seem that $\bar{e}_{A<B}$ would then propagate, as we know $a_{b0} + 2a_{b1}$ is at least 1 due to the domain constraints of the variable $A$ and so the right-hand side of the implication in (16) must be false. However, because the falsity does not follow from a contradiction intrinsic to the PB itself – in isolation it can be satisfied by having both $b_{b0}$ and $b_{b1}$ equal to 0 – unit propagation alone is not strong enough to determine that $\bar{e}_{A<B}$ must follow. A similar problem holds for the other tuple, where the contradiction arises due to two constraints being incompatible rather than just one being falsified.

Fortunately, we can ensure that the desired constraints do follow by RUP by first explicitly deriving some intermediate proof steps. For a sequence of guesses $\mathcal{G}$ and newly derived unsupported assignment $X \neq v$, instead of simply logging the PB constraint (13) as above, we first log

$$\wedge \mathcal{G} \implies \bar{x}_{=v} + \bar{s}_k \geq 1 \qquad (17)$$

for each $k \in \{1, \dots, |T|\}$, i.e. first show by RUP that no tuple selector can be set to true without contradiction, and then derive from these the desired constraint (13). So in our running example, we would log

$$\bar{b}_{=1} + \bar{s}_1 \geq 1; \qquad \bar{b}_{=1} + \bar{s}_2 \geq 1; \qquad \bar{b}_{=1} \geq 1. \qquad (18)$$

This is now sufficient for most cases, however, it is still possible to construct examples where even these constraints in the form of (17) do not follow by RUP. Suppose we have four variables $W, X, Y, Z$, all with domain $\{-2, -1, 0\}$, and then a single smart tuple of the form

$$(W < X), (X < Y), (Y < Z). \qquad (19)$$

Clearly this is unsatisfiable, and so in particular $W = -2$ should be unsupported by the tuple. But the constraint

$$\bar{w}_{=-2} + \bar{s}_1 \geq 1 \qquad (20)$$

does not unit propagate to contradiction. This can be seen by considering the PB encoding of the first smart entry constraint,

$$e_{W<X} \iff -2x_{\text{bneg}} + x_{b0} + 2w_{\text{bneg}} - w_{b0} \geq 1. \qquad (21)$$

The negation of (20) will lead to $e_{W<X}, w_{\text{bneg}}$, and $\bar{w}_{b0}$ being propagated, which would mean the two PB constraints encoding both implication directions in (21) are reduced to

$$-2x_{\text{bneg}} + x_{b0} \geq -1; \quad \text{and} \quad 2x_{\text{bneg}} - x_{b0} \geq -1. \qquad (22)$$

But then no further propagation results from these constraints, which can be seen either by calculating the slack value (see Elffers et al. [8]), or simply by observing that the value sets $\{0, 0\}, \{1, 1\}, \{0, 1\}$ would each satisfy both constraints if assigned to $\{x_{\text{bneg}}, x_{b0}\}$ and so both 0 and 1 are both supported possibilities for each variable. None of the other constraints contain bit variables for $W$ and so it is clear that no propagation will result from these either, and hence no contradiction will be reached.
This weak propagation is not surprising given the context given by Gocht et al. [10], where integer linear inequality constraints need to log all the inferences made when enforcing domain consistency and cannot rely on the disallowed values being implicit. Fortunately for us, these same inferences are being made in the solver as part of the smart table propagator anyway, at the tree filtering stage, and so the desired constraints in the form of (17) and then (13) can still be logged providing we have first logged the inferences for binary constraints explicitly. The only difference in proof logging these and the proof logging of equality and inequality constraints described by Gocht et al. [10] is that each inference will be conditional on a tuple selector. For example, if a variable B has domain with maximum value l and the filtering process for A < B updates the domain copy for A so that A < l, then we would log

$$ \wedge \mathcal{G} \implies \bar{s}_1 + a_{<l} \geq 1; \qquad (23) $$

where the flag $a_{<l}$ can be introduced via reification with corresponding constraints on the bit variables.

We can now be confident we have a complete proof logging procedure. We can show by an inductive argument on the binary representation that the PB constraints that encode the filtering inferences made by simple binary constraints such as (23) will always be RUP, providing they are logged in the same order as they are made by the solver. Then, any unsupported value must be unsupported in every tuple, and hence removed in the filtering process for some smart entry constraint. Given the negation of this removal, i.e. a variable being equal to some unsupported value, if it was removed due to a binary constraint then we have already logged a corresponding contradicting restriction and hence will unit propagate to contradiction. Otherwise, if it was removed due to a unary constraint it will propagate to contradiction immediately as the assignment will contradict a single PB constraint.

## 3.3 A Complete Worked Example for Smart Table Propagation

To demonstrate how this procedure works on a single round of the domain-consistent smart table propagator we will now show all the steps required for the running example. Firstly, we will have the PB constraints encoding that the domains of A, B, C are all {1, 2, 3}, i.e.

$$ 1 \leq a_{b0} + 2a_{b1} + 4a_{b2} \leq 3; \quad 1 \leq b_{b0} + 2b_{b1} + 4b_{b2} \leq 3; \quad 1 \leq c_{b0} + 2c_{b1} + 4c_{b2} \leq 3. \qquad (24) $$

Next we can follow Encoding 1 to encode the smart table constraint. First we have the encodings of the four necessary auxiliary smart entry flags

$$ e_{A<B} \iff b_{b0} + 2b_{b1} + 4b_{b2} - a_{b0} - 2a_{b1} - 4a_{b2} \geq 1; \qquad (25a) $$

$$ e_{A \in \{1,2\}} \implies \bar{a}_3 \geq 1; \qquad (25b) $$

$$ \bar{e}_{A \in \{1,2\}} \implies \bar{a}_1 + \bar{a}_2 \geq 2; \qquad (25c) $$

$$ f_{A>B} \iff a_{b0} + 2a_{b1} + 4a_{b2} - b_{b0} - 2b_{b1} - 4b_{b2} \geq 1; \qquad (25d) $$

$$ f_{A<B} \iff b_{b0} + 2b_{b1} + 4b_{b2} - a_{b0} - 2a_{b1} - 4a_{b2} \geq 1; \qquad (25e) $$

$$ e_{A=B} \iff \bar{f}_{A>B} + \bar{f}_{A<B} \geq 2; \qquad (25f) $$

$$ e_{B \geq C} \iff b_{b0} + 2b_{b1} + 4b_{b2} - c_{b0} - 2c_{b1} - 4c_{b2} \geq 0; \qquad (25g) $$

Then we can use these to encode the two tuples

$$ s_0 \iff e_{A<B} + e_{A \in \{1,2\}} + c_3 \geq 3; \qquad (25h) $$

$$ s_1 \iff e_{A=B} + \bar{a}_1 + e_{B \geq C} \geq 3; \qquad (25i) $$
<figure><img src="image_265.png" /></figure>

(a) The graph for $P_{\sigma_1}$, which consists of two trees.

<figure><img src="image_266.png" /></figure>

(b) The graph for $P_{\sigma_2}$, which consists of a single tree.

■ **Figure 1** The constraint graphs for the two sub-CSPs $P_{\sigma_1}$ and $P_{\sigma_2}$ represented by the smart tuples $\sigma_1$ and $\sigma_2$. Both graphs are acyclic and thus composed of trees.

And then finally we require

$$s_0 + s_1 \geq 1. \qquad (25j)$$

Constraints that define the $x_{=v}$ flags are also included but are omitted here for brevity.

Now, for the domain-consistent propagator itself, we note that constraint networks of the smart tuples do indeed form acyclic graphs and so the procedure is applicable. This can be seen in Figure 1, where unary constraints are represented as edges to a dummy node. The first tuple is clearly made up of two small trees, and the second consists of a single tree.

Initially all variables have {1, 2, 3} in their domains and the set of unsupported values sl is set to:

$$sl = \{A: 1, 2, 3; B: 1, 2, 3; C: 1, 2, 3\} \qquad (26)$$

To find the set of values supported by $P_{\sigma_1}$ we find the values supported by each tree making it up. This involves two filtering passes over copies of the domain, working first from the root outwards, and then back again.

Taking the left tree first, on the first pass 3 is removed from the domain copy for A and 1 from the domain of B due to $A < B$. So we would log

$$\bar{s}_0 + \bar{a}_{<3} \geq 1; \qquad \bar{s}_0 + \bar{b}_{>1} \geq 1; \qquad (27)$$

as these inferences were made due to a binary constraint. No further inferences are made on the second pass. This leaves the values

$$\{A: 1, 2; B: 2, 3\} \qquad (28)$$

supported by the first tree, and so these are removed from sl. Filtering the other tree simply reduces the domain copy for C to {3}, and so $C = 3$ is removed from sl. No inferences need to be logged here as this is a unary constraint. We now have

$$sl = \{A: 3; B: 1; C: 1, 2\} \qquad (29)$$

Moving on to the second tuple, no inference occurs on the first pass due to $A = B$ nor $B \geq C$, but 1 is removed from the (fresh) domain copy for A due to the unary $A \neq 1$ entry. Then on the second pass, 1 is now removed from the domain copy for B due to $A = B$, and so we log

$$\bar{s}_1 + \bar{b}_{=1} \geq 1; \qquad (30)$$

and then no further inference occurs. So the values occurring in some solution for this tree are

$$\{A: 2, 3; B: 2, 3; C: 1, 2, 3\} \qquad (31)$$
and so these can be removed from sl, leaving

$$sl = \{B : 1\}. \qquad (32)$$

Ultimately, $B = 1$ is the only assignment lacking support in the table after this initial execution of the domain-consistent propagator and the only value to be removed from an actual variable domain. We are then able to log the three constraints in (18) as discussed above, and thus we have successfully maintained the required invariant for pseudo-Boolean proof logging of this propagator.

# 4 Proof Logging for Regular Language Membership Constraints

A regular language provides another compact way of representing a constraint extensionally [19]. It is defined as follows.

## ► Definition 2. Regular Language Membership Constraint

Let $\mathbf{X}$ again be a sequence of finite-domain variables $\langle X_1, \dots, X_n \rangle$ and let $M = (Q, \Sigma, \delta, q_0, F)$ denote a deterministic finite automaton (DFA), where

- $Q$ is a set of states;

- $q_0 \in Q$ is the initial state and $F \subseteq Q$ is the set of accepting states;

- $\Sigma$ is a set of symbols having the domain of every variable in $\mathbf{X}$ as a subset; and

- $\delta : (Q \times \Sigma) \to Q$ is the transition function.

A regular language membership constraint Regular($\mathbf{X}, M$) requires that any sequence of values taken by the variables of $\mathbf{X}$ must belong to the regular language recognised by $M$.

Throughout this section we will assume that the DFA is specified explicitly by a trusted source. Verified compilation to create the constraint from a regular expression is left as future work.

## 4.1 PB Encodings for Regular Language Membership Constraints

The Regular constraint is similarly easy to encode using auxiliary PB variables. If we let $M = (Q, \Sigma, \delta, q_0, F)$ denote a deterministic finite automaton (DFA) and we impose a regular constraint using this on $n$ variables, we simply need to define flags $r_{ij}$ that are set to true if and only if the automaton is in the $j_{th}$ state at position before the $(i+1)_{th}$ symbol of the sequence is processed (including an additional set of flags $r_{nj}$ to denote the final state).

## ► Encoding 2. Regular Language Membership PB Encoding

1: **for all** $i \in \{0, \dots, n\}$

2: $r_{i0} + \dots + r_{i|Q|-1} \ge 1$

3: $-r_{i0} - \dots - r_{i|Q|-1} \ge -1$

4: **for all** $i \in \{0, \dots, n-1\}, q \in Q, v \in \Sigma$

5: **if** $\delta(q, v)$ is an allowed transition

6: $r_{iq} \wedge x_{i=v} \implies r_{i+1\delta(q,v)} \ge 1$

7: **else**

8: $\bar{r}_{iq} + \bar{x}_{i=v} \ge 1$

9: $r_{00} \ge 1$

10: $\sum_{f \in F} r_{(n)f} \ge 1$

This is a simpler version of the stronger SAT encoding for Regular given by Bacchus [1].
## 4.2 Justifying Regular Propagation

As with **SmartTable**, the **Regular** propagator enforces domain consistency, and so our goal once again is to log

$$ \wedge \mathcal{G} \implies \bar{x}_{i=v} \geq 1 \qquad (33) $$

where $\mathcal{G}$ is the current sequences of guesses made by the solver, and $X_i = v$ is an assignment that is not supported by the regular language membership constraint. In particular, $X_i = v$ is not supported when there exists no sequence of symbols belonging to the language recognised by $M$ that has $v$ as the $i_{th}$ symbol.

It is trivial to establish that this inference will not immediately follow by unit propagation, which is unsurprising given the relative simplicity of the encoding. Just as we did for **SmartTable**, we will have to log some additional facts during the execution of the propagation algorithm.

We will consider here the original domain-consistent **Regular** propagator by Pesant [19], on a constraint over a sequence of $n$ variables $\mathbf{X} = \langle X_1, \dots, X_n \rangle$ and associated with a DFA $M = (Q, \Sigma, \delta, q_0, F)$. This works by maintaining a layered, directed multigraph that has $n+1$ layers, with the $i_{th}$ layer having $|Q|$ nodes $q_0^i, \dots, q_{|Q|-1}^i$. The edges in this graph are labelled with values in $\Sigma$, and are required to respect several properties that are satisfied from the outset and then maintained through propagation. These are:

1. Edges can only appear between nodes in consecutive layers.

2. An edge labelled with the value $v$ can only be included between nodes $q_k^i$ and nodes $q_l^{i+1}$ if there exists a transition between the states numbered $k$ and $l$ in the DFA for a $v$, with $v$ currently in the domain of $X_i$.

3. Furthermore, every edge included must appear in a path from the node $q_0^0$ (representing the initial state) to a node representing a final state in the last layer.

The graph is constructed at the start of the solving process, and then updated incrementally during propagation to reflect changes in the variable domains and preserve the required properties. It is clear then that once this is achieved, an assignment $X_i = v$ loses support on the regular constraint exactly when there are no edges labelled with $v$ between nodes in the $i_{th}$ and $(i+1)_{th}$ layers.

The strategy for proof logging is then as follows. Every time an edge $(q_k^i, q_l^{i+1})$ labelled with the symbol $v$ is removed by the propagation algorithm (or excluded when the graph is first built), we log the PB constraint

$$ \wedge \mathcal{G} \implies \bar{r}_{ik} + \bar{x}_{i=v} \geq 1, \qquad (34) $$

which can be interpreted as saying: “given these guesses, we can’t both be in state $k$ after we’ve processed the first $i-1$ symbols, and have the $i_{th}$ symbol be $v$”. If we do this, then when an assignment $X_i = v$ loses support we will have logged a constraint in the form of (34) for all $r_{iq}$ where $\delta(q, v)$ is an allowed transition. Thus, the negation of (33) will result in $\bar{r}_{iq}$ being propagated for all $q \in Q$, and then this will contradict one of the constraints specified by the second line of Encoding 2, which says that at least one $r_{iq}$ must be set to 1. So with these constraints in place we would be able to log our desired constraint (33) by RUP.

The problem is therefore reduced to that of deriving the constraints in the form of (34). It turns out that these will either follow immediately by RUP, or require additional justification depending on how the corresponding edge elimination is inferred by the propagation.

Edges are removed by the **Regular** propagator in one of three ways. Firstly, when an assignment $X_i = v$ has already been removed, either because another value has been guessed or because it has been eliminated by a propagator for another CP constraint, all the
edges extending from the $i_{th}$ layer labelled with value $v$ are also removed. In this case the corresponding constraints in (34) can be logged by RUP, as $x_{i=v}$ will immediately contradict the guesses or previous inferences.

Secondly, when a particular node loses all of its outgoing edges, none of its incoming edges can be part of a valid path and so should be removed. These removals occur recursively in a backward pass (Algorithm 3 from Pesant [19]), so if a removed incoming edge was the last outgoing edge of a node in a previous layer then the incoming edges of that previous node are also removed, and so on. In this case too, the constraints in the form of (34) follow by RUP. This can be shown inductively. Suppose all of the outgoing edges for a node $q_l^i$ have been removed and the corresponding constraints in the form of (34) have already been logged. Then for a given incoming edge $(q_k^{i-1}, q_l^i)$ labelled with the value $v$, the negation of the corresponding constraint in the form of (34) will entail $r_{ik}$, and since $\delta(q_k, v)$ must be $q_l$, $r_{il}$ will then unit propagate by one of the constraint specified by line 6 of Encoding 2. Then, since we have already eliminated all possible outgoing edges, we will propagate $\bar{x}_{i=v'}$ for all $v' \in \Sigma$, which will contradict the encoding of the variable $X_i$, as it must take at least one value.

Finally, in the other direction, when a particular node loses all of its incoming edges, similarly none of its outgoing edges can be part of a valid path and so should be removed. These removals also occur recursively, this time in a forward pass (Algorithm 4 from Pesant[19]). In contrast to the previous two cases, this inference requires more justification as the desired constraint will not always follow by RUP, as is the case for (45) in the worked example below. So when a state $q_l^i$ loses all of its incoming edges we first need to log the constraints

$$ \wedge \mathcal{G} \implies \bar{r}_{i-1k} + \bar{r}_{il} \geq 1 \qquad (35) $$

for each $k \in \{1, \dots, k\}$, which intuitively say that no previous state could lead to this state, before logging constraints in the form of (34) for each outgoing edge. These constraints in the form of (35) will all follow by RUP, since for every $v \in \Sigma$ we either have a constraint

$$ \wedge \mathcal{G} \implies \bar{r}_{i-1k} + \bar{x}_{i-1=v} \geq 1; \qquad (36) $$

from those in the form of (34) logged at the previous stage of the recursion; or there is a constraint

$$ \wedge \mathcal{G} \implies \bar{x}_{i-1=v} \geq 1; \qquad (37) $$

logged during earlier search/propagation (when this node was the first layer of the recursion); or else there is a constraint

$$ \bar{r}_{i-1k} + x_{i-1=v} \geq 1; \qquad (38) $$

present in the PB model. In all cases, the negation of (35) will result in $\bar{x}_{i-1=v}$ for every $v \in \Sigma$, giving a contradiction.

Once all the constraints in the form of (35) have been derived, the desired constraint in the form of (34) will definitely follow by RUP, as its negation would then propagate $\bar{r}_{i-1k}$ for each $k \in Q$, contradicting one of the constraints on the second line of Encoding 2.

Putting all of this together, we can emit full justification using RUP for every edge elimination performed by Regular, and thus we are able to derive the constraint in the form of (33) required by the proof logging invariant.