$$ \text{sol3: } U = \begin{pmatrix} -1 & 0 & 0 & 0 & -1 & -1 \\ 0 & 0 & 0 & -1 & 0 & 0 \\ 1 & 0 & -1 & 0 & 0 & 0 \\ 0 & -1 & -1 & 1 & 0 & -1 \\ 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}, \quad V = \begin{pmatrix} 1 & 1 & -1 & 0 & 0 & -1 \\ 0 & 0 & 0 & 0 & -1 & 0 \\ 0 & -1 & 0 & -1 & 0 & 0 \\ 0 & 0 & 0 & -1 & 1 & -1 \\ 0 & 0 & 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 0 & 1 & 1 \end{pmatrix}, \quad W = \begin{pmatrix} 0 & 1 & 0 & 1 & 0 & 1 \\ 0 & 0 & 0 & 0 & 1 & 1 \\ 0 & 1 & 1 & 0 & 0 & 0 \\ 1 & 0 & -1 & 0 & 1 & 0 \\ 1 & 0 & -1 & 0 & 1 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 \end{pmatrix} $$

Figure 3 sol3 is the solution derived from enforcing the sign symmetry constraints on sol2.

Similarly, the same type of sign symmetry-breaking constraints can be applied to the W factor matrix as follows:

$$ w_{1,r} \le 0 \\ w_{k,r} \le \sum_{k'=1}^{k-1} |w_{k',r}| \quad \forall r \in \mathcal{R}, k > 1, k \in \mathcal{W}. $$

The interpretation is as follows. In Figure 1c, consider $c_4 = m_1 - m_2 + m_3 + m_6$, and notice that one can redefine $m_6 = (a_3 - a_1)(b_1 + b_2)$ to become $m_6 = (a_3 - a_1)(-b_1 - b_2)$ and then rewrite $c_4$ as $c_4 = m_1 - m_2 + m_3 - m_6$. Recall that the coefficients of the m terms in an output entry $c_k$ are the entries of row k of factor matrix W. The transformation we just performed produces two equivalent solutions and is an instance of “value symmetry” that is broken by the above constraint set as it forces the first non-zero entry of a column of W to be -1.

## 4.2 Valid Inequalities

Based on the structure of this problem, we can also introduce a series of valid inequalities that could potentially help a CP solver with propagation.

First, for the W matrix, we know that each multiplicative term $m_r$ must be used at least once for sufficiently small R (i.e., for non-trivial cases of the FMM problem where $R \le NMP$). This means that the sum of each column in W must be at least one:

$$ \sum_{k \in \mathcal{W}} |w_{k,r}| \ge 1, \quad \forall r \in \mathcal{R}. $$

Each result term $c_l$ must use at least M terms. This is due to a basic fact in algebraic complexity theory which states that the dot-product of two vectors of size M requires at least M multiplications [17]. This means that the sum of each row of W must be greater or equal to M:

$$ \sum_{r \in \mathcal{R}} |w_{k,r}| \ge M, \quad \forall k \in \mathcal{W}. $$

Each result term $c_l$ must differ in at least two $m_r$ terms; a simple proof by contradiction is omitted for brevity. This can be modelled as follows:

$$ \sum_{r \in \mathcal{R}} |w_{k,r} - w_{k',r}| \ge 2, \quad \forall k \ne k' \in \mathcal{W}. $$

Each term in the A and B matrices must appear in at least one of the multiplicative terms $m_r$. This translates to each row of U and V having at least one non-zero term as shown in the constraints below:

$$ \sum_{r \in \mathcal{R}} |u_{i,r}| \ge 1, \quad \forall i \in \mathcal{U} $$

$$ \sum_{r \in \mathcal{R}} |v_{j,r}| \ge 1, \quad \forall j \in \mathcal{V}. $$
Furthermore, each valid product of two terms from the A and B matrices, e.g., $a_2b_3$ for 2 × 2 matrices, must appear in at least one of the $\mathbb{R}$ multiplication terms. For $a_2b_3$ appears in $c_1$ and $c_2$, see Figure 1a. This can be modelled as follows:

$$ \sum_{r \in \mathcal{R}} |u_{i,r} \cdot v_{j,r}| \geq 1, \quad \forall \text{ valid } i, j. $$

## 4.3 Full CP Model

Finally, the full CP model is presented in Figure 4. The constraints in Equation (2) ensure that the output matches the original multiplication tensor and thus the validity of an assignment as a matrix multiplication algorithm. We enforce permutation symmetry-breaking with Equation (3) and sign symmetry-breaking with Equations (4)–(7). The valid inequalities are modelled through Equations (8)–(13).

$$ \sum_{r \in \mathcal{R}} (u_{i,r} \cdot v_{j,r} \cdot w_{k,r}) = T_{i,j,k}, \quad \forall i \in \mathcal{U}, j \in \mathcal{V}, k \in \mathcal{W} \quad (2) $$

$$ \text{lexicographic-strict}([u_{:,r}; v_{:,r}], [u_{:,r+1}; v_{:,r+1}]), \quad \forall r \in \mathcal{R} \quad (3) $$

$$ u_{1,r} \le 0 \qquad \forall r \in \mathcal{R} \quad (4) $$

$$ u_{i,r} \le \sum_{i'=1}^{i-1} |u_{i',r}|, \quad \forall r \in \mathcal{R}, i > 1, i \in \mathcal{U} \quad (5) $$

$$ w_{1,r} \le 0 \qquad \forall r \in \mathcal{R} \quad (6) $$

$$ w_{k,r} \le \sum_{k'=1}^{k-1} |w_{k',r}| \quad \forall r \in \mathcal{R}, k > 1, k \in \mathcal{W} \quad (7) $$

$$ \sum_{k \in \mathcal{W}} |w_{k,r}| \ge 1, \quad \forall r \in \mathcal{R} \quad (8) $$

$$ \sum_{r \in \mathcal{R}} |w_{k,r}| \ge M, \quad \forall k \in \mathcal{W} \quad (9) $$

$$ \sum_{r \in \mathcal{R}} |w_{k,r} - w_{k',r}| \ge 2, \quad \forall k \ne k' \in \mathcal{W} \quad (10) $$

$$ \sum_{r \in \mathcal{R}} |u_{i,r}| \ge 1, \quad \forall i \in \mathcal{U} \quad (11) $$

$$ \sum_{r \in \mathcal{R}} |v_{j,r}| \ge 1, \quad \forall j \in \mathcal{V} \quad (12) $$

$$ \sum_{r \in \mathcal{R}} |u_{i,r} \cdot v_{j,r}| \ge 1, \quad \forall \text{ valid } i, j \quad (13) $$

Figure 4 Full CP Model with symmetry-breaking constraints and valid inequalities.

## 4.4 Sparsity-based Problem Decomposition

Given that the factor matrices that have been found for known decompositions tend to be sparse, we introduce some inexact inequalities to induce sparsity and trim candidate assignments that have a high likelihood to be infeasible or that are unnecessarily dense. For example, observe that Strassen's solution in Figure 1c leads to many zeros in the factor
14:10 Fast Matrix Multiplication Without Tears: A Constraint Programming Approach

matrices; no m term uses more than 2 out of 4 of the a or b terms, no c term uses more than 4 out of the 7 m terms. It has been observed that as the matrix sizes grow, the best solutions become even sparser.

We first introduce a constraint limiting the number of active (i.e., nonzero) terms in each column r (i.e., multiplication term) of U and V. This constraint is written as:

$$ \sum_{i \in \mathcal{U}} |u_{i,r}| + \sum_{j \in \mathcal{V}} |v_{j,r}| \leq K_1, \quad \forall r \in \mathcal{R}. $$

A similar constraint can be imposed on W, by restricting that each output must use at most $K_2$ multiplication terms. This constraint is written as:

$$ \sum_{r \in \mathcal{R}} |w_{k,r}| \leq K_2, \quad \forall k \in \mathcal{W}. $$

Based on these constraints, $K_1$ has an upper bound of $(NM+MP)$ and $K_2$ is upper bounded by $R$. By observing decompositions for small to medium-scale matrices, we can estimate $K_1$ and $K_2$. For example, for $3 \times 3$ matrices with $R = 23$, we observe that $K_1 = 9$ and $K_2 = 10$ is the safest estimate possible compared to the upper bounds of 18 and 23, respectively, which could restrict the CP search dramatically. Note that one could start with any such estimates of the decomposition parameters $K_1$ and $K_2$, iteratively increasing them if the restricted instances are found to be infeasible by the CP solver, eventually resulting in a complete resolution of the original problem.

## 4.5 Cyclic Invariant Formulation

In contrast to the symmetries of the factor matrices discussed in Section 4.1, there exists well-known cyclic symmetry for the multiplication tensors $T_N$ of square matrices. More precisely, it is known that $T_{i,j,k} = T_{j,k,i} = T_{k,i,j}$. The authors in [1] proposed to leverage this cyclic symmetry property and parameterize FMM algorithms with cyclic invariant factor matrices: $U = [ABCD]$, $V = [ADBC]$, $W = [ACDB]$ with $A \in \{-1, 0, 1\}^{N^2 \times S}$ and $B, C, D \in \{-1, 0, 1\}^{N^2 \times T}$ corresponding to a rank $R = S + 3T$.

Although this parametrization reduces the number of integer variables by a factor of three, helping with the combinatorial nature of the problem, there is no guarantee that the minimal rank decomposition corresponds to solutions that exhibit cyclic symmetry. That being said, Strassen's solution of $R = 7$ for $N = 2$, which is optimal, exhibits such a symmetry ($S \in \{1, 4\}$), as does the best-known rank of 23 for $N = 3$ ($S \in \{2, 5, 11\}$). Performing two steps of Strassen's algorithm for $N = 2$ yields a rank 49 cyclic invariant solution for $N = 4$. It is currently unknown whether a solution of rank less than 49 exists for $N = 4$, let alone one exhibiting cyclic invariance.

We implement Ballard and Benson's cyclic invariant reduction [1] of the FMM problem for square matrices by reducing the decision variables of our CP formulation as required and imposing the invariant structure on the factor matrices.

# 5 Experiments

In this section, we present our experimental results starting by showing how our CP approach can recover the best-known upper bounds on the rank in a small amount of time on multiplication problems ranging from the trivial $(N, M, P) = (1, 1, 1)$ case all the way up to the much harder $(2, 2, 4)$ and $(3, 3, 3)$ cases. We then present results for the infeasible cases
for (2, 2, 2). We used IBM's CP Optimizer (CPO) 22.1.0⁶ to solve our CP models. We ran our experiments on a compute cluster of AMD Ryzen Threadripper 2990WX cores with 128 GB of RAM per node.

## 5.1 Experimental Setup

To ensure the reproducibility and robustness of our results, all our experiments are run with multiple random seeds. This accounts for the often observed performance variability in combinatorial search; this is documented for example in MILP [12]. To that end, we ran each experiment with 10 different seeds. We assigned 8 cores (CPO's Workers parameter) to the solver for each run (except for more compute-intensive experiments in Section 5.4 where we assigned 20 cores) and timed out the experiments after 2 hours.

## 5.2 Evaluation Metrics

We will report the solver runtime and the number of branches during the solution process for completed runs (i.e., runs that returned a feasible solution or a proof of infeasibility). Given that each problem is attempted with multiple random seeds, the shifted geometric mean with a shift of 0.0001⁷, median, minimum, and maximum of the time in seconds and the number of branches will be reported for a complete picture of the results. Runs that terminated due to the time or memory limits will be discussed where applicable.

<table><caption>Table 1 Runtime results for the base CP model on various matrix dimensions. “geo mean” refers to the shifted geometric mean as described in Section 5.2; “med” refers to the median and “min”/“max” to the minimum and maximum, respectively.</caption><thead><tr><th>N</th><th>M</th><th>P</th><th>R</th><th>Time (sec)<br>geo mean (min, med, max)</th><th>Num Branches<br>geo mean (min, med, max)</th></tr></thead><tbody><tr><td>1</td><td>1</td><td>1</td><td>1</td><td>0.00 (0.00, 0.00, 0.01)</td><td>5.05×10<sup>1</sup> (4.50×10<sup>1</sup>, 5.00×10<sup>1</sup>, 5.80×10<sup>1</sup>)</td></tr><tr><td>1</td><td>1</td><td>2</td><td>2</td><td>0.00 (0.00, 0.00, 0.01)</td><td>1.31×10<sup>2</sup> (1.00×10<sup>2</sup>, 1.27×10<sup>2</sup>, 2.13×10<sup>2</sup>)</td></tr><tr><td>1</td><td>2</td><td>1</td><td>2</td><td>0.00 (0.00, 0.01, 0.01)</td><td>1.68×10<sup>2</sup> (1.08×10<sup>2</sup>, 1.69×10<sup>2</sup>, 2.31×10<sup>2</sup>)</td></tr><tr><td>1</td><td>1</td><td>3</td><td>3</td><td>0.00 (0.00, 0.01, 0.01)</td><td>7.09×10<sup>2</sup> (3.47×10<sup>2</sup>, 7.78×10<sup>2</sup>, 1.38×10<sup>3</sup>)</td></tr><tr><td>1</td><td>3</td><td>1</td><td>3</td><td>0.00 (0.00, 0.01, 0.01)</td><td>9.44×10<sup>2</sup> (3.68×10<sup>2</sup>, 9.83×10<sup>2</sup>, 1.82×10<sup>3</sup>)</td></tr><tr><td>1</td><td>2</td><td>2</td><td>4</td><td>0.01 (0.00, 0.01, 0.01)</td><td>4.86×10<sup>3</sup> (1.68×10<sup>3</sup>, 5.27×10<sup>3</sup>, 7.68×10<sup>3</sup>)</td></tr><tr><td>2</td><td>1</td><td>2</td><td>4</td><td>0.01 (0.01, 0.01, 0.02)</td><td>4.36×10<sup>3</sup> (2.63×10<sup>3</sup>, 4.39×10<sup>3</sup>, 7.14×10<sup>3</sup>)</td></tr><tr><td>1</td><td>2</td><td>3</td><td>6</td><td>0.05 (0.02, 0.04, 0.10)</td><td>3.51×10<sup>4</sup> (1.41×10<sup>4</sup>, 2.96×10<sup>4</sup>, 9.83×10<sup>4</sup>)</td></tr><tr><td>1</td><td>3</td><td>2</td><td>6</td><td>0.05 (0.03, 0.06, 0.12)</td><td>3.19×10<sup>4</sup> (1.40×10<sup>4</sup>, 3.19×10<sup>4</sup>, 7.18×10<sup>4</sup>)</td></tr><tr><td>2</td><td>1</td><td>3</td><td>6</td><td>0.05 (0.02, 0.04, 0.11)</td><td>3.79×10<sup>4</sup> (1.45×10<sup>4</sup>, 3.33×10<sup>4</sup>, 8.93×10<sup>4</sup>)</td></tr><tr><td>2</td><td>2</td><td>2</td><td>7</td><td>0.74 (0.28, 0.75, 1.84)</td><td>6.41×10<sup>5</sup> (2.03×10<sup>5</sup>, 6.69×10<sup>5</sup>, 1.70×10<sup>6</sup>)</td></tr><tr><td>1</td><td>3</td><td>3</td><td>9</td><td>0.37 (0.26, 0.36, 0.60)</td><td>2.92×10<sup>5</sup> (1.84×10<sup>5</sup>, 3.09×10<sup>5</sup>, 4.16×10<sup>5</sup>)</td></tr><tr><td>3</td><td>1</td><td>3</td><td>9</td><td>0.42 (0.18, 0.52, 0.61)</td><td>3.10×10<sup>5</sup> (1.64×10<sup>5</sup>, 3.35×10<sup>5</sup>, 4.36×10<sup>5</sup>)</td></tr><tr><td>2</td><td>2</td><td>3</td><td>11</td><td>49.64 (0.98, 71.82, 245.06)</td><td>3.40×10<sup>7</sup> (6.90×10<sup>5</sup>, 4.71×10<sup>7</sup>, 1.74×10<sup>8</sup>)</td></tr><tr><td>2</td><td>3</td><td>2</td><td>11</td><td>26.47 (6.68, 29.41, 133.29)</td><td>1.56×10<sup>7</sup> (3.52×10<sup>6</sup>, 1.39×10<sup>7</sup>, 9.00×10<sup>7</sup>)</td></tr></tbody></table>

<footnote n="6">⁶ https://www.ibm.com/products/ilog-cplex-optimization-studio/cplex-cp-optimizer</footnote>

<footnote n="7">⁷ The shifted geometric mean of a set of n values $t_1, \dots, t_n$ is defined as $(\prod_{i=1}^n [t_i + \text{shift}])^{\frac{1}{n}} - \text{shift}$. Compared to the arithmetic mean, it is less sensitive to large variations in the values.</footnote>
## 5.3 Feasible Cases: Searching for Solutions with the Base CP Model

Table 1 shows the time and number of branches (Num Branches) required by the base CP model (i.e., without symmetry breaking or valid inequalities) to find solution for a range of problems. Our approach was able to find Strassen’s solution for the 2 × 2 matrix multiplication in less than a second whereas the AlphaTensor paper [7] reports a few minutes of model inference to find that solution.

**Performance variability.** In Table 1, we can see that for (2, 2, 3) with $R = 11$, the worst seed took 245 seconds to find a feasible solution compared to 0.98 seconds for the best seed. This drastic difference in time (and ultimately the number of branches) is an indication that minute parameters such as the seed can significantly impact the CP search. For feasible instances, this phenomenon can be seen as a blessing rather than a curse if one has access to multiple cores: the randomness can be exploited by running multiple copies of the solver, terminating as soon as the first successful run is completed. This has been done in MILP [8].

## 5.4 Feasible Cases: Sparsity Constraints and Cyclic Invariance Help

Solving the base CP formulation, with our current time and memory budgets, does not yet yield feasible decompositions for dimensions higher than (2, 2, 3) or (2, 3, 2). However, after increasing both the time and memory limits, we were able to find a solution to the problem for dimension (2, 2, 4) with $R = 14$ in 19.6 hours using the inexact inequalities ($K_1 = 11$ and $K_2 = 7$) developed in Section 4.4. Furthermore, our cyclic invariant formulation (with $S = 5$) with inexact inequalities ($K_1 = 9$ and $K_2 = 10$) was able to find a solution for (3, 3, 3), $R = 23$. More specifically, we ran the cyclic invariant formulation for 10 hours with 5 different seeds and observed that two seeds produced a feasible solution within one hour whereas the other three seeds hit the time limit. Once again, this indicates that performance variability in the CP search is significant for our problem. Additionally, we ran the base CP formulation for (3, 3, 3) without inexact inequalities for 5 seeds which all hit the time limit of 10 hours, demonstrating the benefit of the reduction of variables for the cyclic invariant formulation and sparsity constraints. We have yet to check whether using only inexact inequalities can help the base formulation for (3, 3, 3). A similar result was observed for the (2, 2, 2) case in which the cyclic invariant formulation ($S = 4$) with inexact inequalities ($K_1 = 6$ and $K_2 = 4$) produced an average solution time of 0.05 seconds across 10 seeds whereas the base CP model has an average of 1.46 seconds. Our current implementation is not able to find cyclic invariant solutions for $N = 4$ with $R = 48$, but we have hope that this approach is a promising tool for the search for new cyclic invariant solutions for square matrix multiplication.

## 5.5 Infeasible Cases: The Importance of Symmetry Breaking

Since CP performs an exhaustive search, it can provide a proof of infeasibility if a given rank $R$ is not achievable for certain matrix dimensions. As expected, the runtime to prove infeasibility significantly increases as we approach the known minimum rank; this can be seen in Table 2 for the (2,2,2) case. It is also apparent that the addition of symmetry-breaking constraints helps tremendously when proving infeasibility given that they reduce the search space significantly. More specifically, for $R = 6$ in Table 2, it is not even currently possible to prove infeasibility without symmetry-breaking constraints in 2 hours whereas the CP model with symmetry-breaking constraints (B+S) requires around 7 minutes. These results highlight the importance of symmetry-breaking constraints when looking to prove infeasibility.
**Table 2** Runtime results for the base CP model and variants to prove infeasibility of $R < 7$ for (2,2,2). “geo mean” refers to the shifted geometric mean as described in Section 5.2; “med” refers to the median and “min”/“max” to the minimum and maximum, respectively. Overall, the use of symmetry-breaking constraints (denoted by the letter “S”) on top of the base CP formulation (“B”) is crucial for efficient proofs of infeasibility. “V” refers to the valid inequalities of Section 4.2 which sometimes complement symmetry-breaking but are not always needed for the fastest results.

<table><thead><tr><th>R</th><th>Method</th><th>Time (sec)<br>geo mean (min, med, max)</th><th>Num Branches<br>geo mean (min, med, max)</th></tr></thead><tbody><tr><td rowspan="4">1</td><td>B</td><td>0.01 (0.00, 0.01, 0.02)</td><td>1.08 × 10<sup>3</sup> (1.06 × 10<sup>3</sup>, 1.08 × 10<sup>3</sup>, 1.10 × 10<sup>3</sup>)</td></tr><tr><td>B+S</td><td>0.00 (0.00, 0.01, 0.01)</td><td>1.00 × 10<sup>-4</sup> (0.00, 0.00, 0.00)</td></tr><tr><td>B+V</td><td>0.00 (0.00, 0.00, 0.01)</td><td>1.00 × 10<sup>-4</sup> (0.00, 0.00, 0.00)</td></tr><tr><td><b>B+V+S</b></td><td><b>0.00</b> (0.00, 0.00, 0.00)</td><td>1.00 × 10<sup>-4</sup> (0.00, 0.00, 0.00)</td></tr><tr><td rowspan="4">2</td><td>B</td><td>0.01 (0.00, 0.01, 0.03)</td><td>4.38 × 10<sup>3</sup> (3.72 × 10<sup>3</sup>, 4.41 × 10<sup>3</sup>, 5.30 × 10<sup>3</sup>)</td></tr><tr><td>B+S</td><td>0.01 (0.01, 0.01, 0.02)</td><td>1.08 × 10<sup>3</sup> (1.08 × 10<sup>3</sup>, 1.08 × 10<sup>3</sup>, 1.08 × 10<sup>3</sup>)</td></tr><tr><td><b>B+V</b></td><td><b>0.00</b> (0.00, 0.01, 0.02)</td><td>1.33 × 10<sup>3</sup> (1.31 × 10<sup>3</sup>, 1.32 × 10<sup>3</sup>, 1.34 × 10<sup>3</sup>)</td></tr><tr><td>B+V+S</td><td>0.01 (0.01, 0.01, 0.03)</td><td>1.09 × 10<sup>3</sup> (1.07 × 10<sup>3</sup>, 1.08 × 10<sup>3</sup>, 1.10 × 10<sup>3</sup>)</td></tr><tr><td rowspan="4">3</td><td>B</td><td>0.21 (0.17, 0.21, 0.28)</td><td>1.70 × 10<sup>5</sup> (1.42 × 10<sup>5</sup>, 1.70 × 10<sup>5</sup>, 1.95 × 10<sup>5</sup>)</td></tr><tr><td>B+S</td><td>0.02 (0.01, 0.02, 0.03)</td><td>4.13 × 10<sup>3</sup> (3.06 × 10<sup>3</sup>, 4.30 × 10<sup>3</sup>, 5.32 × 10<sup>3</sup>)</td></tr><tr><td>B+V</td><td>0.20 (0.13, 0.21, 0.25)</td><td>1.38 × 10<sup>5</sup> (1.14 × 10<sup>5</sup>, 1.35 × 10<sup>5</sup>, 1.78 × 10<sup>5</sup>)</td></tr><tr><td><b>B+V+S</b></td><td><b>0.02</b> (0.01, 0.02, 0.03)</td><td>3.61 × 10<sup>3</sup> (2.52 × 10<sup>3</sup>, 3.67 × 10<sup>3</sup>, 4.92 × 10<sup>3</sup>)</td></tr><tr><td rowspan="4">4</td><td>B</td><td>43.79 (31.33, 41.89, 66.93)</td><td>3.85 × 10<sup>7</sup> (3.06 × 10<sup>7</sup>, 3.80 × 10<sup>7</sup>, 5.00 × 10<sup>7</sup>)</td></tr><tr><td><b>B+S</b></td><td><b>0.12</b> (0.07, 0.13, 0.18)</td><td>8.50 × 10<sup>4</sup> (6.94 × 10<sup>4</sup>, 8.54 × 10<sup>4</sup>, 1.03 × 10<sup>5</sup>)</td></tr><tr><td>B+V</td><td>53.49 (39.10, 48.87, 69.35)</td><td>3.70 × 10<sup>7</sup> (3.18 × 10<sup>7</sup>, 3.69 × 10<sup>7</sup>, 4.23 × 10<sup>7</sup>)</td></tr><tr><td>B+V+S</td><td>0.15 (0.11, 0.15, 0.20)</td><td>8.53 × 10<sup>4</sup> (7.34 × 10<sup>4</sup>, 8.10 × 10<sup>4</sup>, 1.06 × 10<sup>5</sup>)</td></tr><tr><td rowspan="4">5</td><td>B</td><td>T.O. (N/A, N/A, N/A)</td><td>6.03 × 10<sup>9</sup> (5.56 × 10<sup>9</sup>, 5.75 × 10<sup>9</sup>, 7.14 × 10<sup>9</sup>)</td></tr><tr><td>B+S</td><td>3.06 (2.28, 3.02, 4.15)</td><td>2.22 × 10<sup>6</sup> (1.89 × 10<sup>6</sup>, 2.19 × 10<sup>6</sup>, 2.67 × 10<sup>6</sup>)</td></tr><tr><td>B+V</td><td>T.O. (N/A, N/A, N/A)</td><td>5.57 × 10<sup>9</sup> (3.97 × 10<sup>9</sup>, 5.83 × 10<sup>9</sup>, 6.16 × 10<sup>9</sup>)</td></tr><tr><td><b>B+V+S</b></td><td><b>2.98</b> (2.56, 2.94, 3.44)</td><td>2.14 × 10<sup>6</sup> (1.91 × 10<sup>6</sup>, 2.12 × 10<sup>6</sup>, 2.56 × 10<sup>6</sup>)</td></tr><tr><td rowspan="4">6</td><td>B</td><td>T.O. (N/A, N/A, N/A)</td><td>5.99 × 10<sup>9</sup> (4.53 × 10<sup>9</sup>, 5.79 × 10<sup>9</sup>, 6.93 × 10<sup>9</sup>)</td></tr><tr><td><b>B+S</b></td><td><b>429.26</b> (333.88, 441.63, 528.61)</td><td>3.28 × 10<sup>8</sup> (2.94 × 10<sup>8</sup>, 3.31 × 10<sup>8</sup>, 3.76 × 10<sup>8</sup>)</td></tr><tr><td>B+V</td><td>T.O. (N/A, N/A, N/A)</td><td>4.67 × 10<sup>9</sup> (3.82 × 10<sup>9</sup>, 4.73 × 10<sup>9</sup>, 5.48 × 10<sup>9</sup>)</td></tr><tr><td>B+V+S</td><td>517.33 (414.07, 522.81, 640.65)</td><td>3.35 × 10<sup>8</sup> (2.97 × 10<sup>8</sup>, 3.28 × 10<sup>8</sup>, 3.95 × 10<sup>8</sup>)</td></tr></tbody></table>

# 6 Conclusion

We have proposed a novel CP approach to solve the fast matrix multiplication problem. We have provided a set of constraints for breaking permutation and sign symmetries as well as a set of valid inequality constraints to help CP prune and propagate more efficiently. We provide a decomposition framework that is beneficial for finding feasible solutions for the largest case we have attempted, i.e., 3 × 3 matrix multiplication. Based on our experimental results, we have been able to solve small instances of this problem within a reasonable amount of time. This is in contrast to some existing search-based approaches (MILP, planning) that seem to struggle. In contrast to the AlphaTensor approach [7], the CP model is far more natural for this combinatorial task and is uniquely positioned to provide proof of infeasibility for some open problems in this space.

While the results of our approach are promising given the limited amount of computing used, there are several limitations that we aim to address in future work. First, our algorithm struggles to scale for larger matrix dimensions or ranks due to the quick increase in the number of variables of the CP model. Secondly, we have found that the base CP model outperforms the addition of symmetry constraints and valid inequalities in the case of feasible
14:14 Fast Matrix Multiplication Without Tears: A Constraint Programming Approach

solutions, likely due to the latter's tendency to prune symmetric solutions early in the tree search. However, we believe that our experiment's small matrix dimensions may have skewed these results and valid inequalities may be crucial for larger sizes. Moving forward, we propose several areas for further exploration and improvement:

- Conduct larger-scale experiments using larger compute clusters to take advantage of the parallelizability of the CP solver's search procedure.

- Analyze the highly structured nature of this problem to develop more valid inequalities that can further reduce the search space of our CP model, including inexact inequalities that may not hold for all matrix multiplication dimensions but help for some cases.

- Explore solver parameter tuning, particularly for branching strategies and other important search-related decisions.

- Further investigate the idea of sparsity-based problem decomposition as a means of improving the scalability and performance of our approach.

References

1 Austin R. Benson and Grey Ballard. A framework for practical parallel fast matrix multiplication. *ACM SIGPLAN Notices*, 50(8):42–53, January 2015. doi:10.1145/2858788.2688513.

2 Markus Blaser. On the complexity of the multiplication of matrices of small formats. *Journal of Complexity*, 19(1):43–60, 2003.

3 Markus Bläser. *Fast Matrix Multiplication*. Number 5 in Graduate Surveys. Theory of Computing Library, 2013. doi:10.4086/toc.gs.2013.005.

4 Roger W Brockett and David Dobkin. On the optimal evaluation of a set of bilinear forms. In *Proceedings of the fifth annual ACM symposium on Theory of computing*, pages 88–95, 1973.

5 Nicolas T Courtois, Gregory V Bard, and Daniel Hulme. A new general-purpose method to multiply 3x3 matrices using only 23 multiplications, 2011.

6 Hans F de Groote. On varieties of optimal algorithms for the computation of bilinear mappings ii. optimal algorithms for 2×2-matrix multiplication. *Theoretical Computer Science*, 7(2):127–148, 1978.

7 Alhussein Fawzi, Matej Balog, Aja Huang, Thomas Hubert, Bernardino Romera-Paredes, Mohammadamin Barekatain, Alexander Novikov, Francisco J R Ruiz, Julian Schrittwieser, and Grzegorz Swirszcz. Discovering faster matrix multiplication algorithms with reinforcement learning. *Nature*, 610(7930):47–53, 2022.

8 Matteo Fischetti and Michele Monaci. Exploiting erraticism in search. *Operations Research*, 62(1):114–122, 2014.

9 Ambros Gleixner, Gregor Hendel, Gerald Gamrath, Tobias Achterberg, Michael Bastubbe, Timo Berthold, Philipp Christophel, Kati Jarck, Thorsten Koch, Jeff Linderoth, et al. Miplib 2017: data-driven compilation of the 6th mixed-integer programming library. *Mathematical Programming Computation*, 13(3):443–490, 2021.

10 Marijn JH Heule, Manuel Kauers, and Martina Seidl. Local search for fast matrix multiplication. In *Theory and Applications of Satisfiability Testing-SAT 2019: 22nd International Conference, SAT 2019, Lisbon, Portugal, July 9–12, 2019, Proceedings 22*, pages 155–163. Springer, 2019.

11 Marijn JH Heule, Manuel Kauers, and Martina Seidl. New ways to multiply 3 × 3-matrices. *Journal of Symbolic Computation*, 104:899–916, 2021.

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

13 Alexey V. Smirnov. The bilinear complexity and practical algorithms for matrix multiplication. *Computational Mathematics and Mathematical Physics*, 53:1781–1795, 2013.

14 Laurent Sorber and Marc Van Barel. A mixed-integer linear program formulation for fast matrix multiplication, 2017.
15 David Speck, Paul Höft, Daniel Gnad, and Jendrik Seipp. Finding matrix multiplication algorithms with classical planning. In Sven Koenig, Roni Stern, and Mauro Vallati, editors, *Proceedings of the Thirty-Third International Conference on Automated Planning and Scheduling (ICAPS 2023)*. AAAI Press, 2023.

16 Volker Strassen. Gaussian elimination is not optimal. *Numerische mathematik*, 13(4):354–356, 1969.

17 Shmuel Winograd. On the number of multiplications necessary to compute certain functions. *Communications on Pure and Applied Mathematics*, 23(2):165–179, 1970.
<figure><img src="image_133.png" /></figure>
# Probabilistic Inference by Projected Weighted Model Counting on Horn Clauses

Alexandre Dubray ✉ ID

Institute of Information and Communication Technologies, Electronics and Applied Mathematics
(ICTEAM), UCLouvain, Belgium

Pierre Schaus ✉ ID

Institute of Information and Communication Technologies, Electronics and Applied Mathematics
(ICTEAM), UCLouvain, Belgium

Siegfried Nijssen ✉ ID

Institute of Information and Communication Technologies, Electronics and Applied Mathematics
(ICTEAM), UCLouvain, Belgium

## Abstract

Weighted model counting, that is, counting the weighted number of satisfying assignments of a propositional formula, is an important tool in probabilistic reasoning. Recently, the use of projected weighted model counting (PWMC) has been proposed as an approach to formulate and answer probabilistic queries. In this work, we propose a new simplified modeling language based on PWMC in which probabilistic inference tasks are modeled using a conjunction of Horn clauses and a particular weighting scheme for the variables. We show that the major problems of inference for Bayesian Networks, network reachability and probabilistic logic programming can be modeled in this language. Subsequently, we propose a new, relatively simple solver that is specifically optimized to solve the PWMC problem for such formulas. Our experiments show that our new solver is competitive with state-of-the-art solvers on the major problems studied.

**2012 ACM Subject Classification** Mathematics of computing → Probabilistic inference problems; Computing methodologies → Probabilistic reasoning; Mathematics of computing → Bayesian networks

**Keywords and phrases** Model Counting, Bayesian Networks, Probabilistic Networks

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

## Supplementary Material

*Software (Source Code): https://github.com/aia-uclouvain/schlandals*
*archived at swh:1:dir:38ff0f7ce1c12353938ed4056877ef497389440c*

# 1 Introduction

Weighted model counters are systems that calculate the weighted number of assignments that satisfy a formula in conjunctive normal form (CNF). Weighted model counting (WMC) is a hard problem: it is #P complete, and hence building efficient systems for this task is a challenge.

An important application of WMC is probabilistic inference, that is, calculating the probability of an observation according to a given probabilistic model. With increasing sizes of probabilistic models, the performance of inference remains important. Model counters have already been used to solve probabilistic inference problems on

- Bayesian networks (BNs) [6, 9, 20];

- Probabilistic networks (PNs) [12, 27]

- Probabilistic logic programs (PLPs) [13].

In each of these studies, approaches were proposed for how to model a probabilistic inference task as a WMC task on a CNF formula. Unfortunately, the resulting models are not always

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

© Alexandre Dubray, Pierre Schaus, and Siegfried Nijssen;
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. 15; pp. 15:1–15:17

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

LIPICS

Leibniz International Proceedings in Informatics

Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl Publishing, Germany
15:2 Probabilistic Inference by Projected Weighted Model Counting on Horn Clauses

simple. For Bayesian networks the CNF formulas are polynomial in size given the size of the BN, but a number of papers have proposed increasingly complex models to make solving efficient [2, 4, 5]; for PNs and PLPs the CNF can even be exponentially larger than the original probabilistic model [27]. A challenge remains how to model inference tasks in a simple manner and obtain good performance at the same time.

A promising solution to this modeling problem may be the use of projected (weighted) model counting (PWCM) [1] as a general approach for probabilistic inference; it was shown that using PWMC a CNF of polynomial size can be used to solve reachability problems on PNs [12]. Where in weighted model counting a sum is calculated over all models of a formula *F*, in PWMC the models are projected on a subset of the variables (the priority variables) and a weight is given only to each resulting projected model. For example, if *F* = (*a* ∨ *b*) ∧ (*b* ∨ *c* ∨ *d*), then a model counter will sum over the 11 models of *F*. On the other hand, if the priority variables are {*a*, *c*}, then a projected model counter sums over 4 projected models (as setting *b* = ⊤ always makes *F* true).

In this work, we study the challenge of efficient and simple probabilistic inference using PWMC in more detail. We show that the aforementioned probabilistic inference tasks can be modeled as PWMC tasks over a simpler form of CNF formulas: CNF formulas of *Horn* clauses with a specific weighting scheme. We will call this task the task of *projected probabilistic Horn model counting* (PPHMC). We will argue that the resulting weighted CNF models are simple and polynomial in size given the original probabilistic models.

Subsequently, we will introduce a new solver, the Schlandals solver, which takes full advantage of the probabilistic weighting scheme and the fact that all clauses in our formulas are Horn clauses. The main intuition behind our solver is that it exploits the well-known fact that the SAT problem over conjunctions of Horn clauses can be solved in linear time. Using this observation, it is able to efficiently find assignments to the non-priority variables that greatly reduce the number of constrained clauses in the input formula (like $b = \top$ in the example above). At the same time, it is still able to exploit optimizations found in other DPLL-style model counters, such as component caching, unit propagation and branching strategies, to run with a bounded use of memory.

An experimental evaluation of our solver shows that it is competitive with state-of-the-art existing tools. We compare our solver with the best performing solvers of the 2022 projected (weighted) model counting competition on two probabilistic inference tasks: inference in Bayesian Networks and reliability estimation in PNs. Our experiments show that our solver is able to solve most of the instances for the BN task, beating other tools when using simple forms of CNF models, and getting similar performance when compared to optimized complex models. On PNs, our solver even outperforms the state-of-the-art, without any optimization of our model. These results open up future possibilities for the use of PPHMC.

2 Related Work

This work is concerned with the task of *Projected Weighted Model Counting* (PWMC). In the classical *Model Counting* problem, one is interested in finding the number of models of a Boolean formula *F*, the assignments that makes the formula true. In this work we focus on the most common setting, in which *F* is a formula in conjunctive normal form (CNF). In *Projected Model Counting*, the goal is to count the assignments to a subset of the variables (the *priority variables*) such that there exists an assignment to the other variables that makes *F* true. If the set of priority variables contains all variables of *F*, the projected and unprojected problems are the same. In the rest of this work, we thus assume that the subset of priority variable is not empty. In the weighted version of these problems, each valid assignment is weighted and the goal is to return the sum of the weights of the models.
More formally, let $F$ be a Boolean formula over a set of variables $V$ with $\mathcal{P} \subseteq V$ the priority variables and $\mathcal{D} = V \setminus \mathcal{P}$ the non-priority variables. In traditional projected weighted model counting each priority variable is given two weights, $w(v)$ and $w(\neg v)$, one for each polarity. We denote by $S_X$ an assignment to the variables in $X \subseteq V$ and $S_X[x]$ the assigned value of $x \in X$. If $X, Y \subset V$ are two disjoint sets of variables, then $S_X \cup S_Y$ is defined as the conjunction of the assignments. We define

$$S_F = \{S_P \mid \exists S_D : F[S_P \cup S_D] = \top\}$$

as the set of assignments to the priority variables that can be extended with an assignment to the non-priority variables and satisfy the formula. The goal of the traditional PWMC task is then to find

$$\sum_{S \in S_F} \left( \prod_{v \in V | S[v] = \top} w(v) \times \prod_{v \in V | S[v] = \bot} w(\neg v) \right).$$

Note that when weighted model counters are used to model probabilistic inference tasks, the weights of assignments are set such that this sum corresponds to a probability. If for every variable $p \in \mathcal{P}$ it holds that $w(p) + w(\neg p) = 1$ the resulting sum will always be $\le 1$. Unfortunately, the models for some problems require that $w(p) + w(\neg p) > 1$, and hence solvers cannot assume that $w(p) + w(\neg p) = 1$.

Various solvers exist to solve the PWMC problem and its unweighted version. In [22], the authors present **Ganak** is a model counter build upon **sharpSAT**, a DPLL-style model counter with component caching [26]. **Ganak** uses probabilistic component caching while ensuring guarantees on the validity of the returned count. Furthermore, the authors propose to use information about the cache in the branching heuristics and show that this is beneficial to model counters. **Ganak** can also be used for projected model counting, by first branching on the priority variables, but it does not use other specialized techniques. Further modifications of **Ganak** have been proposed and in particular it has been shown that integrating tree decomposition in the branching heuristic can have a positive impact [15]. This has led to the development of **SharpSAT-TD** [15], a (weighted) model counter, and **GPMC** [25], which is also able to solve the PWMC problem.

The **projMC** solver [17] uses another approach to solve the PWMC problem<fnref n="1" />. To compute the count over a formula $F$ they first compute a disjunctive decomposition from a model of $F$. They then use the pairwise incompatible parts of the decomposition to simplify $F$, and they recursively solve the new sub-problems.

The aforementioned solvers use a strategy in which a limited amount of memory is used. This is also the focus of our work. An alternative strategy is to use a model compilation in combination with dynamic programming. Recently, Dudek et al. proposed the two-phased **ProCount** solver [11]. In the first phase, the input formula is transformed into a graded project-join tree (by a *planner*). Then, in the second phase, an executor (based on algebraic decision diagrams) is used to compute the count, using a dynamic programming approach.

### 3 Problem Definition

First, let us note that, since this work focuses on probabilistic problems, the weights on the priority variables are used to model probabilities. Hence, we refer to those as *probabilistic variables* while the set of non-priority variables are called *deterministic variables*. In the rest

<footnote n="1">¹ Although the original paper only describes the unweighted problem, a parameter in the solver enables the retrieval of the count as a float, taking into account the weights</footnote>
15:4 Probabilistic Inference by Projected Weighted Model Counting on Horn Clauses

of this paper, we denote with $P(E)$ the probability that an event $E$ occurs. In this work, we introduce a combination of two novelties in how to model problems using PWMC. First, a constraint on the clauses in $F$ is imposed: we only allow clauses to be *Horn* clauses.

► **Definition 1 (Horn clause).** A *Horne clause* $C$ is a formula of the form

$$v_1 \wedge \dots \wedge v_n \Rightarrow v_t$$

where $I = v_1 \wedge \dots \wedge v_n$ is called the implicant of the clause and $h = v_t$ is the head of the clause. Here $v_i \in \mathcal{V}$ is a variable, and $v_t$ is either a variable in $\mathcal{V}$ or $\perp$. If $n = 0$ (the implicant is empty) then the left-hand side reduces to $\top$.

It can be observed that when a Horn clause is written as an implication, as above, all the literals in the clause have the same (positive) polarity. Hence, in order to simplify our discussion and notation, we only talk about variables, and not literals. A Horn clause $C_i$ can be identified uniquely by its implicant $I_i$ and its head $h_i$. In the rest of this paper, we will use the notation $C_i$ and $(I_i, h_i)$ interchangeably. For simplicity of notation, we also denote by $v \in I_i$ the fact that $v \in \mathcal{V}$ is a variable of the implicant of $C_i$.

Horn clauses have been well studied in the literature. An important result is that the SAT problem over a CNF formula of Horn clauses can be solved in linear time [10]; given that the SAT problem in its general form is NP hard, this is a significant simplification.

Secondly, we add the notion of *distributions over the probabilistic variables*. We assume that each probabilistic variable $p \in \mathcal{P}$ belongs to exactly one partition $P_i \subseteq \mathcal{P}$ of the probabilistic variables. We define a *distribution* over each such partition, in the following simple manner: we require that one weight is specified for each probabilistic variable and require that $\sum_{p \in P_i} w(p) = 1$ for all variables in the partition. Subsequently, we calculate the weight of an assignment $S$ to the probabilistic variables as follows. The weight of a partition, given the assignment $S$, is defined as follows:

$$w_{P_i}(S) = \begin{cases} w(p) & \text{if there is exactly one } p \in P_i \text{ for which } S[p] = \top \\ 0 & \text{otherwise,} \end{cases}$$

or, in other words, if exactly one variable in the partition is set to $\top$, the weight of that variable is given to that partition; otherwise, the assignment is invalid. Implicitly, we allow only one variable within a partition to be true at the same time. The weight of the assignment $S$ is the product of the partition weights given $S$ and thus the solution of the PPHMC problem is given by

$$\sum_{S \in S_F} \prod_{i} w_{P_i}(S),$$

where we assume there is at least one partition of probabilistic variables.

## 3.1 Models

As our earlier discussion makes clear, in this work we study a simpler modeling language in which non-Horn clauses are not allowed, and we combine this with a different approach to weighting. In this section, we will show that even though we apply the aforementioned restrictions, a number of different problems can still be modeled in our language.
**► Example 2 (Bayesian Networks).** A Bayesian Network (BN) is a probabilistic model that can be represented by a directed acyclic graph. Random variables are represented by nodes and conditional dependency relationships by edges. Figure 1(a) shows an example of a BN with four nodes. For simplicity, we illustrate the encoding in our language for a network with nodes that have two values, but this generalizes for nodes with more than two values.

In this network, both *B* and *C* depend on *A*, while *D* depends on *B* and *C*. In addition to the network structure, conditional probability tables (CPTs) give, for each node, the probability of its values conditioned by its parent's value. These are called the *parameters* of the network and for conciseness we write $P(x \mid \mathbf{u})$ for the parameter corresponding to the probability that a node *X* takes value *x* given that its parents take values $\mathbf{u} = u_1, \dots, u_n$.

Various CNF encodings have been proposed for Bayesian networks [2, 4, 5, 7]. Even though the author of these works use Horn clauses, the un-projected nature of their target solvers imposes additional, non-Horn, clauses.

Let us present our encoding for BN by first defining the logical variables. For every value *x* of node *X* in the network, we define one deterministic variable $v^x$. For the BN in Figure 1, we have the following variables: $v^{a_0}$, $v^{a_1}$, $v^{b_0}$, $v^{b_1}$, $v^{c_0}$, $v^{c_1}$, $v^{d_0}$, $v^{d_1}$. Moreover for each parameter $P(x \mid \mathbf{u})$ we define a corresponding probabilistic variable $p_{\mathbf{u}}^x$. For the CPT of node B, we introduce four such variables: $p_{a_0}^{b_0}$, $p_{a_0}^{b_1}$, $p_{a_1}^{b_0}$, $p_{a_1}^{b_1}$.

In our approach, we have to define the distributions over the probabilistic variables. The CPTs of the network give a natural partition of the probabilistic variables. That is, for a node *X* in the network, we define one partition for each line of its CPT. Hence, for node *B* there are two partitions: $P_{B_1} = \{p_{a_0}^{b_0}, p_{a_0}^{b_1}\}$ and $P_{B_2} = \{p_{a_1}^{b_0}, p_{a_1}^{b_1}\}$. Next, we define the weight on the variables. Contrary to previous encodings, only weights on probabilistic variables are needed, and we use as weight the parameter they represent: $w(p_{\mathbf{u}}^x) = P(x \mid \mathbf{u})$. For example, we have that $w(p_{a_0}^{b_0}) = 0.6$ and $w(p_{a_1}^{b_1}) = 0.4$.

Finally, we define the clauses: for each parameter $P(x \mid \mathbf{u})$ with $\mathbf{u} = u_1, \dots, u_n$, we introduce one clause $v^{u_1} \wedge \dots \wedge v^{u_n} \wedge p_{\mathbf{u}}^x \Rightarrow v^x$. The clauses for the BN in Figure 1 are shown below.

$$
\begin{array}{l @{\quad} l @{\quad} l @{\quad} l}
p^{a_0} \Rightarrow v^{a_0} & v^{a_0} \wedge p_{a_0}^{c_0} \Rightarrow v^{c_0} & v^{b_1} \wedge v^{c_0} \wedge p_{b_1 c_0}^{d_0} & \Rightarrow v^{d_0} \\
p^{a_1} \Rightarrow v^{a_1} & v^{a_0} \wedge p_{a_0}^{c_1} \Rightarrow v^{c_1} & v^{b_1} \wedge v^{c_0} \wedge p_{b_1 c_0}^{d_1} & \Rightarrow v^{d_1} \\
v^{a_0} \wedge p_{a_0}^{b_0} \Rightarrow v^{b_0} & v^{a_1} \wedge p_{a_1}^{c_0} \Rightarrow v^{c_0} & v^{b_0} \wedge v^{c_1} \wedge p_{b_0 c_1}^{d_0} & \Rightarrow v^{d_0} \\
v^{a_0} \wedge p_{a_0}^{b_1} \Rightarrow v^{b_1} & v^{a_1} \wedge p_{a_1}^{c_1} \Rightarrow v^{c_1} & v^{b_0} \wedge v^{c_1} \wedge p_{b_0 c_1}^{d_1} & \Rightarrow v^{d_1} \\
v^{a_1} \wedge p_{a_1}^{b_0} \Rightarrow v^{b_0} & v^{b_0} \wedge v^{c_0} \wedge p_{b_0 c_0}^{d_0} & v^{b_1} \wedge v^{c_1} \wedge p_{b_1 c_1}^{d_0} & \Rightarrow v^{d_0} \\
v^{a_1} \wedge p_{a_1}^{b_1} \Rightarrow v^{b_1} & v^{b_0} \wedge v^{c_0} \wedge p_{b_0 c_0}^{d_1} & v^{b_1} \wedge v^{c_1} \wedge p_{b_1 c_1}^{d_1} & \Rightarrow v^{d_1}
\end{array}
$$

The clause $v^{a_0} \wedge p_{a_0}^{b_0} \Rightarrow v^{b_0}$, for instance, represents that if *A* has value $a_0$ and we pick the variable $p_{a_0}^{b_0}$ from the distribution $\{p_{a_0}^{b_0}, p_{a_0}^{b_1}\}$, *B* will have value $b_0$; we believe this is a natural and simple representation that directly reflects the BN.

Note that we can satisfy all clauses by setting the *v* variables to true. A common inference problem on Bayesian Networks is that of calculating a probability $P(X = x)$. We can solve this problem by adding a clause $v^{x'} \Rightarrow \perp$ for each value *x'* of *X* such that *x' ≠ x*. Effectively this removes from the sum those assignments in which $X \neq x$.

This encoding differs in a number of ways from the encodings used in WMC [2, 6]. The general idea is similar: for rows of the CPT, clauses are created; probabilistic variables receive weights that represent entries in the CPTs. Compared to earlier encodings, we do not generate clauses to impose that the indicators variables ($v^x$) are mutually exclusive for a node *X*. Our weighting scheme takes care of this. Furthermore, the earlier encodings have
15:6 Probabilistic Inference by Projected Weighted Model Counting on Horn Clauses

<figure><img src="image_136.png" /><figcaption>■ **Figure 1** a) An example of Bayesian Network with four binary variables. In this network *B* and *C* depends on *A* and *D* depends on *B* and *C*. The probability tables are given next to the nodes. b) An example of probabilistic network for reliability problems. The numbers labeled on the edges are their probability of being present.</figcaption></figure>

the parameter variables as the head of the implications in the CNFs, while in our encoding
they are in the implicants; while both representations are equivalent, we believe that in our
representation the structure of the BN is more closely reflected in the clauses.

**► Example 3 (Reliability in Networks).** Reliability in network (RN) problems study the connectivity of nodes in probabilistic graphs. In such graphs, as shown in Figure 1(b), each edge has a probability of being present. In this work, we consider the computation of the probability that two nodes are connected.

More formally, let $G = (V, E)$ be a probabilistic graph and $f_w : E \to [0, 1]$ a weighting function that assign to each edge a probability of being present or not. We denote $s$ the source node, $t$ the target node and $R_t^s$ ($\bar{R}_t^s$) the fact that $t$ is (not) reachable from $s$. The goal is then to compute $P(R_t^s)$.

The encoding of this problem in our language is similar to that in [12], in which the authors propose to compute $P(\bar{R}_t^s)$ and then use the fact that $P(R_t^s) = 1 - P(\bar{R}_t^s)$ to answer the initial query.

Let us first define the logical variables. For each node $X \in V$, we introduce one
deterministic variable $v_X$. For each edge $e \in E$ from $u$ to $v$, we introduce two probabilistic
variables $p_{uv}$ (the edge is present) and $\bar{p}_{uv}$ (the edge is not present). The weighting scheme
of the probabilistic variables uses the weighting function of the edges: $w(p_{uv}) = f_w(e)$ and
$w(\bar{p}_{uv}) = 1 - f_w(e)$. A distribution is defined for each edge, containing these two variables:
$P_e = \{p_{uv}, \bar{p}_{uv}\}$. Since the probabilistic variables for each edge are in their own distribution,
and each distribution contains no other variables, it is easy to see that an assignment to the
probabilistic variables corresponds to a possible instance of the graph.

We will use the deterministic variables $v_X$ to represent whether in the instance implied
by the probabilistic variables, $X$ is reachable from $s$. The clauses, that we define hereafter,
must ensure that an assignment is a model only if $v_t = \perp$ given the choices for the edges.

To ensure that, the clauses use the transitive nature of the connectivity. That is, for each
edge $e$ from $u$ to $v$, a clause $v_u \wedge p_{uv} \Rightarrow v_v$ is created. This can be interpreted as that if
$u$ is reachable from $s$ and $e$ is present, then $v$ is also reachable from $s$. We impose that $s$
is reachable from $s$ by adding the clause $\top \Rightarrow s$ and that $t$ is not reachable by adding the
clause $v_t \Rightarrow \perp$. The clauses for the query $P(\bar{R}_E^A)$ for the graph in Figure 1 are shown below.

$v_A \wedge p_{AB} \Rightarrow v_B$

$v_B \wedge p_{BD} \Rightarrow v_D$

$v_C \wedge p_{CE} \Rightarrow v_E$

$v_E \Rightarrow \perp$

$v_A \wedge p_{AC} \Rightarrow v_C$

$v_C \wedge p_{CD} \Rightarrow v_D$

$v_D \wedge p_{DE} \Rightarrow v_E$

$\top \Rightarrow v_A$
► **Example 4 (Probabilistic programming).** ProbLog is a probabilistic programming language that extends Prolog with probabilistic predicates [8]. It can be used to represent both the aforementioned inference problems. While a full discussion of this language is beyond the scope of this paper, we wish to illustrate how PPHMC can be used in ProbLog. Consider the following example ProbLog program, taken from the ProbLog website<fnref n="2" />:

```
0.4 :: heads.
0.3 :: col(1,red); 0.7 :: col(1,blue).
0.2 :: col(2,red); 0.3 :: col(2,green); 0.5 :: col(2,blue).
win :- heads, col(_,red).
win :- col(1,C), col(2,C).
query win.
```

For such a logic program, ProbLog performs *grounding*, creating a propositional version of the program. This propositional version can be represented as follows in our language:

$$ 
\begin{align*}
p_{\text{heads}}, p_{\text{col}(1,\text{red})} &\Rightarrow v_{\text{win}} & p_{\text{col}(1,\text{blue})}, p_{\text{col}(2,\text{blue})} &\Rightarrow v_{\text{win}} \\
p_{\text{heads}}, p_{\text{col}(2,\text{red})} &\Rightarrow v_{\text{win}} & v_{\text{win}} &\Rightarrow \perp \\
p_{\text{col}(1,\text{red})}, p_{\text{col}(2,\text{red})} &\Rightarrow v_{\text{win}}
\end{align*}
$$

Here we have these partitions: $w(p_{\text{heads}}) = 0.4$, $w(\bar{p}_{\text{heads}}) = 0.6$; $w(p_{\text{col}(1,\text{red})}) = 0.3$, $w(p_{\text{col}(1,\text{blue})}) = 0.7$; $w(p_{\text{col}(2,\text{red})}) = 0.2$, $w(p_{\text{col}(2,\text{green})}) = 0.3$, $w(p_{\text{col}(2,\text{blue})}) = 0.5$. The probability of the query is obtained by PPHMC on this ground version. Given that our previous example showed how reliability problems can be modeled in polynomial space, while the model for this problem is exponential in the grounded ProbLog model, we hypothesize that this is possible for all grounded programs without cycle breaking.

# 4 Algorithm

In the section we present the main algorithms of our solver. At the core, the problem is solved by a backtracking search over the possible assignments for the distributions. When a value is assigned to a variable, we call a new propagator designed specifically for the structure of the clauses. We first present the general algorithm for a search-based solver for the PWMC problem and then the specific propagation algorithm as well as some branching heuristics.

## 4.1 General Approach

The main algorithm is shown in Algorithm 1. In essence, it is similar to other DPLL-based solvers. The computation of the count for a non-empty formula $F$ starts by looking into a cache (line 7) to determine if the formula has already been counted. If not, then it chooses an unfixed variable (lines 8, 10) and assigns it a value. A residual formula is computed after calling a propagation algorithm (line 11) which is then divided into independent components (line 14). The components are then solved independently (line 17) and their count is stored in the cache (line 18). In order to bound the memory consumption of search based solvers, the cache has a limited number of entries. As cache cleaning techniques are not the focus of this work, we simply fully clear the cache when the limit is reached.

<footnote n="2">² https://dtai.cs.kuleuven.be/problog/tutorial/advanced/00_inference.html</footnote>
15:8 Probabilistic Inference by Projected Weighted Model Counting on Horn Clauses

<table><caption>■ Algorithm 1 General PPHMC search algorithm.</caption>
  <tbody>
    <tr>
      <td>1</td>
      <td>Function PWMC(F, P)</td>
    </tr>
    <tr>
      <td></td>
      <td style="padding-left: 2em;">input : A Horn-CNF formula F with probabilistic variables P</td>
    </tr>
    <tr>
      <td></td>
      <td style="padding-left: 2em;">output: The projected (on P) weighted model count of F</td>
    </tr>
    <tr>
      <td>2</td>
      <td style="padding-left: 2em;">C &larr; newCache()</td>
    </tr>
    <tr>
      <td>3</td>
      <td style="padding-left: 2em;">return PWMCr(F, P, C)</td>
    </tr>
    <tr>
      <td>4</td>
      <td></td>
    </tr>
    <tr>
      <td>5</td>
      <td>Function PWMCr(F, P, C)</td>
    </tr>
    <tr>
      <td></td>
      <td style="padding-left: 2em;">input : F, P same as PWMC()</td>
    </tr>
    <tr>
      <td></td>
      <td style="padding-left: 2em;">input : C the cache of previously found counts</td>
    </tr>
    <tr>
      <td></td>
      <td style="padding-left: 2em;">output: Same as PWMC()</td>
    </tr>
    <tr>
      <td>6</td>
      <td style="padding-left: 2em;">if P = &emptyset; then return 1</td>
    </tr>
    <tr>
      <td>7</td>
      <td style="padding-left: 2em;">if F &isin; C then return C[F]</td>
    </tr>
    <tr>
      <td>8</td>
      <td style="padding-left: 2em;">P &larr; a distribution of P such that &exist;v &isin; P | v is not fixed</td>
    </tr>
    <tr>
      <td>9</td>
      <td style="padding-left: 2em;">count &larr; 0</td>
    </tr>
    <tr>
      <td>10</td>
      <td style="padding-left: 2em;">foreach v &isin; P | v is not fixed do</td>
    </tr>
    <tr>
      <td></td>
      <td style="padding-left: 4em;">/* Assign v = &top; and call the propagation algorithm. Returns the residual formula F' and the unconstrained probability of F given F'. */</td>
    </tr>
    <tr>
      <td>11</td>
      <td style="padding-left: 4em;">(F', U<sub>F|F'</sub>) &larr; Propagate(F, P, v)</td>
    </tr>
    <tr>
      <td>12</td>
      <td style="padding-left: 4em;">proba &larr; U<sub>F|F'</sub></td>
    </tr>
    <tr>
      <td>13</td>
      <td style="padding-left: 4em;">if F' is not UNSAT then</td>
    </tr>
    <tr>
      <td>14</td>
      <td style="padding-left: 6em;">Components &larr; all connected components of F'</td>
    </tr>
    <tr>
      <td>15</td>
      <td style="padding-left: 6em;">foreach Comp &isin; Components do</td>
    </tr>
    <tr>
      <td>16</td>
      <td style="padding-left: 8em;">P' &larr; P reduced to the variables in Comp</td>
    </tr>
    <tr>
      <td>17</td>
      <td style="padding-left: 8em;">proba<sub>Comp</sub> &larr; PWMCr (Comp, P', C)</td>
    </tr>
    <tr>
      <td>18</td>
      <td style="padding-left: 8em;">C[Comp] &larr; proba<sub>Comp</sub></td>
    </tr>
    <tr>
      <td>19</td>
      <td style="padding-left: 8em;">proba &larr; proba * proba<sub>Comp</sub></td>
    </tr>
    <tr>
      <td>20</td>
      <td style="padding-left: 6em;">end</td>
    </tr>
    <tr>
      <td>21</td>
      <td style="padding-left: 6em;">count &larr; count + proba</td>
    </tr>
    <tr>
      <td>22</td>
      <td style="padding-left: 6em;">end</td>
    </tr>
    <tr>
      <td>23</td>
      <td style="padding-left: 6em;">end</td>
    </tr>
    <tr>
      <td>24</td>
      <td style="padding-left: 6em;">return count</td>
    </tr>
  </tbody>
</table>

There are a few differences between our solver and other model counters that need to be pointed out. First, let us note that the variables in a distribution are mutually exclusive. Indeed, if a distribution does not have exactly one variable set to $\top$ in an assignment, then its weight is 0 and $\prod_i w_{P_i}(S) = 0$, which does not contribute to the count.

A consequence of this is that, unlike classical DPLL-based counters, the branching decision is made on distributions and not variables. Indeed, since the distributions partition of the probabilistic variables (i.e., the variables on which the number of weighted models must be counted), assigning one variable to $\top$ in each distribution means that all other probabilistic variables are set to $\bot$. Moreover, when fixing a variable $v$ in the selected distribution, only the case of fixing to $\top$ needs to be explored. The case $v = \bot$ is explored when the other variables in the distribution are selected for branching. Notice that when there is no distribution left in $F$, then 1 is returned as the remaining formula is SAT (line 6). Indeed, there are no unit clauses (they are removed during the propagation) and all that remains in $F$ are
Horn clauses with deterministic variables. By setting all remaining variables to $\perp$ we obtain a model of this formula, as the implicant of all clauses evaluates to $\perp$. This would not be possible if $F$ was not solely composed of Horn clauses. In this case, the SAT problem needs to be solved on $F$, which is NP-hard.

The computation of independent components is also slightly different. In traditional model counters, when the input formula $F$ can be decomposed into multiple sub-formulas that do not share any variable, each sub-formula is solved independently and the count of $F$ is the product of the sub-formulas' count. However, since the probabilistic variables are linked in partitions, we must add as condition that two independent sub-formulas cannot share any variable in the same distribution.

Finally, we also devised a new propagation algorithm, explained in the next section. It returns (line 11) a residual formula $F'$ as well as what we call the unconstrained probability of $F$ given $F'$: $U_F|_{F'}$. This probability accounts for the distributions of $F$ that are unconstrained in $F'$. Indeed, if for a distribution $P_i = \{p_i^1, \dots, p_i^m\}$ occurring in $F$, there are no clause in $F'$ that contains any of the $p_i^k \in P_i$, $P_i$ will never be selected by the branching heuristic. Moreover, branching on it will not impact $F'$. Hence, the probability obtained by branching of $P_i$ can be precomputed and is given by $\sum_{p \in P_i} w(p)$. Note that this sum may not be 1 if propagation fixed one of the $p_i^k$ variables to $\perp$ earlier, and $|P_i| > 2$. More generally, for $l$ such distributions $P_{u_1}, \dots, P_{u_l}$ we have $U_F = \prod_{i=1}^l \sum_{p \in P_{u_i} | p \text{ is not fixed}} w(p)$.

## 4.2 Propagation

In this section, we describe the propagation algorithm used by our solver, summarized in Algorithm 2. In brief, we first apply the classical Boolean Unit Propagation (BUP), in which the links between probabilistic variables in partitions are also enforced, until a fixed point is reached. Then we remove from the remaining formula the clauses that do not constrain the count anymore. First, let us detail what the BUP does when a variable $v$ is fixed in $F$:

- If $v = \top$ then

- If $v$ is a probabilistic variable in a distribution $P$, apply the BUP on all $v' \in P, v' \neq v$ with $v' = \perp$.

- Every clause $C = (I, h)$ such that $h = v$ is removed from $F$. Indeed, for every assignment on its remaining variables, it evaluates to $\top$.

- For every clause $C = (I, h)$ such that $v \in I$, replace $I$ by $I' = I \setminus \{v\}$.

- If $v = \perp$ then

- If $v$ is a probabilistic variable in a distribution $P$ and only one variable $v'$ remains unfixed in $P$, apply BUP with $v' = \top$.

- Every clause $C = (I, h)$ such that $v \in I$ is removed from $F$ for the same reason as above.

- The head of every clause $C = (I, h)$ such that $h = v$ is replaced by $\perp$.

There are two cases in which, after a call to the BUP, a variable is forced to take a value, and the BUP algorithm needs to be called again. First, when the last variable of an implicant is removed from it, then the head of the clause must be $\top$. Secondly, when a clause $(I, h)$ has its head set to $\perp$ and there is only one variable left in $I$, then this variable must be set to $\perp$. In Algorithm 2, the call to BUP at line 2 is executed until no such clauses can be found.

However, in the context of projected model counting, a key insight for our work is that further propagation can be done that is not entailed by BUP. The intuition is that the deterministic variables can be used to remove additional clauses from $F$. For example, if our formula includes a clause $\neg a \lor \neg b \lor \neg c$ with $a, b$ being probabilistic variables and $c$ a
15:10 Probabilistic Inference by Projected Weighted Model Counting on Horn Clauses

<figure><img src="image_137.png" /><figcaption>Figure 2 On the left: an instance of a reliability problem on probabilistic networks. In the center: the associated clauses for the query $P(\bar{R}_t^s)$. On the right: the graph of clauses implication: there is one node per clause and a link between $C_i$ and $C_j$ if $h_i \in I_j$.</figcaption></figure>

deterministic variable that does not appear in any other clauses, then setting $c$ to $\perp$ makes the clause evaluate to $\top$, regardless of the choice for $a$ and $b$. Since we are only interested in finding an assignment to the deterministic variables, this assignment does not impact the final count.

Let us show how this works on the clauses in Figure 2, taken from a reliability query in a probabilistic networks, in which we want to compute $P(\bar{R}_t^s)$, the probability that $s$ and $t$ are not connected in the graph on the left. Applying BUP gives $v_s = \top$ and $v_t = \perp$, which removes these variables from $C_1, C_3$ and $C_6$. No further propagation can be done with BUP. However, when looking at the graph on the left, it is clear that only the nodes $B, C$ and $D$ impact the connectivity between $s$ and $t$, but this is not detected by BPU. For instance, let us look at nodes $A$ and the associated clause $C_0$, which contains the distribution $\{e_{AB}, \bar{e}_{AB}\}$. It can be seen that for both choices for the edge $e_{AB}$, setting $v_A = \perp$ (a deterministic variable) reduces the clause $C_0$ to $\top$. Since there are no clauses that have $v_A$ in their head, setting $v_A$ to $\perp$ has no impact on the other variables in $F$. Moreover, since it is a deterministic variable, it does not impact the projected weighted model count. A similar reasoning can be made for $H$ and $C_{11}$ by setting $v_H = \top$, since $v_h$ does not appear in any implicant in $F$.

The intuition behind our propagation is the following. In order for an assignment to not be a model of the input formula, it must generate a clause $\top \Rightarrow \perp$. Some clauses cannot contribute to such a contradiction, by setting a deterministic variable to $\perp$ in its implicant or its head to $\top$. We formalize this intuition next.

First, let us define in which case we cannot know if a clause will have its head set to $\perp$ or its implicant set to $\top$.

**► Definition 5** ($\{\top, \perp\}$-reachability). A clause $C_i = (I_i, h_i) \in F$ is $\perp$-reachable if one of the two following conditions is met:

1. $C_i$ is of the form $I_i \Rightarrow \perp$ or $I_i \Rightarrow p$ with $p \in \mathcal{P}$

2. There exists a clause $C_j = (I_j, h_j) \in F$ such that $C_j$ is $\perp$-reachable and $h_i \in I_j$. Similarly, $C_i$ is $\top$-reachable if one of the two following conditions is met:

1. There exists no deterministic variables in $I_i$

2. There exists a clause $C_j = (I_j, h_j) \in F$ such that $C_j$ is $\top$-reachable and $h_j \in I_j$.

We say that a clause is constrained if it is $\perp$-reachable and $\top$-reachable. The $\{\perp, \top\}$-reachability can be seen for the clauses in Figure 2 on the right, on the implication graph of the clauses. The implication graph of a set of Horn clauses is a graph $G = (V, E)$ such that there is one node per clause and an edge from a clause $C_i = (I_i, h_i)$ to $C_j = (I_j, h_j)$ if $h_i \in I_j$. If a clause is $\top$-reachable ($\perp$-reachable), so are all its descendants (ancestors) in the implication graph.
<table><caption>**Algorithm 2** Propagation algorithm of Schlandals.</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>1</td><td>Function Propagate(F, P, v)</td></tr><tr><td></td><td>input : A boolean formula F with probabilistic variables P</td></tr><tr><td></td><td>input : A variable v set to ⊤</td></tr><tr><td></td><td>output: The residual formula F' and a propagation probability p<sub>prog</sub></td></tr><tr><td></td><td>/* Call the BUP procedure with the initial assignment of ⊤ to v until fix point is reached */</td></tr><tr><td>2</td><td>F' ← BUP(F, v, ⊤)</td></tr><tr><td>3</td><td>foreach C = (I, h) ∈ F' do</td></tr><tr><td>4</td><td>if h = ⊥ or h ∈ P then SetFReachable(F', C)</td></tr><tr><td>5</td><td>if ∃v ∈ I | v ∉ P then SetTReachable(F', C)</td></tr><tr><td>6</td><td>end</td></tr><tr><td>7</td><td>foreach C ∈ F' do</td></tr><tr><td>8</td><td>if C is not ⊤-reachable or C is not ⊥-reachable then</td></tr><tr><td>9</td><td>F' ← F' - C</td></tr><tr><td>10</td><td>end</td></tr><tr><td>11</td><td>end</td></tr><tr><td>12</td><td>U<sub>F|F'</sub> ← 1</td></tr><tr><td>13</td><td>foreach distribution P such that P ∈ F ∧ P ∉ F' do</td></tr><tr><td>14</td><td>U<sub>F</sub> ← U<sub>F</sub> * ∑<sub>p∈P|p</sub> is not fixed w(p)</td></tr><tr><td>15</td><td>end</td></tr><tr><td>16</td><td>return (F', U<sub>F|F'</sub>)</td></tr></tbody></table>

We have the following theorem, which states that if a clause is not constrained, then it can be removed safely (without impacting the count) from F.

**Theorem 6.** Let $F = C_1 \wedge \dots \wedge C_n$ be a formula with $n$ Horn clauses over the variables $V$, $P \subseteq V$ the set of probabilistic variables and $D = V \setminus P$ the set of deterministic variables. Let $C_{u_1}, \dots, C_{u_k}$ be $k$ unconstrained clauses of $F$ with $C_{u_i} = (I_{u_i}, h_{u_i})$.

There exists a subset of $k$ deterministic variables $X = \{x_1, \dots, x_k\} \subseteq D$ with $x_i \in I_{u_i} \cup \{h_{u_i}\}$ and an assignment $S_X$ on $X$ such that

$$S_F = S_{F[S_X]}$$

where $S_F$ denotes the set of models of $F$, projected on $P$, and $F[S_X]$ the formula obtained by applying the BUP algorithm on $F$ with the assignment $S_X$.

We now prove this theorem and give the procedure to find the assignment on the deterministic variables.

**Proof.** Let $G = (V, E)$ be the graph of the implications of $F$ and $C_i$ an unconstrained clause in $F$. We prove that an assignment can be found for one of the deterministic variables of $C_i$ such that it does not impact the count of $F$.

First, let us assume that $C_i$ is not $⊥$-reachable. We denote $G_{C_i}^d = (V_{C_i}^d, E_{C_i}^d)$ the sub-graph of $G$ that contains $C_i$ and all its descendants. By definition there is no clause $C_j \in V_{C_i}^d$ that is $⊥$-reachable, otherwise $C_i$ would be $⊥$-reachable. Hence, all clauses in $V_{C_i}^d$ are unconstrained. Let $X^d = \cup_{C_j \in V_{C_i}^d} \{h_j\}$ be the set of (deterministic) heads in $G_{C_i}^d$. We define the assignment $S_{X^d}$ such that $S_{X^d}[x] = ⊤, \forall x \in X^d$. This removes all the clauses in $V_{C_i}^d$ from $F$ without impacting the values of the other variables in the clauses.
15:12 Probabilistic Inference by Projected Weighted Model Counting on Horn Clauses

Next, if $C_i$ is $\perp$-reachable but unconstrained, then it is not $\top$-reachable. Let $G_{C_i}^p = (V_{C_i}^p, E_{C_i}^p)$ be the sub-graph of $G$ that contains $C_i$ and all its parents. Since $C_i$ is $\perp$-reachable, so is every clause in $V_{C_i}^p$. Thus we have that every clause $C_j \in V_{C_i}^p$ has at least one deterministic node $d_j \in I_j$, otherwise they would be $\top$-reachable. Let $X^p = \cup_{C_j \in V_{C_i}^p} \{d_j\}$ be the set of such nodes. We set $S_{X^p}$ such that $S_{X^p}[d_j] = \perp$ for all clause $d_j \in X^p$ which removes all clauses in $V_{C_i}^p$ from $F$ without constraining the other variables. ◀

<caption>■ Algorithm 3 Procedure to mark the clauses of a formula $F$ as $\perp$-reachable.</caption>

1 **Function SetFReachable(F, C)**
| **input** : A boolean formula $F$, a clause $C = (I, h)$ of $F$
| **if** $C$ is not marked as $\perp$-reachable **then**
| | Mark $C$ as $\perp$-reachable
| 2
| 3 | **foreach** $C' = (I', h') \in F \mid h' \in I$ **do** SetFReachable(F, C')
| 4
| 5 | **end**

This additional propagation is shown in lines 3-11 of Algorithm 2. After the application of BUP until a fix point is reached, the remaining clauses are iterated over. If a clause is of the form $C = (I, \perp)$ or has a head with an unfixed probabilistic variable, then the procedure SetFReachable is called, for which the code is shown in Algorithm 3. This algorithm basically traverses the implication graph of $F$, starting from $C$ and marks every clause it encounters as $\perp$-reachable. Notice that it does not mark multiple times the same clause. Hence, the cost of marking all the $\perp$-reachable clauses is $O(n)$ with $n$ the number of clauses in $F'$. A similar procedure is defined for the $\top$-reachability, but we do not include it for conciseness. After marking the clauses (lines 3-5), every unconstrained clause is removed from the formula obtained after BUP (lines 7-11). Finally, the algorithm computes the unconstrained probability $U_{F|F'}$ (lines 12-15).

It can be noted that in some cases our propagation is similar to *Pure Literal Elimination* (PLE) on the deterministic variables. For instance, in Figure 2, $v_A$ never appears in the head of any clause. Hence, it can be set as $\perp$, which is also the value found by our procedure. However, PLE cannot detect that the nodes E, F and G are not useful for the query. It should also be noted that although our propagation implies that PLE is performed, the procedure only reasons about the clauses, which makes it more efficient.

4.3 Branching Heuristic

At line 8 of Algorithm 1, a distribution is selected within the set of distributions for which no element is set to $\top$. In this section we explain the heuristics implemented in our solver to choose the distribution. We must note that our solver is not a conflict-driven clause learning (CDCL) solver, making it impossible to use the heuristics of such solvers [19, 22]. Instead, we implemented some easy to understand, fast to evaluate, but effective heuristics based on the implication graph of a formula $F$. We provide three heuristics related to the degree of a clause in the implication graph. These heuristics select a clause with i) the lowest in-degree ii) the lowest out-degree and iii) the maximum degree, and then a distribution in this clause. Each time they are called, they re-evaluate the score of each clause in the current sub-formula.
# 5 Results

In order to evaluate the effectiveness of our approach, we compare our solver<fnref n="3" /> with GPMC [25] and ProjMC [17] on Bayesian networks and reliability queries in probabilistic networks. We chose these solvers as they were the best performing on the Projected Model Counting and Projected Weighted Model Counting tracks of the 2022 model counting competition. We did not include the proCount solver [11] as it is a knowledge compilation based solver. Ganak [22] also was not used as we found it was the worst performing on the projected tasks (as also observed in the 2022 competition). For the reliability queries, we also ran ApproxMC [3, 23, 24] as it is known to perform well on this problem, although it should be stressed ApproxMC solves an approximate model counting problem, while we solve an exact model counting problem. For each of the two problems we first present the methodology used to generate the instances and then present the results<fnref n="4" />. For each instance, a timeout of 600 seconds and a memory limit of 15 GB were set.

## 5.1 Bayesian Networks

The Bayesian networks come from the bnlearn R package [21] repository and range from small networks (fewer than 20 nodes) to large (up to hundreds of nodes). We selected all the networks in the repository for which at least one solver did not time out. For each of these networks, the queries are done on the leaves of the networks without any evidence. Hence, for a leaf L of the network that takes values $l_1, \dots, l_n$, we create n instances for the queries $P(L = l_1), \dots, P(L = l_n)$.

The goal of this experiment is to answer the following questions: i) how do the evaluated PWMC solvers perform using our simple encoding? ii) How does our encoding perform compared to other state-of-the-art, but more complex encodings for BNs, designed for weighted model counters? For i) we use a CNF formula in other PWMC solvers similar to our model, where we add additional clauses to enforce the distribution constraints that we have in our weighting scheme. For ii), we follow the nomenclature in [6], and we compare our encoding with ENC1 [7], ENC3 [4], and ENC4 [5]. We also add the encoding recently presented in [2], which we denote ENC4linp. Briefly, ENC1 is the simplest encoding: there is one indicator variable $\lambda^x$ for each value x of node X and one parameter variable $\theta_u^x$ per parameter $P(x \mid u)$ of the network. For the clauses, the indicator variables of the same node are mutually exclusive; a clause $\lambda^{u_1} \wedge \dots \wedge \lambda^{u_n} \wedge \lambda^x$ $\Leftrightarrow$ $\theta_u^x$ is created for each parameter $P(x \mid u)$. The main addition of ENC3 is that, for a given CPT in the network, the same variable $\theta_u^x$ is reused for all the entries having the same probability. In ENC4, The CPTs are simplified using a modified version of Quine/McCluskey algorithm for finding prime implicants, resulting in a better decomposition of the input CNF. Finally, in ENC4linp, the authors propose two novelties for encoding BNs in a CNF: the use of log-encoding for the elementary assignment of a variable and one parameter variable per CPT is kept implicit. For these encodings, GPMC and projMC are run in WMC mode (which is the d4 [16] model counter for the latter).

