Therefore, the slack of LCut($T, u$) is greater than the one of LCut($T, \delta(i)$). The contribution of task $i$ in LCut($T, \delta(i)$) is greater than its contribution in LCut($T, u$). When this contribution is considered during the scheduling of LCut($T, u$), it results in an overload. Therefore, the same overload happens when the contribution of task $i$ is considered during the scheduling of the interval LCut($T, \delta(i)$). Thus, the rule (FTHE-EF) holds for LCut($T, \delta(i)$) and task $i$. ◀

According to Lemmas 4 and 6, the application of rule (FTHE-EF) to the intervals LCut($T, \beta(i)$) and LCut($T, \delta(i)$) is enough to subsume the (extended) edge finding detection rule. The new detection rule is based on the application of the rule (FTHE-EF) on the intervals LCut($T, \beta(i)$) and LCut($T, \delta(i)$) for all $i \in T$. This rule, denoted *Slack-Density Horizontally Elastic Edge Finder* detection, is specified by the formula: for all $i \in T$,

$$ \left\{ \begin{array}{l} \mathrm{ect}_{\mathrm{LCut}(T, \beta(i)) \cup \{i'\}}^{H} > \mathrm{lct}_{\beta(i)} \Rightarrow \mathrm{LCut}(T, \beta(i)) \ll i \\ \mathrm{ect}_{\mathrm{LCut}(T, \delta(i)) \cup \{i'\}}^{H} > \mathrm{lct}_{\delta(i)} \Rightarrow \mathrm{LCut}(T, \delta(i)) \ll i \end{array} \right. \qquad (\mathrm{SDHE-EF}) $$

We denote by (SDHE-EF) the filtering rule resulting from the combination of the detection rule (SDHE-EF) with the adjustment of [15]. This rule subsumes the edge finding rule and its extension combined with the rule (Adj).

**Theorem 7.** *When the detection rule (SDHE-EF) is combined with the adjustment of [15], the resulting filtering rule subsumes the conjunction of rules (EF) and (EEF) combined with the adjustment rule (Adj).*

**Proof.** According to Lemmas 4 and 6, any detection made by (EF) and (EEF) is also done by (SDHE-EF). It is proved in (Theorem 4 of [15]) that the adjustment of [15] is better than the one performed by rule (Adj) after detection made by rules (EF) and (EEF). Therefore, the propagation of (SDHE-EF) combined with the adjustment of [15] subsumes the conjunction of rules (EF) and (EEF) combined with the adjustment rule (Adj). ◀

Back to the CuSP of Figure 1a, after filling the lower part of the resource of capacity 1, the overflow of three units of energy remains. This overflow is scheduled on the upper part of the resource of capacity 1. No energy is consumed in the interval $[1, 2)$ since no task is scheduled in this interval. The scheduling of the upper part of the profile ends at time 4 which corresponds to the same adjustment value made by rule (Adj) with $\Theta = \{b, c\}$. Indeed, $\mathrm{rest}(\Theta, c_d) = 2 > 0$ and $\mathrm{est}_{\Theta} + \mathrm{rest}(\Theta, c_d)/c_d = 4$. It is known from [11] that (GQHE-EF) is a relaxation of (FTHE-EF). Therefore, (GQHE-EF) and (SDHE-EF) are both relaxations of (FTHE-EF) which subsumes the (extended) edge finding rule. In the following theorem, we prove that the rules (GQHE-EF) and (SDHE-EF) are not comparable, i.e., there exists propagation only performed by (GQHE-EF) and propagation only performed by (SDHE-EF).

**Theorem 8.** *The rules (SDHE-EF) and (GQHE-EF) are not comparable.*

**Proof.** Consider the CuSP instance of Figure 2a, where three tasks $\{x, y, z\}$ share a resource of capacity 2. We have $\Omega_{\beta(z)}^{\alpha(z)} = \{x, y\}$ and the rule (SDHE-EF) detects $\{x, y\} \ll z$. This precedence is missed by the rule (GQHE-EF) since, in the profile of $\{x, y\}$, from time 8 back to time 2, we have enough free energy to schedule task $z$.

Converse, we consider the CuSP instance of figure 2b, where four tasks $\{a, b, c, d\}$ share a resource of capacity 2. In the profile of $\{a, b, c\}$, it remains three units of free energy of height 1 when we move from time 8 back to time 2. This free energy is not enough to schedule task $d$ and the rule (GQHE-EF) detects the relation $\{a, b, c\} \ll d$. We have $\Omega_{\beta(d)}^{\alpha(d)} = \{b, c\}$ and the rule (SDHE-EF) misses the relation $\{a, b, c\} \ll d$. ◀
<figure><img src="image_174.png" /><figcaption>(a) CuSP instance of three tasks sharing a resource of capacity 2</figcaption></figure>

<figure><img src="image_175.png" /><figcaption>(b) CuSP instance of four tasks sharing a resource of capacity 2</figcaption></figure>

■ **Figure 2** 2a: (SDHE-EF) detects $\{x, y\} \le z$ while (GQHE-EF) misses the detection; (2b): (GQHE-EF) detects $\{a, b, c\} \le d$ while (SDHE-EF) misses the detection.

The rule (SDHE-EF) is a relaxation of the rule (FTHE-EF), since it restricts the detection to two intervals instead of $n$. It is proved in [11] that the rule (FTHE-EF) is not comparable to the rule (TT-EF). This result remains between (SDHE-EF) and (TT-EF) as it is proved in the following theorem.

► **Theorem 9.** *The rules (SDHE-EF) and (TTEF) are not comparable.*

*Proof.* Consider the CuSP instance of Figure 3a where five tasks $\{a, b, c, d, e\}$ share a resource of capacity 2. We have $\Omega_{\beta(d)}^{\alpha(d)} = \{a, b, c, d\}$ and the rule (SDHE-EF) detects $\{a, b, c, d\} \le e$. This precedence is missed by the rule (TTEF). Indeed, none of the tasks has a fix part and the edge finding rule (EF) fails to detect $\{a, b, c, d\} \le e$.

Conversely, we consider the CuSP instance of figure 3b where six tasks $\{u, v, w, x, y, z\}$ share a resource of capacity 2. When the fix part of task $x$ which intersects the interval $[0, 6]$ is considered, it is detected that $\{u, v, w\} \le z$ and the rule (TTEF) holds. This detection is missed by the rule (SDHE-EF), since the contribution of task $x$ in the interval $LCut(T, u) = \{u, v, w\}$ is not considered. ◀

<figure><img src="image_176.png" /><figcaption>(a) CuSP instance of five tasks sharing a resource of capacity 2</figcaption></figure>

<figure><img src="image_177.png" /><figcaption>(b) CuSP instance of six tasks sharing a resource of capacity 2</figcaption></figure>

■ **Figure 3** 3a: (SDHE-EF) detects $\{a, b, c, d\} \le e$ while (TTEF) misses the detection; 3b: (TTEF) detects $\{u, v, w\} \le z$ while (SDHE-EF) misses the detection.

## 4 Slack-Density Horizontally Elastic Edge Finder Algorithm

The detection algorithm identifies the right bound of the task intervals $\Omega_{\beta(i)}^{\alpha(i)}$ and $\Omega_{\delta(i)}^{\gamma(i)}$ for any task $i \in T$. To do so, the function *ComputesBound()* receives as input the set $T_{est}$ (resp. $T_{lct}$) of tasks sorted in increasing order of est (resp. lct). The global maximum density and minimum slack are initialized at line 2, while the local maximum density and minimum slack are initialized at line 4 of the loop of line 3. The loop of line 5 updates the value of the local maximum density at line 9 and the global value at line 12. Similarly, the loop of line 14 updates the value of the local minimum slack at line 17 and the global value at line 19.

The function *ComputesBound* has a quadratic complexity as it is shown in Proposition 10.

► **Proposition 10.** *ComputesBound* runs in $O(n^2)$ in time.
■ **Algorithm 2** *ComputesBound*(*T*<sub>est</sub>, *T*<sub>lct</sub>) in *O*(*n*<sup>2</sup>) time.

**Input:** *T*<sub>est</sub> and *T*<sub>lct</sub> the sets of tasks sorted in increasing order of est and lct respectively.

**Output:** The bounds $\beta(i)$ and $\delta(i)$ for all task $i \in T$.

1 **forall** $i \in T_{est}$ **do**
2 &nbsp;&nbsp;&nbsp; $\beta(i) \leftarrow -1$, $\delta(i) \leftarrow -1$, $\maxDensity(i) \leftarrow 0$, $\minSlack(i) \leftarrow \infty$
3 **forall** $j \in T_{lct}$ with $lct_j < lct_{j+1}$ **do**
4 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $E \leftarrow 0$, $\maxD \leftarrow 0$, $\minS \leftarrow \infty$
5 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **forall** $k \in T_{est}$ in reverse order of est **do**
6 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **if** $lct_k \le lct_j$ **then**
7 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $E \leftarrow E + e_k$, $\text{density} \leftarrow E/(lct_j - est_k)$
8 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **if** $\text{density} \ge \maxD$ **then**
9 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\maxD \leftarrow \text{density}$
10 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **else**
11 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **if** $\maxD \ge \maxDensity(k)$ **then**
12 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\maxDensity(k) \leftarrow \maxD$, $\delta(k) \leftarrow j$
13 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\text{Energy}(k) \leftarrow E$
14 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **forall** $k \in T_{est}$ **do**
15 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\text{slack} \leftarrow C(lct_j - est_k) - \text{Energy}(k)$
16 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **if** $\text{slack} < \minS$ **then**
17 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\minS \leftarrow \text{slack}$
18 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **if** $lct_k > lct_j \land \minS < \minSlack(k)$ **then**
19 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\minSlack(k) \leftarrow \minS$, $\beta(k) \leftarrow j$

**Proof.** The loop of line 3 of complexity $\mathcal{O}(n)$ calls sequentially the loops of lines 5 and 14 of complexity $\mathcal{O}(n)$ each. Therefore, the overall complexity of *ComputesBound* is $\mathcal{O}(n(n+n)) = \mathcal{O}(n^2)$. ◀

In Algorithm 3, the bounds of the task interval of minimum slack and maximum density are computed at line 1. In the loop of line 2, for any unscheduled task (line 3), it is checked at line 6 whether the task interval of minimum slack precedes the task and the relation is recorded at line 7. When no detection is made by the task interval of minimum slack (line 8), the task interval of maximum density is used to check the relation at line 10 and the relation is recorded at line 11.

The complexity of *Slack-Density Horizontally Elastic Edge Finder* is analyzed in Proposition 11.

► **Proposition 11.** *Slack-Density Horizontally Elastic Edge Finder runs in *O*(n<sup>2</sup>) in time.*

**Proof.** The linear function *ScheduleTasks* is sequentially called twice in the linear loop of line 2. Combined with the quadratic complexity of the function *ComputesBound*, the overall complexity of *Slack-Density Horizontally Elastic Edge Finder* is $\mathcal{O}(n^2+n(n+n)) = \mathcal{O}(n^2)$. ◀

We combine this detection algorithm with the quadratic adjustment algorithm of [15] proposed for the rule (GQHE-EF). Therefore, the complexity of the two-phase algorithm (Detection and Adjustment) is $\mathcal{O}(n^2)$. To our knowledge, this is the first quadratic horizontally elastic edge finder algorithm.
**Algorithm 3** Slack-Density Horizontally Elastic Edge Finder($T_{est}$, $T_{lct}$) in $O(n^2)$ time.

**Input:** $T_{est}$ and $T_{lct}$ the sets of tasks sorted in increasing order of est and lct respectively.

**Output:** Prec: the precedence relation.

1. $\beta, \delta \leftarrow \text{ComputesBound}(T_{est}, T_{lct})$

2. **forall** $i \in T_{lct}$ **do**

<table><tr><td>3</td><td>if $ect_i < lct_i$ **then**</td></tr><tr><td>4</td><td>&nbsp;&nbsp;&nbsp;&nbsp;if $\beta(i) \neq -1$ **then**</td></tr><tr><td>5</td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;ect &larr; ScheduleTask(LCut(T, $\beta(i)$) &cup; {$i'$})</td></tr><tr><td>6</td><td>&nbsp;&nbsp;&nbsp;&nbsp;if $ect > lct_{\beta(i)}$ **then**</td></tr><tr><td>7</td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Prec(i) &larr; $\beta(i)$</td></tr><tr><td>8</td><td>&nbsp;&nbsp;&nbsp;&nbsp;if Prec(i) &ne; -1 &and; $\delta(i) \neq -1$ **then**</td></tr><tr><td>9</td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;ect &larr; ScheduleTask(LCut(T, $\delta(i)$) &cup; {$i'$})</td></tr><tr><td>10</td><td>&nbsp;&nbsp;&nbsp;&nbsp;if $ect > lct_{\delta(i)}$ **then**</td></tr><tr><td>11</td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Prec(i) &larr; $\delta(i)$</td></tr></table>

<figure><img src="image_178.png" /><figcaption>**Figure 4** A CuSP instance of five tasks sharing a resource of capacity 3.</figcaption></figure>

**► Example 12.** Consider the CuSP instance of Figure 4 where five tasks share a resource of capacity 3.

In this example, $\Omega_u^u = \{u, v\}$, $sl(\Omega_u^u) = 5$ while $\Omega_w^u = \{u, v, w, x\}$, $sl(\Omega_w^u) = 6$. Therefore, $\Omega_{\beta(y)}^{\alpha(y)} = \{u, v\}$ and the rules (EF) and (SDHE-EF) detect that $\Omega_{\beta(y)}^{\alpha(y)} \le y$.

## 5 Improvements

To improve our algorithm, we have taken into account the fix part of tasks $T \setminus \text{LCut}(T, \beta(i)) \cup \{i\}$ (resp. $T \setminus \text{LCut}(T, \delta(i)) \cup \{i\}$) which overlap with $\text{LCut}(T, \beta(i))$ (resp. $\text{LCut}(T, \delta(i))$) during the computation of $\text{ect}_{\text{LCut}(T, \beta(i)) \cup \{i'\}}^H$ (resp. $\text{ect}_{\text{LCut}(T, \delta(i)) \cup \{i'\}}^H$) and the adjustment. For a given task $j \in T \setminus \text{LCut}(T, \beta(i)) \cup \{i\}$ with a fix part (i.e., $\text{lst}_j < \text{lct}_{\beta(i)}$), which overlap with $\text{LCut}(T, \beta(i))$ (i.e., $\text{lst}_j < \text{lct}_{\beta(i)}$), a new task denoted $f(j, \beta(i))$ is deduced with the following attributes $\langle \text{lst}_j, \min(\text{ect}_j, \text{lct}_{\beta(i)}), \min(\text{ect}_j, \text{lct}_{\beta(i)} - \text{lst}_j, c_j \rangle$. We denote by $f(T \setminus \text{LCut}(T, \beta(i)) \cup \{i\}, \beta(i))$ the set of deduced fix parts of tasks which overlap with $\text{LCut}(T, \beta(i))$. Analogously, the set of fix part of tasks $f(T \setminus \text{LCut}(T, \delta(i)) \cup \{i\}, \delta(i))$ which overlap with $\text{LCut}(T, \delta(i))$ is considered during the computation of $\text{ect}_{\text{LCut}(T, \delta(i)) \cup \{i'\}}^H$.

To do so, we first extend the initial time points, by adding those corresponding to the $\text{lst}_i$ for all $i \in T$. The number of time points moves from $4n + 1$ to $5n + 1$ and the function ScheduleTasks remains linear. During the initialization of increments, for any task $j \in T \setminus \text{LCut}(T, \beta(i)) \cup \{i\}$ (resp. $j \in T \setminus \text{LCut}(T, \delta(i)) \cup \{i\}$), if tasks $j$ has a fix part which overlap with $\text{LCut}(T, \beta(i))$ (resp. $\text{LCut}(T, \delta(i))$), $\Delta_{max}$ and $\Delta_{req}$ are increased by $c_j$ at $t$. $\text{lst}_j$. If $\text{ect}_j < \text{lct}_{\beta(i)}$ (resp. $\text{ect}_j < \text{lct}_{\delta(i)}$) then $\Delta_{max}$ and $\Delta_{req}$ are decreased by $c_j$.
at $t.ect_j$. If $ect_j \ge \text{lct}_{\beta(i)}$ (resp. $ect_j \ge \text{lct}_{\delta(i)}$) then $\Delta_{max}$ and $\Delta_{req}$ are decreased by $c_j$ at $t.ect_{\beta(i)}$ (resp. $t.ect_{\delta(i)}$). For the adjustment, the fix part tasks $f(T \setminus \text{LCut}(T, \beta(i))) \cup \{i\}, \beta(i))$ (resp. $f(T \setminus \text{LCut}(T, \delta(i))) \cup \{i\}, \delta(i)))$ is considered during the computation of the maximum overflow obtained when the set $\text{LCut}(T, \beta(i))$ (resp. $\text{LCut}(T, \delta(i)))$ is scheduled on a resource of restricted capacity $C - c_i$. The additional fix part of tasks can increase both the detection and the adjustment of the *Slack-Density Horizontally Elastic Edge Finder* algorithm. For example, in the CuSP of Figure 2b, when the fix part of task *a* is considered, the task interval $\Omega_{\beta(d)}^{\alpha(d)} = \{b, c\}$ can detect that $\Omega_{\beta(d)}^{\alpha(d)} \le d$ and the adjustment follows. It is also the case for the CuSP instance of Figure 3b, when the fix part of task *x* is considered.

# 6 Experimental Results

The new algorithm with improvements was compared to the state-of-the-art horizontally elastic edge finder algorithms ([15, 11]) on Resource-Constrained Project Scheduling Problems (RCPSPs). Comparisons are done on instances of benchmark suites of RCPSP of libraries BL [2] and PSPLib [20] (on sets j30, j60 and j90). Starting time of tasks and makespan was used as variables of the model. They were constrained by the precedence graph and resource limitations. Each resource was modeled with a single CUMULATIVE constraint [1]. The lower bound of the makespan variable was updated with the horizontally elastic earliest completion time of the whole set of tasks *T*. The TimeTabling algorithm of [21] was added to the common core model. The optimization is based on depth-first search and restart. Anytime a solution is found, the resolution restarts with an additional constraint, which states that the next makespan must be (strictly) better than the current one. The optimum solution is the last solution found within the time limit.

Three configurations of the global constraint CUMULATIVE was considered. The first one denoted GQHE-EF uses the horizontally elastic edge-finder from [15], while the second one denoted FTHE-EF considers the horizontally elastic edge-finder algorithm of [11]. The last configuration denoted SDHE-EF is based on the Algorithm 3 for detection, combined with the adjustment algorithm of [15], all with improvements of Section 5.

Three strategies of selection of variables and values were used to speed up the solving process. Static heuristic were unscheduled tasks are selected in the chronological order of the indices and assigned to their lower bound value. COS + DomOvWDeg where the Conflict Ordering Search heuristic (COS) [13] is combined with the default search strategy of Choco solver *domOverWDeg* [5]. COS + Smallest where the COS is combined to the smallest heuristic (a variable of smallest lower bound among those not yet instantiated is selected and assigned to its lower bound). The implementation was done in Java using Choco solver 4.10.8 [25]. Any search taking more than 10 minutes was counted as a failure.

In Table 1, the column “solve” indicates the number of instances solved to optimality by each configuration, “common solve” indicates the numbers of instances each configuration solves commonly with the baseline configuration SDHE-EF. The column “back<” (resp. “back>”) indicates the number of instances where each configuration reduces (resp. need) more backtracks than the baseline configuration. The column “back” (resp. “time”) indicates the average number of backtracks (resp. runtime in sec) on common solved instances for each configuration with the baseline configuration.

SDHE-EF always solves more instances than the other configurations, whatever the heuristic selection considers. On common instances solves with static heuristic by SDHE-EF and GQHE-EF (resp. SDHE-EF and FTHE-EF) 100 (resp. 67) instances were solved by SDHE-EF with fewer backtracks (see Table 1). Figures 5a, 5c and 5e compare the number of
<table><caption>**Table 1** Number of instances solved, commonly solved with the baseline configuration SDHE-EF, where the number of backtracks is less (resp. great) than the one of the baseline. The average number of backtracks and runtime (in sec) are also provided for commonly solved instances. The one of the baseline configuration are in brackets.</caption><thead><tr><th></th><th>solve</th><th>common solve</th><th>back&lt;</th><th>back&gt;</th><th>back</th><th>time</th></tr></thead><tbody><tr><td colspan="7" style="text-align: center;">Static</td></tr><tr><td>GQHE-EF</td><td>1037</td><td>1030</td><td>72</td><td>100</td><td><strong>5569</strong> (5939)</td><td>8.294 (<strong>4.666</strong>)</td></tr><tr><td>FTHE-EF</td><td>1031</td><td>1024</td><td>93</td><td>67</td><td><strong>6997</strong> (7810)</td><td>6.552 (<strong>3.695</strong>)</td></tr><tr><td>SDHE-EF</td><td><strong>1044</strong></td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td colspan="7" style="text-align: center;">COS + DomOvWDeg</td></tr><tr><td>GQHE-EF</td><td>1097</td><td>1093</td><td>213</td><td>238</td><td><strong>4675</strong> (4912)</td><td>9.377 (<strong>5.207</strong>)</td></tr><tr><td>FTHE-EF</td><td>1086</td><td>1085</td><td>209</td><td>219</td><td>3304 (<strong>3209</strong>)</td><td>11.788 (<strong>3.479</strong>)</td></tr><tr><td>SDHE-EF</td><td><strong>1121</strong></td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td colspan="7" style="text-align: center;">COS + Smallest</td></tr><tr><td>GQHE-EF</td><td>1053</td><td>1053</td><td>32</td><td>87</td><td>6219 (<strong>5859</strong>)</td><td>7.976 (<strong>4.489</strong>)</td></tr><tr><td>FTHE-EF</td><td>1049</td><td>1049</td><td>42</td><td>69</td><td>2417 (<strong>2360</strong>)</td><td>10.003 (<strong>3.132</strong>)</td></tr><tr><td>SDHE-EF</td><td><strong>1060</strong></td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

backtracks of GQHE-EF and FTHE-EF with the baseline configuration SDHE-EF on each instance commonly solved. The points above the line $y = x$ shown that SDHE-EF records less backtracks than the other configurations. Figures 5b, 5d and 5f represent the number of solved instances as a function of time for each configuration. The running time gains between SDHE-EF and the two other configurations for static heuristic, is not significant enough as it is illustrated in Figure 5b.

When the heuristic COS + DomOvWDeg is considered, the average number of backtracks of SDHE-EF is 3209, which is less than the average number of backtracks of FTHE-EF (3304), on commonly solved instances with baseline configuration (see Table 1). The average number of backtracks of GQHE-EF (4675) is less than the one of SDHE-EF (4912). The average runtime of SDHE-EF is half of the other configurations whatever the heuristic selection considered (see Table 1). The running time gains between SDHE-EF and the two other configurations is more prominent as illustrated in Figure 5d and 5f.

## 7 Conclusion

In this paper, we have proposed a new strengthening of the edge finder rule based on the *Profile* data structure. The minimum slack and the maximum density are used to select the potential task interval for the edge finder rule. The application of the *Profile* on those task intervals results in a new rule named *Slack-Density Horizontally Elastic Edge Finder*. This new rule subsumes the classic (extended) edge finder rule, but is not comparable to *Gingras* and *Quimper's* horizontally elastic edge finder rule [15], and the *TimeTable* edge finding rule [30]. A quadratic detection algorithm for the new rule is combined with the quadratic adjustment algorithm of [15], which results in an overall complexity of $O(n^2)$ in time. Improvements based on the *TimeTable* data structure are obtained by considering fix part of external tasks during the horizontally elastic scheduling of task intervals. Experimental results showed that our new algorithm is competitive with start-of-the-art strengthening of edge finding algorithms, with *Profile* for time and tree search reduction. It is a good trade-off between the speed and the filtering power for the rule (FTHE-EF) of [11].

Future works will be devoted to the utilization of the *Profile* data structure in the energetic reasoning to reduce its complexity and increase its filtering power.
<figure><img src="image_179.png" /><figcaption>(a) # Backtrack Comparison</figcaption></figure>

<figure><img src="image_180.png" /><figcaption>(c) # Backtrack Comparison</figcaption></figure>

<figure><img src="image_181.png" /><figcaption>(e) # Backtrack Comparison</figcaption></figure>

<caption>(f) Runtime(sec) Comparison</caption>

■ **Figure 5** 5a, 5c and 5e: Comparison of the number of backtracks of different configurations of the CUMULATIVE constraint. SDHE-EF is used as the baseline model. 5b, 5d and 5f: Number of solved instances as a function of time for the different configurations of the CUMULATIVE constraint.

References

1. Abderrahmane Aggoun and Nicolas Beldiceanu. Extending CHIP in order to Solve Complex Scheduling and Placement Problems. *Mathl. Comput. Modelling*, 17(7):57–73, 1993. URL: https://hal.archives-ouvertes.fr/hal-00442821.

2. P. Baptiste, C.L. Pape, and W. Nuijten. *Constraint-Based Scheduling: Applying Constraint Programming to Scheduling Problems*. International Series in Operations Research & Management
20:16 Horizontally Elastic Edge Finder Rule Based on Slack and Density

ment Science. Springer US, 2012. URL: https://books.google.com/books?id=qUzhBwAAQBAJ.

3 Philippe Baptiste, Philippe Laborie, Claude Le Pape, and Wim Nuijten. Constraint-based scheduling and planning. In *Handbook of Constraint Programming*, volume 2 of *Foundations of Artificial Intelligence*, pages 761–799. Elsevier, 2006.

4 Nicolas Beldiceanu and Mats Carlsson. A new multi-resource cumulatives constraint with negative heights. In *CP*, volume 2470 of *Lecture Notes in Computer Science*, pages 63–79. Springer, 2002.

5 Frédéric Boussemart, Fred Hemery, Christophe Lecoutre, and Lakhdar Sais. Boosting systematic search by weighting constraints. In *ECAI*, pages 146–150. IOS Press, 2004.

6 Rajkumar Buyya, M. Manzur Murshed, David Abramson, and Srikumar Venugopal. Scheduling parameter sweep applications on global grids: a deadline and budget constrained cost-time optimization algorithm. *Softw. Pract. Exp.*, 35(5):491–512, 2005.

7 Jacques Carlier, Eric Pinson, Abderrahim Sahli, and Antoine Jouglet. An $O(n^2)$ algorithm for time-bound adjustments for the cumulative scheduling problem. *Eur. J. Oper. Res.*, 286(2):468–476, 2020.

8 Yves Caseau and François Laburthe. Improved CLP scheduling with task intervals. In *ICLP*, pages 369–383. MIT Press, 1994.

9 Hamed Fahimi, Yanick Ouellet, and Claude-Guy Quimper. Linear-time filtering algorithms for the disjunctive constraint and a quadratic filtering algorithm for the cumulative not-first not-last. *Constraints An Int. J.*, 23(3):272–293, 2018.

10 Sévérine Fetgo Betmbe and Clémentin Tayou Djamegni. Horizontally Elastic Edge-Finder Algorithm for Cumulative Resource Constraint Revisited. In *CARI 2020, THIES, Senegal, October 2020*. URL: https://hal.archives-ouvertes.fr/hal-02931383.

11 Sévérine Fetgo Betmbe and Clémentin Tayou Djamegni. Horizontally elastic edge-finder algorithm for cumulative resource constraint revisited. *Oper. Res. Forum*, 3(65), 2022. doi: 10.1007/s43069-022-00172-6.

12 M. R. Garey and David S. Johnson. *Computers and Intractability: A Guide to the Theory of NP-Completeness*. W. H. Freeman, 1979.

13 Steven Gay, Renaud Hartert, Christophe Lecoutre, and Pierre Schaus. Conflict ordering search for scheduling problems. In *CP*, volume 9255 of *Lecture Notes in Computer Science*, pages 140–148. Springer, 2015.

14 Steven Gay, Renaud Hartert, and Pierre Schaus. Simple and scalable time-table filtering for the cumulative constraint. In *CP*, volume 9255 of *Lecture Notes in Computer Science*, pages 149–157. Springer, 2015.

15 Vincent Gingras and Claude-Guy Quimper. Generalizing the edge-finder rule for the cumulative constraint. In *IJCAI*, pages 3103–3109. IJCAI/AAAI Press, 2016.

16 Roger Kameugne, Sévérine Betmbe Fetgo, Vincent Gingras, Yanick Ouellet, and Claude-Guy Quimper. Horizontally elastic not-first/not-last filtering algorithm for cumulative resource constraint. In *CPAIOR*, volume 10848 of *Lecture Notes in Computer Science*, pages 316–332. Springer, 2018.

17 Roger Kameugne and Laure Pauline Fotso. A cumulative not-first/not-last filtering algorithm in $o(n2\log(n))$. *Indian Journal of Pure and Applied Mathematics*, 44:95–115, 2013.

18 Roger Kameugne, Laure Pauline Fotso, and Joseph D. Scott. A quadratic extended edge-finding filtering algorithm for cumulative resource constraints. *International Journal of Planning and Scheduling (IJPS)*, 1(4), 2013.

19 Roger Kameugne, Laure Pauline Fotso, Joseph D. Scott, and Youcheu Ngo-Kateu. A quadratic edge-finding filtering algorithm for cumulative resource constraints. *Constraints An Int. J.*, 19(3):243–269, 2014.

20 Rainer Kolisch and Arno Sprecher. Psplib – a project scheduling problem library. *EUROPEAN JOURNAL OF OPERATIONAL RESEARCH*, 96:205–216, 1996.
21 Arnaud Letort, Nicolas Beldiceanu, and Mats Carlsson. A scalable sweep algorithm for the cumulative constraint. In *CP*, volume 7514 of *Lecture Notes in Computer Science*, pages 439–454. Springer, 2012.

22 Luc Mercier and Pascal Van Hentenryck. Edge finding for cumulative scheduling. *INFORMS J. Comput.*, 20(1):143–153, 2008.

23 Pierre Ouellet and Claude-Guy Quimper. Time-table extended-edge-finding for the cumulative constraint. In *CP*, volume 8124 of *Lecture Notes in Computer Science*, pages 562–577. Springer, 2013.

24 Yanick Ouellet and Claude-Guy Quimper. A $o(n \log^2 n)$ checker and $o(n^2 \log n)$ filtering algorithm for the energetic reasoning. In *CPAIOR*, volume 10848 of *Lecture Notes in Computer Science*, pages 477–494. Springer, 2018.

25 Charles Prud'homme, Jean-Guillaume Fages, and Xavier Lorca. *Choco Solver Documentation*. TASC, INRIA Rennes, LINA CNRS UMR 6241, COSLING S.A.S., 2016. URL: http://www.choco-solver.org.

26 Andreas Schutt, Thibaut Feydy, and Peter J. Stuckey. Explaining time-table-edge-finding propagation for the cumulative resource constraint. In *CPAIOR*, volume 7874 of *Lecture Notes in Computer Science*, pages 234–250. Springer, 2013.

27 Andreas Schutt and Armin Wolf. A new $O(n^2 \log n)$ not-first/not-last pruning algorithm for cumulative resource constraints. In *CP*, volume 6308 of *Lecture Notes in Computer Science*, pages 445–459. Springer, 2010.

28 Petr Vilím. *Global Constraints in Scheduling*. PhD thesis, Charles University in Prague, Faculty of Mathematics and Physics, Department of Theoretical Computer Science and Mathematical Logic, KTIML MFF, Universita Karlova, Malostranské náměstí 2/25, 118 00 Praha 1, Czech Republic, August 2007. URL: http://vilim.eu/petr/disertace.pdf.

29 Petr Vilím. Edge finding filtering algorithm for discrete cumulative resources in $O(kn \log n)$. In *CP’09: Proceedings of the 15th international conference on Principles and practice of constraint programming*, pages 802–816, Berlin, Heidelberg, 2009. Springer-Verlag. URL: http://vilim.eu/petr/cp2009.pdf.

30 Petr Vilím. Timetable edge finding filtering algorithm for discrete cumulative resources. In *CPAIOR*, volume 6697 of *Lecture Notes in Computer Science*, pages 230–245. Springer, 2011.
<figure><img src="image_182.png" /></figure>
# Exploring Hydrogen Supply/Demand Networks: Modeller and Domain Expert Views

Matthias Klapperstueck

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

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

Department of Human-Centred Computing, Faculty of IT, Monash University, Clayton, VIC, Australia

Frits de Nijs

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

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

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

Ilankaikone Senthooran

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

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

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

Jack Lee-Kopij

Woodside Energy Ltd., Perth, Australia

Maria Garcia de la Banda

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

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

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

Michael Wybrow

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

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

Department of Human-Centred Computing, Faculty of IT, Monash University, Clayton, VIC, Australia

## Abstract

Energy companies are considering producing renewable fuels such as hydrogen/ammonia. Setting up a production network means deciding where to build production plants, and how to operate them at minimum electricity and transport costs. These decisions are complicated by many factors including the difficulty in obtaining accurate current data (e.g., electricity price and transport costs) for potential supply locations, the accuracy of data predictions (e.g., for demand and costs), and the need for some decisions to be made due to external (not modelled) factors. Thus, decision-makers need access to a user-centric decision system that helps them visualise, explore, interact and compare the many possible solutions of many different scenarios. This paper describes the system we have built to support our energy partner in making such decisions, and shows the advantages of having a graphical user-focused interactive tool, and of using a high-level constraint modelling language (MINIZINC) to implement the underlying model.

**2012 ACM Subject Classification** Theory of computation → Constraint and logic programming; Theory of computation → Integer programming; Human-centered computing → Information visualization

**Keywords and phrases** Facility Location, Hydrogen Supply Chain, Human-Centric Optimisation

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

**Funding** Woodside Energy Ltd.

## 1 Introduction

Moving away from fossil fuels is needed to achieve (net) zero carbon emissions [12]. Part of this move focuses on the production of hydrogen and ammonia as storage carriers for their use in mobile applications and seasonal storage. It is estimated the world production of hydrogen needs to more than double by 2030 [9] and, given around 50% of the current production [5] relies on natural gas, the demand for “green” hydrogen is even higher. This

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

© Matthias Klapperstueck, Frits de Nijs, Ilankaikone Senthooran, Jack Lee-Kopij, Maria Garcia de la Banda, and Michael Wybrow;
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. 21; pp. 21:1–21:18

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

Leibniz International Proceedings in Informatics
LIPICS Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl Publishing, Germany
requires energy companies to expand their production capacity as efficiently as possible. To
achieve this, energy companies need to solve a complex combinatorial optimisation problem
with energy sourcing, operation, transport and demand constraints that can be summarised
as making two interdependent decisions – *where to build* the hydrogen production plants and
*how to operate* the production and supply process – that minimise costs.

This paper describes the system we have built to support our energy partner in solving this problem, which we refer to as the facility location and operation problem. While already deployed at our partner’s servers, the system is in constant evolution to incorporate the new capabilities our partner requests, while they continue using it. Importantly, the system does not just consist of a problem model that is instantiated with input data and solved using a solver, as this did not satisfy our partner’s needs. Instead, we took advantage of the capabilities of constraint programming modelling languages (MINIZINC [15]) and our knowledge of interactive user interfaces to build a system that supports users in obtaining, exploring and comparing solutions in different ways, tailored to our partner’s needs.

In particular, the system provides users with (a) a diverse set of solutions that can be
compared visually and numerically, (b) conflict resolution methods that, upon infeasibility,
identify conflicted constraints and guide users on which to modify, (c) different kinds of data
and model approximations that allow users to explore near-optimal solutions quickly, (d) the
ability to interactively add/remove/modify constraints via the interface to explore different
scenarios, (e) the ability to impose simple robustness constraints that ensure the solution is
resilient to unexpected facility outages, and (f) two different models of the same problem that
are used to cross-validate the correctness of the solutions. The amount of work required to
implement the system was significantly reduced thanks to the use of MINIZINC: its compiler
directly supports (a-b), while its high-level nature significantly simplifies (c-f).

Our industry partner has assessed our code against alternative tools and selected ours
based both on quantitative and qualitative measures (software design & architecture, user
interface, and technology stack), as well as its additional features (diversity and robustness).
This provides some measure of its quality and usefulness.

## 2 Facility location and operation problem overview

Intuitively, the facility location and operation problem solved by our Hydrogen Network Op-
timisation System (or HyNetOS) aims to select the location and configuration of the hydrogen
production plants needed to supply the demand of the given demand locations, in such a
way as to minimise the total cost of building and operating the resulting production/supply

<figure><img src="image_195.png" /><figcaption>■ **Figure 1** Left: solution showing 4 plant locations (squares) and the demand locations each of them supplies (circles in the plant's colour). Right: abstraction of a single plant.</figcaption></figure>

<figure><img src="image_196.png" /></figure>
network over a number of years. The left-hand-side of Figure 1 shows a visualisation of a particular solution, where 4 plant locations (the squares) out of all 9 plant locations given as input are selected by HyNetOS to supply all demand locations (the circles), also given as input. The decision of who supplies who is visualised via colour, with demand locations supplied by the same plant sharing the colour of that plant.

As shown in the right-hand-side of Figure 1, each production plant is itself built using one or more units of three production components – electrolysers, compressors and liquefiers – and two storage systems – for liquid and for compressed gas. The units in which each of these five types of hardware can be delivered, referred to as stock keeping units (SKUs), have different characteristics such as capacity, electricity consumption, and maximum change per hour in production rates, which are given as input. The hydrogen is transported between supply and demand locations using a truck-based transportation network, where the transport costs between each supply-demand location are given via an input table. Note that during the plant sizing, we take into account the impact boil-off has on the final quantity delivered after transport. Boil-off is the decompression due to the boiling of residual hydrogen liquid as a result of increased temperature inside a (nearly) empty truck's storage tank. Further, we also take into account the location specific price of the available electricity sources necessary to operate each facility, including the generation of a facility operation schedule to minimise electricity costs under variable electricity price and availability.

Importantly, the problem is defined across two timescales: *periods* and *hours*. A period is the number of consecutive years during which the demand for hydrogen is assumed to stay constant. Plants are built during some period and must only increase in size in later periods since, currently, demand is assumed never to decrease. In contrast, the electricity prices for some markets are given hourly for an entire year, yielding up to 8760 possible different prices per market. Because of this, we can potentially perform market arbitrage by adjusting production rates to fit the market price for electricity. However, this means the optimisation needs to make power consumption decisions on an hourly resolution.

The problem's objective is to minimise the total cost of the network, which can be broken down into the following interdependent cost elements: *building* (CAPEX) and *operating* (OPEX) the hydrogen production plants; *transporting* the product from supply to demand locations; and *powering* the production plants using the available electricity sources. The latter includes a monthly cost (*demand charge*) some electricity providers add to try to flatten their consumer's load profile, and is proportional to the plant's highest monthly load.

From an optimisation perspective, this problem is challenging due to the large search space created by the high number of demand locations to be supplied (several hundreds), the different electricity sources often available, and the hourly electricity prices the system must consider across an entire year (up to 8760 per source). This is further complicated by decisions on two different timescales: plant construction decisions are taken for each period, while operating decisions are taken hourly. In addition, the dramatic scale difference between CAPEX (millions) and electricity costs (cents per kWh), makes the model numerically unstable. Together, these factors mean that solving the hydrogen facility model to optimality can quickly become out of reach even for commercial mixed-integer solvers such as GUROBI [7]. Furthermore, due to the costs involved, plant construction decisions necessitate the examination of a range of scenarios, requiring the problem to be solved many times.
### 3 Modeller's view

As HyNetOS is developed for industry and there is no ground-truth, we put special attention on reducing modelling errors. To do this, we separately implemented and integrated two models of the problem whose functional equivalence is continuously checked by ensuring solutions to any model instance can be given to the other without drop in objective or infeasibility. This practice increases the redundancy of the programming task, thereby reducing the residual probability of errors [21]. This section presents the input data used by both models (referred to as Ori for the original model, and Alt for the alternative one, both always compiled into a MILP problem by MINIZINC), the constraints implemented by one (Alt) due to space constraints, and the modelling changes that most improved solving time.

#### 3.1 Input data

Both models require the following input data, which MINIZINC refers to as *parameters*:

- $\mathcal{D}$ and $\mathcal{S}$: set of client demand locations and set of candidate supply locations, respectively.

- $\mathcal{K}$: set of products; currently $\mathcal{K} = \{\text{LIQ, GAS}\}$, i.e., liquid and gas hydrogen, respectively.

- $\mathcal{H}$: set of hardware; currently $\mathcal{H} = \{\text{ELEC, COMP, LIQF, GAS\_STR, LIQ\_STR}\}$ corresponding to electrolyser, compressor, liquifier, gas and liquid storage types, respectively.

- $\mathcal{P}$: set of constant demand periods and, for each period $p \in \mathcal{P}$, how many years $p_y \in \mathbb{N}$ it covers within the given plan horizon (typically 20 years).

- $\mathcal{SO}_i$: set of electricity sources (market labels) available at supply location $s_i \in \mathcal{S}$. Elements of set $\mathcal{SO}_i$ correspond to year-long time series of electricity prices, recorded at an hourly resolution. Markets are one of three types: (a) UTILITY, a conventional metered connection with mostly fixed prices, (b) PPA, a fixed-price power purchase agreement with a renewable energy provider (e.g., solar or wind farm), and (c) WHOLESALE, a market with generally the lowest price but exposed to volatile price fluctuations. We distinguish them for their unique features; demand charges typically only apply to UTILITY, while PPA is a zero carbon emissions source, necessary to achieve carbon targets in the future.

- $\mathcal{T}_i \subseteq \{t_1, t_2, \dots, t_\tau\}$: set of time steps considered for a single representative year at supply location $s_i \in \mathcal{S}$. Each time step $t \in \mathcal{T}_i$ has a duration in hours $h_{i,t} \in \{1, \dots, 24\}$, during which the price of electricity is constant. The value of $\tau$ ranges between 365 and 8760, corresponding to the cases where every time step represents 1 day and 1 hour, respectively. Their sum $\sum_{t \in \mathcal{T}_i} h_{i,t}$ must always be 8760.

- $D_{p,j}^k \in \mathbb{R}_{\ge 0}$: daily demand of product $k \in \mathcal{K}$ from location $d_j \in \mathcal{D}$ throughout period $p \in \mathcal{P}$, in tonnes per day.

- $T_{p,i,j}^k \in \mathbb{R}_{> 0}$: transport cost in $ per kg of product $k$ sent from supply location $s_i \in \mathcal{S}$ to demand location $d_j \in \mathcal{D}$ during period $p \in \mathcal{P}$.

- $c_{p,i,t}^{\text{so}} \in \mathbb{R}_{> 0}$: electricity cost of source so $\in \mathcal{SO}_i$ at supply location $s_i \in \mathcal{S}$ during time step $t \in p_y$ of any year of period $p \in \mathcal{P}$, in millions of $ per MW.

- $\hat{c}_{p,i,t}^{\text{so}} \in \mathbb{R}_{\ge 0}$: demand charge of source so $\in \mathcal{SO}_i$ at supply location $s_i \in \mathcal{S}$ during time step $t \in p_y$ of any year of period $p \in \mathcal{P}$ in millions of $ per MW. Common for UTILITY.

- $\text{LB}_i^{\text{so}} \in \mathbb{R}_{\ge 0}$: minimum annual energy usage in MWh required to be allowed to consume energy from source so $\in \mathcal{SO}_i$ at location $s_i \in \mathcal{S}$. In practice, only relevant for PPA.

For each hardware type $h \in \mathcal{H}$, we also need the set of concrete stock keeping units (SKUs) in which $h$ can be delivered, and the following hardware-specific properties:

- $N^h$: number of available SKUs for $h$, from which we build index set $\mathcal{C} = \{1..N^c\}$.

- $C^c \in \mathbb{R}_{> 0}$: capacity of SKU $c \in \mathcal{C}$ in tonnes per day.

- $K_p^c$: ownership cost in millions of $ for one SKU $c \in \mathcal{C}^h$ starting from period $p \in \mathcal{P}$.
- $\nu^h \in [0..1]:$ minimum production (turndown) rate proportional to installed capacity for $h$.

- $e^h \in \mathbb{R}_{\ge 0}$: electricity usage in MWh per tonne of production needed to run $h$.

- $\rho^h \in [0..1]:$ overhead to production and storage for $h \in \{\text{COMP, LIQF, GAS\_STR, LIQ\_STR}\}$ required to compensate for boil-off.

- $\mu \in [0..1]:$ maximum ramping rate of the liquefaction units.

## 3.2 Decision variables, constraints and objective function

Figure 2 shows the objective and constraints in model Alt, where blue is used for decision variables and black for input data, split into the following six groups.

**Network constraints.** As hydrogen demand is assumed to stay constant throughout any period $p \in \mathcal{P}$, we assume the supply network will also stay constant throughout $p$. Thus, we only have to decide once per period $p$ which location $s_i \in \mathcal{S}$ supplies (part of) the daily demand of location $d_j \in \mathcal{D}$ for product $k \in \mathcal{K}$. The network constraints encode these assumptions. Variable $x_{p,i,j}^k$ represents the percentage of the demand $d_j$ of product $k$ supplied by location $s_i$ in period $p$. Constraint (1) ensures that every product of every demand location is serviced in its entirety by the supply locations. In order to satisfy the demand induced by the network, production plants must be built at one or more of the supply locations. Binary variable $b_{p,i}^k$ represents whether a plant for product $k$ is built at location $s_i$ at (or before) period $p$, i.e., whether there is a plant in $s_i$ producing $k$ in period $p$. Constraint (2) ensures $b_{p,i}^k = 1$ whenever any $x_{p,i}^k > 0$, thereby capturing the need for a plant. Currently, hybrid facilities (producing both gas and liquid for offtake) are disallowed via constraint (3); a requirement of our industry partner.

The objective, shown in expression (26), includes as its first component the sum of the total cost of transporting the product across the links selected with non-zero weight $x$, where transporting the entire demand (of $D_{p,j}^k$ tonnes per day) of product $k$ from $s_i \to d_j$ costs $T_{p,i,j}^k$ in period $p$.

**Offtake constraints.** During a period $p \in \mathcal{P}$, the daily demand $D_{p,j}^k$ of product $k \in \mathcal{K}$ at every location $d_j \in \mathcal{D}$ is constant. Thus, we denote as $z_{p,i}^k$ the offtake rate (in tonnes per day) of any supply location $s_i \in \mathcal{S}$ for period $p$ and product $k$, i.e., the rate at which the plant at location $s_i$ must produce $k$ over the long term, using storage to buffer production rate changes in the short term. Constraints (4) and (5) define $z_{p,i}^k$ for $k = \text{LIQ}$ and $k = \text{GAS}$, respectively. Note that the offtake of gas $z_{p,i}^{\text{GAS}}$ includes the material needed to produce liquid *from* that gas through liquefaction, and not just the transport requirements of gas.

**Plant capacity constraints.** Producing $z_{p,i}^k$ of product $k \in \mathcal{K}$ requires a production plant at location $s_i \in \mathcal{S}$ during period $p \in \mathcal{P}$ that is of suitable capacity. This can be modelled using the following five constraints. First, let variable $u_{p,i}^c$ represent how many of SKU $c \in \mathcal{C}^h$ of hardware type $h \in \mathcal{H}$ are newly installed at the start of period $p$ (and thus available for use in $p$) in supply location $s_i$. The maximum (or peak) production capacity of the plant at location $s_i$ in period $p$ is defined by constraint (6) as the sum of the SKUs capacities installed *up to* $p$, and denoted by variable $y_{p,i}^h$. By tying the capacity to the cumulative number of newly installed units, we ensure plants can only grow in size.

Second, the plant must be big enough for its peak capacity to meet the offtake rate $z_{p,i}^k$. This is ensured by constraint (7), which adds the overhead factor $\rho^h$ to the normal offtake to take transportation boil-off into account. The boil-off factor is applied differently to different
<table><caption>Figure 2 Objective and constraints for the Alt model of the HyNetOS decision system.</caption><tbody><tr><td rowspan="3">Network</td><td>Σ<sub>s<sub>i</sub>∈S</sub> x<sup>k</sup><sub>p,i,j</sub> = 1</td><td>∀p, d<sub>j</sub>, k where D<sup>k</sup><sub>p,j</sub> > 0</td><td>(1)</td></tr><tr><td>x<sup>k</sup><sub>p,i,j</sub> ≤ b<sup>k</sup><sub>p,i</sub></td><td>∀p, s<sub>i</sub>, d<sub>j</sub>, k</td><td>(2)</td></tr><tr><td>Σ<sub>k∈K</sub> b<sup>k</sup><sub>p,i</sub> ≤ 1</td><td>∀p, s<sub>i</sub></td><td>(3)</td></tr><tr><td rowspan="2">Offtake</td><td>Σ<sub>d<sub>j</sub>∈D</sub> D<sup>LIQ</sup><sub>p,j</sub>x<sup>LIQ</sup><sub>p,i,j</sub> = z<sup>LIQ</sup><sub>p,i</sub></td><td>∀p, s<sub>i</sub></td><td>(4)</td></tr><tr><td>z<sup>LIQ</sup><sub>p,i</sub> + Σ<sub>d<sub>j</sub>∈D</sub> D<sup>GAS</sup><sub>p,j</sub>x<sup>GAS</sup><sub>p,i,j</sub> = z<sup>GAS</sup><sub>p,i</sub></td><td>∀p, s<sub>i</sub></td><td>(5)</td></tr><tr><td rowspan="4">Plant capacity</td><td>Σ<sub>1≤p'<sub>i</sub>≤p, c∈C<sup>h</sup></sub> C<sup>c</sup>u<sup>c</sup><sub>p'<sub>i</sub></sub> = y<sup>h</sup><sub>p,i</sub></td><td>∀p, s<sub>i</sub>, h</td><td>(6)</td></tr><tr><td>(1 + ρ<sup>h</sup>)z<sup>h→k</sup><sub>p,i</sub> ≤ y<sup>h</sup><sub>p,i</sub></td><td>∀p, s<sub>i</sub>, h ∈ {ELEC, COMP, LIQF}</td><td>(7)</td></tr><tr><td>z<sup>h→k</sup><sub>p,i</sub> ≥ ν<sup>h</sup>y<sup>h</sup><sub>p,i</sub></td><td>∀p, s<sub>i</sub>, h ∈ {ELEC, COMP, LIQF}</td><td>(8)</td></tr><tr><td>ρ<sup>h→k</sup>z<sup>h→k</sup><sub>p,i</sub> + m<sup>h→k</sup>b<sup>h→k</sup><sub>p,i</sub> ≤ y<sup>h</sup><sub>p,i</sub></td><td>∀p, s<sub>i</sub>, h ∈ {GAS_STR, LIQ_STR}</td><td>(9)</td></tr><tr><td rowspan="4">Production</td><td>m<sup>LIQ</sup>z<sup>LIQ</sup><sub>p,i</sub> ≤ y<sup>LIQ_STR</sup><sub>p,i</sub></td><td>∀p, s<sub>i</sub></td><td>(10)</td></tr><tr><td>ν<sup>h</sup>y<sup>h→k</sup><sub>p,i</sub> ≤ ˆp<sup>h→k</sup><sub>p,i</sub></td><td>∀p, s<sub>i</sub>, h ∈ {ELEC, COMP, LIQF}</td><td>(11)</td></tr><tr><td>ˆp<sup>h→k</sup><sub>p,i</sub> + p<sup>h→k</sup><sub>p,i,t</sub> + ρ<sup>h</sup>z<sup>h→k</sup><sub>p,i</sub> ≤ y<sup>h</sup><sub>p,i</sub></td><td>∀p, s<sub>i</sub>, t, h ∈ {ELEC, COMP, LIQF}</td><td>(12)</td></tr><tr><td>p<sup>LIQ</sup><sub>p,i,t+1</sub> - p<sup>LIQ</sup><sub>p,i,t</sub> ≤ μh<sub>i,t</sub>y<sup>LIQF</sup><sub>p,i</sub></td><td>∀p, s<sub>i</sub>, t</td><td>(13)</td></tr><tr><td rowspan="3">Storage</td><td>p<sup>LIQ</sup><sub>p,i,t</sub> - p<sup>LIQ</sup><sub>p,i,t+1</sub> ≤ μh<sub>i,t</sub>y<sup>LIQF</sup><sub>p,i</sub></td><td>∀p, s<sub>i</sub>, t</td><td>(14)</td></tr><tr><td>p<sup>LIQ</sup><sub>p,i,0</sub> - p<sup>LIQ</sup><sub>p,i,|T<sub>i</sub>|</sub> ≤ μh<sub>i,t</sub>y<sup>LIQF</sup><sub>p,i</sub></td><td>∀p, s<sub>i</sub>, t</td><td>(15)</td></tr><tr><td>p<sup>LIQ</sup><sub>p,i,|T<sub>i</sub>|</sub> - p<sup>LIQ</sup><sub>p,i,0</sub> ≤ μh<sub>i,t</sub>y<sup>LIQF</sup><sub>p,i</sub></td><td>∀p, s<sub>i</sub>, t</td><td>(16)</td></tr><tr><td rowspan="4">Electricity usage</td><td>s<sup>h→k</sup><sub>p,i,t</sub> ≤ y<sup>h</sup><sub>p,i</sub></td><td>∀p, s<sub>i</sub>, t, h ∈ {GAS_STR, LIQ_STR}</td><td>(17)</td></tr><tr><td>s<sup>GAS</sup><sub>p,i,t+1</sub> = s<sup>GAS</sup><sub>p,i,t</sub> + h<sub>i,t</sub>(ˆp<sup>GAS</sup><sub>p,i</sub> + p<sup>GAS</sup><sub>p,i,t</sub> - ˆp<sup>LIQ</sup><sub>p,i</sub> - p<sup>LIQ</sup><sub>p,i,t</sub> + p<sup>LIQ</sup><sub>p,i,t+1</sub> - z<sup>GAS</sup><sub>p,i</sub> + z<sup>LIQ</sup><sub>p,i</sub>)</td><td>∀p, i, t</td><td>(18)</td></tr><tr><td>s<sup>LIQ</sup><sub>p,i,t+1</sub> = s<sup>LIQ</sup><sub>p,i,t</sub> + h<sub>i,t</sub>(ˆp<sup>LIQ</sup><sub>p,i</sub> + p<sup>LIQ</sup><sub>p,i,t</sub> + p<sup>LIQ</sup><sub>p,i,t+1</sub> - z<sup>LIQ</sup><sub>p,i</sub>)</td><td>∀p, i, t</td><td>(19)</td></tr><tr><td>s<sup>k</sup><sub>p,i,1</sub> = s<sup>k</sup><sub>p,i,|T<sub>i</sub>|</sub> = 0</td><td>∀p, s<sub>i</sub>, k</td><td>(20)</td></tr><tr><td rowspan="4">Obj.</td><td>24e<sub>p,i,t</sub> = e<sup>ELEC</sup>(ˆp<sup>GAS</sup><sub>p,i</sub> + p<sup>GAS</sup><sub>p,i,t</sub>) + e<sup>COMP</sup>(ˆp<sup>GAS</sup><sub>p,i</sub> + p<sup>GAS</sup><sub>p,i,t</sub> + ρ<sup>GAS</sup>z<sup>GAS</sup><sub>p,i</sub>) + e<sup>LIQF</sup>(ˆp<sup>LIQ</sup><sub>p,i</sub> + p<sup>LIQ</sup><sub>p,i,t</sub> + p<sup>LIQ</sup><sub>p,i,t+1</sub> + ρ<sup>LIQ</sup>z<sup>LIQ</sup><sub>p,i</sub>)</td><td>∀p, i, t</td><td>(21)</td></tr><tr><td>Σ<sub>so∈S</sub> e<sup>SO</sup><sub>p,i,t</sub> = e<sup>SO</sup><sub>p,i,t</sub></td><td>∀p, s<sub>i</sub>, t</td><td>(22)</td></tr><tr><td>ˆc<sup>SO</sup><sub>p,i,t</sub>e<sup>SO</sup><sub>p,i,t</sub> ≤ ˆe<sup>SO</sup><sub>p,i,m[t]</sub></td><td>∀p, s<sub>i</sub>, t, so</td><td>(23)</td></tr><tr><td>(Σ<sub>t∈T<sub>i</sub></sub>(h<sub>i,t</sub>e<sup>SO</sup><sub>p,i,t</sub>) > 0) ⇒ (Σ<sub>t∈T<sub>i</sub></sub>(h<sub>i,t</sub>e<sup>SO</sup><sub>p,i,t</sub>) ≥ LB<sup>SO</sup><sub>i</sub>)</td><td>∀p, s<sub>i</sub>, so</td><td>(24)</td></tr><tr><td></td><td>u<sup>c</sup><sub>p,i</sub> ∈ N<sub>0</sub>; b<sup>k</sup><sub>p,i</sub> ∈ {0, 1}; 0 ≤ x<sup>k</sup><sub>p,i,j</sub> ≤ 1; y, z, p, s, e ∈ ℝ≥0</td><td>∀p, s<sub>i</sub>, t, so</td><td>(25)</td></tr><tr><td></td><td>min Σ<sub>p,i,j,k</sub> (T<sup>k</sup><sub>p,i,j</sub>x<sup>k</sup><sub>p,i,j</sub>) + Σ<sub>p,i,c</sub> (K<sup>c</sup><sub>p</sub>u<sup>c</sup><sub>p,i</sub>) + Σ<sub>p,i,t,so</sub> (c<sup>SO</sup><sub>p,i,t</sub>e<sup>SO</sup><sub>p,i,t</sub>) + Σ<sub>p,i,m,so</sub> (ˆe<sup>SO</sup><sub>p,i,m</sub>)</td><td>∀p, s<sub>i</sub>, t, so</td><td>(26)</td></tr></tbody></table>
types of hardware $h \in \mathcal{H}$ because boil-off comes in uncompressed gas form. Thus, it has to pass through the compressor and the liquefactor but not the electrolyser. We use $h \rightarrow k$ to indicate a mapping from hardware to product. For example, since the electrolyser $h = \text{ELEC}$ produces $k = \text{GAS}$, the electrolyser must be scaled to meet gas demand $z_{p,i}^{h \rightarrow k}$. Third, the plant should not be so big that running at minimum capacity yields more product on average than is demanded. Constraint (8) captures this requirement, where the minimum production rate depends on the minimum turndown capabilities $\nu^h$ of the hardware. Finally, we impose two kinds of minimum storage size constraints. The first constraint ensures two things: 1) that every plant has at least some gas storage for buffering boil-off during transport, which is captured as a fraction $\rho^{h \rightarrow k}$ of daily production rate $z_{p,i}^k$, and 2) that every plant meets a global minimum storage amount $m^{h \rightarrow k}$ if it is built ($z_{p,i}^{h \rightarrow k} > 0$) (constraint (9)).

