35:10 Efficient Enumeration of Fixed Points in Complex Boolean Networks

smaller than the number of fixed points. Hence, the proposed method can have memory
benefits as compared to other ASP-based methods. Third, with a much smaller number of
stable models, the solving time can be much smaller. Although the proposed method needs
to spend time for the post-processing, it can have running-time benefits as compared to
other ASP-based methods. The second and third points shall be analyzed in our experiments
shown in Section 5. Note that even if some of the extreme examples of Section 5 do not
correspond to plausible biological processes, finding common patterns among them could give
biological insights. Since the source method can provide a compact representation of the set
of fixed points, it can be used to find such common patterns. For example, we can know that
a key-value of the hash table corresponds to a set of fixed points that only differ in values of
source nodes. This is a very useful result since different source nodes values usually represent
different (environmental) states of a biological system. Some analyses made possible by this
encoding are: listing all fixed points, counting the number of fixed points, listing/counting all
fixed points under a specific value combination of source nodes, listing the core fixed points
projected on only normal nodes (this is exactly the list of key values of the hash table).

5 Experimental results

We implemented the newly proposed methods as a Python package named fASP<fnref n="3" />. For convenience, we name the method based the conjunctive encoding presented in Section 4.1 as fASP-conj and the method based on the source encoding for the case of source nodes presented in Section 4.2 as fASP-src. To evaluate their performance, we compared them with the four state-of-the-art methods for fixed point enumeration in BNs, including PyBoolNet [28], mpbn [34], AN-ASP [1], and FPCollector [7].

We benchmarked all the compared methods on the BBM repository<fnref n="4" />, a collection of
real-world Boolean models from various sources used in systems biology. BBM consists of
211 models, peaking at 321 variables, 1100 regulations, and 133 source nodes, respectively.
Furthermore, we also included a selection of 13 real-world models that are not covered by the
BBM repository. The BNs of this selection peak at 3158 variables, 43642 regulations, and 237
source nodes, respectively. To our knowledge, the BBM repository along with this selected
set is a highly representative sample of Boolean models currently available in the literature.

To solve the ASP problems, we used the same ASP solver Clingo [18] and the same con-
figuration as that used in PyBoolNet, mpbn, and AN-ASP. Specifically for computing maximal
set-inclusion stable models, we used the configuration -heuristic=Domain
-enum-mod=domRec -dom-mod=3 (subset maximality, equivalent to the deprecated
-dom-pref=32 -heuristic=domain -dom-mod=7 used by PyBoolNet). We ran all the bench-
marks on a machine whose environment is CPU: Intel® Core™ i9-11950H 2.60GHz × 16, 16
GB DDR4 RAM, Ubuntu 20.04.5 LTS. Note that for the methods PyBoolNet, mpbn, AN-ASP,
and fASP-conj, we can control the maximum number of fixed points returned because a
resulting stable model corresponds to a single fixed point. In contrast, FPCollector requires
to compute all fixed points, and fASP-src allows its user to set the maximum number of
resulting stable models but not of resulting fixed points. Since a model can have a huge
number of fixed points due to many source nodes, which might not be biologically plausible,
obtaining a sample of those can prove to be very useful to invalidate the model and lead
to further refinement. Hence, to obtain a relevant, reliable and fair comparison, in our

<footnote n="3">³ The source code and benchmarks are freely accessible at https://github.com/giang-trinh/fASP.</footnote>

<footnote n="4">⁴ https://github.com/sybila/biodivine-boolean-models</footnote>
benchmarks, we searched for both all the fixed points and the first 1000 fixed points for each model. In addition, existing analysis shown in the literature usually revolves then around fixing some source nodes to plausible values and reducing the model accordingly. Although this approach biologically makes sense, it relies on potentially arbitrary decisions, and *hides away* critical modelling choices that were actually not part of the original BN. Hence, we did not fix specific values for source nodes in all the considered models. Finally, we set a time limit of two minutes (resp. ten minutes) for each model with respect to enumerating the first 1000 fixed points (resp. all the fixed points).

## 5.1 BBM repository

With regards to enumerating the first 1000 fixed points, the number of BBM models solved within two minutes of each method is: PyBoolNet (198), mpbn (185), AN-ASP (210), and fASP-conj (211). Note that there are 24 non-locally-monotonic models that mpbn cannot handle. Figure 1 shows the cumulative numbers of models solved by the four compared methods with respect to enumerating the first 1000 fixed points. As it can be observed, the AN-ASP and fASP-conj methods are comparable and they vastly outperform the PyBoolNet and mpbn methods. For every time limit, their numbers of solved models are always greater than those of PyBoolNet and mpbn, especially the difference is large for the time limit of 0.5s. In particular, AN-ASP could handle all but one model and fASP-conj could handle all models within 10s.

<figure><img src="image_330.png" /><figcaption>Figure 1 Cumulative numbers of BBM models solved by the four compared methods with respect to enumerating the first 1000 fixed points.</figcaption></figure>

With regards to enumerating all the fixed points, the number of BBM models solved within ten minutes of each method is: FPCollector (172), PyBoolNet (170), mpbn (162), AN-ASP (179), fASP-conj (181), and fASP-src (195). We can see that fASP-src solved more models than all the other methods. One can note that the models for which fASP-src was the only successful method, have extremely large numbers of fixed points, each time because of many source nodes. This confirms our expectation when proposing fASP-src to deal with the case of many source nodes. Upon closer inspection, Figure 2 depicts the cumulative numbers of models solved by the six studied methods with respect to enumerating all the fixed points. As it can be observed, the AN-ASP and fASP-conj methods are still comparable and they outperform the FPCollector, PyBoolNet, and mpbn methods. Furthermore, for every time limit, the number of solved models using fASP-src is always the highest.

It is worth noting that, most models in the BBM repository have moderate numbers of nodes ($n < 200$) and quite simple Boolean functions. Hence, the difference in performance among the three best methods in terms of BBM models (i.e., AN-ASP, fASP-conj, and fASP-src) is not much exhibited. We shall test more in the following subsections.
<figure><img src="image_331.png" /><figcaption>■ **Figure 2** Cumulative numbers of BBM models solved by the six compared methods with respect to enumerating all the fixed points.</figcaption></figure>

## 5.2 Selected real-world models

Table 1 shows the running time of the four compared methods on the 13 selected BNs with respect to enumerating the first 1000 fixed points. Columns *n*, *s*, and |*A*| denote the number of nodes, the number of source nodes, and the number of fixed points, respectively. “DNF” means that the method did not finish the computation within the time limit of two minutes. “NM” indicates a non-locally-monotonic model. We can easily see that *fASP-conj* is the best method in all models as it is much faster than all the other methods on every model. In particular, it only took 35.24s to enumerate the first 1000 fixed points of the hardest model (i.e., the Cell-Cycle-Control model), whereas none of the other methods could finish the computation. Upon closer inspection, *AN-ASP* is the second best method with running time less than 1s in most models. However, it ran quite slowly on the Insulin model, could not handle the Yeast-Pheromone model, and even got an *Out of Memory* (OOM) error on the Cell-Cycle-Control model before finishing the automata network construction. We note that the above problems all come from the bottleneck of *AN-ASP*, i.e., the high number of transitions of the corresponding AN. Finally, consistent with the observations reported in the previous subsection, *PyBoolNet* and *mpbn* still performed much slower than the *AN-ASP* and *fASP-conj* methods, and there are four non-locally-monotonic models (out of the 13) that *mpbn* cannot handle.

Table 2 shows the running time of the six benchmarked methods on the selected BNs with respect to enumerating all the fixed points. Column |*A*| denotes the number of fixed points, and a “?” denotes the case where none of the competing methods returned the result. We can first see that for 6 of the 13 models, none of the competing methods returned all the fixed points. This is clearly because the numbers of fixed points of these models are extremely large. *FPCollector* only succeeded for two models, with each less than 100 nodes. Of course, it is difficult to handle models of larger size. *PyBoolNet* and *mpbn* only succeeded for three models each. All these models are easy for the other methods. Hereafter, we shall present closer inspection on the three most efficient methods: *AN-ASP*, *fASP-conj*, and *fASP-src*.

First, *fASP-conj* is still much faster than *AN-ASP* for all the models where they both succeeded. Second, *fASP-src* is the sole method that could return all the fixed points of the T-Cell-Co-Receptor model within the time limit. Note that the number of fixed points of this model is huge, and it is apparent that none of the other methods can handle it. Moreover,
<table><caption>**Table 1** Timing comparisons (in seconds) among PyBoolNet (PBN), mpbn, AN-ASP, and fASP-conj (f-conj) on the selected models from the literature for enumerating the first 1000 fixed points.</caption><thead><tr><th></th><th>model</th><th>n</th><th>s</th><th>|A|</th><th>PBN</th><th>mpbn</th><th>AN-ASP</th><th>f-conj</th></tr></thead><tbody><tr><td>1</td><td>Cell-Cycle-Control [36]</td><td>3158</td><td>61</td><td>1000<sup>+</sup></td><td>DNF</td><td>DNF</td><td>OOM</td><td>35.24</td></tr><tr><td>2</td><td>EMT-Mechan [38]</td><td>136</td><td>4</td><td>82</td><td>6.02</td><td>0.17</td><td>0.22</td><td>0.05</td></tr><tr><td>3</td><td>EMT-Mechan-TGFbeta [38]</td><td>150</td><td>6</td><td>478</td><td>7.12</td><td>0.39</td><td>0.27</td><td>0.08</td></tr><tr><td>4</td><td>Alzheimer [2]</td><td>762</td><td>237</td><td>1000<sup>+</sup></td><td>DNF</td><td>NM</td><td>0.55</td><td>0.31</td></tr><tr><td>5</td><td>MAPK [2]</td><td>181</td><td>37</td><td>1000<sup>+</sup></td><td>8.81</td><td>0.91</td><td>0.23</td><td>0.09</td></tr><tr><td>6</td><td>Mast-Cell-Activation [2]</td><td>73</td><td>19</td><td>1000<sup>+</sup></td><td>0.07</td><td>0.30</td><td>0.14</td><td>0.03</td></tr><tr><td>7</td><td>Cholocystokinin [2]</td><td>383</td><td>74</td><td>1000<sup>+</sup></td><td>0.60</td><td>1.47</td><td>0.29</td><td>0.13</td></tr><tr><td>8</td><td>HOG [36]</td><td>43</td><td>5</td><td>1000<sup>+</sup></td><td>11.08</td><td>NM</td><td>0.17</td><td>0.07</td></tr><tr><td>9</td><td>Insulin [36]</td><td>82</td><td>7</td><td>1000<sup>+</sup></td><td>DNF</td><td>DNF</td><td>13.65</td><td>0.21</td></tr><tr><td>10</td><td>Leishmania [17]</td><td>342</td><td>81</td><td>1000<sup>+</sup></td><td>DNF</td><td>1.33</td><td>0.37</td><td>0.14</td></tr><tr><td>11</td><td>Pluripotency [49]</td><td>36</td><td>7</td><td>412</td><td>DNF</td><td>NM</td><td>0.22</td><td>0.02</td></tr><tr><td>12</td><td>T-Cell-Co-Receptor [16]</td><td>206</td><td>39</td><td>1000<sup>+</sup></td><td>DNF</td><td>0.79</td><td>0.28</td><td>0.09</td></tr><tr><td>13</td><td>Yeast-Pheromone [36]</td><td>246</td><td>17</td><td>1000<sup>+</sup></td><td>DNF</td><td>NM</td><td>DNF</td><td>1.01</td></tr></tbody></table>

<table><caption>**Table 2** Timing comparisons (in seconds) among FPCollector (FP), PyBoolNet (PBN), mpbn, AN-ASP, fASP-conj (f-conj), and fASP-src (f-src) on the selected models from the literature with respect to enumerating all the fixed points.</caption><thead><tr><th></th><th>model</th><th>|A|</th><th>FP</th><th>PBN</th><th>mpbn</th><th>AN-ASP</th><th>f-conj</th><th>f-src</th></tr></thead><tbody><tr><td>1</td><td>Cell-Cycle-Control</td><td>?</td><td>DNF</td><td>DNF</td><td>DNF</td><td>OOM</td><td>DNF</td><td>DNF</td></tr><tr><td>2</td><td>EMT-Mechan</td><td>82</td><td>DNF</td><td>5.82</td><td>0.15</td><td>0.30</td><td>0.05</td><td>0.08</td></tr><tr><td>3</td><td>EMT-Mechan-TGFbeta</td><td>478</td><td>DNF</td><td>7.05</td><td>0.38</td><td>0.26</td><td>0.08</td><td>0.09</td></tr><tr><td>4</td><td>Alzheimer</td><td>?</td><td>DNF</td><td>DNF</td><td>NM</td><td>OOM</td><td>OOM</td><td>DNF</td></tr><tr><td>5</td><td>MAPK</td><td>?</td><td>DNF</td><td>DNF</td><td>DNF</td><td>OOM</td><td>OOM</td><td>DNF</td></tr><tr><td>6</td><td>Mast-Cell-Activation</td><td>524288</td><td>5.45</td><td>DNF</td><td>145.87</td><td>11.84</td><td>8.96</td><td>6.34</td></tr><tr><td>7</td><td>Cholocystokinin</td><td>?</td><td>DNF</td><td>OOM</td><td>DNF</td><td>OOM</td><td>OOM</td><td>DNF</td></tr><tr><td>8</td><td>HOG</td><td>25632</td><td>149.20</td><td>17.77</td><td>NM</td><td>0.54</td><td>0.79</td><td>1.81</td></tr><tr><td>9</td><td>Insulin</td><td>563200</td><td>DNF</td><td>DNF</td><td>DNF</td><td>81.93</td><td>36.73</td><td>107.60</td></tr><tr><td>10</td><td>Leishmania</td><td>?</td><td>DNF</td><td>DNF</td><td>DNF</td><td>OOM</td><td>OOM</td><td>DNF</td></tr><tr><td>11</td><td>Pluripotency</td><td>412</td><td>1.30</td><td>DNF</td><td>NM</td><td>0.27</td><td>0.05</td><td>0.02</td></tr><tr><td>12</td><td>T-Cell-Co-Receptor</td><td>441039454208</td><td>DNF</td><td>DNF</td><td>DNF</td><td>OOM</td><td>OOM</td><td>349.15</td></tr><tr><td>13</td><td>Yeast-Pheromone</td><td>?</td><td>DNF</td><td>DNF</td><td>NM</td><td>DNF</td><td>DNF</td><td>DNF</td></tr></tbody></table>

both AN-ASP and fASP-conj met the OOM error, which confirms the memory advantage of fASP-src. For the six other models where it succeeded, fASP-src is comparable to AN-ASP and fASP-conj, with a bit slower running time on average. It is apparent because fASP-src suffers from the overhead of its post-processing based on BDDs. Indeed, we confirmed that in most of these models (also of other models considered in our experiments), the ASP solving time of fASP-src is negligible and most of its running time was spent for the post-processing. Note that the running time of this BDD-based post-processing depends on several factors, such as, the number of resulting stable models, the number of source nodes (the number of BDD variables), and the BDD variable ordering.

We also note that for the six failed models, the difference among the three best methods (i.e., AN-ASP, fASP-conj, and fASP-src) is not clear because they all did not finish or met the OOM error. Hence, we conduct new analysis on these failed models by restricting the maximum number of stable models for fASP-src. We then use the number of fixed points obtained by fASP-src as the maximum number of fixed points for the case of AN-ASP or fASP-conj. The experimental settings are the same as those used for the case of enumerating
35:14 Efficient Enumeration of Fixed Points in Complex Boolean Networks

all fixed points. Table 3 shows the results of this analysis. For fASP-src, we consider two maximum numbers: 100 and 200. For each case, Column |A| denotes the number of fixed points obtained by fASP-src, and the other columns denotes the running time of the corresponding methods to obtain that number of fixed points. We can see that for the four models (Cell-Cycle-Control, Alzheimer, Cholocystokinin, and Leishmania), both AN-ASP and fASP-conj met the OOM error, whereas fASP-src can handle them in reasonable time. For the MAPK and Yeast-Pheromone models, fASP-src is much faster than both fASP-conj and AN-ASP. This is also true for the Mast-Cell-Activation and Pluripotency models in the case of enumerating all fixed points (see Table 2). The above observations confirm an advantage in both memory and run-time of using fASP-src in the case of many source nodes.

<table><caption>Table 3 Timing comparisons (in seconds) among AN-ASP, fASP-conj (f-conj), and fASP-src (f-src) on the six failed models.</caption><thead><tr><th rowspan="2">model</th><th rowspan="2"></th><th colspan="4">100</th><th colspan="4">200</th></tr><tr><th>|A|</th><th>f-src</th><th>AN-ASP</th><th>f-conj</th><th>|A|</th><th>f-src</th><th>AN-ASP</th><th>f-conj</th></tr></thead><tbody><tr><td>1</td><td>Cell-Cycle-Control</td><td>159744</td><td>21.98</td><td>OOM</td><td>OOM</td><td>323584</td><td>23.73</td><td>OOM</td><td>OOM</td></tr><tr><td>2</td><td>Alzheimer</td><td>&gt; 10<sup>32</sup></td><td>10.89</td><td>OOM</td><td>OOM</td><td>&gt; 10<sup>32</sup></td><td>29.93</td><td>OOM</td><td>OOM</td></tr><tr><td>3</td><td>MAPK</td><td>31744</td><td>0.19</td><td>1.74</td><td>1.37</td><td>52480</td><td>0.42</td><td>2.77</td><td>2.11</td></tr><tr><td>4</td><td>Cholocystokinin</td><td>&gt; 10<sup>13</sup></td><td>0.28</td><td>OOM</td><td>OOM</td><td>&gt; 10<sup>13</sup></td><td>0.55</td><td>OOM</td><td>OOM</td></tr><tr><td>5</td><td>Leishmania</td><td>&gt; 10<sup>13</sup></td><td>0.46</td><td>OOM</td><td>OOM</td><td>&gt; 10<sup>13</sup></td><td>0.79</td><td>OOM</td><td>OOM</td></tr><tr><td>6</td><td>Yeast-Pheromone</td><td>1552</td><td>0.91</td><td>DNF</td><td>1.14</td><td>3152</td><td>0.91</td><td>DNF</td><td>1.45</td></tr></tbody></table>

5.3 Pseudo-random models

The results on the 224 real-world models reported in the two previous subsections draw a quite clear picture about the performance of the six compared methods. However, we observed that there is only one model with more than 1000 nodes (i.e., the Cell-Cycle-Control model). Moreover, in most of them, the Boolean functions are quite simple, even sometimes just simple conjunctions/disjunctions of literals. The reason for these facts may be that the modelers were restricted by the limited performance of the tools supported at the time they created the models. This is not the case of the Cell-Cycle-Control model, since its authors only conducted simulations instead of formal analysis [36]. Since in the present work we target large and complex BNs, we set out to test the performance of our proposed methods on larger and more complicated models than the ones available in the literature to date. Specifically, we wanted to test models with 1000 or more nodes and Boolean functions in complicated forms. Such a model is arguably not possible to achieve yet with hand-made modeling, even with a fully or semi-automated inference technique [2], but might be in the near future.

To create a benchmark set of larger and more complex models, we decided to generate pseudo-random models following the generation approach proposed by the research group who created and is maintaining the BBM repository. This generation approach is described in detail in [9] and its implementation is provided at https://github.com/daemontus/artifact_cav2021. In general, it generates Boolean models structurally similar to the real-world models in the BBM repository. To ensure this structural similarity, the generator uses a node-degree distribution sampled from the BBM repository, as opposed to other theoretical random network models. Once the regulators of a node are specified, its Boolean function is generated by randomly choosing between ∧ and ∨ when connecting the positive/negative
literals of the regulators. Note that this option does not cover the full spectrum of possible Boolean functions, but it can make the generated Boolean functions complicated enough for evaluation.

In the end, we created 400 pseudo-random models ranging from 1000 to 5000 variables, 4145 to 63507 regulations, and 127 to 1171 source nodes, respectively. We then tested all the competing methods on these models. We first reported that all the competing methods failed to obtain all the fixed points as they quickly met the OOM error. The reason is that the number of all fixed points is actually too large due to a lot of source nodes (> 100). Hence, we here only searched for the first 1000 fixed points, which might also be more biologically relevant. The time limit for each model was set to two minutes. The other settings are the same as those used in the two previous subsections.

<figure><img src="image_332.png" /><figcaption>Figure 3 Cumulative numbers of pseudo-random models solved by the competing methods with respect to enumerating the first 1000 fixed points.</figcaption></figure>

With regards to enumerating the first 1000 fixed points, the number of pseudo-random models solved within two minutes of each competing method is: PyBoolNet (0), mpbn (338), AN-ASP (13), and fASP-conj (400). PyBoolNet could not handle any model, and it failed at the phase of computing prime implicants in most cases. This is not surprising since the models are large in size and the formulae quite complex. Interestingly, mpbn could handle far more models than AN-ASP. This can be explained by the fact that the number of transitions in the corresponding AN is very large in most models, whereas the size of the DNF for each Boolean function is still moderate. Moreover, the models are all locally-monotonic, which is an assumption of the generator [9]. Upon closer inspection, Figure 3 depicts the cumulative numbers of pseudo-random models solved by the four competing methods with respect to enumerating the first 1000 fixed points. We can see that fASP-conj is the best method as it vastly outperforms the three other methods for every time limit except the time limit of 0.5s where all the methods could handle no model. In particular, it could handle all the 400 pseudo-random models within 50s.

# 6 Conclusion and future work

In this work we have proposed two new ASP-based methods called fASP-conj and fASP-src for efficiently enumerating fixed points of Boolean networks, which are crucial in modeling and analysis of biological systems. fASP-conj is based on conjunctive ASP and fASP-src is a modification of fASP-conj to handle the case of a large number of source nodes. The
35:16 Efficient Enumeration of Fixed Points in Complex Boolean Networks

main advantage of these methods is that they only rely on NNFs of Boolean functions, which
are much more efficient to obtain than other representations used by previous methods (e.g.,
prime-implicants, disjunctive normal forms, automata networks). The main advantage of
fASP-src is that it provides a more compact representation of the results based on BDDs,
which can give both memory and run-time benefits. We have also formally proved the
correctness of the above new methods.

We have then benchmarked their performance against the four other state-of-the-art tools: FPCollector, PyBoolNet, mpbn, and AN-ASP. The experimental results on both real-world and pseudo-random models show that the new methods vastly outperform the state-of-the-art as they can robustly handle various types of large and complex models, whereas the other methods cannot. In particular, they can handle models of up to 5000 nodes with very complicated Boolean update functions. In terms of enumerating the first 1000 fixed points (resp. all fixed points), the experimental results show that fASP-conj is the best (resp. second best) method. fASP-src, which is based on fASP-conj, shows its superiority to all the other methods in enumerating all the fixed points of models with many source nodes.

