MaxSAT solvers (e.g., CASHWMaxSAT-CorePlus [19]) can be directly applied to solving PBO, we mainly focus on evaluating the performance of specialized solvers for PBO. The codes of all competitors were kindly provided by the authors. As for LS-PBO, RoundingSat and PBO-IHS, we employ the default parameters in the corresponding literature, respectively. Our code will be made publicly available. All algorithms are implemented in C++ and compiled by g++ with -O3 option. All the algorithms are run on Intel Xeon Gold 6238 CPU @ 2.10GHz with 512GB RAM under CentOS 7.9.

For the application benchmarks, our proposed algorithms and LS-PBO are both run 20 times whose seed is from 1 to 20 on each instance, whereas the exact solvers are run once on each instance. For the other four standard benchmarks (i.e., PB2016, MIPLIB, CRAFT and KNAP), following the settings of previous works [10, 27], all the algorithms are run only once on each instance. We test the algorithms with a time limit of 3600 seconds.

For the application benchmarks, we use *min* to denote the best solution value found and *avg* to denote the average value over the 20 runs. For all the benchmarks, we report the number of instances where the algorithm finds the best solution value among all algorithms, denoted by *#win*. There are some unsatisfied instances in the PB2016 benchmark. RoundingSat and PBO-IHS can guarantee the optimality of the solutions they obtain and thus can prove some of these unsatisfied instances, whereas DeepOpt-PBO and LS-PBO cannot do it because these two algorithms belong to incomplete algorithms. For the above case, following the similar method from the literature [18], if all the algorithms fail to obtain any feasible solution for such an unsatisfied instance, then *#win* value of all the algorithms for this instance needs to be increased by 1. The bold value indicates the best solution value obtained by all the algorithms. In detail, the bold value of each *avg* column indicates the best average solution when some algorithms obtain the same minimal solution values. For one instance, if only one algorithm finds the best minimal solution value, only its corresponding *min* column should be marked.

## 6.2 Experimental Results

Note that for the application benchmarks, two exact solvers RoundingSat and PBO-IHS can obtain the same best solution as our proposed algorithms for only 4 instances. Thus, for the sake of space, we do not report the detailed results of these two exact solvers. We mainly compare DeepOpt-PBO-v1 and DeepOpt-PBO-v2 with LS-PBO. The experimental results on the application benchmarks are presented in Tables 2–4. According to our results, our proposed algorithms are consistently superior on the MWCB and SAP benchmarks. For the WSNO benchmark, two proposed algorithms and LS-PBO can obtain the same best solution values. Furthermore, LS-PBO can find a minimal average solution for 15 instances, while our proposed algorithms do it for only 6 instances. Our proposed algorithms adopt some perturbation mechanisms, which may make them difficult for the algorithms to steadily obtain a good solution.

Table 5 gives a summary on all the benchmarks. According to the experimental results, DeepOpt-PBO-v2 outperforms other algorithms in terms of obtaining more optimal solution values on PB2016, MIPLIB and CRAFT, except for the benchmark KNAP where PBO-IHS has a better performance. In addition, under a given time limit, we observe that PBO-IHS can prove the optimality of a solution for more instances compared to RoundingSat in all benchmarks. Although our proposed algorithms cannot prove the optimality due to the natural property of local search, DeepOpt-PBO-v2 can find more minimal solution values overall. Because some instances from the above benchmarks have different kinds of problem structures and all the local search algorithms are run only once on each instance, the restart
41:14 Improving Local Search for Pseudo Boolean Optimization

<table><caption>Table 2 Experiment results on MWCB.</caption><thead><tr><th rowspan="2">Instance</th><th colspan="2">DeepOpt-PBO-v1</th><th colspan="2">DeepOpt-PBO-v2</th><th colspan="2">LS-PBO</th></tr><tr><th>min</th><th>avg</th><th>min</th><th>avg</th><th>min</th><th>avg</th></tr></thead><tbody><tr><td>100_200_90</td><td><strong>103363</strong></td><td>104188.35</td><td>103430</td><td>104233</td><td>110450</td><td>111314.25</td></tr><tr><td>1000_250_90</td><td><strong>140223</strong></td><td>141182.2</td><td>140237</td><td>141194.9</td><td>148181</td><td>149828.65</td></tr><tr><td>1200_200_90</td><td><strong>104063</strong></td><td><strong>104572.25</strong></td><td><strong>104063</strong></td><td>104589.1</td><td>110993</td><td>112897.1</td></tr><tr><td>1200_250_90</td><td><strong>141211</strong></td><td>142310.6</td><td>141238</td><td>142343.2</td><td>150212</td><td>152888.7</td></tr><tr><td>1400_200_90</td><td><strong>103948</strong></td><td>104867.35</td><td>104057</td><td>104901.55</td><td>110792</td><td>112975.65</td></tr><tr><td>1400_250_90</td><td><strong>141567</strong></td><td>142675.7</td><td>141746</td><td>142722.15</td><td>150981</td><td>152932.9</td></tr><tr><td>1600_200_90</td><td><strong>119226</strong></td><td><strong>120182.4</strong></td><td><strong>119226</strong></td><td>120215.9</td><td>136944</td><td>143371.9</td></tr><tr><td>1600_250_90</td><td><strong>162651</strong></td><td>163893.3</td><td>162656</td><td>163937.1</td><td>183797</td><td>196547.75</td></tr><tr><td>1800_200_90</td><td><strong>203135</strong></td><td><strong>205888.75</strong></td><td><strong>203135</strong></td><td>205959.2</td><td>219536</td><td>224144.95</td></tr><tr><td>1800_250_90</td><td><strong>253073</strong></td><td>257501.75</td><td>253078</td><td>257715.25</td><td>276336</td><td>283192.95</td></tr><tr><td>2000_200_90</td><td><strong>227294</strong></td><td>229591.95</td><td>227424</td><td>229676.1</td><td>246109</td><td>250958.4</td></tr><tr><td>2000_250_90</td><td><strong>286970</strong></td><td><strong>289889.4</strong></td><td><strong>286970</strong></td><td>290058.1</td><td>309645</td><td>314528.5</td></tr><tr><td>1000_200_95</td><td><strong>113384</strong></td><td>114749.8</td><td>113501</td><td>114855.95</td><td>117064</td><td>117945.7</td></tr><tr><td>1000_250_95</td><td><strong>151882</strong></td><td>153581.05</td><td>152007</td><td>153622</td><td>156543</td><td>157976.1</td></tr><tr><td>1200_200_95</td><td><strong>114585</strong></td><td>115920.2</td><td>114675</td><td>115975</td><td>118045</td><td>119544.4</td></tr><tr><td>1200_250_95</td><td><strong>153104</strong></td><td>155765.8</td><td>153138</td><td>155798.35</td><td>159310</td><td>161454.1</td></tr><tr><td>1400_200_95</td><td><strong>114195</strong></td><td><strong>115231.3</strong></td><td><strong>114195</strong></td><td>115305.95</td><td>118913</td><td>119779.4</td></tr><tr><td>1400_250_95</td><td><strong>155158</strong></td><td>156149.65</td><td>155174</td><td>156180.15</td><td>161658</td><td>162917.95</td></tr><tr><td>1600_200_95</td><td><strong>168271</strong></td><td><strong>171985.9</strong></td><td><strong>168271</strong></td><td>172024.6</td><td>185707</td><td>190763.55</td></tr><tr><td>1600_250_95</td><td><strong>215266</strong></td><td>218013.5</td><td>215316</td><td>218069.55</td><td>236547</td><td>244060.35</td></tr><tr><td>1800_200_95</td><td><strong>238699</strong></td><td>241690.55</td><td>238702</td><td>241779.85</td><td>251744</td><td>256988.75</td></tr><tr><td>1800_250_95</td><td><strong>297981</strong></td><td><strong>302487.3</strong></td><td><strong>297981</strong></td><td>302927.25</td><td>314968</td><td>318961.25</td></tr><tr><td>2000_200_95</td><td><strong>257696</strong></td><td><strong>262062.9</strong></td><td><strong>257696</strong></td><td>262113.6</td><td>272832</td><td>277406.7</td></tr><tr><td>2000_250_95</td><td><strong>324379</strong></td><td><strong>328952.6</strong></td><td><strong>324379</strong></td><td>329032.9</td><td>340859</td><td>346007.95</td></tr></tbody></table>

<table><caption>Table 3 Experiment results on WSNO.</caption><thead><tr><th rowspan="2">Instance</th><th colspan="2">DeepOpt-PBO-v1</th><th colspan="2">DeepOpt-PBO-v2</th><th colspan="2">LS-PBO</th></tr><tr><th>min</th><th>avg</th><th>min</th><th>avg</th><th>min</th><th>avg</th></tr></thead><tbody><tr><td>100_40_4</td><td><strong>210</strong></td><td><strong>210</strong></td><td><strong>210</strong></td><td><strong>210</strong></td><td><strong>210</strong></td><td><strong>210</strong></td></tr><tr><td>150_60_4</td><td><strong>602</strong></td><td>612.85</td><td><strong>602</strong></td><td><strong>602</strong></td><td><strong>602</strong></td><td>602.2</td></tr><tr><td>200_80_4</td><td><strong>715</strong></td><td>719.45</td><td><strong>715</strong></td><td>716.65</td><td><strong>715</strong></td><td><strong>715.1</strong></td></tr><tr><td>250_100_4</td><td><strong>1305</strong></td><td>1401.2</td><td><strong>1305</strong></td><td>1330.05</td><td><strong>1305</strong></td><td><strong>1305</strong></td></tr><tr><td>300_120_4</td><td><strong>1257</strong></td><td>1330.25</td><td><strong>1257</strong></td><td>1373</td><td><strong>1257</strong></td><td><strong>1257.05</strong></td></tr><tr><td>350_140_4</td><td><strong>1737</strong></td><td>1957.4</td><td><strong>1737</strong></td><td>1997.95</td><td><strong>1737</strong></td><td><strong>1744.05</strong></td></tr><tr><td>400_160_4</td><td><strong>2240</strong></td><td>2509.55</td><td>2241</td><td>2598.85</td><td><strong>2240</strong></td><td><strong>2240.5</strong></td></tr><tr><td>450_180_4</td><td><strong>1869</strong></td><td>2780.7</td><td>1878</td><td>2598.7</td><td><strong>1869</strong></td><td><strong>1889.25</strong></td></tr><tr><td>500_200_4</td><td><strong>2577</strong></td><td>3676.8</td><td>2674</td><td>3637.95</td><td><strong>2577</strong></td><td><strong>2616.2</strong></td></tr><tr><td>100_40_6</td><td><strong>140</strong></td><td><strong>140</strong></td><td><strong>140</strong></td><td><strong>140</strong></td><td><strong>140</strong></td><td><strong>140</strong></td></tr><tr><td>150_60_6</td><td><strong>402</strong></td><td><strong>402</strong></td><td><strong>402</strong></td><td><strong>402</strong></td><td><strong>402</strong></td><td><strong>402</strong></td></tr><tr><td>200_80_6</td><td><strong>477</strong></td><td>480</td><td><strong>477</strong></td><td><strong>477.05</strong></td><td><strong>477</strong></td><td><strong>477.7</strong></td></tr><tr><td>250_100_6</td><td><strong>870</strong></td><td>893.35</td><td><strong>870</strong></td><td><strong>870.1</strong></td><td><strong>870</strong></td><td>870.5</td></tr><tr><td>300_120_6</td><td><strong>839</strong></td><td>866.55</td><td><strong>839</strong></td><td>862.15</td><td><strong>839</strong></td><td><strong>839.3</strong></td></tr><tr><td>350_140_6</td><td><strong>1158</strong></td><td>1267.7</td><td><strong>1158</strong></td><td>1288.25</td><td><strong>1158</strong></td><td><strong>1158.85</strong></td></tr><tr><td>400_160_6</td><td><strong>1493</strong></td><td>1671.8</td><td><strong>1493</strong></td><td>1656.5</td><td><strong>1493</strong></td><td><strong>1494.25</strong></td></tr><tr><td>450_180_6</td><td><strong>1246</strong></td><td>1588.45</td><td>1265</td><td>1641.6</td><td><strong>1246</strong></td><td><strong>1247.8</strong></td></tr><tr><td>500_200_6</td><td><strong>1718</strong></td><td>1984.05</td><td><strong>1718</strong></td><td>1927</td><td><strong>1718</strong></td><td><strong>1727.65</strong></td></tr></tbody></table>

strategy plays a key role in the performance of DeepOpt-PBO-v2. To sum up, for three application benchmarks, PB2016, MIPLIB, and CRAFT, our algorithm totally dominates all the competitors, whereas PBO-IHS performs better than other PBO solvers for the KNAP benchmark.

Additionally, we evaluate the performance of all the solvers on seven benchmarks using performance profile [11]. As shown in Figure 2, the plot captures the probability of reaching a fixed quality in a time, at most a factor $\tau$ slower than the optimal algorithm. Specially,
<table><caption>Table 4 Experiment results on SAP.</caption><thead><tr><th rowspan="2">Instance</th><th colspan="2">DeepOpt-PBO-v1</th><th colspan="2">DeepOpt-PBO-v2</th><th colspan="2">LS-PBO</th></tr><tr><th>min</th><th>avg</th><th>min</th><th>avg</th><th>min</th><th>avg</th></tr></thead><tbody><tr><td>100</td><td><strong>579</strong></td><td><strong>579.7</strong></td><td><strong>579</strong></td><td>579.8</td><td>580</td><td>583</td></tr><tr><td>110</td><td><strong>618</strong></td><td><strong>618.75</strong></td><td><strong>618</strong></td><td>618.85</td><td>619</td><td>626.05</td></tr><tr><td>120</td><td><strong>675</strong></td><td><strong>677.15</strong></td><td><strong>675</strong></td><td>677.8</td><td>679</td><td>685.6</td></tr><tr><td>130</td><td><strong>733</strong></td><td>736.9</td><td>734</td><td>738.4</td><td>738</td><td>744.6</td></tr><tr><td>140</td><td><strong>750</strong></td><td><strong>751.9</strong></td><td><strong>750</strong></td><td>753.75</td><td>757</td><td>763.65</td></tr><tr><td>150</td><td><strong>803</strong></td><td><strong>808.4</strong></td><td><strong>803</strong></td><td>810.3</td><td>821</td><td>828.45</td></tr><tr><td>160</td><td><strong>849</strong></td><td><strong>854.85</strong></td><td><strong>849</strong></td><td>855.75</td><td>867</td><td>872.75</td></tr><tr><td>170</td><td><strong>880</strong></td><td><strong>884.4</strong></td><td><strong>880</strong></td><td>885.55</td><td>897</td><td>907.5</td></tr><tr><td>180</td><td><strong>939</strong></td><td>948.45</td><td>941</td><td>951.7</td><td>971</td><td>977.4</td></tr><tr><td>190</td><td><strong>965</strong></td><td>973.2</td><td>970</td><td>977</td><td>996</td><td>1005.25</td></tr><tr><td>200</td><td><strong>1029</strong></td><td><strong>1040</strong></td><td><strong>1029</strong></td><td>1042.55</td><td>1067</td><td>1073.4</td></tr><tr><td>210</td><td><strong>1067</strong></td><td>1074.7</td><td>1068</td><td>1077.9</td><td>1094</td><td>1112.75</td></tr><tr><td>220</td><td><strong>1114</strong></td><td><strong>1127.3</strong></td><td><strong>1114</strong></td><td>1129.6</td><td>1151</td><td>1163.1</td></tr><tr><td>230</td><td><strong>1151</strong></td><td>1166.85</td><td>1162</td><td>1169.7</td><td>1195</td><td>1205.75</td></tr><tr><td>240</td><td><strong>1179</strong></td><td>1192.6</td><td>1183</td><td>1195.5</td><td>1219</td><td>1234.55</td></tr><tr><td>250</td><td><strong>1235</strong></td><td>1243.1</td><td>1237</td><td>1245.5</td><td>1274</td><td>1290.55</td></tr><tr><td>260</td><td><strong>1275</strong></td><td>1286.15</td><td>1277</td><td>1287.45</td><td>1318</td><td>1334.9</td></tr><tr><td>270</td><td><strong>1344</strong></td><td>1353.1</td><td>1348</td><td>1356.1</td><td>1392</td><td>1403.55</td></tr><tr><td>280</td><td><strong>1348</strong></td><td>1370.75</td><td>1354</td><td>1375.35</td><td>1407</td><td>1424.3</td></tr><tr><td>290</td><td><strong>1403</strong></td><td>1416.35</td><td>1405</td><td>1419.5</td><td>1448</td><td>1471.15</td></tr><tr><td>300</td><td><strong>1471</strong></td><td>1491.5</td><td>1480</td><td>1496.1</td><td>1538</td><td>1548.4</td></tr></tbody></table>

<table><caption>Table 5 Summary results of comparing our proposed algorithms to its competitors on all the benchmarks. #inst denotes the number of instances in each benchmark.</caption><thead><tr><th>Benchmark</th><th>#inst</th><th>DeepOpt-PBO-v1</th><th>DeepOpt-PBO-v2</th><th>LS-PBO</th><th>RoundingSat</th><th>PBO-IHS</th></tr><tr><th></th><th></th><th>#win</th><th>#win</th><th>#win</th><th>#win</th><th>#win</th></tr></thead><tbody><tr><td>MWCB</td><td>24</td><td><strong>24</strong></td><td>9</td><td>0</td><td>0</td><td>0</td></tr><tr><td>WSNO</td><td>18</td><td><strong>18</strong></td><td>14</td><td><strong>18</strong></td><td>4</td><td>4</td></tr><tr><td>SAP</td><td>21</td><td><strong>21</strong></td><td>9</td><td>0</td><td>0</td><td>0</td></tr><tr><td>PB2016</td><td>1600</td><td>1141</td><td><strong>1226</strong></td><td>1060</td><td>1195</td><td>1101</td></tr><tr><td>MIPLIB</td><td>267</td><td>163</td><td><strong>168</strong></td><td>135</td><td>125</td><td>130</td></tr><tr><td>CRAFT</td><td>1025</td><td>903</td><td><strong>930</strong></td><td>878</td><td>882</td><td>917</td></tr><tr><td>KNAP</td><td>783</td><td>695</td><td>693</td><td>649</td><td>392</td><td><strong>778</strong></td></tr><tr><td>Total</td><td>3738</td><td>2965</td><td><strong>3049</strong></td><td>2740</td><td>2598</td><td>2930</td></tr></tbody></table>

<table><caption>Table 6 Comparative results for DeepOpt-PBO-v1 and its modified versions with different strategies on the application benchmarks. #better and #worse denote the number of instances where DeepOpt-PBO-v1 obtains better and worse solution values, respectively.</caption><thead><tr><th rowspan="2">Benchmark</th><th rowspan="2">#inst</th><th colspan="2">vs. v1+nohh</th><th colspan="2">vs. v1+nodo</th><th colspan="2">vs. v1+nosf</th></tr><tr><th>#better</th><th>#worse</th><th>#better</th><th>#worse</th><th>#better</th><th>#worse</th></tr></thead><tbody><tr><td>MWCB</td><td>24</td><td><strong>11</strong></td><td>8</td><td><strong>24</strong></td><td>0</td><td><strong>24</strong></td><td>0</td></tr><tr><td>WSNO</td><td>18</td><td><strong>1</strong></td><td>0</td><td><strong>1</strong></td><td>0</td><td>0</td><td>0</td></tr><tr><td>SAP</td><td>21</td><td><strong>11</strong></td><td>5</td><td><strong>18</strong></td><td>0</td><td><strong>20</strong></td><td>0</td></tr><tr><td>Total</td><td>63</td><td><strong>23</strong></td><td>13</td><td><strong>43</strong></td><td>0</td><td><strong>44</strong></td><td>0</td></tr></tbody></table>

when $\tau$ equals 1, we obtain the probability that the algorithm is the fastest. Figure 2 (a) demonstrates the superior performance of DeepOpt-PBO-v1 on the first three application benchmarks. It can be observed that DeepOpt-PBO-v1 has a higher probability of finding the optimal solution for each $\tau$. In Figure 2 (b), LS-PBO performs the best within a short time scale. However, DeepOpt-PBO-v2 surpasses LS-PBO around $\tau = 6$ and maintains its leadership consistently.
41:16 Improving Local Search for Pseudo Boolean Optimization

<figure><img src="image_422.png" /><figcaption>Figure 2 Performance profiles for DeepOpt-PBO and all competitors for reaching the best solution on the three application benchmarks (a) and four standard benchmarks (b).</figcaption></figure>

For three application benchmarks, to determine better solution values, we increase the execution time to 5400 seconds, and once again our algorithms perform better than LS-PBO. For the PB2016 benchmark, when increasing the execution time to 5400 seconds, although the solution values obtained by RoundingSat and PBO-HIS are better than before, our algorithm can still find more best solution values than RoundingSat and PBO-IHS for 29 and 140 instances, respectively.

To verify the effectiveness of the proposed strategies, we compare DeepOpt-PBO-v1 with three alternative versions: 1) v1+nohh utilizes the *age* strategy instead of *hhscore*; 2) v1+nodo does not use the DeepOpt method; 3) v1+nosf does not employ the *SampleFlip* function. The results in Table 6 demonstrate that all proposed strategies are effective. In addition, we randomly select 20 instances for each standard benchmark. In total, 80 instances are picked. We also test the performance of our proposed algorithm and three alternative versions on these picked instances. All the algorithms are run 20 times on each instance. Once again, the results show that our proposed algorithm obviously performs better than three alternative versions.

Here, we give a discussion to compare our algorithm with MIP solvers and MaxSAT solvers. First, although MaxSAT solvers (e.g., CASHWMaxSAT-CorePlus [19]) can be directly applied to solving PBO, Lei et al. [18] observed that several of the instances from some benchmarks, such as PB2016, are too large to admit practical encodings into MaxSAT. Thus, MaxSAT solvers have usually poor performance for PBO. Second, more general solvers (e.g., MIP solvers) could also be used for comparison, but in our work, we mainly focus on evaluating the performance of specialized solvers for PBO. We have also tested the performance of the non-commercial MIP solver SCIP [4] on all the seven benchmarks. In detail, for the KNAP benchmark, the performance of PBO-IHS and SCIP is better than our algorithm. But, for the remaining six benchmarks, our algorithms perform better than SCIP.

7 Conclusion

In this paper, we propose a two-level selection strategy, a novel deep optimization strategy, and a sampling flipping method. Based on the above strategies and some other tricks, we develop two local search algorithms. Experiments show that the proposed algorithms significantly outperform the state-of-the-art PBO algorithms.
References

1. Carlos Ansotegui, Fahiem Bacchus, Matti Järvisalo, and Ruben Martins. MaxSAT evaluation 2017: Solver and benchmark descriptions, 2017.

2. Jeremias Berg, Emilia Oikarinen, Matti Järvisalo, and Kai Puolamäki. Minimum-width confidence bands via constraint optimization. In *Proceedings of the Twenty-Third Principles and Practice of Constraint Programming*, volume 10416, pages 443–459, 2017.

