# Efficient Enumeration of Fixed Points in Complex Boolean Networks Using Answer Set Programming

Van-Giang Trinh, Belaid Benhamou, Sylvain Soliman

► To cite this version:

Van-Giang Trinh, Belaid Benhamou, Sylvain Soliman. Efficient Enumeration of Fixed Points in Complex Boolean Networks Using Answer Set Programming. 29th International Conference on Principles and Practice of Constraint Programming (CP 2023), Aug 2023, Toronto, Canada. pp.35:1–35:19, 10.4230/LIPIcs.CP.2023.35 . hal-04209296

HAL Id: hal-04209296

https://amu.hal.science/hal-04209296v1

Submitted on 22 Sep 2023

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
29th International Conference on
Principles and Practice of
Constraint Programming

CP 2023, August 27–31, 2023, Toronto, Canada

Edited by

Roland H. C. Yap

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

LIPICS

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

TORONTO, CANADA
Editors

Roland H. C. Yap

National University of Singapore, School of Computing, 13 Computing Drive, Singapore
ryap@comp.nus.edu.sg

ACM Classification 2012

Theory of computation → Constraint and logic programming

ISBN 978-3-95977-300-3

Published online and open access by

Schloss Dagstuhl – Leibniz-Zentrum für Informatik GmbH, Dagstuhl Publishing, Saarbrücken/Wadern, Germany. Online available at https://www.dagstuhl.de/dagpub/978-3-95977-300-3.

Publication date

September, 2023

Bibliographic information published by the Deutsche Nationalbibliothek

The Deutsche Nationalbibliothek lists this publication in the Deutsche Nationalbibliografie; detailed bibliographic data are available in the Internet at https://portal.dnb.de.

License

This work is licensed under a Creative Commons Attribution 4.0 International license (CC-BY 4.0):
https://creativecommons.org/licenses/by/4.0/legalcode.

In brief, this license authorizes each and everybody to share (to copy, distribute and transmit) the work under the following conditions, without impairing or restricting the authors' moral rights:

Attribution: The work must be attributed to its authors.

The copyright is retained by the corresponding authors.

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

Digital Object Identifier: 10.4230/LIPIcs.CP.2023.0

ISBN 978-3-95977-300-3

ISSN 1868-8969

https://www.dagstuhl.de/lipics
LIPIcs – Leibniz International Proceedings in Informatics

LIPIcs is a series of high-quality conference proceedings across all fields in informatics. LIPIcs volumes are published according to the principle of Open Access, i.e., they are available online and free of charge.

Editorial Board

- Luca Aceto (Chair, Reykjavik University, IS and Gran Sasso Science Institute, IT)

- Christel Baier (TU Dresden, DE)

- Roberto Di Cosmo (Inria and Université de Paris, FR)

- Faith Ellen (University of Toronto, CA)

- Javier Esparza (TU München, DE)

- Daniel Král' (Masaryk University, Brno, CZ)

- Meena Mahajan (Institute of Mathematical Sciences, Chennai, IN)

- Anca Muscholl (University of Bordeaux, FR)

- Chih-Hao Luke Ong (University of Oxford, GB)

- Phillip Rogaway (University of California, Davis, US)

- Eva Rotenberg (Technical University of Denmark, Lyngby, DK)

- Raimund Seidel (Universität des Saarlandes, Saarbrücken, DE and Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Wadern, DE)

- Pierre Senellart (ENS, Université PSL, Paris, FR)

ISSN 1868-8969

https://www.dagstuhl.de/lipics
<figure><img src="image_5.png" /></figure>
<figure><img src="image_6.png" /></figure>

# Contents

## Preface

Roland H.C. Yap 0:ix-0:x

## Senior Program Committee

## Program Committee

## Additional Reviewers

## Authors

## Invited Talks

### Beyond Optimal Solutions for Real-World Problems

Maria Garcia de la Banda 1:1-1:4

### A Tale of Two Cities: Teaching CP with Story-Telling

Jimmy H.M. Lee 2:1-2:1

### The CP-SAT-LP Solver

Laurent Perron, Frédéric Didier, and Steven Gay 3:1-3:2

### Coupling CP with Deep Learning for Molecular Design and SARS-CoV2 Variants Exploration

Thomas Schiex 4:1-4:3

### CP Solver Design for Maximum CPU Utilization

Petr Vilím 5:1-5:1

## Regular Papers

### Optimization of Short-Term Underground Mine Planning Using Constraint Programming

Younes Aalian, Gilles Pesant, and Michel Gamache 6:1-6:16

### Exploiting Configurations of MaxSAT Solvers

Josep Alòs, Carlos Ansótegui, Josep M. Salvia, and Eduard Torres 7:1-7:16

### Symmetries for Cube-And-Conquer in Finite Model Finding

João Araújo, Choiwah Chow, and Mikoláš Janota 8:1-8:19

### Guiding Backtrack Search by Tracking Variables During Constraint Propagation

Gilles Audemard, Christophe Lecoutre, and Charles Prud'homme 9:1-9:17

### Incremental Constrained Clustering by Minimal Weighted Modification

Aymeric Beauchamp, Thi-Bich-Hanh Dao, Samir Loudni, and Christel Vrain 10:1-10:22

<figure><img src="image_7.png" /></figure>
<table>
  <tbody>
    <tr>
      <td>Simplifying Step-Wise Explanation Sequences<br><i>Ignace Bleukx, Jo Devriendt, Emilio Gamba, Bart Bogaerts, and Tias Guns</i></td>
      <td>11:1-11:20</td>
    </tr>
    <tr>
      <td>Towards More Efficient Local Search for Pseudo-Boolean Optimization<br><i>Yi Chu, Shaowei Cai, Chuan Luo, Zhendong Lei, and Cong Peng</i></td>
      <td>12:1-12:18</td>
    </tr>
    <tr>
      <td>Boosting Decision Diagram-Based Branch-And-Bound by Pre-Solving with Aggregate Dynamic Programming<br><i>Vianney Coppé, Xavier Gillard, and Pierre Schaus</i></td>
      <td>13:1-13:17</td>
    </tr>
    <tr>
      <td>Fast Matrix Multiplication Without Tears: A Constraint Programming Approach<br><i>Arnaud Deza, Chang Liu, Pashootan Vaezipoor, and Elias B. Khalil</i></td>
      <td>14:1-14:15</td>
    </tr>
    <tr>
      <td>Probabilistic Inference by Projected Weighted Model Counting on Horn Clauses<br><i>Alexandre Dubray, Pierre Schaus, and Siegfried Nijssen</i></td>
      <td>15:1-15:17</td>
    </tr>
    <tr>
      <td>A CP Approach for the Liner Shipping Network Design Problem<br><i>Yousra El Ghazi, Djamal Habet, and Cyril Terrioux</i></td>
      <td>16:1-16:21</td>
    </tr>
    <tr>
      <td>Optimization Models for Pickup-And-Delivery Problems with Reconfigurable Capacities<br><i>Arnoosh Golestanian, Giovanni Lo Bianco, Chengyu Tao, and J. Christopher Beck</i></td>
      <td>17:1-17:17</td>
    </tr>
    <tr>
      <td>Preprocessing in SAT-Based Multi-Objective Combinatorial Optimization<br><i>Christoph Jabs, Jeremias Berg, Hannes Ihalainen, and Matti Järvisalo</i></td>
      <td>18:1-18:20</td>
    </tr>
    <tr>
      <td>An Efficient Constraint Programming Approach to Preemptive Job Shop Scheduling<br><i>Carla Juvin, Emmanuel Hebrard, Laurent Houssin, and Pierre Lopez</i></td>
      <td>19:1-19:16</td>
    </tr>
    <tr>
      <td>Horizontally Elastic Edge Finder Rule for Cumulative Constraint Based on Slack and Density<br><i>Roger Kameugne, Sévérine Fetgo Betmbe, Thierry Noulamo, and Clémentin Tayou Djamegni</i></td>
      <td>20:1-20:17</td>
    </tr>
    <tr>
      <td>Exploring Hydrogen Supply/Demand Networks: Modeller and Domain Expert Views<br><i>Matthias Klapperstueck, Frits de Nijs, Ilankaikone Senthooran, Jack Lee-Kopij, Maria Garcia de la Banda, and Michael Wybrow</i></td>
      <td>21:1-21:18</td>
    </tr>
    <tr>
      <td>Binary Constraint Trees and Structured Decomposability<br><i>Petr Kučera</i></td>
      <td>22:1-22:19</td>
    </tr>
    <tr>
      <td>Large Neighborhood Beam Search for Domain-Independent Dynamic Programming<br><i>Ryo Kuroiwa and J. Christopher Beck</i></td>
      <td>23:1-23:22</td>
    </tr>
    <tr>
      <td>MDD Archive for Boosting the Pareto Constraint<br><i>Steve Malalel, Arnaud Malapert, Marie Pelleau, and Jean-Charles Régin</i></td>
      <td>24:1-24:15</td>
    </tr>
    <tr>
      <td>Learning a Generic Value-Selection Heuristic Inside a Constraint Programming Solver<br><i>Tom Marty, Tristan François, Pierre Tessier, Louis Gautier, Louis-Martin Rousseau, and Quentin Cappart</i></td>
      <td>25:1-25:19</td>
    </tr>
  </tbody>
</table>
<table><tr><td>Proof Logging for Smart Extensional Constraints<br><i>Matthew J. McIlree and Ciaran McCreesh</i></td><td>26:1-26:17</td></tr><tr><td>Improving Conflict Analysis in MIP Solvers by Pseudo-Boolean Reasoning<br><i>Gioni Mexi, Timo Berthold, Ambros Gleixner, and Jakob Nordström</i></td><td>27:1-27:19</td></tr><tr><td>Using Canonical Codes to Efficiently Solve the Benzenoid Generation Problem<br>with Constraint Programming<br><i>Xiao Peng and Christine Solnon</i></td><td>28:1-28:17</td></tr><tr><td>Distribution Optimization in Constraint Programming<br><i>Guillaume Perez, Gaël Glorian, Wijnand Suijlen, and Arnaud Lallouet</i></td><td>29:1-29:19</td></tr><tr><td>The p-Dispersion Problem with Distance Constraints<br><i>Nikolaos Ploskas, Kostas Stergiou, and Dimosthenis C. Tsouros</i></td><td>30:1-30:18</td></tr><tr><td>Partially Preemptive Multi Skill/Mode Resource-Constrained Project Scheduling<br>with Generalized Precedence Relations and Calendars<br><i>Guillaume Povéda, Nahum Alvarez, and Christian Artigues</i></td><td>31:1-31:21</td></tr><tr><td>Assembly Line Preliminary Design Optimization for an Aircraft<br><i>Stéphanie Roussel, Thomas Polacsek, and Anouck Chan</i></td><td>32:1-32:19</td></tr><tr><td>SAT-Based Learning of Compact Binary Decision Diagrams for Classification<br><i>Pouya Shati, Eldan Cohen, and Sheila McIlraith</i></td><td>33:1-33:19</td></tr><tr><td>Constraint Programming with External Worst-Case Traversal Time Analysis<br><i>Pierre Talbot, Tingting Hu, and Nicolas Navet</i></td><td>34:1-34:20</td></tr><tr><td>Efficient Enumeration of Fixed Points in Complex Boolean Networks Using<br>Answer Set Programming<br><i>Van-Giang Trinh, Belaid Benhamou, and Sylvain Soliman</i></td><td>35:1-35:19</td></tr><tr><td>Guided Bottom-Up Interactive Constraint Acquisition<br><i>Dimosthenis C. Tsouros, Senne Berden, and Tias Guns</i></td><td>36:1-36:20</td></tr><tr><td>Addressing Problem Drift in UNHCR Fund Allocation<br><i>Sameela Suharshani Wijesundara, Maria Garcia de la Banda, and Guido Tack</i></td><td>37:1-37:18</td></tr><tr><td>From Formal Boosted Tree Explanations to Interpretable Rule Sets<br><i>Jinqiang Yu, Alexey Ignatiev, and Peter J. Stuckey</i></td><td>38:1-38:21</td></tr><tr><td>Searching for Smallest Universal Graphs and Tournaments with SAT<br><i>Tianwei Zhang and Stefan Szeider</i></td><td>39:1-39:20</td></tr><tr><td>FastMapSVM for Predicting CSP Satisfiability<br><i>Kexin Zheng, Ang Li, Han Zhang, and T. K. Satish Kumar</i></td><td>40:1-40:17</td></tr><tr><td>Improving Local Search for Pseudo Boolean Optimization by Fragile Scoring<br>Function and Deep Optimization<br><i>Wenbo Zhou, Yujiao Zhao, Yiyuan Wang, Shaowei Cai, Shimao Wang, Xinyu Wang, and Minghao Yin</i></td><td>41:1-41:18</td></tr></table>
Short Papers

<table>
  <tbody>
    <tr>
      <td>Predict-Then-Optimise Strategies for Water Flow Control<br><i>Vincent Barbosa Vaz, James Bailey, Christopher Leckie, and Peter J. Stuckey</i></td>
      <td>42:1-42:10</td>
    </tr>
    <tr>
      <td>Constraint Programming Models for Depth-Optimal Qubit Assignment and<br>SWAP-Based Routing<br><i>Kyle E. C. Booth</i></td>
      <td>43:1-43:10</td>
    </tr>
    <tr>
      <td>Constraint Model for the Satellite Image Mosaic Selection Problem<br><i>Manuel Combarro Simón, Pierre Talbot, Grégoire Danoy, Jedrzej Musial, Mohammed Alswaitti, and Pascal Bouvry</i></td>
      <td>44:1-44:15</td>
    </tr>
    <tr>
      <td>Partitioning a Map into Homogeneous Contiguous Regions: A Branch-And-Bound<br>Approach Using Decision Diagrams<br><i>Nicolas Golenvaux, Xavier Gillard, Siegfried Nijssen, and Pierre Schaus</i></td>
      <td>45:1-45:10</td>
    </tr>
    <tr>
      <td>Constraint Programming to Improve Hub Utilization in Autonomous Transfer<br>Hub Networks<br><i>Chungjae Lee, Wirattawut Boonbandansook, Vahid Eghbal Akhlaghi, Kevin Dalmeijer, and Pascal Van Hentenryck</i></td>
      <td>46:1-46:11</td>
    </tr>
    <tr>
      <td>A New Approach to Finding 2 x n Partially Spatially Balanced Latin Rectangles<br><i>Renee Mirka, Laura Greenstreet, Marc Grimson, and Carla P. Gomes</i></td>
      <td>47:1-47:11</td>
    </tr>
    <tr>
      <td>Proven Optimally-Balanced Latin Rectangles with SAT<br><i>Vaidyanathan Peruvemba Ramaswamy and Stefan Szeider</i></td>
      <td>48:1-48:10</td>
    </tr>
    <tr>
      <td>Enumerative Level-2 Solution Counting for Quantified Boolean Formulas<br><i>Andreas Plank, Sibylle Möhle, and Martina Seidl</i></td>
      <td>49:1-49:10</td>
    </tr>
  </tbody>
</table>
<figure><img src="image_8.png" /></figure>

# Preface

This volume contains the proceedings of the 29th International Conference on Principles and Practice of Constraint Programming (CP 2023), which was held in Toronto, Canada, August 27–31, 2023. More details of the conference can be found at https://cp2023.a4cp.org/index.html.

Held annually, CP is the premier international conference on constraint programming. CP is concerned with all aspects of computing with constraints, including, but not restricted to: theory, algorithms, environments, languages, models, systems, and applications.

As is customary for CP, papers could be submitted to multiple tracks. CP 2023 had the following tracks:

- Technical Track

- Applications Track

- Machine Learning Track

- Operations Research Track

- Trustworthy Decision Making

Chair: Roland Yap

Chair: Helmut Simonis

Chair: Tias Guns

Chair: Gilles Pesant

Chair: Peter Stuckey

In addition, there was a SAT Fast Track to synchronize with the SAT 2023 conference.

A total of 109 papers (excluding abstracts) were submitted to these tracks. Authors could submit either a full paper with a maximum length of 15 pages (references excluded) or a short paper with a maximum length of 8 pages (references excluded). Each paper was reviewed with a senior Program Committee together with Program Committees and additional reviewers recruited by the Program Committee. The track chairs managed the review process for their respective tracks. All papers had at least three reviews. Authors had the opportunity to answer questions from reviewers in an author response phase. Based on extensive discussion on the papers from reviewers, program and senior program committee, and track chairs taking into account reviews and author responses, a total of 44 papers were accepted. A meta-review was prepared for each paper by a senior program committee member summarizing the decision with suggestions to authors. The Senior Program Committee and Chairs nominated papers for the best paper prizes. A select committee from the Senior Program Committee together with the Program Chair awarded the Best Paper Prize to Mathew J. McIlree and Ciaran McCreesh for “Proof Logging for Smart Extensional Constraints” and the best application paper prize to Matthias Klapperstueck, Frits De Nijs, Ilankaikone Senthooran, Jack Lee-Kopij, Maria Garcia De La Banda and Michael Wybrow for “Exploring Hydrogen Supply/Demand Networks: Modeller and Domain Expert views”.

In addition to the paper tracks, the conference had a number of satellite events. Lars Kotthoff (University of Wyoming) organized the workshops on the first day of the conference with five workshops. The doctoral program, also on the first day, was organized by Xavier Gillard (Université catholique de Louvain) and forms an important part of the conference to give students an environment to present their research with discussions with senior researchers, and networking activities combined with a poster presentation during the conference reception. CP features two tutorials organized by Emir Demirović (TU Delft) on the timely topics of explainable constraint solving and machine learning for solvers. The conference program featured four invited talks given by Maria Garcia de la Banda, Jimmy Lee, Laurent Perron, Thomas Schiex and Petr Vilím. The talks were selected to showcase a range of research in CP and include constraint solver design, designing proteins with applications to the SARS-Cov2 virus, teaching of CP, and making optimisation technology more usable.
Many people have contributed to make the conference a success. The conference would not be possible for the hard work of the authors in submitting high-quality scientific work which forms the basis of the conference proceedings. The program and senior committee together with track chairs had the challenging tasks of selecting papers as well as providing authors with suggestions on paper improvements. A special thanks goes to Andre Cire and Eldan Cohen, the conference and local chairs respectively, both at the University of Toronto, who made CP possible in Toronto. In addition to the above mentioned conference organizers, I would like to thank the following co-organizers. The diversity, equity, and inclusion (DEI) chairs, Maria Andreina Francisco Rodriguez (Uppsala University) and Andrea Rendl (Satalia). Arvind Raghunathan (Mitsubishi Electric Research Laboratories), the sponsorship chair. Anna Latour (National University of Singapore) was instrumental in managing the website and also handled the publicity and social media.

I would also like to thank the Association for Constraint Programming (ACP) which makes the CP conference series possible. The conference is grateful to the ACP president, David Bergman (University of Connecticut) and the ACP conference coordinator, Hélène Verhaeghe (KU Leuven) for their support and feedback.

The conference acknowledges the generous support of all our sponsors:

- ACP

- The Artificial Intelligence Journal (Elsevier)

- Cosling

- Department of Management, University of Toronto Scarborough

- Google

- IBM

- MERL

- The Optimization Firm

- ScheduleOpt

July 2023, Singapore

Roland Yap
<figure><img src="image_9.png" /></figure>

# Organization

## Senior Program Committee

Christian Bessiere

Mats Carlsson

Berthe Choueiry

David Cohen

Maria Garcia De La Banda

Ian Gent

Philip Kilby

Christophe Lecoutre

Jimmy Lee

Kuldeep Meel

Ian Miguel

Michela Milano

Barry O'Sullivan

Justin Pearson

Thomas Schiex

Laurent Simon Labri

Christine Solnon

Peter J. Stuckey

Michael Trick

Neil Yorke-Smith

CNRS, University of Montpellier, France

RISE Research Institutes of Sweden, Sweden

University of Nebraska-Lincoln, USA

Royal Holloway, University of London

Monash University, Australia

University of St Andrews, United Kingdom

Data61 & The Australian National University

CRIL, University of Artois, France

The Chinese University of Hong Kong, Hong Kong

National University of Singapore, Singapore

University of St Andrews, United Kingdom

University of Bologna, Italy

University College Cork, Ireland

Uppsala University, Sweden

Universite Fédérale de Toulouse, ANITI, INRAE, France

Bordeaux Institute of Technology, France

INSA Lyon, France

Monash University, Australia

Carnegie Mellon University, USA

Delft University of Technology, Netherlands

<figure><img src="image_10.png" /></figure>
Program Committee

