26:12 Proof Logging for Smart Extensional Constraints

It should be noted that the same justifications apply when building the graph for the first time, as it can be conceptually viewed as eliminating edges from the complete layered graph with $n + 1$ layers and $|Q|$ nodes in each layer. The deletion occurs in two passes. Firstly it eliminates all outgoing edges from every node in layer 0 that does not correspond to the initial state, recursing for each subsequent layer. Then secondly, it eliminates all incoming edges from every node in layer $n + 1$ that does not correspond to a final state, similarly recursing for each previous layer.

Note finally that stronger encodings of the Regular constraint are possible, and these would not require as much proof logging during search. However, these encodings are not so obviously correct, which is a drawback since the constraint encoding process is not (currently) verified. It is not even clear that such encodings would give more efficient proof verification, since they would involve having a larger set of active PB constraints for the entire verification process.

4.3 A Complete Worked Example for Regular Propagation

Once again, to demonstrate how the proof logging procedure works more concretely, we will now show a worked example, including the initial building of the graph and a round of domain-consistent propagation. We will take the simple Example 1 from Pesant [19], of a regular language membership constraint on five variables $X_1, \dots, X_5$ each with domain $\{0, 1, 2\}$, and using the DFA $M$ as shown in Figure 2a.

For the PB encoding, we will omit the constraints encoding the domains and equals flags and focus on the constraints present in Encoding 2. We can see that we have five states numbered 0...4, and five variables in the sequence, so the model would first define:

$$r_{00} + \dots + r_{04} \geq 1; \quad \dots \quad r_{50} + \dots + r_{54} \geq 1; \quad (39)$$

$$-r_{00} - \dots - r_{04} \geq -1; \quad \dots \quad -r_{50} - \dots - r_{54} \geq -1; \quad (40)$$

saying that we have to have exactly one of the state-position flags set for each position. Next for each position $i$ we would define eight PB constraints corresponding to the eight valid transitions:

$$\begin{align}
\bar{r}_{i0} + \bar{x}_{i=0} + r_{i+11} \geq 1; \quad \bar{r}_{i0} + \bar{x}_{i=2} + r_{i+14} \geq 1; \quad \bar{r}_{i1} + \bar{x}_{i=0} + r_{i+11} \geq 1; \nonumber \\
\bar{r}_{i1} + \bar{x}_{i=1} + r_{i+12} \geq 1; \quad \bar{r}_{i2} + \bar{x}_{i=0} + r_{i+13} \geq 1; \quad \bar{r}_{i2} + \bar{x}_{i=1} + r_{i+12} \geq 1; \tag{41} \\
\bar{r}_{i3} + \bar{x}_{i=0} + r_{i+13} \geq 1; \quad \bar{r}_{i4} + \bar{x}_{i=2} + r_{i+14} \geq 1; \nonumber
\end{align}$$

along with constraints of the form

$$\bar{r}_{iq} + \bar{x}_{i=v} \geq 1; \quad (42)$$

for the remaining seventeen invalid transitions. Finally, we would have constraints saying that we have to start in an initial state and end in a final state:

$$r_{00} \geq 1; \quad r_{53} + r_{54} \geq 1. \quad (43)$$

Propagation depends on us having first built a valid graph, with nodes as shown in Figure 2b. As mentioned, this occurs in two passes, with the forward pass collecting nodes that can be reached from the initial state. This corresponds to removing outgoing edges only reachable from nodes other than the initial state.

So to begin with, we traverse the two edges incident to node $q_0^0$, reaching nodes $q_1^1$ and $q_4^1$ and thereby eliminating all other edges. These are the two edges shown in Figure 2b. Since we are eliminating outgoing edges we are effectively in the third case of the proof logging
<figure><img src="image_267.png" /><figcaption>(a) A DFA $M$ recognising the regular expression `00*11*00**2*`. Double circles indicate accepting states.</figcaption></figure>

<figure><img src="image_268.png" /><figcaption>(c) The same multigraph after the forward pass is complete.</figcaption></figure>

<figure><img src="image_269.png" /><figcaption>(e) After inferring that $X_3 \neq 1$, the multigraph is in an invalid state.</figcaption></figure>

<figure><img src="image_270.png" /><figcaption>(b) The first traversal in the forward pass for building the layered multigraph used in propagation of Regular($X_1, \dots, X_5, M$).</figcaption></figure>

<figure><img src="image_271.png" /><figcaption>(d) The same multigraph after the backwards pass is complete.</figcaption></figure>

<figure><img src="image_272.png" /><figcaption>(f) The same multigraph after its consistency properties have been re-established.</figcaption></figure>

■ **Figure 2** Propagation of the regular expression `00*11*00**2*` on the variables $X_1, \dots, X_5$, showing the initial state, and then the effect of propagating $X_3 \neq 1$.
26:14 Proof Logging for Smart Extensional Constraints

procedure from the previous section, however, as there are no previous states to eliminate we simply need to log constraints in the form of (34) for each of the six eliminated edges:

$$
\begin{align}
\bar{r}_{01} + \bar{x}_{0=0} \geq 1; \quad & \bar{r}_{01} + \bar{x}_{0=1} \geq 1; \quad & \bar{r}_{02} + \bar{x}_{0=0} \geq 1; \nonumber \\
\bar{r}_{02} + \bar{x}_{0=1} \geq 1; \quad & \bar{r}_{03} + \bar{x}_{0=0} \geq 1; \quad & \bar{r}_{04} + \bar{x}_{0=2} \geq 1; \tag{44}
\end{align}
$$

and these will all follow by RUP because we have $r_{00} \ge 1$ in the PB model. Continuing, we recursively collect edges reachable from the initial node, consequently eliminating all outgoing edges for any node that is not reached. The result of this is shown in Figure 2c.

At each layer, we can follow this same proof logging procedure for eliminating outgoing edges. For example, we do not collect $(q_3^2, q_3^3)$ at layer 2, since $q_3^2$ is not reached as part of the forward pass. So we would like to log

$$
\bar{r}_{23} + \bar{x}_{2=0} \geq 1; \tag{45}
$$

for this eliminated edge, but this does not follow by RUP. Despite us having logged constraints
in the form of (34) for each eliminated incoming edge to $q_2^3$, none of these constraints contain
either $r_{23}$ or $x_{2=0}$, as they concern the previous layer, and so no further propagation takes
place from the negation of (45). Hence, as discussed, we first log

$$
\bar{r}_{10} + \bar{r}_{23} \ge 1; \quad \bar{r}_{11} + \bar{r}_{23} \ge 1; \quad \bar{r}_{12} + \bar{r}_{23} \ge 1; \quad \bar{r}_{13} + \bar{r}_{23} \ge 1; \quad \bar{r}_{14} + \bar{r}_{23} \ge 1; \quad (46)
$$

which all do follow by RUP, and then our constraint (45) can be logged.

After this process is complete, the backwards pass takes place, first eliminating incoming edges from any nodes in the last layer corresponding to non-final states, and then recursively eliminating all incoming edges from any node that lost all of its outgoing edges. The complete, correctly initialised state of the graph after this backwards pass is shown in Figure 2d.

Following the proof logging procedure for eliminating incoming edges, in this direction
we simply need to log a constraint in the form of (34) for each, and these will follow by RUP,
as discussed.

At this point, some basic inferences about variable value pairs can already be made, such as $X_0 \neq 1$ and $X_4 \neq 1$. Due to all the constraints that have already been logged, eliminating all possible corresponding edges from consideration, the required proof logging invariant constraint for each will follow by RUP. This concludes the graph initialisation for the Regular propagator, and from this point onwards the structure can be updated incrementally and restored upon backtrack, as described in detail by Pesant [19].

We will now demonstrate one such incremental update, corresponding to an execution of
domain-consistent propagation for Regular. Suppose through the course of the computation,
the assignment $X_3 = 1$ loses support. This means the edges $(q_1^3, q_2^4)$ and $(q_2^3, q_2^4)$ are
immediately removed from the graph, and we log the RUP constraints:

$$
\bar{r}_{13} + \bar{x}_{3=1} \geq 1; \qquad \bar{r}_{23} + \bar{x}_{3=1} \geq 1. \tag{47}
$$

At this point the state of the graph is as shown in Figure 2e, with the two highlighted
nodes having lost all of their outgoing and incoming edges respectively.

The recursive incremental update is then executed for each of these, further removing
the edge $(q_1^2, q_1^3)$ in a backwards pass, and so logging the RUP constraint

$$
\wedge \mathcal{G} \implies \bar{r}_{31} + \bar{x}_{2=0}. \tag{48}
$$
The edge $(q_2^4, q_3^5)$ is also removed, but because this is removed in a forward pass, we first log

$$ \wedge \mathcal{G} \Rightarrow \bar{r}_{30} + \bar{r}_{42} \geq 1; \quad \wedge \mathcal{G} \Rightarrow \bar{r}_{31} + \bar{r}_{42} \geq 1; \quad \dots \quad \wedge \mathcal{G} \Rightarrow \bar{r}_{34} + \bar{r}_{42} \geq 1; \quad (49) $$

before logging the required

$$ \wedge \mathcal{G} \Rightarrow \bar{r}_{42} + \bar{x}_{2=0}. \quad (50) $$

No further nodes lose their last incoming or outgoing edge as a result of this, and so the required properties (and hence domain consistency) have been re-established on the graph. The final consistent state is shown in Figure 2f.

# 5 Implementation and Validation

We have implemented both of these proof logging propagation algorithms inside the open-source Glasgow Constraint Solver [17]. Our implementation adds two new constraints Regular and SmartTable, with specific functionality to produce the required encodings and justifying propagators. They can therefore be used in conjunction with the rest of the constraint types already available in the solver.

We validated the implementations by first modelling some key examples. For SmartTable, these included the representation of lexicographic ordering problems, and problems where at most one variable takes a value, as given by Mairy et al. [16], as well as the illustrative examples from Section 3. For Regular, we implemented both examples 1 and 2 from Pesant [19].

We carried out further experimental validation by generating random (acyclic) smart tables, and random DFAs on sequences of up to five variables and solving the corresponding single-constraint problems. In all tests and examples, the solver produced an OPB model and proof files, and we checked these using the VeriPB proof checker. We also found that although proof logging incurs an obvious performance cost, the observed overheads were not unreasonable, giving a slowdown factor of between 2 and 10. As proofs are written currently written directly to disk, this can be very hardware dependent, and also dependent on how optimised the propagator implementation is. We decided to leave further engineering and optimisation of proof writing to future work.

During development some very subtle bugs in both propagator implementations that had eluded conventional testing were caught immediately by proof logging. For example, in an incremental version of the Regular propagator, an unsound inference was being made due to mixing up variable names, but in such a way that a situation where this unsound inference would actually lead to an incorrect solution was extremely rare (only in specially constructed instances, created once we were made aware of the bug due to proof logging). After correcting bugs such as these, all proofs were certified by VeriPB as being correct.

# 6 Conclusion

We have shown that we can efficiently justify all the reasoning that could possibly be carried out for two families of smart extensional constraints. An interesting observation is that justifying this reasoning required only the RUP rule, and this rule was only to provide hints of algorithmic steps that were already being carried out by the propagators. In effect, we are logging a sequence of “lookahead to see immediate contradiction” steps. We did not require any explicit cutting planes derivations, and although we relied upon pseudo-Boolean
26:16 Proof Logging for Smart Extensional Constraints

constraints to make it simple to express reifications and negations, in principle everything we
did should also be possible in a weaker proof encoding and proof system such as CNF and
DRAT. The only caveat is that we do rely upon strong propagation properties for encoded
integer variables, which would limit approaches based upon Boolean satisfiability to integer
variables with very small domains. This is in contrast to constraints like AllDifferent and
Linear, which cannot be logged efficiently in resolution-based approaches [8, 10].

We expect that other global constraints, even those with complex propagation algorithms
will be similarly feasible to proof log using this technique. Richer smart tables, such as those
with offset or ternary restrictions [4], also fit into the framework given in Section 3.2, although
the justification of the filtering inferences for each restriction may require additional cutting
planes steps. Furthermore, the propagation of “Multi-Valued Decision Diagram”-based
constraints [14] can be viewed as a generalisation of the techniques used for the regular
language membership constraint, and so it seems very plausible that a similar proof logging
methodology as demonstrated in this paper could work here. This bodes well for being
able to provide auditable solving for most global constraints that might occur in a modern
constraint solver.

References

1 Fahemi Bacchus. GAC Via Unit Propagation. In Christian Bessière, editor, *Principles and Practice of Constraint Programming – CP 2007*, Lecture Notes in Computer Science, pages 133–147, Berlin, Heidelberg, 2007. Springer. doi:10.1007/978-3-540-74970-7_12.

2 Seulkee Baek, Mario Carneiro, and Marijn J. H. Heule. A flexible proof format for SAT solver-elaborator communication. In *Tools and Algorithms for the Construction and Analysis of Systems: 27th International Conference, TACAS 2021, Held as Part of the European Joint Conferences on Theory and Practice of Software, ETAPS 2021, Luxembourg City, Luxembourg, March 27 – April 1, 2021, Proceedings, Part I*, pages 59–75, Berlin, Heidelberg, 2021. Springer-Verlag. doi:10.1007/978-3-030-72016-2_4.

3 Bart Bogaerts, Stephan Gocht, Ciaran McCreesh, and Jakob Nordström. Certified Symmetry and Dominance Breaking for Combinatorial Optimisation, March 2022. doi:10.48550/arXiv. 2203.12275.

4 Frédéric Boussemart, Christophe Lecoutre, and Cédric Piette. XCSP3: an integrated format for benchmarking combinatorial constrained problems. *CoRR*, abs/1611.03398, 2016. arXiv: 1611.03398.

5 Chiu Wo Choi, Warwick Harvey, J. H. M. Lee, and Peter J. Stuckey. Finite domain bounds consistency revisited. In *AI 2006: Advances in Artificial Intelligence, 19th Australian Joint Conference on Artificial Intelligence, Hobart, Australia, December 4-8, 2006, Proceedings*, pages 49–58, 2006. doi:10.1007/11941439_9.

6 W. Cook, C.R. Coullard, and Gy. Turán. On the complexity of cutting-plane proofs. *Discrete Applied Mathematics*, 18(1):25–38, 1987. doi:10.1016/0166-218X(87)90039-4.

7 Luís Cruz-Filipe, Marijn J. H. Heule, Warren A. Hunt, Matt Kaufmann, and Peter Schneider-Kamp. Efficient certified RAT verification. In Leonardo de Moura, editor, *Automated Deduction – CADE 26*, pages 220–236, Cham, 2017. Springer International Publishing. doi:10.1007/978-3-319-63046-5_14.

8 Jan Elffers, Stephan Gocht, Ciaran McCreesh, and Jakob Nordström. Justifying All Differences Using Pseudo-Boolean Reasoning. *Proceedings of the AAAI Conference on Artificial Intelligence*, 34(02):1486–1494, April 2020. doi:10.1609/aaai.v34i02.5507.

9 Stephan Gocht. *Certifying Correctness for Combinatorial Algorithms by Using Pseudo-Boolean Reasoning*. PhD thesis, Lund University, Sweden, 2022. URL: https://lup.lub.lu.se/record/3550cb96-83d5-4fc7-9e62-190083a3c10a.
10 Stephan Gocht, Ciaran McCreesh, and Jakob Nordström. An Auditable Constraint Programming Solver. 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 25:1–25:18, Dagstuhl, Germany, 2022. Schloss Dagstuhl – Leibniz-Zentrum für Informatik. doi:10.4230/LIPIcs.CP.2022.25.

11 E. Goldberg and Y. Novikov. Verification of proofs of unsatisfiability for cnf formulas. In *2003 Design, Automation and Test in Europe Conference and Exhibition*, pages 886–891, 2003. doi:10.1109/DATE.2003.1253718.

12 Marijn J. H. Heule, Warren A. Hunt, and Nathan Wetzler. Verifying refutations with extended resolution. In Maria Paola Bonacina, editor, *Automated Deduction – CADE-24*, pages 345–359, Berlin, Heidelberg, 2013. Springer Berlin Heidelberg. doi:10.1007/978-3-642-38574-2_24.

13 Marijn J.H. Heule, Warren A. Hunt, and Nathan Wetzler. Trimming while checking clausal proofs. In *2013 Formal Methods in Computer-Aided Design*, pages 181–188, 2013. doi:10.1109/FMCAD.2013.6679408.

14 Samid Hoda, Willem-Jan van Hoeve, and J. N. Hooker. A systematic approach to MDD-Based constraint programming. In David Cohen, editor, *Principles and Practice of Constraint Programming – CP 2010*, pages 266–280, Berlin, Heidelberg, 2010. Springer Berlin Heidelberg. doi:10.1007/978-3-642-15396-9_23.

15 Mikael Z Lagerkvist and Gilles Pesant. Modeling irregular shape placement problems with regular constraints. In *First workshop on bin packing and placement constraints BPPC'08*, 2008.

16 Jean-Baptiste Mairy, Yves Deville, and Christophe Lecoutre. The Smart Table Constraint. In Laurent Michel, editor, *12th International Conference on Integration of AI and OR Techniques in Constraint Programming (CPAIOR 2015)*, Lecture Notes in Computer Science, pages 271–287, Cham, 2015. Springer International Publishing. doi:10.1007/978-3-319-18008-3_19.

17 Ciaran McCreesh and Matthew McIlree. The Glasgow Constraint Solver. GitHub repository, 2023. URL: https://github.com/ciaranm/glasgow-constraint-solver.

18 Olga Ohrimenko, Peter J. Stuckey, and Michael Codish. Propagation = lazy clause generation. In Christian Bessiere, editor, *Principles and Practice of Constraint Programming - CP 2007, 13th International Conference, CP 2007, Providence, RI, USA, September 23-27, 2007, Proceedings*, volume 4741 of *Lecture Notes in Computer Science*, pages 544–558. Springer, 2007. doi:10.1007/978-3-540-74970-7_39.

19 Gilles Pesant. A Regular Language Membership Constraint for Finite Sequences of Variables. In Mark Wallace, editor, *10th International Conference on Principles and Practice of Constraint Programming (CP 2004)*, Lecture Notes in Computer Science, pages 482–495, Berlin, Heidelberg, 2004. Springer. doi:10.1007/978-3-540-30201-8_36.

20 Hélène Verhaeghe. *The extensional constraint*. PhD thesis, Catholic University of Louvain, Louvain-la-Neuve, Belgium, 2021.

21 Hélène Verhaeghe, Christophe Lecoutre, Yves Deville, and Pierre Schaus. Extending Compact-Table to Basic Smart Tables. In J. Christopher Beck, editor, *Principles and Practice of Constraint Programming*, volume 10416, pages 297–307. Springer International Publishing, Cham, 2017. doi:10.1007/978-3-319-66158-2_19.
<figure><img src="image_273.png" /></figure>
# Improving Conflict Analysis in MIP Solvers by Pseudo-Boolean Reasoning

Gioni Mexi ⬒ ⬒ ID

Zuse Institute Berlin, Germany

Timo Berthold ⬒ ⬒ ID

Fair Isaac Deutschland GmbH, Berlin, Germany

TU Berlin, Germany

Ambros Gleixner ⬒ ⬒ ID

HTW Berlin, Germany

Zuse Institute Berlin, Germany

Jakob Nordström ⬒ ⬒ ID

University of Copenhagen, Denmark

Lund University, Sweden

## Abstract

Conflict analysis has been successfully generalized from Boolean satisfiability (SAT) solving to mixed integer programming (MIP) solvers, but although MIP solvers operate with general linear inequalities, the conflict analysis in MIP has been limited to reasoning with the more restricted class of clausal constraint. This is in contrast to how conflict analysis is performed in so-called pseudo-Boolean solving, where solvers can reason directly with 0-1 integer linear inequalities rather than with clausal constraints extracted from such inequalities.

In this work, we investigate how pseudo-Boolean conflict analysis can be integrated in MIP solving, focusing on 0-1 integer linear programs (0-1 ILPs). Phrased in MIP terminology, conflict analysis can be understood as a sequence of linear combinations and cuts. We leverage this perspective to design a new conflict analysis algorithm based on mixed integer rounding (MIR) cuts, which theoretically dominates the state-of-the-art division-based method in pseudo-Boolean solving.

We also report results from a first proof-of-concept implementation of different pseudo-Boolean conflict analysis methods in the open-source MIP solver SCIP. When evaluated on a large and diverse set of 0-1 ILP instances from MIPLIB 2017, our new MIR-based conflict analysis outperforms both previous pseudo-Boolean methods and the clause-based method used in MIP. Our conclusion is that pseudo-Boolean conflict analysis in MIP is a promising research direction that merits further study, and that it might also make sense to investigate the use of such conflict analysis to generate stronger no-goods in constraint programming.

**2012 ACM Subject Classification** Theory of computation → Discrete optimization; Mathematics of computing → Solvers

**Keywords and phrases** Integer programming, pseudo-Boolean solving, conflict analysis, cutting planes proof system, mixed integer rounding, division, saturation

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

**Funding** The work for this article has been conducted within the Research Campus Modal funded by the German Federal Ministry of Education and Research (BMBF grant numbers 05M14ZAM, 05M20ZBM).

Jakob Nordström: supported by the Swedish Research Council grant 2016-00782 and the Independent Research Fund Denmark grant 9040-00389B.

**Acknowledgements** Part of this work was carried out while some of the authors participated in the extended reunion for the program *Satisfiability: Theory, Practice, and Beyond* at the Simons Institute for the Theory of Computing at UC Berkeley in the spring of 2023. This work has also benefited greatly from discussions during the Dagstuhl Seminar 22411 *Theory and Practice of SAT and Combinatorial Solving*.

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

© Gioni Mexi, Timo Berthold, Ambros Gleixner, and Jakob Nordström;

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. 27; pp. 27:1–27:19

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

LIPICS

Leibniz International Proceedings in Informatics

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

The area of Boolean satisfiability (SAT) solving has witnessed dramatic performance improvements over the last couple of decades, and several techniques from SAT have also inspired developments for other combinatorial optimization paradigms such as SAT-based and (linear) pseudo-Boolean optimization, constraint programming, and mixed integer programming. In particular, conflict analysis as introduced in the works on conflict-driven clause learning (CDCL) [3, 38, 40] ushering in the modern SAT solving revolution has been picked up and generalized in different ways to these more general settings. Interestingly, precursors of this version of conflict analysis and nonchronological backtracking can be traced back all the way to early work in the AI community [47], and related ideas have been used in constraint programming for decades [24, 31]. Our focus in this paper is on conflict analysis in mixed integer programming and pseudo-Boolean optimization, which we proceed to discuss next.

## 1.1 Mixed Integer Programming and Conflict Analysis

The core method of mixed integer programming (MIP) is that a linear programming (LP) relaxation of the problem is fed to an LP solver. If the LP solver finds a solution that assigns real values to integral variables, then either additional cut constraints can be generated that eliminate such solutions, or the problem can be split into subproblems by branching on integer variables, generating new nodes in the search tree. During the solving process infeasible nodes in the search tree are pruned. Unlike in SAT, there can be different reasons for backtracking due to infeasibility of the LP relaxation, node presolving (propagation), or to the current objective value of the relaxed problem being worse than the best solution found so far (branch-and-bound). MIP solvers employ a multitude of further techniques such as symmetry detection, disjoint subtree detection, restarts, et cetera. For a comprehensive description of MIP solving we refer the reader to, e.g., [2].

The use of SAT techniques in MIP solvers has been a fruitful direction of research over the last decades. Specifically, CDCL conflict analysis has proven to be a useful tool to enhance the performance of MIP solvers by learning constraints from infeasibilities detected by propagation or from the LP relaxation [1, 44, 48]. However, SAT and MIP solvers differ fundamentally in how they explore the search space, in that SAT solvers search depth-first, maintaining only the current state of the search, whereas in MIP the search tree is generated in a “best-first” manner based on careful analysis on search statistics such as dual bounds and integrality of LP solutions to subproblems. These differences make it harder for MIP solvers to profit from conflict analysis, and so in contrast to SAT solving, for which this technique is absolutely crucial, in MIP solving it plays more of a supplemental if still highly valuable role.

