<table><caption>**Table 2** CP Solving Time Speedup compared to Omitting Redundant Bounds and Constraints.</caption><thead><tr><th rowspan="2">Redundant Bounds</th><th rowspan="2">Redundant Constraints</th><th colspan="5">Number of Trucks</th></tr><tr><th>50</th><th>100</th><th>150</th><th>200</th><th>250</th></tr></thead><tbody><tr><td>×</td><td>×</td><td>1.0x</td><td>1.0x</td><td>1.0x</td><td>1.0x</td><td>1.0x</td></tr><tr><td>✓</td><td>×</td><td>1.1x</td><td>1.4x</td><td>1.2x</td><td>0.9x</td><td>0.8x</td></tr><tr><td>×</td><td>✓</td><td>3.9x</td><td>6.8x</td><td>7.4x</td><td>3.9x</td><td>5.4x</td></tr><tr><td>✓</td><td>✓</td><td>3.5x</td><td>5.5x</td><td>6.3x</td><td>3.1x</td><td>5.5x</td></tr></tbody></table>

to find optimal schedules, which suggests that the CP computation time is more sensitive to changes in $\Delta$ than in $\sigma$. This is explained by the fact that increasing the flexibility significantly increases the search space. Again, the CP model shows substantial improvement over the initial solution with improvements ranging from 39% to 50% for the case study. When the flexibility increases, the necessary capacity goes down even before CP is applied, but the CP model is able to benefit more from the additional freedom.

## Impact of Redundant Bounds and Constraints

Table 2 shows the impact of the redundant bounds and constraints introduced in Section 2. These results are based on the baseline parameter settings for different numbers of trucks. It can be seen that only adding redundant bounds can both speed up or slow down the solver. For example, the solver becomes 1.4 times faster for $K = 100$ trucks, but a factor 0.8 slower for $K = 250$ trucks. The main benefit comes from adding redundant constraints, which speeds up the solver by up to 7.4 times for $K = 150$ trucks. When the redundant bounds and constraints are combined, performance may be improved further, as is the case for $K = 250$ trucks, but it appears that only including redundant constraints is the most efficient setting for these experiments.

# 5 Conclusion

The Autonomous Transfer Hub Network (ATHN) is one of the most promising ways to adapt self-driving trucks for the freight industry. This paper proposes a framework for optimizing ATHN operations with respect to hub utilization. To accomplish this, a MIP model generates autonomous vehicle routes and an initial schedule, which is then optimized using a CP model to reduce the required hub capacity for loading and unloading trucks. Results from the Ryder case study demonstrate the effectiveness of this approach, with the CP model reducing the total capacity by at least 42% while requiring only minor schedule modifications. This may save $15.2M per year in labor cost and is close to the lowest possible capacity for any initial schedule. Sensitivity analysis on loading/unloading duration and flexibility provides practical insights for ATHN operations and systematically shows the benefit of the CP model. Future work may explore alternative objective functions, such as minimizing parking space. Another interesting direction is to attempt to jointly optimize routes and hub capacity directly.

## References

1. Abderrahmane Aggoun and Nicolas Beldiceanu. Extending chip in order to solve complex scheduling and placement problems. *Mathematical and Computer Modelling*, 17(7):57–73, 1993. doi:10.1016/0895-7177(93)90068-A.
46:10 Constraint Programming to Improve Hub Utilization in ATHNs

2 Natashia Boland, Mike Hewitt, Luke Marshall, and Martin Savelsbergh. The Continuous-Time Service Network Design Problem. *Operations Research*, 65(5):1303-1321, 2017. doi: 10.1287/opre.2017.1624.

3 Natashia Boland and Martin Savelsbergh. Perspectives on integer programming for time-dependent models. *TOP*, 27:147-173, 2019. doi:10.1007/s11750-019-00514-4.

4 Kevin Dalmeijer and Pascal Van Hentenryck. Optimizing Freight Operations for Autonomous Transfer Hub Networks. arXiv:2110.12327, 2021. arXiv:2110.12327.

5 FleetOwner. TuSimple among autonomous truck companies to join Self-Driving Coalition. https://www.fleetowner.com/, 2021. URL: https://www.fleetowner.com/technology/autonomous-vehicles/article/21152006/tusimple-among-autonomous-truck-companies-to-join-selfdriving-coalition.

6 Forbes. Plus Partners With IVECO To Develop Automated Trucks For Global Deployment. https://www.forbes.com/, 2021. URL: https://www.forbes.com/sites/richardbishop1/2021/04/12/plus-partners-with-iveco-to-develop-automated-trucks-for-global-deployment.

7 FreightWaves. Gatik, Isuzu to partner on autonomous truck platform. https://www.freightwaves.com/, 2021. URL: https://www.freightwaves.com/news/gatik-isuzu-to-partner-on-autonomous-truck-platform.

8 Glassdoor. How much does a Truck Mechanic make? https://www.glassdoor.com/Salaries/truck-mechanic-salary-SRCH_K00.14.htm, 2023. Accessed: 2023-04-28.

9 Sönke Hartmann and Dirk Briskorn. An updated survey of variants and extensions of the resource-constrained project scheduling problem. *European Journal of Operational Research*, 297(1):1-14, 2022. doi:10.1016/j.ejor.2021.05.004.

10 Edward He, Natashia Boland, George Nemhauser, and Martin Savelsbergh. An exact algorithm for the service network design problem with hub capacity constraints. *Networks*, 80(4):572-596, 2022. doi:10.1002/net.22128.

11 Heavy Duty Trucking. Daimler's Redundant Chassis for Autonomous-Truck Operation. https://www.truckinginfo.com/, 2021. URL: https://www.truckinginfo.com/10157338/daimlers-redundant-chassis-for-autonomous-truck-operation.

12 Willy Herroelen, Bert De Reyck, and Erik Demeulemeester. Resource-constrained project scheduling: A survey of recent developments. *Computers & Operations Research*, 25(4):279-302, 1998. doi:10.1016/S0305-0548(97)00055-5.

13 Chungjae Lee, Kevin Dalmeijer, and Pascal Van Hentenryck. Optimization Models for Autonomous Transfer Hub Networks. arXiv:2201.06137, 2022. arXiv:2201.06137.

14 Chungjae Lee, Kevin Dalmeijer, Pascal Van Hentenryck, and Peibo Zhang. Optimizing Autonomous Transfer Hub Networks: Quantifying the Potential Impact of Self-Driving Trucks. arXiv:2305.03119, 2023. arXiv:2305.03119.

15 Luke Marshall, Natashia Boland, Martin Savelsbergh, and Mike Hewitt. Interval-Based Dynamic Discretization Discovery for Solving the Continuous-Time Service Network Design Problem. *Transportation Science*, 55(1):29-51, 2021. doi:10.1287/trsc.2020.0994.

16 Clair E. Miller, Albert W. Tucker, and Richard A. Zemlin. Integer Programming Formulation of Traveling Salesman Problems. *Journal of the ACM*, 7(4):326-329, 1960. doi:10.1145/321043.321046.

17 OpenStreetMap. Planet dump retrieved from https://planet.osm.org, 2021. URL: https://www.openstreetmap.org.

18 PACE. Partnership for an Advanced Computing Environment (PACE), 2017. URL: http://www.pace.gatech.edu.

19 Laurent Perron and Vincent Furnon. Google OR-Tools v9.6. https://developers.google.com/optimization/, 2023.

20 A. Alan B. Pritsker, Lawrence J. Waiters, and Wolfe Philip M. Multiproject Scheduling with Limited Resources: A Zero-One Programming Approach. *Management Science*, 60(1):93-108, 1969. doi:10.1287/mnsc.16.1.93.
21 Roland Berger. Shifting up a gear – Automation, electrification and digitalization in the trucking industry. https://www.rolandberger.com/, 2018. URL: https://www.rolandberger.com/publications/publication pdf/roland_berger_trucking_industry.pdf.

22 Ryder System, Inc. and Socially Aware Mobility Lab. The Impact of Autonomous Trucking: A Case-Study of Ryder's Dedicated Transportation Network. *Ryder Newsroom*, 2021. URL: https://newsroom.ryder.com/news/news-details/2021/Ryder-Teams-Up-with-Georgia-Tech-for-Industrys-First-Data-Driven-Study-on-Impact-of-Autonomous-Trucking/.

23 Yannick Oskar Scherr, Mike Hewitt, Bruno Albert Neumann-Saavedra, and Dirk Christian Mattfeld. Dynamic discretization discovery for the service network design problem with mixed autonomous fleets. *Transportation Research Part B: Methodological*, 141:164–195, 2020. doi:10.1016/j.trb.2020.09.009.

24 Andreas Schutt, Thibaut Feydy, Peter J. Stuckey, and Mark G. Wallace. Explaining the cumulative propagator. *Constraints*, 16:250–282, 2011. doi:10.1007/s10601-010-9103-2.

25 Mohsen Shahandasht, Binaya Pudasaini, and Sean Logan McCauley. Autonomous Vehicles and Freight Transportation Analysis. Technical report, The University of Texas at Arlington, 2019. doi:10.13140/RG.2.2.28484.78726.

26 TechCrunch. Aurora and Volvo partner to bring autonomous long-haul trucks to North America. https://www.techcrunch.com/, 2021. URL: https://techcrunch.com/2021/03/30/aurora-and-volvo-partner-to-bring-autonomous-long-haul-trucks-to-north-america/.

27 Transport Topics. Navistar, TuSimple Partner to Launch Self-Driving Trucks in 2024. https://www.ttnews.com/, 2020. URL: https://www.ttnews.com/articles/navistar-tusimple-partner-launch-self-driving-trucks-2024.

28 Steve Viscelli. Driverless? Autonomous Trucks and the Future of the American Trucker. Center for Labor Research and Education, University of California, Berkeley, and Working Partnerships USA, 2018. URL: http://driverlessreport.org/.

29 Duc Minh Vu, Mike Hewitt, Natashia Boland, and Martin Savelsbergh. Dynamic Discretization Discovery for Solving the Time-Dependent Traveling Salesman Problem with Time Windows. *Transportation Science*, 54(3):703–720, 2020. doi:10.1287/trsc.2019.0911.

30 Haotian Wu, Ian Herszterg, Martin Savelsbergh, and Yixiao Huang. Service Network Design for Same-Day Delivery with Hub Capacity Constraints. *Transportation Science*, 57(1):273–287, 2023. doi:10.1287/trsc.2022.1155.
<figure><img src="image_459.png" /></figure>
# A New Approach to Finding 2 x n Partially Spatially Balanced Latin Rectangles

Renee Mirka ☐

Cornell University, Ithaca, NY, USA

Laura Greenstreet ☐

Cornell University, Ithaca, NY, USA

Marc Grimson ☐

Cornell University, Ithaca, NY, USA

Carla P. Gomes ☐

Cornell University, Ithaca, NY, USA

## Abstract

Partially spatially balanced Latin rectangles are combinatorial structures that are important for experimental design. However, it is computationally challenging to find even small optimally balanced rectangles, where previous work has not been able to prove optimality for any rectangle with a dimension above size 11. Here we introduce a graph-based encoding for the 2 × *n* case based on finding the minimum-cost clique of size *n*. This encoding inspires a new mixed-integer programming (MIP) formulation, which finds exact solutions for the 2 × 12 and 2 × 13 cases and provides improved bounds up to *n* = 20. Compared to three other methods, the new formulation establishes the best lower bound in all cases and establishes the best upper bound in five out of seven cases.

**2012 ACM Subject Classification** Applied computing → Operations research

**Keywords and phrases** Spatially balanced Latin squares, partially spatially balanced Latin rectangles, minimum edge weight clique, combinatorial optimization, mixed integer programming, imbalance, cliques

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

**Category** Short Paper

**Funding** This material is based upon work partially supported by the National Science Foundation under Awards CCF-1522054 and CCF-2007009 and by AFOSR DURIP grant FA9550-21-1-0316, AFOSR MURI grant FA9550-18-1-0136, and AFOSR grant FA9550-20-1-0421.

## 1 Introduction

Latin squares and rectangles are combinatorial objects represented by *n* × *n* or *k* × *n* grids with entries assigned from {1, 2, ..., *n*} such that no entry is repeated in a row or column. They are important structures for experimental design. For example, in agronomic field experiments, experts are interested in applying *n* treatments consisting of *n* fertilizers in different orderings, which can be achieved by designing the treatment sequences following a Latin square. Arbitrary Latin squares are not hard to generate. However, geometric imbalance due to some treatments occurring closer together more frequently can bias experimental results [21]. This motivates the use of *spatially balanced* Latin squares and rectangles which require additional structure to capture the notion of distance between any two treatments in the square or rectangle and are more computationally challenging to construct.

A large body of previous work has focused on spatially balanced Latin squares, including introducing streamlining constraints [8], using stochastic optimization [11, 10, 12, 6], and applying local search methods [22]. While spatially balanced Latin squares have been extensively studied [9, 16, 17, 21], there has been much less work on spatially balanced Latin

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

© Renee Mirka, Laura Greenstreet, Marc Grimson, and Carla P. Gomes; 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. 47; pp. 47:1–47:11

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

LIPICS

Leibniz International Proceedings in Informatics

Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl Publishing, Germany
47:2 A New Approach to Finding 2 x n Partially Spatially Balanced Latin Rectangles

rectangles [4], despite rectangular conditions occurring more frequently in practice. For
example, in the case of fertilizer treatments, the number of treatments is often less than the
number of fertilizers resulting in a rectangular structure.

The main previous work investigating spatially balanced Latin rectangles established the non-existence of perfectly balanced Latin rectangles for an infinite family of sizes and shifted focus to constructing *partially* spatially balanced Latin rectangles (PSBLRs) [4]. Diaz et al. introduce and experimentally compare three approaches to generating PSBLRs: an assignment-based mixed-integer program (MIP), a constraint satisfaction program (CP), and a random-restart hill climbing local search. Using these approaches, Diaz et al. were able to find the provably optimal imbalance for rectangles up to a size of 2 × 11 and provide bounds up to 12 × 12.

In this work, we focus on the 2 × n case of constructing PSBLRs, introducing a new graph-based encoding based on a reduction to a minimum-edge weight clique problem. The maximum/minimum edge-weight clique problem (MEWC) is a generalization of the classic max clique problem, where given an input graph G = (V, E) instead of finding the largest complete subgraph the goal is to find the subgraph with the largest/smallest sum of weighted edges. The MEWC problem has many applications including in experimental design [3], molecular biology [20], and materials discovery [1]. Multiple approaches have been developed for MEWC problems including linear and quadratic mixed-integer programs [5, 7, 18], branch-and-cut [5, 15, 19], and heuristic methods [2, 14]. Several benchmarks have been developed for MEWC, though less than half the instances have been solved to optimality [13].