3. Daniel Le Berre and Romain Wallon. On dedicated CDCL strategies for PB solvers. In *Proceedings of the Twenty-Fourth Theory and Applications of Satisfiability Testing*, volume 12831, pages 315–331, 2021.

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

5. Frédéric Boussemart, Fred Hemery, Christophe Lecoutre, and Lakhdar Sais. Boosting systematic search by weighting constraints. In *Proceedings of the Sixteenth European Conference on Artificial Intelligence*, volume 16, pages 146–150, 2004.

6. Shaowei Cai, Jinkun Lin, and Chuan Luo. Finding a small vertex cover in massive sparse graphs: Construct, local search, and preprocess. *Journal of Artificial Intelligence Research*, 59:463–494, 2017.

7. Shaowei Cai and Kaile Su. Comprehensive score: Towards efficient local search for SAT with long clauses. In *Proceedings of the Twenty-Third International Joint Conference on Artificial Intelligence*, pages 489–495, 2013.

8. Jiejiang Chen, Shaowei Cai, Yiyuan Wang, Wenhao Xu, Jia Ji, and Minghao Yin. Improved local search for the minimum weight dominating set problem in massive graphs by using a deep optimization mechanism. *Artificial Intelligence*, 314:103819, 2023.

9. Jo Devriendt, Ambros Gleixner, and Jakob Nordström. Learn to relax: Integrating 0-1 integer linear programming with pseudo-Boolean conflict-driven search. *Constraints*, 26(1):26–55, 2021.

10. Jo Devriendt, Stephan Gocht, Emir Demirović, Jakob Nordström, and Peter J. Stuckey. Cutting to the core of pseudo-Boolean optimization: Combining core-guided search with cutting planes reasoning. In *Proceedings of the Thirty-Fifth AAAI Conference on Artificial Intelligence*, volume 35, pages 3750–3758, 2021.

11. Elizabeth D. Dolan and Jorge J. Moré. Benchmarking optimization software with performance profiles. *Mathematical Programming*, 91:201–213, 2002.

12. Jan Elffers and Jakob Nordström. Divide and conquer: Towards faster pseudo-Boolean solving. In *Proceedings of the Twenty-Seventh International Joint Conference on Artificial Intelligence*, pages 1291–1299, 2018.

13. Jan Elffers and Jakob Nordström. A cardinal improvement to pseudo-Boolean solving. In *Proceedings of the Thirty-Fourth AAAI Conference on Artificial Intelligence*, volume 34, pages 1495–1503, 2020.

14. John Franco and John Martin. A history of satisfiability. In *Handbook of Satisfiability*, volume 185 of *Frontiers in Artificial Intelligence and Applications*, pages 3–74. 2009.

15. Pascal Van Hentenryck and Laurent Michel. *Constraint-Based Local Search*. The MIT press, 2009.

16. Gergely Kovásznai, Balázs Erdélyi, and Csaba Biró. Investigations of graph properties in terms of wireless sensor network optimization. In *Proceedings of the IEEE International Conference on Future IoT Technologies*, pages 1–8, 2018.

17. Gergely Kovásznai, Krisztián Gajdár, and Laura Kovács. Portfolio SAT and SMT solving of cardinality constraints in sensor network optimization. In *Proceedings of the Twenty-First International Symposium on Symbolic and Numeric Algorithms for Scientific Computing*, pages 85–91, 2019.
41:18 Improving Local Search for Pseudo Boolean Optimization

18 Zhendong Lei, Shaowei Cai, Chuan Luo, and Holger H. Hoos. Efficient local search for pseudo Boolean optimization. In Proceedings of the Twenty-Fourth Theory and Applications of Satisfiability Testing, volume 12831, pages 332–348, 2021.

19 Zhendong Lei, Yiyuan Wang, Shiwei Pan, Shaowei Cai, and Minghao Yin. CASHWMaxSAT-CorePlus: Solver description. MaxSAT Evaluation, page 8, 2022.

20 Manuel López-Ibáñez, Jérémie Dubois-Lacoste, Leslie Pérez Cáceres, Mauro Birattari, and Thomas Stützle. The irace package: Iterated racing for automatic algorithm configuration. Operations Research Perspectives, 3:43–58, 2016.

21 Ruben Martins, Vasco Manquinho, and Inês Lynce. Open-WBO: A modular maxsat solver. In Proceedings of the Seventeenth Theory and Applications of Satisfiability Testing, volume 8561, pages 438–445, 2014.

22 David Pisinger. Where are the hard knapsack problems? Computers & Operations Research, 32(9):2271–2284, 2005.

23 Olivier Roussel and Vasco Manquinho. Pseudo-Boolean and cardinality constraints. In Handbook of Satisfiability, volume 336 of Frontiers in Artificial Intelligence and Applications, pages 1087–1129. 2021.

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

25 Paul Shaw. Using constraint programming and local search methods to solve vehicle routing problems. In Proceedings of the Fourth Principles and Practice of Constraint Programming, volume 1520, pages 417–431, 1998.

26 Pavel Smirnov, Jeremias Berg, and Matti Järvisalo. Pseudo-Boolean optimization by implicit hitting sets. In Proceedings of the Twenty-Seventh Principles and Practice of Constraint Programming, volume 210, pages 51:1–51:20, 2021.

27 Pavel Smirnov, Jeremias Berg, and Matti Järvisalo. Improvements to the implicit hitting set approach to pseudo-Boolean optimization. In Proceedings of the Twenty-Fifth International Conference on Theory and Applications of Satisfiability Testing, volume 236, pages 13:1–13:18, 2022.

28 Zhouxing Su, Qingyun Zhang, Zhipeng Lü, Chu-Min Li, Weibo Lin, and Fuda Ma. Weighting-based variable neighborhood search for optimal camera placement. In Proceedings of the Thirty-Fifth AAAI Conference on Artificial Intelligence, volume 35, pages 12400–12408, 2021.

29 Renato Tinós, Michal W. Przewozniczek, and Darrell Whitley. Iterated local search with perturbation based on variables interaction for pseudo-Boolean optimization. In Proceedings of the Genetic and Evolutionary Computation Conference, pages 296–304, 2022.

30 Marc Vinyals, Jan Elffers, Jesús Giráldez-Cru, Stephan Gocht, and Jakob Nordström. In between resolution and cutting planes: A study of proof systems for pseudo-Boolean SAT solving. In Proceedings of the Twenty-First Theory and Applications of Satisfiability Testing, volume 10929, pages 292–310, 2018.

31 Chris Voudouris and Edward Tsang. Partial constraint satisfaction problems and guided local search. In Proceedings of the Practical Application of Constraint Technology, pages 337–356, 1996.

32 Yiyuan Wang, Shaowei Cai, Jiejiang Chen, and Minghao Yin. Sccwalk: An efficient local search algorithm and its improvements for maximum weight clique problem. Artificial Intelligence, 280:103230, 2020.

33 Yiyuan Wang, Dantong Ouyang, Liming Zhang, and Minghao Yin. A novel local search for unicost set covering problem using hyperedge configuration checking and weight diversity. Science China Information Sciences, 60:062103, 2017.

34 Jiongzhi Zheng, Jianrong Zhou, and Kun He. Farsighted probabilistic sampling based local search for (weighted) partial maxsat. CoRR, abs/2108.09988, 2021.
# Predict-Then-Optimise Strategies for Water Flow Control

Vincent Barbosa Vaz 📉 ⬆️ ID

The University of Melbourne, Australia
the Australian Research Council OPTIMA ITTC, Melbourne, Australia

James Bailey 📉 ⬆️ ID

The University of Melbourne, Australia

Christopher Leckie 📉 ⬆️ ID

The University of Melbourne, Australia

Peter J. Stuckey 📉 ⬆️ ID

Monash University, Australia
the Australian Research Council OPTIMA ITTC, Melbourne, Australia

## Abstract

A pressure sewer system is a network of pump stations used to collect and manage sewage from individual properties that cannot be directly connected to the gravity driven sewer network due to the topography of the terrain. We consider a common scenario for a pressure sewer system, where individual sites collect sewage in a local tank, and then pump it into the gravity fed sewage network. Standard control systems simply wait until the local tank reaches (near) capacity and begin pumping out. Unfortunately such simple control usually leads to peaks in sewage flow in the morning and evening, corresponding to peak water usage in the properties. High peak flows require equalization basins or overflow systems, or larger capacity sewage treatment plants. In this paper we investigate combining prediction and optimisation to better manage peak sewage flows. We use simple prediction methods to generate realistic possible future scenarios, and then develop optimisation models to generate pumping plans that try to smooth out flows into the network. The solutions of these models create a policy for pumping out that is specialized to individual properties and which overall is able to substantially reduce peak flows.

**2012 ACM Subject Classification** Applied computing → Operations research; Mathematics of computing → Mixed discrete-continuous optimization

**Keywords and phrases** Water Flow Control, Optimization, Machine Learning

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

**Category** Short Paper

**Funding** This research was partially funded by the Australian Government through the Australian Research Council Industrial Transformation Training Centre in Optimisation Technologies, Integrated Methodologies, and Applications (OPTIMA), with industry partner South East Water, Project ID IC200100009.

## 1 Introduction

A Pressure Sewer System (PSS) is a network of pump stations used to collect and manage sewage from individual properties that cannot be directly connected to the classic gravity sewer network due to the topography of the terrain (gravity limitation). The sewerage gravitates to the pump station and is then pumped through a pressure main to a main sewer and on to wastewater treatment plants. This is illustrated in Figure 1a. A PSS is often composed of intermediate, bigger pump stations between sub parts of the whole network (Figure 1b). Pump stations collect household sewage from a sub part of the network and pump it to the main. Our

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

© Vincent Barbosa Vaz, James Bailey, Christopher Leckie, and Peter J. Stuckey;
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. 42; pp. 42:1–42:10

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

LIPICS

Leibniz International Proceedings in Informatics

Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl Publishing, Germany
<figure><img src="image_425.png" /><figcaption>Figure 1 (a) Pressure sewer unit. (b) Pressure sewer network.</figcaption></figure>

focus in this paper is to optimise the operation of individual pump stations to balance the over-
all load at the treatment plant. In this paper, we only consider an existing network managed by
our industry partner South East Water. The network that we consider is free of intermediate
pump stations. Extending the approach to handle intermediate pump stations would require
adding constraints to ensure the capacity of the pump station is not violated at any time.

We consider a classic scenario for PSS where volumes at the residential level are released without optimised control. When the collection tanks reach their capacity, water is released entirely into the network. This represents a considerable amount of sewage conveyed through the sewer network that must be handled by the treatment plant, resulting in stress on the treatment processes and increased capital costs of upsized pipe and pump networks. This simple control policy does not make good use of the network. Volumes can be retained in pressure sewer tanks at the residential level and selectively released in a way that can optimise the flows to provide network capacity increases for the current network, and improve operations of the downstream treatment plant.

Most of the water usage occurs in the morning and the evening peak loads, when people are home. This results in two, identifiable peaks of activity in the network. This sudden surge in activity represents a challenge for the water treatment processes that follow. Huge amounts of resources need to be deployed at the treatment plant at these times to cope with the volumes to be treated and the associated unpredictability. South East Water would like to flatten the input flow at the treatment plant to achieve greater plant operational efficiency. This is possible by leveraging the buffer capabilities of the collection tanks, assuming that each tank can be controlled individually.

**Current operational strategy.** The water industry has yet to integrate data-driven models and optimisation techniques to facilitate and control their processes in a more efficient and systematic way. Most water companies rely on operators' knowledge and experience to parameterise and control their network. In current operation each tank fills until it reaches its capacity (cut-in setpoint) and then fully (until cut-off setpoint) empties the tank. In a previous approach to tackling the problem of reducing maximum outflows on the network we investigated simply modifying the (cut-in) set points of the tanks to try to reduce homogeneous behaviour but this was not really successful. Without adjustment the tanks quickly reached a steady state where the usual morning and evening peak inflows again resulted in high outflow. Due to the location of the network in a holiday zone the set points needed to be adjusted frequently as usage behaviour fluctuated rapidly over weekends, summer etc.
In this paper we define an optimisation based approach to controlling the maximum outflow, by deciding in which time periods to empty the tanks. In order to flatten out peaks we need to have some idea of the future inflows into the tanks. Because patterns of water usage are quite distinct generating realistic time series of inflows is challenging with machine learning models. Hence we use simple historical sampling to generate plausible future inflows. We show how combining (simplistic) prediction models with an optimisation model to determine when to release sewage into the network from individual properties, can substantially reduce the peak flows in the network.

## 2 Problem Description

We model the problem of controlling a PSS system to reduce peak flows as a MIP. We discretize the control problem by working over a finite time horizon of $n$ discrete time steps, which for our experiments are always 1 hour long. Note that this discretization is fine enough that none of our historical data has examples where a tank is filled from empty in under an hour. Finer discretization would allow some further reduction in peak outflows, but we expect that we capture most of the possible reduction using 1 hour discretization. The parameters for the water flow control model are shown in Table 1. In each time step we decide whether to empty the tank at a particular site. In the default control mechanism, when we decide to pump out of a tank, we empty it. This has the least wear and tear on the pumping and control mechanism.

<table><caption>**Table 1** Parameters for the water flow control model.</caption><thead><tr><td>Description</td><td>Parameter</td></tr></thead><tbody><tr><td>Time horizon</td><td>T = 1..n</td></tr><tr><td>Set of tanks</td><td>S</td></tr><tr><td>Capacity of tank at site s (cut-in setpoint)</td><td>c<sub>s</sub></td></tr><tr><td>Minimum amount of water to be pumped out</td><td>m</td></tr><tr><td>Inflow to tank at site s during time period t</td><td>i<sub>s,t</sub></td></tr><tr><td>Tank level of site s tank at the beginning of the first time period</td><td>l<sub>s,0</sub></td></tr></tbody></table>

### 2.1 Direct Formulation

In this section, we propose a direct MIP formulation for the PSS problem, over a given set of tank sites $S$ and time horizon $T$. The control systems for the tanks have simple functionality, we can (re-)set minimum and maximum tank levels, or initiate a pump out of the tank, but not control exactly how much volume is pumped out of the tank. This leads to the important decision we must make – for each tank $s$ at what times $t$ should it be emptied? $X_{s,t} \in \{0,1\}$. The model makes use of auxiliary variables: $l_{s,t}$ the level of tank $s$ at (the end of) time $t$; and $o_{s,t}$ the volume of water pumped out of tank $s$ during time period $t$. The model is:

$$ \underset{t \in T}{\operatorname{minimize}} \max_{s \in S} \sum_{s \in S} o_{s,t} \\ l_{s,t} = (l_{s,t-1} + i_{s,t}) * (1 - X_{s,t}), \quad \forall s \in S, t \in T \quad (1) $$

$$ l_{s,t} \le c_s, \quad \forall s \in S, t \in T \quad (2) $$

$$ o_{s,t} = (l_{s,t-1} + i_{s,t}) * X_{s,t}, \quad \forall s \in S, t \in T \quad (3) $$

$$ X_{s,t} = 1 \rightarrow o_{s,t} \ge m, \quad \forall s \in S, t \in T \quad (4) $$

$$ X_{s,t} \in \{0,1\}, \quad \forall s \in S, t \in T $$
where Equation (1) computes the level $l_{s,t}$ in each tank $s$ at the end of time period $t$ (previous level plus inflows, unless emptied); Equation (2) ensures that each tank's level remains below tank capacity; Equation (3) computes the outflow $o_{s,t}$ from tank $s$ at time $t$; and Equation (4) ensures that if the tank is emptied there is a minimum volume present (to prevent accelerated wear and tear on the infrastructure). The objective is to minimize maximum outflow across the period considered. Note that each of the constraints, and the objective are linear or easy to linearise.

## 2.2 Packing Formulation

The model above straightforwardly models the problem, but can become challenging to solve as the problem size grows. Next we instead consider the inflow to tank $s$ at time $t$ as a fixed amount of flow, we then decide when this should be pumped out. By precomputation we can then specify simple constraints to enforce that the capacity of tank is not exceeded. This leaves a packing problem, deciding when each chunk of water is pumped into the network.