Figure 3 shows the percentage of instances solved within the time limit.

First, it can be seen that if we use our model in other PWMC solvers, their performance is much worse; these solvers solve a few instances.

<footnote n="3">³ Available at https://github.com/aia-uclouvain/schlandals</footnote>

<footnote n="4">⁴ The instances can be found here https://github.com/AlexandreDubray/pwmc_benchmarks</footnote>
15:14 Probabilistic Inference by Projected Weighted Model Counting on Horn Clauses

<figure><img src="image_138.png" /><figcaption>Figure 3 Percentage of instances solved in 600 seconds for the Bayesian Network data sets.</figcaption></figure>

Compared to the ENC1 model used in traditional WMC, which is arguably the model for
traditional WMC most similar to our model, we can still observe that our approach performs
significantly better, solving 80% of the instances for the minimum in degree heuristic. We
believe it is highly encouraging that our approach works so well for such a simple model.

When evaluating the more complex, optimized models ENC3, ENC4 and ENC4linp
developed for WMC, we can see that these have increasingly better results. This confirms
the benefit of the various improvements made, over the years, to the encodings. Overall, the
GPMC solver is the best performing of all, even when using only the ENC3 encoding. The
biggest gap between the ENC encodings is between ENC1 and ENC3. Indeed, using the
same parameter variables for multiple parameters of a CPT gives a huge drop in the number
of variables, especially for the hardest instances, which contain large CPTs. The benefits of
the Quine/McCluskey reduction are clearly visible for the d4 solver.