For our experiments, we use a straight-forward MIP formulation of MEWC to emphasize the benefits of the reduction instead of advanced techniques developed for MEWC. This new MIP formulation for the problem of finding 2 × n PSBLRs finds optimal solutions for the previously unsolved cases of n = 12 and 13 and provides new bounds for n = 14 − 20. Further, we demonstrate the new clique-based MIP formulation outperforms the previous assignment-based MIP formulation and A* search and finds comparable solutions to local search while providing lower bounds.

2 Preliminaries

► **Definition 1 (Latin Rectangle).** Let $k$ and $n$ be positive integers with $k \le n$ and $R$ be a $k \times n$ matrix. Then $R$ is a Latin rectangle if every entry of $R$ contains a number in $[n] = \{1, 2, \dots, n\}$ and no number is repeated in any row or column.

In order to discuss partially spatially balanced Latin rectangles, we first introduce the notion of distance between two symbols in a Latin rectangle and imbalance:

► **Definition 2 (Imbalance).** For two symbols $u, v \in [n]$ and $i \le k$, the distance between $u$ and $v$ in row $i$, denoted $d_i(u, v)$, is the absolute value of the difference of the indices of the positions of $u$ and $v$ in row $i$. The overall distance between $u$ and $v$ is then $d(u, v) = \sum_{i \le k} d_i(u, v)$.

The spatial imbalance of a Latin rectangle, $R$, is defined by

$$
\mathbb{I}(R) = \sum_{i,j} \left| d(i,j) - \frac{k(n+1)}{3} \right|.
$$

**Definition 3 (Spatially Balanced Latin Rectangle).** We say that a $k \times n$ rectangle $R$ is spatially balanced if all distances $d(u, v)$ are the same. For brevity, we denote this as $SBLR(k, n)$.

Note that in this case, Diaz et al. show that each distance must be exactly $\frac{k(n+1)}{3}$ and thus the imbalance of the rectangle is 0; see [4] for a short proof of this proposition.
**► Proposition 4.** If there exists a solution for SBLR(k, n) then the distance between any pair of symbols is equal to $k(n + 1)/3$.

As a corollary, $k \equiv_3 0$ or $n \equiv_2 2$ are necessary conditions for the existence of SBLR(k, n) as $k(n+1) \equiv_3 0$ must hold. Diaz et al. have further shown the non-existence of SBLRs of size $2 \times n$ for $n \neq 2$ and $3 \times n$ for $n \neq 3$, as well as experimentally demonstrated that no perfectly spatially balanced rectangles with $k \neq n$ exist up to size $7 \times 7$ [4]. In practice, it is useful to minimize the imbalance even if it is not possible to reduce it to zero. This motivates the definition of the object of our study- partially spatially balanced Latin rectangles (PSBLRs):

**► Definition 5 (Partially Spatially Balanced Latin Rectangle).** A $k \times n$ Latin rectangle $R$ is partially spatially balanced if $\mathbb{I}(R)$ is minimized. In particular, $R$ is partially spatially balanced if for any $k \times n$ Latin rectangle $R'$, $\mathbb{I}(R) \le \mathbb{I}(R')$.

### 3 Graph Encoding

In this section, we describe a new graph-based encoding for the problem of constructing PSBLRs. For the $2 \times n$ Latin rectangles we consider, we will always assume the first row is ordered and given by $1 \ 2 \ \dots \ n$, as any other solution can be transformed into a solution in this form through relabelling symbols. This reduces the problem to selecting a single imbalance-minimizing derangement, or permutation with no fixed points, for the second row of the rectangle. It also simplifies the computation of the distance between any two symbols $u$ and $v$. Particularly, for $1 \le u < v \le n$, we have $d(u, v) = v - u + |i_v - i_u|$ where $i_v, i_u$ are the indices of $v, u$ in the second row, respectively. For fixed $u, v$ and $n$, $d(u, v)$ is solely determined by the choices for $i_u$ and $i_v$. We exploit this through the construction of a graph which encodes how these choices affect the imbalance.

The graph encoding is as follows. For a given $n$, construct a complete graph $G = (V, E)$ with $V = \{v_{ij} : 1 \le i, j \le n, i \ne j\}$ and $E = \{(u, v) : u, v \in V\}$. Intuitively, a vertex $v_{ij}$ represents a Latin rectangle where $j$ is the index of $i$ in the second row. We include a cost on each edge ($v_{st}, v_{qr}$) which represents how much is contributed to the imbalance from $d(s, q)$ if $t$ and $r$ are the indices of $s$ and $q$ in the second row, respectively. In particular, for $e = (v_{st}, v_{qr}) \in E$ with $s \neq q$ and $t \neq r$, let

$$ 
\begin{aligned}
c_e &= \left| d(s, q) - \frac{2(n+1)}{3} \right| \\
&= \left| s - q + |i_s - i_q| - \frac{2(n+1)}{3} \right| \\
&= \left| s - q + |t - r| - \frac{2(n+1)}{3} \right|.
\end{aligned}
$$

If $s = q$ or $t = r$, the cost on the edge is $\infty$, since the assignments these vertices represent are not valid for a Latin rectangle. See Figure 1 for an example when $n = 3$.

To complete the problem encoding, we need to discern how to recover a PSBLR from this complete graph. The key observation is that every $n$-clique where all edges have finite cost corresponds to a valid Latin rectangle. The finiteness of the costs ensures that the assignment rules of the Latin rectangle are obeyed. Furthermore, the sum of the costs on the edges of an $n$-clique is exactly equal to the imbalance of its corresponding Latin rectangle. As such, to find a PSBLR, it suffices to find an $n$-clique whose sum of edge costs is minimal. Figure 1 illustrates an optimal clique corresponding to an optimal Latin rectangle for $n = 3$:

$$ \begin{bmatrix} 1 & 2 & 3 \\ 2 & 3 & 1 \end{bmatrix} $$
<figure><img src="image_462.png" /><figcaption>■ **Figure 1** (Left) The graph encoding after including all edges with infinite costs between vertices $v_{sq}$ and $v_{sr}$ for $s = 1, 2, 3$ (thin edges) and between vertices $v_{sq}$ and $v_{rq}$ for $q = 1, 2, 3$ (bold edges). (Center) The graph encoding after including all remaining edges with finite costs (bold edges) between the remaining pairs of vertices. The boxed value is calculated as $|(2 - 1) + (3 - 1) - 8/3| = 1/3$. (Right) A 3-clique (bold edges) in the graph encoding when $n = 3$ corresponding to a $2 \times 3$ PSBLR.</figcaption></figure>

4 MIP Formulation

Now that we have reduced the problem of finding a 2 × n PSBLR to finding a minimum-cost
n-clique, we can model the problem through a new integer program:

$$
\begin{align}
\text{minimize} \quad & \sum_{e \in E} c_e x_e \nonumber \\
\text{subject to} \quad & \sum_{v \in V} z_v = n, \tag{1}
\end{align}
$$

$$
\begin{align*}
x_e &\ge z_v + z_u - 1, && \forall e = (u, v) \in E, \tag{2} \\
x_e &\in \{0, 1\}, && \forall e \in E, \\
z_v &\in \{0, 1\} && \forall v \in V.
\end{align*}
$$

In this model, we have binary variables $z_v$ and $x_e$ for each vertex $v$ and edge $e$ in the graph
representing whether the vertex or edge is included in the *n*-clique. While (1) guarantees that
we select *n* vertices, (2) requires that any edges between two selected vertices are included in
the cost; together these ensure the model selects an *n*-clique.

Note that the size of this model is polynomial in $n$ with $n(n-1)+n(n-1)(n(n-1)-1)/2$ variables and $n(n-1)(n(n-1)-1)/2 + 1$ constraints. Furthermore, by construction, an optimal solution will be a PSBLR. However, the presence of so many binary variables makes finding the optimal solution a cumbersome computational task.

The following observation allows us to partially relax the integer program and significantly
reduce the number of binary variables:

**► Observation 6.** For each edge $e \in E$, $x_e$ will be 0 or 1 if $z_v$ is binary for all $v \in V$. In other words, we can relax the restriction that $x_e$ is binary for all $e \in E$ by replacing it with $x_e \ge 0$ for all $e \in E$ and still guarantee a binary solution.

The relaxed mixed integer program only requires $n(n-1)$ binary variables, one for each vertex, where the remaining edge variables are continuous.
# 5 Methods

We compared the min-cost n-clique MIP formulation, hereafter referred to as the clique formulation, to an assignment-based MIP formulation, A\* search, and a local search approach.

In the assignment based MIP formulation, we again assume the first row is fixed to 1 2 ... n. Thus, we only need to introduce binary variables indicating whether value i occurs in position j of the second row, i, j ∈ {1, 2, ..., n} ≡ [n], i ≠ j. Note that these exactly correspond to the vertex variables in the clique formulation, z_v = z_{v_{ij}} for v ∈ {v_{ij} : i, j ∈ [n], i ≠ j}. As the imbalance is dependent on pairwise positions in the second row, we also introduce binary variables that are 1 if a pair of variables is included in the solution. Note that these correspond exactly to the edge variables in the clique formulation, x_e for e ∈ {(u, v) : u, v ∈ V}, where the pairwise cost corresponds to c_e. Finally, we introduce constraints ensuring each value appears exactly once in the second row (3) and each position in the second row is assigned exactly one value (4):

$$ \text{minimize} \sum_{e \in E} c_e x_e $$

$$ \text{subject to} \quad \sum_{j \in [n], j \neq i} z_{v_{ij}} = 1, \quad \forall i \in [n] \quad (3) $$

$$ \sum_{i \in [n], i \neq j} z_{v_{ij}} = 1, \quad \forall j \in [n] \quad (4) $$

$$ x_e \geq z_v + z_u - 1, \quad \forall e = (u, v) \in E, $$

$$ x_e \in \{0, 1\}, \quad \forall e \in E, $$

$$ z_v \in \{0, 1\} \quad \forall v \in V. $$

Similar to the clique formulation, we can partially relax the integer program by allowing the $x_e$ variables to assume continuous values, reducing the number of binary variables to $n(n-1)$.

Both MIP formulations were implemented in Gurobi and CPLEX using their respective Python APIs. As previous work has established optimal solutions up to $n = 11$, we replicated previous results and further tested instances with $n = 12 - 20$. For each $n$, both models were tested for 6 hours running on a Intel Xeon 6154 processor and allowed 32GB of memory in three configurations: a single instance running on a single thread, a single instance running on 32 threads, and 8 instances running on a total of 32 threads.

For A\* search, we again fix the first row and only consider positions in the second row. Starting with an empty second row, we add all valid placements of the first digit to a priority queue, assigning each node a cost based on the partial imbalance of the filled digits and an admissible heuristic for the imbalance due to the unplaced digits. We repeatedly evaluate the first node in the queue until we reach a node with all digits filled which is guaranteed to be an optimal solution. As a heuristic, for each unplaced digit we calculated the minimum increase to the imbalance that would result from placing it in one of the remaining positions, accounting only for interactions with already placed digits. We then summed these values across all unplaced digits.

A\* search was implemented in C++ and was also tested for 6 hours each on instances of size $n = 12 - 20$ on an Intel Xeon 6154 processor. As A\* search is memory intensive, runs were allocated 128GB of memory, compared to 32GB for the other methods. If A\* search does not complete, the value of the node at the front of the queue can be used as a lower bound on the solution value. However if the algorithm does not run to completion, no feasible solutions are found, resulting in no upper bound on the imbalance.
47:6 A New Approach to Finding 2 x n Partially Spatially Balanced Latin Rectangles

We implemented local search using a random walk from a random initial derangement, where at each step the position of two random digits were swapped, maintaining feasibility. The random walk was implemented in C++ and tested for 6 hours each on instances of size $n = 12 - 20$ on an Intel Xeon 6154 processor and allocated 32GB of memory. Unlike the other methods, a random walk cannot prove optimality or provide a lower bound on the solution.

In addition to testing methods to bound the optimal solution, we performed several tests to better characterize the solution space. First, for $n = 2 - 10$ we brute force computed all solutions to determine the number of optimal solutions. We were inspired to characterize the number of solutions due to the observation that mirroring either the first or second row results in a rectangle with the same imbalance. For example if we have a solution $R$, we can construct a rectangle $R'$ with the second row mirrored by setting $R'_{2i} = n + 1 - R_{2i}$ for $i \in [n]$. This results in the same imbalance:

$$ 
\begin{aligned}
\mathbb{I}(R') &= \sum_{i,j} \left| |R_{1i} - R_{1j}| - |(n+1 - R_{2i}) - (n+1 - R_{2j})| - \frac{2(n+1)}{3} \right| \\
&= \sum_{i,j} \left| |R_{1i} - R_{1j}| - |R_{2i} - R_{2j}| - \frac{2(n+1)}{3} \right| \\
&= \mathbb{I}(R)
\end{aligned}
$$

However, it is not guaranteed that $R'$ is a Latin rectangle, as the rearrangement may not respect the column constraint. Similarly, we can create a rectangle, $R^\dagger$, with the first row mirrored by setting $R^\dagger_{1i} = n + 1 - R_{1i}$ for $i \in [n]$. While this does not result in a rectangle with the first row ordered as 1 2 ... n, we can relabel the variables so the first row goes from 1 2 ... n which results in the second row having the labels $R^\dagger_{2i} = R_{2(n+1-i)}$. While this could potentially allow up to four solutions from a single solution, there is no guarantee that these solutions satisfy the column constraints. As multiple solutions are not guaranteed, we cannot use methods like streamlining. However, the presence of multiple solutions may make it difficult for the MIP formulations to prove optimality.

Further, we sought to characterize the distribution of feasible solutions. For $n = 2 - 10$ we were able to compute the distribution exactly. For $n = 11 - 16$, we randomly sampled a million solutions to approximate the distribution. The distribution of solutions impacts all four methods, where having many solutions with nearly optimal imbalance will benefit a random walk, while it will make it difficult for the MIP formulation and A* search to prove optimality.

6 Results

While CPLEX had better runtimes up to $n = 10$, only Gurobi was able to prove optimality for $n = 12$ and 13 and Gurobi consistently found better bounds for $n = 14 - 20$ (see Tables 3-5). Thus, for the remainder of the results we report the performance of the Guorbi model for both MIP formulations. In general, running the model using a single instance and a single thread performed better than running a single instance with multiple threads or multiple instances, with the only notable exception being for $n = 20$, where running with multiple instances found the smallest upper bound of 816 (See Tables 6-7).

While both MIP formulations were able to prove optimality for $n = 12$, only the clique formulation was able to prove optimality for $n = 13$ within the 6 hours allotted (See Table 8). While local search is not able to prove optimality, it found the optimal solution for both
$n = 12$ and $13$. For all $n = 2-20$, the clique formulation resulted in the largest lower bound (See Table 1). The results were more varied for the upper bound, where the random walk found the lowest upper bound for $n = 14$ and $17$, the clique found the lowest upper bound for $n = 18, 19$, and $20$, and both methods found the same upper bound for $n = 15$ and $16$ (See Table 2). For $n = 16$, the assignment formulation matched the upper bound found by the clique formulation and random walk.

A* was only able to run to completion for $n = 12$. For $n = 13-15$, A* ran out memory, despite being allotted 128GB where all other methods used under 32GB. For $n > 15$, A* was only able to reach configurations with at most seven fixed values, resulting in a weak lower bound on the optimal solution.

<table><caption>Table 1 Lower bound by method for $n = 12-20$. Bold values indicate the best bound, i.e. the largest lower bound, for each $n$. All methods were allotted 6 hours.</caption><thead><tr><th></th><th colspan="10">n</th></tr><tr><th></th><th>12</th><th>13</th><th>14</th><th>15</th><th>16</th><th>17</th><th>18</th><th>19</th><th>20</th></tr></thead><tbody><tr><th>Clique MIP</th><td><strong>168</strong></td><td><strong>218.66</strong></td><td><strong>191</strong></td><td><strong>252.66</strong></td><td><strong>273.66</strong></td><td><strong>252</strong></td><td><strong>402.66</strong></td><td><strong>465.33</strong></td><td><strong>155</strong></td></tr><tr><th>Assignment MIP</th><td><strong>168</strong></td><td><strong>218.66</strong></td><td>131</td><td>141.33</td><td>133</td><td>91</td><td>121</td><td>138.66</td><td>47</td></tr><tr><th>A* search</th><td><strong>168</strong></td><td>191</td><td>166</td><td>169</td><td>164.66</td><td>134</td><td>159.33</td><td>158.66</td><td>106</td></tr><tr><th>Random walk</th><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

<table><caption>Table 2 Upper bound by method for $n = 12-20$. Bold values indicate the best bound, i.e. the smallest upper bound, for each $n$. All methods were allotted of 6 hours.</caption><thead><tr><th></th><th colspan="10">n</th></tr><tr><th></th><th>12</th><th>13</th><th>14</th><th>15</th><th>16</th><th>17</th><th>18</th><th>19</th><th>20</th></tr></thead><tbody><tr><th>Clique MIP</th><td><strong>168</strong></td><td><strong>218.66</strong></td><td>272</td><td><strong>345.33</strong></td><td><strong>427.33</strong></td><td>522</td><td><strong>617.33</strong></td><td><strong>732</strong></td><td><strong>816</strong></td></tr><tr><th>Assignment MIP</th><td><strong>168</strong></td><td><strong>218.66</strong></td><td>272</td><td>346</td><td><strong>427.33</strong></td><td>526</td><td>621.33</td><td><strong>732</strong></td><td>886</td></tr><tr><th>A* search</th><td><strong>168</strong></td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><th>Random walk</th><td><strong>168</strong></td><td><strong>218.66</strong></td><td><strong>268</strong></td><td><strong>345.33</strong></td><td><strong>427.33</strong></td><td><strong>508</strong></td><td>623.33</td><td>742.66</td><td>866</td></tr></tbody></table>

For the exploratory tests, for all $n = 2-10$ there are multiple optimal solutions (See Figure 2). The number of solutions did not show any obvious patterns, such as a monotonic increase or alternation between odd and even digits. For $n = 6-10$ where we were able to compute the distribution of solutions exactly, up to 10% of solutions were within 10% of optimal and as many as 40% of solutions were within 20% of optimal (See Figure 2). For $n = 11-20$ we looked at the distribution over a million random solutions instead of the full solution space, so the best solution likely does not represent the optimal solution. However, the trend of a large proportion of the cumulative distribution being within 20% of the best solution continued, where at the high-end of the range over 40% of solutions were within 10% of the best solution and over 90% were within 20% of the best solution.

# 7 Discussion

Both the clique and assignment formulation were able to prove optimality for the previously unsolved case of $n = 12$ and the clique formulation was further able to prove optimality for $n = 13$. For $n = 14-20$, the clique formulation consistently found the best lower-bounds, outperforming the assignment formulation and A* in all seven cases. While the clique formulation and a random walk were comparable in establishing upper bounds, finding the
<figure><img src="image_463.png" /><figcaption>■ **Figure 2** (Left) Number of optimal solutions by *n*. (Right) Cumulative distribution of solutions as proportion of the best solution. For *n* = 6–10 (green), distributions represent all solutions. For *n* = 11–20 (blue), distributions are over a million random solutions. Dashed lines indicate solutions that are 10% and 20% greater than the best solution respectively.</figcaption></figure>

best upper bound in five and four cases respectively, the clique formulation appeared to perform better for large *n*, establishing the best upper bound for *n* = 18 – 20. Thus, the clique formulation outperformed each of the other three methods individually, as well as performed better than using a combination of other methods, such as using A* to establish lower bounds and a random walk to establish upper bounds.

While the clique formulation outperformed the other methods, no method was able to prove the optimality of a solution above *n* = 13. The exploratory results shed some light on why it is difficult to find even relatively small PSBLRs. First, the potential for mirror solutions as well as the presence of multiple solutions for all *n* = 3 – 10 make it challenging for the MIP formulations to prove optimality, as potentially many branches of the branch-and-bound tree include optimal solutions and must be extensively explored before they can be pruned. Further, while there are likely multiple solutions, the possibility of a mirror rectangle violating a column constraint means we cannot guarantee multiple solutions and use techniques like streamlining constraints to guide the model towards a single optimal solution.

The large portion of the cumulative distribution close to the optimal solution further hampers the MIP formulations, making it difficult to prune branches that include near-optimal solutions. This large number of near-optimal solutions also hampers A* search, where a very tight heuristic is needed to allow for any significant amount of pruning, where experimentally we found that most nodes could generally only be pruned once all but one or two values had been assigned. While a random walk would benefit from multiple optimal solutions, a random walk or other local-search based methods cannot establish optimality.

8 Conclusion

Spatially and partially spatially balanced Latin squares and rectangles are important com-
binatorial structures used for experimental design. While spatially balanced Latin squares
have been extensively studied, relatively little work has considered spatially balanced Latin
rectangles, despite them occurring more frequently in practice.
We introduce a new graph-based encoding for a 2 × *n* Latin rectangle which inspires a new MIP formulation based on the MEWC problem. This new formulation found optimal solutions previously unsolved 2 × 12 and 2 × 13 cases and outperforms an assignment-based MIP formulation, A\* search, and a random-walk based local-search, establishing improved bounds up to *n* = 20. Given the success of our straight-forward MIP implementation for MEWC, a direction for future work is to explore whether more advanced MEWC methods are able to find optimal PSBLRs for *n* > 13.

Further, our exploratory results help characterize what make finding PSBLRs computationally challenging. The potential, but not guarantee, of multiple optimal solutions makes it difficult for mathematical programming methods to establish optimality and prevents the use of standard methods for handling multiple solutions, such as streamlining. The large number of near-optimal solutions makes it difficult for informed search-based methods to prune any significant amount of solutions, requiring near brute-force exploration of the search space. While simple to encode, this problem is quite challenging, making it an ideal benchmark for future search and optimization methods and we encourage further exploration of the problem.

References

1. Luis A Agapito, Marco Fornari, Davide Ceresoli, Andrea Ferretti, Stefano Curtarolo, and Marco Buongiorno Nardelli. Accurate tight-binding hamiltonians for two-dimensional and layered materials. *Physical Review B*, 93(12):125137, 2016.

2. Bahram Alidaee, Fred Glover, Gary Kochenberger, and Haibo Wang. Solving the maximum edge weight clique problem via unconstrained quadratic programming. *European Journal of Operational Research*, 181(2):592–597, 2007.

3. Egon Balas and Chang Sung Yu. Finding a maximum clique in an arbitrary graph. *SIAM Journal on Computing*, 15(4):1054–1068, 1986.

4. Mateo Díaz, Ronan Le Bras, and Carla Gomes. In search of balance: The challenge of generating balanced latin rectangles. In *Proceedings of Integration of AI and OR Techniques in Constraint Programming: 14th International Conference*, pages 68–76. Springer, 2017.

5. G Dijkhuizen and Ulrich Faigle. A cutting-plane approach to the edge-weighted maximal clique problem. *European Journal of Operational Research*, 69(1):121–130, 1993.

6. Stefano Ermon, Carla Gomes, Ashish Sabharwal, and Bart Selman. Low-density parity constraints for hashing-based discrete integration. In *Proceedings of International Conference on Machine Learning*, pages 271–279. PMLR, 2014.

7. U Faigle, R Garbe, K Heerink, and B Spieker. Lp—relaxations for the edge—weighted subclique problem. In *Operations Research'93: Extended Abstracts of the 18th Symposium on Operations Research held at the University of Cologne September 1–3, 1993*, pages 157–160. Springer, 1994.

8. Carla Gomes and Meinolf Sellmann. Streamlined constraint reasoning. In *Proceedings of Principles and Practice of Constraint Programming: 10th International Conference*, pages 274–289. Springer, 2004.

9. Carla Gomes, Meinolf Sellmann, Cindy Van Es, and Harold Van Es. The challenge of generating spatially balanced scientific experiment designs. In *Proceedings of the International Conference on Integration of AI and OR Techniques in Constraint Programming for Combinatorial Optimization Problems*, pages 387–394. Springer, 2004.

10. Carla P Gomes, Joerg Hoffmann, Ashish Sabharwal, and Bart Selman. Short xors for model counting: from theory to practice. In *Proceedings of Theory and Applications of Satisfiability Testing-SAT 2007: 10th International Conference*, pages 100–106. Springer, 2007.
47:10 A New Approach to Finding 2 x n Partially Spatially Balanced Latin Rectangles

11 Carla P Gomes, Ashish Sabharwal, and Bart Selman. Model counting: A new strategy for obtaining good bounds. In *Proceedings of the AAAI Conference on Artificial Intelligence*, volume 10, pages 1597538–1597548, 2006.

12 Carla P Gomes, Willem Jan van Hoeve, Ashish Sabharwal, and Bart Selman. Counting csp solutions using generalized xor constraints. In *Proceedings of the AAAI Conference on Artificial Intelligence*, pages 204–209, 2007.

13 Seyedmohammadhossein Hosseinian, Dalila BMM Fontes, Sergiy Butenko, Marco Buongiorno Nardelli, Marco Fornari, and Stefano Curtarolo. The maximum edge weight clique problem: formulations and solution approaches. *Optimization Methods and Applications: In Honor of Ivan V. Sergienko's 80th Birthday*, pages 217–237, 2017.

14 Seyedmohammadhossein Hosseinian, DBMM Fontes, and Sergiy Butenko. A quadratic approach to the maximum edge weight clique problem. In *XIII Global Optimization Workshop GOW*, volume 16, pages 125–128, 2016.

15 Marcel Hunting, Ulrich Faigle, and Walter Kern. A lagrangian relaxation approach to the edge-weighted clique problem. *European Journal of Operational Research*, 131(1):119–131, 2001.

16 Ronan Le Bras, Carla Gomes, and Bart Selman. From streamlined combinatorial search to efficient constructive procedures. In *Proceedings of the AAAI Conference on Artificial Intelligence*, volume 26-1, pages 499–506, 2012.

17 Ronan Le Bras, Andrew Perrault, and Carla P Gomes. Polynomial time construction for spatially balanced latin squares. Technical report, Cornell University, 2012.

18 Anuj Mehrotra. Cardinality constrained boolean quadratic polytope. *Discrete Applied Mathematics*, 79(1-3):137–154, 1997.

19 Michael M Sørensen. New facets and a branch-and-cut algorithm for the weighted clique problem. *European Journal of Operational Research*, 154(1):57–70, 2004.

20 Etsuji Tomita, Tatsuya Akutsu, and Tsutomu Matsunaga. *Efficient algorithms for finding maximum and maximal cliques: Effective tools for bioinformatics*. IntechOpen, 2011.

21 HM Van Es and CL Van Es. Spatial nature of randomization and its effect on the outcome of field experiments. *Agronomy Journal*, 85(2):420–428, 1993.

22 Pascal Van Hentenryck and Laurent Michel. Differentiable invariants. In *Proceedings of Principles and Practice of Constraint Programming: 12th International Conference*, pages 604–619. Springer, 2006.

A Supplementary Results

■ **Table 3** Runtime in seconds for single instances of the clique formulation run in Gurobi and CPLEX using a single thread for $n = 8 - 13$. The dash indicates that the model was not able to prove optimality in 6 hours (21,600 s). For all $n$ less than 8, both models ran in under one second.

<table><thead><tr><th></th><th colspan="6">n</th></tr><tr><th></th><th>8</th><th>9</th><th>10</th><th>11</th><th>12</th><th>13</th></tr></thead><tbody><tr><th>CPLEX</th><td><strong>1.10</strong></td><td><strong>4.48</strong></td><td><strong>28.76</strong></td><td>675.16</td><td>-</td><td>-</td></tr><tr><th>Gurobi</th><td>1.15</td><td>4.99</td><td>51.13</td><td><strong>171.22</strong></td><td><strong>614.22</strong></td><td><strong>4750.96</strong></td></tr></tbody></table>
<table><caption>Table 4 Comparison of the greatest lower bound found by single instances of the clique formulation run in Gurobi and CPLEX using a single thread for 6 hours for $n = 14 - 20$.</caption><thead><tr><th></th><th colspan="7">n</th></tr><tr><th></th><th>14</th><th>15</th><th>16</th><th>17</th><th>18</th><th>19</th><th>20</th></tr></thead><tbody><tr><th>CPLEX</th><td>132</td><td>119.29</td><td>97.54</td><td>88</td><td>47.79</td><td>40.06</td><td>45.06</td></tr><tr><th>Gurobi</th><td><strong>191</strong></td><td><strong>252.66</strong></td><td><strong>273.66</strong></td><td><strong>252</strong></td><td><strong>402.66</strong></td><td><strong>465.33</strong></td><td><strong>155</strong></td></tr></tbody></table>

<table><caption>Table 5 Comparison of the smallest upper bound found by single instances of the clique formulation run in Gurobi and CPLEX using a single thread for 6 hours for $n = 14 - 20$.</caption><thead><tr><th></th><th colspan="7">n</th></tr><tr><th></th><th>14</th><th>15</th><th>16</th><th>17</th><th>18</th><th>19</th><th>20</th></tr></thead><tbody><tr><th>CPLEX</th><td><strong>272</strong></td><td>354</td><td>437.33</td><td>526</td><td>631.33</td><td>762</td><td>878</td></tr><tr><th>Gurobi</th><td><strong>272</strong></td><td><strong>345.33</strong></td><td><strong>427.33</strong></td><td><strong>522</strong></td><td><strong>617.33</strong></td><td><strong>732</strong></td><td><strong>876</strong></td></tr></tbody></table>

<table><caption>Table 6 Comparison of the greatest lower bound found by three Gurobi configurations run for 6 hours for $n = 14 - 20$. The dashed cell indicates a run that failed due to numerical instability.</caption><thead><tr><th></th><th colspan="7">n</th></tr><tr><th>Instances</th><th>Threads/Inst.</th><th>14</th><th>15</th><th>16</th><th>17</th><th>18</th><th>19</th><th>20</th></tr></thead><tbody><tr><td>1</td><td>1</td><td>191</td><td><strong>252.66</strong></td><td>273.66</td><td><strong>252</strong></td><td><strong>402.66</strong></td><td><strong>465.33</strong></td><td><strong>155</strong></td></tr><tr><td>8</td><td>1</td><td>159.69</td><td>185.28</td><td>203.11</td><td>198.56</td><td>-</td><td>389.50</td><td>92.5</td></tr><tr><td>1</td><td>32</td><td><strong>201.57</strong></td><td>249.79</td><td><strong>273.82</strong></td><td>240.90</td><td>185.75</td><td>371.95</td><td>28</td></tr></tbody></table>

<table><caption>Table 7 Comparison of the smallest upper bound found by three Gurobi configurations run for 6 hours for $n = 14 - 20$. The dashed cell indicates a run that failed due to numerical instability.</caption><thead><tr><th></th><th colspan="7">n</th></tr><tr><th>Instances</th><th>Threads/Inst.</th><th>14</th><th>15</th><th>16</th><th>17</th><th>18</th><th>19</th><th>20</th></tr></thead><tbody><tr><td>1</td><td>1</td><td>272</td><td><strong>345.33</strong></td><td><strong>427.33</strong></td><td><strong>522</strong></td><td><strong>617.33</strong></td><td><strong>732</strong></td><td>886</td></tr><tr><td>8</td><td>1</td><td><strong>270</strong></td><td>358.66</td><td>442.66</td><td>530</td><td>-</td><td>747.33</td><td><strong>816</strong></td></tr><tr><td>1</td><td>32</td><td>276</td><td>347.33</td><td>430.66</td><td><strong>522</strong></td><td>631.33</td><td>766</td><td>900</td></tr></tbody></table>

<table><caption>Table 8 Runtime in seconds for single instances of the clique and assignment formulations run in Gurobi using a single thread and A* search for $n = 8 - 13$. The dash indicates that the model was not able to prove optimality in 6 hours (21,600 s). For all $n$ less than 8, all models finished in under one second. The random walk is not included in this table as the method does not prove optimality.</caption><thead><tr><th></th><th colspan="6">n</th></tr><tr><th></th><th>8</th><th>9</th><th>10</th><th>11</th><th>12</th><th>13</th></tr></thead><tbody><tr><th>Assignment MIP</th><td>1.09</td><td>6.23</td><td>59.02</td><td>678.39</td><td>3790.28</td><td>-</td></tr><tr><th>Clique MIP</th><td>1.15</td><td>4.99</td><td>51.13</td><td>171.22</td><td><strong>614.22</strong></td><td><strong>8573.43</strong></td></tr><tr><th>A*</th><td><strong>0.07</strong></td><td><strong>0.78</strong></td><td><strong>8.72</strong></td><td><strong>98.70</strong></td><td>1299.66</td><td>-</td></tr></tbody></table>
<figure><img src="image_464.png" /></figure>
# Proven Optimally-Balanced Latin Rectangles with SAT

Vaidyanathan Peruvemba Ramaswamy

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

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

Algorithms and Complexity Group, TU Wien, Austria

Stefan Szeider

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

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

Algorithms and Complexity Group, TU Wien, Austria

## Abstract

Motivated by applications from agronomic field experiments, Díaz, Le Bras, and Gomes [CPAIOR 2015] introduced Partially Balanced Latin Rectangles as a generalization of Spatially Balanced Latin Squares. They observed that the generation of Latin rectangles that are optimally balanced is a highly challenging computational problem. They computed, utilizing CSP and MIP encodings, Latin rectangles up to $12 \times 12$, some optimally balanced, some suboptimally balanced.

In this paper, we develop a SAT encoding for generating balanced Latin rectangles. We compare experimentally encoding variants. Our results indicate that SAT encodings perform competitively with the MIP encoding, in some cases better. In some cases we could find Latin rectangles that are more balanced than previously known ones. This finding is significant, as there are many arithmetic constraints involved. The SAT approach offers the advantage that we can certify that Latin rectangles are optimally balanced through DRAT proofs that can be verified independently.

**2012 ACM Subject Classification** Hardware → Theorem proving and SAT solving; Theory of computation → Constraint and logic programming; Mathematics of computing → Combinatoric problems; Software and its engineering → Constraints; Theory of computation → Integer programming; Mathematics of computing → Combinatorial optimization; Mathematics of computing → Enumeration; Mathematics of computing → Solvers; Computing methodologies → Theorem proving algorithms; Applied computing → Agriculture

**Keywords and phrases** combinatorial design, SAT encodings, certified optimality, arithmetic constraints, spatially balanced Latin rectangles

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

**Category** Short Paper

**Supplementary Material** Software (Script, Encodings and Models): https://doi.org/10.5281/zenodo.8151905

**Funding** Supported by the Austrian Science Funds (FWF) within the project P36420, and the Vienna Science and Technology Fund (WWTF) within the project ICT19-065.

**Acknowledgements** This work was carried out in part while the second author visited the Simons Institute for the Theory of Computing, University of Berkeley, within the program *Extended Reunion: Satisfiability*.

## 1 Introduction

A *Latin square* is an $n \times n$ array filled with $n$ different symbols, each occurring exactly once in each row and exactly once in each column. The notion goes back to the 18th century and Leonard Euler, who studied the problem of generating Latin squares. More recently, additional constraints have been added, making the combinatorial design problem even more challenging. An interesting additional constraint asks for a Latin square to be spatially balanced, i.e., any two symbols $u, v$ have the same distance. The distance of $u$ and $v$ is

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

© Vaidyanathan Peruvemba Ramaswamy and Stefan Szeider; 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. 48; pp. 48:1–48:10

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

LIPICS

Leibniz International Proceedings in Informatics

Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl Publishing, Germany
48:2 Proven Optimally-Balanced Latin Rectangles with SAT

defined as the sum of their distances over all the rows. A polynomial-time algorithm is
known that constructs for any given *n* such that 2*n* + 1 is prime a spatially balanced Latin
square [10].

Díaz, Le Bras, and Gomes [3] introduced *spatially balanced Latin rectangles*, which are
arrays with *n* columns and *k* rows filled with *n* symbols, such that no symbol occurs more
than once in any row or column, and again the distance of any pair of symbols is the same.
As it turned out that spatially balanced Latin rectangles only exist for a relatively restricted
set of combinations of *n* and *k*, Díaz et al. defined a notation of *imbalance* and asked for Latin
rectangles where the imbalance is minimal. In a spatially balanced *k* × *n* Latin rectangle,
the distance between any two symbols is *k*(*n* + 1)/3. Therefore, the imbalance of a pair of
symbols *i*, *j* is naturally defined as I(*u*, *v*) = |*d*(*u*, *v*) - *k*(*n* + 1)/3|. The *imbalance* of a Latin
rectangle *L* is then I(*L*) = Σ*u*<*v* I(*u*, *v*) = Σ*u*<*v* |*d*(*u*, *v*) - *k*(*n* + 1)/3|. See Figure 1b for a
sample working out of this quantity. Latin squares are of general importance for the design of
agronomic experiments [9]. More specifically, however, the balanced version is essential to the
design of bias-free agronomic experiments, as they avoid unintentional patterns introduced
due to spatial auto-correlation [13].

Díaz et al. computed Latin rectangles up to 12 × 12, trying to minimize the total imbalance. They used CSP, local search, and MIP methods, with the latter being the most efficient. In some cases, they could obtain optimal imbalance; in other cases, upper bounds. They conclude that finding Latin rectangles with minimum imbalance seems to be a very challenging computational problem and suggest it as an ideal benchmark for different search and optimization approaches.

In this paper, we follow this insight and consider the research question of whether propositional satisfiability (SAT) solvers [4] are competitive in finding optimally balanced Latin rectangles. This question is particularly interesting as the definition of a partially balanced Latin rectangle entails many arithmetic constraints, including addition, subtraction, and absolute values, which adds to the challenge. A SAT approach that is competitive with MIP brings the added value of certification: one can certify that a Latin rectangle is optimally balanced through a DRAT proof that can be checked independently. This provides an additional layer of confidence to the results.

In brief, the SAT encoding developed by us represents the position of each symbol on a row
by a set of propositional variables. The assignment of these variables is then propagated by
means of auxiliary variables and clauses through the rest of the encoding. This propagation
terminates at a set of variables which represent the total imbalance of the Latin rectangle
expressed in unary. Bounding these variables effectively bounds the imbalance. Special care
must be taken to encode the absolute value arithmetic and the nested summations with the
auxiliary variables and clauses.

From our experiments, we observe that the SAT approach is able to replicate 44 out of
the 63 upper bounds shown by the MIP approach and in 14 cases even find a tighter upper
bound. From the infeasibility side as well, the SAT approach is comparable to the MIP
approach and is able to confirm almost all the optimality results with the added bonus of
producing DRAT proofs with some minor overhead.

2 Preliminaries

We denote the set of integers {1, ..., n} by [n]. For positive integers *k* ≤ *n*, a *k* × *n*
*Latin rectangle* is a matrix *L* with *k* rows and *n* columns with elements from [n] such that
*L(i,j) ≠ L(i',j')* whenever *i = i' ∧ j ≠ j' or i ≠ i' ∧ j = j'*.
Consider a $k \times n$ Latin rectangle $L$ and $i, j \in [n]$. The row $i$ distance between $u, v$, if $L(i, c_u) = u$ and $L(i, c_v) = v$, is

$$d_i(u, v) := |c_u - c_v|. \qquad (1)$$

The distance between $u$ and $v$ in $L$ is

$$d(u, v) := \sum_{i \in [k]} d_i(u, v), \qquad (2)$$

and the imbalance of the pair $u, v$ is

$$\mathbb{I}(u, v) := |d(u, v) - k(n + 1)/3|. \qquad (3)$$

The imbalance of a Latin rectangle $L$ is the sum of the pairwise imbalances, i.e.,

$$\mathbb{I}(L) := \sum_{u<v} \mathbb{I}(u, v). \qquad (4)$$

We are interested in finding Latin rectangles with low imbalance. A Latin rectangle $L$ is optimally balanced (or partially spatially balanced as Díaz et al. [3] call it) if $\mathbb{I}(L)$ is minimum and suboptimally balanced otherwise. See Figure 1a for an example.

<figure><img src="image_471.png" /><figcaption>(a) 3 × 5 optimally balanced Latin rectangle with imbalance 6.</figcaption></figure>

<table><caption>(b) Breakdown (for each pair) of the total imbalance, where $Z = k(n + 1)/3 = 6$. For illustration, $d(1, 2) = d_1(1, 2) + d_2(1, 2) + d_3(1, 2)$ which is $1 + 2 + 2 = 5$.</caption><thead><tr><th>u, v</th><th>d(u, v)</th><th>|d(u, v) - Z|</th><th>u, v</th><th>d(u, v)</th><th>|d(u, v) - Z|</th></tr></thead><tbody><tr><td>1, 2</td><td>5</td><td>1</td><td>2, 4</td><td>9</td><td>3</td></tr><tr><td>1, 3</td><td>6</td><td>0</td><td>2, 5</td><td>5</td><td>1</td></tr><tr><td>1, 4</td><td>6</td><td>0</td><td>3, 4</td><td>6</td><td>0</td></tr><tr><td>1, 5</td><td>6</td><td>0</td><td>3, 5</td><td>6</td><td>0</td></tr><tr><td>2, 3</td><td>5</td><td>1</td><td>4, 5</td><td>6</td><td>0</td></tr></tbody></table>

<caption>Figure 1 Example of an optimally balanced 3 × 5 Latin rectangle.</caption>

## 3 Encoding

In this section, we describe the details of the SAT encoding. We describe the clauses that are conjuncted together to form the CNF formula $\mathcal{F}(k, n, t)$. This formula is satisfiable if and only if there exists a $k \times n$ Latin rectangle $L$ such that $\mathbb{I}(L) \le t$. For the sake of convenience, we sometimes shorten $\mathcal{F}(n, k, t)$ to $\mathcal{F}$, and use $\Delta_n$ to denote all ordered pairs $(u, v)$ such that $u < v \in [n]$. We also use $Z := k(n + 1)/3$, assuming that 3 divides $k(n + 1)$. We address the case of indivisibility later. Further, we make extensive use of higher-level cardinality constraints written as $\{v_i | i \in [\ell]\}_{\le k}$, $\{v_i | i \in [\ell]\}_{= k}$, or $\{v_i | i \in [\ell]\}_{\ge k}$ to encode the condition that only at most $k$, exactly $k$, or at least $k$ variables, respectively, from the set $\{v_1, \dots, v_\ell\}$ are set to true. Each higher-level constraint can then be rewritten as a conjunction of several elementary clauses utilizing auxiliary variables.

The core variables in $\mathcal{F}(n, k, t)$ are $p(i, u, j)$ for $u, j \in [n], i \in [k]$ which are true if and only if the column where element $u$ appears in the $i$th row is $j$. The following cardinality constraint encodes the condition that each cell of the matrix contains exactly one element:

$$\bigwedge_{i \in [k], j \in [n]} \{p(i, u, j) | u \in [n]\} = 1$$
48:4 Proven Optimally-Balanced Latin Rectangles with SAT

The following cardinality constraints encode the requirement that each element appears
exactly once within a row:

$$
\bigwedge_{i \in [k]} \{ p(i, u, j) \mid j \in [n] \}_{\ge 1} \wedge \{ p(i, u, j) \mid j \in [n] \}_{\le 1}
$$

Lastly, the following cardinality constraint encodes the requirement that no element appears
twice within the same column:

$$
\bigwedge_{u,j \in [n]} \{ p(i, u, j) \mid i \in [k] \}_{\le 1}
$$

Then, to capture Equation (1), we introduce the auxiliary variables $c(i, u, v, d)$ which are
true if $d_i(u, v) \ge d$ and constraint them using the following three sets of clauses:

$$
\bigwedge_{\substack{i \in [k] \\ (u,v) \in \Delta_n \\ j_1 \neq j_2 \in [n]}} (p(i, u, j_1) \wedge p(i, v, j_2)) \rightarrow c(i, u, v, |j_1 - j_2|)
$$

$$
\bigwedge_{\substack{i \in [k] \\ (u,v) \in \Delta_n \\ j_1 \neq j_2 \in [n]}} (p(i, u, j_1) \wedge p(i, v, j_2)) \rightarrow \neg c(i, u, v, |j_1 - j_2| + 1)
$$

$$
\bigwedge_{\substack{i \in [k] \\ (u,v) \in \Delta_n \\ d \in [n-1]}} c(i, u, v, d) \rightarrow c(i, u, v, d-1)
$$

Then we represent the inner sum $d(u, v) = \sum_i d_i(u, v)$ from Equation (2) using the two
sets of auxiliary variables; $s_g(u, v, p)$ which when falsified enforces the condition $d(u, v) \le p$
and similarly $s_l(u, v, p)$ which when falsified enforces the condition $\sum_i d_i(u, v) \ge p$. This is
achieved using the following conditional cardinality constraints:

$$
\bigwedge_{\substack{i \in [k] \\ (u,v) \in \Delta_n \\ k \le p \le k(n-1)}} \neg s_g(u, v, p) \rightarrow \{c(i, u, v, d) \mid d \in [n-1]\}_{\le p}
$$

$$
\bigwedge_{\substack{i \in [k] \\ (u,v) \in \Delta_n \\ k \le p \le k(n-1)}} \neg s_l(u, v, p) \rightarrow \{c(i, u, v, d) \mid d \in [n-1]\}_{\ge p}
$$

Recall from Equations (3) and (4), that our final goal is to bound the quantity $\mathbb{I}(L) = \sum_{(u,v) \in \Delta_n} |d(u,v) - Z|$, where $Z = k(n+1)/3$. The redundancy in the form of two sets of variables $s_g$ and $s_l$ is to allow us to later deal with the expression $|d(u,v) - Z|$ and help count up from $Z$ in either direction. We thus introduce auxiliary variables which encode the sign of the term $d(u,v) - Z$. $s_0(u,v)$ is true if and only if $d(u,v) = Z$; $s_+(u,v)$ is true if $d(u,v) > Z$, false if $d(u,v) < Z$, and undefined otherwise. The following sets of clauses encode these auxiliary variables:

$$
\bigwedge_{(u,v) \in \Delta_n} s_0(u,v) \rightarrow \neg s_g(u,v,Z) \wedge \neg s_l(u,v,Z)
$$

$$
\bigwedge_{(u,v) \in \Delta_n} (\neg s_0(u,v) \wedge s_+(u,v)) \rightarrow \neg s_l(u,v,Z+1)
$$
$$ \bigwedge_{(u,v) \in \Delta_n} (\neg s_0(u, v) \wedge \neg s_+(u, v)) \rightarrow \neg s_g(u, v, Z-1) $$

Then, using the auxiliary variables $f(u, v, q)$, we represent the condition $|d(u, v) - Z| \le q$. We encode these variables using the following set of clauses:

$$ \bigwedge_{\substack{(u,v) \in \Delta_n \\ 1 \le q \le k(n-1)-Z}} (\neg s_0(u, v) \wedge s_+(u, v) \wedge s_g(u, v, Z+q)) \rightarrow f(u, v, q) $$

$$ \bigwedge_{\substack{(u,v) \in \Delta_n \\ 1 \le q \le Z-k}} (\neg s_0(u, v) \wedge \neg s_+(u, v) \wedge s_l(u, v, Z-q)) \rightarrow f(u, v, q) $$

And finally, we express the bound $\mathbb{I}(L) \le t$ by adding the cardinality constraint:

$$ \{ f(u, v, q) \mid (u, v) \in \Delta_n, 1 \le q \le \max(k(n-1) - Z, Z - k) \}_{\le t} $$

## Symmetry Breaking

We refer to the encoding so far as the *base version* and denote the formula by $\mathcal{F}(n, k, t)$. We optionally add two sets of *symmetry-breaking* clauses which reduce the size of the search space by identifying equivalent solutions. These constraints are identical to those mentioned by Díaz et al. [3], however, in our approach, we selectively enable or disable these constraints depending on whether we are seeking a SAT instance (i.e., upper bound) or an UNSAT instance (i.e., infeasibility). The following set of clauses fixes the first row to be the identity permutation:

$$ \bigwedge_{u \in [n]} p(1, u, u) $$

We refer to the encoding with these additional clauses as $\mathcal{F}'(n, k, t)$. Similarly, the following set of clauses forces the first column to follow a strictly increasing order, i.e., a lexicographic order:

$$ \bigwedge_{(u,v), i \in [k-1]} p(i, v, 1) \rightarrow \neg p(i+1, u, 1) $$

We refer to the encoding with both these sets of additional clauses as $\mathcal{F}''(n, k, t)$. These two sets of symmetry breaking clauses enforce that the Latin rectangle is in its *reduced form* [3].

## Non-integral Instances

In the description of the encoding above, we made the assumption that 3 divides $k(n+1)$. In this case, we do not need to deal with fractional $\frac{1}{3}$ values. We call these the *integral instances* and all other instances *non-integral instances*. To deal with the non-integral instances, we modify our encoding by multiplying the distances by 3. To this end, we use $Z = k(n+1)$ instead of $Z = k(n+1)/3$. We use $3|j_1 - j_2|$ instead of $|j_1 - j_2|$ and we also multiply the bounds for $d$, $p$, and $q$ by 3.

# 4 Experimental Evaluation

To analyze the performance of our proposed encodings, we implemented a method for generating the encodings in Python 3.10.6 with help from the PySAT 0.1.7 library [8] for handling some of the higher-level cardinality constraints [12, 1, 5]. We then ran the SAT solver
kissat [2] with default settings on the generated encodings. We ran all our experiments on a 10-core Intel Xeon E5-2640 v4, 2.40 GHz CPU, with each process having access to 8GB RAM.

kissat, like other state-of-the-art CDCL solvers, is capable of producing a proof in the DRAT (Deletion Resolution Asymmetric Tautology) format for formulas which it claims to be unsatisfiable. Similar to how a mathematical proof of a theorem follows a sequence of smaller lemmas, a DRAT proof also consists of lemmas (in the form of a set of literals) and lemma deletion instructions. Due to the elementary nature of the proof format, it is straightforward for a CDCL solver to emit a DRAT proof with minimal overhead. This proof can then be independently verified using tools such as DRAT-trim [6, 15]. DRAT-trim takes as input the original CNF formula and the DRAT proof produced by the SAT solver and can verify, in time comparable to original proof producing time, that the proof indeed holds.

There are several ways of translating the higher-level cardinality constraints into elementary CNF clauses. Thus, we conducted a preliminary investigation to identify the most promising translations from those provided by the PySAT library<fnref n="1" />. We also tested the three versions of the encoding with different levels of symmetry breaking $F, F', F''$ in combination with the different cardinality encoding types. From this investigation, we found that $F'$ combined with the bitwise cardinality encoding [12] and $F''$ combined with the ladder cardinality encoding [1] were the two best performing encodings for SAT and UNSAT instances, respectively. Hence, we stuck with these two combinations for the main experiments. In this preliminary investigation, we worked with instances which were already known to be either SAT or UNSAT. In the main experiments, typically, we don't know beforehand if the instance is SAT or UNSAT, thus we try both these encodings.

This paper focuses on the obtained bounds and their verification rather than on exact running times; the results reported by Díaz et al. serve as a basic comparison of the runtime performance with other methods. We conducted cursory experiments with a CSP model formulated with MiniZinc [11] for a direct and reproducible comparison on the same hardware. We tested on the solvers Gecode 6.3.0 and Chuffed 0.11.0, the former outperforming the latter. On small rectangles, up to around 5 × 7, the CSP approach is significantly faster than our SAT approach, but its performance quickly deteriorates as the size grows. For example, the CSP approach fails to replicate known bounds for 6 × 8 within 10 hours even though, on the other hand, our approach produced solutions in under an hour.

Apart from the prototyping experiments, we conducted two main experiments rigorously. The first was for showing upper bounds and consequently generating Latin rectangles matching those bounds, and the second was to show lower bounds and consequently generating DRAT proofs of infeasibility. We worked with the same set of instances as Díaz et al. and used a timeout of 24 hours. We used our proposed encodings to replicate their upper bounds, produce verifiable proofs for their lower bounds, and also in some cases improve their upper bounds. The results of these experiments are summarized in Table 1. Table 2 shows some attributes of the proofs generated by our method for a representative set of instances. Figure 2 shows some of the Latin rectangles that witness the new upper bounds shown by us. In the supplementary material, we include the script used to generate the encodings and some generated encodings along with their corresponding models (if satisfiable) or DRAT proofs (if unsatisfiable).

<footnote n="1">¹ List of cardinality encoding types: pysat.card.EncType</footnote>
**Table 1** Table of results. The top row indicates the value of *k* and the left column indicates the value of *n*. Grey cells indicate values claimed optimal by Díaz et al. using CSP/MIP. The cells with bold text indicate upper bounds that we could replicate. Cells with a '*' indicate lower bounds that we could confirm and certify with a DRAT proof (not applicable for rectangles with minimum imbalance 0). Green cells indicate an earlier upper bound (strike-through text) that we could improve.

<table><thead><tr><th></th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th><th>12</th></tr></thead><tbody><tr><th>4</th><td>2.6*</td><td>4*</td><td>5.3*</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><th>5</th><td>8*</td><td>6*</td><td>8*</td><td>0</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><th>6</th><td>16*</td><td>12*</td><td>13.3*</td><td>16*</td><td>0</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><th>7</th><td>28*</td><td>22*</td><td>22.6*</td><td>22.6*</td><td>20*</td><td>18.6*</td><td></td><td></td><td></td><td></td><td></td></tr><tr><th>8</th><td>40*</td><td>36*</td><td>32</td><td>30</td><td>24</td><td>28</td><td>0</td><td></td><td></td><td></td><td></td></tr><tr><th>9</th><td>65.3*</td><td>56</td><td>56.6<br>56.0</td><td>56<br>54</td><td>52<br>48</td><td>66<br>64.6</td><td>60<br>58.6</td><td>0</td><td></td><td></td><td></td></tr><tr><th>10</th><td>92*</td><td>86<br>84</td><td>92<br>91.3</td><td>66.6</td><td>102<br>96</td><td>100<br>99.3</td><td>99.3</td><td>80</td><td>40</td><td></td><td></td></tr><tr><th>11</th><td>124*</td><td>120<br>118</td><td>122<br>118</td><td>122<br>120</td><td>126<br>124</td><td>136<br>132</td><td>132</td><td>128</td><td>110</td><td>0</td><td></td></tr><tr><th>12</th><td>168</td><td>158</td><td>162.6</td><td>170.6</td><td>120</td><td>183</td><td>184.6</td><td>178</td><td>174.6</td><td>147.3</td><td>0</td></tr></tbody></table>

**Table 2** Encoding size for some representative instances along with time required to prove optimality and the size of the generated DRAT proof. All these proofs were generated by running the SAT solver with the `--unsat` preset on $F''$ combined with the `ladder` cardinality encoding. Notice that the encoding size for the non-integral instance $4 \times 6$ is much bigger than that of the integral instance $6 \times 7$.

<table><thead><tr><th>$k \times n$</th><th>#variables</th><th>#clauses</th><th>Unsat Time</th><th>Proof Overhead</th><th>Proof Size</th></tr></thead><tbody><tr><td>3 \times 5</td><td>2364</td><td>5750</td><td>&lt;1 s</td><td>&lt;1 s</td><td>20 KB</td></tr><tr><td>4 \times 6</td><td>73184</td><td>147903</td><td>45 s</td><td>2 s</td><td>11 MB</td></tr><tr><td>6 \times 7</td><td>37809</td><td>84906</td><td>1 hr 12 min</td><td>1 hr 45 min</td><td>1.6 GB</td></tr></tbody></table>
48:8 Proven Optimally-Balanced Latin Rectangles with SAT

<table><caption>(a) 3 × 10 rectangle with imbalance 84 (32 min).</caption><tr><td>1</td><td>2</td><td>3</td><td>4</td><td>5</td><td>6</td><td>7</td><td>8</td><td>9</td><td>10</td></tr><tr><td>4</td><td>10</td><td>2</td><td>8</td><td>7</td><td>5</td><td>1</td><td>6</td><td>3</td><td>9</td></tr><tr><td>7</td><td>3</td><td>8</td><td>1</td><td>4</td><td>9</td><td>2</td><td>10</td><td>6</td><td>5</td></tr></table>

<table><caption>(c) 6 × 10 rectangle with imbalance 96 (40 min).</caption><tr><td>1</td><td>2</td><td>3</td><td>4</td><td>5</td><td>6</td><td>7</td><td>8</td><td>9</td><td>10</td></tr><tr><td>5</td><td>8</td><td>6</td><td>1</td><td>9</td><td>3</td><td>10</td><td>4</td><td>2</td><td>7</td></tr><tr><td>6</td><td>5</td><td>9</td><td>3</td><td>2</td><td>7</td><td>4</td><td>10</td><td>1</td><td>8</td></tr><tr><td>7</td><td>6</td><td>1</td><td>5</td><td>10</td><td>2</td><td>8</td><td>9</td><td>4</td><td>3</td></tr><tr><td>9</td><td>1</td><td>7</td><td>10</td><td>4</td><td>5</td><td>6</td><td>3</td><td>8</td><td>2</td></tr><tr><td>10</td><td>3</td><td>5</td><td>7</td><td>1</td><td>8</td><td>9</td><td>2</td><td>6</td><td>4</td></tr></table>

<table><caption>(b) 4 × 11 rectangle with imbalance 118 (102 min).</caption><tr><td>1</td><td>2</td><td>3</td><td>4</td><td>5</td><td>6</td><td>7</td><td>8</td><td>9</td><td>10</td><td>11</td></tr><tr><td>2</td><td>6</td><td>5</td><td>7</td><td>11</td><td>8</td><td>1</td><td>3</td><td>10</td><td>9</td><td>4</td></tr><tr><td>3</td><td>5</td><td>11</td><td>8</td><td>9</td><td>2</td><td>10</td><td>4</td><td>1</td><td>6</td><td>7</td></tr><tr><td>5</td><td>1</td><td>9</td><td>10</td><td>7</td><td>3</td><td>6</td><td>2</td><td>11</td><td>4</td><td>8</td></tr></table>

<table><caption>(d) 6 × 11 rectangle with imbalance 124 (5.5 hr).</caption><tr><td>1</td><td>2</td><td>3</td><td>4</td><td>5</td><td>6</td><td>7</td><td>8</td><td>9</td><td>10</td><td>11</td></tr><tr><td>2</td><td>11</td><td>7</td><td>6</td><td>4</td><td>5</td><td>1</td><td>9</td><td>8</td><td>3</td><td>10</td></tr><tr><td>3</td><td>7</td><td>9</td><td>5</td><td>10</td><td>2</td><td>6</td><td>1</td><td>11</td><td>4</td><td>8</td></tr><tr><td>5</td><td>1</td><td>11</td><td>10</td><td>7</td><td>3</td><td>8</td><td>6</td><td>4</td><td>2</td><td>9</td></tr><tr><td>6</td><td>10</td><td>1</td><td>3</td><td>8</td><td>7</td><td>2</td><td>11</td><td>5</td><td>9</td><td>4</td></tr><tr><td>7</td><td>4</td><td>10</td><td>1</td><td>9</td><td>11</td><td>3</td><td>2</td><td>6</td><td>8</td><td>5</td></tr></table>

■ **Figure 2** Examples of Latin rectangles witnessing new upper bounds along with the time required by the SAT solver to find these solutions in parentheses. All these rectangles were found by the version of the encoding that combines $\mathcal{F}''$ with the **ladder** cardinality encoding. Note that, no lower bounds are known for these sizes.

5 Conclusion

This work presented a SAT approach to the balanced Latin rectangle problem treating complex arithmetic expressions by translating them into propositional logic. Our approach is competitive with the prior CSP/MIP-based techniques with a slight advantage. On the one hand, the SAT approach can find tighter upper bounds for realizing Latin rectangles. On the other hand, our approach performs comparably to the MIP approach when it comes to proving the infeasibility of particular Latin rectangles. This finding is significant, as propositional logic is often considered inferior to MIP when handling complex arithmetic expressions. A key advantage of our SAT approach is its ability to produce DRAT proofs of optimality that can be independently verified. The proof generation causes only some minor computational overhead. On a more technical level, our work shows how to handle non-trivial nested summations in a SAT encoding by chaining ladders and cardinality counters. We observe that it is sometimes beneficial to construct two different encodings tailored towards SAT and UNSAT instances in problems where one is interested in both upper and lower bounds. For instance, these encodings can use different levels of symmetry-breaking and different types of cardinality encodings.

We see many potential avenues for future work. One could use techniques such as cubing [7] to tackle the instances for which proving optimality is still challenging. The cubing can be performed by either fixing the values of the first column, fixing the imbalance of certain pairs to be 0, or bounding the number of pairs with 0 imbalance. We can enumerate all unique (up to isomorphism) Latin rectangles of a particular size and imbalance by extending the SAT approach to an incremental solving approach. Lastly, another natural direction is to apply Pseudo-Boolean solving to this problem. Pseudo-Boolean solving can also produce cutting planes proofs for verification [14] and might be even better suited to handle the arithmetic constraints involved in this problem.
References

1 Carlos Ansótegui and Felip Manyà. Mapping problems with finite-domain variables to problems with boolean variables. In Holger H. Hoos and David G. Mitchell, editors, *Theory and Applications of Satisfiability Testing, 7th International Conference, SAT 2004, Vancouver, BC, Canada, May 10-13, 2004, Revised Selected Papers*, volume 3542 of *Lecture Notes in Computer Science*, pages 1–15. Springer Verlag, 2004. doi:10.1007/11527695_1.

2 Armin Biere, Katalin Fazekas, Mathias Fleury, and Maximillian Heisinger. CaDiCaL, Kissat, Paracooba, Plingeling and Treengeling entering the SAT Competition 2020. In Tomas Balyo, Nils Froleyks, Marijn Heule, Markus Iser, Matti Järvisalo, and Martin Suda, editors, *Proc. of SAT Competition 2020 - Solver and Benchmark Descriptions*, volume B-2020-1 of *Department of Computer Science Report Series B*, pages 51–53. University of Helsinki, 2020. URL: https://tuhat.helsinki.fi/ws/files/142452772/sc2020_proceedings.pdf.

3 Mateo Díaz, Ronan Le Bras, and Carla P. Gomes. In search of balance: The challenge of generating balanced Latin rectangles. In Domenico Salvagnin and Michele Lombardi, editors, *Integration of AI and OR Techniques in Constraint Programming - 14th International Conference, CPAIOR 2017, Padua, Italy, June 5-8, 2017, Proceedings*, volume 10335 of *Lecture Notes in Computer Science*, pages 68–76. Springer, 2017. doi:10.1007/978-3-319-59776-8_6.

4 Johannes K. Fichte, Daniel Le Berre, Markus Hecher, and Stefan Szeider. The silent (r)evolution of SAT. *Communications of the ACM*, 66(6):64–72, June 2023. doi:10.1145/3560469.

5 Ian P. Gent and Peter Nightingale. A new encoding of alldifferent into SAT. In *International Workshop on Modelling and Reformulating Constraint Satisfaction*, volume 3, pages 95–110, 2004.

6 Marijn Heule, Warren A. Hunt Jr., Matt Kaufmann, and Nathan Wetzler. Efficient, verified checking of propositional proofs. In Mauricio Ayala-Rincón and César A. Muñoz, editors, *Interactive Theorem Proving - 8th International Conference, ITP 2017, Brasília, Brazil, September 26-29, 2017, Proceedings*, volume 10499 of *Lecture Notes in Computer Science*, pages 269–284. Springer Verlag, 2017. doi:10.1007/978-3-319-66107-0_18.

7 Marijn J. H. Heule, Oliver Kullmann, and Victor W. Marek. Solving very hard problems: Cube-and-conquer, a hybrid SAT solving method. In Carles Sierra, editor, *Proceedings of the Twenty-Sixth International Joint Conference on Artificial Intelligence, IJCAI 2017, Melbourne, Australia, August 19-25, 2017*, pages 4864–4868. ijcai.org, 2017. doi:10.24963/ijcai.2017/683.

8 Alexey Ignatiev, Antonio Morgado, and Joao Marques-Silva. PySAT: A Python toolkit for prototyping with SAT oracles. In *SAT*, pages 428–437, 2018. doi:10.1007/978-3-319-94144-8_26.

9 Marcus Jones, Richard Woodward, and Jerry Stoller. Increasing precision in agronomic field trials using latin square designs. *Agronomy Journal*, 107(1):20–24, 2015. doi:10.2134/agronj14.0232.

10 Ronan LeBras, Carla P. Gomes, and Bart Selman. From streamlined combinatorial search to efficient constructive procedures. In Jörg Hoffmann and Bart Selman, editors, *Proceedings of the Twenty-Sixth AAAI Conference on Artificial Intelligence, July 22-26, 2012, Toronto, Ontario, Canada*. AAAI Press, 2012. doi:10.1609/aaai.v26i1.8147.

11 Nicholas Nethercote, Peter J Stuckey, Ralph Becket, Sebastian Brand, Gregory J Duck, and Guido Tack. Minizinc: Towards a standard cp modelling language. In *International Conference on Principles and Practice of Constraint Programming*, volume 4741 of *Lecture Notes in Computer Science*, pages 529–543. Springer, 2007. doi:10.1007/978-3-540-74970-7_38.

12 Steven David Prestwich. CNF encodings. In Armin Biere, Marijn Heule, Hans van Maaren, and Toby Walsh, editors, *Handbook of Satisfiability*, pages 75–97. IOS Press, 2009.

13 H.M. Van Es and C.L. Van Es. Spatial nature of randomization and its effect on the outcome of field experiments. *Agronomy Journal*, 85(2):420–428, 1993.

14 Marc Vinyals, Jan Elffers, Jesús Giráldez-Cru, Stephan Gocht, and Jakob Nordström. In between resolution and cutting planes: A study of proof systems for pseudo-boolean SAT
48:10 Proven Optimally-Balanced Latin Rectangles with SAT

solving. In Olaf Beyersdorff and Christoph M. Wintersteiger, editors, *Theory and Applications of Satisfiability Testing - SAT 2018 - 21st International Conference, SAT 2018. Held as Part of the Federated Logic Conference, FloC 2018, Oxford, UK, July 9-12, 2018, Proceedings*, volume 10929 of *Lecture Notes in Computer Science*, pages 292–310. Springer Verlag, 2018. doi:10.1007/978-3-319-94144-8_18.

15 Nathan Wetzler, Marijn J. H. Heule, and Warren A. Hunt. DRAT-trim: Efficient checking and trimming using expressive clausal proofs. In *Theory and Applications of Satisfiability Testing - SAT 2014*, volume 8561 of *Lecture Notes in Computer Science*, pages 422–429. Springer Verlag, 2014. doi:10.1007/978-3-319-09284-3_31.
# Enumerative Level-2 Solution Counting for Quantified Boolean Formulas

Andreas Plank

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

Institute for Symbolic Artificial Intelligence, JKU Linz, Austria

Sibylle Möhle

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

Max Planck Institute for Informatics, Saarbrücken, Germany

Martina Seidl

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

Institute for Symbolic Artificial Intelligence, JKU Linz, Austria

## Abstract

We lift the problem of enumerative solution counting to quantified Boolean formulas (QBFs) at the second level. In contrast to the well-explored model counting problem for SAT (#SAT), where models are simply assignments to the Boolean variables of a formula, we are now dealing with tree (counter-)models reflecting the dependencies between the variables of the first and the second quantifier block. It turns out that enumerative counting on the second level does not give the complete model count. We present the – to the best of our knowledge – first approach of counting tree (counter-)models together with a counting tool that exploits state-of-the-art QBF technology. We provide several kinds of benchmarks for testing our implementation and illustrate in several case studies that solution counting provides valuable insights into QBF encodings.

**2012 ACM Subject Classification** Theory of computation → Automated reasoning

**Keywords and phrases** QBF, Second-Level Model Counting

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

**Category** Short Paper

**Supplementary Material** Software: https://qcount.pages.sai.jku.at/12count

**Funding** This work has been supported by the LIT AI Lab funded by the State of Upper Austria.

# 1 Introduction

Over the last decade, solvers for quantified Boolean formulas (QBFs) have become appealing tools to handle PSPACE-hard problems as found in many applications from, for example, artificial intelligence and formal verification (see [21] for a survey). Much progress has been made in the theory and practice of solving [3, 18], partly relying on generalizations of techniques from SAT, but partly being enabled by new genuine QBF techniques. However, some aspects of QBFs are hardly explored yet. One of these is counting the number of solutions of a formula. The problem of counting the number of solutions of a given QBF is also known as #QBF [13]. In contrast to #SAT [9], the counting problem of propositional logic, which is used in many application domains including probabilistic reasoning [7, 19], verification of neural networks [1, 16] and the analysis of software vulnerability [5, 24], #QBF has mainly been studied theoretically [13, 10, 2].

Only recently, some work has been presented dealing with counting the solutions of a QBF at the outer-level [22]. Given a true QBF with first quantifier block $∃X$, this work is concerned with the problem of counting the number of assignments to the variables in $X$ leading to a QBF that is true. The authors also consider the dual problem, i.e., given a false QBF starting with quantifier block $∀Y$, how many assignments of the variables $Y$ result in a false QBF? In order to answer these questions, they presented an enumerative approach that

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

© Andreas Plank, Sibylle Möhle, and Martina Seidl;
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. 49; pp. 49:1–49:10

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

LIPICS

Leibniz International Proceedings in Informatics

Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl Publishing, Germany
directly resulted from enumerative propositional model counting [6, 8]. Therefore, a solver is
used to obtain one solution, which is then excluded from the formula to obtain a different
solution. This process is repeated until no further solution is found.

In this work, we go one step further by dealing with solutions of true QBFs that start with quantifier prefix $\forall X \exists Y$ and false QBFs that start with quantifier prefix $\exists X \forall Y$, i.e., in the first case, we count models and in the second case, we count counter-models. For both cases, the solutions are not simply propositional assignments as in outer-level counting, but tree models and counter-models capturing the dependencies between the variables of the first and second quantifier block. While tree (counter-)models are a very convenient way to describe QBF solutions, in practice, QBF solvers represent solutions by Boolean functions. We will use both views on QBF solutions to investigate to what extend enumeration-based counting can be lifted to the second level. Therefore, solutions are excluded by conjunctively/disjunctively adding (negated) functions to the formula until no further solutions are found. It turns out, that enumeration-based solution counting at the second level does not lead to the full model count, but still gives valuable insights about the solutions of a formula. We implemented the approach in a first prototype to evaluate how it works in practice. To this end, we build on state-of-the-art QBF solving technology like incremental solving and function extraction.

This paper is structured as follows. First, we introduce the necessary preliminaries in the next section. Then we present the level-2 counting problem by an example in Section 3 before we discuss enumeration-based solution counting in Section 4. In Section 5, evaluate our implementation on three different types of benchmarks, before we conclude in Section 6.

2 Preliminaries

We consider quantified Boolean formulas of the form $\Pi.\phi$, where $\Pi = Q_1X_1...Q_n$ is called quantifier prefix (with $Q_i \in \{\forall, \exists\}, Q_i \neq Q_{i+1}$ for $i \in \{1, ..., n-1\}$ and $X_1, ..., X_n$ are pairwise disjoint, non-empty sets of Boolean variables). The matrix $\phi$ is a propositional formula over variables $X_i$ with standard Boolean connectives $\neg, \vee, \wedge, \to, \leftrightarrow$ and $\oplus$ (XOR). The prefix $\Pi$ induces an ordering on the variables: $x_i <_{\Pi} x_j$ if $x_i \in X_i, x_j \in X_j$ and $i < j$. If prefix $\Pi$ is clear from the context, we just write $x_i < x_j$. In this paper, we consider only closed formulas, i.e., every variable that occurs in matrix $\phi$ also occurs in the prefix $\Pi$.

A QBF $\Pi.\phi$ is in *prenex conjunctive normal form* (PCNF), if $\phi$ is a conjunction of clauses. A *clause* is a disjunction of literals and a *literal* is a variable or a negated variable. If $l$ is a literal, then $var(l) = x$ if $l = x$ or $l = \neg x$. As usual, $\bar{l} = x$ if $l = \neg x$ and $\bar{l} = \neg x$ otherwise. For a QBF $\varphi = Q_1X_1...Q_nX_n.\phi$, $var(\varphi) = X_1 \cup ... \cup X_n$. An *assignment* $\sigma$ of $\varphi$ is a set of literals over (a subset of) $var(\varphi)$ such that there is no $l \in \sigma$ with $\bar{l} \in \sigma$. For an assignment $\sigma$, $var(\sigma) = \{var(l) | l \in \sigma\}$. If $var(\sigma) = var(\varphi)$, then $\sigma$ is a full assignment, else it is a partial assignment. A (partial) X-assignment is an assignment over a (sub-)set of variables $X$.

Given a propositional formula $\phi$ and an assignment $\sigma$, then $\phi_\sigma$ denotes the formula
obtained when setting all variables $x \in var(\sigma)$ to true if $x \in \sigma$ and to false if $\neg x \in \sigma$,
respectively. Based on this notation, the semantics of a QBF is defined as follows: $\forall x \Pi.\phi$
is true iff $\Pi.\phi_{\{x\}}$ and $\Pi.\phi_{\{\neg x\}}$ are true. A QBF $\exists x \Pi.\phi$ is true iff $\Pi.\phi_{\{x\}}$
or $\Pi.\phi_{\{\neg x\}}$ is true.
For example, the QBF $\forall x \exists y.(x \leftrightarrow y)$ is true, while the QBF $\exists y \forall x.(x \leftrightarrow y)$ is false. A *model*
for a QBF $\varphi = \Pi.\phi$ with $|var(\varphi)| = m$ is a tree of height $m+1$ such that every node at level
$k \in \{1, \dots, m\}$ is labeled with a variable $x_k$ in the order of the prefix, i.e., if variable $x_j$ is at
level $j$ and variable $x_k$ is at level $k$ with $j < k$ then $x_j \le \Pi x_k$. A node at level $k$ has one
child if $x_k \in X_i$ and $Q_i = \exists$, and two children if $Q_i = \forall$. In the graphical representation, a
dashed edge indicates that the parent variable is set to false and a solid edge indicates that
it is set to true. Examples are shown in Figure 1. The path from the root to the leaves gives
a full assignment under which $\phi$ evaluates to true.
<figure><img src="image_477.png" /><figcaption>Figure 1 The left tree shows the full assignment tree for the given QBF. Tree models (16 in total) are subtrees of the assignment tree such that universal nodes have two children, existential nodes have one child and all leaves are $\top$. An example of a tree model is shown on the right.</figcaption></figure>

A counter-model is defined dually, with the difference that universal nodes have one child node, existential nodes have two child nodes, and the matrix evaluates to false under the respective assignments. An alternative representation of tree models that is practically more relevant are sets of Skolem functions, so-called Skolem sets. A Skolem set $F$ of a true QBF $\Phi$ is a set of Boolean functions such that for each existential variable $y$ of $\Phi$, there is a function $f_y(x_1, \dots, x_n) \in F$ where $x_1, \dots, x_n$ are the universal variables of $\Phi$ with $x_i < y$. The elements of a Skolem set are called Skolem functions. If we take the matrix of $\Phi$ and replace all the existential variables by their Skolem functions, we obtain a valid propositional formula over the universal variables of $\Phi$. A Skolem set for the true QBF $\forall x \exists y. (x \leftrightarrow y)$ is $\{f_y = x\}$. If we replace $y$ by $f_y$, we obtain the valid formula $(x \leftrightarrow x)$. If $\sigma_X$ is an assignment to a set of universal variables $X$, then $f_y(\sigma_X)$ denotes the value of $y$ w.r.t. the assignment of the variables in $X$. For false QBFs and their tree counter-models, functions for the universal variables are defined dually. Such functions are called Herbrand functions and are collected in Herbrand sets. If we replace all universal variables by their Herbrand functions of a Herbrand set, we obtain an unsatisfiable formula.

## 3 Counting at Level Two

In the following, we introduce the QBF level-2 counting problem by an example. Consider the true QBF

$$ \Phi = \forall x, y \exists a, b. (\neg x \lor a) \land (\neg y \lor b). $$

On the left of Figure 1 we see the full assignment tree of $\varphi = \Pi.\phi$. Each complete path from the root to a leaf node represents a full assignment to the variables in $\varphi$. A tree model of $\varphi$ is a subtree of this assignment tree such that nodes $x$ and $y$ have two children, nodes $a$ and $b$ have one child each, and the leaves are labeled by $\top$. An example of a tree model is shown on the right of Figure 1. We are concerned with the question of counting the number of tree models of a given true QBF. Alternatively, we could also ask the question how many different Skolem sets a true QBF has. Here it is important that syntactically different, but semantically equivalent Skolem sets are not counted multiple times. For example the Skolem functions $f_a = \top$ and $f_a = (x \lor \neg x)$ of our example are semantically equivalent even though they are syntactically different.

A direct approach to count tree models is to iterate over all assignments $\sigma$ of the universal variables in the first quantifier block of a true QBF $\varphi = \Pi.\phi$. With a propositional model counter we can determine the number of satisfying assignments of $\phi_\sigma$. If $S$ is the set of all assignments of the universal variables in the first quantifier block, then the total number
49:4 Enumerative Solution Counting for QBFs

<figure><img src="image_478.png" /><figcaption>Figure 2 Tree models of $\forall x\forall y\exists a\exists b\exists c.(\neg x \lor a \lor c) \land (\neg y \lor b \lor \neg c)$.</figcaption></figure>

of tree models is $\Pi_{\sigma \in S}SAT(\phi_{\sigma})$. Counting counter-models is defined dually, with the
difference that the role of existential and universal quantifiers is exchanged, i.e., the formulas
need to have a prefix starting with $\exists X\forall Y$. While this direct approach results in a complete
solution counter, it is very inefficient to iterate over all assignments of the variables in the
outermost quantifier block and solve a propositional counting problem for each assignment.
In the following section, we therefore investigate if an enumeration-based approach relying
on recent QBF solving technology is possible.

4 Enumerative Model Counting for 2QBF

In this section, we formalize enumeration-based solution counting for the second level. For
technical simplicity, we focus on true 2QBFs (formulas with quantifier prefix $\forall X \exists Y$). To
this end, we lift the idea of blocking clauses to blocking functions. Recall that a blocking
clause in SAT is the negation of a model $\sigma$ of a formula $\phi$. If $\phi$ is enriched with $\neg \sigma$, then $\sigma$
is excluded from the solution space. For QBFs, we introduce the notion of blocking Skolem
set as follows.

**Definition 1.** Let $\Phi = \forall X \exists Y. \phi$ be a true QBF and let $F$ be a Skolem set of $\Phi$. Then $\neg \phi_F$ is a blocking Skolem set of $\Phi$ where $\phi_F = \bigwedge_{f_y \in F} (y \leftrightarrow f_y)$.

In the following example, we now use blocking Skolem sets to exclude solutions.

**Example 2.** Consider the true QBF $\Phi = \forall x \forall y \exists a \exists b \exists c. (\neg x \lor a \lor c) \land (\neg y \lor b \lor \neg c)$. We now incrementally add blocking Skolem sets $\neg \psi_{F_1}$ until the formula becomes false. The resulting models are shown in Figure 2.

1. One solution of this formula is tree $T_1$, which is represented by Skolem set $F_1 = \{f_a(x, y) = \top, f_b(x, y) = \top, f_c(x, y) = \neg y\}$. To exclude $F_1$, we add $\neg \psi_{F_1}$ to $\Phi$ and obtain $\Phi_1 = \forall x \forall y \exists a \exists b \exists c. (\neg x \lor a \lor c) \land (\neg y \lor b \lor \neg c) \land \neg \psi_{F_1}$. Now $T_1$ is no solution of $\Phi_1$.

2. A solution of $\Phi_1$ is $T_2$ with Skolem set $F_2 = \{f_a(x, y) = \bot, f_b(x, y) = (x \leftrightarrow y), f_c(x, y) = x\}$. To exclude $F_2$ as well, we get $\Phi_2 = \forall x \forall y \exists a \exists b \exists c. (\neg x \lor a \lor c) \land (\neg y \lor b \lor \neg c) \land \neg \psi_{F_1} \land \neg \psi_{F_2}$.

3. Next, we get tree model $T_3$ with Skolem set $F_3 = \{f_a(x, y) = \top, f_b(x, y) = (x \oplus y), f_c(x, y) = \neg x\}$. We exclude $F_3$ from $\Phi_2$ as before and obtain $\Phi_3$.

4. The next tree model we find is $T_4$ with Skolem set $F_4 = \{f_a(x, y) = x, f_b(x, y) = y, f_c(x, y) = y\}$.

5. Finally, the QBF $\Phi_4 = \forall x \forall y \exists a \exists b \exists c. \phi \land \neg \psi_{F_1} \land \neg \psi_{F_2} \land \neg \psi_{F_3} \land \neg \psi_{F_4}$ is false.

We found four different tree models, each with four branches. Hereby Tseitin transform-
ation is used to add the blocking function to the PCNF formula before each solver call.
From these four models, we can assemble more models by combining the red branches of
Figure 2. To calculate the full model count, all possible combinations of the paths need to be considered, resulting in $4^4$ models in total (there are four choices for the branch {$\bar{x}, \bar{y}$}, four choices for the branch {$x, \bar{y}$} and so on). When we take a closer look, however, the count is not correct, as also {$\bar{x}, \bar{y}, a, b, \bar{c}$} is a model of the matrix of $\Phi$, but it has not been considered.

As the example above shows, blocking Skolem sets can indeed be used to exclude solutions, but they are sometimes too restrictive. To describe what is happening when adding blocking Skolem sets to a formula, we need to introduce the following definitions.

**Definition 3.** Let $\Phi = \forall X \exists Y. \phi$ be a true QBF and let $T_1$ and $T_2$ be tree models of $\Phi$. Then $T_1$ and $T_2$ are disjoint if there are no complete paths $\sigma_1$ of $T_1$ and $\sigma_2$ of $T_2$ with $\sigma_1 = \sigma_2$.

A complete path describes a full assignment (i.e., an assignment of all variables), such that the according branch of the tree model evaluates to true. The trees of Figure 2 are disjoint, because all their paths are different. If we want to characterize the models of a QBF in terms of Skolem sets, we get the following definition.

**Definition 4.** Let $\Phi = \forall X \exists Y. \phi$ be a true QBF and let $F$ and $G$ be Skolem sets of $\Phi$. Then $F$ and $G$ are called disjoint if for each full assignment $\sigma_X$ of the universal variables $X$ there exists an existential variable $y \in Y$ such that $f_y(\sigma_X) \neq g_y(\sigma_X)$ with $f_y \in F$ and $g_y \in G$.

The notion of disjoint models is rather strong as it forbids to have to have common paths in the tree model representation. We therefore also introduce the notion of different models requiring that at least one path of two tree models is different. This notion is transferred to Skolem sets as follows.

**Definition 5.** Let $\Phi = \forall X \exists Y. \phi$ be a true QBF. Furthermore, let $F$ and $G$ be Skolem sets of $\Phi$. Then $F$ and $G$ are called different if there exists a full assignment $\sigma_X$ of the universal variables $X$ such that there is an existential variable $y \in Y$ with $f_y(\sigma_X) \neq g_y(\sigma_X)$, $f_y \in F$, and $g_y \in G$.

Obviously, any two disjoint models are different. The other direction does not hold. To count all models of a QBF, we need to count the different models. As the following lemma shows, using blocking Skolem sets excludes disjoint models only.

**Lemma 6.** Let $\Phi = \forall X \exists Y. \phi$ be a QBF and $F$ be a Skolem set of $\Phi$. If $G$ is a Skolem set of $\Phi' = \forall X \exists Y. \phi \wedge \neg \phi_F$, where $\neg \phi_F$ is a blocking Skolem set as defined above, then $F$ and $G$ are disjoint.

**Proof.** Assume that $F$ and $G$ are not disjoint. Then there is an assignment $\sigma_X$ such that for all $y \in Y$ it holds that $f_y(\sigma_X) = g_y(\sigma_X)$ for $f_y \in F, g_y \in G$. Now we extend the assignment $\sigma_X$ to an assignment over $X \cup Y$ as follows: $\sigma = \sigma_X \cup \{y | f_y(\sigma_X) = \top, f_y \in F\} \cup \{\neg y | f_y(\sigma_X) = \bot, f_y \in F\} = \sigma_X \cup \{y | g_y(\sigma_X) = \top, g_y \in G\} \cup \{\neg y | g_y(\sigma_X) = \bot, g_y \in G\}$. As $G$ is a Skolem set of $\Phi'$, $\sigma$ has to satisfy $\neg \phi_F$. But by its construction, $\sigma$ also satisfies $\phi_F$, leading to a contradiction. Hence, $F$ and $G$ have to be disjoint. ◀

As the following lemma shows, we can use the enumerative approach to calculate the maximum number of pairwise disjoint Skolem sets.

**Proposition 7.** Let $\Phi = \forall X \exists Y. \phi$ be a true QBF and let $F_1, \dots, F_m$ be pairwise disjoint Skolem sets of $\Phi$ such that

$$ \Phi' = \forall X \exists Y. \phi' = \forall X \exists Y. (\phi \wedge \neg \phi_{F_1} \wedge \dots \wedge \neg \phi_{F_m}) $$

is false. Then $m$ is the maximum number of pairwise disjoint Skolem sets.
**Proof.** Since $\Phi'$ is false, there is at least one assignment $\sigma_X$ such that $\phi'$ is unsatisfiable under $\sigma_X$. Assume there is a Skolem set $F_{m+1}$ that is pairwise disjoint with $F_1, \dots, F_m$. Let $\sigma = \sigma_X \cup \{x \mid f_x(\sigma_X) = \top, f_x \in F\} \cup \{\neg x \mid f_x(\sigma_X) = \bot, f_x \in F\}$. Since $F_{m+1}$ is a Skolem set of $\Phi$, $\phi$ is satisfied by $\sigma$. Further, $\sigma$ satisfies $\neg\phi_{F_i}$ with $1 \le i \le m$, because $F_i$ and $F_{m+1}$ are disjoint. Hence, $\phi'$ is satisfied by $\sigma$. This contradicts the assumption that $\phi'$ is unsatisfiable under $\sigma_X$. ◀

The maximum number of pairwise disjoint models is determined by the assignments of the variables in the outermost universal quantifier block that, when the universal variables are assigned accordingly, lead to the fewest propositional models.

**Lemma 8.** Let $\Phi = \forall X \exists Y. \phi$ be a true QBF and let $S$ be the set of all full assignments of universal variables $X$. Then the maximum number of pairwise disjoint models of $\Phi$ is $\min(\{\#SAT(\phi_\sigma) \mid \sigma \in S\})$.

The proof of the lemma above follows directly from the construction of disjoint models. For example, the QBF $\forall x, y \exists a, b. (\neg x \lor a) \land (\neg y \lor b)$ with the assignment tree shown in Figure 1 has only one disjoint model, because for the assignment $\sigma = \{x, y\}$, there is only one assignment to $\{a, b\}$ that satisfies the matrix under $\sigma$. In contrast, the QBF $\forall x \forall y \exists a \exists b \exists c. (\neg x \lor a \lor c) \land (\neg y \lor b \lor \neg c)$ from Example 2 has four disjoint models, because of the assignment $\{x, y\}$. In practical encodings, those assignments of the universal variables determine the number of disjoint models for which the assignment of the outermost variables not immediately falsifies the formula. The number of disjoint models gives us a lower bound for the full model count.

**Proposition 9.** Given a true QBF $\Phi = \forall X \exists Y. \phi$ with $|X| = n$ and $m$ pairwise disjoint tree models. Then $\Phi$ has at least $m^{2^n}$ different tree models.

If we know that a true QBF $\Phi = \forall X \exists Y. \phi$ with $|X| = n$ has $m$ pairwise disjoint models, we can calculate the full model count as follows. Let $\Phi' = \forall X \exists Y. \phi'$ with $\phi' = \phi \land \neg \phi_{F_1} \land \dots \land \neg \phi_{F_m}$ be the false QBF obtained by enriching $\Phi$ with the $m$ blocking Skolem sets. Furthermore, let $S$ be the set of assignments of variables $X$ such that $\phi'_\sigma$ is true $\forall \sigma \in S$. Then the full model count is

$$ \prod_{\sigma \in S} \#SAT(\phi'_{\sigma}) * m^{2^n - |S|} $$

While the enumerative approach also works for true QBFs with prefix $\forall X \exists Y \forall Z \Pi. \phi$ for counting disjoint partial Skolem sets that consider only functions of the set $Y$, it is not possible to obtain a full model count of partial Skolem sets in the $Y$ variables by using a propositional model counter. Here, the approach needs to be applied recursively.

The enumerative approach for counting disjoint models also directly transfers to an enumerative approach for counting disjoint counter-models of false QBFs starting with the prefix $\exists X \forall Y \exists Z$.<fnref n="1" />

**Proposition 10.** Let $\Phi = \exists X \forall Y \exists Z. \phi$ be a false QBF and let $H_1, \dots, H_m$ be disjoint Herbrand sets such that

$$ \Phi' = \forall X \exists Y. (\phi \lor \phi_{H_1} \lor \dots \lor \phi_{H_m}) $$

is true with $\phi_H = \bigwedge_{h_y \in H} (y \leftrightarrow h_y)$. Then $m$ is the maximum number of pairwise disjoint Herbrand sets over the variables from $Y$.

<footnote n="1">¹ Note that we include the last quantifier block to deal with formulas in PCNF. Otherwise, the universal variables could be immediately removed.</footnote>
<figure><img src="image_479.png" /><figcaption>Figure 3 Runtime related to number of disjoint models for randomly generated formulas (left) and Unique-SAT/QBFEval benchmarks (right).</figcaption></figure>

# 5 Evaluation

We implemented the previously described approach in the tool qCounter.<fnref n="2" /> As backend solver our tool relies on the QBF solver DepQBF 6.03 [15] which is able to produce Q-resolution proofs for true and false formulas from which Skolem and Herbrand functions can be extracted. We used the QBF certification framework QBFCert [17] to obtain the functions in the Aiger format.³ We implemented a small tool to convert the blocking functions of the variables from the second quantifier block. Therefore, the Skolem functions needed to be negated and appended conjunctively to the current QBF by using the incremental interface of DepQBF. To disjunctively add the Herbrand functions, an extra Tseitin transformation step is necessary to obtain a formula in PCNF. While this transformation introduces auxiliary variables we can avoid an exponential increase of the formula using this technique[23, 12]. For this we also used the incremental interface of DepQBF with its ability to add temporary clauses. To calculate the full model count as described above, we employed the propositional model counter Ganak [20] together with the SAT solver Lingeling [4].

Currently, there exist no standard benchmark sets which we could use to evaluate the correctness and performance of our implementation. To this end, we propose three different benchmark sets: (1) randomly generated formulas, (2) QBF encodings of the unique-SAT problem (does a given propositional formula have exactly one solution?) and (3) benchmarks from past QBF Evaluations with a suitable quantifier structure. Whereas the benchmarks from past QBF Evaluations with a suitable quantifier structure. Whereas the benchmarks from (1) and (2) are constructed in such a way that the number of models is known, this is not the case for the benchmarks from (3). Details follow below. All experiments were performed on a cluster of dual-socket AMD EPYC 7313 @ 16 × 3.7GHz machines with 4GB memory limit and 1800 seconds as timeout.

<footnote n="2">² The tool, the benchmarks and the log-files are publicly available at https://qcount.pages.sai.jku.at/12count
³ http://fmv.jku.at/aiger/</footnote>
## 5.1 Randomly Generated Formulas

We provide a generator that returns true 2-QBFs in PCNF with quantifier structure $\forall X \exists Y$ such that the number of (disjoint) models is known. The parameters $n$ (size of $X$), $m$ (size of $Y$) have to be provided. For producing the clauses of the matrix, the generator iterates over all possible assignments of the universal variables $X$ and excludes propositional models by randomly assigning the existential variables $Y$ and building blocking clauses. For example, for a prefix $\forall x_1, x_2 \exists y_1, y_2$ the clause $(x_1 \lor \bar{x}_2 \lor y_1 \lor y_2)$ prohibits in any QBF model that in the branch where $x_1$ is set to false and $x_2$ is set to true, both $y_1$ and $y_2$ are both set to false. The number of disjoint models is determined by the branch with the smallest number of propositional models (see also Lemma 8).

We generated 3950 formulas with $2 \le |X|, |Y| \le 11$ having up to 100 disjoint models each. For 2608 of these formulas, the number of disjoint models could be determined. The results are shown in the left plot of Figure 3. We cluster the results according to the size of $X$, indicating that not only the number of disjoint models impacts the runtime, but also the size of the first quantifier block.

## 5.2 Unique-SAT Encodings

The question whether a given propositional formula $\phi$ has exactly one model can be encoded by a QBF $\exists X \forall \psi$ as shown in [11]. In this QBF, $X$ contains the variables of $\phi$, $Y$ is a copy of the variables of $\phi$, and $\psi$ is constructed in such a way that the QBF is true iff $\phi$ has exactly one model. The prefix fits for counting the counter-models at level-2 in case the QBF is false. The number of QBF counter-models can be determined based on the number of the models of $\phi$. If $\phi$ is unsatisfiable and has $n$ variables, then the number of disjoint counter-models equals the number of different counter-models and is $2^{2^n}$. If $\phi$ has $m$ models, then it has $(m-1)$ disjoint models and $(2^n)^{2^n-m}(m-1)^m$ different models.

We generated 497 Unique-SAT formulas based on true propositional formulas with 26 to 240 variables and 61 and 643 clauses by using the random generator from [14]. We considered only propositional formulas with more than one model in order to count disjoint counter-models of the resulting false QBFs. The performance of our tool is shown in the right plot of Figure 3. The number of disjoint counter-models enumerated by our tool could be quickly validated with a propositional model counter. In this way, we could check the results of our tool for false formulas.

## 5.3 Benchmarks from QBF Evaluations

As a final set of benchmarks, we considered formulas from the main tracks of QBFEval 2022 and QBFEval 2008.<fnref n="4" /> From the QBFEval 2022 set, we identified 18 true formulas and 102 false formulas with a suitable prefix structure that could be solved by plain DepQBF within a time limit of 1800 seconds. In a similar way, we selected two true formulas and 683 false formulas from the QBFEval 2008 set. For the 2022 formulas, we could determine the number of disjoint models of 3 formulas and the number of disjoint counter-models of 53 formulas. For the 2008 formulas, we could determine the number of disjoint models for both formulas and the number of disjoint counter-models of 285 formulas. Interestingly, many of the formulas have only one disjoint model indicating that the search space is strongly

<footnote n="4">⁴ http://www.qbflib.org</footnote>
restricted for certain assignments of the variables in the outermost quantifier block. On
the other hand, for some of the formulas we could find up to 82 disjoint (counter-)models.
Details are shown in the right plot of Figure 3.

# 6 Conclusion

We considered the problem of counting level-2 (counter-)models for QBFs. It turned out that enumerating Skolem/Herbrand functions does not provide the full count of solutions as it is the case for propositional model counting and counting QBF models at the outer level. We characterized the subset of solutions which can be counted in an enumerative manner. This subset often connects to interesting solutions of a QBF encoding. We also provided the first practical implementation of an enumerative method for level-2 solution counting. Our approach cannot only be used for counting, but also for explicitly enumerating solutions which might also have some applications in better understanding and debugging QBF encodings. We provided several sets of benchmarks of true and false instances to evaluate our implementation.

We consider this work as an important first step to full solution counting, i.e., for QBFs
with an arbitrary number of quantifier blocks. While a recursive application of our approach
seems possible, it has to be expected that this approach will be inefficient. Hence, alternative
ways need to be explored like a tight integration of solving and counting or exploiting
approaches that process the prefix in a reverse order as it is for example done in certain
preprocessing techniques.

References

1 Teodora Baluta, Shiqi Shen, Shweta Shinde, Kuldeep S. Meel, and Prateek Saxena. Quantitative verification of neural networks and its security applications. In *Proc. of the 2019 ACM SIGSAC Conf. on Computer and Communications Security*, pages 1249–1264. ACM, 2019.

2 Michael Bauland, Elmar Böhler, Nadia Creignou, Steffen Reith, Henning Schnoor, and Heribert Vollmer. Quantified constraints: The complexity of decision and counting for bounded alternation. *Electron. Colloquium Comput. Complex.*, TR05-024, 2005.

3 Olaf Beyersdorff, Janota Mikolás, Florian Lonsing, and Martina Seidl. Quantified Boolean Formulas. In *Handbook of Satisfiability*, volume 336, pages 1177–1221. IOS Press, 2021.

4 Armin Biere. CaDiCaL, Lingeling, Plingeling, Treengeling and YalSAT Entering the SAT Competition 2018. In *Proc. of SAT Competition 2018 - Solver and Benchmark Descriptions*, volume B-2018-1 of *Department of Computer Science Series of Publications B*, pages 13–14. University of Helsinki, 2018.

5 Fabrizio Biondi, Michael A. Enescu, Annelie Heuser, Axel Legay, Kuldeep S. Meel, and Jean Quilbeuf. Scalable approximation of quantitative information flow in programs. In *Proc. of Int. Conf. on Verification, Model Checking, and Abstract Interpretation*, volume 10747 of *LNCS*, pages 71–93. Springer, 2018.

6 Elazar Birnbaum and Eliezer L. Lozinskii. The good old Davis-Putnam procedure helps counting models. *J. Artif. Intell. Res.*, 10:457–477, 1999.

7 Supratik Chakraborty, Kuldeep S. Meel, and Moshe Y. Vardi. Algorithmic improvements in approximate counting for probabilistic inference: From linear to logarithmic SAT calls. In *Proc. of Int. Joint Conf. on Artificial Intelligence*, pages 3569–3576. IJCAI/AAAI Press, 2016.

8 Olivier Dubois. Counting the number of solutions for instances of satisfiability. *Theor. Comput. Sci.*, 81(1):49–64, 1991.

9 Carla P. Gomes, Ashish Sabharwal, and Bart Selman. Model counting. In *Handbook of Satisfiability*, pages 993–1014. IOS Press, 2021.
49:10 Enumerative Solution Counting for QBFs

10 Lane A. Hemaspaandra and Heribert Vollmer. The satanic notations: counting classes beyond #P and other definitional adventures. *SIGACT News*, 26(1):2–13, 1995.

11 Hans Kleine Büning and Theodor Lettmann. *Propositional logic - deduction and algorithms*. Cambridge University Press, 1999.

12 Elias Kuiter, Sebastian Krieter, Chico Sundermann, Thomas Thüm, and Gunter Saake. Tseitin or Not Tseitin? The Impact of CNF Transformations on Feature-Model Analyses. In *Proc. of the 37th IEEE/ACM Int. Conf. on Automated Software Engineering*. ACM, 2023.

13 Richard E. Ladner. Polynomial space counting problems. *SIAM J. Comput.*, 18(6):1087–1097, 1989.

14 Massimo Lauria, Jan Elffers, Jakob Nordström, and Marc Vinyals. CNFgen: A Generator of Crafted Benchmarks. In *Proc. of the 20th Int. Conf. on Theory and Applications of Satisfiability Testing*, volume 10491 of LNCS, pages 464–473. Springer, 2017.

15 Florian Lonsing and Uwe Egly. DepQBF 6.0: A search-based QBF solver beyond traditional QCDCL. In *Proc. of the 26th Conf. on Automated Deduction*, volume 10395 of LNCS, pages 371–384. Springer, 2017.

16 Nina Narodytska, Aditya A. Shrotri, Kuldeep S. Meel, Alexey Ignatiev, and João Marques-Silva. Assessing heuristic machine learning explanations with model counting. In *Proc. of the Int. Conf. on Theory and Applications of Satisfiability Testing*, volume 11628 of LNCS, pages 267–278. Springer, 2019.

17 Aina Niemetz, Mathias Preiner, Martina Seidl, and Armin Biere. Resolution-based certificate extraction for QBF - (tool presentation). In *Proc. of the 15th Int. Conference on Theory and Applications of Satisfiability Testing*, volume 7317 of LNCS, pages 430–435. Springer, 2012.

18 Luca Pulina and Martina Seidl. The 2016 and 2017 QBF solvers evaluations (QBFEVAL’16 and QBFEVAL’17). *Artif. Intell.*, 274:224–248, 2019.

19 Tian Sang, Paul Beame, and Henry A. Kautz. Performing Bayesian inference by weighted model counting. In *Proc. of the 20th Nat. Conf. on Artificial Intelligence*, pages 475–482. AAAI Press / The MIT Press, 2005.

20 Shubham Sharma, Subhajit Roy, Mate Soos, and Kuldeep S. Meel. Ganak: A scalable probabilistic exact model counter. In *Proc. of Int. Joint Conf. on Artificial Intelligence*, pages 1169–1176. Int. Joint Conf. on Artificial Intelligence Organization, 2019.

21 Ankit Shukla, Armin Biere, Luca Pulina, and Martina Seidl. A survey on applications of quantified Boolean formulas. In *Proc. of the Int. Conf. on Tools with Artificial Intelligence*, pages 78–84. IEEE, 2019.

22 Ankit Shukla, Sibylle Möhle, Manuel Kauers, and Martina Seidl. Outercount: A first-level solution-counter for quantified boolean formulas. In *Proc. of the 15th Int. Conf on Intelligent Computer Mathematics*, volume 13467 of LNCS, pages 272–284. Springer, 2022.

23 G. S. Tseitin. *On the Complexity of Derivation in Propositional Calculus*, pages 466–483. Springer, 1983.

24 Ziqiao Zhou, Zhiyun Qian, Michael K. Reiter, and Yinqian Zhang. Static evaluation of noninterference using approximate model counting. In *Proc. of IEEE Symposium on Security and Privacy*, pages 514–528. IEEE Computer Society, 2018.