Although the setting is different, the graph-based conflict analysis [1] used to learn from infeasibilities in MIP is very similar to the classic SAT approach. First, a partial assignment is extracted that consists of branching decisions and implications that led to the infeasibility. If the LP relaxation is infeasible, the information which bound changes led to infeasibility is gathered from the non-zero duals of the LP. Next, a directed acyclic graph is constructed that encodes information about the conflict, in that source nodes correspond to branching decisions, non-source nodes encode implications, and the sink node represents the infeasibility. Each cut in this graph that separates the source nodes from the sink is a valid constraint. It is important to note that all implications correspond to clausal constraints, and so this conflict analysis operates not on the linear constraints of the problem but on clauses extracted from these linear constraints. (There are also methods that can learn general linear constraints
from infeasibilities, one notable example being dual-proof analysis [48], but this technique is
limited to conflicts arising from infeasibility of the LP relaxation and does not analyze or
strengthen the partial assignment that led to infeasibility.)

## 1.2 Pseudo-Boolean Solving and Conflict Analysis

Pseudo-Boolean (PB) solving is another approach specific to integer linear programs with only binary variables, or 0-1 ILPs, which are referred to as (linear) pseudo-Boolean formulas in the PB solving literature. While MIP solvers find real-valued solutions and try to push such solutions closer and closer to integrality, PB solvers follow the SAT approach of considering only Boolean assignments and trying to extend partial assignments to more and more variables without violating any constraints. Just as in SAT, this search is performed in a depth-first manner.

Some PB solvers stick very closely to SAT in that they immediately translate the 0-1 ILP into conjunctive normal form (CNF) using auxiliary variables and then run a standard CDCL SAT solver [21, 39, 43]. Another approach, which is what is of interest in the context of this work, is to extend the solvers to reason natively with 0-1 linear inequalities [12, 46, 34, 23]. Such conflict-driven pseudo-Boolean solvers have the potential to run exponentially faster than CDCL-based solvers, since their conflict analysis method is exponentially stronger than that used in CDCL SAT solvers.

Since it is crucial for our work to understand the differences between conflict analysis in MIP and PB solvers, let us try to provide a somewhat simplified exposition of PB solving in a language that is meant to convey a MIP perspective (and where what follows below is heavily indebted to [19]). During the search phase, the pseudo-Boolean solver always first tries to extend the current partial solution with any variable assignments that are propagated by some linear inequality. When no further propagations are possible, the solver chooses some unassigned variable and makes a *decision* to assign this variable 0 or 1, after which it again turns to propagation. This cycle of decisions and propagations repeats until either a satisfying assignment is found or some 0-1 linear inequality *C* is violated. In the latter case, the solver switches to the conflict analysis phase, which works as follows:

1. The linear inequality *R* responsible for propagating the last variable *x* in *C* to the “wrong value” from the point of view of *C* is identified; this inequality *R* is referred to as the *reason constraint* for *x*.

2. A *division* or *saturation* rule is applied to *R* to generate a modified inequality *R*\* that propagates *x* tightly to its assigned value even when considered over the reals.

3. A new linear constraint *D* is computed as the smallest integer linear combination of *R*\* and *C* for which the variable *x* cancels and is eliminated. It is not too hard to show that it follows from the description above that this constraint *D* is violated by the current partial assignment of the solver with the value of *x* removed, and we can set *C* := *D* and go to step 1 again.

This continues until a termination criterion analogous to the *unique implication point (UIP)* notion used in SAT solving leads to *D* being declared as the learned constraint. At this point, the solver undoes further assignments in reverse chronological order until *D* is no longer violated, and then switches back to the search phase. We refer the reader to the chapter [11] for a more detailed description of conflict-driven pseudo-Boolean solving (and to the handbook [8] for an in-depth treatment of SAT and related topics in general).

In contrast to MIP conflict analysis, the algorithm described above is not phrased in terms of the conflict graph, but focuses on the syntactic *resolution* method [9, 17, 16, 42] employed in CDCL conflict analysis and harnesses the observation by Hooker [29, 30] that
resolution can be understood as a cut rule and extended to 0–1 integer linear inequalities. The conflict-graph-based analysis in MIP does not operate on the reason constraints *R* as described above, but instead on disjunctive clauses extracted from these constraints. It is not hard to prove formally (appealing to [4, 14, 28]) that this incurs an exponential loss in reasoning power compared to performing derivations on the linear constraints themselves.

In practice, however, it seems fair to say that current pseudo-Boolean solvers do not quite deliver on this promise of exponential gains in performance. Although there are specific problem domains where PB solvers outperform even commercial MIP solvers [35, 45], evaluations over larger sets of benchmarks [5, 19, 20] have demonstrated that the open-source MIP solver SCIP [7] tends to be clearly more effective in solving pseudo-Boolean optimization problems, and is also quite competitive for decision problems. This is especially so for some decision problems that are in some sense close to LP-infeasibility – such problems are almost trivial for MIP solvers, but can be extremely challenging for pseudo-Boolean solvers [22].

## 1.3 Questions Studied in This Work and Our Contributions

Since mixed integer programming solvers and pseudo-Boolean solvers approach 0–1 integer linear problems from quite different angles, and seem to have complementary performance profiles, it is natural to ask whether techniques from one of the paradigms can be used to improve solvers based on the other paradigm.

Some MIP-inspired approaches have been integrated with success in SAT and PB solvers, perhaps most recently in [19], where the PB solver ROUNDINGSAT [23] makes careful use of the LP solver SOPLEX [7] to detect infeasibility of LP relaxations and generate cut constraints (though this paper also raises many questions that would seem to merit further study). However, in the other direction we are not aware of any work trying to harness state-of-the-art techniques from pseudo-Boolean solving to improve the performance of MIP solvers.

In this work, we consider how the clausal conflict analysis in MIP solvers can be replaced by pseudo-Boolean reasoning, focusing on 0–1 integer linear programs. A key difference between the clausal and pseudo-Boolean conflict analysis methods is that in the latter algorithm the linear reason constraint *R* propagating a variable assignment might need to be modified, or *reduced*, to another constraint *R*\* that propagates tightly also over the reals (which is already guaranteed to hold if *R* is a clausal constraint). Viewed from a MIP perspective, this *reduction* step deriving *R*\* from *R* can be seen to be an application of one of two specific cut rules, where saturation-based reduction as in [34] corresponds to coefficient tightening and division-based reduction as in [23] uses Chvátal-Gomory cuts.

This observation raises the question of whether more general cuts could also used to obtain other, and potentially more powerful, reduction methods for pseudo-Boolean conflict analysis. The answer turns out to be yes, and, in particular, we introduce a new reduction algorithm utilizing mixed integer rounding (MIR) cuts [27, 37]. A theoretical comparison of the MIR-based reduction rule with the reduction methods currently used in PB solvers show that MIR-based reduction dominates the division-based method that is considered to be state of the art in pseudo-Boolean solving, while saturation-based reduction and MIR-based reduction appear to be incomparable.

We implement pseudo-Boolean conflict analysis for 0–1 ILPs in the MIP solver SCIP, including all three reduction methods discussed above, and compare these different flavours of PB conflict analysis with each other as well as with clausal MIP conflict analysis on a large benchmark set consisting of pure 0–1 ILP instances from MIPLIB 2017. We find that the MIR-based pseudo-Boolean conflict analysis has the best performance, beating not only the
conflict analysis methods in the PB literature but also the standard clausal conflict analysis in SCIP. Interestingly, the new method is better measured not only in terms of number of nodes in the search tree, but also in terms of the number of instances solved, even though we only provide a proof-of-concept implementation lacking many of the optimizations that would be included in a full integration of this method into the SCIP codebase. Although our experimental data cannot provide conclusive evidence as to what causes this improved performance, we observe that the constraints learned from pseudo-Boolean conflict analysis seem more useful in that they take part more actively in propagations than constraints obtained by clausal conflict analysis.

## 1.4 Organization of This Paper

After reviewing preliminaries in Section 2, we give a detailed description of clausal and pseudo-Boolean conflict analysis for 0-1 integer linear programs in Section 3, including a discussion of the reduction methods found in the PB literature and our new version using mixed integer rounding cuts, and study how the different reduction rules compare in theory. In Section 4 we present our experimental results. We conclude the paper in Section 5 by summarizing our work and discussing direction for future research.

## 2 Preliminaries and Notation

Let $n \in \mathbb{Z}_{>0}$, and $\mathcal{N} := [1, \dots, n]$. We let $x_i$ denote Boolean (i.e., $\{0, 1\}$-valued) variables and $\ell_i$ denote literals, which can be either $x_i$ or its negation $\bar{x}_i = 1 - x_i$. A pseudo-Boolean constraint is a 0-1 integer linear inequality

$$ \sum_{i \in \mathcal{N}} a_i \ell_i \geq b, \qquad (1) $$

where we can assume without loss of generality that $a_i \in \mathbb{Z}_{\ge 0}$ for all $i \in \mathcal{N}$ and $b \in \mathbb{Z}_{\ge 0}$ (so-called normalized form). We can convert “$\le$”-constraints with 0-1 variables to “$\ge$”-constraints by multiplying the constraint by $-1$ and normalizing, i.e., replacing the variables by literals. Moreover, equalities “$=$” can be viewed as syntactic sugar for two opposing inequalities, which can also be transformed into normalized pseudo-Boolean format. In particular, every pure 0-1 integer linear program can be transformed to a normalized pseudo-Boolean representation. Note that in Section 3 we develop our theory and algorithms using normalized PB constraints for simplicity of exposition. However, in our actual implementation and experiments (described in Section 4), we directly operate on general linear constraints.

A (partial) assignment $\rho$ is a (partial) map from variables to $\{0, 1\}$, which is extended to literals by respecting the meaning of negation. We call a literal $\ell_i$ falsified or false if $\rho(\ell_i) = 0$ and satisfied or true if $\rho(\ell_i) = 1$. If $\rho$ is undefined for a literal, we call the literal unassigned or free. A constraint is satisfied under some partial assignment $\rho$ if the respective inequality holds, independently of which values the unassigned literals take, and is falsified if no assignment to the unassigned literals can make the inequality true.