The nature of our modeling language is such that the optimizations of ENC3, ENC4 and ENC4linp cannot be directly applied to our models. Unlike the ENC encodings, we require that one variable per distribution is set to $\top$. However, we hypothesize that similar optimizations can also be developed for our solver in future work.

## 5.2 Connectivity in Probabilistic Networks

For this problem, the data sets used come from two sources. First, we used the power grid network of Europe and USA as extracted by the GridKit tool [18, 28]. The extracted graphs represent the electric power system in these geographical areas. In order to have various sizes of instances, both networks are divided by country (Europe) or by state (USA). Then for each subnetwork, five random pairs of nodes are selected and the probability that they are connected is computed. Notice that since the edges in the power-grid network have no orientation (the graph is undirected), we transformed it into a directed graph by replacing each edge between *u* and *v* by one edge from *u* to *v* and one edge from *v* to *u*. We assume that each edge has a probability of 0.125 to be down, as done in [12].

Secondly, we used (oriented) graphs representing water distribution systems from the WNTR Python package [14]. We considered as sources (sinks) nodes having no parents (children). Then, we compute $P(\bar{R}_i^s)$ between each source-target pair $(s, t)$ such that there exists a path between $s$ and $t$. As for the power grid network data set, we assume a fix probability of 0.125 that each link is down.

Unlike for Bayesian Networks, our encoding and the one proposed in [12] are almost
identical and can be reused for all the evaluated solvers. The only difference is that we do not
create two variables per edge for the encoding passed to GPMC and projMC. Figure 4 shows
the percentage of instances solved in the given amount of time for both of these data sets.
<figure><img src="image_139.png" /><figcaption>Figure 4 Percentage of instances solved in 600 seconds for the power-grid networks (left) and the water networks (right).</figcaption></figure>

First, let us note that GPMC performs the least well on both data sets. While it is able to solve up to 25% of the instances for the power grid network, it is unable to solve the largest instances of the water networks. For the latter data set, it quickly reaches the memory limit (15 GB) and stagnates until the time-out. On the other hand, it can be seen that projMC and Schlandals solve roughly the same number of instances on the power grid networks, while Schlandals solves 10% more instances on the water networks. On these types of data sets, the min-in degree and min-out degree heuristics work best for our solver. It can be seen that for the water networks, Schlandals outperforms projMC.

Interestingly, on the power grid data sets, our solver is still able to match the performance of projMC even though our propagation is not used to its full capacity. Indeed, since the original graphs are not directed, they were transformed by adding two directed edges for each undirected edge. As a result, either all clauses in a component are constrained, or no clause is constrained. Indeed, by doubling the edges in the original graph, the edges in the implication graph of the resulting CNF are also doubled. Hence, every clause which is a descendant of a $\perp$-reachable clause is also one of its ancestors and a similar argument can be made for $\top$-reachable clauses. Hence, our propagation algorithm has less pruning power. On the other hand, for the water network data set, our solver uses the full strength of our propagation and is more efficient than projMC.

Let us briefly comment on the performance of approxMC [12]. Again, it should be noted that approxMC is an approximate model counter that provides guarantees $(\delta - \epsilon)$-guarantees and does not consider weights on the literals. Hence, the problem it solves is quite different from the other solvers compared in this work: if the count returned by approxMC is $C$, and the exact count $C^*$, it ensures that $\frac{C^*}{1+\epsilon} < C < C^*(1+\epsilon)$ with a confidence of $1 - \delta$ $(\delta, \epsilon \in [0, 1])$. In Figure 4 the results are shown for $\epsilon = 0.1$ and $\delta = 0.05$. This is an example of configuration in which the solver is quite confident in its solution, and hence the results are more comparable to that of the other solvers. It can be seen that it performs very well on the power grid network instances, but poorly on the water supply network ones. Using $\epsilon = 0.8$ and $\delta = 0.2$, as reported in [12], the solver is able to solve all instances on both data sets, but it is much less confident in its solution. Overall the performance of approxMC heavily depends on the acceptable margin of error.
# 6 Conclusion

Weighted model counters have become an essential tool in probabilistic reasoning, but the CNF models have grown more and more complex. A step towards simplicity for some problems has been taken by the introduction of projected weighted model counting. In this work, we propose the Projected Probabilistic Horn model counting (PPHMC) problem, in which only Horn clauses are allowed, making the modeling language simpler. We have shown that with an appropriate weighting scheme it is possible to model important probabilistic problems as PPHMC problems. In particular, we provide an encoding for Bayesian networks, probabilistic networks and probabilistic logic programs. We also introduced a new tool, the Schlandals solver, specifically designed for our language. Our experiments show that our solver is competitive with state-of-the-art solvers and opens the path to further work on PPHMC.

As we have seen for the Bayesian Networks, the encoding can have a great impact on the performance. An interesting line of work would be to investigate how the optimizations developed during the past twenty years, for the encoding of Bayesian Networks into a logical formula, can be applied with our weighting schema, and how our solver can be integrated in probabilistic logic programming systems. On the other hand, a lot of work has been done to make model counters efficient. In the future, integrating such techniques (probabilistic caching, tree decomposition, symmetry breaking, advanced branching heuristics, etc.) with the specificity of PPHMC can also increase the performance of our solver. Finally, this paper has been focused on PWMC within bounded memory. However, the core of our solver can be reused in a knowledge compiler or an approximate solver.

# References

1 Rehan Abdul Aziz, Geoffrey Chu, Christian Muise, and Peter Stuckey. #∃ SAT: Projected model counting. In *International Conference on Theory and Applications of Satisfiability Testing 2015*. Springer, 2015.

2 Anicet Bart, Frédéric Koriche, Jean-Marie Lagniez, and Pierre Marquis. An improved CNF encoding scheme for probabilistic inference. In *Proceedings of the Twenty-second European Conference on Artificial Intelligence*, 2016.

3 Supratik Chakraborty, Kuldeep S. Meel, and Moshe Y. Vardi. Algorithmic improvements in approximate counting for probabilistic inference: From linear to logarithmic sat calls. In *IJCAI*, 2016.

4 Mark Chavira and Adnan Darwiche. Compiling Bayesian networks with local structure. In *IJCAI*, volume 5, 2005.

5 Mark Chavira and Adnan Darwiche. Encoding CNFs to empower component analysis. In *Theory and Applications of Satisfiability Testing-SAT 2006: 9th International Conference, Seattle, WA, USA, August 12-15, 2006. Proceedings 9*. Springer, 2006.

6 Mark Chavira and Adnan Darwiche. On probabilistic inference by weighted model counting. *Artificial Intelligence*, 172(6-7), 2008.

7 Adnan Darwiche. A logical approach to factoring belief networks. *KR*, 2, 2002.

8 Luc De Raedt, Angelika Kimmig, and Hannu Toivonen. Problog: A probabilistic prolog and its application in link discovery. In *IJCAI*, volume 7. Hyderabad, 2007.

9 Paulius Dilkas and Vaishak Belle. Weighted model counting with conditional weights for Bayesian networks. In *Uncertainty in Artificial Intelligence*. PMLR, 2021.

10 William F. Dowling and Jean H. Gallier. Linear-time algorithms for testing the satisfiability of propositional Horn formulae. *The Journal of Logic Programming*, 1(3), 1984.

11 Jeffrey M. Dudek, Vu H. N. Phan, and Moshe Y. Vardi. ProCount: Weighted Projected Model Counting with Graded Project-Join Trees. In Chu-Min Li and Felip Manyà, editors, *Theory of Probabilistic Models*, 2019.
and Applications of Satisfiability Testing - SAT 2021, Lecture Notes in Computer Science,
Cham, 2021. Springer International Publishing. doi:10.1007/978-3-030-80223-3_11.

12 Leonardo Duenas-Osorio, Kuldeep Meel, Roger Paredes, and Moshe Vardi. Counting-based reliability estimation for power-transmission grids. In *Proceedings of the AAAI Conference on Artificial Intelligence*, volume 31, 2017.

13 Daan Fierens, Guy Van den Broeck, Joris Renkens, Dimitar Shterionov, Bernd Gutmann, Ingo Thon, Gerda Janssens, and Luc De Raedt. Inference and learning in probabilistic logic programs using weighted Boolean formulas. *Theory and Practice of Logic Programming*, 15(3), 2015.

14 Katherine A. Klise, Michael Bynum, Dylan Moriarty, and Regan Murray. A software framework for assessing the resilience of drinking water systems to disasters with an example earthquake case study. *Environmental modelling & software*, 95, 2017.

15 Tuukka Korhonen and Matti Järvisalo. Integrating tree decompositions into decision heuristics of propositional model counters. In *27th International Conference on Principles and Practice of Constraint Programming (CP 2021)*. Schloss Dagstuhl-Leibniz-Zentrum für Informatik, 2021.

16 Jean-Marie Lagniez and Pierre Marquis. An Improved Decision-DNNF Compiler. In *IJCAI*, volume 17, 2017.

17 Jean-Marie Lagniez and Pierre Marquis. A recursive algorithm for projected model counting. In *Proceedings of the AAAI Conference on Artificial Intelligence*, volume 33, 2019.

18 Wided Medjroubi, Ulf Philipp Müller, Malte Scharf, Carsten Matke, and David Kleinhans. Open data in power grid modelling: New approaches towards transparent grid models. *Energy Reports*, 3, 2017.

19 Tian Sang, Paul Beame, and Henry Kautz. Heuristics for fast exact model counting. In *Theory and Applications of Satisfiability Testing: 8th International Conference, SAT 2005, St Andrews, UK, June 19-23, 2005. Proceedings 8*. Springer, 2005.

20 Tian Sang, Paul Beame, and Henry Kautz. Solving Bayesian networks by weighted model counting. In *Proceedings of the Twentieth National Conference on Artificial Intelligence (AAAI-05)*, volume 1. AAAI Press, 2005.

21 Marco Scutari. Learning Bayesian networks with the bnlearn R package. *arXiv preprint arXiv:0908.3817*, 2009. *arXiv:0908.3817*.

22 Shubham Sharma, Subhajit Roy, Mate Soos, and Kuldeep S. Meel. GANAK: A Scalable Probabilistic Exact Model Counter. In *IJCAI*, volume 19, 2019.

23 Mate Soos, Stephan Gocht, and Kuldeep S. Meel. Tinted, detached, and lazy cnf-xor solving and its applications to counting and sampling. In *International Conference on Computer Aided Verification*, 2020.

24 Mate Soos and Kuldeep S. Meel. Bird: engineering an efficient cnf-xor sat solver and its applications to approximate model counting. In *AAAI*, 2019.

25 Ryosuke Suzuki, Kenji Hashimoto, and Masahiko Sakai. Improvement of projected model-counting solver with component decomposition using SAT solving in components. Technical report, JSAI Technical Report, SIG-FPAI-506-07, 2017.

26 Marc Thurley. sharpSAT-counting models with advanced component caching and implicit BCP. *SAT*, 4121, 2006.

27 Jonas Vlasselaer, Angelika Kimmig, Anton Dries, Wannes Meert, and Luc De Raedt. Knowledge compilation and weighted model counting for inference in probabilistic logic programs. In *Proceedings of the First Workshop on Beyond NP*. AAAI Press, 2016.

28 Bart Wiegmans. Gridkit: European And North-American Extracts, March 2016. doi: 10.5281/ZENODO.47317.
<figure><img src="image_140.png" /></figure>
# A CP Approach for the Liner Shipping Network Design Problem