<table>
  <tr>
    <td>Carlos Ansotegui</td>
    <td>Universitat de Lleida, Spain</td>
  </tr>
  <tr>
    <td>Christopher Beck</td>
    <td>University of Toronto, Canada</td>
  </tr>
  <tr>
    <td>Nicolas Beldiceanu</td>
    <td>IMT Atlantique (LS2N), France</td>
  </tr>
  <tr>
    <td>David Bergman</td>
    <td>University of Connecticut, USA</td>
  </tr>
  <tr>
    <td>Armin Biere</td>
    <td>University of Freiburg, Germany</td>
  </tr>
  <tr>
    <td>Nikolaj Bjorner</td>
    <td>Microsoft, USA</td>
  </tr>
  <tr>
    <td>Ken Brown</td>
    <td>University College Cork</td>
  </tr>
  <tr>
    <td>Quentin Cappart</td>
    <td>Ecole Polytechnique de Montréal, Canada</td>
  </tr>
  <tr>
    <td>Martin Cooper</td>
    <td>IRIT - Universite Paul Sabatier, France</td>
  </tr>
  <tr>
    <td>Simon de Givry</td>
    <td>INRA - MIAT, France</td>
  </tr>
  <tr>
    <td>Jean-Guillaume Fages</td>
    <td>COSLING, France</td>
  </tr>
  <tr>
    <td>Pierre Flener</td>
    <td>Uppsala University, Sweden</td>
  </tr>
  <tr>
    <td>Emmanuel Hebrard</td>
    <td>LAAS, CNRS, France</td>
  </tr>
  <tr>
    <td>Alexey Ignatiev</td>
    <td>Monash University, Australia</td>
  </tr>
  <tr>
    <td>George Katsirelos</td>
    <td>MIA Paris, INRAE, AgroParisTech, France</td>
  </tr>
  <tr>
    <td>Zeynep Kiziltan</td>
    <td>University of Bologna, Italy</td>
  </tr>
  <tr>
    <td>T. K. Satish Kumar</td>
    <td>University of Southern California, USA</td>
  </tr>
  <tr>
    <td>Mikael Lagerkvist</td>
    <td>Optischedule, Sweden</td>
  </tr>
  <tr>
    <td>Nadjib Lazaar</td>
    <td>UM2-LIRMM, France</td>
  </tr>
  <tr>
    <td>Chu-Min Li</td>
    <td>Université de Picardie Jules Verne, France</td>
  </tr>
  <tr>
    <td>Michele Lombardi</td>
    <td>University of Bologna, Italy</td>
  </tr>
  <tr>
    <td>Ciaran McCreech</td>
    <td>University of Glasgow, United Kingdom</td>
  </tr>
  <tr>
    <td>Laurent Michel</td>
    <td>University of Connecticut, USA</td>
  </tr>
  <tr>
    <td>Nysret Musliu</td>
    <td>TU Wien, Austria</td>
  </tr>
  <tr>
    <td>Peter Nightingale</td>
    <td>University of York, United Kingdom</td>
  </tr>
  <tr>
    <td>Jakob Nordstrom</td>
    <td>University of Copenhagen & Lund University</td>
  </tr>
  <tr>
    <td>Marie Pelleau</td>
    <td>Université Côte d'Azur, France</td>
  </tr>
  <tr>
    <td>Guillaume Perez</td>
    <td>University of Nice-Sophia Antipolis/I3S, France</td>
  </tr>
  <tr>
    <td>Laurent Perron</td>
    <td>Google France, France</td>
  </tr>
  <tr>
    <td>Steve Prestwich</td>
    <td>Insight Centre for Data Analytics, Ireland</td>
  </tr>
  <tr>
    <td>Patrick Prosser</td>
    <td>Glasgow University, United Kingdom</td>
  </tr>
  <tr>
    <td>Charles Prud'Homme</td>
    <td>IMT Atlantique, LS2N, France</td>
  </tr>
  <tr>
    <td>Claude-Guy Quimper</td>
    <td>Laval University, Canada</td>
  </tr>
  <tr>
    <td>Jean-Charles Regin</td>
    <td>University Nice-Sophia Antipolis/I3S/CNRS, France</td>
  </tr>
  <tr>
    <td>Pierre Schaus</td>
    <td>UCLouvain, Belgium</td>
  </tr>
  <tr>
    <td>Paul Shaw</td>
    <td>IBM, France</td>
  </tr>
  <tr>
    <td>Mohamed Siala</td>
    <td>INSA Toulouse & LAAS-CNRS, France</td>
  </tr>
  <tr>
    <td>Mate Soos</td>
    <td>National University of Singapore, Singapore</td>
  </tr>
  <tr>
    <td>Kostas Stergiou</td>
    <td>University of Western Macedonia, Greece</td>
  </tr>
  <tr>
    <td>Guido Tack</td>
    <td>Monash University, Australia</td>
  </tr>
  <tr>
    <td>Cyril Terrioux</td>
    <td>LIS - UMR CNRS 7020 - Aix-Marseille Université, France</td>
  </tr>
  <tr>
    <td>Gilles Trombettoni</td>
    <td>LIRMM, University of Montpellier, France</td>
  </tr>
  <tr>
    <td>Charlotte Truchet</td>
    <td>Université de Nantes, France</td>
  </tr>
  <tr>
    <td>Willem-Jan Van Hoeve</td>
    <td>Carnegie Mellon University, USA</td>
  </tr>
  <tr>
    <td>Hélène Verhaeghe</td>
    <td>KU Leuven, Belgium</td>
  </tr>
  <tr>
    <td>Petr Vilím</td>
    <td>ScheduleOpt</td>
  </tr>
  <tr>
    <td>Mark Wallace</td>
    <td>Monash University, Australia</td>
  </tr>
  <tr>
    <td>Ruiwei Wang</td>
    <td>National University of Singapore, Singapore</td>
  </tr>
  <tr>
    <td>Armin Wolf</td>
    <td>Fraunhofer Institute, Germany</td>
  </tr>
  <tr>
    <td>Lebbah Yahia</td>
    <td>University of Oran 1, Algeria</td>
  </tr>
  <tr>
    <td>Neng-Fa Zhou</td>
    <td>CUNY Brooklyn College & Graduate Center, USA</td>
  </tr>
</table>
Additional Reviewers

Josep Alòs

Valentin Antuori

Noureddine Aribi

Federico Baldo

Nassim Belmecheri

Hendrik Bierlee

Andrea Borghesi

Touati Chahira

Vianney Coppé

Timothy Curry

Toby O. Davies

Guillaume Derval

Julien Ferry

Romain Fontaine

Andrea Formisano

Matteo Francobaldi

Xavier Gillard

Ramiz Gindullin

Luca Giuliani

Arthur Gontier

Cristian Grozea

Matthias Horn

Victor Jung

Anthony Karahalios

Amina Kemmar

Lucas Kletzander

Stepan Kochemazov

Tanguy Lapègue

Shuolin Li

Matthew J. McIlree

Claude Michel

Florian Mischek

Eleonora Misino

Pierre Montalbano

Antonio Morgado

Bertrand Neveu

Andy Oertel

Jolan Philippe

Philippe Refalo

András Z. Salamon

Noah Schutte

Jean-Baptiste Sciau

Dimosthenis C. Tsouros

Marc Vinyals

Damien T. Wojtowicz

Oleg Zaikin

Yuanlin Zhang
<figure><img src="image_11.png" /></figure>
List of Authors

Younes Aalian (6)
Département de mathématiques et de génie industriel, Polytechnique Montréal, Québec, Canada

Vahid Eghbal Akhlaghi (46)
H. Milton Stewart School of Industrial and Systems Engineering, Georgia Institute of Technology, Atlanta, GA, USA

Mohammed Alswaitti (44)
University of Luxembourg, Luxembourg; Interdisciplinary Centre for Security, Reliability and Trust (SnT), Luxembourg

Nahum Alvarez (31)
Airbus (AI Research), Toulouse, France

Josep Alòs (7)
Logic & Optimization Group (LOG), University of Lleida, Spain

Carlos Ansótegui (7)
Logic & Optimization Group (LOG), University of Lleida, Spain

João Araújo (8)
Universidade Nova de Lisboa, Lisbon, Portugal

Christian Artigues (31)
LAAS-CNRS, Université de Toulouse, CNRS, Toulouse, France

Gilles Audemard (9)
CRIL, Univ. Artois & CNRS, France

James Bailey (42)
The University of Melbourne, Australia

Vincent Barbosa Vaz (42)
The University of Melbourne, Australia; the Australian Research Council OPTIMA ITTC, Melbourne, Australia

Aymeric Beauchamp (10)
University of Orléans, INSA Centre Val de Loire, LIFO EA 4022, France

J. Christopher Beck (17, 23)
Department of Mechanical and Industrial Engineering, University of Toronto, Canada

Belaid Benhamou (35)
LIS, Aix-Marseille University, Marseille, France

Senne Berden (36)
KU Leuven, Belgium

Jeremias Berg (18)
HIIT, Department of Computer Science, University of Helsinki, Finland

Timo Berthold (27)
Fair Isaac Deutschland GmbH, Berlin, Germany; TU Berlin, Germany

Sévérine Fetgo Betmbe (20)
Faculty of Sciences, Department of Mathematics and Computer Science, University of Dschang, Cameroon

Giovanni Lo Bianco (17)
Department of Mechanical and Industrial Engineering, University of Toronto, Canada

Ignace Bleukx (11)
DTAI, KU Leuven, Belgium

Bart Bogaerts (11)
Artificial Intelligence Lab, VUB, Brussels, Belgium

Wirattawut Boonbandansook (46)
H. Milton Stewart School of Industrial and Systems Engineering, Georgia Institute of Technology, Atlanta, GA, USA

Kyle E. C. Booth (43)
Amazon Quantum Solutions Lab, Seattle, WA, USA

Pascal Bouvry (44)
University of Luxembourg, Luxembourg; Interdisciplinary Centre for Security, Reliability and Trust (SnT), Luxembourg

Shaowei Cai (12, 41)
State Key Laboratory of Computer Science, Institute of Software, Chinese Academy of Sciences, Beijing, China

Quentin Cappart (25)
Polytechnique Montréal, Montreal, Canada

Anouck Chan (32)
ONERA, ONERA DTIS, Toulouse, Université de Toulouse, France

Choiwah Chow (8)
Universidade Aberta, Lisbon, Portugal

Yi Chu (12)
Institute of Software, Chinese Academy of Sciences, Beijing, China
Eldan Cohen (33)
Department of Mechanical and Industrial
Engineering, University of Toronto, Canada

Manuel Combarro Simón (44)
University of Luxembourg, Luxembourg;
Interdisciplinary Centre for Security, Reliability
and Trust (SnT), Luxembourg

Vianney Coppé (13)
UCLouvain, Louvain-la-Neuve, Belgium

Kevin Dalmeijer (46)
H. Milton Stewart School of Industrial and
Systems Engineering, Georgia Institute of
Technology, Atlanta, GA, USA

Grégoire Danoy (44)
University of Luxembourg, Luxembourg;
Interdisciplinary Centre for Security, Reliability
and Trust (SnT), Luxembourg

Thi-Bich-Hanh Dao (10)
University of Orléans, INSA Centre Val de Loire,
LIFO EA 4022, France

Frits de Nijs (21)
Department of Data Science and AI, Faculty of
IT, Monash University, Clayton, VIC, Australia;
ARC Industrial Training and Transformation
Centre OPTIMA, Carlton, VIC, Australia

Jo Devriendt (11)
DTAI, KU Leuven, Belgium

Arnaud Deza (14)
Department of Mechanical and Industrial
Engineering, University of Toronto, Canada

Frédéric Didier (3)
Google, Paris, France

Clémentin Tayou Djamegni (20)
Faculty of Sciences, Department of Mathematics
and Computer Science, University of Dschang,
Cameroon; UIT Fotso Victor of Bandjoun,
Department of Computer Engineering,
University of Dschang, Cameroon

Alexandre Dubray (15)
Institute of Information and Communication
Technologies, Electonics and Applied
Mathematics (ICTEAM), UCLouvain, Belgium

Yousra El Ghazi (16)
Aix Marseille Univ, Université de Toulon, CNRS,
LIS, Marseille, France

Tristan François (25)
Ecole Polytechnique, Palaiseau, France

Michel Gamache (6)
Département de mathématiques et de génie
industriel, Polytechnique Montréal, Québec,
Canada

Emilio Gamba (11)
Data Analytics Lab, VUB, Brussels, Belgium;
DTAI, KU Leuven, Belgium

Maria Garcia de la Banda (1, 21, 37)
Department of Data Science and AI, Faculty of
IT, Monash University, Clayton, Victoria,
Australia; ARC Industrial Training and
Transformation Centre OPTIMA, Carlton,
Victoria, Australia

Louis Gautier (25)
Ecole Polytechnique, Palaiseau, France

Steven Gay (3)
Google, Paris, France

Xavier Gillard (13, 45)
UCLouvain, Louvain-la-Neuve, Belgium

Ambros Gleixner (27)
HTW Berlin, Germany; Zuse Institute Berlin,
Germany

Gaël Glorian (29)
Huawei Technologies Ltd, CSI Paris,
Boulogne-Billancourt, France

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

Arnoosh Golestanian (17)
Department of Mechanical and Industrial
Engineering, University of Toronto, Canada

Carla P. Gomes (47)
Cornell University, Ithaca, NY, USA

Laura Greenstreet (47)
Cornell University, Ithaca, NY, USA

Marc Grimson (47)
Cornell University, Ithaca, NY, USA

Tias Guns (11, 36)
DTAI, KU Leuven, Belgium; Data Analytics
Lab, VUB, Brussels, Belgium

Djamal Habet (16)
Aix Marseille Univ, Université de Toulon, CNRS,
LIS, Marseille, France

Emmanuel Hebrard (19)
LAAS-CNRS, Université de Toulouse, France
Laurent Houssin ID (19)
ISAE-SUPAERO, Université de Toulouse,
France

Tingting Hu (34)
University of Luxembourg, Luxembourg

Alexey Ignatiev ID (38)
Department of Data Science and AI, Monash
University, Clayton, Australia

Hannes Ihalainen ID (18)
HIIT, Department of Computer Science,
University of Helsinki, Finland

Peter J. Stuckey ID (42)
Monash University, Australia; the Australian
Research Council OPTIMA ITTC, Melbourne,
Australia

Christoph Jabs ID (18)
HIIT, Department of Computer Science,
University of Helsinki, Finland

Mikoláš Janota ID (8)
Czech Technical University in Prague, Czech
Republic

Carla Juvin (19)
LAAS-CNRS, Université de Toulouse, France

Matti Järvisalo ID (18)
HIIT, Department of Computer Science,
University of Helsinki, Finland

Roger Kameugne ID (20)
Faculty of Sciences, Department of Mathematics
and Computer Science, University of Maroua,
Cameroon

Elias B. Khalil (14)
Department of Mechanical and Industrial
Engineering, University of Toronto, Canada

Matthias Klapperstueck ID (21)
Department of Human-Centred Computing,
Faculty of IT, Monash University, Clayton, VIC,
Australia

T. K. Satish Kumar (40)
University of Southern California, Los Angeles,
CA, USA

Ryo Kuroiwa ID (23)
Department of Mechanical and Industrial
Engineering, University of Toronto, Canada

Petr Kučera ID (22)
Department of Theoretical Computer Science
and Mathematical Logic, Faculty of
Mathematics and Physics, Charles University,
Prague, Czech Republic

Arnaud Lallouet ID (29)
Huawei Technologies Ltd, CSI Paris,
Boulogne-Billancourt, France

Christopher Leckie ID (42)
The University of Melbourne, Australia

Christophe Lecoutre (9)
CRIL, Univ. Artois & CNRS, France

Chungjae Lee ID (46)
H. Milton Stewart School of Industrial and
Systems Engineering, Georgia Institute of
Technology, Atlanta, GA, USA

Jimmy H.M. Lee (2)
Department of Computer Science and
Engineering, The Chinese University of Hong
Kong, China

Jack Lee-Kopij (21)
Woodside Energy Ltd., Perth, Australia

Zhendong Lei ID (12)
State Key Laboratory of Computer Science,
Institute of Software, Chinese Academy of
Sciences, Beijing, China

Ang Li (40)
University of Southern California, Los Angeles,
CA, USA

Chang Liu (14)
Department of Mechanical and Industrial
Engineering, University of Toronto, Canada

Pierre Lopez ID (19)
LAAS-CNRS, Université de Toulouse, France

Samir Loudni ID (10)
TASC (LS2N-CNRS), IMT Atlantique, France,
GREYC, University of Caen Normandy, France

Chuan Luo ID (12)
School of Software, Beihang University, Beijing,
China

Steve Malalel (24)
Université Côte d'Azur, CNRS, I3S, Nice,
France

Arnaud Malapert (24)
Université Côte d'Azur, CNRS, I3S, Nice,
France
Tom Marty (25)
Polytechnique Montréal, Montreal, Canada;
Ecole Polytechnique, Palaiseau, France

Ciaran McCreesh (26)
University of Glasgow, UK

Sheila McIlraith (33)
Department of Computer Science, University of
Toronto, Canada; Vector Institute, Toronto,
Canada

Matthew J. McIlree (26)
University of Glasgow, UK

Gioni Mexi (27)
Zuse Institute Berlin, Germany

Renee Mirka (47)
Cornell University, Ithaca, NY, USA

Jedrzej Musial (44)
Poznan University of Technology, Poland

Sibylle Möhle (49)
Max Planck Institute for Informatics,
Saarbrücken, Germany

Nicolas Navet (34)
University of Luxembourg, Luxembourg

Siegfried Nijssen (15, 45)
Institute of Information and Communication
Technologies, Electonics and Applied
Mathematics (ICTEAM), UCLouvain, Belgium

Jakob Nordström (27)
University of Copenhagen, Denmark; Lund
University, Sweden

Thierry Noulamo (20)
UIT Fotso Victor of Bandjoun, Department of
Computer Engineering, University of Dschang,
Cameroon

Marie Pelleau (24)
Université Côte d'Azur, CNRS, I3S, Nice,
France

Cong Peng (12)
Finovation in CCBFT, Beijing, China

Xiao Peng (28)
Univ Lyon, INSA Lyon, Inria, CITI, EA3720,
69621 Villeurbanne, France

Guillaume Perez (29)
Huawei Technologies Ltd, CSI Paris,
Boulogne-Billancourt, France

Laurent Perron (3)
Google, Paris, France

Vaidyanathan Peruvemba Ramaswamy (48)
Algorithms and Complexity Group, TU Wien,
Austria

Gilles Pesant (6)
Département de génie informatique et génie
logiciel, Polytechnique Montréal, Québec,
Canada

Andreas Plank (49)
Institute for Symbolic Artificial Intelligence,
JKU Linz, Austria

Nikolaos Ploskas (30)
University of Western Macedonia, Kozani,
Greece

Thomas Polacsek (32)
ONERA, ONERA DTIS, Toulouse, Université
de Toulouse, France

Guillaume Povéda (31)
Airbus (AI Research), Toulouse, France

Charles Prud'homme (9)
TASC, IMT-Atlantique, LS2N-CNRS, France

Louis-Martin Rousseau (25)
Polytechnique Montréal, Montreal, Canada

Stéphanie Roussel (32)
ONERA, ONERA DTIS, Toulouse, Université
de Toulouse, France

Jean-Charles Régin (24)
Université Côte d'Azur, CNRS, I3S, Nice,
France

Josep M. Salvia (7)
Logic & Optimization Group (LOG), University
of Lleida, Spain

Pierre Schaus (13, 15, 45)
UCLouvain, Louvain-la-Neuve, Belgium

Thomas Schiex (4)
Université Fédérale de Toulouse, ANITI,
INRAE, UR 875, 31326 Toulouse, France

Martina Seidl (49)
Institute for Symbolic Artificial Intelligence,
JKU Linz, Austria

Ilankaikone Senthooran (21)
Department of Data Science and AI, Faculty of
IT, Monash University, Clayton, VIC, Australia;
ARC Industrial Training and Transformation
Centre OPTIMA, Carlton, VIC, Australia
Pouya Shati (33)
Department of Computer Science, University of
Toronto, Canada; Vector Institute, Toronto,
Canada

Sylvain Soliman ID (35)
Lifeware team, Inria Saclay, Palaiseau, France

Christine Solnon (28)
Univ Lyon, INSA Lyon, Inria, CITI, EA3720,
69621 Villeurbanne, France

Kostas Stergiou ID (30)
University of Western Macedonia, Kozani,
Greece

Peter J. Stuckey ID (38)
Department of Data Science and AI, Monash
University, Clayton, Australia; Australian
Research Council OPTIMA ITTC, Clayton,
Australia

Wijnand Suijlen ID (29)
Huawei Technologies Ltd, CSI Paris,
Boulogne-Billancourt, France

Stefan Szeider ID (39, 48)
Algorithms and Complexity Group, TU Wien,
Austria

Guido Tack ID (37)
Department of Data Science and AI, Faculty of
IT, Monash University, Clayton, Australia; ARC
Industrial Training and Transformation Centre
OPTIMA, Clayton, Australia

Pierre Talbot ID (34, 44)
University of Luxembourg, Luxembourg;
Interdisciplinary Centre for Security, Reliability
and Trust (SnT), Luxembourg

Chengyu Tao ID (17)
Department of Mechanical and Industrial
Engineering, University of Toronto, Canada

Cyril Terrioux ID (16)
Aix Marseille Univ, Université de Toulon, CNRS,
LIS, Marseille, France

Pierre Tessier (25)
Ecole Polytechnique, Palaiseau, France

Eduard Torres ID (7)
Logic & Optimization Group (LOG), University
of Lleida, Spain

Van-Giang Trinh ID (35)
LIS, Aix-Marseille University, Marseille, France

Dimosthenis C. Tsouros ID (30, 36)
KU Leuven, Belgium

Pashootan Vaezipoor (14)
Department of Computer Science, University of
Toronto, Canada

Pascal Van Hentenryck ID (46)
H. Milton Stewart School of Industrial and
Systems Engineering, Georgia Institute of
Technology, Atlanta, GA, USA

Petr Vilím ID (5)
ScheduleOpt, Nový Knín, Czech Republic

Christel Vrain ID (10)
University of Orléans, INSA Centre Val de Loire,
LIFO EA 4022, France

Shimao Wang (41)
School of Information Science and Technology,
Northeast Normal University, Changchun, China

Xinyu Wang (41)
School of Information Science and Technology,
Northeast Normal University, Changchun, China

Yiyuan Wang ID (41)
School of Information Science and Technology,
Northeast Normal University, Changchun,
China; Key Laboratory of Applied Statistics of
MOE, Northeast Normal University, Changchun,
China

Sameela Suharshani Wijesundara ID (37)
Department of Data Science and AI, Faculty of
IT, Monash University, Clayton, Australia; ARC
Industrial Training and Transformation Centre
OPTIMA, Clayton, Australia

Michael Wybrow ID (21)
Department of Human-Centred Computing,
Faculty of IT, Monash University, Clayton, VIC,
Australia

Minghao Yin ID (41)
School of Information Science and Technology,
Northeast Normal University, Changchun,
China; Key Laboratory of Applied Statistics of
MOE, Northeast Normal University, Changchun,
China

Jinqiang Yu ID (38)
Department of Data Science and AI, Monash
University, Clayton, Australia; Australian
Research Council OPTIMA ITTC, Clayton,
Australia

Han Zhang (40)
University of Southern California, Los Angeles,
CA, USA
0:xx

Authors

Tianwei Zhang (39)

Algorithms and Complexity Group, TU Wien,
Austria

Yujiao Zhao (41)

School of Information Science and Technology,
Northeast Normal University, Changchun, China

Kexin Zheng (40)

University of Southern California, Los Angeles,
CA, USA

Wenbo Zhou (41)

School of Information Science and Technology,
Northeast Normal University, Changchun,
China; Key Laboratory of Applied Statistics of
MOE, Northeast Normal University, Changchun,
China; Key Laboratory of Symbolic
Computation and Knowledge Engineering of
MOE, Jilin University, Changchun, China
# Beyond Optimal Solutions for Real-World Problems

Maria Garcia de la Banda

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

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

## Abstract

Combinatorial optimisation technology has come a long way. We now have mature high-level modelling languages in which to specify a model of the particular problem of interest [18, 7, 24, 6]; robust complete solvers in each major constraint paradigm, including Constraint Programming (CP) [1, 19], MaxSAT [5, 11], and Mixed Integer Programming (MIP) [2, 3]; effective incomplete search techniques that can easily be combined with complete solvers to speed up the search such as Large Neighbourhood Search [23]; and enough general knowledge about modelling techniques to understand the need for our models to incorporate components such as global constraints [25], symmetry constraints [8], and more. All this has significantly reduced the amount of knowledge required to apply this technology successfully to the many different combinatorial optimisation problems that permeate our society.

And yet, not many organisations use such advanced optimisation technology; instead, they often rely on the solutions provided by problem-specific algorithms that are implemented in traditional imperative languages and lack any of the above advances. Further, while advanced optimisation technology is particularly suitable for the kind of complex human-in-the-loop decision-making problems that occur in critical sectors of our society, including health, transport, energy, disaster management, environment and finance, these decisions are often still made by people with little or no technological support. In this extended abstract I argue that to change this state of affairs, our research focus needs to change from improving the technology on its own, to improving it so that users can better trust, use, and maintain the optimisation systems that we develop with it. The rest of this extended abstract discusses my personal experiences and opinion on these three points.

