32:10 Assembly Line Preliminary Design Optimization for an Aircraft

<figure><img src="image_304.png" /><figcaption>(a) Solution $sol_{3machines}$, with one machine for each skill, (b) Solution $sol_{4machines}$, with two machines for skill $s_1$, a takt-time equal to 4 and a leadtime equal to 8.</figcaption></figure>

<caption>Figure 2 Two solution examples that each has two workstations.</caption>

## 4.2 Constraint Programming Encoding

We present here a Constraint Programming encoding for the ALPDP problem. For this encoding, we suppose that we are given two additional input parameters: an upper bound for the number of workstations in the assembly line, denoted $nW$, and an upper bound $H$ for the leadtime.

We base our CP encoding on data structures and functions that are available in the Optimization Programming Language (OPL - [28]). In particular, we use interval variables and state functions. An interval variable encompasses a start date variable and a duration variable. It can be optional, i.e. it might not be present in the schedule, and it can be temporally bounded. State functions allow to express values that a function should take over given temporal intervals.

We consider the following decision variables:

- for each activity $t \in T$, $itv_t$ is an interval variable in $[0, H]$ that represented the execution of activity $a$. For atomic activities $a \in A$, the duration of the interval is fixed to $dur_a$;

- for each atomic activity $a \in A$, for each workstation $w \in [1..nW]$, $itv_{a,w}$ is an optional interval variable in $[0, H]$ with a duration fixed to $dur_a$. The interval $itv_{a,w}$ is present if and only if $a$ starts being performed in workstation $w$;

- for each workstation $w \in [1..nW]$, $itv_w$ is an optional interval variable in $[0, H]$. $itv_w$ is present if and only if $w$ is a used workstation, i.e. a workstation in which some activities are performed;

- $takt$ is an integer variable in $[1, H]$;

- for each skill $s \in S$, for each workstation $w \in [1..nW]$, $skUsed_{s,w}$ is a boolean variable that indicates whether at least one machine with skill $s$ is assigned to workstation $w$ and $nMachines_{s,w}$ is an integer variable that represents the number of machines with skill $s$ that are assigned to $w$. An upper bound for that variable is the number of atomic activities.

- for each zone $z \in Z$, $stateOcc_z$ is a state function that represents the occupation state of zone $z$. Such a state is 1 when $z$ is occupied by an activity being performed and is equal to 0 if $z$ is neutralized by an activity.

With those decision variables, the ALPDP constraints can be encoded as follows.
## Assembly Line Consistency Constraints.

$$ \forall w \in [1..nW], \quad \text{sizeOf}(itv_w, 0) \leq \text{takt} \qquad (1) $$

$$ \forall w \in [1..nW], \quad \text{sizeOf}(itv_w, H) \geq \text{takt} \qquad (2) $$

$$ \text{startOf}(itv_1, 0) = 0 \qquad (3) $$

$$ \forall w \in [2..nW], \quad \text{endAtStart}(itv_{w-1}, itv_w) \qquad (4) $$

$$ \forall w \in [2..nW], \quad \text{presenceOf}(itv_{w-1}) \geq \text{presenceOf}(itv_w) \qquad (5) $$

Constraints (1) and (2) together ensure that the duration of a used workstation is equal to the takt-time. The second parameter in the `sizeOf` function defines the value of an interval variable that is not present in the final schedule. Through Constraint (3), the first workstation starts at time 0. Constraint (4) prevents temporal holes between successive workstations. Finally, Constraint (5) ensures that no unused workstation is placed in between used workstations.

## Activities and Workstation Constraints.

$$ \forall c \in \mathcal{C}, \quad \text{span}(itv_c, \{itv_a|a \in \mathcal{A}_c\}) \qquad (6) $$

$$ \forall (t, t') \in \mathcal{P}, \quad \text{endBeforeStart}(itv_t, itv_{t'}) \qquad (7) $$

$$ \forall a \in \mathcal{T}, \quad \text{alternative}(itv_a, \{itv_{a,w}|w \in [1..nW]\}) \qquad (8) $$

$$ \forall a \in \mathcal{A}, \forall w \in [1..nW], \quad \text{startBeforeStart}(itv_w, itv_{a,w}) \qquad (9) $$

$$ \forall a \in \mathcal{A} \text{ s.t. } S_a \neq \emptyset, \forall w \in [1..nW], \quad \text{endBeforeEnd}(itv_{a,w}, itv_w) \qquad (10) $$

$$ \forall a \in \mathcal{T}, \forall w \in [1..nW], \quad \text{presenceOf}(itv_{a,w}) \leq \text{presenceOf}(itv_w) \qquad (11) $$

$$ \forall w \in [1..nW], \quad \text{presenceOf}(itv_w) \leq \sum_{a \in \mathcal{A}} \text{presenceOf}(itv_{a,w}) \qquad (12) $$

$$ \forall a \in \mathcal{A}, \quad \text{endOf}(itv_a) \leq \max_{w \in [1..nW]} \text{endOf}(itv_w, 0) \qquad (13) $$

Constraint (6) makes interval variables associated to composite activities span over interval variables of their children. Constraint (7) enforces precedence constraints associated with precedence relation $\mathcal{P}$. Constraint (8) uses the alternative constraint in CP Optimizer and guarantees that there exists exactly one workstation $w$ such that $itv_{a,w}$ coincides with interval $itv_a$ (only for atomic activities). Constraint (9) ensures the consistency between start dates of workstation $w$ and of activity $a$ if the latter starts in workstation $w$. Constraint (10) ensures a similar consistency for end dates but only for activities that require machines. Constraints (11) and (12) guarantees the consistency between presence of activities on workstations intervals and presence of workstation intervals. Constraint (13) ensures that workstations on which atomic activities finish are used.
32:12 Assembly Line Preliminary Design Optimization for an Aircraft

Machines Constraints.

$$
\forall s \in \mathcal{S}, \forall w \in [1..nW], \quad \sum_{a \in \mathcal{A}, s \in \mathcal{S}_a} \mathrm{pulse}(itv_{a,w}, 1) \le n\mathrm{Machines}_{s,w} \quad (14)
$$

$$
\forall s \in \mathcal{S}, \forall w \in [1..nW], \quad n\text{Machines}_{s,w} \geq \text{skUsed}_{s,w} \quad (15)
$$

$$
\forall s \in \mathcal{S}, \forall w \in [1..nW], \quad n\text{Machines}_{s,w} \le |\mathcal{A}| \cdot \text{skUsed}_{s,w} \quad (16)
$$