Boolean network models of biological systems usually contain many source nodes, which
might be hard to avoid in the modeling process [2]. Currently, there are many such models
that fASP-src cannot handle. Hence, improving fASP-src is necessary. Note that, in some
cases, the number of auxiliary atoms in the core encoding of fASP-conj and fASP-src can
be reduced. Such optimization will be looked into in the future. Furthermore, we also plan
to extend the methods proposed in this present paper to those for computing trap spaces of
Boolean networks, which are more general than fixed points and useful approximations for
complex attractors in Boolean networks. It is crucial because the state-of-the-art methods
for the trap space computation are all unable to robustly handle large and complex models,
for instance, the models used in our experiments here.

References

1 Emna Ben Abdallah, Maxime Folschette, Olivier F. Roux, and Morgan Magnin. ASP-based method for the enumeration of attractors in non-deterministic synchronous and asynchronous multi-valued networks. *Algorithms Mol. Biol.*, 12(1):20:1–20:23, 2017. doi: 10.1186/s13015-017-0111-2.

2 Sara Sadat Aghamiri, Vidisha Singh, Aurélien Naldi, Tomás Helikar, Sylvain Soliman, Anna Niarakis, and Jinbo Xu. Automated inference of Boolean models from molecular interaction maps using CaSQ. *Bioinform.*, 36(16):4473–4482, 2020. doi:10.1093/bioinformatics/btaa484.

3 Tatsuya Akutsu, Satoru Kuhara, Osamu Maruyama, and Satoru Miyano. A system for identifying genetic networks from gene expression patterns produced by gene disruptions and overexpressions. *Genome Informatics*, 9:151–160, 1998. doi:10.11234/gi1990.9.151.

4 Tatsuya Akutsu, Avraham A. Melkman, Takeyuki Tamura, and Masaki Yamamoto. Determining a singleton attractor of a Boolean network with nested canalyzing functions. *J. Comput. Biol.*, 18(10):1275–1290, 2011. doi:10.1089/cmb.2010.0281.

5 Tatsuya Akutsu, Yang Zhao, Morihiro Hayashida, and Takeyuki Tamura. Integer programming-based approach to attractor detection and control of Boolean networks. *IEICE Trans. Inf. Syst.*, 95-D(12):2960–2970, 2012. doi:10.1587/transinf.E95.D.2960.

6 Julio Aracena. Maximum number of fixed points in regulatory Boolean networks. *Bull. Math. Biol.*, 70:1398–1409, 2008. doi:10.1007/s11538-008-9304-7.

7 Julio Aracena, Luis Cabrera-Crot, and Lilian Salinas. Finding the fixed points of a Boolean network from a positive feedback vertex set. *Bioinform.*, 37(8):1148–1155, 2021. doi:10.1093/bioinformatics/btaa922.
8 Ferhat Ay, Günhan Gülsoy, and Tamer Kahveci. Finding steady states of large scale regulatory networks through partitioning. In *International Workshop on Genomic Signal Processing and Statistics*, pages 1–4. IEEE, 2010. doi:10.1109/GENSIPS.2010.5719669.

9 Nikola Benes, Lubos Brim, Samuel Pastva, and David Safránek. Computing bottom SCCs symbolically using transition guided reduction. In *International Conference on Computer Aided Verification*, pages 505–528. Springer, 2021. doi:10.1007/978-3-030-81685-8_24.

10 Nikolaos Berntenis and Martin Ebeling. Detection of attractors of large Boolean networks via exhaustive enumeration of appropriate subspaces of the state space. *BMC Bioinform.*, 14(1):1–10, 2013. doi:10.1186/1471-2105-14-361.

11 R. E. Bryant. Graph-based algorithms for Boolean function manipulation. *IEEE Trans. Comput.*, 35(8):677–691, 1986.

12 Stéphanie Chevalier, Vincent Noël, Laurence Calzone, Andrei Yu. Zinovyev, and Loïc Paulevé. Synthesis and simulation of ensembles of Boolean networks for cell fate decision. In *International Conference on Computational Methods in Systems Biology*, pages 193–209. Springer, 2020. doi:10.1007/978-3-030-60327-4_11.

13 Karla Fabiola Corral-Jara, Camille Chauvin, Wassim Abou-Jaoudé, Maximilien Grandclaudon, Aurélien Naldi, Vassili Soumelis, and Denis Thieffry. Interplay between SMAD2 and STAT5A is a critical determinant of IL-17A/IL-17F differential expression. *Mol. Biomed.*, 2(1):9, 2021. doi:10.1186/s43556-021-00034-3.

14 Steve Dworschak, Susanne Grell, Victoria J. Nikiforova, Torsten Schaub, and Joachim Selbig. Modeling biological networks by action languages via answer set programming. *Constraints An Int. J.*, 13(1-2):21–65, 2008. doi:10.1007/s10601-007-9031-y.

15 Swann Floc'Hlay, Maria Dolores Molina, Céline Hernandez, Emmanuel Haillot, Morgane Thomas-Chollier, Thierry Lepage, and Denis Thieffry. Deciphering and modelling the TGF-β signalling interplays specifying the dorsal-ventral axis of the sea urchin embryo. *Dev.*, 148(2):dev189944, 2021. doi:10.1242/dev.189944.

16 Piyali Ganguli, Saikat Chowdhury, Rupa Bhowmick, and Ram Rup Sarkar. Temporal protein expression pattern in intracellular signalling cascade during T-cell activation: A computational study. *J. Biosci.*, 40(4):769–789, September 2015. doi:10.1007/s12038-015-9561-1.

17 Piyali Ganguli, Saikat Chowdhury, Shomeek Chowdhury, and Ram Rup Sarkar. Identification of Th1/Th2 regulatory switch to promote healing response during leishmaniasis: a computational approach. *EURASIP J. Bioinform. Syst. Biol.*, 2015:13, 2015. doi:10.1186/s13637-015-0032-7.

18 Martin Gebser, Benjamin Kaufmann, Roland Kaminski, Max Ostrowski, Torsten Schaub, and Marius Schneider. Potassco: The Potsdam answer set solving collection. *AI Commun.*, 24(2):107–124, 2011. doi:10.3233/AIC-2011-0491.

19 Michael Gelfond and Vladimir Lifschitz. The stable model semantics for logic programming. In *International Conference and Symposium on Logic Programming*, pages 1070–1080. MIT Press, 1988.

20 Carlos Gershenson. Classification of random Boolean networks. In *Artificial Life VIII, Proceedings of the Eighth International Conference on Artificial Life*, pages 1–8, 2003.

21 Trinh Van Giang, Tatsuya Akutsu, and Kunihiko Hiraishi. An FVS-based approach to attractor detection in asynchronous random Boolean networks. *IEEE ACM Trans. Comput. Biol. Bioinform.*, 19(2):806–818, 2022. doi:10.1109/TCBB.2020.3028862.

22 Leon Glass and Stuart A Kauffman. The logical analysis of continuous, non-linear biochemical control networks. *J. Theor. Biol.*, 39(1):103–129, 1973. doi:10.1016/0022-5193(73)90208-7.

23 Eric Goles and Lilian Salinas. Sequential operator for filtering cycles in Boolean networks. *Adv. Appl. Math.*, 45(3):346–358, 2010. doi:10.1016/j.aam.2010.03.002.

24 A Gonzalez Gonzalez, Aurélien Naldi, Lucas Sanchez, Denis Thieffry, and Claudine Chaouiya. GINsim: a software suite for the qualitative modelling, simulation and analysis of regulatory networks. *Biosyst.*, 84(2):91–100, 2006. doi:10.1016/j.biosystems.2005.10.003.
35:18 Efficient Enumeration of Fixed Points in Complex Boolean Networks

25 Changki Hong, Jeewon Hwang, Kwang-Hyun Cho, and Insik Shin. An efficient steady-state analysis method for large Boolean networks with high maximum node connectivity. *PLoS One*, 10(12):e0145734, December 2015. doi:10.1371/journal.pone.0145734.

26 Katsumi Inoue. Logic programming for Boolean networks. In *International Joint Conference on Artificial Intelligence*, pages 924–930. IJCAI/AAAI, 2011. doi:10.5591/978-1-57735-516-8/ IJCAI11-160.

27 Roland Kaminski, Torsten Schaub, Anne Siegel, and Santiago Videla. Minimal intervention strategies in logical signaling networks with ASP. *Theory Pract. Log. Program.*, 13(4-5):675–690, 2013. doi:10.1017/S1471068413000422.

28 Hannes Klarner, Adam Streck, and Heike Siebert. PyBoolNet: a python package for the generation, analysis and visualization of Boolean networks. *Bioinform.*, 33(5):770–772, 2017. doi:10.1093/bioinformatics/btw682.

29 Tomoya Mori and Tatsuya Akutsu. Attractor detection and enumeration algorithms for Boolean networks. *Comput. Struct. Biotechnol. J.*, 20:2512–2520, 2022. doi:10.1016/j.csbj. 2022.05.027.

30 Mushthofa Mushthofa, Gustavo Torres, Yves Van de Peer, Kathleen Marchal, and Martine De Cock. ASP-G: an ASP-based method for finding attractors in genetic regulatory networks. *Bioinform.*, 30(21):3086–3092, 2014. doi:10.1093/bioinformatics/btu481.

31 Aurélien Naldi. BioLQM: a Java toolkit for the manipulation and conversion of logical qualitative models of biological networks. *Front. Physiol.*, 9:1605, 2018. doi:10.3389/fphys. 2018.01605.

32 Aurélien Naldi, Pedro T. Monteiro, Christoph Müssel, Hans A. Kestler, Denis Thieffry, Ioannis Xenarios, Julio Saez-Rodriguez, Tomás Helikar, and Claudine Chaouiya. Cooperative development of logical modelling standards and tools with CoLoMoTo. *Bioinform.*, 31(7):1154–1159, 2015. doi:10.1093/bioinformatics/btv013.

33 Karen J Nuñez-Reza, Aurélien Naldi, Arantza Sánchez-Jiménez, Ana V Leon-Apodaca, M Angélica Santana, Morgane Thomas-Chollier, Denis Thieffry, and Alejandra Medina-Rivera. Logical modelling of in vitro differentiation of human monocytes into dendritic cells unravels novel transcriptional regulatory interactions. *Interface Focus*, 11(4):20200061, 2021. doi:10.1098/rsfs.2020.0061.

34 Loïc Paulevé, Juraj Kolčák, Thomas Chatain, and Stefan Haar. Reconciling qualitative, abstract, and scalable modeling of biological networks. *Nat. Commun.*, 11(1), August 2020. doi:10.1038/s41467-020-18112-5.

35 Alexandre Rocca, Nicolas Mobilia, Eric Fanchon, Tony Ribeiro, Laurent Trilling, and Katsumi Inoue. ASP for construction and validation of regulatory biological networks. *Logical Modeling of Biological Systems*, pages 167–206, 2014.

36 Jesper C Romers and Marcus Krantz. rxncon 2.0: a language for executable molecular systems biology. *BioRxiv*, page 107136, 2017. doi:10.1101/107136.

37 Torsten Schaub and Sven Thiele. Metabolic network expansion with answer set programming. In *International Conference on Logic Programming*, pages 312–326. Springer, 2009. doi: 10.1007/978-3-642-02846-5_27.

38 Emmalee Sullivan, Marlayna Harris, Arnav Bhatnagar, Eric Guberman, Ian Zonfa, and Erzsébet Ravasz Regan. Boolean modeling of mechanosensitive Epithelial to Mesenchymal Transition and its reversal. *bioRxiv*, September 2022. doi:10.1101/2022.08.29.505701.

39 René Thomas. Boolean formalisation of genetic control circuits. *J. Theor. Biol.*, 42:565–583, 1973. doi:10.1016/0022-5193(73)90247-6.

40 René Thomas. Regulatory networks seen as asynchronous automata: a logical description. *J. Theor. Biol.*, 153(1):1–23, 1991. doi:10.1016/S0022-5193(05)80350-9.

41 René Thomas and Richard d'Ari. *Biological feedback*. CRC press, 1990.

42 Van-Giang Trinh, Kunihiko Hiraishi, and Belaid Benhamou. Computing attractors of large-scale asynchronous Boolean networks using minimal trap spaces. In *ACM International Conference*
on *Bioinformatics*, Computational Biology and Health Informatics, pages 13:1–13:10. ACM, 2022. doi:10.1145/3535508.3545520.

43 Eirini Tsirvouli, Felicity Ashcroft, Berit Johansen, and Martin Kuiper. Logical and experimental modeling of cytokine and eicosanoid signaling in psoriatic keratinocytes. *iScience*, 24(12):103451, 2021. doi:10.1016/j.isci.2021.103451.

44 Leslie G Valiant. The complexity of enumeration and reliability problems. *SIAM J. Comput.*, 8(3):410–421, 1979.

45 Alan Veliz-Cuba, Boris Aguilar, Franziska Hinkelmann, and Reinhard C. Laubenbacher. Steady state analysis of Boolean molecular network models via model reduction and computational algebra. *BMC Bioinform.*, 15(1):1–8, 2014. doi:10.1186/1471-2105-15-221.

46 Santiago Videla, Carito Guziolowski, Federica Eduati, Sven Thiele, Martin Gebser, Jacques Nicolas, Julio Saez-Rodriguez, Torsten Schaub, and Anne Siegel. Learning Boolean logic models of signaling networks with ASP. *Theor. Comput. Sci.*, 599:79–101, 2015. doi: 10.1016/j.tcs.2014.06.022.

47 Santiago Videla, Julio Saez-Rodriguez, Carito Guziolowski, and Anne Siegel. caspo: a toolbox for automated reasoning on the response of logical signaling networks families. *Bioinform.*, 33(6):947–950, 2017. doi:10.1093/bioinformatics/btw738.

48 Rui-Sheng Wang, Assieh Saadatpour, and Reka Albert. Boolean modeling in systems biology: an overview of methodology and applications. *Phys. Biol.*, 9(5):055001, 2012. doi:10.1088/1478-3975/9/5/055001.

49 Ayako Yachie-Kinoshita, Kento Onishi, Joel Ostblom, Matthew A Langley, Eszter Posfai, Janet Rossant, and Peter W Zandstra. Modeling signaling-dependent pluripotency with Boolean logic to predict cell fate transitions. *Mol. Syst. Biol.*, 14(1):e7952, 2018. doi: 10.15252/msb.20177952.
<figure><img src="image_333.png" /></figure>
# Guided Bottom-Up Interactive Constraint Acquisition

Dimosthenis C. Tsouros ✉️ ID

KU Leuven, Belgium

Senne Berden ✉️ ID

KU Leuven, Belgium

Tias Guns ✉️ ID

KU Leuven, Belgium

## Abstract

Constraint Acquisition (CA) systems can be used to assist in the modeling of constraint satisfaction problems. In (inter)active CA, the system is given a set of candidate constraints and posts queries to the user with the goal of finding the right constraints among the candidates. Current interactive CA algorithms suffer from at least two major bottlenecks. First, in order to converge, they require a large number of queries to be asked to the user. Second, they cannot handle large sets of candidate constraints, since these lead to large waiting times for the user. For this reason, the user must have fairly precise knowledge about what constraints the system should consider. In this paper, we alleviate these bottlenecks by presenting two novel methods that improve the efficiency of CA. First, we introduce a bottom-up approach named GROWACQ that reduces the maximum waiting time for the user and allows the system to handle much larger sets of candidate constraints. It also reduces the total number of queries for problems in which the target constraint network is not sparse. Second, we propose a probability-based method to guide query generation and show that it can significantly reduce the number of queries required to converge. We also propose a new technique that allows the use of openly accessible CP solvers in query generation, removing the dependency of existing methods on less well-maintained custom solvers that are not publicly available. Experimental results show that our proposed methods outperform state-of-the-art CA methods, reducing the number of queries by up to 60%. Our methods work well even in cases where the set of candidate constraints is 50 times larger than the ones commonly used in the literature.

**2012 ACM Subject Classification** Computing methodologies → Discrete space search; Computing methodologies → Supervised learning by classification; Theory of computation → Constraint and logic programming; Computing methodologies → Active learning settings

**Keywords and phrases** Constraint acquisition, Constraint learning, Active learning, Modelling

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

**Supplementary Material** Software: https://github.com/Dimosthenis/ActiveConLearn

**Funding** This research received funding from the European Research Council (ERC) under the EU Horizon 2020 research and innovation programme (Grant No 101002802, CHAT-Opt) and from the European Union's Horizon 2020 research and innovation programme under grant agreement No 101070149, project Tuples

# 1 Introduction and related work

Constraint programming (CP) is considered one of the main paradigms for solving combinatorial problems, with many successful applications in a variety of domains. However, there are still challenges to be faced in order for CP technology to become even more widely used. One of the most important challenges is to ease the modeling process. The current assumption in CP is that the user first models the problem and that a solver is then used to

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

© Dimosthenis C. Tsouros, Senne Berden, and Tias Guns; 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. 36; pp. 36:1–36:20

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

LIPICS

Leibniz International Proceedings in Informatics

Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl Publishing, Germany
solve it. However, modeling is a non-trivial task. Expressing a combinatorial problem as
a set of constraints over decision variables is not straightforward and requires substantial
expertise [12]. As a result, modeling is considered a major bottleneck for the widespread
adoption of CP [12, 13, 14].

This obstacle has led to research into a very different approach to modeling: that of
*learning* the constraint problem from data, as opposed to manually constructing it. This is the
focus of the research area of *constraint acquisition (CA)*, in which CP meets machine learning.
In CA, the model of a constraint problem is acquired (i.e., learned) (semi-)automatically from
a set of examples of solutions, and possibly non-solutions. CA methods can be categorized
as *active* or *passive* on the basis of whether a user provides feedback during learning or not.

In *passive acquisition*, a dataset of examples of solutions and non-solutions is provided by the user upfront. Based on these examples, the system learns a set of constraints modeling the problem [4, 5, 7, 9, 11, 17, 18, 19, 21]. Approaches vary in the types of constraints they are able to learn and the methodologies they employ: CONACQ.1 is a version space algorithm for learning fixed-arity constraints [7, 9, 11], ModelSeeker learns global constraints that are taken from a predefined constraint catalog [4], and COUNT-CP is a generate-and-aggregate approach that can learn expressive first-order constraints [18]. None of these approaches are robust to errors in the labeled data. To this end, SEQACQ and BAYESACQ were introduced, being robust to noise in the training set. In SEQACQ, a statistical approach based on sequential analysis is used [22], while in BAYESACQ, a naive Bayes classifier is trained, from which a constraint network is then derived [23].

In contrast to passive learning, *active* or *interactive acquisition* systems learn the con-
straints through interaction with the user, by asking queries. The main type of query used is
the *membership query*, which asks the user to classify a given example (i.e., an assignment
to the variables of the problem) as a solution or a non-solution. An early work in active
CA is the Matchmaker agent [15], where users, when they answer a membership query
negatively, also have to provide a violated constraint. In order to lower the expertise level
required from the user, Bessiere et al. later proposed CONACQ.2 [10, 11] – an active version
of CONACQ.1 that uses membership queries and does not require the user to provide any
violated constraints. In [24], CONACQ.2 was in turn extended to also accept arguments
regarding *why* examples should be rejected or accepted.

As the number of membership queries needed can be exponentially large for these
methods [11], a new family of interactive algorithms was proposed that use *partial queries*
instead [3, 8, 20, 25, 26, 27, 28, 29]. A partial query asks the user to classify a partial
assignment to the variables. Using partial queries, CA systems are able to converge faster.
QUACQ was the first system to use partial queries [6, 8], and was later extended into
MULTIACQ [3]. MQUACQ was later introduced to reduce the number of queries needed per
learned constraint [25, 29], and MQUACQ-2 further improved the performance by exploiting
the structure of the constraints already learned [27].

Despite these advancements in active CA, there are still significant obstacles for the
technology to become usable in practice. One of the main limitations is that it typically still
requires asking a large number of queries to the user in order to find all constraints. In addition,
existing systems cannot handle large sets of candidate constraints in reasonable run times, and
thus require significant expertise from the user in limiting the constraints the system should
consider (and thus the size of the candidate set) upfront. Finally, query generation – a highly
important part of the CA process – currently requires the use of customized solvers that are
not publicly available and are not as well-maintained as conventional solvers. Without the
use of such customized solvers, current active CA algorithms can lead to very high query
generation times or are sometimes unable to converge to the correct set of constraints when
time limits are imposed [1, 25].
We focus on the above limitations, and contribute the following improvements:

- We present a novel query generation method named PQ-GEN that allows conventional constraint solvers to be used by CA algorithms while also ensuring convergence, removing the dependency on customized solvers.

- We propose a bottom-up learning approach named GROWACQ that uses any other CA algorithm to learn the constraints of an increasingly large problem. It starts learning with only a subset of variables and an associated subset of candidate constraints, and incrementally grows this set of variables and constraints. This allows it to handle significantly larger sets of candidate constraints and reduces the maximum waiting time for the user.

- Finally, we introduce a better way to guide the query generation process, with the goal of generating queries that learn the set of constraints faster. We propose an objective function for query generation that uses probabilistic estimates of whether constraints are likely to hold or not. We demonstrate the potential of this method by using a simple counting-based approach as probabilistic estimator.

The rest of the paper is structured as follows. Some background on CA is given in Section 2. Sections 3–5 present our proposed methods. An experimental evaluation is given in Section 6. Finally, Section 7 concludes the paper.

# 2 Background

We now introduce some basic notions regarding constraint satisfaction problems and interactive constraint acquisition.

## 2.1 Constraint satisfaction problems

A constraint satisfaction problem (CSP) is a triple $P = (X, D, C)$, consisting of:

- a set of $n$ variables $X = \{x_1, x_2, ..., x_n\}$, representing the entities of the problem,

- a set of $n$ domains $D = \{D_1, D_2, ..., D_n\}$, where $D_i \subset \mathbb{Z}$ is the finite set of values for $x_i$,

- a constraint set (also called constraint network) $C = \{c_1, c_2, ..., c_t\}$.

A constraint $c$ is a pair $(rel(c), var(c))$, where $var(c) \subseteq X$ is the scope of the constraint and $rel(c)$ is a relation over the domains of the variables in $var(c)$ that specifies (implicitly or explicitly) what assignments are allowed. $|var(c)|$ is called the *arity* of the constraint. The constraint set $C[Y]$, where $Y \subseteq X$, denotes the set of constraints from $C$ whose scope is a subset of $Y$. The set of solutions of a constraint set $C$ is denoted by $sol(C)$. A *redundant* or *implied* constraint $c \in C$ is a constraint in $C$ such that $sol(C) = sol(C \setminus \{c\})$.

A (partial) assignment $e_Y$ is an assignment over a set of variables $Y \subseteq X$. $e_Y$ is rejected by a constraint $c$ iff $var(c) \subseteq Y$ and the projection $e_{var(c)}$ of $e_Y$ on the variables in the scope $var(c)$, is not in $rel(c)$, that is, is not allowed by the constraint. $\kappa_C(e_Y)$ represents the subset of constraints from $C[Y]$ that reject $e_Y$, i.e., $\kappa_C(e_Y) = \{c \mid c \in C[Y] \wedge e_{var(c)} \notin rel(c)\}$.