## Trust

I highlight trust (which focuses on the user's point of view) rather than trustworthiness (which is a characteristic of the software itself) because I think it is the former rather than the latter that is at stake for the adoption of optimisation technology.

One of the biggest hurdles I have found for trust in the context of optimisation systems is for the domain experts to (feel like they) *understand the underlying model*. While many users will never do (or have to), I believe it is key for domain experts to have a high-level understanding of the constraints in the model, since their (dis)trust will likely spread through the organisation, impacting the adoption of the system. Thanks to the use of high-level modelling languages in CP, our group has achieved this [13] by documenting the constraints in a language the user knows (mathematics) and linking each constraint to the particular part of the model that implements it (via comments). While domain experts do not completely understand the model, the similarity between the format they understand (mathematics) and the model constraint has helped them verify our perception of their problem and improved their trust in the model. However, more needs to be done in this direction via the development of formal techniques. For example, our group is exploring the use of *domain-specific languages* [10] as a bridge between domain experts and modellers that helps both trust and maintenance (see later). This [27] and other approaches need to be explored.

A very significant source of trust for our domain experts (and of trustworthiness for the software) has been the development of two different models implemented by two different people for the same problem [13]. While this can be seen as a prohibitively expensive exercise, it did not take that long once the first model was mature, is a good way to onboard new optimisation team members, and has helped up detect not only bugs but also differences in the interpretation of domain expert information. For optimisation problems where it is not possible to verify the optimality (or even correctness) of the solution, we see such *redundant modelling* as the only solution for now. Interestingly, a significant

<figure><img src="image_13.png" /></figure>
1:2 Beyond Optimal Solutions

step forward in obtaining the trust of our domain experts has been the *generation of an optimality gap* whenever an optimal solution could not be found due to time constraints. While explaining this concept took time, once understood it has boosted their trust, particularly when tackling problems where the solution is not easy verifiable or when approximated models/data are used (needed for speed, see later). This makes it difficult to work with CP and SAT solvers, as they usually lack tight lower bounds. Finally, trust is often developed through the use of the system, which I discuss below.

Use

Usability is known to be key for the deployment of software systems. By “system” in our context, I refer to the combination of the problem model(s), the associated solver(s) and, importantly, the User Interface (UI) that often integrates them and is fundamental to their success. In addition to the traditional usability characteristics of software systems, I believe an optimisation system requires particular care in the following areas. *Interaction*, i.e., the system must allow users to interact with the UI not only to provide and modify the input data, but also to modify the constraints (at the very least by turning some on/off) as well as explore and compare solutions, as argued in [17, 15]. Incremental compilers and solvers would significantly help in making this easier, as well as generic ways for the UIs to communicate with them. *Conflict resolution*, that is, ensuring the system can not only detect infeasible instances, but also support users in understanding the data/constraints that cause infeasibility and how to modify the instance to make it feasible. Any interactive optimisation system that has users, will likely have conflicts. Thus, it is mandatory for CP to improve its conflict resolution technology which, while existent [16, 14, 22], is not widespread and it is often still problem-dependent, overwhelming (in the number of constraints shown to the user) and slow. Without it, users will be “stumped” when (rather than if) infeasibility is reached. *Solution diversity*, that is, supporting users in obtaining a diverse set of (close-to-optimal) solutions, where diversity is measured by a user-provided metric modelled somehow. While some solver-independent technology has been developed and implemented for this [9, 20, 12], it should be easier to use and more widespread. Further, it requires sophisticated solution comparison capabilities and, importantly, for optimal solutions to be found in seconds rather than hours. This brings me to speed, an area where CP solvers are falling behind. Most of our research group applications now use MIP solvers due to the need for floats (which precludes us from using learning solvers such as Chuffed [4]), but also to the lack of effective warm-start processes that are available in MIP solvers. Interestingly, data and model approximations have been proved to achieve orders of magnitude speedups with small reductions in optimality [13]. Developing generic (i.e., problem independent) accurate approximations would be extremely useful for complex decision systems. Other areas where I think generic CP methods are worth investigating more include dealing with *uncertainty* and *online* problems, ensuring solution *fairness* (even if it is over time), and studying *predict + optimise* approaches.

Maintain

I know very few papers devoted to the issue of maintenance in optimisation technology. While this
may be due to my lack of knowledge, I suspect it is also due to the limited adoption of optimisation
technology. While the issues in this area are again common to other software systems, I believe the
solutions for CP require special attention. For example, the issue of changes in user requirements
(that our research group calls *problem drift*) seems particularly prevalent in decision-making systems,
as such problems can evolve rapidly due to unforeseen circumstances. This can make optimisation
systems obsolete faster than expected. Our research group has proposed to tackle problem drift
by developing a *requirements model* implemented in the above-mentioned MDSLs and created
by both domain experts and modellers that, when modified re-generates parts of the model to
support the modifications [27]. This and other approaches such as the creation of reusable models
components [21, 26], or instantiatable classes for common problem domains, are worth investigating.
**2012 ACM Subject Classification** Theory of computation → Constraint and logic programming; Theory of computation → Integer programming; Human-centered computing → Information visual-ization

**Keywords and phrases** Combinatorial optimisation systems, usability, trust, maintenance

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

**Category** Invited Talk

**Acknowledgements** Any good idea here is the result of working in many different projects with my colleagues at the Optimisation and at the Data Visualisation and Immersive Analytics research groups in the Faculty of IT, Monash University. I have learned a lot from them.

References

1 GECODE - An open, free, efficient constraint solving toolkit. URL: https://www.gecode.org/.

2 Gurobi Optimizer Reference Manual. URL: https://www.gurobi.com/documentation/9.5/refman/index.html.

3 IBM ILOG CPLEX optimizer. URL: https://www.ibm.com/products/ilog-cplex-optimization-studio/cplex-optimizer.

4 Geoffrey Chu. Improving combinatorial optimization - extended abstract. In Francesca Rossi, editor, *23rd International Joint Conference on Artificial Intelligence, IJCAI 2013*, pages 3116–3120. IJCAI/AAAI, 2013. URL: http://www.aaai.org/ocs/index.php/IJCAI/IJCAI13/paper/view/6687.

5 Jessica Davies and Fahemi Bacchus. Solving MAXSAT by solving a sequence of simpler SAT instances. In Jimmy Ho-Man Lee, editor, *17th International Conference on Principles and Practice of Constraint Programming - CP 2011*, Lecture Notes in Computer Science, pages 225–239. Springer, 2011. doi:10.1007/978-3-642-23786-7_19.

6 Robert Fourer, David M Gay, and Brian W Kernighan. A modeling language for mathematical programming. *Management Science*, 36(5):519–554, 1990.

7 Alan M. Frisch, Warwick Harvey, Chris Jefferson, Bernadette Martínez-Hernández, and Ian Miguel. Essence: A constraint language for specifying combinatorial problems. *Constraints*, 13(3):268–306, September 2008. doi:10.1007/s10601-008-9047-y.

8 Ian P. Gent, Karen E. Petrie, and Jean-François Puget. Symmetry in constraint programming. In Francesca Rossi, Peter van Beek, and Toby Walsh, editors, *Handbook of Constraint Programming*, volume 2 of *Foundations of Artificial Intelligence*, pages 329–376. Elsevier, 2006. doi:10.1016/S1574-6526(06)80014-3.

9 Emmanuel Hebrard, Brahim Hnich, Barry O'Sullivan, and Toby Walsh. Finding diverse and similar solutions in constraint programming. In Manuela M. Veloso and Subbarao Kambhampati, editors, *20th National Conference on Artificial Intelligence, AAAI 2005*, pages 372–377. AAAI Press / The MIT Press, 2005. URL: http://www.aaai.org/Library/AAAI/2005/aaai05-059.php.

10 Paul Hudak. Domain-specific languages. *Handbook of programming languages*, 3(39-60):21, 1997.

11 Alexey Ignatiev, António Morgado, and João Marques-Silva. RC2: an efficient maxsat solver. *J. Satisf. Boolean Model. Comput.*, 11(1):53–64, 2019. doi:10.3233/SAT190116.

12 Linnea Ingmar, Maria Garcia de la Banda, Peter J. Stuckey, and Guido Tack. Modelling diversity of solutions. In *34th Conference on Artificial Intelligence, AAAI 2020*, pages 1528–1535. AAAI Press, 2020. URL: https://ojs.aaai.org/index.php/AAAI/article/view/5512.

13 Matthias Klapperstueck, Frits de Nijs, Ilankaikone Senthooran, Jack Lee-Kopij, Maria Garcia de la Banda, and Michael Wybrow. Exploring hydrogen supply/demand networks: Modeller and domain expert views. In Roland Yap, editor, *29th International Conference on Principles and Practice of Constraint Programming - CP23*, LIPIcs. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, to appear, 2023.
1:4 Beyond Optimal Solutions

14 Niklas Lauffer and Ufuk Topcu. Human-understandable explanations of infeasibility for resource-constrained scheduling problems. In *ICAPS 2019 Workshop XAIP*, 2019.

15 Jie Liu, Tim Dwyer, Guido Tack, Samuel Gratzl, and Kim Marriott. Supporting the problem-solving loop: Designing highly interactive optimisation systems. *IEEE Trans. Vis. Comput. Graph.*, 27(2):1764–1774, 2021. doi:10.1109/TVCG.2020.3030364.

16 João Marques-Silva and Alessandro Previti. On computing preferred muses and mcses. In *Carsten Sinz and Uwe Egly, editors, 17th International Conference on Theory and Applications of Satisfiability Testing - SAT 2014*, Lecture Notes in Computer Science, pages 58–74. Springer, 2014. doi:10.1007/978-3-319-09284-3_6.

17 David Meignan, Sigrid Knust, Jean-Marc Frayret, Gilles Pesant, and Nicolas Gaud. A review and taxonomy of interactive optimization methods in operations research. *ACM Trans. Interact. Intell. Syst.*, 5(3):17:1–17:43, 2015. doi:10.1145/2808234.

18 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, *13th International Conference on Principles and Practice of Constraint Programming - CP 2007*, Lecture Notes in Computer Science, pages 529–543. Springer, 2007. doi:10.1007/978-3-540-74970-7_38.

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

20 Thierry Petit and Andrew C. Trapp. Finding diverse solutions of high quality to constraint optimization problems. In Qiang Yang and Michael J. Wooldridge, editors, *24th International Joint Conference on Artificial Intelligence, IJCAI 2015*, pages 260–267. AAAI Press, 2015. URL: http://ijcai.org/Abstract/15/043.

21 Sophia Saller and Jana Koehler. Easy, adaptable and high-quality modelling with domain-specific constraint patterns. *CoRR*, abs/2206.02479, 2022. doi:10.48550/arXiv.2206.02479.

22 Ilankaikone Senthooran, Matthias Klapperstueck, Gleb Belov, Tobias Czauderna, Kevin Leo, Mark Wallace, Michael Wybrow, and Maria Garcia de la Banda. Human-centred feasibility restoration in practice. *Constraints*, to appear, 2023.

23 Paul Shaw. Using constraint programming and local search methods to solve vehicle routing problems. In Michael J. Maher and Jean-François Puget, editors, *4th International Conference on Principles and Practice of Constraint Programming - CP98*, Lecture Notes in Computer Science, pages 417–431. Springer, 1998. doi:10.1007/3-540-49481-2_30.

24 Pascal Van Hentenryck. *The OPL optimization programming language*. MIT Press, Cambridge, MA, USA, 1999.

25 Willem-Jan van Hoeve and Irit Katriel. Global constraints. In Francesca Rossi, Peter van Beek, and Toby Walsh, editors, *Handbook of Constraint Programming*, volume 2 of *Foundations of Artificial Intelligence*, pages 169–208. Elsevier, 2006. doi:10.1016/S1574-6526(06)80010-6.

26 Toby Walsh. Constraint patterns. In Francesca Rossi, editor, *9th International Conference on Principles and Practice of Constraint Programming - CP 2003*, Lecture Notes in Computer Science, pages 53–64. Springer, 2003. doi:10.1007/978-3-540-45193-8_4.

27 Sameela Suharshani Wijesundara, Maria Garcia de la Banda, and Guido Tack. Addressing problem drift in UNHCR fund allocation. In Roland Yap, editor, *29th International Conference on Principles and Practice of Constraint Programming - CP23*, LIPIcs. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, to appear, 2023.
# A Tale of Two Cities: Teaching CP with Story-Telling

Jimmy H.M. Lee

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

Department of Computer Science and Engineering, The Chinese University of Hong Kong, China

## Abstract

This presentation is all about story-telling. It tells the story, the pedagogical innovations and experience of the co-development of three MOOCs on the subject of “Modeling and Solving Discrete Optimization Problems” by The Chinese University of Hong Kong (CUHK) and the University of Melbourne, each with unique culture and tradition. The MOOCs feature the Fable-based Learning approach, which is a form of problem-based learning encapsulated in a story plot. Each MOOC video begins with an animation that follows a story adapted from a Chinese classic. The heroes of the story encounter various optimization problems requiring technical assistance from two professors from modern time via a magical tablet granted to the heroes by a genie old man. The animation thus sets the stage for lecturing modeling and solving techniques. The new pedagogy provides a movie-like immersive experience to the learners, and aims at increasing learners’ motivation and interests as well as situating them in a coherent learning context. In addition to scriptwriting, animation production and embedding the teaching materials in the story plot, another challenge of the project is the remote distance between the two institutions as well as the need to produce all teaching materials in both (Mandarin) Chinese and English to cater for different geographic learning needs. The project and production spanned across 2016 and 2017. The MOOCs have been running recurrently on Coursera since January, 2017. We present learner statistics and feedback, and discuss our experience and preliminary observations of adopting the online materials in a Flipped Classroom setting at CUHK.

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

**Keywords and phrases** Constraint Programming, MOOCs, Fable-based Learning

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

**Category** Invited Talk

<figure><img src="image_15.png" /></figure>
<figure><img src="image_16.png" /></figure>
# The CP-SAT-LP Solver

Laurent Perron

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

Google, Paris, France

Frédéric Didier

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

Google, Paris, France

Steven Gay

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

Google, Paris, France

## Abstract

The CP-SAT-LP solver is developed by the Operations Research team at Google and is part of the OR-Tools [8] open-source optimization suite. It is an implementation of a purely integral Constraint Programming solver on top of a SAT solver using Lazy Clause Generation [11]. It draws its inspiration from the chuffed solver [4], and from the CP 2013 plenary by Peter Stuckey on Lazy Clause Generation [12].

The CP-SAT-LP solver improves upon the chuffed solver [4] in two main directions. First, it uses a simplex alongside the SAT engine. Second, it implements and relies upon a portfolio of diverse workers for its search part.

The use of the simplex brings the obvious advantages of a linear relaxation on the linear part of the full model. It also started the integration of MIP technology into CP-SAT-LP. This is a huge endeavour, as MIP solvers are mature and complex. It includes presolve – which was already a part of CP-SAT –, dual reductions, specific branching rules, cuts, reduced cost fixing, and more advanced techniques. It also allows to integrate tightly the research from the Scheduling on MIP community [3, 1, 9] along with the most advanced scheduling algorithms [13]. This has enabled breakthroughs in solving and proving hard scheduling instances of the Job-Shop problems [5] and Resource Constraint Project Scheduling Problems [6, 2].

Using a portfolio of different workers makes it easier to try new ideas and to incorporate orthogonal techniques with little complication, except controlling the explosion of potential workers. These workers can be categorized along multiple criteria like finding primal solutions – either using complete solvers, Local Search [7] or Large Neighborhood Search [10] –, improving dual bounds, trying to reduce the problem with the help of continuous probing. This diversity of behaviors has increased the robustness of the solver, while the continuous sharing of information between workers has produced massive speedups when running multiple workers in parallel.

All in all, CP-SAT-LP is a state-of-the-art solver, with unsurpassed performance in the Constraint Programming community, breakthrough results on Scheduling benchmarks (with the closure of many open problems), and competitive results with the best MIP solvers (on purely integral problems).

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

**Keywords and phrases** Constraint Programming, Operations Research, Sat Solver

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

**Category** Invited Talk

## References

1. David Applegate and William Cook. A computational study of the job-shop scheduling problem. *ORSA Journal on Computing*, 3(2):149–156, 1991. doi:10.1287/ijoc.3.2.149.

2. Christian Artigues, Sophie Demassey, and Emmanuel Neron. *Resource-Constrained Project Scheduling: Models, Algorithms, Extensions and Applications*. ISTE/Wiley, 2008. URL: https://hal.science/hal-00482946.

3. Egon Balas. Disjunctive programming and a hierarchy of relaxations for discrete optimization problems. *SIAM Journal on Algebraic Discrete Methods*, 6(3):466–486, 1985.

<figure><img src="image_20.png" /></figure>
3:2 The CP-SAT-LP Solver

4 Geoffrey Chu, Peter J. Stuckey, Andreas Schutt, Thorsten Ehlers, Graeme Gange, and Kathryn Francis. the chuffed solver, June 2023. URL: https://github.com/chuffed/chuffed.

5 Junwen Ding, Zhipeng Lü, Chu-Min Li, Liji Shen, Liping Xu, and Fred Glover. A two-individual based evolutionary algorithm for the flexible job shop scheduling problem. In *Proceedings of the AAAI Conference on Artificial Intelligence*, volume 33, pages 2262–2271, 2019.

6 Rainer Kolisch and Arno Sprecher. Psplib - a project scheduling problem library: Or software - orsep operations research software exchange program. *European Journal of Operational Research*, 96(1):205–216, 1997. doi:10.1016/S0377-2217(96)00170-1.

7 Bjørnar Luteberget and Giorgio Sartor. Feasibility jump: an lp-free lagrangian mip heuristic. *Mathematical Programming Computation*, 15, March 2023. doi:10.1007/s12532-023-00234-8.

8 Laurent Perron and Vincent Furnon. Or-tools, March 2023. URL: https://developers.google.com/optimization/.

9 Maurice Queyranne. Structure of a simple scheduling polyhedron. *Math. Program.*, 58:263–285, 1993. doi:10.1007/BF01581271.

10 Paul Shaw. Using constraint programming and local search methods to solve vehicle routing problems. In Michael J. Maher and Jean-François Puget, editors, *Principles and Practice of Constraint Programming - CP98, 4th International Conference, Pisa, Italy, October 26-30, 1998, Proceedings*, volume 1520 of *Lecture Notes in Computer Science*, pages 417–431. Springer, 1998. doi:10.1007/3-540-49481-2_30.

11 Peter J. Stuckey. Lazy clause generation: Combining the power of sat and cp (and mip?) solving. In Andrea Lodi, Michela Milano, and Paolo Toth, editors, *Integration of AI and OR Techniques in Constraint Programming for Combinatorial Optimization Problems*, pages 5–9, Berlin, Heidelberg, 2010. Springer Berlin Heidelberg.

12 Peter J. Stuckey. Those who cannot remember the past are condemned to repeat it. In Christian Schulte, editor, *Principles and Practice of Constraint Programming*, pages 5–6, Berlin, Heidelberg, 2013. Springer Berlin Heidelberg.

13 Petr Vilím. Timetable edge finding filtering algorithm for discrete cumulative resources. In Tobias Achterberg and J. Christopher Beck, editors, *Integration of AI and OR Techniques in Constraint Programming for Combinatorial Optimization Problems*, pages 230–245, Berlin, Heidelberg, 2011. Springer Berlin Heidelberg.
# Coupling CP with Deep Learning for Molecular Design and SARS-CoV2 Variants Exploration

Thomas Schiex

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

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

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

Université Fédérale de Toulouse, ANITI, INRAE, UR 875, 31326 Toulouse, France

## Abstract

The use of discrete optimization, including Constraint Programming, for designing objects that we completely understand is quite usual. In this talk, I'll show how designing specific biomolecules (proteins) raises new challenges, requiring solving problems that combine precise design targets, approximate laws, and design rules that can be deep-learned from data.

**2012 ACM Subject Classification** Computing methodologies → Artificial intelligence; Computing methodologies → Machine learning; Theory of computation → Constraint and logic programming; Computing methodologies → Learning in probabilistic graphical models

**Keywords and phrases** graphical models, deep learning, constraint programming, cost function networks, random Markov fields, decision-focused learning, protein design

Digital Object Identifier 10.4230/LIPIcs.CP.2023.4

Category Invited Talk

## 1 Introduction

Proteins are biomolecules that support most mechanisms in living organisms, from viruses to human beings. They already have major commercial applications in green chemistry (as enzymes) but also in the health domain (e.g., the anti-CoViD Regeneron™ antibodies are proteins). Most commercially used proteins are either natural proteins or engineered versions of natural proteins. To go beyond the repertoire of natural proteins, it is important to be able to reliably and efficiently design new proteins, with new capacities [9]. Proteins are defined by their amino acid sequence, a discrete object defined over an alphabet of 20 characters. Once the sequence of a protein is fixed, it can be encoded into a suitable microbe, enabling the cheap manufacturing of these complex microscopic assemblies.

Optimization is often used to design objects such as schedules, assignments, time-tables or packing, which we completely understand. Instead, proteins are tiny physical objects that live in the realm of quantum physics. Their behavior is hard to formally, precisely and concisely capture. Designing new proteins therefore requires to combine knowledge, expressed as approximate laws of physics, with targeted design constraints and criteria, in the context of large sets of data of past successful designs (natural proteins) that also embody the many hidden laws which a successfully expressed protein must satisfy.

## 2 Designing proteins and SARS-CoV2 variants with CP

In this talk, we will see how Cost Function Networks (CFNs), a weighted variant of Constraint Networks/CP) can help us design new proteins [1, 6]. Alone, CFNs can already capture both logical information (constraints) and numerical information, enabling the simultaneous representation of approximate laws of physics and design targets. By solving suitable instances of Weighted Constraint Satisfaction Problems, one can already produce protein sequences that can be tested in silico (with e.g., AlphaFold2 [7]), characterized experimentally, and lead to successful designs [8].

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

© Thomas Schiex; 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. 4; pp. 4:1–4:3
Leibniz International Proceedings in Informatics

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

LIPICS Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl Publishing, Germany
## 4:2 Coupling CP with Deep Learning for Molecular Design

By leveraging the exhaustive enumeration capabilities of exact discrete solvers, it becomes possible to tackle previously unsolvable questions. To infect us, the SARS-Cov2 virus relies on its own spike protein, designed by evolution to be stable and efficiently bind to the human ACE2 receptor. Using a protein structure produced in the early months of 2020, we exhaustively enumerated SARS-CoV2 variants that would, in theory, bind to ACE2 and kept those that remained sufficiently stable. After a drastic selection among tenths of millions of predicted variants, 59 sequences were tested experimentally for affinity, infectivity, and resistance to antibodies, resulting in a list of non-yet-existing infectious therapeutic-antibodies-resistant variants that could be used to design vaccines proactively [3].

### 3 Learning how to play the Protein Design and Sudoku games