$$
\forall (s, s') \in \mathcal{E} \text{ s.t. } s \neq s', \forall w \in [1..nW], \quad skUsed_{s,w} + skUsed_{s',w} \le 1 \quad (17)
$$

$$
\forall s \in \mathcal{S} \text{ s.t. } (s, s) \in \mathcal{E}, \forall w \in [1..nW], \quad n\text{Machines}_{s,w} \le 1 \quad (18)
$$

$$
\forall s \in \mathcal{S} \text{ s.t. } \exists a \in \mathcal{A} \text{ with } s \in \mathcal{S}_a, \quad \sum_{w=1}^{nW} skUsed_{s,w} \ge 1 \quad (19)
$$

$$
\sum_{w=1}^{nW} \sum_{s \in \mathcal{S}} n\text{Machines}_{s,w} \geq M_{\text{UB}} \quad (20)
$$

$$
\text{with } M_{\text{UB}} = |\{s \in \mathcal{S} | \exists a \in \mathcal{A}, s \in \mathcal{S}_a\}|
$$

Through Constraint (14), we build the consumption profile of machines with skill *s* by using the *pulse* primitive and we ensure that this consumption does not exceed the number of machines with skill *s* in the workstation. Constraints (15) and (16) express the relationship between *nMachines* and *skUsed* variables: *skUsed*<sub>s,w</sub> is equal to 0 if and only if the number of machines is equal to 0. Constraint (17) forbids to use excluded skills in the same workstation. Constraint (18) handles the specific case in which exclusion targets the same skill. Constraints (19) and (20) give an upper bound for the total number of machines.

Zones Constraints.

$$
\forall z \in \mathcal{Z} \text{ s.t. } \mathrm{capa}_z = 1, \quad \mathrm{noOverlap}(\{itv_a|\mathrm{occ}_{a,z} > 0\}) \quad (21)
$$

$$
\forall z \in \mathcal{Z} \text{ s.t. } \mathrm{capa}_z > 1, \quad \sum_{\substack{a \in \mathcal{A} \\ \mathrm{occ}_{a,z} > 0}} \mathrm{pulse}(itv_a, \mathrm{occ}_{a,z}) \le \mathrm{capa}_z \quad (22)
$$

$$
\forall a \in \mathcal{A}, \forall z \in \mathcal{Z} \text{ s.t. } \mathrm{occ}_{a,z} > 0, \quad \mathrm{alwaysEqual}(\mathrm{stateOcc}_z, itv_a, 1) \quad (23)
$$

$$
\forall a \in \mathcal{A}, \forall z \in \mathcal{Z}_a^{\text{neutr}}, \quad \text{alwaysEqual}(\text{stateOcc}_z, itv_a, 0) \quad (24)
$$

Constraints (21) and (22) express that the capacity of zones is never exceeded (con-
straints (21) considers the specific case when zone have a capacity equal to 1). Finally,
Constraints (23) and (24) enforce the state functions to be equal to 0 or 1 depending of
occupation and neutralization of activities. As state function have one value for each instant,
this prevents activities occupying a zone and activities neutralizing it to overlap.

We consider the following four criteria:

$$
\text{minimize} \quad takt \tag{25}
$$

$$
\text{minimize} \quad \max_{w \in [1..nW]} \text{endOf}(itv_w, 0) \tag{26}
$$

$$
\text{minimize} \quad \sum_{w \in [1..nW]} \text{presenceOf}(itv_w) \qquad (27)
$$

$$
\text{minimize} \quad \sum_{w \in [1..nW]} \sum_{s \in \mathcal{S}} n\text{Machines}_{s,w} \qquad (28)
$$

Criterion (25), (26), (27) and (28) respectively encode the minimization of the rate, the
leadtime, the number of workstations in the assemby line and the total number of machines
used. Note that, because of the equation takt · |W| = Tmax, if the takt-time and the leadtime
are fixed, so is the number of workstations. This means that one of the three first criteria is
redundant. Therefore, in the following, we only consider the takt-time, the leadtime and the
number of machines minimization.
# 5 Algorithmic approach

In this section, we describe how we have used the previously defined Constraint Programming encoding to explore the Pareto solutions set.

The algorithm follows an epsilon-constraint based strategy ([4]), which requires bounds for each criteria. In our problem, there are three criteria and we therefore have to find in which order we explore them, how to fix the bounds and the number of steps on each part of the front. As the resulting procedure is not trivial, we present it in this section. We believe that this exploration algorithm could be a starting point for other problems with three or more criteria. In fact, to the best of our knowledge, the literature mostly contains epsilon-constraint based strategies for problems with two criteria.

In order to compute bounds for the criteria values, we use a lexicographic objective function. Then, for a given number of machines in the assembly line, we compute a takt-time step that allows to explore solutions. From the set of all computed solutions, we extract the Pareto ones.

The algorithm is formally described in Listing 1. In the pseudo-code, a run of the Constraint Programming encoding is represented through the function **cpSolve** that requires two parameters. The first one is the set of constraints to satisfy and the second one is the objective function to optimize. Such an objective function can be simple (one criterion) or a lexicographic order over several criteria (denoted *lex(...)*). In the pseudo-code, we consider that we have a data structure for storing solution. If *sol* is a solution, then *sol.machines*, *sol.takt* and *sol*.T<sub>max</sub> respectively return the number of machines, the takt-time and the leadtime values associated with *sol*.

The algorithm's inputs are all constraints ((1) to (24)) of an ALPDP instance (denoted *P*) and a maximum number of points in the Pareto front for each number of machines. We start by computing three solutions *sol*<sup>*M*</sup>, *sol*<sup>*T*</sup> and *sol*<sup>*L*</sup> that respectively minimize the number of machines, the takt and the leadtime values (lines 2-4). Solution *sol*<sup>*M*</sup> provides a lower bound for the number of machines. The three solutions are added to the set of solutions through the function *updateFront* that maintains the set of Pareto solutions (line 5). While minimizing the takt and the leadtime, we also minimize the number of machines as the last criterion of a lexicographic objective function, which allows to compute an upper bound of the number of machines (line 6).

Next, we go through each value of the total number of machines, starting from the upper-bound. For each value, we consider the associated constraint $c_{MACHINES}$ that limits the number of machines in the assembly line (line 8). Then, we compute a lower bound and an upper bound for the takt-time, respectively through $sol_1$ and $sol_2$ (lines 9-10), that are stored in variables *takt*<sub>1</sub> and *takt*<sub>2</sub>. If there is a Pareto front to explore (difference between both takt-times strictly greater than 1), we compute a step, denoted $\delta$, for decreasing the considered Stakt-time according to the number of points in input (line 14). Starting from the takt upper-bound, we consider a constraint $c_{TAKT}$ that fixes the takt-time and look for a solution satisfying that constraint (line 18). If it exists, the takt-value is decreased from $\delta$ (line 22) and this step is iterated until reaching the lower bound. Otherwise, there is no solution with a smaller takt value for that number of machines and the loop is stopped. The algorithm finally returns the Pareto front.

Note that the number of workstations criteria is never used in the algorithm. In fact, if two of the three criteria *takt*, *nWorkstations* and $T_{max}$ are fixed, so is the third. As presented in the Experiments Section, results tend to show that **cpSolve** is less efficient when it comes to minimize the number of workstations. Therefore, we use the two other criteria.
■ **Algorithm 1** Epsilon approach for ALPDP.

1: **function epsilonParetoALPDP(P, nMaxPointsPerMachine)**
2: &nbsp;&nbsp;&nbsp;&nbsp;sol<sup>M</sup> &larr; **cpSolve**(P, MACHINES)
3: &nbsp;&nbsp;&nbsp;&nbsp;sol<sup>T</sup> &larr; **cpSolve**(P, lex(TAKT, LEADTIME, MACHINES))
4: &nbsp;&nbsp;&nbsp;&nbsp;sol<sup>L</sup> &larr; **cpSolve**(P, lex(LEADTIME, TAKT, MACHINES))
5: &nbsp;&nbsp;&nbsp;&nbsp;front &larr; updateFront( {}, {sol<sup>M</sup>, sol<sup>R</sup>, sol<sup>T</sup>} )
6: &nbsp;&nbsp;&nbsp;&nbsp;m &larr; max( sol<sup>L</sup> . machines, sol<sup>T</sup> . machines)
7: &nbsp;&nbsp;&nbsp;&nbsp;**while** m &ge; sol<sup>M</sup> . machines **do**
8: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;c<sub>MACHINES</sub> &larr; CONSTRAINT( &sum;<sub>w&isin;[1..nW]</sub> &sum;<sub>s&isin;S</sub> nMachines<sub>s,w</sub> &le; m )
9: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;sol<sub>1</sub> &larr; **cpSolve**(P &cup; {c<sub>MACHINES</sub>}, lex(TAKT, LEADTIME))
10: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;sol<sub>2</sub> &larr; **cpSolve**(P &cup; {c<sub>MACHINES</sub>}, lex(LEADTIME, TAKT))
11: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;front &larr; updateFront( front, {sol<sub>1</sub>, sol<sub>2</sub>} )
12: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;takt<sub>1</sub> &larr; sol<sub>1</sub> . takt; takt<sub>2</sub> &larr; sol<sub>2</sub> . takt
13: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** takt<sub>2</sub> &minus; takt<sub>1</sub> &gt; 1 **then**
14: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&delta; = &lfloor; takt<sub>2</sub> &minus; takt<sub>1</sub> &rfloor;
15: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;t &larr; takt<sub>2</sub> &minus; &delta;
16: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;sol &larr; sol<sub>1</sub>
17: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**while** t &gt; takt<sub>1</sub> &and; sol &ne; nil **do**
18: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;c<sub>TAKT</sub> &larr; CONSTRAINT( takt = t )
19: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;sol &larr; **cpSolve**(P &cup; {c<sub>MACHINES</sub>, c<sub>TAKT</sub>}, LEADTIME)
20: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** sol &ne; nil **then**
21: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;front &larr; updateFront( front, {sol} )
22: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;t &larr; t &minus; &delta;
23: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;m &larr; m &minus; 1
24: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**return** front

Functions *updateFront* and *dominates* are quite straightforward and therefore not detailed in the paper. Note that the Pareto front does not take into account the number of workstations.

6 Experiments

This section is dedicated to the experiments presentation. We first describe the benchmarks that we have used. We present two sets of experiments. First, we have solved the CP encoding with each criteria alone. This allows us to customize the **cpSolve** procedure in terms of heuristics and time out. Then, we present and analyze results associated with the Pareto solutions computation.

6.1 Benchmark Description

We have worked on a real industrial use-case dealing with the airframe assembly line of an aircraft. The airframe is the physical structure that supports all the other components of the aircraft, including avionics, passenger and cargo compartments. The materials used in its construction can include aluminium alloys, composites, titanium and other high-strength materials. The materials are cut and shaped into the various components using tools such as saws and drills. The assembly process for the airframe components typically involves careful alignment and attachment using specialised tools and techniques.
We were given three aircraft designs candidates, denoted *Design 1*, *Design 2* and *Design 3*. Features of instances are detailed in Table 1. Designs 1 and 2 have the same granularity level in terms of macro-tasks and have to be compared with each other. Design 3 is a detailed version of Design 1 with more fine-grain tasks that allows to study the scalability of the approach. These benchmarks are available at https://github.com/stephroussel/assemblyLine.

<table><caption>**Table 1** Instances features: number of zones, skills, incompatible skills, tasks and atomic tasks, precedence, maximum number of workstations and time horizon.</caption><thead><tr><th>Instance</th><th>|Z|</th><th>|S|</th><th>|E|</th><th>|T|</th><th>|A|</th><th>|P|</th><th>nW</th><th>H</th></tr></thead><tbody><tr><td>Design 1</td><td>48</td><td>5</td><td>6</td><td>176</td><td>153</td><td>186</td><td>20</td><td>40000</td></tr><tr><td>Design 2</td><td>48</td><td>5</td><td>6</td><td>187</td><td>187</td><td>279</td><td>20</td><td>40000</td></tr><tr><td>Design 3</td><td>48</td><td>5</td><td>6</td><td>628</td><td>461</td><td>417</td><td>20</td><td>40000</td></tr></tbody></table>

In this use case, the number of skills is the same in each instance. In fact, all the designs require the same types of machines but these machines are not used the same way.

## 6.2 Single Objective Experiments

**Experiments setup.** We use IBM CP Optimizer 22.1.1 through the Java API with a timeout equal to 5 minutes. Experiments were all run on a 20-core Intel(R) Xeon(R) CPU E5-2660 v3 @ 2.60GHz, 62GB RAM.

In order to test the CP encoding, we have first experimented it on the instances by considering one unique criterion for every run. The objective was to improve the solver setup and specifically the search strategy.

In CP Optimizer, the default search strategy is generally efficient. However, in our problem, we noticed a real performance improvement when asking the solver to instantiate first the *takt* variable. In fact, it allows to fix the duration of workstations, which is a key element for the remaining variables. Once this variable is fixed, we have tried other strategies such as deciding the presence of workstations or the number of machines in each workstation but it globally downgraded the performance compared to the solver default strategy. Therefore, we only present in this paper results associated with the *takt* variable instantiation strategy.

Results for the mono-objective solving are presented in Table 2 for each design and each criteria. The solver finds the optimal solution in several cases, even for the largest instance. The optimal solution is either found in less than a minute or reached the time out. Note that the solver was unable to establish the optimality of the solution with respect to the number of workstations. On that point, it would be possible to compute a lower bound using the skill exclusion relationships. We could, for instance, find the size of the largest clique in the exclusion skills graph.

<table><caption>**Table 2** Mono-criteria CP solving with search phase customization on the *takt* variable. For each criteria, value and time (in seconds) are given. The symbol * denotes that the solution found is optimal.</caption><thead><tr><th rowspan="2">Instance<br>Criterion</th><th colspan="2">takt</th><th colspan="2">nMachines</th><th colspan="2">T<sub>max</sub></th><th colspan="2">nWorkstations</th></tr><tr><th>Value</th><th>Time</th><th>Value</th><th>Time</th><th>Value</th><th>Time</th><th>Value</th><th>Time</th></tr></thead><tbody><tr><td>Design 1</td><td>550*</td><td>25</td><td>5*</td><td>69</td><td>3606*</td><td>59</td><td>2</td><td>300</td></tr><tr><td>Design 2</td><td>550*</td><td>32</td><td>5*</td><td>28</td><td>3536</td><td>300</td><td>2</td><td>300</td></tr><tr><td>Design 3</td><td>340*</td><td>35</td><td>8</td><td>300</td><td>3456</td><td>300</td><td>2</td><td>300</td></tr></tbody></table>
6.3 Multi Objective Results

**Experiments setup.** We did not give a global timeout to the whole algorithm but instead a 1 minute timeout to each call to **cpSolve**. Note that we have let CP Optimizer handle directly the lexicographic criteria. In fact, preliminary experiments showed that there was no real performance improvement when decomposing the resolution into optimizing the first criterion and fix it, then optimizing the next one and so on.

<figure><img src="image_305.png" /><figcaption>Figure 3 Pareto solutions obtained for Designs 1 and 2.</figcaption></figure>

Figures 3a and 3b show the Pareto front obtained respectively for Designs 1 and 2. The leadtime is on the horizontal axis, the takt-time on the vertical axis, the number of machines is represented through colors and the number of stations are the lines on each figure. Each small cross corresponds a to dominated solution. Note that for Design 2, we have not represented a Pareto solution that had 18 stations and a smaller takt-time.

As expected, for each design, the more machines on the assembly line, the better takt-time and leadtime values it is possible to get. Each figure shows the possible trade-offs that have to be made for each design. For instance, in Design 1, with 2 workstations, the best takt-time is equal to 18 hours. When adding 1 more workstation, such a value falls to 12 hours but requires many more machines in order to maintain the leadtime value. With one additional workstation and the same leadtime, the takt-time is less than 10 hours. We also observe that it is not worth considering more than 8 workstations for that design. Similar remarks can be made for Design 2.

Figure 4a presents Pareto solutions for Designs 1 and 2 when the assembly line is equipped with 5 machines. It shows that Design 1 (in blue on the Figure) allows to reach lower takt-time values. If one more machine is available (Figure 4b), then solutions associated with Design 2 dominate solutions of Design 1.

<figure><img src="image_306.png" /><figcaption>Figure 4 Comparison of Design 1 and 2 with 5 and 6 machines.</figcaption></figure>
These results have been presented to industrial partners. They have appreciated the possibility to visualize the Pareto front, which allows them to foresee the trade offs that could be made not only within the assembly line (number of workstations, number of machines, etc) but also between the aircraft design and the associated assembly line. Following results presented on Figure 4, they were surprised that Design 2 seems more promising than Design 1 only in the presence of an additional machine. Indeed, before this study, they had the false intuition than Design 2 would dominate Design 1 even with 5 machines. Such a positive feedback from end-users experts shows the added value of computing and exploring the Pareto front in this project.

# 7 Conclusion

In this paper, we have proposed a Constraint Programming based approach for supporting manufacturer in the early assembly design phase. We have formally defined the associated problem and have developed an algorithm for exploring the Pareto front. Each solution in this front is a trade off between the takt-time, the leadtime, the number of machines equipping the assembly line and the number of workstations. Such tools allow the aircraft manufacturer not only to assess the performance of aircraft candidates with respect to the assembly line performances but also to compare candidates with each other.

There are several directions for future works. First, the CP encoding could be improved by computing some lower and upper bounds for the targeted criteria. In fact, the addition of redundant constraints for boosting the solving might change the way the Pareto exploration should be performed and should be studied more deeply. Then, the Pareto front exploration could benefit from some recent works on computing representative Pareto solutions, such as [27]. It could also be possible to compare the CP approach and the Pareto exploration with evolutionary algorithms in terms of quality of results. We would also like to port the model to other CP solvers in order to test them. While most of the constraints could easily be written in a more classical language such as PyCSP3 ([17]), the neutralization constraints modeled by a state function in OPL might be a little trickier to encode in order to stay compact.

The instances we have used here come from a real assembly line. In order to test the approach more broadly, it would be possible to modify these instances by adding random resources or changing the precedence between tasks. Similarly, we could also adapt existing instances of the literature.

Finally, a last perspective is to consider uncertainty in the tasks duration. To do so, it might be worth considering coupling CP solvers with learning approaches that are particularly suited for managing uncertainty.

# References

1. Hacı Mehmet Alakaş. General resource-constrained assembly line balancing problem: conjunction normal form based constraint programming models. *Soft Computing*, 25(8):6101–6111, 2021.

2. Hacı Mehmet Alakaş, Mehmet Pınarbaşı, and Mustafa Yüzkırmızı. Constraint programming model for resource-constrained assembly line balancing problem. *Soft Computing*, 24:5367–5375, 2020.

3. Dmitry Arkhipov, Olga Battaïa, Julien Cegarra, and Alexander Lazarev. Operator assignment problem in aircraft assembly lines: a new planning approach taking into account economic and ergonomic constraints. *Procedia CIRP*, 76:63–66, 2018.
32:18 Assembly Line Preliminary Design Optimization for an Aircraft

4 Francisco Ballestín and Rosa Blanco. Theoretical and practical fundamentals for multi-objective optimisation in resource-constrained project scheduling problems. *Computers & Operations Research*, 38(1):51-62, 2011. Project Management and Scheduling. doi:10.1016/j.cor.2010.02.004.

5 Olga Battaïa and Alexandre Dolgui. A taxonomy of line balancing problems and their solution approaches. *International Journal of Production Economics*, 142(2):259-277, 2013. Anticipation of risks impacts and industrial performance evaluation in distributed organizations life cycles. doi:10.1016/j.ijpe.2012.10.020.

6 Alexander Biele and Lars Mönch. Hybrid approaches to optimize mixed-model assembly lines in low-volume manufacturing. *Journal of Heuristics*, 24(1):49-81, 2018.

7 Tamara Borreguero, Alvaro García, and Miguel Ortega. Scheduling in the aeronautical industry using a mixed integer linear problem formulation. *Procedia engineering*, 132:982-989, 2015.

8 Nils Boysen, Philipp Schulze, and Armin Scholl. Assembly line balancing: What happened in the last fifteen years? *European Journal of Operational Research*, 301(3):797-814, 2022. doi:10.1016/j.ejor.2021.11.043.

9 Peter Brucker, Andreas Drexl, Rolf Möhring, Klaus Neumann, and Erwin Pesch. Resource-constrained project scheduling: Notation, classification, models, and methods. *European journal of operational research*, 112(1):3-41, 1999.

10 Jens Buergin, Sina Helming, Jan Andreas, Philippe Blaettchen, Yannick Schweizer, Frank Bitte, Benjamin Haefner, and Gisela Lanza. Local order scheduling for mixed-model assembly lines in the aircraft manufacturing industry. *Production Engineering*, 12:759-767, 2018.

11 Yossi Bukchin and Tal Raviv. Constraint programming for solving various assembly line balancing problems. *Omega*, 78:57-68, 2018. doi:10.1016/j.omega.2017.06.008.

12 Anouck Chan, Anthony Fernandes Pires, Thomas Polacsek, and Stéphanie Roussel. The aircraft and its manufacturing system: From early requirements to global design. In Xavier Franch, Geert Poels, Frederik Gailly, and Monique Snoeck, editors, *Advanced Information Systems Engineering - 34th International Conference, CAiSE 2022, Leuven, Belgium, June 6-10, 2022, Proceedings*, volume 13295 of *Lecture Notes in Computer Science*, pages 164-179. Springer, 2022. doi:10.1007/978-3-031-07472-1_10.

13 Hicham Chehade, Alexandre Dolgui, Frédéric Dugardin, Lina Makdessian, and Farouk Yalaoui. Multi-objective approach for production line equipment selection. *Management and Production Engineering Review*, 3(1):4-17, 2012.

14 G. Heike, M. Ramulu, E. Sorenson, P Shanahan, and Kamran Moinzadeh. Mixed model assembly alternatives for low-volume manufacturing: the case of the aerospace industry. *International Journal of Production Economics*, 72(2):103-120, 2001.

15 Damla Kizilay and Zeynel Abidin Çil. Constraint programming approach for multi-objective two-sided assembly line balancing problem with multi-operator stations. *Engineering Optimization*, 53(8):1315-1330, 2021. doi:10.1080/0305215X.2020.1786081.

16 Philippe Laborie, Jérôme Rogerie, Paul Shaw, and Petr Vilím. IBM ILOG CP optimizer for scheduling: 20+ years of scheduling with constraints at IBM/ILOG. *Constraints*, 23:210-250, 2018.

17 Christophe Lecoutre and Nicolas Szczepanski. PyCSP3: Modeling Combinatorial Constrained Problems in Python. Technical report, arXiv, December 2021. URL: https://hal-univ-artois.archives-ouvertes.fr/hal-03701203.

18 Fernando Mas, José Luis Menéndez, Manuel Oliva, Javier Servan, Rebeca Arista, and Carmelo Del Valle. Design within complex environments: Collaborative engineering in the aerospace industry. In *Information System Development: Improving Enterprise Communication*, pages 197-205. Springer, 2014.

19 Jonathan Oesterle and Lionel Amodeo. Efficient multi-objective optimization method for the mixed-model-line assembly line design problem. *Procedia CIRP*, 17:82-87, 2014. Variety Management in Manufacturing. doi:10.1016/j.procir.2014.01.038.
20 Jonathan Oesterle, Lionel Amodeo, and Farouk Yalaoui. A comparative study of multi-objective algorithms for the assembly line balancing and equipment selection problem under consideration of product design alternatives. *Journal of intelligent Manufacturing*, 30:1021–1046, 2019.

21 Tenda Okimoto, Yongjoon Joe, Atsushi Iwasaki, Toshihiro Matsui, Katsutoshi Hirayama, and Makoto Yokoo. Interactive algorithm for multi-objective constraint optimization. In Michela Milano, editor, *Principles and Practice of Constraint Programming*, pages 561–576, Berlin, Heidelberg, 2012. Springer Berlin Heidelberg.

22 Thomas Polacsek, Stéphanie Roussel, Cédric Pralet, and Claude Cuiller. Design for efficient production, A model-based approach. In Manuel Kolp, Jean Vanderdonckt, Monique Snoeck, and Yves Wautelet, editors, *13th International Conference on Research Challenges in Information Science, RCIS 2019, Brussels, Belgium, May 29-31, 2019*, pages 1–6. IEEE, 2019. doi:10.1109/RCIS.2019.8877088.

23 Cédric Pralet, Stéphanie Roussel, Thomas Polacsek, François Bouissière, Claude Cuiller, Pierre-Eric Dereux, Stéphane Kersuzan, and Marc Lelay. A scheduling tool for bridging the gap between aircraft design and aircraft manufacturing. *Proceedings of the International Conference on Automated Planning and Scheduling*, 28(1):347–355, June 2018. doi:10.1609/icaps.v28i1.13910.

24 Brahim Rekiek, Alain Delchambre, Alexandre Dolgui, and Antoneta Bratcu. Assembly line design: A survey. *IFAC Proceedings Volumes*, 35(1):155–166, 2002. 15th IFAC World Congress. doi:10.3182/20020721-6-ES-1901.01647.

25 Emma Rollon and Javier Larrosa. Multi-objective propagation in constraint programming. *Frontiers in Artificial Intelligence and Applications*, 141:128, 2006.

26 Tamara Borreguero Sanchidrián, Tom Portoleau, Christian Artigues, Alvaro García Sánchez, Miguel Ortega Mier, and Pierre Lopez. Exact and heuristic methods for an aeronautical assembly line time-constrained scheduling problem with multiple modes and a resource leveling objective. *hal-03344445*, 2021.

27 Nicolas Schwind and Demirović Emir. Representative solutions for bi-objective optimisation. In *Proceedings of the AAAI Conference on Artificial Intelligence*, volume 34, pages 1436–1443, 2020.

28 Pascal Van Hentenryck, Laurent Michel, Laurent Perron, and Jean-Charles Régin. Constraint programming in opl. In *PPDP*, volume 99, pages 98–116. Springer, 1999.

29 Yeo Keun Kim, Yong Ju Kim, and Yeongho Kim. Genetic algorithms for assembly line balancing with various objectives. *Computers & Industrial Engineering*, 30(3):397–409, 1996. IE in Korea. doi:10.1016/0360-8352(96)00009-5.
<figure><img src="image_307.png" /></figure>
# SAT-Based Learning of Compact Binary Decision Diagrams for Classification

Pouya Shati

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

Department of Computer Science, University of Toronto, Canada
Vector Institute, Toronto, Canada

Eldan Cohen

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

Department of Mechanical and Industrial Engineering, University of Toronto, Canada

Sheila McIlraith

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

Department of Computer Science, University of Toronto, Canada
Vector Institute, Toronto, Canada

## Abstract

Decision trees are a popular classification model in machine learning due to their interpretability and performance. However, the number of splits in decision trees grow exponentially with their depth which can incur a higher computational cost, increase data fragmentation, hinder interpretability, and restrict their applicability to memory-constrained hardware. In constrast, binary decision diagrams (BDD) utilize the same split across each level, leading to a linear number of splits in total. Recent work has considered optimal binary decision diagrams (BDD) as compact and accurate classification models, but has only focused on binary datasets and has not explicitly optimized the compactness of the resulting diagrams. In this work, we present a SAT-based encoding for a multi-terminal variant of BDDs (MTBDDs) that incorporates a state-of-the-art direct encoding of numerical features. We then develop and evaluate different approaches to explicitly optimize the compactness of the diagrams. In one family of approaches, we learn a tree BDD first and model the size of the diagram the tree will be reduced to as a secondary objective, in a one-stage or two-stage optimization scheme. Alternatively, we directly learn diagrams that support multi-dimensional splits for improved expressiveness. Our experiments show that direct encoding of numerical features leads to better performance. Furthermore, we show that exact optimization of size leads to more compact solutions while maintaining higher accuracy. Finally, our experiments show that multi-dimensional splits are a viable approach to achieving higher expressiveness with a lower computational cost.

**2012 ACM Subject Classification** Mathematics of computing → Combinatorial optimization; Computing methodologies → Machine learning

**Keywords and phrases** Binary Decision Diagram, Classification, Compactness, Numeric Data, MaxSAT

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

**Supplementary Material** Software (Source code): https://github.com/PouyaShati/BDD

**Funding** We gratefully acknowledge funding from Natural Sciences and Engineering Research Council of Canada (NSERC) and the CIFAR AI Chairs program (Vector Institute).

## 1 Introduction

Classifiers are complete functions for assigning labels to datapoints, learned from a limited set of supervised training data. A classifier is trained to focus on informative aspects of the input and extract patterns from the data to make decisions. In complex black-box classifiers such as deep and convoluted neural networks, it is often challenging to understand the influential features of the data and the inner-workings of the decision-making [36, 38]. In contrast, interpretable classifiers such as decision trees [1, 26, 33, 34] and diagrams [16, 24] are concise and easy to understand [10]. The simplicity in interpretable solutions makes them the

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

© Pouya Shati, Eldan Cohen, and Sheila McIlraith;
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. 33; pp. 33:1–33:19

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

LIPICS

Leibniz International Proceedings in Informatics

Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl Publishing, Germany
perfect candidates for when we need to formally analyze or explain the classifier's behavior
[27, 19, 20]. Surprisingly, the benefit of interpretability does not come with significant cost
to accuracy in many applications [31, 25].

Binary decision diagrams (BDD) are graph representations of functions with binary inputs
and are widely used in logical synthesis and formal verification methods [7, 2, 28, 23]. While
equivalent to truth tables in purpose, BDDs are more compact since redundant sub-tables
can be eliminated and merged [23]. A multi-terminal BDD (MTBDD) is a BDD variant that
supports multiple outputs [11].

Decision trees are the most common form of interpretable classifiers (e.g., [7, 29, 30, 3,
37, 18, 5, 33, 34]). However, the number of splits double at each level of a decision tree,
making it challenging to interpret the solution as depth increases [15]. The large number
of splits in decision trees can hinder the learning performance given that the search space
grows exponentially with each level. Further, deep splits that only affect a small number of
datapoints can cause data fragmentation and overfitting [35, 21, 14]. Lastly, the exponential
number of splits in decision trees restricts their applicability to memory-constrained hardware
[35]. The sequential nature of decision-making in BDDs resembles that of decision trees, but
unlike decision trees, the same split is used across a level, leading to a linear number of splits.

In this paper, we choose BDDs as our interpretable classifiers in order to emphasize compactness to an even greater degree compared to popular interpretable approaches, e.g. ones based on decision trees. Alongside the primary objective of accuracy, the size of a BDD classifier is also encoded and optimized. Our encoding for learning BDDs is inspired by the encoding of numerical branchings in Shati et al. [33, 34], which allows us to directly learn splits over numeric features without explicitly binarizing them. We then expand the notion of learning splits to multiple dimensions in order to learn solutions from a limited but distinctively interpretable family of diagrams.

The main contributions of this paper are as follows:

1. We present a novel SAT encoding of maximum accuracy BDD classifiers for numerical datasets. Our model represents a non-reduced BDD which can be reduced according to its sequence of terminals.

2. We extend our encoding by modelling the size of the final reduced BDD. The size encoding can be used as a secondary weighted objective, or it can be optimized in a second stage.

3. We present a variant of our encoding which supports direct learning of diagrams through multi-dimensional splits. Multi-dimensional BDDs enable learning more efficiently as they provide more expressive solutions within the same number of splits.

4. We run extensive experiments which demonstrate that our base encoding outperforms the state of the art in runtime and accuracy. Additionally, we show that explicitly optimizing size leads to significantly more compact solutions while maintaining high accuracy. Finally, we show that multi-dimensional BDDs scale better than BDDs due to their expressiveness and the size of their encoding.

2 Related Work

Decision tree classifiers are traditionally constructed via local search and heuristics [7, 29, 30].
Recent advances in tools and techniques for exact optimization have enabled approaches for
finding optimal decision trees via branch-and-bound search [1], SAT-based encodings [33, 18,
3], Constraint Programming [37], or Mixed Integer Programming [5]. Exact optimization
produces solutions that are optimal in accuracy, size, or both [33]. The size of a decision tree is
often measured by its depth or number of nodes, which does not take into account the amount
of redundancy between nodes. Decision diagrams address redundancy directly by merging equivalent nodes. Approaches to learn optimal decision diagrams via exact optimization usually require a pre-determined skeleton as input. For example, Florio et al. [14] requires user-provided width for each level, limiting solutions but guaranteeing compactness and improving performance. Oblivious decision diagrams encourage compactness even further by requiring all splits of the same level to be the same, leading to only a linear number of splits instead of exponential as in ordinary decision trees.

Binary Decision Diagrams (BDD) are oblivious decision diagrams that can be reduced and ordered, making them ideal candidates to represent boolean functions. While BDDs are commonly used for hardware synthesis [8, 2, 28, 23], there has been a recent focus on utilizing BDD classifiers as interpretable solutions [16, 9]. For example, Hu et al. [16] uses MaxSAT to learn a maximum accuracy tree BDD which will then be reduced into a diagram. On the other hand, Cabodi et al. [9] considers a SAT-based approach for learning BDDs of minimum size that correctly classify all the training data.

# 3 Background

## 3.1 Binary Decision Diagrams

In this section, we first define standard Binary Decision Diagrams which can be used for binary classification. To support multi-class classification, we extend the BDD definition to its multi-terminal variant.

**► Definition 1 (Binary Decision Diagram).** *Given a boolean function operating on boolean input, a Binary Decision Diagram (BDD) is a graph representation of the function. A BDD is a rooted, directed, and acyclic graph with sink nodes as terminals ($N_T$) and the rest as decision nodes ($N_D$). Each decision node is labelled with a split and has two outgoing edges, corresponding to the two (1/0) possible values of the boolean input. The terminal nodes are assigned output values 0 or 1. The size of a BDD is the number of its decision nodes.*

**► Definition 2 (Multi-Terminal Binary Decision Diagram).** *Given a function operating on binary features with constant range $[1..k]$, a Multi-Terminal Binary Decision Diagram (MTBDD) is a graph representation of the function. An MTBDD is a BDD with its terminals supporting $k$ output values instead of 2.*

Throughout the paper, we will simply refer to MTBDDs as BDDs to be more concise. We adopt standard BDD terminology with minor exceptions. As a concise way of addressing the nodes, we borrow from the decision tree literature to view the outgoing edges as parent-child relationships. Specifically, we address the node connected via the 1 (0) outgoing edge as the left (right) child. Moreover, we refer to the boolean inputs assigned to decision nodes as splits rather than features or variables, to avoid confusion with the features of our data and the variables of our encoding. Finally, considering that the features in our data are numerical rather than binary, we overload the definition of splits to represent pairs of numerical features and a valid threshold.

**► Definition 3 (Split).** *Given a finite set of numerical features $F$, a split is a binarization of a numerical feature through pairing with a threshold value ($f, \alpha$) $\in F \times \mathbb{R}$. A split assigns a value of 1 to each point $x$ that comes before the threshold ($x[f] \le \alpha$). Similarly, it assigns a value of 0 to each point $x$ that comes after the threshold ($x[f] > \alpha$).*
<figure><img src="image_313.png" /><figcaption>Figure 1 Two BDDs representing the same boolean function. Solid (dotted) lines correspond to value 1 (0).</figcaption></figure>

In order to calculate the value of a function for a given input using its BDD representation, we start from the root, move to the left (right) child if the input has value 1 (0) for the current split, and assign an output when a terminal is reached. A node is said to contain an input if the node is on the input's path from the root to a terminal.

The same function with binary splits can be represented with multiple BDDs. In order to uniquely represent a function, we define the BDD properties of being *ordered* and *reduced*. The ordered property enforces the same sequence of splits along each path, making BDDs practically equivalent to Oblivious Read-Once Decision Graphs (OODGs). The reduced property requires all equivalent parts of the BDD to be merged and all of the redundant parts to be removed. Given a function and an ordering of splits, there only exists one ordered and reduced BDD representing the function.

**Definition 4 (Ordered BDD and Split Sequence).** A BDD is said to be ordered if the splits observed along every path from the root to a terminal respect a singular total ordering, called its split sequence.

**Definition 5 (Node Equivalency and Redundancy).** Given a binary decision diagram $\mathcal{T}$ and two of its nodes $t_1, t_2 \in \mathcal{N}$, $t_1$ and $t_2$ are equivalent if they are both decision nodes with the same split and the sub-graph of $t_1$ and its descendants is isomorphic to the sub-graph of $t_2$ and its descendants, or if they are both terminal nodes with the same label. Furthermore, a decision node $t \in \mathcal{N}_D$ is redundant if its left and right children are equivalent.

**Definition 6 (Reduced BDD).** A BDD is said to be reduced if all of its equivalent nodes are merged and all of its redundant nodes are replaced with their children.

We assume a BDD to be ordered and reduced unless noted otherwise. Figure 1 depicts two BDDs representing the same boolean function, one ordered and reduced and the other one not. It has been shown in the literature that independent of the merging and elimination process, the final result of reducing a BDD is always the same. We formally state the uniqueness of reduced BDDs in Proposition 7 and refer the reader to the BDD literature for more details [23].<fnref n="1" />

**Proposition 7.** Every function operating on splits and a given split sequence is represented by exactly one ordered and reduced BDD up to renaming.

<footnote n="1">¹ Note that when employing standard BDD terminology, redundant and equivalent nodes are defined through binary sequences called *beads*. Beads correspond to nodes of the uniquely reduced BDD described in Proposition 7.</footnote>
## 3.2 Weighted Partial MaxSAT

In this section, we describe the MaxSAT paradigm which we use to learn binary decision diagrams. A SAT formula is a conjunction of clauses, where each clause is a disjunction of literals, and each literal is either a Boolean variable or its negation. In a weighted partial MaxSAT instance, clauses are categorized into **hard** and **weighted soft** clauses. The goal is then to find a truth assignment to variables which satisfies all of the hard clauses and maximizes the total weight of the satisfied soft clauses.

# 4 MaxSAT Encoding of BDD Classifiers

In this section, we propose a MaxSAT encoding to find a BDD classifier with maximum accuracy. The inputs to our problem are: a training dataset $X$ over a set of numerical features $F$ and labels $K$, a ground-truth label $y_i$ for each point $x_i \in X$, and a maximum number of splits $s_{max}$. The outputs are $s_{max}$ splits $\theta \in \{F \times \mathbb{R}\}^{s_{max}}$ with terminal labelling $\gamma: \{0,1\}^{s_{max}} \to K$. Similar to previous work of Hu et al. [16], we learn an ordered yet non-reduced (tree) BDD first. We then focus on merging and replacing nodes towards a reduced BDD. Inspired by the numerical branching in Shati et al. [33, 34], we encode splits (Definition 3) without the need for prior binarization of data which was shown to lead to significant performance degradation in decision trees. We instead directly encode how each split directs each point, while making sure the order of values is respected given a selected feature.

## 4.1 Variables

The following set of binary variables represents different aspects of modeling splits, labels, and accuracy in our BDD encoding.

- $a_{s,j}$: The feature chosen at split $s$ is or comes before $j$.

- $d_{s,i}$: Point $x_i$ is directed to the left child at split $s$.

- $c_{t,l}$: Output label $l$ is assigned to terminal node $t$.

- $o_i$: Point $x_i$ is classified correctly.

## 4.2 Clauses

We propose the following sets of Conjunctive Normal Form (CNF) clauses for modelling a non-reduced ordered BDD. The clauses in Eqs. (1-2) guarantee that one feature is selected at each split by enforcing the ordered encoding of $a_{s,j}$ variables. The clauses in Eqs. (3-5) guarantee that for each split and chosen feature, the points directed to the left (resp. right) have a lesser (resp. greater) value compared to a threshold. The clauses in Eqs. (6-7) guarantee that one output label is chosen for each terminal node. The clauses in Eq. (8) guarantee that a point can only be considered classified if it ends up in a terminal node with the same label as its ground truth. We use $O_j(X)$ to denote the set of all consecutive pairs when the members of $X$ are sorted based on their $j$ values, $O_j^-(X)$ to denote the subset of pairs with equal $j$ values, i.e., $O_j(X) \cap \{(i_1, i_2) \mid x_{i_1}[j] = x_{i_2}[j]\}$, and $\#_j^1$ to denote the index of the point with the smallest $j$ value. Furthermore, we set $a_{s,|F|}$ to the false constant and define $A_R(t)$ ($A_L(t)$) as the set of right (left) ancestors of terminal $t$.
## 33:6 SAT-Based Learning of Compact Binary Decision Diagrams for Classification

$$ (a_{s,j}, \neg a_{s,j+1}) \qquad s < s_{max}, j \in F \quad (1) $$

$$ (a_{s,0}) \qquad s < s_{max} \quad (2) $$

$$ (\neg a_{s,j}, a_{s,j+1}, d_{s,i_1}, \neg d_{s,i_2}) \qquad s < s_{max}, j \in F, (i_1, i_2) \in O_j(X) \quad (3) $$

$$ (\neg a_{s,j}, a_{s,j+1}, \neg d_{s,i_1}, d_{s,i_2}) \qquad s < s_{max}, j \in F, (i_1, i_2) \in O_j^-(X) \quad (4) $$

$$ (\neg a_{s,j}, a_{s,j+1}, d_{s,\#j}) \qquad s < s_{max}, j \in F \quad (5) $$

$$ (\neg c_{t,l_1}, \neg c_{t,l_2}) \qquad t \in \mathcal{N}_T, l_1, l_2 \in K \quad (6) $$

$$ (\bigvee_{l \in L} c_{t,l}) \qquad t \in \mathcal{N}_T \quad (7) $$

$$ (\bigvee_{s \in A_L(t)} \neg d_{s,i}, \bigvee_{s \in A_R(t)} d_{s,i}, c_{t,y_i}, \neg o_i) \qquad t \in \mathcal{N}_T, x_i \in X \quad (8) $$

In order to model the objective, we include the soft clauses in Eq. (9) with unit weights, which represent correctly classifying as many training points as possible.

$$ (o_i) \qquad x_i \in X \quad (9) $$

## 4.3 Decoding

An assignment to the variables in Section 4.1 which is satisfying with regard to the hard clauses in Section 4.2 is decoded into a reduced BDD in two steps. First, the assignment is decoded into a non-reduced BDD. Specifically, the labels of terminals are decoded from the selected labels ($c_{t,l}$) and the sequence of splits are decoded from the pairings of selected features ($a_{s,j}$ values) and thresholds ($d_{s,i}$ values). Second, the resulting BDD is reduced by merging equivalent nodes and eliminating redundant nodes.

# 5 Size Optimization

In the encoding presented in Section 4, not all terminals are guaranteed to contain training points. Thus, we may end up with empty terminals, which we can relabel without affecting the training accuracy. However, the labels of such terminals can affect the prediction on unseen datapoints as well as the size of the final reduced BDD. In Hu et al. [16], the labels of empty terminals are decided arbitrarily by the solver and they investigate the impact of two post-processing relabelling heuristics on the testing accuracy (i.e., accuracy on unseen data). The first heuristic assigns the majority label of the terminal's first non-empty ancestor. The second heuristic finds and merges equivalent nodes in a greedy top-down search. Hu et al.'s experiments showed that the heuristics do not have significant impact on the testing accuracy. However, as we demonstrate in Figure 2a and Figure 2b, such heuristic approaches can increase the size of the final reduced BDDs.

In this section, we propose to use exact optimization for solving the same problem as in Section 4, with additional compactness considerations via deciding the labels of empty terminals. Specifically, we aim to model the size of the BDD after the merging of its equivalent nodes and the removal of its redundant nodes. The variables and clauses introduced in this section can either be added to the encoding in Section 4 as a secondary objective (i.e. 1-stage approach), or they can be used as a separate post-processing stage alongside variables $c_{t,l}$ and the clauses in Eqs. (6-7) with the labels of non-empty terminals fixed (i.e. 2-stage approach).
<figure><img src="image_314.png" /><figcaption>■ **Figure 2** Examples for heuristics of assigning labels to empty terminals where size is increased. Black labels represent non-empty terminals, red labels are chosen by the heuristics, and green labels are the original labels which lead to a more compact solution.</figcaption></figure>

## 5.1 Variables

The following set of binary variables represents uniqueness and repetition aspects of terminal labels, required for modelling the size of the ultimately reduced BDD.

- $\sigma_{t_1,t_2}$: Terminals $t_1$ and $t_2$ have been assigned different output labels.

- $b_{t,\Delta}$: The sequence of $\Delta$ labels starting from terminal node $t$ (inclusive) cannot be divided into two equal sub-sequences.

- $r_{t_1,t_2,\Delta}$: The sequence of $\Delta$ labels starting from terminal node $t_1$ is equal to the sequence of $\Delta$ labels starting from terminal node $t_2$ (both inclusive).

Note that the terminal order which is referred to in the variable definitions, is the order of appearance in the depth-first search of the tree with the left outgoing edges prioritized. We say a sequence of terminals correspond to a decision node, if they are the sequence of terminals in the decision node's sub-graph. Moreover, we say a single terminal coincides with a decision node if it marks the beginning the node's sequence of terminals.

Note that variables $\sigma_{t_1,t_2}$, $b_{t,\Delta}$, and $r_{t_1,t_2,\Delta}$ are not defined for all possible index combinations. Specifically, $\Delta$ always refers to the length of a sequence of terminals corresponding to a decision node (a power of 2), and can uniquely determine the node, if paired with a coinciding terminal. Consequently, $b_{t,\Delta}$ is only used to represent that the decision node at level $s_{max} - \log_2(\Delta) + 1$ coinciding with $t$ is not redundant. Moreover, $r_{t_1,t_2,\Delta}$ is only used to represent that the two decision nodes at level $s_{max} - \log_2(\Delta) + 1$ coinciding with $t_1$ and $t_2$ are equivalent. Lastly, $\sigma_{t_1,t_2}$ is only used when it affects the equivalency of two decision nodes.

## 5.2 Clauses

We first define two sets to help with the clause construction. The set $P(s_{max})$ contains all possible lengths of terminal sequences corresponding to decision nodes except for the root (i.e., all $2^{p'}$ where $p' < s_{max}$). The set $G(s_{max})$ contains all triples $(t_1, t_2, \Delta)$ where comparing the labels of terminals $t_1$ and $t_2$ has impact on a decision node at level $s_{max} - \log_2(\Delta) + 1$ being redundant.

$$ 
\begin{aligned}
G(1) &= \{(0, 1, 2)\} \\
G(p) &= G(p-1) \cup \{(t_1 + 2^{p-1}, t_2 + 2^{p-1}, \Delta) | (t_1, t_2, \Delta) \in G(p-1)\} \\
&\quad \cup \{(t, t + 2^{p-1}, 2^p) | 0 \le t < 2^{p-1}\}
\end{aligned}
$$

We propose the following sets of CNF clauses for modelling the size of a reduced BDD given its terminal labels. The clauses in Eqs. (10-12) guarantee that $\sigma_{t_1,t_2}$ variables correctly represent the difference in labels. The clauses in Eq. (13) guarantee that each decision node
33:8 SAT-Based Learning of Compact Binary Decision Diagrams for Classification

can only be redundant if all of the corresponding terminal pairs have the same label. The clauses in Eqs. (14-15) guarantee that a pair of decision nodes can only be equivalent if their corresponding sequence terminals have the same labels.

$$ (\neg c_{t_1,l}, \neg c_{t_2,l}, \neg \sigma_{t_1,t_2}) \qquad (t_1, t_2, \Delta) \in G(s_{max}), l \in K \quad (10) $$

$$ (\neg c_{t_1,l}, c_{t_2,l}, \sigma_{t_1,t_2}) \qquad (t_1, t_2, \Delta) \in G(s_{max}), l \in K \quad (11) $$

$$ (c_{t_1,l}, \neg c_{t_2,l}, \sigma_{t_1,t_2}) \qquad (t_1, t_2, \Delta) \in G(s_{max}), l \in K \quad (12) $$

$$ (b_{\Delta[t_1/\Delta],\Delta}, \neg \sigma_{t_1,t_2}) \qquad (t_1, t_2, \Delta) \in G(s_{max}) \quad (13) $$

$$ (\neg r_{t_1\Delta, t_2\Delta, \Delta}, \neg c_{t_1\Delta+\delta, l}, c_{t_2\Delta+\delta, l}) \qquad \Delta \in P(s_{max}), t_1 < t_2 < 2^{s_{max}}/\Delta, \delta < \Delta, l \in K \quad (14) $$

$$ (\neg r_{t_1\Delta, t_2\Delta, \Delta}, c_{t_1\Delta+\delta, l}, \neg c_{t_2\Delta+\delta, l}) \qquad \Delta \in P(s_{max}), t_1 < t_2 < 2^{s_{max}}/\Delta, \delta < \Delta, l \in K \quad (15) $$

To model our objective, we include the soft clauses in Eq. (16) which represent each decision node to be either redundant or equivalent to a previous one. Maximizing the number of redundant of equivalent decision nodes will consequently minimize the size of our reduced BDD. Note that if size is being considered as a secondary objective to accuracy, the clauses in Eq. (16) can be weighted against the clauses in Eq. (9) proportionately. Otherwise, if size is our second-stage objective, we can use the clauses in Eq. (16) with unit weights.

$$ (\bigvee_{0 \le t_2 < t} r_{t_2\Delta, t\Delta, \Delta}, \neg b_{t\Delta, \Delta}) \qquad \Delta \in P(s_{max}), t < 2^{s_{max}}/\Delta \quad (16) $$

5.3 Decoding

We first show the soundness of the size encoding.

► **Proposition 8.** Given a sequence of terminal labels $c_{t,l}$, the encoding in Sections (5.1,5.2) has an optimal objective value equal to the reduced size of an ordered tree BDD with the given terminals.

**Proof Sketch.** Given the alternative interpretation of variables described in Section 5.2, we can conclude that the objective clauses in Eq. (16) aim to minimize the number of decision nodes that are not redundant and not equivalent to any of the nodes that come before them in the same level. Considering that equivalent decision node pairs can only appear in the same level, we can restate the objective as minimizing the number of decision nodes that remain when all of equivalent ones are merged and redundant ones are eliminated. According to Proposition 7, this objective will lead to a unique reduced BDD in its most compact form, proving that our encoding correctly models and minimizes the size of our solution. ◀

To decode the solution, we simply need to do as we did in Section 4.3, since the additional variables for size encoding are not involved in structuring the BDD. Once we decode and reduce the solution as before, we will end up with one decision node for every unsatisfied soft clause in Eq. (16).

6 MaxSAT Encoding of Multi-Dimensional BDDs

In Section 4 and Section 5, we consider learning BDDs by finding trees and reducing them to diagrams. However, the clauses in Eq.(8) are exponential in the number of splits and multiplied by the size of the dataset. Learning diagrams directly can help us consider more expressive solutions without increasing the number of splits, resulting in smaller encodings. In this section, we look at a family of diagrams that contain splits over multiple features at each level, i.e., multi-dimensional splits which themselves are BDDs with two labels.
**► Definition 9 (Directional Inner BDD).** Given a finite set of numerical features $F$ and a dimension $D$, a directional inner BDD $\mathcal{T} \leftrightarrow$ is a BDD operating on $D$ splits $\theta \leftrightarrow \in \{F \times \mathbb{R}\}^D$ with 0 and 1 as labels.

A multi-dimensional BDD operates on multi-dimensional splits. A multi-dimensional split uses a directional inner BDD to direct points towards left (label 1) and right (label 0), rather than using a single split.

**► Problem 10 (Multi-dimensional BDD Learning Problem).** Given a finite set of numerical features $F$, a finite set of labels $K$, a set of training points $x_i \in X$ with corresponding labels $y_i \in K$, a sequence of dimensions $\mathcal{D} = \{D_0, D_1, D_{s_{max}-1}\}$, and a number of multi-dimensional splits $s_{max}$, the goal is to find a sequence of directional inner BDDs $\theta^M = \{\mathcal{T}_0 \leftrightarrow, \mathcal{T}_1 \leftrightarrow, \dots, \mathcal{T}_{s_{max}-1} \leftrightarrow\}$ of respective dimensions $D_0, D_1, \dots, D_{s_{max}-1}$ with terminal labelling $\gamma^M : \{0, 1\}^{s_{max}} \to K$ to construct a BDD with the directional inner BDDs as splits. The objective for learning Multi-dimensional BDDs is high accuracy.

A multi-dimensional split is a generalization of an ordinary split, which makes a multi-dimensional BDD more expressive compared to an ordinary BDD with the same number of splits. Previous works have considered other representations that are designed to be more expressive than BDDs, such as sentential decision diagrams [12] or read-k-times branching programs [6, 22]. Next, we formally compare the expressiveness of BDDs against their multi-dimensional variants in the other direction.

**► Theorem 11.** Consider a multi-dimensional BDD $\mathcal{T}^M$ operating on $s_{max}$ directional inner BDDs $\mathcal{T}_s \leftrightarrow$ with split sequences $\theta_s \leftrightarrow$ and terminal labellings $\gamma_s \leftrightarrow$, which itself has a terminal labelling $\gamma^M$. The binary function that $\mathcal{T}^M$ represents is the same as the binary function represented by a BDD $\mathcal{T}$ with split sequence $\theta(\#(s,h)) = \theta_s^\leftrightarrow(h)$ and terminal labelling $\gamma(t) = \gamma^M(\gamma_0^\leftrightarrow(t_0)\gamma_1^\leftrightarrow(t_1)\dots\gamma_{s_{max}-1}^\leftrightarrow(t_{s_{max}-1}))$ where $\#(s,h)$ produces a complete ordering of $S_H$ by concatenating the split sequences of each directional inner BDD, and $t_0, t_1, \dots, t_{s_{max}-1}$ are sub-sequences of $t$ divided based on the sequence of dimensions.

**Proof Sketch.** The BDD using multi-dimensional splits can be transformed into a tree BDD with ordinary splits. In the transformation, every split of dimension $D$ will correspond to $D$ levels where nodes are expanded $2^D$-fold. With each expansion, we add annotations to the nodes specifying the label which they were assigned to by the corresponding directional inner BDD. Once all multi-dimensional splits are added, the final level of nodes are considered as terminals and are labelled according to how $\gamma^M$ labels their annotations. Note that the resulting tree can be turned into a diagram by merging all of the nodes that have equal annotations and replacing all of redundant nodes with their children. ◀

We can conclude from Theorem 11, that multi-dimensional BDDs operating on dimensions $\{D_0, D_1, \dots, D_{s_{max}-1}\}$ with $D_{total} = \sum_{i<s_{max}} D_i$ are less expressive than BDDs operating on $D_{total}$ splits. Combining with the aforementioned fact that multi-dimensional BDDs operating on $\{D_0, D_1, \dots, D_{s_{max}-1}\}$ dimensions are trivially more expressive than BDDs operating on $s_{max}$ splits, we have shown a lower and upper bound number of BDD splits when discussing the expressiveness of a multi-dimensional BDD.

Figure 3 shows how a multi-dimensional BDD can operate on two directional inner BDDs and how it can be unfolded to operate on ordinary splits.
<figure><img src="image_315.png" /><figcaption>Figure 3 Multi-dimensional BDD $\mathcal{T}^M$ operating on directional inner BDDs $\mathcal{T}_1 \leftrightarrow$ and $\mathcal{T}_2 \leftrightarrow$ and its unfolded version.</figcaption></figure>

## 6.1 Variables

In order to encode a multi-dimensional BDD, we use the same set of variables as Section 4.1 but remove $a_{s,j}$ variables since multi-dimensional splits require more features to be selected at each split. Furthermore, we add the following variables to encode the inner-workings of each directional inner BDD.

- $\hat{a}_{(s,h),j}$: The feature chosen at split $h$ of directional inner BDD $s$ is or comes before $j$.

- $\hat{d}_{(s,h),i}$: Point $x_i$ is directed to left at split $h$ of directional inner BDD $s$.

- $\hat{c}_{s,t}$: Terminal $t$ of directional inner BDD $s$ is assigned the label 1.

## 6.2 Clauses

We discard clauses Eqs. (1-6) from the original encoding since splits are learned differently than before. We keep the clauses Eqs. (6-9) however, since we still need the labels and the appearance of points at terminals to be modelled correctly. Lastly, we add CNF clauses for modelling multi-dimensional splits, in order to complete our encoding of a multi-dimensional BDD.

The clauses in Eqs. (17-18) guarantee that one feature is selected at each split of each directional inner BDD by enforcing the ordered encoding of $\hat{a}_{(s,h),j}$ variables. The clauses in Eqs. (19-21) guarantee that each split of each directional inner BDD conforms to a pairing of selected feature and threshold. The clauses in Eqs. (22-23) guarantee that the direction of a point in a multi-dimensional split matches the label of its containing leaf in the corresponding directional inner BDD. The clauses in Eq. (24) guarantee the leftmost leaf in each directional inner BDD to be assigned the 0 label (analogous to Eq. (5), see Shati et al. [33, 34] for further discussion). The set $S_H$ contains all pairs of directional inner BDDs and their corresponding splits $S_H = \{(s, h) | s < s_{max}, h < D_s\}$, $N_T^s$ contains the terminals of directional inner BDD $s$, $A_R(s, t)$ ($A_L(s, t)$) contains the right (left) ancestors of terminal $t$ within directional inner BDD $s$, and $\hat{a}_{(s,h),|F|}$ is set to the false constant.

$$ (\hat{a}_{(s,h),j}, \neg\hat{a}_{(s,h),j+1}) \qquad (s, h) \in S_H, j \in F \quad (17) $$

$$ (\hat{a}_{(s,h),0}) \qquad (s, h) \in S_H \quad (18) $$

$$ (\neg\hat{a}_{(s,h),j}, \hat{a}_{(s,h),j+1}, \hat{d}_{(s,h),i_1}, \neg\hat{d}_{(s,h),i_2}) \qquad (s, h) \in S_H, j \in F, (i_1, i_2) \in O_j(X) \quad (19) $$

$$ (\neg\hat{a}_{(s,h),j}, \hat{a}_{(s,h),j+1}, \neg\hat{d}_{(s,h),i_1}, \hat{d}_{(s,h),i_2}) \qquad (s, h) \in S_H, j \in F, (i_1, i_2) \in O_j^-(X) \quad (20) $$

$$ (\neg\hat{a}_{(s,h),j}, \hat{a}_{(s,h),j+1}, \hat{d}_{(s,h),\#_j}) \qquad (s, h) \in S_H, j \in F \quad (21) $$

$$ (\bigvee_{h \in A_R(s,t)} \hat{d}_{(s,h),i}, \bigvee_{h \in A_L(s,t)} \neg\hat{d}_{(s,h),i}, d_{s,i}, \neg\hat{c}_{s,t}) \qquad s \in S, x_i \in X, t \in \mathcal{N}_T^s \quad (22) $$
$$ ( \bigvee_{h \in A_R(s,t)} \hat{d}_{(s,h),i}, \bigvee_{h \in A_L(s,t)} \neg \hat{d}_{(s,h),i}, \neg d_{s,i}, \hat{c}_{s,t} ) \qquad s \in S, x_i \in X, t \in \mathcal{N}_T^s \quad (23) $$

$$ (\hat{c}_{s,0}) \qquad s \in S \quad (24) $$

## 6.3 Decoding

A satisfying assignment to the variables in Section 6.1 with regard to the clauses in Section 6.2, can be decoded into a tree BDD operating on multi-dimensional splits. In order to further transform the results into a BDD operating on 1-dimensional splits, we unfold the multi-dimensional splits and multiply the terminals according to the proof sketch for Theorem 11. Note that the number of splits in the resulting BDD is equal to the sum of dimensions from the multi-dimensional BDD.

The size optimization presented in Section 5 can also be utilized as a second stage after the solution is decoded. Specifically, we consider the multiplied sequence of terminal labels as input and treat empty terminals as before. In order to respect the structure of a multi-dimensional BDD, we also need to have additional clauses guaranteeing that multiplied instances of the same original terminal have the same label. Note that the size encoding cannot be employed in a 1-stage approach since the structure of the BDD is not yet determined at first.

# 7 Experiments

In this section, we perform studies to experimentally analyze the performance of our tree (Section 4), size (Section 5), and multi-dimensional (Section 6) encodings for learning BDDs. We compare the performance of our approach against state-of-the-art BDD learning baseline and investigate the trade-off between size, accuracy, and performance in 1-stage, 2-stage, and multi-dimensional approaches to learning compact diagrams. Throughout the experiments, we seek to find out whether compactness can be achieved without significant compromise, and investigate its impact on testing accuracy. Furthermore, we aim to understand if the added expressiveness of multi-dimensional BDDs allows us to achieve high quality solutions with lower number of splits.

## 7.1 Setup

We use the Java programming language to produce encodings and the Loandra MaxSAT solver [4] to solve each instance. We set the solver timeout limit to 15 minutes and use the best found solution in case of a timeout. Our experiments are run on an AMD EPYC 7502 32-core processor and 256GB of RAM.

## 7.2 Baseline

We compare our approach against the recent work on SAT-based learning of BDDs. Hu et al. [16] aims to find BDD classifiers similar to our approach. However, unlike our approach, they only support binary classification and require pre-processing of each numeric feature into a set of binary features. Furthermore, they do not explicitly model the size of the reduced BDD.

Note that Hu et al. [16, 17] have compared optimal BDDs against optimal decision trees and heuristic decision trees and found that BDDs are competitive in terms of testing accuracy and have smaller size. Our approach aims to improve the performance of Hu et al. and learn more compact BDD classifiers with comparable accuracy. We therefore compare our approach to Hu et al. in terms of runtime, accuracy, and size of BDDs.
## 7.3 Datasets

We run our experiments over a range of datasets from the UCI repository [13] covering different number of labels, both numerical and binary features, and different dataset sizes.

## 7.4 Results

**Results on Comparing Our Base Encoding Against Hu et al. [16]**

In our first set of results, we compare the performance of our approach in terms of runtime and solution quality against Hu et al. for different numbers of splits. Note that while our approach learns multi-terminal BDDs, Hu et al.'s approach learns standard BDDs and is unable to support datasets with more than two labels. Furthermore, since the objective in Hu et al. does not include compactness considerations, we use our base encoding that only optimizes accuracy as well (Section 4). As both approaches explore the same space of (feasible and) optimal BDD solutions, we focus our comparison on optimization performance (i.e., training accuracy and runtime). In contrast, in the next two sets of experiments, we will also evaluate the testing accuracy over 5-fold cross-validation.

Based on the results presented in Table 1, we see that our approach is able to achieve a higher than or equal to Hu et al. [16] accuracy in all but one case. Furthermore, the runtimes of non-timeout cases show that our approach can also prove optimality much faster. As we expected given our direct encoding of non-binary values, the improvement is most noticeable in datasets with highly numerical features, namely Banknote and Ionosphere.

**Results on 1-Stage and 2-Stage Size Optimization**

Next, we evaluate the encoding presented in Section 5 to minimize the size of our learned BDDs. We perform size optimization in 1-stage and 2-stage approaches. In the 1-stage approach, we use different values for the weight of the size objective against the accuracy objective. Given a weight $\beta$ and $s_{max}$ splits, the combined size and accuracy objective aims to find a solution $f$ with maximum $\beta(acc(f)) - (2^{s_{max}} - 1)|N_D^f|$, where $N_D^f$ is the set of decision nodes in the reduced version of $f$. We use $\beta = \infty^2$ to denote the 2-stage approach, $\beta = \infty^1$ to denote that accuracy is completely prioritized over size in the 1-stage approach, and Def. to denote the approach without any size optimization. Note that even in the 1-stage approach, running the second stage can further optimize the size if the first stage have yielded a sub-optimal solution (timeout). However, we opt against mixing the two approaches for the sake of clarity in comparison between the two. We use 5-fold cross validation and report the average value across folds to understand the effects of compactness on generalization and testing accuracy.

The results are presented in Figure 4. As expected, we see increase in training accuracy and size as we shift the priority to accuracy by changing the $\beta$ value from 1 to $\infty^2$. However, the improvement in accuracy stagnates while the size continues to grow, indicating that a large enough $\beta$ value can act as complete prioritization of accuracy. Interestingly, testing accuracy stagnates sooner and suffers from a larger variability as $\beta$ increases, demonstrating a limited but positive effect of compactness on testing accuracy and generalization.

The solutions for lower $\beta$ values in the 1-stage approaches are significantly smaller than the $\beta = \infty^2$ case, highlighting the importance of adding size considerations while the solution is being learned. However, optimizing size as a second stage still provides a significant size reduction compared to the case with no size optimization. For a detailed report of the results per dataset, we refer the reader to Table 2 in Appendix A.1.
<table><caption>Table 1 Results for learning max-accuracy BDDs.</caption><thead><tr><th rowspan="2">Dataset</th><th rowspan="2">Splits</th><th colspan="2">Accuracy (%)</th><th colspan="2">Time (s)</th></tr><tr><th>Ours</th><th>Hu et al. [16]</th><th>Ours</th><th>Hu et al. [16]</th></tr></thead><tbody><tr><td rowspan="3">Banknote</td><td>4</td><td>96.8</td><td>96.8</td><td>843.01</td><td>TO</td></tr><tr><td>5</td><td>97.6</td><td>90.4</td><td>TO</td><td>TO</td></tr><tr><td>6</td><td>98.5</td><td>91.8</td><td>TO</td><td>TO</td></tr><tr><td rowspan="3">Breast</td><td>4</td><td>89.7</td><td>89.7</td><td>TO</td><td>TO</td></tr><tr><td>5</td><td>92.2</td><td>91.4</td><td>TO</td><td>TO</td></tr><tr><td>6</td><td>94.8</td><td>94.8</td><td>TO</td><td>TO</td></tr><tr><td rowspan="3">Cryotherapy</td><td>4</td><td>97.8</td><td>97.8</td><td>1.23</td><td>2.43</td></tr><tr><td>5</td><td>98.9</td><td>98.9</td><td>3.97</td><td>9.41</td></tr><tr><td>6</td><td>100</td><td>100</td><td>0.44</td><td>1.64</td></tr><tr><td rowspan="3">Immunotherapy</td><td>4</td><td>95.6</td><td>95.6</td><td>8.18</td><td>26.65</td></tr><tr><td>5</td><td>96.7</td><td>96.7</td><td>74.08</td><td>290.99</td></tr><tr><td>6</td><td>97.8</td><td>97.8</td><td>433.35</td><td>TO</td></tr><tr><td rowspan="3">Ionosphere</td><td>4</td><td>94.9</td><td>90.6</td><td>TO</td><td>TO</td></tr><tr><td>5</td><td>95.2</td><td>85.8</td><td>TO</td><td>TO</td></tr><tr><td>6</td><td>96.6</td><td>90.6</td><td>TO</td><td>TO</td></tr><tr><td rowspan="3">Iris</td><td>4</td><td>98.7</td><td>-</td><td>0.67</td><td>-</td></tr><tr><td>5</td><td>99.3</td><td>-</td><td>0.62</td><td>-</td></tr><tr><td>6</td><td>100</td><td>-</td><td>0.43</td><td>-</td></tr><tr><td rowspan="3">User</td><td>4</td><td>94.2</td><td>-</td><td>54.57</td><td>-</td></tr><tr><td>5</td><td>95.7</td><td>-</td><td>828.34</td><td>-</td></tr><tr><td>6</td><td>97.7</td><td>-</td><td>TO</td><td>-</td></tr><tr><td rowspan="3">Vertebral</td><td>4</td><td>88.1</td><td>87.7</td><td>TO</td><td>TO</td></tr><tr><td>5</td><td>89.7</td><td>90</td><td>TO</td><td>TO</td></tr><tr><td>6</td><td>91</td><td>90.3</td><td>TO</td><td>TO</td></tr><tr><td rowspan="3">Wine</td><td>4</td><td>99.4</td><td>-</td><td>29.29</td><td>-</td></tr><tr><td>5</td><td>100</td><td>-</td><td>2.07</td><td>-</td></tr><tr><td>6</td><td>100</td><td>-</td><td>1.14</td><td>-</td></tr><tr><td rowspan="3">Car</td><td>4</td><td>92.5</td><td>92.5</td><td>316.88</td><td>TO</td></tr><tr><td>5</td><td>92.9</td><td>92.9</td><td>TO</td><td>TO</td></tr><tr><td>6</td><td>95.4</td><td>95.4</td><td>TO</td><td>TO</td></tr><tr><td rowspan="3">Monk2</td><td>4</td><td>74.6</td><td>74.6</td><td>85.49</td><td>473.43</td></tr><tr><td>5</td><td>84.6</td><td>84.6</td><td>185.3</td><td>416.37</td></tr><tr><td>6</td><td>100</td><td>100</td><td>0.35</td><td>1.09</td></tr></tbody></table>

## Results on Learning Multi-Dimensional BDDs

In our next set of experiments, we evaluate our approach for directly learning multi-dimensional BDDs. Our goal is to understand whether the added expressiveness of multi-dimensional BDDs allows us to find high quality solutions with a lower number of splits. Furthermore, we aim to study the effects of dimension division for multi-dimensional BDDs
<figure><img src="image_316.png" /><figcaption>Figure 4 Distributional result of size, training accuracy, and testing accuracy across datasets for different $\beta$ values.</figcaption></figure>

by considering one balanced diagram with three multi-dimensional splits (2-2-2) and one unbalanced diagram with one 3-dimensional split followed by three 1-dimensional splits (3-1-1-1). Finally, we have used the size encoding in Section 5 to optimally decide the labels of empty terminals towards compactness in a second stage. Note that the size of a multi-dimensional BDD is considered to be its number of decision nodes ($|N_D^f|$) after it is unfolded and reduced.

<figure><img src="image_317.png" /><figcaption>Figure 5 Distributional result of size, accuracy, number of clauses, and average clause length for BDDs of two number of splits and Multi-dimensional BDDs of two dimensions.</figcaption></figure>

The results of the third set of experiments presented in Figure 5 show the two multi-dimensional approaches achieve training accuracy that is higher than a one-dimensional diagram with 5 splits and lower than a one-dimensional diagram with 6 splits. Since the sum of dimensions equal 6 both case, the upper bound is expected according to Theorem 11. However, the comparison against 5 splits, shows that the multi-dimensional approaches
perform better than the theoretical lower bound (resp. 3 and 4 splits) in practice. The same comparison is observed in testing accuracy and size. However, the 3-1-1-1 approach is able to achieve smaller variability in testing accuracy while still being close to the best approach in size. Finally, we see a lower number of clauses and a significantly shorter average clause length for our multi-dimensional approaches against ordinary BDDs of similar expressiveness. We refer the reader to Table 3 in Appendix A.1 for the complete results of this experiment with a larger set of dimension sequences considered.

Next, we run experiments for a significantly larger dataset, namely the Adult dataset ($|X| = 32561$, $|F| = 105$, $|K| = 2$), to investigate the difference in encoding size and performance between one-dimensional and multi-dimensional approaches on large datasets. Given the scale of the tasks, we increase the timeout limit to 60 minutes. We compare, one-dimensional BDDs against ones with 2, 3, and 4-dimensional splits according to their total number of dimensions in Figure 6 and Figure 7. Figure 7 depicts the exponential growth of encoding size for the one-dimensional approach. The training accuracy presented in Figure 6 shows that the exponential size causes the ordinary approach to decrease in quality and finally cease to produce any solutions due to memory overflow. We observe that by using two dimensional splits, we maintain a higher accuracy over a significantly larger number of dimensions. Three-dimensional and four-dimensional splits prove to be more challenging compared to two-dimensional splits. However, we are still able to find solutions for higher total number of dimensions compared to the one-dimensional splits. In this experiment, we avoid 5-fold cross-validation as we focus on the optimization performance of the two methods.

<figure><img src="image_318.png" /><figcaption>Figure 6 Training accuracy for the Adult dataset for different dimensions and number of total splits.</figcaption></figure>

<figure><img src="image_319.png" /><figcaption>Figure 7 Number and average length of clauses of the encoding for the Adult dataset for different dimensions and number of total splits.</figcaption></figure>

# 8 Conclusion

In this paper, we present a novel MaxSAT encoding for learning Binary Decision Diagrams. Our BDD encoding represents a tree which can be reduced to an equivalent diagram. We extend our encoding with optimization for the size of the reduced diagram. The size objective can be balanced against accuracy in a 1-stage approach or optimized as a second stage. Furthermore, we present a variant of our encoding using multi-dimensional splits, which are inner BDDs themselves. Our experiments show that we outperform the state-of-the-art SAT-based BDD learning baseline due to our direct encoding of numerical splits. We further show that our 1-stage and 2-stage size optimization approaches lead to significantly more compact solutions while maintaining testing accuracy. Finally, we show that the expressiveness of our multi-dimensional BDDs allows us to produce high quality solutions in smaller number of splits, mitigating the exponential growth in the size of the encoding.
33:16 SAT-Based Learning of Compact Binary Decision Diagrams for Classification

Our work, can be extended in a number of ways. Our size encoding for the 2-stage optimization can be extended by allowing more than empty terminal labels to be changed, e.g., feature ordering or splits. Moreover, nested directional BDDs can be added to our multi-dimensional BDD encoding to improve expressiveness and avoid exponentiation even further. Other interesting directions for future work involve investigating different strategies for balancing the two objectives, e.g., producing Pareto optimal solutions, or a more comprehensive analysis of parameters including but not limited to the dimension sequences. Finally, investigating the impact of more expressive BDDs such as free BDDs, which are not constrained to be ordered [32], is also an interesting direction for future research.

References

1 Gaël Aglin, Siegfried Nijssen, and Pierre Schaus. Learning optimal decision trees using caching branch-and-bound search. In AAAI Conference on Artificial Intelligence (AAAI), pages 3146–3153, 2020.

2 Sheldon B. Akers. Binary decision diagrams. IEEE Transactions on computers, 27(06):509–516, 1978.

3 Florent Avellaneda. Efficient inference of optimal decision trees. In AAAI Conference on Artificial Intelligence (AAAI), pages 3195–3202, 2020.

4 Jeremias Berg, Emir Demirović, and Peter J Stuckey. Core-boosted linear search for incomplete MaxSAT. In International Conference on Integration of Constraint Programming, Artificial Intelligence, and Operations Research (CPAIOR), pages 39–56. Springer, 2019.

5 Dimitris Bertsimas and Jack Dunn. Optimal classification trees. Machine Learning, 106(7):1039–1082, 2017.

6 Allan Borodin, Alexander Razborov, and Roman Smolensky. On lower bounds for read-k-times branching programs. Computational Complexity, 3:1–18, 1993.

7 Leo Breiman, Jerome H Friedman, Richard A Olshen, and Charles J Stone. Classification and regression trees. Wadsworth & Brooks/Cole Advanced Books & Software, 1984.

8 Randal E Bryant. Graph-based algorithms for boolean function manipulation. Computers, IEEE Transactions on, 100(8):677–691, 1986.

9 Gianpiero Cabodi, Paolo E Camurati, Alexey Ignatiev, Joao Marques-Silva, Marco Palena, and Paolo Pasini. Optimizing binary decision diagrams for interpretable machine learning classification. In 2021 Design, Automation & Test in Europe Conference & Exhibition (DATE), pages 1122–1125. IEEE, 2021.

10 Diogo V Carvalho, Eduardo M Pereira, and Jaime S Cardoso. Machine learning interpretability: A survey on methods and metrics. Electronics, 8(8):832, 2019.

11 Edmund M Clarke, Masahiro Fujita, and Xudong Zhao. Multi-terminal binary decision diagrams and hybrid decision diagrams. Representations of discrete functions, pages 93–108, 1996.

12 Adnan Darwiche. Sdd: A new canonical representation of propositional knowledge bases. In Twenty-Second International Joint Conference on Artificial Intelligence, 2011.

13 Dheeru Dua and Casey Graff. UCI machine learning repository, 2017. URL: http://archive.ics.uci.edu/ml.

14 Alexandre M Florio, Pedro Martins, Maximilian Schiffer, Thiago Serra, and Thibaut Vidal. Optimal decision diagrams for classification. arXiv preprint arXiv:2205.14500, 2022.

15 Oktay Günlük, Jayant Kalagnanam, Minhan Li, Matt Menickelly, and Katya Scheinberg. Optimal decision trees for categorical data via integer programming. Journal of Global Optimization, pages 1–28, 2021.

16 Hao Hu, Marie-José Huguet, and Mohamed Siala. Optimizing binary decision diagrams with maxsat for classification. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 36(4), pages 3767–3775, 2022.
17 Hao Hu, Marie-José Huguet, and Mohamed Siala. Optimizing binary decision diagrams with maxsat for classification, 2022. arXiv:2203.11386.

18 Hao Hu, Mohamed Siala, Emmanuel Hébrard, and Marie-José Huguet. Learning optimal decision trees with MaxSAT and its integration in AdaBoost. In *International Joint Conference on Artificial Intelligence and Pacific Rim International Conference on Artificial Intelligence (IJCAI-PRICAI)*, 2020.

19 Alexey Ignatiev and Joao Marques-Silva. SAT-based rigorous explanations for decision lists. In *Theory and Applications of Satisfiability Testing-SAT 2021: 24th International Conference, Barcelona, Spain, July 5-9, 2021, Proceedings 24*, pages 251-269. Springer, 2021.

20 Alexey Ignatiev, Filipe Pereira, Nina Narodytska, and Joao Marques-Silva. A SAT-based approach to learn explainable decision sets. In *Automated Reasoning: 9th International Joint Conference, IJCAR 2018, Held as Part of the Federated Logic Conference, FloC 2018, Oxford, UK, July 14-17, 2018, Proceedings 9*, pages 627-645. Springer, 2018.

21 Dmitry Ignatov and Andrey Ignatov. Decision stream: Cultivating deep decision trees. In *2017 IEEE 29th international conference on tools with artificial intelligence (ictai)*, pages 905-912.

22 Stasys Jukna. A note on read-times branching programs. *RAIRO-Theoretical Informatics and Applications*, 29(1):75-83, 1995.

23 Donald E Knuth. *The Art of Computer Programming, Volume 4, Fascicle 1: Bitwise Tricks & Techniques; Binary Decision Diagrams*. Addison-Wesley Professional, 2009.

24 Ron Kohavi. Bottom-up induction of oblivious read-once decision graphs. In *Machine Learning: ECML-94: European Conference on Machine Learning Catania, Italy, April 6-8, 1994 Proceedings 7*, pages 154-169. Springer, 1994.

25 Benjamin Letham, Cynthia Rudin, Tyler H McCormick, and David Madigan. Interpretable classifiers using rules and bayesian analysis: Building a better stroke prediction model, 2015.

26 Breiman LI, Jerome Friedman, RA Olshen, and C.J. Stone. Classification and regression trees (CART). *Biometrics*, 40:358, September 1984. doi:10.2307/2530946.

27 Joao Marques-Silva and Alexey Ignatiev. Delivering trustworthy AI through formal XAI. In *Proceedings of the AAAI Conference on Artificial Intelligence*, volume 36(11), pages 12342-12350, 2022.

28 Bernard ME Moret. Decision trees and diagrams. *ACM Computing Surveys (CSUR)*, 14(4):593-623, 1982.

29 J Ross Quinlan. Induction of decision trees. *Machine learning*, 1(1):81-106, 1986.

30 J Ross Quinlan. *C4.5: programs for machine learning*. Elsevier, 2014.

31 Cynthia Rudin. Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead. nat mach intell 1: 206-215. DOI: https://doi.org/10.1038/s42256-019-0048-x, 2019.

32 Petr Savický and Ingo Wegener. Efficient algorithms for the transformation between different types of binary decision diagrams. *Acta Informatica*, 34(4):245-256, 1997.

33 Pouya Shati, Eldan Cohen, and Sheila McIlraith. SAT-based approach for learning optimal decision trees with non-binary features. In *27th International Conference on Principles and Practice of Constraint Programming (CP 2021)*. Schloss Dagstuhl-Leibniz-Zentrum für Informatik, 2021.

34 Pouya Shati, Eldan Cohen, and Sheila A. McIlraith. SAT-based optimal classification trees for non-binary data. *Constraints*, July 2023. doi:10.1007/s10601-023-09348-1.

35 Jamie Shotton, Toby Sharp, Pushmeet Kohli, Sebastian Nowozin, John Winn, and Antonio Criminisi. Decision jungles: Compact and rich models for classification. *Advances in neural information processing systems*, 26, 2013.

36 Hazem Torfah, Shetal Shah, Supratik Chakraborty, S Akshay, and Sanjit A Seshia. Synthesizing pareto-optimal interpretations for black-box models. In *2021 Formal Methods in Computer Aided Design (FMCAD)*, pages 153-162. IEEE, 2021.
37 Hélene Verhaeghe, Siegfried Nijssen, Gilles Pesant, Claude-Guy Quimper, and Pierre Schaus. Learning optimal decision trees using constraint programming. *Constraints*, 25(3):226–250, 2020.

38 Yu Zhang, Peter Tiňo, Aleš Leonardis, and Ke Tang. A survey on neural network interpretability. *IEEE Transactions on Emerging Topics in Computational Intelligence*, 5(5):726–742, 2021.

# A Appendix

## A.1 Additional Experimental Results

<table><caption>Table 2 Average 5-fold results for learning BDDs with the combined objective of accuracy and compactness.</caption><thead><tr><th rowspan="2">D.S.</th><th rowspan="2">Splits</th><th colspan="5">Size</th><th colspan="5">Training Acc. (%)</th><th colspan="5">Testing Acc. (%)</th></tr><tr><th>1</th><th>3</th><th>10</th><th>∞<sup>1</sup></th><th>∞<sup>2</sup></th><th>1</th><th>3</th><th>10</th><th>∞<sup>1</sup></th><th>∞<sup>2</sup></th><th>1</th><th>3</th><th>10</th><th>∞<sup>1</sup></th><th>∞<sup>2</sup></th></tr></thead><tbody><tr><td rowspan="3">Banknote</td><td>4</td><td>1</td><td>3.6</td><td>4</td><td>4</td><td>4.4</td><td>85.3</td><td>95.8</td><td>96.8</td><td>96.8</td><td>96.8</td><td>85.2</td><td>95</td><td>96.6</td><td>96.6</td><td>96.4</td></tr><tr><td>5</td><td>3</td><td>4</td><td>5</td><td>11</td><td>6</td><td>94</td><td>96.8</td><td>97.6</td><td>97.5</td><td>97.6</td><td>93.7</td><td>96.4</td><td>96.9</td><td>96.8</td><td>96.7</td></tr><tr><td>6</td><td>4</td><td>6</td><td>5.8</td><td>19.8</td><td>6.4</td><td>96.8</td><td>98.6</td><td>98.4</td><td>98.6</td><td>98.6</td><td>96.5</td><td>97.4</td><td>97.4</td><td>97.2</td><td>97.2</td></tr><tr><td rowspan="3">Breast</td><td>4</td><td>1</td><td>4</td><td>5.2</td><td>6.4</td><td>6.8</td><td>72.4</td><td>88.8</td><td>90.7</td><td>90.7</td><td>90.7</td><td>71.5</td><td>75.9</td><td>75.9</td><td>78.4</td><td>77.6</td></tr><tr><td>5</td><td>3.8</td><td>5.4</td><td>11.4</td><td>11.8</td><td>11</td><td>88.8</td><td>91.8</td><td>94</td><td>94</td><td>94</td><td>74.1</td><td>70.7</td><td>74.1</td><td>72.4</td><td>72.5</td></tr><tr><td>6</td><td>5.6</td><td>9.4</td><td>22.2</td><td>22.4</td><td>16</td><td>92.9</td><td>96.1</td><td>97</td><td>97.2</td><td>97.2</td><td>75</td><td>71.6</td><td>63</td><td>71.6</td><td>65.5</td></tr><tr><td rowspan="3">Cryo.</td><td>4</td><td>1</td><td>3</td><td>4.8</td><td>4.8</td><td>6</td><td>86.1</td><td>94.7</td><td>97.8</td><td>97.8</td><td>97.8</td><td>75.6</td><td>88.9</td><td>92.2</td><td>92.2</td><td>94.4</td></tr><tr><td>5</td><td>2.4</td><td>5.2</td><td>5.8</td><td>5.8</td><td>8</td><td>93.1</td><td>99.2</td><td>99.4</td><td>99.4</td><td>99.4</td><td>83.3</td><td>94.4</td><td>91.1</td><td>91.1</td><td>90</td></tr><tr><td>6</td><td>4.6</td><td>5.8</td><td>5.8</td><td>5.8</td><td>10.4</td><td>98.3</td><td>100</td><td>100</td><td>100</td><td>100</td><td>95.6</td><td>87.8</td><td>90</td><td>90</td><td>82.2</td></tr><tr><td rowspan="3">Immuno.</td><td>4</td><td>1</td><td>2.2</td><td>5.4</td><td>5.4</td><td>6.8</td><td>87.2</td><td>91.4</td><td>95.6</td><td>95.6</td><td>95.6</td><td>82.2</td><td>83.3</td><td>83.3</td><td>83.3</td><td>85.6</td></tr><tr><td>5</td><td>1.6</td><td>4.4</td><td>5.4</td><td>5.4</td><td>8.4</td><td>89.7</td><td>95.3</td><td>96.7</td><td>96.7</td><td>96.7</td><td>81.1</td><td>81.1</td><td>78.9</td><td>81.1</td><td>73.3</td></tr><tr><td>6</td><td>3.8</td><td>7.2</td><td>10</td><td>10</td><td>14</td><td>94.7</td><td>98.3</td><td>99.2</td><td>99.2</td><td>99.2</td><td>80</td><td>76.7</td><td>73.3</td><td>72.2</td><td>75.6</td></tr><tr><td rowspan="3">Iono.</td><td>4</td><td>1.8</td><td>2.4</td><td>4</td><td>5.6</td><td>4</td><td>90</td><td>92</td><td>95.2</td><td>94.9</td><td>95.2</td><td>87.2</td><td>90.9</td><td>90.9</td><td>90.3</td><td>92.3</td></tr><tr><td>5</td><td>2</td><td>3.6</td><td>5.2</td><td>10.2</td><td>8.4</td><td>91.2</td><td>94.4</td><td>95.9</td><td>95.7</td><td>95.9</td><td>89.7</td><td>91.7</td><td>89.5</td><td>85.8</td><td>89.2</td></tr><tr><td>6</td><td>2.4</td><td>6</td><td>7</td><td>17.2</td><td>12.8</td><td>92</td><td>96.3</td><td>96.9</td><td>97.2</td><td>97.1</td><td>90.9</td><td>88.6</td><td>88</td><td>86.9</td><td>86.3</td></tr><tr><td rowspan="3">Iris</td><td>4</td><td>2</td><td>2.4</td><td>3.8</td><td>3.8</td><td>4.2</td><td>96.3</td><td>97.3</td><td>98.8</td><td>98.8</td><td>98.8</td><td>94</td><td>94</td><td>94.7</td><td>94.7</td><td>94</td></tr><tr><td>5</td><td>2</td><td>2.8</td><td>5.2</td><td>5.2</td><td>6</td><td>96.3</td><td>98</td><td>99.5</td><td>99.5</td><td>99.5</td><td>94</td><td>95.3</td><td>96.7</td><td>96</td><td>96</td></tr><tr><td>6</td><td>2.8</td><td>5.4</td><td>5.8</td><td>5.8</td><td>7</td><td>98</td><td>99.8</td><td>100</td><td>100</td><td>100</td><td>95.3</td><td>95.3</td><td>94.7</td><td>96</td><td>94.7</td></tr><tr><td rowspan="3">User</td><td>4</td><td>2.2</td><td>4.4</td><td>6</td><td>6</td><td>6</td><td>82.8</td><td>91.5</td><td>94.3</td><td>94.3</td><td>94.3</td><td>79.8</td><td>87.6</td><td>93.4</td><td>93.4</td><td>93.4</td></tr><tr><td>5</td><td>4</td><td>6</td><td>6.2</td><td>9.2</td><td>7</td><td>90.5</td><td>95.9</td><td>96</td><td>96</td><td>96</td><td>86</td><td>93</td><td>91.9</td><td>93.8</td><td>92.2</td></tr><tr><td>6</td><td>6</td><td>7</td><td>8</td><td>31.6</td><td>9.2</td><td>96.2</td><td>97.2</td><td>97.8</td><td>97.8</td><td>97.8</td><td>93</td><td>93</td><td>95.7</td><td>95.7</td><td>95.4</td></tr><tr><td rowspan="3">Vertebral</td><td>4</td><td>1.2</td><td>2.4</td><td>4</td><td>5.4</td><td>4.6</td><td>81.8</td><td>86.8</td><td>89.8</td><td>89.8</td><td>89.8</td><td>76.8</td><td>76.5</td><td>79</td><td>78.1</td><td>77.7</td></tr><tr><td>5</td><td>2.2</td><td>4.2</td><td>5.2</td><td>10.6</td><td>8.6</td><td>86.1</td><td>89.4</td><td>89.8</td><td>90.8</td><td>91</td><td>74.8</td><td>81.6</td><td>81.9</td><td>79.7</td><td>76.1</td></tr><tr><td>6</td><td>3.4</td><td>6.2</td><td>10.2</td><td>20.6</td><td>15</td><td>88.5</td><td>91.6</td><td>92.4</td><td>92.8</td><td>92.7</td><td>80.3</td><td>80.3</td><td>79</td><td>76.8</td><td>78.7</td></tr><tr><td rowspan="3">Wine</td><td>4</td><td>2.2</td><td>3</td><td>4</td><td>5</td><td>6.4</td><td>93.8</td><td>98.5</td><td>99.6</td><td>99.9</td><td>99.9</td><td>87.7</td><td>96.6</td><td>93.3</td><td>93.3</td><td>93.9</td></tr><tr><td>5</td><td>3</td><td>3.6</td><td>4.6</td><td>4.6</td><td>9.4</td><td>98.5</td><td>99.3</td><td>100</td><td>100</td><td>100</td><td>95</td><td>95</td><td>96.7</td><td>95.5</td><td>92.7</td></tr><tr><td>6</td><td>3</td><td>4.6</td><td>4.6</td><td>4.6</td><td>10</td><td>98.5</td><td>100</td><td>100</td><td>100</td><td>100</td><td>95</td><td>93.3</td><td>92.2</td><td>93.3</td><td>93.3</td></tr><tr><td rowspan="3">Car</td><td>4</td><td>2</td><td>4</td><td>4</td><td>4</td><td>5.4</td><td>85.5</td><td>92.5</td><td>92.5</td><td>92.5</td><td>92.5</td><td>85.5</td><td>92.5</td><td>92.5</td><td>92.5</td><td>92.5</td></tr><tr><td>5</td><td>3.2</td><td>4</td><td>4.6</td><td>6.8</td><td>7.4</td><td>90</td><td>92.5</td><td>92.9</td><td>93.1</td><td>93.1</td><td>88.6</td><td>92.5</td><td>91.7</td><td>92.1</td><td>92.1</td></tr><tr><td>6</td><td>4</td><td>6.6</td><td>6.8</td><td>12.2</td><td>9.6</td><td>89</td><td>95</td><td>95</td><td>95.4</td><td>95.4</td><td>87.6</td><td>93.8</td><td>93.9</td><td>94.7</td><td>94.7</td></tr><tr><td rowspan="3">Monk2</td><td>4</td><td>0</td><td>3.4</td><td>5.4</td><td>5.4</td><td>6.6</td><td>62.1</td><td>71.7</td><td>75.4</td><td>75.4</td><td>75.4</td><td>62.1</td><td>63.3</td><td>64.5</td><td>65.6</td><td>63.9</td></tr><tr><td>5</td><td>0.6</td><td>7.8</td><td>9</td><td>9</td><td>13</td><td>64.6</td><td>84.3</td><td>84.9</td><td>84.9</td><td>84.9</td><td>59.7</td><td>79.3</td><td>76.4</td><td>77.6</td><td>74.6</td></tr><tr><td>6</td><td>7.6</td><td>11.8</td><td>11.8</td><td>11.8</td><td>13.4</td><td>86.2</td><td>100</td><td>100</td><td>100</td><td>100</td><td>81.7</td><td>98.2</td><td>98.8</td><td>98.2</td><td>99.4</td></tr></tbody></table>
<table><caption>**Table 3** Results for learning max-accuracy multi-dimensional BDDs.</caption><thead><tr><th rowspan="2">Dataset</th><th rowspan="2">Dimensions</th><th colspan="2">Accuracy (%)</th><th colspan="2">Size</th><th rowspan="2">Time (s)</th></tr><tr><th>Training</th><th>Testing</th><th>Stage 1</th><th>Stage 2</th></tr></thead><tbody><tr><td rowspan="4">Banknote</td><td>2-2-2</td><td><strong>98.6</strong></td><td>97.5</td><td>7.4</td><td><strong>6.6</strong></td><td>TO</td></tr><tr><td>3-3</td><td>98.4</td><td>97.5</td><td><strong>6.6</strong></td><td><strong>6.6</strong></td><td>TO</td></tr><tr><td>1-2-3</td><td><strong>98.6</strong></td><td>96.9</td><td>8.2</td><td>6.8</td><td>TO</td></tr><tr><td>3-1-1-1</td><td><strong>98.6</strong></td><td><strong>97.7</strong></td><td>8.8</td><td>6.8</td><td>TO</td></tr><tr><td rowspan="4">Breast</td><td>2-2-2</td><td><strong>94.6</strong></td><td><strong>74.1</strong></td><td><strong>8.6</strong></td><td><strong>8.6</strong></td><td>TO</td></tr><tr><td>3-3</td><td>92.9</td><td>64.6</td><td>8.8</td><td>8.8</td><td>TO</td></tr><tr><td>1-2-3</td><td>94.2</td><td>68.1</td><td>9.6</td><td>9.6</td><td>TO</td></tr><tr><td>3-1-1-1</td><td>94.2</td><td>72.4</td><td>11</td><td>9.8</td><td>TO</td></tr><tr><td rowspan="4">Cryotherapy</td><td>2-2-2</td><td>99.7</td><td>86.7</td><td>10.6</td><td>10.2</td><td><strong>6.18</strong></td></tr><tr><td>3-3</td><td>99.4</td><td><strong>92.2</strong></td><td>10.6</td><td>10</td><td>25.29</td></tr><tr><td>1-2-3</td><td>99.7</td><td>90</td><td>10.6</td><td>9.8</td><td>13.66</td></tr><tr><td>3-1-1-1</td><td><strong>100</strong></td><td>86.7</td><td><strong>10</strong></td><td><strong>9.4</strong></td><td>6.45</td></tr><tr><td rowspan="4">Immunotherapy</td><td>2-2-2</td><td>97.2</td><td>75.6</td><td><strong>9.2</strong></td><td><strong>8.6</strong></td><td>515.83</td></tr><tr><td>3-3</td><td>96.9</td><td>76.7</td><td>10.2</td><td>10.2</td><td>750.07</td></tr><tr><td>1-2-3</td><td>97.2</td><td>76.7</td><td>11.6</td><td>11.2</td><td>767.69</td></tr><tr><td>3-1-1-1</td><td><strong>98.1</strong></td><td><strong>84.4</strong></td><td>9.8</td><td>9.6</td><td><strong>489.51</strong></td></tr><tr><td rowspan="4">Ionosphere</td><td>2-2-2</td><td>96.9</td><td>90</td><td><strong>6.6</strong></td><td><strong>6.6</strong></td><td>TO</td></tr><tr><td>3-3</td><td>95.4</td><td>88.3</td><td>8.8</td><td>8.8</td><td>TO</td></tr><tr><td>1-2-3</td><td><strong>97</strong></td><td>90.3</td><td>7.4</td><td>7.4</td><td>TO</td></tr><tr><td>3-1-1-1</td><td>96.7</td><td><strong>91.5</strong></td><td>8.2</td><td>8.2</td><td>TO</td></tr><tr><td rowspan="4">Iris</td><td>2-2-2</td><td>99.5</td><td><strong>95.3</strong></td><td><strong>10.4</strong></td><td>8.4</td><td>1.17</td></tr><tr><td>3-3</td><td>99.5</td><td>94.7</td><td>11.2</td><td>11.2</td><td>0.67</td></tr><tr><td>1-2-3</td><td>99.8</td><td>94</td><td>14.6</td><td><strong>8</strong></td><td>1.2</td></tr><tr><td>3-1-1-1</td><td><strong>100</strong></td><td>93.3</td><td>12.4</td><td>8.8</td><td><strong>0.44</strong></td></tr><tr><td rowspan="4">User</td><td>2-2-2</td><td>95.6</td><td>91.1</td><td><strong>12</strong></td><td>11.6</td><td>TO</td></tr><tr><td>3-3</td><td>93</td><td>88.7</td><td>12.2</td><td>12.2</td><td><strong>143.96</strong></td></tr><tr><td>1-2-3</td><td><strong>97.7</strong></td><td><strong>96.1</strong></td><td>18.2</td><td>18.2</td><td>494.02</td></tr><tr><td>3-1-1-1</td><td><strong>97.7</strong></td><td>95.4</td><td>13.6</td><td><strong>9.8</strong></td><td>783.39</td></tr><tr><td rowspan="4">Vertebral</td><td>2-2-2</td><td>91.3</td><td>79.4</td><td>7.6</td><td><strong>7.2</strong></td><td>TO</td></tr><tr><td>3-3</td><td><strong>91.8</strong></td><td><strong>80</strong></td><td><strong>7.4</strong></td><td>7.4</td><td>TO</td></tr><tr><td>1-2-3</td><td>91.1</td><td>79.4</td><td>9.4</td><td>8.2</td><td>TO</td></tr><tr><td>3-1-1-1</td><td>90.9</td><td>77.4</td><td>10.6</td><td>8.6</td><td>TO</td></tr><tr><td rowspan="4">Wine</td><td>2-2-2</td><td><strong>100</strong></td><td>93.9</td><td><strong>9</strong></td><td>8.8</td><td>3.43</td></tr><tr><td>3-3</td><td><strong>100</strong></td><td>90.5</td><td>9.6</td><td>9.6</td><td><strong>1.03</strong></td></tr><tr><td>1-2-3</td><td><strong>100</strong></td><td>93.3</td><td>14.2</td><td>11.6</td><td>5.29</td></tr><tr><td>3-1-1-1</td><td><strong>100</strong></td><td><strong>95</strong></td><td><strong>9</strong></td><td><strong>7.6</strong></td><td>9.03</td></tr><tr><td rowspan="4">Car</td><td>2-2-2</td><td>94</td><td>92.5</td><td>6.4</td><td>6.4</td><td>TO</td></tr><tr><td>3-3</td><td>93</td><td>92.3</td><td>6.4</td><td>6.4</td><td>TO</td></tr><tr><td>1-2-3</td><td>92.8</td><td>92.6</td><td><strong>6.2</strong></td><td><strong>5.8</strong></td><td>TO</td></tr><tr><td>3-1-1-1</td><td><strong>94.1</strong></td><td><strong>93.2</strong></td><td>6.8</td><td>6.8</td><td>TO</td></tr><tr><td rowspan="4">Monk2</td><td>2-2-2</td><td>87.9</td><td>79.3</td><td>12.2</td><td>12.2</td><td>TO</td></tr><tr><td>3-3</td><td>89.6</td><td>80.5</td><td><strong>10.6</strong></td><td><strong>10.6</strong></td><td>TO</td></tr><tr><td>1-2-3</td><td>90.1</td><td>74.6</td><td>13</td><td>13</td><td>TO</td></tr><tr><td>3-1-1-1</td><td><strong>92.9</strong></td><td><strong>85.8</strong></td><td>11.8</td><td>11.8</td><td><strong>695.58</strong></td></tr></tbody></table>
<figure><img src="image_320.png" /></figure>
# Constraint Programming with External Worst-Case Traversal Time Analysis

Pierre Talbot ✉️ ID

University of Luxembourg, Luxembourg
Interdisciplinary Centre for Security, Reliability and Trust (SnT), Luxembourg

Tingting Hu ☒

University of Luxembourg, Luxembourg

Nicolas Navet ✉️ ID

University of Luxembourg, Luxembourg

## Abstract

The allocation of software functions to processors under compute capacity and network links constraints is an important optimization problem in the field of embedded distributed systems. We present a hybrid approach to solve the allocation problem combining a constraint solver and a worst-case traversal time (WCTT) analysis that verifies the network timing constraints. The WCTT analysis is implemented as an industrial black-box program, which makes a tight integration with constraint solving challenging. We contribute to a new multi-objective constraint solving algorithm for integrating external under-approximating functions, such as the WCTT analysis, with constraint solving, and prove its correctness. We apply this new algorithm to the allocation problem in the context of automotive service-oriented architectures based on Ethernet networks, and provide a new dataset of realistic instances to evaluate our approach.

**2012 ACM Subject Classification** Theory of computation → Constraint and logic programming; Computer systems organization → Real-time systems; Networks → Network performance evaluation

**Keywords and phrases** Constraint programming; external function, multi-objective optimization, network analysis, worst-case traversal time analysis, abstract interpretation

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

**Supplementary Material Software:** https://github.com/ptal/automotive-network-cp/tree/cp2023, archived at swh:1:dir:e0fe246fdd3c6654293c9f5183cae7782b540d07

**Funding** Pierre Talbot: This work is supported by the Luxembourg National Research Fund (FNR) – COMOC Project, ref. C21/IS/16101289.

**Acknowledgements** We are grateful to the reviewers for their detailed comments.

## 1 Introduction

The hardware architecture of automobiles consists of dozens of interconnected electronic control units (ECUs). An important optimization problem in the field of distributed embedded system, called the *deployment problem*, is to allocate software functions to the ECUs without overloading their compute capacity and overloading the network communication links. Worst-case traversal time analysis (WCTT) is critical to ensure the communications among software functions meet hard deadlines. In most works on the deployment problem, the network considered is a controller area network (CAN), whose operating principles are relatively simple and for which an exact WCTT analysis is available [7, 39]. Therefore, the constraint model can specify both the allocation problem and the WCTT analysis. However, the newest automotive electrical-electronic (E/E) architectures rely on high-speed switched Ethernet networks. WCTT analysis already exists for Ethernet networks but is much harder to model

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

© Pierre Talbot, Tingting Hu, and Nicolas Navet;
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. 34; pp. 34:1–34:20

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

LIPICS

Leibniz International Proceedings in Informatics

Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl Publishing, Germany
as a constraint problem, especially considering the wealth of complicated quality-of-service
mechanisms available in the time-sensitive networking standards (TSN, see [22]) that are
used on top of standard Ethernet.

Constraint programming is a declarative paradigm for solving combinatorial problems. In general, the solvers are designed to work with a closed-world assumption, that is, they do not interact with external entities during solving. However, in practice, there are often parts of the model that are difficult to express as constraints, and are already programmed in another language. This is the case of the WCTT analysis. Industrial constraint solvers such as IBM ILOG CP OPTIMIZER<fnref n="1" /> and LOCAL SOLVER<fnref n="2" /> both propose *black-box expressions* to plug external functions in the model. However, these extensions are “last resort solutions” as they do not come with a semantics, and there is no guarantee on when the function is called, and how it is used within the solver. We contribute to a rigorous approach to this problem when the external function is under-approximating, i.e., it only produces valid solutions but not necessarily all.

In this work, we integrate a constraint solver for the allocation problem and a WCTT
analyser for the timing constraints. We propose a general framework, based on abstract inter-
pretation [6], for integrating external under-approximating functions (here, the WCTT ana-
lyser) in a constraint solving algorithm, and prove its correctness. The under-approximating
external function validates each solution produced by the constraint solver. Moreover, when
the external function can explain its failure, we dynamically add a new constraint to the
constraint model, approximating the reason of the failure of the external function, to improve
the quality of the subsequent solutions. In the following, we call these constraints *conflicts*<fnref n="3" />.
Because the WCTT analysis is a black-box function, deriving useful *over-approximating
conflicts* – which do not remove solutions from the problem, but might accept non-solutions –
can be difficult, or even impossible depending on the information provided by the analyser.
Our main contribution is to propose *CUSOLVE MO* a multi-objective constraint solving al-
gorithm that is *over-approximating* even if the generated conflicts are not over-approximating.
This algorithm extends the well-known multi-objective constraint programming algorithm of
Gavanelli [13] which has been frequently used in constraint optimization [19, 31, 14]. Further,
our framework can be used on top of any constraint solvers. Finally, we contribute to a new
set of benchmarks for the deployment problem and evaluate our solving algorithms on them.

2 Service Deployment Problem

For the sake of conciseness, we present the *service deployment problem* in mathematical
notation. In Appendix A, we give the constraint model in the MINIZINC constraint modelling
language [25]. The MINIZINC model is very close from the mathematical definition given
here and does not contain any particular modelling trick.

Let $\langle H, L, hc, lc \rangle$ be a weighted graph where $H$ is a set of hardware units connected by communication links $L \subseteq H \times H$. Moreover, each unit $h \in H$ has a compute capacity $hc(h)$ and each link $\ell \in L$ has a link capacity $lc(\ell)$. This graph represents a network of connected heterogeneous hardware units such as processors and switches.

<footnote n="1">¹ https://www.ibm.com/docs/en/icos/20.1.0?topic=2010-cp-optimizer-black-box-expressions</footnote>

<footnote n="2">² https://www.localsolver.com/docs/last/modelingfeatures/externalfunctions.html</footnote>

<footnote n="3">³ We avoid using the terminology of *nogood* because as we will see later, these conflicting constraints might not always preserve all solutions of the problem (over-approximating).</footnote>
Let $\langle S, Com, sc, cc \rangle$ be a weighted graph where $S$ is a set of software functions that we call services and $Com \subseteq S \times S$ is the set of communications between the services. Each service $s \in S$ consumes a certain amount of computational power $sc(s)$ and for any communication $c \in Com, cc(c)$ represents the network utilization due to this communication.

The core of the service deployment problem is to find a deployment function $d: S \to H$ allocating each service on a processor. We illustrate this problem in Figure 1 where the software graph (Figure 1b) must be deployed on the hardware graph (Figure 1a) while satisfying a number of constraints – several solutions to this particular instance are shown on Figure 4. The first constraint is on the compute capacity of each processor:

$$ \forall h \in H, \quad \sum_{s \in d^{-1}(h)} sc(s) \le hc(h) $$

It guarantees that the sum of the computational power required by all services allocated on processor $h$ does not exceed the compute capacity of $h$.

The second constraint is on the communication network:

$$ \forall \ell \in L, \quad \sum_{c \in Com} com(c, \ell) \le lc(\ell) $$

where the function $com(c, \ell)$ returns the cost on the link $\ell$ of communication $c$, and is defined by:

$$ com(c, \ell) = \begin{cases} cc(c) & \text{iff } \ell \in path(d(x), d(y)), c = (x, y) \\ 0 & \text{otherwise} \end{cases} $$

This constraint guarantees that the maximal capacity of a network link $\ell$ is never exceeded by all communications $c \in Com$ deployed on processors communicating through this link $\ell$. The function $path(h_1, h_2)$ returns a path in the hardware graph between two hardware units $h_1$ and $h_2$. In our implementation, this function represents the routing table and is given by the user as a parameter of the model. Shortest path is the standard routing strategy in automotive networks, that we use in our experiments. Follow-up work may consider the routing table as a decision variable of the model.

Finally, we note that additional constraints may be considered in similar deployment problems [15, 11]. For instance, a service might need to be allocated on a specific processor (locality constraint) or on the same processor than another service (co-location constraint). For brevity and because these constraints can be taken into account in standard ways, we choose to focus on the core problem presented above.

## 2.1 Multi-Objective Optimization

In automotive applications, we usually want to find a deployment function $d$ optimizing various objectives such as reliability [21], extensibility and cost reduction [11, 17]. We focus on two new extensibility objectives and a well-known cost reduction objective. Typically, once the services are deployed on the processors, they cannot be moved to other processors. This poses challenges when updating the system with new services. Therefore, an important goal is that the deployment function favours extensibility, that is the ability to add further services over the lifetime of the vehicle. This requirement is captured by two extensibility objectives as follows:

$$ \min \max_{h \in H} \sum_{s \in d^{-1}(h)} sc(s) $$
<figure><img src="image_323.png" /><figcaption>(a) Hardware graph with 3 ECUs and 1 switch. The values of *hc* are displayed in white and those of *lc* in blue.</figcaption></figure>

<figure><img src="image_324.png" /><figcaption>(b) Software graph with 5 services and 5 communications. The values of *sc* are displayed in red and those of *cc* in blue.</figcaption></figure>

Figure 1 An example of the software deployment problem.

which minimizes the maximum utilization rate of a processor, and

$$ \min \max_{\ell \in L} \left( \sum_{c \in Com} \frac{com(c, \ell)}{lc(\ell)} \right) $$

which minimizes the maximum utilization rate of a network link. The maximum value is considered as the corresponding hardware resource will be the bottleneck of the system. Indeed, we want to avoid a processor to be fully occupied in case a new service requires to be placed on this processor in a system update, e.g., due to the proximity of a sensor.

At the same time, we want to minimize the number of processors in order to reduce the costs:

$$ \min |d(S)| $$

Of course, this objective directly conflicts with the first one. The tradeoffs between extensibility and cost of the architecture are exposed to the system expert through a Pareto front of solutions. The final decision will involve many factors that are external to the model, such as re-use of existing processors and safety concerns.

# 3 Constraint Programming with External Function

## 3.1 Constraint Programming

A constraint satisfaction problem (CSP) is a tuple $(X, D, C)$ where $X$ is a set of variables, $D = D_1 \times \dots \times D_n$ the sets of values taken by each variable $x_i \in X$, and $C$ a set of relations over variables, called *constraints*. An assignment is a function $asn(x_i) = v_i$ from variable $x_i \in X$ to values where $v_i \in D_i$. Let $asn$ be an assignment and $c \in C$ a constraint defined on the variables $x_1, \dots, x_n$. Then the constraint $c$ is *satisfied* when $c(asn(x_1), \dots, asn(x_n))$ holds, or for short $c(asn)$. An assignment $asn$ is a *solution* when each constraint is satisfied. We write $ASN$ the set of all assignments. We call the *concrete domain* the powerset lattice $D^b = \langle P(ASN), \supseteq \rangle$ where the least element is the set of all assignments $ASN$ and the greatest element is the empty set (no solution). The set of solutions of a CSP $P = (X, D, C)$ is an element of $D^b$ computed by the following function:

$$ sol(P) := \{asn \in ASN | \forall c \in C, c(asn)\} $$
A multi-objective constraint optimization problem $M = (X, D, C, \leq)$ extends the previous definition with a partial order relation $\leq$: $ASN|_{OBJ} \times ASN|_{OBJ}$ where $ASN|_{OBJ}$ is the restriction of the assignments<fnref n="4" /> to a set of objective variables $OBJ \subseteq X$. In contrast to a single-objective optimization problem, there can be several solutions such that none is better than the others, which is why we need a partial order. Let $a, b \in ASN|_{OBJ}$. In the case of maximization over integer variables, we can define $a \leq b \Leftrightarrow \forall x \in OBJ, a(x) \leq b(x)$, that is, all the objectives of $b$ are greater or equal to the ones of $a$. We say that $b$ dominates $a$ when $a \leq b$. We write $a \succ b$ for $a \neq b \land a \geq b$ and we have $a \geq b \Leftrightarrow b \leq a$ and $a \succ b \Leftrightarrow b \prec a$. Importantly, due to the partial order, the fact that $a$ does not dominate $b$ ($a \not\geq b$) does not imply that $a$ is dominated by $b$ ($a \prec b$).

For the sake of clarity, we overload $\leq$ to work on arbitrary assignments. For any $a, b \in ASN$, we have $a \leq b \Leftrightarrow a|_{OBJ} \leq b|_{OBJ}$, and similarly for $\prec, \geq$ and $\succ$. In that case, $\leq$ is a preorder since two assignments with the same objective values might not be equal: $\leq$ lacks antisymmetry. This is not an issue since antisymmetry can be recovered by considering classes of equivalent assignments, but we do not need this construction here. The solution function for a multi-objective constraint optimization problem is then defined as:

$$sol(X, D, C, \leq) := \{a \in sol(X, D, C) | \forall b \in sol(X, D, C), b \not\leq a\}$$

Multi-objective optimization in the context of constraint programming is described in greater length in, e.g., [13, 31, 14].

In the following, we will also need to merge and generate constraints from the Pareto front. The type of a Pareto front is a set of assignments $PF := P(ASN)$. We define the operator $\sqcup: PF \times PF \rightarrow PF$ merging two Pareto fronts as $A \sqcup B := \{c \in A \cup B | \forall d \in A \cup B, d \not\leq c\}$. An equivalent definition of the solutions set is possible using $\sqcup$:

$$sol(X, D, C, \leq) := \sqcup \{\{a\} | a \in sol(X, D, C)\}$$

where $\sqcup\{s_1, \dots, s_n\} := s_1 \sqcup \dots \sqcup s_n$.

Finally, we define the function $opt: ASN \rightarrow C$ which returns a constraint ensuring that for all solutions $a \in sol(X, D, C)$, no solution in $b \in sol(X, D, C \land opt(a))$ is dominated by $a$, i.e. $a \not\geq b$. This function is defined by:

$$opt(a) := \bigvee_{x \in OBJ} x < a(x)$$

It generates a constraint requiring at least one of the objective variables to be strictly better than the one obtained in $a$. This approach to multi-objective optimization was pioneered by Gavanelli [13].

## 3.2 Abstract Constraint Programming

In general, the set $sol(P)$ might not be efficiently computed on the concrete domain. In constraint reasoning by abstract interpretation [10, 27, 33, 35], they design an abstract solving function $sol_o^*(P)$ which over-approximates the solution set, i.e., $sol_o^*(P) \supset sol(P)$. Dually, we can also design an under-approximating solving function such that $sol_u^*(P) \subseteq sol(P)$. Over-approximation contains all solutions but might contains non-solution assignments as well, while under-approximation only contains solutions but not necessarily all solutions. For

<footnote n="4">⁴ Formally, $ASN|_{OBJ} := \{asn|_{OBJ} | asn \in ASN\}$ where $asn|_{OBJ}(x) = asn(x)$ for all $x \in OBJ$.</footnote>
instance, discrete constraint programming solvers are both under- and over-approximating [10, 35], and continuous constraint programming solvers are over-approximating [27]. Incomplete discrete solvers, such as those based on local search, can be viewed as under-approximating solving functions.

## 3.3 Abstract Constraint Model

We can also use the abstraction framework at the level of the constraint model. In industry, some elements of a constraint model might already be available and tested, and it is usually not practical to spend time redeveloping those parts as a constraint problem. Sometimes, the problem is just too difficult to be expressed as a constraint model in a reasonable amount of time; this is the case of the WCTT analysis for instance. In these cases, the problem *P* is never explicitly written as a constraint model. Instead, we can rely on an over-approximating model *O* of *P*, such that *sol*(*O*) ⊇ *sol*(*P*). We often have an idea of some constraints that must be satisfied in any solution of the model but we do not necessarily know them all. This model *O* can be solved by an over-approximating function *sol*<sub>*o*</sub><sup>‡</sup>(*O*) ⊇ *sol*(*O*) – although *O* simplifies *P*, it might still not be efficiently computable. If *O* is unsatisfiable (*sol*<sub>*o*</sub><sup>‡</sup>(*O*) = {}), then the problem *P* is unsatisfiable as well since only *sol*(*P*) = {} satisfies *sol*<sub>*o*</sub><sup>‡</sup>(*O*) ⊇ *sol*(*P*). In the following, we denote *OSOLVE*(*O*) ∈ *sol*<sub>*o*</sub><sup>‡</sup>(*O*), the solving algorithm computing a single solution of *O* and returning {} if *O* is unsatisfiable. Its definition in terms of abstract interpretation can be found in [1, 10, 27]. Dually, we can also propose an under-approximating model *U* of *P* and its solving function *sol*<sub>*u*</sub><sup>‡</sup>(*U*) such that *sol*<sub>*u*</sub><sup>‡</sup>(*U*) ⊆ *sol*(*U*) ⊆ *sol*(*P*). If *U* is satisfiable, then the problem *P* is satisfiable as well. An under-approximation makes additional assumptions about the reality, and therefore might discard solutions of *P*. Therefore, the real problem *P* is framed between an over-approximating model *O* and an under-approximating model *U*, which is summarized by *sol*<sub>*o*</sub><sup>‡</sup>(*O*) ⊇ *sol*(*O*) ⊇ *sol*(*P*) ⊇ *sol*(*U*) ⊇ *sol*<sub>*u*</sub><sup>‡</sup>(*U*).

## 3.4 Under-Approximating External Function

We must go one step further for our abstract framework to be useful in practice. If we cannot explicitly list the constraints of the concrete problem *P*, it seems unlikely that we could list *more constraints* in an under-approximating model *U*. Nevertheless, when given a solution to *O*, it can often be validated by existing code developed by domain experts. For instance, worst-case analysis such as WCTT and feasibility tests fall in this category. They conservatively analyse the network architecture, and discard some solutions that would be valid but could not be proven valid by the analysis. In practice, worst-case analysis are not directly working with constraints and domains, and thus do not explicitly define an under-approximating constraint model *U*, but they work on assignments.

We formally define the analysis as an under-approximating function *uf*: *ASN* → *C*. The function *uf* returns a conflict constraint when the assignment generated by *OSOLVE* is not in *sol*(*U*), or *true* if it is in *sol*(*U*). Actually, we define the solutions of the under-approximating model *U* as the set of all solutions accepted by *uf*:

$$ sol(U) := uf^{-1}(true) = \{asn \in ASN | uf(asn) = true\} $$

A general conflict automatically available to all functions *uf*, is the logical negation of the assignment (**NA**): $\neg asn := \neg(x_1 = asn(x_1) \land \dots \land x_n = asn(x_n)) \Leftrightarrow x_1 \neq asn(x_1) \lor \dots \lor x_n \neq asn(x_n)$. However, it is a weak conflict since it only prevents *OSOLVE* from returning to this assignment, without providing additional pruning. A conflict is *over-approximating* if it does not remove valid solutions: for all assignments *asn*, we have *sol*(*U*) ⊆ *sol*(*O* ∧ *uf*(*asn*)).
**function** `USOLVE(O, ufo)`
  $S \leftarrow \{\}$
  `asn` $\leftarrow$ `OSOLVE(O)`
**while** `asn` $\neq \{\}$ **do**
  **if** `ufo(asn)` = `true` **then**
    $S \leftarrow S \cup \{asn\}$
    $O \leftarrow O \land \neg asn$
  **else**
    $O \leftarrow O \land ufo(asn)$
**end if**
`asn` $\leftarrow$ `OSOLVE(O)`
**end while**
**return** $S$
**end function**

(a) Find all satisfiable solutions.

**function** `USOLVE_MO(O, ufo, ∪, opt)`
  $F \leftarrow \{\}$
  `asn` $\leftarrow$ `OSOLVE(O)`
**while** `asn` $\neq \{\}$ **do**
  **if** `ufo(asn)` = `true` **then**
    $F \leftarrow F \cup \{asn\}$
    $O \leftarrow O \land opt(asn)$
  **else**
    $O \leftarrow O \land ufo(asn)$
**end if**
`asn` $\leftarrow$ `OSOLVE(O)`
**end while**
**return** $F$
**end function**

(b) Multi-objective version of `USOLVE`.

■ **Figure 2** Constraint solving with external under-approximating function producing over-approximating conflicts.

We say that a conflict *c* is *sound* if it implies **NA**; in other terms, it excludes the current assignment: `asn` $\notin$ `sol(c)`. **NA** is a sound over-approximating conflict. We denote by `ufo` an under-approximating external function returning sound over-approximating conflicts.

Let *O* be an over-approximating model and *ufo* a sound under-approximating external function. The function *USOLVE* presented in Algorithm 2a constructs the solutions set *S* of the under-approximating model *U*. Constraint programming helps us navigating in the under-approximated solution space of the external function more efficiently. Without it, we would need to call *ufo* on many more unsatisfiable assignments, since those would not be removed by a constraint solver. The next proposition shows that *USOLVE* computes an under-approximation of *P*.

**Proposition 1.** *USOLVE is a sound under-approximating function, that is, $USOLVE(O, ufo) = ufo^{-1}(true) \subseteq sol(P)$.*

*Proof.* Let $O_i$ and $S_i$ be the variables *O* and *S* at the *i*th iteration of the loop where $O_0 = O$ and $S_0 = \{\}$. We must show that at the final iteration *n*, we have $ufo^{-1}(true) = S_n$. We proceed inductively by defining $O_{i+1}$ and $S_{i+1}$ as follows:

1. If `asn` is a solution to `ufo`: $O_{i+1} = O_i \land \neg asn$ and $S_{i+1} = S \cup \{asn\}$. In that case we have $sol(O_i) \cup S_i = sol(O_{i+1}) \cup S_{i+1}$.

2. If `asn` is not a solution to `ufo`: $O_{i+1} = O_i \land ufo(asn)$ and $S_{i+1} = S$. In that case we have $sol(O_i) \cup S_i \supseteq sol(O_{i+1}) \cup S_{i+1}$. Since the conflict must be over-approximating, no assignment removed from $O_i$ is in $ufo^{-1}(true)$ and therefore $ufo^{-1}(true) \subseteq sol(O_{i+1}) \cup S_{i+1}$.

The final iteration is necessarily with $O_n = \{\}$, hence we must have $ufo^{-1}(true) \subseteq S_n$. Since we only add in $S_n$ the assignments `asn` such that `ufo(asn) = true`, we also have $ufo^{-1}(true) \supseteq S_n$. ◀

The extension to multi-objective optimization is a small modification of `USOLVE`. The set *F* represents the Pareto front of the problem. To compute the Pareto front, we introduce the algorithm `USOLVE_MO` in Figure 2b, and we highlight the differences with `USOLVE` in green.
**► Proposition 2.** `USOLVE MO` is a sound under-approximating function. Moreover, we have:

1. `USOLVE MO(O, ufo, ⋃, opt)` = `sol(U, ⩽)`,

2. `USOLVE MO(O, ufo, ⋃, opt)` ⊆ `USOLVE(O, ufo)`, and

3. `USOLVE MO(O, ufo, ⋃, opt)` = `{}` ⇔ `USOLVE(O, ufo)` = `{}`.

*Proof.* We prove each statement in turn:

1. Each time we reach a solution `asn`, we add the constraint `opt(asn)` to `O`. By definition, all solutions removed by this constraint are dominated by `asn`, and therefore it cannot remove solutions from `sol(U, ⩽)`.

2. Notice that `opt(asn) ⇒ ¬asn`, and therefore it can only prune more solutions in comparison to `USOLVE`.

3. Before reaching the first solution, the algorithm behaves in the same way than `USOLVE`. ◀

The risk when navigating in an under-approximating solution space is to find no solution at all. Therefore, it is useful to notice that the multi-objective algorithm returns an empty set of solution only if `USOLVE` does as well (by Proposition 2(3)).

## 4 Worst-Case Traversal Time Analysis

We focus on the worst-case traversal time (WCTT) analysis, which verifies if the network communications among the deployed services meet timing constraints, i.e., deadline constraints in this work. In an automotive network, we must ensure the deadlines of network packets are met, which is crucial for safety reasons (e.g., a message sent to an airbag arrives on time) and other non-functional requirements (e.g., the speakers must be synchronized when playing music). WCTT analysis is a formal method which provides upper bounds on the worst-case delay of every packet sent in the network. It is therefore an under-approximating external function because all deployments passing this analysis will also fulfill the real-time constraints in reality. But some deployments, that in fact meets all timing constraints, will not pass the WCTT analysis because it only gives an upper-bound on the delay: it is a sufficient but not necessary condition.

The WCTT analysis for Ethernet networks, based on *network calculus* [18], is mathematically complicated (see for instance [28]). We think it would take tremendous efforts to model WCTT analysis as a constraint problem if the goal is to develop an implementation that is sufficiently accurate to be used on real-world problems, given the complexity of the WCTT analysis after 30 years of research. As an illustration, the network calculus engine from the company RTAW we use in the paper has been developed for 15 years and implements state-of-the-art techniques such as [34, 3]. We take a more pragmatic approach where we reuse an existing WCTT analyser, and integrate it in our framework as an under-approximating external function.

Let us first describe the output of a WCTT analysis. For each communication $(x, y) \in Com$, the WCTT analysis outputs the worst-case end-to-end delay of a packet traversing the network from $d(x)$ to $d(y)$. A negative delay means the deadline for that communication cannot be met and thus the deployment $d$ is unsatisfiable from the point of view of the analysis. From an unsatisfiable assignment, we can think of various conflicts such as forcing the network load of the problematic link to be smaller, or forbidding to allocate the services $x$ or $y$ on their current processors. Unfortunately, conflicts that are intuitive are often not over-approximating. To complicate the finding of over-approximating conflicts, the WCTT analysis is non-monotonic w.r.t. the network load. Indeed, it can happen that increasing the load of an unsatisfiable network turns it into a satisfiable network according to the WCTT analysis. This is a well-known phenomenon referred to as *timing anomaly* (see, e.g., [23]), that may occur with non-preemptive scheduling as in communication networks.
function CUSOLVE MO(F, O, C, uf, , opt)

$asn \leftarrow OSOLVE(O \land C)$

if asn ≠ {} then

$co \leftarrow uf(asn)$

if co = true then

$F \leftarrow F \sqcup \{asn\}$

$O \leftarrow O \land opt(asn)$

$CUSOLVE\_MO(F, O, C, uf, \sqcup, opt)$

else

$CUSOLVE\_MO(F, O, C \land co, uf, \sqcup, opt)$

$CUSOLVE\_MO(F, O, C \land \neg co \land \neg asn, uf, \sqcup, opt)$

end if

end if

return F

end function

Figure 3 Multi-objective constraint solving with under-approximating external function returning conflicts that are not over-approximating.

Although the conflicts mentioned above are not over-approximating, they can nevertheless be useful as heuristics to find a solution faster. In Figure 3, we extend `USOLVE MO` in the case where `uf` is returning sound conflicts that are not over-approximating. The intuition is that the conflict is viewed as a branching decision, and thus backtracked when the sub-problem has been fully explored. Doing so, we do not lose the completeness of our solving algorithm. We provide an example unrolling this algorithm in Section 4.2.

We note that if the conflict `co` is over-approximating, then `O ∧ ¬co` is necessarily unsatisfiable. In that case, `CUSOLVE MO` remains correct but the second recursive call `CUSOLVE MO(F, O, C ∧ ¬co, uf, , opt)` is unnecessary.

► **Proposition 3.** CUSOLVE MO(Ø, O, Ø, uf, , opt) = sol(U, ≤)

*Proof.* The difference with Proposition 2, is that *uf* does not necessarily produce over-approximating conflicts. However, given any constraint problem *O* and any constraint *co*, we always have *sol*(*O* ∧ *co*) ∪ *sol*(*O* ∧ ¬*co*) = *sol*(*O*). Therefore, the same solution space is eventually explored, and the same Pareto front is found. However, the conflicting assignment *asn* might be reexplored in the right branch. Indeed, *sol*(*¬co*) = *ASN* \ *sol*(*co*) and therefore *asn* ∈ *sol*(*¬co*). By forbiding revisiting this assignment in both the left and right branches, we guarantee progress and thus termination of the algorithm.

◀

## 4.1 Conflicts for WCTT

Let asn be the current assignment for which WCTT detected that the communication (x, y) ∈ Com does not meet its deadline. Among the possible conflicts, we experiment with the following ones in the next section:

- Forbid the source service *x* to be allocated on its current hardware unit or the one of *y*:

$$d(x) \notin \{asn(x), asn(y)\} \qquad (\text{FS})$$

- Forbid the target service *y* to be allocated on its current hardware unit or the one of *x*:

$$d(y) \notin \{asn(x), asn(y)\} \qquad (\text{FT})$$
34:10 Constraint Programming with External Worst-Case Traversal Time Analysis

Decrease the number of hops in the network between x and y:

$$ |path(d(x), d(y))| < |path(asn(x), asn(y))| \qquad (\mathbf{DH}) $$

We also test two conflicts that are global to the network, thus do not consider a specific communication. It is based on the observation that a communication can fail to meet its deadline because of other communications. We let the variable load<sub>l</sub> = &Sigma;<sub>c&isin;Com</sub> com(c, l) to be the network load of the link l in the current assignment asn.

Decrease the load of at least one network link:

$$ \bigvee_{l \in L} \sum_{c \in Com} com(c, l) \leq load_{l} \qquad (\mathbf{D1L}) $$

Decrease the load of the most occupied network link, with m = maxarg<sub>l&isin;L</sub> (load<sub>l</sub>/lc(l)):

$$ \sum_{c \in Com} com(c, m) \leq load_{m} \qquad (\mathbf{DML}) $$

Taking the conjunction of any two conflicts will further prune the search tree, while taking their disjunction will create a weaker conflict. In the experiments, we test the conjunction of FS and FT that we name FST. In our investigations, we found that, in general, the disjunction of conflicts did not help to reach a (better) solution faster.

## 4.2 An Example of the Algorithm cusolve\_mo

We unroll the algorithm CUSOLVE\_MO with the DH conflict strategy on the software deployment problem given in Figure 1. Initially, the Pareto front F and the set of conflicts C are empty. The first call to OSOLVE returns a solution to the problem as depicted in the orange box labelled A. Because one of the network link is used at 100% capacity, we suppose the communication between S<sub>1</sub> and S<sub>4</sub> fails to meet its deadline, hence the WCTT analysis fails on this solution. The conflict |path(d(S<sub>1</sub>), d(S<sub>4</sub>))| < 1, reducing the number of hops between S<sub>1</sub> and S<sub>4</sub>, is added to the conflicts set C. In this case, this conflict forces both services to be allocated on the same processor.

The model is solved again with this new conflict and OSOLVE returns a solution as depicted in the box B. This time the WCTT analysis succeeds, and the solution is added to the Pareto front F. The objectives are (100, 80, 2) where 100 is the maximum utilization rate among all processors, 80 is the maximum utilization among all network links and 2 is the number of cores used.

As long as the WCTT analysis succeeds, the OSOLVE procedure is iteratively called with the updated Pareto front. In the box C, we have a solution (90, 90, 3) incomparable to (100, 80, 2), hence the Pareto front now contains both. In the box D, we find the solution (90, 80, 3) which dominates the previous one. Afterwards, the OSOLVE procedure finds the problem unsatisfiable, which means there is no solution better than the ones found previously. However, we have previously added a conflict which was not over-approximating, and therefore we might have missed solutions of the problem. Therefore, we need to backtrack and explore the problem with the negation of the conflict. In box E, the model is solved again with the latest Pareto front, and a new non-dominated solution (80, 90, 3) is found. This is repeated and a better solution (70, 90, 3) is found in box F. As long as the WCTT analysis succeeds, the solving procedure continues, and when it fails we branch as we did in the first node.

This example demonstrates that conflicts are heuristics which are used in a way that do not prevent to find all solutions. This is also why the non-monotonicity of the WCTT analysis is not an issue: the entire solution space is eventually explored.
<figure><img src="image_325.png" /><figcaption>Figure 4 An example unrolling CUSOLVE_MO on a small network. The orange boxes represent the solutions found by OSOLVE. F is the current Pareto front and is automatically added to the model before solving; note that the Pareto front is preserved on backtracking.</figcaption></figure>
# 5 Implementation and Experiments

The code of the MINIZINC model, the data of the instances and the implementation of the algorithms used can be found online at https://github.com/ptal/automotive-network-cp/tree/cp2023.

## 5.1 Experimental Setting

We run all the experiments on an AMD Epyc ROME 7H12 processor (64 cores, 280W). The constraint programming solver implementing the OSOLVE solving function is GECODE 6.3.0 [32] in parallel mode with 8 cores and 16 threads. We use a first-fail variable selection strategy (the variable with the smallest domain is chosen first) and a random value selection, as it shows better performance than the free search strategy of GECODE. We also tried CHUFFED 0.10.4 [26], a hybrid solver between SAT and constraint programming, but it did not outperform Gecode on our problem. Due to the lack of open-source alternative, the WCTT analysis is performed by the proprietary software RTAW-PEGASE-4.3.7 [29], which implements state-of-the-art network calculus algorithms [18, 2, 4]. The WCTT analysis takes on average 1.5 seconds to run, and this time remains stable across instances. The algorithms presented in this paper – USOLVE, USOLVE MO and CUSOLVE MO – are implemented in Python using MINIZINC PYTHON 0.9.0 [9]. In addition we provide OSOLVE MO which implements the multi-objective optimization solving procedure of [13] – it is the same than USOLVE MO but without the external function filtering. Although multi-objective optimization is very important in practice, it is not natively available in every constraint programming solver (for instance in GECODE, CHUFFED or ORTOOLS). Similarly to [14], our approach does not require to modify the constraint solver, but the solver state is lost between two calls to OSOLVE which might be less efficient – although the issue is mitigated since we use a random search strategy. It can be seen as a restart strategy triggered on every solution. It is also similar to what is done in MINISEARCH [30] to design search strategies generically across solvers.

## 5.2 Dataset Description

The instances are derived from a realistic automotive Ethernet network, shown in Figure 5, consisting of 19 network devices (14 ECUs and 5 switches) provided by the company RealTime-at-Work<fnref n="5" />. The experiments consider 5 problem instances of 50 services, 5 instances of 75 services and 8 instances of 100 services. The numbers of communications vary among the instances, but are between 125% and 135% of the number of services. An information missing in the network description is the CPU usage for each service. For each of the 18 instances, we generated 10 versions where the sum of all computational requirements is 20%, 40%, 60%, 80% and 90% of the total computational capacity of all ECUs with a uniform distribution among services. To summarize, an instance named I5_75-14-u60 has 75 services allocated on 14 processors and using 60% of the total computational power, and I5 denotes the fifth instance with 75 services. In total, we have a new dataset of 90 MiniZinc instances for the deployment problem. In the following, we present experimental results for a subset of these instances (I{1,2,5}_{50,75,100}-14-u{20,40,60,80,90}), totalizing 45 instances. We set a timeout on the constraint solver of 30 minutes for each instance, and unlimited time for the WCTT analysis.

<footnote n="5">⁵ https://www.realtimeatwork.com/</footnote>
<figure><img src="image_326.png" /><figcaption>Figure 5 Realistic automotive Ethernet network used in the experiments.</figcaption></figure>

<figure><img src="image_327.png" /><figcaption>Figure 6 Cumulated hypervolume score for each experiment over all instances.</figcaption></figure>

## 5.3 Evaluation of cusolve\_mo

We evaluate the algorithm CUSOLVE\_MO on NA and the seven strategies presented in Section 4.1. For a single assignment, it is possible that several communications cannot meet their deadlines, and thus several conflicts are generated. We write $FS_\lor$ when these conflicts are combined disjunctively and $FS_\land$ when they are combined conjunctively. It only impacts the conflicts that are local to a communication ($FS$, $FT$, $DH$ and $FST$), thus we have 11 conflicts in total.

Our main comparison metrics is the hypervolume of the Pareto front which is standard in multi-objective optimization. For all 45 instances, none of the algorithms tested could find the optimum within the time limit. We give a general picture of the situation in Figures 6 and 7. Overall, the decreasing hops strategy is the best, and finds the best hypervolumes on 19 of the 45 instances. We also witness a smaller number of conflicts, which means that $DH$ is effective to search the state-space of $sol(U)$. When considering the score, forbidding the
<figure><img src="image_328.png" /><figcaption>Figure 7 Number of times each experiment computed the best hypervolume.</figcaption></figure>

source and target services on a particular ECU are usually not better than simply using the **NA** conflict. These strategies are still superior regarding the number of times they find the best hypervolume.

In addition, we evaluate `CUSOLVE_MO` against a more straightforward *two-steps* algorithm `OSOLVE_MO THEN UF` (denoted by `MO UF`) where the Pareto front is first fully generated, and then filtered by the `uf` function. To improve this method, we keep all intermediate solutions when building the Pareto front in a set `S`. During the filtering step, if a solution `a` is discarded by the function `uf`, we remove `a` from the Pareto front and reconstruct it with $\bigcup\{\{b\} \mid b \in S \setminus \{a\}\}$. It does not make this algorithm over-approximating as it can still discard assignments accepted by `uf`, but it improves the filtered Pareto front.

As shown in Figure 6, `MO UF` ranks fourth, and therefore is a good approach to solve the deployment problem when we do not seek (or cannot find) the true optimal solution. Interestingly, for 18 instances over the 45, the hypervolume before and after filtering is the same, which means that all solutions of the Pareto front were valid w.r.t. `uf`. This is particularly true with 50 services where 14/15 instances have the same hypervolume before and after filtering. This result is explained by noticing that adding more services has a higher impact on the network load, and thus the WCTT analysis fails more often. Over the 30 instances with 75 and 100 services, there are 24 instances that have a filtered hypervolume within 3% of the unfiltered hypervolume. It is not always the case as for the instances `I1_100-14-u60` and `I1_100-14-u80`, the filtered hypervolume is respectively 57% and 73% of the unfiltered hypervolume. An advantage of the offline filtering proposed by `MO UF` is to call `uf` an order of magnitude less than with `DH`. Over all instances, `MO UF` calls `uf` 1180 times while `DH` calls `uf` 12245 times. Therefore, depending on the time taken by the external function and the number of conflicts, `MO UF` can be better – especially for low number of conflicts and long evaluation time.

From an implementation perspective, using `MINIZINC PYTHON` allowed us to implement an algorithm generic across solvers, but it incurs a cost. Besides losing the state of the solver between calls, we must call the `MINIZINC` (source of the model) to `FLATZINC` (simpler format supported by solvers) translator, and it takes on average around 40% of the total solving time. An improvement to `MINIZINC PYTHON` would be to directly add `FLATZINC` constraints to avoid recompiling the `MINIZINC` model each time.
# 6 Related Work

## 6.1 CAN Networks

The deployment problem has been extensively studied due to its importance in distributed real-time embedded systems. It was pioneered in [38] for tasks allocation on *controller area network* (CAN). In CAN network, the hardware units are all connected on a broadcast bus, and therefore the hardware network is fully connected. The main difference with our work is that we consider a more general switch-based network. In the context of real-time and critical systems, such as those found in the automotive industry, it is crucial to ensure the network will not be overloaded by communication, and when required, that the network packet deadlines are met. The WCTT analysis on switch-based networks is more complicated and under-approximating (it does not give an exact upper bound), while it is an exact analysis for CAN network [37, 7, 39].

Due to its simpler nature, schedulability analysis, such as WCTT, over CAN networks has been directly incorporated in the constraint model before. The work of Hladik et al. [15] is the first to model the deployment problem over CAN network using constraint programming. They model the schedulability analysis as a global constraint. Alternatively, they also use a method inspired by logic-based Benders decomposition (LBBD) [16] to separate the allocation problem solved using constraint programming and the schedulability analysis solved by an ad-hoc algorithm. It differs from our approach mainly because there is no notion of approximation, and the conjunction of both parts models the problem exactly. Moreover, they consider only the satisfiability of the problem, and they do not seek to optimize one or more objectives.

Other techniques were proposed to solve the deployment over CAN networks with multi-objective optimization, for instance, evolutionary optimization [21], ant colony system with constraint propagation and without searching [36], mixed integer linear programming (MIP) [24] and satisfiability modulo theories (SMT) [11]. These methods are either incomplete (no proof of optimality or unsatisfiability) and thus under-approximating, or they are complete (MIP and SMT) which is only possible because they model a simpler problem (CAN network).

## 6.2 Switched Networks

To the best of our knowledge, Kugele et al. [17] are the first to configure and analyse an application distributed over a switched network using a SMT solver. Similarly to `OSOLVE MO`, they generate a Pareto front and then verify the produced solutions. However, the verification is performed using simulation, which does not give a formal worst-case guarantee. Therefore, their solving method is over-approximating and the obtained solutions are not guaranteed to be valid. Besides, they do not provide the constraint model and only tested their algorithm on a small network of 3 ECUs and 25 services.

## 6.3 Other Applications

Campeanu et al. [5] study the deployment problem in heterogeneous architectures (CPU, GPU and FGPA), but with communication still happening over a CAN network. *Satisfiability Modulo Discrete Event Simulation* [20] combines a SAT solver with discrete event simulation (DES) for a railway construction planning problem. The combination of both techniques share similarities with `CUSOLVE MO` since the DES simulator is encapsulated as a theory and provide conflict to the SAT solver – but, like us, only on full assignments. However, simulation is an over-approximating technique and therefore the global method remains over-approximating. Moreover, the algorithms are specialized to the railway construction problem and no general algorithm or correctness proof is given.
## 6.4 Online and Dynamic Constraint Programming

Our work is related to online constraint programming [12] – also called dynamic constraint programming [8] – as in both approaches the model is incrementally refined. A difference is that online constraint programming is primarily designed when the solutions generated are used in real-time, and variables impacting the past decisions cannot be modified in the subsequent solving steps. We do not have such real-time requirements since the new data are obtained from an offline analysis. Moreover, in [12], they propose to internalize the dynamic part of the problem inside the model. Here, we purposely delegated a part of the model to an external function, which would have been prohibitively complicated to model otherwise.

# 7 Conclusion

We study the deployment problem, an important problem in the field of distributed real-time embedded systems, and more specifically in the automotive industry. This problem has a task allocation part, which is efficiently solved by constraint programming, and a scheduling network analysis part, which is efficiently solved by a WCTT analysis. The integration of both techniques is difficult since both parts are black-box functions. We propose the algorithm CUSOLVE\_MO which combines both parts in a loosely coupled manner, thus making our framework reusable on other similar problems. Our approach is based on abstract interpretation, a formal method allowing us to prove properties of our algorithms. Finally, we evaluated our approach on a new dataset for the deployment problem – since none existed before – and conclude that the cooperation scheme proposed by CUSOLVE\_MO works better than solving each part in sequence.

# References

1 Krzysztof R. Apt. The essence of constraint propagation. *Theoretical computer science*, 221(1-2):179-210, 1999. doi:10.1016/S0304-3975(99)00032-8.

2 Anne Bouillard and Éric Thierry. An algorithmic toolbox for network calculus. *Discrete Event Dynamic Systems*, 18(1):3-49, 2008. doi:10.1007/s10626-007-0028-x.

3 Marc Boyer and Hugo Daigmorte. Improved service curve for element with known transmission rate. *IEEE Networking Letters*, 5(1):46-49, 2023. doi:10.1109/LNET.2022.3150649.

4 Marc Boyer, Jörn Migge, and Nicolas Navet. An efficient and simple class of functions to model arrival curve of packetised flows. In *Proceedings of the 1st International Workshop on Worst-Case Traversal Time*, WCTT '11, pages 43-50, New York, NY, USA, 2011. Association for Computing Machinery. doi:10.1145/2071589.2071595.

5 Gabriel Campeanu, Jan Carlson, and Severine Sentilles. Component Allocation Optimization for Heterogeneous CPU-GPU Embedded Systems. In *2014 40th EUROMICRO Conference on Software Engineering and Advanced Applications*, pages 229-236, Verona, Italy, 2014. IEEE. doi:10.1109/SEAA.2014.29.

6 Patrick Cousot and Radhia Cousot. Abstract Interpretation: A Unified Lattice Model for Static Analysis of Programs by Construction or Approximation of Fixpoints. In *Proceedings of the 4th ACM SIGACT-SIGPLAN Symposium on Principles of Programming Languages*, POPL '77, pages 238-252, New York, NY, USA, 1977. Association for Computing Machinery. doi:10.1145/512950.512973.

7 Robert I. Davis, Alan Burns, Reinder J. Bril, and Johan J. Lukkien. Controller Area Network (CAN) Schedulability Analysis: Refuted, Revisited and Revised. *Real-Time Systems*, 35(3):239-272, 2007. doi:10.1007/s11241-007-9012-7.
8 Rina Dechter and Avi Dechter. Belief Maintenance in Dynamic Constraint Networks. In *Proceedings of the Seventh AAAI National Conference on Artificial Intelligence*, AAAI'88, pages 37-42. AAAI Press, 1988.

9 Jip J. Dekker. MiniZinc Python, 2023. URL: https://github.com/MiniZinc/minizinc-python.

10 Vijay D'Silva, Leopold Haller, and Daniel Kroening. Abstract satisfaction. In *Proceedings of the 41st ACM SIGPLAN-SIGACT Symposium on Principles of Programming Languages - POPL '14*, pages 139-150, San Diego, California, USA, 2014. ACM Press. doi:10.1145/2535838.2535868.

11 Johannes Eder, Sebastian Voss, Andreas Bayha, Alexandru Ipatiov, and Maged Khalil. Hardware architecture exploration: automatic exploration of distributed automotive hardware architectures. *Software and Systems Modeling*, 19(4):911-934, 2020. doi:10.1007/s10270-020-00786-6.

12 Alexander Ek, Maria Garcia de la Banda, Andreas Schutt, Peter J. Stuckey, and Guido Tack. Modelling and Solving Online Optimisation Problems. *Proceedings of the AAAI Conference on Artificial Intelligence*, 34(02):1477-1485, 2020. doi:10.1609/aaai.v34i02.5506.

13 Marco Gavanelli. An algorithm for multi-criteria optimization in CSPs. In *ECAI 2002: 15th European Conference on Artificial Intelligence, July 21-26, 2002, Lyon France: Including Prestigious Applications of Intelligent Systems (PAIS 2002): Proceedings*, volume 77, page 136. IOS Press, 2002.

14 Tias Guns, Peter J. Stuckey, and Guido Tack. Solution Dominance over Constraint Satisfaction Problems, 2018. arXiv:1812.09207 [cs]. URL: http://arxiv.org/abs/1812.09207.

15 Pierre-Emmanuel Hladik, Hadrien Cambazard, Anne-Marie Déplanche, and Narendra Jussien. Solving a real-time allocation problem with constraint programming. *Journal of Systems and Software*, 81(1):132-149, 2008. doi:10.1016/j.jss.2007.02.032.

16 J.N. Hooker and G. Ottosson. Logic-based Benders decomposition. *Mathematical Programming*, 96(1):33-60, 2003. doi:10.1007/s10107-003-0375-9.

17 Stefan Kugele, Philipp Obergfell, and Eric Sax. Model-based resource analysis and synthesis of service-oriented automotive software architectures. *Software and Systems Modeling*, 20(6):1945-1975, December 2021. doi:10.1007/s10270-021-00896-9.

18 Jean-Yves Le Boudec and Patrick Thiran. Network Calculus. In *Network Calculus: A Theory of Deterministic Queuing Systems for the Internet*, pages 3-81. Springer Berlin Heidelberg, Berlin, Heidelberg, 2001. doi:10.1007/3-540-45318-0_1.

19 Martin Lukasiewycz, Michael Glaß, Christian Haubelt, and Jürgen Teich. Solving Multi-objective Pseudo-Boolean Problems. In *Theory and Applications of Satisfiability Testing - SAT 2007*, pages 56-69. Springer Berlin Heidelberg, Berlin, Heidelberg, 2007. doi:10.1007/978-3-540-72788-0_9.

20 Bjørnar Luteberget, Koen Claessen, Christian Johansen, and Martin Steffen. SAT modulo discrete event simulation applied to railway design capacity analysis. *Formal Methods in System Design*, 57(2):211-245, August 2021. doi:10.1007/s10703-021-00368-2.

21 Irene Moser and Sanaz Mostaghim. The automotive deployment problem: A practical application for constrained multiobjective evolutionary optimisation. In *IEEE Congress on Evolutionary Computation*, pages 1-8, Barcelona, Spain, July 2010. IEEE. doi:10.1109/CEC.2010.5585991.

22 Ahmed Nasrallah, Akhilesh S. Thyagaturu, Ziyad Alharbi, Cuixiang Wang, Xing Shao, Martin Reisslein, and Hesham ElBakoury. Ultra-Low Latency (ULL) Networks: The IEEE TSN and IETF DetNet Standards and Related 5G ULL Research. *IEEE Communications Surveys & Tutorials*, 21(1):88-145, 2019. doi:10.1109/COMST.2018.2869350.

23 Mitra Nasri, Sanjoy Baruah, Gerhard Fohler, and Mehdi Kargahi. On the Optimality of RM and EDF for Non-Preemptive Real-Time Harmonic Tasks. In *Proceedings of the 22nd International Conference on Real-Time Networks and Systems - RTNS '14*, pages 331-340, Versaille, France, 2014. ACM Press. doi:10.1145/2659787.2659806.
34:18 Constraint Programming with External Worst-Case Traversal Time Analysis

24 Asef Nazari, Dhananjay Thiruvady, Aldeida Aleti, and Irene Moser. A mixed integer linear programming model for reliability optimisation in the component deployment problem. *Journal of the Operational Research Society*, 67(8):1050–1060, August 2016. doi:10.1057/jors.2015.119.

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

26 Olga Ohrimenko, Peter J. Stuckey, and Michael Codish. Propagation via Lazy Clause Generation. *Constraints*, 14(3):357–391, September 2009. doi:10.1007/s10601-008-9064-x.

27 Marie Pelleau, Antoine Miné, Charlotte Truchet, and Frédéric Benhamou. A Constraint Solver Based on Abstract Domains. In *Verification, Model Checking, and Abstract Interpretation*, pages 434–454. Springer, 2013. doi:10.1007/978-3-642-35873-9_26.

28 R. Queck. Analysis of Ethernet AVB for automotive networks using Network Calculus. In *2012 IEEE International Conference on Vehicular Electronics and Safety (ICVES 2012)*, pages 61–67, July 2012. doi:10.1109/ICVES.2012.6294261.

29 RealTime-at-Work. Rtaw-pegase, 2022. URL: https://www.realtimeatwork.com/rtaw-pegase/.

30 Andrea Rendl, Tias Guns, Peter J. Stuckey, and Guido Tack. MiniSearch: a solver-independent meta-search language for MiniZinc. In *Principles and Practice of Constraint Programming*, pages 376–392. Springer, 2015. doi:10.1007/978-3-319-23219-5_27.

31 Pierre Schaus and Renaud Hartert. Multi-Objective Large Neighborhood Search. In *Principles and Practice of Constraint Programming*, volume 8124, pages 611–627. Springer Berlin Heidelberg, Berlin, Heidelberg, 2013. doi:10.1007/978-3-642-40627-0_46.

32 Christian Schulte, Guido Tack, and Mikael Lagerkvist. *Modeling and Programming with Gecode*, 2020.

33 Joseph Scott. *Other Things Besides Number: Abstraction, Constraint Propagation, and String Variable Types*. PhD thesis, Acta Universitatis Upsaliensis, Uppsala, 2016. OCLC: 943721122.

34 Seyed Mohammadhossein Tabatabaei, Marc Boyer, Jean-Yves Le Boudec, and Jörn Migge. Efficient and accurate handling of periodic flows in time-sensitive networks. In *2023 IEEE 29th Real-Time and Embedded Technology and Applications Symposium (RTAS)*, pages 303–315, 2023. doi:10.1109/RTAS58335.2023.00031.

35 Pierre Talbot, Éric Monfroy, and Charlotte Truchet. Modular Constraint Solver Cooperation via Abstract Interpretation. *Theory and Practice of Logic Programming*, 20(6):848–863, 2020. doi:10.1017/S1471068420000162.

36 Dhananjay Thiruvady, I. Moser, Aldeida Aleti, and Asef Nazari. Constraint Programming and Ant Colony System for the Component Deployment Problem. *Procedia Computer Science*, 29:1937–1947, 2014. doi:10.1016/j.procs.2014.05.178.

37 Tindell, Hansson, and Wellings. Analysing real-time communications: controller area network (CAN). In *Proceedings Real-Time Systems Symposium REAL-94*, pages 259–263, San Juan, Puerto Rico, 1994. IEEE Comput. Soc. Press. doi:10.1109/REAL.1994.342710.

38 K. W. Tindell, A. Burns, and A. J. Wellings. Allocating hard real-time tasks: An NP-Hard problem made easy. *Real-Time Systems*, 4(2):145–165, June 1992. doi:10.1007/BF00365407.

39 P. M. Yomsi, D. Bertrand, N. Navet, and R. Davis. Controller Area Network (CAN): Response Time Analysis with Offsets. In *9th IEEE International Workshop on Factory Communication Systems*, pages 43–52, United States, May 2012. IEEE. doi:10.1109/WFCS.2012.6242539.

A MiniZinc Model

We describe the full MINIZINC constraint model implementing the mathematical model given in Section 2. We first give the parameters of the model with the corresponding mathematical notations in blue comments:
% I. The hardware graph $\langle H, L, hc, lc \rangle$.

int: locations;
set of int: LOCATIONS = 1..locations; % H
int: numLinks;
set of int: NUM_LINKS = 1..numLinks; % L
array [LOCATIONS] of int: cpu_capacity; % hc
array [NUM_LINKS] of int: capacity; % lc

% II. The software graph $\langle S, Com, sc, cc \rangle$.

% Com is implicitly represented by the adjacency matrix coms where
coms[s_i][s_j] = 0 if the services s_i and s_j do not communicate.

int: services;
set of int: SERVICES = 1..services; % S
array [SERVICES] of int: services_cpu_usage; % sc
array [SERVICES, SERVICES] of int: coms; % cc

% III. The path function

% shortest_path[h_i, h_j] contains all the edges belonging to the shortest path
between h_i and h_j.

% Interestingly, we do not need to know the order of the edges on the
shortest path, thus we can use a set.

array [LOCATIONS, LOCATIONS] of set of NUM_LINKS: shortest_path;

% IV. Not part of the mathematical specification: this is to display the
solutions with locations and services names instead of indexes.

array [LOCATIONS] of string: locations2names;

array [SERVICES] of string: services2names;

The decision variable is the function $d: S \to H$ which is modelled as a MiniZinc array:

array [SERVICES] of var LOCATIONS: services2locs; % d: S \to H

The constraints are defined using intermediate arrays of variables to simplify their definitions.

% I. CPU load constraint.

% $\forall h \in H, \sum_{s \in d^{-1}(h)} sc(s) \le hc(h)$

array [LOCATIONS] of var int: cpu_usage;
constraintforall (1 in LOCATIONS)
(cpu_usage[1] =
  sum (s in SERVICES)
  (services_cpu_usage[s] * (services2locs[s] == 1)))
;
constraintforall (1 in LOCATIONS)(cpu_usage[1] >= 0 /\ cpu_usage[1] <= cpu_capacity[1]);

% II. Network load constraint.

% $\forall l \in L, \sum_{c \in Com} com(c, l) \le lc(l)$

array [NUM_LINKS] of var int: slack;
constraintforall (link in NUM_LINKS) (
  slack[link] = capacity[link] -
  sum (s1, s2 in SERVICES) (
    coms[s1, s2] * (link in shortest_path[services2locs[s1],
      services2locs[s2]])
  ))
);

% Then we ensure the slack is always greater or equal to 0.

constraintforall (link in NUM_LINKS)(slack[link] >= 0 /\ slack[link] <= capacity[link]);
# Efficient Enumeration of Fixed Points in Complex Boolean Networks Using Answer Set Programming

Van-Giang Trinh<fnref n="1" /> ⬒ ID

LIS, Aix-Marseille University, Marseille, France

Belaid Benhamou ⬒

LIS, Aix-Marseille University, Marseille, France

Sylvain Soliman<fnref n="1" /> ⬒ ID

Lifeware team, Inria Saclay, Palaiseau, France

## Abstract

Boolean Networks (BNs) are an efficient modeling formalism with applications in various research fields such as mathematics, computer science, and more recently systems biology. One crucial problem in the BN research is to enumerate all fixed points, which has been proven crucial in the analysis and control of biological systems. Indeed, in that field, BNs originated from the pioneering work of R. Thomas on gene regulation and from the start were characterized by their asymptotic behavior: complex attractors and fixed points. The former being notably more difficult to compute exactly, and specific to certain biological systems, the computation of stable states (fixed points) has been the standard way to analyze those BNs for years. However, with the increase in model size and complexity of Boolean update functions, the existing methods for this problem show their limitations. To our knowledge, the most efficient state-of-the-art methods for the fixed point enumeration problem rely on Answer Set Programming (ASP). Motivated by these facts, in this work we propose two new efficient ASP-based methods to solve this problem. We evaluate them on both real-world and pseudo-random models, showing that they vastly outperform four state-of-the-art methods as well as can handle very large and complex models.

**2012 ACM Subject Classification** Computing methodologies → Logic programming and answer set programming; Applied computing → Computational biology; Applied computing → Systems biology

**Keywords and phrases** Computational systems biology, Boolean network, Fixed point, Answer set programming

Digital Object Identifier 10.4230/LIPIcs.CP.2023.35

**Supplementary Material** Software (Source Code): https://github.com/giang-trinh/fASP
archived at swh:1:dir:86a458f51d92dda88a71499b7e6bea7178377da0

**Funding** Van-Giang Trinh: This work is supported by Institut Carnot STAR, Marseille, France.

## 1 Introduction

In molecular biology, the regulation of the transcription of a gene is the process by which a cell modulates the conversion of its DNA into RNA. Transcription is a vital process in all living organisms and leads to orchestrating the whole gene activity. Its regulation can take many different forms, mostly affecting the binding of the RNA polymerase on the DNA.

The lack of precise quantitative information about transcriptional regulation and the sigmoid nature of its kinetics led, about fifty years ago, to the idea to represent models of gene regulation as discrete event systems. Those gene regulation networks use thresholds or equivalently logical functions to represent the different regulations [22, 39, 41, 40]. Over the years, Boolean Network (BN) modelling has proven that it can bring powerful analyses and

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

<footnote n="1">¹ Corresponding authors</footnote>
corresponding insight to the many cases where enough quantitative biological data is not
available [48], even for modelling post-transcriptional mechanisms. This is even more true
for very large models where such data is frequently missing and led to a constant increase
in size of logical models à la Thomas [2] and more and more complex logical formulae to
describe the dynamics of those models.

Besides simulation, the analysis of such models is mostly based on attractor computation, since those correspond roughly to observable biological phenotypes [48]. An attractor of a BN is a minimal set of states from which the dynamics of this BN cannot escape once entered [39, 48]. An attractor of size one is called a stable state or *fixed point*. Otherwise, it is called a cyclic attractor or complex attractor. To date, the analysis of the set of fixed points of a BN remains a very useful tool in understanding the behavior of those complex biological models. This is not only due to the fact that in some cases the full computation of complex attractors remains intractable, but also because for many biological systems, the expected long-term behavior is not cyclic (as in the Cell Cycle, or Circadian rhythms for instance) but rather a stabilization to an observable *phenotype* (cell differentiation, apoptosis, proliferation, signal transduction, protein transcription, etc.). See for instance [33, 15, 13, 43] for some recent publications using stable states as main validation. It is also worth noting that the fixed point computation is the crucial starting point for several state-of-the-art methods for computing complex attractors of BNs [21, 42].

Answer Set Programming (ASP) [19] has been widely applied in the field of computa-
tional systems biology [46] because of its declarative characteristics as well as strong tools'
support [18]. Very early, ASP has been used to model biological networks [14, 37]. Since BNs
have become a popular modeling formalism in systems biology, it is naturally that ASP has
been quickly applied to modeling and analysis of BNs. One of the first connections between
ASP and BNs is the theoretical work by [26], but nowadays we can find in the literature many
references showing the successful application of ASP to model and reason over biological
systems modeled as BNs. The notable use of ASP in the analysis of BNs ranges from
enumerating fixed points [28, 1, 34], enumerating or approximating attractors [30, 28, 1, 34],
and inferring BNs from biological data [35, 46, 47, 12], to controlling BNs [27, 47].

There is a rich history of research on enumerating fixed points of BNs since this modeling formalism was proposed [22, 39]. The fixed point enumeration problem has attracted researchers from various communities and many methods have been proposed [29]. We can classify the existing methods into the following main approaches: algorithmic [29], structure-based [10, 45, 25, 7], Boolean resolution-based [24, 32, 31], integer linear programming-based [5], and ASP-based [28, 1, 34]. A more detailed summary of the existing methods shall be given in Section 3. Note however that with the increase in model size and complexity of Boolean update functions, the existing methods for this problem show their limitations [29]. One reason is that they require an intermediate representation of the original BN that may be computationally expensive or even intractable to obtain, e.g., prime implicants [28], transition-based representations [1], disjunctive normal forms [34].

Inspired by the above elements along with the fact that the most recent and most efficient
fixed point enumeration methods all rely on ASP, in this work we propose two new ASP-based
methods for efficiently enumerating all fixed points of a BN. The first method is based on
conjunctive ASP, and the second method is a modification of the first one to handle the case
of a large number of source nodes. If a BN has many source nodes, its number of fixed points
may be extremely large, leading to both long running time and high memory consumption.
The main advantage of the two proposed methods is that they rely on negative normal forms
of Boolean functions whose computation is more efficient than that of other intermediate
representations used by the previous methods. After some preliminaries on BNs and the problem of enumerating their fixed points, we present the two new ASP-based methods and benchmark them against four other state-of-the-art tools. The experimental results on both real-world and pseudo-random models show that they vastly outperform the state-of-the-art and can handle very large and complex models.

# 2 Preliminaries

We start by recalling some classical definitions.

## 2.1 Boolean networks

**► Definition 1 (Boolean network).** A Boolean Network (BN) [40] is a pair $\mathcal{N} = (V, F)$ where:

- $V = \{v_1, \dots, v_n\}$ is the set of nodes. We use $v_i$ to denote both the node $v_i$ and its associated Boolean variable.

- $F = \{f_1, \dots, f_n\}$ is the set of Boolean update functions. Each function $f_i$ is associated with node $v_i$ and satisfies $f_i: \mathbb{B}^{|IN(v_i)|} \to \mathbb{B}$ where $\mathbb{B} = \{0, 1\}$ and $IN(v_i)$ denotes the set of input nodes of $v_i$. If $v_j \in IN(v_i)$, we say that there is a regulation between $v_j$ and $v_i$, and $v_j$ is a regulator of $v_i$. Note that a node $v_i \in V$ is called a source node if and only if $f_i$ is an identity function on Boolean variable $v_i$ (i.e., $f_i = v_i$).

**► Example 2.** We give a BN $\mathcal{N} = (V, F)$, where $V = \{v_1, v_2\}$ and $F = \{f_1, f_2\}$ with $f_1 = (v_1 \land v_2) \lor (\neg v_1 \land \neg v_2)$, $f_2 = (v_1 \land v_2) \lor (\neg v_1 \land \neg v_2)$.

**► Definition 3 (local monotonicity).** A Boolean function is locally-monotonic if it can be represented by a formula in Disjunctive Normal Form (DNF) in which all occurrences of any given literal are either negated or non-negated [34].

A BN is said to be locally-monotonic if all its update functions are locally-monotonic. Otherwise, this model is said to be non-locally-monotonic.

The BN of Example 2 is non-locally-monotonic.

## 2.2 Fixed points

A state $x \in \mathbb{B}^n$ is as a mapping $x: V \mapsto \mathbb{B}$ that assigns either 0 (inactive) or 1 (active) to each node. We also write $x_i$ to denote $x(v_i)$ for short and for simplicity we write $f_i(x)$ even when $IN(v_i) \subsetneq V$, i.e., $IN(v_i)$ does not contain some nodes of $V$.

**► Definition 4 (fixed point).** A fixed point of $\mathcal{N}$ is a state $s$ such that $s_i = f_i(s)$ for every $v_i \in V$.

The state space of the BN of Example 2 includes four states: 00, 01, 10, and 11. However, this BN has only one fixed point: 11.

### Complexity

Note that the fixed points do not depend on the choice of *update scheme* of the BN, they are the same for synchronous, asynchronous or even generalized updates [20]. These *update schemes* are what precisely defines a transition relation on states from the update functions. In general they allow one (asynchronous), all (synchronous) or any number (generalized) of nodes to change their value $v_i$ to their update value specified by $f_i$. However, if running a simulation in the synchronous update is a feasible way to find *the only reachable fixed point* starting from a completely known initial state, this does not scale up to big networks with many source nodes with unknown values, or to other update schemes.
From a theoretical view point, the problems of detecting a fixed point and enumerating all fixed points of a general BN have been shown to be respectively NP-hard and #P-hard [3]. In fact, for general BNs, there is no existing method that works faster than $k \times 2^n$ for any $k \ge 1$ [29].

3 Related work

The recent review of Mori and Akutsu [29] shows quite well that because of the above complexity result, a certain amount of work has been done on restricted versions of the problem, limiting the type of BN to simpler regulations like [6, 23], or simpler Boolean formulae like for instance nested canalizing functions [4]. However, when working with real-world models, built by biologists, such restrictions are often impossible to enforce.

Hence, various methods exploiting the structure of a BN have been proposed, using feedback vertex sets [3, 7], subspaces [10], graph-reductions for low connectivity [45], network decomposition [8, 25], etc. Unfortunately, these still do not scale to the size of the most recent BNs (above 1000 nodes) with average connectivity and complex logical formulae.

Other methods with broader generality are also common, using classical Boolean resolution techniques, like BDD or SAT. That is the case for instance of the *BioLQM* library [31] that is at the core of the GINsim Boolean modelling tool [24] and of the CoLoMoTo Docker images [32]. Since integer linear programming is another useful method to efficiently solve Boolean constraints, it has been applied to addressing the fixed point enumeration [5]. The evaluation in [5] shows that this method can handle well models of up to 200 nodes with small average connectivity.

However, the most recent and most efficient fixed point enumeration methods all rely on ASP [19]. This is probably due to the fact that it links the efficiency of SAT for the Boolean constraint solving, having adapted and implemented some techniques like lazy clause generation to the point of winning certain categories of the SAT competition, and the ease of enumeration of all solutions, which is crucial here, in a declarative language.

More precisely, while there is indeed a direct encoding of the fixed-point problem into SAT [29], it creates two issues. First a SAT solver needs to convert the original Boolean formula into a CNF. It is of course possible to use a polynomial transformation like our conjunctive ASP encoding (see Section 4) or Tseitin's transformation, but this introduces auxiliary variables. This in turn leads to enumerating models that encode no fixed points or other redundant models that encode the same fixed points. A step to eliminate spurious and redundant SAT models is therefore necessary to guarantee the correctness and this would add complexity to the SAT/CP approach. In contrast, the ASP approach can avoid the above issue because of the stable model semantics, i.e., only searching for minimal Herbrand models, since the set of Herbrand models one-to-one corresponds to the set of SAT/CP models, whereas the set of minimal ones one-to-one corresponds directly to the set of fixed points.

One of the first connections between ASP and BNs is the theoretical work by [26], but nowadays ASP is used for many different BN analyses, from computing fixed point as we will show, to trap-spaces [28] which are an approximation of complex attractors, and even for representing sets of BNs [12].

By constraining the number of ground atoms in a stable model, the trap-space computation method [28] was adapted to compute fixed points. Note however that this method still requires to compute prime implicants of a Boolean function, and the number of prime implicants may be exponential in the number of inputs of this function. Moreover, the computation of prime implicants of a Boolean function is also a computationally demanding task, and gets intractable when the number of source nodes exceeds 10.
It is worth noting that Paulevé et al. [34] have proposed a new method for computing minimal trap spaces/fixed points that can avoid computing prime implicants. This method has been implemented in the tool mpbn<fnref n="2" /> demonstrated in [34] for handling medium-sized models from the literature and very large synthetic models (up to 100,000 nodes) with respect to minimal trap spaces. Although there is no benchmark designed for the fixed point computation in [34], mpbn should handle well very large models with respect to fixed points. However, there are two drawbacks limiting the applicability of mpbn. First, it requires that the original BN is locally-monotonic. The class of locally-monotonic BNs is too small as compared to the class of all possible BNs, since a BN is non-locally-monotonic if just one of its Boolean functions is non-locally-monotonic (see Definition 3). Moreover, we also found many non-locally-monotonic Boolean models in the literature (see Section 5 for some of them). Second, it requires a DNF of a Boolean function. Note that obtaining a single DNF may be exponential in the size of the Boolean function (i.e., the number of inputs |IN| of this Boolean function).

Not using the concept of trap spaces, the method by [1] characterizes fixed points of a Boolean network as dead configurations (or deadlocks) of its corresponding Automata Network (AN). ANs are formal models similar to Petri nets with transitions representing the updates of the whole system. A transition includes the current configuration, the next configuration, and the condition for enabling this transition. A configuration is said to be dead if and only if there is no transition whose enabling condition is satisfied by this configuration. Then the above characterization is encoded as an ASP. This method has been reported to be able to handle well large-scale models [1]. However, its bottleneck lies in the construction of the corresponding AN, which in general requires to obtain two DNFs for each Boolean variable, one for the update function $f_i$ and one for its negation $\neg f_i$.

# 4 Answer set programming-based methods

We will now describe the two new ASP-based encodings that we propose for the fixed point enumeration problem.

## 4.1 Conjunctive encoding

Let $\mathcal{N} = (V, F)$ be a BN. We intend to build an ASP encoding for $\mathcal{N}$ such that a stable model of the encoded ASP (say $\mathcal{L}$) is equivalent to a fixed point of $\mathcal{N}$. First, for each node $v_i$, we introduce two atoms $p_i$ and $n_i$. The translation from a stable model $A$ of $\mathcal{L}$ to a state $x$ of $\mathcal{N}$ is that for every $v_i \in V$, $x_i = 1$ if and only if $p_i \in A$, and $x_i = 0$ if and only if $n_i \in A$. The below ASP rules ensure that a stable model of $\mathcal{L}$ corresponds to a state of $\mathcal{N}$:

$$ :- p_i, n_i. \qquad (1) $$

meaning $false \leftarrow p_i \wedge n_i$, and

$$ p_i, n_i. \qquad (2) $$

meaning $p_i \lor n_i \Leftarrow true$, for every $v_i \in V$. Recall that state $x$ is a fixed point of $\mathcal{N}$ if and only if the relation $x_i = f_i(x)$ holds for all $v_i \in V$. This relation can be seen as the conjunction of $x_i \leftarrow f_i(x)$ and $\neg x_i \leftarrow \neg f_i(x)$, which can be characterized by $v_i \leftarrow f_i$ and $\neg v_i \leftarrow \neg f_i$, respectively. Hereafter, we show how to encode the two parts for every $v_i \in V$ as conjunctive ASP rules.

<footnote n="2">² https://github.com/bnediction/mpbn</footnote>
35:6 Efficient Enumeration of Fixed Points in Complex Boolean Networks

For the former part, to avoid the presence of negation, we convert $f_i$ into its Negative Normal Form (NNF). The NNF is obtained by recursively applying De Morgan laws until all negations that remain are on literals. We now associate ASP rules to $f_i$ as follows:

$$ \gamma(v_i) :- \gamma(\text{NNF}(f_i)). $$

where we define function $\gamma$ as

$$ \begin{array}{l} \gamma(v_i) = p_i \\ \gamma(\neg v_i) = n_i \\ \gamma(\bigwedge_{1 \le j \le J} \alpha_j) = \gamma(\alpha_1), \dots, \gamma(\alpha_J) \\ \gamma(\bigvee_{1 \le j \le J} \alpha_j) = aux_k \text{ where } aux_k \text{ is a new atom and for each } j \text{ add the rule } aux_k :- \gamma(\alpha_j). \end{array} $$

Note that $k$ is here a global counter starting from 1 and will be increased by 1 after a new atom is created. For the latter part, we similarly apply the above process with $\neg v_i$ and $\neg f_i$ instead of respectively $v_i$ and $f_i$. Note that it also requires to convert $\neg f_i$ into an NNF first.

Listing 1 shows the encoded ASP of the BN shown in Example 2 following the above encoding. Atoms $p_1$ and $n_1$ (resp. $p_2$ and $n_2$) correspond to node $v_1$ (resp. $v_2$). Lines 1 and 2 represent the rules shown in Equation (1) and Equation (2), respectively. The rules for the part $v_1 \leftarrow f_1$ (resp. $\neg v_1 \leftarrow \neg f_1$) are presented in Lines 4–5 (resp. Lines 6–8). Similarly, Lines 10–14 represent the rules for node $v_2$. Line 16 indicates that we omit auxiliary atoms in the resulting stable models. This ASP has only one stable model: $\{p_1, p_2\}$, which corresponds to the sole fixed point (i.e., 11) of $\mathcal{N}$.

<table><caption>Listing 1 Conjunctive ASP encoding for the BN shown in Example 2.</caption><tr><td>:- p1, n1.</td><td>:- p2, n2.</td><td>1</td></tr><tr><td>p1, n1.</td><td>p2, n2.</td><td>2</td></tr><tr><td>p1 :- aux1.</td><td></td><td>3</td></tr><tr><td>aux1 :- p1, p2.</td><td>aux1 :- n1, n2.</td><td>4</td></tr><tr><td>n1 :- aux2, aux3.</td><td></td><td>5</td></tr><tr><td>aux2 :- n1.</td><td>aux2 :- n2.</td><td>6</td></tr><tr><td>aux3 :- p1.</td><td>aux3 :- p2.</td><td>7</td></tr><tr><td>p2 :- aux4.</td><td></td><td>8</td></tr><tr><td>aux4 :- p1, p2.</td><td>aux4 :- n1, n2.</td><td>9</td></tr><tr><td>n2 :- aux5, aux6.</td><td></td><td>10</td></tr><tr><td>aux5 :- n1.</td><td>aux5 :- n2.</td><td>11</td></tr><tr><td>aux6 :- p1.</td><td>aux6 :- p2.</td><td>12</td></tr><tr><td>#show p1/0. #show n1/0.</td><td>#show p2/0. #show n2/0.</td><td>13</td></tr><tr><td>#show p1/0. #show n1/0.</td><td>#show p2/0. #show n2/0.</td><td>14</td></tr><tr><td>#show p1/0. #show n1/0.</td><td>#show p2/0. #show n2/0.</td><td>15</td></tr><tr><td>#show p1/0. #show n1/0.</td><td>#show p2/0. #show n2/0.</td><td>16</td></tr></table>

We here discuss the advantages of the above ASP encoding. First, the ASP $\mathcal{L}$ has no negation besides Equation (1) that may hinder the efficiency of ASP solvers. Note also that obtaining the NNF of a Boolean function is linear in its size and thus quite efficient. Second, except the rules of Equation (2), all the rules in $\mathcal{L}$ are conjunctive. This is the reason we name the above encoding as the *conjunctive* encoding.

The next result shows that our ASP encoding is sound and complete with respect to the fixed points of a BN.
**► Proposition 5.** The set of stable models of $\mathcal{L}$ one-to-one corresponds to the set of fixed points of $\mathcal{N}$.

**Proof.** Note that the only negations in the conjunctive encoding come from the state constraints of Equation (1), i.e., $:-p_i, n_i$. Since this equation is coupled with Equation (2), all stable models will contain exactly one of $p_i$ or $n_i$, i.e., they will correspond with the state $x$ where $x_i$ is true or false depending on the atom in the stable model.

The remainder of the ASP has no negation, no disjunction in heads and is logically equivalent to $\forall v_i \in V, x_i = f_i(x)$ via the introduction of some existentially quantified $aux_j$. Hence there is for each state $x$ at most a single stable model, equal to the smallest Herbrand model containing the $p_i$ or $n_i$ corresponding to $x$, the necessary $aux_j$, and that satisfies $x = f(x)$ by construction.

All stable models of $\mathcal{L}$ will therefore correspond one to one with states $x$ described by the $p_i$ or $n_i$ that are true, and such that $\forall v_i \in V, x_i = f_i(x)$, hence they correspond one to one with fixed points of $\mathcal{N}$.

## 4.2 Source encoding

Recall that the number of fixed points of a BN may be extremely large if it has many source nodes (see Definition 1). Specifically, that number may be exponential in the number of source nodes. In the conjunctive encoding as well as those of the state-of-the-art methods [28, 1, 34], a resulting stable model always corresponds to a single fixed point. Hence, having many source nodes is actually a bottleneck for these methods. To overcome this issue, we propose a new encoding based on the conjunctive encoding of Section 4.1.

Let $\mathcal{L}_c$ be the encoded ASP of the BN following the conjunctive encoding. Let $V^s$ be the set of source nodes of the BN. Our main idea is to group two stable models $A_1$ and $A_2$ of $\mathcal{L}_c$ into a stable model $A$ if they only differ in the atoms corresponding to a source node. More specifically, if there is a source node $v_i$ such that $p_i \in A_1, n_i \in A_2$, and $A_1 \setminus \{p_i\} = A_2 \setminus \{n_i\}$, then we can group $A_1$ and $A_2$ into a stable model $A$ such that $A = A_1 \cup \{n_i\} = A_2 \cup \{p_i\}$. For example, let $A_1$ and $A_2$ be the stable models respectively corresponding to fixed points 01 and 11 of the BN shown in Example 7. Herein, $A_1 = \{n_1, p_2\}$ and $A_2 = \{p_1, p_2\}$. They can be grouped into stable model $A = \{p_1, n_1, p_2\}$. Now, we add $A$ to the set of stable models of $\mathcal{L}_c$, and then repeat the grouping process until there is no new stable model. Note that this process introduces more stable models than before, e.g., we need to consider all $A, A_1$, and $A_2$. However, the new stable model covers all the fixed points represented by the two stable models constituting it. Hence, we just need to consider the maximal set-inclusion stable models. We adjust the conjunctive encoding to make the above approach fully automated in the ASP solver. Since the new encoding aims to handle the case of many source nodes, we name it the *source encoding*.

Similar to the conjunctive encoding, for each node $v_i \in V$, we introduce two atoms $p_i$ and $n_i$. For each node in $V$, we associate to this node the ASP rules identical to whose of the conjunctive encoding. For each $v_i \in V^s$, we remove from the encoded ASP (say $\mathcal{L}_s$) the rule of Equation (1), i.e., $:-p_i, n_i$. By releasing this condition, $\mathcal{L}_s$ can have Herbrand models that contain both $p_i$ and $n_i$, $v_i \in V^s$. To make such Herbrand models to be stable models of $\mathcal{L}_s$, we add to $\mathcal{L}_s$ the choice rules $\{p_i\}$. and $\{n_i\}$. for all $v_i \in V^s$. A choice rule $\{p_i\}$. is equivalent to the rule $p_i$:- not $p_i$. where not denotes the default negation. Finally, we add to $\mathcal{L}_s$ the rules \#show $p_i/0$ and \#show $n_i/0$ for all $v_i \in V$, which indicate that we omit auxiliary atoms in the resulting stable models.
35:8 Efficient Enumeration of Fixed Points in Complex Boolean Networks

Note that a stable model of $\mathcal{L}_s$ may correspond to multiple fixed points of the BN. Given a stable model $A$ of $\mathcal{L}_s$, the set $F$ of fixed points represented by $A$ is specified as follows. For each node $v_i \in V^s$, it can receive value 0 if $n_i \in A$ and $p_i \notin A$, value 1 if $p_i \in A$ and $n_i \notin A$, both if $n_i \in A$ and $p_i \in A$. For each node $v_i \in V \setminus V^s$, it can receive value 0 if $n_i \in A$ and value 1 if $p_i \in A$. Then, $F$ is equivalent to the set of all possible value combinations of all nodes.

The next result shows that our ASP encoding is correct.

**Proposition 6.** The set of maximal set-inclusion stable models of $\mathcal{L}_s$ with respect to the shown atoms exactly covers all fixed points of the BN.

**Proof.** First, we see that $\mathcal{L}_s$ still contains all Herbrand models of $\mathcal{L}_c$. However, by releasing the condition of Equation (1) for all source nodes, $p_i$ and $n_i$ can both appear in a Herbrand model of $\mathcal{L}_s$ if $v_i \in V^s$. Assume that $A_1$ and $A_2$ are two stable models of $\mathcal{L}_s$ such that $p_i \in A_1$, $n_i \in A_2$, and $A_1 \setminus \{p_i\} = A_2 \setminus \{n_i\} = B$ for a node $v_i \in V^s$. The ASP rules corresponding to node $v_i$ are tautology. In all the remaining rules, $p_i$ and $n_i$ never appear in the left hand side. Hence, $A = B \cup \{p_i, n_i\}$ is also a Herbrand model of $\mathcal{L}_s$. By introducing the choice rules for all source nodes, such Herbrand models will be stable models of $\mathcal{L}_s$. Hence, the set of all stable models of $\mathcal{L}_s$ is equivalent to the set of stable models obtained by the grouping approach on the set of stable models of $\mathcal{L}_c$. All fixed points represented by a stable model of $\mathcal{L}_s$ are also covered by a maximal set-inclusion stable model of $\mathcal{L}_s$ with respect to the shown atoms (corresponding to nodes in the BN). Hence, the set of all maximal set-inclusion stable models of $\mathcal{L}_s$ exactly covers all fixed points of the BN. ◀

For illustration, consider the BN shown in Example 7. Listing 2 shows the encoded ASP of this BN following the above encoding. Atoms $p_1$ and $n_1$ (resp. $p_2$ and $n_2$) correspond to node $v_1$ (resp. node $v_2$). Line 1 represents the rule of Equation (2). Note that the rule $:-p_1, n_1$ is removed because $v_1$ is a source node. Instead, two choice rules $\{p_1\}$ and $\{n_1\}$ are added in Line 2. Line 3 represents the rules shown in Equation (1) and Equation (2) of node $v_2$. The rules for the parts $v_1 \leftarrow f_1$ and $\neg v_1 \leftarrow \neg f_1$ are presented in Line 5. Similarly, Lines 7–9 represent the rules for node $v_2$. Line 11 indicates that we omit auxiliary atoms in the resulting stable models. The encoded ASP has four stable models including: $\{n_1, n_2\}$ (corresponding to fixed point 00), $\{n_1, p_2\}$ (corresponding to fixed point 01), $\{p_1, p_2\}$ (corresponding to fixed point 11), and $\{p_1, n_1, p_2\}$ (corresponding to fixed points 01 and 11). From these results, we can see that the encoded ASP has two maximal set-inclusion stable models $(\{n_1, n_2\}$ and $\{p_1, n_1, p_2\})$, which cover all the fixed points of the BN.

**Example 7.** We give a BN $\mathcal{N} = (V, F)$, where $V = \{v_1, v_2\}$ and $F = \{f_1, f_2\}$ with $f_1 = v_1, f_2 = v_1 \lor v_2$. $v_1$ is a source node of $\mathcal{N}$. $\mathcal{N}$ has three fixed points: 00, 01, 11.

<table><caption>Listing 2 Source ASP encoding for the BN shown in Example 7.</caption><tr><td>p1, n1.</td><td></td><td>1</td></tr><tr><td>{p1}.</td><td>{n1}.</td><td>2</td></tr><tr><td>:- p2, n2.</td><td>p2, n2.</td><td>3</td></tr><tr><td>p1 :- p1.</td><td>n1 :- n1.</td><td>4</td></tr><tr><td>p2 :- aux1.</td><td></td><td>5</td></tr><tr><td>aux1 :- p1.</td><td>aux1 :- p2.</td><td>6</td></tr><tr><td>n2 :- n1, n2.</td><td></td><td>7</td></tr><tr><td>#show p1/0.</td><td>#show n1/0.</td><td>8</td></tr><tr><td>#show p2/0.</td><td>#show n2/0.</td><td>9</td></tr><tr><td>#show n2/0.</td><td></td><td>10</td></tr><tr><td></td><td></td><td>11</td></tr></table>
## 4.3 Post-processing

The set of maximal set-inclusion stable models of the encoded ASP $\mathcal{L}_s$ with respect to the shown atoms can be seen as a meta result from which we can easily retrieve the fixed points. Note that a stable model can be group-able with multiple ones, and thus the sets of fixed points of two maximal set-inclusion stable models of $\mathcal{L}_s$ can intersect. In other words, a fixed point of the BN may be included in two different maximal set-inclusion stable models of $\mathcal{L}_s$. Hence, it is not straightforward to obtain the number of fixed points, which by itself is a common difficult problem related to #SAT (model-counting, see [44]). Moreover, many more precise analysis questions can be answered, but not directly from the above meta result. For example, given a state $x^n$ on $V \setminus V^s$, return the set of states $x^s$ on $V^s$ such that $(x^n, x^s)$ is a fixed point of the BN. With this motivation, we propose a post-processing step as follows.

We maintain a hash table (denoted by $H$) with as *key* a state $x^n$ on $V \setminus V^s$ and as *associated value* the set of states $x^s$ on $V^s$ such that $(x^n, x^s)$ is a fixed point of the BN. For each stable model $A$ of the meta result, we extract $x^n$ from it by simply checking either $p_i$ or $n_i$ belongs to $A$ for all $v_i \in V \setminus V^s$. For each $v_i \in V^s$, $x_i$ can receive value 0 if $n_i \in A$ and value 1 if $p_i \in A$. Then we get the set of states on $V^s$ (denoted by $S^s$) as the combinations of all possible values of all $v_i \in V^s$. If $x^n$ is not a key of $H$, we just add the pair $(x^n, S^s)$ to $H$. Otherwise, we replace the current associated value of $x^n$ in $H$ by the union of it and $S^s$ because the current associated value and $S^s$ may intersect. When there are many nodes in $V^s$, $S^s$ may contain a large number of states, even exponential in the number of nodes in $V^s$. Binary Decision Diagrams (BDDs) [11] are an efficient data structure for representing a set of states as well as performing set operations. Hence, we store $S^s$ as a BDD where each node $v_i \in V^s$ corresponds to a BDD variable.

Now, let us show how to answer some analysis questions from the hash table $H$. First, for the example analysis question mentioned in the beginning of Subsection 4.3, if $x^n$ is not a key of $H$, we return the empty set. Otherwise, the set of all states on $V^s$ can be easily retrieved from the BDD as the value of $x^n$ in $H$. Specifically, we list all satisfying valuations of this BDD. Analogously, by traversing all items in $H$, we can also answer the question, given a combination of values on source nodes, to return the set of fixed points of the BN restricted by this combination.

Second, we can efficiently compute the number of fixed points in the BN based on $H$. For a given BDD, we can efficiently compute its number of satisfying valuations. Such procedure is linear in the number of nodes in this BDD [11]. Since any two keys in $H$ are distinct, the sets of fixed points on $V$ corresponding to them are also distinct. Hence, the number of fixed points of the BN is equivalent to the sum of all numbers of satisfying valuations for all BDDs in $H$. This is in contrast with most model-counting problems, where approximate methods are usually necessary to minimize the number of calls to a SAT solver.

Note that an alternative approach might be to represent the hash table as a sole BDD where every node of the BN corresponds to a BDD variable. As such, the new BDD will have $|V|$ variables, whereas each BDD in $H$ has $|V^s|$ variables. Recall that the size of a BDD may be exponential in its number of variables [11]. Hence, splitting a BDD into multiple BDDs with smaller numbers of variables is a good strategy in most cases, especially when $|V|$ is much larger than $|V^s|$. Indeed, in our experiments shown in Section 5, we observed that the alternative approach using a sole BDD gave poorer performance than the approach using a hash table in most cases.

To conclude this section, we discuss the advantages of the source encoding. First, it inherits all the advantages of the conjunctive encoding, which can be seen as its core part. Second, the number of maximal set-inclusion stable models of the encoded ASP may be much