A complete assignment $e$ that is accepted by all the constraints in $C$ is a *solution* to $C$, i.e., $e \in sol(C)$. A partial assignment $e_Y$ is called a *partial solution* to $C$ iff it is accepted by all the constraints in $C[Y]$. Note that a partial solution to $C$ may not be extendable to a complete one, due to constraints not in $C[Y]$.
36:4 Guided Bottom-Up Interactive Constraint Acquisition

■ **Algorithm 1** Constraint Acquisition through partial queries.

**Input:** X, D, B, C<sub>in</sub> (X: the set of variables, D: the set of domains, B: the bias, C<sub>in</sub>: an optional set of known constraints)

**Output:** C<sub>L</sub>: a learned constraint network

1: C<sub>L</sub> ← C<sub>in</sub>

2: **while** True **do**

3: Generate an *e* accepted by C<sub>L</sub> and rejected by B

4: **if** *e* = nil **then** return C<sub>L</sub> ▷ Stopping condition

5: **if** ASK(*e*) = Yes **then** ▷ Ask (partial) membership query *e*

6: Remove the constraints rejecting *e*, namely κ<sub>B</sub>(*e*), from B

7: **else**

8: Find one (or more) minimal scopes S in *e* for which |κ<sub>B</sub>(*e*<sub>S</sub>)| ≥ 1 and ASK(*e*<sub>S</sub>) = No

9: Find all {c ∈ C<sub>T</sub> | var(c) = S} through partial queries; add to C<sub>L</sub>, remove from B

2.2 Active constraint acquisition with partial membership queries

In CA, the pair (X, D) is called the *vocabulary* of the problem at hand and is common knowledge shared by the user and the system. Besides the vocabulary, the learner is also given a *language* Γ consisting of *fixed-arity* constraint relations. Using the vocabulary (X, D) and the constraint language Γ, the system generates the *constraint bias B*, which is the set of all possible candidate constraints for the problem.

Let $C_T$, the target constraint network, be an unknown set of constraints such that for every assignment $e$ over $X$ it holds that $e \in \text{sol}(C_T)$ iff $e$ is a solution to the problem the user has in mind. The goal of CA is to learn a constraint set $C_L \subseteq B$ that is equivalent to $C_T$. Like other works, we assume that the bias $B$ can represent $C_T$, i.e., there exists a $C \subseteq B$ s.t. $\text{sol}(C) = \text{sol}(C_T)$.

In active CA, the system interacts with the user while learning the constraints. A *membership query* [2] in this setting is a question ASK(*e*<sub>X</sub>), asking the user whether a complete assignment *e*<sub>X</sub> is a solution to the problem that the user has in mind. A *partial query* ASK(*e*<sub>Y</sub>), with *Y* ⊂ *X*, asks the user to determine if *e*<sub>Y</sub>, which is an assignment in *D*<sup>Y</sup>, is a partial solution with respect to *C*<sub>T</sub>[*Y*]. We use the notation *c* ∈ *C*<sub>T</sub> iff ∀*e* ∈ *D*<sup>Y</sup> with *var*(*c*) ⊆ *Y* ⊂ *X*, ASK(*e*<sub>Y</sub>) = True ⇒ *e*<sub>var(*c*)</sub> ∈ *sol*(*c*).

While in passive acquisition there are methods that can handle noisy answers [22, 23], this is not the case for active acquisition. For this reason, in this work, we follow the assumption that the user answers all queries correctly.

A query ASK(e<sub>Y</sub>) is called *irredundant* iff the answer is not implied by any information already available to the system. That is, the query is *irredundant* iff e<sub>Y</sub> is rejected by at least one constraint from the bias B and is not rejected by the network C<sub>L</sub> learned thus far. The first condition captures that κ<sub>B</sub>(e<sub>Y</sub>) cannot be empty, since if κ<sub>B</sub>(e<sub>Y</sub>) would be empty, the answer to the query ASK(e<sub>Y</sub>) would have to be “yes”, based on the assumption that C<sub>T</sub> is representable by the constraints in B. The second condition captures that e<sub>Y</sub> should not be rejected by any constraint in the learned network C<sub>L</sub>, since otherwise the user would certainly answer “no” to the query.

Algorithm 1 presents the generic process followed by active CA methods with partial queries. The learned set $C_L$ is first initialized either to the empty set or to a set of constraints given by the user that is known to be part of $C_T$ (i.e., $C_{in} \subset C_T$) (line 1). Then the main
loop of the acquisition process begins, where, in every iteration, the system first generates an irredundant query (line 3) and posts it to the user (line 5). If the query is answered positively, then the candidate constraints from B that violate it are removed (line 6). Otherwise, the system has to find one or more constraints from $C_T$ that violate the query. This is done in two steps. First, queries are asked to find the scope of a constraint in $\kappa_{C_T}(e)$ (line 8). Then, queries are asked to find all constraints $c \in C_T$ with that scope (line 9).

The acquisition process has converged on the learned network $C_L \subseteq B$ iff $C_L$ agrees with the set of all labeled examples E, and for every other network $C \subseteq B$ that agrees with E, it holds that $sol(C) = sol(C_L)$. This is proved if no query could be generated at line 3, as in this case, all remaining constraints in B (if any) are redundant. If the first condition is true but the second condition has not been proved when the acquisition process finishes, premature convergence has occurred. This can happen when the query generation at line 3 returns $e = nil$, but without having proved that an irredundant query does not exist (e.g., because of a time limit).

Existing algorithms like QUACQ [6, 8], MQUACQ [25, 29] and MQUACQ-2 [27] follow this template, but differ mainly in how they implement lines 3, 8 and 9, and hence how many constraints they are able to learn in each iteration. Examples of functions used to locate the scope of a constraint (line 8) are FindScope [6, 8] or the more efficient FindScope-2 [25]. To learn the constraints in the scope found (line 9), the FindC function is typically used [6, 8].

# 3 Using conventional solvers for query generation

Query generation (line 3 of Algorithm 1) is one of the most important parts of the CA process. It aims to find an irredundant membership query (i.e., a (partial) assignment that does not violate $C_L$ but violates at least one $c \in B$) that will be asked to the user. Thus, it can be formalized as follows:

$$ \text{find } e_Y \text{ s.t. } e_Y \in sol(C_L[Y] \wedge \bigvee_{c_i \in B[Y]} \neg c_i), $$

which can be formulated as a CSP with variables Y and constraints $C_L[Y] \wedge \bigvee_{c_i \in B[Y]} \neg c_i$.

## 3.1 Problems when using conventional solvers

In principle, this CSP could be solved using any conventional CP solver. However, this can lead to issues for the following two reasons.

**A large bias.** At the start of the acquisition process, the set of candidate constraints B can be very large. This makes the propagation of the constraint $\bigvee_{c_i \in B[Y]} \neg c_i$ time-consuming, and severely slows down the query-generation process.

**Indirectly implied constraints.** At the end of the acquisition process, only constraints that are implied by $C_L$ remain in B, if any. In this case, it will be impossible to generate a query that does not violate $C_L$ and violates at least one constraint from B. However, propagation is often unable to prove such implications when they are indirect and involve multiple variables and constraints. For this reason, solvers internally end up enumerating all possible variable assignments satisfying $C_L$ and checking if the constraint $\bigvee_{c_i \in B[Y]} \neg c_i$ can be satisfied. This can be very time-consuming, and a time limit is usually imposed on query generation, leading to premature convergence.
36:6 Guided Bottom-Up Interactive Constraint Acquisition

Algorithm 2 PQ-GEN: Projection-based Query Generation.

Input: $C_L, B, l, t$ ($B$: the set of candidate constraints (bias), $C_L$: set of known constraints, $l$: size limit, $t$: time limit)

Output: $e$: the query generated

1: timer.start()

2: $Y \leftarrow \bigcup_{c \in B} \text{var}(c)$

3: **if** $|B| > l$ **then**

4: &nbsp;&nbsp;&nbsp;&nbsp;$e \leftarrow \text{solve}(C_L[Y])$

5: &nbsp;&nbsp;&nbsp;&nbsp;**if** $\exists c \in B: e \notin \text{sol}(c)$ **then**

6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**return** $e$

7: &nbsp;&nbsp;&nbsp;&nbsp;$e \leftarrow \text{solve}(C_L[Y] \wedge \bigvee_{c_i \in B[Y]} \neg c_i)$

8: &nbsp;&nbsp;&nbsp;&nbsp;**if** $\text{timer.end()} < t$ **then**

9: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$e' \leftarrow \text{solve}(C_L[Y] \wedge \bigvee_{c_i \in B[Y]} \neg c_i)$, maximize: $\text{obj}$, time limit:$t - \text{timer.end()}$)

10: &nbsp;&nbsp;&nbsp;&nbsp;**if** $e' \neq \text{nil}$ **then**

11: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**return** $e'$

12: **return** $e$