Because the laws of physics and modeling assumptions used in such approaches lead to approximate results, it becomes crucial to exploit the massive amount of data that has been produced by experimentalists in terms of natural protein structures and sequences. This raises the exciting question of learning CFNs describing the “quality” of sequences for a given protein structure to eventually learn how to design proteins. This problem is reminiscent of learning “how to reason” or “how to play Sudoku” which has been addressed by various recent decision-focused learning architectures. By leveraging a usual probabilistic interpretation of CFNs, we recently proposed a simple scalable learning architecture [4] that combines Deep Learning with an exact CFN solver (toulbar2 [2]) to learn how to design proteins (or how to play Sudoku) which outperforms existing architectures in terms of training time, data-efficiency and accuracy. Because solving the WCSP is NP-hard, powerful polynomial time relaxations then become handy [5].

#### References

1 David Allouche, Isabelle André, Sophie Barbe, Jessica Davies, Simon de Givry, George Katsirelos, Barry O'Sullivan, Steve Prestwich, Thomas Schiex, and Seydou Traoré. Computational protein design as an optimization problem. *Artificial Intelligence*, 212:59-79, 2014.

2 David Allouche, Simon De Givry, George Katsirelos, Thomas Schiex, and Matthias Zytnicki. Anytime hybrid best-first search with tree decomposition for weighted CSP. In *Principles and Practice of Constraint Programming: 21st International Conference, CP 2015, Cork, Ireland, August 31-September 4, 2015, Proceedings 21*, pages 12-29. Springer, 2015.

3 Mireia Solà Colom, Jelena Vucinic, Jared Adolf-Bryfogle, James W Bowman, Sébastien Verel, Isabelle Moczygemba, Thomas Schiex, David Simoncini, and Christopher D Bahl. Deep evolutionary forecasting identifies highly-mutated SARS-CoV-2 variants via functional sequence-landscape enumeration. *Research Square*, pages rs-3, 2022.

4 M. Defresne, S. Barbe, and T. Schiex. Scalable coupling of deep learning with logical reasoning. In *Proc. of the 32nd IJCAI*, Macau, A.S.R., China, 2023.

5 Valentin Durante, George Katsirelos, and Thomas Schiex. Efficient low rank convex bounds for pairwise discrete graphical models. In *International Conference on Machine Learning*, pages 5726-5741. PMLR, 2022.

6 Mark A Hallen and Bruce R Donald. Protein design by provable algorithms. *Communications of the ACM*, 62(10):76-84, 2019.

7 John Jumper, Richard Evans, Alexander Pritzel, Tim Green, Michael Figurnov, Olaf Ronneberger, Kathryn Tunyasuvunakool, Russ Bates, Augustin Žídek, Anna Potapenko, et al. Highly accurate protein structure prediction with alphafold. *Nature*, 596(7873):583-589, 2021.
8 Hiroki Noguchi, Christine Addy, David Simoncini, Staf Wouters, Bram Mylemans, Luc Van Meervelt, Thomas Schiex, Kam YJ Zhang, Jeremy RH Tame, and Arnout RD Voet. Computational design of symmetrical eight-bladed β-propeller proteins. *IUCrJ*, 6(1):46–55, 2019.

9 Ilan Samish, editor. *Computational Protein Design*. Humana New York, NY, 2017. doi: 10.1007/978-1-4939-6637-0.
<figure><img src="image_26.png" /></figure>
# CP Solver Design for Maximum CPU Utilization

## Petr Vilím

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

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

ScheduleOpt, Nový Knín, Czech Republic

## Abstract

In this talk, I explain how to improve the performance of a solver without focusing on algorithms, search, propagation or parallelism. Performance is achieved instead with better CPU utilization, efficient code and more precise design of the solver itself.

In the words of Fedor G. Pikus [1], the time of “performance taking care of itself” is over. In today’s hardware the number of cores is increasing while the CPU clock speed has reached a plateau. Main memory access is slow in comparison to the CPU. And despite multiple memory cache levels, the CPU can easily become idle waiting for data from the memory, slowing down the computation considerably. Unfortunately, those trends are probably not going to change in the near future.

For those reasons we are witnessing revived interest in efficient code and performance-centered software design, especially in areas where the performance is critical: computer games, compilers, internet browsers, language interpreters (e.g. JavaScript or Python), etc.

The good news is that many of the tricks used in the above-mentioned areas, can be used in constraint programming as well. The bad news is that the performance has to be taken into account from the very beginning of the design. It is not possible to add it easily later. Sometimes, better performance can be achieved only by radical shifts in the design such as from object-oriented to data-oriented programming.

The design of a CP solver is not an exception in this regard. Without the efficient core of the CP solver, it is not possible to write truly efficient propagation or search algorithms. On the other hand, all algorithms in the solver must take the design of the solver into account and leverage it.

In this talk, I will describe what I consider the most important aspects of the design of *ScheduleOpt Optal* solver. I will concentrate on the performance, but I will also mention other aspects such as ease of use, maintainability, and testing.

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

**Keywords and phrases** Constraint Programming, Software Design, Efficient Code

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

**Category** Invited Talk

## References

1. F.G. Pikus. *The Art of Writing Efficient Programs: An advanced programmer's guide to efficient hardware utilization and compiler optimizations using C++ examples*. Packt Publishing, 2021.

<figure><img src="image_29.png" /></figure>
<figure><img src="image_30.png" /></figure>
# Optimization of Short-Term Underground Mine Planning Using Constraint Programming

Younes Aalian

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

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

Département de mathématiques et de génie industriel, Polytechnique Montréal, Québec, Canada

Gilles Pesant

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

Département de génie informatique et génie logiciel, Polytechnique Montréal, Québec, Canada

Michel Gamache

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

Département de mathématiques et de génie industriel, Polytechnique Montréal, Québec, Canada

## Abstract

Short-term underground mine planning problems are often difficult to solve due to the large number of activities and diverse machine types to be scheduled, as well as multiple operational constraints. This paper presents a Constraint Programming (CP) model to optimize short-term scheduling for the Meliadine underground gold mine in Nunavut, Canada, taking into consideration operational constraints and the daily development and production targets of the mine plan. To evaluate the efficacy of the developed CP short-term planning model, we compare schedules generated by the CP model with the ones created manually by the mine planner for two real data sets. Results demonstrate that the CP model outperforms the manual approach by generating more efficient schedules with lower makespans.

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

**Keywords and phrases** Mine planning, Constraint Programming, Short-term planning, Underground mine, Scheduling

Digital Object Identifier 10.4230/LIPIcs.CP.2023.6

# 1 Introduction

The mining industry is an important component of Canada's economic vitality. In 2019, its economic contribution was estimated at $ 109 billion, or 5 % of Canada's GDP [10]. Mining projects involve a variety of operations that handle significant amounts of material and require substantial investment. Even small reductions in costs or increases in ore yield can have a considerable economic impact. These projects can generate significant profits when they are managed efficiently. Furthermore, the mining industry is evolving and transitioning towards automated mining. With the advent of new communication and data collection tools, mining operation data is becoming more easily accessible. This creates opportunities to develop new optimization tools that can use the available data to enhance the operational efficiency in mines.

The model presented in this study is designed for an underground gold mine. The price of gold is set by the market and the same for all mining companies. Among other things, 47 % of the gold produced in Canada is purchased by the London Bullion Market, which trades gold worldwide. The only options for gold mines to increase their profits is to reduce their operating costs. One way to reduce operating costs is to make better use of available resources. Minimizing the makespan indirectly reduces operating costs by doing more activities with the same equipment and reducing downtime.

Short-term planning in underground mines plays a crucial role in ensuring the profitability and success of mining operations. It involves allocating resources to activities and determining the sequence and start time of activities during each work shift over a planning horizon

<figure><img src="image_35.png" /></figure>
## 6:2 Short-Term Underground Mine Planning Using CP

ranging from one to two weeks [1, 3]. Currently, scheduling decisions in underground mines are typically made manually based on the planner's experience. Planning has been done manually for several reasons. First, communication systems in the underground mines were virtually non-existent. As a result, the exchange of information between the planning teams was essentially done between shifts. In addition, the management systems are not yet standardized in the mines, which means that information on geology, equipment maintenance and production management are found in different systems and the transfer of one system to another is not always trivial. However, manual planning is prone to errors and often results in infeasible schedules with low accuracy and efficiency. Therefore, developing a decision tool to optimize short-term scheduling in underground mines can help achieve high-quality schedules, improve mine productivity, and reduce reliance on the planner's experience, while ensuring technical and safety requirements are met [17]. In this paper, a Constraint Programming (CP) model is presented for the short-term scheduling of activities at the Meliadine underground gold mine located in Nunavut, Canada. The model considers both operational constraints and the mine's development and production targets to generate more practical and reliable schedules.

### 1.1 Why CP?

Previous research has shown that Constraint Programming is an effective and efficient method for solving scheduling problems across various industries, including planning, scheduling, transportation, and automated systems [9]. CP uses a wide variety of variable types, functions, and global constraints to offer modeling at a high level of abstraction, making it a more flexible and intuitive approach than other model-based methods such as Mixed Integer Programming (MIP) [4]. Consequently CP models are more concise and require fewer decision variables and constraints which makes them an attractive tool for addressing large-scale scheduling problems. In the context of underground mining, the use of CP functions (as described in detail in Section 3) makes it easier to model operational constraints in the short-term scheduling problem, resulting in a more compact and efficient model.

### 1.2 Plan of the Paper

Section 2 describes the problem we address, Section 3 introduces the CP model we developed to solve it, and Section 4 discusses the outcomes of implementing the presented model on two actual data sets. Section 5 highlights the advantages of using CP for this short-term underground mine planning problem. Section 6 presents an overview of related computational approaches in the literature. Finally Section 7 concludes the paper.

## 2 Problem Description

Underground mining operations involve two primary categories of activities: development and production. In order to access economically valuable ore deposits, development activities are conducted in waste rocks that lack financial value. Production activities take place in economically significant rocks located in areas referred to as stopes [5]. Mining activities occur in a cycle at one of several sites that serve as a workplace to perform these activities. Figures 1a and 1b show the development and production cycles with activities arranged in a sequence-dependent order. Table 1 provides the description of activities in the cycle, along with the required machine type. There are several machines available for each type of activity.
Each machine can be viewed as a renewable unitary resource, limited to performing one activity at a time. Short-term scheduling for underground mines includes assigning activities in the cycle to eligible machines and determining the start and end times for each activity [1].

<figure><img src="image_36.png" /><figcaption>Figure 1 Typical development (a) and production (b) cycles in underground mining.</figcaption></figure>

<table><caption>Table 1 Activities and the required machine type in the cycle.</caption><thead><tr><td>Activity</td><td>Machine</td><td>Description</td></tr></thead><tbody><tr><td>Drilling</td><td>Drilling rigs</td><td>Drilling blast holes in the rock face</td></tr><tr><td>Charging</td><td>Anfo loader</td><td>Charging drilled holes with explosives</td></tr><tr><td>Loading</td><td>Scooptram</td><td>Removing broken rocks after blasting</td></tr><tr><td>Bolting</td><td>Bolter</td><td>Stabilizing drifts by installing bolts into the rock mass</td></tr><tr><td>Cleaning</td><td>Scooptram</td><td>Removing small rocks from the site (the gallery)</td></tr><tr><td>Cabling</td><td>Cabling machine</td><td>Reinforcing stope by installing steel cables into rock mass</td></tr><tr><td>Slot raising</td><td>Raise borer</td><td>Creating a vertical or inclined hole into the rock</td></tr></tbody></table>

There are several underground mining methods for extracting deep mineral deposits. The Meliadine mine uses the long-hole stoping method, which is one of the most commonly-used underground mining techniques that involve extracting a significant amount of material from each stope (Figure 2). This method is particularly suitable for large-scale and steeply dipping ore deposits with preferably tabular shapes. The long-hole stoping method begins with the development of main shafts or declines for transportation and ventilation purposes. Next, drill drives are excavated to access the intended location of the ore body and to create stopes. In each stope, production holes are drilled and charged with explosives. Once the blasting is completed, the fragmented rock is accessed through draw points developed at the bottom level of the stope. Scoop trams and trucks are used to collect the broken ore and transport it to the surface or other underground locations via drifts or ramps. In the final stage, the evacuated space in the stope is filled with a mixture of waste rocks and concrete to provide sufficient stability for the subsequent adjacent opening stopes [16].

At Meliadine work is organized into a succession of day and night shifts, each lasting 55 time units. Blasting activities are performed only during designated blast windows. A blast window corresponds to the period between shifts in the morning, during which resources cannot be used by operators due to safety regulations. The team rotation takes approximately 1 hour and 30 minutes, and the blast window requires roughly 4 hours and 30 minutes, including the time needed for team rotation and gas clearance (18 working hours for both
<figure><img src="image_37.png" /><figcaption>Figure 2 Typical representation of the long-hole stoping operation [8].</figcaption></figure>

day and night shifts + 1 hour 30 minutes for team rotation at the end of day shift + 4 hours
30 minutes for blast window at the end of night shift = 24 hours). The shift at the end of
which team rotation occurs is referred to as the day shift, while the subsequent shift, which
includes the blasting window at the end, is known as the night shift. Figure 3 illustrates the
shift organization in the studied underground mine. Mining activities are preemptive as they
can be interrupted at the end of each shift and continue in the next shift.

Shift Organization

<figure><img src="image_38.png" /><figcaption>■ **Figure 3** Timeline of alternating day and night shifts including time to rotate the teams and to perform blasting (above). Its representation in the CP model (below).</figcaption></figure>

Short-term planning at the Meliadine mine incorporates several key performance indicators (KPI) such as progress of development rounds in the drift, total length of production holes drilled in the stope, and total amount of material mucked from the stope to meet the medium-term planning goals. The KPI values vary monthly and are updated every three months by the medium-term mine planner. Development and production constraints will be introduced in our CP model to consider the defined KPIs in short-term scheduling.
# 3 How our Problem is Modeled in CP

An optimization model is developed using Constraint Programming (CP) for short-term underground mine scheduling, taking into account operational requirements of underground mining operations. Additional constraints are introduced to ensure that the mine planning development and production targets are met and that practical and reliable short-term schedules are generated. In other words, the produced schedule determines the detailed execution of mining activities in underground operations considering the required daily rates of development and production. It is important to note that the same model can be used for both development and production activities in underground mining, which ensures consistency and accuracy in short-term scheduling.

CP Optimizer (CPO) from IBM ILOG Optimization Studio [9] was used to create the model presented in this article. In this CPO model, interval variables are used to represent activities, each with several related optional interval variables depicting the choice of resource. Optional interval variables include a Boolean status reflecting the fact that the corresponding activity is present or absent from the solution (i.e. not considered by the constraints). The ordering of resources can be represented by a set of interval variables, known as a sequence variable. This sequence variable is used in the scheduling model to prevent activities in the sequence from overlapping in time. More formally, an interval variable *a* is defined by a start time *s* and an end time *e*, which are non-negative integer values, such that $a \in \{[s, e) | s, e \in \mathbb{N}, e \geq s\}$. Optional interval variable *b* is presented such that $b \in \{\emptyset\} \cup \{[s, e) | s, e \in \mathbb{N}, e \geq s\}$. Additionally the developed CPO model uses various functions and constraints that are described as follows [9]:

- **endOf**: A function that provides the end value of an interval variable if it exists, or else returns zero.

- **alternative**: This constraint ensures that if a given interval variable is present, then only one related optional interval variable is chosen with the same start and end values.

- **noOverlap**: This constraint is used to ensure that a set of interval variables defined by a sequence variable do not overlap, while maintaining a minimum distance between them as specified by a transition distance matrix.

- **endBeforeStart**: This constraint guarantees that if two interval variables are present, then the first ends before the second starts, with an optional minimum delay between them.

- **forbidExtent**: This constraint makes sure that an interval variable cannot overlap with a forbidden region where the value of the step function is zero. As a result, the interval variable must either end before the forbidden region or start after it.

- **stepAtEnd**: This step function returns an elementary cumulative function with a non-negative integer value at the end of an interval variable. Such functions model a known function of time, such as the resources used during a particular time period, by returning a non-negative integer value (height of the elementary function) within the range of the interval variable and zero outside of it.

- **cumulFunction**: This expression models a known function of time, such as the cumulative amount of resources used by an activity during a specific time period.

- **alwaysIn**: This constraint restricts the potential values of a cumulative function to a specific range during a time interval.
## 6:6 Short-Term Underground Mine Planning Using CP

Tables 2 and 3 present lists of sets, parameters, and variables used in the CP model, along with their corresponding descriptions.

<table><caption>Table 2 Sets and parameters of the CP model.</caption><thead><tr><th>Set</th><th>Description</th></tr></thead><tbody><tr><td>J</td><td>Index set of activities</td></tr><tr><td>M</td><td>Index set of all available equipment</td></tr><tr><td>M<sub>j</sub></td><td>Index set of eligible machines to perform activity j</td></tr><tr><td>A<sub>j</sub></td><td>Index set of activities that must occur after activity j</td></tr><tr><td>B</td><td>Index set of blast activities</td></tr><tr><td>T</td><td>Index set of time windows (starting at 1)</td></tr><tr><th>Parameter</th><th>Description</th></tr><tr><td>p<sub>j</sub></td><td>Processing time of activity j</td></tr><tr><td>D</td><td>Matrix of transition time between sites where the value of its element is equal to 0 for the same site and greater than 0 otherwise</td></tr><tr><td>d<sub>j</sub></td><td>Development (meter) of activity j</td></tr><tr><td>h<sub>j</sub></td><td>Production hole drilling (meter) of activity j</td></tr><tr><td>o<sub>j</sub></td><td>Stope ore mucking (ton) of activity j</td></tr><tr><td>s<sub>t</sub></td><td>Starting time of time window t</td></tr><tr><td>e<sub>t</sub></td><td>Ending time of time window t</td></tr><tr><td>d&#773;</td><td>Lower bound for daily development</td></tr><tr><td>d&#773;&#773;</td><td>Upper bound for daily development</td></tr><tr><td>h&#773;</td><td>Lower bound for daily hole drilling</td></tr><tr><td>h&#773;&#773;</td><td>Upper bound for daily hole drilling</td></tr><tr><td>o&#773;</td><td>Lower bound for daily ore mucking</td></tr><tr><td>o&#773;&#773;</td><td>Upper bound for daily ore mucking</td></tr><tr><td>Blast\_calendar</td><td>The time periods during which only blasting activities are permitted (all activities except blasting are forbidden to be performed during these periods</td></tr></tbody></table>

<table><caption>Table 3 Decision variables of the CP model.</caption><thead><tr><th>Variable</th><th>Description</th></tr></thead><tbody><tr><td>Y<sub>j</sub></td><td>Interval variable for activity j</td></tr><tr><td>X<sub>jm</sub></td><td>Optional interval variable to perform activity j using machine m</td></tr><tr><td>S<sub>m</sub></td><td>Sequence variable for machine m (S<sub>m</sub> = {X<sub>jm</sub> | j ∈ J})</td></tr><tr><td>Q<sup>d</sup></td><td>Integer variable for total development in drifts</td></tr><tr><td>Q<sup>h</sup></td><td>Integer variable for total amount of production hole drilling in stopes</td></tr><tr><td>Q<sup>o</sup></td><td>Integer variable for total amount of ore material mucked from stopes</td></tr></tbody></table>
The CP model is given as (1)-(11):

**Objective function**

$$ \text{Minimize} \quad \max_{j \in J} (\text{endOf}(Y_j)) \qquad (1) $$

**Constraints**

$$ \text{alternative}(Y_j, X_{jm} | m \in M_j) \qquad \forall j \in J \qquad (2) $$

$$ \text{noOverlap}(S_m, D) \qquad \forall m \in M \qquad (3) $$

$$ \text{endBeforeStart}(Y_j, Y_i) \qquad \forall j \in J, i \in A_j \qquad (4) $$

$$ \text{forbidExtent}(Y_j, \text{Blast\_calendar}) \qquad \forall j \in J \setminus B \qquad (5) $$

$$ \text{cumulFunction}(Q^d) = \sum_{j \in J} \text{stepAtEnd}(Y_j, d_j) \qquad (6) $$

$$ \text{alwaysIn}(Q^d, s_t, e_t, t \times \underline{d}, t \times \bar{d}) \qquad \forall t \in T \qquad (7) $$

$$ \text{cumulFunction}(Q^h) = \sum_{j \in J} \text{stepAtEnd}(Y_j, h_j) \qquad (8) $$

$$ \text{alwaysIn}(Q^h, s_t, e_t, t \times \underline{h}, t \times \bar{h}) \qquad \forall t \in T \qquad (9) $$

$$ \text{cumulFunction}(Q^o) = \sum_{j \in J} \text{stepAtEnd}(Y_j, o_j) \qquad (10) $$

$$ \text{alwaysIn}(Q^o, s_t, e_t, t \times \underline{o}, t \times \bar{o}) \qquad \forall t \in T \qquad (11) $$

Objective (1) of the CP model is to minimize the makespan. Constraint (2) ensures that only one optional variable is chosen for an interval variable i.e. only one machine (with the appropriate type) is used to perform a given activity. Constraint (3) prevents machines from being used simultaneously, meaning that each machine can only be assigned to one activity at a time. Constraint (4) takes into account the order in which activities must be performed at a site, with most activities having only one predecessor and some having none. It is important to note that the site where each activity must be carried out is predefined in the input data. Therefore, all activities can be executed in their respective predetermined sites.

Constraint (5) is used to ensure that only blasting activities occur during designated blast windows. In order to model the blasting constraint in the CP model, the day and night work shifts are compressed into a 110-time unit period (each shift consists of 55 time units), where each time unit represents 10 minutes in the real world. This compression allows for blasting activities with a length of zero time units to occur only at the end of compressed periods, every 110 time units (see Figure 3). Multiple blasting activities can be performed at the same time during each blasting window.

Constraints (6) and (7) are introduced to ensure that the progress of development rounds each day (measured in meters per day) is maintained within specific limits based on the defined development target. To model these development constraints, we define time windows $[s_t, e_t)$ each representing a day in the schedule. For each time window, we establish cumulative upper and lower bounds (based on daily bounds) for total development in drifts (in meters) that must be achieved. Next, we use the `cumulFunction` to model the cumulative amount of development per meter and apply the `alwaysIn` constraint to ensure that the cumulative function stays within the target value bounds for each time window. The function `stepAtEnd(i, j)` returns an elementary cumulative function with a step of height j (a non-negative integer value) at the end of interval variable i. The presented development
6:8 Short-Term Underground Mine Planning Using CP

constraints aim to achieve the desired daily development target in the generated schedule. Furthermore, by using the cumulFunction in this constraint, the model is able to flexibly compensate in the following days for any shortfall in achieving the daily development goal (see e.g. Figure 6a). This feature of the constraints closely resembles what is taken into account in actual short-term underground mine planning, making the model more practical for real-world operations.