The slack of a PB constraint $C: \sum_{i \in \mathcal{N}} a_i \ell_i \geq b$ under a partial assignment $\rho$ is defined as $\text{slack}(C, \rho) := \sum_{\{i \in \mathcal{N}: \rho(i) \neq 0\}} a_i - b$. With this definition, $C$ is falsified under $\rho$ if and only if $\text{slack}(C, \rho) < 0$. For example the constraint $C: 2\bar{x}_1 + 2x_2 + 3x_3 \geq 4$ is falsified under the partial assignment $\rho = \{x_1 = 1, x_2 = 0\}$ since $\text{slack}(C, \rho) = -1 < 0$. For a non-falsified constraint $C$ and an unassigned literal $\ell_i$ with coefficient $a_i$, the constraint propagates $\ell_i$ if and only if $\text{slack}(C, \rho) < a_i$. For instance, the same constraint $C: 2\bar{x}_1 + 2x_2 + 3x_3 \geq 4$
propagates both variables $x_2$ and $x_3$ to 1 under the partial assignment $\rho = \{x_1 = 1\}$ since slack($C, \rho$) = 1 is strictly smaller than the coefficients of each of the variables. A constraint propagates the assignment of a free variable tightly if the slack under the current partial assignment is 0. For any two pseudo-Boolean constraints $C$ and $C'$ and partial assignment $\rho$ it holds that the slack is subadditive, i.e., slack($C + C'$, $\rho$) $\le$ slack($C, \rho$) + slack($C'$, $\rho$). The decision level of a literal $\ell_i$ under a partial assignment $\rho$ is the number of decisions prior to the fixing of $\ell_i$. Note that the first fixing in every decision level is a decision literal.

# 3 Conflict Analysis Algorithms

For simplicity, in this section we present all algorithms in a pseudo-Boolean framework, where all coefficients and constants are integral, and the proofs of correctness that we provide also make crucial use of this fact. It is important to note that this is not the case in the actual implementation in SCIP, which operates with real-valued coefficients and constants. In fact, one of the challenges in implementing pseudo-Boolean conflict analysis in a MIP framework is that careful thought is required to rephrase the algorithms in such a way that they can deal with real-valued data but are still correct. Next, we describe the details of conflict analysis algorithms used in PB solvers and the different techniques that we consider in this paper.

## 3.1 Clausal Conflict Analysis

To explain the idea of conflict analysis, we first consider the case where all constraints are clauses. Conflict analysis begins at the stage where a conflict clause $C_{\text{confl}}$ is falsified by the current partial assignment $\rho$. Let $\ell_r$ be the literal in $C_{\text{confl}}$ that was last propagated to false, and let $C_{\text{reason}}$ be the reason clause in chronological order that is responsible for the propagation, i.e., we have $C_{\text{confl}} = C' \lor \ell_r$ and $C_{\text{reason}} = C'' \lor \ell_r$. Using the resolution rule, we can derive the so-called *resolvent* $C' \lor C''$ as a new learned clause $C_{\text{learn}}$.

Note that, even after removing $\ell_r$ from the partial assignment $\rho$, both $C'$ and $C''$ remain falsified: $C'$ because $C_{\text{confl}} = C' \lor \ell_r$ and $\ell_r$ were false, and $C''$ because $C_{\text{reason}} = C'' \lor \ell_r$ propagated. This is the key invariant of the algorithm: At any point during the algorithm the resolvent is falsified by the remaining partial assignment $\rho$.

Hence, we can replace the conflict clause by the resolvent and continue this process. At each step either a propagating literal is removed from $\rho$ or the learned clause is empty (at which point unsatisfiability is proven) or the last fixed literal is a decision literal. In the third case, we have reached a *first unique implication point* (FUIP) and conflict analysis terminates, with the final resolvent being the learned clause $C_{\text{learn}}$. With $C_{\text{learn}}$ added, propagation on the previous decision level will prevent the last infeasible decisions to happen as the search continues.

It is straightforward to apply this algorithm to problems with 0-1 linear constraints. Suppose $\sum_{i \in N} a_i \ell_i \ge b$ is the initial conflict constraint falsified under $\rho$, then $\bigvee_{i: a_i > 0 \land \rho(\ell_i) = 0} \ell_j$ can be used as initial conflict clause. Analogously, we can extract at each step a reason clause from the linear constraint that propagated the last literal and perform resolution. After terminating at an FUIP, the learned clause can be added as linear constraint to the solver.

## 3.2 PB Conflict Analysis

As in the clausal version, the main idea of PB conflict analysis is also to find a new constraint that explains the infeasibility of the current subproblem under a falsifying partial assignment. Algorithm 1 shows the base algorithm for all variants of PB conflict analysis considered in this
paper, using the first unique implication point (FUIP) learning scheme. It is initialized with a falsifying partial assignment $\rho$ and a conflicting constraint $C_{\text{confl}}$ under $\rho$. First, the learned conflict constraint $C_{\text{learn}}$ is set equal to the conflict constraint $C_{\text{confl}}$. In each iteration, we extract the latest literal $\ell_r$ from $\rho$. If the literal assignment was due to propagation of a constraint and the negated literal $\bar{\ell}_r$ occurs in $C_{\text{confl}}$, then we extract the reason constraint $C_{\text{reason}}$ that propagated $\ell_r$. In line 6 we “reduce” the reason constraint such that the resolvent of $C_{\text{learn}}$ and the reduced reason $C_{\text{reason}}$ (Line 7) that cancel the last literal $\ell_r$ is still falsified under the remaining partial assignment $\rho$. The conflict constraint is set to the resolvent and we continue until we reach an FUIP ($C_{\text{learn}}$ is asserting) or we prove $C_{\text{learn}}$ makes the problem infeasible. We have reached an FUIP if $C_{\text{learn}}$ would propagate some literal after removing at least all literal assignments in the current decision level from $\rho$. We have shown that the problem is infeasible if $C_{\text{learn}}$ is falsified under an empty partial assignment $\rho$. At this point, the learned constraint can be added to the constraint database of our problem to prevent the solver from exploring the same search space again.

<table><caption>**Algorithm 1** Pseudo-Boolean Conflict Analysis Algorithm.</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>Input</td><td>: conflict constraint $C_{\text{confl}}$, falsifying partial assignment $\rho$</td></tr><tr><td>Output</td><td>: learned conflict constraint $C_{\text{learn}}$</td></tr><tr><td>1</td><td>$C_{\text{learn}} \leftarrow C_{\text{confl}}$</td></tr><tr><td>2</td><td>while $C_{\text{learn}}$ not asserting and $C_{\text{learn}} \neq \perp$ do</td></tr><tr><td>3</td><td>&nbsp;&nbsp;&nbsp;&nbsp;$\ell_r \leftarrow$ literal last assigned on $\rho$</td></tr><tr><td>4</td><td>&nbsp;&nbsp;&nbsp;&nbsp;if $\ell_r$ propagated and $\bar{\ell}_r$ occurs in $C_{\text{learn}}$ then</td></tr><tr><td>5</td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$C_{\text{reason}} \leftarrow \text{reason}(\ell_r, \rho)$</td></tr><tr><td>6</td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$C_{\text{reason}} \leftarrow \text{reduce}(C_{\text{reason}}, C_{\text{learn}}, \ell_r, \rho)$</td></tr><tr><td>7</td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$C_{\text{learn}} \leftarrow \text{resolve}(C_{\text{learn}}, C_{\text{reason}}, \ell_r)$</td></tr><tr><td>8</td><td>&nbsp;&nbsp;&nbsp;&nbsp;$\rho \leftarrow \rho \setminus \{\ell_r\}$</td></tr><tr><td>9</td><td>return $C_{\text{learn}}$</td></tr></tbody></table>

The key invariant of Algorithm 1 is that in each iteration the resolvent $C_{\text{learn}}$ remains falsified. In the clausal version this holds even without the reduction step in line 6. However, for general linear constraints this is not the case, as shown by the next example.

**► Example 1.** Consider the two PB constraints $C_{\text{reason}} = x_1 + x_2 + 2x_3 \ge 2$ and $C_{\text{confl}} = x_1 + 2\bar{x}_3 + x_4 + x_5 \ge 3$ and the partial assignment $\rho = \{x_1 = 0, x_3 = 1\}$ where $x_1 = 0$ is a decision, and $x_3 = 1$ is propagated by $C_{\text{reason}}$. Under $\rho$ the constraint $C_{\text{confl}}$ is falsified. Applying generalized resolution to cancel $x_3$ yields the constraint $2x_1 + x_2 + x_4 + x_5 \ge 3$ which is not falsified under $\rho$.

In the following sections, we present three different reduction techniques for Algorithm 1 that operate directly on PB constraints. The main idea is to apply valid operations on the reason constraint to reduce the slack and ensure that the resolvent will have negative slack. The two main ingredients of the reduction techniques are *weakening* and *cutting planes* and are applied to the reason constraint until the invariant is fulfilled.

Weakening a literal in a PB constraint simply sets it to 1. For example, weakening a constraint $C : x_1 + x_2 + 2x_3 \ge 2$ on $x_1$ yields $x_2 + 2x_3 \ge 1$. Weakening is a valid operation since it simply adds a multiple of the valid bound constraint $\bar{x}_1 \ge 0$ to $C$. Note that weakening entails a loss of information. However, as we will see, it is a necessary operation to reduce the slack of the reason constraint. Note that whenever weakening is applied on non-falsified literal, it does not change the slack of the constraint. See Section 3.7 for more details on weakening.
27:8 Improving Conflict Analysis in MIP Solvers by Pseudo-Boolean Reasoning

Our main focus in this paper, however, is the second necessary ingredient of the reduction algorithm: cutting planes (cuts). Cuts are applied to the “weakened” version of the reason constraint in order to reduce its slack. We first present two well-documented cuts from existing literature, namely *Saturation* (Section 3.3) and *Division* (Section 3.4). Both ensure the reduction of the slack of the reason constraint to 0 at least after weakening all non-falsified literals in the original reason constraint. In Section 3.5, we introduce a new cut based on the *Mixed Integer Rounding (MIR)* procedure and prove that it has the same property. In Section 3.6 we show that the reduction algorithm using MIR always returns an equally strong or stronger reason constraint than the reduction using Division.

3.3 Saturation-based Reduction

First, we present the *Saturation* cut. Then, we provide details about the *Saturation-based Reduction* algorithm and demonstrate how the reduction ensures that the key invariant of conflict analysis holds.

► **Definition 2 (Saturation Cut).** Let $C: \sum_{i \in \mathcal{N}} a_i \ell_i \ge b$. The **Saturation Cut** of $C$ is given by the constraint

$$ \sum_{i \in \mathcal{N}} \min\{a_i, b\} \ell_i \ge b. $$

Saturation is a valid cut known as coefficient tightening cut in the MIP literature [10] and does not entail a loss of information. Algorithm 2 is used to reduce the reason constraint $C_{reason}$ before applying generalized resolution. Similar to the implementation in [12], in each iteration, the algorithm picks a non-falsified literal in the reason constraint different from the literal we are resolving on and weakens it. Then it applies the Saturation cut to the resulting constraint. The algorithm terminates when the slack of the resolvent becomes negative.

<table><caption>■ **Algorithm 2** Saturation-based Reduction Algorithm.</caption><tr><td>Input</td><td>: conflict constraint $C_{confl}$, reason constraint $C_{reason}$, literal to resolve $\ell_r$, partial assignment $\rho$</td></tr><tr><td>Output</td><td>: reduced reason $C_{reason}$</td></tr><tr><td>1</td><td>while slack((resolve($C_{reason}$, $C_{confl}$, $\ell_r$), $\rho$) &ge; 0 do</td></tr><tr><td>2</td><td>&nbsp;&nbsp;&nbsp;&nbsp;$\ell_j$ &leftarrow; non falsified literal in $C_{reason} \setminus \{\ell_r\}$</td></tr><tr><td>3</td><td>&nbsp;&nbsp;&nbsp;&nbsp;$C_{reason}$ &leftarrow; <b>weaken</b>($C_{reason}$, $\ell_j$)</td></tr><tr><td>4</td><td>&nbsp;&nbsp;&nbsp;&nbsp;$C_{reason}$ &leftarrow; <b>saturate</b>($C_{reason}$)</td></tr><tr><td>5</td><td>return $C_{reason}$</td></tr></table>

For completeness, we prove the following well-known fact that demonstrates that the slack of the reason constraint will be reduced to 0 at the latest after weakening the last non-falsified literal and applying the Saturation cut. Since the slack is subadditive, the resolvent's slack becomes negative and the resolvent is thus falsified.

► **Lemma 3.** Let $\rho$ be a partial assignment, and $C_{reason}: \sum_{i \in \mathcal{N}} a_i \ell_i \ge b$ a constraint propagating a literal $\ell_r$ to 1. Further, assume that slack($C_{reason}$, $\rho$) > 0. Then, after weakening all non-falsified literals in $C_{reason}$ (except for $\ell_r$) and applying Saturation on $C_{reason}$, the slack of the reduced reason constraint is 0.
**Proof.** First, we rewrite the constraint $C_{\text{reason}}$ as

$$ \sum_{j:\rho(j)=0} a_j \ell_j + \sum_{i \neq r:\rho(i) \neq 0} a_i \ell_i + a_r \ell_r \geq b. $$

Since slack($C_{\text{reason}}$, $\rho$) := $\sum_{i \neq r:\rho(i) \neq 0} a_i + a_r - b > 0$, it holds that

$$ a_r > b - \sum_{i \neq r:\rho(i) \neq 0} a_i. \qquad (2) $$

After weakening all literals from $\{i \neq r : \rho(i) \neq 0\}$ the constraint $C_{\text{reason}}$ becomes

$$ \sum_{j:\rho(j)=0} a_j \ell_j + a_r \ell_r \geq \tilde{b} := b - \sum_{i \neq r:\rho(i) \neq 0} a_i. \qquad (3) $$

Applying Saturation on (3) sets $a_r$ to $\tilde{b}$ because of (2). Therefore the slack of the reduced reason constraint becomes $\tilde{b} - \tilde{b} = 0$. ◀

## 3.4 Division-based Reduction

A very competitive alternative to Saturation in the reduction algorithm is based on Division cuts. Division is also a valid cut known as Chvátal-Gomory cut in the MIP literature [13].

**Definition 4 (Division Cut).** Let $C: \sum_{i \in \mathcal{N}} a_i \ell_i \geq b$. The Division Cut of $C$ with divisor $d \in \mathbb{Z}_{>0}$ is given by the constraint

$$ \sum_{i \in \mathcal{N}} \lfloor \frac{a_i}{d} \rfloor \ell_i \geq \lfloor \frac{b}{d} \rfloor. $$

To see why this procedure is valid, we can think of it as three steps: dividing by $d$ maintains the validity of the constraint; rounding up coefficients on the left-hand side relaxes the constraint and is hence valid; the validity of rounding up the right-hand side follows from the integrality of the left-hand side coefficients and literals.

In the Division-based reduction algorithm, the divisor $d$ used is the coefficient of the literal $\ell_r$ we are resolving on. As proven in [23], it suffices to weaken non-falsified variables with a coefficient that is not a multiple of $a_r$, i.e., from the index set $W := \{i \in \mathcal{N} : \rho(i) \neq 0 \text{ and } a_i \nmid a_r\}$. After weakening all literals in $W$ and applying Division on $C_{\text{reason}}$, the slack of the reduced reason constraint is 0, which for completeness we include in Lemma 6 below.

## 3.5 MIR-based Reduction

Next, we define a new cut for the reduction algorithm based on the *Mixed Integer Rounding* formula [37], which is a generalization of Gomory's mixed integer cuts [27].

**Definition 5 (Mixed Integer Rounding Cut).** Let $C: \sum_{i \in \mathcal{N}} a_i \ell_i \geq b$. The Mixed Integer Rounding (MIR) Cut of $C$ with divisor $d \in \mathbb{Z}_{>0}$ is given by the constraint

$$ \sum_{i \in I_1} \lfloor \frac{a_i}{d} \rfloor \ell_i + \sum_{i \in I_2} \left( \lfloor \frac{a_i}{d} \rfloor + \frac{f(a_i/d)}{f(b/d)} \right) \ell_i \geq \lfloor \frac{b}{d} \rfloor, \qquad (4) $$

where

$$ I_1 = \{i \in \mathcal{N} : f(a_i/d) \geq f(b/d) \text{ or } f(a_i/d) \in \mathbb{Z}\}, $$
27:10 Improving Conflict Analysis in MIP Solvers by Pseudo-Boolean Reasoning

$$
I_2 = \{i' \in \mathcal{N} : f(a_{i'}/d) < f(b/d) \text{ and } f(a_{i'}/d) \notin \mathbb{Z}\},
$$

and $f(\cdot) = -\lfloor \cdot \rfloor$. To obtain a normalized version of the MIR cut, we multiply both sides of
the constraint by $(b \mod d)$.

The proof that applying MIR to a constraint is a valid procedure can be found in [37]. Similar
to the Division-based reduction, it suffices to weaken non-falsified variables with a coefficient
that is not a multiple of $a_r$ before applying MIR in order to reduce the slack of the reason
constraint to at most 0. This is shown in the following lemma.

**Lemma 6.** Let $\rho$ be a partial assignment and $C_{\text{reason}} : \sum_{i \in \mathcal{N}} a_i \ell_i \ge b$ a constraint propagating a literal $\ell_r$ to 1. Then, after weakening all non-falsified literal in $W := \{i \in \mathcal{N} : \rho(i) \ne 0 \text{ and } a_r \nmid a_i\}$ and applying Division or MIR on $C_{\text{reason}}$ with $d = a_r$, the slack of the reduced reason is at most 0.

**Proof.** After weakening all literals in W, the constraint $C_{\text{reason}}$ becomes

$$
a_r \ell_r + \sum_{j \in \mathcal{N} \setminus W} a_j \ell_j \geq \tilde{b} := b - \sum_{i \in W} a_i. \quad (5)
$$

Its slack is

$$
\text{slack}(C_{\text{reason}}, \rho) = a_r + \sum_{\substack{j \in \mathcal{N} \setminus W: \\ \rho(j) \neq 0}} a_j - \tilde{b} = a_r + \sum_{\substack{j \in \mathcal{N}: \\ \rho(j) \neq 0, a_r | a_j}} a_j - \tilde{b}.
$$

Since weakening does not affect the slack, we have slack($C_{\text{reason}}, \rho) < a_r.

1. After applying the Division cut to (5) with $d = a_r$, the slack becomes

$$
\text{slack}(C_{\text{reason}}, \rho) = 1 + \sum_{\substack{j \in \mathcal{N}: \\ \rho(j) \neq 0, a_r | a_j}} \left[ \frac{a_j}{a_r} \right] - \left[ \frac{\tilde{b}}{a_r} \right] \le 1 + \sum_{\substack{j \in \mathcal{N}: \\ \rho(j) \neq 0, a_r | a_j}} \frac{a_j}{a_r} - \frac{\tilde{b}}{a_r} < \frac{a_r}{a_r} = 1. \quad (6)
$$

Because $C_{\text{reason}}$ contains only integer coefficients after applying the division rule, its slack is
integer; hence, it must be at most 0.

2. Applying the MIR cut to (5) with $d = a_r$ results in the same slack as in (6). This is because all left-hand side coefficients in the slack computation are divisible by $d$, hence they fall into the index set $I_1$ and are transformed the same way as by the Division cut. ◀

3.6 Dominance Relationships

In this section, we would like to discuss briefly known dominance relationships between the
different reduction techniques. The ultimate goal is to find a reduction technique that yields
the strongest possible reason constraint to use in the resolution step of conflict analysis. The
following lemma states the well-known fact that constraints from Saturation-based reduction
are always at least as strong as the resolvents created during clausal conflict analysis as
described in Section 3.1.

**Lemma 7.** Let $\rho$ be a partial assignment and $C_{\text{reason}} : \sum_{i \in \mathcal{N}} a_i \ell_i \ge b$ be a PB constraint which propagates literal $\ell_r$ to 1. Let $C'_{\text{reason}}$ and $C''_{\text{reason}}$ be the constraints obtained by clausal and Saturation-based reduction, respectively. Then $C'''_{\text{reason}}$ implies $C'_{\text{reason}}$.

**Proof.** Under the current partial assignment $\rho$, the disjunctive clause reason is given by
$\ell_r \lor_{j:\rho(\ell_j)=0} \ell_j$, which can be linearized as

$$
C'_{\text{reason}} : \ell_r + \sum_{j : \rho(\ell_j) = 0} \ell_j \geq 1.
$$
Now let W be the set of all non-falsified literals, except $l_r$. After weakening all literals in W and applying Saturation, we obtain the constraint

$$C''_{\text{reason}} : \min\{a_r, b - \sum_{i \in W} a_i\}l_r + \sum_{j: \rho(l_j)=0} \min\{a_j, b - \sum_{i \in W} a_i\}l_j \geq b - \sum_{i \in W} a_i.$$

As in the proof of Lemma 3, it holds that $\min\{a_r, b - \sum_{i \in W} a_i\} = b - \sum_{i \in W} a_i$. Now, after scaling $C''_{\text{reason}}$ by $b - \sum_{i \in W} a_i$ we see that $C''_{\text{reason}}$ has the same right-hand side as $C'_{\text{reason}}$, but smaller or equal coefficients on the left-hand side. ◀

In [26] the authors show that using Division instead of Saturation can be exponentially stronger, and that a single Saturation step can be simulated by an exponential number of Division steps.

The dominance of MIR cuts over Chvátal-Gomory cuts is a well-known fact in the MIP literature. The following lemma shows essentially the same result as in [15], but in the context of conflict analysis for pseudo-Boolean problems.

**Lemma 8.** Let $\rho$, $C_{\text{reason}}$, $l_r$ be given as in Lemma 7. Let $C'_{\text{reason}}$ and $C''_{\text{reason}}$ be the constraints obtained by Division-based and MIR-based reduction, respectively. Then $C''_{\text{reason}}$ implies $C'_{\text{reason}}$.

**Proof.** Let $C'_{\text{reason}}$, $C''_{\text{reason}}$ be constraints as in Definition 4 and 5, respectively, with divisor $d = a_r$. The constraints have the same right-hand side and the same coefficients for all literals $l_i$ with $i \in I_1$. For $i \in I_2$ the coefficient of literal $l_i$ in $C'_{\text{reason}}$ is given by $\lfloor \frac{a_i}{a_r} \rfloor$ and in $C''_{\text{reason}}$ by $\lfloor \frac{a_i}{a_r} \rfloor + \frac{f(a_i/a_r)}{f(b/a_r)}$. The coefficients of the literals $l_i$ in $C'_{\text{reason}}$ are always greater than or equal to the coefficients in $C''_{\text{reason}}$, since by definition of the set $I_2$ it holds that $f(a_i/a_r)/f(b/a_r) < 1$. Therefore $C''_{\text{reason}}$ implies $C'_{\text{reason}}$. ◀

As an example, consider the partial assignment $\rho = \{x_1 = 0, x_2 = 0, x_3 = 1\}$ and the constraint $C_{\text{reason}} : 2x_1 + 6x_2 + 10x_3 \geq 8$ which propagates variable $x_3$ to 1. Then the Division cut with divisor 10 is $x_1 + x_2 + x_3 \geq 1$. The MIR cut with the same divisor is $\frac{0.2}{0.8}x_1 + \frac{0.6}{0.8}x_2 + x_3 \geq 1$. Multiplying with 8 mod 10 = 8 gives the normalized MIR cut $2x_1 + 6x_2 + 8x_3 \geq 8$. The normalized MIR cut is stronger than the Division cut, which can be easily seen after scaling the Division cut by 8.

## 3.7 Practical Aspects of Weakening

While the evaluation of different weakening strategies is not the focus of this paper, we would like to discuss briefly some practical aspects of weakening literals. In our implementation we consider the following iterative weakening strategy: weaken free literals first followed by implied literals. We stop as soon as the resolvent is falsified under the remaining partial assignment. Intuitively, this order is motivated by the fact that free literals are not relevant for the propagation of literals in the reason constraint and do not affect the falsification of the conflict constraint.

However, the optimal order in which to weaken literals is not yet fully understood, and remains an open research question. Possible approaches include weakening literals in order of increasing or decreasing coefficient size. In [33] the authors conducted experiments with various weakening techniques, including partial weakening of literals and applying weakening on the conflict constraint, but the results did not yield a conclusive “best” weakening strategy.

A simple alternative is to weaken literals in a single sweep. For all three reduction algorithms, we can weaken the entire candidate set of literals as stated in Lemma 3 and Lemma 6 at once. Weakening literals all at once leads to a faster reduction algorithm
since repeated slack computations are avoided and only one cut is applied in each iteration. However, this may result in the constraint being less informative due to unnecessary weakening of literals.

4 Experiments

It is well known in the SAT and PB communities that efficient conflict-driven search requires substantial amounts of very careful engineering. In this first work, our focus has been on importing and adapting the pseudo-Boolean conflict analysis to a MIP setting – which is a nontrivial task in its own right – leaving further optimizations as future work.

All techniques from Section 3 have been implemented in the open source MIP solver
SCIP 8.0.3 [7] and we conducted extensive experiments to compare the different reduction
techniques in isolation. Obtaining accurate performance results for MIP solvers requires a
carefully designed experimental setup since even small changes to algorithms or the input
data can have a large impact on the behavior and the performance of the solver. This is
a well-known fact in the MIP literature known as performance variability [36]. To lessen
the effects of performance variability and obtain a fair comparison of the different reduction
techniques in the context of MIP solving, we use a fairly large and diverse testset of instances
and different permutations of each instance, see, e.g., [25]. Our experiments were carried out
on all pure 0–1 models from the MIPLIB 2017 collection [25]. After removing numerically
unstable models (with the tag “ numerics ”) our testset consists of 195 instances permuted by
5 different random seeds, giving a total of 975 measurements per run. For the remainder of
this paper, we will refer to the combination of a model and a permutation as an instance.
All experiments are conducted on a cluster with Intel Xeon Gold 6338 CPUs with a limit of
16GB of RAM.

It's worth noting that SCIP, along with its underlying LP solver, is based on floating-point arithmetic. Implementing a Pseudo-Boolean Optimization solver using a limited-precision LP-based branch-and-cut framework comes with some technical challenges which are discussed, e.g., in [5, 6]. From a theoretical standpoint, switching between reals and integers (rather than between limited and arbitrary precision) is straightforward:

All the algorithms presented in Section 3 can be naturally extended to the case of 0–1
constraints with coefficients that are real numbers instead of nonnegative integers. The
Chvátal-Gomory procedure, MIR cutting, and coefficient tightening algorithm were originally
designed for MIP with real coefficients.

However, in practice, floating-point arithmetic may cause numerical issues due to imprecise
representations of real numbers and cancellation effects. To mitigate the risk of numeric
instability, many components of SCIP, such as MIR-cut generation, utilize double-double
precision arithmetic [18], which could be also employed in conflict analysis. Currently, for
constraints generated in conflict analysis, we use the following standard techniques:

- We terminate conflict analysis if the coefficients of the constraints span too many orders of magnitude. Specifically, if the quotient of the largest to smallest coefficient is large (in our implementation, $10^6$), we stop conflict analysis.

- We remove variables from the conflict constraint if their coefficients are too small (in our implementation, $10^{-9}$), thereby relaxing the constraint slightly.

The latter threshold is a common default value for the zero tolerance in MIP solvers, and
the former is a common modeling recommendation for MIP.
**Table 1** Average percentage of true or unassigned literals that should be weakened to preserve the conflict analysis invariant. This experiment is conducted on the test set with 3 random seeds.

<table><thead><tr><td>Setting</td><td>avg(%) literals weakened</td></tr></thead><tbody><tr><td>Division</td><td>98.0</td></tr><tr><td>Saturation</td><td>99.7</td></tr><tr><td>MIR</td><td>97.3</td></tr></tbody></table>

## 4.1 Pre-Experiment: Weaken-All-At-Once vs. Weaken-Iteratively

As noted earlier, the weakening rule can be applied iteratively or in a single sweep. In preliminary experiments, we noticed that in almost all cases, most unassigned or true literals must be weakened to achieve the conflict analysis invariant that the resolved constraint has a negative slack. Table 1 summarizes this finding for different reduction techniques: Over all instances and all conflict analysis calls, an average between 97.3% (MIR) and 99.7% (Saturation) of all literals had to be weakened. Furthermore, for most instances both weakening variants did not lead to different execution paths.

In this case, weakening all literals at once avoids the overhead of iterative use of cuts and expensive slack computations. Consequently, we decided to always weaken all literals at once and apply the cut rule on the reason side only once for the remaining experiments presented in this paper.

## 4.2 Main Experiments: Comparing Different Reasoning Techniques

In the following, we compare all different reduction techniques from Section 3 to SCIP without any conflict analysis.

In our comparisons, we report for each technique the number of optimally solved instances, as well as the shifted geometric means of the number of processed nodes and the CPU time in seconds. The shifted geometric mean, a standard performance aggregator in the MIP literature, of the values $t_1, \dots, t_n$ is defined as

$$ \left( \prod_{i=1}^{n} (t_i + s) \right)^{1/n} - s, \qquad (7) $$

for some $s > 0$. We set the shift $s$ to 1 second for the CPU time and to 100 nodes for the number of nodes. Our base of comparison is SCIP without conflict analysis ("No Conflicts"). We report absolute values for the shifted means, and also quotients comparing them to our base setting. A factor below 1 means that a setting was faster (or needed less nodes), and a factor greater than 1 means that it was detrimental.

In Table 2 we report the results of our experiments. The table is split in four parts. We show results for "all" instances, as well as for three subsets of instances: (i) instances that are "affected" by conflict analysis, hence where the execution path of at least one setting differs from the others, (ii) "[100, limit]" instances, which take at least 100 seconds to solve to optimality or hit the time limit and (iii) "all-optimal", which are instances solved by all settings. Note that the number of nodes can only be fairly compared on the "all-optimal" subset, since the number of nodes when hitting a time limit is hard to interpret and hard to aggregate with the same statistics on instances that are solved to optimality.

The variant of SCIP with clausal conflict analysis is referred to as "Clausal-CA". For a fair comparison of the different strategies, we disabled the upgrading of constraints to specialized types, i.e., all generated conflicts are treated as linear constraints, and further
27:14 Improving Conflict Analysis in MIP Solvers by Pseudo-Boolean Reasoning

<table><caption>Table 2 Main results.</caption><thead><tr><th></th><th>Setting</th><th>solved</th><th>time(s)</th><th># nodes</th><th>time quot</th><th>nodes quot</th></tr></thead><tbody><tr><td rowspan="5">all(975)</td><td>No Conflicts</td><td>394</td><td>656.75</td><td>784</td><td>1.0</td><td>1.0</td></tr><tr><td>Clausal-CA</td><td>405</td><td>603.55</td><td>682</td><td>0.92</td><td>0.87</td></tr><tr><td>Division</td><td>419</td><td>601.4</td><td>683</td><td>0.92</td><td>0.87</td></tr><tr><td>MIR</td><td>420</td><td>599.37</td><td>677</td><td>0.91</td><td>0.86</td></tr><tr><td>Saturation</td><td>418</td><td>599.76</td><td>692</td><td>0.91</td><td>0.88</td></tr><tr><td rowspan="5">affected(295)</td><td>No Conflicts</td><td>259</td><td>160.46</td><td>1096</td><td>1.0</td><td>1.0</td></tr><tr><td>Clausal-CA</td><td>270</td><td>122.64</td><td>776</td><td>0.76</td><td>0.71</td></tr><tr><td>Division</td><td>284</td><td>119.24</td><td>707</td><td>0.74</td><td>0.65</td></tr><tr><td>MIR</td><td>285</td><td>118.29</td><td>700</td><td>0.74</td><td>0.64</td></tr><tr><td>Saturation</td><td>283</td><td>118.09</td><td>735</td><td>0.74</td><td>0.67</td></tr><tr><td rowspan="5">[100, limit](218)</td><td>No Conflicts</td><td>182</td><td>667.14</td><td>2056</td><td>1.0</td><td>1.0</td></tr><tr><td>Clausal-CA</td><td>193</td><td>486.45</td><td>1466</td><td>0.73</td><td>0.71</td></tr><tr><td>Division</td><td>207</td><td>486.23</td><td>1345</td><td>0.73</td><td>0.65</td></tr><tr><td>MIR</td><td>208</td><td>485.26</td><td>1336</td><td>0.73</td><td>0.65</td></tr><tr><td>Saturation</td><td>206</td><td>491.98</td><td>1428</td><td>0.74</td><td>0.69</td></tr><tr><td rowspan="5">all-optimal(374)</td><td>No Conflicts</td><td>374</td><td>46.16</td><td>320</td><td>1.0</td><td>1.0</td></tr><tr><td>Clausal-CA</td><td>374</td><td>40.58</td><td>259</td><td>0.88</td><td>0.81</td></tr><tr><td>Division</td><td>374</td><td>40.75</td><td>244</td><td>0.88</td><td>0.77</td></tr><tr><td>MIR</td><td>374</td><td>40.40</td><td>241</td><td>0.88</td><td>0.75</td></tr><tr><td>Saturation</td><td>374</td><td>40.24</td><td>246</td><td>0.87</td><td>0.77</td></tr></tbody></table>

only generated one conflict per call. Conversely, we accept PB reasoning conflicts only if the number of nonzeros is less than 15% of the original problem variables, as in the default clausal implementation in SCIP. Our preliminary experiments confirmed that in our implementation, it is indeed detrimental to accept too-long conflicts. We did, however, add a fallback strategy, of applying weakening on the conflict constraint if the constraints are too long. This happens for about 9% of the instances.

We observe that all conflict analysis variants solved more instances than SCIP without conflict analysis, needed significantly less nodes on the all-optimal set, and less time on all four instance sets. Note that on average, the time spent in conflict analysis is only about 0.1% of the total run time. The three PB conflict analysis variants could solve more instances than the clausal variant, and needed significantly less nodes. The difference in time was less pronounced.

The performance of the PB conflict analysis variants is quite similar in all three cases. Nevertheless, MIR-based reduction could solve the most instances and needed the least nodes on the all-optimal set. When looking at the seemingly identical time-wise performance in more detail, it turns out that MIR also slightly improves on the other settings in this measure. There are 104 instances for which the path differs between Saturation-based resolution and MIR-based resolution and MIR was on average 1.1% faster on those. There are 86 instances for which the path differs between Division-based resolution and MIR-based resolution and MIR was on average 3.6% faster on those. Consequently, we decided to concentrate on MIR-based resolution for our next statistic.

Ultimately, the purpose of conflict constraints is to restrict the future search space by propagating literal assignments and pruning the search tree. Hence we analyzed how many conflicts each of the methods generates in shifted geometric mean, how large these conflicts are on average, and how many of them lead to propagations down the road. Table 3 shows
<table><caption>**Table 3** Shifted geometric mean of number of conflicts, average percentage of conflict constraints that propagate at least once and average length of learned conflicts.</caption><thead><tr><td>Setting</td><td>mean # conflicts</td><td>avg % prop. conflicts</td><td>avg # literals</td></tr></thead><tbody><tr><td>Clausal-CA</td><td>290.77</td><td>34.54</td><td>82.45</td></tr><tr><td>MIR</td><td>169.61</td><td>58.54</td><td>80.20</td></tr></tbody></table>

the results on the set of all instances that have a search tree of at least 100 nodes (to get a decent chance of conflict generation and propagation) and for which at least one conflict was generated with one of the methods. We consider only instances where the two settings have the same execution path. We observe that our MIR-based conflict analysis generated about a third less conflicts, but at the same time, they are much more likely to propagate: for the classic clausal conflict analysis of SCIP, about a third of the generated conflicts are used for propagation later on, while for our MIR-based variant, slightly more than half (58.54%) of all conflicts propagate at least once. At the same time, MIR-based conflicts are about the same size as clausal conflicts.

At first glance, this might appear as a contradiction, given that, as a rule of thumb, shorter conflicts tend to propagate more often and one might expect similar-sized conflicts to be similarly likely to propagate. Note, however, that the conflicts are of a quite different nature in the two cases. On the one hand, clausal conflicts are always logic clauses that only propagate when all but one literal are assigned. On the other hand, MIR-based conflicts are general pseudo-Boolean constraints, which might propagate some assignments (of literals with large coefficients) even when a majority of literals are still unassigned. This goes nicely together with the above observation that the reduction in the number of nodes is more pronounced than the reduction in runtime. As a final remark, integrating PB conflict analysis in a production-grade MIP solver would require substantially more work, but should also be expected to provide substantial further improvements measured in wallclock time.

# 5 Conclusion

In this work, we study how to integrate pseudo-Boolean conflict analysis for 0-1 integer linear programs into a MIP solving framework. In contrast to standard MIP conflict analysis, the pseudo-Boolean method operates directly on the linear constraints, rather than on clauses extracted from these constraints, and this makes it exponentially stronger in terms of reasoning power. Viewing PB conflict analysis from a MIP perspective is also helpful since it provides a view of the algorithm as a sequence of linear combinations and cuts, and we use this to strengthen the pseudo-Boolean conflict analysis further by developing a new conflict analysis method using the powerful mixed integer rounding (MIR) cuts.

We have made a first proof-of-concept implementation of our new pseudo-Boolean conflict analysis method, as well as methods from the PB literature based on saturation [34] and division [23], in the open-source MIP solver SCIP, and have run experiments on 0-1 ILP instances from MIPLIB 2017 comparing the different methods with each other and with standard clause-based MIP conflict analysis. We find that solving 0-1 ILPs with MIR-based pseudo-Boolean conflict analysis performs better than other methods, not only in the sense that it reduces the size of the search tree, but also in that our implementation can beat the highly optimized MIP conflict analysis currently used in SCIP in terms of actual running time. In our opinion, this demonstrates convincingly that pseudo-Boolean conflict analysis in MIP is a research direction that should be worth pursuing further, and that similar proof-of-concept studies could also be relevant to investigate for other combinatorial solving paradigms such as constraint programming.
27:16 Improving Conflict Analysis in MIP Solvers by Pseudo-Boolean Reasoning

As already noted above, an obvious direction of future work is to provide a more carefully
engineered version of pseudo-Boolean conflict analysis that could deliver more fully on the
potential for improved performance identified by our experiments. In addition to optimizing
the existing code, however, it would be valuable to develop a better understanding of how
and why the conflict analysis works and of ways in which the reasoning could be improved.

Pseudo-Boolean conflict analysis alternates between weakening constraints (to eliminate seemingly less relevant variables) and strenghtening them by applying cut rules (to get tighter propagation on the variables that remain). The interplay between these two operations is quite poorly understood even for pseudo-Boolean solvers, and so both PB solvers and MIP solvers could gain from a careful study of how to strike the right balance. Since PB conflict analysis can be performed with several different reduction methods, and since different reduction methods can be employed independently in consecutive steps in one and the same conflict analysis, it would also be good to be able to assess the quality of constraints derived during conflict analysis, so as to select the most promising candidate at each step to pass on to the next step in the conflict analysis.

Arguably the most interesting research question, though, is whether pseudo-Boolean
conflict analysis could be extended beyond 0-1 ILPs to 0-1 mixed linear problems, and/or
to general integer linear programs. It is worth noting that the latter has been attempted
in [32, 41], but so far with quite limited success. It is clear that the algorithms presented
in this paper cannot work for 0-1 mixed LPs or general ILPs if generalized in the obvious,
naive way, and so additional, new ideas will be needed.

References

1 Tobias Achterberg. Conflict analysis in mixed integer programming. Discrete Optimization, 4(1):4–20, March 2007.

2 Tobias Achterberg. *Constraint Integer Programming*. PhD thesis, Technische Universität Berlin, 2007. Available at https://opus4.kobv.de/opus4-zib/files/1112/Achterberg_C constraint-IntegerProgramming.pdf.

3 Roberto J. Bayardo Jr. and Robert Schrag. Using CSP look-back techniques to solve real-world SAT instances. In *Proceedings of the 14th National Conference on Artificial Intelligence (AAAI '97)*, pages 203–208, July 1997.

4 Paul Beame, Henry Kautz, and Ashish Sabharwal. Towards understanding and harnessing the potential of clause learning. *Journal of Artificial Intelligence Research*, 22:319–351, December 2004. Preliminary version in *IJCAI '03*.

5 Timo Berthold, Stefan Heinz, and Marc Pfetsch. Solving Pseudo-Boolean problems with SCIP, 2008.

6 Timo Berthold, Stefan Heinz, and Marc E Pfetsch. Nonlinear Pseudo-Boolean optimization: relaxation or propagation? In *Theory and Applications of Satisfiability Testing-SAT 2009: 12th International Conference, SAT 2009, Swansea, UK, June 30-July 3, 2009*. Proceedings 12, pages 441–446. Springer, 2009.

7 Ksenia Bestuzheva, Mathieu Besançon, Wei-Kun Chen, Antonia Chmiela, Tim Donkiewicz, Jasper van Doornmalen, Leon Eifler, Oliver Gaul, Gerald Gamrath, Ambros Gleixner, et al. The SCIP Optimization Suite 8.0. *arXiv preprint arXiv:2112.08872*, 2021.

8 Armin Biere, Marijn J. H. Heule, Hans van Maaren, and Toby Walsh, editors. *Handbook of Satisfiability*, volume 336 of *Frontiers in Artificial Intelligence and Applications*. IOS Press, 2nd edition, February 2021.

9 Archie Blake. *Canonical Expressions in Boolean Algebra*. PhD thesis, University of Chicago, 1937.

10 AL Brearley, Gautam Mitra, and H Paul Williams. Analysis of mathematical programming problems prior to applying the simplex algorithm. *Mathematical programming*, 8:54–83, 1975.
11 Samuel R. Buss and Jakob Nordström. Proof complexity and SAT solving. In Armin Biere, Marijn J. H. Heule, Hans van Maaren, and Toby Walsh, editors, *Handbook of Satisfiability, volume 336 of Frontiers in Artificial Intelligence and Applications*, chapter 7, pages 233–350. IOS Press, 2nd edition, February 2021.

12 Donald Chai and Andreas Kuehlmann. A fast pseudo-Boolean constraint solver. *IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems*, 24(3):305–317, March 2005. Preliminary version in *DAC '03*.

13 Vasek Chvátal. Edmonds polytopes and a hierarchy of combinatorial problems. *Discrete mathematics*, 4(4):305–337, 1973.

14 William Cook, Collette Rene Coullard, and György Turán. On the complexity of cutting-plane proofs. *Discrete Applied Mathematics*, 18(1):25–38, November 1987.

15 Gérard Cornuéjols and Yanjun Li. Elementary closures for integer programs. *Operations Research Letters*, 28(1):1–8, 2001. doi:10.1016/S0167-6377(00)00067-5.

16 Martin Davis, George Logemann, and Donald Loveland. A machine program for theorem proving. *Communications of the ACM*, 5(7):394–397, July 1962.

17 Martin Davis and Hilary Putnam. A computing procedure for quantification theory. *Journal of the ACM*, 7(3):201–215, 1960.

18 Theodorus Jozef Dekker. A floating-point technique for extending the available precision. *Numerische Mathematik*, 18(3):224–242, 1971.

19 Jo Devriendt, Ambros Gleixner, and Jakob Nordström. Learn to relax: Integrating 0-1 integer linear programming with pseudo-Boolean conflict-driven search. *Constraints*, 26(1-4):26–55, October 2021. Preliminary version in *CPAIOR '20*.

20 Jo Devriendt, Stephan Gocht, Emir Demirović, Jakob Nordström, and Peter Stuckey. Cutting to the core of pseudo-Boolean optimization: Combining core-guided search with cutting planes reasoning. In *Proceedings of the 35th AAAI Conference on Artificial Intelligence (AAAI '21)*, pages 3750–3758, February 2021.

21 Niklas Eén and Niklas Sörensson. Translating pseudo-Boolean constraints into SAT. *Journal on Satisfiability, Boolean Modeling and Computation*, 2(1-4):1–26, March 2006.

22 Jan Elffers, Jesús Giráldez-Cru, Jakob Nordström, and Marc Vinyals. Using combinatorial benchmarks to probe the reasoning power of pseudo-Boolean solvers. In *Proceedings of the 21st International Conference on Theory and Applications of Satisfiability Testing (SAT '18)*, volume 10929 of *Lecture Notes in Computer Science*, pages 75–93. Springer, July 2018.

23 Jan Elffers and Jakob Nordström. Divide and conquer: Towards faster pseudo-Boolean solving. In *Proceedings of the 27th International Joint Conference on Artificial Intelligence (IJCAI '18)*, pages 1291–1299, July 2018.

24 Matthew L Ginsberg. Dynamic backtracking. *Journal of artificial intelligence research*, 1:25–46, 1993.

25 Ambros Gleixner, Gregor Hendel, Gerald Gamrath, Tobias Achterberg, Michael Bas-tubebe, Timo Berthold, Philipp M. Christophel, Kati Jarck, Thorsten Koch, Jeff Linderoth, Marco Lübbecke, Hans D. Mittelmann, Derya Ozyurt, Ted K. Ralphs, Domenico Sal-vagnin, and Yuji Shinano. MIPLIB 2017: Data-Driven Compilation of the 6th Mixed-Integer Programming Library. *Mathematical Programming Computation*, 13:443–490, 2021. doi:10.1007/s12532-020-00194-3.

26 Stephan Gocht, Jakob Nordström, and Amir Yehudayoff. On division versus saturation in pseudo-Boolean solving. In *Proceedings of the 28th International Joint Conference on Artificial Intelligence (IJCAI '19)*, pages 1711–1718, August 2019.

27 Ralph E. Gomory. An algorithm for the mixed integer problem. Technical Report P-1885, The RAND Corporation, June 1960.

28 Armin Haken. The intractability of resolution. *Theoretical Computer Science*, 39(2-3):297–308, August 1985.

29 John N. Hooker. Generalized resolution and cutting planes. *Annals of Operations Research*, 12(1):217–239, December 1988.
27:18 Improving Conflict Analysis in MIP Solvers by Pseudo-Boolean Reasoning

30 John N. Hooker. Generalized resolution for 0-1 linear inequalities. *Annals of Mathematics and Artificial Intelligence*, 6(1):271-286, March 1992.

31 Yuejun Jiang, Thomas Richards, and Barry Richards. Nogood backmarking with min-conflict repair in constraint satisfaction and optimization. In *Principles and Practice of Constraint Programming: Second International Workshop, PPCP'94 Rosario, Orcas Island, WA, USA, May 2-4, 1994 Proceedings*, pages 21-39. Springer, 1994.

32 Dejan Jovanovic and Leonardo de Moura. Cutting to the chase solving linear integer arithmetic. *Journal of Automated Reasoning*, 51(1):79-108, June 2013. Preliminary version in CADE-23.

33 Daniel Le Berre, Pierre Marquis, and Romain Wallon. On weakening strategies for PB solvers. In *Proceedings of the 23rd International Conference on Theory and Applications of Satisfiability Testing (SAT '20)*, volume 12178 of *Lecture Notes in Computer Science*, pages 322-331. Springer, July 2020.

34 Daniel Le Berre and Anne Parrain. The Sat4j library, release 2.2. *Journal on Satisfiability, Boolean Modeling and Computation*, 7:59-64, July 2010.

35 Vincent Liew, Paul Beame, Jo Devriendt, Jan Elffers, and Jakob Nordström. Verifying properties of bit-vector multiplication using cutting planes reasoning. In *Proceedings of the 20th Conference on Formal Methods in Computer-Aided Design (FMCAD '20)*, pages 194-204, September 2020.

36 Andrea Lodi and Andrea Tramontani. Performance variability in mixed-integer programming. In *Theory driven by influential applications*, pages 1-12. INFORMS, 2013.

37 Hugues Marchand and Laurence A Wolsey. Aggregation and mixed integer rounding to solve mips. *Operations research*, 49(3):363-371, 2001.

38 João P. Marques-Silva and Karem A. Sakallah. GRASP: A search algorithm for propositional satisfiability. *IEEE Transactions on Computers*, 48(5):506-521, May 1999. Preliminary version in *ICCAD '96*.

39 Ruben Martins, Vasco M. Manquinho, and Inês Lynce. Open-WBO: A modular MaxSAT solver. In *Proceedings of the 17th International Conference on Theory and Applications of Satisfiability Testing (SAT '14)*, volume 8561 of *Lecture Notes in Computer Science*, pages 438-445. Springer, July 2014.

40 Matthew W. Moskewicz, Conor F. Madigan, Ying Zhao, Lintao Zhang, and Sharad Malik. Chaff: Engineering an efficient SAT solver. In *Proceedings of the 38th Design Automation Conference (DAC '01)*, pages 530-535, June 2001.

41 Robert Nieuwenhuis. The IntSat method for integer linear programming. In *Proceedings of the 20th International Conference on Principles and Practice of Constraint Programming (CP '14)*, volume 8656 of *Lecture Notes in Computer Science*, pages 574-589. Springer, September 2014.

42 John Alan Robinson. A machine-oriented logic based on the resolution principle. *Journal of the ACM*, 12(1):23-41, January 1965.

43 Masahiko Sakai and Hidetomo Nabeshima. Construction of an ROBDD for a PB-constraint in band form and related techniques for PB-solvers. *IEICE Transactions on Information and Systems*, 98-D(6):1121-1127, June 2015.

44 Tuomas Sandholm and Robert Shields. Nogood learning for mixed integer programming. In *Workshop on Hybrid Methods and Branching Rules in Combinatorial Optimization, Montréal*, volume 20, pages 21-22, 2006.

45 Buser Say, Jo Devriendt, Jakob Nordström, and Peter Stuckey. Theoretical and experimental results for planning with learned binarized neural network transition models. In *Proceedings of the 26th International Conference on Principles and Practice of Constraint Programming (CP '20)*, volume 12333 of *Lecture Notes in Computer Science*, pages 917-934. Springer, September 2020.

46 Hossein M. Sheini and Karem A. Sakallah. Pueblo: A hybrid pseudo-Boolean SAT solver. *Journal on Satisfiability, Boolean Modeling and Computation*, 2(1-4):165-189, March 2006. Preliminary version in DATE '05.
47 Richard M Stallman and Gerald J Sussman. Forward reasoning and dependency-directed backtracking in a system for computer-aided circuit analysis. *Artificial intelligence*, 9(2):135-196, 1977.

48 Jakob Witzig, Timo Berthold, and Stefan Heinz. Computational aspects of infeasibility analysis in mixed integer programming. Technical Report 19-54, ZIB, Takustr. 7, 14195 Berlin, 2019.
<figure><img src="image_276.png" /></figure>
# Using Canonical Codes to Efficiently Solve the Benzenoid Generation Problem with Constraint Programming

Xiao Peng

Univ Lyon, INSA Lyon, Inria, CITI, EA3720, 69621 Villeurbanne, France

Christine Solnon

Univ Lyon, INSA Lyon, Inria, CITI, EA3720, 69621 Villeurbanne, France

**Abstract**

The Benzenoid Generation Problem (BGP) aims at generating all benzenoid molecules that satisfy some given properties. This problem has important applications in chemistry, and Carissan et al (2021) have shown us that Constraint Programming (CP) is well suited for modelling this problem because properties defined by chemists are easy to express by means of constraints. Benzenoids are described by hexagon graphs and a key point for an efficient enumeration of these graphs is to be invariant to rotations and symmetries. In this paper, we introduce canonical codes that uniquely characterise hexagon graphs while being invariant to rotations and symmetries. We show that these codes may be defined by means of constraints. We also introduce a global constraint for ensuring that codes are canonical, and a global constraint for ensuring that a pattern is included in a code. We experimentally compare our new CP model with the CP-based approach of Carissan et al (2021), and we show that it has better scale-up properties.

**2012 ACM Subject Classification** Mathematics of computing → Graph enumeration; Computing methodologies → Artificial intelligence

**Keywords and phrases** Benzenoid Generation Problem, Canonical Code, Hexagon Graph

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

**Acknowledgements** We want to thank authors of [3] who helped us reproduce their results.

## 1 Introduction

Benzenoids are hydrocarbon molecules whose carbon atoms are forming cycles of size 6, i.e., hexagons. An important chemistry problem concerns the generation of all benzenoids that satisfy some given properties [5]. This problem is called the Benzenoid Generation Problem (BGP) in [2]. As benzenoids are regular tillings with hexagonal faces, they may be represented by hexagon graphs such that a vertex is associated with every hexagonal face and an edge with every pair of vertices corresponding to adjacent faces [1, 2]. For example, we display in Figure 1 a benzenoid composed of five hexagonal faces and its associated hexagon graph. These hexagon graphs are always connected.

Different benzenoids may be represented with isomorphic hexagon graphs. For example, let us consider the hexagon graph displayed in Figure 1. The subgraph induced by *c*, *d*, and *e* is isomorphic to the subgraph induced by *a*, *c*, and *d* whereas their associated molecules are different because hexagons *c*, *d*, and *e* are not aligned whereas hexagons *a*, *c*, and *d* are aligned. To overcome this problem, we take into account edge directions, considering the six directions defined in Figure 1. In this case, the subgraph induced by *c*, *d*, and *e* is no longer isomorphic to the subgraph induced by *a*, *c*, and *d* because edges (*c*, *d*) and (*d*, *e*) have different directions (0 and 1), whereas edges (*a*, *c*) and (*c*, *d*) have the same direction (0).

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

© Xiao Peng and Christine Solnon;
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. 28; pp. 28:1–28:17

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

LIPICS

Leibniz International Proceedings in Informatics
Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl Publishing, Germany
<figure><img src="image_279.png" /><figcaption>■ **Figure 1** Example of benzenoid *B* and its associated hexagon graph *G*, and symmetrical graph $\bar{G}$. Right: Edge directions from a vertex *v* to each of its 6 possible neighbours.</figcaption></figure>

To solve the BGP, we have to enumerate hexagon graphs that satisfy some given properties and these properties are usually easily expressed by means of constraints (such as, for example, the absence of cliques of order three, or the inclusion of some patterns). Two main approaches may be used:

- The dedicated approach of [1] uses canonical codes to represent hexagon graphs. These codes are invariant to rotations and symmetries and they are used to decide whether two graphs are isomorphic or not in linear time. However, these codes cannot represent benzenoids with holes (called coronoids). Furthermore, this approach is not declarative and it does not allow one to easily add constraints on the graphs to be enumerated.

- The Constraint Programming (CP)-based approach of [2, 3] basically searches for all connected subgraphs within an initial hexagon graph, and it is implemented in Choco [11] using graph variables. Using CP allows one to easily add constraints and this approach is able to generate all kinds of benzenoids, including coronoids. However, it is less efficient than the dedicated approach.

In this paper, we introduce an approach which is both efficient and declarative: like [1], it is based on canonical codes but these canonical codes can represent all benzenoids (including coronoids); like [2, 3], it is based on CP so that one may easily add constraints to specify structural properties. In Section 2, we introduce our new canonical code and study some of its properties that are used to efficiently generate canonical codes with CP. In Section 3, we extend canonical codes to the case where benzenoids are constrained to contain some given patterns. In Section 4, we introduce CP models for solving BGPs. In Section 5, we report experimental results and compare our approach with the CP-based approach of [2, 3].

**Notations**

Given two integer values *i* and *j*, we note [*i*, *j*] the set of all integer values ranging from *i* to *j*. Given a hexagon graph *G*, we note $\bar{G}$ the symmetrical graph obtained by mirroring *G* with respect to the *x*-axis, as displayed in Figure 1. Given a hexagon graph *G* = (*V*, *E*) and a subset of vertices *S* ⊆ *V*, we note $G_{↓S}$ the subgraph of *G* induced by *S*, i.e., $G_{↓S} = (S, E \cap S \times S)$.

2 Representation of Hexagon Graphs with Canonical Codes

A key point for an efficient enumeration of hexagon graphs is to be invariant to rotations and symmetries. For example, the graph displayed in Figure 1 is isomorphic to any other graph obtained by rotating it of $k \times 60^\circ$ with $k \in \mathbb{N}$ and/or mirroring it with respect to the *x*-axis. In [3], constraints are added to break these symmetries and avoid generating several times the same graph. In this paper, we propose to use another approach based on canonical codes, i.e., integer sequences that uniquely characterise graphs and that are invariant to rotations and symmetries. We introduce our canonical code in Section 2.1. We study its properties in Section 2.2. We compare it with related work in Section 2.3.
■ **Algorithm 1** BFS<sub>G</sub>(v<sub>0</sub>, v<sub>1</sub>).

**Input:** A hexagon graph $G$ and an initial edge $(v_0, v_1)$ of $G$
**Output:** The code associated with a BFS of $G$ started from $(v_0, v_1)$
1 &nbsp;&nbsp;rotate $G$ so that the direction of edge $(v_0, v_1)$ is equal to 0
2 &nbsp;&nbsp;for each vertex $v_i$ of $G$ do initialise $num[v_i]$ to $-1$;
3 &nbsp;&nbsp;set $num[v_0]$ to 0 and initialise a counter $c$ to 1
4 &nbsp;&nbsp;let $q$ be an empty FIFO queue; add $v_0$ in $q$
5 &nbsp;&nbsp;initialise code to an empty sequence
6 &nbsp;&nbsp;while $q$ is not empty do
7 &nbsp;&nbsp;&nbsp;&nbsp;remove from $q$ its oldest vertex $v_i$
8 &nbsp;&nbsp;&nbsp;&nbsp;for each edge $(v_i, v_j)$ of $G$ taken by order of increasing direction do
9 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;if $num[v_j] < 0$ then
10 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;add $v_j$ in $q$
11 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;set $num[v_j]$ to $c$ and increment $c$
12 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;add $num[v_i]$ and the direction of edge $(v_i, v_j)$ at the end of code
13 &nbsp;&nbsp;return code

<figure><img src="image_280.png" /><figcaption>■ **Figure 2** Examples of runs of BFS<sub>G</sub>(v<sub>0</sub>, v<sub>1</sub>) when G is the graph of Figure 1 and (v<sub>0</sub>, v<sub>1</sub>) is equal to (a, c) for Ex1, (c, a) for Ex2, (c, b) for Ex3, and (e, d) for Ex4. For each vertex v<sub>i</sub> ∈ {a, b, c, d, e}, we display v<sub>i</sub> : num[v<sub>i</sub>] in circles. Each edge (v<sub>i</sub>, v<sub>j</sub>) used lines 10-12 is displayed in bold and its direction is displayed on top of it (the initial edge (v<sub>0</sub>, v<sub>1</sub>) is displayed in green).</figcaption></figure>

## 2.1 Definition of Canonical Codes

Given a hexagon graph $G$ with $n + 1$ vertices, a code is a sequence of $n$ couples of integer values which is associated with a Breadth First Search (BFS) of $G$ starting from a given initial edge $(v_0, v_1)$, as described in Algorithm 1. The initial edge is used to define the orientation of the graph (line 1): we always consider the orientation such that $v_1$ is on the right of $v_0$. This allows us to be invariant to rotations.

► **Example 1.** We display in Figure 2 different runs of BFS on the hexagon graph of Figure 1, starting from different edges. In Ex1, we start from edge (a, c) which already has direction 0 so that we do not have to rotate the graph. When starting from edge (c, a) (resp. (c, b) and (e, d)), we have to rotate the graph of 180° (resp. 240° and 60°).

BFS assigns a different number $num[v_i]$ to each vertex $v_i$ corresponding to the order vertices are discovered: $v_0$ is numbered with 0 (line 3) and when a vertex $v_j$ is reached for the first time it is numbered with the next non-assigned value $c$ (line 11). Each time the search discovers a new vertex $v_j$ using an edge $(v_i, v_j)$, $num[v_i]$ and $d(v_i, v_j)$ are added at the end of the code, where $d(v_i, v_j)$ is the direction of edge $(v_i, v_j)$ (line 12).
<figure><img src="image_281.png" /><figcaption>■ **Figure 3** Run of BFS$_{\bar{G}}(c, a)$ where $\bar{G}$ is the graph symmetrical to the graph $G$ of Figure 1.</figcaption></figure>

**Example 2.** In Figure 2 (Ex1), we illustrate a run of BFS$_G(a, c)$ when $G$ is the graph of Figure 1. In this case, the code is 00 01 10 31: the first couple is 00 because vertex 1 (c) has been reached from vertex 0 (a) and $d(a, c) = 0$; the second couple is 01 because vertex 2 (b) has been reached from vertex 0 (a) and $d(a, b) = 1$; the third couple is 10 because vertex 3 (d) has been reached from vertex 1 (c) and $d(c, d) = 1$; the fourth couple is 31 because vertex 4 (e) has been reached from vertex 3 (d) and $d(d, e) = 1$.

Note that the loop lines 8-12 considers edges outgoing from $v_i$ by order of increasing directions. This ensures that we always compute the same code provided that the direction of the first edge ($v_0, v_1$) is fixed to 0.

Given a code $p_1d_1 \ p_2d_2 \ \dots \ p_n d_n$ computed from a graph with $n + 1$ vertices, we can draw this graph, starting from vertex number 0: for each $i \in [1, n]$, $p_i$ gives the number of the predecessor of vertex number $i$ and $d_i$ gives the direction of the edge $(p_i, i)$. Once all vertices have been drawn, missing edges can be added as every vertex is connected to all its neighbours. For example, in Ex1, we add an edge between vertices 1 and 2 because they have neighbour positions.

There exist different possible codes for a given graph. Each code corresponds to a different BFS starting from a different initial edge ($v_0, v_1$). We define a total order on the set of all possible codes that may be associated with a given graph by considering a lexicographic order. Among all the possible codes for a graph, the smallest one according to this order is called the canonical code of this graph and it is unique.

So far, our canonical code is invariant to rotations but not to symmetries. To become invariant to symmetries, we must also compute all codes starting from all possible edges while considering the symmetrical graph $\bar{G}$. This leads us to the following definition.

**Definition 3 (Canonical code $cc(G)$).** The canonical code of a hexagon graph $G = (V, E)$ is: $cc(G) = \min\{BFS_G(v_0, v_1), BFS_{\bar{G}}(v_0, v_1) : (v_0, v_1) \in E\}$ when considering a lexicographic order to compare codes.

**Example 4.** The smallest code that may be computed for the graph $G$ of Figure 1 is BFS$_G(c, b)$, displayed in Figure 2. However, if we consider the symmetrical graph $\bar{G}$, BFS$_{\bar{G}}(c, a)$ is smaller than BFS$_G(c, b)$ (see Figure 3). This code is the smallest one for both $G$ and $\bar{G}$, i.e., $cc(G) = 00 01 03 32$.

## 2.2 Properties of Canonical Codes

The next two theorems are used to efficiently enumerate canonical codes with CP in Section 4.

**Theorem 5.** Given a code $p_1d_1 \ p_2d_2 \ \dots \ p_n d_n$, we have:

**Property a:** $\forall i \in [1, n-1], p_i \le p_{i+1}$;

**Property b:** $\forall i \in [1, n-1], (p_i = p_{i+1}) \Rightarrow (d_i < d_{i+1})$.

**Proof.** Property a is a consequence of the fact that (i) $q$ is a FIFO queue, (ii) vertices are added in $q$ by order of increasing number, and (iii) when a vertex $v_i$ is removed from $q$ we add to code all couples $p_i d_i$ such that $p_i = \text{num}[v_i]$ and the vertex at direction $d_i$ from $v_i$ is
not yet numbered. Hence, all couples added to code during one iteration of the loop lines 6-12 have the same value for $p_i$. Property b is a straightforward consequence of the fact that the loop lines 8-12 considers edges by order of increasing directions. ◀

**▶ Theorem 6.** Let $c = p_1d_1 p_2d_2 \dots p_nd_n$ be a canonical code. For each $i \in [1, n-1]$, the prefix $c_i = p_1d_1 p_2d_2 \dots pimd_i$ of c is also a canonical code.

**Proof.** Let $G$ and $G_i$ be the graphs associated with $c$ and $c_i$, where each vertex is numbered with $num[v_i]$. As $c_i$ is a prefix of $c$, $G_i$ is the subgraph of $G$ induced by the vertices numbered from 0 to $i$. Let us assume that $c_i$ is not canonical, i.e., there exists another code $c'_i = p'_1d'_1 p'_2d'_2 \dots p'_id'_i$ that is smaller than $c_i$ while representing the same subgraph $G_i$. Let $(v_0, v_1)$ be the edge of $G_i$ such that $c'_i = \text{BFS}_{G_i}(v_0, v_1)$ (resp. $c'_i = \text{BFS}_{\overline{G_i}}(v_0, v_1)$ if $c'_i$ is computed in the symmetrical graph). Let $c' = \text{BFS}_G(v_0, v_1)$ (resp. $c' = \text{BFS}_{\overline{G}}(v_0, v_1)$) be the code obtained by performing a search of $G$ (resp. $\overline{G}$) that starts from the same initial edge as the one used to obtain $c'_i$. Let $j$ be the smallest vertex number for which $c'$ and $c'_i$ have different values, i.e., $\forall k < j$, vertex number $k$ has the same edges in both $G$ and $G_i$, whereas vertex number $j$ has more edges in $G$ than in $G_i$. More precisely, let $S$ (resp. $S_i$) be the set of edges outgoing from $j$ in $G$ (resp. $G_i$). We have $S_i \subset S$ as $G_i$ is a subgraph of $G$. We have to distinguish two different cases.

**Case 1:** the largest edge direction in $S_i$ is equal to the largest edge direction in $S$. In this case, $c'$ is smaller than $c'_i$ because the sequence of couples associated with the successors of $j$ in $c'$ is necessarily smaller than the sequence of couples associated with the successors of $j$ in $c'_i$ (e.g., if edge directions in $S_i$ are $\{1, 5\}$ whereas edge directions in $S$ are $\{1, 2, 5\}$, then $j1 j2 j5 < j1 j5$).

**Case 2:** the largest edge direction in $S_i$ is smaller than the largest edge direction in $S$. In this case, either $c'_i$ is a prefix of $c'$ (if $j$ is the greatest vertex with an outgoing edge), or $c'$ is smaller than $c'_i$ (e.g., if edge directions in $S_i$ are $\{1, 2\}$ whereas edge directions in $S$ are $\{1, 2, 5\}$, $j1 j2 j5 < j1 j2 k$ because $k > j$, as stated in Property a).

In both cases, this implies that $c' < c$ which contradicts the fact that $c$ is canonical. ◀

**Computational complexity.** The canonical code of a graph with $n$ vertices is computed in $O(n^2)$: BFS is run twice for each edge and the number of edges is in $O(n)$ as a vertex cannot have more than six neighbours; the time complexity of BFS is $O(n)$ as the loop lines 6-12 is iterated $n$ times and the loop lines 8-12 is iterated at most six times.

We may speed-up the computation by avoiding some calls to BFS. To this aim, we display in Figure 4 the 13 different possible patterns for the neighbourhood of a vertex $v$ (up to rotations). For each pattern $P_i$ such that $v$ has $k_i$ successors, every code computed by BFS with $v_0 = v$ starts with a prefix $0d_1 0d_2 \dots 0d_{k_i}$, and we give below the smallest possible prefix composed of $2k_i$ integers, denoted $\pi_i$, the set $S$ of edges $(v, v_k)$ of $G$ such that $\text{BFS}_G(v, v_k)$ starts with $\pi_i$, and the set $\overline{S}$ of edges $(v, v_k)$ of $G$ such that $\text{BFS}_{\overline{G}}(v, v_k)$ starts with $\pi_i$ (these edges are highlighted in Figure 4).

<table><tr><td>π<sub>6</sub> = 00 01 02 03 04 05</td><td>S = {(v, v<sub>i</sub>) : i ∈ [1, 6]}</td><td>S&#773; = {(v, v<sub>i</sub>) : i ∈ [1, 6]}</td></tr><tr><td>π<sub>5</sub> = 00 01 02 03 04</td><td>S = {(v, v<sub>1</sub>)}</td><td>S&#773; = {(v, v<sub>5</sub>)}</td></tr><tr><td>π<sub>4</sub> = 00 01 02 03</td><td>S = {(v, v<sub>1</sub>)}</td><td>S&#773; = {(v, v<sub>4</sub>)}</td></tr><tr><td>π<sub>3+1</sub> = 00 01 02 04</td><td>S = {(v, v<sub>1</sub>)}</td><td>S&#773; = {(v, v<sub>3</sub>)}</td></tr><tr><td>π<sub>2+2</sub> = 00 01 03 04</td><td>S = {(v, v<sub>1</sub>), (v, v<sub>3</sub>)}</td><td>S&#773; = {(v, v<sub>2</sub>), (v, v<sub>4</sub>)}</td></tr></table>
<figure><img src="image_282.png" /><figcaption>Figure 4 Different neighbourhood patterns for a vertex v. We highlight in blue (resp. green and red) every edge (v, v_k) such that BFS_G(v, v_k) (resp. BFS_Ḡ(v, v_k) and both BFS_G(v, v_k) and BFS_Ḡ(v, v_k)) starts with the smallest possible prefix.</figcaption></figure>

<table><tr><td>π<sub>3</sub> = 00 01 02</td><td>S = {(v, v<sub>1</sub>)}</td><td>S&#773; = {(v, v<sub>3</sub>)}</td></tr><tr><td>π<sub>2+1a</sub> = 00 01 03</td><td>S = {(v, v<sub>1</sub>)}</td><td>S&#773; = ∅</td></tr><tr><td>π<sub>2+1b</sub> = 00 01 03</td><td>S = ∅</td><td>S&#773; = {(v, v<sub>2</sub>)}</td></tr><tr><td>π<sub>1+1+1</sub> = 00 03 05</td><td>S = {(v, v<sub>1</sub>), (v, v<sub>2</sub>), (v, v<sub>3</sub>)}</td><td>S&#773; = {(v, v<sub>1</sub>), (v, v<sub>2</sub>), (v, v<sub>3</sub>)}</td></tr><tr><td>π<sub>2</sub> = 00 01</td><td>S = {(v, v<sub>1</sub>)}</td><td>S&#773; = {(v, v<sub>2</sub>)}</td></tr><tr><td>π<sub>1+1a</sub> = 00 02</td><td>S = {(v, v<sub>1</sub>)}</td><td>S&#773; = {(v, v<sub>2</sub>)}</td></tr><tr><td>π<sub>1+1b</sub> = 00 03</td><td>S = {(v, v<sub>1</sub>), (v, v<sub>2</sub>)}</td><td>S&#773; = {(v, v<sub>1</sub>), (v, v<sub>2</sub>)}</td></tr><tr><td>π<sub>1</sub> = 00</td><td>S = {(v, v<sub>1</sub>)}</td><td>S&#773; = {(v, v<sub>1</sub>)}</td></tr></table>

**Example 7.** Let us consider the graph G of Figure 1. The pattern of vertex a in G is P<sub>2</sub> (with a = v, c = v<sub>1</sub>, and b = v<sub>2</sub>). Therefore, the smallest possible code computed with v<sub>0</sub> = a starts with the prefix π<sub>2</sub> = 00 01, and the two codes that start with this prefix are BFS<sub>G</sub>(a, c) and BFS<sub>Ḡ</sub>(a, b). Also, the pattern of vertex c in G is P<sub>2+1b</sub> (with c = v, a = v<sub>2</sub>, b = v<sub>1</sub>, and d = v<sub>3</sub>). Therefore, the smallest possible code computed with v<sub>0</sub> = c starts with the prefix π<sub>2+1a</sub> = 00 01 03, and the only code that starts with this prefix is BFS<sub>Ḡ</sub>(c, a).

To further reduce the number of searches, we define a total order among patterns.

**Definition 8 (Order between patterns).** Let P<sub>i</sub> and P<sub>j</sub> be two patterns (as listed in Figure 4), and let π<sub>i</sub>1 and π<sub>j</sub>1 be the sequences obtained by adding “1” at the end of π<sub>i</sub> and π<sub>j</sub>, respectively. We define P<sub>i</sub> < P<sub>j</sub> (resp. P<sub>i</sub> = P<sub>j</sub>) if π<sub>i</sub>1 is lexicographically smaller than (resp. equal to) π<sub>j</sub>1. We obtain the following order:

$$P_6 < P_5 < P_4 < P_{3+1} < P_3 < P_{2+2} < P_{2+1a} = P_{2+1b} < P_2 < P_{1+1+1} < P_{1+1a} < P_{1+1b} < P_1$$

Given two vertices v<sub>i</sub> and v<sub>j</sub>, if the pattern of v<sub>i</sub> is smaller than the pattern of v<sub>j</sub>, then we know that every code computed from v<sub>i</sub> is smaller than every code computed from v<sub>j</sub>. Hence, to compute the canonical code of a hexagon graph G, we first search for the smallest pattern P<sub>k</sub> in G. Then, for each vertex v<sub>i</sub> such that the pattern of v<sub>i</sub> in G is equal to P<sub>k</sub>, we compute all codes from the starting edges highlighted in Figure 4. Finally, we return the smallest computed code.

**Example 9.** Let us consider the graph G of Figure 1. The pattern of vertex a (resp. b, c, d, and e) is P<sub>2</sub> (resp. P<sub>2</sub>, P<sub>2+1b</sub>, P<sub>1+1a</sub>, and P<sub>1</sub>). The smallest pattern is P<sub>2+1b</sub>. Therefore, the canonical code is obtained by running BFS<sub>G</sub>(c, a), and it is useless to run BFS<sub>G</sub> or BFS<sub>Ḡ</sub> from any other edge.
Note that this does not change the time complexity as in the worst case we may have $O(n)$ vertices with a smallest pattern. For example, when we have an horizontal chain of vertices, all vertices have pattern $P_{1+1b}$, except the two endpoints which have pattern $P_1$. In this case, we must run BFS$_G$ and BFS$_{\bar{G}}$ for every edge outgoing from vertices with pattern $P_{1+1b}$ (i.e., all vertices but the two endpoints).

## 2.3 Comparison with other canonical codes

In the general case, building the canonical code of a graph is a problem which is not known to be in $\mathcal{P}$ nor to be $\mathcal{NP}$-complete (it is isomorphic-complete). There exist rather efficient algorithms such as Nauty [9], but these algorithms have an exponential time complexity in the worst case. Canonical codes are widely used in graph mining tools such as gSpan [12] or Gaston [10]. When graphs are embedded in a 2D space, canonical codes may be computed in polynomial time [8]. These canonical codes may be simplified when considering grid graphs such that all faces are squares [6]. Canonical codes defined in [12, 10, 6] share some similarities with our canonical codes as they are computed by performing graph traversals and assigning numbers to vertices according to the order of discovery. However, in [12, 10, 6], traversals are Depth First Searches (DFSs). Performing BFSs instead of DFSs allows us to avoid some calls to BFS (as explained in Section 2.2). It also allows us to exploit properties a and b to define codes by means of constraints (see Section 4).

In [1], canonical codes are introduced for benzenoid structures. These codes only describe boundary cycles and assume that there is no hole within these cycles, thus preventing the representation of coronoids. Our codes fully describe hexagon graphs, including vertices inside the boundary cycle and, therefore, we can represent coronoids.

# 3 Canonical Codes in Case of Required Patterns

Many BGPs involve enumerating hexagon graphs that contain some given patterns, where patterns both specify mandatory and forbidden vertices [3]. More precisely, a pattern is defined by a couple $(P, F)$ such that $P$ is a hexagon graph that specifies mandatory vertices, and $F$ is a set of couples that specify forbidden vertex positions, i.e., for each couple $(p_i, d_i) \in F$, vertex $p_i$ must not have an outgoing edge in direction $d_i$.

**Example 10.** The pattern $(P, F)$ displayed in Figure 5 has four mandatory vertices in grey and one forbidden vertex in red. $(P, F)$ occurs twice in $G$ as $P$ is isomorphic to $G_{\downarrow\{d,a,b,c\}}$ (resp. $G_{\downarrow\{g,d,a,f\}}$) and $a$ has no outgoing edge in direction 1 (resp. $a$ has no outgoing edge in direction 1 when rotating $G$ so that $G_{\downarrow\{g,d,a,f\}}$ has the same orientation as $P$).

To efficiently generate graphs that contain a pattern, the idea is to constrain canonical codes to start with this pattern so that all generated graphs contain it by construction. We introduce this new canonical code in Section 3.1, and we study its properties in Section 3.2.

### 3.1 Definition of P-canonical codes

To ensure that codes start with the canonical code of a given pattern $(P, F)$, we perform BFSs of $G$ and $\bar{G}$ from a given P-sequence that corresponds to an occurrence of $(P, F)$ in $G$ or $\bar{G}$. More precisely, this P-sequence is defined as follows.

**Definition 11 (P-sequence of $(G, P, F)$).** Let $G$ be a hexagon graph with $n$ vertices, and $(P, F)$ be a pattern with $1 < k < n$ mandatory vertices. A P-sequence of $(G, P, F)$ is a sequence $s = \langle v_0, v_1, \dots, v_{k-1} \rangle$ of $k$ vertices of $G$ such that
<figure><img src="image_283.png" /><figcaption>■ **Figure 5** Left: Example of hexagon graph $G$ and pattern $(P, F)$. Right: To compute $cc_{P,F}(G)$, we run P-BFS$_G$ (resp. P-BFS$_{\bar{G}}$) with the P-sequences $\langle a, b, d, c \rangle$ and $\langle a, d, f, g \rangle$ (resp. $\langle b, a, c, d \rangle$ and $\langle d, a, g, f \rangle$). The P-canonical code of $G$ is $cc_{P,F}(G) = 00\ 02\ 11\ 03\ 04\ 56$.</figcaption></figure>

(i) the subgraph of $G$ induced by $\{v_0, \dots, v_{k-1}\}$, i.e., $G_{\downarrow s}$, is isomorphic to $P$;

(ii) $\forall (p_i, d_i) \in F$, $G_{\downarrow s}$ has no outgoing edge from the vertex corresponding to $p_i$ in direction $d_i$ when rotating $G$ so that $G_{\downarrow s}$ and $P$ have the same orientation;

(iii) the code returned by $BFS_{G_{\downarrow s}}(v_0, v_1)$ is canonical;

(iv) for each $i \in [0, k-1]$ the number assigned by $BFS_{G_{\downarrow s}}(v_0, v_1)$ to $v_i$ is equal to $i$, i.e., $\text{num}[v_i] = i$ at line 13 of Algorithm 1.

Each P-sequence corresponds to an occurrence of $(P, F)$ in $G$, and $(P, F)$ may occur more than once in $G$. To compute all P-sequences of $(G, P, F)$, we iterate on each edge of $G$ and try to build a P-sequence that starts with this edge using Algorithm 2: given the canonical code $c$ of $P$, an edge $(v_0, v_1)$ of $G$, and a set $F$ of forbidden positions, $\text{seq}(G, v_0, v_1, c, F)$ returns the P-sequence of $(G, P, F)$ that starts with $\langle v_0, v_1 \rangle$ if it exists, and it returns null otherwise. To build the P-sequence, it uses the canonical code $c$ to associate vertices of $G$ to vertex numbers used in the canonical code: for each vertex number $i \in [0, k-1]$, $\text{vertex}[i]$ is the vertex of $G$ that corresponds to $i$, and for each couple $p_i d_i$ in the canonical code, we ensure that vertex $i$ in $P$ corresponds to the vertex of $G$ that is at direction $d_i$ from $\text{vertex}[p_i]$, if it exists (line 5); if it does not exist, then there is no P-sequence that starts from $(v_0, v_1)$ and we return null (line 4). When a P-sequence is completed, we check that there are no vertices at forbidden positions (line 6).

To be invariant to symmetries, we must also compute all P-sequences of $(\bar{G}, P, F)$ by running $\text{seq}(\bar{G}, v_0, v_1, c, F)$ for each edge $(v_0, v_1)$ of $\bar{G}$.

► **Example 12.** Let us consider the graph $G$ and the pattern $(P, F)$ of Figure 5. The canonical code of $P$ is 00 02 11, and it may be computed either by $BFS_P(0, 1)$ or by $BFS_{\bar{P}}(1, 0)$ because $P$ is symmetrical. Also, this pattern is isomorphic to two subgraphs of $G$. Hence, the four runs of seq that return a P-sequence corresponding to an occurrence of $(P, F)$ are:
■ **Algorithm 2** seq(G, v₀, v₁, c, F).

**Input:** A hexagon graph G, an edge (v₀, v₁) of G, the canonical code
c = p₁d₁ ... pₖ₋₁dₖ₋₁ of a pattern P, and the set
F = {(p'₁, d''₁), ..., (p'ₕ, d''ₕ)} of forbidden vertices

**Output:** A P-sequence of (G, P, F) starting with ⟨v₀, v₁⟩, or null if such a P-sequence does not exist

1 rotate G so that the direction of edge (v₀, v₁) is equal to 0

2 vertex[0] ← v₀; vertex[1] ← v₁

3 **for** i ∈ [2, k − 1] **do**

4   **if** vertex[pᵢ] does not have an outgoing edge in direction dᵢ **then** return null;

5   vertex[i] ← vertex at direction dᵢ from vertex[pᵢ]

6 **if** ∃(p'_i, d''_i) ∈ F s.t. vertex[p'_i] has an outgoing edge in direction d''_i **then** return null;

7 **return** ⟨vertex[0], vertex[1], ..., vertex[k − 1]⟩

■ seq(G, a, b, 00 02 11, {(0, 1)}) = ⟨a, b, d, c⟩ corresponding to 0 = a, 1 = b, 2 = d, 3 = c;

■ seq(G, b, a, 00 02 11, {(0, 1)}) = ⟨b, a, c, d⟩ corresponding to 0 = b, 1 = a, 2 = c, 3 = c;

■ seq(G, a, d, 00 02 11, {(0, 1)}) = ⟨a, d, f, g⟩ corresponding to 0 = a, 1 = d, 2 = f, 3 = g;

■ seq(G, d, a, 00 02 11, {(0, 1)}) = ⟨d, a, g, f⟩ corresponding to 0 = d, 1 = a, 2 = g, 3 = f.

All other runs of seq return null.

In Figure 6, the graph G has only one P-sequence for pattern P (i.e., seq(G, 0, 1, cc(P), F) = ⟨0, 1, 2, 3⟩). seq(G, 2, 4, cc(P), F) = null because G already has a vertex at direction 4 from 2 when rotating G so that the direction of edge (2, 4) is equal to 0.

Given a hexagon graph G and a P-sequence s = ⟨v₀, v₁, ..., vₖ₋₁⟩, we define an algorithm, called P-BFS<sub>G</sub>(s), which performs a BFS of G that starts from s, and which is obtained from Algorithm 1 by replacing lines 3 and 4 by the following lines:

display the canonical code of G<sub>↓s</sub>

let q be an empty FIFO queue

**for** i ∈ [0, k − 1] **do** set num[vᵢ] to i, and add vᵢ in q

initialise the counter c to k

These lines ensure that the code returned by P-BFS starts with the canonical code
of the pattern associated with s so that this pattern is included in the returned code by
construction. The numbers assigned to the vertices of s correspond to those assigned by
BFS when computing the canonical code of the pattern and we add the vertices of s in q to
ensure that these vertices are treated first.

The P-canonical code of (G, P, F) is the smallest code obtained from any P-sequence.

► **Definition 13 (P-canonical code (cc<sub>P,F</sub>(G)).)** The P-canonical code of G = (V, E) with respect to a pattern (P, F) is cc<sub>P,F</sub>(G) = min {P-BFS<sub>G</sub>(s) : s ∈ S} ∪ {P-BFS<sub>Ḡ</sub>(s) : s ∈ S̄} where

■ S = {seq(G, v₀, v₁, cc(P), F) : (v₀, v₁) ∈ E ∧ seq(G, v₀, v₁, c, F) ≠ null},

■ $\bar{S} = \{ \overline{\text{seq}}(\bar{G}, v_0, v_1, cc(P), F) : (v_0, v_1) \in E \land \text{seq}(\bar{G}, v_0, v_1, c, F) \neq \text{null} \}$.

► **Example 14.** Let us consider the graph *G* and the pattern (*P*, *F*) of Figure 5. As seen in Example 12, there are four different candidate P-sequences *s* from which we may compute a code. We display in Figure 5 the runs of P-BFS from these four P-sequences.
<figure><img src="image_284.png" /><figcaption>■ **Figure 6** Example of P-canonical code with a non P-canonical prefix. Left: Pattern $(P, F)$. Middle: Graph $G$. Right: $G_{\downarrow\{0,1,2,3,4\}}$. The prefix `00 01 03 21` of $cc_{P,F}(G)$ corresponds to $G_{\downarrow\{0,1,2,3,4\}}$ and it is not P-canonical as $cc_{P,F}(G_{\downarrow\{0,1,2,3,4\}})$ is smaller.</figcaption></figure>

## 3.2 Properties of P-canonical codes

The time complexity for computing $cc_{P,F}(G)$ when $G$ has $n$ vertices, $P$ has $k$ vertices (with $k < n$), and $F$ has $f$ couples is $O(n(n+f))$:

- the canonical code of $P$ is computed in $O(k^2)$;

- Algorithm 2 is in $O(k+f)$, and it is called $O(n)$ times to build all P-sequences;

- P-BFS is in $O(n)$ and it is called once per P-sequence, i.e., $O(n)$ times.

Properties of Theorem 5 are no longer valid as illustrated in Figure 5: the code returned by $P-BFS_G(\langle a, d, f, g \rangle)$ does not satisfy Property a as $p_3 = 1$ and $p_4 = 0$; the code returned by $P-BFS_G(\langle d, a, g, f \rangle)$ does not satisfy Property b as $p_3 = p_4 = 1$ whereas $d_3 = 1$ and $d_4 = 0$. This comes from the fact that the canonical code of $P$ is inserted at the beginning of the sequence, before starting the search. If a vertex of $G$ corresponding to a pattern vertex has some extra edges outgoing from it in $G$ but not in $P$, then the integer couple associated with this edge is added after $cc(P)$ and it may violate Properties a or b. However, we can easily show that these properties are satisfied in the second part of the code, corresponding to the vertices that do not correspond to pattern vertices, i.e., given a code $p_1d_1 \dots p_nd_n$ returned by $P-BFS_G(\langle v_0, \dots, v_{k-1} \rangle)$, we have:

**Property a':** $\forall i \in [k, n-1], p_i \le p_{i+1}$;

**Property b':** $\forall i \in [k, n-1], (p_i = p_{i+1}) \Rightarrow (d_i < d_{i+1})$.

When there are no forbidden vertices (i.e., $F = \emptyset$), Theorem 6 is still valid, i.e., every prefix of a canonical code is also canonical (the proof is similar to the proof of Theorem 6 and it is omitted due to space limits).

However, this theorem is no longer valid whenever $F \neq \emptyset$, i.e., the prefix of a P-canonical code may not be P-canonical, as shown in the following example.

► **Example 15.** In Figure 6, the occurrence of $(P, F)$ in $G$ corresponds to $G_{\downarrow\{0,1,2,3\}}$. The prefix `00 01 03 21` of $cc_{P,F}(G)$ is not P-canonical: this prefix corresponds to $G_{\downarrow\{0,1,2,3,4\}}$ and $(P, F)$ occurs twice in it (once in $G_{\downarrow\{0,1,2,3\}}$ with the forbidden vertex displayed in brown, and once in $G_{\downarrow\{0,1,2,4\}}$ with the forbidden vertex displayed in red). If we extend $G_{\downarrow\{0,1,2,3,4\}}$ by adding a vertex on the brown position, we obtain the graph $G$ that contains only one occurrence of $(P, F)$ and whose P-canonical code starts with a prefix greater than `00 01 03 15`. If we extend $G_{\downarrow\{0,1,2,3,4\}}$ by adding a vertex on the red position, we obtain another graph (not isomorphic to $G$) that also contains only one occurrence of $(P, F)$.
This example shows us that a non P-canonical code may be extended to a P-canonical code. This comes from the fact that a pattern may have more occurrences in the graph associated with a prefix of a code *c* than in the graph associated with *c*, due to forbidden vertices. However, when the pattern is symmetrical (including its forbidden vertices), we may discard some non P-canonical codes. More precisely, we define below the dominated codes that may be discarded because they cannot be extended to P-canonical codes.

**► Definition 16.** (Dominated code) Let (P, F) be a symmetrical pattern, G = (V, E) be a hexagon graph that contains at least one occurrence of (P, F), and W ⊂ V be a subset of vertices such that G↓W is isomorphic to P. Let C be the set of codes computed from a P-sequence s = (v₀, v₁, ..., vₖ) such that (v₀, v₁) is an edge of G↓W. A code c ∈ C is dominated if there exists another code c' ∈ C such that c' < c.

**► Theorem 17.** For each P-canonical code *c* and each prefix *c'* of *c*, *c'* is not dominated.

*Proof.* This is a straightforward consequence of the fact that all P-sequences that start from an edge (*v*₀, *v*₁) of *G*↓*W* correspond to different automorphisms of (P, F) and are interchangeable.

**► Example 18.** In Figure 5, (P, F) is symmetrical. The dominated codes are P-BFS<sub>G</sub>(⟨b, a, c, d⟩) (when W = {a, b, c, d}) and P-BFS<sub>G</sub>(⟨d, a, g, f⟩) (when W = {a, d, f, g}. However, we cannot discard P-BFS<sub>G</sub>(⟨a, b, d, c⟩) (though it is larger than P-BFS<sub>G</sub>(⟨a, d, f, g⟩) as it may lead to a P-canonical code if a vertex is added on the forbidden position for W = {a, d, f, g}.

# 4 CP Models for Enumerating Hexagon Graphs

As hexagon graphs may be uniquely represented with canonical codes, we solve BGPs by enumerating canonical codes, and we propose to use CP to achieve this enumeration as this allows us to easily add constraints on the structures. In Section 4.1, we introduce a first CP model for enumerating consistent codes. In Section 4.2, we introduce a global constraint for ensuring that codes are canonical. In Section 4.3, we introduce a CP model for enumerating all structures that contain some given patterns.

## 4.1 CP Model for Enumerating Consistent Codes

We say that a code is consistent if there exist a hexagon graph *G* and an edge (*v*₀, *v*₁) of *G* such that BFS<sub>*G*</sub>(*v*₀, *v*₁) or BFS<sub>*G*</sub>(*v*₀, *v*₁) return this code. Obviously, the code must satisfy Properties a and b listed in Theorem 5. However, this is not enough because we must also ensure that all vertices have different positions in the plane. To this aim, we associate 2D coordinates with vertices. We assume that all edges have a length equal to two and that vertex 0 has coordinates (0, 0). For every other vertex *i*, its coordinates depend on the coordinates of its predecessor *p*ᵢ and the direction *d*ᵢ of edge (*p*ᵢ, *i*): if the coordinates of *p*ᵢ are (*x*, *y*) and *d*ᵢ = 0 (resp. 1, 2, 3, 4, and 5), then vertex *i* is at coordinates (*x* + 2, *y*) (resp. (*x* + 1, *y* + √3), (*x* − 1, *y* + √3), (*x* − 2, *y*), (*x* − 1, *y* − √3), and (*x* + 1, *y* − √3)).

**► Example 19.** Let us consider the code 00 02 24. Vertex 1 is at coordinates (2, 0) (because the edge that reaches 1 from 0 has direction 0), vertex 2 is at coordinates (1, √3) (because the edge that reaches 2 from 0 has direction 1), and vertex 3 is at coordinates (0, 0) (because the edge that reaches 3 from 2 has direction 4). As vertices 0 and 3 have the same coordinates, this code is not consistent.
28:12 Using Canonical Codes to Efficiently Solve the BGP with CP

For each vertex $i \in [1, n]$, our CP model uses four integer variables $x_i, y_i, p_i, d_i$:

$x_i$ is the abscissa of $i$, and its domain is $D(x_i) = [-2n, 2n]$;

$y_i$ is the ordinate of $i$ divided by $\sqrt{3}$, and its domain is $D(y_i) = [-n, n]$;

$p_i$ is the vertex which has discovered $i$, and its domain is $D(p_i) = [0, i - 1]$ ($p_i < i$ because $i$ cannot be discovered by a vertex $j$ which has not yet been discovered);

$d_i$ is the direction of edge $(p_i, i)$, and its domain is $D(d_i) = [0, 5]$.

For vertex 0, we define $x_0 = 0$, $y_0 = 0$, $p_0 = -1$, and $d_0 = -1$. The code associated with these variables is $p_1d_1 \dots p_n d_n$

To ensure the consistency of the generated codes, we post the following constraints:

$C_1$: all vertices have different coordinates, i.e., $\forall i, j \in [0, n], i \neq j \Rightarrow x_i \neq x_j \lor y_i \neq y_j$;

$C_2$: Properties a and b (defined in Theorem 5) are satisfied, i.e.,

$$ \forall i \in [1, n-1], p_i \le p_{i+1} \land p_i = p_{i+1} \Rightarrow d_i < d_{i+1}; $$

$C_3$: For each edge $(p_i, i)$, coordinates of $i$ and $p_i$ are consistent with direction $d_i$, i.e.,

$$ \forall i \in [1, n], (d_i, x_i, x_{p_i}, y_i, y_{p_i}) \in T \text{ where } T \text{ is the table defined as follows:} $$

$$ \begin{aligned} T = & \{(0, x + 2, x, y, y) : x \in [-2n, 2n - 2], y \in [-n, n]\} \\ & \cup \{(1, x + 1, x, y + 1, y) : x \in [-2n, 2n - 1], y \in [-n, n - 1]\} \\ & \cup \{(2, x - 1, x, y + 1, y) : x \in [-2n + 1, 2n], y \in [-n, n - 1]\} \\ & \cup \{(3, x - 2, x, y, y) : x \in [-2n + 1, 2n], y \in [-n, n]\} \\ & \cup \{(4, x - 1, x, y - 1, y) : x \in [-2n + 1, 2n], y \in [-n + 1, n]\} \\ & \cup \{(5, x + 1, x, y - 1, y) : x \in [-2n, 2n - 1], y \in [-n + 1, n]\}. \end{aligned} $$

## 4.2 Global Constraint for Ensuring Canonicity

We must ensure that codes are canonical to become invariant to rotations and symmetries. To this aim, we introduce a new global constraint defined below.

**► Definition 20.** Given an integer value $n \ge 1$ and a tuple of four integer variables $(p_i, d_i, x_i, y_i)$ for each $i \in [1, n]$, the constraint canonical($\{(p_i, d_i, x_i, y_i) : i \in [1, n]\}$) is satisfied if the code $p_1d_1 \dots p_n d_n$ is canonical.

To propagate this constraint, we maintain a vector *pat* such that, for each vertex *i*, *pat[i]* is equal to the current pattern of *i* (with respect to the edges that have been added so far). We also maintain a matrix *M* such that for each vertex *i* and each direction *d* ∈ [0, 5], *M[i][d]* is either equal to the vertex at direction *d* from *i*, if such a vertex exists, or to -1 otherwise. *pat* and *M* are updated each time a new couple (*p*ᵢ, *d*ᵢ) is instantiated. This may be done in constant time by exploiting vertex coordinates.

We trigger a failure whenever there is a vertex $i > 0$ such that *pat[i] < pat[0]*, as we have seen in Section 2.2 that the pattern of *i* cannot be smaller than the pattern of 0 (according to the order of Definition 8).

Also, we trigger a failure whenever assigned variables define a code prefix which is not canonical, because Theorem 6 tells us that a non-canonical prefix cannot be extended to a canonical code. More precisely, when all variables that occur in a code prefix $p_1d_1 \dots p_k d_k$ with $k \le n$ are assigned, we check that this prefix is canonical and, if it is not the case, we trigger a failure. To check whether a prefix is canonical or not, we try to build smaller codes by running BFS from other edges than (0, 1). We limit the number of BFSs to be performed by exploiting patterns, as explained in Section 2.2.
Finally, we ensure that the degree of every vertex is upper bounded by a value *g* which depends on *pat[0]*: *g* = 5 when *pat[0]* = *P*<sub>5</sub>; *g* = 4 when *pat[0]* ∈ {*P*<sub>4</sub>, *P*<sub>3+1</sub>, *P*<sub>3</sub>, *P*<sub>2+2</sub>}; *g* = 3 when *pat[0]* ∈ {*P*<sub>2+1a</sub>, *P*<sub>2+1b</sub>, *P*<sub>2</sub>, *P*<sub>1+1+1</sub>}; *g* = 2 when *pat[0]* = *P*<sub>1+1</sub>; and *g* = 1 when *pat[0]* = *P*<sub>1</sub>. When the degree of a vertex *i* reaches *g*, we remove *i* from the domain of every non-assigned variable *p*<sub>*k*</sub> in order to prevent *i* from having more outgoing edges than *g*.

4.3 CP Model for Enumerating Graphs with a Given Pattern

Given a pattern (*P*, *F*) with *m* mandatory vertices, the CP model for enumerating all graphs with *n* + 1 > *m* vertices that contain (*P*, *F*) uses four integer variables *x*<sub>*i*</sub>, *y*<sub>*i*</sub>, *p*<sub>*i*</sub>, and *d*<sub>*i*</sub> for each vertex *i* ∈ [0, *n*], like in Section 4.1. However, for each pattern vertex *i* ∈ [0, *m* − 1], we assign the variables *x*<sub>*i*</sub>, *y*<sub>*i*</sub>, *p*<sub>*i*</sub>, and *d*<sub>*i*</sub> according to the canonical code *cc*(*P*).

**Example 21.** For the pattern *P* displayed in Figure 6, we set $x_0 = y_0 = 0$, $p_1 = p_2 = p_3 = 0$, $d_1 = 0$, $d_2 = 1$, $d_3 = 3$, $x_1 = 2$, $y_1 = 0$, $x_2 = 1$, $y_2 = 1$, $x_3 = -2$, and $y_3 = 0$.

Like in Section 4.1, we post constraint $C_1$ to ensure that all vertices have different coordinates. We modify constraint $C_2$ as Properties of Theorem 5 are satisfied only for vertices that are not in the pattern, i.e., $\forall i \in [m, n-1]$, $p_i \le p_{i+1} \land p_i = p_{i+1} \Rightarrow d_i < d_{i+1}$. We post constraint $C_3$ but with a modified table $T'$ to ensure that no vertex is at a forbidden position. More precisely, we compute the set $X_F$ of coordinates of forbidden positions in $F$ and define $T'$ from table $T$ as follows: $T' = \{(d, x, x', y, y') \in T : (x, y) \notin X_F \land (x', y') \notin X_F\}$.

Finally, we post a global constraint to ensure that the generated codes are P-canonical.
This constraint is defined as follows.

**Definition 22.** Given the canonical code *cc*(*P*) of a pattern *P* with *m* vertices, an integer value *n* > *m* and a tuple of four integer variables (*p*<sub>*i*</sub>, *d*<sub>*i*</sub>, *x*<sub>*i*</sub>, *y*<sub>*i*</sub>) for each *i* ∈ [1, *n*], the constraint *P-canonical*<sub>*cc*(*P*)</sub>({(*p*<sub>*i*</sub>, *d*<sub>*i*</sub>, *x*<sub>*i*</sub>, *y*<sub>*i*</sub>) : *i* ∈ [1, *n*]}) is satisfied if the code $p_1 d_1 \dots p_n d_n$ is *P-canonical*.

To propagate this constraint, we cannot trigger a failure whenever assigned variables
define a code prefix which is not P-canonical as we have seen in Section 3.2 that Theorem 6
is no longer valid. However, we exploit Theorem 17 to trigger a failure whenever the current
prefix is dominated. Also when all $p_i$ and $d_i$ variables are assigned, we check that the code is
P-canonical and trigger a failure whenever this is not the case.

Some BGPs involve enumerating graphs that contain two given patterns (*P*, *F*) and (*P'*, *F'*). Without loss of generality, we assume that #P ≥ #P' where #G denotes the number of vertices of a graph *G*. In this case, we post the constraint *P-canonical*<sub>*cc*(*P*)</sub> to ensure that all generated codes contain pattern (*P*, *F*) by construction. To ensure that the graph *G* associated with the generated code also contains pattern (*P'*, *F'*), we post a global constraint that ensures that $G_{\downarrow[k,\#G]}$ contains an occurrence of (*P'*, *F'*), where *k* = #P whenever the two patterns must be disjoint, whereas *k* = 0 otherwise. This constraint is defined below.

**Definition 23.** Given a pattern (*P*, *F*), two integer values *k* and *n* such that *k*+#*P* ≤ *n*, and a tuple of two integer variables (*p*<sub>*i*</sub>, *d*<sub>*i*</sub>) for each *i* ∈ [0, *n*], the constraint *subgraph*<sub>*P*,*k*</sub>({(*p*<sub>*i*</sub>, *d*<sub>*i*</sub>) : *i* ∈ [0, *n*]}) is satisfied if there exists *S* ⊆ [*k*, *n*] such that $G_{\downarrow S}$ is isomorphic to *P* and there is no vertex at a forbidden position.

We have implemented a very basic propagator for this constraint: it is not propagated dur-
ing the search, and we only check that it is entailed when all variables have been assigned. This
is done by searching for an edge (*v*₀, *v*₁) of $G_{\downarrow[k,n]}$ such that $\text{seq}(G_{\downarrow[k,n]}, v_0, v_1, \text{cc}(P'), F') \neq \text{null}$. If we do not find such an edge, then the constraint is not satisfied; otherwise we have
found an occurrence of $(P', F')$ in $G_{\downarrow[k,n]}$.
28:14 Using Canonical Codes to Efficiently Solve the BGP with CP

<table><caption>Table 1 Results of BENZAI and CCODE for BGP1 and BGP2. For each $n \in [2, 10]$ and each problem, we report the number of solutions (#sol) and the time in seconds of BENZAI and CCODE.</caption><thead><tr><th rowspan="2"></th><th colspan="3">BGP1</th><th colspan="3">BGP2</th></tr><tr><th>#sol</th><th>BENZAI</th><th>CCODE</th><th>#sol</th><th>BENZAI</th><th>CCODE</th></tr></thead><tbody><tr><td>$n = 2$</td><td>1</td><td>0.00</td><td>0.01</td><td>1</td><td>0.01</td><td>0.01</td></tr><tr><td>$n = 3$</td><td>3</td><td>0.00</td><td>0.01</td><td>2</td><td>0.01</td><td>0.01</td></tr><tr><td>$n = 4$</td><td>7</td><td>0.05</td><td>0.03</td><td>5</td><td>0.01</td><td>0.02</td></tr><tr><td>$n = 5$</td><td>22</td><td>0.05</td><td>0.11</td><td>12</td><td>0.01</td><td>0.07</td></tr><tr><td>$n = 6$</td><td>81</td><td>0.14</td><td>0.22</td><td>36</td><td>0.07</td><td>0.15</td></tr><tr><td>$n = 7$</td><td>331</td><td>0.67</td><td>0.30</td><td>118</td><td>0.20</td><td>0.28</td></tr><tr><td>$n = 8$</td><td>1436</td><td>19.58</td><td>1.08</td><td>412</td><td>3.88</td><td>0.40</td></tr><tr><td>$n = 9$</td><td>6510</td><td>392.24</td><td>4.66</td><td>1492</td><td>19.78</td><td>2.21</td></tr><tr><td>$n = 10$</td><td>30129</td><td>22874.00</td><td>25.16</td><td>5587</td><td>732.75</td><td>7.90</td></tr></tbody></table>

5 Experimental Evaluation

In this section, we experimentally evaluate our approach on four different BGPs (see [2, 3, 4] for more details):

BGP1 enumerates hexagon graphs without single vertex holes (where a single vertex hole is a position without vertex which is surrounded by a cycle of 6 vertices);

BGP2 enumerates catacondensed graphs, i.e., hexagon graphs without single vertex holes and without cliques of order 3;

BGP3 enumerates hexagon graphs without single vertex holes and that contain a given pattern (taken in the list of eight patterns used in [3] and recalled in Appendix A);

BGP4 enumerates hexagon graphs without single vertex holes and that contain two given patterns (taken in the list of eight patterns used in [3] and recalled in Appendix A) so that the two patterns are not overlapping. (In [3], a variant is considered where patterns may share vertices; we do not display results for this variant as the conclusions are very similar.)

For each problem, the number $n$ of vertices is given, and we report results for $n \in [2, 10]$.

We consider the CP models described in Section 4, implemented in Choco [11], and available at https://gitlab.inria.fr/xipeng/bgp. To forbid single vertex holes, we exploit *pat* and *M* data structures. To forbid cliques of order 3, we ensure that, for each vertex *i*, *pat[i]* ∈ {P<sub>1+1+1</sub>, P<sub>1+1a</sub>, P<sub>1+1b</sub>, P<sub>1</sub>} as all other patterns contain a clique of order 3.

Our approach, denoted CCODE, is compared with the CP-based approach of [3], denoted BENZAI, and we used the Choco implementation available at https://github.com/benzAI-team/BenzAI. We do not compare our approach with the dedicated approach of [1] as it cannot enumerate hexagon graphs with holes (for BGP1, this approach finds less solutions than our approach), and it is not possible to add extra constraints in a declarative way so that it cannot be used to solve BGP2, BGP3, or BGP4 without modifying the code, or performing a post-processing to remove solutions that do not satisfy the extra constraints.

All experiments are carried out on an Intel Core Xeon E5-2623v3 of 3.0GHz×16 with 32GB of RAM.

In Table 1, we report CPU times of BENZAI and CCODE on BGP1 and BGP2. For BGP1, BENZAI is faster when $n \le 3$, but CCODE is faster when $n \ge 4$ and when $n = 10$ it is 824 times as fast. For BGP2, BENZAI is faster when $n \le 7$, but CCODE is faster when $n \ge 8$ and when $n = 10$ it is 93 times as fast. In Table 2, we report CPU times of BENZAI
<table><caption>**Table 2** Results of BENZAI and CCODE for BGP3 and BGP4. For each $n \in [2, 10]$ and each problem, we report the number of instances (#i) and the average, maximum and minimum time in seconds over these instances for BENZAI and CCODE. We report '-' when all runs exceed 3600s.</caption><thead><tr><th rowspan="2">n</th><th colspan="6">BGP3</th><th colspan="6">BGP4</th></tr><tr><th>#i</th><th colspan="3">BENZAI</th><th colspan="3">CCODE</th><th>#i</th><th colspan="3">BENZAI</th><th colspan="3">CCODE</th></tr><tr><th></th><th></th><th>avg</th><th>max</th><th>min</th><th>avg</th><th>max</th><th>min</th><th></th><th>avg</th><th>max</th><th>min</th><th>avg</th><th>max</th><th>min</th></tr></thead><tbody><tr><td>2</td><td>2</td><td>0.02</td><td>0.02</td><td>0.02</td><td>0.01</td><td>0.01</td><td>0.00</td><td>0</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>3</td><td>6</td><td>0.02</td><td>0.03</td><td>0.01</td><td>0.00</td><td>0.01</td><td>0.00</td><td>0</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>4</td><td>7</td><td>0.03</td><td>0.04</td><td>0.02</td><td>0.00</td><td>0.02</td><td>0.00</td><td>3</td><td>0.07</td><td>0.16</td><td>0.02</td><td>0.01</td><td>0.03</td><td>0.00</td></tr><tr><td>5</td><td>8</td><td>0.05</td><td>0.07</td><td>0.01</td><td>0.01</td><td>0.07</td><td>0.00</td><td>11</td><td>0.06</td><td>0.20</td><td>0.04</td><td>0.01</td><td>0.06</td><td>0.00</td></tr><tr><td>6</td><td>8</td><td>0.38</td><td>0.43</td><td>0.23</td><td>0.03</td><td>0.15</td><td>0.00</td><td>23</td><td>1.41</td><td>2.51</td><td>0.93</td><td>0.02</td><td>0.16</td><td>0.01</td></tr><tr><td>7</td><td>8</td><td>0.95</td><td>1.06</td><td>0.79</td><td>0.13</td><td>0.35</td><td>0.01</td><td>29</td><td>3.78</td><td>5.42</td><td>2.30</td><td>0.10</td><td>0.50</td><td>0.01</td></tr><tr><td>8</td><td>8</td><td>20.33</td><td>23.84</td><td>13.94</td><td>0.73</td><td>1.61</td><td>0.04</td><td>34</td><td>65.06</td><td>83.21</td><td>50.28</td><td>0.49</td><td>1.69</td><td>0.03</td></tr><tr><td>9</td><td>8</td><td>306.30</td><td>399.25</td><td>75.05</td><td>3.94</td><td>10.52</td><td>0.23</td><td>35</td><td>400.51</td><td>574.15</td><td>226.87</td><td>2.94</td><td>10.85</td><td>0.20</td></tr><tr><td>10</td><td>8</td><td>-</td><td>-</td><td>-</td><td>34.60</td><td>69.75</td><td>2.20</td><td>36</td><td>-</td><td>-</td><td>-</td><td>24.19</td><td>71.03</td><td>1.61</td></tr></tbody></table>

and CCODE on BGP3 and BGP4. CCODE is always faster. When $n = 9$, it is 78 (resp. 136) times as fast as BENZAI for BGP3 (resp. BGP4). When $n = 10$, BENZAI is not able to solve any of the 8 (resp. 36) instances of BGP3 (resp. BGP4) within one hour whereas CCODE never exceeds 71s.

# 6 Conclusion

We have introduced a new canonical code for representing hexagon graphs while being invariant to rotations and symmetries. This canonical code is the smallest code associated with a BFS of the graph, and we have shown how to define codes by means of constraints, thus allowing us to generate all consistent codes with CP. We have shown that every prefix of a canonical code is canonical and this property allows us to trigger a failure whenever the current code prefix is not canonical.

We have also shown how to efficiently enumerate all codes that start with a given prefix, thus ensuring that a given pattern is included in the generated graphs. However, in this case it may happen that the prefix of a canonical code is no longer canonical, due to the fact that patterns may specify forbidden vertices. In this case, we exhibit a weaker property which allows us to filter some codes that cannot be extended to a canonical code.

We have experimentally compared our approach with the CP-based approach of [3], and we have shown that it scales much better, especially when there are a lot of solutions to enumerate.

As our approach is based on CP, new properties that must be satisfied by the generated graphs are defined in a declarative way by means of constraints. To this aim, the user has access to four variables for each vertex $i$, i.e., $p_i$, $d_i$, $x_i$, and $y_i$. In some cases, it may be convenient to be able to add constraints on vertex patterns. This is the case, for example, for forbidding cliques of order 3 (in BGP2). In our first implementation, which is a proof of concept, the pattern of a vertex is maintained with a backtrackable data structure (using `IStateInt` Choco objects). As future work, we plan to introduce integer variables that represent patterns in order to allow the user to easily define new constraints on patterns.

We also plan to improve the propagation of the global constraint *subgraph* in order to detect some inconsistencies earlier: we could adapt Algorithm 2 (seq) in order to search for the largest sequence in polynomial time, and trigger a failure whenever the size of this sequence plus the number of non assigned vertices is smaller than $n$.
Finally, we would like to extend our approach to other kinds of graphs such as grid graphs (where vertices are associated with square faces instead of hexagonal faces), for example. We hope this will pave the way for new applications of CP such as, for example, applications that search for patterns in cellular automata as in [7].

References

1 Gunnar Brinkmann, Gilles Caporossi, and Pierre Hansen. A constructive enumeration of fusenes and benzenoids. *J. Algorithms*, 45(2):155-166, 2002.

2 Yannick Carissan, Denis Hagebaum-Reignier, Nicolas Prcovic, Cyril Terrioux, and Adrien Varet. Using constraint programming to generate benzenoid structures in theoretical chemistry. In Helmut Simonis, editor, *Principles and Practice of Constraint Programming - 26th International Conference, CP*, volume 12333 of *Lecture Notes in Computer Science*, pages 690-706. Springer, 2020.

3 Yannick Carissan, Denis Hagebaum-Reignier, Nicolas Prcovic, Cyril Terrioux, and Adrien Varet. Exhaustive generation of benzenoid structures sharing common patterns. In Laurent D. Michel, editor, *27th International Conference on Principles and Practice of Constraint Programming, CP*, volume 210 of *LIPIcs*, pages 19:1-19:18. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2021.

4 Yannick Carissan, Denis Hagebaum-Reignier, Nicolas Prcovic, Cyril Terrioux, and Adrien Varet. How constraint programming can help chemists to generate benzenoid structures and assess the local aromaticity of benzenoids. *Constraints An Int. J.*, 27(3):192-248, 2022.

5 S. J. Cyvin, J. Brunvoll, and B. N. Cyvin. Search for concealed non-kekulean benzenoids and coronoids. *Journal of Chemical Information and Computer Sciences*, 29(4):236-244, 1989.

6 Romain Deville, Élisa Fromont, Baptiste Jeudy, and Christine Solnon. Grima: A grid mining algorithm for bag-of-grid-based classification. In Antonio Robles-Kelly, Marco Loog, Battista Biggio, Francisco Escolano, and Richard C. Wilson, editors, *Structural, Syntactic, and Statistical Pattern Recognition - Joint IAPR International Workshop, S+SSPR, Proceedings*, volume 10029 of *Lecture Notes in Computer Science*, pages 132-142, 2016.

7 Romain Deville, Élisa Fromont, Baptiste Jeudy, and Christine Solnon. Mining frequent patterns in 2d+t grid graphs for cellular automata analysis. In Pasquale Foggia, Cheng-Lin Liu, and Mario Vento, editors, *Graph-Based Representations in Pattern Recognition - 11th IAPR-TC-15 International Workshop, GbRPR 2017, Anacapri, Italy, May 16-18, 2017, Proceedings*, volume 10310 of *Lecture Notes in Computer Science*, pages 177-186, 2017.

8 Stéphane Gosselin, Guillaume Damiand, and Christine Solnon. Efficient search of combinatorial maps using signatures. *Theor. Comput. Sci.*, 412(15):1392-1405, 2011.

9 Brendan D. McKay and Adolfo Piperno. Practical graph isomorphism, II. *J. Symb. Comput.*, 60:94-112, 2014.

10 Siegfried Nijssen and Joost N. Kok. The gaston tool for frequent subgraph mining. In *Proceedings of the 2nd International Workshop on Graph-Based Tools, GraBaTs 2004, Rome, Italy, October 2, 2004*, volume 127 of *Electronic Notes in Theoretical Computer Science*, pages 77-87. Elsevier, 2004. doi:10.1016/j.entcs.2004.12.039.

11 C. Prud'homme, J.-G. Fages, and X. Lorca. *Choco Documentation*. TASC, INRIA Rennes, LINA CNRS UMR 6241, COSLING S.A.S., 2016. URL: http://www.choco-solver.org.

12 Xifeng Yan and Jiawei Han. gspan: graph-based substructure pattern mining. In *2002 IEEE International Conference on Data Mining, 2002. Proceedings.*, pages 721-724, 2002. doi:10.1109/ICDM.2002.1184038.
A Patterns used in BGP3 and BGP4

The eight patterns used in BGP3 and BGP4 are displayed below (mandatory and forbidden vertices are displayed in black and red, respectively).

<figure><img src="image_285.png" /></figure>
<figure><img src="image_286.png" /></figure>
# Distribution Optimization in Constraint Programming

Guillaume Perez

Huawei Technologies Ltd, CSI Paris, Boulogne-Billancourt, France

Gaël Glorian

Huawei Technologies Ltd, CSI Paris, Boulogne-Billancourt, France

Wijnand Suijlen

Huawei Technologies Ltd, CSI Paris, Boulogne-Billancourt, France

Arnaud Lallouet

Huawei Technologies Ltd, CSI Paris, Boulogne-Billancourt, France

## Abstract

Stochastic Constraint Programming introduces stochastic variables following a probability distribution to model uncertainty. In the classical setting, probability distributions are given and constant. We propose a framework in which random variables are given a set of possible distributions and only one should be selected. A solution is obtained when all variable distributions are assigned, and all decision variables are assigned too. In such a setting, a constraint on random variables limits the possible distributions its random variables may take. We generalize the notion of *chance* as the probability of satisfaction of a constraint, called *probabilization*, given variable distributions. Probabilization can be seen as a generalization of reification in a random setting whose result is a random variable. We define minimal arithmetic to work with stochastic variables having a variable distribution. Using the introduced representation, our framework can in theory save an exponential number of decisions, and represents problems that were previously not representable with finite integer domains. Finally, we model and solve two industrial problems that require this extension – virtual network configuration and assignment of chemical delivery – and show improvement in terms of quality of solution and speed.

**2012 ACM Subject Classification** Theory of computation → Constraint and logic programming; Mathematics of computing → Solvers; Computing methodologies → Probabilistic reasoning

**Keywords and phrases** Constraint Programming, Optimization, Stochastic Optimization

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

## 1 Introduction

Stochastic optimization and chance-constrained programming [6, 41] are classes of problems in which uncertainty is present. In such a setting, both decision variables and random variables are present. Usually the probability distributions of the random variables are known and constant. The goal is to optimize a given objective function on these variable sets and to satisfy a set of constraints with sufficient probability. Optimal production planning, optimal power flow, textile manufacturing, vehicle sharing, and parcel delivery services are just a few of the many industrial areas where stochastic optimization is required [28, 29, 34, 55, 21].

Learning probabilistic distributions, or distribution learning, is a machine learning framework that consists of learning the probability distribution that could generate a given set of samples [22, 24]. In the usual settings, the input is a set of samples drawn from an unknown distribution and the goal is to uncover this unknown distribution. Since, many methods have been proposed, for example using assumptions on the class of probability distribution to be a mixture of Gaussian or a Poisson law etc. [10, 9], or directly using neural networks [1].

<figure><img src="image_287.png" /></figure>
Consider the virtual network functions design problem [11, 13, 15, 47]. The main part of this problem consists of selecting the settings of different nodes of a network function, such that the global computation latency is robust. Robustness here implies that the total latency (i.e. the sum of the latency of each node) is smaller than a given value, with a probability of at least $\gamma$. Locally, for each node the latency distribution is given by a random variable whose distribution is conditioned by the settings of the node. Different settings will lead to different probability distributions for the nodes.

In this paper, we propose to work on *probability distribution optimization*. In this setting, the input is a stochastic constraint optimization problem, and the goal is to find both the assignment of the decision variables and the distribution of the random variables. This can be seen as a generalization of distribution learning in the sense that it is not restricted to the fitting constraints.

Constraint Programming (CP) is an expressive optimization framework often used to solve combinatorial problems such as scheduling. In CP, optimization under chance, confidence, probability, or statistical constraints is a prolific research area [29, 39, 33, 27]. The early and impacting works on stochastic constraint programming defined the basics [14, 53]. Then, optimization methods for chance constraints, Markov, or sampling probability distribution constraints etc. have been proposed [48, 18, 46, 39]. These works have focused either on the multi-stage framework or on one global constraint to extend the optimization process to a stochastic context. But they all consider random variables with fixed distributions.

In this paper, the stochastic CP framework is extended to handle distribution optimization. It is another set of stochastic problems where random variables are given a domain of possible distributions, from which only *one* should be selected. The abstract object *probability distribution variable* is introduced for modeling purposes. It represents the variability of its sample space and its probability. A random variable is *assigned* when its associated probability distribution is known and fixed. *Distribution variables* are now one of the many variables of the problem to solve. A solution is found when all the decision variables are assigned, and when the distribution of all the random variables is fixed. In addition, we propose the definition of *distribution constraints*, which are constraints involving distribution variables. In the *hard* case, a distribution constraint restricts the possible distributions of the random variables it involves. Using the introduced representation, our framework can in theory save an exponential number of decisions, and represents problems that previously could not be represented with a finite integer domain.

Then, we propose to focus on a particular case of constraints namely *relational constraints*. These constraints represent the usual CP constraints, as they are defined by the set of allowed tuples. For these constraints, we propose two new consistency levels, namely *P* consistency, for probability consistency, and *Ω* consistency, which is a probability distribution encoding of the usual arc consistency. A direct implication is that most existing stochastic constraints (confidence, sampling, PMF, etc.) that consider the distribution as data of the problem, can be upgraded to deal with variable distributions. That is why in this paper we extend the Confidence/Chance of relational constraints [49, 30, 36]. Furthermore, the notion of *probabilization* of a constraint is proposed. It is a generalization of the chance concept [49], closely related to the reification of constraint, but for probability purpose. The *probabilization* returns a Boolean random variable associated to a distribution variable representing the probability of satisfaction of the constraint.

Distribution variables are abstract objects that can be implemented in diverse ways. We propose an implementation of distribution variables using a direct encoding as a probability mass function. For each possible value of a random variable, we define its probability as a continuous variable. Then, we propose filtering rules defining the minimal arithmetic allowing to model distribution optimization problems.
Finally, we exploit this implementation to solve two distribution optimization problems in the experimental section. The first problem is the stochastic design of virtual network functions, where a simple, yet efficient, model allows to find optimal designs. The second problem we model is a chemical delivery application. It shows how using *variable distributions* allows to have smaller search trees and higher quality solutions.

## 2 Related Work

Chance constrained optimization [6, 31, 33] is a prolific research area. They are classes of problems in which uncertainty is present. In such setting, both decision variables and random variables are present. Usually the probability distributions of the random variables are known and constant. The first definition was given for discrete distributions and piece-wise linear functions with linear inequalities involving random variables. These inequalities had to be maintained at a given level of probability. In these problems, the random variables can be defined using a joint probability distribution, making them conditional on each other [8], which usually makes the problem harder to solve. This framework has been used for the case where the distribution of random variables was conditioned on storage levels or stream flows [20]. Moreover, joint chance constrained optimization problem are problems that contain multiple uncertain constraints on multivariate random variables. They are jointly required to be satisfied with probability exceeding a threshold [54]. However, it is uncommon to have the probability distribution defined as conditioned by the decision variables, or even constrained. The work proposed in this paper is part of chance constrained optimization and aims to use a constraint programming solver as a modeling and solving framework for problems where finding a distribution that satisfies constraints is also part of the problem.

The probabilistic graphical model community proposed Bayesian networks and influence diagram framework [35, 26] to represent the interactions between decision variables (agents) and the probability distributions of random variables. Such influence diagrams, while they tend to grow exponentially, could be powerful tools to represent the constrained distribution during the optimization process. Several methods have been defined to solve chance constrained optimization. In general, non-linear programming solvers are used [28] and dedicated models are defined. For example, dynamic programming coupled with influence diagram has been used to optimize agent decisions to maximize utility functions [50], in a similar way as reinforcement learning. Furthermore, it is not unusual to use sampling methods to approximate chance optimization [3, 33].

In CP, the *chance* constraint is defined as a policy under uncertainty, guaranteeing robustness of the assignment [53, 49]. More precisely, consider a problem where first a set $X_1$ of decision variables are assigned. Next, some random variables $Y_1$ reveal their value, then values are selected for the set of variables $X_2$ and so on for as many iterations as it is required. The solution of such problem is a tree, and no longer a tuple, for which the specified fraction of the scenarios is satisfied. Since then, this work has been extended multiple times, [17, 18] proposed generic algorithms reusing existing filtering for the global chance constraint. More recently, [30] proposed to restrict the chance constraint to 2 stages optimization, as it is a widely used approach. They proposed a filtering for the conjunction of binary inequalities. Finally, [36] proposed to use MDDs as a new generic filtering for the confidence constraint. All of these are different from the work proposed in this paper. On the one hand, policies are not considered, the proposed framework is restricted to 2-stages similarly to [30, 36]. On the other hand, the proposed framework is the first introducing variable distribution inside of constraint programming. Note that it is one of the possible extensions to stochastic CP
as mentioned in the *Extensions* section of [53]. Finally, a recent trend called randomness optimization [23, 2, 15] is used to optimize the randomness of the solution of a decision problem. For example randomness of solution is important for the design of systems to prevent reverse engineering. Random solution as a constraint is already part of the constraint programming framework [42].

# 3 Distribution Optimization

## 3.1 Preliminaries

A sample space $\Omega$ is a set representing all the possible outcomes of an experiment. In the general case, the probability of an outcome $v \in \Omega$ is noted $\text{Pr}[v]$. An event is a set of outcomes from the sample space and an event space $\mathcal{F}$ is a set of events. Events are often useful for characterising particular subsets of outcomes. For example, the latency of a process might be a real number, but the probability will be usually defined by segments of time (i.e. intervals). Finally, $P$ is a probability function that maps each element $e$ of $\mathcal{F}$ to a probability $P(e) \in [0, 1]$. These three elements form a probability space, denoted by the triplet $\psi = (\Omega, \mathcal{F}, P)$. A random variable $y$ follows a probability space $\psi$ if its distribution is defined by the probability space $\psi$. This will be noted $y \sim \psi$. When $\Omega$ is discrete or countable, a probability mass function (PMF) $\text{pmf}: \Omega \to \mathbb{R}$ assigns a probability $\text{Pr}[v]$ to each value $v \in \Omega$. Let $Y = (y_1, ..., y_r)$ be a vector of $r$ independent random variables. By definition, the probability of a sample of independent variables (i.e. a tuple) is given by the product of the probabilities of the selected values. More precisely, the probability of a tuple $t = (a_1, ..., a_r)$ is defined by: $\text{Pr}[Y = t] = \prod_{i=1}^{r} \text{Pr}[y_i = a_i]$. When a subset $Y$ is composed of variables that are not independent, it is possible to replace them by a single random variable $Y^*$ representing their aggregation [7]. The domain of outcomes of this variable is the Cartesian product of the outcomes of the variables of $Y$. So we can assume in this paper that all random variables are independent.

In the rest of this paper, $X = \{x_1, ..., x_r\}$ denotes a set of variables (usually integer). $Y = \{y_1, ..., y_k\}$ denotes a set of random variables. Given a constraint $C$, $T_C(v_1, ..., v_n) = 1$ (resp. = 0) implies that constraint $C$ is satisfied by the tuple $(v_1, ..., v_n)$ (resp. unsatisfied). $D(x)$ denotes the current domain of variable $x$. We denote by $\bar{D}(x)$ (resp. $D(x)$) the upper bound (resp. lower bound) of variable $x$. We denote by $v \in D(x)$ the fact that a value $v$ belongs to the domain of a variable $x$. In a similar way, let $\forall t \in D^\times(X) = \prod_{i=1}^n D(X_i)$ denotes all the tuples in the Cartesian product of the current domain of variables in $X$.

## 3.2 Distribution Variables

In this section, we extend the CP framework with the definition of probability distribution variables.

**► Definition 1.** A probability distribution variable $r$ is a variable defining a probability space. The domain $D(r)$ of a distribution variable $r$ is a set of probability spaces.

Recall that a probability space defines its sample space $\Omega$, its event space $\mathcal{F}$ and the probability $P$ of each event in it. *Distribution variables* are now one of the many variables of the problem, similar to any other decision variable. A probability distribution variable is assigned when the probability of each value in its event space is known and fixed. Let $r$ be a random distribution variable. $\psi \in D(r)$ denotes the membership of a probability space $\psi$ to the current domain of the distribution variable $r$. We denote by $y \sim r$ the
fact that the random variable $y$ will be drawn following the distribution variable $r$, and by $y \sim \psi$ that random variable $y$ follows the distribution defined by the probability space $\psi$ (a constant). Let $y_1$ and $y_2$ be two random variables such that $y_1 \sim r$, $y_2 \sim r$. Such notation implies that they are independent and identically distributed (IID). We denote $D(y)$ or $D(r)$ the domain of the probability distribution variable of random variable $y \sim r$. Let $\forall \Psi \in D^\times(Y)$ denotes all the tuples of probability spaces in the Cartesian product of the current domain of distribution variables of $Y$. Let $\Omega_y$ be the set of all the outcomes such that there exists a probability space in $D(y)$ where the outcome has a non-zero probability. Let $\forall t_Y \in D_\Omega^\times(Y)$ denote all the tuples of outcomes (i.e. value) in the Cartesian product of the current sample space $\Omega_y$ of random variables in $y \in Y$. We note by $\overline{Pr[Y=t]}$ (resp. $\overline{Pr[Y=t]}$) the upper bound (resp. lower bound) probability of tuple $t$ to be drawn by random variables in $Y$. We have $\overline{Pr[Y=t]} = \max\{Pr[Y \sim \Psi = t] | \Psi \in D^\times(Y)\}$ and $\overline{Pr[Y=t]} = \min\{Pr[Y \sim \Psi = t] | \Psi \in D^\times(Y)\}$. All along the paper, $\Psi$ represents a vector of probability spaces, $\psi$ a probability space, and $\Psi[\{y\}]$ the probability space associated to variable $y$ in $\Psi$. Finally, $\mathbb{P}$ is the space of all the possible probability spaces (i.e. $\forall \psi \in \mathbb{P}$).

It is interesting to note that integer variables can be represented by the proposed distribution variables.

**Proposition 2.** *A integer variable always has an equivalent distribution variable.*

*Proof.* Let $x$ be an integer variable with domain $D(x)$. Let $r$ be a distribution variable such that $D(r) = \{(\{v\}, \{\emptyset, \{v\}\}, \text{Pr}[v]=1) | \forall v \in D(x)\}$. In other words, for each value $v$ in the domain of $x$, there is a probability space whose sample space $\Omega = \{v\}$ is restricted to value $v$. There is a one-to-one correspondence between the values in the domains of $r$ and $x$, hence $x$ and $r$ are equivalent. Note that the converse is not true as the number of probability spaces in a distribution variable may not be countable. In the rest of this paper, without loss of generality, we assume that all the integer variables have been replaced by equivalent distribution variables for simplicity of notation. ◀

## 3.3 Distribution Constraints

A random distribution variable being a modeling object of CP, it can be constrained. In this section, we propose a basic definition of distribution constraints. We use this definition to derive a generalization of arc consistency to distribution variables. Finally, we propose to focus on a sub-part of the constraint space, namely relational constraints, that generalizes usual CP constraints.

Distribution constraints are all the constraints involving random variables with variable distribution. A distribution constraint defines the possible probability spaces of the random variables it involves. Integer constraints can always be defined by the set of satisfying tuples. This basic notion is extended to distribution constraints:

**Definition 3.** Let $C$ be a constraint defined on random variables $Y$. $C$ is defined by a set of tuples $T_R$ such that $\forall \Psi \in T_R$, $\Psi$ is a tuple of valid probability spaces for the probability distribution of random variables in $Y$.

**Example 4.** Table of PMFs The following table represents a constraint on two variables: integer variable $x$ and random variable $y \sim r$.

<table><thead><tr><th>Pr[x = 0]</th><th>Pr[x = 1]</th><th>Pr[y = 1]</th><th>Pr[y = 2]</th><th>Pr[y = 3]</th><th>Pr[y = 4]</th><th>Pr[y = 5]</th><th>Pr[y = 6]</th></tr></thead><tbody><tr><td>1</td><td>0</td><td><sup>1</sup>&frasl;<sub>6</sub></td><td><sup>1</sup>&frasl;<sub>6</sub></td><td><sup>1</sup>&frasl;<sub>6</sub></td><td><sup>1</sup>&frasl;<sub>6</sub></td><td><sup>1</sup>&frasl;<sub>6</sub></td><td><sup>1</sup>&frasl;<sub>6</sub></td></tr><tr><td>0</td><td>1</td><td><sup>1</sup>&frasl;<sub>3</sub></td><td><sup>1</sup>&frasl;<sub>3</sub></td><td><sup>1</sup>&frasl;<sub>12</sub></td><td><sup>1</sup>&frasl;<sub>12</sub></td><td><sup>1</sup>&frasl;<sub>12</sub></td><td><sup>1</sup>&frasl;<sub>12</sub></td></tr></tbody></table>
29:6 Distribution Optimization in CP

Such a table restricts the distribution variables of x and y to two probability spaces, defined by probability mass functions. x is a integer variable encoded as a probability distribution. The possible probabilities are either 0 or 1. y is a true random variable. Depending on the assigned value of x, random variable y will behave differently. Note that PMFs are not the only way to encode probability spaces.

The notion of consistency is of utmost important in constraint programming [45]. We propose to directly map the generalized arc consistency to distribution variables. The principal difference is that the domain of a distribution variable is a set of probability spaces, and not a set of integers as for integer variables for example. Let C be a constraint defined on random variables Y defined by a set of tuples of probability spaces $T_R$. The consistency properties are defined by:

- A probability space $\psi \in D(y)$ of $y \in Y$ is consistent with C if it exists a support $\Psi \in T_R$ such that $\Psi[\{y\}] = \psi$ and $\Psi \in D^\times(Y)$.

- A distribution variable of $y \in Y$ is consistent if each probability space in its domain is consistent.

- A constraint C is generalized arc consistent if all its variables are consistent.

3.4 Relational Constraints

Relational constraints are a subset of distribution constraints. They consider the sampled values of the random variables for the satisfiability check of the constraint. In other words, the constraint can always be defined by the possible assignments for the integer variables and for the random variables. They are composed of all the usual constraints of the CP framework but generalized to variable distributions.

**► Definition 5.** Let C be a relational constraint defined on random variables Y. C is defined by a table of tuples $T_C$, such that $\forall t \in T_C$ is a valid sample that can be drawn for random variables in Y.

Note that, for each table $T_C$, an infinite number of probability spaces are valid. Indeed, consider the binary constraint defined by the tuples [(0, 0), (0, 2)] on variables x and y. This constraint restricts the sampled value of x to be 0, and restricts the probability distribution of y to any distribution such that $Pr[y = 0] + Pr[y = 2] = 1$. This implies that such a definition can be highly compressing, as an infinite number of distributions satisfy this equality. On the other hand, such a definition is not expressive enough to restrict to all subsets of probability spaces. For example, it is impossible to define $E(y) = 3$ (expectation), or simply the constraint defined by the table from the example in section 3.3. That is the reason why relational constraints are a subset of distribution constraints.

**► Definition 6 (Relational Constraint).** Let $T_C(t) = 1$ if $t \in T_C$ and 0 otherwise be a predicate function for constraint C. We propose to use the following definition for the relational constraint:

$$ \int_{t \in D_{\Omega}^{\times}(Y)} Pr[Y = t]T_C(t) = 1 \quad ; \quad \int_{t \in D_{\Omega}^{\times}(Y)} Pr[Y = t] \neg T_C(t) = 0 \qquad (1) $$

The most intuitive description of such constraint is that it does not imply that all valid tuples must be reached, but that all reachable tuples must be valid.

The notion of consistency of a probability space in the case of a relational constraint can be redefined using its definition. Given a relational constraint C defined on random variables Y by the relational table $T_C$.
**► Definition 7.** A probability space $\psi \in D(y)$ of $y \in Y$ is consistent with $C$ if and only if:

$$ \exists \Psi \in D^\times(Y), \Psi[\{y\}] = \psi \wedge \forall t_Y \in D_\Omega^\times(Y), \Pr[Y \sim \Psi = t_Y] \le T_C(t) \quad (2) $$

It is a direct rewriting of the consistency of distribution constraints using the relational definition (equation (1)). Then, the domain of probability variable $y \in Y$ is consistent with constraint $C$ if all the probability spaces in it are consistent with $C$. Finally, a constraint $C$ is consistent if all the variables in it are consistent.

Let $\psi$ be a probability space and $Y_y^\psi$ be a random distribution variable vector where variable $y$ is restricted to probability space $\psi$. A weaker consistency can be defined for relational constraints and distribution variables. This consistency is close to bound-consistency. It considers the bounds of the probability of events given the current domain of a distribution variable.

**► Definition 8.** A probability space $\psi \in D(y)$ of $y \in Y$ is P-bound consistent with $C$ if:

$$ \forall t \in D_\Omega^\times(Y_y^\psi), \Pr[Y_y^\psi = t] \le T_C(t) \quad (3) $$

Then, the domain of probability variable $y \in Y$ is P-bound consistent with constraint $C$ if all the probability spaces in it are P-bound consistent with $C$. Finally, a constraint $C$ is P-bound consistent if all the variables in it are P-bound consistent: $\forall t \in D_\Omega^\times(Y), \Pr[Y = t] \le T_C(t)$. Now that the P-bound consistency is defined, we propose to focus on the case $T_C(t) = 0$ as it is the only one able to invalidate probability spaces. Let $Y_{\setminus\{y\}}$ be the vector of variables $Y$ without variable $y$. Let $t_{\setminus\{y\}}$ be the tuple $t$ without the value at the position of variable $y$. Propagating the following logical implication is enough to enforce P-bound consistency constraint $C$:

$$ \forall y \in Y, \forall t \in D_\Omega^\times(Y), \Pr[Y_{\setminus\{y\}} = t_{\setminus\{y\}}] > T_C(t) \implies \Pr[y = t_y] = 0 \quad (4) $$

This pruning rule implies that all the probability spaces that could lead to an invalid tuple should be removed. It will later be derived to extract filtering rules for the implementation of this framework.

**► Definition 9.** Given a relational constraint $C$ defined on random variables $Y$ by the relational table $T_C$. A probability space $\psi \in D(y)$ of $y \in Y$ is $\Omega$-consistent with $C$ if:

$$ \exists t \in D_\Omega^\times(Y_y^\psi), T_C(t) \quad (5) $$

$\Omega$-consistency ensures that it exists at least one tuple with non-zero probability to satisfy the predicate given the other domains. $\Omega$-consistency, when all the variables are integer variables encoded as distribution variables, is equivalent to global arc consistency as it finds a support for each value in the domain of each variable.

Consider a total order on the outcomes of $\Omega_y$. Let $\overline{\Omega_y}$ (resp. $\overline{\Omega_y}$) be the largest (resp. smallest) element of $\Omega_y$ with respect to the total order. Let $D_\Omega^\times(Y)$ denote the Cartesian product of the intervals $[\overline{\Omega_y}, \overline{\Omega_y}]$, $\forall y \in Y$.

**► Definition 10.** A probability space $\psi \in D(y)$ of $y \in Y$ is $\Omega$-bound-consistent with $C$ if:

$$ \exists t \in D_\Omega^\times(Y_y^\psi), T_C(t) \quad (6) $$

When all the variables are integer variables encoded as distribution variables and the total order is defined by the operator $<$, $\Omega$-bound-consistency is equivalent to bound-consistency.

Consider the Confidence constraint that restricts the decision variables to values such that the probability of satisfaction is greater than a given threshold [49, 30, 36]. We propose to generalize it by considering variable distributions for random variables.
**► Definition 11.** Given a relational constraint $C(Y)$, with $Y$ a vector of random variables with variable distributions. Given a confidence threshold $\gamma$. The generalized confidence is:

$$ \int_{t \in D_{\Omega}^{\times}(Y)} \Pr[Y = t] \mathbb{T}_{C}(t) > \gamma \qquad (7) $$

Previously, propagation of chance constraint was only removing values from integer variables that could not lead to a robust enough solution. The new propagation of this constraint will impact both the integer variables (here encoded as distributions) and the distribution variables as probability spaces that are not robust enough will be removed too.

## 3.5 Views on Random Variables

We propose to make a clear distinction between the three equality operators: “$y_1 = y_2$”, “$y_1 \simeq y_2$”, and “$y_1 \equiv y_2 + y_3$”. First, $y_1 = y_2$ is the classic equality constraint, it implies that the value of the random variables must be equal, but they are considered different random variables. Second, the constraint $y_1 \simeq y_2$ ensures that $y_1 \sim r_1$ and $y_2 \sim r_2$ are equal in law. More precisely it ensures that $r_1 = r_2$. This implies that their random assignment might be different, but the constraint is satisfied if they follow the same distribution. From a probability point of view, they are independent. Finally, $y_1 \equiv y_2 + y_3$, where $y_1$ is the actual summation of the two random variables $y_2$ and $y_3$. From a probability point of view, they are dependent.

Views are a useful abstraction in constraint programming [47, 51]. They prevent the creation of intermediate variables representing some function on a variable ($x' = f(x)$). Views can be adapted to random variables directly. Consider for example the affine view $y_1 \equiv ay + b$. Iterating on the values of $y_1$ is done by iterating on the values of $y$ and applying the affine transformation. In addition, $\Pr[y = i] = \Pr[y_1 = ai + b]$. In the rest of this paper, notations such as $-y$ represent a view on random variable $y$.

## 3.6 Probability Valuation of Constraints

In constraint programming, the reification of a constraint is very useful while modeling a problem. The reification variable is usually a Boolean variable denoting the satisfiability of the related constraints. Other valuations of constraints have been proposed, such as soft constraints [52, 25] and cost-version of constraints [44]. In this paper, we introduce the probability valuation of a constraint. It is a generalization of the chance from [18] for the case of variable distributions as the result is a Boolean random variable.

**► Definition 12.** Let $C$ be a relational constraint defined on the random variables $Y$. Let $p_C$ be a Boolean random variable. The probability valuation $p_C \equiv C(Y)$ is defined by:

$$ \Pr[p_C = 1] = \int_{t \in D_{\Omega}^{\times}(Y)} \Pr[Y = t] \mathbb{T}_{C}(t) \qquad (8) $$

Note that in the general case processing the probability of satisfaction of a constraint $C$ is hard because the table $T_C$ of the constraint may be untractable.

**► Proposition 13.** The probabilization of a constraint is a generalization of the reification.

**Proof.** Let $C$ be a relational constraint defined on the random variables $Y$. Let $p_C \in [0, 1]$ be a Boolean random variable such that $\Pr(p_C = 1) = \Pr[C(Y)]$. The two possible values for the reification of a constraint are satisfied and unsatisfied. If $\Pr(p_C = 1) = 1$ the constraint is always satisfied, if $\Pr(p_C = 1) = 0$, the constraint is always unsatisfied. For all the other values, $p_C$ represents the probability of satisfaction of the constraint. ◀
Note that the probabilization can be seen as a random variable representing the projection of the distribution of the random variables on the constraint. This implies that $p_C$ and the random variables in $C$ are dependents. In the rest of this paper, $p_C$ will be used to denote $\Pr[p_C = 1]$ when no ambiguity is present.

## 3.7 Multiple Stochastic Relational Constraints

In chance-constrained optimization, applications can be made of several independent parts [54], and the joint probability must be satisfied. Those are problems for which constraints of the problem can be split into independent ones. Let **1** be the vector of ones. More formally, these problems can be reformulated as:

$$ \arg\min_{Y \in \mathbb{D}} F(Y) \qquad (9) $$

$$ \text{subject to} \quad p_1 = C_1(Y_1) \qquad (10) $$

$$ \dots \qquad (11) $$

$$ p_n = C_n(Y_n) \qquad (12) $$

$$ \Pr[(p_1, \dots, p_n) = \mathbf{1}] \geq \gamma \qquad (13) $$

Where $Y_i \cap Y_j = \emptyset, \forall i, j$, and $\cap$ being the set intersection operator. Using the probability valuation of constraints (10-12), the global chance or confidence constraint is (13). This implies that the hard part of the chance constraining will be located into the probability valuation of the constraints.

# 4 Implementation

## 4.1 Distribution Implementation

In the rest of this paper, a possible implementation for the discrete random variables is proposed. Then, propagation algorithms are proposed using the proposed implementation. Implementing a variable distribution random variable inside a solver might be challenging, as probability distributions can be expressed in very diverse ways [26]. Indeed, as for set or sequence variables [11], the choice of implementation will impact the design of propagation algorithms, and the capability of representation. For known probabilistic distributions, a straightforward implementation would be the parameters that describe it. For example, it would be sufficient to have two continuous variables $\mu$ and $\sigma$ for a normal distribution. Indeed, with such a definition, the propagation could use the closed form of the cumulative distribution function and probability density functions. In this paper, for the implementation of the distribution, we assume that random variables are independent and discrete. We propose to represent distribution variables using probability mass function variables PMF.

**Definition 14.** Let $y$ be a discrete independent random variable with sample space $\Omega_y = \{v_1, \dots, v_d\}$. The PMF variable of $y$ is represented as a vector $d^y = (d_v^y)_{v \in \Omega_y}$ of continuous variables where $D(d_v^y) \subseteq [0, 1]$ represents the probability $\Pr[y = v]$ and $\sum_i d_v^y = 1$.

Note that using continuous or floating point variables inside of CP solver is a known topic [5, 43]. They are usually represented by intervals. For evident reasons, integer variables are considered encoded as usual.

**Example 15.** Consider a Constraint Satisfaction Problem (CSP) with one decision variable $x \in \{0, 1\}$ and one random variable $y \in [1, 6]$. The CSP contains two constraints: one distribution constraint and one relational constraint. The first one defines the possible
distributions of $y$ and is defined by the table from the example in section 3.3. The second constraint is the constraint $\text{Pr}[y > 2] \ge \gamma$ with $\gamma = 0.6$. A solution of such a CSP is an assignment of both $x$ and the $d_i^y$ variables that satisfies both constraints.

Consider the propagation of the $y > 2$ constraint. For each of the variables $d_3^y$ to $d_6^y$, the lower bound is set to $\gamma - 3\frac{1}{6} = 0.1$. Indeed, for each of the possible values, the cumulative probability must be at least equal to $\gamma$. Once the propagation of constraint $y > 2$ is done, the distribution constraint can filter out value 1 from $x$ as the associated distribution is no longer valid ($\underline{d_4^y} = 0.1 > \frac{1}{12}$). Values of $d^y$ are assigned to the tuple $(\frac{1}{6}, \dots, \frac{1}{6})$, and a solution is found.

**Example 16 (Probability valuation).** Consider again the problem of the previous example. Once $y$ is assigned, the probability valuation of constraint $y > 2$ is $p_{y>2} = \sum_{i=3}^6 d_i^y = \frac{4}{6}$.

## 4.2 Filtering Algorithms for Relational Constraints on Random variables

It is stated in [28] that the major challenge towards solving chance constrained optimization problems lies in the computation of the probability and its derivatives of satisfying inequality constraints. In this section, probability filtering is provided for several constraints. Assumptions made by these filtering algorithms are that all random variables are independent and implemented using the distribution variables proposed in this paper. Given $p_C \equiv C(X, Y)$ the probability valuation of a constraint, the filtering algorithm should filter inconsistent values from $X$, from the distributions of $Y$ and $p_C$.

Let $C$ be a relational constraint defined on the random variables $Y$. We propose here a few filtering rules deriving from equations (4) and (8).

**Example 17 (Unary Equal).** Let $p_{y=c} \equiv (y = c)$ be the probability valuation of the equal constraint. In such settings, $\text{Pr}[p_{y=c} = 1] = d_c^y$. In addition for all the other values, another filtering can be defined by $\forall v \ne c, d_v^y \le 1 - \text{Pr}[p_{y=c} = 1]$. Note that this filtering should be directly done by the domain definition of the probability distribution.

**Example 18 (Unary Greater Than).** Let $p_{y>c} \equiv (y > c)$ be the probability valuation of the unary greater-than constraint. We propose to decompose equation (8) into two equations on the bounds of the random variable: $p_{y>c} \ge \sum_{v_i=c+1}^{\Omega(y)} d_{v_i}^y$; $p_{y>c} \le \sum_{v_i=c+1}^{\Omega(y)} d_{v_i}^y$. This rewriting is usual in CP, and linear propagation can be applied. Moreover, analogously to the reification, filtering can also be defined for values lower than, or equal to $c$. Indeed, as $p_{y>c}$ is the probability of satisfaction of the constraint, $p_{y\le c} = 1 - p_{y>c}$ is the probability of violating the constraint. This implies that the two filtering rules above should be used to propagate $d_v^y$, $\forall v \le c$.

**Example 19 (Binary Equal).** Let $p_{y_1=y_2} = \text{Pr}[y_1 = y_2]$ be the probability valuation of the equal constraint between two random variables. A filtering can be extracted from equation: $p_{y_1=y_2} = \sum_i d_i^{y_1} d_i^{y_2}$.

**Example 20 (Constraints $y \approx y_1 \text{ op } y_2$).** $y \approx y_1 + y_2$ is the definition of the distribution of variable $y$ from the distributions of variables $y_1$ and $y_2$. It is a probability-based constraint. In the general case, the probability distribution of the sum of two independent variables is given by the convolution: $\text{Pr}[y = k] = \sum_{i=-\infty}^{\infty} \text{Pr}[y_1 = k - i] \text{Pr}[y_2 = i]$. This convolution can be used to extract propagation rules for the ternary constraint. For each value $v_k$ of the possible values of $y$, a constraint of the form of $d_{v_k}^y = \sum_i d_{v_{k-i}}^{y_1} d_{v_i}^{y_2}$ is posted. Note that the same equation can be defined for the multiplication, subtraction, division, and modulo constraints.
**► Example 21 (Constraints $y = y_1 \text{ op } y_2$).** For example, $y = y_1 + y_2$ is the test of equality between the sampled values of $y$ and the sum of the sampled values of $y_1$ and $y_2$. This constraint can be expressed using the two constraints above. First, let $y_{1+2} \equiv y_1 + y_2$ be the distribution of the addition of $y_1$ and $y_2$. Second, let $p_{y=y_1+y_2}$ be the probability valuation of $y = y_{1+2}$. Finally, $p_{y=y_1+y_2} = \Pr[y = y_{1+2}]$.

## 4.3 Filtering Relational Constraints on Integer and Random Variables

Consider a relational constraint $C$ on decision variables $X$ and on random variables $Y$.

**► Proposition 22.** Set the variable $d_{v_i}^{y_i}$ from $y_i \in Y$ for constraint $C$ to 0 if:

$$ \forall t_X \in D^\times(X), \exists t_Y \in D_\Omega^\times(Y), t_{Y_i} = v_i \wedge \prod_{\forall v_j \in t_{Y_i} \setminus \{v_i\}} d_{v_j}^{y_j} > \mathbb{T}_C(t_X, t_Y) \quad (14) $$

This proposition is a direct rewriting of equation equations (4) and (8). In addition to the classical filtering of values in the decision variables if they do not have a valid tuple, another filtering should be defined.

**► Proposition 23.** Filter value $v_i$ from the domain of $x_i \in X$ for constraint $C$ if:

$$ \exists t_Y \in D_\Omega^\times(Y), \forall t_X = (\dots, v_i, \dots) \in D^\times(X), \prod_{\forall v_j \in t_Y} d_{v_j}^{y_j} > \mathbb{T}_C(t_X, t_Y) \quad (15) $$

Let $p_C(X,Y) \equiv C(X,Y)$ be the probability valuation of the constraint. A direct simple filtering can be extracted from equation (8). In addition, the filtering should be propagated back to the decision variables, for example:

**► Proposition 24.** Filter value $v_i$ from the domain of $x_i \in X$ for constraint $C$ if:

$$ \forall t_x = (\dots, v_i, \dots) \in D^\times(X), \Pr[C(t_x, Y)] \neq p_C \quad (16) $$

Here, equation (8) is applied to each tuple containing the value $v_i$. If none of them satisfies the probability valuation, then the value can be safely removed. Recall that $p_C$ represents $\Pr[p_C = 1]$. Moreover, in general, the $\neq$ will be checked for $>$ and $<$ to filter the bounds of the variables. The exact rule depends on the expression of the satisfaction of the constraint in terms of distribution variables, as shown in the following examples.

**► Example 25 (Lower Than).** Let $p_{x<y} \equiv (x < y)$ be the probability valuation of the lower-than constraint. This constraint is one of the most used in chance optimization as it can be used to model the deviation of scheduled amounts [8]. Analogously to the equal constraint, the pruning of variable $x$ can be defined by the following rule: $x < \max\{v_j \in x \mid \sum_{v_i=v_j} \overline{\Omega(y)} D(d_{v_i}^y) > D(p_{x<y})\}$ First, the opposite rule can also be defined by using $\overline{D(p_{x<y})}$. Then, note that this pruning rule, when both the $d_v^y$ and $p_{x<y}$ are fixed, is equivalent to the pruning of [30]. In addition, as the authors proposed, a closed formula may be used to directly extract the bounds of $x$. Now that $x$ is propagated, the probability distribution and valuation should be filtered. This can be done by using the following sum: $p_{x<y} = \sum_{v_i=x+1}^{\overline{\Omega(y)}} d_{v_i}^y$.

**► Example 26 (Equal).** Let $p_{x=y} \equiv (x = y)$ be the probability valuation of the equal constraint. First the pruning of variable $x$ can be defined using the following rule:

$$ \forall v \in \Omega(y), \overline{D(d_v^y)} < \overline{D(p_{x=y})} \lor \overline{D(d_v^y)} < \overline{D(p_{x=y})} \implies x \neq v \quad (17) $$
Then, the pruning of the probability valuation $p_{x=y}$ can be defined using:

$$p_{x=y} \geq \min(\{\overline{D(d_v^y)}|\forall v \in x\}); p_{x=y} \leq \max(\{\overline{D(d_v^y)}|\forall v \in x\}) \quad (18)$$

Finally, the pruning of the distribution of $y$ is an adaptation of the usual propagation:

$$x = v \implies d_v^y = p_{x=y} \quad (19)$$

## 4.4 Rewriting of constraints

In the previous sections, several constraints and propagation algorithms have been proposed. In this section, it is shown how to combine these constraints to model most use cases. Consider the not equal constraint. Even though no dedicated algorithm has been defined, propagation can be based on other ones. Let $p_{x \neq y} = (x \neq y)$ be the probability valuation of the not-equal constraint, $p_{x \neq y} = 1 - p_{x=y}$. This is analogous to the reification, indeed, $(x_1 = x_2) = \neg(x_1 \neq x_2)$. Several more constraints can be reformulated, for example the greater-less-than relationship: $p_{x \leq y} = p_{x-1 < y}$; $p_{x > y} = p_{-x < -y}$ (alternatively, $p_{x \geq y} = p_{-x-1 < -y}$). The in-not in relationship: $p_{y \notin X} = 1 - p_{y \in X}$. And the many logical relationships such as $p_{c_1 \land c_2} = p_{c_1} p_{c_2}$; $p_{c_1 \lor c_2} = 1 - (1 - p_{c_1})(1 - p_{c_2})$; and $p_{c_1 \oplus c_2} = p_{c_1} (1 - p_{c_2}) + (1 - p_{c_1}) p_{c_2}$. With $\oplus$ the xor logical operator. Note that some of these reformulations are direct encodings of the probability rules, such as for the and and or constraints.

# 5 Application: Virtual Network Design

In virtual network functions design, the main performance criterion is the total latency of the virtual network [12, 13, 16, 48]. In such a graph, each node is either a machine or a virtual machine, and has to process blocks of data. The processing time of a data block by a node is following a random distribution, which is conditioned by its configuration (both hardware and software). The global latency is the sum of the processing times of the nodes of the network. The goal of the virtual network design problem is to set the configurations for all the nodes such that the global latency is below $L$ milliseconds with probability $\gamma$, while the cost is minimized and other configuration constraints are satisfied.

<figure><img src="image_288.png" /><figcaption>Figure 1 Example of virtual network graph.</figcaption></figure>

Consider the network from Figure 1. This network contains 4 nodes ($c_1, ..., c_4$). In this graph, when a data block needs to be sent, it will go through $c_1$, then $c_2$, then $c_3$ and finally $c_4$. Note that in practice, the depth of the graph is often shallow. In our industrial application the number of nodes was 3. Moreover, each node contains several settings (cpu redundancy, queue maximum read times, queue batch pkt num, dpe, etc.) The latency of a node is influenced by all these settings. In addition, each node serves a particular role, and some operations must be done in at least one node. All of these requirements are contained in the constraints on the settings variables. The influence of the setting variables on the random variables' distribution is encoded as a table $T_s$ for each node.
## 5.1 Model

Let $N = (n_1, ..., n_k)$ be the nodes to configure. Let $y_i \sim D^{y_i}$ be the random variable representing the latency of node $n_i$. Let $S = (S^1, ..., S^k)$ with $S^i = s_1^i, ..., s_m^i$ being the settings variables associated to each node. Let $L$ be the maximum latency and $\gamma$ the minimal probability. Let $acc_i \in (acc_2, ..., acc_k)$ be the accumulated sum of the $i$ first random variables. Let $c_i$ be the cost associated with node $n_i$.

$$ \arg \min_{S \in \mathbb{Z}, Y \in \mathbb{P}} \sum_{i=1}^{k} c_i \qquad (20) $$

$$ \text{such that} \quad \text{Table}(S^i, c_i, y_i, T_s), \quad \forall i \in [1, k] \qquad (21) $$

$$ acc_1 \approx y_1 \qquad (22) $$

$$ acc_i \approx acc_{i-1} + y_i, \quad \forall i \in [2, k] \qquad (23) $$

$$ \Pr[acc_k < L] \geq \gamma \qquad (24) $$

$$ \text{ValidSettings}(S) \qquad (25) $$

## 5.2 Data

A large number of configurations have been heavily sampled to extract the distribution of probability. Once those were built, they have been used to extrapolate to unknown configurations. The latency is given in micro-seconds, with values in the set (50, 100, ..., 1450, 1500). Each node may be defined by around 100 configurations (valid assignment of the settings). The number of nodes varies between 2 to 5. Two algorithms are compared. First, *Variable* is the work proposed by this paper, where distribution variables are used. The second one is *Fixed*, it is an adaptation of the algorithm from [30, 36]. In *Fixed*, the distributions are considered unknown until all the setting variables used to define them are instantiated. It can be seen as a form of generate and test. The search strategy is the same for both algorithms. It starts by assigning the possible settings. Note that this is the only fair comparison with *Fixed* as the setting variables are the one influencing the possible distributions. In the general case, better strategies could be defined for *Variable*. Time out is set to 30 minutes. The integer CP solver used for all the experiments of this paper is our internal solver. The continuous part is solved by call to IBEX [4], analogously to CHOCO-solver [43].

## 5.3 Results

First, consider the results from Table 1. As we can see, finding a solution for the fixed model is not too hard, but as the size increases, it is not able to find the optimal one, or to prove it. In contrast, for the variable model, the optimal solution is always found and proved. The main reason is that invalid distributions are removed by the latency constraint before any decision leading to them is taken.

<table><caption>Table 1 For each instance set, the values are \#SAT (\#OPT). The last row is the mean time in seconds for proving optimality.</caption><thead><tr><th># Nodes</th><th>2</th><th>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>Fixed (baseline)</td><td>30 (30)</td><td>30 (22)</td><td>30 (0)</td><td>30 (0)</td></tr><tr><td>Variable</td><td>30 (30)</td><td>30 (30)</td><td>30 (30)</td><td>30 (30)</td></tr><tr><td>Variable avg time</td><td>0.57 s</td><td>10.2 s</td><td>125.4 s</td><td>752.3 s</td></tr></tbody></table>
Moreover, consider the plots from Figure 2. Instances are sorted by size. On the middle, the comparison on the objective value is made. As we can see, the objective value of the *fixed* method is optimal for small instances, even if no proof can be reached. Then, for larger instances, it is often quite far from the optimal (objective $\in [0, 3]$). In addition, on the left and right plots, it is seen that most of the time, the *fixed* method reaches the timeout. Note that the time increases exponentially with the number of nodes for both methods, but not the number of decisions for the *variable* method.

<figure><img src="image_289.png" /><figcaption>Figure 2 Instances are ordered w.r.t. size. (left) Number of decision (log base 10). (middle) Best objective value comparison. (right) Running time (log base 10).</figcaption></figure>

In conclusion, this experiment shows that the *variable* method proposed in this paper is better both in terms of quality of solution and in terms of time. In practice, solving the industrial instances of size 3 takes 17.75 seconds on average for the optimal solution.

## 6 Application: Chemical Deliveries

The following problem is part of a pipeline of chemicals processing and is composed of two assignment problems. First, chemicals are received by the factory every day. Once received, they must be stored into containers. These containers are restricted to some type of products. Each container already contains a known amount of products and has a maximum storage capacity. In addition, for some chemicals the total amount delivered is larger than the remaining storage capacity of the containers, given their currently stored quantities. In practice, this is not an issue as the chemicals are also processed, hence emptying the containers. The second problem concerns product assignment. Indeed, the chemicals are used to manufacture bio-sourced bases for perfumes. Those are later used by home-perfume makers to create reeds, candles, etc. A dozen of teams are spread over 4 buildings. Each building/team is specialized in a set of types of product, yet there are overlaps. Every day, in addition to the incoming deliveries, the factory must produce a given amount of several products. Selecting which team should work on which product is part of the problem to solve. The probability of emptying a container depends on the products associated to the teams in the building. A solution is an assignment of products to teams and deliveries to containers such that the stored quantity in each container minus the quantity that will be used does not exceed the maximum capacity with a high confidence.

More formally, this is an assignment problem. Each product $p \in P$ must be produced by a team $t \in T$. Each team $t \in T$ has a work capacity $W_t$, and work in a building $b_t \in B$. Each product $w$ requires a given amount of work $w_p$, and is compatible to a subset of the teams $T_p \in T$. Each delivery $j \in D$ stores a quantity $q_j$ in a container $c_i \in C$. Each container $c_i$ has a maximum capacity $C_i$. Each container $c_i$ will be emptied in parallel of a quantity $y_i$ unknown in advance. The exact quantity $y_i$ is unknown, but its distribution is influenced by the products processed in its building.