The second constraint ensures enough liquid storage is built to sustain liquid hydrogen offtake for a minimum number of days $m^{\text{LIQ}}$ (constraint (10)).

The construction of any SKU $c \in \mathcal{C}^h$ of hardware $h \in \mathcal{H}$ incurs CAPEX and OPEX costs. These are aggregated into the period-specific component cost $K_p^c$ (in millions of $ per unit) and accumulated via the number of SKUs newly installed in $p$ (given by $u_{p,i}^c$), yielding the objective term $K_p^c u_{p,i}^c$. The sum of these objective terms forms the second component of the objective, shown in expression (26).

**Production constraints.** For the plant operation, we model the average operating costs through a representative year per period. This means we create a production schedule for one year in each (constant demand) period, which then repeats for however many years the period is long. The production rate of hardware $h \in \mathcal{H}$ for product $k \in \mathcal{K}$ at supply location $s_i \in \mathcal{S}$ during period $p \in \mathcal{P}$ is allowed to change during each of the time steps $t \in \mathcal{T}_i$ defined for that location. The production rate is measured in tonnes per day and modelled via two terms: $\hat{p}_{p,i}^{h \rightarrow k} + p_{p,i,t}^{h \rightarrow k}$. Here, $\hat{p}^{h \rightarrow k}$ is the constant baseline production amount, and $p_{p,i,t}^{h \rightarrow k}$ the component that is variable in time step $t$. Constraint (11) ensures the baseline exceeds the minimum turndown production rate. Constraint (12) ensures the flexible component never exceeds the installed capacity, while also keeping sufficient headroom for the boil-off fraction $\rho$ that has to be re-compressed and re-liquefied. Constraints (13) and (14) deal with $h = \text{LIQF}$ having a slow ramping speed. They constrain the ramp up and down rate, respectively, between two consecutive production rates (where the consecutive steps wrap around the year, via constraints (15) and (16)), such that the change in rate does not exceed the ramping capability $\mu$ of liquefaction (0.1 of total capacity, i.e., 10%) by the number of hours $h_{i,t}$ for which the ramping is maintained.