Constraints (8) and (9) model the production drilling constraint to ensure that the amount of production holes drilled in the stope per day (measured in meters per day) is restricted within certain bounds, defined based on the production drilling objectives. Furthermore, Constraints (10) and (11) are used to model the stope mucking constraint, which ensures that the amount of ore material mucked from stopes each day (measured in tons per day) is maintained within specific limits determined based on the production target. These constraints (Constraints (8)-(11)) aim to achieve the production plan in the produced schedule by controlling the daily amount of production holes drilled and ore mucked. Similar to the development constraints, the production constraints also allow for making up for shortfalls in meeting the daily production goal. An activity can perform either development or production depending on the type of cycle. If the activity is part of a development cycle, it can have development ($d_j$), and if it is involved in a production cycle, it can have either production hole drilling ($h_j$) or stope mucking ($o_j$). The development cycle includes activities that are performed in waste rocks lacking financial value to access economically valuable deposits, while the production cycle is conducted in valuable rocks to extract ore material from the stope.

In the presented short-term mine planning model, operational development and production targets (KPIs) are considered to achieve the tactical decisions made at the medium-term planning level. Specifically, tactical decisions in underground mine planning are typically associated with defining the extraction sequence over a planning horizon of one to three months [7].

4 Implementation and Results

The experiments were conducted on a computer featuring an Intel(R) Core(TM) i7-9750H CPU @ 2.60GHz and 16 GB of RAM. The CP models were solved using the Constraint Programming Optimizer in IBM ILOG CPLEX Optimization Studio version 12.8.0.

The model was tested on two real data sets collected from the Meliadine underground gold mine in Nunavut, Canada. Both data sets involve scheduling activities for a roughly one-week planning horizon. The first data set (Instance 1) relates to development operations, which consist of 15 machines and 291 activities to be performed across 18 sites. The total advancement achieved by all available development rounds in this instance is equal to 188 meters. Specifically, each round (cycle) results in approximately 4 meters of advancement in the development drift. The second data set (Instance 2) concerns production operations and includes 27 machines, 185 activities, and 27 sites. In this instance, a total of 1500 meters of production holes have been drilled across all accessible stopes, resulting in the extraction of 27,000 tons of ore material. The available resources are categorized into different equipment types, and the number of each type is reported in Table 4. Although the developed CP model takes into account both development and production activities, there was no data available (in the mine) that included both activities together. Therefore, we applied our model separately to two different datasets: one for development and another for production.
<table><caption>Table 4 Number of machines per equipment type for Instances 1 and 2.</caption><thead><tr><td>Equipment type</td><td>Instance 1</td><td>Instance 2</td></tr></thead><tbody><tr><td>Scooptram</td><td>2</td><td>7</td></tr><tr><td>Bolter</td><td>6</td><td>-</td></tr><tr><td>Scooptram clean face</td><td>1</td><td>-</td></tr><tr><td>Jumbo</td><td>3</td><td>-</td></tr><tr><td>Anfo loader</td><td>3</td><td>3</td></tr><tr><td>Truck</td><td>-</td><td>7</td></tr><tr><td>Raise borer</td><td>-</td><td>1</td></tr><tr><td>Production drill rig</td><td>-</td><td>4</td></tr><tr><td>Cabling machine</td><td>-</td><td>5</td></tr></tbody></table>

The results obtained by implementing the CP model on Instances 1 and 2 are presented in the following subsections.

## 4.1 Instance 1

Table 5 presents the results of schedules generated for Instance 1 with different daily development upper bounds ($\bar{d}$) in the development constraint. All the models in the table are solved to optimality in a short amount of time. As the primary objective of the scheduling model is to minimize the makespan, the lower bound on daily development specified in the development constraint is readily satisfied. Therefore to evaluate the effect of different development targets on the resulting schedule, we only modify the upper bound value for the total amount of development to be accomplished per day.

<table><caption>Table 5 Results of different CP models on Instance 1.</caption><thead><tr><td>Model</td><td>Development upper bound ($\bar{d}$)</td><td>Makespan</td><td>Solving time</td></tr></thead><tbody><tr><td>1</td><td>24</td><td>882</td><td>12 sec</td></tr><tr><td>2</td><td>28</td><td>772</td><td>13 sec</td></tr><tr><td>3</td><td>32</td><td>678</td><td>12 sec</td></tr><tr><td>4</td><td>36</td><td>635</td><td>13 sec</td></tr><tr><td>5</td><td>40</td><td>600</td><td>11 sec</td></tr><tr><td>6</td><td>44</td><td>600</td><td>10 sec</td></tr><tr><td>7</td><td>∞</td><td>600</td><td>10 sec</td></tr></tbody></table>

As can be seen from Table 5, increasing the upper bound in the development constraint results in lower makespans in the produced schedule. Furthermore, the schedule makespan remains unchanged for bounds greater than 40. Therefore, $\bar{d} = 40$ can be considered a suitable daily development target for generating a short-term schedule on Instance 1. Interestingly, this value coincides with the daily development target employed by the human planner at the mine – our model confirms this empirical choice. Figure 4 displays the location-based Gantt chart for the short-term schedule generated on Instance 1 with $\bar{d} = 40$.

The daily and cumulative development resulting from the schedule produced using Model 3 on Instance 1 with $\bar{d} = 32$ are displayed in Figures 5a and 5b. As seen in Figure 5a, the maximum daily development limit of 30 meters is respected, resulting in a total cumulative development of 188 meters in six days (Figure 5b).
6:10 Short-Term Underground Mine Planning Using CP

<figure><img src="image_39.png" /><figcaption>Figure 4 Location-based Gantt chart for the generated schedule on Instance 1.</figcaption></figure>

<figure><img src="image_40.png" /><figcaption>Figure 5 Daily (a) and cumulative (b) development in Model 3 ($\bar{d} = 32$) on Instance 1.</figcaption></figure>

Figures 6a and 6b show the daily and cumulative development obtained from Model 6 on Instance 1 with $\bar{d} = 44$. According to Figure 6a, 36 meters of development are achieved on Day 2, which is lower than the maximum daily target of 44 meters. However, this shortfall is made up on Day 3 by completing 48 meters of development, above the maximum daily target. In other words, 48 meters of development are completed on Day 3 to compensate for the shortfall on Day 2. After Day 3, it is not possible to meet the maximum daily goal due to the limited number of drifts available. As shown in Figure 6b, the total cumulative development of 188 meters is reached in five days. This feature of the development constraints in the CP model can be practical for short-term planning in underground mines, where operational restrictions or a relatively small number of accessible drifts (sites) prevent the achievement of the development target on certain days.

Figure 7 shows the comparison of the average utilization rate of several machine types in schedules produced using CP models on Instance 1 with different $\bar{d}$ for the development constraint. The utilization rate of a machine is the total amount of time units during which the machine was actively operating at the site relative to the total amount of time for which it was available for use. As can be seen from this figure, increasing $\bar{d}$ results in a higher
<figure><img src="image_41.png" /><figcaption>Figure 6 Daily (a) and cumulative (b) development in Model 6 ($\bar{d} = 44$) on Instance 1.</figcaption></figure>

average utilization rate of machines in the schedule. This is due to the fact that larger $\bar{d}$ values lead to more compact schedules with lower makespans, which in turn, reduces waiting time for machines.

<figure><img src="image_42.png" /><figcaption>Figure 7 Average utilization rate of machines in schedules with different development bounds ($\bar{d}$) on Instance 1.</figcaption></figure>

## 4.2 Instance 2

Table 6 displays the makespan of schedules generated by implementing different models on Instance 2, with distinct upper bounds for the daily production drilling ($\bar{h}$) and stope mucking ($\bar{o}$) in production constraints. According to Table 6, reducing the upper bound values in production constraints leads to longer makespans in the generated schedule. Specifically, for the production drilling constraint, the suitable $\bar{h}$ value is 400 meters, as it leads to the lowest makespan value that remains unchanged for larger upper bounds. Similarly, for the stope mucking constraint, $\bar{o} = 6,000$ appears to be an appropriate stope mucking target for the schedule generated on Instance 2. Figure 8 shows the location-based Gantt chart for the created schedule on Instance 2 with $\bar{o} = 6,000$.

Figures 9a and 9b present the daily and cumulative production drilling rates in the schedule generated using Model 5 with $\bar{h} = 500$ on Instance 2. Figure 9a demonstrates that the drilling rate exceeds the daily limit by reaching 600 meters on Day 2 to compensate for the shortfall on Day 1. As depicted in Figure 9b, the total production drilling of 1500 meters is achieved within four days.
6:12 Short-Term Underground Mine Planning Using CP

<table><caption>Table 6 Results of different CP models on Instance 2.</caption><thead><tr><th>Model</th><th>Production drilling upper bound (h&#772;)</th><th>Stope mucking upper bound (o&#772;)</th><th>Makespan</th><th>Solving time</th></tr></thead><tbody><tr><td>1</td><td>&infin;</td><td>&infin;</td><td>730</td><td>16 sec</td></tr><tr><td>2</td><td>200</td><td>-</td><td>1060</td><td>17 sec</td></tr><tr><td>3</td><td>300</td><td>-</td><td>840</td><td>16 sec</td></tr><tr><td>4</td><td>400</td><td>-</td><td>730</td><td>16 sec</td></tr><tr><td>5</td><td>500</td><td>-</td><td>730</td><td>16 sec</td></tr><tr><td>6</td><td>-</td><td>4000</td><td>881</td><td>17 sec</td></tr><tr><td>7</td><td>-</td><td>5000</td><td>771</td><td>16 sec</td></tr><tr><td>8</td><td>-</td><td>6000</td><td>730</td><td>17 sec</td></tr><tr><td>9</td><td>-</td><td>7000</td><td>730</td><td>17 sec</td></tr></tbody></table>

<figure><img src="image_43.png" /><figcaption>Figure 8 Location-based Gantt chart for the generated schedule on Instance 2.</figcaption></figure>

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

(a)

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

(b)

<caption>Figure 9 Daily (a) and cumulative (b) drilling in Model 5 ($\bar{h} = 500$) on Instance 2.</caption>

Figures 10 and 11 compare the average utilization rate of several machine types in schedules produced using CP models on Instance 2, with different values for $\bar{h}$ and $\bar{o}$, respectively. According to these figures, the average utilization rate of machines increases for schedules with higher upper bounds in production constraints.
<figure><img src="image_46.png" /><figcaption>Figure 10 Average utilization rates of machines in schedules with different production drilling bounds ($\bar{h}$) on Instance 2.</figcaption></figure>

<figure><img src="image_47.png" /><figcaption>Figure 11 Average utilization rates of machines in schedules with different stope mucking bounds ($\bar{o}$) on Instance 2.</figcaption></figure>

# 5 Added Value of CP

Constraint Programming allowed us to efficiently address short-term underground mine planning by quickly producing optimal schedules minimizing the makespan. Moreover the model identifies $\bar{d} = 40$ as the appropriate daily development target, which aligns with the value selected by the mine planner and confirms the practice of setting the development upper bound at 40. Additionally, the model can explore what-if scenarios by varying parameter values, such as the impact of changing daily development or production targets in the model on machine utilization rates in the generated schedule. These results demonstrate the practicality and efficiency of using the CP model for short-term scheduling in real-world underground mining operations.

## 5.1 Comparison of CP model and manual approach

In order to demonstrate the effectiveness of our optimization model, we compared the short-term schedules produced by the CP model with those manually created by the mine planner for the same instance. Since a detailed schedule of activities with similar time fidelity to the schedule produced using the CP model was not provided in the studied mine, we only compared the schedule makespan. In particular, we compared the number of shifts required to complete all activities in the generated short-term schedule using the CP model and manual approach for both Instances 1 and 2, as shown in Table 7. The development KPI considered for scheduling activities in Instance 1 is 40 meters per day (m/day). In Instance
## 6:14 Short-Term Underground Mine Planning Using CP

2, the production drilling KPI is 400 m/day, and the stope mucking KPI is 6,000 tons per day. As previously mentioned, each day consists of two working shifts, where each shift is equivalent to 55 time units in the CP model.

<table><caption>**Table 7** Comparison of schedule makespan between CP model and manual approach for Instances 1 and 2.</caption><thead><tr><td>Instance</td><td>CP model</td><td>Manual approach</td></tr></thead><tbody><tr><td>1</td><td>11 Shifts</td><td>14 Shifts</td></tr><tr><td>2</td><td>14 Shifts</td><td>16 Shifts</td></tr></tbody></table>

Table 7 shows that the CP approach outperforms the manual scheduling method on both Instances 1 and 2 by creating more compact schedules with lower makespans (based on the number of shifts) while satisfying the daily development and production targets (KPIs). Additionally, the CP models are quickly solved to optimality, making it an efficient tool for mine planners to rapidly generate updated short-term schedules whenever changes occur in the underground mine plan.

The results of this study demonstrate advantages of the developed CP model for optimizing short-term planning in underground mines and reducing the reliance on manual scheduling, which is highly dependent on the planner's experience. Moreover, the CP model can be easily adjusted to accommodate or exclude additional activity types and related constraints based on the specific requirements of underground mining operations.

## 6 Literature review

Short-term underground mine planning models are often difficult to solve (NP-hard) due to various operational constraints to consider and to the large number of variables involved. However there has been notable research interest in developing new mathematical models and algorithms to optimize short-term scheduling in underground mines.

Nehring et al. (2010) designed a MIP model to optimize the short-term scheduling and allocation of loader-trucks in sublevel stoping mines. The model allows for the reallocation of equipment in response to changes in underground operations. The proposed model was applied to a copper mine, demonstrating satisfactory results in terms of tonnage deviations from predetermined amounts throughout the planning period [11]. O'Sullivan and Newman (2015) introduced an Integer Programming (IP) model for scheduling activities in an Irish lead and zinc underground mine to maximize the discounted amount of produced metal. Both exact and heuristic solutions were used to reduce the number of variables in the model. Additionally, an optimization-based decomposition heuristic was developed to generate feasible schedules in less computation time for complicated problem instances [12]. Song et al. (2015) developed a decision support tool to determine the scheduling of activities in underground mines. The tool was tested on a real mine dataset in Finland and significantly decreased the makespan compared to manual scheduling methods, thereby improving operational performance. However, the proposed method did not take into account uncertainty related to unexpected activities in underground operations [15].

Schulze and Zimmermann (2017) introduced a solution approach for short-term production scheduling in underground mining. The developed approach assigns staff and machines to mining activities while considering operational constraints with the goal of minimizing deviations from targeted production in a potash mine. The method was tested on various instances and demonstrated superior performance when compared to manual scheduling [13]. Seifi et al. (2019) proposed a two-stage solution approach for scheduling machines and staff in an underground potash mine in Germany. The first step involves solving the relaxation
of the MIP model, and in the second step, a heuristic algorithm is used to modify the
solutions obtained from the relaxation model to achieve feasible schedules. The experiments
conducted on real-word datasets show that the developed approach outperforms the heuristic
procedure presented by Schulze and Zimmermann (2017)[14]. Wang et al. (2020) utilized
a genetic algorithm (GA) for optimizing the scheduling of equipment used in underground
mining. A Non-Linear Programming (NLP) model is presented with a significant number of
decision variables associated with multiple mining sites and equipment types [17]. A MIP
model was presented by Campeau and Gamache (2020) to optimize short-term planning in
underground mines. The goal was to maximize material extraction while ensuring a minimum
ore production rate to keep the mill active. The model considers operational and resource
constraints to generate feasible schedules. When applied to a gold mine data set, the model
produced an optimal short-term schedule [5]. Campeau et al. (2022) introduced a novel
MIP model to address short- and medium-term planning in underground mines. The model
integrated continuous variables for time discretization, resulting in realistic schedules. The
effectiveness of the model was demonstrated by applying it to a dataset from a Canadian
gold mine, which produced promising results [7].

Over the last few years, several CP approaches have been proposed to tackle the short-
term underground mine planning problem. A model using CP was suggested by Astrand et al.
(2018) for scheduling a mobile fleet in underground operations, which was tested on data from
an actual underground mine [2]. Astrand et al. (2020) extended the previously developed CP
model by incorporating the time it takes for mobile machinery to travel between different sites
in an underground cut-and-fill mine. They also proposed a revised CP model with compressed
blasting time and post-processed solutions to obtain schedules for the primary problem.
In order to improve the quality of schedules and reduce computation time, a specialized
neighborhood definition was implemented in a Large Neighborhood Search (LNS) algorithm.
The effectiveness of this algorithm was assessed using several instances of an underground
mine in Sweden. The outcome showed that the suggested method successfully enhanced the
initial feasible solution and generated high-quality schedules [3]. Campeau and Gamache
(2022) presented a CP model for short- and medium-term planning in underground mining.
They evaluated the model's ability to address long-term production planning objectives by
testing it on five data sets from a Canadian underground gold mine, considering a planning
horizon of up to one year. The outcomes revealed that the CP model was superior to the
equivalent MIP model in terms of computational efficiency and application [6].

These previous CP approaches for short-term underground mine planning exhibit a
limited ability to incorporate daily mine planning development and operational goals during
the short-term scheduling process. To overcome this limitation, this paper introduced a
CP model for the short-term scheduling of activities in underground mining that takes into
consideration operational constraints and the development and production targets of the
mine plan to generate more practical and reliable schedules.

7 Conclusion

This paper presented a CP model that takes into account various operational constraints
and daily development and production targets for short-term scheduling optimization in
underground mines. The model was tested on two data sets from the Meliadine gold mine
using the long-hole stoping mining method. We conducted a comparative analysis of the
schedules generated by our CP model and those created manually by the mine planner.
The experiments showed that the CP model outperforms the manual approach, resulting
in more efficient schedules with lower makespans. Results highlight the potential benefits of implementing the CP model in actual underground mining operations to improve both development and production through optimized short-term mine planning. Underground mines are somewhat unpredictable environments which may affect how long an activity actually takes. For future work, it could be beneficial to incorporate uncertainty in activity durations which would improve the robustness of the short-term schedule.

References

1 Max Åstrand, Mikael Johansson, and Jenny Greberg. Underground mine scheduling modelled as a flow shop: a review of relevant work and future challenges. *Journal of the Southern African Institute of Mining and Metallurgy*, 118(12):1265–1276, 2018.

2 Max Åstrand, Mikael Johansson, and Alessandro Zanarini. Fleet scheduling in underground mines using constraint programming. In *Integration of Constraint Programming, Artificial Intelligence, and Operations Research: 15th International Conference, CPAIOR 2018, Delft, The Netherlands, June 26–29, 2018, Proceedings 15*, pages 605–613. Springer, 2018.

3 Max Åstrand, Mikael Johansson, and Alessandro Zanarini. Underground mine scheduling of mobile machines using constraint programming and large neighborhood search. *Computers & Operations Research*, 123:105036, 2020.

4 Philippe Baptiste, Claude Le Pape, and Wim Nuijten. *Constraint-based scheduling: applying constraint programming to scheduling problems*, volume 39. Springer Science & Business Media, 2001.

5 Louis-Pierre Campeau and Michel Gamache. Short-term planning optimization model for underground mines. *Computers & Operations Research*, 115:104642, 2020.

6 Louis-Pierre Campeau and Michel Gamache. Short-and medium-term optimization of underground mine planning using constraint programming. *Constraints*, 27(4):414–431, 2022.

7 Louis-Pierre Campeau, Michel Gamache, and Rafael Martinelli. Integrated optimisation of short-and medium-term planning in underground mines. *International Journal of Mining, Reclamation and Environment*, 36(4):235–253, 2022.

8 Atlas Copco. Mining methods in underground mining. *Atlas Copco: Nacka, Sweden*, 2007.

9 Philippe Laborie, Jérôme Rogerie, Paul Shaw, and Petr Vilím. Ibm ilog cp optimizer for scheduling: 20+ years of scheduling with constraints at ibm/ilog. *Constraints*, 23:210–250, 2018.

10 B. Marshall. Facts and figures : The state of canada's mining industry. Technical report, The Mining Association of Canada, 2020.

11 Micah Nehring, Erkan Topal, and Peter Knights. Dynamic short term production scheduling and machine allocation in underground mining using mathematical programming. *Mining Technology*, 119(4):212–220, 2010.

12 Dónal O'Sullivan and Alexandra Newman. Optimization-based heuristics for underground mine scheduling. *European Journal of Operational Research*, 241(1):248–259, 2015.

13 Marco Schulze and Jürgen Zimmermann. Staff and machine shift scheduling in a german potash mine. *Journal of Scheduling*, 20:635–656, 2017.

14 Cinna Seifi, Marco Schulze, and Jürgen Zimmermann. A two-stage solution approach for a shift scheduling problem with a simultaneous assignment of machines and workers. In *Mining Goes Digital*, pages 377–385. CRC Press, 2019.

15 Zhen Song, Håkan Schunnesson, Mikael Rinne, and John Sturgul. Intelligent scheduling for underground mobile mining equipment. *PloS one*, 10(6):e0131003, 2015.

16 Farzad Sotoudeh, Micah Nehring, Mehmet Kizil, Peter Knights, and Amin Mousavi. Production scheduling optimisation for sublevel stoping mines using mathematical programming: A review of literature and future directions. *Resources Policy*, 68:101809, 2020.

17 Hao Wang, Victor Tenorio, Guoqing Li, Jie Hou, and Nailian Hu. Optimization of trackless equipment scheduling in underground mines using genetic algorithms. *Mining, Metallurgy & Exploration*, 37:1531–1544, 2020.
# Exploiting Configurations of MaxSAT Solvers

Josep Alòs

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

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

Logic & Optimization Group (LOG), University of Lleida, Spain

Carlos Ansótegui

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

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

Logic & Optimization Group (LOG), University of Lleida, Spain

Josep M. Salvia

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

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

Logic & Optimization Group (LOG), University of Lleida, Spain

Eduard Torres

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

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

Logic & Optimization Group (LOG), University of Lleida, Spain

## Abstract

In this paper, we describe how we can effectively exploit alternative parameter configurations to a MaxSAT solver. We describe how these configurations can be computed in the context of MaxSAT. In particular, we experimentally show how to easily combine configurations of a non-competitive solver to obtain a better solving approach.

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

**Keywords and phrases** maximum satisfiability, maxsat evaluation, automatic configuration

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

**Funding** This work was supported by MCIN/AEI/10.13039/501100011033 (Grant: PID2019-109137GB-C21), Agencia de Gestio d'Ajuts Universitaris i de Recerca (AGAUR), Departament d'Empresa i Coneixement de la Generalitat de Catalunya (Grant: 2022 FI_B 00010)

**Acknowledgements** We want to thank Alexander Nadel for sharing the solver TT-Open-WBO with the configurable parameters exposed.

## 1 Introduction

Since 2006, the MaxSAT Evaluation (MSE) [5] has been held annually with the primary objective of advancing MaxSAT technology and assessing its current state-of-the-art. The evaluation consists of multiple solvers being tested on various benchmarks across different evaluation tracks. This event has undeniably spurred the MaxSAT community to create more cutting-edge solvers and enhance their competitiveness.

It is not surprising that solver performance depends on several factors, including the power of the algorithm implemented by the solver, proper configuration of solver parameters to unleash its full potential, and implementation issues. Therefore, we must interpret the MaxSAT Evaluation ranking results carefully and derive conclusions according to the goal of our analysis. For example, a similar or weaker algorithm could outperform other approaches thanks to better implementation of data structures or a previous tuning process of its input parameters.

From an industrial point of view, we mainly care about obtaining an effective solving approach that is ready for deployment for a particular problem subject to available resources (computing power, environment restrictions, licenses available, etc). From a research point of view, we are more interested in identifying the potential of new solving approaches that lead to further promising research avenues.

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

© Josep Alòs, Carlos Ansótegui, Josep M. Salvia, and Eduard Torres;
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. 7; pp. 7:1–7:16

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

LIPICS

Leibniz International Proceedings in Informatics

Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl Publishing, Germany
## 7:2 Exploiting Configurations of MaxSAT Solvers

Our aim is to satisfy both industrial and research perspectives by identifying the best possible solving approach that can be achieved from a single solver while adhering to certain restrictions. In particular, we treat the solver as a black box, meaning that we cannot access its source code, nor do we have any domain knowledge of the problem to be solved, meaning that we cannot utilize any specific structure feature.

Despite these constraints, our approach enables us to unleash the hidden potential of the solver and avoid incorrect rankings of better algorithms that have not been appropriately configured or restarted. Additionally, our study emphasizes the importance of being cautious when interpreting rankings based on the MaxSAT Evaluation, as we mentioned previously.

In this paper, we first show how to effectively configure MaxSAT solvers using Automatic Configuration (AC) tools (tuners), specifically GGA [4] and SMAC [14]. Then, we show that we can take advantage of not only the best configuration returned by the tuners but also a selection of the configurations seen by the tuner during the AC process. With these configurations, we can then build a simple portfolio that runs in parallel these configurations if enough computational resources are available.

We also demonstrate how to create a sequential portfolio that schedules the execution of different parametrizations of a single MaxSAT solver on a given number of cores within a specified timeout. This approach can be thought of as a restarting strategy, where a different configuration of the solver parameters is selected at each restart.

It is worth mentioning that all these approaches are agnostic of the structure of the instances. Otherwise, we should explore extending other approaches available in the literature such as ISAC++ [13].

Finally, we integrate all these building processes in the OptiLog framework [1]. With the new APIs, the user can provide an input MaxSAT solver and its parameters through the *BlackBox Module*, and OptiLog automatically generates a new solving approach for a given number of cores.

We conducted an extensive experimental investigation on the Weighted Incomplete track of the MaxSAT Evaluation 2022, with a particular focus on the highly configurable MaxSAT solver Loandra [6]. In this track, Loandra ranked sixth when restricted to a timeout of 60 seconds. Our approach involves the construction of parallel and sequential portfolios based solely on Loandra, which significantly improves its performance.

# 2 Preliminaries

MaxSAT is the optimization variant of the SAT decision problem. While for SAT the goal is to find an assignment to the Boolean variables (solution) that satisfies all the clauses in the input CNF formula, in MaxSAT we look for a solution that satisfies the maximum possible number of clauses. Since some of these clauses can be falsified we refer to them as *soft* clauses. Within the MaxSAT community, it is typical to reformulate the problem from a minimization perspective aiming to find a solution that falsifies the minimum possible number of soft clauses.

There are several variants of the MaxSAT problem. We can add weights to the soft clauses that represent the cost of falsifying the clause. In this case, we want to look for a solution that minimizes the aggregated cost of the *Weighted* soft clauses. Additionally, we can have *hard* clauses, i.e., clauses that cannot be falsified by the solution.

MaxSAT solvers have experimented a great success in the last decade. Among these solvers, we find complete (or exact) solvers and incomplete solvers. Complete solvers provide optimal solutions while incomplete solvers report solutions as good as possible, but are not
required to guarantee their optimality. These solvers can either refine a lower bound (lb) on the cost of the optimal solutions or an upper bound (ub), or both. In particular, incomplete solvers iteratively report (whenever possible) a better (smaller) upper bound on the optimal solution.

# 3 The MaxSAT Evaluation

The MaxSAT Evaluation 2022 was structured into three tracks: main track complete (unweighted and weighted variants), main track incomplete (unweighted and weighted variants), and the special incremental MaxSAT track. In this paper, we focus on the incomplete track for weighted MaxSAT instances with a timeout of 60 seconds.

The term *incomplete* refers to the type of MaxSAT solvers which are not required to be exact, i.e., they do not need to certify the optimum. Their goal is to report the best possible solution within a given timeout. The term *weighted* refers to the variant of MaxSAT instances. The weighted MaxSAT variant allows integer weights for the soft clauses plus the hard clauses.

We consider the timeout of 60 seconds useful for our study since it is a realistic scenario of industrial applications where we require a suboptimal solution in a short time window and because our automatic configuration process, given the computational resources we have available, can be restricted to two days (see Section 5).

The MaxSAT Evaluation 2022 incomplete (weighted) track involved 197 MaxSAT instances and 10 incomplete solvers: DT-HyWalk [18], noSAT-MaxSAT [15], NuWLS-c [7], Exact [8, 11], Loandra [6], Open-WBO-inc (two variants) [12], and TT-Open-WBO-Inc (three variants) [16].

Each solver *s* was ranked according to the scoring function *score(s)* shown in Equation 1.

$$ score(s) = \frac{\sum_{i=1}^{n} score(s, i)}{n} \qquad (1) $$

Given a set of *n* instances, the *score(s)* of a MaxSAT solver *s* is the average of the scores for each instance, computed by *score(s, i)* in Equation 2.

$$ score(s, i) = \frac{1 + \text{best-known ub for instance } i}{1 + \text{ub for } i \text{ found by } s} \qquad (2) $$

The *score(s, i)* function computes the ratio between the *best-known upper bound* of an instance *i* and the bound reported by the solver *s* on the same instance. Assuming that (*best-known ub*) ≤ *ub* (which is the case for the MaxSAT Evaluation), the computed value ranges between 0 and 1, where higher values correspond to better upper bounds.

The competition has some specific rules about what is and is not allowed in the implementation of the solvers. In particular, the solvers are not allowed to employ triggers to modify their behavior, which is deemed to be specific to particular instances. However, solvers can concatenate the usage of different solving techniques.

In the most recent competition, the MaxSAT solvers used a variety of strategies and solvers. Some of these solvers are outlined below, along with the various approaches they employ in order to find improved solutions.

1. **NuWLS-c:** This solver adopts two solvers, the NuWLS solver, which is an improvement of SATLike, and the integration of TT-Open-WBO-Inc.

2. **TT-Open-WBO-Inc:** This solver uses four different strategies, including SATLike for inprocessing, a modified version of Mrs. Beaver for unweighted instances, BMO-clustering for weighted instances, and Polosat, a SAT-based local search method. This MaxSAT solver has three different distributions:
7:4 Exploiting Configurations of MaxSAT Solvers

■ (g) Which incorporates the Glucose 4.1 SAT Solver.

- (i) Which incorporates the new Intel SAT Solver.

- (is) Which incorporates the new Intel SAT Solver and is tuned for short invocations.

3. **DT-HyWalk:** This solver employs three distinct strategies, including a direct call to a SatSolver, the SATLike solver for local search, and the use of another MaxSAT solver, TT-Open-WBO-Inc.

4. **Loandra:** This solver utilizes two core algorithms, namely a Core-Based algorithm and a Linear algorithm.

Table 1, column “MSE”, shows the results of the MaxSAT Evaluation 2022 for the top six solvers at the incomplete weighted track with 60 seconds timeout. As we can see, the MaxSAT solver Loandra was not competitive within this category. In this paper, we propose an approach that is agnostic of the structure of instances and only allows the usage of alternative configurations of the same input MaxSAT solver. We experimentally show the goodness of our approach on the MaxSAT solver Loandra.

**3.1 Reproducing the MaxSAT Evaluation for the Incomplete track**

All of our executions of the MaxSAT solvers are run on a computation cluster composed of nodes with two AMD 7402 processors (each with 24 cores at 2.8 GHz) and 21 GB of RAM per core, managed by Sun Grid Engine (SGE). All the experiments are managed using the *Running Module* of the OptiLog framework.

Each execution is given 60s of CPU time and 32 GB of memory. As the memory requirements exceed the memory per core available, two slots are reserved and an affinity mask is set by SGE to restrict the execution to only one of the two cores. In contrast to the MaxSAT evaluation, each solver was evaluated with 50 different random seeds and we report results on the average score, and in some of the experiments, we also show the minimum and maximum scores, and the standard deviation.

In the course of developing our experiments, we detected two problems with some executions of the solvers: 1) some executions report a bound that does not correspond to the real cost of the solution reported, and 2) some executions report a solution that does not satisfy the *hard* clauses.