In order to limit the large runtimes in a more advanced way than by simply imposing a time bound $t$, Addi et al. proposed a method using conventional solvers named TQ-GEN [1]. It iteratively tries to solve the query generation problem, by gradually reducing the number of variables taken into account by a proportion $\alpha \in ]0, 1[$, until a query can be generated within a small time limit $\tau$. This is repeated until either an irredundant query is generated, or a global time bound $t$ is reached, leading to premature convergence. However, choosing the right hyperparameters for $t$ and $\tau$ is problem-specific [1] and requires tuning, and thus more interaction with the user.

3.2 Customized solvers

To avoid premature convergence, a CP solver can be customized to store partial assignments that satisfy every $c \in C_L[Y]$ and violate at least one $c \in B[Y]$ during the search. Given an objective function, such as maximizing the number of assigned variables, in every non-failing node of the search tree it will check the above property and, if fulfilled, store the best-scoring partial assignment.

As these customized solvers are guaranteed to find valid partial solutions, their use will never lead to premature convergence. In addition, finding a partial query to return is not time-consuming (especially when combined with specialized search heuristics [25]), even when the bias is large. However, such custom solvers are not publicly available and are typically not based on the latest version of state-of-the-art solvers. This also means that the corresponding active CA methods are heavily tied to those particular customized solvers.

3.3 Projection-based Query Generation

We now introduce a method named *Projection-based Query Generation (PQ-GEN)* that makes it possible to use state-of-the-art conventional solvers for query generation, without premature convergence. Our proposed method is shown in Algorithm 2.

**Avoiding indirectly implied constraints.** A key observation we make is that when generating a query on line 3, it might be that $\bigcup_{c \in B} \text{var}(c) \subset X$, that is, some variables have no more candidate constraints in $B$. These have become irrelevant, as both lines 6 and 8 are only
concerned with $\kappa_B(e)$, which will not include these variables. So, to generate an irredundant query, it is sufficient to consider only the variables in $B$. This is not only faster, but also avoids indirectly implied constraints, as these are indirect through variables not used in $B$.

Thus, our proposed query generator projects the variables down to $Y \subseteq X$, with $Y = \bigcup_{c \in B} \text{var}(c)$, thereby simplifying the problem to finding an assignment over $Y \subseteq X$. This will inherently result in a partial assignment when $Y$ is a strict subset of $X$, without requiring a custom solver. Thus, we first compute the set of variables $Y$ relevant to the query (line 2), and project $C_L$ down to those variables (on lines 4 and 7). The solver then has to prove that there exists a query that satisfies $C_L[Y]$ and violates at least one constraint from $B$.

**Dealing with large biases.** As mentioned above, having a large bias $B$ can severely slow down the solver during query generation because propagating the $\bigvee_{c_i \in B[Y]} \neg c_i$ constraint takes a long time. However, we observe that when $B$ contains many constraints, the property that a query $e$ violates at least one of these is usually satisfied without needing to enforce this. Hence, we propose not using this constraint when the bias is larger than some threshold (lines 3 to 6 in Algorithm 2). If in a post-hoc check, it turns out that the generated query violates at least one $c \in B$, it is directly returned (line 6). Otherwise, we again generated a query, this time with the constraint enforcing that there must exist a constraint in $B$ that is violated (line 7).

**Optimizing the query.** The above ensures that we will always find a valid query. However, much better queries – according to some objective function – can often be found. This would take additional time, but is safe because, since a valid query has already been found, the optimization can always safely be interrupted. Given a time limit, we can hence call an optimization solver for the remaining time after a first valid query has been found (lines 8-11).

As expressed in Proposition 1, Algorithm 2 is correct.

**► Proposition 1.** Given a bias $B$, with an unknown target network $C_T$ being representable by $B$, and a learned constraint set $C_L$, if `nil` is returned by Algorithm 2, then the system has converged on $C_T[X]$.

*Proof.* When `nil` is returned by Algorithm 2, it means that $\not\exists e \in \text{sol}(C_L[Y] \wedge \bigvee_{c_i \in B[Y]} \neg c_i)$, with $Y = \bigcup_{c \in B} \text{var}(c)$, i.e., $\not\exists e \in \text{sol}(C_L[Y] \wedge \bigvee_{c_i \in B[Y]} \neg c_i)$. In order to prove convergence over all of $X$, we must have $\not\exists e \in \text{sol}(C_L[X] \wedge \bigvee_{c_i \in B[X]} \neg c_i)$. We will now show that when $Y = \bigcup_{c \in B} \text{var}(c)$, it means that

$$ \not\exists e \in \text{sol}(C_L[Y] \wedge \bigvee_{c_i \in B[Y]} \neg c_i) \implies \not\exists e \in \text{sol}(C_L[X] \wedge \bigvee_{c_i \in B[X]} \neg c_i) $$

Assume that Algorithm 2 returns `nil`, i.e., that no assignment exists in a $Y \subset X$ that is accepted by $C_L[Y]$ and rejected by $B[Y]$. This means that all the constraints in $B[Y]$ are proved to be implied by the constraints in $C_L[Y]$. Thus, the remaining constraints in $B$, that are not proved to be redundant, are the constraints $c \in B \setminus B[Y]$. When we know that $Y = \bigcup_{c \in B} \text{var}(c)$ it means that $B[Y] = B$, so $B \setminus B[Y] = \emptyset$. As a result, in this case, all the constraints in $B$ are proved to be implied. Hence, no assignment that is accepted by $C_L$ and rejected by $B$ exists in $X$. ◀

# 4 Bottom-up Constraint Acquisition

We start by observing that all current active CA algorithms always consider either the full set of variables $X$, or a large subset $Y \subseteq X$, in their top-level loop (lines 2-9 in Algorithm 1). This generally leads to complete or almost-complete queries getting generated (line 3 of
■ Algorithm 3 Growing Acquisition.

**Input:** Γ, X, D, C<sub>in</sub> (Γ: the language, X: the set of variables, D: the set of domains, C<sub>in</sub>: an optional set of known constraints)
**Output:** C<sub>L</sub>: a constraint network

1: C<sub>L</sub> ← ∅

2: Y ← ∅

3: while |Y| ≤ |X| do

4: x ← x ∈ (X \ Y)

5: Y ← Y ∪ {x}

6: B ← {c | rel(c) ∈ Γ ∧ var(c) ⊆ Y ∧ x ∈ var(c)}

7: C<sub>L</sub> ← Acq(Y, D<sup>Y</sup>, B, C<sub>L</sub> ∪ C<sub>in</sub>[Y])

8: return C<sub>L</sub>

Algorithm 1). However, larger queries are generally harder to answer than smaller queries [25]. Also, a large initial query leads to many additional queries getting posed in the scope-finding method on line 8. That is because the worst-case complexity of the best scope-finding methods, in terms of the number of queries required, is $Θ(\log(|Y|))$, where $Y \subseteq X$ is the set of variables considered [25].

Additionally, by directly considering the whole set of variables, the CA algorithm has to represent and operate on the entire set of candidate constraints (i.e., the bias B) at once. The bias is used in many parts of the acquisition process. Hence, the memory requirements and the run time of the acquisition process increase significantly as the bias grows, either because the problems contain more variables or because the language Γ given to the system includes a larger number of relations. This means that, in practice, state-of-the-art active CA methods are only applicable to problems with not too many variables or problems for which the user already has relatively precise knowledge about what constraints the system should consider (which corresponds to the bias being small).

To improve on this, we propose a novel *meta-algorithm* named GROWACQ (Algorithm 3). The key idea is to call a CA algorithm on an increasingly large subset of the variables $Y \subseteq X$, each time using only a relevant unexplored subset of the bias. GROWACQ begins with $Y = \emptyset$ (line 2) and gradually incorporates more variables (lines 3-5). Once a new variable $x_i \in X$ has been added to $Y$, the new problem becomes to find the new $C_T[Y]$. However, as $C_T[Y \setminus \{x_i\}]$ was already found in the previous iterations, the set of constraints to seek is actually $C_T[Y] \setminus C_T[Y \setminus \{x_i\}]$. To find $C_T[Y] \setminus C_T[Y \setminus \{x_i\}]$, any existing active CA algorithm can be used. We represent this with the function Acq (line 7). In every iteration, only a part of the bias B is needed, namely $B[Y] \setminus B[Y \setminus \{x_i\}]$, and as shown in Lemma 1, the bias constructed at line 6 is equivalent to $B[Y] \setminus B[Y \setminus \{x_i\}]$.

**► Lemma 1.** Let $Y_i$ be the set of variables $Y$ in iteration $i$ after line 5 of Algorithm 3 and $B_i = \{c \mid \text{rel}(c) \in \Gamma \wedge \text{var}(c) \subseteq Y_i \wedge x_i \in \text{var}(c)\}$ be the bias constructed at line 6 in iteration $i$. It holds that $B_i = B[Y_i] \setminus B[Y_{i-1}]$.

**Proof.** At line 6 of Algorithm 3, the bias B is constructed. For each iteration $i$, it is constructed as $B_i = \{c \mid \text{rel}(c) \in \Gamma \wedge \text{var}(c) \subseteq Y_i \wedge x_i \in \text{var}(c)\}$. For a set of variables $Y_i$, the full bias, which includes all candidate constraints, is $B[Y_i] = \{c \mid \text{rel}(c) \in \Gamma \wedge \text{var}(c) \subseteq Y_i\}$. For the previous iteration, as $Y_{i-1} = Y_i \setminus \{x_i\}$, we know that $B[Y_{i-1}] = \{c \mid \text{rel}(c) \in \Gamma \wedge \text{var}(c) \subseteq Y_i \setminus \{x_i\}\}$. Thus, the additional constraints that are in $B[Y_i]$ and not in $B[Y_{i-1}]$
are the ones with a scope $var(c) \subseteq Y_i$ for which $x_i \in var(c)$:

$$ 
\begin{aligned}
B[Y_i] \setminus B[Y_{i-1}] &= \{c \mid rel(c) \in \Gamma \land var(c) \subseteq Y_i\} \setminus \{c \mid rel(c) \in \Gamma \land var(c) \subseteq Y_i \setminus \{x_i\}\} \\
&= \{c \mid rel(c) \in \Gamma \land var(c) \subseteq Y_i \land x_i \in var(c)\} = B_i
\end{aligned}
$$

Hence, it holds that $B_i = B[Y_i] \setminus B[Y_{i-1}]$ for $B_i = \{c \mid rel(c) \in \Gamma \land var(c) \subseteq Y_i \land x_i \in var(c)\}$. ◀

This bottom-up approach alleviates the problems described above, i.e., starting from large initial queries and having to represent the whole bias from the beginning, in two ways. First, it naturally leads to partial queries of increasing size in the first step of the “inner” CA system (Algorithm 1 line 5). This is valuable since smaller queries are generally easier for the user to answer [25], and also a smaller initial query leads to a lower worst-case number of additional queries to locate scopes. Second, since the algorithm only stores and uses a small part of the bias at a time (line 6 of Algorithm 3), it is able to handle significantly larger biases than the state-of-the-art. Not representing the whole bias in every iteration does not affect the algorithm’s correctness, as we state in Proposition 2.

**Proposition 2.** Given a bias $B$ built from a language $\Gamma$, with bounded arity constraints, and a target network $C_T$ representable by $B$, $GROWACQ$ is correct (i.e., will learn a constraint set $C_L$ that is equivalent to $C_T$), as long as a correct (i.e., sound and complete) CA algorithm is used in line 7.

**Proof.** (Sketch)

Let us now prove that if any correct algorithm is used in line 7 of Algorithm 3 – like QUACQ, MQUACQ or MQUACQ-2 – GROWACQ remains correct. We will subscript sets with the number of the iteration that they occur in to distinguish between the iterations. Even though the full bias $B$ is never constructed and never kept in memory all at once in GROWACQ, we will still refer to it in this proof and denote it with $B$, i.e., $B = \{c \mid rel(c) \in \Gamma \land var(c) \subseteq X\}$. When we instead write $B_i$, we refer to the part of the bias that is constructed and used in iteration $i$ (line 6 of Algorithm 3), which is $B[Y_i] \setminus B[Y_{i-1}]$ (Lemma 1).

**Soundness.** $GROWACQ$ adds constraints to $C_L$ only at line 7 of Algorithm 3. At that line, only constraints returned from the inner interactive CA algorithm are added to $C_L$. Since the assumption is that a sound algorithm is used in the Acq function, $GROWACQ$ is sound.

**Completeness.** We prove that $GROWACQ$ is complete by proving by induction that, after each iteration $i$, $C_L$ is equivalent to $C_T[Y_i]$, meaning that after the last iteration, $C_L$ is equivalent to $C_T[X]$. $GROWACQ$ starts with $Y_1 = \emptyset$, so both $C_T[Y_1]$ and $B_1$ are empty. The first iteration where the algorithm has to actually learn any constraints will be the one where $Y$ grows large enough so that $C_T[Y] \neq \emptyset$. Assume that this happens at iteration $k$. In this case, $C_T[Y_k]$ will be representable by $B_k$, because $B_k = B[Y_k] \setminus B[Y_{k-1}]$ and we know that $C_T[Y_{k-1}] = \emptyset$. Since $C_T[Y_k]$ is representable by $B_k$, it will be successfully learned in line 7, as long as a complete interactive CA algorithm is used.

Assuming now that $C_L = C_T[Y_n]$ holds at the end of the $n$-th iteration, let us now prove that $C_L = C_T[Y_{n+1}]$ will hold at the end of the $n+1$-th iteration. From the assumption that $C_L = C_T[Y_n]$, it follows that $(B[Y_n] \setminus C_L) \cap C_T = \emptyset$. As a result, $B_{n+1}$, being equal to $B[Y_{n+1}] \setminus B[Y_n]$ does not exclude any constraint from $C_T[Y_{n+1}]$ that has not already been learned. From this, it follows that $(C_T[Y_{n+1}] \setminus C_L) \subseteq B_{n+1}$, and thus this set of constraints will be learned in line 7 as long as a complete interactive CA algorithm is used. Hence, $GROWACQ$ is complete. ◀
## 5 Guided query generation

We now turn our attention to the objective function used at line 9 of Algorithm 2. Since when GROWACQ is used, the size of B used in every iteration is reduced, query generation is now often fast, leaving sufficient room for using optimization to find a good query.

The objective function used in existing query generation systems [6, 25] tries to maximize the number of constraints from B that are violated by the generated query e. The motivation is that this can potentially help shrink the bias faster. The objective function is

$$e = \arg \max_{e} \sum_{c \in B} [e \notin sol(\{c\})]$$

where [·] is the Iverson bracket which converts True/False into 1/0.

However, looking only at the number of violated constraints in B does not fully capture what a good query is:

- We want queries that lead to a positive answer to violate many constraints from the bias B, as these can then all be removed from B, shrinking it faster.

- On the other hand, we want queries that lead to a negative answer to violate a small number of constraints from B, as it allows the CA system to find the conflicting constraint faster.

Based on this, in order to generate good queries regardless of the user's answer, we want query generation to minimize the violation of constraints that are in the unknown target set $C_T$, seeking a query to which the user's answer will be “yes”. At the same time, we want to maximize the violation of constraints in B that are not in $C_T$, so that positive answers can shrink the bias faster (the first bullet point above). Note that we also have the constraint ensuring that at least one constraint from B has to be violated. This means that when $B \setminus C_T = \emptyset$, we want a minimum number of constraints in $C_T$ that we have not already learned to be violated. This leads to negative queries that violate a small number of constraints in B (the second bullet point above).

Assume we have access to an oracle $\mathcal{O}$ that tells us whether a constraint c belongs to the unknown target set or not: $\mathcal{O}(c) = (c \in C_T)$. Using this oracle we can formulate an objective function for query generation, using the reasoning above, as follows:

$$\sum_{c \in B} [e \notin sol(\{c\})] \cdot (1 - |\Gamma| \cdot [\mathcal{O}(c)])$$

On the one hand, every time that the oracle returns False for a constraint from the bias that is violated by e, the objective function is increased by 1, thereby maximizing the violation of these constraints. Conversely, for constraints where $\mathcal{O}$ returns True, we aim to minimize the violations, which requires a reduction in the objective value for each such violated constraint. However, it is possible that violating a set of constraints C (where $\forall c_i \in C | \mathcal{O}(c_i) = False$) may imply the violation of a constraint $c_j$ with $\mathcal{O}(c_j) = True$. In such cases, if the reduction in the objective value for violating $c_j$ is not large enough, the system will violate both C and $c_j$, maximizing the objective. To address this issue, we introduce a “penalty” of $|\Gamma|$, which is equal to the upper bound of the number of constraints in each scope. This ensures that the system prioritizes satisfying a constraint with $\mathcal{O}(c_j) = True$, over violating other constraints from B.

**Modeling the oracle.** Observe how the current objective of maximizing violations corresponds to using a model of the oracle M that always answers False, i.e., that assumes that none of the candidate constraints belong to $C_T$. On the other hand, if we used an oracle M
that always answers True, then the query generation would try to violate as few constraints as possible. However, the $\bigvee_{c_i \in B[Y]} \neg c_i$ constraint would still need to be satisfied, in the extreme case leading every query to violate exactly one constraint from B. Based on this observation, we propose to model the oracle using the following model M, which tries to determine for every constraint c whether violating or satisfying c would lead to the least amount of queries later on in the algorithm.

$$M(c) = \left( \frac{1}{P[c \in C_T]} \leq \log(|Y|) \right)$$

On the one hand, in the extreme case, the constraints for which $M(c)$ answers True will be violated one by one in the later queries (once most of the constraints for which $M(c)$ answers False have been dealt with). Let $P[c \in C_T]$ be a probabilistic estimate of whether c is part of $C_T$. Then, if the generated queries would violate the constraints with that probability one by one, we would in expectation need $1/P[c \in C_T]$ queries to find a constraint from $C_T$. For example, for a set of constraints that each has a probability of 25%, 1 in every 4 queries is expected to lead to a $c \in C_T$ being learned.

On the other hand, for each constraint $c \in C_T$ for which $M(c)$ answers False, a scope-finding procedure is needed to locate the violated constraint. The most efficient functions commonly used to do it (i.e., FindScope [6] or FindScope-2 [25]) have been shown to require $\Theta(\log(|Y|))$ queries to find a violated constraint $c \in C_T$ in the worst case, where Y is the number of variables considered in query generation. As a result, we estimate the number of queries needed in this case as $k \cdot \log(|Y|)$, with k a constant. We found $k = 1$ to work well in practice.

**Probability estimation.** To compute the probability $P(c \in C_T)$ of a constraint $c \in B$, we use a simple approach, considering only information from the relations $rel(c)$ of the constraints. More specifically, to compute $P(c \in C_T)$, we count the number of times a constraint with relation $rel(c)$ has been added to $C_L$, and divide it by the total number of times that such a constraint has been removed from B. Much more advanced estimation techniques, including machine learning methods, can be used for more accurate estimation. We leave this for future work.

# 6 Experimental evaluation

In this section, we empirically answer the following research questions:

(Q1) Does using PQ-GEN with conventional solvers avoid premature convergence, and how do CA systems perform when they use it?

(Q2) Does GROWACQ (using MQUACQ-2) perform better than using MQUACQ-2 directly?

(Q3) How does our probability-guided query generation objective function perform compared to the one used in current CA systems?

(Q4) How does the combination of our methods perform?

(Q5) How do our methods perform on problems with a huge bias B?

## 6.1 Benchmarks

We used the following benchmarks:

**Jigsaw Sudoku.** The Jigsaw Sudoku is a variant of Sudoku in which the 3 × 3 boxes are replaced by irregular shapes. It consists of 81 variables with domains of size 9. The target network consists of 811 binary $\neq$ constraints, on rows, columns, and shapes. The bias B was constructed using the language $\Gamma = \{\geq, \leq, <, >, \neq, =\}$ and contains 19 440 binary constraints.
36:12 Guided Bottom-Up Interactive Constraint Acquisition

**Murder.** The Murder puzzle problem consists of 20 variables with domains of size 5. The target network contains 4 cliques of 10 ≠ constraints and 12 additional binary constraints. The bias was initialized with 760 constraints based on the language $\Gamma = \{\ge, \le, <, >, \ne, =\}$.

**Random.** We used a problem with 100 variables and domains of size 5. We generated a random target network with 495 ≠ constraints. The bias was initialized with 19800 constraints, using the language $\Gamma = \{\ge, \le, <, >, \ne, =\}$.

**Golomb rulers.** The problem is to find a ruler where the distance between any two marks is different from that between any other two marks. We built a simplified version of a Golomb ruler with 8 marks, with the target network consisting only of quaternary constraints.<fnref n="1" /> The bias, consisting of 238 binary and quaternary constraints, was created with the language $\Gamma = \{\ge, \le, <, >, \ne, =, |x_i - x_j| \ne |x_k - x_l|\}$.

**Job-shop scheduling.** The job-shop scheduling problem involves scheduling a number of jobs, consisting of several tasks, across a number of machines, over a certain time horizon. The decision variables are the start and end times of each task. There is a total order over each job's tasks, expressed by binary precedence constraints. There are also constraints capturing the duration of the tasks and that tasks should not overlap on the same machine. The language $\Gamma = \{\ge, \le, <, >, \ne, =, x_i + c = x_k\}$ was used, with $c$ being a constant from 0 up to the maximal duration of the jobs. We used a problem instance containing 10 jobs, 3 machines (i.e., $|X| = 60$) and a time horizon of 15 steps, leading to a bias containing 14160 constraints.

## 6.2 Experimental setup

Let us now give some details about the experimental settings:

- All the experiments were conducted on a system carrying an Intel(R) Core(TM) i9-11900H, 2.50 GHz clock speed, with 16 GB of RAM.

- We measure the total number of queries $\#q$, the average time of the query generation process $\bar{T}_{gen}$ (line 3 of Algorithm 1), the average waiting time $\bar{T}$ per query for the user, and the total time needed (to converge) $T_{total}$. All times are presented in seconds. The difference between $\bar{T}_{gen}$ and $\bar{T}$ is that the latter takes into account also the queries posed on lines 8-9 of Algorithm 1, which are very fast to compute.

- We evaluate our methods in comparison with the state-of-the-art method MQuAcQ-2 [27].

- All methods and benchmarks were implemented in Python<fnref n="2" /> using the CPMpy constraint programming and modeling library [16], except for the experiments using custom solvers.<fnref n="3" />

- The results presented in each benchmark, for each algorithm, are the means of 10 runs.

We now discuss the results of our experimental evaluation, based on the questions we posed at the beginning of the section.

## 6.3 [Q1] Performance of PQ-Gen

Both PQ-GEN, our projection-based query generation approach, and TQ-GEN [1] (discussed in Section 3.1) involve hyperparameters that affect their performance. Thus, we first performed a hyperparameter sensitivity analysis to assess their performance under different

<footnote n="1">¹ The ternary constraints derived when $i=k$ or $j=l$ in $|x_i - x_j| \neq |x_k - x_l|$ were excluded, as also done in the literature [25, 27]</footnote>

<footnote n="2">² Our code is available online in: https://github.com/Dimosts/ActiveConLearn</footnote>

<footnote n="3">³ For the custom solver based query generators from [6, 25], we obtained the implementations (in C++) through personal communication with the authors.</footnote>
configurations. In tandem with TQ-GEN, we also used the adjust function described in [1]. We used the JSudoku benchmark for this comparison. For TQ-GEN, we fixed the hyperparameter $\alpha$ to 0.8 as recommended in [1], and used $\tau = \{0.05, 0.1, 0.2, 0.3\}$ and $t = \{0.5, 1, 1.5, 2\}$. For PQ-GEN hyperparameters, we used $l = \{3000, 5000, 7500, 10000\}$ and $t = \{0.5, 1, 1.5, 2\}$. Thus, we examined 16 different configurations for each. A summary of the results are shown in Table 1.<fnref n="4" />

<table><caption>**Table 1** A summary of the performance of TQ-GEN and PQ-GEN with different configurations.</caption><thead><tr><th>Problem</th><th>Conv</th><th>#q</th><th>T<sub>max</sub></th><th>T<sub>total</sub></th></tr></thead><tbody><tr><td>MQuAcq-2 with TQ-GEN [1]</td><td>32%</td><td>7555</td><td>20.66</td><td>2371.40</td></tr><tr><td>MQuAcq-2 with PQ-GEN (ours)</td><td>100%</td><td>6551</td><td>4.42</td><td>728.25</td></tr></tbody></table>

Confirming our analysis, with our PQ-GEN there is never a case of premature convergence, no matter what hyperparameters are used. On the other hand, when TQ-GEN is used, the system fails to converge in the majority of cases, and specific hyperparameter values have to be chosen to ensure convergence. In addition, our PQ-GEN shows much better performance both in terms of the number of queries needed and, especially, in terms of runtime.

In more detail, we compared our projection-based query generation (PQ-GEN) with a baseline where we run a conventional CP solver to directly solve the query generation problem, using a one-hour time limit, as well as with query generation methods from the literature, i.e., TQ-GEN and the custom solver based query generators from [6, 27]. For PQ-GEN and TQ-GEN we used the best configuration found in the previous experiment. That is, we run PQ-GEN with $l = 5000$ and $t = 1$ and TQ-GEN with $\tau = 0.2$ and $t = 2$. We used benchmarks that are similar to the ones used in [6, 27]. For consistency, we used the same state-of-the-art query generation objective function across all methods that accept one, i.e., our PQ-GEN and the custom solvers, which tries to maximize the number of violated constraints from B. The results are shown in Table 2.

We can observe that convergence was reached in all cases, except for the baseline in which a conventional solver was used directly using a time limit. Our method, PQ-GEN, and the baseline show similar performance in terms of the number of queries needed, while being much better than TQ-GEN in JSudoku and Random, where B is larger, especially when considering time performance.

On the other hand, when custom solvers are used, we can see that the time performance has improved and the number of queries has decreased. This happens because the custom solver can return a partial assignment of any size, trying only to maximize the value of the objective function used, and utilizing heuristics from the literature, while when our PQ-GEN is used, the query generated has to be a solution in a specific (sub)set of variables, which takes more time to compute. As a result, we observed that custom solvers often return queries that violate more constraints from B, which helps MQUAcq-2 shrink the bias faster in terms of the number of queries needed.

## 6.4 [Q2 - Q4] Evaluating GrowAcq and guided query generation

Hereafter, we continue our experiments using PQ-GEN, as the motivation was to investigate techniques that work with any solver. As using PQ-GEN allows us to use conventional solvers, able to run on any given benchmark, in contrast to the custom solvers from [6, 27], where specific constraint relations are implemented, from now on we will use all of the benchmarks mentioned in Section 6.1.

<footnote n="4">⁴ More details regarding this experiment can be found in Appendix A.</footnote>
36:14 Guided Bottom-Up Interactive Constraint Acquisition

<table><caption>**Table 2** Comparing PQ-GEN with state-of-the-art query generators.</caption><thead><tr><th>Method</th><th>Problem</th><th>#q</th><th>T&#773;<sub>gen</sub></th><th>T&#773;</th><th>T<sub>max</sub></th><th>T<sub>total</sub></th><th>Convergence</th></tr></thead><tbody><tr><td colspan="8" style="text-align:center; font-weight:bold;">Using conventional solvers</td></tr><tr><td rowspan="3">MQUACQ-2 <sub>BASELINE</sub></td><td>JSudoku</td><td>6 337</td><td>2.05</td><td>0.21</td><td>11.67</td><td>-</td><td>0%</td></tr><tr><td>Murder</td><td>347</td><td>0.09</td><td>0.01</td><td>0.31</td><td>4.69</td><td>100%</td></tr><tr><td>Random</td><td>5 694</td><td>2.90</td><td>0.05</td><td>20.78</td><td>294.94</td><td>100%</td></tr><tr><td rowspan="3">MQUACQ-2 <sub>TQ-GEN [1]</sub></td><td>JSudoku</td><td>7 153</td><td>0.26</td><td>0.13</td><td>8.46</td><td>919.34</td><td>100%</td></tr><tr><td>Murder</td><td>394</td><td>0.04</td><td>0.01</td><td>0.32</td><td>5.22</td><td>100%</td></tr><tr><td>Random</td><td>5 787</td><td>5.26</td><td>1.23</td><td>17.61</td><td>7100.44</td><td>100%</td></tr><tr><td rowspan="3">MQUACQ-2 <sub>PQ-GEN (OURS)</sub></td><td>JSudoku</td><td>6 458</td><td>0.77</td><td>0.10</td><td>3.19</td><td>666.91</td><td>100%</td></tr><tr><td>Murder</td><td>370</td><td>0.66</td><td>0.03</td><td>1.10</td><td>12.28</td><td>100%</td></tr><tr><td>Random</td><td>5 708</td><td>0.60</td><td>0.04</td><td>2.22</td><td>233.62</td><td>100%</td></tr><tr><td colspan="8" style="text-align:center; font-weight:bold;">Using custom solvers</td></tr><tr><td rowspan="3">MQUACQ-2 <sub>GenerateQuery.cutoff [6]</sub></td><td>JSudoku</td><td>5 321</td><td>0.99</td><td>0.06</td><td>2.04</td><td>336.42</td><td>100%</td></tr><tr><td>Murder</td><td>421</td><td>0.76</td><td>0.13</td><td>1.02</td><td>53.21</td><td>100%</td></tr><tr><td>Random</td><td>5 349</td><td>0.94</td><td>0.04</td><td>3.24</td><td>198.65</td><td>100%</td></tr><tr><td rowspan="3">MQUACQ-2 <sub>max<sub>B</sub> [25]</sub></td><td>JSudoku</td><td>5 042</td><td>1.00</td><td>0.06</td><td>2.34</td><td>277.36</td><td>100%</td></tr><tr><td>Murder</td><td>325</td><td>0.86</td><td>0.04</td><td>1.01</td><td>12.90</td><td>100%</td></tr><tr><td>Random</td><td>5 012</td><td>0.94</td><td>0.02</td><td>1.52</td><td>95.21</td><td>100%</td></tr></tbody></table>

**[Q2] Using GrowAcq within MQuAcq-2.** We now evaluate the performance of GROWACQ, our proposed bottom-up CA approach. To evaluate it, we used MQUACQ-2, as the inner CA algorithm within GROWACQ (line 7 of Algorithm 3) and compared this to using MQUACQ-2, directly on the full-sized problem. Table 3, top two blocks, presents the results.

We can observe that the usage of GROWACQ results in a reduction of the number of queries in JSudoku, Murder, and Random, while a slight increase can be seen in Golomb. In Job-shop, the increase in the number of queries is somewhat larger (25%). This is the case because the target constraint network in this benchmark is sparse, with most of the iterations of GROWACQ in a $Y \subset X$ not learning any constraint from $C_T$ and only shrinking the bias. So, when the full-sized problem is looked at directly when MQUACQ-2 is used, the bias $B$ can shrink with fewer queries. On the other hand, when the target network is not sparse, there is a decrease in the number of queries of up to 19%, due to the fact that the system can locate the scopes of the constraints faster, starting from a $Y \subset X$ every time. Based on the above observations, we can see that using GROWACQ leads to learning constraints in a lower amount of queries, but on the other hand, needs more queries to shrink the bias.

Finally, although the total time is almost the same in most problems, and slightly increased in JSudoku and Golomb, the average time per query has not noticeably increased, while the maximum time the user has to wait between two queries has decreased significantly (up to 88% in the Job-shop benchmark), due to the overall reduction in the time needed in query generation in almost all problems (as indicated by the $T_{qgen}$ column). As the (maximum) waiting time for the user is of paramount importance for interactive settings, we can see that GROWACQ improves this aspect of time performance of interactive CA systems.

**[Q3] Guided query generation.** In order to evaluate the performance of our proposed objective function for guiding query generation, we compare it with the use of the most popular objective function used in state-of-the-art CA systems, i.e., maximizing violations of constraints from $B$. The objective functions are utilized in line 9 of Algorithm 2. For this comparison, GROWACQ is used, again with MQUACQ-2 as the inner acquisition algorithm at
<table><caption>Table 3 Evaluation of GROWACQ and the proposed approach for guiding query generation.</caption><thead><tr><th>Problem</th><th>#q</th><th>T&#773;<sub>gen</sub></th><th>T&#773;</th><th>T<sub>max</sub></th><th>T<sub>total</sub></th></tr></thead><tbody><tr><td colspan="6">MQUACQ-2</td></tr><tr><td>JSudoku</td><td>6 458</td><td>0.77</td><td>0.10</td><td>3.19</td><td>666.91</td></tr><tr><td>Murder</td><td>370</td><td>0.66</td><td>0.03</td><td>1.10</td><td>12.28</td></tr><tr><td>Random</td><td>5 708</td><td>0.60</td><td>0.04</td><td>2.22</td><td>233.62</td></tr><tr><td>Golomb</td><td>233</td><td>0.96</td><td>0.22</td><td>1.22</td><td>50.89</td></tr><tr><td>Job-shop</td><td>590</td><td>1.03</td><td>0.10</td><td>5.36</td><td>56.23</td></tr><tr><td colspan="6">GROWACQ + MQUACQ-2</td></tr><tr><td>Sudoku</td><td>5 863</td><td>0.15</td><td>0.12</td><td>1.98</td><td>721.15</td></tr><tr><td>Murder</td><td>357</td><td>0.04</td><td>0.02</td><td>0.11</td><td>7.97</td></tr><tr><td>Random</td><td>4 804</td><td>0.14</td><td>0.05</td><td>1.30</td><td>230.36</td></tr><tr><td>Golomb</td><td>270</td><td>0.80</td><td>0.29</td><td>1.30</td><td>78.47</td></tr><tr><td>Job-shop</td><td>786</td><td>0.13</td><td>0.06</td><td>0.66</td><td>48.18</td></tr><tr><td colspan="6">GROWACQ + MQUACQ-2 <sub>guided</sub></td></tr><tr><td>JSudoku</td><td>3 963</td><td>0.15</td><td>0.24</td><td>1.96</td><td>963.42</td></tr><tr><td>Murder</td><td>250</td><td>0.04</td><td>0.04</td><td>0.27</td><td>9.07</td></tr><tr><td>Random</td><td>4 820</td><td>0.14</td><td>0.05</td><td>1.19</td><td>229.47</td></tr><tr><td>Golomb</td><td>100</td><td>0.16</td><td>0.27</td><td>0.95</td><td>27.44</td></tr><tr><td>Job-shop</td><td>776</td><td>0.13</td><td>0.06</td><td>0.64</td><td>47.53</td></tr></tbody></table>

line 7 of Algorithm 3. The results using the guided query generation can be seen in Table 3, bottom-two blocks, comparing GROWACQ + MQUACQ-2 against GROWACQ + MQUACQ-2 guided.

We can see that, when using our probability-based guidance for query generation, the number of queries has significantly decreased in JSudoku, Murder, and Golomb, while it has remained nearly the same in Random and Job-Shop. In the latter cases, the number of queries has not decreased because these are under-constrained problems, and thus the probability derived from the constraints' relations was small. This led to maximizing the violations of all constraints in *B* (i.e., the same behavior as with the existing objective). On the other hand, in the problems that do not have a sparse constraint network, where using the simple counting method to compute the probabilities of the constraints could effectively guide the acquisition system, the decrease observed in the number of queries is substantial (32% in JSudoku, 30% in Murder, and 64% in Golomb). However, as violating constraints one-by-one leads to more queries *generated* at line 3 of Algorithm 1, yet fewer queries at lines 8-9, which are very fast to compute, there is a small increase in the total time on JSudoku.

**[Q4] Combination of our methods.** Comparing the combination of our methods (i.e., GROWACQ + MQUACQ-2 guided) with MQUACQ-2 (Table 3), we can see that combining our bottom-up approach with guiding the query generation greatly outperforms MQUACQ-2 in terms of the number of queries needed to achieve convergence on most of the benchmarks. The number of queries has decreased on all benchmarks except Job-shop, where, because of its sparse target network, we need 23% more queries, as GROWACQ increases the number of queries to converge in under-constrained problems, due to the reasons described in section 6.4, while guiding the query generation does not improve it, as the probabilities estimated are always low. In the rest of the problems, we observe a total decrease of 16% in Random, 39% in JSudoku, 32% in Murder, and up to 60% in Golomb.
These results demonstrate the effectiveness of the proposed methods in reducing the number of queries needed for CA algorithms, which is crucial in interactive scenarios.

## 6.5 [Q5] Dealing with larger biases

To answer this question, we evaluated GROWACQ and the combination of our methods on larger instances of the Job-shop benchmark, using the same language as before. We used two instances: one with 15 jobs, 11 machines, and 40 steps (denoted as JS-15-11), which resulted in a bias consisting of 542 850 constraints, and one with 19 jobs, 12 machines, and again 40 steps (denoted as JS-19-12), resulting in a bias of 1 037 400 constraints. The results are presented in Table 4.

On the one hand, GROWACQ needs more queries to converge (like on the smaller Job-Shop instance) because the constraint network of this problem is sparse. Yet the total time needed to converge is one order of magnitude lower than in MQUACQ-2, being 24.4 times faster in the instance with a bias size of 0.5 million constraints and 25.6 times faster in the instance with $|B| > 1M$. In addition, the maximum waiting time has drastically decreased by using GROWACQ (and the combination GROWACQ and guiding query generation), from 5 499 seconds to only 3 (resp. 8) seconds in JS-15-11 and from more than 20 371 seconds to only 7 (resp. 6) seconds in JS-19-12. Importantly, the average waiting time is more than 30 times lower when using GROWACQ. Note that, as in the smaller job-shop instance, guiding does not lead to improvement in terms of the number of queries. However, it does not noticeably worsen the time performance of the system.

Hence, the experiments confirm that the proposed methodology can efficiently handle significantly larger sets of candidate constraints than the state of the art, up to 50 times larger than the ones commonly used in the literature [6, 8, 27, 29].

<table><caption>**Table 4** Experimental results on instances with a large bias.</caption><thead><tr><th>Problem</th><th>|B|</th><th>#q</th><th>T&#773;<sub>qgen</sub></th><th>T&#773;</th><th>T<sub>max</sub></th><th>T<sub>total</sub></th></tr></thead><tbody><tr><td colspan="7" style="text-align:center; font-weight:bold;">MQUACQ-2</td></tr><tr><td>JS-15-11</td><td>≈ 0.5M</td><td>5 456</td><td>66.12</td><td>6.25</td><td>5 499.76</td><td>34 085.73</td></tr><tr><td>JS-19-12</td><td>≈ 1M</td><td>8 012</td><td>80.99</td><td>9.75</td><td>20 371.74</td><td>78 124.41</td></tr><tr><td colspan="7" style="text-align:center; font-weight:bold;">GROWACQ + MQUACQ-2</td></tr><tr><td>JS-15-11</td><td>≈ 0.5M</td><td>7 015</td><td>0.44</td><td>0.20</td><td>2.84</td><td>1 422.93</td></tr><tr><td>JS-19-12</td><td>≈ 1M</td><td>10 309</td><td>0.62</td><td>0.29</td><td>6.92</td><td>2 984.77</td></tr><tr><td colspan="7" style="text-align:center; font-weight:bold;">GROWACQ + MQUACQ-2<sub>guided</sub></td></tr><tr><td>JS-15-11</td><td>≈ 0.5M</td><td>7 062</td><td>0.44</td><td>0.20</td><td>7.88</td><td>1 399.63</td></tr><tr><td>JS-19-12</td><td>≈ 1M</td><td>10 219</td><td>0.64</td><td>0.30</td><td>6.19</td><td>3 054.68</td></tr></tbody></table>

## 7 Conclusions

Some of the most important limitations of interactive CA methods are the large number of queries needed to converge, as well as the size of the candidate constraint set that they can handle efficiently. In this work, we presented novel methods to alleviate these issues, improving the efficiency of CA systems. We proposed a bottom-up approach, which allows the system to handle significantly larger biases, reducing the maximum waiting time for the user, and also reducing the total number of queries needed when the target constraint network is not sparse. We also introduced a probabilistic method to guide query generation,
further reducing the number of posted queries when our simple counting method could guide the acquisition system to learn constraints more efficiently. In addition, we presented a new query generation technique, named PQ-GEN, that allows the use of conventional CP solvers, removing the dependency of existing methods on customized solvers to converge. Our experimental evaluation showed that our proposed methods outperform state-of-the-art systems in terms of the number of queries in problems with non-sparse constraint networks, reducing this number up to 60%. In addition, the experiments show that GROWACQ can handle up to 50 times larger biases than the ones commonly used in the literature, allowing CA to tackle increasingly large and complex problems. The biggest avenue for future work is to further investigate additional ways to reduce the number of queries needed, e.g., by using guidance in all parts of the acquisition process (not just the query generation), and with more advanced probabilistic models. Another important avenue is to consider the setting in which user answers can be noisy as has been investigated for passive systems.

References

1. Hajar Ait Addi, Christian Bessiere, Redouane Ezzahir, and Nadjib Lazaar. Time-bounded query generator for constraint acquisition. In *International Conference on the Integration of Constraint Programming, Artificial Intelligence, and Operations Research*, pages 1–17. Springer, 2018.

2. Dana Angluin. Queries and concept learning. *Machine learning*, 2(4):319–342, 1988.

3. Robin Arcangioli, Christian Bessiere, and Nadjib Lazaar. Multiple constraint acquisition. In *IJCAI: International Joint Conference on Artificial Intelligence*, pages 698–704, 2016.

4. Nicolas Beldiceanu and Helmut Simonis. A model seeker: Extracting global constraint models from positive examples. In *Principles and practice of constraint programming*, pages 141–157. Springer, 2012.

5. Senne Berden, Mohit Kumar, Samuel Kolb, and Tias Guns. Learning max-sat models from examples using genetic algorithms and knowledge compilation. In *28th International Conference on Principles and Practice of Constraint Programming (CP 2022)*, 2022.

6. Christian Bessiere, Clement Carbonnel, Anton Dries, Emmanuel Hebrard, George Katsirelos, Nina Narodytska, Claude-Guy Quimper, Kostas Stergiou, Dimosthenis C Tsouros, and Toby Walsh. Learning constraints through partial queries. *Artificial Intelligence*, 319:103896, 2023.

7. Christian Bessiere, Remi Coletta, Eugene C Freuder, and Barry O'Sullivan. Leveraging the learning power of examples in automated constraint acquisition. In *International Conference on Principles and Practice of Constraint Programming*, pages 123–137. Springer, 2004.

8. Christian Bessiere, Remi Coletta, Emmanuel Hebrard, George Katsirelos, Nadjib Lazaar, Nina Narodytska, Claude-Guy Quimper, Toby Walsh, et al. Constraint acquisition via partial queries. In *IJCAI*, volume 13, pages 475–481, 2013.

9. Christian Bessiere, Remi Coletta, Frédéric Koriche, and Barry O'Sullivan. A sat-based version space algorithm for acquiring constraint satisfaction problems. In *European Conference on Machine Learning*, pages 23–34. Springer, 2005.

10. Christian Bessiere, Remi Coletta, Barry O'Sullivan, Mathias Paulin, et al. Query-driven constraint acquisition. In *IJCAI*, volume 7, pages 50–55, 2007.

11. Christian Bessiere, Frédéric Koriche, Nadjib Lazaar, and Barry O'Sullivan. Constraint acquisition. *Artificial Intelligence*, 244:315–342, 2017.

12. Eugene C Freuder. Modeling: the final frontier. In *The First International Conference on The Practical Application of Constraint Technologies and Logic Programming (PACLP)*, London, pages 15–21, 1999.

13. Eugene C Freuder. Progress towards the holy grail. *Constraints*, 23(2):158–171, 2018.

14. Eugene C Freuder and Barry O'Sullivan. Grand challenges for constraint programming. *Constraints*, 19(2):150–162, 2014.
15 Eugene C Freuder and Richard J Wallace. Suggestion strategies for constraint-based match-maker agents. In *International Conference on Principles and Practice of Constraint Programming*, pages 192–204. Springer, 1998.

16 Tias Guns. Increasing modeling language convenience with a universal n-dimensional array, cppy as python-embedded example. In *Proceedings of the 18th workshop on Constraint Modelling and Reformulation at CP (Modref 2019)*, volume 19, 2019.

17 Mohit Kumar et al. Acquiring integer programs from data. In *Proceedings of the Twenty-Eighth International Joint Conference on Artificial Intelligence*. IJCAI, 2019.

18 Mohit Kumar, Samuel Kolb, and Tias Guns. Learning constraint programming models from data using generate-and-aggregate. In *28th International Conference on Principles and Practice of Constraint Programming (CP 2022)*. Schloss Dagstuhl-Leibniz-Zentrum für Informatik, 2022.

19 Arnaud Lallouet, Matthieu Lopez, Lionel Martin, and Christel Vrain. On learning constraint problems. In *Tools with Artificial Intelligence (ICTAI), 2010 22nd IEEE International Conference on*, volume 1, pages 45–52. IEEE, 2010.

20 Nadjib Lazaar. Parallel constraint acquisition. In *Proceedings of the AAAI Conference on Artificial Intelligence*, volume 35, pages 3860–3867, 2021.

21 Michele Lombardi, Michela Milano, and Andrea Bartolini. Empirical decision model learning. *Artificial Intelligence*, 244:343–367, 2017.

22 Steven D Prestwich. Robust constraint acquisition by sequential analysis. *Frontiers in Artificial Intelligence and Applications*, 325:355–362, 2020.

23 Steven D Prestwich, Eugene C Freuder, Barry O'Sullivan, and David Browne. Classifier-based constraint acquisition. *Annals of Mathematics and Artificial Intelligence*, pages 1–20, 2021.

24 Kostyantyn Shchekotykhin and Gerhard Friedrich. Argumentation based constraint acquisition. In *Data Mining, 2009. ICDM'09. Ninth IEEE International Conference on*, pages 476–482. IEEE, 2009.

25 Dimosthenis C Tsouros and Kostas Stergiou. Efficient multiple constraint acquisition. *Constraints*, 25(3):180–225, 2020.

26 Dimosthenis C Tsouros and Kostas Stergiou. Learning max-csps via active constraint acquisition. In *27th International Conference on Principles and Practice of Constraint Programming (CP 2021)*. Schloss Dagstuhl-Leibniz-Zentrum für Informatik, 2021.

27 Dimosthenis C Tsouros, Kostas Stergiou, and Christian Bessiere. Structure-driven multiple constraint acquisition. In *International Conference on Principles and Practice of Constraint Programming*, pages 709–725. Springer, 2019.

28 Dimosthenis C Tsouros, Kostas Stergiou, and Christian Bessiere. Omissions in constraint acquisition. In *International Conference on Principles and Practice of Constraint Programming*, pages 935–951. Springer, 2020.

29 Dimosthenis C. Tsouros, Kostas Stergiou, and Panagiotis G. Sarigiannidis. Efficient methods for constraint acquisition. In *24th International Conference on Principles and Practice of Constraint Programming*, 2018.
A Hyperparameter evaluation for PQ-Gen and TQ-Gen

<figure><img src="image_336.png" /><figcaption>(a)</figcaption></figure>

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

Figure 1 Performance of PQ-GEN with different hyperparameters values, in terms of: a) the number of queries posted and b) the time (s) needed (in brackets we show the maximum waiting time for the user).