**Storage constraints.** Variable production is buffered via storage: overproduction causes the storage level to increase, while underproduction causes the constant offtake to drain the storage. At every time step $t \in \mathcal{T}_i$ of supply location $s_i \in \mathcal{S}$ in period $p \in \mathcal{P}$, we track the storage level of product $k \in \mathcal{K}$ in tonnes via variable $s_{p,i,t}^k$. Constraint (17) ensures we never store more than the plants' installed storage capacity. Changes in storage levels are captured by constraints (18) for $k = \text{GAS}$, and (19) for $k = \text{LIQ}$. In both, the storage in $t+1$ is an offset from the storage in $t$, plus hourly production (baseline plus flexible), minus offtake (either by the liquefaction unit, or by the demand locations constant offtake factor). In addition, we anchor the storage to a baseline with a wrap-around constraint in (20). Adding a baseline of 0 at each year's end removes some of the symmetries in storage schedule solutions and ensures the storage level is implementable across period changes.
**Electricity cost constraints.** The last set of constraints connects the daily hydrogen production to the electricity cost required to produce it. To this end, at every time step $t \in T_i$ for the plant of supply location $s_i \in S$ in period $p \in P$, we represent the power consumption in MW of the plant via variable $e_{p,i,t}$. This variable is defined by constraint (21) as the sum of the power consumed by the plant's hardware $h \in \{ELEC, COMP, LIQF\}$ to produce the required hydrogen in $t$ of $p$, i.e., as the electricity usage $e^h$ multiplied by the production rate at that time which, as before, is formed by a baseline $\hat{p}_{p,i,t}^k$ component, plus a flexible and possibly a boiloff one. Since our production rates are daily, we downscale the power consumption by factor 24 to the per-hour rate. Constraint (22) distributes the power consumption among sources via variables $e_{p,i,t}^{so}$, which represent the power consumed from each electricity source $so \in SO$ available at that location. These variables are used to derive the cost of this power, which is accumulated in the third sum of the objective. Constraint (23) defines variable $\hat{e}_{p,i,m[t]}^{so}$ representing the monthly demand charge for the month $m[t] \in \{1, \dots, 12\}$ where time step $t$ falls. It is the maximum demand charge (cost $c_{p,i,t}^{so}$ of each MW consumed $e_{p,i,t}^{so}$) that can be obtained due to the power consumed at each $t$ of that month. This accounts for the last sum of the objective. Finally, constraint (24) ensures electricity source $so \in SO$ is selected only if its minimum annual energy usage $LB_i^{so}$ is met.

## 3.3 Effective modelling and solving in practice

One of the primary benefits of using a *Constraint Modelling Language* such as MINIZINC is that it takes care of the details of mapping models to solvers efficiently (e.g., linearisation). This is important because writing good (i.e. efficiently solvable) models can otherwise be a process of significant trial-and-error, as much as art form as a science. While there are important strategies for good modelling [22], the “last mile” of good practice involves rules of thumb, such as minimizing the number of integer variables in favour of floats, or avoiding equalities if possible. MINIZINC takes care of these during compilation, allowing the modeller to focus on semantic changes. Nevertheless, because of its high level nature, it is equally well-suited to rapid prototyping of model changes. This is why for our industry partner’s quantitative evaluation of the system, we were able to use MINIZINC to try to find model improvements that reduce the total runtime. To do this we took a reference implementation of the model and evaluated the runtime of several model changes. The progression of these changes is shown as a box plot in Figure 3 and discussed below.

**1. Split power equality.** The power equality (21) involves variables $e_{p,i,t}^{so}$ that are directly minimized in the objective. As such, it is tempting to remove the equality constraint by replacing it with $\ge$, such that the model will always demand at least as much power as needed to produce the hydrogen. However, this change allows overconsumption of power to meet the PPA minimum requirement (24). Therefore, we split the power equality constraint into a greater than and a *conditionally applied* less than part, which is only applied if the source has a minimum usage requirement through a conditional constraint:

forall (i) if (MIN_USE[i]) then $\sum_{so \in SO} e_{p,i,t}^{so} = e_{p,i,t}$ else $\sum_{so \in SO} e_{p,i,t}^{so} \ge e_{p,i,t}$ endif;

**2. Tighter constraints.** Plant production boolean $b_{p,i}^k$ is tied to many floating point supply indicator variables. Instead of defining them separately, we can define the state activity indicator with a single constraint, which will be active if any demand is supplied, i.e.:

$$ \sum_{d_j \in D} x_{p,i,j}^k \le |D| b_{p,i}^k, \quad \forall p, s_i, k $$
<figure><img src="image_197.png" /><figcaption>Figure 3 Response of model solve time as a result of model changes; mean time annotated.</figcaption></figure>

**3. Direct objective formulation.** In our reference implementation, the objective terms are captured by intermediate variables for each of the sums that make up the objective. We can rewrite the objective directly in terms of the base variables, allowing the compiler to group terms together, and furthermore helping the MIP solver to prove optimality faster.

We evaluated the time-to-optimality with each model change over 10 runs with different seeds for the solver, to average out effects of solver-internal randomized branching choices and heuristics, and capture statistically meaningful averages on the total runtime. We used a Linux machine with AMD Ryzen 9 3950X 16-Core Processor (3.5-4.7 GHz) and the input data file A shown in Table 1, which has 9 supply locations, 100 demand locations, 3 SKUS per hardware type, and an average period length of 642 tiers (min 365, max 1338). Figure 3 tracks the distributions of total runtime as the changes are implemented, for the GUROBI 9.5.1 and HIGHS 1.5.0 [8] solvers. We also evaluated the combination of 2 and 3 which individually seemed to produce improvements. Nevertheless, the combination of compilation and pre-solving means that most changes have no significant effect. The biggest improvement comes from adjusting how GUROBI scales the objective terms, by changing the SCALEFLAG parameter. We believe this is due to the nature of the hydrogen facility location problem, where many candidate integer assignments are inherently feasible: we can always redistribute the supply network and adjust production rates. The optimisation is thus primarily guided by the (large) objective. We hypothesise that a core-guided search [6], which always assumes constraints can be satisfied and iteratively tightens the model upon detecting an infeasibility, would be a better strategy to solve this kind of model. However, due to the prevalence of floating point data, no core-guided solver available to us can be applied to this problem.

**4 From modellers to users**

To integrate our models into a useful decision system we extended them with functionality that its users (engineers and business experts) could access via the system’s user interface (see Section 5). This includes minor changes to allow many decisions to be set to true/false by the user to explore different scenarios, as well as the more significant ones discussed here.

**4.1 Run-time versus accuracy – approximations and warm-starts**

The large number of variables involved in the hourly operational decisions, can make solving our models to optimality too time consuming for the user. Further, approximate solutions may be sufficient whenever users are interactively exploring what-if scenarios. Our system gives users several ways to control the time versus optimality trade-off, by means of the following two approximation options and warm-start strategy.
**Constraint approximation – annual constant production.** This approximation (referred to as App) allows each plant to run at a constant daily production rate sufficient to exactly meet the daily demand that the plant supplies. This means no intra-day or seasonal storage is required and plants can simply build the minimum storage capacity. As a result, in this version of the model all the production and storage constraints (11–20) are removed, and the electricity usage constraints are reformulated to operate on the total annual energy demand under a constant production. As such, the (hourly) energy cost coefficients $c_{p,i,t}^{\text{so}}$ are rescaled to the total annual cost for a given offtake rate in tonnes per day:

$$c_{p,i}^{\text{so}} = (e^{\text{E}} + e^{\text{C}} + e^{\text{L}}) \sum_{t \in \mathcal{T}_i} \left( c_{p,i,t}^{\text{so}} \cdot \frac{h_{i,t}}{24} \right)$$

**Data approximation – price-tier grouping.** In practice, many hours in a year often have identical prices for electricity. Since the only driver for operational decisions is a reduction in electricity price, we can group consecutive time steps with equal price. We use the utility market to inform the grouping operation, as utility prices often stay constant for hours, e.g. a daytime tariff and a cheaper night-time tariff from 11pm to 6am. For the other markets in the same compressed time step, we take the median price as the constant price. This approximation (referred to as Tie) reduces the cardinality $\tau$ of the set $\mathcal{T}_i$, thereby significantly reducing the total number of floating point variables in the instance (see Table 1).

**Warm-starting strategy.** Both App and Tie approximations are valid, i.e., if their solutions are fed to the non-approximated model (referred to as Exa and Hou, respectively) they also yield a solution: we can always set the production schedule at a finer time granularity to be equal to the constant-production (rate) assumptions at the coarser level. In addition, they can be combined to further speed up the search at the cost of accuracy. This observation informs our three-stage warm-starting strategy:

1. Solve using an approximate version of the model (i.e., App-Hou or App-Tie) to optimality;

2. Assign the integer decisions of that solution to the equivalent variables using an exact version of the model (i.e., Exa-Hou or Exa-Tie, resp) and resolve;

3. Warm-start that last used model with the solution from step 2 and resolve to optimality. The final step arrives at the same objective as a traditional execution (referred to as cold-start) using the same model; however, it will find a good quality initial solution much faster.

<table><caption>Table 1 Size, solving time and objective value of the instances obtained with four data files.</caption><thead><tr><th colspan="5">Input data size</th><th colspan="5">Instance size for GUROBI solver</th></tr><tr><th>File</th><th>S</th><th>D</th><th>N</th><th>M</th><th>t/day</th><th>Model</th><th>Float</th><th>Int</th><th>Cons.</th><th>Time(s)</th><th>Ratio</th></tr></thead><tbody><tr><td rowspan="4">A</td><td rowspan="4">9</td><td rowspan="4">100</td><td rowspan="4">9</td><td rowspan="4">18</td><td rowspan="4">[40]</td><td>Exa-Hou</td><td>597.7 k</td><td>144</td><td>835.6 k</td><td>5114.85</td><td>1.000</td></tr><tr><td>Exa-Tie</td><td>96.6 k</td><td>144</td><td>65.3 k</td><td>34.96</td><td>1.020</td></tr><tr><td>App-Hou</td><td>2.0 k</td><td>144</td><td>3.4 k</td><td>0.22</td><td>1.055</td></tr><tr><td>App-Tie</td><td>2.0 k</td><td>144</td><td>3.4 k</td><td>0.19</td><td>1.068</td></tr><tr><td rowspan="4">B</td><td rowspan="4">7</td><td rowspan="4">100</td><td rowspan="4">9</td><td rowspan="4">14</td><td rowspan="4">[278]</td><td>Exa-Hou</td><td>457.1 k</td><td>112</td><td>642.2 k</td><td>20982.06</td><td>1.000</td></tr><tr><td>Exa-Tie</td><td>74.7 k</td><td>112</td><td>56.5 k</td><td>38.06</td><td>1.008</td></tr><tr><td>App-Hou</td><td>1.6 k</td><td>112</td><td>2.7 k</td><td>0.19</td><td>1.017</td></tr><tr><td>App-Tie</td><td>1.6 k</td><td>112</td><td>2.7 k</td><td>0.17</td><td>1.016</td></tr><tr><td rowspan="4">C</td><td rowspan="4">9</td><td rowspan="4">100</td><td rowspan="4">9</td><td rowspan="4">18</td><td rowspan="4">[278]</td><td>Exa-Hou</td><td>597.7 k</td><td>144</td><td>835.6 k</td><td>22605.64</td><td>1.000</td></tr><tr><td>Exa-Tie</td><td>96.6 k</td><td>144</td><td>65.3 k</td><td>84.38</td><td>1.009</td></tr><tr><td>App-Hou</td><td>2.0 k</td><td>144</td><td>3.4 k</td><td>0.31</td><td>1.017</td></tr><tr><td>App-Tie</td><td>2.0 k</td><td>144</td><td>3.4 k</td><td>0.31</td><td>1.018</td></tr><tr><td rowspan="4">D</td><td rowspan="4">9</td><td rowspan="4">709</td><td rowspan="4">9</td><td rowspan="4">18</td><td rowspan="4">[283, 293]</td><td>Exa-Hou</td><td>1217.4 k</td><td>360</td><td>1706.7 k</td><td>—</td><td>—</td></tr><tr><td>Exa-Tie</td><td>215.2 k</td><td>360</td><td>166.0 k</td><td>1309.47</td><td>&gt; 1.016</td></tr><tr><td>App-Hou</td><td>26.1 k</td><td>360</td><td>42.2 k</td><td>7.81</td><td>&gt; 1.027</td></tr><tr><td>App-Tie</td><td>26.1 k</td><td>360</td><td>42.2 k</td><td>5.37</td><td>&gt; 1.027</td></tr></tbody></table>
<figure><img src="image_198.png" /><figcaption>Figure 4 Quality of incumbent solution (ratio over optimum) as a function of log runtime obtained by the warm- and cold-start strategies on each instance. Stage solutions are annotated with points.</figcaption></figure>

**Experimental snapshot.** Table 1 shows a snapshot of the trade-offs obtained by the above methods instantiated with four data files (A-D). For each file, it shows the number of supply locations $|S|$; demand locations $|D|$; SKUs over all hardware types, $N = \sum_h N^h$; markets $M = \sum_i |SO_i|$, and total demand per period in tonnes/day. It then shows the number (in thousands) of floating-point variables, integer variables, and constraints in the instance resulting from compiling each data file with either an Exa-Hou, Exa-Tie, App-Hou or App-Tie model for the GUROBI solver. The last two columns show the time (where – indicates a timeout) to find an optimal solution and prove optimality when executing each instance in cold-start mode using MINIZINC 2.7.6 and GUROBI 10.0.0.2, and the ratio between the best objective value found using the most accurate Exa-Hou model and the others.

Figure 4 compares the improvement of feasible solutions over time obtained with a cold-start run of the Exa-Hou model instantiated with the data files of Table 1, and with our warm-start strategy, using MINIZINC 2.7.6 and GUROBI 10.0.0.2. Typically, by the time the cold-start finds its first feasible solution (15+ minutes), the warm-start is in its third stage and up to 5% gap to optimality.

## 4.2 Diversity of solutions

As is common in optimisation, every instance of our problem has many optimal (and even more close-to-optimal) solutions. It is thus useful for decision-makers to obtain a number of close-to-optimal alternative solutions that are meaningfully different in terms of the components being optimised (e.g., transport vs electricity cost) or some major decisions (e.g., plants locations or capacities). Ingmar et al. [11] proposed strategies to find *N* diverse solutions of a model instance where the model includes the required user-defined diversity measures. One of these strategies iterates at most *N* – 1 times from an optimal solution looking for a new solution that is as close to optimality as desired by the user, and is maximally diverse w.r.t. all previously found solutions.