Yousra El Ghazi [ID

Aix Marseille Univ, Université de Toulon, CNRS, LIS, Marseille, France

Djamal Habet [ID

Aix Marseille Univ, Université de Toulon, CNRS, LIS, Marseille, France

Cyril Terrioux [ID

Aix Marseille Univ, Université de Toulon, CNRS, LIS, Marseille, France

## Abstract

The liner shipping network design problem consists, for a shipowner, in determining, on the one hand, which maritime lines (in the form of rotations serving a set of ports) to open, and, on the other hand, the assignment of ships (container ships) with the adapted sizes for the different lines to carry all the container flows. In this paper, we propose a modeling of this problem using constraint programming. Then, we present a preliminary study of its solving using a state-of-the-art solver, namely the OR-Tools CP-SAT solver.

**2012 ACM Subject Classification** Computing methodologies → Artificial intelligence

**Keywords and phrases** Constraint optimization problem, modeling, solving, industrial application

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

**Supplementary Material** Dataset: https://pageperso.lis-lab.fr/cyril.terrioux/LSNDP/instances.zip

**Funding** This work has been funded by Bpifrance under the PIA project Transformation Numérique du Transport Maritime (TNTM).

**Acknowledgements** The authors would like to thank the anonymous reviewers for their valuable comments and suggestions.

## 1 Introduction

Nowadays, maritime transport plays a major role in world trade. According to the International Maritime Organization (IMO), more than 80% of international trade is carried out by sea. The transport of containerized commodities constitutes the major part of this trade. It relies on more than 5,000 container ships that serve more than 500 ports worldwide. In this context, many combinatorial optimization problems [6, 9, 28] may arise with non-negligible economic and ecological impacts given their scale.

In this paper, we focus on the Liner Shipping Network Design Problem (LSNDP [9]). A shipping line, also called a service, is defined by a cyclic route (called a rotation) that visits a given set of ports in a given order and at regular times (see, for example, Figure 1). Generally, each port is thus visited by a vessel of the line at a weekly or biweekly frequency. All the vessels on a line are assumed to be homogeneous in terms of their main features (loading capacity, speed, fuel consumption, engines, ...). Operating a weekly line with a rotation lasting *k* weeks requires *k* vessels of the same type. Given a set of ports, a fleet of container ships, and a container flow (defined by a set of triples consisting of the original port of the commodities, their destination port, and the number of containers they represent), the LSNDP problem consists, for a shipowner, in determining, on the one hand, which shipping lines to open, and, on the other hand, which ships to operate on each line in order to carry as

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

© Yousra El Ghazi, Djamal Habet, and Cyril Terrioux;
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. 16; pp. 16:1–16:21

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

LIPICS

Leibniz International Proceedings in Informatics

Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl Publishing, Germany
<figure><img src="image_143.png" /><figcaption>Figure 1 Example of a line connecting Asia and Europe.</figcaption></figure>

many commodities as possible while ensuring a weekly frequency of visits to each port and optimizing costs. It is classified as NP-hard [7]. To give an idea of its difficulty, we can note that, taken separately, each of its two subproblems is already an NP-hard problem [7, 13]. Moreover, from a practical viewpoint, its solving by exact methods is currently limited to instances with a dozen ports at most. Although recent, this problem has been the subject of many works, especially in the last ten years. Most of them are from operational research. Note that, in the literature, different variants of the LSNDP problem are studied according to the assumptions and properties taken into account (transit time, transshipment, constant or variable speed from one rotation to another or from one leg (i.e. a trip between two consecutive ports in the rotation) to another, type of service, possibility of refusing some commodities, ...).

Different approaches have been considered, often based on integer (mixed) linear programming (e.g. [18, 19, 21, 22, 23, 27, 30]). They are mainly based on two types of formulations. The first type of formulation is service-oriented. The set of possible services is calculated upstream and provided as input to the model. The latter is then limited to selecting the services to be kept among the candidates. The main disadvantage of this type of formulation is that the number of possible services grows exponentially with the number of ports, which limits its practical interest in the context of a solving performed with complete methods. On the other hand, it can be interesting in the context of incomplete methods, because one can then consider only a subset of the possible services. In practice, the solution proposed in [2] and based on a tabu method coupled with column generation allowed handling instances up to 120 ports. The second type of formulation is based on the selection of the arcs of the graph representing the possible links between each pair of ports. A service is then defined by the arcs that compose it, and the same arc can be used to define several services. From a practical viewpoint, such modeling coupled with complete methods [18, 23, 22] allows handling instances with up to a dozen ports [9].

Other approaches (e.g. [1, 2]) are based on two-step solving. Since the LSNDP problem consists of two subproblems, they process each subproblem separately. For example, the approaches presented in [2, 5, 13, 12] solve, in the first step, the problem of creating services and, in the second step, consider the vessel assignment and the management of the commodity flow based on the services found by the first step. In [17], the first phase is devoted to the management of the flow, while the second defines the services. Generally, the solving is done in several passes, the first phase then benefiting from some feedbacks from the second phase of the previous pass. Of course, this type of approach corresponds to incomplete methods. In practice, these approaches provide satisfactory results for instances with up to 120 ports [2]. They are often based on matheuristics (e.g. [5, 13, 12]) or the Variable Neighborhood Search algorithm (VNS) such as [17].
Beyond that, there are many related problems to the LSNDP problem. For example, the Vehicle Routing Problem (VRP [15]) and its variants have strong similarities with the LSNDP problem. In particular, the routes are circuits and, for some variants, the vehicle load or transit times can be taken into account. In maritime transport, the Liner Shipping Fleet Repositioning Problem (LSFRP [28]) consists in moving container ships from one service to another while taking into account the commodities to be transported, the empty containers to be relocated and maximizing the difference between the revenues and the costs generated. Among the approaches studied to solve this problem, we can underline the interest in using constraint programming (CP) put forward in [14].

While the VRP and LSFRP problems (and other shipping-related problems such as [16, 24]) have been studied from a CP perspective, this does not seem to be the case for the LSNDP problem. In this paper, we propose a model to handle a relatively general version of the LSNDP problem. Our model considers variable speeds from one trip to another and takes into account transshipments and transit times. Although developed in partnership with one of the world's leading container shipping companies, the model presented takes into account a relatively general version of the LSDNP problem. It can, of course, be adapted to specific needs, taking advantage of the flexibility of constraint programming. One of the aims of this work is to study the interest of a CP-based approach to modeling and solving such problems.

This paper is organized as follows. Section 2 introduces the notions necessary to understand the paper. Then, in Section 3, we propose a CP model for the LSNDP problem. Finally, we present some experimental results, in Section 4, before discussing related work in Section 5 and concluding in Section 6.

# 2 Preliminary

## 2.1 Constraint Programming

An instance $P$ of the Constraint Optimization Problem (COP) can be defined as a 4-tuple $(X, D, C, f)$ where $X = \{x_1, \dots, x_n\}$ is the set of variables, $D = \{D_{x_1}, \dots, D_{x_n}\}$ is the set of domains, the domain $D_{x_i}$ being related to the variable $x_i$, $C = \{c_1, \dots, c_e\}$ represents the set of the constraints which define the interactions between the variables and describe the allowed combinations of values and $f$ specifies the criterion to be optimized. Solving a COP instance $P = (X, D, C, f)$ amounts to finding an assignment of all variables of $X$ satisfying all constraints of $C$ and optimizing the criterion given by $f$. This is an NP-hard problem.

One of the advantages of constraint programming lies in the existence of specialized constraints (the global constraints) which will make easier the modeling of problems, but also, their solving thanks to their dedicated filtering algorithms. In the following, we will exploit the following global constraints (where $\odot$ denotes a relational operator among $\leq$, $<$, $=, \neq, >$ or $\geq$):

- **Alldiff-except**($Y, v$) [3, 10] which ensures that the values of the variables of $Y$ are pairwise distinct, except in the case where they are equal to $v$,

- **Circuit**($Y, \ell$) [4] which imposes that the values of the variables of $Y$ form a circuit of length $\ell$ (in the sense of graph theory), each variable $y_i$ having for value $i$ if it does not take part in the circuit, and $j$ (with $i \neq j$) if $j$ is the successor of $i$ in the circuit,

- **Count**($Y, V$) $\odot$ $k$ [4, 8] which ensures that the number of variables of $Y$ whose value belongs to $V$ satisfies the condition imposed by the relation $\odot$ with respect to $k$,

- **Elt**($Y, i$) = $k$ [29] which ensures that the $i$th value of $Y$ (using a 0-based indexing) is equal to $k$ ($Y$ can be here an ordered set of variables or values),
16:4 A CP Approach for the Liner Shipping Network Design Problem

- $\text{Elt}_m(Y, i, j) = k$ [4, 3] which ensures the same property as $\text{Elt}$, but, for an ordered set $Y$ of variables or values organized in the form of a two-dimensional matrix,

- $\text{Maximum}(Y) = k$ which ensures that the greatest value of $Y$ is equal to $k$ ($Y$ can here be a set of variables or expressions),

- $\text{Sum}(Y, \Lambda) \odot k$ which imposes that the sum of the values of $Y$ weighted by the coefficients of $\Lambda$ satisfies the condition imposed by the relation $\odot$ with respect to $k$. In the following, this constraint will be represented in the more explicit form $\sum_i \lambda_i \cdot y_i \odot k$.

2.2 The Liner Shipping Network Design Problem

Liner shipping involves the use of standardized vessels that will reliably move cargo between ports according to a pre-determined route and schedule. It is often compared to scheduled passenger service, such as a train or bus service because it operates on a fixed schedule and provides regular and predictable service for shippers and receivers of commodities. A shipping line, also called a service, is defined by a cyclic route (called a rotation) that serves a set of ports in a specific order and on a regular schedule. Figure 1 describes the example of a line connecting Asia and Europe.

In this paper, we consider only the transportation of commodities in their containerized form, as this mode of transportation constitutes the bulk of freight transportation in terms of quantity and value. Thus, a customer wishing to move commodities from a port of loading (POL) to a port of destination (POD) needs to place them inside one or more containers. Containers have the advantage, for the shipowner, that their dimensions are standardized, thus facilitating their handling and placement on board of specialized vessels such as container ships. There are mainly two sizes of containers: 20-foot containers (about 6.1 m) and 40-foot containers (12.2 m) with a height of 8.6 feet (2.6 m) and a width of 8 feet (2.4 m). The majority of containers transported are of one of these two sizes. Also, the storage space of the vessels is divided into 40-foot unit spaces on which it is possible to stack both 40-foot and 20-foot units. The Twenty-foot Equivalent Unit (TEU) is the unit generally used to count a number of containers. For example, a 40-foot container counts as 2 TEUs.

From the carrier's viewpoint, each commodity $k$ is seen as a quantity $q(k)$ of containers (expressed in TEUs) to be transported from the port of origin $pol(k)$ to the port of destination $pod(k)$ in exchange for a revenue $rev(k)$ per TEU (expressed in dollars). This revenue may be zero in the case of empty containers. Some commodities may have a maximal transit time $tt_{max}(k)$ that must be respected. This time is the maximum time allowed for their transport. Generally, such commodities are transported within the framework of premium offers proposed by the carriers. It should be noted that a batch of containers sent by a customer from one port to another cannot be divided into several sub-batches. Finally, cargo can be transported from its port of origin to its port of destination via the successive use of different lines. The operation of unloading a commodity from one line and loading it on another is called "transshipment". It may require the commodities to be stored for several days in the transshipment port until the vessel of the next line arrives and loads them on board. This can result in costs (see transshipment costs below) and longer travel times.

Concerning vessels, container ships are grouped by type of vessel with identical or similar features. Thus, each class $v$ is characterized by its capacity $\kappa(v)$ (i.e. the maximum number of containers (in TEU) that can be transported), its daily charter rate $tc(v)$ (corresponding to the daily cost of using the vessel), its interval of possible speeds [$\nu_{min}(v), \nu_{max}(v)$] (in knots), its hourly consumption $\text{cons}(v, \nu)$ of fuel for the main engine (in tons per hour), for each type of fuel $\text{fuel}(v)$ and each possible speed $\nu$. Regarding consumption, other parameters
that could have an impact such as wind strength, sea currents, draft or load on board are ignored. These parameters can be variable in time and difficult to anticipate, as the lines are generally defined on a yearly scale.

Each port $p$ also has its own features, namely its productivity $prod(p, v)$ (i.e., the number of containers loaded or unloaded per hour for vessels of type $v$), its waiting time $wt(p, v)$ (time at anchor before entering the port), its maneuvering time $man^{in}(p, v)$ and $man^{out}(p, v)$ to enter and leave the port respectively, its call charges $pc(p, v)$ (in dollars), and its transshipment cost $ts(p)$ (in dollars). The times are given in hours and depend on the type $v$ of vessels, as do the call charges. A canal $c$ (e.g. the Suez or Panama Canal) is characterized by a waiting time $wt(c, v)$, a traversal time $trav(c)$ (in hours), and a traversal cost $pc(c, v)$ (in dollars).

The number of vessels operating on a service is determined by the length of the rotation and the frequency of departures. Indeed, a rotation must guarantee a regular frequency of visits to each port it serves. This rotation frequency is generally weekly or biweekly. For a weekly frequency, the duration of the rotation must be a multiple of seven days. The number of vessels deployed per rotation must then be equal to the number of weeks in that rotation. For example, the line shown in Figure 1 has a duration of 91 days or 13 weeks. It is therefore operated with 13 vessels.

The liner shipping network design problem (LSNDP) can be defined as follows: Given a set of ports, a set of vessels divided by type (each type $v$ having $nb(v)$ vessels) and a set of commodities to be transported, define a set of rotations having a weekly frequency and determine the vessels operating them to transport the commodities while respecting, if necessary, the transit times and maximizing the profit. The profit is defined as the difference between the revenues generated by the commodities transported and all the costs generated by this transport (fuel costs, vessel operating costs, port call and canal costs, transshipment costs, ...). To calculate fuel costs, for each type of fuel $f$, we have the price $fp(f)$ per ton of fuel (expressed in dollars).

While the primary purpose of this problem is to design maritime transportation networks, it can also be used to assist in decision-making. For example, it can be used to simulate situations such as traffic jams to enter certain ports and determine whether or not it is relevant to adapt existing rotations. It can also be used to consider changes in the flow of containers to be transported, to evaluate the interest in taking market share in certain commodity flows or to anticipate the construction of new ships.

# 3 Model

## 3.1 Modeling Choices

In our model, we adopt the usual assumptions of the literature. In particular, we assume that all container ships of a given type have identical features and that the frequency of services is weekly. Furthermore, we choose to treat canals (such as the Suez and Panama canals) in the same way as ports. The time it takes to cross a canal replaces the time it takes to load/unload a ship in a port. As a result, the notion of rotation now also takes canals into account. Since a rotation can pass several times through the same canal, but not through the same port, we consider, in our model, two instances of each canal so that a canal can be used both on the “outbound” and on the “return”. For example, we can see that the line represented in Figure 1 passes twice through the Suez Canal, once on the “outbound” (blue route) and once on the “return” (red route). Note that creating more than two instances of the same canal is of little interest because a solution passing more than twice through the same canal has little chance of being optimal.
16:6 A CP Approach for the Liner Shipping Network Design Problem

<figure><img src="image_144.png" /><figcaption>Figure 2 A rotation involving five ports.</figcaption></figure>

$$ \text{Circuit}(\{s_p^r | p \in \mathcal{P} \cup \mathcal{C}, \ell_r\}) \qquad r \in \mathcal{R} \qquad (R.1) $$

$$ v_r > 0 \Rightarrow \sum_{p \in \mathcal{P}} (s_p^r \neq p) \geq 3 \qquad r \in \mathcal{R} \qquad (R.2) $$

$$ v_r > 0 \iff \ell_r \geq 3 \qquad r \in \mathcal{R} \qquad (R.3) $$

<caption>Figure 3 Constraints related to rotations and routes.</caption>

Our model takes, as inputs, all the information about the ports and canals, the flow of commodities, the types of vessels, and the distances between pairs of ports/canals. It also relies on the maximum number $r_{max}$ of rotations to be defined, the time horizon $h_{max}$ (in hours), that is the maximum duration to achieve a rotation, and the maximum number $ts_{max}$ of transshipments allowed per commodity. One of the particularities of this model is that the main operations will be time-stamped to handle rotation or transit times of the commodities as precisely as possible. Moreover, speeds can be different from one leg to another. Finally, it takes into account the possibility of refusing to transport a commodity in the network if it is not profitable or impossible.

Thereafter, given the large number of variables, we define the variables progressively when needed. The set of ports is denoted $\mathcal{P} = \{0, 1, \dots, |\mathcal{P}| - 1\}$, the set of canals $\mathcal{C} = \{c, c+|\mathcal{C}_0|, \dots, |\mathcal{P}|, \dots, |\mathcal{C}_0| - 1\}$ the set of canals before duplication), the set of type of vessels available $\mathcal{V} = \{1, 2, \dots, |\mathcal{V}|\}$ and that of the commodities $\mathcal{K} = \{1, 2, \dots, |\mathcal{K}|\}$. We note respectively $\mathcal{I}$ and $\mathcal{I}^+$ the index sets $[0, ts_{max}]$ and $[0, ts_{max} + 1]$. Let $\mathcal{R} = \{1, \dots, r_{max}\}$ be the set of the indices of the possible rotations.

## 3.2 Definition of Rotations

Our model does not necessarily use all possible $r_{max}$ rotations. Therefore, we consider a variable $v_r$ per rotation $r$. Its value is an integer between 1 and $|\mathcal{V}|$ representing the type of vessels exploited if the rotation is used, 0 otherwise. Each rotation $r$ must correspond to a circuit. To define such circuits, we introduce a variable $s_p^r$ per port/canal $p$ and rotation $r$. Its value is $p$ if the port/canal $p$ is not involved in the circuit, the successor of port/canal $p$ otherwise. Figure 2 illustrates this for a rotation involving five ports. We also introduce a variable $\ell_r$ specifying the length of the circuit associated with rotation $r$. Thus, for each rotation $r$, the existence of a circuit can be guaranteed by constraint (R.1) (see Figure 3). Note that this constraint avoids the existence of subtours, a property that is not generally easy to guarantee. For instance, in MIP, avoiding subtours requires adding non-linear constraints that must then be linearized. Then, thanks to constraint (R.2), a circuit must involve at least three ports (this is a business rule generally desired by carriers), and so, cannot involve only canals. Finally, constraint (R.3) ensures that rotation $r$ is used if and only if the associated circuit has a length at least equal to three.
<figure><img src="image_145.png" /><figcaption>Figure 4 Transport of a commodity *k* from its origin port *pol*(*k*) to its destination port *pod*(*k*).</figcaption></figure>

## 3.3 Cargo Flow

Our model allows for the possibility of not carrying a commodity *k* if it is not possible or not profitable. To do this, we define a Boolean variable $\alpha_k$ that is true if commodity *k* is accepted in the network, false otherwise. Taking charge of commodity *k* means that it is loaded at the original port *pol*(*k*) and unloaded at the destination one *pod*(*k*), possibly passing through intermediate ports. In our model, we consider only the intermediate ports where the commodity will be transshipped as shown in Figure 4. Also, we introduce a variable $rot_{k,i}$ per commodity *k* and step *i* to represent the *i*th rotation used to transport the commodity *k*. $rot_{k,i}$ has the value *r* if commodity *k* is carried thanks to rotation *r* during the *i*th step, -1 if this step is not needed. Similarly, the variable $port_{k,i}$ represents the port where the commodity *k* enters its *i*th rotation. By so doing, commodity *k* enters in its *i*th rotation at port $port_{k,i}$ and leaves it at port $port_{k,i+1}$ (i.e. the port in which it starts its next rotation if $port_{k,i+1}$ differs from $pod(k)$). These variables have the value of the corresponding port *p* ($p \in \mathcal{P}$) if the *i*th rotation is used, -1 otherwise. More precisely, the domain of $port_{k,i}$ is $\{-1, pol(k)\}$ if $i = 0$, $\{-1, pod(k)\}$ if $i = ts_{max} + 1$, $\{-1\} \cup \mathcal{P} - \{pol(k)\}$ otherwise. For each commodity *k*, a variable *nbr*<sub>*k*</sub> specifies the number of rotations used (between 0 and $ts_{max} + 1$).

We can now define the associated constraints (see Figure 5). First, constraint (F.1) specifies that commodity *k* is accepted in the network if and only if there is at least one rotation that carries it. Of course, the port of departure of an accepted commodity *k* is necessarily its port of origin *pol*(*k*) (constraint (F.2)). The last port used is necessarily the destination port *pod*(*k*). To ensure this, we introduce a variable *pod*<sub>*k*</sub> per commodity *k* that can take two values: either *pod*(*k*) if commodity *k* is accepted, or -1 otherwise. Constraint (F.3) guarantees that *pod*<sub>*k*</sub> has the relevant value depending on the value of $\alpha_k$ while constraint (F.4) ensures that the last used port is consistent with *pod*<sub>*k*</sub>. Naturally, no port or rotation can be used beyond the destination port (constraints (F.5)-(F.7)). Moreover, a commodity cannot transit several times through the same port or the same rotation, what is ensured by constraints (F.8) and (F.9). To facilitate the expression of some constraints about the path followed by the commodities, we introduce Boolean variables $from_{k,p}$ (resp. $to_{k,p}$) which are true if commodity *k* enters (resp. leaves) the rotation *r* at port *p*, false otherwise. These variables are directly related to the previous ones as seen in constraints (F.10) and (F.11). These variables are also used to link the cargo flow to the definition of routes. Indeed, if commodity *k* is (un)loaded at a port *p* for a rotation *r*, it implies that the port *p* is used in this rotation (constraints (F.12) and (F.13)). Conversely, if a port *p* is used in rotation *r*, then there is at least one commodity that is (un)loaded in that port for that rotation (see constraints (F.14) and (F.15)). Finally, if a rotation is not used, no commodities can transit through it, and vice versa (constraint (F.16)). Note that this constraint is redundant, but in most cases, it allows finding some conflicts earlier.

## 3.4 Properties of Rotations and Vessels

In some MIP models in the literature (e.g. [18, 9]), each service is associated with a predefined type of vessel. While this choice facilitates the consideration of the specificities of each type of vessel, it leads to handling many rotations, few of which will be used in the end. Moreover,
16:8 A CP Approach for the Liner Shipping Network Design Problem

$$ \alpha_k = 1 \iff nbr_k > 0 \qquad k \in \mathcal{K} \qquad (\text{F.1}) $$

$$ \alpha_k = 1 \iff port_{k,1} = pol(k) \qquad k \in \mathcal{K} \qquad (\text{F.2}) $$

$$ \alpha_k = 1 \iff pod_k = pod(k) \qquad k \in \mathcal{K} \qquad (\text{F.3}) $$

$$ \text{Elt}(\{port_{k,i+1} | i \in [1, ts_{max}]\}, nbr_k) = pod_k \qquad k \in \mathcal{K} \qquad (\text{F.4}) $$

$$ i \ge nbr_k \iff rot_{k,i} = -1 \qquad k \in \mathcal{K}, i \in \mathcal{I} \qquad (\text{F.5}) $$

$$ nbr_k < i \Rightarrow port_{k,i} = -1 \qquad k \in \mathcal{K}, i \in \mathcal{I}^{+} \qquad (\text{F.6}) $$

$$ port_{k,i} = -1 \Rightarrow nbr_k \le i \qquad k \in \mathcal{K}, i \in \mathcal{I}^{+} \qquad (\text{F.7}) $$

$$ \text{Alldiff-except}(\{port_{k,i} | i \in \mathcal{I}^{+}\}, -1) \qquad k \in \mathcal{K} \qquad (\text{F.8}) $$

$$ \text{Alldiff-except}(\{rot_{k,i} | i \in \mathcal{I}\}, -1) \qquad k \in \mathcal{K} \qquad (\text{F.9}) $$

$$ from_{k,p}^{r} = \sum_{i \in \mathcal{I}} (port_{k,i} = p) \cdot (rot_{k,i} = r) \qquad k \in \mathcal{K}, p \in \mathcal{P}, r \in \mathcal{R} \qquad (\text{F.10}) $$

$$ to_{k,p}^{r} = \sum_{i \in \mathcal{I}} (port_{k,i+1} = p) \cdot (rot_{k,i} = r) \qquad k \in \mathcal{K}, p \in \mathcal{P}, r \in \mathcal{R} \qquad (\text{F.11}) $$

$$ from_{k,p}^{r} = 1 \Rightarrow s_{p}^{r} \neq p \qquad k \in \mathcal{K}, p \in \mathcal{P}, r \in \mathcal{R} \qquad (\text{F.12}) $$

$$ to_{k,p}^{r} = 1 \Rightarrow s_{p}^{r} \neq p \qquad k \in \mathcal{K}, p \in \mathcal{P}, r \in \mathcal{R} \qquad (\text{F.13}) $$

$$ s_{p}^{r} \neq p \Rightarrow \text{Count}(\{port_{k,i} | k \in \mathcal{K}, i \in \mathcal{I}^{+}\}, \{p\}) \ge 1 \qquad p \in \mathcal{P}, r \in \mathcal{R} \qquad (\text{F.14}) $$

$$ s_{p}^{r} \neq p \Rightarrow \text{Count}(\{rot_{k,i} | k \in \mathcal{K}, i \in \mathcal{I}\}, \{r\}) \ge 1 \qquad p \in \mathcal{P}, r \in \mathcal{R} \qquad (\text{F.15}) $$

$$ v_r = 0 \iff \text{Count}(\{rot_{k,i} | k \in \mathcal{K}, i \in \mathcal{I}^{+}\}, \{r\}) = 0 \qquad r \in \mathcal{R} \qquad (\text{F.16}) $$

<caption>Figure 5 Constraints related to cargo flow.</caption>

trying a new vessel type for a rotation requires the solver to change the assignment of a lot of variables. In our model, we choose to let the solver decide on the type of vessels associated with each rotation. Therefore, it is necessary to ensure that the type of vessels chosen for a rotation matches the features of the rotation. This requires the introduction of a certain number of variables whose values will then be fixed using Elt constraints. The variable $\kappa_r$ represents the maximum capacity (expressed in TEUs) of commodities that can be transported via rotation $r$. Its value is 0 if the rotation is not used, the capacity of the type of vessel used otherwise. The variables $v_{min}^r$ and $v_{max}^r$ specify the minimum and maximum speeds of the rotation $r$ (0 if the rotation is not used). The variable $fp_r$ expresses the price per ton of fuel for rotation $r$ (0 if the rotation is not operated). For each rotation $r$, we post the constraints (P.1)–(P.4) (see Figure 6). Likewise, some information (call costs, waiting time, ...) related to ports or canals also depends on the type of vessels associated with rotation $r$. For each rotation $r$, we then introduce the variables $wt_{p,r}$, $man_{p,r}^{in}$ and $man_{p,r}^{out}$ which specify respectively the waiting time of port/canal $p$ and the maneuvering time to enter and leave port $p$. The cost of calling at port $p$ for rotation $r$ is represented by the variable $pc_{p,r}$ while the productivity for port $p$ and rotation $r$ is expressed by the variable $prod_{p,r}$. Constraints (P.5)–(P.9) ensure the consistency of these features. Finally, we consider the variable $tc_r$ which, for each rotation $r$, specifies the daily cost of using the vessels associated with the rotation and its related constraint (P.10).
$$ \mathrm{Elt}(\{0\} \cup \{\kappa(v) \mid v \in \mathcal{V}\}, v_r) = \kappa_r \qquad r \in \mathcal{R} \qquad (\mathrm{P}.1) $$

$$ \mathrm{Elt}(\{0\} \cup \{\nu_{\min}(v) \mid v \in \mathcal{V}\}, v_r) = \nu_{\min}^r \qquad r \in \mathcal{R} \qquad (\mathrm{P}.2) $$

$$ \mathrm{Elt}(\{0\} \cup \{\nu_{\max}(v) \mid v \in \mathcal{V}\}, v_r) = \nu_{\max}^r \qquad r \in \mathcal{R} \qquad (\mathrm{P}.3) $$

$$ \mathrm{Elt}(\{0\} \cup \{fp(fuel(v)) \mid v \in \mathcal{V}\}, v_r) = fp_r \qquad r \in \mathcal{R} \qquad (\mathrm{P}.4) $$

$$ \mathrm{Elt}(\{0\} \cup \{wt(p, v) \mid v \in \mathcal{V}\}, v_r) = wt_{p,r} \qquad r \in \mathcal{R}, p \in \mathcal{P} \cup \mathcal{C} \qquad (\mathrm{P}.5) $$

$$ \mathrm{Elt}(\{0\} \cup \{man^{in}(p, v) \mid v \in \mathcal{V}\}, v_r) = man_{p,r}^{in} \qquad r \in \mathcal{R}, p \in \mathcal{P} \cup \mathcal{C} \qquad (\mathrm{P}.6) $$

$$ \mathrm{Elt}(\{0\} \cup \{man^{out}(p, v) \mid v \in \mathcal{V}\}, v_r) = man_{p,r}^{out} \qquad r \in \mathcal{R}, p \in \mathcal{P} \cup \mathcal{C} \qquad (\mathrm{P}.7) $$

$$ \mathrm{Elt}(\{0\} \cup \{pc(p, v) \mid v \in \mathcal{V}\}, v_r) = pc_{p,r} \qquad r \in \mathcal{R}, p \in \mathcal{P} \cup \mathcal{C} \qquad (\mathrm{P}.8) $$

$$ \mathrm{Elt}(\{1\} \cup \{prod(p, v) \mid v \in \mathcal{V}\}, v_r) = prod_{p,r} \qquad r \in \mathcal{R}, p \in \mathcal{P} \cup \mathcal{C} \qquad (\mathrm{P}.9) $$

$$ \mathrm{Elt}(\{0\} \cup \{tc(v) \mid v \in \mathcal{V}\}, v_r) = tc_r \qquad r \in \mathcal{R} \qquad (\mathrm{P}.10) $$

<caption>Figure 6 Constraint related to properties of rotations and vessels.</caption>

$$ from_{k,p}^r = 1 \Rightarrow leave_{k,p}^r = 1 \qquad k \in \mathcal{K}, r \in \mathcal{R}, p \in \mathcal{P} \qquad (\mathrm{L}.1) $$

$$ to_{k,p}^r = 1 \Rightarrow leave_{k,p}^r = 0 \qquad k \in \mathcal{K}, r \in \mathcal{R}, p \in \mathcal{P} \qquad (\mathrm{L}.2) $$

$$ (s_p^r = p' \land from_{k,p'}^r = 0 \land to_{k,p'}^r = 0) \Rightarrow leave_{k,p}^r = leave_{k,p'}^r \qquad k \in \mathcal{K}, r \in \mathcal{R}, p, p' \in \mathcal{P} \qquad (\mathrm{L}.3) $$

$$ s_p^r = p \Rightarrow leave_{k,p}^r = 0 \qquad k \in \mathcal{K}, r \in \mathcal{R}, p \in \mathcal{P} \qquad (\mathrm{L}.4) $$

$$ \sum_{k \in \mathcal{K}} q(k) \cdot leave_{k,p}^r \le \kappa_r \qquad r \in \mathcal{R}, p \in \mathcal{P} \qquad (\mathrm{L}.5) $$

<caption>Figure 7 Constraints related to loads.</caption>

## 3.5 Vessel Load

We need to ensure that vessels do not leave each port loaded beyond their maximum capacity. This requires knowing, for each rotation, which commodities it carries at the exit of each port. To do this, we use a Boolean variable $leave_{k,p}^r$ per commodity $k$, port $p$, and rotation $r$. This variable is true if commodity $k$ leaves port $p$ via rotation $r$, false otherwise. Constraints (L.1) and (L.2) (see Figure 7) deal with the cases when the commodities are respectively loaded in rotation $r$ and unloaded from rotation $r$ while constraint (L.3) guarantees the transitivity all along the trip. Finally, the constraint (L.4) corresponds to the case where a port $p$ does not appear in rotation $r$. Constraint (L.5) then allows ensuring that, for each port $p$, the load, when leaving the port, does not exceed the maximum capacity $\kappa_r$ of rotation $r$.

## 3.6 Timestamp and Transit Times

### 3.6.1 Duration of Port Operations and Canal Traversal

To express the duration of loading/unloading operations in a port or the duration of traversing a canal, we introduce a variable $t_p^r$ per port/canal $p$ and rotation $r$. In the case of a canal the value of this variable is defined as equal to the duration of the traversal if the canal is used, 0 otherwise (see constraint (T.1) in Figure 8). For a port $p$, two cases are possible. If the port
16:10 A CP Approach for the Liner Shipping Network Design Problem

$$t_p^r = \text{trav}(p) \cdot (s_p^r \neq p) \qquad r \in \mathcal{R}, p \in \mathcal{C} \qquad (\text{T.1})$$

$$s_p^r = p \iff t_p^r = 0 \qquad r \in \mathcal{R}, p \in \mathcal{P} \qquad (\text{T.2})$$

$$teu_{p,r} = \sum_{k \in \mathcal{K}} (from_{k,p}^r + to_{k,p}^r) \cdot q(k) \qquad r \in \mathcal{R}, p \in \mathcal{P} \qquad (\text{T.3})$$

$$s_p^r \neq p \Rightarrow t_p^r = \left\lceil \frac{\mu.teu_{p,r}}{prod_{p,r}} \right\rceil \qquad r \in \mathcal{R}, p \in \mathcal{P} \qquad (\text{T.4})$$

$$v_r = 0 \iff dep_r = -1 \qquad r \in \mathcal{R} \qquad (\text{T.5})$$

$$v_r > 0 \Rightarrow \text{Maximum}(\{p.(s_p^r \neq p) | p \in \mathcal{P}\}) = dep_r \qquad r \in \mathcal{R} \qquad (\text{T.6})$$

$$dep_r = p \Rightarrow time_{p,r}^{in} = 0 \qquad r \in \mathcal{R}, p \in \mathcal{P} \qquad (\text{T.7})$$

$$time_{p,r}^{out} = time_{p,r}^{in} + t_p^r \qquad r \in \mathcal{R}, p \in \mathcal{P} \cup \mathcal{C} \qquad (\text{T.8})$$

$$s_p^r = p' \Rightarrow st_p^r = \left\lceil \frac{\delta(p, p')}{\nu_p^r} \right\rceil \qquad r \in \mathcal{R}, p, p' \in \mathcal{P} \cup \mathcal{C} \qquad (\text{T.9})$$