<figure><img src="image_338.png" /><figcaption>(a)</figcaption></figure>

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

Figure 2 Performance of TQ-GEN with different parameters, in terms of: a) the convergence rate (in brackets we show the number of queries posted when it converged) and b) the time (s) needed (in brackets we show the maximum waiting time for the user).

Both PQ-GEN, our projection-based query generation approach, and TQ-GEN [1] (discussed in Section 3.1) involve hyperparameters that affect their performance. As mentioned in Section 6.3, we performed a sensitivity analysis of the performance with respect to the hyperparameter configuration used of PQ-GEN and TQ-GEN [1]. In this comparison, both query generation methods were used within the state-of-the-art active CA method MQUACQ-2. We used the JSudoku benchmark for this comparison, as from the benchmarks considered in this paper, this is shown to be the hardest one to reach convergence on (see Table 2).

In more detail, we varied the hyperparameters of both PQ-GEN and TQ-GEN to assess their performance under different configurations. While we fixed the hyperparameter $\alpha$ of [1] to 0.8 as recommended in [1], we had to try different values for the time-related hyperparameters, $\tau$ and $t$. We did not use the values proposed by the authors in [1] because we use a different solver and system, and this can affect significantly the time performance.
In our evaluation we used $\tau = [0.05s, 0.1s, 0.2s, 0.3s]$ and $t = [0.5s, 1s, 1.5s, 2s]$ for TQ-GEN. We also used the adjust function described in [1], as it has been shown to improve its performance. For PQ-GEN hyperparameters, we used $l = \{3000, 5000, 7500, 10000\}$ and $t = \{0.5, 1, 1.5, 2\}$. Thus, we examined 16 different configurations for each. The results of our experiments are presented in Figures 1 and 2, respectively, for PQ-GEN and TQ-GEN.

Focusing on Figure 1, we can see that the performance of PQ-GEN is stable across all configurations, both in terms of the number of queries and time performance, having also converged in all cases. Let us now shift our focus to Figure 2 and the performance of TQ-GEN. The first observation is that in the majority of the cases, MQUACQ-2 failed to converge when using TQ-GEN as the query generator. Only when the time limit was set to 2s, we see at least one run achieving convergence for all values of $\tau$. In addition, the performance of MQUACQ-2 using TQ-GEN is highly sensitive to changes in hyperparameter values, particularly with respect to time.

Overall, comparing the results of PQ-GEN and TQ-GEN, we observe that PQ-GEN exhibits superior performance in terms of convergence rate, fully overcoming the issue of premature convergence. PQ-GEN also requires a lower number of queries to reach convergence and offers improved time performance, resulting in reduced waiting times for the user.
# Addressing Problem Drift in UNHCR Fund Allocation

Sameela Suharshani Wijesundara

Department of Data Science and AI, Faculty of IT, Monash University, Clayton, Australia
ARC Industrial Training and Transformation Centre OPTIMA, Clayton, Australia

Maria Garcia de la Banda

Department of Data Science and AI, Faculty of IT, Monash University, Clayton, Australia
ARC Industrial Training and Transformation Centre OPTIMA, Clayton, Australia

Guido Tack

Department of Data Science and AI, Faculty of IT, Monash University, Clayton, Australia
ARC Industrial Training and Transformation Centre OPTIMA, Clayton, Australia

## Abstract

Optimisation models are concise mathematical representations of real-world problems, usually developed by modelling experts in consultation with domain experts. Typically, domain experts are only indirectly involved in the problem modelling process, providing information and feedback, and thus perceive the deployed model as a black box. Unfortunately, real-world problems “drift” over time, where changes in the input data parameters and/or requirements cause the developed model to fail. This requires modelling experts to revisit and update deployed models. This paper identifies the issue of problem drift in optimisation problems using as case study a model we developed for the United Nations High Commissioner for Refugees (UNHCR) to help them allocate funds to different crises. We describe the initial model and the challenges due to problem drift that occurred over the following years. We then use this case study to explore techniques for mitigating problem drift by including domain experts in the modelling process via techniques such as domain specific languages.

**2012 ACM Subject Classification** Theory of computation → Constraint and logic programming; Software and its engineering → Constraint and logic languages

**Keywords and phrases** Fund Allocation, Problem Drift, Domain Specific Languages, MiniZinc

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

## 1 Introduction

Most software engineering projects suffer from one or several forms of *drift*, where a change in requirements during the lifetime of the project causes an existing, deployed software to degrade, to fail, or in the worst case, to produce incorrect solutions. Software projects that include optimisation components are no exception. In this paper, we present a real-world problem as a case study in *problem drift*, and discuss how we can mitigate its effects by enabling the problem owners to become part of the modelling and maintenance workflow.

Our case study is the United Nations High Commissioner for Refugees (UNHCR) fund allocation. UNHCR works in 135 countries and territories to provide crucial assistance to the millions of people forced to flee as a result of conflicts, crises, persecution or natural disasters; recently estimated as 112.6 million [19]. To achieve this, every year UNHCR proposes an annual *budget* for each of the many *projects* it wants to tackle and releases a Global Appeal requesting donations to cover it. For 2023 alone the Global Appeal was more than US$10 billion [18]. The large amount of donations, or *funding*, received as a result need to be distributed among the budgeted projects according to different priorities, while ensuring the allocation respects any constraints set by the donors. For example, it is common for funding to be “earmarked” for a specific crisis, or a geographical area such as a region or a country. These allocation constraints often make it difficult for the funds to be optimally allocated.

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

© Sameela Suharshani Wijesundara, Maria Garcia de la Banda, and Guido Tack;
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. 37; pp. 37:1–37:18

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

LIPICS

Leibniz International Proceedings in Informatics

Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl Publishing, Germany
As a result, excess funds would be carried forward from one year to the next rather than
being directed to UNHCR's activities such as life-saving assistance. The aim of UNHCR is
to maximise the funds available for activities.

This *fund allocation problem* can be represented as a combinatorial optimisation problem. As such, it can be modelled in a high-level constraint modelling language (e.g., [5, 7, 14, 22]) and solved with state-of-the-art optimisation solving technology (e.g., [1, 17, 23]). In 2016, we developed a model for this problem that would have been able to provide US$400 million more in allocations than their previous method, leaving only US$100 million in unallocated funds. As it is common when developing such models, UNHCR experts were only indirectly involved in the modelling process: they provided information and feedback but did not fully understand the implemented model as they lacked the required expertise. This is unfortunate since, as it often happens with real-world problems, the problem requirements *drifted* over time, experiencing changes in the format and kind of the input data, as well as in the rules for allocating funds. After a while, the drift grew to the point where the developed model became unusable. In addition, personnel changes led to a loss of organisational memory, severing the contact between the organisation and us, the developers. As a result, UNHCR never implemented a full system using the model and continued to use their original greedy allocation algorithm, which they were able to maintain and update using in-house expertise.

This paper makes **three contributions**. The first is identifying the issue of problem drift in optimisation models and highlighting its importance to the community as an area worth investigating. The second contribution is our modelling of the UNHCR fund allocation problem: an initial model developed in the first phase of our work in 2016, and a new model developed in a second phase to incorporate the changes that caused the problem to drift. Our third and final contribution results from an analysis of this problem drift, and a subsequent exploration of techniques to mitigate this issue. We present an approach to include domain experts in the modelling process by adapting known software engineering techniques such as domain specific languages. This phase of our work is, so far, exploratory. However, we believe our approach can increase the domain experts' trust in the developed solution, by making them an integral part of the development process, and can also enable the domain experts to adapt the model over time to keep up with problem drift.

## 2 Background

Combinatorial optimisation problems require finding a combination of choices (i.e., a **solution**) that satisfies a set of requirements and optimises the quality of the solutions. Modern approaches for solving these problems first specify a *model* of the problem that formally describes its choices, in terms of **variables** representing the decisions and their **domains** representing the choices available; its requirements, in terms of **constraints** defined over the variables; and the quality of the solutions in terms of an **objective function** also defined over the variables. Models are often specified using **parameters** so that input data can be provided independently of the model. This allows models to be used to solve any **instance** of the problem by instantiating the model parameters with the concrete input data and translating these instances into input for a **solver**, which then produces solutions to the problem. At a high level, this modelling and solving approach can be depicted as the step-wise process shown in Figure 1.

- *Step 1* formalises the real world problem via a model. This is currently a manual process where modelling experts communicate with domain experts to capture all relevant requirements, often requiring several iterations to progressively improve the formalisation (*refinement process*). The result is a model of the problem requirements written in a high-level modelling language such as MiniZinc [14], Essence [7], OPL [22] or AMPL [5].
<figure><img src="image_342.png" /><figcaption>Figure 1 Step-wise illustration of the combinatorial problem modelling and solving process.</figcaption></figure>

- *Step 2* takes the model along with concrete input data formatted according to the model's parameters and representing an instance of the problem, and compiles this instance into a format suitable for the solver selected by the modeller. Each of the modelling languages mentioned in *Step 1* has its own compilation methods and lower level languages such as FlatZinc for MiniZinc [14] and Essence' for Essence [7].

- *Step 3* solves the compiled instance using the selected solver. Solvers use state-of-the-art algorithms and heuristics that are efficient in tackling the unique challenges inherent in the kind of constraints they support. There are many different solving technology paradigms including MIP [23] (for Mixed Integer Programming), CP [17] (for Constraint Programming), and SAT [1] (for Boolean Satisfiability Problems).

*Steps 2 and 3* are extensively automated thanks to many years of research on modelling language and solver design. However, *Step 1* is a manual process and, in practice, considerably more complex than depicted in Figure 1, involving multiple refinement and feedback loops.

# 3 Phase 1: The 2016 UNHCR Fund Allocation Problem and Model

## 3.1 Problem Description

In 2016 we worked with UNHCR on an optimisation model for their fund allocation. According to the rules at the time, part of or all of the donations pledged by an organisation (a fund) could be either *earmarked* or not, indicating whether the donor put some restrictions on the kind of projects the fund could cover or not, respectively. The kind of restrictions supported by the algorithm used at the time by UNHCR were based on the following fund and project attributes:

- **Level**: can take one of three values: *Region*, *Sub-region*, or *Country*, indicating the geographic scale of a project, e.g., a continent such as Europe, a part of a given continent such as Eastern Europe, or a country such as Ukraine. For a fund, the level is used to restrict its allocation, and can have an additional value *Global* to indicate that the fund is not restricted to a geographical area.

- **Region**: the region of the world where the project is located or to which the fund is allocated. Only relevant when the level is *Region*; disregarded otherwise.

- **Sub-region**: the sub-region of the world where the project is located or to which the fund is allocated. Only relevant when the level is *Sub-region*; disregarded otherwise.
37:4 Addressing Problem Drift in UNHCR Fund Allocation

■ Area of Budgetary Control (ABC): the name of the UNHCR office which could be
responsible either for a single country or a collection of countries each with fewer activities.

- Pillar: the broad target area of support. **Pillar 1** stands for refugees, **Pillar 2** for stateless people, **Pillar 3** for reintegration, and **Pillar 4** for internally displaced people. Pillars can be combined, e.g., **Pillars 1-2**. A fund that can be spent in any pillar is given the value **All Pillars**.

- Situation: the particular issue a project is tackling or a fund can be spent on, e.g., **South Sudan Situation**. A fund that can be spent on any situation is given the value **Country/Regular Program**.

The input data is given by an Excel sheet where each row is referred to as an *item*. The first six columns correspond to an item's pillar, situation, level, region, sub-region, and ABC. The next two columns contain the *budget* and *income*. An item with an empty budget column is a fund, and with a non-empty budget column is a project. Note that the budget of some projects can already be partially covered by money from other sources (e.g., bank interest). In such cases both the budget and income columns will have non-zero amounts.

The objective is to move money from fund items to project items to minimise the total amount of money left unspent in the funds, while ensuring the allocation constraints are satisfied. This means, ensuring money does not move from item A to item B if they have:

1. Different pillars, unless the pillar of item A is the generic All Pillars.

2. Different regions if the level of item A is Region.

3. Different sub-regions if the level of item A is Sub-region.

4. Different ABCs if the level of item A is Country.

5. Different situations, unless the situation of item A is Country/Regular Programme, its level is Country, its pillar is All Pillars and item A is in surplus for its country.

The first four preclude a fund earmarked for a non-generic pillar, region, sub-region, or country to be spent on a different one. The last one ensures that funds for generic country situations are spent first on their country. There are also common-sense rules ensuring, for example, that funds are not overspent and projects do not get more income than budgeted.

3.2 Model

Modelling this problem in MiniZinc was easy. We will not reproduce the full model here, but
we will give enough details for later discussions.

**Parameters.** The Excel input was fed into 9 MiniZinc parameters: integer *n* representing the number of rows and used to build set *ITEM* of 1..*n*, and eight one-dimensional arrays indexed by *ITEM* representing the values of the eight columns for each *ITEM*, i.e., its pillar, situation, level, region, sub-region, ABC, budget and income. The following MiniZinc code shows the definition of some of these parameters together with two enumerated data types (*ensems*) whose definitions (strings from the input data) are the values for *REGION* and *SUBREGION*, an auxiliary parameter *maxinc* set to the maximum income, and an array *excess*, where *excess[i]* is the initial amount of money either available in fund *i*, or needed by project *i*:

set of int: ITEM = 1..n;
enum REGION;
enum SUBREGION;
array[ITEM] of REGION: region;
array[ITEM] of SUBREGION: subregion;
array[ITEM] of int:income;
array[ITEM] of int:budget;
int: maxinc = max(income);
array[ITEM] of int: excess = [ income[i] - budget[i] | i in ITEM ];

While `ensems` were not used in the model, they are more intuitive and used here for brevity.
**Variables.** The decision variables are stored in a two-dimensional array `movement`, where `movement[i,j]` represents the amount of money moved from item `i` to item `j`, and is constrained to be between 0 and `maxinc`. Note that if `i` is a project, the value of `movement[i,j]` will be later set to 0. This array is defined as follows:

```javascript
array[ITEM,ITEM] of var 0..maxinc: movement;
```

In addition, several auxiliary variables are defined based on `movement`. For example, `out[i]` and `iin[i]` represent, respectively, the amount of money moved out of item `i` or into it. They are defined as follows:

```javascript
array[ITEM] of var 0..maxinc: out =
  [ sum(j in ITEM)(movement[i,j]) | i in ITEM ];
array[ITEM] of var 0..maxinc: iin =
  [ sum(j in ITEM)(movement[j,i]) | i in ITEM ];
```

**Constraints.** Are quite simple and modelled using the above parameters and variables. For example, the common-sense constraint ensuring funds are not overspent is modelled as:

```javascript
constraintforall(i in ITEM where excess[i] > 0)(out[i] <= excess[i]);
constraintforall(i in ITEM where excess[i] <= 0)
  (forall(j in ITEM)(movement[i,j] = 0);
```

ensuring that if item `i` is a fund (`excess[i]>0`), the amount of money that moves out of `i` is not higher than its excess, and otherwise no money moves out of `i`. Similarly, constraint 2 which correctly earmarks funds based on region is modelled as follows:

```javascript
constraintforall(i,j in ITEM
  where region[i] != region[j] /\ level[i] = Region)
  (movement[i,j] = 0);
```

ensuring no money moves from `i` to `j` if level of `i` is `Region` and their regions do not match. Note that constraints 3 and 4 follow exactly the same structure.

Our last example is constraint 5; the most complex. First we define auxiliary parameters `surplus_to_country fuels` as either a surplus if `i` is a generic country fund (i.e., one whose pillar is `all_pillars`, level `Country` and situation `country fuels`), or 0, otherwise. The surplus is defined as the fund's income minus the money needed by that country. We also define Boolean parameter `has_sufficiency` to true iff the associated surplus parameter is positive.

```javascript
array[ITEM] of -maxinc..maxinc: surplus_to_country fuels =
  [ if situation[i] != country fuels
    \ level[i] != Country
    \ pillar[i] != all_pillars
  then 0
  else income[i] -
    sum(j in ITEM where abc[i] = abc[j] /
        situation[j] = country fuels)
    (budget[j] - income[j])
  elif
  | i in ITEM ];
array[ITEM] of bool: has_sufficiency_to_country fuels =
  [ surplus_to_country fuels[i] > 0 | i in ITEM ];
```

These auxiliary parameters are then used to model part of constraint 5 as follows:
37:6 Addressing Problem Drift in UNHCR Fund Allocation

```pascal
constraintforall(i in ITEM
  where situation[i] = countryжи regular programme
  /\ level[i] = Country
  /\ pillar[i] = all_pillars
  /\ has_surplus_to_country府 regular_situation[i])
  (sum(j in ITEM where situation[i] != situation[j]) (movement[i,j])
  <= surplus_to_country府 regular_situation[i]);
```

If there is no surplus, the other part (not shown) sets `movement[i,j]=0` for every j in the same country but different situation.

**Objective.** Minimises the excess of funds and deficit of projects after allocation by means of the intermediate variable *difference*, which captures the excess for fund items (positive number) and the deficit for projects (negative number). The sum *obj* of the absolute value *absdiff* of these values is then computed and minimised as follows.

```pascal
array[ITEM] of var -maxinc .. maxinc: difference =
  [ excess[i] - out[i] + inn[i] | i in ITEM ];
array[ITEM] of var 0..maxinc: absdiff;
```

```pascal
constraintforall(i in ITEM)(absdiff[i] >= difference[i] /\
  absdiff[i] >= -difference[i]);
```

```pascal
var int: obj = sum(i in ITEM)(absdiff[i]);
solve minimize obj;
```

This model solves the above fund allocation problem in less than 10s for the sample data we were given. With the 2016 sample data, the model was able to allocate US$400 million more than the UNHCR greedy algorithm used at the time. UNHCR domain experts saw great value in the results produced by the optimisation model, and used it occasionally to provide supplementary input into their decision making process.

# 4 Phase 2: The 2022 UNHCR Fund Allocation Problem and Model

## 4.1 Problem Changes

When we contacted UNHCR again several years later, we learned that problem drift had made our model obsolete. After a number of discussions that spanned several months, we identified the following main changes to the problem.

**Data format.** Instead of the data being provided as a table where each row was an item representing both funds and projects, it is now provided as two separate tables (that we will refer to as *entities*), one for funds and one for projects.

**Attributes.** Some attributes changed names (e.g., fund income was now called *Balance*), others disappeared (e.g., funds no longer had budget as attribute, since they were always empty), and new ones appeared. The complete list of attributes for funds is: *Unique ID*, *Level*, *Domain*, *Region*, *Sub-Region*, *Country*, *Cost Center*, *Pillar*, *Situation*, *Account Status*, *Account Type*, *Balance*. Projects have the same attributes except that instead of *Balance*, they have a *Requirement* attribute that is equivalent to previously subtracting income from budget. Also, while previously the smallest geographical operational area provided in the sample data was the value of attribute ABC, countries are now further divided into smaller geographical areas defined by the attribute *Cost Centres*. More importantly, funds and
projects are now given priority levels via two different attributes. One is *Account Type*, which can take values *OL* (*operational level*) and *AOL* (*above operational level*), the latter having lesser priority [20]. The other is *Domain* which can take several values such as, in decreasing priority order, Headquarters (for leadership, management, policy guidance, etc), Global (for global programmes undertaken at the headquarters but of direct benefit to field operations), and Field (for field operations, which still takes 88% of the funding according to the 2022 Global Report [19]).

**Constraints.** Are now defined in terms of these new attributes. In particular, the constraints for Cost Centre and Domain are similar to those for other attributes.

**Objective.** While the main objective remained unchanged, it now also prioritised the allocation of funds based on their Account Type and Domain, as well as funds carried over from previous years and those with the highest number of constraints.

## 4.2 Model Changes

**Parameters.** Changes to the parameters of the model were mainly caused by UNHCR's decision to split *ITEM* into two entities, *FUNDS* and *PROJECTS*, and change their associated attributes. We modelled this similarly to before but using a different set of indexed arrays per entity, as follows:

```pascal
set of int: PROJECTS = 1..numProjects;
set of int: FUNDS = 1..numFunds;

array[PROJECTS] of LEVEL: p_level;
array[PROJECTS] of DOMAIN: p_domain;
array[PROJECTS] of REGION: p_region;
array[PROJECTS] of ACCOUNT_TYPE: p_account_type;
array[PROJECTS] of int: p_requirement;

array[FUNDS] of LEVEL: f_level;
array[FUNDS] of DOMAIN: f_domain;
array[FUNDS] of REGION: f_region;
array[FUNDS] of ACCOUNT_TYPE: f_account_type;
array[FUNDS] of int: f_balance;
```

**Variables.** The main decision variables are almost identical to those in the first model, except for the index sets now being known to be funds for the first dimension and projects for the second, and the maximum value being represented by fund parameter *f_balance*.

```pascal
array[FUNDS,PROJECTS] of var 0..max(f_balance) : movement;
```

Most auxiliary variables also undergo very superficial changes, such as `out[i]` and `iin[i]`:

```pascal
array[FUNDS] of var int : out =
  [ sum(p in PROJECTS)(movement[f,p]) | f in FUNDS ];
array[PROJECTS] of int : iin =
  [ sum(f in FUNDS)(movement[f,p]) | p in PROJECTS ];
```

**Constraints.** Some constraints encoding the same rules only require relatively superficial changes. For example, the constraint that earmarks funds based on region looks very similar:
37:8 Addressing Problem Drift in UNHCR Fund Allocation

```pascal
constraintforall(f in FUNDS, p in PROJECTS
  where f_region[f] != p_region[p] /\ f_level[f] == L3)
  (movement[f,p] = 0);
```

if one knows L3 is equivalent to the old Region. Others change more, e.g., ensuring funds are
not overspent is much simpler now that funds and projects are clearly separated:

```pascal
constraintforall(f in FUNDS)(out[f] <= f_balance[f]);
```

In addition, we need new constraints to correctly earmark funds based on new attributes,
such as cost centre, which look very similar to those implemented for Phase 1 attributes:

```pascal
constraintforall(f in FUNDS, p in PROJECTS
  where f_cost_centre[f] != p_cost_centre[p] /\ f_level[f] == L6)
  (movement[f,p] = 0);
```

**Objective.** As described above, the old objective maximising the overall fund allocation was now extended to prioritise carry-over funds from the previous round, tightly earmarked funds (funds with the tightest restrictions), and funds allocated according to their Domain and Account Type. The new objective function was implemented as a weighted sum of these different aspects in order to achieve the required prioritisation.

4.3 Implementation and Adoption

We implemented the changes to the model as described above over the course of several months,
in consultation with the UNHCR domain expert. Compared to the in-house algorithm, the
new model is faster (it solves the sample data in under 3 minutes) and more effective in
terms of total allocation and priority order. We were also able to experiment with extensions
of the model, such as achieving fairness and balance in the allocation, which is easier to do
using optimisation technology compared to greedy algorithms.

Unfortunately, the new model was not adopted by UNHCR. However, as a result of
our work with them, they proposed the use of optimisation techniques with the new re-
implementation of their enterprise resource planning (ERP) system. The viability of this is
currently under investigation.

5 Phase 3: Addressing Problem Drift via Interactive Modelling DSLs

The transition from Phase 1 to Phase 2 of the UNHCR project encountered the usual
challenges that exist for optimisation models even after they are successfully integrated
into real-world workflows, which is an already challenging task. First, the problem had
drifted between the two phases and the end users, i.e., the domain experts, did not have
sufficient in-house expertise to adjust the implementation to tackle the drift. As a result,
they continued to use their original method which was less effective but much more familiar.
And second, due to subsequent personnel changes in UNHCR, Phase 2 required not just
updating the implemented model, but on-boarding a new domain expert who had no reason
to be invested in a new method other than our promises of a better future. While we were
lucky enough for our new domain expert to trust us and become invested in developing a
better approach, it is easier for organisations to simply continue using their old systems and
move on.
In this section, we will present our *initial exploration* of how domain experts could be integrated better into the modelling process. Our proposal is based on our experience with Phases 1 and 2 of the UNHCR problem. Our goal is to develop a new *framework* for integrating domain experts into the modelling process in such a way they can (a) understand and verify how their user requirements have been mapped to the model, thus increasing the model's *transparency*; and (b) modify the model to cope with some level of problem drift, thus increasing the model's *flexibility*. Additionally, such an approach may help organisations obtain a better optimisation model, achieve higher levels of trust in the model, and develop more in-house expertise to manage the life cycle of the model. The rest of this section discusses this framework.

## 5.1 Related Work

The field of *Interactive Optimisation* has a long tradition of introducing domain experts into the development process early on. Meignan et al. [12] present a comprehensive overview of interactive optimisation approaches from operations research studies where domain experts are involved via a User Interface (UI) to guide the *solving* process. It also presents the high level components of an interactive optimisation system that includes a *preference model*, that is, an intermediate model between the UI and the optimisation model that captures user preferences when guiding the solver towards the preferred solution. This approach can be useful for users to find their preferred solutions, particularly for multi-objective optimisation problems with a large set of solutions. However, it can be achieved without domain experts understanding much about the model or being able to modify anything except its search and/or objective function.

Liu et al. [10] introduce a theoretical framework for interactive optimisation problem solving, called the *Problem Solving Loop* (left of Figure 2) and based on the *Sense Making Loop* [16] from the field of Visual Analytics [9]. The Problem Solving Loop aims at directly engaging the domain experts in the optimisation process. In doing this, it captures the high-level user goals and tasks in two major sub-loops: the *model-defining loop* and the *optimisation loop*, which correspond to the modelling design phase and the decision making phase, respectively.

<figure><img src="image_343.png" /><figcaption>Figure 2 Problem solving loop presented by Liu et al. [10] and updated Model Defining Loop.</figcaption></figure>
However, Liu et al. focus mostly on the latter, which covers the utilisation of the model by the domain experts from the moment they execute a *Problem Instance* of the *General Optimisation Model*, to the moment they make a final decision based on a *List of Ranked Solutions* (see left of Figure 2). Their model-defining loop does capture the fact that, in practice, problems are not converted to models in a singular step; rather, this involves multiple iterations of discussion among different stakeholders. However, they mostly ignore this sub-loop because they argue that its tasks, which lead to the creation of the *General Optimisation Model* and its *Problem Instance*, are usually done by modelling experts. While this is certainly the case currently, it does not mean it is appropriate, as we discuss below.

## 5.2 An Updated Model-Defining Loop

Our experience in several real-world projects indicates that to increase the model's transparency and flexibility, we must deeply involve domain experts in the model-defining loop. Further, we need to do so by giving them a formal role and a formal language with which to communicate not only with the modelling experts but with the model itself. This will enable domain experts to verify the modelling experts' perception of the problem, as well as to directly refine and modify the model, both as part of its initial definition and its ongoing maintenance. Thus, it should help address problem drift. To this end, we focus on extending this loop with the addition of a *Requirements Model* as depicted on the right of Figure 2. The Requirements Model captures user requirements in a *high-level, expressive* intermediate formal language at the early stages of the project. This language will have to be different for each application area, in order to allow both the modelling experts and the domain experts to express the parts of the model they are responsible for. Thus, we propose for the Requirements Model to be expressed in a *Domain Specific Language (DSL)* [13, 21, 8]. DSLs offer powerful expressivity while being tailored to a specific application domain, and have proven successful in many different application areas [6].

We call the DSLs that are used for the Requirements Model *Modelling DSLs, or MDSLs*. Our proposal is for the Requirements Model, expressed in an application-specific MDSL, to be compiled into a General Optimisation Model, which can then be instantiated and solved in the usual way. It is the responsibility of the MDSL designers to develop a compiler (often referred to as a *generator* in the DSL literature [15]) that generates efficient and correct code for the General Optimisation Model.

An important design decision is the level of abstraction and complexity required of an MDSL. On the one hand, we could make it as expressive as a General Purpose Modelling Language (GPML) such as MiniZinc, but this would defeat the purpose of simplicity and easy communication with the domain experts. On the other hand, significantly restricting its expressivity could severely limit its usefulness. For this reason, we propose to split the General Optimisation Model into two parts. A *fixed part*, which is implemented by the modelling expert; and a *flexible part*, which is generated (compiled) from an MDSL specification.

In the following, we will discuss which parts of an optimisation model such an MDSL should be able to express, and at what level of abstraction the MDSL should be designed. We will then see a concrete example of an MDSL for the UNHCR fund allocation problem.
## 5.3 MDSL expressivity

When designing an MDSL, it is important to select an appropriate level of expressivity, striking a balance between ease of use (especially for non-experts) and usefulness. While there are guidelines for general DSL design [8, 21, 13], we will focus here on the specific aspects related to Modelling DSLs. Let us look at the basic parts of an optimisation model and discuss how those parts are amenable to being expressed within MDSLs.

### 5.3.1 Variables

The choice of variables for an optimisation problem is one of the fundamental modelling decisions. It has a profound impact on the way the constraints are modelled and on the performance of solving algorithms – a fact that domain experts would usually not be aware of. Therefore, while it is important to define variables in a way that is easy to understand by the domain experts, we argue that the choice itself should not be left to domain experts. The modelling experts should define the variables, which are then made available to domain experts to be used (but not modified) via the MDSL.

This results in an architecture where the General Optimisation Model in our updated model-defining loop is split into two parts: a *fixed* part and a *flexible* one. The modelling experts design the former, which includes the choice of decision variables and is expressed in a General Purpose Modelling Language. This fixed part of the model cannot be modified through the MDSL, and any updates require input from a modelling expert. In addition, the modelling and domain experts collaborate on the *flexible model part*, which captures the Requirements Model and is defined using an MDSL. The two parts together are compiled into the General Optimisation Model, which can then be instantiated with data and solved.

Figure 3 shows a more detailed view of the new model-defining loop that includes this split into a fixed core and a flexible MDSL model.

### 5.3.2 Parameters

Parameters capture the core objects and data that the model is concerned with. We can split parameters into two kinds: *entities* and *attributes*. *Entities* define the main objects in a problem. For example, in Phase 1 of the UNHCR model, the main entity was an *ITEM*, while Phase 2 had two entities, *PROJECT* and *FUND*. *Attributes* define the *properties* of the entities. For example, in our Phase 1 model the *ITEM* entity had attributes such as *income* and *budget*, while in Phase 2 we saw *requirement* for *PROJECT* and *balance* for *FUND*.

Entities and their attributes commonly appear in the definitions of the constraints and/or the objective function in the model. For example, the balance of a fund will appear in the rules that govern how that fund can be used. The transition from Phase 1 to Phase 2 in the UNHCR problem saw the addition of new attributes (such as cost centre) and the removal of others (such as budget), which required changes to the corresponding funding rules. An MDSL can accommodate attribute changes by adding support for defining attributes and then using them consistently in the constraints and objective function. To do this we can build on the numerous DSLs commonly used in practice for this purpose, such as UML diagramming languages [2], class diagrams used in the object oriented paradigm and entity-relationship diagramming languages [3] used in database design. Note that this part of the MDSL is not application specific and can be used across many domains.

Accommodating changes in entities is more complex, as user decisions often relate to entities. For example, in the UNHCR problem we had to decide how much funding to move from each *fund* to each *project*. Thus, entity changes often require changes in the variables.
Based on our rationale above for leaving the choice of variables to the modelling experts, we do not propose yet to add support into MDSLs for defining new entities. However, a more abstract version of entity might be needed to avoid simple changes in input format.

### 5.3.3 Constraints

Supporting a close collaboration between modelling and domain experts during the definition and refinement of constraints is key for ensuring transparency. Supporting domain experts in modifying them is key for flexibility. As a result, any MDSL will need to be designed such that constraints can be specified and modified in a way that is natural for the domain experts. The Requirements Model, expressed in the MDSL, will serve as a means of communication and documentation between the modelling and domain experts. Therefore, the main challenge when applying our approach is designing an application-specific MDSL with sufficient expressivity to capture current and likely future constraints, yet sufficient simplicity to be understandable by domain experts. A prototype MDSL for the UNHCR problem is presented in Section 5.4.

### 5.3.4 Objective function

Changes to the objective function are also common, both during the modelling process and during the life cycle of a model [11]. To support this, MDSLs can provide basic arithmetic expressions involving all model variables and parameters that can then be used to define and/or modify the objective. This however is not enough for defining extra variables that may be needed to modify the objective. Since we do not want domain experts to introduce new variables using the MDSL, modelling experts will have to foresee and pre-define any auxiliary variables that may be exposed via the MDSL. For example, if the auxiliary variables for `total_allocation` and `minimum_allocation` are pre-modelled, domain experts will subsequently be able to formulate new objective functions using them.

The combination of basic arithmetic expressions and pre-defined variables also allows MDSLs to support domain experts in experimenting with multi-objective solving by combining different objectives in a weighted sum. These weights can then be changed by domain experts (including setting them to 0 to remove certain terms from the objective). Domain experts can even add new terms into the objective if required.

### 5.3.5 Proposed framework

Given the above discussion, we designed our framework to expand the functionality of the updated model-defining loop suggested in Section 5.2 with a new workflow depicted in Figure 3. As discussed earlier, the Requirements Model is used to *extract and simplify* the change requests and to capture them using an MDSL. Domain experts use it to *understand* how constraints are mapped to the model and to build *trust* in the model. The Requirements Model is then *compiled* to generate GPML code to be added to the *General Optimisation Model*. An interesting extension would be to support bidirectional compilation, where GPML code can be reflected back into the Requirements Model. This would in turn allow the model to be *explained* in terms of the MDSL in the Requirements Model.

### 5.3.6 Natural vs Formal Languages for MDSLs

The main function of the MDSL is to serve as the interface between the domain expert and the optimisation model. Therefore, it has to strike a balance between being easy to understand by the domain expert, easy to implement, and easy to compile into a general
optimisation model. Being easy to understand suggests designing a language closer to the
natural language a domain expert would use to describe the problem in a conversation; being
easy to implement and to compile suggests one closer to the mathematical model that is
going to be the result of the compilation process. Unfortunately, the vague and ambiguous
nature of natural language can be problematic for an optimisation model, as its aim is to
capture a problem specification precisely and unambiguously.

Let us illustrate this with an example from our UNHCR model. The following are
statements from our domain experts, communicating to us their rules for transferring money:

1. Can transfer money from a fund to a project if the pillar of the fund is “all pillars” or the pillars of the fund and the projects match

2. Can transfer money from a fund to a project if the situation of the fund is “regular situation” or the situation of the fund and the projects match

3. Can transfer money from a fund to a project if the level of fund is equal to 2 and the domains of the fund and project match

4. Can transfer money from a fund to a project if the level of the fund is equal to 3 and the regions of the fund and project match

These statements have several issues. First, while they are expressed as *Can transfer*, expressing “possibility” in an optimisation model is not easy. For example, our model’s *movement[i,j]* variable represents whether there is a transfer (*movement[i,j] > 0*) or not (*movement[i,j]=0*). We can represent this rule as *If movement[i,j] > 0, then the pillar of the fund must be “all pillars”, or the pillars of the fund and the projects must match*. While this captures the correct meaning, it may not be immediately intuitive for the domain expert (or anyone else). An alternative would be to state the contrapositive: *If the pillar of the fund is not “all pillars” and the pillars of the fund and the projects do not match, then we must not transfer funds, so movement[i,j]=0*. This may be more intuitive, but it requires negating all conditions in the rules, which might be difficult to implement.

The second issue is the ambiguity regarding whether the statements are logically connected
by “and” or by “or”, or if they are supposed to be mutually exclusive. This is where the
statements in contra-positive form have an advantage, as each of them stands on its own and
it is clearer they should be connected by an “and”. However, even then we still have a third
source of ambiguity, as we do not know whether the statements are exhaustive or not. In

<figure><img src="image_344.png" /><figcaption>**Figure 3** Further detailed model-defining loop incorporating MDSLs.</figcaption></figure>
other words: if none of the rules holds, should we set `movement[i,j] = 0` or not? While it turns out the answer is “no”, the opposite could have been true, which is more difficult to achieve in an optimisation model.

