Lecture 10: Solvers and How to Help Them

So far we have written models in pulp and called .solve(...), and the answer came back within a fraction of a second. For larger ILO models this is no longer true: the solver can run for minutes or hours, and it is useful to know what it is doing in the meantime. This notebook first compares the solvers that are available, then shows how to read the log that a solver writes while it works, and what you can do when a solve takes too long. The applications of this lecture are Multi-Period Inventory Planning and Robust Regression.

Learning outcomes

On completion of this notebook, you will be able to:

• compare open-source and commercial solvers, and explain where their speed comes from;

• read a solver log: find the incumbent, the best bound and the gap, and tell whether finding a good solution or proving its optimality is the hard part;

• help a slow solve with a time limit or gap tolerance, a warm start, or a tighter formulation.

Solvers¶

The solver is the software that does the optimizing; pulp builds the model and passes it to a solver. Roughly, there are two groups:

Quality¶

Solvers are compared on public sets of benchmark problems, such as the benchmarks that Hans Mittelmann maintains at https://plato.asu.edu/bench.html: how many of the problems a solver solves within a time limit, and how long it takes on average. On hard ILO problems the commercial solvers solve more problems, and solve them faster, than the open-source ones. For the models in this course, CBC is more than fast enough. The commercial solvers cost money outside academia, which is why not every company uses the fastest one.

Speedup¶

For LO, the classic method is the simplex algorithm (corner to corner). Since the 1980s, interior-point methods move through the interior of the feasible region instead and solve LO in provably polynomial time; modern solvers offer both. For ILO, branch and bound (see Integer Optimization) wraps around an LO solver, with many improvements on top: preprocessing that simplifies the model before the search starts, extra constraints (cuts) that tighten the LO relaxations, and heuristics that look for good integer solutions. Together with faster computers, this has made ILO solving many orders of magnitude faster. Linderoth (2017) estimates that from 1988 to 2017 the ILO algorithms became about 150,000 times faster and the computers about 17,000 times, so that a typical ILO that would have taken 124 years to solve in 1988 is solved in one second in 2017. Koch et al. (2022) ran the solvers of 2001 and of 2020 on the same benchmark problems: over these two decades the computers became about 20 times faster and the ILO algorithms about 50 times, a total speedup of about 1000.

pulp can call any installed solver without changing the model; only the .solve(...) line changes. To see what is available here:

Without installing anything, you can also export the model to a standard .mps file and submit it to the free NEOS Server, which hosts many solvers including commercial ones. If the variable and constraint names might leak information about your data, pulp can anonymize them on export with rename=1.

Calling Another Solver¶

The list above has two solvers: CBC, which comes with pulp, and HiGHS, an open-source solver for LO and ILO (Huangfu & Hall, 2018). In Colab, both are available without installing anything. On your own computer, pip install "pulp[highs]==3.3.2" installs pulp together with HiGHS (the Python package highspy). To solve a model with HiGHS, only the solver in the .solve(...) call changes:

product_mix.solve(pulp.HiGHS(msg=False))

The options msg, timeLimit and gapRel work the same as for PULP_CBC_CMD. pulp also has pulp.HiGHS_CMD, which calls a separate HiGHS program instead of the Python package. That program is not installed by pip or in Colab, so use pulp.HiGHS.

Reading the Solve Log¶

An Example That Takes Some Effort¶

The small examples of the previous lectures are solved before there is anything to watch. To see a solver at work, we take the single-machine scheduling model from Machine Scheduling and give it 10 jobs instead of 3. The job data (durations, release dates, due dates and tardiness costs) is drawn at random with a fixed seed, so every run gets the same instance. The model code is the same as in Lecture 9, placed in a function so that we can build it again for other instances.

With msg=True, CBC writes its log to the place where Python itself runs. In a notebook (also in Colab) that is often not the notebook, so the log does not always appear below the cell. Passing logPath writes the log to a file instead, which we can then read and print. The full log is long; the lines that matter most are printed below.

Incumbent, Best Bound and Gap¶

Recall from branch and bound that the solver keeps two numbers during the search. For a minimization problem:

• The incumbent is the best integer solution found so far. Its objective value, the incumbent value, is an upper bound on the optimal value: the optimum is at least as good. In the CBC log, every line Integer solution of ... found reports a new incumbent, and best solution in a progress line is the incumbent value. (Solver logs often label the incumbent value just “incumbent”; Gurobi’s log, for example, has a column Incumbent.)

• The best bound is a lower bound on the optimal value: no integer solution can be better. It is the smallest LO relaxation value among the nodes of the tree that are still open, and it can only go up during the search. CBC calls it best possible.

The gap is the distance between the two, relative to the incumbent value:

A gap of 0 means that the incumbent is proven optimal, and that is when the solver stops. (Solvers differ slightly in what they divide by; CBC’s summary line at the end of a stopped run divides by the best bound.)

In the log above, CBC finds an incumbent right away, with a heuristic called the feasibility pump. This happens at the root node, the first node of the branch-and-bound tree: the LO relaxation of the whole model, before any variable is branched on. The line Cbc0010I After 0 nodes is a progress line at the root node, with the incumbent value (best solution) and the best bound (best possible). There, the best bound is still below 10, while the incumbent value is almost 80. CBC writes such a progress line after every 1000 nodes of the tree (see CbcModel.cpp in the CBC source code); this solve needs fewer than 1000 nodes, so there is only the one at the root. The later incumbents are found during the search, and the last one is the optimum. After that, CBC keeps working until the bound has come up to the incumbent value: the line Search completed is the moment optimality is proven. Figure 1 shows both numbers over time.

Finding Versus Proving¶

A log like this one tells you which part of the work is hard for the solver. Three patterns come back again and again.