$$(v_r > 0 \land s_p^r \neq p) \Rightarrow \nu_p^r \geq \nu_{min}^r \qquad r \in \mathcal{R}, p \in \mathcal{P} \cup \mathcal{C} \qquad (\text{T.10})$$

$$v_r > 0 \Rightarrow \nu_p^r \leq \nu_{max}^r \qquad r \in \mathcal{R}, p \in \mathcal{P} \cup \mathcal{C} \qquad (\text{T.11})$$

$$s_p^r = p \iff st_p^r = 0 \qquad r \in \mathcal{R}, p \in \mathcal{P} \cup \mathcal{C} \qquad (\text{T.12})$$

$$s_p^r = p \Rightarrow time_{p,r}^{in} = 0 \qquad r \in \mathcal{R}, p \in \mathcal{P} \cup \mathcal{C} \qquad (\text{T.13})$$

$$s_p^r = p \iff \nu_p^r = 0 \qquad r \in \mathcal{R}, p \in \mathcal{P} \cup \mathcal{C} \qquad (\text{T.14})$$

$$(s_p^r = p' \land p' \neq dep_r) \Rightarrow time_{p',r}^{in} = time_{p,r}^{out} + man_{p,r}^{out} + st_p^r + wt_{p',r} + man_{p',r}^{in} \qquad r \in \mathcal{R}, p, p' \in \mathcal{P} \cup \mathcal{C} \qquad (\text{T.15})$$

$$(from_{k,p}^r \land to_{k,p'}^r) \Rightarrow (time_{p,r}^{in} < time_{p',r}^{in} \lor time_{p',r}^{in} < time_{p,r}^{in} < time_{p',r}^{in} + T_r) \qquad r \in \mathcal{R}, p, p' \in \mathcal{P} \cup \mathcal{C}, k \in \mathcal{K} \qquad (\text{T.16})$$

<caption>Figure 8 Constraints related to timestamps.</caption>

is not used in the rotation *r*, the variable $t_p^r$ is 0 (see constraint (T.2)). Otherwise, its value depends on the number of TEUs loaded and unloaded in the port *p* for the rotation *r*. So, we consider the variable *teu*<sub>*p*,r*</sub> which indicates the number of TEUs loaded and unloaded at port *p* for rotation *r*. Its value can be defined thanks to constraint (T.3). We can express the duration of the operations thanks to constraint (T.4). A crane movement allows moving a container whatever its size. To take into account the existence of 20-foot and 40-foot containers among the commodities to be handled, the parameter $\mu$ makes it possible to calculate the number of containers to be handled and thus the number of crane movements necessary from the number of containers expressed in TEUs.

3.6.2 Call Timestamps

In order to establish the schedule for each call, it is necessary to designate a port as the departure port in each rotation. To do this, we consider a variable *dep*<sub>*r*</sub> per rotation *r* which has, for value, a port *p* if the rotation is used, -1 otherwise. The choice of the starting port being purely arbitrary, we choose the one with the largest index. This can be achieved thanks to constraints (T.5) and (T.6). Then, in our model, we consider two key moments: the moment when the vessel arrives at the berth (resp. enters the canal) and the moment when it leaves the berth (resp. leaves the canal). For each rotation *r* and each port/canal *p*, these two moments are represented respectively by the variables *time*<sub>*p*,r*</sub><sup>in</sup> and *time*<sub>*p*,r*</sub><sup>out</sup> which take their values in [0, *h*<sub>max</sub>]. For each rotation, we consider that time 0 coincides with the time of arrival at the berth in the departure port thanks to constraint (T.7). The time
of leaving a port or a canal depends only on the time of arrival and the duration of the operations in the port or the traversal of the canal (constraint (T.8)). Then, to determine the arrival time at a port or canal as a function of the exit time from the port/canal ahead of it in the rotation, we need to define the travel time as a function of the vessel's speed. The variables $st_p^r$ and $v_p^r$ represent respectively the travel time from the port/canal $p$ to its successor (if any) in rotation $r$ and the speed (expressed in knots) used on this leg. The two variables are correlated by constraint (T.9). Of course, the speeds used must be consistent with the capabilities of the vessels operating the rotation (constraints (T.10) and (T.11)). We can now define the time of arrival at the port/canal $p'$ from its predecessor $p$ in rotation $r$ thanks to constraint (T.15). Note that for canals, we consider that the variables $man_{p,r}^{in}$ and $man_{p,r}^{out}$ are 0. This allows us to avoid defining the previous constraint according to the different possibilities of port/canal successions. In the case where a port $p$ is not operated in a rotation $r$, we set the values of the variables $st_p^r$, $v_p^r$ and $time_{p,r}^{in}$ to the value 0 (constraints (T.12)-(T.14)).

Finally, if a commodity $k$ is loaded in rotation $r$ at port $p$ and unloaded at port $p'$, this imposes that the arrival at port $p$ takes place before the arrival at port $p'$ if the trip between $p$ and $p'$ does not pass through the departure port of rotation $r$. If this path passes through the departure port, then the arrival at port $p'$ will occur in the next rotation, and the arrival at port $p$ is between the two visits to port $p'$. For example, if we consider the rotation in Figure 2 (which lasts 840 hours) and a commodity sent from port $p_2$ to port $p_4$, a vessel operating this rotation enters port $p_2$ at hour 50 (in blue) and arrives in port $p_4$ at hour 600. In this case, we have $time_{p_2,r}^{in} < time_{p_4,r}^{in}$. On the other hand, if we consider a commodity going from $p_4$ to $p_3$, the vessel enters port $p_3$ at hour 275 before visiting $p_4$. This commodity will then be delivered only at the next passage of the vessel at time 1,115. We then have $time_{p_3,r}^{in} < time_{p_4,r}^{in} < time_{p_3,r}^{out} + T_r$. This is ensured by constraint (T.16).

### 3.6.3 Transit Times

To accurately consider the transit time of commodities, we need to know the key moments in their transportation, namely when they are loaded on board a rotation or unloaded. To simplify the model, we consider that a commodity is loaded on board a rotation when the rotation leaves the port and that it is unloaded when the rotation arrives at the port. These two times are represented by the variables $ctime_{i,k}^{in}$ and $ctime_{i,k}^{out}$ respectively. Constraints (T.17) and (T.18) (see Figure 9) ensure the correspondence between the key times of the rotations and the ones of the commodities.

The time spent by the commodity $k$ in its $i$th rotation is represented by the variable $\delta_{i,k}$. It corresponds naturally to the difference between the exit time and the entry time. However, we must take into account the particular case where the journey passes through the port of departure. In this case, the commodities are unloaded at the next rotation. For example, if we consider the previous example, a commodity sent from port $p_2$ to port $p_4$ leaves port $p_2$ at hour 60 (in red) and arrives in port $p_4$ at hour 600 (in blue). This gives a travel time of 540 hours. On the other hand, a commodity shipped from port $p_4$ to port $p_3$ leaves port $p_4$ at time 620 and arrives at port $p_3$ at time 1,115 and thus takes 495 hours to reach its destination. Constraints (T.19) and (T.20) deal respectively with the first and second cases.

In the case where a transshipment takes place, the time that the commodities spend on the quay between the two rotations must be taken into account. Given the weekly frequency of the rotations, this time can be of the order of a week at most. To consider it more precisely, we introduce a variable $\Delta_{i,k}$ per commodity $k$ and $i$th rotation used. The value of this variable is related to the weekly frequency of rotations. For example, consider a commodity
16:12 A CP Approach for the Liner Shipping Network Design Problem

$$
\mathrm{Elt}_{\mathrm{m}}(\{\mathrm{time}_{p,r}^{\mathrm{out}} \mid p \in \mathcal{P}, r \in \mathbb{R}\}, \mathrm{port}_{k,i}, \mathrm{rot}_{k,i}) = \mathrm{ctime}_{i,k}^{\mathrm{in}} \qquad k \in \mathcal{K}, i \in \mathcal{I} \qquad (\mathrm{T}.17)
$$

$$
\mathrm{Elt}_{\mathrm{m}}(\{\mathrm{time}_{p,r}^{\mathrm{in}} | p \in \mathcal{P}, r \in \mathbb{R}\}, \mathrm{port}_{k,i+1}, \mathrm{rot}_{k,i}) = \mathrm{ctime}_{i,k}^{\mathrm{out}} \quad k \in \mathcal{K}, i \in \mathcal{I} \quad (\mathrm{T}.18)
$$

$$
\mathrm{ctime}_{i,k}^{\mathrm{in}} \le \mathrm{ctime}_{i,k}^{\mathrm{out}} \Rightarrow \delta_{i,k} = \mathrm{ctime}_{i,k}^{\mathrm{out}} - \mathrm{ctime}_{i,k}^{\mathrm{in}} \qquad k \in \mathcal{K}, i \in \mathcal{I}^{+} \qquad (\mathrm{T}.19)
$$

$$
(\mathrm{rot}_{k,i} = r \land \mathrm{ctime}_{i,k}^{\mathrm{in}} > \mathrm{ctime}_{i,k}^{\mathrm{out}}) \Rightarrow \delta_{i,k} = \mathrm{ctime}_{i,k}^{\mathrm{out}} - \mathrm{ctime}_{i,k}^{\mathrm{in}} + T_r \quad k \in \mathcal{K}, i \in \mathcal{I}^{+} \tag{T.20}
$$

$$
i + 1 < nbr_k \Rightarrow \Delta_{i,k} = (ctime_{i+1,k}^{in} - ctime_{i,k}^{out})\%168 \qquad k \in \mathcal{K}, i \in \mathcal{I} \qquad (T.21)
$$

$$
i + 1 \geq nbr_k \Rightarrow \Delta_{i,k} = 0 \qquad k \in \mathcal{K}, i \in \mathcal{I} \qquad (T.22)
$$

$$
\sum_{i \in \mathcal{I}^{+}} \delta_{i,k} + \sum_{i \in \mathcal{I}} \Delta_{i,k} \leq tt_{max}(k) \qquad k \in \mathcal{K} \qquad (T.23)
$$

Figure 9 Constraints related to transit times.

k that arrives at a port *p* at hour 200 (according to its *ctime*<sub>*i*,*k*</sub><sup>*out*</sup> value) on a rotation *r* and
leaves it at hour 2,000 (according to its *ctime*<sub>*i+1*,*k*</sub><sup>*in*</sup> value) via a rotation *r'*. The weekly
frequency of the rotations *r* and *r'* implies that, in practice, commodity *k* leaves the port
at hour 320. Indeed, a vessel of rotation *r'* leaves the port at hours 152 (i.e. 2,000 modulo
(7 × 24)), 320, 488, ... If the arrival in the port is later than the departure from the port
(according to the values *ctime*<sub>*i*,*k*</sub><sup>*out*</sup> and *ctime*<sub>*i+1*,*k*</sub><sup>*in*</sup>), it means that the commodity must wait
for the vessel of the next week and so we have to add 168 hours to the considered difference.
Constraint (T.21) takes into account these two cases. The *i*th rotations that are not used for
commodity *k* have a *Δ*<sub>*i*,*k*</sub> variable whose value is zero (constraint (T.22)). The transit time
of commodity *k* can then be guaranteed by constraint (T.23). Note that not all commodities
have a transit time constraint. Also, if commodity *k* does not have a maximum transit time
constraint, the variables and constraints presented here are not considered for it.

3.7 Vessel Availability

Since the frequency of the rotations is weekly, each port is visited by one vessel operating the
rotation each week. The number of vessels needed is therefore the total time of the rotation
divided by the duration of one week. If the variables $n_r$ and $T_r$ represent respectively the
number of vessels used by the rotation $r$ and the total time of the rotation $r$, we can impose
constraint (A.1) (see Figure 10). The time of the rotation is, of course, zero if the rotation
is not used (constraint (A.2)). Otherwise, since each rotation starts at time 0, the total
time is the arrival time at the departure port from the last port of the rotation (constraint
(A.3)). Finally, we guarantee that, for each type of vessel, the number of vessels used does
not exceed the number of available vessels thanks to constraint (A.4).

3.8 Objective Function

Briefly, the objective function is the difference between the revenues generated by accepting
commodities into the network and the total costs of transporting them (fuel, vessel operations,
port calls, ...). The fuel cost depends on the fuel price and the fuel consumption of each trip
made. For this latter, we consider a variable cons<sub>p</sub><sup>r</sup> per port and rotation that specifies the
amount of fuel consumed per hour by the rotation r for the trip made between the port/canal
p and its successor. In the absence of successors, the variable cons<sub>p</sub><sup>r</sup> has, of course, the value
$$n_r = \left\lfloor \frac{T_r}{7 \times 24} \right\rfloor \qquad r \in \mathcal{R} \qquad (A.1)$$

$$v_r = 0 \iff T_r = 0 \qquad r \in \mathcal{R} \qquad (A.2)$$

$$\text{dep}_r = p' \wedge s_p^r = p' \Rightarrow T_r = \text{time}_{p,r}^{\text{out}} + \text{man}_{p,r}^{\text{out}} + \text{st}_p^r + \text{wt}_{p',r} + \text{man}_{p',r}^{\text{in}} \qquad r \in \mathcal{R}, p, p' \in \mathcal{P} \qquad (A.3)$$

$$\sum_{r \in \mathcal{R}} n_r (v_r = v) \le nb(v) \qquad v \in \mathcal{V} \qquad (A.4)$$

Figure 10 Constraints related to vessel availability.

$$\text{Elt}_m(\{\text{cons}(v, \nu) \mid v \in \{0\} \cup \mathcal{V}, \nu \in \{0\} \cup [\nu_{\min}(v), \nu_{\max}(v)]\}, v_r, v_p^r) = \text{cons}_p^r \qquad r \in \mathcal{R}, p \in \mathcal{P} \cup \mathcal{C} \qquad (O.1)$$

$$\text{teu}_p^{\text{ts}} = \sum_{k \in \mathcal{K} | p \neq \text{pod}(k)} \text{to}_{k,p}^r \cdot q(k) \qquad p \in \mathcal{P} \qquad (O.2)$$

Figure 11 Constraints related to the objective function.

0. The quantity consumed here depends only on the type of vessel used and the speed. In constraint (O.1) (see Figure 11), we assume that cons(v, ν) is 0 if v or ν is 0.

The costs associated with transshipment depend on the port and the quantity of commodities transshipped. So we need to represent the quantity of commodities transshipped at each port. To do this, we introduce a variable teu<sub>p</sub><sup>ts</sup> per port. The commodities k transshipped at port p are those that are unloaded at port p (i.e., those for which to<sub>k,p</sub><sup>r</sup> is 1) and for which port p is not their destination port (see constraint (O.2)). We can now express our objective function based on revenues (R), fuel costs (C), canal and port call costs (E), vessel operating costs (X), and transshipment costs (T):

$$\max \quad \sum_{k \in \mathcal{K}} rev(k) \cdot q(k) \cdot \alpha_k \qquad (R)$$

$$- \sum_{r \in \mathcal{R}} \sum_{p \in \mathcal{P} \cup \mathcal{C}} f p_r \cdot \text{cons}_p^r \cdot \text{st}_p^r \qquad (C)$$

$$- \sum_{r \in \mathcal{R}} \sum_{p \in \mathcal{P} \cup \mathcal{C}} p c_{pr} \cdot (s_p^r \neq p) \qquad (E)$$

$$- 7 \sum_{r \in \mathcal{R}} t c_r \cdot n_r \qquad (X)$$

$$- \sum_{p \in \mathcal{P}} ts(p) \cdot \text{teu}_p^{\text{ts}} \qquad (T)$$

## 3.9 Additional Constraints

Given the size of the search space, it may be desirable to avoid certain symmetries as much as possible. Starting each rotation at time 0 (see constraint (T.7)) allows for avoiding any translation on the time axis. However, other symmetries may exist. For example, the rotations are interchangeable. To avoid this, we can ensure that the first rotations are used in priority and these rotations are sorted in decreasing order of their duration thanks to constraint (S.1) (see Figure 12). One of the main practical difficulties of the LSNDP problem lies in the enumeration of the different possible circuits. For a circuit of length
16:14 A CP Approach for the Liner Shipping Network Design Problem

$$T_1 \ge T_2 \ge \dots \ge T_r \qquad (S.1)$$

$$\sum_{r \in \mathcal{R}} (s_p^r \neq p) \le \max(|K_p^{pol}|, |K_p^{pod}|) \qquad p \in \mathcal{P}_2 \qquad (S.2)$$

$$v_r = v \Rightarrow s_p^r = p \qquad r \in \mathcal{R}, p \in \mathcal{P}_3 \qquad (S.3)$$

$$teu_p^{ts} = 0 \qquad p \in \mathcal{P}_4 \qquad (S.4)$$

■ Figure 12 Possible additional constraints where $\mathcal{P}_i$ denotes the set of ports impacted by the constraint (S.i).

$\ell$, since it is possible to go from one port to any other, the solver may have to consider a non-negligible part of the $\ell$! possible permutations. Since, in addition, several rotations are usually considered simultaneously, this can quickly become very time-consuming. To reduce the number of rotations and thus of circuits to consider, we introduce constraint (S.2). Let $K_p^{pol}$ and $K_p^{pod}$ denote respectively the set of commodities for which $p$ is the origin port and one for which $p$ is the destination port. This constraint ensures that the number of rotations that uses a port $p$ does not exceed the maximum between $|K_p^{pol}|$ and $|K_p^{pod}|$. In a way, it eliminates some solutions in which the call in the port would only be used for transshipments. Generally, calling a port for only transshipments is not wished by shipping companies, except for some particular ports (e.g. hubs). Thanks to the flexibility of CP, we can add this constraint depending on the needs of the shipowner.

On the other hand, some ports cannot handle certain types of vessels. For example, the port of Dutch Harbor in Alaska is not deep enough. It can therefore only handle small container ships. Thus, if $v$ vessels cannot berth at port $p$, we impose constraint (S.3) for each rotation $r$. Similarly, some ports do not have enough space to store containers. It is therefore impossible to carry out transshipments there. For such ports, we can then exploit constraint (S.4) to prohibit any transshipment.

# 4 Experiments

## 4.1 Experimental Protocol and Implementation Details

The LINER-LIB benchmark<fnref n="1" /> [7] is the reference for experiments on the LSNDP problem. It consists of seven instances with 12 to 197 ports, thus allowing the evaluation of both complete and incomplete methods. In order to have instances of a reasonable and varied size, we have produced sub-instances from instances of the LINER-LIB benchmark. To do this, from an instance, we select $n$ ports in the following way. The first selected port is the one that handles the most commodities. The next $n-1$ ports are the ones that exchange the most commodities with the already selected ports. For our test set, we considered the smallest instance (Baltic) of the LINER-LIB set and 40 instances produced from the instances Baltic, EuropeAsia, Mediterranean and WAF. The number of ports varies from 3 to respectively 11, 10, 10 and 17. Moreover, as the LINER-LIB benchmark does not take into account productivity, waiting or maneuvering times, we generate randomly these values. Note that this partial random generation introduces no bias, since these values have a negligible impact on the solving efficiency. For the following experiments, the number of

<footnote n="1">¹ https://github.com/blof/LINERLIB/</footnote>
rotations $r_{max}$ is fixed at 4 and that of transhipments at 1 per commodity (a higher value not being desired by the experts) while the maximum duration $h_{max}$ is set to 12 weeks. The chosen values of $r_{max}$ and $h_{max}$ seem reasonable to us taking into account the commodities to be transported and the distances to be covered. In practice, the optimal solutions we found require less time and fewer rotations than our choice of parameters allows. For $\mu$, we take the value 0.54 given by the experts. The instances we consider are available at https://pageperso.lis-lab.fr/cyril.terrioux/LSNDP/instances.zip.

The presented model is implemented in the OR-Tools CP-SAT solver (version 9.6.2534 [26]) via its Python interface. This choice is first guided by the solver efficiency, since the OR-Tools CP-SAT solver won several gold medals during the past MiniZinc Challenges [25]. Moreover, lazy clause generation [20] provides good results for the LSFRP problem [14]. Finally, another advantage is the possibility of exploiting a certain form of parallelism. Hence, for the solving, we run from 1 to 16 threads. When a single thread is run, it corresponds to the CP-SAT solver. Except for the number of threads, all the parameters are the default ones. The experiments are being conducted on servers with Intel Xeon Gold 5218R processors running at 2.1 GHz and 192 GB of memory with a time limit of two hours. When exploiting several threads, each instance is solved 10 times and the reported runtime is the average time. The solving step involving $t$ threads is denoted $\times t$. We apply it to our model $M$, but also, to two derived versions denoted $M-1,2$ and $M-2$. Model $M-1,2$ does not consider constraints (S.1) and (S.2) while model $M-2$ uses constraints (S.1), but not constraints (S.2).

## 4.2 Results

First, we compare our model with its two derived versions from the efficiency viewpoint (see Table 1² and Figure 13). Clearly, model $M$ is the most efficient. Indeed, the addition of constraints (S.1) and (S.2) allows us to solve optimally more instances (24 instances for model $M$ against 14 for models $M-1,2$ and $M-2$ when using a single thread) while reducing significantly the runtime. As there often exists an arc between each pair of ports/canals, the circuit constraint admits a huge number of allowed tuples. In practice, the number of allowed tuples studied by the solver is mainly restricted by the load and transit time constraints or the objective function. In the latter case, as the objective function considers all the rotations, it may take some time for the solver to realize that a rotation is not suitable. So finding an optimal solution may require exploring a huge number of feasible solutions. Constraints (S.1) and (S.2) allow us to reduce this number significantly, as we can see in our results. Then, as shown in Table 1, exploiting several threads allows improving the efficiency, but mostly by reducing the runtime. For instance, using 16 threads instead of a single one leads to reducing the runtime by a factor of 5 on average and up to 20 at best.

Figure 13 indicates that the runtime increases exponentially with the number of ports. However, other parameters affect the runtime like the number of commodities to be processed or the type of instances. For example, our model performs well on Baltic and WAF instances, which corresponds to feedering instances (i.e. a collection of small services that ensures the transport of commodities between some main ports and satellite ones). In contrast, it turns out to be less efficient for instances like EuropeAsia ones that connect the more important ports of two commercial areas. Then, our model finds interesting solutions even if it does not accept all the commodities or visit all the ports (see Table 2). One explanation is related to

<footnote n="2">² See Appendix A for the other instances.</footnote>
16:16 A CP Approach for the Liner Shipping Network Design Problem

<table><caption>Table 1 Runtimes in seconds for some representative instances (a runtime different from 7,200 s corresponds to an instance solved optimally).</caption><thead><tr><th>Instance</th><th>M-1,2 x 1</th><th>M-2 x 1</th><th>M x 1</th><th>M x 2</th><th>M x 4</th><th>M x 8</th><th>M x 16</th></tr></thead><tbody><tr><td>Baltic</td><td>7,200</td><td>7,200</td><td>3,502</td><td>3,385</td><td>386</td><td>372</td><td>343</td></tr><tr><td>Baltic_sub7</td><td>1,483</td><td>301</td><td>6.7</td><td>4.8</td><td>3.0</td><td>3.3</td><td>3.2</td></tr><tr><td>EuropeAsia_sub7</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td></tr><tr><td>Mediterranean_sub7</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,001</td><td>6,234</td></tr><tr><td>WAF_sub7</td><td>7,200</td><td>7,200</td><td>27.4</td><td>19.0</td><td>6.6</td><td>5.6</td><td>5.7</td></tr><tr><td>WAF_sub17</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,021</td><td>5,117</td></tr></tbody></table>

<table><caption>Table 2 Information of some instances and solutions (the value in k$ of the best solution found, the number of visited ports, accepted commodities and used rotations).</caption><thead><tr><th colspan="3">Instance</th><th colspan="4">Solution</th></tr><tr><th>Name</th><th>|P|</th><th>|K|</th><th>Cost</th><th>#ports</th><th>#comm.</th><th>#rot</th></tr></thead><tbody><tr><td>Baltic</td><td>12</td><td>22</td><td>4,752</td><td>10</td><td>16</td><td>3</td></tr><tr><td>Baltic_sub7</td><td>7</td><td>12</td><td>2,508</td><td>6</td><td>9</td><td>2</td></tr><tr><td>EuropeAsia_sub7</td><td>7</td><td>42</td><td>4,228</td><td>6</td><td>23</td><td>3</td></tr><tr><td>Mediterranean_sub7</td><td>7</td><td>26</td><td>225</td><td>6</td><td>15</td><td>2</td></tr><tr><td>WAF_sub7</td><td>7</td><td>12</td><td>5,823</td><td>7</td><td>12</td><td>3</td></tr><tr><td>WAF_sub17</td><td>17</td><td>32</td><td>11,952</td><td>10</td><td>16</td><td>4</td></tr></tbody></table>