To address these issues we conduct a validation step executed after the solver exhausts the 60s of CPU time. In particular:

For 1), we trust the cost we compute from the solution reported, ignoring the bound reported by the solver.

For 2), we consider the solver was not able to find any solution at all.

This validation step is also conducted during the automatic configuration process when we evaluate a particular configuration of the solver on a given instance (see Section 5).

The score for each solver is computed using the MaxSAT evaluation rules. In particular, it is important to define which is the set of best-known upper bounds that we use to compute the score. Table 1 shows in column “MSE 2022”, the scores reported in the MaxSAT Evaluation 2022.

The rest of the columns present the results of the experimentation we conducted (in our cluster) using different sets of *best-known upper bounds*. “*VBS*<sub>*b*</sub>” uses the upper bounds found by the Virtual Best Solver of the solvers we executed, “*MSE*<sub>*b*</sub>” uses the set of best-known upper bounds provided by the MaxSAT Evaluation, and “*LRUNS*<sub>*b*</sub>” (Long Runs) uses a set of new *best-known upper bounds* we computed by running Loandra and NuWLS-c (both with the default parameters) with a timeout of 12 hours. We recall that the score presented
is the average score on 50 seeds in contrast to the MSE results where only 1 seed is used. As we can observe in Table 1, the relative ranking of the solvers is preserved although we can observe variations in the scores reported. We will use in the rest of the paper the best bounds from $MSE_b + VBS_b + LRUNS_b$.

<table><caption>**Table 1** Results of the MaxSAT Evaluation 2022 on the incomplete weighted track (60 seconds timeout) and reproduction of the Evaluation in our system with different sets of Best-Known Upper Bounds.</caption><thead><tr><th></th><th>MSE 2022</th><th colspan="3">MSE 2022 on our system</th></tr><tr><th>Best-known UBs</th><th>MSE<sub>b</sub></th><th>VBS<sub>b</sub></th><th>VBS<sub>b</sub> + MSE<sub>b</sub></th><th>VBS<sub>b</sub> + MSE<sub>b</sub> + LRUNS<sub>b</sub></th></tr></thead><tbody><tr><th>Solvers</th><td></td><td></td><td></td><td></td></tr><tr><th>NuWLS-c</th><td>0.759</td><td>0.7831</td><td>0.7590</td><td>0.7524</td></tr><tr><th>DT-Hywalk</th><td>0.732</td><td>0.7625</td><td>0.7414</td><td>0.7351</td></tr><tr><th>TT-Open-WBO-inc (g)</th><td>0.728</td><td>0.7412</td><td>0.7221</td><td>0.7164</td></tr><tr><th>TT-Open-WBO-inc (is)</th><td>0.726</td><td>0.7354</td><td>0.7201</td><td>0.7141</td></tr><tr><th>TT-Open-WBO-inc (i)</th><td>0.720</td><td>0.7354</td><td>0.7178</td><td>0.7118</td></tr><tr><th>Loandra</th><td>0.693</td><td>0.7107</td><td>0.7003</td><td>0.6953</td></tr></tbody></table>

# 4 Automatic Configurators (AC)

In this section, we review the Automatic Configuration Problem and two state-of-the-art automatic configuration algorithms or tuners.

## 4.1 The Automatic Configuration Problem

Given a target algorithm $A$ with parameters $\{p_1, \dots, p_n\}$ of domain $d(p_i)$. We define the parameter space $\Theta$ of $A$ as the subset $d(p_1) \times \dots \times d(p_n)$ of valid parameter combinations. Depending on the parameter, $d(p_i)$ can be categorical, a discrete domain of fixed values with no predefined order, or numerical, which represent integer or real values. Then, we define the *Automatic Algorithm Configuration* (AAC) problem as the optimization problem that consists of exploring $\Theta$ to find a configuration $\theta \in \Theta$ of $A$, which given a set of problem instances $\Pi$, minimizes a cost metric $\hat{c}: \Theta \times \Pi \to \mathbb{R}$, without exceeding a configuration budget $B$.

It is common for $A$ to be a black box (target algorithm), meaning it accepts some inputs (the parameters) and provides some output (e.g., $\hat{c}$), but we cannot see its internal functionality. This allows AAC to generalize to any type of algorithm but makes it more challenging for algorithm tuners since they cannot use $A$ to infer additional information about $\Theta$. In practice, $A$ is implemented as a binary file that outputs its results in a format suitable for its domain but may not be suitable for the AAC tool. Moreover, it may be necessary to limit the resources that $A$ can use to solve an instance, such as memory or CPU time. The standard way of addressing these issues in AAC tools is for the user to replace $A$ with a wrapper script that handles these and any other necessary aspects. Figure 1 describes the automatic configuration process where the tuner is a solver for the AAC problem.
<figure><img src="image_58.png" /><figcaption>Figure 1 Visualization of the Automatic Configuration process.</figcaption></figure>

## 4.2 The GGA Automatic Configurator

The Gender-Based Genetic Automatic Algorithm Configuration (GGA) is a genetic algorithm that was introduced in [4] to search for high-quality configurations. It was one of the pioneering algorithms that supported continuous parameters and introduced the novel concept of gender to apply diverse selection pressures to the population's individuals.

**Algorithm 1 GGA.**

**Input:** Target Algorithm $A$, Parameter Space $\Theta$, Instances $\Pi$, Performance Metric $\hat{c}$, # MiniTournaments $N$, Configuration Budget $B$
1: **function** GGA($A, \Theta, \Pi, \hat{c}, N, B$)
2: &nbsp;&nbsp;&nbsp;&nbsp;pop $\leftarrow$ initPopulation($\Theta$)
3: &nbsp;&nbsp;&nbsp;&nbsp;$j = 0$
4: &nbsp;&nbsp;&nbsp;&nbsp;**while** $B$ not exhausted **and** threshold not achieved **do**
5: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$j = j + 1$
6: &nbsp;&nbsp;&nbsp;&nbsp;$\Pi_j \leftarrow$ selectInstances($\Pi, j$)
7: &nbsp;&nbsp;&nbsp;&nbsp;$<w_1, ..., w_N>$ $\leftarrow$ runMiniTournaments($A, \text{pop.comp}, \Pi_j, \hat{c}, \text{pop.comp}/N$)
8: &nbsp;&nbsp;&nbsp;&nbsp;$\text{offspring} \leftarrow$ applyCrossoverAndMutate($\text{pop.noncomp}, <w_1, ..., w_N>, \Theta$)
9: &nbsp;&nbsp;&nbsp;&nbsp;pop $\leftarrow$ agingAndDeath($w_1, \text{pop}) \cup \text{offspring}$
&nbsp;&nbsp;&nbsp;&nbsp;**return** $w_1$

Algorithm 1 shows the pseudocode of the GGA algorithm, which takes as input the target algorithm $A$, its parameter space $\Theta$, a set of training instances $\Pi$, a performance metric $\hat{c}$ to optimize (e.g., time, accuracy, quality within a fixed timeout, etc), the number $N$ of GGA mini-tournaments (which will be explained shortly), and a configuration time budget $B$.

GGA starts by initializing a *population* (*pop*) of configurations (named *genomes*) as a subset of $\Theta$ in line 2. This population is partitioned into a *competitive* group (*pop.comp*, which is directly evaluated on the target algorithm) and *non-competitive* group (*pop.noncomp*, which simply acts as a source of diversity).

The algorithm proceeds in a main loop that finishes when GGA reaches the configuration budget $B$ or a threshold on the performance (line 4). At each iteration (which we call *generation*), GGA selects a subset of the instances $\Pi_j$ to evaluate the genomes in line 6¹. Then, in line 7, GGA evaluates the competitive genomes of the population over the selected instances $\Pi_j$ using a parallel racing scheme called *mini-tournament*. This procedure returns

<footnote n="1">¹ There are different policies that can be applied to select the instances at each generation, see [4].</footnote>
a set of N winners, $<w_1, ..., w_N>$, which will be the only competitive genomes that will generate new offspring in this generation (line 8). Finally, GGA applies an ageing policy in line 9 that is used to prevent population growth. The only exception is the overall best competitive genome ($w_1$), which survives as long as it performs better than the other mini-tournament winners. At the end of the main loop, GGA returns the best competitive genome $w_1$ of the last generation.

For more details on the GGA algorithm, we refer the reader to [4].

## 4.3 The SMAC Automatic Configurator

Sequential Model-Based Algorithm Configuration (SMAC) is an automatic configuration algorithm based on Bayesian optimization [10, 14]. In Bayesian optimization, we use a few evaluations of the target algorithm to train a *surrogate model* that predicts the performance of the algorithm for a given configuration. This fast-to-evaluate surrogate model is used to search for promising new configurations that will be executed on the training instances.

■ Algorithm 2 SMAC.

**Input:** Target Algorithm A, Parameter Space $\Theta$, Instances $\Pi$, Performance Metric $\hat{c}$, Configuration Budget $B$
1: **function** SMAC(A, $\Theta$, $\Pi$, $\hat{c}$, B)
2: &nbsp;&nbsp;&nbsp;&nbsp;[$R, \theta_{inc}$] $\leftarrow$ initialize($\Theta$, $\Pi$)
3: &nbsp;&nbsp;&nbsp;&nbsp;**while** B not exhausted **do**
4: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[$M, t_{fit}$] $\leftarrow$ fitModel(R)
5: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[$\tilde{\Theta}_{new}, t_{select}$] $\leftarrow$ selectConfigurations($M, \theta_{inc}, \Theta$)
6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[$R, \theta_{inc}$] $\leftarrow$ intensify(A, $\tilde{\Theta}_{new}, \theta_{inc}, R, \Pi, \hat{c}$)
&nbsp;&nbsp;&nbsp;&nbsp;**return** $\theta_{inc}$

Algorithm 2 shows the pseudocode of SMAC. This algorithm receives as input the target algorithm A, its parameter space $\Theta$, a set of training instances $\Pi$, a performance metric $\hat{c}$ to optimize and a configuration time budget B. First, SMAC initializes a *best candidate* configuration $\theta_{inc}$ and the history of conducted evaluations of different (*configuration, instance*) pairs R (which might be empty) in line 2.

As in GGA (see Section 4.2), SMAC has a main loop defined in line 3 that proceeds until the configuration budget $B$ is reached. At each iteration, it fits a surrogate model $M$ using the information in $R$ in line 4. Then, it uses $M$ to select a new set of promising candidate configurations $\tilde{\Theta}_{new}$ in line 5. Finally, it evaluates $\tilde{\Theta}_{new}$ and $\theta_{inc}$ on instances from $\Pi$ to determine the next best candidate $\theta_{inc}$, according to $\hat{c}$ in line 6. Similar to the GGA algorithm, SMAC returns the best candidate configuration $\theta_{inc}$.

## 4.4 Support for Tuning into the OptiLog framework

In this section, we present an excerpt of the code that uses the OptiLog framework to generate the configuration environment (from now on Tuning Scenario) of the solver Loandra for GGA and SMAC tuners.

1 \# example\_ac.py

2

3 from optilog.blackbox import *

4 from optilog running import ParsingInfo

5 from optilog.tuning import *

6
## 7:8 Exploiting Configurations of MaxSAT Solvers

```python
class LoandraBB(SystemBlackBox):
    config = {
        "weight-strategy": Int(0, 2, default=2),
        "preprocess": Bool(default=True),
        (...)
    }
    (...)
```

### Listing 1 Sample code to wrap the solver Loandra into an OptiLog BlackBox.

Listing 1 defines a custom BlackBox class named LoandraBB that inherits from System-BlackBox. This class represents the binary that we want to optimize. The `config` dictionary defines the parameters of this binary that can be tuned by the optimization algorithm. We show the parameters "weight-strategy" and "preprocess", with their respective types and default values. In particular, we need 18 lines of code to wrap Loandra, with 40 additional lines defining the parameters.

```python
1 # scenario_gga.py
2
3 from optilog.tuning.configurators import *
4 from example_ac import LoandraBB
5
6 if __name__ == "__main__":
7     configurator = GGAConfigurator(
8         LoandraBB(),
9         input_data="/path/to/instances/*",
10        ...)
11
12 configurator(generatescenario("./scenario"))
```

### Listing 2 Sample code to create a Tuning Scenario for the solver defined in Listing 1.

Listing 2 is a definition of a Tuning Scenario for the solver Loandra. It imports the custom `LoandraBB` class defined in Listing 1, and sets up a tuner (in this case GGA) to optimize the parameters of `LoandraBB`. The `input_data` parameter specifies the path to the instances used during the optimization process. Lastly, the `generatescenario` method is called with the desired output path for the scenario that is being created.

Similar code to Listing 2 could be used to define a Tuning Scenario to be used with the SMAC AC tool, as OptiLog supports both GGA and SMAC.

## 5 Configuring MaxSAT Solvers

Although the AC tools (tuners) presented in the previous section have also parameters that impact the effectiveness of the configuration process, *tuning the tuner* is out of reach in this paper and we focus on providing a good cost function to be used during the tuning process.

Ideally, we would use the `score(s, i)` function from the MaxSAT Evaluation (Equation 2). Notice though that in the MaxSAT Evaluation we are trying to maximize this scoring function, whereas tuners minimize a cost (see Section 4). Therefore, we have to convert the `score` function to a `cost` function. Additionally, it is not guaranteed that the bounds found are equal or worse than the previously *best-known upper bounds* (see Section 3), and we cannot update the *best-known upper bounds* sets during the tuning process (otherwise previous results computed in the same tuning process would not be comparable).
We define the cost<sub>ac</sub>(s, i) function, shown in Equation 3, as follows. First, we split the function in two cases: 1) the reported bound is worse (or equal) than the previous best-known upper bound, and 2) the bound reported is better.

For 1), we compute $1 - score(s, i)$ to obtain a value in the range [0, 1), where better bounds are closer to 0.

For 2), notice that we are breaking the assumption (best-known ub) $\le ub$, which may lead to unbounded values that tend to $\infty$. To restrict the values to the range (-1, 0), we use the inverse of the score(s, i) function, and then subtract 1.

The cost<sub>ac</sub>(s, i) function returns values between (-1, 1), where better bounds correspond to values closer to -1.

$$ cost_{ac}(s, i) = \begin{cases} 1 - score(s, i), & \text{if ub for } i \text{ found by } s \ge \text{best-known ub for } i \\ \frac{1}{score(s, i)} - 1, & \text{otherwise} \end{cases} \quad (3) $$

As stated in Section 3.1, some executions might report a bound that does not match the reported solution. Thus, we integrate a validation step (see Figure 2) that certifies the real cost of the solution returned by the solver and reports it to the tuner.