The first pattern is the one above, and it is the most common one for ILO models: the optimum (or a solution close to it) is found early, and most of the time goes into proving that nothing better exists. The reason is a weak best bound. In this model the LO relaxation may set each $y_{ij}$ to $\tfrac12$, so every no-overlap constraint is relaxed by $M/2$ and the jobs may overlap. We can check this by solving only the LO relaxation, which mip=False does:

The LO relaxation allows zero tardiness, while the optimum has a total weighted tardiness of 75. Branch and bound has to close this whole gap by branching.

The second pattern appears when an instance is too large to solve to optimality in the time we have. Then the incumbent improves step by step over time, and the question is when it is good enough to stop. With 12 jobs, CBC does not finish within 30 seconds. We stop it after 5 seconds with a time limit, timeLimit=5:

Two things stand out. First, pulp.LpStatus says Optimal, even though CBC stopped on the time limit. The solution is feasible, but not proven optimal: pulp’s sol_status tells the difference (Solution Found rather than Optimal Solution Found), and the log says Stopped on time limit. Always check one of these after a solve with a limit. Second, the best bound has hardly moved from its root value, so the gap is still very large. The incumbent may well be close to optimal, but the solver cannot show it (Figure 2).

The third pattern is a solver that has a hard time finding any integer solution at all: the log shows many nodes but no Integer solution line, and the incumbent in the progress lines stays at a huge placeholder value (CBC prints 1e+50). Without an incumbent, the solver can prune nothing by bound, and a time limit leaves you with nothing. This happens with many equality constraints or tight hard deadlines, where few combinations of integer values are feasible. The remedy is to give the solver a feasible solution to start from, which is the warm start below.

Helping a Solver¶

Once the log shows where the difficulty is, there are three things you can do about it.

Accept a Good Solution: Time Limits and Gap Tolerances¶

Data such as forecast demand or processing times is rarely known up to 1%, so a solution that is proven to be within 1% of optimal is usually good enough. A gap tolerance tells the solver to stop as soon as the gap is below a given fraction: gapRel=0.01 in PULP_CBC_CMD stops at a gap of 1%. A time limit (timeLimit, in seconds) stops the solver after a fixed time, whatever the gap. Both are options in the same call:

schedule.solve(pulp.PULP_CBC_CMD(msg=False, gapRel=0.01, timeLimit=60))

A gap tolerance only helps if the best bound comes close to the incumbent before the end. In the scheduling examples it does not: the bound stays low until the very end, so gapRel=0.01 would not stop the search any earlier, and only a time limit saves time. The price is that you then do not know how good the solution is.

Warm Starting¶

In a warm start, you give the solver a feasible solution to begin with. It then has an incumbent from the start: this helps with the third pattern above, and it lets the solver prune nodes by bound from the first node on. Good starting solutions come from a simple rule of thumb (a heuristic, see Heuristics), or from the solution of an earlier, almost identical problem, such as yesterday’s schedule when one job has changed.

For the scheduling instance, a simple rule is earliest due date first: process the jobs in order of their due dates, each as soon as the machine is free and the job is released. In pulp, .setInitialValue(...) sets a variable’s starting value, and warmStart=True passes these values to CBC. The starting solution must give a value to every variable, including the binary order variables.

The line Cbc0045I MIPStart provided solution shows that CBC accepted the starting solution as its first incumbent. The search then finds the same optimum as before, in about the same time: a warm start helps to find good solutions, not to prove optimality, and here finding was never the hard part. Its value shows in problems of the third pattern, and in loops that re-solve a large ILO many times with small data changes, such as the simulation-optimization loop in Simulation Optimization.

A Tighter Formulation¶

The most effective help is often a better model. The weak bound of the scheduling model comes from its big-M constraints, and a large $M$ makes the LO relaxation weaker. That is why Machine Scheduling chooses $M$ as small as possible while keeping the model correct. For this instance even the smallest valid $M$ leaves an LO relaxation value of 0, and other formulations of the same problem, with more variables, give much stronger bounds. Choosing between formulations goes beyond this course, but the log tells you when it is worth trying: when the best bound hardly moves.

A Note on Modeling Tools¶

Before Python libraries like pulp, models were usually written in an algebraic modeling language (AML), such as AMPL, GAMS or AIMMS: a language for writing a model in near-mathematical notation, kept separate from the data, and handed to whichever solver you choose. pulp gives the same model and data separation (Separating Data from the Model) inside Python, together with everything Python offers for preparing data and analyzing solutions. AMLs are still in use, mostly for large models that have been maintained for years, but for new projects a Python library is now the usual choice.

References¶

• Koole, G. (2019). An Introduction to Business Analytics. §6.6 “Modeling Tools”. AMPL and spreadsheet material replaced with pulp.

• CBC source code, Cbc/src/CbcModel.cpp (version 2.10.3, the version included in PuLP 3.3.2). 2.10.3/Cbc/src/CbcModel.cpp

• Huangfu, Q., & Hall, J. A. J. (2018). Parallelizing the dual revised simplex method. Mathematical Programming Computation, 10(1), 119-142. HiGHS: https://highs.dev

• Koch, T., Berthold, T., Pedersen, J., & Vanaret, C. (2022). Progress in mathematical programming solvers from 2001 to 2020. EURO Journal on Computational Optimization, 10, 100031. Koch et al. (2022)

• Linderoth, J. (2017). Talk at FOCAPO 2017 (Foundations of Computer-Aided Process Operations). Numbers as shown in the Lecture 10 slides.

• Mittelmann, H. Benchmarks for optimization software. https://plato.asu.edu/bench.html

• PuLP issue #448, Warm Start feature not working for CBC with PuLP on Windows. coin-or/pulp#448

• PuLP documentation: https://coin-or.github.io/pulp/