the way the LINER-LIB benchmark was built (namely by aggregating data from different shipowners without ensuring that the considered fleet can handle all the commodities).

Finally, our approach manages to optimally solve some instances of up to 17 ports. While the literature reports solving of up to a dozen ports, it is difficult to compare accurately with existing exact methods: implementations are generally not available and each work treats the problem with a different point of view, in particular regarding the working hypotheses or the cost function to optimize (see Section 5).

<figure><img src="image_146.png" /><figcaption>Figure 13 Runtime (in s) for Baltic (a) and WAF (b) instances w.r.t. the number of ports.</figcaption></figure>

5 Related Work

The Liner Shipping Fleet Repositioning Problem (LSFRP [28]) aims to adapt the network in order to take into account some evolution of the customer needs (e.g. seasonality, port congestion, increase or decrease of the demands, ...). Moving container ships from one service to another while considering commodity transport is a complex and expensive task
for shipping companies. In this context, the CP approach presented in [14] turns out to be more efficient than MIP ones. In the case of LSFRP, services are already defined while LSNDP aims to design them. Moreover, the number of vessels and commodities to handle may be reduced for LSFRP. Finally, an LSFRP instance corresponds to a one-shot task while an LSNDP one leads to a schedule for several weeks or months.

LSNDP is close to Vehicle Routing Problem (VRP [15]) and its variants like the pickup and delivery problem (PDP [11]). Indeed, in both cases, vehicles carry commodities from one location to another. The first difference is that generally, for VRP and PDP, a commodity is carried by a single vehicle whereas, for LSNDP, the transport can be achieved by several vessels operating different rotations thanks to transshipments. Taking into account transshipments makes the problem more difficult. This requires additional variables and constraints, while significantly increasing the number of possible routes for commodities (and therefore the number of feasible solutions to study). Moreover, in general, VRP and PDP aim to transport all the commodities while, for LSNDP, some commodities may be rejected. Finally, the objective function is often more complex for LSNDP than for VRP and PDP. For instance, the variety of considered costs (e.g. call cost, transshipment cost, fuel cost, charter rate, ...) is more important.

Regarding the exact solving of LSNDP, unlike our model, the models proposed in [18, 23, 22] use a constant speed for each leg and do not handle transit time constraints. Those of [18, 23] cannot reject a commodity. In contrast, [18] takes into account the empty container repositioning while [23, 22] consider the transshipment costs. The objective functions consist in minimizing costs [18, 23] or maximizing the profit [22]. Regarding the type of service, the three models consider a more general form than ours. For instance, they can exploit butterfly services (i.e. services that can call several times in the same port). However, such an extension could be taken into account in our model by duplicating ports as we do for canals and relaxing some constraints like constraints (F.9). Note that the experimentations achieved in [22] rely on the LINER-LIB benchmark, but the proposed approach does not succeed in solving optimally the Baltic instance.

# 6 Conclusions and Perspectives

In this paper, we have proposed a first CP model to solve the LSNDP problem. The first practical results are very encouraging with optimally solved instances with up to 17 ports and show the interest in a CP approach. In the future, this model will have to be extended to better handle some kinds of instances and take into account other forms of services (e.g. butterflies) or the constraints imposed by the International Maritime Organization (IMO) concerning the gas emissions of ships (e.g. related to carbon intensity indicator). Another extension would be to differentiate containers by type (full, empty, reefers, ... ). In particular, the repositioning of empty containers is an important issue, while the transport of reefers is highly profitable and raises specific questions. Moreover, to facilitate scaling up, the use of incomplete methods will need to be explored more deeply.

# References

1 R. Agarwal and Ö. Ergun. Ship Scheduling and Network Design for Cargo Routing in Liner Shipping. *Transportation Science*, 42(2):175–196, 2008. doi:10.1287/trsc.1070.0205.

2 J. F. Álvarez. Joint routing and deployment of a fleet of container vessels. *Maritime Economics & Logistics*, 11(2):186–208, 2009. doi:10.1057/mel.2009.5.
16:18 A CP Approach for the Liner Shipping Network Design Problem

3 N. Beldiceanu, M. Carlsson, and J.-X. Rampon. Global constraint catalog. Technical Report T2012:03, TASC-SICS-LINA, 2014.

4 N. Beldiceanu and E. Contejean. Introducing global constraints in CHIP. *Mathematical and Computer Modelling*, 20(12):97-123, 1994.

5 B. D. Brouer, G. Desaulniers, C. V. Karsten, and D Pisinger. A Matheuristic for the Liner Shipping Network Design Problem with Transit Time Restrictions. In *Computational Logistics - 6th International Conference, ICCL*, volume 9335 of *Lecture Notes in Computer Science*, pages 195-208. Springer, 2015. doi:10.1007/978-3-319-24264-4_14.

6 B. D. Brouer, C. V. Karsten, and D. Pisinger. Optimization in liner shipping. *Annals of Operations Research*, 271(1):205-236, 2018. doi:10.1007/s10479-018-3023-8.

7 B. D. Brouer, J. F. Álvarez, C. E. M. Plum, D. Pisinger, and M. M. Sigurd. A base integer programming model and benchmark suite for liner-shipping network design. *Transportation Science*, 48(2):281-312, 2014. doi:10.1287/trsc.2013.0471.

8 M. Carlsson, G. Ottosson, and B. Carlson. An open-ended finite domain constraint solver. In *Programming Languages: Implementations, Logics, and Programs (PLILP)*, volume 1292 of *Lecture Notes in Computer Science*, pages 191-206, 1997. doi:10.1007/BFb0033845.

9 M. Christiansen, E. O. Hellsten, D. Pisinger, D. Sacramento, and C. Vilhelmsen. Liner shipping network design. *European Journal of Operational Research*, 286, 2020. doi:10.1016/j.ejor.2019.09.057.

10 R. Cymer. Dulmage-Mendelsohn Canonical Decomposition as a generic pruning technique. *Constraints*, 17(3):234-272, 2012.

11 Y. Dumas, J. Desrosiers, and F. Soumis. The pickup and delivery problem with time windows. *European Journal of Operational Research*, 54(1):7-22, 1991. doi:10.1016/0377-2217(91)90319-Q.

12 C. V. Karsten, B. D. Brouer, G. Desaulniers, and D. Pisinger. Time constrained liner shipping network design. *Transportation Research. Part E: Logistics and Transportation Review*, 105:152-162, 2017. doi:10.1016/j.tre.2016.03.010.

13 C. V. Karsten, D. Pisinger, S. Røpke, and B. D. Brouer. The time constrained multi-commodity network flow problem and its application to liner shipping network design. *Transportation Research Part E: Logistics and Transportation Review*, 76:122-138, 2015. doi:10.1016/j.tre.2015.01.005.

14 E. Kelareva, K. Tierney, and P. Kilby. CP methods for scheduling and routing with time-dependent task costs. *EURO Journal on Computational Optimization*, 2(3):147-194, 2014. doi:10.1007/s13675-014-0022-7.

15 P. Kilby and P. Shaw. Vehicle routing. In F. Rossi, P. van Beek, and T. Walsh, editors, *Handbook of Constraint Programming*, volume 2 of *Foundations of Artificial Intelligence*, pages 801-836. Elsevier, 2006. doi:10.1016/S1574-6526(06)80027-1.

16 D. Kizilay, P. Van Hentenryck, and D. Türsel Eliiyi. Constraint programming models for integrated container terminal operations. *European Journal of Operational Research*, 286(3):945-962, 2020. doi:10.1016/j.ejor.2020.04.025.

17 A. Krogsgaard, D. Pisinger, and J. Thorsen. A flow-first route-next heuristic for liner shipping network design. *Networks*, 72(3):358-381, 2018. doi:10.1002/net.21819.

18 Q. Meng and S. Wang. Liner shipping service network design with empty container repositioning. *Transportation Research Part E: Logistics and Transportation Review*, 47(5):695-708, 2011. doi:10.1016/j.tre.2011.02.004.

19 Q. Meng and S. Wang. Optimal operating strategy for a long-haul liner service route. *European Journal of Operational Research*, 215(1):105-114, 2011. doi:10.1016/j.ejor.2011.05.057.

20 O. Ohrimenko, P. J. Stuckey, and M. Codish. Propagation via lazy clause generation. *Constraints An International Journal*, 14(3):357-391, 2009. doi:10.1007/s10601-008-9064-x.

21 C. E. M. Plum, D. Pisinger, J. J. Salazar-González, and M. M. Sigurd. Single liner shipping service design. *Computers & Operations Research*, 45:1-6, 2014. doi:10.1016/j.cor.2013.11.018.
22 C. E. M. Plum, D. Pisinger, and M. M. Sigurd. A service flow model for the liner shipping network design problem. *European Journal of Operational Research*, 235(2):378-386, 2014. doi:10.1016/j.ejor.2013.10.057.

23 L. B. Reinhardt and D. Pisinger. A branch and cut algorithm for the container shipping network design problem. *Flexible Services and Manufacturing Journal*, 24(3):349-374, 2012. doi:10.1007/s10696-011-9105-4.

24 D. Sacramento, C. Solnon, and D. Pisinger. Constraint Programming and Local Search Heuristic: a Matheuristic Approach for Routing and Scheduling Feeder Vessels in Multi-terminal Ports. *Annals of Operations Research*, 1(4), 2020. doi:10.1007/s43069-020-00036-x.

25 P. J. Stuckey, T. Feydy, A. Schutt, G. Tack, and J. Fischer. The MiniZinc Challenge 2008-2013. *AI Magazine*, 35(2):55-60, 2014.

26 OR-Tools team. OR-Tools CP-SAT solver, 2023. URL: https://developers.google.com/optimization/cp/cp_solver.

27 K. Thun, H. Andersson, and M. Christiansen. Analyzing complex service structures in liner shipping network design. *Flexible Services and Manufacturing Journal*, 29(3-4):535-552, 2017. doi:10.1007/s10696-016-9262-6.

28 K. Tierney. *Optimizing Liner Shipping Fleet Repositioning Plans*, volume 57 of *Operations research / computer science interfaces series*. Springer, 2015. doi:10.1007/978-3-319-17665-9.

29 P. Van Hentenryck and J.-P. Carillon. Generality versus Specificity: An Experience with AI and OR Techniques. In *Proceedings of the National Conference on Artificial Intelligence*, pages 660–664, 1988.

30 S. Wang and Q. Meng. Liner shipping network design with deadlines. *Computers & Operations Research*, 41:140-149, 2014. doi:10.1016/j.cor.2013.08.014.

## A Detailed Results

Tables 3 and 4 provide the same information as Tables 1 and 2 but for all the instances we consider.
<table><caption>Table 3 Runtimes in seconds for all the considered instances (a runtime different from 7,200 s corresponds to an instance solved optimally).</caption><thead><tr><th>Instance</th><th>M-1,2 x 1</th><th>M-2 x 1</th><th>M x 1</th><th>M x 2</th><th>M x 4</th><th>M x 8</th><th>M x 16</th></tr></thead><tbody><tr><td>Baltic_sub3</td><td>0.4</td><td>0.2</td><td>0.2</td><td>0.2</td><td>0.2</td><td>0.3</td><td>0.3</td></tr><tr><td>Baltic_sub4</td><td>15.8</td><td>4.5</td><td>0.6</td><td>0.4</td><td>0.4</td><td>0.4</td><td>0.5</td></tr><tr><td>Baltic_sub5</td><td>134</td><td>16.9</td><td>1.6</td><td>2.8</td><td>2.3</td><td>1.6</td><td>1.2</td></tr><tr><td>Baltic_sub6</td><td>256</td><td>120</td><td>2.9</td><td>1.9</td><td>1.3</td><td>1.4</td><td>1.4</td></tr><tr><td>Baltic_sub7</td><td>1,483</td><td>301</td><td>6.7</td><td>4.8</td><td>3.0</td><td>3.3</td><td>3.2</td></tr><tr><td>Baltic_sub8</td><td>7,200</td><td>7,200</td><td>29.2</td><td>14.2</td><td>5.8</td><td>5.2</td><td>5.1</td></tr><tr><td>Baltic_sub9</td><td>7,200</td><td>7,200</td><td>87.7</td><td>55.8</td><td>12.9</td><td>9.8</td><td>10.6</td></tr><tr><td>Baltic_sub10</td><td>7,200</td><td>7,200</td><td>158</td><td>126</td><td>22.9</td><td>16.9</td><td>15.2</td></tr><tr><td>Baltic_sub11</td><td>7,200</td><td>7,200</td><td>637</td><td>504</td><td>49.6</td><td>42.6</td><td>38.7</td></tr><tr><td>Baltic</td><td>7,200</td><td>7,200</td><td>3,502</td><td>3,385</td><td>386</td><td>372</td><td>343</td></tr><tr><td>EuropeAsia_sub3</td><td>2.5</td><td>1.0</td><td>0.9</td><td>1.0</td><td>1.2</td><td>1.3</td><td>1.2</td></tr><tr><td>EuropeAsia_sub4</td><td>157</td><td>244</td><td>134</td><td>100</td><td>37.3</td><td>18.1</td><td>10.5</td></tr><tr><td>EuropeAsia_sub5</td><td>7,200</td><td>7,200</td><td>7,200</td><td>3,736</td><td>4,101</td><td>4,459</td><td>2,356</td></tr><tr><td>EuropeAsia_sub6</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td></tr><tr><td>EuropeAsia_sub7</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td></tr><tr><td>EuropeAsia_sub8</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td></tr><tr><td>EuropeAsia_sub9</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td></tr><tr><td>EuropeAsia_sub10</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td></tr><tr><td>Mediterranean_sub3</td><td>0.4</td><td>0.2</td><td>0.3</td><td>0.3</td><td>0.3</td><td>0.3</td><td>0.3</td></tr><tr><td>Mediterranean_sub4</td><td>2.8</td><td>0.8</td><td>0.8</td><td>0.8</td><td>0.8</td><td>0.9</td><td>1.0</td></tr><tr><td>Mediterranean_sub5</td><td>85.0</td><td>28.9</td><td>26.5</td><td>22.4</td><td>5.2</td><td>3.8</td><td>3.8</td></tr><tr><td>Mediterranean_sub6</td><td>621</td><td>229</td><td>281</td><td>307</td><td>30.6</td><td>25.2</td><td>21.5</td></tr><tr><td>Mediterranean_sub7</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,001</td><td>6,234</td></tr><tr><td>Mediterranean_sub8</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td></tr><tr><td>Mediterranean_sub9</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td></tr><tr><td>Mediterranean_sub10</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td></tr><tr><td>WAF_sub3</td><td>0.9</td><td>0.7</td><td>0.6</td><td>0.3</td><td>0.3</td><td>0.4</td><td>0.3</td></tr><tr><td>WAF_sub4</td><td>401</td><td>42.8</td><td>0.8</td><td>0.6</td><td>0.7</td><td>0.6</td><td>1.1</td></tr><tr><td>WAF_sub5</td><td>896</td><td>248</td><td>1.8</td><td>1.8</td><td>1.6</td><td>1.6</td><td>1.4</td></tr><tr><td>WAF_sub6</td><td>7,200</td><td>7,200</td><td>7.5</td><td>4.4</td><td>3.9</td><td>3.4</td><td>3.5</td></tr><tr><td>WAF_sub7</td><td>7,200</td><td>7,200</td><td>27.4</td><td>19.0</td><td>6.6</td><td>5.6</td><td>5.7</td></tr><tr><td>WAF_sub8</td><td>7,200</td><td>7,200</td><td>31.3</td><td>19.9</td><td>6.7</td><td>6.1</td><td>6.1</td></tr><tr><td>WAF_sub9</td><td>7,200</td><td>7,200</td><td>17.0</td><td>18.0</td><td>7.5</td><td>7.0</td><td>7.0</td></tr><tr><td>WAF_sub10</td><td>7,200</td><td>7,200</td><td>39.4</td><td>35.7</td><td>16.0</td><td>12.3</td><td>12.7</td></tr><tr><td>WAF_sub11</td><td>7,200</td><td>7,200</td><td>249</td><td>187</td><td>32.2</td><td>32.6</td><td>29.3</td></tr><tr><td>WAF_sub12</td><td>7,200</td><td>7,200</td><td>1,093</td><td>600</td><td>73.4</td><td>66.7</td><td>64.9</td></tr><tr><td>WAF_sub13</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>815</td><td>792</td><td>499</td></tr><tr><td>WAF_sub14</td><td>7,200</td><td>7,200</td><td>4,662</td><td>3,779</td><td>323</td><td>262</td><td>263</td></tr><tr><td>WAF_sub15</td><td>7,200</td><td>7,200</td><td>7,034</td><td>5,365</td><td>542</td><td>494</td><td>349</td></tr><tr><td>WAF_sub16</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>4,182</td><td>3,106</td><td>2,154</td></tr><tr><td>WAF_sub17</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,200</td><td>7,021</td><td>5,117</td></tr></tbody></table>
<table><caption>Table 4 Some information about instances and solutions. For solutions, we provide the value (in k$) of the best solution found, the number of visited ports, one of accepted commodities and one of used rotations.</caption><thead><tr><th colspan="3">Instance</th><th rowspan="2">Cost</th><th colspan="3">Solution</th></tr><tr><th>Name</th><th>|P|</th><th>|K|</th><th>#ports</th><th>#comm.</th><th>#rot</th></tr></thead><tbody><tr><td>Baltic_sub3</td><td>3</td><td>4</td><td>1,876</td><td>3</td><td>4</td><td>1</td></tr><tr><td>Baltic_sub4</td><td>4</td><td>6</td><td>1,895</td><td>3</td><td>4</td><td>1</td></tr><tr><td>Baltic_sub5</td><td>5</td><td>8</td><td>2,074</td><td>5</td><td>8</td><td>2</td></tr><tr><td>Baltic_sub6</td><td>6</td><td>10</td><td>2,074</td><td>5</td><td>8</td><td>2</td></tr><tr><td>Baltic_sub7</td><td>7</td><td>12</td><td>2,508</td><td>6</td><td>9</td><td>2</td></tr><tr><td>Baltic_sub8</td><td>8</td><td>14</td><td>3,322</td><td>7</td><td>10</td><td>3</td></tr><tr><td>Baltic_sub9</td><td>9</td><td>16</td><td>3,733</td><td>8</td><td>11</td><td>3</td></tr><tr><td>Baltic_sub10</td><td>10</td><td>18</td><td>4,187</td><td>9</td><td>13</td><td>3</td></tr><tr><td>Baltic_sub11</td><td>11</td><td>20</td><td>4,345</td><td>10</td><td>15</td><td>3</td></tr><tr><td>Baltic_sub12</td><td>12</td><td>22</td><td>4,752</td><td>10</td><td>16</td><td>3</td></tr><tr><td>EuropeAsia_sub3</td><td>3</td><td>6</td><td>616</td><td>3</td><td>4</td><td>1</td></tr><tr><td>EuropeAsia_sub4</td><td>4</td><td>12</td><td>616</td><td>3</td><td>4</td><td>1</td></tr><tr><td>EuropeAsia_sub5</td><td>5</td><td>20</td><td>1,463</td><td>4</td><td>8</td><td>1</td></tr><tr><td>EuropeAsia_sub6</td><td>6</td><td>30</td><td>1,463</td><td>4</td><td>8</td><td>1</td></tr><tr><td>EuropeAsia_sub7</td><td>7</td><td>42</td><td>4,228</td><td>6</td><td>23</td><td>3</td></tr><tr><td>EuropeAsia_sub8</td><td>8</td><td>56</td><td>4,425</td><td>6</td><td>19</td><td>2</td></tr><tr><td>EuropeAsia_sub9</td><td>9</td><td>72</td><td>4,425</td><td>6</td><td>19</td><td>2</td></tr><tr><td>EuropeAsia_sub10</td><td>10</td><td>89</td><td>2,935</td><td>7</td><td>24</td><td>2</td></tr><tr><td>Mediterranean_sub3</td><td>3</td><td>5</td><td>177</td><td>3</td><td>5</td><td>1</td></tr><tr><td>Mediterranean_sub4</td><td>4</td><td>9</td><td>177</td><td>3</td><td>5</td><td>1</td></tr><tr><td>Mediterranean_sub5</td><td>5</td><td>14</td><td>196</td><td>4</td><td>8</td><td>1</td></tr><tr><td>Mediterranean_sub6</td><td>6</td><td>20</td><td>196</td><td>4</td><td>8</td><td>1</td></tr><tr><td>Mediterranean_sub7</td><td>7</td><td>26</td><td>225</td><td>6</td><td>15</td><td>2</td></tr><tr><td>Mediterranean_sub8</td><td>8</td><td>34</td><td>302</td><td>8</td><td>30</td><td>2</td></tr><tr><td>Mediterranean_sub9</td><td>9</td><td>43</td><td>509</td><td>8</td><td>30</td><td>2</td></tr><tr><td>Mediterranean_sub10</td><td>10</td><td>53</td><td>605</td><td>9</td><td>38</td><td>2</td></tr><tr><td>WAF_sub3</td><td>3</td><td>4</td><td>1,293</td><td>3</td><td>4</td><td>1</td></tr><tr><td>WAF_sub4</td><td>4</td><td>6</td><td>1,308</td><td>4</td><td>6</td><td>3</td></tr><tr><td>WAF_sub5</td><td>5</td><td>8</td><td>2,329</td><td>5</td><td>8</td><td>3</td></tr><tr><td>WAF_sub6</td><td>6</td><td>10</td><td>2,911</td><td>6</td><td>10</td><td>3</td></tr><tr><td>WAF_sub7</td><td>7</td><td>12</td><td>5,823</td><td>7</td><td>12</td><td>3</td></tr><tr><td>WAF_sub8</td><td>8</td><td>14</td><td>5,823</td><td>7</td><td>12</td><td>3</td></tr><tr><td>WAF_sub9</td><td>9</td><td>16</td><td>5,823</td><td>7</td><td>12</td><td>3</td></tr><tr><td>WAF_sub10</td><td>10</td><td>18</td><td>7,543</td><td>8</td><td>14</td><td>3</td></tr><tr><td>WAF_sub11</td><td>11</td><td>20</td><td>8,831</td><td>9</td><td>15</td><td>4</td></tr><tr><td>WAF_sub12</td><td>12</td><td>22</td><td>9,764</td><td>10</td><td>16</td><td>4</td></tr><tr><td>WAF_sub13</td><td>13</td><td>24</td><td>11,282</td><td>11</td><td>18</td><td>4</td></tr><tr><td>WAF_sub14</td><td>14</td><td>26</td><td>11,952</td><td>10</td><td>16</td><td>4</td></tr><tr><td>WAF_sub15</td><td>15</td><td>28</td><td>11,952</td><td>10</td><td>16</td><td>4</td></tr><tr><td>WAF_sub16</td><td>16</td><td>30</td><td>11,952</td><td>10</td><td>16</td><td>4</td></tr><tr><td>WAF_sub17</td><td>17</td><td>32</td><td>11,952</td><td>10</td><td>16</td><td>4</td></tr></tbody></table>
<figure><img src="image_147.png" /></figure>
Optimization Models for Pickup-And-Delivery
Problems with Reconfigurable Capacities

Arnoosh Golestanian ⬒ ID

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

Giovanni Lo Bianco ⬒ ID

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

Chengyu Tao ⬒ ID

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

J. Christopher Beck ⬒ ID

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

Abstract

When a transportation service accommodates both people and goods, operators sometimes opt for vehicles that can be dynamically reconfigured for different demands. Motivated by air service in remote communities in Canada's north, we define a pickup-and-delivery problem in which aircraft can add or remove seats during a multi-stop trip to accommodate varying demands. Given the demand for people and cargo as well as a seat inventory at each location, the problem consists in finding a tour that picks up and delivers all demand while potentially reconfiguring the vehicle capacity at each location by adding or removing seats. We develop a total of six models using three different approaches: constraint programming, mixed integer programming, and domain-independent dynamic programming. Our numerical experiments indicate that domain-independent dynamic programming is able to substantially outperform the other technologies on both solution quality and run-time on a set of randomly generated instances spanning the size of real problems in northern Canada.

2012 ACM Subject Classification Computing methodologies → Modeling methodologies

Keywords and phrases Pickup and delivery, Dial-a-ride problem, Optimization

Digital Object Identifier 10.4230/LIPIcs.CP.2023.17

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

# 1 Introduction

Pickup-and-delivery problems involve using vehicles to transport goods and/or passengers from a set of origins to a set of destinations on a given transportation network [1]. A typical pickup-and-delivery problem such as the Pickup and Delivery Traveling Salesperson Problem (PD-TSP) includes a one or more vehicles, requests with different pickup and delivery locations, and an objective to find a minimum-cost tour (or set of routes) that visit(s) each pickup location before its corresponding delivery location [4]. There has been substantial research literature on pickup and delivery problems over the past several years (e.g., [19, 21]) motivated, in part, by global efforts to reduce transportation-related carbon emissions [16]. Many variations of such problems have been proposed and studied in the operations research literature. For example, some problems include handling costs when an item is loaded or unloaded depending on the position of the item in the vehicle [24] and some include subsets of requests that cannot be in a vehicle at the same time [5].