Regarding the tuning environment, all experiments are conducted on the same computation cluster. Each tuning process is given a wall-time tuning budget of 48 hours, a memory limit of 32G per worker, and is allowed to use up to 50 parallel workers unless otherwise specified. Each configuration instance is given a CPU time limit of 60 seconds for the solver, and then a validation step is executed. As training instances for the tuning process, we will use the 197 instances from the MaxSAT Evaluation 2021 (incomplete weighted track, 60 seconds) and we will test the best configuration returned by GGA and SMAC (see Section 4) on the 151 instances from the MaxSAT Evaluation 2022.

GGA allows for the selection of how many instances are used in each generation. Incrementally increasing the training set across generations till including the whole set of available instances is often recommended, as it facilitates discarding bad configurations with less effort, therefore more generations can be reached within the tuning budget. However, as discussed in Section 6, prioritizing the evaluation of more configurations on the whole training set within the same tuning budget may be preferable over having more generations. GGA also can preserve a set of elite configurations that are run at every generation. We define as an elite the default configuration of Loandra. Finally, we use the PyDGGA [2] (version 1.7.0) distribution of GGA which has support for distributed execution. The following non-default parameters were used for GGA: cost tolerance set to 0, population set to 100, generations set to 300, and minimum generations set to 50.

Regarding SMAC, although it can be executed in parallel, it does not report an overall winner in contrast to GGA. Instead, it reports as many winners as computation cores were used since it basically runs several sequential SMACs in parallel. Thus, after SMAC completes, we have to take all the winners from each SMAC sequential execution, which may not have been evaluated on all training instances, and perform the missing evaluations. Then, the winner with the best performance on the training set is selected to be the overall winner. We use the SMAC3 [14] (version 1.4.0) implementation of the algorithm.

As we have described earlier, the cost function used during the tuning process is not strictly the minimization version of the score we maximize according to the MSE 2022 (see Section 3). Therefore, one may argue that it would be better to return a winner for the training set with respect to the score function computed by the MSE. This is easy to do if the tuner provides the logs of each evaluation so, in the case of incomplete MaxSAT solvers, we can retrieve the best bound found by the solver on a given instance.
7:10 Exploiting Configurations of MaxSAT Solvers

Therefore, we add a selection phase (see Figure 2) after the tuning phase that recomputes the scores (according to the MSE) of the configurations traversed by the tool during the tuning process.

<figure><img src="image_59.png" /><figcaption>Figure 2 Visualization of the Automatic Configuration process extended with the validator and a selection phase over the explored configurations.</figcaption></figure>

In particular, for SMAC, we compute the MSE score of the 50 winners reported by the tuner on the training set and select the one with the highest score. Even though we have access to the logs of the evaluations of SMAC and could use those scores to select the winning configuration, we need to make sure that all the configurations are evaluated with all the instances.

For GGA, we order the configurations first by their ranking in a generation (according to the cost function in Equation 5), and within the same rank, we order by the most recent generation. Then, we select the first 50 distinct configurations. We look into their logs, recompute their MSE score according to Section 3, and report the winner<fnref n="2" />.

<table><caption>Table 2 Comparison using GGA and SMAC to tune the Loandra solver (using $VBS_b + MSE_b + LRUNS_b$ bounds).</caption><thead><tr><th></th><th>Mean</th><th>Median</th><th>Min</th><th>Max</th><th>Std</th></tr></thead><tbody><tr><td>NuWLS-c</td><td>0.7524</td><td>0.7522</td><td>0.7484</td><td>0.7560</td><td>0.0017</td></tr><tr><td>Loa (GGA, all-i)</td><td>0.7393</td><td>0.7391</td><td>0.7313</td><td>0.7475</td><td>0.0037</td></tr><tr><td>Loa (GGA, incremental)</td><td>0.7353</td><td>0.7354</td><td>0.7275</td><td>0.7433</td><td>0.0038</td></tr><tr><td>DT-Hywalk</td><td>0.7351</td><td>0.7355</td><td>0.7288</td><td>0.7415</td><td>0.0030</td></tr><tr><td>Loa (SMAC)</td><td>0.7237</td><td>0.7234</td><td>0.7149</td><td>0.7355</td><td>0.0048</td></tr><tr><td>TT-Open-WBO-inc (g)</td><td>0.7164</td><td>0.7165</td><td>0.7128</td><td>0.7194</td><td>0.0015</td></tr><tr><td>TT-Open-WBO-inc (i)</td><td>0.7141</td><td>0.7142</td><td>0.7093</td><td>0.7188</td><td>0.0020</td></tr><tr><td>TT-Open-WBO-inc (is)</td><td>0.7118</td><td>0.7117</td><td>0.7098</td><td>0.7145</td><td>0.0008</td></tr><tr><td>Loandra</td><td>0.6953</td><td>0.6957</td><td>0.6872</td><td>0.7036</td><td>0.0037</td></tr></tbody></table>

Table 2 shows the result of the best configurations provided by GGA using all the instances from the first generation (“Loa (GGA, all-i)”), or adding them incrementally at each step (“Loa (GGA, incremental)”) and the best configuration provided by SMAC (see

<footnote n="2">² In our experiments, the winner reported by GGA was the same configuration as the best one found in the selection phase.</footnote>
"Loa (SMAC)") after the additional selections process described in the paragraph above. For the incremental approach of GGA, we use 20% of instances at the first generation and instruct GGA to use all the instances on generation 25. Those values were selected based on preliminary experiments taking into account the number of generations that GGA can do in the given time. It is clear by the results that the usage of all the instances from the beginning benefits GGA, allowing it to lift Loandra from the sixth position to the second one. In the next section, we will focus on the variant of GGA "Loa (GGA, all-i)".

We also ran the variants “Loa (GGA, all-i)” and “Loandra” on another set of benchmarks, the MSE 2020 instances. The default parameters variant achieves a score of 0.755, and the tuned version a score of 0.777.

# 6 Exploiting Configurations Discarded by the Tuner

As it has been shown in the literature [9, 17], from the most pragmatic point of view, we can obtain an efficient parallel approach by just running the same non-deterministic solver with different seeds in parallel, or we can also run in parallel different configurations of the same solver.

In case resources are limited, we can also schedule the execution of different configurations of the same solver. In this section, we concrete and study these different approaches. We use OptiLog [1] to generate all the portfolios, as we explain in Section 6.3.

## 6.1 Parallel Portfolios of seeds and configurations

As we have already explained, tuners report the best configuration they have found. However, many other potentially good configurations are also explored and discarded during the automatic configuration process with respect to their performance on the particular training set. These configurations may exhibit good performance in different kinds of instances. As observed in [3] on SAT benchmarks, superior performance can be achieved by combining these complementary configurations.

The first approach we explored is the parallel execution of *N* different random seeds over a given MaxSAT solver. This approach can be applied to both the *default* MaxSAT solver and the best configuration obtained in the tuner.

Another approach is to extract *N* configurations of a MaxSAT solver from the ones traversed by the tuner and execute them in parallel. There are many strategies that we could follow to extract these configurations from the tuner. In particular, we use the set of configurations considered during the *selection phase* (after the *tuning phase*) process as explained in Section 5. Notice that we do not analyze any structure of the instances and we only incorporate configurations of the same solver.

Table 3 shows the results of the parallel portfolios that we explained. We tested parallel portfolios with 25, 30, 35, 40, 45, and 50 parallel executions. Each row shows the results of a parallel portfolio (rows marked with (Seeds) refer to a parallel portfolio of seeds, whereas the row marked with (Configs) refer to a parallel portfolio of configurations). We show the score as computed in the MaxSAT evaluation using the $VBS_b + MSE_b + LRUNS_b$ upper bounds and the rank of each portfolio with respect to the others.

As we can observe, the portfolio over different seeds for the default Loandra ("(Seeds) Loandra" in Table 3) is not competitive while the portfolio of different seeds for the best configuration of Loandra computed by GGA (column "(Seeds) Loa (GGA, all-i)") already outperforms NuWLS-c. Additionally, a portfolio of the best configurations provided by the selection phase (column "(Configs) Loa (GGA, all-i)") systematically outperforms the rest of the approaches. These observations hold almost for any number of parallel executions.
7:12 Exploiting Configurations of MaxSAT Solvers