This strategy is now implemented<fnref n="1" /> in MINIZINC and used in our system by simply importing from our models a set of diversity measures our users can choose from. The measure most used currently is the Manhattan distance between the vectors of the Boolean variables $b_{p,i}^k$ indicating there is a plant at supply location $i \in S$ in period $p \in P$. Figures 5 and 7(D) show the results of one of the implemented diversity metrics, i.e. maximum diverse set of supply locations. To highlight the differences between diverse solutions our

<footnote n="1">¹ It will be released as part of the MiniZinc Python package [3] before the end of the year.</footnote>
system allows users to compare solutions in terms of their supply/demand network, the configuration of the plants they selected, their cost components, etc. Such comparisons have already been useful in identifying, for example, supply locations that appear in all optimal and close-to-optimal solutions due to their favourable cost coefficients (e.g. green location in Figure 5) and are thus prime candidates for construction.

<figure><img src="image_199.png" /><figcaption>Figure 5 Three maximally diverse solutions within a maximum optimality gap of 0.02, using as distance measure the Manhattan distance of the vector of supply locations built by each solution.</figcaption></figure>

## 4.3 Robustness against plant shutdown

In practice, optimal solutions may lack robustness against uncertain events. There are many sources of uncertainty in our problem, from changes to input data (e.g., costs, demand or hardware technology) to loss of production due to plant shutdowns. Our system is not yet required to deal with the former, as our users are still determining how to build reliable future scenarios and/or probability distributions. It can however be used to find solutions that, in the event of any plant failure, guarantee the supply of a user-defined minimum percentage of the total demand of the network (denoted as new parameter GS ∈ [0..1], set to 0 as default). This allows users to compare optimal but non-robust solutions against robust but non-optimal ones. We achieved this by adding constraint (27) to the models to ensure that for every period p ∈ P, supply location sᵢ ∈ S, and product k ∈ K, the user-defined percentage GS of the total demand ∑_{dⱼ∈D} D_{p,j}^{k} for product k in period p is covered by the sum of the peak production capacity of the plants in other supply locations sᵢ' ∈ S \ {sᵢ}.

$$ \text{GS} \sum_{d_j \in D} D_{p,j}^k \le \sum_{s_{i'} \in S \setminus \{s_i\}} y_{p,i'}^h \quad \forall p, s_i, k, h \in \{\text{ELEC, COMP, LIQF}\} \qquad (27) $$

## 5 User's view

Our system (shown in Figure 6) connects two main components – optimisation and visualisation – via an application program interface (API). The optimisation calls MINIZINC with the model, input data and configuration to find and return its solutions. The visualisation calls the optimisation via the API and visualises the returned solutions on an interactive web-based dashboard application, referred to as the User Interface (UI). It was implemented using Plotly Dash [18], as requested by our industry partner.

Users experience our system mostly via its powerful UI, part of which appears in Figure 7. The UI allows domain-experts and non-expert users (e.g., decision makers) to load different input data, interactively turn some model constraints on/off, obtain solutions, compare them and explore them in detail. This not only helps them make decisions but also gain confidence and trust in the system's output. In addition, domain-experts and external tools can access the system's functionality independently of the UI via the API and the command line interface we implemented in Python. This allows them to automate the solving of instances, i.e. explore many scenarios with different sets of parameters.
<figure><img src="image_200.png" /><figcaption>Figure 6 The two components integrated by our system and connected by an API.</figcaption></figure>

## 5.1 Configuring an execution run

As described above, our system supports a wide range of execution modes (e.g., warm-starts, diversity, and robustness), which users can control via the UI. Fig. 7(A) shows one of the panes provided for this. In it, users can select a pre-configured scenario (list of supply and demand locations, their demand, and electricity markets) or upload a new one; modify the scenario by turning on/off the allowed markets (Utility, PPA, Wholesale); select a model profile (e.g., Ori-Exa-Tie) from a list of available ones; determine whether to run in diversity, robustness or traditional mode; and modify the default configuration values for each of those. In addition, users can select one of the available solvers; set a timeout and number of threads; and set various solver specific parameters. Once users finish configuring the run, solving can be triggered by pressing a button. If solving time is expected to be long, users can leave the dashboard and load the solution later using the session manager, which keeps track of each execution and reports whether a solution was found, no solution was found due to infeasibility, or the solver is still running.

## 5.2 Visualising and interacting with solutions

The UI displays solutions via the two tabs shown in Figures 7(B) and (C). Figure 7(B) focuses on overview information and the supply/demand network, for which it combines different sub-windows (or *cards*). Figure 7(C) focuses on the operational parts of one plant.

**Overview and Network Tab.** This tab (Figure 7(B)) comprises cards (frames) that contain controller and visualisation elements. With the controller element a specific solution from a set of multiple solutions (for example, solutions from diversity, robustness, or different set of model parameters) and time period can be loaded and visualised. Below the controller a stacked barchart summarises the $/kg cost (i.e. each cost component) for all active supply locations for all periods. The map in the centre prominently shows the supply/demand network, with markers representing supply locations and coloured circles for each demand location, with the colour matching the supply location they are linked to. For each supply location a radial chart with 4 segments (north, west, south east) shows the amount of hydrogen delivery for each direction, and an area overlay highlights the covered area using semi-transparent rectangles drawn between each linked supply/demand locations, indicating compactness, spread or amount of overlap of supply areas. Supply locations can be selected and its details viewed in a separate card. To the right of the map it shows details about the delivery to each demand location; below it shows the plant configuration represented as a flow diagram with details about installed hardware capacity and number of SKUs, including a stacked bar-chart for each type of hardware showing utilisation to installed capacity. For a
<figure><img src="image_201.png" /><figcaption>Figure 7 UI’s four main tabs: A) input tab, which allows detailed configuration of the model and data before execution; B) supply network as an interactive map with details of selected supply location; C) operational details of each supply location; and D) comparison view of multiple solutions.</figcaption></figure>

selected supply location, interactive elements allow changes to maximum and minimum daily production amount, disabling that particular supply location, or force a non-active location to be active. A re-solve of an instance containing the new data can then be triggered.

**Operation Tab.** For a supply location selected from the map card, this tab (Figure 7(C)) shows further details with focus on the operation side. The two cards in the middle show for each hour of the year (8760 hours) component usage and gas storage levels on the left, and electricity source utilisation (Utility, PPA, Wholesale) on the right. Other cards show the overall hydrolyser utilisation and waterfall charts of the investment costs. On the top left the same controller card is shown to select specific solutions and time periods.

## 5.3 Comparing solutions

Fig 7(D) shows part of the solution comparison pane for three solutions (one per column). For each solution, the pane shows a summary of the configuration of the run that created it, a breakdown of its objective value, and many other useful information, including the two kinds shown in the figure: a map of the plants it built and their hardware as a box-matrix plot (one plant per row). Each box represents one SKU of a particular type of hardware (identified by its colour) and capacity (identified by its size). This enables easy visual comparison of significant changes in plant locations and in their associated hardware. If a location appears in all solutions it is highlighted with a light green colour. Users can use the loading manager to compare any set of solutions obtained, for example, by selecting different model profiles, available sources, and diversity definitions; by modifying the robustness parameter, or by forcing a supply location to produce a given product.
## 5.4 Detecting and resolving conflict

It is not uncommon for users to create infeasible instances while exploring solutions by, for example, setting the minimum daily production for a plant to a higher value than its maximum production. This is a significant roadblock for the usability of any optimisation system. To address it, we use the conflict resolution method of [17], which allows users to obtain a visual representation of the conflicts that is easy to understand, select which of the conflicts to relax for a solution to be found, and obtain a solution to this relaxed instance that quantifies the minimum changes needed to restore feasibility for the original instance.

To implement this method the model needs to be changed to a) provide information about any constraint that can cause conflict, and b) soften these constraints to quantify the minimum changes needed. Both can be achieved in MINIZINC: a) by adding an annotation to each constraint that gives meaningful names to the constraint and its objects, and that tracks its values; b) by adding a slack variable to each constraint that quantifies the minimum changes (see [17] for details). The method also requires computing Minimum Unsatisfiable Sets (MUSes), where a MUS is a set of infeasible constraints in the instance such that removal of any one of the constraints makes the set feasible, and/or Minimum Correction Sets (MCS), where an MCS is a minimum set of constraints that if eliminated from the infeasible instance makes it feasible. Note that any minimum set that intersects all MUSes is an MCS. We have currently implemented one of the pathways proposed in [17] which, upon failure, uses FINDMUS (a MUS enumerator available in MINIZINC) to obtain all MUSes, shows them to the user in different formats, and allows them to manually select an MCS, thus ensuring there is a solution if the constraints in the MCS are relaxed. Currently, users must directly modify these constraints to obtain a feasible instance. We are in the process of implementing the automatic relaxation and resolution process of [17].

Figure 8 shows two of the formats proposed in [17] to show MUSes, which we have implemented. The pane on the left shows the MUS-graph, a compact way of showing the conflicting constraints (left of graph), the MUSes they appear in (centre), and the objects they involve (right). Each circle on the left of the graph represents one conflicting constraint, whose colour connects it to the name shown in the legend on the top right of the graph. For the example shown in the figure these include *minimum production* and *maximum production*. Constraints in the same box appear in the same MUSes. The links connect each group of constraints to the MUSes they appear in, and each MUS to the objects in any of its constraints, which in this case are particular supply and demand locations. If users select a constraint on the left of the graph, the constraint and all the MUSes in which it appears are highlighted with a red frame. Users know they have selected an MCS if the selected constraints highlight all MUSes. A geo-network view of the conflicts is given by showing them on a map (Figure 8 right) of the supply/demand locations. We overlay conflicts specific to a location with a frame in the colour of the conflict (e.g., light blue to show that this location is involved in a maximum production conflict), and with a coloured link for those involving a supply and demand location. Conflict overlays disappear when any of its constraints is selected in the MUS-graph, thus showing no conflict overlays once an MCS is selected.

# 6 Literature review

The problem of hydrogen production facility location and supply chain optimisation is widely studied; see Riera, Lima and Knio [16] for a recent survey of the field, both in terms of modelling the problem itself, and of optimisation strategies used. Our key takeaway from this survey is that there are myriads of different contexts in which to study the problem, from the
<figure><img src="image_202.png" /><figcaption>Figure 8 Two conflict visualisations: MUS-graph on the left, Geo-network on the right.</figcaption></figure>

“micro” perspective of a single plant, to the “macro” decisions around the energy delivery system for entire countries (including hydrogen as one of many types of renewable fuels). As such, it is difficult to identify any two papers that study exactly the same mathematical model (e.g., for optimisation benchmarking purposes).

Nevertheless, several papers stand out as closely related to the problem studied here: Ingason, Ingolfsson and Jensson [10] study a electrolysis-based facility location problem for Iceland, although they ignore the operational scheduling of the plant thanks to the assumption of constant-rate renewable sources of electricity (hydro and geothermal). Likewise, Štádlerová et al. [20] study a facility location problem in Norway including uncertainty in the demand and operational considerations such as minimum turn-down rates, although not at an hourly resolution. Demirhan et al. [4] study a multi-fuel, multi-chemical (hydrogen, methanol and ammonia) facility location problem with hourly resolution on the operation of the production plant to capture variability in renewable energy generation from wind and solar. Finally, Kim and Kim [13] also study a joint facility location and hourly operation problem for green hydrogen, although they consider only one time period for facility location without planning the facility’s expansion pathway. All these papers, however, focus mostly on the modelling expert’s view and ignore the domain expert’s one.

# 7 Conclusions and Future Work

This paper describes a Hydrogen Network Optimisation System (HyNetOS) designed to support energy companies in solving the complex combinatorial optimisation problem of producing and supplying hydrogen at minimum cost. HyNetOS integrates two supply network and facility operation models implemented in the high-level constraint modelling language MINIZINC, which natively supplies a range of tools to support the *human* decision-making process, such as finding a diverse sets of near-optimal solutions to present possible alternatives, and conflict resolution technology to explain infeasible instances. Further, HyNetOS supports our industry partner’s decision-making process by means of a powerful and interactive user interface that allows them to view, change, and compare solutions, and rapidly iterate on “what-if” scenarios with efficiently solvable approximate versions of our models. To increase confidence in the quality of the model, we applied a strategy of redundancy and rapid prototyping of model changes to identify the most efficient formulation, and provide robustness measures, to produce resilient solutions against hydrogen production failures.
The HyNetOS system was quantitatively and qualitatively evaluated against an alternative implementation using a direct modelling approach. It was selected as the preferred system based mainly on its significantly higher scores in qualitative criteria such as maintainability, extensibility, user interface, code quality and technology stack, as well as its additional features (diversity and robustness).

While deployed, the system is in constant evolution with many avenues for future work. One of the most pressing and complex is extending the system to produce robust solutions against uncertainties in electricity prices and availability of resources using techniques such as stochastic programming [19], sensitivity analysis [2] and Predict+Optimise [14]. Others include the integration of extra functionality, such as the production of ammonia or the consideration of carbon intensity, the integration of a simulation system for detailed operational modelling on a high resolution time scale, and the integration of a plant layout optimisation system such as [1] to generate optimal hydrogen facility layouts.

References

1. Gleb Belov, Tobias Czauderna, Maria Garcia de la Banda, Matthias Klapperstueck, Ilankaikone Senthooran, Mitch Smith, Michael Wybrow, and Mark Wallace. Process Plant Layout Optimization: Equipment Allocation. In John Hooker, editor, *Principles and Practice of Constraint Programming*, pages 473–489. Springer, 2018. doi:10.1007/978-3-319-98334-9_31.

2. M. W. Dawande and J. N. Hooker. Inference-based sensitivity analysis for mixed integer/linear programming. *Operations Research*, 48(4):505–660, 2000. doi:10.1287/opre.48.4.623.12420.

3. Jip J. Dekker. MiniZinc Python, 2023. URL: https://minizinc-python.readthedocs.io/en/latest/.

4. C. Doga Demirhan, William W. Tso, Joseph B. Powell, and Efstratios N. Pistikopoulos. A multi-scale energy systems engineering approach towards integrated multi-product network optimization. *Applied Energy*, 281:116020, 2021. doi:10.1016/j.apenergy.2020.116020.

5. Ibrahim Dincer and Canan Acar. Review and evaluation of hydrogen production methods for better sustainability. *International Journal of Hydrogen Energy*, 40(34):11094–11111, 2015.

6. Graeme Gange, Jeremias Berg, Emir Demirović, and Peter J. Stuckey. Core-guided and core-boosted search for CP. In Emmanuel Hebrard and Nysret Musliu, editors, *Integration of Constraint Programming, Artificial Intelligence, and Operations Research*, pages 205–221, Cham, 2020. Springer International Publishing. doi:10.1007/978-3-030-58942-4_14.

7. Gurobi Optimization, LLC. Gurobi Optimizer Reference Manual, 2023. URL: https://www.gurobi.com.

8. Q. Huangfu and J. A. J. Hall. Parallelizing the dual revised simplex method. *Mathematical Programming Computation*, 10(1):119–142, March 2018. doi:10.1007/s12532-017-0130-5.

9. IEA. Net zero by 2050. Technical report, IEA, Paris, May 2021. URL: https://www.iea.org/reports/net-zero-by-2050.

10. Helgi Thor Ingason, Hjalti Pall Ingolfsson, and Pall Jensson. Optimizing site selection for hydrogen production in Iceland. *International Journal of Hydrogen Energy*, 33(14):3632–3643, 2008. TMS07: Symposium on Materials in Clean Power Systems. doi:10.1016/j.ijhydene.2008.04.046.

11. Linnea Ingmar, Maria Garcia de la Banda, Peter J Stuckey, and Guido Tack. Modelling diversity of solutions. In *Proceedings of the AAAI Conference on Artificial Intelligence*, pages 1528–1535, 2020.

12. IPCC. Climate change 2022: Impacts, adaptation, and vulnerability. Technical report, Cambridge University Press, Cambridge, UK and New York, NY, USA, 2022. Contribution of Working Group II to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change. doi:10.1017/9781009325844.
21:18 Exploring Hydrogen Supply/Demand Networks

13 Minsoo Kim and Jiyong Kim. Optimization model for the design and analysis of an integrated renewable hydrogen supply (IRHS) system: Application to Korea's hydrogen economy. *International Journal of Hydrogen Energy*, 41(38):16613-16626, 2016. doi: 10.1016/j.ijhydene.2016.07.079.

14 Jayanta Mandi, Emir Demiroviç, Peter J. Stuckey, and Tias Guns. Smart predict-and-optimize for hard combinatorial optimization problems. In *Proceedings of the AAAI Conference on Artificial Intelligence*, pages 1603–1610, 2020. doi:10.1609/aaai.v34i02.5521.

15 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: 13th International Conference, CP 2007, Providence, RI, USA, September 23-27, 2007. Proceedings 13*, pages 529–543. Springer, 2007.

16 Jefferson A. Riera, Ricardo M. Lima, and Omar M. Knio. A review of hydrogen production and supply chain modeling and optimization. *International Journal of Hydrogen Energy*, 48(37):13731–13755, 2023. doi:10.1016/j.ijhydene.2022.12.242.

17 Ilankaikone Senthooran, Matthias Klapperstueck, Gleb Belov, Tobias Czauderna, Kevin Leo, Mark Wallace, Michael Wybrow, and Maria Garcia de la Banda. Human-Centred Feasibility Restoration in Practice. *Constraints*, 2023. (in publication). doi:10.1007/s10601-023-09344-5.

18 Shammamah Hossain. Visualization of Bioinformatics Data with Dash Bio. In Chris Calloway, David Lippa, Dillon Niederhut, and David Shupe, editors, *Proceedings of the 18th Python in Science Conference*, pages 126–133, 2019. doi:10.25080/Majora-7ddc1dd1-012.

19 Lawrence V. Snyder. Facility location under uncertainty: a review. *IIE Transactions*, 38(7):547–564, 2006. doi:10.1080/07408170500216480.

20 Šárka Štádlerová, Trygve Magnus Aglen, Andreas Hofstad, and Peter Schütz. Locating hydrogen production in Norway under uncertainty. In Jesica de Armas, Helena Ramalhinho, and Stefan Voß, editors, *Computational Logistics*, pages 306–321, Cham, 2022. Springer International Publishing. doi:10.1007/978-3-031-16579-5_21.

21 Mark van den Brand and Jan Friso Groote. Software engineering: Redundancy is key. *Science of Computer Programming*, 97:75–81, 2015. Special Issue on New Ideas and Emerging Results in Understanding Software. doi:10.1016/j.scico.2013.11.020.

22 H. Paul Williams. *Model building in mathematical programming*. John Wiley & Sons, 5th edition, 2013.
# Binary Constraint Trees and Structured Decomposability

Petr Kučera

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

Department of Theoretical Computer Science and Mathematical Logic, Faculty of Mathematics and Physics, Charles University, Prague, Czech Republic

## Abstract

A binary constraint tree (BCT, Wang and Yap 2022) is a normalized binary CSP whose constraint graph is a tree. A BCT constraint is a constraint represented with a BCT where some of the variables may be hidden (i.e. existentially quantified and used only for internal representation). Structured decomposable negation normal forms (SDNNF) were introduced by Pipatsrisawat and Darwiche (2008) as a restriction of decomposable negation normal forms (DNNF). Both DNNFs and SDNNFs were studied in the area of knowledge compilation. In this paper we show that the BCT constraints are polynomially equivalent to SDNNFs. In particular, a BCT constraint can be represented with an SDNNF of polynomial size and, on the other hand, a constraint that can be represented with an SDNNF, can be represented as a BCT constraint of polynomial size. This generalizes the result of Wang and Yap (2022) that shows that a multivalued decision diagram (MDD) can be represented with a BCT. Moreover, our result provides a full characterization of binary constraint trees using a language that is well studied in the area of knowledge compilation. It was shown by Wang and Yap (2023) that a CSP on $n$ variables of domain sizes bounded by $d$ that has treewidth $k$ can be encoded as a BCT on $O(n)$ variables with domain sizes $O(d^{k+1})$. We provide an alternative reduction for the case of binary CSPs. This allows us to compile any binary CSP to an SDNNF of size that is parameterized by $d$ and $k$.

**2012 ACM Subject Classification** Theory of computation → Automated reasoning; Theory of computation → Constraint and logic programming

**Keywords and phrases** Binary CSP, Binary Constraint Tree, Structured Decomposability, Structured DNNF, Polynomial Equivalence

Digital Object Identifier 10.4230/LIPIcs.CP.2023.22

## 1 Introduction

Constraint satisfaction problems (CSPs) offer an expressive and natural way of formulating problems. A CSP is a problem of checking satisfiability of a conjunction of constraints on variables with finite domains. Constraints can be represented in various ways, which include tables (see e.g. [2, 17, 18]) or multivalued decision diagrams (MDD, see e.g. [1, 5, 6]).

The representation using binary constraint trees was introduced in [27]. A BCT constraint is a constraint $c$ defined on a set of variables $x$ that is represented with a normalized binary CSP $P$ whose constraint graph is a tree. The CSP $P$ itself is defined on a set of variables $z$ which may include some *hidden* variables in addition to all the *original* variables from $x$. BCTs have a nice property that an arc consistency propagator can be used to check their consistency [12]. Any CSP can be turned into a binary one with an encoding such as dual encoding [11], hidden variable encoding [22], double encoding [24], or bipartite encoding [25].

Decomposable negation normal forms (DNNFs) were introduced in [7] as a tractable language for knowledge representation. Structured DNNFs (SDNNF) were introduced in [20]. The definition of SDNNFs is based on the notion of a *v-tree* which is a rooted binary tree whose leaves are in one-to-one correspondence with the constraint variables (both original and hidden). The conjunction gates in an SDNNF are then required to respect a particular v-tree (see definitions 5 and 6 in Section 2.3 for more details). The structural requirements imposed

<figure><img src="image_204.png" /></figure>
on SDNNFs allow for instance a polynomial time construction of an SDNNF representing the conjunction of two SDNNFs that respect the same v-tree – something that is not possible for DNNFs without any structural requirements. Although DNNFs were introduced as a representation of functions on boolean variables, they were also considered as a representation of constraints on variables with finite domains [1, 14].

Both BCT constraints and SDNNFs have a structure based on a tree. We use this
similarity to show the main result of this paper: BCT constraints and SDNNF constraints are
polynomially equivalent. In particular, we show a polynomial time transformation of a BCT
constraint into an SDNNF and also a polynomial time transformation of an SDNNF into a
BCT constraint. The polynomial equivalence of BCTs with SDNNFs offers a characterization
of BCTs by a language of SDNNFs which has been extensively studied in the area of knowledge
compilation [3, 20, 21, 23]. Our result also generalizes the previous construction of a BCT
constraint for an MDD described in [27]. It was shown in [26] that BCTs are strictly more
succinct than MDDs. This also follows from a combination of our result with the fact that
SDNNFs are strictly more succinct than MDDs by [20].

Recently, [28] studied BCTs from the perspective of knowledge compilation together with several other languages that are being used to represent ad-hoc constraints. The authors studied BCTs with respect to the queries and transformations considered in the knowledge compilation map [10] and showed that BCTs allow answering consistency, clausal entailment and model enumeration queries in polynomial time which is (unsurprisingly) the same as in the case of structured DNNFs [20]. The authors of [28] also studied BCTs with respect to transformations. If the input BCTs or SDNNFs are required to have the same tree structure, then they allow polynomial-time bounded conjunction, unbounded disjunction, forgetting any number of variables, and conditioning [20, 28]. Interestingly, [20] only considers the case of combining SDNNFs that respect the same v-tree while [28] also considers the case of combining BCTs that are not required to have the same tree structure. In this case BCTs do not allow polynomial time bounded conjunction, they do not allow an unbounded disjunction, and the case of bounded disjunction is unresolved in [28]. We believe that our result might help to resolve the case of bounded disjunction for BCTs, because it might be easier to reason about a disjunction of two SDNNFs than BCTs. It is also worth mentioning that according to [20], AOMDDs introduced in [19] are strictly less succinct than SDNNFs and thus also strictly less succinct than BCTs. This already answers one of the questions posed in [28].

Our transformation of a BCT constraint into an SDNNF leads to a smooth SDNNF. It
is thus possible to use a domain consistency propagator for smooth DNNFs described in [14]
as a domain consistency propagator for BCT constraints. An encoding of BCT constraints
with propagation complete CNF formulas was described in [26]. Various CNF encodings
of DNNF theories were considered in [1] and a propagation complete encoding of smooth
DNNFs was introduced in [16]. Our result thus offers an alternative way of reducing BCT
constraints to a CNF encoding.

If CNF $\varphi$ has treewidth $k$, then it can be compiled to an SDNNF of size that is parame-
terized by $k$ by the construction described in [21]. In particular, if $\varphi$ has $n$ variables and $m$
clauses, then an equivalent SDNNF can be constructed in time $O(nm2^k)$. We can obtain
a similar result also for binary CSPs, but we have to take into account also the domain
sizes. It was shown in [28] that if $P$ is a CSP on $n$ variables of domain size $d$ that has
treewidth $k$, then it can be encoded as a BCT with $O(n)$ variables with domain size $d^{k+1}$.
The construction in [28] uses the encoding described in [11]. In addition, if $P$ is a binary CSP,
its consistency can be checked in time $O(nd^{k+1})$ by [13]. To have a complete compilation
procedure of a binary CSP into an SDNNF, we provide a direct reduction of a binary CSP
to a BCT parameterized by the treewidth and the domain size. The bound we obtain is the same as in [28], so our result is not really new in this sense, but we obtain a slightly better bound on the size of an SDNNF constructed for the given binary CSP than if we would simply combine the bound of [28] with our construction of an SDNNF.