In this paper, we propose and study a novel variation of PD-TSP: requests can include both goods (cargo) and passengers and the vehicle has a capacity that can be adjusted en-route depending on the request and equipment stored at locations in the network. The

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

© Arnoosh Golestanian, Giovanni Lo Bianco, Chengyu Tao, and J. Christopher Beck;
licensed under Creative Commons License CC-BY 4.0

29th International Conference on Principles and Practice of Constraint Programming (CP 2023).

Editor: Roland H. C. Yap; Article No. 17; pp. 17:1–17:17

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

LIPICS

Leibniz International Proceedings in Informatics

Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl Publishing, Germany
17:2 Optimization Models for PDPs with Reconfigurable Capacities

problem is motivated by a real transportation problem faced by air services in northern
Canada. Since many communities in this region are reachable only by air during some parts
of the year, their access to basic human needs such as fresh food and healthcare services is
limited. The need for air transportation combined with the relatively small populations and
lack of resources led northern air services to adopt the practice of transporting both cargo
and passengers on the same flights. The vehicles are aircraft with removable seats, allowing
staff to either remove passenger seats and store them at airports to transport more cargo or
add additional seats, previously stored at airports, to carry more passengers. The problem,
which we call the Pickup-and-Delivery with Seat Replacement Problem (PD-SRP), therefore
requires finding the shortest tour delivering all goods and passengers from their origins to
destinations without exceeding aircraft capacity but allowing seats to be removed from or
added to aircraft at each location, subject to seat availability and total aircraft capacity.

To solve the PD-SRP, we developed three types of optimization models: one Constraint Programming (CP) model, three Mixed Integer Programming (MIP) models, and two Domain-Independent Dynamic Programming (DIDP) [14] models. We compare their performance on randomly generated instances based on the size of the problem in Canada's north, demonstrating that both of the DIDP models outperform the CP and MIP models in terms of the number of instances solved and proved optimal, solution quality, and solve time.

2 Related Works

Reconfigurable capacity is a general term in the transportation literature, typically indicating that vehicle capacity can be changed at some cost and/or limited by some constraints [22, 23]. Other terms such as multi-compartment vehicle or multi-purpose vehicle are used to convey a similar meaning [20, 8]. We review the vehicle routing and dial-a-ride problems literature for studies that considered adjustable vehicles.

*Vehicle Routing Problems (VRP)*: The Vehicle Routing Problem and its many variations have been studied extensively over the past 50 years [18]. The idea of adjusting the vehicle to handle different types of demand has been studied in multi-compartment vehicle routing problems [20]. For example, Henke et al. [9] studied how to split the capacity of a truck into different compartments to maintain the separation of different colors of recycled glass. Similarly, for grocery distribution, different temperature-sensitive products can be transported on the same truck with multiple compartments [11]. In both of these problems, a vehicle's capacity configuration is fixed for its entire route and cannot be modified during the trip.

*Dial-a-Ride Problems (DARP)*: In the Dial-a-Ride Problem a transportation request takes the form of pickup and delivery location pair and the service provider utilizes its fleet of vehicles to fulfill the requests while minimizing a cost function that typically includes some travel distance component [10]. Some variants include a reconfigurable vehicle capacity to serve the needs of different users: those who use seats or those who use wheelchairs [23]. Some of the vehicle seats can be folded and stored inside the vehicle to make room for passengers in wheelchairs. Unlike this problem, the seats of the vehicle in the PD-SRP cannot be stored on the aircraft without occupying cargo capacity and are instead detached and stored at the airports.

Hatzenbühler et al. [8] studied a multi-purpose pickup and delivery problem that can deliver passengers or cargo by exchanging the module of the vehicle at a depot or special service site. Each vehicle includes a removable module and a fixed platform such that changing modules modifies the ability of the vehicle from only carrying cargo to only carrying passengers and vice versa. Compared to problems with conventional solo-purpose vehicles,
requests can be served with a fewer vehicles but at the expense of adding extra service sites and visits. We can view the core multi-purpose pickup and delivery problem as a special case of PD-SRP where the seat exchange decisions must be all-or-none: either all seats are removed to maximize cargo space or all seats are installed to maximize passenger capacity.

### 3 Problem Definition

In PD-SRP, we are given $n$ requests, each potentially requiring the transportation of cargo and passenger demands. Let $V = P \cup D$ where $P = \{v_1, \dots, v_n\}$ is the set of pickup locations and $D = \{v_{n+1}, \dots, v_{2n}\}$ is the set of delivery locations. We assumed that cargo is shipped in unit-sized boxes, each having the same weight and volume. Although, in reality cargo is shipped in various shapes and weights, incorporating four-dimensional packing (i.e., combining volume and weight) would substantially complicate the problem. Therefore, similar to approximations done in practice by airlines (e.g., standard weight per passenger), we opted for this simplification.

Each request $i$ includes picking up $\hat{q}_i$ boxes of cargo and $\hat{\pi}_i$ passengers from location $v_i$ and delivering them to location $v_{n+i}$. Thus, the demand of the corresponding delivery location has an equal magnitude negative value (i.e., $-\hat{q}_i = \hat{q}_{i+n}, -\hat{\pi}_i = \hat{\pi}_{i+n}, \forall i \in P$). Note that this representation can model more complex patterns (e.g., requests that share pickup or delivery locations but not both) by copying locations for each unique pickup-delivery pair.

When an aircraft is at its maximum seat capacity, it has $\hat{S}$ seats and can carry $\hat{C}$ boxes of cargo. By removing a seat, $L$ boxes of cargo capacity are added to the aircraft. Therefore, the maximum cargo capacity when removing all the seats is $K = \hat{S}L + \hat{C}$. Each location $i$ starts with $S_i^0$ stored seats and therefore the aircraft can add at most $\min(\hat{S}, S_i^0)$ seats or remove at most $\hat{S}$ seats when visiting location $i$. There is no maximum number of seats that can be stored at a given location.

In order to represent the problem as a path, two nodes are assigned to the depot: $v_0$ is the start node and $v_{2n+1}$ is the end node. For modeling purposes we define sets $V_{N+1} = V \cup \{v_{2n+1}\}, V_0 = V \cup \{v_0\}$ and $V_{0,N+1} = V \cup \{v_0, v_{2n+1}\}$. Therefore, the problem is defined on graph $G = (V_{0,N}, A)$ where $A = \{(i,j) | i, j \in V_{0,N+1}, i \neq j\}$ with each arc having an associated distance, $d_{ij}$. The vehicle is initially at the depot $v_0$ with a cargo and passenger capacity of $\tilde{C}_0$ and $\tilde{S}_0$ where $\tilde{C}_0 = K - L\tilde{S}_0$ and $\tilde{S}_0 \le \hat{S}$, respectively, and must finish the trip at depot $v_{2n+1}$. We assume that the start and end nodes are not the pickup or delivery location of any requests. Again, this assumption is not limiting as such requests can be represented by adding extra nodes at the same location as the start and end nodes.

In PD-SRP we aim to minimize the travel distance while deciding how many seats to add or remove at each location to fulfill all the requests while respecting capacities. The PD-SRP is NP-hard because if we fix the seat decisions and set all the demands to zero, the problem can be reduced to TSP which is known to be NP-hard [13].

An instance of this problem can be seen in Figure 1. The optimal tour is shown in pink, and the seat icon near each vertex represents the number of seats stored at the corresponding base. The optimal tour for this instance is $(v_0, v_2, v_4, v_1, v_3, v_5)$ with two seats left at $v_1$.

### 4 Methods

We develop six models for the PD-SRP using constraint programming (CP), mixed integer programming (MIP), and domain independent dynamic programming (DIDP). One of the MIP models solves a restricted version of the PD-SRP and is used to warm-start the CP model and the two other MIP models. In this section, we describe each of the models in detail.
<figure><img src="image_150.png" /><figcaption>Figure 1 Example of an PD-SRP instance with 2 requests. The optimal tour is shown by dotted pink edges. In the aircraft configuration, white seats and boxes show the current passenger and cargo capacity, respectively. The colored seats and boxes show the corresponding cargo and passenger requests that are picked up.</figcaption></figure>

## 4.1 A Constraint Programming Model (CP)

Our CP model equates distance and time and, thus, uses a one-machine scheduling approach where jobs correspond to the visits and the setup times between two consecutive jobs correspond to the distance between two locations. The model uses $|V_{0,2n+1}|$ interval variables $x_i$ that represent visits to each location, and a sequence variable, $\pi$, that constrains interval variables to form a sequence with an extra end node representing the return to the depot. The size of the interval variable is 0 because there is no service time associated with the visits. For every location $i \in \{0, ..., 2n\}$, variable $s_i$ is introduced to represent the number of seats that are added or removed. The formulation of the CP model is presented in Figure 2. Note that CP model is written in CP Optimizer language.

The objective function is the minimization of the total distance traveled by the aircraft. EndOf($x_{2n+1}$) corresponds to the end-point of the last interval variable in the sequence variable $\pi$: the time (i.e., total distance travelled) when the aircraft returns to the depot. Constraint (1a) ensures that each pair of consecutive interval variables is scheduled with a transition time equal to at least the required travel distance between the two corresponding locations. Constraint (1b) enforces that the pickup location of each request is visited before the delivery location. Constraint (1c) specifies that the aircraft begins and ends at the start and end depot locations.

We used three cumulative functions to represent the following values that are potentially changed by each interval variable (aircraft visits): available cargo space, number of empty seats, and the total number of seats. In particular, cumulative functions (1d) and (1f) are used to represent the available passenger and cargo space as the trip proceeds. $H$ represents the number of empty seats in the aircraft (i.e., the available passenger space) and $C$ represents
$$ \min \text{ EndOf}(x_{2n+1}) \qquad (\mathbb{CP}) $$

$$ \text{s.t.} \ NoOverlap(\pi, \{d_{i,j} : (i, j) \in A\}) \qquad (1a) $$

$$ \text{EndBeforeStart}(x_i, x_{n+i}) \qquad \forall i \in \{1, \dots, n\} \qquad (1b) $$

$$ \text{First}(\pi, x_0), \text{Last}(\pi, x_{2n+1}) \qquad (1c) $$

$$ C = \text{StepAt}(x_0, \tilde{C}_0) + \sum_{i=0}^{2n} \text{StepAtStart}(x_i, -\hat{q}_i - L \cdot s_i) \qquad (1d) $$

$$ C \geq 0 \qquad (1e) $$

$$ H = \text{StepAt}(x_0, \tilde{S}_0) + \sum_{i=0}^{2n} \text{StepAtStart}(x_i, -\hat{\pi}_i + s_i) \qquad (1f) $$

$$ H \geq 0 \qquad (1g) $$

$$ S = \text{StepAt}(x_0, \tilde{S}_0) + \sum_{i=0}^{2n} \text{StepAtStart}(x_i, s_i) \qquad (1h) $$

$$ S \leq \hat{S} \qquad (1i) $$

$$ x_i : \text{intervalVar}(0) \qquad \forall i \in V_{0,2n+1} \qquad (1j) $$

$$ s_i : \text{integerVar}(-\hat{S}, \min(\hat{S}, S_i^0)) \qquad \forall i \in V_0 \qquad (1k) $$

$$ \pi : \text{sequenceVar}(x_0, \dots, x_{2n+1}) \qquad (1l) $$

<caption>Figure 2 The CP Model for the PD-SRP.</caption>

the available cargo space. Before the trip starts, $K = H + C$ and, if there are $\hat{S}_0$ seats in the aircraft at the start, $H = L \cdot \hat{S}_0$. The expression $\text{StepAtStart}(var, \text{impact})$ specifies the change (increment or decrement) to the cumulative function at the start of an interval variable. The available cargo space $C$ will decrease as cargo and seats are picked up, therefore we use $\text{StepAtStart}(x_i, -\hat{q}_i - L \cdot s_i)$ to represent the changes to available cargo space at each location $i \in \{1, \dots, 2n\}$. The available passenger space will decrease when cargo is picked up, while increasing when adding seats as represented by $\text{StepAtStart}(x_i, -\hat{\pi}_i + s_i)$ at each location $i \in \{0, \dots, 2n\}$. The cumulative function $S$ is introduced in constraint (1h) to describe the change of the total number of seats in the aircraft. $S$ will change with the number of seats being added or removed as represented by $s_i$. Constraint (1i) restricts the total number of seats by the maximum seat capacity $\hat{S}$. In constraint (1k), the domain of $s_i$ is $[-\hat{S}, S_i^0]$ reflecting the range of the number of seats that the aircraft can remove or add at location $i$.

It should be noted that every interval variable contributes to the cumulative constraint, which means that these bounds are maintained throughout the sequence. Therefore, we do not need to have a separate cumulative function for every location of the tour.

## 4.2 Mixed Integer Programming Models

In this section, we describe three MIP models motivated by existing models for pickup and delivery problems. The first two models exactly represent the PD-SRP and therefore admit optimal solutions. The final model is a restriction of the PD-SRP problem that can be used to quickly find a feasible solution and, therefore, an upper bound for the PD-SRP. In our experiments, we investigate the use of this restricted model to warm-start the CP model and two other MIP models.
17:6 Optimization Models for PDPs with Reconfigurable Capacities

$$
\min \sum_{i \in V_0} \sum_{j \in V_{N+1}, i \neq j} d_{ij} x_{ij} \qquad (\mathit{MIP}_{loc})
$$

$$
\sum_{j \in V_{N+1}} x_{ij} = 1 \qquad i \in V_0 \qquad (2a)
$$

$$
\sum_{j \in V_0, j \neq i} x_{ji} - \sum_{j \in V_{N+1}, i \neq j} x_{ij} = 0 \qquad i \in V \qquad (2b)
$$

$$
\tau_i + x_{ij} - |V|(1 - x_{ij}) \le \tau_j \qquad i \in V_0, j \in V_{N+1}, i \neq j \qquad (2c)
$$

$$
1 \le \tau_i \le |V| \qquad i \in V \qquad (2d)
$$

$$
\tau_i + 1 \le \tau_{n+i} \qquad i \in P \qquad (2e)
$$

$$
y_i + \pi_i \le \hat{S} \qquad i \in V_0 \qquad (2f)
$$

$$
y_i + \pi_i + s_i \le \hat{S} \qquad i \in V_0 \qquad (2g)
$$

$$
\pi_i + \hat{\pi}_i \le \hat{S} \qquad i \in V_0 \qquad (2h)
$$

$$
u_i + q_i + Ly_i + L\pi_i \le K \qquad i \in V_0 \qquad (2i)
$$

$$
u_j \le u_i - Ls_i - \hat{q}_i x_{ij} + (2K)(1 - x_{ij}) \quad i \in V, j \in V_{N+1}, i \neq j \quad (2j)
$$

$$
u_j \geq u_i - Ls_i - \hat{q}_i x_{ij} - (2K)(1 - x_{ij}) \qquad i \in V, j \in V_{N+1}, i \neq j \qquad (2k)
$$

$$
y_j \le y_i + s_i - \hat{\pi}_i x_{ij} + 2\hat{S}(1 - x_{ij}) \qquad i \in V, j \in V_{N+1}, i \neq j \qquad (2l)
$$

$$
y_j \geq y_i + s_i - \hat{\pi}_i x_{ij} - 2\hat{S}(1 - x_{ij}) \qquad i \in V, j \in V_{N+1}, i \neq j \quad (2m)
$$

$$
\pi_j \ge \pi_i + \hat{\pi}_i x_{ij} - \hat{S}(1 - x_{ij}) \qquad i \in V, j \in V_{N+1}, i \neq j \quad (2n)
$$

$$
q_j \geq q_i + \hat{q}_i x_{ij} - K(1 - x_{ij}) \qquad i \in V, j \in V_{N+1}, i \neq j \qquad (2o)
$$

$$
-y_i \le s_i \le \min(\hat{S}, S_i^0) \qquad i \in V \qquad (2p)
$$

$$
\tau_0 = \pi_0 = q_0 = 0, u_0 = \tilde{C}_0, y_0 = \tilde{S}_0 \qquad (2q)
$$

$$
x_{ij} \in \{0, 1\} \qquad i \in V_0, j \in V_{N+1}, i \neq j \qquad (2r)
$$

$$
u_i, y_i, \pi_i, q_i, \tau_i \in \mathbb{R}^{0+}, s_i \in \mathbb{R} \qquad i \in V_{0,N+1} \qquad (2s)
$$

<caption>Figure 3 The MIP<sub>loc</sub> Model for the PD-SRP.</caption>

4.2.1 Two-indexed Location-Based MIP (MIP<sub>loc</sub>)

We propose a two-indexed location-based MIP model for PD-SRP (MIP<sub>loc</sub>) based on a model
for an existing pickup and delivery variant [7]. In MIP<sub>loc</sub>, x<sub>ij</sub> is a binary variable that is 1 if
arc (i, j) ∈ A is traveled and is 0 otherwise. Non-negative continuous variables τ<sub>i</sub>, u<sub>i</sub>, and
y<sub>i</sub> represent the distance, available cargo space, and empty seats, respectively, on arrival at
vertex i ∈ V<sub>0,N+1</sub>. As above, let variable s<sub>i</sub> be the number of seats that are added (s<sub>i</sub> > 0) or
removed (s<sub>i</sub> < 0) at location i. Finally, let π<sub>i</sub> and q<sub>i</sub> be the number of passengers and boxes
of cargo on the aircraft on arrival at vertex i ∈ V<sub>0,N+1</sub>, respectively. The MIP<sub>loc</sub> model is
shown in Figure 3.

The objective function minimizes the total distance traveled. Constraint (2a) ensures that
each customer is visited exactly once while constraint (2b) forces an arrival and departure
at each non-depot vertex. Constraints (2c) and (2d) prevent the formation of the subtours,
using Miller-Tucker-Zemlin (MTZ) constraints [17]. Constraint (2e) forces the aircraft to
visit the pickup location of each commodity before the delivery location. Constraints (2f) and
(2g) respectively ensure that the total number of seats before and after adding or removing
new seats does not exceed the passenger capacity. Similarly, constraint (2h) ensures that the total number of passengers on the aircraft after fulfilling the demand of vertex *i* does not exceed the passenger capacity. Constraint (2i) enforces the relationship between $u_i$ and $y_i$. Note that the left hand side of the constraint restricts the picked-up cargo and passengers to not exceed the aircraft capacity. Constraints (2j) and (2k) define the upper bound and lower bound on the available cargo space, respectively. Similarly, constraints (2l) and (2m) set the upper and lower bounds on the number of empty seats. Constraint (2n) ensures that the passenger demand is met at each location, while constraint (2o) does the same thing for the cargo demand. Constraint (2p) restricts the number of the seats that can be added based on their availability. The lower bound on the number of removed seats, when $s_i < 0$, is always the number of seats on the aircraft at the arrival of location *i*. Lastly, constraints (2q) - (2s) specify binary and continuous variable domains.

### 4.2.2 Three-indexed Rank-Based MIP ($MIP_{rank}$)

The three-indexed ranked-based MIP model for PD-SRP ($MIP_{rank}$) is adapted from a model for the multi-commodity pickup and delivery traveling salesperson problem [3]. In $MIP_{rank}$, $z_{i,j}^t$ is a binary variable indicating that aircraft goes directly from location *i* to location *j* and location *i* is at position *t* of the tour, for $i, j \in V_{0,N+1}$, $i \neq j$, $t \in \{0, ..., 2n+1\}$. Binary variable $y_{i,t}$ is 1 if location *i* is visited at position *t* of the tour, $i \in V_{0,N+1}$, $t \in \{0, ..., 2n+1\}$ and 0 otherwise. Variable $s_t$ is the number of seats added or removed at the *t*-th position of the tour, for $t \in \{0, ..., 2n+1\}$, with a negative value corresponding to the number of seats removed. Variables $w_t$ and $u_t$ represent the empty seats and available cargo space on arrival at *t*-th position of the tour, for $t \in \{0, ..., 2n+1\}$. Finally, let $\pi_t$ and $q_t$ be the number of passengers and boxes of cargo on arrival at *t*-th position of the tour. The $MIP_{rank}$ is presented in Figure 4.

The objective function minimizes the total travel distance. Constraints (3a) and (3b) ensure that tour positions are assigned to exactly one location and that each location is visited exactly once, respectively. Constraint (3c) calculates the number of empty seats just before visit *t*, where $\sum_{i=1}^n y_{i,t-1}\hat{n}_i$ is the number of passengers picked up at position *t* - 1 of the tour. Similarly, constraint (3d) calculates the available cargo space just before visit *t*, where $\sum_{i=1}^n y_{i,t}\hat{q}_i$ is the amount of cargo picked up at position *t* of the tour. Constraint (3e) states that each commodity is picked up before it is delivered. Constraint (3f) enforces the relationship between $w_t$ and $u_t$. The left hand side of the constraint enforces that the picked-up cargo and passengers do not exceed the available aircraft capacity. Constraint (3g) ensures that there is always $\hat{C}$ space available for cargo on the aircraft. From (3f) and (3g) we can conclude that $\pi_t + w_t \le \hat{S}$: the total number of seats does not exceed the passenger capacity. Constraint (3h) ensures the feasibility of the number of seats to be added or removed. Constraints (3i) and (3j) calculate the number of passengers and boxes of cargo at each position of the tour, respectively. Constraints (3k) and (3l) enforce the relationship between *y* and *z* variables and, together with (3e) and (3m), prevent the formation of subtours in a MTZ fashion. Lastly, constraints (3n) - (3q) specify the domains of the variables.

### 4.2.3 Upper bound MIP Model ($MIP_{UB}$)

Our preliminary experiments suggested that the CP and MIP models presented above struggled to find good feasible solutions. We, therefore, investigate the use of a third MIP model, designed to quickly find an upper bound on the PD-SRP by solving a restriction of the full problem. Such a model provides a heuristic solution as well as a potential warm-start solution for the complete models.
$$
\begin{array}{ll}
\min & \displaystyle \sum_{t=0}^{2n} \sum_{i \in V_{0,N+1}} \sum_{j \in V_{0,N+1}, j \neq i} d_{i,j} z_{i,j}^t \\
& \qquad (\text{MIP}_{\text{rank}}) \\
& \displaystyle \sum_{i \in V_0} y_{i,t} = 1 & t \in \{0, \dots, 2n\} \qquad (3a) \\
& \displaystyle \sum_{t=1}^{2n} y_{i,t} = 1 & i \in V \qquad (3b) \\
& w_t = w_{t-1} + s_{t-1} - \sum_{i=1}^{n} y_{i,t-1} \hat{\pi}_i & t \in \{1, \dots, 2n+1\} \qquad (3c) \\
& u_t = u_{t-1} - L s_{t-1} - \sum_{i=1}^{n} y_{i,t-1} \hat{q}_i & t \in \{1, \dots, 2n+1\} \qquad (3d) \\
& \displaystyle \sum_{t=1}^{n} t y_{i,t} - \sum_{t=1}^{n} t y_{n+i,t} \le -1 & i \in P \qquad (3e) \\
& q_t + u_t + L w_t + L \pi_t \le K & t \in \{0, \dots, 2n+1\} \qquad (3f) \\
& q_t + u_t \ge \hat{C} & t \in \{0, \dots, 2n+1\} \qquad (3g) \\
& -w_t \le s_t \le \min (\hat{S}, \sum_{i=1}^{n} S_i^0 y_{i,t}) & t \in \{0, \dots, 2n+1\} \qquad (3h) \\
& \pi_t = \pi_{t-1} + \sum_{i=1}^{n} y_{i,t-1} \hat{\pi}_i & t \in \{1, \dots, 2n+1\} \qquad (3i) \\
& q_t = q_{t-1} + \sum_{i=1}^{n} y_{i,t-1} \hat{q}_i & t \in \{1, \dots, 2n+1\} \qquad (3j) \\
& y_{i,t} - \sum_{j=0}^{n} z_{i,j}^t = 0 & i \in V_0, t \in \{0, \dots, 2n\} \qquad (3k) \\
& y_{j,t} - \sum_{i=0}^{n} z_{i,j}^{t-1} = 0 & j \in V_{0,N+1}, t \in \{1, \dots, 2n+1\} \qquad (3l) \\
& y_{0,0} = y_{2n+1,2n+1} = 1, y_{0,t} = 0 & t \in \{1, \dots, 2n\} \qquad (3m) \\
& s_t \le \hat{S}, w_t \le \hat{S}, u_t \le K & t \in \{0, \dots, 2n+1\} \qquad (3n) \\
& u_0 = \tilde{C}_0, w_0 = \tilde{S}_0, \pi_0 = q_0 = 0 & t \in \{0, \dots, 2n+1\} \qquad (3o) \\
& y_{i,t} \in \{0,1\}, z_{i,j}^t \in \{0,1\} & i, j \in V_{0,N+1}, t \in \{0, \dots, 2n+1\} \qquad (3p) \\
& u_t, w_t, \pi_t, q_t \in \mathbb{R}^{0+}, s_t \in \mathbb{R} & t \in \{0, \dots, 2n+1\} \qquad (3q)
\end{array}
$$

<caption>Figure 4 A Three-Indexed MIP Model for the PD-SRP.</caption>