A final point to note is that rules 3 and 4, when interpreted as mentioned above, in fact contradict each other: Rule 3 would mean that if the level is not 2, or the domains do not match, no movement is possible, while rule 4 states the same for level 3 and the regions. Clearly, the intended interpretation here is that if `movement[i,j] > 0` and the level is 2, then the domains must match.

These examples make it clear that, while a language closer to natural language seems convenient for domain experts, it is not precise enough to capture a unique logical meaning. As capturing the constraints in a mathematically precise format is crucial for optimisation modelling, care must be taken for the MDSLs to strike the right balance between being easy to understand and mathematically precise. Given the narrow scope of the language, domain experts will require some training to use an MDSL. However, we anticipate this to be only a minor hurdle when compared to learning a complete modelling language such as MiniZinc.

With the current fast advance of AI-based large-language-model discourse systems, it would be interesting to explore how this technology could be used to bridge the gap between precise mathematical MDSL and natural language. However, it is unclear how the inherent ambiguity in natural language formulations could be avoided.

## 5.4 An MDSL for UNHCR fund allocation problem

Our design goal for the MDSL grammar of the new UNHCR problem was to support domain experts in formalising any new required changes to the model in terms of the given set of parameters and variables. Compared to the full MiniZinc language, the MDSL grammar should be significantly simpler, covering only a small sub-set of MiniZinc's capabilities specific to this problem domain. Most of the constraints in the UNHCR problem are expressed in the form of ad-hoc *rules* that consist of *conditions* and *consequences* in the form of *if-then-else* statements. This level of modelling is close to the actual constraints, while still being relatively easy to understand even for someone without any programming experience.

We first define the *abstract* syntax for the MDSL in the form of the following grammar:

```
<rule>     ::= <logical-expr> | <iteration> | <if-then-else>
<logical-expr> ::= <expr> <compare-op> <expr>
                  | "not" <logical-expr>
                  | <logical-expr> "and" <logical-expr>
                  | <logical-expr> "or" <logical-expr>
<iteration> ::= "forall" (<identifier>, ...) <logical-expr>
<if-then-else> ::= "if" <logical-expr> "then" <logical-expr> ("else" <logical-expr>)
<compare-op>  ::= "=" | "!=" | "<" | "<=" | ">" | ">="
<expr>       ::= <named-expr> | <number> | <expr> <arithmetic-op> <expr>
<named-expr> ::= <identifier> | <identifier> "of" (<identifier>, ...)
                  | <identifier> "from" <identifier> "to" <identifier>
<arithmetic-op> ::= "+" | "-" | "*" | "/"
```

The resulting abstract grammar is much simpler than that of MiniZinc. In particular, if-then-else expressions can only have Boolean type and cannot be nested; the grammar only supports one kind of iteration (`forall`); and it restricts the way variables and parameters are addressed. These restrictions were identified as sufficient for UNHCR problem.
While the grammar defines the kind of constraints that can be expressed in the MDSL, it does not define the concrete representation of the constraints the domain experts would use. Thus, the grammar can be implemented in several different ways. We now show a concrete text-based language as well as a simple graphical interface below to illustrate the possibilities. However, we do not claim that those are the most intuitive or best suited for this purpose.

A concrete text-based representation of the above abstract MDSL grammar could use syntax that is close to a more natural way of stating constraints, while still being concise and unambiguous. The following example shows a rule as defined in MiniZinc and in a concrete version of the abstract MDSL grammar:

### MiniZinc constraint

```c
constraintforall(f in FUNDS, p in PROJECTS where
    f_level[f] == L4 /\ f_subregion[f] != p_subregion[p])
    (movement[f,p] = 0);
```

### Constraint represented using a concrete version of the abstract MDSL grammar

```c
For Every FUND, PROJECT :
    IF [(level OF FUND )== (L4)) and
        ((subregion OF FUND )!= (subregion OF PROJECT ))]
    THEN
        [(movement FROM FUND TO PROJECT )== (0)]
```

Note how the concrete MDSL grammar simplifies the access to entity attributes, by using notation such as `level OF FUND` instead of having to introduce a new identifier `f` and corresponding expressions `f_level[f]`, as it is done in the MiniZinc model. For two-dimensional arrays like `movement`, it defines an explicit syntax `movement FROM FUND TO PROJECT` that avoids the ambiguity of the direction of movement present in the MiniZinc model. Many other concrete versions of the abstract MSDSL grammar are also possible. For example, it would be easy to make the concrete grammar more verbose by writing the iteration as `For Every FUND and PROJECT`. Importantly, the MDSL grammar does not have hard-wired

<figure><img src="image_345.png" /><figcaption>Figure 4 User Interface for forming intermediate formalisation.</figcaption></figure>
attribute names. If the set of attributes of an entity is later extended, the MDSL can easily
be extend accordingly. Another important choice is that of operator precedences. We opted
for full explicit bracketing, in order to avoid ambiguity. However, there are many other
approaches, such as indentation-based grouping of expressions that belong together. It is
beyond the scope of this paper to explore the advantages and disadvantages of different
choices in concrete syntax.

In addition to the text-based representation, a graphical representation of the MDSL specifications could be used. Figure 4 shows a prototype browser-based interface that is still quite close to the text version. It has two main sections: The *Rule Base* section displays the already defined rules, and allows users to turn them on and off as needed. The *Rule Constructor* section allows users to define new rules using drop down list boxes containing various components of the language grammar, such as parameters (entity attributes) and variables. A graphical representation like this (or an integration of the text-based language into a graphical development environment) may have several advantages over pure text-based languages. For instance, precedence between different expressions can be indicated graphically, which may help avoid mistakes in the specification. Furthermore, the graphical interface can restrict the choices to those that are valid in a particular situation.

6 Conclusion and Future Work

Many organisations need to solve optimisation problems as part of their core operation, but they still perceive the introduction of technologies like model-based optimisation in their decision support systems as a significant risk. They may appreciate that this technology can solve difficult problems, that it can produce better outcomes than simple algorithms or manual solutions, and that it can incorporate additional factors such as uncertainty, fairness or diversity of solutions. However, the lack of in-house expertise for implementation and maintenance and the perceived black-box nature of the technology are often hurdles that are difficult to overcome. The change of input data formats or requirements after an optimisation solution has been deployed, which we refer to as *problem drift*, adds further risks to this process.

This paper presents the UNHCR fund allocation problem, a combinatorial optimisation problem that deals with the distribution of donated funds to humanitarian projects according to ad-hoc rules set by the donors and operational requirements of UNHCR. While the initial model we implemented in 2016 worked well for the sample data, the problem drift that occurred over the next few years prevented the organisation from fully adopting it, instead continuing to use their in-house algorithm. Phase 2 of our project focused on adapting the initial optimisation model to the changed requirements.

Based on this experience, this paper identifies problem drift in optimisation models and highlights its importance to the community as an area worth investigating. It then proposes to tackle it via an extended framework for model development that incorporates a *Requirements Model* into the Model Defining Loop defined by Liu et al. [10]. The requirements model is based on an application specific Modelling Domain Specific Language (MDSL). We have explored guidelines for the trade-off between general expressivity and ease of use of MDSLs. Our recommendations include that certain parts of a model such as the choice of decision variables need to be left to the optimisation expert, while other parts, such as entity attributes, constraints and objective function, can be exposed in the MDSL. We believe that the introduction of a Requirements Model and MDSL into the modelling process can facilitate collaboration between modelling and domain experts, achieve higher levels of trust with the end users of optimisation technology, as well as equip them with the skills and tools to address certain forms of problem drift in-house. sp An important area for future
work is to evaluate the feasibility, suitability and effectiveness of our proposed framework for addressing problem drift in real-world applications. We acknowledge that the development of an application-specific MDSL together with a model may seem prohibitively expensive. Future work will therefore include exploring the use of template MDSLs for larger problem classes to enable code reuse across applications. Furthermore, we will investigate how to speed up and streamline the MDSL development process using tools such as Language Workbenches [6, 4]. We will also study graphical and hybrid text/graphical MDSLs, which can draw on techniques from research areas such as interactive optimisation and visual analytics. Finally, it would be interesting to explore the use of AI-based dialogue systems such as large language models as a user-facing representation of an MDSL, where the potential ambiguity of natural language becomes a challenging research topic.

References

1 A. Biere, M. Heule, and H. van Maaren. *Handbook of Satisfiability*. IOS Press, January 2009.

2 Grady Booch, James E. Rumbaugh, and Ivar Jacobson. *The unified modeling language user guide - covers UML 2.0, Second Edition*. Addison Wesley object technology series. Addison-Wesley, 2005.

3 Peter P. Chen. The entity-relationship model - toward a unified view of data. *ACM Trans. Database Syst.*, 1(1):9–36, 1976. doi:10.1145/320434.320440.

4 Sebastian Erdweg, Tijs van der Storm, Markus Völter, Meinte Boersma, Remi Bosman, William R. Cook, Albert Gerritsen, Angelo Hulshout, Steven Kelly, Alex Loh, Gabriël D. P. Konat, Pedro J. Molina, Martin Palatnik, Risto Pohjonen, Eugen Schindler, Klemens Schindler, Riccardo Solmi, Vlad A. Vergu, Eelco Visser, Kevin van der Vlist, Guido H. Wachsmuth, and Jimi van der Woning. The State of the Art in Language Workbenches. In David Hutchison, Takeo Kanade, Josef Kittler, Jon M. Kleinberg, Friedemann Mattern, John C. Mitchell, Moni Naor, Oscar Nierstrasz, C. Pandu Rangan, Bernhard Steffen, Madhu Sudan, Demetri Terzopoulos, Doug Tygar, Moshe Y. Vardi, Gerhard Weikum, Martin Erwig, Richard F. Paige, and Eric Van Wyk, editors, *Software Language Engineering*, volume 8225, pages 197–217. Springer International Publishing, Cham, 2013. Series Title: Lecture Notes in Computer Science. doi:10.1007/978-3-319-02654-1_11.

5 Robert Fourer, David M Gay, and Brian W Kernighan. A modeling language for mathematical programming. *Management Science*, 36(5):519–554, 1990.

6 Martin Fowler. *Domain-Specific Languages*. Pearson Education, September 2010.

7 Alan M. Frisch, Warwick Harvey, Chris Jefferson, Bernadette Martínez-Hernández, and Ian Miguel. Essence: A constraint language for specifying combinatorial problems. *Constraints*, 13(3):268–306, September 2008. doi:10.1007/s10601-008-9047-y.

8 Paul Hudak. Domain Specific Languages. *Handbook of programming languages*, 3(39-60):23, 1997.

9 Youn-ah Kang and John Stasko. Examining the Use of a Visual Analytics System for Sensemaking Tasks: Case Studies with Domain Experts. *IEEE Transactions on Visualization and Computer Graphics*, 18(12):2869–2878, December 2012. doi:10.1109/TVCG.2012.224.

10 Jie Liu, Tim Dwyer, Kim Marriott, Jeremy Millar, and Annette Haworth. Understanding the Relationship Between Interactive Optimisation and Visual Analytics in the Context of Prostate Brachytherapy. *IEEE Transactions on Visualization and Computer Graphics*, 24(1):319–329, January 2018. doi:10.1109/TVCG.2017.2744418.

11 Jie Liu, Tim Dwyer, Guido Tack, Samuel Gratzl, and Kim Marriott. Supporting the Problem-Solving Loop: Designing Highly Interactive Optimisation Systems. *IEEE Transactions on Visualization and Computer Graphics*, 27(2):1764–1774, February 2021. doi:10.1109/TVCG.2020.3030364.
37:18 Addressing Problem Drift in UNHCR Fund Allocation

12 David Meignan, Sigrid Knust, Jean-Marc Frayret, Gilles Pesant, and Nicolas Gaud. A Review and Taxonomy of Interactive Optimization Methods in Operations Research. *ACM Transactions on Interactive Intelligent Systems*, 5(3):1–43, October 2015. doi:10.1145/2808234.

13 Marjan Mernik, Jan Heering, and Anthony M Sloane. When and how to develop domain-specific languages. *ACM computing surveys (CSUR)*, 37(4):316–344, 2005.

14 Nicholas Nethercote, Peter J. Stuckey, Ralph Becket, Sebastian Brand, Gregory J. Duck, and Guido Tack. MiniZinc: Towards a Standard CP Modelling Language. In Christian Bessière, editor, *Principles and Practice of Constraint Programming – CP 2007*, Lecture Notes in Computer Science, pages 529–543, Berlin, Heidelberg, 2007. Springer. doi:10.1007/978-3-540-74970-7_38.

15 Václav Pech.JetBrains MPS: Why Modern Language Workbenches Matter. In Antonio Bucchiarone, Antonio Cicchetti, Federico Ciccozzi, and Alfonso Pierantonio, editors, *Domain- Specific Languages in Practice: with睛MPS*, pages 1–22. Springer International Publishing, Cham, 2021. doi:10.1007/978-3-030-73758-0_1.

16 P. Pirolli and S. Card. The sensemaking process and leverage points for analyst technology as identified through cognitive task analysis. In *Proceedings of international conference on intelligence analysis*, volume 5. McLean, VA, USA, 2005.

17 Francesca Rossi, Peter van Beek, and Toby Walsh. *Handbook of Constraint Programming*. Elsevier, August 2006.

18 UNHCR Global Appeal 2023. Accessed on 5 May, 2023. URL: https://reporting.unhcr.org/globalappeal2023.

19 UNHCR Global Report 2022. Accessed on 5 May, 2023. URL: https://reporting.unhcr.org/global-report-2022.

20 UNHCR operations plan in emergencies. Accessed on 5 May, 2023. URL: https://emergency.unhcr.org/support-response/planning-and-programming/unhcr-operations-plan-emergencies.

21 Arie van Deursen, Paul Klint, and Joost Visser. Domain-specific languages: an annotated bibliography. *ACM SIGPLAN Notices*, 35(6):26–36, June 2000. doi:10.1145/352029.352035.

22 Pascal Van Hentenryck. *The OPL optimization programming language*. MIT Press, Cambridge, MA, USA, 1999.

23 Laurence A. Wolsey. *Integer Programming*. John Wiley & Sons, September 2020.
# From Formal Boosted Tree Explanations to Interpretable Rule Sets

Jinqiang Yu

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

Department of Data Science and AI, Monash University, Clayton, Australia

Australian Research Council OPTIMA ITTC, Clayton, Australia

Alexey Ignatiev

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

Department of Data Science and AI, Monash University, Clayton, Australia

Peter J. Stuckey

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

Department of Data Science and AI, Monash University, Clayton, Australia

Australian Research Council OPTIMA ITTC, Clayton, Australia

## Abstract

The rapid rise of Artificial Intelligence (AI) and Machine Learning (ML) has invoked the need for *explainable AI* (XAI). One of the most prominent approaches to XAI is to train rule-based ML models, e.g. decision trees, lists and sets, that are deemed interpretable due to their transparent nature. Recent years have witnessed a large body of work in the area of constraints- and reasoning-based approaches to the inference of interpretable models, in particular decision sets (DSes). Despite being shown to outperform heuristic approaches in terms of accuracy, most of them suffer from scalability issues and often fail to handle large training data, in which case no solution is offered. Motivated by this limitation and the success of gradient boosted trees, we propose a novel anytime approach to producing DSes that are both accurate and interpretable. The approach makes use of the concept of a generalized formal explanation and builds on the recent advances in formal explainability of gradient boosted trees. Experimental results obtained on a wide range of datasets, demonstrate that our approach produces DSes that more accurate than those of the state-of-the-art algorithms and comparable with them in terms of explanation size.

**2012 ACM Subject Classification** Computing methodologies → Machine learning

**Keywords and phrases** Decision set, interpretable model, gradient boosted tree, BT compilation

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

**Supplementary Material** Software (Source Code): https://github.com/jinqiang-yu/cpl/
archived at swh:1:dir:40fde451a732a518f78caa9ad372a7c267446836

**Funding** This research was partially funded by the Australian Government through the Australian Research Council Industrial Transformation Training Centre in Optimisation Technologies, Integrated Methodologies, and Applications (OPTIMA), Project ID IC200100009.

## 1 Introduction

Rapid development of Artificial Intelligence (AI) and Machine Learning (ML) have revolu-
tionized all aspects of human lives in recent years [30, 1]. However, decisions made by most
widely used ML models are hard for humans to understand hence the interest in the theory
and practice of *Explainable AI* (XAI) rises.

One major approach to XAI is to compute *post-hoc* explanations for ML predictions
to answer a “why” question [34, 44], i.e. why the prediction is made. Although heuristic
approaches to post-hoc explanations prevail [34, 44, 43], they suffer from a number of
weaknesses [21, 16, 49, 52]. Formal methods [48, 20, 37] provide alternative approaches
to explanations that avoid these weaknesses. Another alternative approach to XAI is to
compute *interpretable* ML models, i.e. logic-based models, including decision trees [40],

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

© Jinqiang Yu, Alexey Ignatiev, and Peter J. Stuckey;

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. 38; pp. 38:1–38:21

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

LIPICS

Leibniz International Proceedings in Informatics

Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl Publishing, Germany
38:2 From Formal Boosted Tree Explanations to Interpretable Rule Sets

decision lists [46], and decision sets [29]. These models enable decision makers to obtain
succinct explanations from the models directly. In this paper, we focus on the decision
set (DS) models.

Decisions sets are particularly easy to explain: the rule that fired is an explanation of the decision. This led to an upsurge in interest of decision sets that are both interpretable and accurate. Recent work [50] uses propositional satisfiability (SAT) to generate minimum-size decision sets that are perfectly accurate on the training data, and demonstrates that decision sets that completely agree with the training data outperform others in terms of accuracy. A more scalable maximum satisfiability (MaxSAT) approach [18] to this problem was then proposed. Unfortunately, both of these methods are unable to provide any decision information if a dataset is not completely solved.

Motivated by these works and their limitations, this paper aims at making a bridge between formal post-hoc explainability and interpretable DS models. In particular, the paper focuses on developing a novel anytime approach to computing decision sets that are both interpretable and accurate, by compiling a gradient boosted tree model into a decision set on demand with the use of formal explanations. This is done with the use of the recent approach [17] to compute abductive explanations for gradient boosted trees using maximum satisfiability (MaxSAT). Furthermore, the paper proposes a range of post-hoc model reduction heuristics aiming at enhancing interpretability of the result models, done with MaxSAT and integer linear programming (ILP). The experimental results show that compared with other state-of-the-art methods, decision sets generated by the proposed approach are more accurate, and comparable with the competition in terms of interpretability.

2 Preliminaries

**SAT and MaxSAT.** The standard definitions for propositional satisfiability (SAT) and maximum satisfiability (MaxSAT) solving are assumed [3]. A propositional formula $\phi$ is said to be in conjunctive normal form (CNF) if it is a conjunction of clauses. A clause is a disjunction of literals, where a literal is either a Boolean variable $b$ or its negation $\neg b$. A truth assignment $\mu$ is a mapping from the set of variables to $\{0,1\}$. A clause is said to be satisfied by truth assignment $\mu$ if one of the literals in the clause is assigned value 1; otherwise, the clause is falsified. If all clauses in formula $\phi$ are satisfied by assignment $\mu$, $\phi$ is satisfied; otherwise, assignment $\mu$ falsifies $\phi$. A CNF formula $\phi$ is unsatisfiable if there exists no assignment satisfying $\phi$.

In the context of unsatisfiable formulas, the MaxSAT problem consists in finding a truth
assignment that maximizes the number of satisfied clauses. Hereinafter, we use a variant
of MaxSAT called Partial Weighted MaxSAT [3, Chapters 23 and 24]. The formula $\phi$ in
this variant is represented as a conjunction of *hard* clauses $\mathcal{H}$, which must be satisfied, and
*soft* clauses $\mathcal{S}$ where each of them is associated with a weight representing a preference to
satisfy them, i.e. $\phi = \mathcal{H} \wedge \mathcal{S}$. Partial Weighted MaxSAT problems aim at finding a truth
assignment $\mu$ that satisfies all hard clauses and maximizes the total weight of satisfied soft
clauses.

**Classification Problems.** We consider classification problems with a set of classes<fnref n="1" /> $\mathcal{K} = \{1, \dots, k\}$, and a set of features $\mathcal{F} = \{1, \dots, m\}$. The value of each feature $i \in \mathcal{F}$ is taken from its corresponding (numeric) domain $D_i$. As a result, the entire feature space is defined as

<footnote n="1">¹ Non-integer class labels can be mapped to a set {1, ..., |K|}.</footnote>
<table><caption>(a) Decision set.</caption><tr><td>IF "petal.length" &lt; 2.60</td><td>THEN class = "setosa"</td></tr><tr><td>IF 2.60 ≤ "petal.length" &lt; 4.95 ∧ "petal.width" &lt; 1.75</td><td>THEN class = "versicolor"</td></tr><tr><td>IF "petal.length" ≥ 2.60 ∧ "petal.width" ≥ 1.75</td><td>THEN class = "virginica"</td></tr><tr><td>IF "petal.length" ≥ 4.95</td><td>THEN class = "virginica"</td></tr></table>

<figure><img src="image_351.png" /><figcaption>(b) BT model [5] consisting of 2 trees per class, each of depth $\le 2$, adopted from [17].</figcaption></figure>

**Figure 1** Example DS and BT models computed on the well-known Iris classification dataset.

$\mathbb{F} \triangleq \prod_{i=1}^{m} D_i$. A concrete point represented by $\mathbf{v} = (v_1, \dots, v_m) \in \mathbb{F}$, s.t. each $v_i$ is a constant value taken by feature $i \in \mathcal{F}$, together with its corresponding class $c \in \mathcal{K}$, represented by a pair $(\mathbf{v}, c)$, indicate a *data instance* or *example*. With a slight abuse of notation and whenever convenient, a data point $\mathbf{v} \in \mathbb{F}$ is also referred to as an instance. Finally, $\mathbf{x} = (x_1, \dots, x_m)$ denotes a vector of feature variables $x_i \in D_i$, $i \in \mathcal{F}$, used for reasoning over points in $\mathbb{F}$.

A classifier defines a *classification function* $\tau: \mathbb{F} \to \mathcal{K}$. The objective of classification problems is to learn a function $\tau$ to generalize well on unseen data given a training dataset $\mathcal{E} = \{e_1, e_2, \dots, e_n\}$, where each instance $e_d \in \mathcal{E}$ is a pair of $(\mathbf{v}_d, c_d)$. Classification problems are conventionally posed as an optimization problem, i.e. either to minimize the complexity of $\tau$, or maximize its accuracy, or both.

**Rules, Decision Sets and Gradient Boosted Trees.** Multiple ways exist to learn classifiers given data $\mathcal{E}$. This paper focuses on arguably one of the most interpretable models, i.e. decision sets, trained by *compiling* gradient boosted trees.