The paper is organized as follows. We introduce the necessary notation in Section 2, including the definitions of BCTs and structured DNNFs. In Section 3, we show that a BCT constraint can be represented with an SDNNF. The transformation of an SDNNF to a BCT is described in Section 4. A transformation of a binary CSP with bounded treewidth into an SDNNF is described in Section 5. Section 6 closes the paper with a few concluding remarks.

# 2 Definitions

In this section, we shall recall the necessary notation and notions used in the paper.

## 2.1 BCT Constraint

We use a notation adapted from [27] where binary constraint trees were introduced.

A CSP $P$ is a pair $(\mathbf{x}, C)$ where $\mathbf{x}$ is a set of variables and $C$ is a set of constraints. Each variable $x$ has a finite domain denoted $dom(x)$. A *literal* on a variable $x$ is a variable value assignment $(x, a)$. A *tuple* over a set of variables $\{x_{i_1}, x_{i_2}, \dots, x_{i_r}\}$ is a *set of literals* $\{(x_{i_1}, a_1), (x_{i_2}, a_2), \dots, (x_{i_r}, a_r)\}$. Each constraint $c_j$ has a constraint scope $scp(c_j) \subseteq \mathbf{x}$ and a relation $rel(c_j)$ defined by a set of tuples over $scp(c_j)$. A constraint $c$ is a *binary constraint* if $|scp(c)| = 2$ and it is a *unary constraint* if $|scp(c)| = 1$. A CSP $P$ is called a *binary CSP* if it consists of binary and unary constraints. A binary CSP is *normalized* if its constraints have pairwise different scopes. Given any set of variables $\mathbf{z}$ and literals $\tau$, we use $\tau[\mathbf{z}] = \{(x, a) \in \tau \mid x \in \mathbf{z}\}$ to denote a subset of $\tau$, while $T[\mathbf{z}] = \{\tau[\mathbf{z}] \mid \tau \in T\}$ is the *projection* of a set of tuples $T$ on $\mathbf{z}$. A tuple $\tau$ over $\mathbf{x}$ is a solution of $P$ if $\tau[scp(c)] \in rel(c)$ for all constraints $c \in C$ and $a \in dom(x)$ for all $(x, a) \in \tau$. We use $sol(\mathbf{x}, C)$ (or $sol(P)$) to denote the set of all solutions of $P$. We also say that $P$ is *satisfied* by its solution and that a solution of $P$ satisfies all constraints in $C$. A *support* of a value $a \in dom(x)$ in a constraint $c$ is a tuple $\tau \in rel(c)$ such that $(x, a) \in \tau$ and $b \in dom(y)$ for all $(y, b) \in \tau$.

**Definition 1** ([27]). A Binary Constraint Tree (BCT) is a normalized binary CSP whose constraint graph is a tree. A BCT constraint $c$ is a pair $(\mathbf{x}, P)$ such that $P = (\mathbf{z}, C)$ is a BCT, $scp(c) = \mathbf{x} \subseteq \mathbf{z}$, and $rel(c) = sol(\mathbf{z}, C)[\mathbf{x}]$. A tree binary encoding (TBE) of a constraint $c^*$ is a BCT $P = (\mathbf{z}, C)$ such that the BCT constraint $(scp(c^*), P)$ has the same constraint relation as $c^*$ where the variables in $scp(c^*)$ and $\mathbf{z} \setminus scp(c^*)$ are called the original and hidden variables, respectively.

**Example 2.** Let us consider a BCT constraint $c^* = (\mathbf{x}, P)$ on three variables $\mathbf{x} = \{x_1, x_2, x_3\}$ where $P = (\mathbf{z}, C)$ is a BCT described as follows. We have one hidden variable $y$ in $P$, i.e. $\mathbf{z} = \{x_1, x_2, x_3, y\}$, and three constraints $C = \{c_1, c_2, c_3\}$ with $scp(c_i) = \{x_i, y\}$, $i = 1, 2, 3$. The domain of all variables (original and hidden) is $\{1, 2, 3\}$. For $i = 1, 2, 3$, we set $rel(c_i) = \{( (x_i, a), (y, b) ) \mid a \neq b \}$. That is, $c_i$ enforces that $y$ has a different value from $x_i$ in any solution to $c^*$. Altogether, $c^*$ is equal to the negation of the *alldifferent* constraint over the variables $x_1, x_2, x_3$.

A general construction of a BCT representing the negation of the *alldifferent* constraint over variables $x_1, \dots, x_r$ with domains $D = \{1, \dots, r\}$ was described in [27] where the authors also noted that the size of the MDD representing the constraint is exponential in $r$.
<figure><img src="image_205.png" /><figcaption>(a) The binary constraint tree of P.</figcaption></figure>

<figure><img src="image_206.png" /><figcaption>(b) Relation of constraint $c_1$ ($z_1 \le z_3$).</figcaption></figure>

<figure><img src="image_207.png" /><figcaption>(c) Relation of constraint $c_2$ ($z_2 \le z_3$).</figcaption></figure>

<figure><img src="image_208.png" /><figcaption>(d) Relation of constraint $c_3$ ($z_3 < z_4$).</figcaption></figure>

■ **Figure 1** A binary constraint tree $P = (\mathbf{z}, C)$ from Example 3.

To demonstrate the techniques described in the paper, we shall consider a simple constraint that is a bit less symmetrical than the negation of the *alldifferent* constraint on three variables.

► **Example 3.** Figure 1 shows a BCT $P = (\mathbf{z}, C)$ that is defined on variables $\mathbf{z} = (z_1, z_2, z_3, z_4)$ with domains $dom(z_i) = \{0, 1, 2\}$ for all $i = 1, \dots, 4$. $C$ consists of three constraints. Constraints $c_1$ and $c_2$ represent inequalities $z_1 \le z_3$ and $z_2 \le z_3$ respectively and its relation is shown in figures 1b and 1c. Constraint $c_3$ represents inequality $z_3 < z_4$ and its relation is shown in Figure 1d. Note that literals $(z_4, 0)$ and $(z_3, 2)$ do not have support in $c_3$.

Let us now consider a constraint $c^*$ with scope $scp(c^*) = \{x_1, x_2, x_3\}$ where $dom(x_i) = \{0, 1, 2\}$ for $i = 1, 2, 3$ and the set of tuples $rel(c^*)$ represents inequality $max(x_1, x_2) < x_3$. If we identify variables $z_1$, $z_2$, and $z_4$ with $x_1$, $x_2$, and $x_3$ respectively, then $P$ is a tree binary encoding of $c^*$ in which $z_1$, $z_2$, and $z_4$ are original variables and $z_3$ is a hidden variable.

Note that the hidden variable $z_3$ is not actually needed for the constraint representation. We keep it to demonstrate how a hidden variable can be later forgotten in an SDNNF. We may also observe that literals $(z_4, 0)$ and $(z_3, 2)$ do not have a support in constraint $c_3$. We shall see later how this situation is dealt with during the construction of an SDNNF representing $c^*$.

## 2.2 DNNF

The notion of a DNNF was introduced in [7] as a restriction of NNF. We consider a multivalued variant that was used for instance in [14, 16]. This form is suitable for using DNNFs to represent constraints.

Consider a set of variables $\mathbf{x} = \{x_1, \dots, x_n\}$ with a finite domain $dom(\mathbf{x})$ for each $\mathbf{x} \in \mathbf{x}$. A sentence in *negation normal form* (NNF) $D$ is a rooted DAG with vertices $V$, root $\rho \in V$, the set of edges $E$, and the set of leaves $L \subseteq V$. The inner vertices (also called *gates*) are labeled with logical connectives $\wedge$ or $\vee$. Each edge $(v, u)$ in $D$ connects an inner vertex $v$ labeled $\wedge$ or $\vee$ with one of its inputs $u$. The edge is directed from $v$ to $u$, so the inputs of a
vertex are its successors (or child vertices). The leaves are labeled with literals of variables **x**, i.e. each leaf is labeled with a literal (*x*<sub>*i*</sub>, *a*) on a variable *x*<sub>*i*</sub> ∈ **x** and a value *a* ∈ dom(*x*<sub>*i*</sub>). We assume that each literal (*x*<sub>*i*</sub>, *a*) is used as a label of at most one leaf. Some of the literals may be missing in *D*, however, we assume that for each *i* = 1, ..., *n* at least one leaf of *D* is labeled with a literal on variable *x*<sub>*i*</sub>. For technical reasons, we also allow to label leaves with constants 0 or 1. If a NNF is nonempty, constants can always be propagated and these constants are thus needed only on a NNF without any variables.

If all the constraint variables $x_i \in \mathbf{x}$ are boolean (i.e. $dom(x_i) = \{0,1\}$) and we identify
literal $(x_i, 1)$ with the propositional literal $x_i$ and literal $(x_i, 0)$ with the propositional literal
$
\neg x_i$, then we obtain the usual definition of a NNF for representing a boolean function.

Assume that *c* is a constraint with the scope *scp*(*c*) = **x** and that *D* is a NNF defined on the variables **x**. We say that *D* represents constraint *c* if for every tuple *τ* over variables **x** we have that *τ* ∈ *rel*(*c*) if and only if *D* evaluates to true on the tuple *τ*. Evaluating *D* on *τ* is done in a straightforward manner, we simply set the leaves (*x*<sub>*i*</sub>, *a*) ∈ *τ* to true and the remaining leaves to false, then we use the usual semantic of the circuit *D* to get the value on this assignment.

Following [14], we define the decomposability and smoothness properties with respect to constraint variables $x_1, \dots, x_n$. For a vertex $v \in V$, let us denote $var(v) \subseteq \mathbf{x}$ the set of variables in the subcircuit of $D$ rooted at $v$. More precisely, a variable $x_i \in \mathbf{x}$ belongs to $var(v)$ if and only if there is a directed path from $v$ to a leaf labeled with a literal $(x_i, a)$ for a value $a \in \text{dom}(x_i)$. We have by assumption that $var(a) = \mathbf{x}$.

► **Definition 4.** We define the following structural restrictions of NNFs.

- We say that NNF *D* is decomposable (DNNF), if for every vertex *v* = *v*₁ ∧ ... ∧ *v*ₖ the sets of variables var(*v*₁), ..., var(*v*ₖ) are pairwise disjoint.

- We say that DNNF *D* is smooth if for every vertex *v* = *v*₁ ∨ ... ∨ *v*ₖ we have var(*v*) = var(*v*₁) = ... = var(*v*ₖ).

Assume that *D* is a DNNF representing constraint *c* with scope *scp*(*c*) = **z**. Let **x** ⊆ **z**
and **y** = **z** \ **x**. By *forgetting* variables **y** in *D* we mean the construction of a DNNF *D*'
that represents the constraint *c*' which is a projection of *c* on variables **x**. In particular,
*scp*(*c*') = **x** and *rel*(*c*') = *rel*(*c*)[**x**]. Forgetting can be done efficiently on a DNNF, we simply
replace every literal (*y*, *a*) with constant 1 for all *y* ∈ **y** and *a* ∈ dom(*y*) [8].

2.3 Structured DNNF

Structured DNNFs were introduced in [20]. Structured decomposability is based on the
notion of a v-tree defined as follows.

**Definition 5 ([20]).** A v-tree for a set of variables **x** is a full, rooted binary tree whose leaves are in one-to-one correspondence with the variables in **x**.

Given a node *t* of a v-tree *T*, we denote *var*(*t*) the set of variables associated with the
leaves in the subtree of *T* rooted at *t*. We also denote *var*(*T*) = *var*(*σ*) where *σ* is the root of
*T*. For a non-leaf node *t*, we use *t*<sup>*l*</sup> (*t*<sup>*r*</sup>) to denote the left (right) child node of *t*. For the rest
of this paper, we will assume that each conjunction in a DNNF has exactly two non-constant
inputs, while a disjunction can have any number of inputs. This is a technical assumption
used in [20] mainly to simplify the definition of a SDNNF, since with this assumption, it is
enough to consider only binary v-trees. Note also that we can make this assumption without
loss of generality due to the associativity and commutativity of conjunction.
<figure><img src="image_209.png" /><figcaption>■ **Figure 2** An example of an SDNNF $D_P$ and the corresponding v-tree $T$. This particular SDNNF is the result of our construction on the BCT $P$ from Example 3. The labels of the nodes mark the steps of our construction.</figcaption></figure>

► **Definition 6** ([20]). A DNNF $D$ respects a v-tree $T$ if for every conjunction $v = v_1 \wedge v_2$ in $D$, there is a node $t$ in $T$, such that $var(v_1) \subseteq var(t^l)$ and $var(v_2) \subseteq var(t^r)$.

Let $v$ be a node in a DNNF $D$ that respects a v-tree $T$. The decomposition node ($d$-node) of $v$ is defined as the deepest node $d$ in $T$ such that $var(v) \subseteq var(d)$.

► **Definition 7** ([20]). A DNNF that respects a given v-tree $T$ is denoted as $DNNF_T$. Moreover, the language of structured DNNFs (SDNNF) consists of all $DNNF_T$ for any v-tree $T$.

Given a $DNNF_T D$, we can construct an equivalent smooth $DNNF_T D'$ in quadratic time [23]. It means that we can always assume that the input $DNNF_T$ is smooth.

► **Example 8.** Figure 2 shows an example of a smooth $DNNF_T$. In particular, $DNNF_T D_P$ on Figure 2a respects v-tree $T$ on Figure 2b. $DNNF_T D_P$ represents the BCT constraint $(z, P)$ where $P$ is the BCT from Example 3.

For the construction of an SDNNF representing a BCT constraint, we will need the following operation for composing two v-trees. Assume $T_1$ and $T_2$ are two v-trees on disjoint sets of variables, i.e. $var(T_1) \cap var(T_2) = \emptyset$. Then $T = T_1 \circ T_2$ denotes the v-tree with a newly added root $\sigma$ whose left child node $\sigma^l$ is set to the root of $T_1$ and the right child node $\sigma^r$ is set to the root of $T_2$. It follows that $var(T) = var(T_1) \cup var(T_2)$.
### 3 Compiling a BCT Constraint into a Structured DNNF

We shall show in this section that we can construct an SDNNF representing a given BCT constraint in polynomial time.

**► Theorem 9.** Let $c^* = (\mathbf{x}, P)$ be a BCT constraint where $P = (\mathbf{z}, C)$ is a BCT. Then there is a smooth SDNNF $D$ representing $c^*$ with $O(m\mathbf{z})$ nodes and $O(m\mathbf{z}^2)$ edges where $m = |\mathbf{z}|$ and $d = \max_{z_i \in \mathbf{z}} |\mathrm{dom}(z_i)|$.

We will describe the construction of $D$ in the rest of this section, thus proving Theorem 9. Assume that $n = |\mathbf{x}|$ and $\mathbf{z} = \{z_1, \dots, z_m\}$ where $\mathbf{x} \subseteq \mathbf{z}$ and $m \ge n$ is the number of all variables. We assume that $|\mathrm{dom}(z_i)| \le d$ for every $i = 1, \dots, m$.

The construction proceeds in two steps. First, we describe a construction of a SDNNF $D_P$ representing $P$. A SDNNF $D$ that represents $c^*$ then originates from $D_P$ by forgetting the hidden variables $\mathbf{y} = \mathbf{z} \setminus \mathbf{x}$. This step can be done efficiently by [20].

Let $G$ be the constraint graph of $P$. $G$ is a tree with the set of nodes $\mathbf{z}$, each edge corresponds to a single constraint from $C$. Let $G^+$ denote a directed tree that originates from $G$ by picking an arbitrary node as a root and directing all edges from the root towards the leaves. Let us assume that the nodes $z_1, \dots, z_m$ are ordered in a reverse topological order with respect to $G^+$. It means that $z_m$ is the root and if $(z_i, z_j)$ is an edge in $G^+$, then $i > j$. See Figure 1a for an example.

For every $i = 1, \dots, m$, let us consider the subtree $G_i$ of $G^+$ rooted at $z_i$. Let $C_i \subseteq C$ denote the set of constraints corresponding to the edges of $G_i$. Denote $\mathbf{z}_i = \bigcup_{c \in C_i} \mathrm{scp}(c)$. In this way, we have defined BCT $P_i = (\mathbf{z}_i, C_i)$. For every value $a \in \mathrm{dom}(z_i)$, we also define BCT $P_{i,a}$ as a restriction of $P_i$ to the solutions that contain literal $(z_i, a)$. This can be best understood as adding a unary constraint with scope $z_i$ and a single relation $(z_i, a)$ to $C_i$. Another way of looking at it is restricting the relation of every constraint $c \in C_i$ with $z_i \in \mathrm{scp}(c)$ to the tuples containing $(z_i, a)$ and setting $\mathrm{dom}(z_i)$ to $\{a\}$.

The algorithm proceeds for every $i = 1, \dots, m$ in order and constructs for every $a \in \mathrm{dom}(z_i)$ a SDNNF $D_{i,a}$ representing BCT constraint $(\mathbf{z}_i, P_{i,a})$ and a v-tree $T_i$ that is respected by $D_{i,a}$ using the following steps:

(A1) If $z_i$ is a leaf (no edges leave $z_i$ in $G^+$), then $P_i$ has no constraints and $P_{i,a}$ has the domain of $z_i$ restricted to the single value $a$. DNNF $D_{i,a}$ is a single node labeled with literal $(z_i, a)$ and v-tree $T_i$ is a single node labeled with variable $z_i$.

(A2) Assume that $z_i$ is not a leaf and it has $k$ outgoing edges $(z_i, z_{i_1}), \dots, (z_i, z_{i_k})$ associated with constraints $c_1, \dots, c_k \in C$. For every $p = 1, \dots, k$ we have that $i_p < i$, because the nodes are processed in a reverse topological order, and thus we have already constructed $D_{i_p,b}$ and $T_{i_p}$ for each $b \in \mathrm{dom}(z_{i_p})$. Let us now construct $D_{i,a}$ for $a \in \mathrm{dom}(z_i)$ in the following two steps.

(A2a) For every $p = 1, \dots, k$, define $E_{i,a}^{i_p}$ as follows:

$$E_{i,a}^{i_p} = \bigvee_{\{(z_i,a),(z_{i_p},b)\} \in \mathrm{rel}(c_p)} D_{i_p,b}.$$
(A2b) We construct $D_{i,j}$ as a conjunction of literal $(z_i, a)$ with all DNNFs $E_{i,a}^{ip}$ for $p = 1, \dots, k$. However, since we only allow conjunctions with two inputs, we construct $D_{i,a}$ in $k+1$ steps as follows.

$$D_{i,a}^{(1)} = E_{i,a}^{i_1} \qquad (1)$$

$$D_{i,a}^{(p)} = D_{i,a}^{(p-1)} \wedge E_{i,a}^{ip} \quad \text{for } p = 2, \dots, k \qquad (2)$$

$$D_{i,a} = (z_i, a) \wedge D_{i,a}^{(k)} \qquad (3)$$

In addition, we define $T_i$ as follows:

$$T_i^{(1)} = T_{i_1} \qquad (4)$$

$$T_i^{(p)} = T_i^{(p-1)} \circ T_{i_p} \quad \text{for } p = 2, \dots, k \qquad (5)$$

$$T_i = z_i \circ T_i^{(k)} \qquad (6)$$

Variable $z_i$ is identified with a tree consisting of a single leaf labeled with $z_i$ in step (6).

Once we have constructed $D_{m,a}$ for every $a \in \text{dom}(z_m)$, we compose them to obtain $D_P$ that respects v-tree $T = T_m$ as follows:

$$D_P = \bigvee_{a \in \text{dom}(z_m)} D_{m,a}. \qquad (7)$$

Intuitively, $E_{i,a}^{ip}$ represents the fact that $(z_i, a)$ has a support $((z_i, a), (z_{i_p}, b))$ in $c_p$. Moreover, the constraints below $z_{i_p}$ in $G^+$ can be satisfied with the value of $z_{i_p}$ set to $b$. SDNNF $D_{i,a}$ represents the models of all constraints that correspond to the edges in the subtree of $G^+$ rooted at $z_i$ and that contain literal $(z_i, a)$. If $(z_i, a)$ does not have a support in $c_p$ for some $p = 1, \dots, k$, then the empty disjunction $E_{i,a}^{ip}$ is equal to constant 0 and so is $D_{i,a}$. However, it is possible that $E_{i,a}^{ip}$ is inconsistent even if $(z_i, a)$ has a support $((z_i, a), (z_{i_p}, b))$ in $c_p$, but $D_{i_p,b}$ is inconsistent for every such $b$.

**Example 10.** Figure 2 shows the result of the construction when applied to the BCT $P$ from Example 3. Note that there is no leaf labeled with literal $(z_4, 0)$, because $(z_4, 0)$ has no support in constraint $c_3$. For the same reason, there is no leaf labeled with literal $(z_3, 2)$. Consequently, there are no leaves labeled with literals $(z_1, 2)$ and $(z_2, 2)$. It is worth noting that value 2 would be removed from the domains of variables $z_1, z_2$, and $z_3$ and value 0 would be removed from the domain of $z_4$ when enforcing arc consistency. In this way, enforcing arc consistency is part of the construction.

If variable $z_3$ would be forgotten from $D_P$, we would obtain a SDNNF representing the constraint $c^*$ from Example 3. This would amount to replacing leaves labeled with literals $(z_3, 0)$ and $(z_3, 1)$ with constant 1. In this case, it just means removing these leaves altogether. Afterwards, we could simplify the SDNNF by removing the trivial gates with a single input, the result of this simplification can be seen in Figure 3.

We will now show that $D_P$ is a smooth SDNNF of polynomial size that represents $P$. We will start by showing that $D_P$ is a smooth SDNNF.

**Lemma 11.** $D_P$ is a smooth SDNNF that respects v-tree $T_m$.

**Proof.** Let us first show that $D_P$ is an SDNNF that respects v-tree $T_m$. We will proceed by induction on $i = 1, \dots, m$. If $z_i$ is a leaf of $G^+$, then by step (A1), $D_{i,a}$ consists of a single leaf node for every $a \in \text{dom}(z_i)$. It follows that $D_{i,a}$ respects $T_i$ which also consists of a single leaf node.
<figure><img src="image_210.png" /><figcaption>■ **Figure 3** Simplified SDNNF for the BCT constraint $c^*$ from Example 3 representing constraint $\max(x_1, x_2) < x_3$. The SDNNF originated from the SDNNF in Figure 2a by forgetting $z_3$, identifying $x_1 = z_1, x_2 = z_2, x_3 = z_4$, and removing gates with a single input.</figcaption></figure>

Let us now assume that $z_i$ is not a leaf, in particular $i > 1$. Let us assume that $z_i$ has $k$ outgoing edges to nodes $z_{i_1}$ to $z_{i_k}$. Assume a value $a \in \text{dom}(z_i)$. By induction hypothesis, for every $p = 1, \dots, k$ and every $b \in \text{dom}(z_{i_p})$ we have that $D_{i_p, b}$ is a SDNNF respecting v-tree $T_{i_p}$. It follows that $E_{i,a}^{i_p}$ constructed in step (A2a) is a SDNNF respecting $T_{i_p}$. We have that $\text{var}(E_{i,a}^{i_p}) = \mathbf{z}_{i_p}$. Since the subtrees rooted at nodes $z_{i_1}, \dots, z_{i_k}$ are pairwise disjoint, the same is true for sets $\mathbf{z}_{i_1}, \dots, \mathbf{z}_{i_p}$. Moreover, variable $z_i$ is not in any of these sets. Therefore, $D_{i,a}$ constructed in step (A2b) is a DNNF. The construction of the tree $T_i$ proceeds in a way similar to the construction of $D_{i,a}$, and thus $D_{i,a}$ respects $T_i$. In particular, the node of $T_i$ introduced in (5) is the d-node of the conjunction (2) for the same value of $p$, and the node introduced in (6) is the d-node of the conjunction (3).

$D_P$ is constructed in step (7) as a disjunction of $D_{m,a}$, $a \in \text{dom}(z_m)$. As each of these SDNNFs respects $T_m$, the same is true for $D_P$.

Let us now show the smoothness. Assume that $D_{i,a}$ is nontrivial, i.e. it is not just a single leaf labeled with 0. We show by induction that then $D_{i,a}$ is a SDNNF that depends on all variables in $\mathbf{z}_i$. This is true for the leaves. If $z_i$ is not a leaf, $D_{i,a}$ is constructed in step (A2b). By induction hypothesis used on each $D_{i_p, b}$ we get that $E_{i,a}^{i_p}$ is a smooth disjunction that depends on all variables in $\mathbf{z}_{i_p}$. Thus also $D_{i,a}$ depends on all variables in $\mathbf{z}_i = \{z_i\} \cup \bigcup_{p=1}^k \mathbf{z}_{i_p}$. It follows that also the disjunction introduced in the final step (7) is smooth. ◀

Now, let us estimate the size of $D_P$.