The tank inflows are aggregated by hours and constitute the items of the problem. We require that the inflows are pumped out in order of inflow. The inflow at time $t$, $i_{s,t}$, must be pumped out, no earlier than $t$, and not later than latest($s,t$) = min($\{t' - 1 | t \le t' \le n, \sum_{t \le i \le t'} i_{s,i} \ge c_s\} \cup \{n\}$) which would mean the tank capacity was exceeded, since no later flows can be pumped out before $i_{s,t}$. Note that we treat the starting tank level $l_{s,0}$ as an inflow at time 0, $i_{s,0} = l_{s,0}$. We can also define the last inflow that must be pumped out in the considered time period, last($s$) = min{$t | \sum_{t \le i \le n} i_{s,i} \le c_s$} - 1. And for each tank and time define below($s,t$) = max{$t' | t \le t' < \text{latest}(s,t)$, $\sum_{t \le i \le t'} i_{s,i} < m$} to be the latest time $t'$ such that if the tank is emptied at time $t$ the sum of inflows after it up to $t'$ is too small to empty.

The new decision variables $p_{s,t,t'}$ determine the time $t'$ when the inflow to $s$ at time $t$ is pumped out (including the original tank volume $t = 0$). The model is defined by:

$$ \begin{array}{ll} \underset{t \in T}{\operatorname{max}} & \sum_{s \in S, t' \le t} p_{s,t',t} \times i_{s,t} \\ & \sum_{t \le i \le \text{latest}(s,t)} p_{s,t,i} = 1 & s \in S, t \in 0..\text{last}(s) \quad (5) \\ & p_{s,0,0} = 0 & s \in S \end{array} $$

$$ \sum_{t \le i \le t'} p_{s,t,i} \ge p_{s,t+1,t'}, \quad \forall s \in S, t \in 1..n, t' \in t+1..\text{latest}(s,t) \quad (7) $$

$$ X_{s,t'} = 1 \rightarrow \sum_{t \le i \le t'} p_{s,t,i} \ge 1, \quad s \in S, t \in T, t' \in t..\text{latest}(s,t) \quad (8) $$

$$ \sum_{t' \le t \le \text{latest}(s,t')} p_{s,t',t} \ge 1 \rightarrow X_{s,t} = 1, \quad s \in S, t \in T \quad (9) $$

$$ X_{s,t} = 1 \rightarrow \sum_{t < i \le \text{below}(s,t)} X_{s,i} \le 0, \quad \forall s \in S, t \in T \quad (10) $$

$$ X_{s,t} \in \{0,1\}, \quad \forall s \in S, t \in T $$

$$ p_{s,t,t'} \in \{0,1\}, \quad \forall s \in S, t \in 0..n, t' \in t..\text{latest}(s,t) $$

where Equation (5) enforces we pump out each inflow (before last($s$)) exactly once; Equation (6) enforces that nothing is pumped out at time 0; Equation (7) enforces that the inflow to tank $s$ at time $t$ is pumped out no later than the time the inflow at time $t+1$ is pumped out, i.e. the inflows are pumped out in order; Equation (8) connects the emptied variables to the
pump out variables by requiring that if tank *s* is emptied at time *t*' then each inflow before *t*'
is pumped out by time *t*'; Equation (9) connects them in the other direction requiring that if
any inflow is pumped out of tank *s* at time *t* then *tank *s* is emptied at time *t*; and Equation
(10) enforces that if *tank *s* is emptied at time *t* then it is not emptied again until at least *m*
units of flow have entered the tank. The objective minimizes maximum outflow, computed
from the pumped out variables. Again the entire model is easy to linearise.

3 Predicted Water Usage Generation

Online optimisation problems can be augmented with a predictor that informs the model on
future instances. Simple predictors can sample the inputs or sample the distribution of the
inputs. More complex Machine Learning based predictors can learn from the distribution
of the inputs as the online problem is being solved. Research [5, 1, 2, 3] has demonstrated
that sampling the distribution of the inputs or providing estimates can significantly improve
the quality of the solution. In practice, not all optimisation problems have access to the
distribution of the inputs.

We wish to generate water usage predictions for each site. The collected water usage
comes from diverse households exhibiting different behavioural patterns. Care must be taken
when building a predictor to capture the seasonality and cycles in the data. Because of these
properties, building an individual predictor for each site and time instance is unlikely to
produce realistic distributions of water usage.

**3.1 Historical Sampling**

We use *historical sampling* as described by Bent and Van Hentenryck [4]. This generates
samples from past subsequences in the historical data. Despite its simplicity, historical
sampling captures the structure of the sequence while providing fast outcomes compared to
ML techniques that require training during the online algorithm. We adapt the algorithm to
sample historical data while preserving the structural periodic information of our samples. In
particular, we “retrieve” a prediction sequence on inflows for tank *s* for times *T* = 1..*n* from
the historical data sequence *S* of inflows for tank *s* including *d̄* days of data, by randomly
selecting starting position *t*′ which is the start of some day (since our experiments always
commence from the first hour of a day) in *S*, where *S*[*l*..*u*] returns the slice of sequence *S*
from index *l* to *u* inclusive.

<caption>**Algorithm 1 Historical sampling.**</caption>

1 historicalSample(S, n)
2   d̄ ← |S| div 24
3   t' ← 24 × RANDOM([0, d̄ − 1]) + 1
4   return S[t'..t' + n − 1]

4 Experimental Evaluation

In this section, we present the results of using proposed models on a set of benchmark
instances. The models are implemented in MiniZinc [8], a high-level and solver-independent
modelling language, allowing for fast experimentation across existing solvers (OR-tools,
Gurobi, CPLEX) without compromising on efficiency. All experiments were run on the same
<figure><img src="image_426.png" /><figcaption>Figure 2 (a) Original level data available from the SCADA system, and (b) the reconstructed level (from where we compute inflow data).</figcaption></figure>

machine which has an Apple M1 Pro 3.22 GHz CPU with 10 cores and 16 GB of RAM. All approaches were given a time limit of one hour per instance. The solver used was Gurobi version 9.5.2.

## 4.1 Benchmark Instances

The data is provided by our industry partner, South East Water, and corresponds to pressure sewer readings. This catchment has been selected as it is of reasonable scale and is free of infiltration. The data is collected, through a series of scripts, from the SCADA server and corresponds to 3 years (2019-2021) of historical readings from approximately 4200 individual households. A range of attributes can be extracted from the readings; we focus on the water level and pump activation. The network is free of intermediate pump stations.

In order to make different decisions, we first need to rollback any previous decisions to obtain the original system inputs. The tanks are not equipped with water meters to measure the inflows but have level sensors. From the water level, we can derive the inflows to be the difference between two positive consecutive readings. This is illustrated in Figures 2(a) and 2(b). We generated inflow data for each site for each hour of the day in the periods used for instance generation.

We cluster similar sites using the k-means algorithm where distance is defined as the sum of absolute differences over their inflow data. We observe the average inflow for each cluster. We determined four identifiable water usage profiles amongst the households: two diurnal/bimodal (with a morning and afternoon peak) and two uni-modal usages (with just a morning peak), at a time translation respectively. We combined the bimodal and uni-modal sites respectively. The average inflow for each cluster is shown in Figure 4 in the supplementary material.

We created 6 problem instances from our industry partner. These instances represent different levels of complexity. To ensure they are different we choose different kinds of flows. For each instance we pick $|S|$ sites to use, either from the unimodal clusters, the bimodal clusters, or the complete set of clusters. We choose a number of hours $n$ to solve over and a uniform capacity $C$ for each tank. For each instance we create 30 scenarios, corresponding to different date ranges for the actual flow, and generate different historical sampling predictions for each site in each of the scenarios. Details of the problem instances are shown in Table 2. We consider scenarios of 24 and 48 hours length, the 48 hour instances allow water to be kept overnight in the tanks in order to smooth the outflows. We also briefly explored longer scenarios of 1 week but they did not lead to much greater peak outflow reductions than 48 hour scenarios.
**Table 2** Statistics of the 6 difference problem instances giving: the kinds water usages: unimodal (only using tanks that have a single peak in usage), bimodal (only using tanks that have bimodal water usage) and complete (using all types of tanks); number of tanks $|S|$; number of periods $n$; and the peak capacity of each tank $C$.

<table><thead><tr><th>Inst.</th><th>Types of usage</th><th>|S|</th><th>n</th><th>C</th></tr></thead><tbody><tr><td><b>I1</b></td><td>unimodal</td><td>1250</td><td>24</td><td>500</td></tr><tr><td><b>I2</b></td><td>unimodal</td><td>1250</td><td>48</td><td>500</td></tr><tr><td><b>I3</b></td><td>bimodal</td><td>1250</td><td>24</td><td>500</td></tr><tr><td><b>I4</b></td><td>bimodal</td><td>1250</td><td>48</td><td>500</td></tr><tr><td><b>I5</b></td><td>complete</td><td>2500</td><td>24</td><td>500</td></tr><tr><td><b>I6</b></td><td>complete</td><td>2500</td><td>48</td><td>500</td></tr></tbody></table>

## 4.2 Alternative approaches to controlling outflow

Because we only have a prediction of the future, the decisions made by the optimisation models may not be implementable with the actual inflows. Thus we implement out decisions as a “policy” for the tank to follow. If $X_{s,t} = 1$ then tank $s$ will empty only if there is sufficient volume in the tank ($\ge m$), and if $X_{s,t} = 0$ then the tank will still empty if the level would reach capacity $c_s$. This means that the decisions always lead to operation of the tank within specification. We denote this approach as HS (historical sampling).

The current operational approach and baseline is that a tank $s$ only empties when it reaches capacity. We can understand this as a set of decisions where $X_{s,t} = 0, \forall s \in S, t \in T$, since emptying only happens when capacity is reached. We denote this policy $[0, 0]$.

An alternate baseline strategy is to set $X_{s,t} = 1, \forall s \in S, t \in T$, this guarantees to empty each tank as soon as it has more than the minimal capacity. While unattractive in practice, since it induces significant wear and tear on the pumps, this may reduce peak outflows. We denote this policy $[1, 1]$.

We also consider a random policy by counting the average proportion of pump out periods $prop_s$ for each site $s$ using the baseline $[0, 0]$ policy. We draw a random number in $[0..1)$ for each site, and try to pump out if it is below $prop_s$. The random policy is defined as $X_{s,t} = 1 \rightarrow RANDOM([0..1]) \le prop_s, \forall s \in S, t \in T$.

Only using a single historical scenario is not robust, although since we are sampling for many tanks, certainly some of the variance of the problem is considered. We can make more robust plans by sampling for each tank $s$ instead $k$ scenarios, and computing the plan that leads to the minimum average maximum outflow across the scenarios (the deterministic equivalent). But the models are already slow, and this would substantially increases solve time. Instead we construct an artificial worst case scenario $w$, by in each time period defining the $i_{s,t}^w = \max\{i_{s,t}^i | i \in 1..k\}$, that is, the maximum inflow over all the $k$ scenarios. This has the same size, and hence solving difficulty, as a single scenario. We denote this approach as WHS (worst case historical sampling) where we use $k = 7$.

## 4.3 Results

The direct formulation while asymptotically smaller, $O(|S||T|)$), is challenging to solve for the size of problem we tackle. On instances with 2500 sites it struggles to find solutions quickly. The packing formulation which is ostensibly $O(|S||T|^2)$ but since $latest(s,t) - t$ is either small, or there are many time periods where $i_{s,t} = 0$ which do not require any pumped out decision variables, the actual size grows as $O(|S||T|)$), and the constraints are much simpler and hence faster to solve. For the remainder of the paper, we only report results for the packing formulation.
42:8 Predict-Then-Optimise Strategies for Water Flow Control

<table><caption>**Table 3** (Maximum mean outflow/standard deviation of outflow) for each problem instance and method across 30 scenarios.</caption><thead><tr><th></th><th>[0,0]</th><th>[1,1]</th><th>random</th><th>clairvoyant</th><th>HS</th><th>WHS</th></tr></thead><tbody><tr><th>I1</th><td>(13,793/223)</td><td>(14,444/170)</td><td>(10,858/204)</td><td>(8,596/169)</td><td><strong>(8,998/170)</strong></td><td>(9,514/173)</td></tr><tr><th>I2</th><td>(13,793/223)</td><td>(14,444/170)</td><td>(11,408/209)</td><td>(8,110/206)</td><td>(8,702/205)</td><td><strong>(8,592/202)</strong></td></tr><tr><th>I3</th><td>(12,176/254)</td><td>(12,425/217)</td><td>(11,156/236)</td><td>(10,859/212)</td><td>(11,387/221)</td><td><strong>(10,877/213)</strong></td></tr><tr><th>I4</th><td>(12,176/254)</td><td>(12,425/217)</td><td>(11,042/239)</td><td>(9,717/272)</td><td>(10,382/241)</td><td>(10,865/237)</td></tr><tr><th>I5</th><td>(24,438/303)</td><td>(25,165/273)</td><td>(19,578/284)</td><td>(18,446/261)</td><td>(18,793/276)</td><td>(20,309/270)</td></tr><tr><th>I6</th><td>(24,438/303)</td><td>(25,165/273)</td><td>(19,225/285)</td><td>(17,035/281)</td><td>(17,314/290)</td><td>(17,368/287)</td></tr></tbody></table>

Our proposed methods are compared against the current operational strategy [0,0] and an alternate base line [1,1]. In order to see how close to optimal we get we also compare against the *clairvoyant* approach which is running the optimisation model with the actual inflow data for the tested time period. This computes the minimal maximum outflow possible.

Figure 3 shows the behavior of the models running on instance I5. The clairvoyant approach illustrates that there is a significant reduction in peak outflow available if we make wise emptying decisions. The historical sampling approach actually gets quite close to the best solution on average but it is clear that the variance is large, often well over the clairvoyant solution. The worst case approach pays some penalty, it is unable to reduce the peak flow as well, but still is not that far from the clairvoyant solution, and its standard deviation is much smaller.

<figure><img src="image_427.png" /><figcaption>**Figure 3** Mean total outflows in each hour of the day for different approaches applied to instance I5. The shaded regions show the 25% - 75% confidence interval, across the 30 scenarios.</figcaption></figure>

Table 3 gives the summary results across the 6 instances. First note that the current baseline [0,0] is much better at reducing mean peak outflow than the alternative [1,1] of always pumping out when possible, but the standard deviation of the second method is much lower. The clairvoyant method shows that there is considerable reduction in peak available compared to the current baseline. The historical sampling optimisation approach HS is able to capture much of the available reduction in peak outflow and while the standard
deviation is larger than the clairvoyant approach it is not that much larger and considerably better than the current baseline. The worse case WHS approach also beats the baseline, and for some instances can be better than HS, its main strength is that usually reduces the standard deviation compared to HS. The random policy performs well for 24 period instances comparatively to 48 period instances. The HS and WHS approaches consistently perform better at reducing the peak flow when considering larger periods. The random policy is able to use the morning buffer capability of the tanks that is overlooked by the other approach for 24 period instances. For larger periods instances, the HS and WHS approaches systematically beat the random policy. This suggests a potential performance gain for CP based approach for larger period instances.

# 5 Conclusion

In this work, we introduced a novel practical problem from the water industry, controlling a pressure sewer system to reduce peak outflow. We proposed a direct MIP formulation and a packing formulation to solve practical scenarios. We provide an extensive experimental study on challenging and realistic instances of considerable size. The results show an optimisation model can significantly reduce the peak outflow of the system compared to the current operational approach.

So far we have only considered the most simple prediction approach, we plan to investigate forecasting models such as LSTM and LightGBM, to see whether they can produce realistic future flows. Ideally we would also extend the forecast to take into account parameters that affect the likely inflows, such as day of the week, season, and weather (the PSS we study is in a holiday zone, so inflow patterns change significantly on weekends, during summer, and when the weather is hot). It would be interesting to investigate Predict+Optimise approaches [6, 7] applied to the problem, but seeing that the predictions are for individual tanks and the objective results from considering all tanks simultaneously this appears challenging. While the simple optimisation approach we use here works we also plan to investigate decomposition approaches such as Benders or column generation, since the tanks are only weakly coupled by the objective. As future work, we plan to take into account more of the real features of the network such as the distance of tanks from the sewer treatment plant, the full topology of the network, and the inclusion of intermediate pump stations.

# References

1. Russell Bent and Pascal Van Hentenryck. Regrets only! Online stochastic optimization under time constraints. In *Proceedings of the 19th National Conference on Artificial Intelligence*, AAAI'04, pages 501-506. AAAI Press, 2004.

2. Russell Bent and Pascal Van Hentenryck. Scenario-based planning for partially dynamic vehicle routing with stochastic customers. *Operations Research*, 52:977-987, 2004. doi: 10.1287/opre.1040.0124.

3. Russell Bent and Pascal Van Hentenryck. The value of consensus in online stochastic scheduling. In *Proceedings of the 14th International Conference on Automated Planning and Scheduling*, ICAPS 2004, pages 219-226, 2004.

4. Russell Bent and Pascal Van Hentenryck. Online stochastic optimization without distributions. In *ICAPS 2005 - Proceedings of the 15th International Conference on Automated Planning and Scheduling*, pages 171-180, 2005.

5. Hyeong Soo Chang, Robert Givan, and Edwin K. P. Chong. On-line scheduling via sampling. In *Proceedings of the Fifth International Conference on Artificial Intelligence Planning Systems*, AIPS'00, pages 62-71. AAAI Press, 2000.
42:10 Predict-Then-Optimise Strategies for Water Flow Control

6 Adam Elmachtoub and Paul Grigas. Smart "predict, then optimize". Management Science, 68(1):9-26. 2022.

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

8 Nicholas Nethercote, Peter J. Stuckey, Ralph Becket, Sebastian Brand, Gregory J. Duck, and Guido Tack. Minizinc: Towards a standard CP modelling language. In Christian Bessière, editor, *Principles and Practice of Constraint Programming – CP 2007*, pages 529–543, Berlin, Heidelberg, 2007. Springer Berlin Heidelberg. doi:10.1007/978-3-540-74970-7_38.

A Clusters

<figure><img src="image_428.png" /><figcaption>Figure 4 Inflows (averaged) for the two types of flow: unimodal and bimodal. Cluster 1 shows a diurnal water usage (morning and evening peak). Cluster 2 shows a unique morning peak. Each consists of 2 sub clusters where the peaks are time shifted.</figcaption></figure>
# Constraint Programming Models for Depth-Optimal Qubit Assignment and SWAP-Based Routing

Kyle E. C. Booth

Amazon Quantum Solutions Lab, Seattle, WA, USA

## Abstract

Due to the limited connectivity of gate model quantum devices, logical quantum circuits must be compiled to target hardware before they can be executed. Often, this process involves the insertion of SWAP gates into the logical circuit, usually increasing the depth of the circuit, achieved by solving a so-called qubit assignment and routing problem. Recently, a number of integer linear programming (ILP) models have been proposed for solving the qubit assignment and routing problem to proven optimality. These models encode the objective function and constraints of the problem, and leverage the use of automated solver technology to find hardware-compliant quantum circuits. In this work, we propose constraint programming (CP) models for this problem and compare their performance against ILP for circuit depth minimization for both linear and two-dimensional grid lattice device topologies on a set of randomly generated instances. Our empirical analysis indicates that the proposed CP approaches outperform the ILP models both in terms of solution quality and runtime.

**2012 ACM Subject Classification** Computing methodologies → Planning and scheduling; Theory of computation → Constraint and logic programming

**Keywords and phrases** Qubit routing, quantum computing, constraint programming, combinatorial optimization

Digital Object Identifier 10.4230/LIPIcs.CP.2023.43

Category Short Paper

## 1 Introduction

The quantum circuit model of computation specifies quantum algorithms as sequences of logical quantum gates [12]. Quantum circuit compilation is the process of compiling a logical quantum circuit to a target quantum device such that the compiled circuit adheres to device-specific connectivity constraints. Commonly, this process involves the insertion of SWAP gates into the original circuit enabling qubits to move to neighboring locations on the hardware. Determining the initial qubit allocation as well as when and where these gates should be inserted has been studied under the name *qubit routing* [6].

Figure 1 provides an illustration of an input logical quantum circuit, a quantum device topology (represented as an undirected graph that specifies device connectivity), and a valid compiled circuit that inserted a single SWAP gate. Due to the effects of quantum decoherence, particularly in superconducting devices, it is often beneficial to solve the qubit assignment and routing problem while minimizing circuit depth (where the depth is the number of layers of parallelized gates). For this reason, additional gates should be used sparingly and with high parallelization, when possible.

A variety of techniques have been proposed for solving this problem in the literature, including both exact and heuristic approaches. Exact methods typically involve the use of search-based solvers leveraging smart inference techniques that, given enough time, will find and prove the optimal solution [1, 3, 8, 13]. Alternatively, heuristic methods (which have,

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

© Kyle E. C. Booth;

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. 43; pp. 43:1–43:10

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

LIPICS

Leibniz International Proceedings in Informatics

Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl Publishing, Germany
<figure><img src="image_431.png" /><figcaption>■ **Figure 1** Qubit assignment and routing problem example specifying the input logical circuit (a), target device topology (b) and a compiled circuit (c). The input circuit has a depth of 2 (gates $(q_1, q_3)$ and $(q_2, q_4)$ can be executed in parallel), while the compiled circuit has a depth of 3.</figcaption></figure>

until recently, been the focus of previous work) sacrifice completeness in favor of rapidly producing high quality circuits [8, 15, 6]; these methods tend to scale more effectively to larger problem instances as well.

In this short paper, we focus on the development of exact model-based approaches that outperform those from the literature. Specifically, we introduce and report initial results for two new CP models for the qubit assignment and routing problem. Our models have a relatively simple implementation, and leverage constraints supported by a variety of CP solvers. The first model is designed for linear array quantum device topologies, while the second can be used to solve problems involving more general architectures. We conduct an empirical analysis of the ILP models proposed by Boccia et al. [3] and Nannicini et al. [10] using circuit depth minimization as our objective function. We demonstrate, through empirical evaluation using different solvers, that our CP models outperform these existing ILP approaches in terms of solution quality and runtime for depth minimization.

The remainder of this paper is organized as follows. In Section 2 we define the specific variant of the qubit assignment and routing problem that we consider in this paper. In Section 3 we present our new CP models for the studied problem. In Section 4 we conduct an empirical assessment of the presented models on both linear array and 2D lattice architectures of varying size. Finally, in Section 5 we provide concluding remarks.

## 2 Preliminaries and definitions

We follow previously used notation with minor alterations to ease the presentation of the problem definition [10]. The input to the qubit assignment and routing problem consists of: i) a sequence of quantum gate groups (e.g., Figure 1a), and a hardware graph (e.g., Figure 1b). The hardware graph, $H = (V, E)$, consists of nodes $i \in V$, where each node represents a physical qubit on the quantum computer. The graph has edges $e \in E$, where each edge $e = \{i, j\}$ dictates pairs of physical qubits that can execute two-qubit gates (i.e., qubits that are neighboring each other on the architecture). Often, as in previous work [10], it is useful to define a directed graph $A$ such that each undirected edge in $E$ corresponds to two directed edges in $A$, $(i, j)$ and $(j, i)$. In this paper we study both linear and general hardware graphs (e.g., lattices).

The sequence of quantum gate groups is represented by $G = (G^1, G^2, \dots, G^{|L|})$ where $L$ is the set of layers in the original circuit and $G^\ell$ indicates the set of gates that must be executed in the $\ell^{th}$ layer. In this paper, we consider two-qubit quantum gates, $g = \{p, q\}$, involving a pair of logical qubits $p, q \in Q$ such that $G^\ell = \{g_1^\ell, g_2^\ell, \dots\}$. All of the two-qubit
gates executed in a given layer involve a set of non-overlapping pairs of qubits. This way, the gates in a layer are parallelized. To account for inserted SWAP gates, we augment the original circuit with a number of “auxiliary layers” between each of the logical circuit layers. We denote the total set of layers as $L'$, and use $L_{SWAP}$ to specify the auxiliary layers only.

Qubit assignment involves producing an initial mapping of logical qubits, $Q$, to physical qubits, $V$, on the quantum hardware. For simplicity, we refer to logical qubits and physical qubits as qubits and nodes from now on, respectively. We also assume that $|Q| = |V|$.<fnref n="1" />

Qubit routing is the process of moving qubits around the hardware architecture such that the gates at each layer can be executed in a way that satisfies the connectivity of the device. In this work, we consider routing accomplished via the use of SWAP gates (although we note there are other routing techniques [2]), two-qubit gates which exchange the position of two neighboring qubits. These SWAP gates are inserted into the auxiliary layers between the layers of the idealized quantum circuit and (often) increase its overall depth.

The objective of the qubit assignment and routing problem studied in this paper is to: i) assign qubits to nodes on the hardware, and ii) route the positions of the qubits such that the gates at each layer can be executed while satisfying the connectivity of the device. We seek to minimize the number of auxiliary layers utilized, effectively resulting in circuits with lower depth. This problem is known to be NP-Complete [5].

Following previous work, our assumptions on the formulation of the problem include:

- Qubits involved in a two-qubit gate can also swap positions. This is also called “merging” SWAPs with adjacent two qubit gates [7].

- Two-qubit gates are undirectional (i.e., the cost of executing a two-qubit gate in the forward configuration is the same as in the reverse configuration).

- We restrict our attention to circuits involving two-qubit gates, since single qubit gates can be merged with two-qubit gates and thus do not need to be considered.

# 3 Constraint programming models

For each of our proposed models, the main decision variable is $x_{ql} \in \{1, ..., |V|\}$ which represents the node location of logical qubit $q \in Q$ at layer $l \in L'$. Since ILP is less flexible than CP (e.g., due to linearity restrictions) the ILP models from the literature [10, 3, 9] must introduce additional variables to properly model the problem. In our models, we exploit the structure of the problem such that it is sufficient to constrain the values of $x_{ql}$ from one layer to the next and avoid additional variable overhead.

## 3.1 Linear array architectures

The first CP model we present is applicable to linear array architectures (such as the one presented in Figure 1b). The linear array is represented as a graph with nodes $V = \{1, 2, ..., |V|\}$ and edges $E = \{(i, i+1) : i \in V \setminus |V|\}$. This model is motivated by previous work in ILP [9], and employs absolute value constraints to ensure that qubit movement between layers is valid, and to ensure that qubits involved in a gate are at adjacent locations on the device.<fnref n="2" />

<footnote n="1">¹ Note that, in the case where $|Q| < |V|$, we can just introduce and route auxiliary qubits.</footnote>

<footnote n="2">² We note that the work in [9] does not minimize circuit depth, but rather the number of SWAPs inserted, and so we do not compare to it in this paper (non-trivial changes to the model must be made in order to express a depth minimization objective).</footnote>
43:4 CP Models for Depth-Optimal Qubit Assignment and Routing

The first constraint in our model ensures that, at each layer, each logical qubit is located at one node on the architecture, and all of the locations of the logical qubits are different. We use the alldifferent global constraint [13] to accomplish this:

$$ \text{alldifferent}(\{x_{1\ell}, \dots, x_{|Q|\ell}\}), \forall \ell \in L' \qquad (1) $$

The next set of constraints ensure that the gates specified in each layer are executed while adhering to the connectivity constraints of the target hardware. These constraints are expressed as follows:

$$ |x_{p\ell} - x_{q\ell}| = 1, \forall \ell \in L, g = (p, q) \in G^{\ell} \qquad (2) $$

Intuitively, these constraints specify that the locations of logical qubits $(p, q)$ involved in two-qubit gate $g$ must be neighboring. Note that, since this model is for a linear architecture (as defined above), this constraint is enough to ensure the gates in a layer can be executed; more complex architectures require more sophisticated modeling (as described in the next section).

A similar set of constraints is used to constrain the movement of qubits from one layer to the next. These are expressed as follows:

$$ |x_{q\ell} - x_{q\ell+1}| \le 1, \forall \ell \in \{1, \dots, |L'| - 1\}, q \in Q \qquad (3) $$

The above constraints dictate that, from one layer to the next, a given qubit cannot move more than one location away from its current location (preventing, for example, two simultaneous SWAP gates involving the same qubit). Additional constraints must be added to ensure that pairs of qubits, where qubit is involved in a gate at a given layer and the other is not (recall that the qubit pairs involved in a gate are permitted to SWAP), cannot SWAP from that layer to the next. This constraint is expressed as follows:

$$ |x_{p\ell} - x_{q\ell+1}| \le 1, \forall \ell \in L, g = (p, q) \in G^{\ell} \qquad (4) $$

Constraint (4) ensures that the qubit pair involved in a gate are still neighboring at the next layer. Note that this permits the qubit pair to exchange positions, but it does not allow either qubit to swap with another qubit not involved in the gate.

Finally, we encode the objective function of the optimization which is to minimize the number of auxiliary layers added to the circuit to support SWAP insertions; this is effectively the same as circuit depth, and follows previous work [10]. To express this objective, we introduce a binary decision variable $z_{\ell} \in \{0, 1\}$ for each layer, $\ell \in L_{SWAP}$. We constrain the variable such that it takes on a value of 1 whenever a SWAP gate is inserted into an auxiliary layer:

$$ z_{\ell} \ge |x_{q\ell} - x_{q\ell+1}|, \forall q \in Q, \ell \in L_{SWAP}. \qquad (5) $$

Constraint (5) tracks each time a qubit changes locations from one layer to the next, and the location change is not due to a (non-SWAP) gate operation. Finally, the objective function is expressed as:

$$ \min \sum_{\ell \in L_{SWAP}} z_{\ell} \qquad (6) $$

Note that if all qubit routing operations can be achieved by merged SWAPs, the optimal solution will be zero as no auxiliary SWAP gates need to be introduced.
The linear architecture CP model has $O(|Q||L'|)$ variables, each with a domain of $O(|V|)$, and $O(|Q||L'|)$ constraints. While the model is designed according to the assumptions stated in Section 2, it can be altered fairly easily to accommodate problem variations. For example, if qubits involved in a two-qubit gate cannot also SWAP places, we could simply constrain the locations of the qubits involved in the gate as necessary; though this would undoubtedly increase the depth of the produced circuits.

## 3.2 General architectures

When the problem is extended to more general architectures, the linear architecture model is no longer valid. This is because in the linear model, there is a clear mapping between the locations of the logical qubits and the chip connectivity; any qubits whose positions are within one of each other could have a gate applied between them. In the more general case, however, this is no longer true.<fnref n="3" />

Our model for general architectures is similar to that for linear array architectures, however, instead of using absolute value constraints we use **table** constraints. The **table** constraint specifies the list of tuples (solutions) to which a vector of variables can be fixed. For example, the constraint **table((y₁, y₂), {(1, 2), (2, 3)})** specifies that, in a solution, the variables $y₁$ and $y₂$ can be assigned values (1, 2) or (2, 3), respectively.

For our model, we utilize the **table** constraint to conveniently encode the connectivity of the hardware device, using $A$ as the list of tuples representing the locations of neighboring pairs of qubits. Our model for general architectures is given as follows:

$$ \min \sum_{\ell \in L_{SWAP}} z_{\ell} \qquad (7) $$

$$ \text{alldifferent}(\{x_{1\ell}, \dots, x_{|Q|\ell}\}) \qquad \forall \ell \in L' \qquad (8) $$

$$ \text{table}((x_{p\ell}, x_{q\ell}), A) \qquad \forall \ell \in L, g = (p, q) \in G^{\ell} \qquad (9) $$

$$ \text{table}((x_{q\ell}, x_{q\ell+1}), A \cup \{(i, i) : i \in V\}) \qquad \forall q \in Q, \forall \ell \in \{1, \dots, |L'| - 1\} \qquad (10) $$

$$ x_{q\ell} = x_{q\ell+1} \lor x_{q\ell} = x_{p\ell+1} \qquad \forall \ell \in L, g = (p, q) \in G^{\ell} \qquad (11) $$

$$ z_{\ell}|Q| \geq |x_{q\ell} - x_{q\ell+1}| \qquad \forall \ell \in L_{SWAP}, q \in Q \qquad (12) $$

The first constraint in the general architecture model, Constraint (8), is the same as Constraint (1) in the linear array model, leveraging the **alldifferent** global constraint.

Constraint (9) uses the **table** constraint to enforce that, at each gate layer, the qubits involved in each gate are at neighboring locations on the architecture (i.e., at locations represented by one of the tuples in arc set $A$).

Constraint (10) uses the **table** constraint to dictate the permissible movement of qubits from one layer to the next. Specifically, this constraint requires a qubit to stay in the same location (represented by the value $(i, i)$), or move to a neighboring location on the hardware. Constraint (11) adds some additional restriction regarding the movement of qubits involved in a two-qubit gate, similar to Constraint (4) in the linear array model. Finally, Constraint (12) links the objective function to the main variables, acting as a flag each time a qubit changes positions from one layer to the next (not including original circuit layers).

As with the linear array model, the general architecture model has $O(|Q||L'|)$ variables and constraints.

<footnote n="3">³ In a simple square topology with four nodes, $v₁$ and $v₄$ can be neighboring, but the absolute value of the difference of their node labels is not equal to one.</footnote>
4 Experiments

We present an empirical analysis of our models on a benchmark set of randomly generated
instances. Following previous work [10] we generate “square” circuits (similar to those used
for quantum volume testing) where |Q| = |L| (e.g., 4 × 4 indicates a circuit with four qubits
and four layers). For a given problem size, we generate a random permutation of the qubits
{1, 2..., |Q|} for each ℓ ∈ L. Based on this permutation, we group neighboring pairs of
qubits into the ⌊|Q|/2⌋ two-qubit gates for that layer (e.g., permutation 3-1-2-4 would yield
gates (3, 1) and (2, 4)). By increasing the size of Q and L, we test the scalability of our
models (e.g., a 10 × 10 instance will involve an initial circuit with 50 two-qubit gates). For
each problem size, we generate ten problem instances.

The target devices used for experimentation include a linear array and two-dimensional
lattice grids of increasing size. We generate the devices such that the number of hardware
vertices, |V|, is equal to the number of logical qubits, |Q|.<fnref n="4" /> The linear array device topologies
are straightforward. The lattice grids are generated for |Q| ∈ {4, 6, 8, 9, 10}, where Q ∈ {4, 9}
corresponds to square lattices, and the remainder are rectangular lattices.

The CP models are implemented with the CP-SAT solver in OR-Tools (v9.3) [11] using the Python interface. We include the symmetry breaking constraints used in previous ILP models [10]. The absolute value expressions in the CP models were implemented using auxiliary variables and the **AddAbsEquality** constraint, while the **table** constraint was implemented using the **AddAllowedAssignments** constraint. Finally, in the general architecture CP model, the dijunction in Constraint (11) was implemented using auxiliary boolean variables and the **OnlyEnforceIf** enforcement literal. Following previous work [10], we set the number of auxiliary layers between each original layer to four.

To thoroughly assess the previously proposed ILP models, we implement them in SCIP [1], Xpress (v8.13.5), and using the CP-SAT solver in OR-Tools.<fnref n="5" /> For the SCIP experiments, we use the OR-Tools modeling interface and select SCIP as the backend solver. All experiments are run with default search and inference settings on a machine with a 2.6 GHz 6-Core Intel i7 processor and 16GB of RAM. Additionally, some modifications to the model presented by Boccia et al. [3] are made to enable fair comparison. Specifically, auxiliary variables and constraints were added to permit the merged SWAP functionality.

**4.1 Linear array architectures**

Our first set of experiments involve linear array architectures. For these experiments, both
of the CP models we propose can be used to solve the problem. The results are visualized
in Figure 2. Our proposed models are denoted “CP-linear” (for the linear array model)
and “CP-general” (for the general model). We split the instances into two classes: the first
involves randomly generated square circuits of size 4 × 4, 5 × 5, 6 × 6, and 7 × 7 while the
second involves square circuits of sizes 8 × 8, 9 × 9, 10 × 10, and 11 × 11.

In Figure 2a, we illustrate the performance of all of the implementations on Class 1
instances with a solver time limit set to 10 seconds. The figure details the number of
instances solved to proven optimality vs. runtime. We can see that all of the methods
implemented with the CP-SAT solver in OR-Tools are able to quickly find and prove the
optimal solution to all 30 problems in less than two seconds. Conversely, the ILP models

<footnote n="4">⁴ Our approaches are not restricted to this case; auxiliary qubits can be used to handle the situation where |V| > |Q|.</footnote>

<footnote n="5">⁵ Since the ILP models do not contain continuous variables, this is straightforward.</footnote>
<figure><img src="image_432.png" /><figcaption>(a) Class 1 instances ($|Q| \in \{4, 5, 6, 7\}$).</figcaption></figure>

<figure><img src="image_433.png" /><figcaption>(b) Class 2 instances ($|Q| \in \{8, 9, 10, 11\}$).</figcaption></figure>

■ **Figure 2** Empirical results: CP models against ILP models from the literature for linear array device topologies. Number of instances solved to proven optimality over time. Time limit of 10 seconds (Class 1) and five minutes (Class 2).

implemented in SCIP and Xpress are much slower at finding and proving optimal solutions. For these approaches, the SCIP implementations seem to perform the worst, and Xpress provides a modest improvement. The experiments indicate that ILP (and ILP-based solvers) may not be the best candidate approach for this problem; this is in-line with previously reported results that showed ILP (implemented with an ILP solver) struggled to solve square circuit instances to proven optimality in short runtimes beyond six qubits [10]. As such, we elect to only investigate the models solved with OR-Tools CP-SAT for larger problem instances.

Figure 2b illustrates the performance of the CP-based implementations on medium-sized instances with a solver time limit of 5 minutes. From the figure, it is evident that our proposed CP models outperform the ILP models even when using the OR-Tools CP-SAT solver for all methods. Further, the linear array model (using absolute value constraints) performs slightly better than the more general model. Both models are able to find and prove optimality for all problem instances in less than 50 seconds, whereas the implementations of the ILP models require significantly more time.

Recall that, due to the insertion of SWAP gates as a result of our optimization, the number of layers in the final compiled circuits almost always increases. In Table 1 we detail the average circuit depth for each of the problem sizes, obtained by our depth-optimal CP methods. From the table, we can see that the inclusion of SWAP gates often doubles the depth of the circuit.

## 4.2 General architectures

Our second set of experiments involve 2D grid lattice topologies, and are visualized in Figure 3. We use the same solvers and settings as in the linear array experiments (recall that we cannot run the proposed linear array model for these lattice topology problems). For these tests, Class 1 instances are 4 × 4 and 6 × 6 (to permit grid lattice construction), while Class 2 instances are 8 × 8, 9 × 9, and 10 × 10.
<figure><img src="image_434.png" /><figcaption>(a) Class 1 instances (4 × 4, 6 × 6).</figcaption></figure>

<figure><img src="image_435.png" /><figcaption>(b) Class 2 instances (8 × 8, 9 × 9, 10 × 10).</figcaption></figure>

Figure 3 Empirical results: CP models against ILP models from the literature for 2D grid lattice device topologies. Number of instances solved to proven optimality over time. Time limit of 10 seconds (Class 1) and five minutes (Class 2).

In Table 2 we summarize the average depths of the optimally compiled circuits for each of the problem sizes. Comparing to the linear array results in Table 2, it is immediately apparent that the increased connectivity offered by 2D lattice topologies results in significantly shorter circuits. In the 10 × 10 case, for example, only 5.4 layers were added (on average) to permit SWAP operations, versus the 17.0 layers added (on average) in the linear case.

In terms of model/solver runtime performance, from Figure 3a we see a similar trend to the linear topology case: the ILP methods implemented in SCIP and Xpress struggle to find and prove optimal solutions within the runtime limit, while the OR-Tools implementations rapidly solve these problems (in less than one second). Figure 3b illustrates the performance of the OR-Tools implementations for the larger class of instances. As visualized in the figure, our proposed CP approach is able to find provably optimal solutions to all instances significantly faster than the ILP models from the literature.

## 5 Conclusions

In this paper we propose CP models for depth-optimal qubit assignment and SWAP-based routing. Our first model is specific to linear array topologies, while our second model is applicable to more general architectures (e.g., grid lattices). We conduct a series of experiments on randomly generated circuits, and demonstrate that the CP-based approaches

Table 1 Empirical results: CP model solution circuit depths by class, linear array device topologies.

<table><thead><tr><th>Class 1</th><th>Compiled circuit depth (avg.)</th><th>Class 2</th><th>Compiled circuit depth (avg.)</th></tr></thead><tbody><tr><td>4 × 4</td><td>6.1</td><td>8 × 8</td><td>18.0</td></tr><tr><td>5 × 5</td><td>6.2</td><td>9 × 9</td><td>20.9</td></tr><tr><td>6 × 6</td><td>10.9</td><td>10 × 10</td><td>27.0</td></tr><tr><td>7 × 7</td><td>12.4</td><td>11 × 11</td><td>30.7</td></tr></tbody></table>
<table><caption>Table 2 Empirical results: CP model solution circuit depths by class, 2D grid lattice device topologies.</caption><thead><tr><th>Class 1</th><th>Compiled circuit depth (avg.)</th><th>Class 2</th><th>Compiled circuit depth (avg.)</th></tr></thead><tbody><tr><td>4 × 4</td><td>4.0</td><td>8 × 8</td><td>10.3</td></tr><tr><td>6 × 6</td><td>6.5</td><td>9 × 9</td><td>12.2</td></tr><tr><td>—</td><td>—</td><td>10 × 10</td><td>15.4</td></tr></tbody></table>

provide superior performance over their ILP counterparts. Our results suggest that CP is a promising technology for producing provably depth-optimal circuits when qubit routing is accomplished via SWAP gate insertion.

References

1. Tobias Achterberg. SCIP: solving constraint integer programs. *Mathematical Programming Computation*, 1(1):1–41, 2009.

2. Aniruddha Bapat, Andrew M Childs, Alexey V Gorshkov, Samuel King, Eddie Schoute, and Hrishee Shastri. Quantum routing with fast reversals. *Quantum*, 5:533, 2021.

3. Maurizio Boccia, Adriano Masone, Antonio Sforza, and Claudio Sterle. SWAP minimization in nearest neighbour quantum circuits: An ILP formulation. In *Advances in Optimization and Decision Science for Society, Services and Enterprises*, pages 255–265. Springer, 2019.

4. Kyle E. C. Booth, Minh Do, J Beck, Eleanor Rieffel, Davide Venturelli, and Jeremy Frank. Comparing and integrating constraint programming and temporal planning for quantum circuit compilation. In *Proceedings of the International Conference on Automated Planning and Scheduling*, volume 28, pages 366–374, 2018.

5. Adi Botea, Akihiro Kishimoto, and Radu Marinescu. On the complexity of quantum circuit compilation. In *Proceedings of the International Symposium on Combinatorial Search*, volume 9, pages 138–142, 2018.

6. Alexander Cowtan, Silas Dilkes, Ross Duncan, Alexandre Krajenbrink, Will Simmons, and Seyon Sivarajah. On the qubit routing problem. *arXiv preprint arXiv:1902.08091*, 2019.

7. Petar Jurcevic, Ali Javadi-Abhari, Lev S Bishop, Isaac Lauer, Daniela F Bogorin, Markus Brink, Lauren Capelluto, Oktay Günlük, Toshinari Itoko, Naoki Kanazawa, et al. Demonstration of quantum volume 64 on a superconducting quantum computing system. *Quantum Science and Technology*, 6(2):025020, 2021.

8. Gushu Li, Yufei Ding, and Yuan Xie. Tackling the qubit mapping problem for NISQ-era quantum devices. In *Proceedings of the Twenty-Fourth International Conference on Architectural Support for Programming Languages and Operating Systems*, pages 1001–1014, 2019.

9. Jesse Mulderij, Karen I Aardal, Irina Chiscop, and Frank Phillipson. A polynomial size model with implicit SWAP gate counting for exact qubit reordering. *arXiv preprint arXiv:2009.08748*, 2020.

10. Giacomo Nannicini, Lev S Bishop, Oktay Günlük, and Petar Jurcevic. Optimal qubit assignment and routing via integer programming. *ACM Transactions on Quantum Computing*, 2021.

11. Laurent Perron. Operations research and constraint programming at Google. In *Principles and Practice of Constraint Programming-CP 2011: 17th International Conference, CP 2011, Perugia, Italy, September 12-16, 2011. Proceedings 17*, pages 2–2. Springer, 2011.

12. Eleanor G Rieffel and Wolfgang H Polak. *Quantum computing: A gentle introduction*. MIT Press, 2011.

13. Willem-Jan Van Hoeve. The AllDifferent constraint: A survey. *arXiv preprint cs/0105015*, 2001.
43:10 CP Models for Depth-Optimal Qubit Assignment and Routing

14 Davide Venturelli, Minh Do, Eleanor Rieffel, and Jeremy Frank. Compiling quantum circuits to realistic hardware architectures using temporal planners. *Quantum Science and Technology*, 3(2):025004, 2018.

15 Alwin Zulehner, Alexandru Paler, and Robert Wille. An efficient methodology for mapping quantum circuits to the IBM QX architectures. *IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems*, 38(7):1226–1236, 2018.
# Constraint Model for the Satellite Image Mosaic Selection Problem

Manuel Combarro Simón

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

Pierre Talbot

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

Grégoire Danoy

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

Jedrzej Musial

Poznan University of Technology, Poland

Mohammed Alswaitti

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

Pascal Bouvry

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

## Abstract

Satellite imagery solutions are widely used to study and monitor different regions of the Earth. However, a single satellite image can cover only a limited area. In cases where a larger area of interest is studied, several images must be stitched together to create a single larger image, called a mosaic, that can cover the area. Today, with the increasing number of satellite images available for commercial use, selecting the images to build the mosaic is challenging, especially when the user wants to optimize one or more parameters, such as the total cost and the cloud coverage percentage in the mosaic. More precisely, for this problem the input is an area of interest, several satellite images intersecting the area, a list of requirements relative to the image and the mosaic, such as cloud coverage percentage, image resolution, and a list of objectives to optimize. We contribute to the constraint and mixed integer linear programming formulation of this new problem, which we call the *satellite image mosaic selection problem*, which is a multi-objective extension of the polygon cover problem. We propose a dataset of realistic and challenging instances, where the images were captured by the satellite constellations SPOT, Pléiades and Pléiades Neo. We evaluate and compare the two proposed models and show their efficiency for large instances, up to 200 images.

2012 ACM Subject Classification Computing methodologies → Discrete space search

Keywords and phrases constraint modeling, satellite imaging, set covering, polygon covering

Digital Object Identifier 10.4230/LIPIcs.CP.2023.44

Category Short Paper

Supplementary Material Software (Source Code): https://github.com/mancs20/mosaic_imagecombination/tree/cp2023, archived at swh:1:dir:c776d54cd26c5685fd8a0a8ddfc2131771884b32

Funding Manuel Combarro Simón: This work is partially funded by the Luxembourg National Research Fund (FNR) – ASTRAL Project, ref. 17043604, and by the joint research programme UL/SnT-ILNAS on Technical Standardization for Trustworthy ICT, Aerospace, and Construction.

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

© Manuel Combarro Simón, Pierre Talbot, Grégoire Danoy, Jedrzej Musial, Mohammed Alswaitti, and Pascal Bouvry;
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. 44; pp. 44:1–44:15

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

LIPICS

Leibniz International Proceedings in Informatics
Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl Publishing, Germany
<figure><img src="image_438.png" /><figcaption>Figure 1 Optimization of the subset of images to make a mosaic of the Tokyo Bay region.</figcaption></figure>

Pierre Talbot: This work is supported by the FNR - COMOC Project, ref. C21/IS/16101289.

Jedrzej Musial: This work is funded by the FNR - PolLux program under the SERENITY Project, ref. C22/IS/17395419.

1 Introduction

The space industry is continuously growing and is no longer an exclusive market for military and government applications. According to the most recent report of the European Union Agency for the Space Programme (EUSPA) [13], the global market for navigation systems and earth observation (EO) had revenues of around €200 billion in 2022 and is expected to reach €500 billion by 2031. As access to space has become cheaper, an increasing number of private companies have entered the space business. Due to advances in satellite design and high-resolution remote sensors, the number of satellite launches dedicated to EO in 2021 was greater than the sum of launches between 2012 and 2016 [37]. In 2020 more than 100 terabytes of satellite imagery was generated per day [28].

There are several EO-based applications that analyze a vast area of interest (AOI) that can only be covered by combining several adjacent images into a larger one, called a mosaic. Mosaics are crucial for applications such as crop classification [21, 14], environmental monitoring [30, 15], and urban development analysis [38, 32]. The mosaicking of satellite images is a complicated process that presents challenges, such as color balancing [39] and image stitching [27].

In this work, we focus on the combinatorial problem of selecting the images to create the mosaic by optimizing one or several criteria. This problem is an extended version of the NP-hard problem of finding the minimum axis-parallel rectangle cover of a rectilinear polygon without holes [11], where the axis-parallel rectangles can be seen as the satellite images, and the rectilinear polygon as the AOI. In our problem, a cover is the subset of images that can be used to generate a mosaic. In Figure 1 a particular example of this problem is shown, where the objective is to build a mosaic using the smaller number of images. There are 30 images to choose from, and the optimization algorithm finds an optimal subset of four images.

In this paper, we present a multi-objective approach for this problem that seeks to optimize four popular parameters of satellite images for mosaic generation: cloud cover, incidence angle [1], resolution and cost of the images. In general, there might not be an objective that is more important than the other, which is why we propose a multi-objective approach, instead of a linear aggregation or lexicographic ordering of the objectives.

As the number of available satellite images has increased significantly, it is becoming more challenging to select the optimal combination of images to build a mosaic. The number of images covering one place can reach hundreds. This is even more difficult if the user is interested in optimizing several parameters. Without a computational approach for this,
users have to select by hand the images they want for the cover. Providing a Pareto front from which users can choose a cover is crucial to save money and time, considering that high-resolution satellite images are expensive.

In this paper, we propose a constraint programming (CP) model and a mixed integer linear programming (MILP) model to solve the problem presented. However, directly modeling this geometric problem with constraints is challenging, as we would need to encode geometric operations such as union and intersection of polygons. Instead, we preprocess each instance by computing a discretization where the intersections of all images are first computed (Section 2.1). We obtain a set of non-overlapping polygons where each polygon is simply an integer and the geometric characteristics can be ignored. This problem is a multi-objective extension of the well-known set covering problem (Section 3).

A unique aspect of this problem is to minimize the cloud coverage of the mosaic, which, in contrast to other objectives, does not have the same value throughout the image. While any part of the image has the same resolution, not all parts of the images have the same amount of clouds, except for images with 0 or 100% of cloud coverage. Because of this particularity of the problem, it is possible to reduce the cloud coverage percentage in the final mosaic by choosing a specific combination of images in such a way that cloudy regions of an image are overlapped by non-cloudy regions of other images. To the best of our knowledge, there is no work taking this into account to reduce the cloud coverage in the final mosaic. We call this problem the satellite image mosaic selection problem (SIMS) (Section 2).

The main contribution of this paper is to introduce the SIMS problem and present a CP model, as well as MILP model that can successfully find solutions to real instances of up to 200 images (Section 4). The constraint and mixed integer programming approaches are part of a larger framework where the images are automatically retrieved from different marketplaces and the solutions found by the solver can be visualized.

Although SIMS can be expressed as a linear problem (Section A), we choose to rely on constraint programming for two reasons: to ease the formulation of the model – in particular, it is convenient to use set variables – and because this problem aims to be extended for new requirements, which can be non-linear. The flexibility of the model is of utmost importance in this work, which is why constraint programming is our main choice.

## 2 Satellite Image Mosaic Selection Problem

The input for the SIMS problem is an area of interest (AOI) on Earth and a set of satellite images that intersect it. Each image has a cost and a list of parameters including the resolution, incidence angle and cloud coverage. The AOI is represented as a simple closed polygon without holes, and the images are represented as quadrilaterals. For both the AOI and satellite images, the corner coordinates are provided. With that information, the AOI and the images can be represented in the plane.

As clouds are not usually even distributed in the images, having images with a certain cloud coverage percentage does not guarantee that the final mosaic has less than that cloud coverage. Depending on the cloud distribution in the images, the final mosaic can have a lower or higher percentage of cloud coverage as depicted in Figure 2. This is not the case for the other objectives because they have a unique value along the image. For example, if all images have a determined resolution, the final mosaic will have the same resolution.
## 44:4 Satellite Image Mosaic Selection Problem

<figure><img src="image_439.png" /><figcaption>Figure 2 A cloudy region of an image can be covered by a non-cloudy region of another image, impacting the cloud coverage percentage of the final mosaic.</figcaption></figure>

<table><caption>Table 1 Number of intersections for the instances covering the Tokyo Bay region.</caption><thead><tr><td>Images</td><td>Intersections</td></tr></thead><tbody><tr><td>30</td><td>298</td></tr><tr><td>50</td><td>806</td></tr><tr><td>100</td><td>3278</td></tr><tr><td>150</td><td>8079</td></tr><tr><td>200</td><td>14855</td></tr></tbody></table>

### 2.1 Preprocessing of the Problem

To make the discretization, we first remove the parts of the images that are outside the AOI, and then we find all the polygons resulting from the intersection of the images, we find the polygons using the GEOS library [17]. The universe is partitioned into a set of polygon elements, and each of them is assigned to its corresponding images. In Figure 3, we show an example of this process, where we generate 254 polygons from 30 images. In Table 1, we show the number of intersections for different cardinality of the images set to cover the Tokyo Bay region.

The following step is to detect the clouds in the images and add them to the universe and to the correspondig sets. The objective of doing this is to know whether a region of the final mosaic, that is represented by one element of the universe, is free of clouds or not. We consider that a region of the AOI is free of clouds if there is at least one image in the cover in which that region does not have clouds.

<figure><img src="image_440.png" /><figcaption>Figure 3 298 polygons are obtained after preprocessing of 30 images. First the area of the images outside the AOI is removed and then the polygons resulting from the images intersections are found.</figcaption></figure>
<figure><img src="image_441.png" /><figcaption>(a) Before the cloud integration the universe had 3 elements and 7 after the integration.</figcaption></figure>

<figure><img src="image_442.png" /><figcaption>(b) Regions represented by element 5 and 6 are free of clouds in the final mosaic.</figcaption></figure>

■ **Figure 4** Possible cases for the integration of the detected clouds to the universe and the corresponding sets. Set I and II are represented by the rectangles with red and blue borders respectively.

A cloudy region is an element in the universe that is present in all the sets that cover that region. For the set that represents the image in which the cloud occurred, we make a distinction and we say that the element is cloudy for that set. In this way, a set is composed of cloudy elements and non-cloudy elements.

In real applications, clouds can be detected using cloud detectors [36, 19, 22]. As this is a problem orthogonal to our work, we do not detect the clouds, but instead randomly allocate them in the parts of the image. For each image, we have metadata indicating the cloud coverage percentage of the image. With that information and knowing the elements that belong to the image, we randomly set one of the elements as cloudy. We repeat this operation until the cloud coverage percentage of the image is achieved.

In Figure 4a, all possible scenarios of how clouds are converted to elements of the universe are shown for the general case. To facilitate the understanding of this process, only two overlapping images are shown, but the procedure is the same when more than two images overlap. Images I and II are represented as rectangles with red and blue borders, respectively. The clouds in image I are colored red, and the ones in image II are colored blue. Initially, both sets have in common its intersection, element 3 ($I = \{1, 3\}$ and $II = \{2, 3\}$). We can see that both clouds of I are partially covered by II; in one case is because one part of the cloud is in the intersection, and in the other case is because one part of the cloud is overlapped by a cloud from II, so it can not be completely covered by a non-cloudy region of II. From the previous, we can see that three elements are created, 4, 5 and 7. Element 4 are the clouds that are only present in I, element 5 is the cloudy region of I that is not cloudy in II, i.e. covered by II, and element 7 is a cloudy region of I that is also cloudy in II. Element 6 is similar to element 5, is a cloudy area in II that is not cloudy in I. When all the clouds are detected and incorporate to the universe, the original three elements 1, 2 and 3 are modified as follows: element 1 is the non-cloudy region of I that is not overlapped by any other image, element 2 is equivalent to element 1 but for image II, and element 3 is the non-cloudy area of the intersection between I and II. Finally, the universe has seven elements, and the sets are $I = \{1, 3, 4, 5, 6, 7\}$, $II = \{2, 3, 5, 6, 7\}$. In Figure 4b, we show a resulting mosaic after covering the clouds. Importantly, taking into account the clouds in this way increases the cardinality of the universe, but the problem itself does not change, which is why we took a simpler approach to randomly assign clouds to each element.
3 Constraint Model

Let $U = \{k_1, \dots, k_n\} \subset \mathbb{N}$ be a set of $n$ parts of the area of interest, called the universe. The set $U$ is a polygon partition of the area of interest, i.e. two parts do not overlap and their union is exactly the area of interest. Each satellite image is represented by a collection $P_i \subset U$ of parts. We write $I = \{P_1, \dots, P_m\}$ the set of all $m$ satellite images. The goal is to find a subset $T \subset \{1, \dots, m\}$ of images that covers the area of interest. The parameters of the model are $U$ and $I$ while $T$ is the main decision variable. The set covering constraint is captured by the following:

$$ \bigcup_{i \in T} P_i = U \qquad (1) $$

A trivial solution to this constraint is to take all the images, but we usually consider an optimization version where the cardinality of $T$ is minimized. In our case, each image $i \in \{1, \dots, m\}$ has a cost $W_i \in \mathbb{N}$ that we seek to minimize:

$$ \min \sum_{i \in T} W_i \qquad (2) $$

Along with Equation 1, this problem is called weighted set cover. Depending on the user requirements, we can consider other objectives such as resolution and incidence angle. For each part $k \in U$, we have its area $A_k \in \mathbb{N}$. And for each image $i \in \{1, \dots, m\}$, we have its resolution $R_i \in \mathbb{N}$ and its incidence angle $F_i \in \mathbb{N}$. We seek to minimize (the resolution is given in how many $cm^2$ represents a pixel, the less the better) the best resolution obtained for each part:

$$ \min \sum_{k \in U} \min\{R_i \mid i \in T, k \in P_i\} \qquad (3) $$

For the incidence angle, we seek to minimize the maximal angle, although other choices would be possible such as minimizing the average.

$$ \min \{\max \{F_i \mid i \in T\}\} \qquad (4) $$

A more challenging aspect of this problem is to minimize the area covered by clouds. To achieve that, we consider that each part is either cloudy or not. We leave the cloud detection and the splitting of the image into cloudy and non-cloudy parts to a preprocessing step. Let $C_i \subset P_i$ the cloudy parts of the image $i$. For each part $k \in U$, we define $D_k := \{i \in \{1, \dots, m\} \mid k \in P_i \setminus C_i\}$ the set of all images containing a non-cloudy view of the part $k$. For each part $k \in U$, we can now define the Boolean variable $V_k$ to be true when the part $k$ is cloudy in the cover:

$$ V_k \Leftrightarrow \bigwedge_{i \in D_k} i \notin T \qquad (5) $$

We can now minimize the area covered by clouds:

$$ \min \sum_{k \in U} V_k * A_k \qquad (6) $$

This objective can also be turned into a constraint if the user only wants covers with a certain cloud coverage threshold.
The model introduced is actually linear, as shown in Appendix A, and can be solved by mixed integer programming solvers. We simply represent the set $T$ by $m$ 0-1 variables $\{x_1, ..., x_m\}$ such that $x_i = 1$ if we take the image and $x_i = 0$ otherwise. We also use this representation for constraint programming solvers, because it is not possible to represent the set covering constraint otherwise – this is due to $T$ having a non-fixed cardinality.

## 3.1 Search Strategy Based On Greedy Algorithm

A well-known greedy algorithm for the set covering problem consists in taking the images covering the most uncovered parts of the universe first [9]. We model this heuristic as a search strategy within the MiniZinc constraint model. This has the advantage of always producing a solution that is at least as good as the greedy heuristics – since it is the first solution found. To achieve that, we reuse an existing search strategy provided by MiniZinc. A second advantage is that our search strategy can be reused with any constraint solver compatible with MiniZinc. We select the variable using the `anti_first_fail` strategy – the variable with the largest domain is selected first – and we take the highest value in its domain (`indomain_max`). The trick is to model a set of variables $\{G_1, ..., G_m\}$ such that $G_i \in \{0, ..., |P_i|\}$ is equal to the number of parts covered by the image $i$. Actually, in any solution, we have $G_i = |P_i|$ since the whole universe must be covered. What is interesting is the value of $G_i$ in partial assignments during the search. The difference between the upper and lower bounds $max(G_i) - min(G_i)$ is the number of parts that are currently uncovered by the partial assignment, and that can be covered by the image $i$. Since the anti-first-fail strategy selects the largest domain first, it effectively implements the greedy heuristics. We model $G_i$ as follows:

$$G_i = \sum_{k \in P_i} (\bigvee_{i \in T} k \in P_i) \qquad (7)$$

We note that the new variables $G_i$ are fully defined with the parts, and therefore once the main decision variable $T$ is assigned, the variables $G_i$ must be assigned as well.

## 3.2 Multi-Objective Constraint Optimization Algorithm

The multi-objective constraint programming algorithm used in this work was pionneered by Gavanelli [16] and has been frequently used in constraint optimization [23, 33, 18]. The main idea is to run a satisfaction constraint solver iteratively and add new constraints representing the Pareto front to ensure the next solution is not dominated by any point in the current Pareto front. To illustrate this algorithm, suppose a biobjective maximization problem where $x$ and $y$ are the two variables to optimize. We run the constraint solver which returns a first satisfiable solution where $x = 10$ and $y = 5$. At that point the Pareto front is $\{(10, 5)\}$. We add to the model the constraint $x > 10 \lor y > 5$ which guarantees that the next solution will not be dominated by the current points in the Pareto front. The solver might then find the solution $x = 2$ and $y = 6$ which is incomparable to the previous solution and is added to the Pareto front $\{(10, 5), (2, 6)\}$. The constraint generated from the Pareto front is now $(x > 10 \lor y > 5) \land (x > 2 \lor y > 6)$. This process continues until the solver finds an unsatisfiable solution, in which case we are guaranteed to have found the optimal Pareto front.
# 4 Evaluation

## 4.1 Dataset Description

For each experiment, there is an AOI and a number of images that cover the AOI. The objective is to find the Pareto front, where each point in the front represents a subset of images that must cover the AOI and optimize four objectives: cost, resolution, incidence angle and cloud coverage.

To carry out this research, we developed a framework capable of retrieving image metadata from different satellite marketplaces, preprocessing it (discretization and cloud integration to the universe), calling a CP or a MILP solver, and visualizing the solutions from the Pareto front.

We selected five AOIs from around the world: Mexico City (Mexico), Rio de Janeiro (Brazil), Paris (France), Lagos (Nigeria), and Tokyo Bay (Japan). For each AOI, we obtained all available images that were captured from 01-01-2021 to 01-01-2023 by the following satellite constellations SPOT [4], Pléiades [2] and Pléiades Neo [3]. We opted for those satellite constellations as they have all the metadata used in the experiments; other satellite constellations lacked some parameters such as cloud coverage or incidence angle.

Five instances were generated for all AOIs, except Lagos. Each of these instances differs from each other by the number of images given to cover the AOIs. The number of images for the instances were 30, 50, 100, 150 and 200. For Lagos, the total number of images available for the specified date range was 145, so the number of images for the Lagos instances were 30, 50, 100 and 145.

## 4.2 Experimental Setup

Each instance was solved using the CP model and the MILP model. We run the CP model with two solvers, OR-tools [31] and Gecode 6.3.0 [34]. For each of these solvers, we run the experiments twice; one with the default solver search strategy, and the other one with the greedy search strategy proposed in Section 3.1. The MILP model was implemented using the Gurobi solver [20]. We will refer to these five approaches as OR-tools default, OR-tools greedy, Gecode default, Gecode greedy and Gurobi.

For the MILP model, the algorithm used to obtain the exact Pareto front was SAUG-MENCON [40] which is based on the AUGMECON [25, 26] algorithm and on the well-known $\epsilon$-constraint method. In the $\epsilon$-constraint methods, one objective is optimized, and the others are added as constraints to the model. The right-hand side of the objective constraints gradually changes from the less restrictive values of the objectives to the most restricted ones. This process continues until all combinations of values for the constraint objectives have been explored. The SAUGMENCON method introduces two acceleration mechanisms to improve the computational efficiency of the front generation.

The experiments were run on an AMD Epyc ROME 7H12 processor (64 cores, 280W). All the solvers were configured to run in parallel with 8 cores and 16 threads. The running time for each experiment was 1 hour.

## 4.3 Experimental Results

To compare the results, we used the hypervolume of the Pareto front, which is a standard metric for comparing fronts in multiobjective optimization. For each instance, we score the strategies, calculating how worst they are compared to the best. For example, a score of 1 means that the strategy has the same hypervolume value as the best one, and a score of 0.5 means that the hypervolume of the strategy is half of the best hypervolume for that instance.
<figure><img src="image_443.png" /><figcaption>Figure 5 Number of times each approach had the best hypervolume.</figcaption></figure>

<figure><img src="image_444.png" /><figcaption>(a) Considering the 24 instances.</figcaption></figure>

<figure><img src="image_445.png" /><figcaption>(b) Considering the 21 instances where all strategies found a solution.</figcaption></figure>

<caption>Figure 6 Average score for each approach.</caption>

In Table 2 from Appendix C we can see the hypervolume values for each strategy for all instances.

For only one instance, the complete Pareto front was found in the running time, the strategies that found the complete front were *OR-tools default*, *OR-tools greedy* and *Gurobi*. For the rest of the instances, the whole Pareto front was not found, and the hypervolume corresponds to the partial front found during the running time. For 3 of the 4 instances with 200 images, the CP strategies could not find any point of the front; the entire running time was employed by the FlatZinc submodule of Minizinc, to flatten the model with the data file. However, Gurobi for 3 of these instances could find one point of the Pareto front. For the rest of the instances, all the strategies could find at least one point of the Pareto front, and generally the CP approaches obtained superior results compared to MILP.

As we can see in Figures 5 and 6a, for these experiments, the best approach was *OR-tools default*, being the best for 13 out of 24 instances and with a score average of 0.806. The second and third best strategies were *OR-tools greedy* and *Gurobi*, with a similar performance. The fourth and fifth places were occupied by Gecode default and greedy, being really close.

If we do not consider the 3 instances in which the CP strategies could not find a solution, the score averages change; see Figure 6b. *OR-tools default* has an average score very close to 1, and *OR-tools greedy* has a much better average score than *Gurobi*, which for these instances is the worst strategy in average. The difference in the average score between both Gecode strategies remains very close.

It is interesting to note that *Gurobi* showed excellent performance for small instances comprising 30 and 50 images. However, its performance was not as impressive for larger instances. This could be related to the way solutions are discovered and added to the front. For future research, it could be interesting to compare the hypervolume anytime behavior for different approaches used for CP and MILP to get the exact Pareto front.
## 5 Related Work

Geometric set covering problems can be divided into two categories based on the requirements of the covering shapes. In one category, the covering shapes do not have a fixed position in the plane, for example, covering a polygonal region with the minimum amount of fixed sized rectangles [24] or with a set of known rectangles that can freely move on the plane [35]. The other category is where the covering shapes have a fixed size and position, for example, covering a polygonal region with discs with a fixed size and position on the plane [10]. SIMS, belongs to the second category as the satellite images represent a fixed region on the Earth.

Considering the cloud coverage percentage in the final mosaic makes SIMS problem different from polygon cover and other geometrical cover problems, where the covering shapes only have to cover the polygon or the universe of points in the space. In this problem, the covering shapes, besides covering the polygon, should also cover certain regions (clouds) that are present in the shapes. Interestingly, this can be seen as solving two weighted set covering problems; in one, the AOI must be covered and in the other, the clouds.

The main approaches to solving geometric set covering problems are local search [6, 12, 5] and linear programming (LP) [7, 8]. In most of the papers, opposite to SIMS, the universe is a set of points instead of a region. In [10], they provide an exact algorithm for the case where the universe is a set of regions, and the covering objects are discs. The algorithm is effective when the minimum number of discs to cover the space is low.

In [29], set covering problem is tackled using constraint programming. There, the authors propose a way to prune the domain of possible solution using a lower and upper bound for the objective value. The lower bound consists of determining the minimum number of sets that can cover the space. This is equivalent to answering the following NP-complete problem: does a cover of the universe exist with *K* sets. They propose a new strategy to get an approximation of this lower bound and compare it against two other well-known lower-bound values: the value of the LP relaxation problem and a greedy algorithm. The proposed prune strategy is good for problems where the size of the sets is small, for bigger subsets, they recommend alternating between the LP relaxation and the greedy algorithm.

## 6 Conclusion

In this paper, we introduce a novel geometrical NP-hard problem, SIMS, inspired by the selection of satellite images for mosaic generation. CP and MILP models are provided for this problem, together with a search strategy for the CP model, based on the well-known greedy algorithm used for set covering problems. In the experiments performed, the CP solved with OR-tools got the best result, evidencing the power of this solver. Our proposed search strategy could not outperform the default search strategies, but in the case of the Gecode solver, it produced similar results. Generally, the CP model outperformed the MILP model. This could be related to the method used to generate the Pareto front. For future work, it will be interesting to compare different approaches to generate the exact Pareto front for the CP and MILP models, based on the metric anytime behavior for the hypervolume. We also plan to propose heuristics to tackle larger instances and to evaluate their performance against the proposed CP and MILP models.

## References

1. Airbus. Incidence angle. URL: https://www.intelligence-airbusds.com/en/8719-angle-conversion.
2 Airbus. Pléiades. URL: https://www.intelligence-airbusds.com/en/8692-pleiades.

3 Airbus. Pléiades Neo. URL: https://www.airbus.com/en/products-services/space/earth-observation/earth-observation-portfolio/pleiades-neo.

4 Airbus. SPOT. URL: https://www.intelligence-airbusds.com/en/8693-spot-67.

5 Pradeesha Ashok, Aniket Basu Roy, and Sathish Govindarajan. Local search strikes again: PTAS for variants of geometric covering and packing. *Journal of Combinatorial Optimization*, 39(2):618-635, June 2019. doi:10.1007/s10878-019-00432-y.

6 Nikhil Bansal and Kirk Pruhs. Weighted geometric set multi-cover via quasi-uniform sampling. In *Algorithms – ESA 2012*, pages 145–156. Springer Berlin Heidelberg, 2012. doi:10.1007/978-3-642-33090-2_14.

7 Timothy M. Chan and Qizheng He. Faster approximation algorithms for geometric set cover. In *Leibniz International Proceedings in Informatics*, LIPIcs, volume 164. Schloss Dagstuhl–Leibniz-Zentrum fur Informatik GmbH, Dagstuhl Publishing, June 2020. doi:10.4230/LIPIcs.SoCG.2020.27.

8 Chandra Chekuri, Sariel Har-Peled, and Kent Quanrud. Fast LP-based approximations for geometric packing and covering problems. In *Proceedings of the Fourteenth Annual ACM-SIAM Symposium on Discrete Algorithms*, pages 1019–1038. Society for Industrial and Applied Mathematics, January 2020. doi:10.1137/1.9781611975994.62.

9 V. Chvatal. A greedy heuristic for the set-covering problem. *Mathematics of Operations Research*, 4(3):233–235, August 1979. doi:10.1287/moor.4.3.233.

10 Claudio Contardo and Alain Hertz. An exact algorithm for a class of geometric set-cover problems. *Discrete Applied Mathematics*, 300:25–35, 2021. doi:10.1016/j.dam.2021.05.005.

11 Joseph C. Culberson and Robert A. Reckhow. Covering polygons is hard. *Annual Symposium on Foundations of Computer Science (Proceedings)*, pages 601–611, 1988. doi:10.1109/SFCS.1988.21976.

12 Minati De and Abhiruk Lahiri. Geometric Dominating-Set and Set-Cover via Local-Search. *Computational Geometry*, 113:102007, August 2023. doi:10.1016/j.comgeo.2023.102007.

13 European Union Agency for the Space Programme. EUSPA EO and GNSS Market Report. Technical Report 1, European Union Agency for the Space Programme, 2022. URL: https://www.euspa.europa.eu/sites/default/files/uploads/euspa_market_report_2022.pdf.

14 Shilan Felegari, Alireza Sharifi, Kamran Moravej, Muhammad Amin, Ahmad Golchin, Anselme Muzirafuti, Aqil Tariq, and Na Zhao. Integration of Sentinel 1 and Sentinel 2 Satellite Images for Crop Mapping. *Applied Sciences*, 11:10104, 2021. doi:10.3390/app112110104.

15 Neil Flood, Fiona Watson, and Lisa Collett. Using a U-net convolutional neural network to map woody vegetation extent from high resolution satellite imagery across Queensland, Australia. *International Journal of Applied Earth Observation and Geoinformation*, 82:101897, October 2019. doi:10.1016/J.JAG.2019.101897.

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

17 GEOS contributors. *GEOS coordinate transformation software library*. Open Source Geospatial Foundation, 2021. URL: https://libgeos.org/.

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

19 Yanan Guo, Xiaoqun Cao, Bainian Liu, and Mei Gao. Cloud detection for satellite imagery using attention-based u-net convolutional neural network. *Symmetry*, 12(6):1056, June 2020. doi:10.3390/sym12061056.

20 Gurobi Optimization, LLC. Gurobi Optimizer Reference Manual, 2023. URL: https://www.gurobi.com.
44:12 Satellite Image Mosaic Selection Problem

21 Ola Hall, Sigrun Dahlin, Håkan Marstorp, Maria Francisca Archila Bustos, Ingrid Öborn, and Magnus Jirström. Classification of maize in complex smallholder farming systems using UAV imagery. *Drones*, 2(3):1–8, 2018. doi:10.3390/drones2030022.

22 Jacob Høxbroe Jeppesen, Rune Hylsberg Jacobsen, Fadil Inceoglu, and Thomas Skjødeberg Toftegaard. A cloud detection algorithm for satellite imagery based on deep learning. *Remote Sensing of Environment*, 229:247–259, August 2019. doi:10.1016/j.rse.2019.03.039.

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

24 Sina Sharif Mansouri, George Georgoulas, Thomas Gustafsson, and George Nikolakopoulos. On the covering of a polygonal region with fixed size rectangles with an application towards aerial inspection. In *2017 25th Mediterranean Conference on Control and Automation (MED)*. IEEE, July 2017. doi:10.1109/med.2017.7984284.

25 George Mavrotas. Effective implementation of the ϵ-constraint method in multi-objective mathematical programming problems. *Applied Mathematics and Computation*, 213(2):455–465, July 2009. doi:10.1016/j.amc.2009.03.037.

26 George Mavrotas and Kostas Florios. An improved version of the augmented ϵ-constraint method (AUGMECON2) for finding the exact pareto set in multi-objective integer programming problems. *Applied Mathematics and Computation*, 219(18):9652–9669, May 2013. doi:10.1016/j.amc.2013.03.002.

27 V. Megha and K. K. Rajkumar. Automatic Satellite Image Stitching Based on Speeded Up Robust Feature. In *Proceedings - 2021 1st IEEE International Conference on Artificial Intelligence and Machine Vision, AIMV 2021*, pages 1–6, Gandhinagar, India, 2021. Institute of Electrical and Electronics Engineers Inc. doi:10.1109/AIMV53313.2021.9670954.

28 Doug Mohney. Terabytes From Space: Satellite Imaging is Filling Data Centers, 2020. URL: https://datacenterfrontier.com/terabytes-from-space-satellite-imaging-is-filling-data-centers/.

29 Sébastien Mouthuy, Yves Deville, and Grégoire Dooms. Global constraint for the set covering problem. *Journées Francophones de Programmation par Contraintes*, pages 183–192, 2007.

30 Konstantinos G. Nikolakopoulos, Paraskevi Lampropoulou, Elias Fakiris, Dimitris Sardelianos, and George Papatheodorou. Synergistic use of UAV and USV data and petrographic analyses for the investigation of beachrock formations: A case study from Syros Island, Aegean sea, Greece. *Minerals*, 8(11):534, November 2018. doi:10.3390/min8110534.

31 Laurent Perron and Vincent Furnon. Or-tools. URL: https://developers.google.com/optimization/.

32 Simone Piaggesi, Laetitia Gauvin, Michele Tizzoni, Natalia Adler, Stefaan Verhulst, Andrew Young, Rahmannan Price, Leo Ferres, Ciro Cattuto, and André Panisson. Predicting city poverty using satellite imagery. In *IEEE Computer Society Conference on Computer Vision and Pattern Recognition Workshops*, volume 2019-June, pages 90–96, 2019.

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

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

35 Y. G. Stoyan, T. Romanova, G. Scheithauer, and A. Krivulya. Covering a polygonal region by rectangles. *Computational Optimization and Applications 2009 48:3*, 48(3):675–695, May 2009. doi:10.1007/S10589-009-9258-1.

36 Lin Sun, Xu Yang, Shangfeng Jia, Chen Jia, Quan Wang, Xinyan Liu, Jing Wei, and Xueying Zhou. Satellite data cloud detection using deep learning supported by hyperspectral data. *International Journal of Remote Sensing*, 41(4):1349–1371, September 2019. doi:10.1080/01431161.2019.1667548.
37 Union of Concerned Scientists. UCS Satellite Database. URL: https://www.ucsusa.org/resources/satellite-database.

38 Haibo Wang, Xueshuang Gong, Bingbing Wang, Chao Deng, and Qiong Cao. Urban development analysis using built-up area maps based on multiple high-resolution satellite data. *International Journal of Applied Earth Observation and Geoinformation*, 103:102500, December 2021. doi:10.1016/J.JAG.2021.102500.

39 Lei Yu, Yongjun Zhang, Mingwei Sun, and Xinyu Zhu. Colour balancing of satellite imagery based on a colour reference library. *International Journal of Remote Sensing*, 37(24):5763–5785, December 2016. doi:10.1080/01431161.2016.1249306.

40 Weihua Zhang and Marc Reimann. A simple augmented $\epsilon$-constraint method for multi-objective mathematical integer programming problems. *European Journal of Operational Research*, 234(1):15–24, April 2014. doi:10.1016/j.ejor.2013.09.001.

A Linear Programming Model

For the mixed integer linear programming, we use the same nomenclature as for the constraint model. We just add the necessary variables to linearize the model.

To linearize the cover constraint (1) it is necessary to associate each image $P_i$ with a decision variable $x_i$ that is equal to 1 if the image $i$ is selected, otherwise it is 0. We rewrite the constraint as follows:

$$ \sum_{i:k \in P_i} x_i \geq 1, \text{ for all } k \in U \qquad (8) $$

The previous constraint guarantees that all the parts are covered by at least one image.

To linearize the cost constraint (2) it is necessary to associate each image $P_i$ with an auxiliary variable $w_i$ that represents the cost of the image. The linear constraint can be written as:

$$ \min \sum_{P_i \in I} x_i w_i \qquad (9) $$

The constraints (8) and (9) are the classical constraints used for set covering problems.

The resolution objective is a min-min problem, where the objective is to minimize the sum of the min resolution of each part. The min resolution of a part is the minimum resolution of the images that contain them and belong to a cover. We need to add an auxiliary decision variable $r_k$ representing the best resolution of the part $k$ and a big constant $B$, bigger than the maximum image resolution. Also, we need to add an auxiliary binary decision variables $z_{k_j}$ for each image $P_j$ that contains $k$. For each part $k$, we define $L_k := \{i \in \{1, \dots, m\} | k \in P_i\}$ as the set of all images containing the part $k$. For each part $k$ we can now define a constraint for the values that can take the variables $z_{k_j}$.

$$ \sum_{k=1}^{|L_k|} z_{k_j} = |L_k| - 1 \qquad (10) $$

The constraint expressed above states that only one of the $z_{k_j}$ variables can be 0, the rest have to be 1. We define the minimum resolution of a part as $r_k$. With the following two constraints, we can linearize (3).

$$ r_k \geq (x_j R_j + B(1 - x_i)) - 2Bz_{k_j} \quad \text{for all } j \in L_k \qquad (11) $$

$$ \min \sum_{k \in U} r_k \qquad (12) $$
## 44:14 Satellite Image Mosaic Selection Problem

The first term on the right-hand side of (11) affects how the part $k$ perceives the resolution of the images to which it belongs. If the image is in a cover, the resolution is equal to $R_j$. If the image is not in a cover, then the resolution is equal to the big constant $B$. As we minimize $r_k$, the lower this first term, the better. This term forced the images with lower resolution to be in a cover. The second term on the right-hand side of (11) is used to force $r_k$ equal to the minimum value of the resolution of the images that contain the part $k$ and are in the cover. When $z_{kj} = 1$ the right-hand side is negative and when 0 the value is positive, and it is the value that $r_k$ takes.

To linearize the incidence angle objective, we need to minimize an auxiliary variable $max_f$ that represents the maximum incidence angle of the images in the cover.

$$ \min max_f \geq x_i F_i, \text{ for all } i = 1, \dots, m \qquad (13) $$

To minimize the area of the clouds, we can model this as a partial set cover problem, where the universe $C = \{1, \dots, c\}$ is formed by all the clouds, and the sets are the images that can cover the clouds. For example, if we have the following set $P_{2c} = \{c_1, c_2, c_5\}$ it means that image 2 can cover clouds 1, 2 and 5, i.e. parts 1, 2 and 5 are not cloudy in image 2. For each cloud $c_i$ we have a variable $y_i$ that is 1 if the cloud is covered or 0 otherwise and $A_c$ indicating the area of the cloud. To maximize the covering of the cloudy areas, we will minimize the following expression:

$$ \min - \sum_{c \in C} y_c A_c \qquad (14) $$

Subject to the following constraint, which forces $y_c$ to be 0 if any of the images that cover it is selected to cover the AOI.

$$ \sum_{i:c \in P_{ic}} x_i \geq y_c, \text{ for all } c \in C \qquad (15) $$

## B MiniZinc Model

We describe the full MINIZINC constraint model implementing the mathematical model given in Section 2.

```c
int: num_images;
int: universe;
int: max_cloud_area;

set of int: IMAGES = 1..num_images;
set of int: UNIVERSE = 1..universe;

array[IMAGES] of set of int: images;
array[IMAGES] of set of int: clouds;
array[IMAGES] of int: costs;
array[UNIVERSE] of int: areas;
array[IMAGES] of int: resolution;
array[IMAGES] of int: incidence_angle;

array[IMAGES] of var bool: taken;

% Which images have a universe 'u' without cloud?
% That is, uclear[u] = {i1, i2, ...} means that the images numbered i1, i2, ... contains 'u' without clouds.
array[UNIVERSE] of set of int: uclear = [{ i | i in IMAGES where not (u
in clouds[i]) /\ u in images[i] } | u in UNIVERSE];
```
```pascal
% Set covering constraint.
constraintforall(u in UNIVERSE) (
  exists(i in IMAGES)(taken[i] /\ u in images[i]));

% cloudy[u] is true iff no image containing a version of 'u' without
clouds is taken.
array[UNIVERSE] of var bool: cloudy;
array[UNIVERSE] of var int: num_clear_images;
constraintforall(u in UNIVERSE) (
  num_clear_images[u] = sum(i in uclear[u])(taken[i])
);
constraintforall(u in UNIVERSE)(cloudy[u] = (num_clear_images[u] == 0));

var int: cloudy_area = sum(u in UNIVERSE)(cloudy[u] * areas[u]);

var int: total_cost = sum(i in IMAGES)(costs[i] * taken[i]);
var int: max_resolve = sum(u in UNIVERSE)(min(i in IMAGES where u in
  images[i] /\ taken[i])(resolution[i]));
var int: max_incidence = max(i in IMAGES)(taken[i] * incidence_angle[i]);

array[1..4] of var int:wide;
constraintwide objects[1] = total_cost;
constraintwide objects[2] = cloudy_area;
constraintwide objects[3] = max_resolve;
constraintwide objects[4] = max_incidence;
```

## C Experimental results detailed

<table><caption>Table 2 Hypervolume values for all the experiments.</caption><thead><tr><th>Instance</th><th>OR-tools default</th><th>OR-tools greedy</th><th>Gecode default</th><th>Gecode greedy</th><th>Gurobi</th></tr></thead><tbody><tr><td>lagos_nigeria_30</td><td><strong>5.46E+33</strong></td><td>5.21E+33</td><td>4.89E+32</td><td>4.91E+33</td><td>4.87E+33</td></tr><tr><td>mexico.city_30</td><td>4.83E+32</td><td>4.82E+33</td><td>4.62E+33</td><td>4.59E+32</td><td><strong>4.83E+33</strong></td></tr><tr><td>paris_30</td><td><strong>1.95E+34</strong></td><td><strong>1.95E+34</strong></td><td>1.59E+34</td><td>1.58E+34</td><td>1.95E+33</td></tr><tr><td>rio_de_janeiro_30</td><td><strong>5.76E+33</strong></td><td>5.65E+33</td><td>5.65E+33</td><td>5.65E+33</td><td><strong>5.76E+33</strong></td></tr><tr><td>tokyo_bay_30</td><td><strong>6.35E+33</strong></td><td>6.33E+33</td><td>6.05E+33</td><td>6.05E+33</td><td>6.30E+33</td></tr><tr><td>lagos_nigeria_50</td><td><strong>3.26E+34</strong></td><td>3.15E+34</td><td>2.19E+34</td><td>2.18E+34</td><td>2.77E+34</td></tr><tr><td>mexico.city_50</td><td>2.88E+33</td><td><strong>2.92E+34</strong></td><td>2.63E+34</td><td>2.57E+34</td><td>2.71E+34</td></tr><tr><td>paris_50</td><td><strong>9.02E+34</strong></td><td>8.85E+34</td><td>7.12E+34</td><td>7.24E+34</td><td>7.85E+34</td></tr><tr><td>rio_de_janeiro_50</td><td>3.94E+33</td><td><strong>3.87E+34</strong></td><td>3.05E+34</td><td>3.04E+34</td><td>3.83E+34</td></tr><tr><td>tokyo_bay_50</td><td><strong>3.32E+34</strong></td><td>3.22E+34</td><td>1.65E+34</td><td>1.78E+34</td><td>1.98E+34</td></tr><tr><td>lagos_nigeria_100</td><td>1.86E+35</td><td><strong>1.87E+35</strong></td><td>1.68E+35</td><td>1.68E+35</td><td>8.55E+34</td></tr><tr><td>mexico.city_100</td><td>1.98E+35</td><td>1.97E+35</td><td><strong>2.11E+35</strong></td><td>2.03E+35</td><td>1.58E+35</td></tr><tr><td>paris_100</td><td><strong>4.73E+35</strong></td><td>3.13E+34</td><td>4.32E+35</td><td>2.91E+35</td><td>4.32E+35</td></tr><tr><td>rio_de_janeiro_100</td><td><strong>4.76E+35</strong></td><td>4.14E+35</td><td>2.28E+35</td><td>1.83E+35</td><td>3.39E+35</td></tr><tr><td>tokyo_bay_100</td><td>1.23E+35</td><td>1.23E+35</td><td>1.89E+35</td><td>1.89E+35</td><td><strong>2.77E+35</strong></td></tr><tr><td>lagos_nigeria_145</td><td>2.98E+36</td><td><strong>3.78E+36</strong></td><td>2.32E+36</td><td>2.32E+36</td><td>9.30E+35</td></tr><tr><td>mexico.city_150</td><td>1.52E+36</td><td>1.11E+36</td><td>2.27E+36</td><td><strong>2.28E+36</strong></td><td>3.42E+35</td></tr><tr><td>paris_150</td><td>2.60E+36</td><td>2.83E+36</td><td>1.11E+36</td><td>1.11E+36</td><td><strong>2.87E+36</strong></td></tr><tr><td>rio_de_janeiro_150</td><td><strong>6.59E+35</strong></td><td>4.81E+34</td><td>4.00E+35</td><td>4.00E+35</td><td>4.61E+34</td></tr><tr><td>tokyo_bay_150</td><td><strong>1.16E+36</strong></td><td>8.68E+35</td><td>8.36E+35</td><td>8.36E+35</td><td>1.12E+36</td></tr><tr><td>mexico.city_200</td><td>2.66E+36</td><td>2.05E+35</td><td><strong>4.28E+36</strong></td><td><strong>4.28E+36</strong></td><td>1.06E+36</td></tr><tr><td>paris_200</td><td>0.00E+00</td><td>0.00E+00</td><td>0.00E+00</td><td>0.00E+00</td><td><strong>1.60E+36</strong></td></tr><tr><td>rio_de_janeiro_200</td><td><strong>0.00E+00</strong></td><td><strong>0.00E+00</strong></td><td><strong>0.00E+00</strong></td><td><strong>0.00E+00</strong></td><td><strong>0.00E+00</strong></td></tr><tr><td>tokyo_bay_200</td><td>0.00E+00</td><td>0.00E+00</td><td>0.00E+00</td><td>0.00E+00</td><td><strong>2.46E+35</strong></td></tr></tbody></table>
<figure><img src="image_446.png" /></figure>
# Partitioning a Map into Homogeneous Contiguous Regions: A Branch-And-Bound Approach Using Decision Diagrams

Nicolas Golenvaux

Institute for Information and Communication Technologies, Electronics and Applied Mathematics (ICTEAM), Université catholique de Louvain, Louvain-la-Neuve, Belgium

Xavier Gillard

Institute for Information and Communication Technologies, Electronics and Applied Mathematics (ICTEAM), Université catholique de Louvain, Louvain-la-Neuve, Belgium

Siegfried Nijssen

Institute for Information and Communication Technologies, Electronics and Applied Mathematics (ICTEAM), Université catholique de Louvain, Louvain-la-Neuve, Belgium

Pierre Schaus

Institute for Information and Communication Technologies, Electronics and Applied Mathematics (ICTEAM), Université catholique de Louvain, Louvain-la-Neuve, Belgium

## Abstract

Regionalization is a crucial spatial analysis technique used for partitioning a map divided into zones into *k* continuous areas, optimizing the similarity of zone attributes within each area. This technique has a variety of applications in fields like urban planning, environmental management, and geographic information systems. The REDCAP algorithm is a well-known approach for addressing the regionalization problem. It consists of two main steps: first, it generates a spatially contiguous tree (SCT) representing the neighborhood structure of the set of spatial objects using a contiguity-constrained hierarchical clustering method. Second, it greedily removes *k* - 1 edges from the SCT to create *k* regions. While this approach has proven to be effective, it may not always produce the most optimal solutions. We propose an alternative method for the second step, an exact dynamic programming (DP) formulation for the k-1 edges removal problem. This DP is solved using a multi-valued decision diagram (MDD)-based branch and bound solver leading to a more optimal solution. We compared our proposed method with the REDCAP state-of-the-art technique on real data and synthetic ones, using different instances of the regionalization problem and different supervised and unsupervised metrics. Our results indicate that our approach provides higher quality partitions than those produced by REDCAP at acceptable computational costs. This suggests that our method could be a viable alternative for addressing the regionalization problem in various applications.

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

**Keywords and phrases** Regionalization, Redcap, Skater, Multivalued Decision Diagrams

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

**Category** Short Paper

## 1 Introduction

Spatial analysis plays a crucial role in comprehending and managing intricate spatial relationships [13]. A fundamental challenge in spatial analysis involves determining homogeneous regions based on similarity computed from shared attributes.

Given a geographical map divided into zones that partition the space, each zone is associated with a set of attributes (e.g., population density, land use, socio-economic factors, etc.) as represented on Figure 1a where the colors represent the attributes. The regionalization

<figure><img src="image_447.png" /></figure>
problem studied in this paper, involves grouping these zones into *k* contiguous areas (also referred to as regions), optimizing the similarity of attributes within each area. The contiguity constraint requires that each grouping forms a single, connected area without any isolated parts or exclaves.

Solving the regionalization problem can be computationally complex, even for a moderate number of zones and areas, as it often entails searching for an optimal solution within a large combinatorial search space. The state-of-the-art method called REDCAP [12] involves solving the problem in a two-step approach as illustrated on Figure 1. In the first step, a spatially contiguous (spanning) tree is created using a hierarchical clustering strategy. In the second step, *k* − 1 edges are greedily deleted from the tree. The final contiguous areas are defined by the *k* remaining connected components of the tree, ensuring both attribute similarity and contiguity.

<figure><img src="image_448.png" /><figcaption>■ **Figure 1** Steps for solving a regionalization problem to cluster four areas: (a) Input of the problem where colors represent the attributes of each zone. (b) Create a spatially contiguous tree that connects all the zones. (c) Remove three edges from the tree. (d) The output areas are then formed.</figcaption></figure>

2 Related Work

The regionalization problem has been a well-studied combinatorial problem since the 1970s
[5]. Regionalization methods can be broadly classified as spatially implicit or spatially explicit
models, depending on how they represent the spatial contiguity constraints of the formed
regions [9]. Implicit methods initially apply traditional or non-spatial clustering methods to
obtain a preliminary solution, which is then adjusted to enforce spatial constraints [16, 17].
Conversely, explicit models enforce spatial contiguity constraints from the outset [9].

Exact methods provide an optimal guarantee for the solutions. However, they are considered computationally intensive and still limited to small problems [8], meaning they are suited for situations with a low number of zones and regions. In contrast, heuristic approaches, which are more scalable, do not guarantee optimal solutions.

Among heuristic approaches with spatially explicit constraints, tree-based methods such as SKATER [14, 2] and REDCAP [12] are widely used and have been demonstrated to generate near-optimal partitions with acceptable computational costs [1, 6]. Both SKATER and REDCAP employ a two-step approach that first constructs a spatially contiguous tree connecting all the zones, then greedily splits it to create the desired number of regions.

Certain spatially explicit regionalization methods, such as those described by [7] and [18],
do not require a predetermined number of regions. Instead, these methods aim to identify
underlying regions while imposing constraints related to these regions.

It is worth noting that the regionalization problem can be viewed as a variant of the
optimal graph partitioning problem [4]. The main difference is that, in graph partitioning,
the dissimilarities between non-adjacent pairs of nodes are not considered by the objective
function.
# 3 Proposed Approach

We start by stating formally the problem. Let us denote by $V = \{1, 2, \dots, n\}$ the zones of our map or more generally the zones. Let $G = (V, E)$ represent the contiguity graph of the map, where edges $(i, j) \in E$ exist if and only if zones $i$ and $j$ share a common border on the map. $G$ must be connected. Let $P = \{V_1, V_2, \dots, V_k\}$ represent a partition (areas) of the zones $V$ into $k$ regions, with $1 \le k \le n$. A partition $P$ is considered feasible if all areas are disjoint, cover the original set of zones, and the induced subgraphs $G(V_u)$ are connected for all $u \in \{1, 2, \dots, k\}$. Let $(a_{i,1}, a_{i,2}, \dots, a_{i,m})$ denote the $m$ numerical attributes of zone $i$. The quadratic distance or dissimilarity between zones based on the set of numerical attributes is $d_{i,j} = \sum_{l=1}^{m} (a_{i,l} - a_{j,l})^2$. The heterogeneity of a region $V_u$ is defined as $h(V_u) = \frac{1}{|V_u|} \sum_{i<j \mid i,j \in V_u} d_{i,j}$. It can be shown that the heterogeneity can equivalently be computed as the sum of squared distances to the mean of attributes $h(V_u) = \sum_{i \in V_u} \sum_{l=1}^{m} (a_{i,l} - \bar{a}_l)^2$ with $\bar{a}_l = \frac{1}{|V_u|} \sum_{i \in V_u} a_{i,l}$. The regionalization problem is to find a feasible $k$-partition minimizing the overall heterogeneity $H(P) = \sum_{u=1}^{k} h(V_u)$.

## 3.1 The REDCAP Two-Step Approach

REDCAP [12] is a state-of-the-art method for solving the regionalization problem. This approach consists of two consecutive steps. First, it identifies a spanning tree, $T$, of graph $G$ (also referred to as a spatially contiguous tree (SCT) in this context) using a hierarchical clustering approach. Second, it identifies $k-1$ edges that partition the tree into a forest of $k$ subtrees, each of which constitutes the final cluster of the regionalization problem.

**Step1.** Starting initially with a set of $n$ clusters $C$, each one containing one of the zones $C = \{c_1 = \{1\}, c_2 = \{2\}, \dots, c_n = \{n\}\}$, the hierarchical clustering approach merges at each step, the two closest contiguous clusters $c_I$ and $c_J$ until one single cluster regroups all the zones. Two clusters $c_I$ and $c_J$ are considered contiguous if there is an edge in the connectivity graph $G$ linking two zones from each cluster $c_I$ and $c_J$. The distance between two clusters $c_I$ and $c_J$ denoted as $D(c_I, c_J)$ and can be computed with different variants. The variant that generally yields the best results is called full-order complete linkage (Full-Order-CLK) defined as $D(c_I, c_J) = \max_{i \in c_I, j \in c_J} d_{i,j}$. Initially empty, one edge is thus added to the SCT $T$ each time two clusters $c_I$ and $c_J$ are merged. This edge $e \in E$ is the one of the original connectivity graph with minimal cost i.e. $\argmin_{(i,j) \in E} d_{i,j}$. At the end of the procedure, $T$ contains $n-1$ edges connecting all the nodes. Overall, the computational complexity for building $T$ using the aforementioned method is $O(n^2 \log n)$.

**Step2.** The second step of REDCAP [12] to obtain $k$ homogeneous regions is to identify $k-1$ edges to remove from $T$ as illustrated in Figure 1c. The remaining components form the final $k$ regions. Due to the inherent complexity of finding an optimal solution for the second-step tree partitioning problem [14], this problem is solved in REDCAP using a greedy heuristic. Let us denote by $F = \{T_1, T_2, \dots, T_k\}$ the spanning forest obtained after the removal of $k-1$ edges from $T$. The set of nodes and edges of each tree $T_u$ are denoted by $V_u$ and $E_u$.

At each iteration, one edge is taken out, splitting one tree of the forest into two trees. Notice that a subtree can possibly contain a single node in case a leaf-edge is removed. The edge that results in the greatest decrease in heterogeneity (or in other words, the highest *homogeneity gain*) is chosen to be eliminated.
For a tree $T_u = (V_u, E_u)$, the homogeneity gain $h_g(e)$ obtained by the removal of an edge $e \in E_u$, dividing the tree into two trees $T_{u1}$ and $T_{u2}$ is defined as $h_g(e) = h(V_u) - h(V_{u1}) - h(V_{u2})$ where $V_{u1}, V_{u2}$ are the nodes in the corresponding sub-trees $T_{u1}$ and $T_{u2}$. The complexity of this greedy algorithm is $O(k \cdot n^2)$ but as $k$ is usually much smaller than $n$, it can be ignored. We propose to replace this second-step greedy algorithm for the tree partitioning problem by an exact formulation using dynamic programming and MDD-based optimization as explained in the next section.

## 3.2 Edge Removal using MDD-based Optimization

We express the problem of the optimal removal of $k-1$ edges from the spanning tree $T = (V_T, E_T)$ to minimize the heterogeneity as a Dynamic Programming Problem. We use a sequence of $k-1$ decision variables $x_i$ representing the edges that are successively removed from the SCT. The domain of each of these variables is the set of edges $E_T$ of the tree. The search-space can be described as Layered Transition Diagram, also called Multivalued Decision Diagrams (MDD) [3]. Let us denote by $F_i$ the set of possible forests obtained by removing exactly $i$ edges from $T$. This set of forests $\cup_{0 \le i \le k-1} F_i$ constitutes the state space and the corresponding nodes of the MDD. A state (forest) is denoted by $f = (V_f, E_f)$. Since we always remove edges on the transitions, but not nodes, the set of nodes of each forest of the state space remains the one of the original contiguity graph $V_f = V, \forall f \in F_i, \forall 0 \le i \le k-1$. Let us now describe the important nodes, the transition function and cost functions for the MDD:

- The set of state-spaces $F = \{F_0, \dots, F_{k-1}\}$ forms the layers, where $F_i$ corresponds to all the states formed by removing exactly $i$ edges from $T$.

- The root of the MDD is denoted as $r$ and corresponds to the state $f^0 \in F_0$, with $f^0 = T$, and its initial value is $v_r = -h(V_T)$, representing the heterogeneity of the entire original map.

- The terminal states are denoted by $t$ and regroup every state $f^{k-1} \in F_{k-1}$.

- The set $\tau$ of transition functions s.t. $\tau_i: F_i \times E \rightarrow F_{i+1}$ for $i = 0, \dots, k-2$ taking the system from one state $f^i$ to the next state $f^{i+1}$ based on the edge removed.

- The set $c$ of transition cost functions $c_i: F_i \times E \rightarrow \mathbb{R}$ s.t. $c_i(f, e)$ is the homogeneity gain $h_g(e)$ of making the decision for $x_i$ to remove the edge $e$ from the forest $f$ at the level $i$. The objective function is then to maximize $v_r + \sum_{i=0}^{k-2} c_i(f^i, x_i)$ so that $f^{i+1} = \tau_i(f^i, x_i)$ and $x_i \in E_T, \forall i \in \{0, \dots, k-2\}; f^i \in F_i, \forall i \in \{0, \dots, k-1\}$. The optimal solution can be obtained by searching the longest path from the root $r$ to one of the terminal nodes $t$.

## 3.3 Branch-and-Bound with MDD

The number of states in the MDD augments rapidly with $k$ and $n$ ($n-1$ choose $k-1$). For such a situation, Bergman et al. [3] have introduced a branch-and-bound (BnB) framework to explore the state space of the MDD without generating it completely upfront and keeping the memory requirement limited. In BnB based on MDDs, relaxed and restricted MDDs, obtained by limiting the width of the MDDs, are used to efficiently explore and prune the solution space. A relaxed MDD is obtained by state-merging. It is an over-approximation of the solution space, where some infeasible solutions might be included. The optimal path of a relaxed MDD provides an upper-bound, which can be used to prune the search space. A restricted MDD, on the other hand, is an under-approximation of the solution space. It is obtained by discarding the less promising states. Some feasible solutions might be excluded but it can nevertheless be used to get lower bounds (similarly to a beam-search). MDD based
BnB enqueues nodes of the original MDD in the queue. When a node is popped, a restricted and a relaxed MDD are compiled from this node to hopefully improve the incumbent solution and prune the search by upper-bounding. This combination of dynamic compilation of restricted and relaxed MDD allows for a more effective exploration of the solution space of the MDD and helps to find the optimal solution in a computationally efficient manner. The nodes can also be pruned by computing a (cheap) upper-bound [11]. We describe next the state-merging procedure and the cheap upper-bound for the optimal edge removal problem.

### 3.3.1 State Merging

As opposed to restricted MDD, a relaxed MDD encodes a superset of the solutions of the original MDD and thus leads to an upper-bound. Its construction is limited to a given width by applying a problem-specific merge operator. In the context of the edge-removal problem, the merge operator applied to two nodes at a same level simply consists in taking the union of the edges in the two forests: merge($f_A = (V, E_{f_A})$, $f_B = (V, E_{f_B})$) = (V, $E_{f_A} \cup E_{f_B}$). The costs of the arcs leading to the merged nodes remain unchanged. Notice that in a relaxed MDD, it is no longer true that the forests $f^i$ at level $i$ have exactly $i+1$ components. One can be convinced that this merging definition correctly includes a superset of the possible paths. It also guarantees an upper-bound on the optimal homogeneity gain that would be obtained without compression. This is a direct consequence of the following property.

**Lemma 1.** For a tree $T = (V, E)$ and a super-tree of $T$ denoted $T' = (V', E')$ with $V \subseteq V'$ and $E \subseteq E'$, let $e$ be an edge present in both $E$ and $E'$. The homogeneity gain of this edge removal in $T$ denoted $h_g(e)$ is lower than $h'_g(e)$ i.e. when this edge is removed from $T'$.

**Proof.** Assuming $e$ connects the two sub-trees $T_1 = (V_1, E_1)$ and $T_2 = (V_2, E_2)$ of $T$ and $T'_1 = (V'_1, E'_1)$ and $T'_2 = (V'_2, E'_2)$ of $T'$. The homogeneity gain in $T$ is $h_g(e) = h(V) - h(V_1) - h(V_2)$ which can equivalently be computed as $h_g(e) = \sum_{i \in V_1} \sum_{j \in V_2} d_{i,j}$. Similarly $h'_g(e) = \sum_{i \in V'_1} \sum_{j \in V'_2} d_{i,j}$. Therefore $h'_g(e) - h_g(e) = \sum_{i \in V'_1 \setminus V_1} \sum_{j \in V'_2 \setminus V'_1} d_{i,j} \ge 0$. ◀

### 3.3.2 Cheap Upper-Bound

To efficiently compute an upper-bound for a state forest $f^i$ at level $i$, one can assume that the $k-i-1$ remaining edges that will be removed induce sub-trees that are perfectly homogeneous (null heterogeneity). The upper-bound on the total homogeneity gain starting from $f^i$ can then be calculated as: $\sum_{l=1}^{\min(k-i-1, |f^i|)} h(V_l^*)$ where $V_l^*$, ..., $V_{k-i-1}^*$ are the $k-i-1$ trees of $f^i$ having the highest heterogeneity value. Notice that the count of trees in $f^i$ may be fewer than $k-i-1$. Should this occur, the formula accounts for the heterogeneity of all present trees in the state.

## 4 Experiments

To assess the effectiveness of our MDD-based method for partitioning spatially contiguous trees, we carried out experiments that compare our approach to REDCAP using both real-life datasets and synthetic datasets with known ground-truth regions.

**Real-life datasets.** We use a set of 5 real-life regionalization datasets varying in size, geometry and number of attributes. Since comparing the absolute value of heterogeneity is meaningless, we use the *Rescaled Overall Heterogeneity H<sub>r</sub>*. That corresponds to the ratio between the overall heterogeneity of the MDD approach and the one of REDCAP.
- **Economic and demographic indicators in NUTS zones:** We collected economic and demographic data on the NUTS areas of Europe to create several datasets that can be used as instances of the regionalization problem. Using the Eurostat database, we gathered data on the density, median age, average GDP per inhabitant, and migration rate for each NUTS-1, NUTS-2, and NUTS-3 European zones in 2019. We constructed one regionalization dataset for each level of NUTS classification, which we refer to as *Ecodemo NUTS1*, *Ecodemo NUTS2*, and *Ecodemo NUTS3*. After removing the unconnected NUTS zones, we are left with 94 zones for the *Ecodemo NUTS1* dataset, 236 for the *Ecodemo NUTS2* dataset, and 1,155 for the *Ecodemo NUTS3* dataset, each having four attributes.

- **Education in Belgium:** We collected data on the level of education in each municipality in Belgium for the year 2017 from the StatBel Open Data. Using this data, we created a regionalization dataset named *Education BE*, where the Belgian municipalities serve as the zones, and their attributes include the share of low, medium, and highly educated inhabitants living in their respective territories. This dataset comprises a total of 563 zones, each having three attributes.

- **USA Ecoregions:** Ecoregions are geographic regions of ecological systems based on vegetation, climate conditions, and land cover [15]. They are frequently used in conservation ecology for planning urban and agricultural development while preserving biodiversity. We use the same dataset as [1] to evaluate our regionalization model. The dataset gathers climatic and land-cover measurements from 1994 in 186 zones of the USA territory. Once the isolated zones are removed, our *Ecoregions USA* dataset includes 172 zones, each one having 15 different ecological attributes.

**Synthetic Maps.** We evaluate the performance of regionalization models in recovering the original regions on synthetic maps using supervised-learning metrics. Our synthetic maps are generated following the methodology proposed in [1], and are parameterized by the number of zones (cells of a square grid), number of regions, region fuzziness, and region geometries.

We created 3 different classes of synthetic maps obtained for different settings of the parameters of the generation process. Each class describes a different level of complexity for regionalization methods. The family A regroups maps with 100 zones divided in 5 regions with simple concentric geometries and well-delimited attribute values. The family B consists of synthetic maps of 400 zones distributed in 10 regions, with more complex and different geometries, and less pronounced frontiers between regions due to their more similar mean values. Finally, the family C comprises maps with 900 zones divided in 20 regions having more complex geometries and more diverse sizes but similar fuzziness than family B. An example of map for each family is represented on Figure 2.

In the case of synthetic datasets, we have access to the original partition, or the ground-truth. Consequently, the optimization problem can be reframed as a machine learning problem, with the goal of recovering the original partition. By comparing the assigned regions with the ground truth, we can calculate various machine learning metrics to assess the performance of the regionalization models. We use the pairwise comparison to evaluate the regionalization algorithms. Each pair of zones is labeled as positive if they belong in the same original region and as negative if they come from different ones. One can then compute the classical True Positive (TP), True Negative (TN), False Positive (FP), False Negative (FN) rates and evaluate the *precision*, *recall* and $F_1$ score metrics to evaluate the regionalization methods on the synthetic maps. In addition to these 3 supervised metrics, we also calculate the ratio between the overall heterogeneity of the partition generated by the regionalization algorithm with the one of ground-truth partition. We name this metric the *True Overall Heterogeneity Ratio* $H_T$.
<figure><img src="image_449.png" /></figure>

Figure 2 Examples of ground-truth regions for each synthetic map class. From left to right, we have an example of the regions' geometry of a synthetic map from family A, family B and family C. The zones belonging to the same region are colored in the same color.

Table 1 Comparison of MDD and REDCAP approaches on real-life regionalization problems. The table presents the rescaled overall heterogeneity for each dataset and partition size ($k \in 5, 10, 15, 20$).

<table><thead><tr><th rowspan="2">Dataset</th><th colspan="4">Number of regions k</th></tr><tr><th>5</th><th>10</th><th>15</th><th>20</th></tr></thead><tbody><tr><td>Ecodemo NUTS1</td><td>1.00</td><td>0.98</td><td>0.96</td><td>0.965</td></tr><tr><td>Ecodemo NUTS2</td><td>0.989</td><td>0.97</td><td>0.958</td><td>0.935</td></tr><tr><td>Ecodemo NUTS3</td><td>0.997</td><td>1.00</td><td>1.00</td><td>1.00</td></tr><tr><td>Education BE</td><td>0.956</td><td>0.861</td><td>0.887</td><td>0.916</td></tr><tr><td>Ecoregions USA</td><td>1.00</td><td>0.997</td><td>0.963</td><td>0.953</td></tr></tbody></table>

# 5 Results

For each regionalization problem, we construct the corresponding SCT using REDCAP's hierarchical clustering. In the second step, we compare the greedy edge removal of REDCAP with the one computed using MDD-based optimization. We employ the DDO solver [10] with a width of 50 for both restricted and relaxed DDs and set a timeout of 100 seconds to identify the final $k$ regions. Before the regionalization process, we normalize the attributes of the zones within a range between 0 and 1 using a MinMax scaler.

## 5.1 Real-life Instances

We assessed the performance of the MDD approach and REDCAP in generating 5, 10, 15, and 20 regions for each dataset described. The comparison between the two methods on the real-life regionalization problems is presented in Table 1. The table provides insights into the quality of the methods' partitions by displaying the rescaled overall heterogeneity for each dataset and for the partition sizes $k = 5, 10, 15, 20$.

Regardless of the requested number of regions, our experimental results show that the MDD approach consistently matches or outperforms REDCAP in terms of overall heterogeneity for each regionalization dataset. This indicates that the partitions generated by the MDD approach are of higher quality. However, the degree of difference varies depending on the dataset. For instance, in the case of the Education BE dataset, the MDD approach reduces the heterogeneity of the partition by almost 15% for 10 regions compared to REDCAP. In contrast, for the Ecodemo NUTS3 dataset, the two algorithms produce similar partitions for all values of $k$. Although the MDD approach only achieved optimality for the Ecodemo NUTS1 dataset when $k = 5$, it outperformed REDCAP for all datasets.
45:8 Partitioning a Map into Homogeneous Regions: A BnB Approach Using DD

Table 2 Number of seconds taken by the MDD approach to find the first solution for each real-life regionalization problem. The time at which it discovers the best solution is presented in parentheses in the cases where it is not equivalent to the first solution founded.

<table><thead><tr><th rowspan="2">Dataset</th><th colspan="4">Number of regions k</th></tr><tr><th>5</th><th>10</th><th>15</th><th>20</th></tr></thead><tbody><tr><td>Ecodemo NUTS1</td><td>&lt;1</td><td>&lt;1</td><td>&lt;1</td><td>&lt;1</td></tr><tr><td>Ecodemo NUTS2</td><td>&lt;1</td><td>&lt;1</td><td>2</td><td>2 (28)</td></tr><tr><td>Ecodemo NUTS3</td><td>6</td><td>13</td><td>19</td><td>21</td></tr><tr><td>Education BE</td><td>2 (18)</td><td>3 (85)</td><td>5 (97)</td><td>6</td></tr><tr><td>Ecoregions USA</td><td>&lt;1</td><td>&lt;1</td><td>&lt;1 (21)</td><td>2 (24)</td></tr></tbody></table>

Table 3 Comparison of our MDD approach with REDCAP on synthetic datasets

<table><thead><tr><th rowspan="2">Metric</th><th colspan="2">A</th><th colspan="2">B</th><th colspan="2">C</th></tr><tr><th>MDD</th><th>REDCAP</th><th>MDD</th><th>REDCAP</th><th>MDD</th><th>REDCAP</th></tr></thead><tbody><tr><td>Precision</td><td><strong>0,988</strong></td><td>0,944</td><td><strong>0,966</strong></td><td>0,912</td><td><strong>0,956</strong></td><td>0,954</td></tr><tr><td>Recall</td><td><strong>0,989</strong></td><td>0,962</td><td><strong>0,949</strong></td><td>0,891</td><td><strong>0,951</strong></td><td>0,908</td></tr><tr><td>F<sub>1</sub> Score</td><td><strong>0,988</strong></td><td>0,952</td><td><strong>0,956</strong></td><td>0,903</td><td><strong>0,954</strong></td><td>0,928</td></tr><tr><td>H<sub>T</sub></td><td><strong>0,99</strong></td><td>1,19</td><td><strong>0,978</strong></td><td>1,108</td><td><strong>1,015</strong></td><td>1,078</td></tr></tbody></table>

Regarding performance, there is a significant difference between the MDD approach and REDCAP. REDCAP takes, on average, between 0.01 and 0.3 seconds to produce a partition depending on the dataset and the number of regions requested, while the MDD approach is allowed to use up to 100 seconds to obtain the best possible solution. However, for the most part, the first partition discovered by the MDD approach is also the best one obtained within the 100-seconds timeframe. Table 2 presents the amount of time the MDD approach searched before finding the first valid partition for each regionalization problem. It is noteworthy that this first partition found by the MDD approach has always a lower or equal overall heterogeneity than the partition found by REDCAP. Additionally, if this first solution is not the best one found within the 100-seconds timeframe, Table 2 reports in parentheses the time taken by the MDD approach to discover the partition with the lowest overall heterogeneity.

## 5.2 Synthetic Maps

We evaluated both the REDCAP and MDD approaches on various synthetic maps of different sizes and complexities (families A, B, and C). For each family, we generated and assessed 20 maps using both methods. Table 3 displays the mean Precision, Recall, F1 Score, and True Overall Heterogeneity Ratio for each method, computed using the ground-truth values. The MDD approach was able to prove the optimal solution for the edge removal problem only for instances belonging to family A.

We can see that using the MDD-based approach for the second step of REDCAP improves the solution in all metrics for the three synthetic map families. Comparing the results of the MDD approach, it can be seen that the precision and recall are lower for families B and C. Both of these families share a characteristic in that their region boundaries are more fuzzy than those of family A. This suggests that our method encounters more difficulty in recovering the initial partition when the attributes between two neighboring regions are more similar, i.e. when the delimitations between regions are less pronounced. The number of zones, the number of regions and their geometric complexity seem to have a lesser impact on the capacity of our model to recover the original regions.
Moreover, for families A and B, the MDD approach generates regions with a lower overall heterogeneity than the ground-truth partition on average. This implies that the original partition is not always the optimal one in terms of overall heterogeneity. Thus, it can be concluded that for the first two families, the regions generated by the MDD approach deviate from the ground-truth ones simply because it discovered partitions with lower overall heterogeneity than the original ones.

# 6 Conclusion

In this paper, we have proposed a novel approach for the regionalization problem, an essential clustering task in a variety of spatial analysis domains. We have improved upon the second step of the well-established REDCAP algorithm by introducing an exact dynamic programming formulation for the edge removal problem solved using a multi-valued decision diagram (MDD)-based branch and bound solver. We have provided comprehensive experiments on both real-life datasets and synthetic ones to illustrate the efficacy of our method. Our comparison with the REDCAP algorithm using a wide range of supervised and unsupervised metrics demonstrated that our approach consistently produces partitions of higher quality.

# References

1 Orhun Aydin, Mark Janikas, Renato Assunção, and Ting-Hwan Lee. A quantitative comparison of regionalization methods. *Int. J. Geogr. Inf. Sci.*, 35:1–29, 2021. doi:10.1080/13658816.2021.1905819.

2 Orhun Aydin, Mark V Janikas, Renato Assunção, and Ting-Hwan Lee. Skater-con: Unsupervised regionalization via stochastic tree partitioning within a consensus framework using random spanning trees. In *Proceedings of the 2nd ACM SIGSPATIAL international workshop on AI for geographic knowledge discovery*, pages 33–42, 2018.

3 David Bergman, Andre A. Cire, Willem-Jan van Hoeve, and J. N. Hooker. Discrete optimization with decision diagrams. *INFORMS Journal on Computing*, 28(1):47–66, 2016. doi:10.1287/ijoc.2015.0648.

4 Aydın Buluç, Henning Meyerhenke, Ilya Safro, Peter Sanders, and Christian Schulz. *Recent advances in graph partitioning*. Springer, 2016.

5 Andrew D Cliff and Peter Haggett. On the efficiency of alternative aggregations in region-building problems. *Environment and Planning A*, 2(3):285–294, 1970.

6 Diep Dao and Jean-Claude Thill. Detecting Attribute-Based Homogeneous Patches Using Spatial Clustering: A Comparison Test. In *Information Fusion and Geographical Information Systems*, pages 37–54. January 2018. doi:10.1007/978-3-319-59539-9_4.

7 Juan C. Duque, Luc Anselin, and Sergio J. Rey. THE MAX-P-REGIONS PROBLEM*. *Journal of Regional Science*, 52(3):397–419, August 2012. doi:10.1111/j.1467-9787.2011.00743.x.

8 Juan C Duque and Richard L Church. A new heuristic model for designing analytical regions. In *North American Meeting of the International Regional Science Association, Seattle*, 2004.

9 Juan Carlos Duque, Raúl Ramos, and Jordi Suriñach. Supervised Regionalization Methods: A Survey. *International Regional Science Review*, 30(3):195–220, July 2007. Publisher: SAGE Publications Inc. doi:10.1177/0160017607301605.

10 X. Gillard, P. Schaus, and V. Coppé. Ddo, a generic and efficient framework for mdd-based optimization. International Joint Conference on Artificial Intelligence (IJCAI-20); DEMO track, 2020.

11 Xavier Gillard, Coppé Vianney, Pierre Schaus, and Ciré André. Improving the filtering of branch-and-bound mdd solvers. In *CPAIOR*, 2021.

12 D. Guo. Regionalization with dynamically constrained agglomerative clustering and partitioning (REDCAP). *International Journal of Geographical Information Science*, 22(7):801–823,
45:10 Partitioning a Map into Homogeneous Regions: A BnB Approach Using DD

July 2008. Publisher: Taylor & Francis _eprint: https://doi.org/10.1080/13658810701674970.
doi:10.1080/13658810701674970.

13 Paul A Longley, Michael F Goodchild, David J Maguire, and David W Rhind. *Geographic information science and systems*. John Wiley & Sons, 2015.

14 Assuncao Martins, Marcos Neves, Gilberto Câmara, and Domingos Da Costa Freitas. Efficient Regionalization Techniques for Socio-Economic Geographical Units Using Minimum Spanning Trees. *International Journal of Geographical Information Science*, 20:797–811, August 2006.
doi:10.1080/13658810600665111.

15 James M Omernik, Shannen S Chapman, Richard A Lillie, Robert T Dumke, et al. Ecoregions of wisconsin. *Transactions of the Wisconsin Academy of Sciences, Arts and Letters*, 88:77–103, 2000.

16 Stan Openshaw. A regionalisation program for large data sets. *Computer Applications*, 3(4):136–147, 1973.

17 Stan Openshaw. Classifying and regionalizing census data. *Census users' handbook*, pages 239–270, 1995.

18 Ran Wei, Sergio Rey, and Elijah Knaap. Efficient regionalization for spatially explicit neighborhood delineation. *International Journal of Geographical Information Science*, 35(1):135–151, 2021. doi:10.1080/13658816.2020.1759806.
# Constraint Programming to Improve Hub Utilization in Autonomous Transfer Hub Networks

Chungjae Lee ⬒ ⬒ ID

H. Milton Stewart School of Industrial and Systems Engineering, Georgia Institute of Technology, Atlanta, GA, USA

Wirattawut Boonbandansook ⬒ ID

H. Milton Stewart School of Industrial and Systems Engineering, Georgia Institute of Technology, Atlanta, GA, USA

Vahid Eghbal Akhlaghi ⬒ ⬒ ID

H. Milton Stewart School of Industrial and Systems Engineering, Georgia Institute of Technology, Atlanta, GA, USA

Kevin Dalmeijer<fnref n="1" /> ⬒ ⬒ ID

H. Milton Stewart School of Industrial and Systems Engineering, Georgia Institute of Technology, Atlanta, GA, USA

Pascal Van Hentenryck ⬒ ⬒ ID

H. Milton Stewart School of Industrial and Systems Engineering, Georgia Institute of Technology, Atlanta, GA, USA

## Abstract

The Autonomous Transfer Hub Network (ATHN) is one of the most promising ways to adapt self-driving trucks for the freight industry. These networks use autonomous trucks for the middle mile, while human drivers perform the first and last miles. This paper extends previous work on optimizing ATHN operations by including transfer hub capacities, which are crucial for labor planning and policy design. It presents a Constraint Programming (CP) model that shifts an initial schedule produced by a Mixed Integer Program to minimize the hub capacities. The scalability of the CP model is demonstrated on a case study at the scale of the United States, based on data provided by Ryder System, Inc. The CP model efficiently finds optimal solutions and lowers the necessary total hub capacity by 42%, saving $15.2M in annual labor costs. The results also show that the reduced capacity is close to a theoretical (optimistic) lower bound.

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

**Keywords and phrases** Constraint Programming, Autonomous Trucking, Transfer Hub Network

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

**Category** Short Paper

**Funding** This research was partly funded through a gift from Ryder and partly supported by the NSF AI Institute for Advances in Optimization (Award 2112533).

**Acknowledgements** Special thanks to the Ryder team for their invaluable support, expertise, and insights.

## 1 Introduction

It is widely believed that autonomous trucks will revolutionize the freight industry, and many companies have started exploring its potential [5, 6, 7, 11, 26, 27]. Some of the major players describe the *transfer hub business model* to be the most likely implementation for

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

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

<footnote n="1">¹ Corresponding author</footnote>
autonomous trucking [21, 25, 28]. An Autonomous Transfer Hub Network (ATHN) is a
network of autonomous truck ports (*transfer hubs*) that leverages the strengths of humans
and automation in their most effective roles. Autonomous trucks handle the monotonous
middle mile to transport goods between the transfer hubs, while humans handle the complex
first and last miles through local cities and deal with customer contacts. According to
Roland Berger [21], implementing a transfer hub model can lower operational costs by 22%
to 40%. A case study in the Southeast of the United States by Ryder System, Inc. and
the Socially Aware Mobility lab [22] reports savings from 27% to 40%, supporting earlier
estimates. The authors model optimizing ATHN operations as a scheduling problem, and a
Constraint Programming (CP) model is used to minimize the empty miles [4]. Subsequent
work by [13] presents a column-generation approach and a bespoke network-flow model to
quickly find high-quality solutions. These findings are combined by [14] into a flow-based
Mixed Integer Programming (MIP) framework that can solve large-scale instances over a
long time horizon (e.g., a month) optimally in reasonable time. These capabilities also allow
for detailed analyses of the benefits and costs of ATHNs.

This paper builds on previous work by introducing a framework for optimizing both ATHN operations and transfer hub capacity utilization. Capacity planning plays a critical role in an actual implementation of the ATHN for operational scheduling. Furthermore, hub capacity has an impact on the number of essential personnel required, which is a crucial factor to consider in labor planning and policy design. Prior research on ATHN operations primarily focused on routing and scheduling, but neglected the capacity considerations involved. Including capacity constraints in a scalable way is not obvious. The CP method in [4] could be extended to incorporate hub capacity constraints by modeling a Resource-Constrained Project Scheduling Problem [9, 12, 20] at every hub and adding cumulative constraints [1, 24] for hub capacity. However, the poor performance reported in [4] makes it unlikely this method will provide good solutions on the national level. The MIP framework from [14] does handle large scale systems, but does not support a cumulative constraint. Alternatives typically require sophisticated modeling techniques and solution methods such as Time-expanded Networks or Dynamic Discretization Discovery [2, 3, 10, 15, 23, 29, 30], which are not obvious to scale either.

The method presented in this work optimizes ATHN operations and improves hub
utilization even for large-scale systems by combining the strengths of MIP and CP. The MIP
model is used to generate routes and an initial schedule, while the CP model is used to shift
the schedule and minimize the required hub capacities. A case study based on real data is
conducted to demonstrate the effectiveness of the new methodology on an ATHN system
spanning the United States for a four-week horizon. The proposed CP model efficiently
finds optimal solutions and lowers the necessary total hub capacity by 42%. This reduction
in capacity may save $15.2M per year in labor cost. Furthermore, it is shown that this
is close to the best possible savings for any initial schedule. This paper also includes a
sensitivity analysis and provides operational insights for future implementation of the ATHN
framework. The remainder of this paper is organized as follows. Section 2 presents the MIP
and CP methodology. Section 3 describes the data and experimental settings used in the
case study. Results and sensitivity analysis are presented in Section 4. Finally, Section 5
provides conclusions and suggests directions for future research.

2 Methodology

This paper uses the methods by [14] to design an ATHN and to optimize the routes of the autonomous vehicles with a MIP. The resulting solution minimizes the cost of the system, but does not take into account the necessary capacity at each of the hubs, which may lead
to low hub utilization. To address this issue, this paper introduces a CP model to shift the schedule in such a way that the original time windows remain satisfied, and the necessary hub capacities are minimized.

## 2.1 ATHN Design and Operations

The input data for the design and optimization of ATHN consists of a set of loads with origins, destinations, and release times. Following [14], K-means clustering is used to determine hub locations $H$, and loads are assigned to hubs according to a hub-assignment rule that minimizes the total driving distance (with autonomous miles discounted by a factor $\gamma$).

To optimize the ATHN operations, first define a set of tasks $T = \{1, 2, ...\}$, where each task $t \in T$ corresponds to moving a load from origin hub $h_t^+ \in H$ to destination hub $h_t^- \in H$ with an autonomous truck. Each task $t \in T$ is associated with a desired pickup time $p_t$ and a flexibility $\Delta \ge 0$, resulting in a pickup time window of $[p_t - \Delta, p_t + \Delta]$. A task graph $G = (V, A)$ is introduced to find the optimal sequence of tasks for every vehicle. The set $V = T \cup \{0\} \cup \{|T| + 1\}$ consists of vertices that correspond to the tasks, together with a source node 0 and a sink node $|T| + 1$. Choosing an arc $a \in A$ indicates that the corresponding tasks are performed sequentially by the same vehicle. An arc between two tasks $t, t' \in T$ represents loading at $h_t^+$, moving freight from $h_t^+$ to $h_t^-$, unloading at $h_t^-$, and relocating from $h_t^-$ to $h_t^+$ to be ready for the next task. Each arc is associated with a corresponding time $\tau_a$. The time for loading or unloading is given by a parameter $\sigma$, and OpenStreetMap times are used for driving and relocation [17].

The cost is calculated in two parts: $d_t$ is the direct cost, which represents the cost of serving task $t \in T$ directly with a conventional truck without using any of the hubs. This cost is taken to be the total distance for delivery and empty return. Note that the current non-autonomous system corresponds to using only direct trips. The second component is a cost differential $c_a$ associated with each arc $a \in A$. For arc $a = (t, t')$ this represents the difference in cost to switch task $t$ from conventional to autonomous delivery, including the relocation cost from $h_t^-$ to $h_{t'}^+$. That is, $d_t$ is the direct trip cost, and $d_t + c_{tt'}$ is the cost to serve task $t \in T$ autonomously and relocate to the start of task $t'$. The autonomous middle miles are discounted by a factor $\alpha \in [0, 1]$ and it is assumed that a fraction $\beta \in [0, 1)$ of the implied first/last miles are empty. Additional details on how $\tau_a$ and $c_a$ are calculated are provided by [14].

Let $y_a$ be a binary variable that indicates that arc $a \in A$ is selected, and let $x_t$ be the start time of task $t \in T$. For convenience, let $\delta_v^+$ and $\delta_v^-$ denote the out-arcs and in-arcs of vertex $v \in V$, respectively. For a given number of autonomous trucks $K$, MIP (1) asks for a set of at most $K$ routes from source to sink that cover different tasks. Objective (1a) minimizes the total cost. If a task is not covered it means it is served by conventional means and the cost is that of a direct trip $d_t$. If the task is covered, the cost differential $c_{tt'}$ is added to calculate the autonomous cost. Constraints (1b) enforce that every task is performed at most once. Constraints (1c) are flow conservation constraints, and Constraint (1d) limits the number of vehicles to $K$. Constraints (1e) are Miller-Tucker-Zemlin constraints [16] that ensure sufficient time passes between subsequent tasks, where $M$ is a sufficiently large constant. These constraints also eliminate cycles. Finally, the variables and their domains are given by Equations (1f) and (1g). The MIP Model (1) is solved with a blackbox solver after applying the acceleration techniques detailed in [14]. The arc-flows ($y$-variables) are translated into a set of routes, which are represented as sequences of tasks, by tracing the flows from the source node 0 to the sink node $|T| + 1$. Note that for given optimal routes, the start times ($x$-variables) are typically not unique. For consistent analysis, the start times are shifted to as early as possible in post processing.
## 46:4 Constraint Programming to Improve Hub Utilization in ATHNs

$$ \min \sum_{t \in T} d_t + \sum_{a \in A} c_a y_a, \qquad (1a) $$

$$ \text{s.t.} \quad \sum_{a \in \delta_t^+} y_a \le 1 \quad \forall t \in T, \qquad (1b) $$

$$ \sum_{a \in \delta_t^+} y_a = \sum_{a \in \delta_t^-} y_a \quad \forall t \in T, \qquad (1c) $$

$$ \sum_{a \in \delta_0^+} y_a \le K, \qquad (1d) $$

$$ x_{t'} \ge x_t + \tau_{tt'} - M(1 - y_{tt'}) \quad \forall t, t' \in T, (t, t') \in A, \qquad (1e) $$

$$ x_t \in [p_t - \Delta, p_t + \Delta] \quad \forall t \in T, \qquad (1f) $$

$$ y_a \in \mathbb{B} \quad \forall a \in A. \qquad (1g) $$

### Figure 1 Mixed Integer Program for Optimizing ATHN Operations.

## 2.2 Minimizing Required Hub Capacities

The routes and schedule obtained from the MIP do not take hub utilization into account. Therefore, it may happen that many trucks are loading and unloading at the same hub at the same time. To address this issue, a CP model is introduced to shift the schedule to minimize the necessary loading/unloading capacity at the hubs, while satisfying the original time windows and maintaining the same route (sequence of tasks) for every vehicle.

The CP model is based on $\mathcal{K}$; the set of routes obtained from MIP (1). Every route $k \in \mathcal{K}$ is split into an sequence of jobs $\mathcal{J}_k = \{1, 2, \dots\}$. For notational convenience, the jobs are numbered sequentially by the order in which they are performed, rather than using the original task numbers. Figure 2 provides a visualization, which will serve as a running example. Each job $j \in \mathcal{J}_k$ has up to four properties:

- type($j$) $\in \{\mathbf{l}, \mathbf{d}, \mathbf{u}, \mathbf{r}, \mathbf{p}\}$; type of job: load, drive, unload, relocate, park, respectively.

- duration($j$) (only for types $\mathbf{l}, \mathbf{d}, \mathbf{u}, \mathbf{r}$); duration of this job.

- task($j$) $\in T$ (only for types $\mathbf{l}, \mathbf{d}, \mathbf{u}$); task associated with this job.

- hub($j$) $\in H$ (only for types $\mathbf{l}, \mathbf{u}, \mathbf{p}$); hub associated with this job.

Each route is modeled with the same repeating sequence of loading at the origin hub ($\mathbf{l}$), waiting before driving ($\mathbf{p}$), driving ($\mathbf{d}$), waiting at the destination hub before unloading ($\mathbf{p}$), unloading ($\mathbf{u}$), waiting before relocating to the next task if any ($\mathbf{p}$), relocating ($\mathbf{r}$), and waiting at the origin hub of the next task until loading ($\mathbf{p}$). Note that jobs of type $\mathbf{l}, \mathbf{d}, \mathbf{u}, \mathbf{r}$ have a fixed duration, while the duration of $\mathbf{p}$ jobs is flexible and can be zero. When no relocation is necessary (the previous destination is equal to the next origin), the $\mathbf{r}$ job is still defined with duration zero for convenience. Every $\mathbf{l}, \mathbf{d}, \mathbf{u}$ job is trivially associated with an original task $t \in T$. Jobs $\mathbf{l}, \mathbf{u}, \mathbf{p}$ for which the vehicle is standing still are associated with a hub as described above.

### CP Model

The CP Model (2) is based on variables $S_j^k$ that indicate the start time of job $j \in \mathcal{J}_k$ in route $k \in \mathcal{K}$, and variables $C_h$ that indicate the necessary loading/unloading capacity at hub $h \in H$. For convenience, $S_{\mathcal{J}_k|+1}^k$ is defined to represent the time at which the final job
<figure><img src="image_452.png" /><figcaption>Figure 2 Jobs $\mathcal{J}_k$ for Route $k \in \mathcal{K}$ (jobs $j, j' \in \mathcal{J}_k$ are used as examples in the main text).</figcaption></figure>

$$ \min \sum_{h \in H} C_h, \qquad (2a) $$

$$ \text{s.t.} \quad I_j^k = \text{Interval}([S_j^k, S_{j+1}^k]) \quad \forall k \in \mathcal{K}, j \in \mathcal{J}_k, \qquad (2b) $$

$$ \text{Cumulative} ([I_j^k | k \in \mathcal{K}, j \in \mathcal{J}_k, \text{type}(j) \in \{\mathbf{l}, \mathbf{u}\}, \text{hub}(j) = \mathbf{h}], C_h) \quad \forall \mathbf{h} \in \mathcal{H}, \qquad (2c) $$

$$ S_{j+1}^k = S_j^k + \text{duration}(j) \quad \forall k \in \mathcal{K}, j \in \mathcal{J}_k, \text{type}(j) \in \{\mathbf{l}, \mathbf{d}, \mathbf{u}, \mathbf{r}\}, \qquad (2d) $$

$$ \text{dom}(S_j^k) = [\underline{s}_j^k, \bar{s}_j^k] \quad \forall k \in \mathcal{K}, j \in \mathcal{J}_k. \qquad (2e) $$

<caption>Figure 3 Constraint Programming Model for Minimizing Required Hub Capacities.</caption>

u is completed (see Figure (2)). Objective (2a) minimizes the total necessary hub capacity. To calculate the hub capacity, the model first defines an interval variable $I_j^k$ for every job (Equation (2b)), which implicitly enforces $S_j^k \le S_{j+1}^k$. Equation (2c) defines a cumulative constraint for every hub $h \in H$ to collect the intervals of the jobs with type l and u at that hub, and to assign the necessary hub capacity to $C_h$. Note that this cumulative constraint has a variable as the capacity. It is also worth mentioning that Constraints (2c) span all the vehicles. Constraints (2d) ensure the correct job duration when a duration is defined (p jobs are flexible). Finally, Equation (2e) defines the domains of the S-variables, where constants $\underline{s}_j^k, \bar{s}_j^k$ remain to be defined. This paper focuses on loading/unloading capacity, but note that CP Model (2) is easily adapted to other objectives, such as minimizing parking space.

To ensure that the time flexibility $\Delta$ is respected, the domains of the S-variables need to be defined accordingly. For job $j \in \mathcal{J}_k$ of route $k \in \mathcal{K}$ and type($j$) = l, the domain is defined around the desired pickup time of the task to match the MIP: $\text{dom}(S_j^k) = [p_{\text{task}(j)} - \Delta, p_{\text{task}(j)} + \Delta]$. This domain is translated to the following d and u jobs to make sure that the flexibility is not exceeded until the task is complete. With slight abuse of notation, this gives the translated domains $\text{dom}(S_j^k) = \text{dom}(S_{j-2}^k) + \text{duration}(j-2)$ for jobs of type($j$) $\in \{\mathbf{d}, \mathbf{u}\}$ (see Figure 2).

## Redundant Bounds and Constraints

Non-trivial bounds for jobs of type r and p are not strictly necessary, but it is straightforward to derive the following:

$$ \text{dom}(S_j^k) = \begin{cases} [\underline{s}_{j-2}^k + \text{duration}(j-2), \bar{s}_{j+2}^k - \text{duration}(j)] & \text{if type}(j) = \mathbf{r}, \\ [\underline{s}_{j-1}^k + \text{duration}(j-1), \bar{s}_{j+1}^k] & \text{if type}(j) = \mathbf{p}. \end{cases} \qquad (3) $$

Relocation can only start after unloading is completed, and has to start to be in time for the next loading (see job $j'$ in Figure 2). In a similar way, parking can only start after the previous job is completed, and parking for duration zero is possible until the upper bound of the next job.
<table><caption>Table 1 Baseline Parameter Values for the Case Study.</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>n</td><td>6842 loads</td></tr><tr><td>|H|</td><td>100 transfer hubs</td></tr><tr><td>γ</td><td>40% discount for autonomous mileage during hub-assignment</td></tr><tr><td>β</td><td>25% first/last-mile inefficiency</td></tr><tr><td>α</td><td>25% discount for autonomous mileage</td></tr><tr><td>Δ</td><td>1 hour pickup-time flexibility</td></tr><tr><td>σ</td><td>30 minutes autonomous truck loading/unloading time</td></tr><tr><td>K</td><td>100 autonomous trucks</td></tr></tbody></table>

For the current Objective (2a), the following observation is used to reduce the search space: To minimize the loading/unloading capacity, it does not matter when relocation takes place between the jobs u and 1. It is therefore possible without loss of generality to impose that relocation starts immediately after unloading:

$$S_j^k = S_{j-1}^k \quad \forall k \in \mathcal{K}, j \in \mathcal{J}_k, \text{type}(j) = \mathbf{r}. \qquad (4)$$

## 3 Case Study

This paper conducts a realistic case study based on data from Ryder System, Inc. (Ryder), one of the largest transportation and logistics companies in North America. Ryder has provided a dataset that is representative for its dedicated transportation business in the US, reducing the scope to orders that are strong candidates for automation. Following [14], the case study focuses on orders that are *challenging* in the sense that they would currently induce an empty return trip. Every order represents a load with an origin, a destination, and a scheduled release time. The hubs are chosen based on data from October to December 2019, while the experiments are based on 6842 loads in the first four weeks of October.

### Experimental Settings

The parameter settings are taken from [14] and summarized in Table 1. All methods from Section 2 are implemented in Python 3.9. The MIP models are solved with Gurobi 9.5.2, and the CP models are solved with CP-SAT 9.6 [19]. All time-related data is rounded to the nearest minute to facilitate the integer domains that are required by CP-SAT. The experiments are conducted on a Linux machine with dual Intel Xeon Gold 6226 CPUs on the PACE Pheonix cluster [18]. Each experiment is assigned to use at most 24 cores and 192GB of RAM. If Gurobi runs out of memory, the experiment is repeated on a machine with 384GB of RAM, and if that fails, the number of cores is halved until the solver terminates successfully. CP-SAT did not encounter memory issues at any point. The Gurobi time limit is set to three hours, except for the sensitivity analysis for flexibility Δ, which is given 12 hours to obtain better solutions and is warm started with a MIP start. The CP-SAT solver is warm started with the solution obtained by the MIP, which is provided as a hint.

## 4 Results

Figure 4 summarizes the results of using the CP model to minimize the necessary hub capacity for loading and unloading trucks. “Before CP” shows the required capacity if the MIP solution were implemented immediately, while “After CP” shows the results after applying the CP model. The ATHN problem was solved to optimality for all instances except
<figure><img src="image_453.png" /><figcaption>Figure 4 Capacity Before and After CP.</figcaption></figure>

<figure><img src="image_454.png" /><figcaption>Figure 5 Loads Served Autonomously.</figcaption></figure>

<figure><img src="image_455.png" /><figcaption>(a) Before CP Optimization.</figcaption></figure>

<figure><img src="image_456.png" /><figcaption>(b) After CP Optimization.</figcaption></figure>

<caption>Figure 6 Capacity Reduction for the 100 Truck Case (circle area proportional to hub capacity).</caption>

when $K = 50$, which remained at a 0.05% optimality gap within the given time limit. For the CP problem, all instances were solved to optimality in under 30 seconds. The figure shows that the CP model is highly effective in reducing the total hub capacity for all instances, reducing the necessary capacity by 31% up to 50%. For the base case of $K = 100$ the CP model can reduce capacity by 42%. Note that this may correspond to significant monetary savings: If each unit of loading/unloading capacity requires the assistance of a mechanic around the clock (three shifts of $57,557 per year [8]), the CP model reduces the annual labor cost by $15.2 million. Figure 4 also shows that the resulting capacity is close to the lower bound for any ATHN solution, which is obtained by removing the time constraints between subsequent tasks, i.e., treat every task as if it is the only task on the route.

As the number of trucks increases, the ATHN starts serving more loads autonomously, as shown by Figure 5. Without the CP model, this leads to a substantial increase in hub capacity as more loads are added to the system. The CP model completely mitigates this effect, and is able to maintain an almost stable hub capacity as the workload increases. This reveals a suprising robustness to accomodate new orders that the CP model is able to exploit. When making investment decisions and hiring personnel to operate the hubs, this is a very desirable property. The maps in Figure 6 demonstrate that the necessary capacity is reduced throughout the system, and the peaks in the South and Northeast have decreased significantly. The capacity at the hub near Louisville, Kentucky for example was reduced from 7 to 2.

Interestingly, only a small portion of the jobs needs to be rescheduled to obtain the substantial savings in hub capacity. E.g., in the 100-truck case only 13% of loadings were rescheduled. Figure 7 provides histograms for the absolute size of the shift for the loading and unloading times that were shifted. It can be seen that a majority of these jobs were
<figure><img src="image_457.png" /><figcaption>(a) Histogram of Shifted Loading Times.</figcaption></figure>

<caption>(b) Histogram of Shifted Unloading Times.</caption>

Figure 7 Shifted Schedule Times for the 100 Truck Case.

<figure><img src="image_458.png" /><figcaption>(a) Impact of Loading/Unloading Time $\sigma$.</figcaption></figure>

<caption>(b) Impact of Flexibility $\Delta$.</caption>

Figure 8 Sensitivity Analysis Loading/Unloading Time and Flexibility for the 100 Truck Case.

moved by no more than 30 minutes, which shows that the ATHN solutions are sufficiently flexible to be adjusted without propagating delays through the schedule. The fact that the CP model is easy to solve and only makes small modifications also makes it an attractive tool at the operational level: if order details are changed or delays are encountered, the CP model can quickly be re-solved to avoid causing overlap in loading and unloading at the hubs.

### Sensitivity Analysis

The loading/unloading time $\sigma$ and the flexibility $\Delta$ are two significant factors that affect the ATHN. Figure 8 summarizes how these parameters impact the total hub capacity and the CP model's performance. The ATHN is solved for loading/unloading time from zero to two hours, and for no flexibility up to one day of flexibility. Note that loading/unloading time up to two hours may be realistic if the vehicle is inspected every time it leaves or enters a hub.

For $\sigma$, the ATHN solver found optimal solutions for all instances and the CP problem found optimal schedules within 15 seconds. As $\sigma$ increases, the total hub capacity increases because jobs overlap more frequently. The CP model is not able to fully compensate for this effect, but can still improve the hub capacity by at least 28%.

For $\Delta$, optimal solutions were found for $\Delta \in \{0h, 1h\}$, a near-optimal solution within 0.05% optimality gap was found for $\Delta = 2h$, while $\Delta \in \{12h, 24h\}$ remained at a 3% optimality gap. The CP solver took 7 minutes for $\Delta = 12h$ and 13 minutes for $\Delta = 24h$.