A *decision rule* is in the form of “IF antecedent THEN prediction”, where the antecedent is a set of feature literals. Informally, a rule is said to classify an instance $\mathbf{v} \in \mathbb{F}$ as class $c \in \mathcal{K}$ if its antecedent is *compatible* with $\mathbf{v}$ (or *matches* $\mathbf{v}$) and its prediction is $c$. A *decision set* (DS) is an unordered set of decision rules $\mathcal{R}$. An instance $(\mathbf{v}, c) \in \mathcal{E}$ is misclassified by a DS if either there exists no rule in $\mathcal{R}$ matching $\mathbf{v}$, or there exists a rule classifying $\mathbf{v}$ as a class $c' \in \mathcal{K}$ s.t. $c' \neq c$.

A *gradient boosted tree* (BT) is a tree ensemble $\mathcal{T}$ defining sets of decision trees $T_c \in \mathcal{T}$ for each class $c \in [|\mathcal{K}|]$, where $T_c$ comprises $N \in \mathbb{N}_{>0}$ trees $t_{kz+c}$, $z \in \{0, \dots, N-1\}$, $k = |\mathcal{K}|$. Given an instance $\mathbf{v} \in \mathbb{F}$, its class is obtained by computing the sum of scores assigned by trees for each class $w(\mathbf{v}, c) = \sum_{t \in T_c} t(\mathbf{v})$ and assigning the class which has the maximum score, i.e. $\argmax_{c \in [|\mathcal{K}|]} w(\mathbf{v}, c)$. Whenever convenient, $\mathbf{n} \in t$ denotes a non-terminal node, where $t \in \mathcal{T}$ represents an arbitrary decision tree. Moreover, each such $\mathbf{n}$ indicates a feature condition in the form of $x_i < d$, where feature $i \in \mathcal{F}$ and *splitting threshold* $d \in D_i$.
<table><caption>Table 1 Several instances extracted from Iris dataset.</caption><thead><tr><th>#</th><th>sepal.length</th><th>sepal.width</th><th>petal.length</th><th>petal.width</th><th>class</th></tr></thead><tbody><tr><td>e<sub>1</sub></td><td>5.1</td><td>3.5</td><td>1.4</td><td>0.2</td><td>setosa</td></tr><tr><td>e<sub>2</sub></td><td>7.7</td><td>2.6</td><td>6.9</td><td>2.3</td><td>virginica</td></tr><tr><td>e<sub>3</sub></td><td>5.6</td><td>2.5</td><td>3.9</td><td>1.1</td><td>versicolor</td></tr><tr><td>e<sub>4</sub></td><td>6.2</td><td>2.8</td><td>4.8</td><td>1.8</td><td>virginica</td></tr><tr><td>e<sub>5</sub></td><td>5.6</td><td>2.8</td><td>4.9</td><td>2.0</td><td>virginica</td></tr></tbody></table>

**► Example 1.** Figure 1 shows DS and BT models trained on the Iris dataset, which has 4 numeric features and 3 classes: “setosa”, “versicolor”, and “virginica”. Observe that instance $\mathbf{v}_1 \in e_1$ shown in Table 1 is classified as “setosa” by the first rule of the DS. In the BT model, each class $c \in [3]$ is represented by 2 trees $t_{3z+c}$, $z \in \{0,1\}$. Thus, it also classifies $\mathbf{v}_1$ as “setosa”, since the score of this class $w(\mathbf{v}_1, 1) = t_1 + t_4 = 0.71928$ is higher than the score of “versicolor” $w(\mathbf{v}_1, 2) = t_2 + t_5 = -0.40253$ and the score of “virginica” $w(\mathbf{v}_1, 3) = t_3 + t_6 = -0.41324$. ▼

**Interpretability and Explanations.** Interpretability is not formally defined as it is considered to be a subjective concept [33]. In this paper interpretability is defined as the overall succinctness of the information offered by an ML model to justify a provided prediction. Moreover, following earlier work [48, 20], we equate explanations for ML models with *abductive explanations* (AXps), which are subset-minimal sets of features sufficient to explain a given prediction. Concretely, given an instance $\mathbf{v} \in \mathbb{F}$ and a prediction $c = \tau(\mathbf{v}) \in \mathcal{K}$, an AXp is a subset-minimal set of features $\mathcal{X} \subseteq \mathcal{F}$ such that

$$ \forall (\mathbf{x} \in \mathbb{F}). \left[ \bigwedge_{i \in \mathcal{X}} (x_i = v_i) \right] \rightarrow (\tau(\mathbf{x}) = c) \qquad (1) $$

**► Example 2.** Consider the setup of Example 1. Given instance $\mathbf{v}_1$, observe that for any instance with “petal.length” = 1.4, the BT is guaranteed to predict “setosa” independently of the values of other features, since the weights for “setosa” and “versicolor” are 0.71928 and -0.40253 respectively as before, and the maximal weight for “virginica” is 0.39408 - 0.08968 = 0.30440. Thus, the (only) AXp $\mathcal{X}$ for the prediction for $e_1$ made by the BT model is {“petal.length”}. ▼

**Explanations in BTs.** Formal reasoning has been recently applied to computing AXps for BT models, with the key difficulty being how to effectively reason about the aggregation over a large number of trees in a BT model. Recent work applied satisfiability modulo theory (SMT) [21] or mixed integer linear programming (MILP) solvers [42, 27] to directly address the linear summations arising in the BT encoding. Hereinafter, we build on the recent MaxSAT approach [17], which maps the aggregation reasoning to a set of MaxSAT queries to avoid a costly encoding of the linear constraints into CNF. Also, [17] demonstrates how a MaxSAT query can be made such that (1) holds if and only if the optimal value of the constructed objective function is negative.<fnref n="2" /> In general, assuming that each feature $i \in \mathcal{F}$ is numeric (continuous), the approach orders the set of splitting thresholds $\{d_{i1}, ..., d_{ih_i}\}$ in a BT $\mathcal{T}$ for each feature $i$, where $h_i$ is the total number of thresholds of feature $i$ in $\mathcal{T}$ and $d_{ij} \in \mathcal{D}_i$ for $j \in [h_i]$. Given an instance $\mathbf{v} = (v_1, ..., v_m) \in \mathbb{F}$, the above approach

<footnote n="2">² The reader is referred to [17] for the details.</footnote>
associates each value $v_i$ with a single interval $I'_i$ from the set of disjoint intervals $D_i = \{ I_{i1} \equiv [\min(D_i), d_{i1}), I_{i2} \equiv [d_{i1}, d_{i2}), \dots, I_{ih_i+1} \equiv [d_{ih_i}, \max(D_i)] \}$. Thus, AXp extraction boils down to finding a subset-minimal subset $\mathcal{X} \in \mathcal{F}$ s.t.

$$ \forall (\mathbf{x} \in \mathbb{F}). \left[ \bigwedge_{i \in \mathcal{X}} x_i \in I'_i \right] \rightarrow (\tau(\mathbf{x}) = c) \qquad (2) $$

**► Example 3.** Recall Example 2 and assume “petal.length” and “petal.width” have indices 3 and 4. Note that the sets of splitting thresholds for feature “petal.length” {$d_{31} = 2.60, d_{32} = 4.75, d_{33} = 4.95$} and for feature “petal.width” {$d_{41} = 1.45, d_{42} = 1.65, d_{43} = 1.75$}. Let $min(D_3) = -\infty$ and $min(D_4) = 0.1$. Then we can associate the values of features 3 and 4 in our instance $v_1 \in e_1$ with intervals $I_{31} \equiv (-\infty, 2.60)$ and $I_{41} \equiv [0.1, 1.45)$. Hence by (2), the AXp shown in Example 2 can in fact be seen as a rule $\langle \text{IF } \text{"petal.length"} < 2.60 \text{ THEN } \text{class} = \text{"setosa"} \rangle$. ▼

## 3 Related Work

Interpretable decision sets are logic-based ML models that can be traced back to the 70s and 80s [39, 15, 4, 45]. To the best of our knowledge, [6] proposed the first approach to decision sets, which were introduced as the variant of decision lists [45, 7]. The first method making use of logic and optimization to synthesize a disjunction of rules that match a given dataset was proposed in [26]. Recent work [29] argued that decision sets are more interpretable than the other logic-based models, i.e. decision lists and decision trees. This work uses smooth local search to generate a set of rules first and heuristically minimizes a linear combination of criteria afterwards, e.g. the size of a rule, their maximum number, overlap or error.

Since then a number of works proposed the use of logic reasoning and optimization procedures to train DS models [22, 36, 12, 50, 18] claiming to significantly outperform the approach of [29] in terms of accuracy and performance. Among those, the works closest to ours are [22, 50, 18]. They proposed SAT-based approaches to computing smallest-size decision sets that *perfectly* agree with the training data by minimizing either the number of rules [22, 18] or the number of literals [50, 18] used in the model. Additionally, [50] is capable of computing *sparse* decisions sets that trade off training accuracy for model size. Despite the dramatic performance increase achieved in [18], all the approaches above suffer from scalability issues.

Post-hoc explainability is one of the major approaches to XAI. Besides a plethora of heuristic sampling-based methods to post-hoc explainability [43, 34, 44], a formal reasoning based approach to computing abductive explanations [48, 20] stands out. AXps can be related with prime implicants of the decision function (hence an alternative name *prime implicant explanations*, *PI-explanations*) associated with ML predictions and are guaranteed to capture the semantics of the ML models in the entire feature space. Although hard to compute in general, AXps were shown to be effectively computable for BT models by an incremental MaxSAT-based approach [17].

Our work aims at making a bridge between interpretable DS models and AXp computation by exploiting the latter for training the former. Given a BT model, it focuses on generating decision rules that agree with the BT. Each rule represents an AXp for the prediction made by the BT model, resulting in a DS model in a way *guided* by the original BT model. The approach is shown to outperform the prior logic-based approaches to DS inference in terms of test accuracy and performance. Note that despite prior attempts to train sparse models guided by tree ensembles [38], to our best knowledge, none of the existing works have applied formal post-hoc explanations to compile interpretable models.
Finally, our approach can be related to the existing line of work on *knowledge distillation* [11, 13], where an interpretable model is trained to approximate a hard-to-interpret black-box model, which is often seen as teacher-to-student knowledge transfer. Note that in contrast to knowledge distillation, our approach is able to *compile* a BT into an *equivalent* DS if we consider the entire feature space, as shown below.

# 4 Decision Sets by Boosted Tree Compilation

Based on [17], this section details a MaxSAT-based approach to compiling a BT into a DS where each rule in the DS is equivalent to a prime implicant of the BT classification function.

## 4.1 Rule Extraction

Recall that an AXp, as defined in (1) and (2), can be seen as an *if-then* rule. Given a hard-to-interpret BT model, the AXp extraction approach of [17] can be modified to compute an interpretable DS consisting of a *set* of AXps for the BT. However, when the features are continuous (numeric), this potential approach suffers from the following issue. Recall that an AXp $\mathcal{X} \in \mathcal{F}$ indicates a set of *concrete* feature values that are sufficient to explain a prediction $c = \tau(\mathbf{v})$ for a certain instance $\mathbf{v} \in \mathbb{F}$. Although this same AXp can explain other instances compatible with it, its applicability in general is at the mercy of expressivity of the feature literals used in the AXp, i.e. equality literals and succinct interval membership in the case of (1) and (2), respectively. Motivated by this limitation, we propose to compute AXps over the literals intrinsic to the BT model aiming at getting feature intervals that are as general as possible, as detailed below.<fnref n="3" />

In contrast to the work of [17], which associates each feature value $v_i \in D_i$ with a single *narrowest* interval $I'_i$ covering the value, we exploit all the splitting points used by the BT for feature $i$ and identify all of the corresponding literals satisfied by the feature value $v_i$. Note that the original MaxSAT encoding [17] introduces a single Boolean variable $o_{ij}$ for each literal $x_i < d_{ij}$ with $d_{ij}$ being a j'th threshold used in the BT for feature $i$, s.t. $o_{ij} = 1$ iff $x_i < d_{ij}$ holds true. This way, each positive $o_{ij}$ represents an upper bound on the value of $x_i$ while each negative $\neg o_{ij}$ represents a lower bound on $x_i$.

**► Example 4.** Feature 3 ("petal.length") from Example 3 has 3 thresholds: $d_{31} = 2.60$, $d_{32} = 4.75$, $d_{33} = 4.95$. Boolean variables $o_{31}$, $o_{32}$, and $o_{33}$ are set to true iff $x_3 < 2.60$, $x_3 < 4.75$, and $x_3 < 4.95$, respectively. Let feature 3 take value 3.9 in the instance we want to explain. Observe how we can immediately assign literals $\neg o_{31}$, $o_{32}$, and $o_{33}$ to true. ▼

Next, given an instance $\mathbf{v} = (v_1, \dots, v_m) \in \mathbb{F}$, let us construct a complete conjunction $\bigwedge_{i \in \mathcal{F}, j \in [h_i]} \tilde{o}_{ij}$ of literals $\tilde{o}_{ij}$ s.t. $\tilde{o}_{ij}$ is to be replaced by $o_{ij}$ if $v_i < d_{ij}$ and replaced by $\neg o_{ij}$ otherwise. By construction, this conjunction holds true for instance $\mathbf{v}$. Now, given this conjunction of literals, we can apply the existing approach of [17] to extract a subset-minimal explanation $\mathcal{Y} \subseteq \{\tilde{o}_{ij} \mid i \in \mathcal{F}, j \in [h_i]\}$ for instance $\mathbf{v}$ over literals $\tilde{o}_{ij}$ s.t.

$$ \forall (\mathbf{x} \in \mathbb{F}). \quad [\bigwedge_{l \in \mathcal{Y}} l] \quad \rightarrow (\tau(\mathbf{x}) = c) \qquad (3) $$

Such an explanation $\mathcal{Y}$ may (or may not) define either a lower bound on feature $i$, an upper bound, or both, aiming to construct the *most general* interval for each feature $i \in \mathcal{Y}$. Hence, we informally refer to such explanations as *generalized AXps* or simply *rules* (hereinafter, we use both interchangeably).

<footnote n="3">³ An alternative to our approach is *inflation* of abductive explanations, which is discussed in [23, 24]. Given an AXp, it aims at extending the set of values covered by each feature literal in the AXp while the AXp condition (1) still holds.</footnote>
■ **Algorithm 1** Deletion-based Rule Extraction.

**Function:** RuleExtract($\mathfrak{T}$, $\mathbf{v}$, $c$, $\mathcal{E}$)

**Input:** $\mathfrak{T}$: BT defining $\tau(\mathbf{x})$, $\mathbf{v}$: Instance, $c$: Prediction, i.e. $c = \tau(\mathbf{v})$, $\mathcal{E}$: Training data

**Output:** $\mathcal{Y}$: Subset-minimal rule

1: $\langle \mathcal{H}, \mathcal{S} \rangle \leftarrow \text{Encode}(\mathfrak{T})$

2: $\mathcal{Y} \leftarrow \text{Init}(\mathfrak{T}, \mathbf{v})$

3: $\mathcal{Y} \leftarrow \text{Sort}(\mathcal{Y}, \mathcal{E})$

4: **for** $l \in \mathcal{Y}$ **do**

5: &nbsp;&nbsp;&nbsp;&nbsp;**if** $\text{EntCheck}(\langle \mathcal{H}, \mathcal{S} \rangle, c, \mathcal{Y} \setminus \{l\})$ **then**

6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**Y** $\leftarrow \mathcal{Y} \setminus \{l\}$

7: **return** $\mathcal{Y}$

► **Example 5.** Consider instance $\mathbf{v}_3$ predicted as *"versicolor"* by the BT (observe that $\mathbf{v}_3 = 3.9$ and $\mathbf{v}_4 = 1.1$) and recall the thresholds for features 3 and 4 discussed in Example 3. We can compute a generalized AXp $\mathcal{Y} = \{\neg o_{31}, o_{33}, o_{43}\}$ representing the second rule of the DS shown in Figure 1a. The original approach of [17] would instead compute an AXp defining the narrowest intervals for features 3 and 4, representing a rule: $\langle \text{IF } 2.60 \le \text{"petal.length"} < 4.75 \land \text{"petal.width"} < 1.45 \ \text{THEN } \text{class} = \text{"versicolor"} \rangle$, which is far less general than $\mathcal{Y}$ . ▼

A possible rule extraction procedure is outlined in Algorithm 1. (Please ignore line 3 for now; feature sorting is described in Section 4.2). The input BT model $\mathfrak{T}$ is encoded into MaxSAT by applying the approach of [17]. Given an instance $\mathbf{v} \in \mathbb{F}$, the initial set of literals $\mathcal{Y} = \{\tilde{o}_{ij} | i \in \mathcal{F}, j \in [h_i]\}$ is created. Note that any feature $i \in \mathcal{F}$ unused in the BT $\mathfrak{T}$ is excluded from $\mathcal{Y}$. The rest of the procedure implements the standard deletion-based AXp extraction [20], i.e. it iterates through all literals in $\mathcal{Y}$ one by one, and checks which of the them can be safely removed such that entailment (3) still holds.

► **Example 6.** Consider our running example model and instance $\mathbf{v}_2 \in e_2$ from Table 1 predicted as *"virginica"* by the BT $\mathfrak{T}$. Given the thresholds for features 3 and 4 in Example 3, set $\mathcal{Y}$ is initialized to $\{\neg o_{31}, \neg o_{32}, \neg o_{33}, \neg o_{41}, \neg o_{42}, \neg o_{43}\}$. The other two features are excluded from $\mathcal{Y}$ since they are irrelevant to the classification function in $\mathfrak{T}$. Applying Algorithm 1 results in extracting a subset-minimal generalized AXp $\mathcal{Y} = \{\neg o_{33}\}$, which represents the rule $\langle \text{IF } \text{petal.length} \ge 4.95 \ \text{THEN } \text{class} = \text{"virginica"} \rangle$. ▼

► **Remark 7.** Algorithm 1 relies on deciding whether formula (3) holds for each feature in explanation $\mathcal{Y}$. Here, this is done by means of a series of incremental core-guided MaxSAT oracle calls [19, 17]. One may wonder whether or not incomplete *anytime* MaxSAT solving [31, 35, 2, 32] can be applied in this setting. Although this may look plausible at first glance, time-restricted anytime MaxSAT algorithms can only *over-approximate* exact MaxSAT solutions while (3) holds *if and only if* the exact value of the objective function is negative. Therefore, an over-approximation of a MaxSAT solution is *never able* to prove the validity of (3) and so none of the features being tested can be discarded in the case of incomplete MaxSAT algorithms, which defies the purpose of Algorithm 1.

## 4.2 Boosted Tree Compilation

As mentioned above, generalized AXps can be seen as general decision rules that can be applied to an enormous number of instances. Therefore, it makes little sense to extract such rules for each instance in the feature space $\mathbb{F}$. Instead, one can devise an on-demand
■ **Algorithm 2** Compile a BT into a DS.

**Function:** Compile($\mathfrak{T}$, $\tau$, $\mathcal{C}$)

**Input:** $\mathfrak{T}$: BT defining $\tau(\mathbf{x})$, $\tau$: Classification function in $\mathfrak{T}$, $\mathcal{C}$: Coverage set

**Output:** $\mathcal{R}$: Set of Rules

1: $\mathcal{R} \leftarrow \emptyset$

2: $C_u \leftarrow \mathcal{C}$

3: **while** $C_u \neq \emptyset$ **do**

4: &nbsp;&nbsp;&nbsp;&nbsp;**v** $\leftarrow$ GetInst($C_u$)

5: &nbsp;&nbsp;&nbsp;&nbsp;$\mathcal{Y} \leftarrow$ RuleExtract($\mathfrak{T}$, **v**, $c = \tau(\mathbf{v})$, $C_u$)

6: &nbsp;&nbsp;&nbsp;&nbsp;$C_c \leftarrow$ GetCover($\mathcal{Y}$, $C_u$)

7: &nbsp;&nbsp;&nbsp;&nbsp;$C_u \leftarrow C_u \setminus C_c$

8: &nbsp;&nbsp;&nbsp;&nbsp;$\mathcal{R} \leftarrow \mathcal{R} \cup \mathcal{Y}$

9: **return** $\mathcal{R}$

compilation process, i.e. given a yet uncovered instance $\mathbf{v} \in \mathbb{F}$, we can apply Algorithm 1 to extract a rule covering $\mathbf{v}$ (and some other instances). Clearly, exhaustive compilation of a BT, i.e. if the target is to cover all the instances in $\mathbb{F}$ with generalized AXps of the BT, is computationally expensive given that AXp extraction for tree ensembles is hard for $\mathbb{D}^P$ [25]. This can also lead to the large size of the resulting DSes making them hard to interpret. In practice, local compilation aiming at capturing the behavior of the BT on the training data only, is sufficient to generate a DS, which is both accurate and interpretable.

The proposed approach to compiling a BT $\mathfrak{T}$ into a DS $\mathcal{R}$ is shown in Algorithm 2. We initialize the set $C_u$ of currently uncovered instances to be equal to $\mathcal{C}$, i.e. the set of examples we wish to cover. The algorithm represents a loop generating rules until the set of computed rules $\mathcal{R}$ covers all instances in coverage set data $\mathcal{C}$, i.e. until there is no uncovered instances in $\mathcal{C}$. Each iteration of the algorithm selects an instance $\mathbf{v}$ from $C_u$. Afterwards, a generalized AXp $\mathcal{Y}$ for the prediction $c = \tau(\mathbf{v})$ by the BT $\mathfrak{T}$ (recall that $\mathfrak{T}$ is meant to compute classification function $\tau(\mathbf{x})$) is extracted by invoking Algorithm 1. The iteration proceeds by updating the set of rules $\mathcal{R}$ and the set of uncovered instances $C_u$. The algorithm terminates when all the instances in the coverage set $\mathcal{C}$ are covered and returns a compiled DS $\mathcal{R}$.

**Proposition 8.** Let $\mathfrak{T}$ be a BT and $\mathcal{R}$ be a DS returned by Algorithm 2 for $\mathfrak{T}$. Then $\mathcal{R} \equiv \mathfrak{T}$ with respect to $\mathcal{C}$.

We consider two usages of the algorithm: for exhaustive compilation the coverage set $\mathcal{C} = \mathbb{F}$ is all possible feature combinations (in practice we model this coverage set implicitly, rather than in its explicit exponential sized form), and for training set compilation where $\mathcal{C} = \mathcal{E}$ is the training set. Based on the properties of prime implicants, Proposition 8 states that as a generalized AXp $\mathcal{Y} \in \mathcal{R}$ is a formal explanation for a prediction made by BT $\mathfrak{T}$, a compiled DS captures the semantics of the original model $\mathfrak{T}$ on coverage set $\mathcal{C}$, assuming everything else is a don't care. Furthermore, if the process is applied subject to coverage set $\mathcal{C} = \mathbb{F}$, i.e. when we target the entire feature space $\mathbb{F}$, then $\mathcal{R}$ and $\mathfrak{T}$ behave identically, i.e. they compute the same classification function $\tau(\mathbf{x})$.

**Corollary 9.** Let Algorithm 2 return a DS $\mathcal{R}$ for a BT $\mathfrak{T}$. Then there is no instance in feature space $\mathbb{F}$ covered by two distinct rules $\mathcal{Y}_1, \mathcal{Y}_2 \in \mathcal{R}$ predicting inconsistent classes $c_1 \neq c_2$.