► **Lemma 12.** SDNNF $D_P$ has $O(md)$ nodes and $O(md + s)$ edges where $s = \sum_{c \in C} |\text{rel}(c)|$.
**Proof.** For every $i = 1, \dots, m$, we add one disjunction for each value $a \in \text{dom}(z_i)$ in step (A2a) and $k$ conjunctions in step (A2b) assuming $z_i$ is not a leaf. One more disjunction gate is added in the final step (7). Altogether, we thus add $O(md)$ gates to $D_P$. Each conjunction gate has at most two inputs, thus $O(md)$ edges are leaving the conjunction gates. For every constraint $c_p$ with scope $\{z_i, z_{i_p}\}$, every tuple $((z_i, a), (z_{i_p}, b))$ adds one edge leaving disjunction gate $E_{i,a}^{i_p}$. The total number of edges leaving the disjunction gates added in step (A2a) is thus at most $s$. The disjunction gate added in the final step has at most $|\text{dom}(z_i)| \le d$ inputs. Altogether, we have $O(md + s)$ edges in $D_P$. ◀

It remains to show that $D_P$ represents $P$.

► Lemma 13. SDNNF $D_P$ represents $P$.

**Proof.** We shall first show by induction on $i$ that $D_{i,a}$ represents $P_{i,a}$ for every $i = 1, \dots, m$ and $a \in \text{dom}(z_i)$. This is true for leaves (and thus also for $i = 1$), because $P_{i,a}$ does not have any constraints, it depends only on $z_i$, and the domain of $z_i$ is restricted to the value $a$. A single node labeled with literal $(z_i, a)$ added in step (A1) is thus a correct representation of $P_{i,a}$ in this case.

Let us now consider a variable $z_i$ with outgoing edges $(z_i, z_{i_1}), \dots, (z_i, z_{i_k})$ associated with constraints $c_1, \dots, c_k \in C$ where $k \ge 1$. $C_i$ is thus a disjoint union of $C_{i_p}, p = 1, \dots, k$ with $\{c_1, \dots, c_k\}$. Let us also consider a value $a \in \text{dom}(z_i)$. Let us assume by induction hypothesis that each $D_{i_p, b}$ represents $P_{i_p, b}$ for every $b \in \text{dom}(z_{i_p})$. Recall that the scope of $P_{i,a}$ is the set $\mathbf{z}_i$ of variables in the subtree of $G^+$ rooted at $z_i$. Let $\tau$ be a tuple of variables $\mathbf{z}_i$ and let us fix some $a \in \text{dom}(z_i)$.

Let us first assume that $\tau \in \text{sol}(P_{i,a})$. It follows that $(z_i, a) \in \tau$. We will show that $\tau$ satisfies $D_{i,a}$. Let us consider literals $(z_{i_1}, b_1), \dots, (z_{i_k}, b_k) \in \tau$. For every $p = 1, \dots, k$ we have that $\tau$ satisfies $P_{i_p}$. It satisfies $P_{i_p, b_p}$ as well since $(z_{i_p}, b_p) \in \tau$. By induction hypothesis, circuit $D_{i_p, b_p}$ represents $P_{i_p, b_p}$ and thus it evaluates to true on $\tau$. By definition of $P_{i,a}$, $\tau$ satisfies $P_i$ and thus also constraint $c_p$. It follows that $\{(z_{i_p}, b_p), (z_i, a)\} \in \text{rel}(c_p)$ and thus, by step (A2a), also $E_{i,a}^{i_p}$ evaluates to true. Since this holds for every $p = 1, \dots, k$ and $(z_i, a) \in \tau$, we have by step (A2b) that $D_{i,a}$ evaluates to true on $\tau$.

Let us now assume that $D_{i,a}$ evaluates to true on $\tau$ and let us show that $\tau \in \text{sol}(P_{i,a})$. We have $(z_i, a) \in \tau$ by (3), it remains to show that $\tau \in \text{sol}(P_i)$. Let $p \in \{1, \dots, k\}$ be arbitrary. We have by (1) to (3) that $E_{i,a}^{i_p}$ evaluates to true. By (A2a), we have that $D_{i_p, b_p}$ evaluates to true on $\tau$ for some $\{(z_i, a), (z_{i_p}, b_p)\} \in \text{rel}(c_p)$. Induction hypothesis implies $\tau[\mathbf{z}_{i_p}] \in \text{sol}(P_{i_p, b_p})$ and thus also $(z_{i_p}, b_p) \in \tau$. Together with $(z_i, a) \in \tau$ we obtain that $c_p$ is satisfied by $\tau$. In addition, all constraints in $P_{i_p}$ are satisfied by $\tau$. Since this holds for every $p = 1, \dots, k$, we get that all constraints of $P_i$ are satisfied and thus $\tau \in \text{sol}(P_i)$.

Let us now show that $D_P$ represents $P$. If $\tau$ is a tuple satisfying $P$ and $(z_m, a) \in \tau$, then $\tau$ satisfies $P_{m,a}$. It follows that $D_{m,a}$ evaluates to true and that $D_P$ evaluates to true as well. If, on the other hand, $D_P$ evaluates to true on $\tau$, then $D_{m,a}$ evaluates to true for some $a \in \text{dom}(z_m)$. Thus $\tau$ is a solution of both $P_{m,a}$ and $P = P_m$. ◀

Theorem 9 now follows from the above propositions.

**Proof of Theorem 9.** The SDNNF $D$ representing $c^*$ originates from $D_P$ by forgetting variables $\mathbf{y} = \mathbf{z} \setminus \mathbf{x}$. This step can be performed in polynomial time by [20] by replacing the literals on variables from $\mathbf{y}$ with constants 1 and then propagating these constants. Note that in $D_P$, this is equivalent to removing the corresponding leaves which were added in (3). In particular, this step preserves smoothness which is ensured for $D_P$ by Lemma 11. The size bound on $D$ follows from Lemma 12 using the fact that $|C| \le m$ and $|\text{rel}(c)| \le d^2$ for every $c \in C$. ◀
# 4 Compiling an SDNNF to a BCT Constraint

In this section, we shall show the following theorem.

**► Theorem 14.** Let $c^*$ be a constraint represented by a smooth SDNNF. Then there is a tree binary encoding of $c^*$ of polynomial size.

It is useful to look at a DNNF in terms of certificates [4], also called minimal satisfying subtrees in [16]. A certificate for a satisfying assignment is simply a minimal satisfied sub-DNNF that contains the output gate. Due to decomposability, the certificates of a DNNF are trees whose leaves are some of the leaves of the DNNF. In addition, no two leaves of a certificate are labeled with a literal of the same variable. Assume $D$ is a smooth DNNF on variables $\mathbf{x} = (x_1, \dots, x_n)$ and $S$ its certificate. We shall also assume that $D$ has no leaves labeled with constants 0 or 1 (as these can always be propagated). Then for each $i = 1, \dots, n$, we have that $S$ contains exactly one leaf associated with a literal of variable $x_i$ (see also [16] for more details). The leaves of $S$ thus determine a tuple $\tau$ on which $D$ evaluates to true. We say in this case that the leaves of $S$ are associated with the literals in tuple $\tau$. The certificates are thus in one-to-one correspondence with the satisfying assignments of $D$.

Let us now fix a constraint $c^*$ with $scp(c^*) = \mathbf{x}$. Let us assume that $c^*$ is represented by a smooth SDNNF $D$ that respects a v-tree $T$ and let $\rho$ denote the root of $D$. In particular, $var(T) = var(\rho) = \mathbf{x}$. Then the certificates of $D$ also respect $T$. In fact, we will show that if $S$ is a certificate of $D$, then the conjunction gates of $S$ are in one-to-one correspondence with the inner nodes of $T$. This property lies at the basis of our construction of a tree binary encoding $P = (\mathbf{z}, C)$ of $c^*$. The idea is to introduce a hidden variable for each inner node $t$ of $T$ with the domain being the $\wedge$-gates whose d-node is $t$. The constraints make sure that the models of $P$ are in one-to-one correspondence with the certificates of $D$.

We will construct a BCT $P = (\mathbf{z}, C)$ satisfying $\mathbf{x} \subseteq \mathbf{z}$ and $rel(c^*) = sol(P)[\mathbf{x}]$. For each inner node $t$ of $T$, we introduce a hidden variable $y_t$. The set of all these hidden variables will be denoted as $\mathbf{y}$. We then define $\mathbf{z} = \mathbf{x} \cup \mathbf{y}$. The domain of an original variable $x_i \in \mathbf{x}$ is $dom(x_i)$ as given by the constraint $c^*$. For a hidden variable $y_t \in \mathbf{y}$, we set $dom(y_t) = \Lambda(t)$ where $\Lambda(t)$ denotes the set of conjunction gates of $D$ that have d-node $t$.

The constraints of $C$ correspond to the edges of $T$. In particular, for each edge $(t, t')$ of $T$ where $t'$ is a child node of $t$, we add a constraint $c_{t,t'}$ to $C$ whose definition differs slightly depending on whether $t'$ is a leaf or an inner node of $T$. A sequence of vertices $v_0, \dots, v_k$ of $D$ is called an $\vee$-path if it is a path, nodes $v_1, \dots, v_{k-1}$ are $\vee$-gates, $v_0$ is a conjunction gate and $v_k$ is either a conjunction gate, or a leaf node.

(C1) Assume $t'$ is a leaf labeled with variable $x_i$. Then $scp(c_{t,t'}) = \{y_t, x_i\}$ and $rel(c_{t,t'})$ is defined as a set of tuples $\{(y_t, v), (x_i, a)\}$ such that $D$ contains a $\vee$-path from $v$ to the leaf labeled with literal $(x_i, a)$.

(C2) Assume $t'$ is an inner node of $T$. Then $scp(c_{t,t'}) = \{y_t, y_{t'}\}$ and $rel(c_{t,t'})$ is defined as a set of tuples $\{(y_t, v), (y_{t'}, v')\}$ such that $D$ contains a $\vee$-path from $v$ to $v'$.

**► Example 15.** Figure 4 shows the result of the application of the construction to the SDNNF from Figure 2. Recall that the SDNNF itself was constructed as a representation of the BCT $P$ from Example 3. BCT in Figure 4 differs from $P$ in that it has three hidden variables $y_1, y_2$, and $y_3$. Note that $y_1$ and $y_2$ are basically equivalent to $z_3$ and $y_3$ is equivalent to $z_4$. The auxiliary variables can thus be easily eliminated by which we obtain the constraints of $P$.
<figure><img src="image_211.png" /><figcaption>Figure 4 Example of the construction of the tree binary encoding of a constraint represented by the SDNNF $D$ in Figure 2.</figcaption></figure>

The size of $P$ defined in this way is clearly polynomial in the size of $D$. The rest of this section is devoted to showing the correctness of the construction. We will start with a few technical propositions on the structure of the certificates of $D$.

**Lemma 16.** Assume $v$ is a node of $D$ with d-node $t$ in $T$. Then $var(v) = var(t)$. Moreover, 
1. if $v = v_1 \wedge v_2$, then $t^l$ is the d-node of $v_1$ and $t^r$ is the d-node of $v_2$, and 
2. if $v = v_1 \vee \cdots \vee v_k$, then $t$ is the d-node of all input nodes $v_1, \ldots, v_k$.

**Proof.** We will proceed by the induction on the structure of $D$. If $v = \rho$ is the root of $D$, then $var(v) = \mathbf{x}$. It follows that its d-node $t$ is the root of $T$ and thus $var(t) = \mathbf{x} = var(v)$.

Let us assume that $v = v_1 \wedge v_2$ and that $var(v) = var(t)$. Since $D$ does not contain leaves labeled with constants, we have that both $var(v_1)$ and $var(v_2)$ are nonempty and thus $var(v_i) \subseteq var(v)$ for $i = 1, 2$ and $var(v) = var(v_1) \cup var(v_2)$. Let $t_i$ be the d-node of $v_i$, $i = 1, 2$. By the definition of structured DNNFs, both $t_1$ and $t_2$ are descendants of $t$ in $T$ and thus $var(t_i) \subseteq var(t)$, $i = 1, 2$. By the definition of d-nodes, we also have that $var(v_i) \subseteq var(t_i)$, $i = 1, 2$. It follows that $var(t) = var(v) = var(v_1) \cup var(v_2) \subseteq var(t_1) \cup var(t_2) \subseteq var(t)$ and thus $var(t) = var(t_1) \cup var(t_2)$. The only possibility is that both $t_1$ and $t_2$ are the child nodes of $t$ and thus $t_1 = t^l$, $t_2 = t^r$, and $var(v_i) = var(t_i)$, $i = 1, 2$.

Assume that $v = v_1 \vee \cdots \vee v_k$ and $var(v) = var(t)$. By smoothness we get that $var(t) = var(v) = var(v_1) = \cdots = var(v_k)$. It follows that $t$ is the d-node of all input nodes $v_1, \ldots, v_k$.

We have shown that if $var(v) = var(t)$ and $v'$ is a child node of $v$ with d-node $t'$, then $var(v') = var(t')$ which also holds for the leaves. $\quad$

**Lemma 17.** Assume that $v_0, \ldots, v_k$ is a $\vee$-path with $k > 0$. Assume that $t$ is the d-node of $v_0$ and $t'$ is the d-node of $v_k$. Then $v_1, \ldots, v_{k-1}$ have d-node $t'$ and $t'$ is a child node of $t$.

**Proof.** By smoothness, $var(v_1) = var(v_2) = \cdots = var(v_k)$ and thus all gates $v_1, \ldots, v_{k-1}$ have the same d-node as $v_k$. In particular, $t'$ is the d-node of $v_1$ which is an input to the conjunction gate $v_0$. By Lemma 16 we have that $t'$ is a child node of $t$. $\quad$

Based on Lemma 17, we can show the following proposition.

**Lemma 18.** Assume $S$ is a certificate and $t$ is an inner node of $T$. Then $S$ contains exactly one conjunction gate $v$ from $\Lambda(t)$.
**Proof.** Consider a variable $x_i \in \text{var}(t)$ and the path $v_0, \dots, v_k$ in $S$ that leads from the root $v_0 = \rho$ to the leaf $v_k$ of $S$ labeled with a literal on $x_i$. It follows that $\text{var}(v_j) \subseteq \text{var}(v_{j-1})$ for every $j = 1, \dots, k$. Moreover $\text{var}(v_0) = \mathbf{x}$ and $\text{var}(v_k) = \{x_i\}$. Let $v_{i_1}, \dots, v_{i_p}$ be the subsequence of $v_0, \dots, v_k$ formed only by conjunction nodes. Then by Lemma 17 we get that $t_{i_1}, \dots, t_{i_p}, t_k$ form a path in $T$ from the root to the leaf labeled with $x_i$. For some index $i_j$ we thus have that $t_{i_j} = t$ and it follows that $v_{i_j}$ is a conjunction gate in $S$ with d-node $t$.

Let us now assume that $S$ contains two $\wedge$-gates $v_1$ and $v_2$ with the same d-node $t$, thus $\text{var}(v_1) = \text{var}(v_2) = \text{var}(t)$ by Lemma 16. However, Lemma 16 also implies that there is no path from $v_1$ to $v_2$ or from $v_2$ to $v_1$. If we take the paths from the root $\rho$ to $v_1$ and $v_2$ in $S$, they have to split in a $\wedge$-gate $v$ (by minimality of $S$), but then $v$ is not decomposable.

It follows that $v_{i_j}$ is the only conjunction gate in $S$ that belongs to $\Lambda(t)$. $\quad$

Note that each literal on a variable from $z = \mathbf{x} \cup \mathbf{y}$ is associated with a node of $D$. In particular, for $x_i \in \mathbf{x}$, literal $(x_i, a)$ is associated with the leaf of $D$ labeled with $(x_i, a)$. To this end, we need to assume that every such literal has a leaf labeled with it. However, if $D$ does not contain any leaf associated with literal $(x_i, a)$, then this literal does not have a support in $c^*$ and $a$ can be removed from $\text{dom}(x_i)$. We may thus assume that no such value is in $\text{dom}(x_i)$. For an inner node $t$ of $T$, a literal $(y_t, v)$ is associated with the node $v \in \Lambda(t)$.

**Lemma 19.** Let $\tau \in \text{sol}(P)$ be a tuple that is a solution to $P$. Then $\tau[\mathbf{x}] \in \text{rel}(c^*)$.

**Proof.** Since $D$ represents $c^*$, it is enough to show that there is a certificate $S$ of $D$ whose leaves are associated with the literals in $\tau[\mathbf{x}]$.

Tuple $\tau$ associates a node $v$ of $D$ with every node $t$ of $T$. We proceed by induction on the structure of $T$ to describe a certificate $S_t$ for the sub-DNNF of $D$ rooted at $v$.

Assume first that $t$ is a leaf of $T$ labeled with variable $x_i$. Consider the literal $(x_i, a) \in \tau$ and set the certificate $S_t$ to a single node labeled with this literal.

Assume now that $t$ is an inner node of $T$. Since $t$ is an inner node of $T$, we have that $(y_t, v) \in \tau$ for some $v \in \Lambda(t)$. Tuple $\tau$ also contains literals associated with $t^l$ and $t^r$. These literals associate a nodes $v_l$ and $v_r$ of $D$ with $t^l$ and $t^r$ respectively. By induction hypothesis, we have constructed certificate $S_l$ and $S_r$ for the sub-DNNFs rooted at $v_l$ and $v_r$ respectively. Since $\tau$ satisfies constraints $c_{t,t^l}$ and $c_{t,t^r}$, $D$ contains a $\vee$-paths from $v$ to $v_l$ and from $v$ to $v^r$. Certificate $S_t$ for the sub-DNNF rooted at $v$ is then constructed as a union of $S_l, S_r$, node $v$ and the $\vee$-paths from $v$ to $v_l$ and $v_r$.

Let $\sigma$ be the root of $T$ and let us assume that $v$ is the node of $D$ associated with $\sigma$ by $\tau$. Let $S_\sigma$ be the certificate of the sub-DNNF rooted at $v$. If $v = \rho$ is the root of $D$, then $S = S_\sigma$ is a certificate of $D$. Otherwise, $D$ contains a path from $\rho$ to $v$ that consists only of $\vee$-gates and we construct $S$ by combining this path with $S_\sigma$. $\quad$

**Lemma 20.** For every $\tau^* \in \text{rel}(c^*)$, there is $\tau \in \text{sol}(P)$ satisfying $\tau^* = \tau[\mathbf{x}]$.

**Proof.** Since $\tau^* \in \text{rel}(c^*)$, there is a certificate $S$ of $D$ whose leaves are associated with the literals from $\tau^*$. By Lemma 18, the certificate $S$ contains exactly one conjunction gate $v_t \in \Lambda(t)$ for each inner node $t$. We form $\tau$ by adding literals $(y_t, v_t)$ to $\tau^*$ for all internal nodes $t$ of $T$. Let us check that $\tau$ satisfies all constraints of $P$. Let $c_{t,t'}$ be a constraint of $P$ where $t'$ is a child node of $t$ in $T$. By Lemma 16, one of the child nodes of $v_t$ in $D$ has d-node $t'$, let us denote it $v_1$. Since $v_t$ is a conjunction gate, $v_1$ must belong to $S$. If $v_1$ is a disjunction, then by Lemma 16, its child nodes have d-node $t'$, too. If we follow the path in $S$ from $v_1$ to a leaf or to the next conjunction gate, we get a $\vee$-path that ends in the node $v_k$ whose d-node is still $t'$ and $v_k$ is either a leaf or a conjunction gate.
If $t'$ is a leaf of $T$ labeled with variable $x_i$, we must have that $v_k$ is a leaf of $S$ labeled with a literal $(x_i, a)$ for some $a \in \text{dom}(x_i)$, it follows that $(x_i, a) \in \tau^* \subseteq \tau$ and $\{(y_t, v_t), (x_i, a)\} \in \text{rel}(c_{t,t'})$, constraint $c_{t,t'}$ is thus satisfied by $\tau$.

If $t'$ is an inner node, then $v_k$ is a conjunction gate $v_{t'}$ associated with $t'$ in $S$. It follows that $(y_{t'}, v_{t'}) \in \tau$ and $\{(y_t, v_t), (y_{t'}, v_{t'})\} \in \text{rel}(c_{t,t'})$, constraint $c_{t,t'}$ is thus satisfied by $\tau$. $\quad$

Theorem 14 now follows by the above construction from the following proposition.

**Theorem 21.** *P* = (z, C) is a tree binary encoding of c* of polynomial size.

*Proof.* C consists of $O(n)$ constraints and the total size of the domains of variables in z is bounded by the number of the nodes in D. Lemmas 19 and 20 imply that P is a TBE of $c^*$. $\quad$

# 5 Binary Constraint Graphs With Bounded Treewidth

In this section, we shall extend the construction from Section 3 to BCG constraints that naturally generalize BCT constraints.

**Definition 22.** A BCG constraint $c$ is a pair $(\mathbf{x}, P)$ such that $P = (\mathbf{z}, C)$ is a normalized binary CSP, $\text{scp}(c) = \mathbf{x} \subseteq \mathbf{z}$ and $\text{rel}(c) = \text{sol}(\mathbf{z}, C)[\mathbf{x}]$.

The construction we describe is parameterized by the treewidth of the underlying constraint graph and the domain size. The treewidth of a graph is defined using a tree decomposition.

Given an undirected graph $G = (V, E)$, its tree decomposition is defined as a pair $(T, \chi)$ where $T = (V_T, E_T)$ is a tree and $\chi: V_T \rightarrow \mathcal{P}(V)$ is a function that assigns each vertex $t \in V_T$ a subset of $V$ called a *bag* that satisfies the following conditions:

(d1) $V = \bigcup_{t \in V_T} \chi(t)$.

(d2) For each edge $\{u, v\} \in E$ there is a node $t \in V_T$ such that $\{u, v\} \subseteq \chi(t)$.

(d3) If a node $v$ is contained in two bags $\chi(t_1)$ and $\chi(t_2)$, then $v \in \chi(t)$ for every node $t$ on the path connecting $t_1$ with $t_2$.

The width of the tree decomposition is defined as $\max_{t \in V_T} |\chi(t)| - 1$. The treewidth $\text{tw}(G)$ of $G$ is the minimum width among all possible tree decompositions of $G$. It should be noted that any tree decomposition of a graph $G$ on $n$ vertices can be transformed into a tree decomposition with the same width and $O(n)$ nodes [15].

We are now ready to formulate the main result of this section.

**Theorem 23.** Assume that $c^* = (\mathbf{x}, P)$ is a BCG constraint defined by a normalized binary CSP $P = (\mathbf{z}, C)$. Denote $G$ the constraint graph of $P$. Denote $m = |\mathbf{z}|$ and $d = \max_{z_i \in \mathbf{z}} |\text{dom}(z_i)|$. Then there is an SDNNF $D$ representing $c^*$ with $O(m^{\text{tw}(G)+1})$ nodes and $O(m^{2^{\text{tw}(G)+1}})$ edges.

The proof of Theorem 23 is based on the following proposition.

**Theorem 24.** Assume that $c^* = (\mathbf{x}, P)$ is a BCG constraint defined by a normalized binary CSP $P = (\mathbf{z}, C)$. Denote $G$ the constraint graph of $P$. Denote $m = |\mathbf{z}|$ and $d = \max_{z_i \in \mathbf{z}} |\text{dom}(z_i)|$. Then $c^*$ has a tree binary encoding $P' = (\mathbf{z}', C')$ with $|\mathbf{z}'| = O(m)$ and $|\text{dom}(z_i')| \le d^{\text{tw}(G)+1}$ for every $z_i' \in \mathbf{z}'$.

Note that Theorem 24 actually follows from Proposition 4 in [28] which is based on the encoding described in [11]. If we would simply combine the bound given by Theorem 24 with the bound given by Theorem 9, we would get an SDNNF $D$ representing $c^*$ with
<figure><img src="image_212.png" /><figcaption>**Figure 5** Example of the construction, see the description in Example 25.</figcaption></figure>

$O(md^{tw(G)+1})$ nodes and $O(md^{2(tw(G)+1)})$ edges. We provide a specific construction that proves Theorem 24 and that can be combined with Lemma 12 to prove a slightly better bound stated in Theorem 23. The construction we describe below is similar to the dual encoding described in [11].

Let us fix a BCG constraint $c^* = (\mathbf{x}, P)$ where $P = (\mathbf{z}, C)$ is a normalized binary CSP with constraint graph $G$. Let us assume that $\mathbf{z} = (z_1, \dots, z_m)$ and $d = \max_{i=1}^m |\text{dom}(z_i)|$. Let us also fix a tree decomposition $(T, \chi)$ of $G$. Let us assume that $V_T = (t_1, \dots, t_N)$ for some $N = O(m)$. We shall describe a BCT $P' = (z', C')$ which is a TBE of $c^*$. First, let us define the variables in $z'$. We associate a new variable $v_i$ with every $t_i$, $i = 1, \dots, N$. Then we set $z' = \mathbf{z} \cup \mathbf{v}$ where $\mathbf{v} = (v_1, \dots, v_N)$. The domains of variables in $\mathbf{z}$ are given by $c^*$. For every $v_i$, $i = 1, \dots, N$, we set the domain as follows. Let us consider the set of constraints defined on variables from $\chi(t_i)$ as $C_i = \{c \in C \mid \text{scp}(c) \subseteq \chi(t_i)\}$. Then the domain of $v_i$ is defined as the set of solutions to CSP $(\chi(t_i), C_i)$, i.e. $\text{dom}(v_i) = \text{sol}(\chi(t_i), C_i)$.

Let us now define the constraints in $C'$.

(T1) For every edge $\{t_i, t_j\} \in E_T$ we add a constraint $c'_{i,j}$ into $C'$ with $\text{scp}(c'_{i,j}) = \{v_i, v_j\}$. The constraint relation $\text{rel}(c'_{i,j})$ consists of pairs $\{(v_i, \tau_1), (v_j, \tau_2)\}$ where $\tau_1 \in \text{dom}(v_i)$, $\tau_2 \in \text{dom}(v_j)$, and $\tau_1[\chi(t_i) \cap \chi(t_j)] = \tau_2[\chi(t_i) \cap \chi(t_j)]$.

(T2) For every $z_i$, $i = 1, \dots, m$, we pick a representative node $t_{r_i} \in V_T$ satisfying $z_i \in \chi(t_{r_i})$. We then add a constraint $c'_i$ into $C'$ with $\text{scp}(c'_i) = \{z_i, t_{r_i}\}$. The set of tuples $\text{rel}(c'_i)$ consists of pairs $\{(z_i, a), (v_{r_i}, \tau)\}$ where $a \in \text{dom}(z_i)$, $\tau \in \text{dom}(v_{r_i})$, and $(z_i, a) \in \tau$.

**Example 25.** Let us consider a binary CSP $P = (\mathbf{z}, C)$ with $\mathbf{z} = \{z_1, \dots, z_5\}$ whose constraint graph $G$ is depicted in Figure 5a. We shall use $c_{i,j} \in C$ to denote the constraint with scope $\{z_i, z_j\}$. Figure 5b shows a tree decomposition $T$ of the graph with the contents of the bags inside the rectangles. The structure of the tree binary encoding $P'$ of $P$ is then shown in Figure 5c. The domain of variable $v_1$ consists of tuples $\tau$ on variables $z_1, z_2$, and $z_3$ satisfying constraints $c_{1,2}, c_{2,3}$, and $c_{1,3}$. Assume a tuple $\sigma' \in \text{sol}(P')$. Constraint $c'_2$ makes sure that if $(z_2, a) \in \sigma'$, then $\sigma'$ contains $(v_2, \tau)$ satisfying $(z_2, a) \in \tau$. Similarly for other variables. Constraints $c'_{1,2}$ and $c'_{2,3}$ extend this property also to nodes $v_2$ and $v_3$. The tuples assigned to variables $v_1, v_2$, and $v_3$ are thus consistent with each other and also with constraints $C$. We thus have that $\sigma'[\mathbf{z}] \in \text{sol}(P)$.
The proof of the correctness of the above construction and thus also the proof of Theo-
rem 24 is moved to the appendix. Here, we will describe its application for proving the main
result of this section.

**Proof of Theorem 23.** Using Theorem 24, we obtain a TBE encoding $P' = (z', C')$ with $O(m)$ variables and the domain sizes bounded by $d^{tw(G)+1}$. We can then apply Theorem 9 to obtain an SDNNF $D$ that represents $c^*$. $D$ has $O(md^{tw(G)+1})$ nodes. By Lemma 12, the number of edges of $D$ is bounded by $O(md^{tw(G)+1} + s)$ where $s = \sum_{c' \in C'} |rel(c')|$. Since $|C'| = O(m)$, it is enough to show that $|rel(c')| \le d^{2 \tw(G)+1}$ for every $c' \in C'$.

Assume first a constraint $c'_{i,j}$ added in step (T1). We may assume that $G$ is connected (otherwise we process each connected component of $G$ separately) and therefore $\chi(t_i) \cap \chi(t_j) \neq \emptyset$. The number of pairs of tuples $\tau_1$ and $\tau_2$ that satisfy $\tau_1[\chi(t_i) \cap \chi(t_j)] = \tau_2[\chi(t_i) \cap \chi(t_j)]$ is thus at most $d^{tw(G)+1} \cdot d^{tw(G)} = d^{2 \tw(G)+1}$. Therefore $|rel(c'_{i,j})| \le d^{2 \tw(G)+1}$.