<table><caption>Table 3 Score and rank (#) for each parallel portfolio, given N parallel processes (using $VBS_b + MSE_b + LRUNS_b$ bounds).</caption><thead><tr><th rowspan="2">N</th><th colspan="2">25</th><th colspan="2">30</th><th colspan="2">35</th><th colspan="2">40</th><th colspan="2">45</th><th colspan="2">50</th></tr><tr><th>#</th><th>score</th><th>#</th><th>score</th><th>#</th><th>score</th><th>#</th><th>score</th><th>#</th><th>score</th><th>#</th><th>score</th></tr></thead><tbody><tr><td>(Configs) Loa (GGA, all-i)</td><td>1</td><td>0.813592</td><td>1</td><td>0.818663</td><td>1</td><td>0.823986</td><td>1</td><td>0.825208</td><td>1</td><td>0.826448</td><td>1</td><td>0.827105</td></tr><tr><td>(Seeds) Loa (GGA, all-i)</td><td>2</td><td>0.806889</td><td>2</td><td>0.808078</td><td>2</td><td>0.809922</td><td>2</td><td>0.812589</td><td>2</td><td>0.813007</td><td>2</td><td>0.815263</td></tr><tr><td>(Seeds) NuWLS-c</td><td>3</td><td>0.768409</td><td>3</td><td>0.769302</td><td>3</td><td>0.769785</td><td>3</td><td>0.770293</td><td>3</td><td>0.771219</td><td>3</td><td>0.771263</td></tr><tr><td>(Seeds) DT-Hywalk</td><td>4</td><td>0.759271</td><td>4</td><td>0.764688</td><td>4</td><td>0.765736</td><td>4</td><td>0.765764</td><td>4</td><td>0.765862</td><td>4</td><td>0.766397</td></tr><tr><td>(Seeds) TT-Open-WBO-inc (g)</td><td>5</td><td>0.728901</td><td>5</td><td>0.728903</td><td>5</td><td>0.729766</td><td>5</td><td>0.729902</td><td>5</td><td>0.730165</td><td>5</td><td>0.730235</td></tr><tr><td>(Seeds) TT-Open-WBO-inc (i)</td><td>6</td><td>0.725971</td><td>6</td><td>0.726171</td><td>6</td><td>0.726253</td><td>6</td><td>0.726315</td><td>6</td><td>0.726510</td><td>6</td><td>0.727468</td></tr><tr><td>(Seeds) Loandra</td><td>8</td><td>0.717788</td><td>8</td><td>0.722663</td><td>7</td><td>0.723954</td><td>7</td><td>0.723996</td><td>7</td><td>0.724521</td><td>7</td><td>0.724648</td></tr><tr><td>(Seeds) TT-Open-WBO-inc (is)</td><td>7</td><td>0.722672</td><td>7</td><td>0.722704</td><td>8</td><td>0.722918</td><td>8</td><td>0.723174</td><td>8</td><td>0.723324</td><td>8</td><td>0.723489</td></tr></tbody></table>

## 6.2 Sequential Portfolios of configurations

In some settings, we will not have enough resources to run a parallel portfolio as described in Section 6.1. Potentially, we can have just one computation core available. In this case, we can schedule the sequential execution of different configurations of Loandra within the given timeout.

Let us describe how we construct this sequential portfolio. We assume we have a sequence of solvers (or configurations of a solver) (S) that iteratively report better solutions, a time budget (TO), and a maximum time budget a solver can exhaust between two consecutive reported solutions (MTBS). The solvers are executed according to their order in the sequence until the time consumed globally by all the solvers exceeds TO.

Each solver is run as follows: first, we wait for the first solution reported by the solver. Once this first solution is reported, we start a timer of MTBS seconds. If the solver reports a new solution before this timer expires, we reset the timer and wait for a new solution. This is repeated until the solver is unable to report a new solution before the timer is consumed. At that point in time, the solver is stopped and the next one in the ordered list of solvers is executed. Note that, at any point in this process, a solver can also be stopped if the global time budget of TO seconds gets exhausted. A special case is the last solver of the sequence, which is allowed to run until the time budget expires (i.e. it is not stopped even if it took more than MTBS seconds to find a new solution). Obviously, we keep track of the best overall solution seen so far.

To identify which sequence of solvers S and MTBS value the portfolio should use, we carry out a simulation of sequential portfolios with the configurations provided by the selection phase (see Section 5) and their respective logs on the training instances computed during the tuning phase. In particular, we explore all sequences of up to size 3 and MTBS values of {2, 3, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50} seconds. Once we identify the best virtual sequential portfolio for the training instances, we simulate again the execution of this virtual sequential portfolio on the test set. In Table 4 we present the results of this simulation.

To implement this virtual sequential portfolio we would need to take into account an additional thread that keeps track of the evolution of the solvers in the sequence, which may decrease the overall performance. Therefore, we see this virtual sequential portfolio as a restarting policy that MaxSAT developers could integrate into their solvers, with the added benefit that they may be able to reuse information computed by each solver in the sequence.

Table 4 shows the results of the virtual sequential portfolios (rows prefixed with “Virtual portfolio”), compared to the results that obtained the solvers from the competition with the default parameters, and with the best approach obtained using a tuner (“Loa (GGA, all-i)”). As in the MSE we run each solver with the same seed, except for NuWLS-c for which we also
<table><caption>Table 4 Score of the *virtual* sequential portfolio compared with the single-execution approach (using $VBS_b + MSE_b + LRUNS_b$ bounds).</caption><thead><tr><th></th><th>Score</th></tr></thead><tbody><tr><td>Virtual sequential portfolio (N=2) - Loa (gga, all-i)</td><td>0.7642</td></tr><tr><td>Virtual sequential portfolio (N=3) - Loa (gga, all-i)</td><td>0.7642</td></tr><tr><td>NuWLS-c (max score on 50 seeds)</td><td>0.7560</td></tr><tr><td>NuWLS-c</td><td>0.7554</td></tr><tr><td>Loa (gga, all-i)</td><td>0.7513</td></tr><tr><td>DT-Hywalk</td><td>0.7432</td></tr><tr><td>TT-Open-WBO-inc (g)</td><td>0.7214</td></tr><tr><td>TT-Open-WBO-inc (i)</td><td>0.7180</td></tr><tr><td>TT-Open-WBO-inc (is)</td><td>0.7180</td></tr><tr><td>Loandra</td><td>0.6965</td></tr></tbody></table>

report on the best score value from 50 seeds. The *N* value shown in the *virtual* sequential portfolios rows indicates the length of the solvers' sequence. The portfolios are built on top of the configurations obtained after the selection phase with ("Loa (GGA, all-i)").

We notice that *virtual* sequential portfolios do perform better than NuWLS-c, and a selection of two configurations suffices to that end. Interestingly, if we build the *virtual* sequential portfolio on the test instances from the MSE 2022, then we get a better portfolio using three configurations that achieves a score of 0.7689, however, we cannot predict this portfolio based on the analysis we perform on the training instances from the MSE 2021.

Additionally, we conducted experiments to analyze the potential of combining solvers. In particular, we used the solver TT-Open-WBO<fnref n="3" />. We tuned this solver, and generated a *virtual* sequential portfolio combining the best configurations of TT-Open-WBO and Loandra. The *virtual* sequential portfolio (based on the analysis of the performance on the MSE 2021) obtained a score of 0.779 on the MSE 2022, which is the best score obtained for single-core evaluations. In comparison, the individual performance of Loandra and TT-Open-WBO after tuning is 0.747 and 0.751 respectively.

## 6.3 OptiLog Portfolio Generator

To facilitate the creation of the parallel and *virtual* sequential portfolios, we added support to compute them using the OptiLog framework.

```python
1 from optilog.portfolio import get_parallel_portfolio
2
3 get_parallel_portfolio(
4     ggaScenario='./tuning-scenario',
5     n_solvers=10,
6     save_to='./parallel-portfolio'
7)
```

<caption>Listing 3 Computing a parallel portfolio with OptiLog.</caption>

<footnote n="3">³ A version provided by the author of the solver with the parameters exposed.</footnote>
7:14 Exploiting Configurations of MaxSAT Solvers

Listing 3 shows how we can generate and save a parallel portfolio with OptiLog. This portfolio is built by selecting *N* configurations as explained in Section 5, thus requiring a Tuning Scenario (generated with OptiLog as seen in Section 6.3). The function `get_parallel_portfolio` receives as parameters the Tuning Scenario that contains the results of the tuning process (`ggaScenario`), the number of solvers that will compose the parallel portfolio (`nSolvers`), and the directory where the scripts to launch each individual solver that composes the portfolio will be saved (`save_to`).

```python
1 from optilog.portfolio import get_sequential_portfolio
2
3 get_sequential_portfolio(
4     pathScenario="./running-scenario",
5     nSolvers=2,
6     solution matchers=[r"^o\s(\d+)",
7     save_to="./sequential-portfolio",
8     score_fn=maxsat_score_fn,
9     max_time-between_solutions=[5, 15, 25, 35]
10 )
```

<caption>Listings 4 Computing a sequential portfolio with OptiLog.</caption>

Listing 4 shows how we are generating a sequential portfolio with the results of a Running Scenario. Note that to generate the *virtual* sequential portfolio we require the full trace of the solvers (in particular for the incomplete MaxSAT case, we need the evolution of the best bound over time), so we cannot build it from a Tuning Scenario directly. The parameters `ggaScenario`, `save_to`, and `nSolvers` mean the same as in the function to compute a parallel portfolio. Additionally, we have to specify the following parameters: `score_fn` is used to transform the lines matched by `solution matchers` to a score that the portfolio will try to maximize (in this example the score function is `score(s)` defined in Section 3), and `max_time-between_solutions` contains the possible values that the portfolio can choose from when selecting the parameter *MTBS*.

7 Conclusions

Given a target solver, we have presented an approach to easily generate a potentially much better solving approach. To this end, we exploit a set of alternative configurations of the same target solver coming from the residues of a tuning process. It is important to notice that we do not exploit any structure feature of the input problem or instance since in some domains these features are not easy to compute. In particular, we have shown how from a MaxSAT solver with a low ranking in one of the tracks of the MSE 2022 we can obtain a more competitive approach.

Our sequential portfolio generation approach can be seen as a first attempt to come up with effective restarting policies for MaxSAT solvers, something that has not been studied in depth in the literature.

Finally, the approach described has been integrated into the OptiLog framework avoiding the tedious process of setting up tuning environments and generating portfolios. Moreover, the API is general enough to be applied not only to MaxSAT solvers but to other solving approaches.
References

1. Josep Alòs, Carlos Ansótegui, Josep M. Salvia, and Eduard Torres. OptiLog V2: Model Solve, Tune and Run. In Kuldeep S. Meel and Ofer Strichman, editors, 25th International Conference on Theory and Applications of Satisfiability Testing (SAT 2022), volume 236 of Leibniz International Proceedings in Informatics (LIPIcs), pages 25:1–25:16, Dagstuhl, Germany, 2022. Schloss Dagstuhl – Leibniz-Zentrum für Informatik. ISSN: 1868-8969. doi: 10.4230/LIPIcs.SAT.2022.25.

2. Carlos Ansótegui, Josep Pon, and Meinolf Sellmann. Boosting evolutionary algorithm configuration. Annals of Mathematics and Artificial Intelligence, 2021. doi:10.1007/s10472-020-09726-y.

3. Carlos Ansótegui, Josep Pon, and Meinolf Sellmann. Boosting evolutionary algorithm configuration. Annals of Mathematics and Artificial Intelligence, 90(7-9):715–734, 2022. Publisher: Springer.

4. Carlos Ansótegui, Meinolf Sellmann, and Kevin Tierney. A Gender-based Genetic Algorithm for the Automatic Configuration of Algorithms. In Proceedings of the 15th International Conference on Principles and Practice of Constraint Programming, CP’09, pages 142–157. Springer-Verlag, 2009. tex.acmid: 1789011 tex.numpages: 16 tex.year: 2009 event-place: Lisbon, Portugal. URL: http://dl.acm.org/citation.cfm?id=1788994.1789011.

5. Fahiem Bacchus, Jeremias Berg, Matti Järvisalo, Ruben Martins, and Andreas Niskanen. MaxSAT Evaluation 2022: Solver and Benchmark Descriptions. Department of Computer Science Series of Publications B, B-2022-2, 2022. Accepted: 2022-08-25T10:09:01Z Publisher: Department of Computer Science, University of Helsinki. URL: https://helda.helsinki.fi/handle/10138/347396.

6. Jeremias Berg, Emir Demirovic, and Peter Stuckey. Core-boosted linear search for incomplete maxsat. Integration of Constraint Programming, Artificial Intelligence, and Operations Research: 16th International Conference, CPAIOR 2019, 2019.

7. Yi Chu, Shaowei Cai, Zhendong Lei, and Xiang He. Nuwls-c: Solver description. MaxSAT Evaluation 2022, page 28, 2022.

8. Jan Elffers and Jakob Nordström. Divide and conquer: Towards faster pseudo-boolean solving. Proceedings of the Twenty-Seventh International Joint Conference on Artificial Intelligence, pages 1291–1299, 2018. URL: https://www.ijcai.org/Proceedings/2018/180.

9. Youssef Hamadi, Said Jabbour, and Lakhdar Sais. ManySAT: a parallel SAT solver. JSAT, 6:245–262, June 2009. doi:10.3233/SAT190070.

10. Frank Hutter, Holger H Hoos, and Kevin Leyton-Brown. Sequential Model-Based Optimization for General Algorithm Configuration (extended version). International Conference on Learning and Intelligent Optimization, 2011.

11. Jo Devriendt. Exact Solver Repository, April 2023. URL: https://gitlab.com/JoD/exact.

12. Saurabh Joshi, Prateek Kumar, Sukrut Rao, and Ruben Martins. Open-wbo-inc: Approximation strategies for incomplete weighted maxsat. J. Satisf. Boolean Model. Comput., 11(1):73–97, 2019. doi:10.3233/SAT190118.

13. Serdar Kadioglu, Yuri Malitsky, Meinolf Sellmann, and Kevin Tierney. ISAC – Instance-Specific Algorithm Configuration. In ECAI 2010, pages 751–756. IOS Press, 2010. doi: 10.3233/978-1-60750-606-5-751.

14. Marius Lindauer, Katharina Eggensperger, Matthias Feurer, André Biedenkapp, Difan Deng, Carolin Benjamins, Tim Ruhkopf, René Sass, and Frank Hutter. SMAC3: A Versatile Bayesian Optimization Package for Hyperparameter Optimization. In ArXiv: 2109.09831, 2021. URL: https://arxiv.org/abs/2109.09831.

15. Ole Lübke and Sibylle Schupp. nosat-maxsat. In MaxSAT Evaluation 2022, pages 29–30. Department of Computer Science, University of Helsinki, 2022.

16. Alexander Nadel. Polarity and Variable Selection Heuristics for SAT-Based Anytime MaxSAT: System Description. Journal on Satisfiability, Boolean Modeling and Computation, 12(1):17–22, September 2020. doi:10.3233/SAT-200126.
7:16 Exploiting Configurations of MaxSAT Solvers

17 Olivier Roussel. Description of ppfolio 2012. In *Proceedings of SAT Challenge 2012: Solver and Benchmark Descriptions*, 2012.

18 Jiongzhi Zheng, Kun He, Zhuo Chen, Jianrong Zhou, and Chu-Min Li. Decision tree based hybrid walking strategies. *MaxSAT Evaluation 2022*, page 24, 2022.
# Symmetries for Cube-And-Conquer in Finite Model Finding

João Araújo

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

ID

Universidade Nova de Lisboa, Lisbon, Portugal

Choiwah Chow

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

Universidade Aberta, Lisbon, Portugal

Mikoláš Janota

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

ID

Czech Technical University in Prague, Czech Republic

## Abstract

The cube-and-conquer paradigm enables massive parallelization of SAT solvers, which has proven to be crucial in solving highly combinatorial problems. In this paper, we apply the paradigm in the context of finite model finding, where we show that isomorphic cubes can be discarded since they lead to isomorphic models. However, we are faced with the complication that a well-known technique, the Least Number Heuristic (LNH), already exists in finite model finders to effectively prune (some) isomorphic models from the search. Therefore, it needs to be shown that isomorphic cubes still can be discarded when the LNH is used. The presented ideas are incorporated into the finite model finder Mace4, where we demonstrate significant improvements in model enumeration.

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

**Keywords and phrases** finite model enumeration, cube-and-conquer, symmetry-breaking, parallel algorithm, least number heuristic

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

**Supplementary Material Software (Source Code):** https://github.com/ChoiwahChow/public/tree/main/CP_2023_Supplement_2.zip
archived at swh:1:cnt:69e6b145a05d29726512c0605cfb027c17c98dd3

**Funding** João Araújo: Fundação para a Ciência e a Tecnologia, through the projects UIDB/00297-2020 (CMA), PTDC/MAT-PUR/31174/2017, UIDB/04621/2020 and UIDP/04621/2020.

Mikoláš Janota: The results were supported by the Ministry of Education, Youth and Sports within the dedicated program ERC CZ under the project POSTMAN no. LL1902. This article is part of the RICAIP project that has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 857306.

## 1 Introduction

An important tool that working algebraists need in their research is libraries of the algebras they are interested in. These libraries allow them to get intuitions, test or refute hypotheses and conjectures, and gain insights into the properties of the algebras (see examples on p. 2891 of [30]). Many libraries of algebraic models of small orders, such as the smallsemi package [14] for semigroups and the loops package [36] for quasigroups, are available in the GAP [16] system. A lot more such libraries are needed, but they often take an inordinate amount of time and computing resources to generate.

First-order logic (FOL) has been the most popular language to define algebras. There are two major resource-intensive steps in generating non-isomorphic models from FOL [27]. The first step is to generate models according to the laws specified by the FOL formula. This step often generates a huge number of isomorphic models. For example, given the first-order

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

© João Araújo, Choiwah Chow, and Mikoláš Janota;
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. 8; pp. 8:1–8:19

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

LIPICS

Leibniz International Proceedings in Informatics

Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl Publishing, Germany
## 8:2 Symmetries for Cube-And-Conquer in Finite Model Finding

formula for semigroups, which is $(x * y) * z = x * (y * z)$, Mace4 [35] generates 1,021,120,198 models of order 7, out of which only 1,627,672 ($\approx$ 0.16%) [44] are pairwise non-isomorphic. The second step is to eliminate the isomorphic models generated in the first step. In this paper, we propose a novel efficient and scalable parallel algorithm that not only speeds up the first step but also generates fewer isomorphic models. Suppressing the generation of isomorphic models in the first step reduces the workloads of both the first and the second steps. Not only does it make the whole process much faster, but the required computing resources (disk space, etc.) are also reduced.

While modern-day general-purpose computers are predominantly multi-core, harnessing parallelism for combinatorial search is surprisingly difficult. Consequently, there are few parallel algorithms in constraints programming in general, and in finite model enumeration in particular. Indeed, in satisfiability modulo theories (SMT), even negative results are concluded for cube-and-conquer [23]. A recent literature review concludes that “there is little overall guidance that can be given on how best to exploit multi-core computers to speed up constraint solving” [18]. We aim to help close this gap by devising new parallel algorithms for finite model enumeration.

This paper advances finite model enumerators toward the following two objectives:

1. Mathematicians can use the tool to quickly generate all models (up to isomorphism) of the classes of algebras of their interests on their multi-core computers.

2. The tool can also take advantage of massively parallel computing architectures to pre-generate models (up to isomorphism) of the classes of algebras of general interest.

We find inspiration in the well-established cube-and-conquer approach introduced for SAT [20]. In SAT this means splitting the search space by mutually exclusive conjunctions of propositional literals (cubes). In the context of finite model finding, the structure is richer – a decision of the solver corresponds to inserting a point into the graph of one of the considered functions, e.g., $f(0, 1) = 3$. We comment on cube-and-conquer in more detail in Section 6.

We show that a cube can be excluded from the search if it is isomorphic to an existing one. Effectively, this is breaking symmetries in the search space. However, the task is nontrivial because finite model finders already contain a technique, called the least number heuristic (LNH), to exclude some isomorphic models. The LNH<fnref n="1" /> enables the solver to consider only certain values from the co-domain for a given decision point. Therefore, we show that isomorphic cubes can be pruned in the presence of the LNH. Like so, we can take advantage of the two powerful and complementary techniques and ultimately suppress the generation of a large number of isomorphic models.

This paper’s contributions are the following:

1. Devise a low runtime overhead parallel algorithm based on the cube-and-conquer approach for finite model enumeration. This scalable algorithm divides finite model enumeration into many independent non-overlapping search jobs to make full use of the available resources.

2. We show that isomorphic cubes can be discarded without losing isomorphic models even in the presence of the well established symmetry breaking technique already present in finite model finders – the least number heuristic (LNH).

3. We extend the model finder Mace4 with the proposed techniques and evaluate it on a large number of problems, where significant speed-up is observed.

<footnote n="1">¹ Despite the technique being called a heuristic, it does not sacrifice the completeness of the solver.</footnote>
## 2 Preliminaries

Familiarity with the general notions of abstract algebra such as groups, semigroups, and quasigroups is assumed, and so is general knowledge about functions and isomorphisms. A good reference is Chapters 2 and 5 of [9].

In this paper, the domain of the search space is denoted by the set $D = \{0, \dots, n-1\}$, where $n \ge 2$, That is, we exclude the trivial case of searching on domains of size 1.

Let $\pi$ denote an arbitrary permutation on $D$, $\pi_{id}$ denote the identity permutation, and $\pi_{(a,b)}$ denote the permutation cycle $(a, b)$. For example, $\pi_{(0,1)}$ is the permutation cycle $(0, 1)$.

### 2.1 Finite Model Enumeration

For a signature $\Sigma$ and a FOL formula $\mathcal{F}$ on $\Sigma$, a traditional finite model finder first expands the FOL formula to its ground representation by its domain elements in $D$, then searches for models by backtracking to exhaustively explore the search space [49]. The domain elements in $D$ are seen as special constants not appearing in the original $\mathcal{F}$, c.f. [40].

Following the terminology of [49], a *value assignment (VA) clause* is a term $f(a_1, \dots, a_k) = v$, where $f$ is a $k$-ary function symbol in $\Sigma$ and $a_j, v \in D$. We refer to the term $f(a_1, \dots, a_k)$ as the *cell term* (or simply *cell*) since conceptually the search fills the operation table of $f$.

To search for finite models in $\mathcal{F}$, the finite model finder employs a *cell selection* strategy to pick cell terms successively, without duplicates, to assign values from $D$ to form VA clauses. If a newly formed VA clause causes any failure in the axioms in $\mathcal{F}$, then a new value will be tried for that cell term. If no value can be assigned to that cell term without failing the axioms in $\mathcal{F}$, then the model finder backtracks to the previous cell term to try to assign another value to it. When all cell terms in $\mathcal{F}$ are assigned values without violating the axioms in $\mathcal{F}$, a model, as represented by its VA clauses, is found. After a model is found, the process can continue with backtracking to find more models.

A set of models can be partitioned into equivalence classes by isomorphisms. Intuitively, a model can be transformed into any other model in the same equivalence class by renaming its domain elements. Two models are said to be isomorphic to each other if an isomorphism exists from one model to the other.

The search space can be organized as a search tree in which nodes are VA clauses and edges join successive nodes with cell terms in the search order. The root node is an empty VA clause. The cell term in each node is selected by the cell selection strategy. A *search path* in a search tree is a path from the root to a node in the search tree. It can be represented by a sequence of VA clauses $\langle t_0 = v_0; t_1 = v_1; \dots \rangle$, where $t_i$ is the cell term in the $i^{th}$ position of the sequence and $v_i \in D$, and $t_i \neq t_j$ when $i \neq j$. Furthermore, a search path will be terminated at the first VA clause that results in a violation of any axiom of $\mathcal{F}$.

If the length of a search path is the same as the total number of cell terms in $\mathcal{F}$, then it is a complete search path and its VA clauses represent a model. Otherwise, it is an incomplete search path representing *partial assignments* of cell values in $\mathcal{F}$.

The backtracking algorithm in its simplest form is to try every possible value assignment for every cell. For example, to search an FOL formula $\mathcal{F}$ with just one binary operation, there are $n^{n^2}$ possible combinations ($n^2$ cells with $n$ possible values each). Even the very small domain size of 4 gives over 4 billion combinations of cell values. However, in practice, the number of viable combinations to check is much smaller than the theoretically maximum number because of the constraints imposed by $\mathcal{F}$. Furthermore, a finite model finder may infer new VA clauses from existing ones by *propagation*.
**Example 1.** Suppose the FOL formula contains only the equation $f(x, y) = f(y, x)$, that is, the operation *f* is commutative. After the assignment $f(0, 1) = 0$, the finite model finder can infer $f(1, 0) = 0$. This is referred to as positive propagation.

On the other hand, if the FOL formula contains the inequality $f(x,y) \neq f(y,x)$, then after the assignment $f(0,1) = 0$, the finite model finder can exclude 0 from the list of possible values for the cell $f(1,0)$. This is referred to as negative propagation. ◀

## 2.2 Least Number Heuristic

The least number heuristic (LNH) [4, 50, 51] is a very effective symmetry-breaking algorithm widely implemented in model finders/enumerators such as Mace4. The main idea of the LNH is that all domain elements that have not yet appeared in any VA clauses and the current cell term in the search are indistinguishable to each other and therefore only one of them, say, the smallest one, needs to be considered in a cell value assignment.

To ease discussions of the LNH, we introduce the notation $\text{Vals}(P)$ to denote the set of all domain elements appearing in $P$, where $P$ can be a search path, a VA clause, or a cell term.

**Example 2.** For the cell term $f(1, 1)$: $\text{Vals}(f(1, 1)) = \{1\}$. For the VA clause $f(1, 1) = 0$: $\text{Vals}(f(1, 1) = 0) = \{0, 1\}$. For the partial search path $S = \langle f(0, 0) = 0; f(1, 1) = 0 \rangle$: $\text{Vals}(S) = \{0, 1\}$. ◀

The LNH can now be stated precisely: In adding a VA clause, $t = v$, to extend a search path $S$, the possible choices of $v$ allowed under the LNH are $\text{Vals}(S) \cup \text{Vals}(t) \cup \{s\}$ where $s$ is the smallest domain element in $D \setminus (\text{Vals}(S) \cup \text{Vals}(t))$, and they are $D$ if $\text{Vals}(S) \cup \text{Vals}(t) = D$. Strictly speaking, it is not necessary to set $s$ to be the smallest domain element not seen so far, it could as well be the biggest one, for example. But the rule to set $s$ must be unambiguous - only one value is consistently picked by the rule each time. In this paper, we always set $s$ to be the smallest domain element not seen so far.

Furthermore, we say a search path is *LNH-compliant* if it respects the LNH restrictions on the choices of values assigned to its VA clauses.

**Example 3.** Suppose the domain size, $|D|$, is 4. Then the complete search path $\langle f(1) = 0; f(0) = 3; f(3) = 1; f(2) = 1 \rangle$ is not LNH-compliant.

For the first VA clause in the search path, $S = \emptyset$ and $t = f(1)$. So, $\text{Vals}(S) \cup \text{Vals}(t) = \emptyset \cup \{1\} = \{1\}$, and therefore $D \setminus (\text{Vals}(S) \cup \text{Vals}(t)) = \{0, 2, 3\}$. Thus, $s = \min(\{0, 2, 3\}) = 0$. The LNH limits the choices of the value for $f(1)$ to $\text{Vals}(S) \cup \text{Vals}(t) \cup \{s\} = \{0, 1\}$. So the first VA clause $f(1) = 0$ is LNH-compliant. However, for the second VA clause in the search path, $S = \{f(1) = 0\}$ and $t = f(0)$. So, $\text{Vals}(S) \cup \text{Vals}(t) = \{0, 1\} \cup \{0\} = \{0, 1\}$, and therefore $D \setminus (\text{Vals}(S) \cup \text{Vals}(t)) = \{2, 3\}$. Thus, $s = \min(\{2, 3\}) = 2$. The LNH limits the choices of the value for $f(0)$ to $\text{Vals}(S) \cup \text{Vals}(t) \cup \{s\} = \{0, 1, 2\}$, so $f(0) = 3$ is not allowed under the LNH. Therefore, the whole search path is not LNH-compliant. ◀

The LNH does not impose any restrictions on the order of the cell terms in the search path<fnref n="2" />. It speeds up the search by limiting the choices of the values for the cell terms. Therefore, its effectiveness decreases with the increase in the length of the search path as more domain elements are used when more VA clauses are added to the search path.

<footnote n="2">² In practice, a number called the maximal designated number (mdn) is often used to partition the domain into 2 subsets so that $\{0, ..., \text{mdn}\}$ are domain elements already seen, and $\{\text{mdn} + 1, ..., n - 1\}$ are domain elements not seen so far [49]. In this case, cell selection strategies that keep the mdn small are preferred because the search tree will be kept narrower.</footnote>
**► Example 4.** The concentric cell selection strategy is a simple cell selection strategy to minimize the growth of choices of values in the finite model search with the LNH. This strategy picks the cell $f(a_0, \dots, a_{k-1})$ with the least $r = \max(a_0, \dots, a_{k-1})$ from all available cells. Any fixed tie-breaker can be used in case of a tie. For example, one of the possible orders of the cells by this cell selection strategy for a binary operation is $f(a_0, a_1) < f(b_0, b_1)$ if $a_0 = a_1 \lor a_0 + a_1 < b_0 + b_1 \lor (a_0 + a_1 = b_0 + b_1 \land a_0 < b_0)$. This gives the sequence $f(0, 0), f(1, 1), f(0, 1), f(1, 0), f(2, 2), f(0, 2), f(2, 0), f(2, 1), f(1, 2), f(3, 3) \dots$ . ◀

## 2.3 Cube

A cube is a prefix of a search path, and as such, it can be specified by a sequence of VA clauses. Permutations and isomorphisms can be applied to a cube by applying them to its VA clauses. Specifically, if $\pi$ is a permutation on $D$ and $B$ is a cube, then $\pi(B) := \{f(\pi(a_1), \dots, \pi(a_k)) = \pi(v) \mid f(a_1, \dots, a_k) = v\}$ is a VA clause in $B$. Observe that $\pi_{id}(B)$ is the (unordered) set of all individual VA clauses in the cube $B$.

Note that predicates in an FOL formula can be implemented as functions with two values, $T$ (true) and $F$ (false), which do not affect the LNH because they are not domain elements. For convenience, we consider $I(T) = T$ and $I(F) = F$ for any isomorphism $I$ so that the same terminology is used for both relations and functions.

Cubes are said to be isomorphic if their VA clauses are isomorphic. In particular, two cubes $B_0$ and $B_1$ are isomorphic if there is a permutation $\pi$ on $D$ such that $\pi(B_0) = \pi_{id}(B_1)$.

**► Example 5.** If $B_0 = \langle f(0) = 0; g(0,0) = 0; f(1) = 0; g(1,1) = 0 \rangle$ and $B_1 = \langle f(0) = 1; g(0,0) = 1; f(1) = 1; g(1,1) = 1 \rangle$, then $B_0$ and $B_1$ are isomorphic because $\pi_{(0,1)}(B_0) = \{f(1) = 1, g(1,1) = 1, f(0) = 1, g(0,0) = 1\} = \pi_{id}(B_1)$. ◀

# 3 Isomorphic Cubes Redundancy

The main objective of this section is to show that isomorphic cubes can be removed from the search. More formally, if cubes $B_0$ and $B_1$ are isomorphic, then it is sufficient to explore assignments extending $B_0$ and ignore all assignments extending $B_1$. We need to prove that any model lost by discarding $B_1$ must necessarily be isomorphic to some model obtained from extending $B_0$ under the LNH. This statement is intuitive, but the proof requires some care as effectively, we are dealing with a combination of two symmetry-breaking techniques: LNH and isomorphic cube pruning, under an arbitrary search strategy.

As a motivational example, consider the cube $\langle f(0,0) = 0 \rangle$, which states that $f$ is idempotent in $0$. But because $0$ does not appear in the original FOL formula, intuitively, the constant $0$ cannot play a special role in the formula. Consequently, this cube searches all interpretations of $f$ that have at least one idempotent. For instance, the cube $\langle f(1,1) = 1 \rangle$ will search the same interpretations, up to isomorphism. Now, we need to show this property formally and that it holds when the solver searches with the LNH restriction.

The key idea of the proof is that given a model $B_1$ with VA clauses $A$, any cube that is isomorphic to a subset of $A$ can be gradually extended to be a model isomorphic to $B_1$. Each extension step of the cube must uphold the following properties: (1) The cube is isomorphic to some subset of $A$. (2) The cube is LNH-compliant. The extension step is illustrated in Figure 1. We are given a cube $B_0$ that is isomorphic to an $A_0 \subseteq A$. When the finite model finder decides on some empty cell $t$, we need to show that it is possible to find a value according to the LNH such that the extended cube is isomorphic to some subset of $A$ containing $A_0$.
<figure><img src="image_65.png" /><figcaption>Figure 1 Extension of a cube according to the VA clauses A.</figcaption></figure>

**Notation 1.** For a mapping $R$ from $D$ to $D$ and a value $d \in D$ we write $\mathcal{E}_R^d$ for a mapping that maps $d$ to $R(d)$ if $d \in \text{dom}(R)$ and otherwise maps $d$ to $\min(D \setminus \text{rang}(R))$. We further write $\mathcal{E}_R^{d_1, \dots, d_k}$ for successive extensions by $d_1, \dots, d_k$, i.e. $\mathcal{E}_R^{d_1, d_2} = \mathcal{E}_{\mathcal{E}_R^{d_1}}^{d_2}$ etc. ◀

**Example 6.** Suppose $D = \{1, 2, 3\}$ and $R : \{1\} \to \{2\}$ s.t. $R(1) = 2$ and $R^{-1}(2) = 1$ (so, $R$ is a bijection). Then $\mathcal{E}_R^2$ s.t. $R(2) = \mathcal{E}_R^2(2) = 1$ is a valid extension of $R$ because $\min(D \setminus \{2\}) = 1$. Furthermore, $\mathcal{E}_R^{2,1}$ s.t. $\mathcal{E}_R^{2,1}(1) = 2$ is a valid (but trivial) extension of $\mathcal{E}_R^2$. ◀

**Lemma 7.** If $R$ is a bijection between some $D_0, D_1 \subseteq D$ and $d \in D$ then $\mathcal{E}_R^d$ is well-defined and also a bijection.

**Proof.** If $d \in D_0$, then $\mathcal{E}_R^d = R$ and there is nothing to proof. If $d \in D \setminus D_0$, then by definition, $\mathcal{E}_R^d = R \cup \{(d, p)\}$ for some $p \in D \setminus D_1$. Since $R$ is a bijection from $D_0$ to $D_1$, $d \notin \text{dom}(R)$, and $p \notin \text{rang}(R)$, so $\mathcal{E}_R^d$ is well-defined, one-one, and onto. That is, it is a bijection. ◀

**Notation 2.** $B \oplus \langle t = u \rangle$ is the new cube formed by extending the cube $B$ with the VA clause $t = u$. ◀

The following lemma is the core of our proof. We have a cube $B$ isomorphic to some partial assignment $A_0$ and now we need to prove that for any model $A$ completing $A_0$ and any search strategy, we are able to extend $B$ while observing the LNH. Then, the lemma is used to prove that isomorphic cubes can be discarded by induction on cube length (Theorem 9).

**Lemma 8.** Let $B$ be an LNH-compliant cube and $A$ a model s.t. $B$ is isomorphic to some $A_0 \subseteq A$. Then for any cell term $t$ not appearing in $B$, there exists a value $u$ and a VA clause $t' = u' \in A \setminus A_0$, s.t. $B \oplus \langle t = u \rangle$ is LNH-compliant and isomorphic to $A_0 \cup \{t' = u'\}$.

**Proof.** Let $R$ be an isomorphism mapping $B$ to $A_0$ and let $t$ be a cell term $f(a_1, \dots, a_k)$. Define $R_1$ as $\mathcal{E}_R^{a_1, \dots, a_k}$, and let $t'$ denote the cell term $f(R_1(a_1), \dots, R_1(a_k))$, i.e., map the cell that the solver searches on into a cell in the prescribed model $A$.

Since $A$ is a model, there must exist a value $u' \in D$ with $(t' = u') \in A$, i.e. $u'$ can be found by a lookup of $t'$ in $A$. Since $t$ is not a cell term in $B$ and $R_1$ is a bijection, so $t'$ is not a cell term in $A_0$ and must therefore be in $A \setminus A_0$. Thus, $t' = u'$ is a VA clause in $A \setminus A_0$.

To obtain $u$ (a value for cell $t$), define $R_2$ as $\mathcal{E}_{R_1^{u'}}$, i.e. map $u'$ back into the search by extending the inverse. Then, set $u = R_2(u')$. By Lemma 7, $R_2$ is bijection and it is therefore an isomorphism from $A_0 \cup \{t' = u'\}$ to $B \oplus \langle t = u \rangle$. Finally, by definition of $R_2$, $u$ either already appears in $B$ or otherwise is the smallest domain element not in $B$. Therefore, the extension of the cube $B$ by the VA clause $t = u$ is LNH-compliant. ◀

**Theorem 9.** Suppose we are searching under the LNH with any cell selection strategy on a signature $\Sigma$ and a FOL formula, $\mathcal{F}$, on $\Sigma$. If $B_0$ and $B_1$, of length $l \ge 0$, are isomorphic cubes, and if $M_1$ is a model obtained by completing (not necessarily under the LNH) the search path in $B_1$, then $B_0$ can be extended by a search path $S$ under the said LNH and cell selection strategy to a model $M_0$ which is isomorphic to $M_1$.