Assume now a constraint $c'_i$ added in step (T2). The number of tuples $\tau$ satisfying that $(z_i, a) \in \tau$ for one particular $a \in \text{dom}(z_i)$ is at most $d^{\text{tw}(G)}$ and thus $|\text{rel}(c'_{i,j})| \le d^{\text{tw}(G)+1} \le d^{2 \text{tw}(G)+1}$. ◇

Note that the size estimate in Theorem 23 is only an upper bound and the real size of $P'$
and the SDNNF $D$ depends on the particular tree decomposition and, in particular, on how
much the bags intersect. Therefore, there is a space for optimization in a practical setting.

6 Conclusion

As the main result of our paper, we have shown that binary constraint trees are polynomially
equivalent to structured DNNF circuits. We would like to note that for a given BCT *P*
the construction in Section 3 leads to a deterministic SDNNF *D*<sub>*P*</sub> (thanks to rule 3 in
step (A2b)). This means that for every pair of distinct children *v*<sub>1</sub> and *v*<sub>2</sub> of a disjunction
node, the sub-NNFs rooted at *v*<sub>1</sub> and *v*<sub>2</sub> do not share any models (see [9, 20]). This property
allows for instance model counting on *D*<sub>*P*</sub>. However, forgetting the hidden variables from *D*<sub>*P*</sub>
does not preserve determinism in general [10] and thus the actual result of the construction
is not a deterministic SDNNF. Introducing hidden variables is thus an important part
of the construction described in Section 4 since SDNNFs are strictly more succinct than
deterministic SDNNFs [20].

Several rules for reducing the number of hidden variables in a BCT constraint were
described in [27], it would be interesting to investigate the effect of these rules on a SDNNF
that is compiled into a BCT constraint, reduction rules are applied to it and then it is
compiled back to a SDNNF. When compiling the BCT constraint back to a SDNNF, we
can pick an arbitrary node of the constraint tree as a root which allows us to change the
structure of the SDNNF to a different orientation of the v-tree. This, for instance, extends
the applicability of the conjoin operation described in [20] to conjoining two SDNNFs whose
v-trees differ, but their undirected versions are the same.

References

1. Ignasi Abío, Graeme Gange, Valentin Mayer-Eichberger, and Peter J. Stuckey. On CNF encodings of decision diagrams. In Claude-Guy Quimper, editor, *Integration of AI and OR Techniques in Constraint Programming*, pages 1–17, Cham, 2016. Springer International Publishing.

2 Christian Bessiere and Jean-Charles Régin. Arc consistency for general constraint networks: preliminary results. In *IJCAI (1)*, pages 398–404, 1997.
3 Simone Bova, Florent Capelli, Stefan Mengel, and Friedrich Slivovsky. On compiling CNFs into structured deterministic DNNFs. In *Marijn Heule and Sean Weaver, editors, Theory and Applications of Satisfiability Testing - SAT 2015*, pages 199–214, Cham, 2015. Springer International Publishing.

4 Simone Bova, Florent Capelli, Stefan Mengel, and Friedrich Slivovsky. Knowledge compilation meets communication complexity. In *Proceedings of the Twenty-Fifth International Joint Conference on Artificial Intelligence, IJCAI’16*, pages 1008–1014. AAAI Press, 2016. URL: http://dl.acm.org/citation.cfm?id=3060621.3060761.

5 Kenil C. Cheng and Roland H. Yap. Maintaining generalized arc consistency on ad hoc r-ary constraints. In *Proceedings of the 14th International Conference on Principles and Practice of Constraint Programming*, CP ’08, pages 509–523, Berlin, Heidelberg, 2008. Springer-Verlag. doi:10.1007/978-3-540-85958-1_34.

6 Kenil C. Cheng and Roland H. Yap. An mdd-based generalized arc consistency algorithm for positive and negative table constraints and some global constraints. *Constraints*, 15(2):265–304, April 2010. doi:10.1007/s10601-009-9087-y.

7 Adnan Darwiche. Compiling knowledge into decomposable negation normal form. In *Proceedings of the 16th International Joint Conference on Artificial Intelligence - Volume 1, IJCAI’99*, pages 284–289, San Francisco, CA, USA, 1999. Morgan Kaufmann Publishers Inc. URL: http://dl.acm.org/citation.cfm?id=1624218.1624260.

8 Adnan Darwiche. Decomposable negation normal form. *J. ACM*, 48(4):608–647, July 2001. doi:10.1145/502090.502091.

9 Adnan Darwiche. A compiler for deterministic, decomposable negation normal form. In *Eighteenth National Conference on Artificial Intelligence*, pages 627–634, Menlo Park, CA, USA, 2002. American Association for Artificial Intelligence. URL: http://dl.acm.org/citation.cfm?id=777092.777189.

10 Adnan Darwiche and Pierre Marquis. A knowledge compilation map. *Journal of Artificial Intelligence Research*, 17:229–264, 2002.

11 Rina Dechter and Judea Pearl. Tree clustering for constraint networks. *Artificial Intelligence*, 38(3):353–366, April 1989. doi:10.1016/0004-3702(89)90037-4.

12 Eugene C. Freuder. A sufficient condition for backtrack-bounded search. *Journal of the ACM*, 32(4):755–761, October 1985. doi:10.1145/4221.4225.

13 Eugene C. Freuder. Complexity of K-tree structured constraint satisfaction problems. In *Proceedings of the eighth National conference on Artificial intelligence - Volume 1*, AAAI’90, pages 4–9, Boston, Massachusetts, 1990. AAAI Press.

14 Graeme Gange and Peter J. Stuckey. Explaining propagators for s-DNNF circuits. In Nicolas Beldiceanu, Narendra Jussien, and Éric Pinson, editors, *Integration of AI and OR Techniques in Constraint Programming for Combinatorial Optimization Problems*, pages 195–210. Springer Berlin Heidelberg, 2012.

15 Ton Kloks. *Treewidth, Computations and Approximations*, volume 842 of *Lecture Notes in Computer Science*. Springer, 1994. doi:10.1007/BFb0045375.

16 Petr Kučera and Petr Savický. Propagation complete encodings of smooth DNNF theories. *Constraints*, 27(3):327–359, July 2022. doi:10.1007/s10601-022-09331-2.

17 Christophe Lecoutre. STR2: optimized simple tabular reduction for table constraints. *Constraints*, 16(4):341–371, October 2011. doi:10.1007/s10601-011-9107-6.

18 Christophe Lecoutre, Chavalit Likitvivatanavong, and Roland H. C. Yap. STR3: A path-optimal filtering algorithm for table constraints. *Artificial Intelligence*, 220:1–27, March 2015. doi:10.1016/j.artint.2014.12.002.

19 Robert Mateescu and Rina Dechter. Compiling Constraint Networks into AND/OR Multi-valued Decision Diagrams (AOMDDs). In Frédéric Benhamou, editor, *Principles and Practice of Constraint Programming - CP 2006*, Lecture Notes in Computer Science, pages 329–343, Berlin, Heidelberg, 2006. Springer. doi:10.1007/11889205_25.
22:18 Binary Constraint Trees and Structured Decomposability

20 Knot Pipatsrisawat and Adnan Darwiche. New compilation languages based on structured decomposability. In *Proceedings of the 23rd National Conference on Artificial Intelligence - Volume 1*, AAAI'08, pages 517–522. AAAI Press, 2008. URL: http://dl.acm.org/citation.cfm?id=1619995.1620079.

21 Knot Pipatsrisawat and Adnan Darwiche. Top-down algorithms for constructing structured DNNF: Theoretical and practical implications. In *Proceedings of the 2010 Conference on ECAI 2010: 19th European Conference on Artificial Intelligence*, pages 3–8, Amsterdam, The Netherlands, The Netherlands, 2010. IOS Press. URL: http://dl.acm.org/citation.cfm?id=1860967.1860970.

22 Francesca Rossi, Charles J. Petrie, and Vasant Dhar. On the equivalence of constraint satisfaction problems. In *Proceedings of the 9th European Conference on Artificial Intelligence*, ECAI'90, pages 550–556, USA, 1990. Pitman Publishing, Inc.

23 Andy Shih, Guy Van den Broeck, Paul Beame, and Antoine Amarilli. Smoothing structured decomposable circuits. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, *Advances in Neural Information Processing Systems*, volume 32. Curran Associates, Inc., 2019. URL: https://proceedings.neurips.cc/paper/2019/file/940392f5f32a7ade1cc201767cf83e31-Paper.pdf.

24 Kostas Stergiou and Toby Walsh. Encodings of non-binary constraint satisfaction problems. In *Proceedings of the Sixteenth National Conference on Artificial Intelligence and the Eleventh Innovative Applications of Artificial Intelligence Conference Innovative Applications of Artificial Intelligence*, AAAI '99/IAAI '99, pages 163–168, USA, 1999. American Association for Artificial Intelligence.

25 Ruiwei Wang and Roland H. C. Yap. Bipartite encoding: A new binary encoding for solving non-binary CSPs. In *Proceedings of the Twenty-Ninth International Joint Conference on Artificial Intelligence*, IJCAI'20, Yokohama, Yokohama, Japan, 2021.

26 Ruiwei Wang and Roland H. C. Yap. CNF Encodings of Binary Constraint Trees. In Christine Solnon, editor, *28th International Conference on Principles and Practice of Constraint Programming (CP 2022)*, volume 235 of *Leibniz International Proceedings in Informatics (LIPIcs)*, pages 40:1–40:19, Dagstuhl, Germany, 2022. Schloss Dagstuhl – Leibniz-Zentrum für Informatik. doi:10.4230/LIPIcs.CP.2022.40.

27 Ruiwei Wang and Roland H.C. Yap. Encoding multi-valued decision diagram constraints as binary constraint trees. *Proceedings of the AAAI Conference on Artificial Intelligence*, 36(4):3850–3858, June 2022. doi:10.1609/aaai.v36i4.20300.

28 Ruiwei Wang and Roland H.C. Yap. The expressive power of ad-hoc constraints for modelling CSPs. *Proceedings of the AAAI Conference on Artificial Intelligence*, 37(4):4104–4114, June 2023. doi:10.1609/aaai.v37i4.25526.

A Proof of Theorem 24

In this section, we shall prove the correctness of the construction described in Section 5. We shall also prove Theorem 24 that states the properties of the construction. We use the same notation that was used in Section 5. In particular, we assume a fixed BCG constraint $c^* = (\mathbf{x}, P)$ where $P = (\mathbf{z}, C)$ is a normalized binary CSP with constraint graph $G$. We assume that that BCT $P' = (\mathbf{z}', C')$ is the result of the construction from Section 5. We shall show that $P'$ is a TBE of $c^*$. The construction of $C'$ implies the following property.

**Lemma 26.** Let $p \in \{1, \dots, m\}$ be arbitrary and let $t_i \in V_T$ be such that $z_p \in \chi(t)$. Assume that $\sigma' \in \text{sol}(P')$ and assume that $(z_p, a), (v_i, \tau) \in \sigma'$. Then $(z_p, a) \in \tau$.

**Proof.** Let $t_{r_p}$ be the representative node picked for $z_p$ in step (T2) and consider the path $t_{r_p} = t_{j_1}, t_{j_2}, \dots, t_{j_k} = t_i$ in $T$ connecting $t_{r_p}$ with $t_i$. For every $q = 1, \dots, k$ we have that $z_p \in \chi(t_{r_p}) \cap \chi(t_i)$ and thus $z_p \in \chi(t_{j_q})$ by condition (d3). Denote $\tau_q \in \text{dom}(v_{j_q})$ the tuple for which $(v_{j_q}, \tau_q) \in \sigma'$. We shall show by induction on $q$ that $(z_p, a) \in \tau_q$ for every $q = 1, \dots, k$. Since $\tau_k = \tau$, we then have that $(z_p, a) \in \tau$.
Since $t_1 = t_{r_p}$, we have that $(z_p, a) \in \tau_1$ by constraint $c'_p$ added to $C'$ in step (T2).
Assume now that $q > 1$ and consider constraint $c'_{j_{q-1},j_q}$ added in step (T1). The induction
hypothesis implies that $(z_p, a) \in \tau_{q-1}$, and by definition of $\mathrm{rel}(c'_{j_{q-1},j_q})$ we also have that
$(z_p, a) \in \tau_q$.
$\quad$

We are now ready to prove the correctness of the construction.

**Lemma 27.** (*x*, *P*') is a tree binary encoding of BCG constraint *c*\* = (*x*, *P*).

**Proof.** The constraint graph of $P'$ is a tree that originates from $T$ by adding leaves corresponding to the constraints $c'_i$ added in step (T2). We shall show that $\mathrm{sol}(P')[\mathbf{z}] = \mathrm{sol}(P)$. Then $\mathrm{rel}(c^*) = \mathrm{sol}(\mathbf{z}, C)[\mathbf{x}] = \mathrm{sol}(\mathbf{z}', C')[\mathbf{x}]$ and the proposition follows.

Assume first that we have a solution $\sigma' \in \mathrm{sol}(P')$. Denote $\sigma = \sigma'[\mathbf{z}]$ and let us show
that $\sigma$ satisfies all constraints of $P$. Let $c \in C$ be a constraint with $\mathrm{scp}(c) = \{z_p, z_q\}$. We
have $(z_p, a) \in \sigma$ and $(z_q, b) \in \sigma$ for some $a \in \mathrm{dom}(z_p)$ and $b \in \mathrm{dom}(z_q)$. By condition (d2)
we have that $\mathrm{scp}(c) \subseteq \chi(t_i)$ for some $t_i \in V_T$. Consider literal $(v_i, \tau) \in \sigma'$. By Lemma 26,
we have that $(z_p, a) \in \tau$ and $(z_q, b) \in \tau$. Since $\tau \in \mathrm{dom}(v_i)$, we have that $\{(z_p, a), (z_q, b)\} =
\tau[\mathrm{scp}(c)] \in \mathrm{rel}(c)$. Since this holds for every constraint $c \in C$, we get that $\sigma \in \mathrm{sol}(P)$.

Assume now that we have a solution $\sigma \in \mathrm{sol}(P)$. Let us now define a tuple $\sigma' =
\sigma \cup \{(v_i, \sigma[\chi(t_i)]) \mid i = 1, \dots, N\}$. Since $\sigma \in \mathrm{sol}(P)$, we have that $\sigma[\chi(t_i)] \in \mathrm{rel}(C_i)$ and thus
$\sigma[\chi(t_i)] \in \mathrm{dom}(v_i)$. Tuple $\sigma'$ is thus correctly defined. It also satisfies all constraints (T1)
and (T2) and thus $\sigma' \in \mathrm{sol}(P')$.
$\quad$

**Proof of Theorem 24.** Assume that $P' = (\mathbf{z}', C')$ is constructed as above. Then $(\mathbf{x}, P')$ is a tree binary encoding of $c^* = (\mathbf{x}, P)$ by Lemma 27. We may assume by [15] that $|V_T| = O(m)$ and thus $|\mathbf{z}'| = m + |V_T| = O(m)$. For every variable $z_i \in \mathbf{z}$ we have $|\mathrm{dom}(z_i)| \le d$ by assumption. For every variable $v_i \in \mathbf{z}' \setminus \mathbf{z}$ we have that $|\chi(t_i)| \le \mathrm{tw}(G) + 1$ and thus $|\mathrm{dom}(v_i)| \le d^{\mathrm{tw}(G)+1}$ by the definition of $\mathrm{dom}(v_i)$.
$\quad$
<figure><img src="image_213.png" /></figure>
# Large Neighborhood Beam Search for Domain-Independent Dynamic Programming

Ryo Kuroiwa

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

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

J. Christopher Beck

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

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

## Abstract

Large neighborhood search (LNS) is an algorithmic framework that removes a part of a solution and performs search in the induced search space to find a better solution. While LNS shows strong performance in constraint programming, little work has combined LNS with state space search. We propose large neighborhood beam search (LNBS), a combination of LNS and state space search. Given a solution path, LNBS removes a partial path between two states and then performs beam search to find a better partial path. We apply LNBS to domain-independent dynamic programming (DIDP), a recently proposed generic framework for combinatorial optimization based on dynamic programming. We empirically show that LNBS finds better quality solutions than a state-of-the-art DIDP solver in five out of nine benchmark problem types with a total of 8570 problem instances. In particular, LNBS shows a significant improvement over the existing state-of-the-art DIDP solver in routing and scheduling problems.

**2012 ACM Subject Classification** Computing methodologies → Discrete space search

**Keywords and phrases** Large Neighborhood Search, Dynamic Programming, State Space Search, Combinatorial Optimization

Digital Object Identifier 10.4230/LIPIcs.CP.2023.23

Supplementary Material Software (Source Code): https://github.com/domain-independent-dp/didp-rs/releases/tag/lnbs-cp23

archived at swh:1:dir:b531b883c4c9eda4e0452fe801a5b50f6a3970cc

**Funding** This research is supported by the Natural Sciences and Engineering Research Council of Canada.

## 1 Introduction

In constraint programming (CP), large neighborhood search (LNS) [34] achieves strong performance in solving combinatorial optimization problems such as routing [24] and scheduling problems [27]. LNS is an algorithmic framework that removes a part of a solution and then performs search in the induced partial search space (neighborhood) to find a better solution. Typically, LNS uses tree search to find a better solution in a partial search space, where each search node represents a partial assignment of decision variables, and a solution of a problem corresponds to a leaf node, where all variables are assigned values.

Dynamic programming (DP) is a powerful method for multiple combinatorial optimization problems [14, 15], and the hybridization of CP, decision diagrams, and DP is a topic of active research [1, 22, 23, 5, 28, 19, 17]. Recently, domain-independent dynamic programming (DIDP), a model-based paradigm for combinatorial optimization based on DP, has been proposed [25]. In DIDP, a model of a problem is represented by a state transition system. A solution corresponds to a path in a state space graph, where each vertex represents a state and each edge represents a transition between two states. The current state-of-the-art DIDP

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

© Ryo Kuroiwa and J. Christopher Beck; 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. 23; pp. 23:1–23:22

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

LIPICS

Leibniz International Proceedings in Informatics

Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl Publishing, Germany
## 23:2 Large Neighborhood Beam Search for Domain-Independent Dynamic Programming

solver is complete anytime beam search (CABS) [26], which is based on beam search, an algorithm that searches for a path in a state space graph by maintaining a fixed number of states at a time.

In this paper, we propose large neighborhood beam search (LNBS), a combination of LNS and beam search in a state space graph. LNBS tries to improve a solution path by removing a partial path between two states and then performing beam search to find a better partial path. While LNBS has the freedom to select a neighborhood (i.e., a partial path to remove), we propose a strategy that dynamically adjusts the size of a neighborhood based on a multi-armed bandit problem. With our strategy, LNBS is complete, i.e., it finds and proves an optimal solution given enough time, but, of course, it is aimed at problems where its solution quality is more important than proved optimality. We implement LNBS for DIDP and empirically evaluate its performance. The experimental results show that LNBS outperforms CABS in five out of nine benchmark problem types in terms of solution quality. In addition, LNBS performs better than a commercial CP solver, which uses LNS [27], in seven problems while CABS is better than CP in six problems. Since LNBS performs particularly well in routing and scheduling problems, we also investigate the reason for this performance and gain insight from empirical analysis.

## 2 Background

We first introduce domain-independent dynamic programming and complete anytime beam search. Then, we present large neighborhood search (LNS). We also describe LNS with decision diagrams [18], a recently proposed method for combinatorial optimization that can be considered a combination of LNS and state space search.

### 2.1 Domain-Independent Dynamic Programming

A combinatorial optimization problem is to find a set of discrete decisions, e.g., a permutation, to minimize or maximize an objective function. In dynamic programming (DP), a problem is recursively formulated by decomposing it into subproblems, represented by *states*, and the optimal objective value of each subproblem is represented by the *value function*, which maps a state to a real number.

Domain-independent dynamic programming (DIDP) is a model-based paradigm for combinatorial optimization based on DP [25]. In DIDP, a DP formulation of a combinatorial optimization problem is defined by a state-transition system in Dynamic Programming Description Language (DyPDL). A DyPDL model is a seven-tuple $\langle V, S^0, K, T, B, C, h \rangle$ consisting of *state variables* $V$, the *target state* $S^0$, *constants* $K$, *transitions* $T$, *base cases* $B$, *state constraints* $C$, and the *dual bound* $h$.

A state variable $v \in V$ has a type of *set*, *element*, or *numeric*. Each set and element variable $v_i$ is associated with a set of objects $N_i = \{0, ..., n_i - 1\}$. The domain of a set variable $v_i$ is $2^{N_i}$, and the domain of an element variable $v_i$ is $N_i$. A numeric variable takes a real value. A constant in $K$ is a value independent of the state variables.

A state is a complete value assignment to the state variables, and the target state $S^0$ is a state. We denote the value of a state variable $v_i$ in a state $S$ by $S[v_i]$. The value of a state $S$ is represented by the value function $V(S)$, and the objective of a DyPDL model is to compute the value of the target state, $V(S^0)$. We can define dominance between states based on state variables. If a state $S$ dominates another state $S'$, denoted by $S' \preceq_S S$, then $V(S)$ is equal or better (less/greater for minimization/maximization) than $V(S')$. For the details of dominance in DyPDL, please refer to previous work [25].
A transition $\tau \in \mathcal{T}$ is a four-tuple $\langle \text{eff}_\tau, \text{cost}_\tau, \text{pre}_\tau, \text{forced}_\tau \rangle$ defining how a state $S$ is transformed to a successor state (a subproblem) $S[\llap{\tau}]$, and how $V(S)$ is computed based on the value of $V(S[\llap{\tau}])$. The set of effects $\text{eff}_\tau = \{(v_i, e_{v_i}) \mid v_i \in \mathcal{V}\}$ defines how each state variable $v_i$ is updated by an expression $e_{v_i}$. The expression $e_{v_i}$ consists of predefined operations on state variables and constants and returns $S[\llap{\tau}][v_i]$ given $S$. For example, for a set variable $U$, $U \setminus \{i\}$ is an expression representing removing an element $i$ from $U$, which can be used as $e_U$ and results in $S[\llap{\tau}][U] = S[U] \setminus \{i\}$. The cost expression $\text{cost}_\tau$ is an expression that takes $V(S[\llap{\tau}])$ in addition to $S$ and returns a real number $\text{cost}_\tau(V(S[\llap{\tau}]), S)$. The preconditions $\text{pre}_\tau$ are conditions on the state variables, i.e., expressions returning a binary value $\top$ or $\bot$. The preconditions define when the transition is applicable. For example, if $\text{pre}_\tau = \{i \in U\}$, then $\tau$ is applicable in $S$ iff $i \in S[U]$. The flag $\text{forced}_\tau$ is a boolean value indicating whether the transition is a forced transition. Let $\mathcal{T}(S)$ be the set of applicable transitions in state $S$. If forced transitions are applicable, the first defined one $\tau$ is selected, and all other forced and non-forced transitions are ignored, i.e., $\mathcal{T}(S) = \{\tau\}$. The value of $V(S)$ is computed by taking the best $\text{cost}_\tau(V(S[\llap{\tau}]), S)$ over all $\tau \in \mathcal{T}(S)$.

Base cases $\mathcal{B}$ are sets of conditions. For any $B \in \mathcal{B}$, if a state $S$ satisfies all conditions in $B$, denoted by $S \models B$, then it is called a base state, and $V(S)$ is defined non-recursively by an expression $e_B(S)$. State constraints $\mathcal{C}$ are conditions that must be satisfied by all states. If one of the state constraints $c \in \mathcal{C}$ is violated, denoted by $S \nvDash \mathcal{C}$, then $V(S) = \infty$ ($V(S) = -\infty$) for minimization (maximization). The dual bound $h(S)$ is a lower (upper) bound on $V(S)$ for minimization (maximization).

Overall, if the problem is minimization, the DP formulation is defined as follows.

$$ \text{compute } V(S^0) \tag{1} $$

$$ \text{s.t. } V(S) = \begin{cases} \min_{\tau \in \mathcal{T}(S)} \text{cost}_{\tau}(V(S[\llap{\tau}]), S) & \text{if } S \models \mathcal{C} \land \forall B \in \mathcal{B}, S \nvDash B \\ e_B(S) & \text{if } S \models \mathcal{C} \land \exists B \in \mathcal{B}, S \models B \\ \infty & \text{if } S \nvDash \mathcal{C} \end{cases} \tag{2} $$

$$ V(S) \leq V(S') \qquad \text{if } S' \preceq S \tag{3} $$

$$ V(S) \geq h(S). \tag{4} $$

The first line states that the optimal objective value is $V(S^0)$. Equation (2) recursively defines the value function $V$. Inequalities (3) and (4) are bounds on the value function. For maximization, we replace min with max in the first line of Equation (2) and swap $\le$ and $\ge$ in Inequalities (3) and (4). A solution for the DP formulation is a sequence of transitions that transforms the target state $S^0$ into a base state. Concretely, for a sequence of transitions $x = \langle x_1, ..., x_n \rangle$, let $S^{i+1} = S^i[\llap{\tau}x_{i+1}]$ for $i = 0, ..., n-1$. Then, $x$ is a solution if $x_{i+1} \in \mathcal{T}(S^i)$, $S^i \models \mathcal{C}$, and $\forall B \in \mathcal{B}, S^i \nvDash B$ for $i = 0, ..., n-1$, $S^n \models \mathcal{C}$, and $\exists B \in \mathcal{B}, S^n \nvDash B$. The solution is an optimal solution if $\text{cost}_{\tau_i}(V(S^{i+1}), S^i) = V(S^i)$ for $i = 0, ..., n-1$ in addition.

## 2.2 Complete Anytime Beam Search for DIDP

Previous research has shown that a subset of DyPDL models can be solved by cost-algebraic heuristic search [11], a generalized version of the shortest path algorithm [25]. Multiple DIDP solvers using cost-algebraic heuristic search algorithms have been proposed [25, 26]. These solvers perform state space search, which finds a path from the target state to a base state in a state space graph, a directed graph where each vertex is a state. In the state space graph, an edge from state $S$ to $S[\llap{\tau}]$ exists if $\tau \in \mathcal{T}(S)$. Following the previous work, we assume that $\text{cost}_\tau(V(S[\llap{\tau}]), S)$ is expressed as $w_\tau(S) \times V(S[\llap{\tau}])$ where $w_\tau$ is an expression returning
■ **Algorithm 1** Beam Search for DyPDL.

1: **function** BEAMSEARCH($\bar{f}$, $b$)
2:     $g(S^0) \leftarrow 0$, $f(S^0) \leftarrow h(S^0)$, $x(S^0) \leftarrow \emptyset$.
3:     $O \leftarrow \{S^0\}$, $\bar{x} \leftarrow \text{NULL}$, **complete** $\leftarrow \top$.
4:     **while** $O \neq \emptyset$ and $\bar{x} = \text{NULL}$ **do**
5:         $G \leftarrow \emptyset$.
6:         **for all** $S \in O$ **do**
7:             **if** $\exists B \in \mathcal{B}$ such that $S \models B$ **then**
8:                 **if** $g(S) \times e_B(S) < \bar{f}$ **then**
9:                     $\bar{f} \leftarrow g(S) \times e_B(S)$, $\bar{x} \leftarrow x(S)$.
10:                **else**
11:                    **for all** $\tau \in \mathcal{T}(S) : S[\tau] \models C$ **do**
12:                        **if** $\exists S' \in G$ such that $S[\tau] \preceq S'$ and $g(S) \times w_\tau(S) \geq g(S')$ **then**
13:                            $x(S[\tau]) \leftarrow \langle x(S); \tau \rangle$.
14:                            $g(S[\tau]) \leftarrow g(S) \times w_\tau(S)$, $f(S[\tau]) \leftarrow g(S[\tau]) \times h(S[\tau])$.
15:                        **if** $\exists S' \in G$ such that $S' \preceq S[\tau]$ and $g(S[\tau]) \leq g(S')$ **then**
16:                            $G \leftarrow G \setminus \{S'\}$.
17:                        **if** $f(S[\tau]) < \bar{f}$ **then**
18:                            $G \leftarrow G \cup \{S[\tau]\}$.
19:                        $O \leftarrow \{S \in G | f(S) < \bar{f}\}$.
20:                    **if** $|O| > b$ **then**
21:                        $O \leftarrow$ the best $b$ states in $G$ minimizing $f$.
22:                        **complete** $\leftarrow \top$.
23:                    **if** $O \neq \emptyset$ **then**
24:                        **complete** $\leftarrow \top$.
25:                    **return** $\bar{x}$, **complete**.

a real value and $\times$ is a binary operator satisfying a cost-algebra. In particular, we focus on nonnegative $w_\tau$ and binary operators $+$ or $\max$, i.e., $\text{cost}_\tau(V(S[\tau]), S) = w_\tau(S) + V(S[\tau])$ or $\text{cost}_\tau(V(S[\tau]), S) = \max\{w_\tau(S), V(S[\tau])\}$. The weight of the edge $(S, S[\tau])$ is defined as $w_\tau(S)$, and the cost of a path from the target state $S^0$, which corresponds to a sequence of the transitions $x = (x_1, ..., x_n)$, is $\text{cost}_x(S^0) = \times_{i=0}^{n-1} w_{x_{i+1}}(S^i)$ where $S^i = [x_{i+1}]$ for $i = 0, ..., n-1$. As each edge weight is nonnegative, the cost of a path is nonnegative and non-decreasing in length. In this paper, we focus on minimization while DyPDL and cost-algebraic heuristic search can handle both minimization and maximization.

The state-of-the-art cost-algebraic heuristic search solver is complete anytime beam search (CABS) [36, 26]. CABS performs beam search, which searches at most $b$ states in the *open list* $O$ at each layer of the state space graph. We show the pseudo-code of beam search in Algorithm 1. In addition to a DyPDL model and $b$, beam search takes the primal bound $\bar{f}$ as an input, which is the best-known objective value and could be infinity. With each state $S$, the best path $x(S)$ from $S^0$ and its cost $g(S) = \text{cost}_{x(S)}(S^0)$ (the *g-value*) are maintained in lines 13 and 14. Starting from $O = \{S^0\}$, beam search processes a state $S$ in $O$. If $S$ is a base state, and the best path to $S$ has a better cost than $\bar{f}$ in line 8, then $\bar{f}$ and the solution $\bar{x}$ are updated. Otherwise, $S[\tau]$ is added to the candidate set $G$ for each transition $\tau \in \mathcal{T}(S)$, which is called the *expansion* of $S$, and we say that $S$ is expanded (lines 11–18). When expanding $S \in O$, if there exists a state $S' \in G$ that dominates $S[\tau]$ and $g(S') \le g(S) \times w_\tau(S)$, then $S[\tau]$ is not added to $G$. In addition to $g(S[\tau])$, the priority
$f(S[[\tau]]) = g(S[[\tau]]) \times h(S[[\tau]])$ (the f-value) is computed in line 14, which is a lower bound on the optimal path cost from $S^0$ to a base state via $S[[\tau]]$. If $f(S[[\tau]]) \ge \bar{f}$, the corresponding path does not lead to an improved solution, so $S[[\tau]]$ is not added to $G$ in line 18. After expanding all states, $O$ is updated to the best $b$ states in $G$ according to $f$ in lines 20-21. This procedure is repeated until a solution whose cost is better than the given primal bound is found, or no successor states are generated. The variable complete maintains whether the search is complete. If states in $G$ are pruned due to the beam width $b$ (line 22), or $O$ is not empty when a solution is found (line 24), there may exist a better solution, so the search is not complete. If complete = $\top$, then $\bar{x}$ is the optimal solution if it is not NULL, or the model is infeasible if $\bar{x} = NULL$. CABS performs a sequence of beam search with exponentially increasing beam width $b = 1, 2, 4, ...$ using the best objective value found so far as the primal bound $\bar{f}$ until the optimality of the best solution or the infeasibility is proved.

## 2.3 Large Neighborhood Search

Large neighborhood search (LNS) iteratively removes a part of a solution and solves the resulting subproblem (neighborhood) [34]. A solution for a CP problem is represented as a complete value assignment to decision variables. Given a solution, LNS removes a subset of the value assignments and solves the subproblem where the remaining variables are fixed to the values assigned in the original solution. Typically, a tree search algorithm is used. In a tree search algorithm, a search node is a partial value assignment to decision variables, successor nodes are generated by assigning a value to an unassigned variable, and a solution corresponds to a leaf node, where all variables are assigned values.

### 2.3.1 Large Neighborhood Search with Decision Diagrams

For combinatorial optimization, recent work proposed LNS with decision diagrams (DD-LNS) [18]. Although DD-LNS was not explicitly framed as state space search, we interpret it as a state space search algorithm. While DD-LNS is independent of DIDP, it also uses the DP formulation of a problem as input while assuming that the solution has $n$ transitions. Given a sequence of transitions $\langle x_1, ..., x_n \rangle$, DD-LNS keeps the first $d$ transitions, $\langle x_1, ..., x_d \rangle$, and searches for the remaining $n - d$ transitions. To find such a sequence, DD-LNS constructs a decision diagram (DD), a directed graph where nodes are partitioned into layers. In the constructed DD, each vertex corresponds to a state, and each edge corresponds to a transition, so it is a state space graph. The first layer in the DD contains only the node corresponding to $S^d$, where $S^i = S^{i-1}[x_i]$ for $i = 1, ..., d$. DD-LNS iteratively constructs a layer by applying transitions to the states in the current layer until reaching layer $n - d + 1$.

Because constructing an exact DD is intractable, DD-LNS constructs a restricted DD, which keeps a subset of states in the exact DD. In each layer, states satisfying certain conditions are kept, some of which are selected randomly. Let the number of such states be $K$. If $K$ is smaller than a parameter $W$, from the remaining states, DD-LNS also keeps the best $W - K$ states that minimize a priority function, the rough lower bound (RLB). The RLB of a state is a lower bound on the cost of a solution path via that state, which is the same as the f-value. If the RLB of a state is larger than the best solution cost, the state is removed from the DD as it does not lead to a better solution. Therefore, the procedure of constructing a restricted DD can be considered beam search with randomization. Indeed, the authors acknowledged that DD-LNS is a hybridization of LNS and beam search [18]. DD-LNS decreases $d$ by 1 if a better solution is not found with $d$, starting from $d = n - 2$ and restarting from $d = n - 2$ if $d = 0$. When $d = 0$ and the restricted DD keeps all states except
<figure><img src="image_218.png" /><figcaption>■ **Figure 1** Partial state space graph induced by the prefix $\langle 1 \rangle$ and the suffix $\langle 4 \rangle$ (highlighted in red) in Example 1. The current partial path is highlighted in blue and an alternative partial path is highlighted in green. Dashed transitions conflict with the suffix (explained in Section 3.2).</figcaption></figure>

for those removed based on RLB, then it is the exact DD, so DD-LNS proves the optimality of the solution. Since DD-LNS can be interpreted as a state space search algorithm and is designed for combinatorial optimization, we implement it for DIDP and experimentally compare it with our method.

## 3 Large Neighborhood Beam Search

We start with a simple idea of LNS for state space search: given a solution path, $x = \langle x_1, \dots, x_n \rangle$ which connects states $\langle S^0, \dots, S^n \rangle$, we remove a partial path $\langle x_i, \dots, x_{i+d-1} \rangle$ and search for a better partial path from $S^{i-1}$ to $S^{i+d-1}$. If we find a better solution $\langle x_1, \dots, x_{i-1}, x'_i, \dots, x'_{i+d'-1}, x_{i+d}, \dots, x_n \rangle$, we repeat this procedure with the new solution. While the overview of the algorithm is simple, there are design choices on how to select a partial path to remove and how to search for a better partial path. The novelty of our method compared to existing methods arises from such choices in addition to the fact that it is used for DIDP. First, we describe the modifications of beam search for DyPDL to search for a partial path from $S^{i-1}$ to $S^{i+d-1}$. Then, we propose strategies to select a partial path to remove.

### 3.1 Beam Search for DyPDL in a Partial State Space Graph

We want to find a path from $S^{i-1}$ to $S^{i+d-1}$ instead of from $S^0$ to a base state. We could modify line 8 in Algorithm 1 so that it checks if $S = S^{i+d-1}$ instead of $\exists B \in \mathcal{B}, S \models B$. However, in DyPDL, it may not be desirable as shown in the following example.

**Example 1.** Consider the following DP formulation, where $U \subseteq \{0, 1, 2, 3, 4\}$ is a set variable, $k \in \{0, 1, 2, 3, 4\}$ is an element variable, and $c_{ij}$ for $l, j \in \{0, 1, 2, 3, 4\}$ is a constant.

compute $V(\{1, 2, 3, 4\}, 0)$

$$V(U, k) = \begin{cases} \min_{j \in U} c_{kj} + V(U \setminus \{j\}, j) & \text{if } U \neq \emptyset \\ 0 & \text{if } U = \emptyset. \end{cases}$$

Each transition in the DyPDL model has precondition $j \in U$ and effect $(U, U \setminus \{j\})$ for some $j$, and so we denote each transition by $j$. Each solution corresponds to a permutation of the transitions 1, 2, 3, 4. A solution $\langle 1, 2, 3, 4 \rangle$ connects a sequence of states $\langle (\{1, 2, 3, 4\}, 0), (\{2, 3, 4\}, 1), (\{3, 4\}, 2), (\{4\}, 3), (\emptyset, 4) \rangle$. Consider removing $\langle 2, 3 \rangle$ from the solution. We visualize the partial state space graph in Figure 1. An algorithm tries to find a path from $(\{2, 3, 4\}, 1)$ to $(\{4\}, 3)$. The original one, $\langle 2, 3 \rangle$, is only the path. However, a partial path $\langle 3, 2 \rangle$ from $(\{2, 3, 4\}, 1)$ to $(\{4\}, 2)$ also results in a valid solution $\langle 1, 3, 2, 4 \rangle$.
<caption>**Algorithm 2 Beam Search in a partial state space graph.**</caption>

1: **function** BEAMSEARCHFORPARTIALPATH($\bar{f}$, $b$, **prefix**, **suffix**)
2:     $\hat{S} \leftarrow S^0[[\text{prefix}]]$, $g(\hat{S}) \leftarrow \text{cost}_{\text{prefix}}(S^0)$, $f(\hat{S}) \leftarrow g(\hat{S}) \times h(\hat{S})$, $x(\hat{S}) \leftarrow \text{prefix}$.
3:     $O \leftarrow \{\hat{S}\}$, $\bar{x} \leftarrow \text{NULL}$, **complete** $\leftarrow \top$.
4:     **while** $O \neq \emptyset$ and $\bar{x} = \text{NULL}$ **do**
5:         $G \leftarrow \emptyset$. ▷ A set of states in the next layer.
6:         **for all** $S \in O$ **do**
7:             $S' \leftarrow S$, **success** $\leftarrow \bot$.
8:             **if** $\exists B \in \mathcal{B}$ and $S' \models B$ **then**
9:                 **success** $\leftarrow \top$. ▷ A base state, success.
10:             **else**
11:                 **for** $\tau \leftarrow \text{suffix}_1, \dots, \text{suffix}_{n-i-d+1}$ **do** ▷ Rollout of the suffix.
12:                 **if** $S' \not\models \text{pre}_\tau$ **then**
13:                     **break**. ▷ Preconditions are not satisfied, fail.
14:                 $S' \leftarrow S'[[\tau]]$, $g(S') \leftarrow g(S') \times w_\tau(S')$, $x(S') \leftarrow \langle x(S'); \tau \rangle$.
15:                 **if** $S' \not\models \mathcal{C}$ **then**
16:                     **break**. ▷ State constraints are not satisfied, fail.
17:                 **if** $\exists B \in \mathcal{B}$ such that $S' \models B$ **then**
18:                     **success** $\leftarrow \top$. ▷ A base state, success.
19:                 **break**.
20:             **if** **success** **then**
21:                 **if** $g(S') \times e_B(S') < \bar{f}$ **then** ▷ A better solution.
22:                 $\bar{f} \leftarrow g(S') \times e_B(S')$, $\bar{x} \leftarrow x(S')$.
23:             **else**
24:                 **for all** $\tau \in \mathcal{T}(S) : S[[\tau]] \models \mathcal{C}$ **do**
25:                     **if** $\nexists S' \in G$ such that $S[[\tau]] \preceq S'$ and $g(S) \times w_\tau(S) \ge g(S')$ **then**
26:                         $x(S[[\tau]]) \leftarrow \langle x(S); \tau \rangle$.
27:                         $g(S[[\tau]]) \leftarrow g(S) \times w_\tau(S)$, $f(S[[\tau]]) \leftarrow g(S[[\tau]]) \times h(S[[\tau]])$.
28:                     **if** $\exists S' \in G$, $S' \preceq S[[\tau]] \wedge g(S[[\tau]]) \le g(S')$ **then**
29:                         $G \leftarrow G \setminus \{S'\}$. ▷ Remove a dominated state.
30:                     **if** $f(S[[\tau]]) < \bar{f}$ **then** ▷ Pruning by the primal bound.
31:                         $G \leftarrow G \cup \{S[[\tau]]\}$.
32:     $O \leftarrow \{S \in G \mid f(S) < \bar{f}\}$.
33:     **if** $|G| > b$ **then**
34:         $O \leftarrow$ the best $b$ states in $G$ minimizing $f$.
35:         **complete** $\leftarrow \bot$.
36:     **if** $O \neq \emptyset$ **then**
37:         **complete** $\leftarrow \bot$.
38:     **return** $\bar{x}$, **complete**.

Considering the above example, instead of focusing on a partial path to a state, we focus on a partial path to a *suffix* of the solution path. Given a solution path $\langle x_1, \dots, x_n \rangle$, if we remove a partial path $\langle x_i, \dots, x_{i+d-1} \rangle$, then $\langle x_1, \dots, x_{i-1} \rangle$ is the prefix, and $\langle x_{i+d}, \dots, x_n \rangle$ is the suffix. For a partial path $\langle x'_i, \dots, x'_{i+d'-1} \rangle$, we want to check if $\langle x_1, \dots, x_{i-1}, x'_i, \dots, x'_{i+d'-1}, x_{i+d}, \dots, x_n \rangle$ is a valid solution. Therefore, for a state $S$ found by a search algorithm, we perform a rollout of the suffix from $S$ and check if each of resulting states satisfies the state constraints and a
■ **Algorithm 3** Large Neighborhood Beam Search (LNBS).

**Input:** initial feasible solution $\bar{x}$.

**Output:** solution $\bar{x}$ and if the optimality or infeasibility is proved.

1. **while** the time limit is not reached **do**

2. $n \leftarrow |\bar{x}|, \bar{f} \leftarrow \text{cost}_{\bar{x}}(S^0)$.

3. Select $d$ such that $2 \le d \le n$.

4. Select $i$ such that $1 \le i \le n - d + 1$.

5. Select beam width $b$.

6. **prefix** $\leftarrow \langle \bar{x}_1, \dots, \bar{x}_{i-1} \rangle$, **suffix** $\leftarrow \langle \bar{x}_{i+d}, \dots, \bar{x}_n \rangle$.

7. $x$, **complete** $\leftarrow \text{BEAMSEARCHFORPARTIALPATH}(\bar{f}, b$, **prefix**, **suffix**)

8. **if** $x \neq \text{NULL}$ **then**

9. $\bar{x} \leftarrow x$.

10. **if** $i = 0 \land d = n \land \text{complete}$ **then**

11. **return** $\bar{x}$, $\top$.

12. **return** $\bar{x}$, $\bot$.

base case. We show the modified version of beam search in Algorithm 2. This algorithm takes a prefix prefix and a suffix suffix as input. In line 2, we denote the state resulting from applying the prefix to the target state by $\hat{S} = S^0 \llbracket \text{prefix} \rrbracket$ and initialize the open list $O$ with $\hat{S}$ in line 3. In lines 8–22, the algorithm performs a rollout of the suffix from $S$ and checks if it results in a better solution. Other parts are the same as Algorithm 1.

We use this modified version of beam search in large neighborhood beam search (LNBS) as shown in Algorithm 3. In line 2, $|\bar{x}|$ denotes the length of the current solution $\bar{x}$. In lines 3–5, LNBS selects parameters $d$, $i$, and $b$. In line 7, LNBS performs beam search in the neighborhood. If an improving solution is found, LNBS updates the current solution $\bar{x}$ in line 9. If the searched neighborhood is the original search space, i.e., $i = 1$ and $d = n$, and beam search proves the optimality or infeasibility, LNBS terminates in line 11. Therefore, if it is guaranteed to select $i = 1$ and $d = n$ with sufficiently large $b$ given enough time, LNBS is guaranteed to find the optimal solution or prove the infeasibility, i.e., it is complete. CABS can be considered a configuration of LNBS, where $i = 1$, $d = n$, and $b$ increases exponentially. DD-LNS can also be considered a configuration of LNBS, where $i$ ranges from $n-2$ to 1, $d$ is $n-i+1$, and $b$ is fixed to $W$ while beam search is extended with the randomization mechanism. We will describe the strategies that we use to select $d$, $i$, and $b$ below.

## 3.2 Removing Conflicting Transitions

In Example 1, consider finding a partial path from a prefix $\langle 1 \rangle$, which results in state $\hat{S} = (\{2, 3, 4\}, 1)$, to a suffix $\langle 4 \rangle$ using beam search. In $\hat{S}$, three transitions 2, 3, and 4 are applicable. However, applying transition 4 does not lead to a feasible solution because it is already used in the suffix and cannot be applied twice: it requires $4 \in U$ and removes 4 from $U$, but no other transition adds 4 to $U$, so applying 4 makes the suffix inapplicable. Generalizing this example, if we know that a transition $\tau$ makes a transition $\tau'$ in the suffix inapplicable, then we can ignore $\tau$ when searching for a partial path. In particular, we focus on the effects of $\tau$ that add/remove an element to/from a set variable and the preconditions of $\tau'$ that require the element to be/not to be in that set variable.

**Proposition 2.** Suppose that a DyPDL model $\langle \mathcal{V}, S^0, \mathcal{K}, \mathcal{T}, \mathcal{B}, \mathcal{C}, h \rangle$ has a set variable $U \in \mathcal{V}$ whose domain is $2^N$, where $N$ is a set of objects. There does not exist a solution $\langle x_1, \dots, x_n \rangle$ such that any pair of $\tau = x_i$ and $\tau' = x_j$ for $1 \le i < j \le n$ satisfy either of the following