Use initial solution for MIP (and LP?) with APPSI HiGHS interface

[FGlazov]

Summary

I wish to use initial solutions to the MIP and LP HiGHS solver using the APPSI HighsPy interface.

Rationale

It is not uncommon to be able to create decent initial solutions to a MIP problem, for example, by using a heuristic. The MIP solver can then use the initial solution to speed up the solving process - potentially giving you an optimal or better solution, if you're time constrained in solving the MIP model.

Being able to pass solutions to a MIP solver also allows for better naive implementations of various LP/MIP solving techniques, such as a Bender's Decomposition or a Rolling Horizon Heuristic.

This can also be helpful for LP problems, I'm not sure if a similar thing can be done for the QP solver.

Description

The current solve functions found in APPSI Highspy interface look like this:

    def _solve(self, timer: HierarchicalTimer):
        config = self.config
        options = self.highs_options

        ostreams = [
            LogStream(
                level=self.config.log_level, logger=self.config.solver_output_logger
            )
        ]
        if self.config.stream_solver:
            ostreams.append(sys.stdout)

        with TeeStream(*ostreams) as t:
            with capture_output(output=t.STDOUT, capture_fd=True):
                self._solver_model.setOptionValue('log_to_console', True)
                if config.logfile != '':
                    self._solver_model.setOptionValue('log_file', config.logfile)

                if config.time_limit is not None:
                    self._solver_model.setOptionValue('time_limit', config.time_limit)
                if config.mip_gap is not None:
                    self._solver_model.setOptionValue('mip_rel_gap', config.mip_gap)

                for key, option in options.items():
                    self._solver_model.setOptionValue(key, option)
                timer.start('optimize')
                self._solver_model.run()
                timer.stop('optimize')

        return self._postsolve(timer)

    def solve(self, model, timer: HierarchicalTimer = None) -> Results:
        StaleFlagManager.mark_all_as_stale()
        if self._last_results_object is not None:
            self._last_results_object.solution_loader.invalidate()
        if timer is None:
            timer = HierarchicalTimer()
        if model is not self._model:
            timer.start('set_instance')
            self.set_instance(model)
            timer.stop('set_instance')
        else:
            timer.start('update')
            self.update(timer=timer)
            timer.stop('update')
        res = self._solve(timer)
        self._last_results_object = res
        if self.config.report_timing:
            logger.info('\n' + str(timer))
        return res

I expect there to be a call to self._solver_model.setSolution filled with the values I've passed to my Pyomo model in order to pass an initial solution for an LP/MIP problem. However, I was not able to find this call anywhere in the methods I've listed above, nor in any of the methods called inside solve and _solve.

Why I think solver_model.setSolution is the correct function to use here:

1. setSolution exists in self._solver_model, at least for recent versions: https://github.com/ERGO-Code/HiGHS/blob/master/src/highspy/_core/__init__.pyi#L1002C9-L1002C20
2. The documentation here suggests that this function corresponds, in name at least, to the C++ function that sets the initial solution (both for MIPs and LPs): https://ergo-code.github.io/HiGHS/dev/guide/further/#MIP
3. I've found this comment from Julian Hall (one of the maintainers of HiGHS) which suggest that this C++ binding was created in Python a year ago: example for LP warm start? ERGO-Code/HiGHS#617 (comment)

I am also open to other suggestions on how to pass initial solutions to HiGHs. For example, the CLI also supports passing a solution file: https://ergo-code.github.io/HiGHS/dev/executable/ - So an alternative could be to build an old-school CLI-based Pyomo solver, that passes the initial solution using the --read_solution_file option. But I think this is a rather crappy alternative, because this sort of implementation does not provide the benefits APPSI does (i.e. a persistent solver).

Additional information

I'm interested in this topic because I want to implement a Rolling Horizon Heursitic using HiGHS and Pyomo. My code looks like this:

def rolling_horizon_solve(model: pyo.Model, solver: Highs, config: OptConfig):
    original_time_limit = solver.config.time_limit
    solver.config.time_limit = config.MILP_ROLLING_HORIZON_TIME_LIMIT

    all_dates = sorted(model.dates)

    logger.info(f"Starting rolling horizon solve for {len(all_dates)} dates.")

    for i, rolling_horizon_date in enumerate(all_dates):
        logger.info(f"Rolling horizon solve for date {rolling_horizon_date}.")

        __fix_variables_for_date(model, all_dates, rolling_horizon_date)
        solver.solve(model)
        __unfix_variables(model)

        # TODO: Error handling if no solution found!
        # (Escape from loop to master problem solve)

        solver.load_vars()

    solver.config.time_limit = original_time_limit

    logger.info("Completed rolling horizon solve. Starting master problem solve.")

    # Master problem solve
    return solver.solve(model)

Where __fix_variables_for_date fixes past variables to their current values, and future ones to zero. Essentially, many calls to model.example_var_set.fix(0.0) (for future variables) and my_var.fix(my_var.value) (for past variables).

In any case, this heuristic is able to find a decent initial solution rather quickly (<10 minutes). What I want to do is to feed this initial solution back into my master problem, which runs for about 6 hours in total. The solution I find using this heuristic is about as good as the solution by cold starting my master problem and letting it compute for 2 hours. So if I'm able to pass the initial solution from the heuristic into my master problem, I'd get a little under 2 hours of more headroom to find a better solution to my master problem.

The thing is: The master problem solve does not benefit from the rolling horizon heuristic I did above it in my python snippet. Indeed, the solution logs for my master problem solve, with the rolling horizon heuristc activated, look like this for a simplified dev-model I ran for 200 seconds:

Coefficient ranges:
  Matrix [8e-03, 1e+02]
  Cost   [7e-06, 1e+02]
  Bound  [9e-02, 1e+04]
  RHS    [3e-02, 1e+03]
Presolving model
33964 rows, 42545 cols, 93058 nonzeros  0s
23554 rows, 32120 cols, 84902 nonzeros  0s
17549 rows, 29674 cols, 75776 nonzeros  0s
17166 rows, 27582 cols, 69505 nonzeros  0s

Solving MIP model with:
   17166 rows
   27582 cols (216 binary, 9708 integer, 1264 implied int., 16394 continuous)
   69505 nonzeros

Src: B => Branching; C => Central rounding; F => Feasibility pump; H => Heuristic; L => Sub-MIP;
     P => Empty MIP; R => Randomized rounding; S => Solve LP; T => Evaluate node; U => Unbounded;
     z => Trivial zero; l => Trivial lower; u => Trivial upper; p => Trivial point

        Nodes      |    B&B Tree     |            Objective Bounds              |  Dynamic Constraints |       Work      
Src  Proc. InQueue |  Leaves   Expl. | BestBound       BestSol              Gap |   Cuts   InLp Confl. | LpIters     Time

         0       0         0   0.00%   1685.05445      inf                  inf        0      0      0         0     0.2s
         0       0         0   0.00%   13022.035354    inf                  inf        0      0      3      9818     0.7s
         0       0         0   0.00%   18757.634001    inf                  inf    12843   2343    353     18405     8.7s
         0       0         0   0.00%   21149.900474    inf                  inf    10732   2377   2101    103396    43.3s

Symmetry detection completed in 0.0s
No symmetry present

        37       0         1   0.00%   21149.900474    inf                  inf    10749   2377   2121    164113    70.9s
       251     179         2   0.00%   21149.900474    inf                  inf    10760   2377   2360    231301   108.0s
       454     344         5   0.01%   21451.831279    inf                  inf    11741   2750   2450    245490   117.6s
       578     449         6   0.01%   21451.831279    inf                  inf    12338   3304   2517    257624   124.6s
       692     556         7   0.01%   21451.831279    inf                  inf    12548   3385   2580    274398   133.6s
       745     555         8   0.01%   21451.831279    inf                  inf    12911   3512   2593    291077   143.9s
       894     728        11   0.01%   21451.831279    inf                  inf    12913   3512   2676    302730   151.1s
       961     727        12   0.01%   21451.831279    inf                  inf    13501   3605   2739    435927   197.6s
       983     810        13   0.01%   21451.831279    inf                  inf    13502   3605   2746    437340   200.0s

Solving report
  Status            Time limit reached
  Primal bound      inf
  Dual bound        21451.8312792
  Gap               inf
  P-D integral      0
  Solution status   -
  Timing            200.03 (total)
                    0.00 (presolve)
                    0.00 (solve)
                    0.00 (postsolve)
  Max sub-MIP depth 2
  Nodes             983
  Repair LPs        0 (0 feasible; 0 iterations)
  LP iterations     437340 (total)
                    148725 (strong br.)
                    21937 (separation)
                    191302 (heuristics)

Whereas the solution logs for the last iteration of the rolling horizon heuristic look like this:

Coefficient ranges:
  Matrix [8e-03, 1e+02]
  Cost   [7e-06, 1e+02]
  Bound  [9e-02, 1e+04]
  RHS    [3e-02, 1e+03]
Presolving model
30994 rows, 35919 cols, 77100 nonzeros  0s
20416 rows, 26216 cols, 64779 nonzeros  0s
17305 rows, 20801 cols, 55271 nonzeros  0s
12128 rows, 19736 cols, 43220 nonzeros  0s
10736 rows, 17155 cols, 37234 nonzeros  0s

Solving MIP model with:
   10736 rows
   17155 cols (51 binary, 5588 integer, 1255 implied int., 10261 continuous)
   37234 nonzeros

Src: B => Branching; C => Central rounding; F => Feasibility pump; H => Heuristic; L => Sub-MIP;
     P => Empty MIP; R => Randomized rounding; S => Solve LP; T => Evaluate node; U => Unbounded;
     z => Trivial zero; l => Trivial lower; u => Trivial upper; p => Trivial point

        Nodes      |    B&B Tree     |            Objective Bounds              |  Dynamic Constraints |       Work      
Src  Proc. InQueue |  Leaves   Expl. | BestBound       BestSol              Gap |   Cuts   InLp Confl. | LpIters     Time

         0       0         0   0.00%   12201.033907    inf                  inf        0      0      0         0     0.2s
         0       0         0   0.00%   40658.399534    inf                  inf        0      0      3      4281     0.4s
 C       0       0         0   0.00%   47334.451323    84815.714809      44.19%     5755   1105    172      8601     1.8s
 L       0       0         0   0.00%   49247.175821    50209.5156         1.92%     8960   1255    172     10864     7.3s

0.9% inactive integer columns, restarting
Model after restart has 10639 rows, 17001 cols (48 bin., 5542 int., 1251 impl., 10160 cont.), and 36922 nonzeros

         0       0         0   0.00%   49299.615591    50209.5156         1.81%      975      0      0     21174     8.7s
         0       0         0   0.00%   49300.365256    50209.5156         1.81%      975    919      3     26000     8.9s
 L       0       0         0   0.00%   49702.558372    49876.730482       0.35%     2993   1091      3     27104    15.0s
         1       0         1 100.00%   49702.558372    49876.730482       0.35%     2993   1091      3     33356    15.1s

Solving report
  Status            Optimal
  Primal bound      49876.7304815
  Dual bound        49702.5583724
  Gap               0.349% (tolerance: 1%)
  P-D integral      2.3908187846
  Solution status   feasible
                    49876.7304815 (objective)
                    0 (bound viol.)
                    2.50022225146e-13 (int. viol.)
                    0 (row viol.)
  Timing            15.07 (total)
                    0.00 (presolve)
                    0.00 (solve)
                    0.00 (postsolve)
  Max sub-MIP depth 7
  Nodes             1
  Repair LPs        0 (0 feasible; 0 iterations)
  LP iterations     33356 (total)
                    0 (strong br.)
                    7776 (separation)
                    16473 (heuristics)

What I expect to happen is that the master problem solve starts off with the initial solution provided by the last iteration of the rolling horizon heuristic, and thus start with a primal bound of 49876.730482 - the dual bound should get reset and recomputed, since the underlying MIP has less constraints.

[FGlazov]

I think I've done enough research here that I might be able to do implement this feature on my own into the repository if nobody else wants to grab this issue.

What would be nice to know beforehand is if I misunderstood the documentation and if there's actually a way to pass initial solutions to the HiGHS-APPSI interface, so I don't stress myself out doing some unnecessary changes :) - I'd also be very thankful if someone else does the changes.

[mrmundt]

Hi @FGlazov - this isn't an answer to your request, but I wanted you to be aware. We are actually working on adding a new interface to HiGHS (and other solvers) using our future functionality, so I wouldn't recommend taking time to make changes to the APPSI code.

[kutlay]

I also think that this is a very needed feature, actually the only blocker for me right now to fully adapt Pyomo as a HiGHS user. I would be happy to contribute in implementation and testing of this feature as well.

[sk-surya]

I also need this feature badly. While the new interface is being developed, is there a workaround to do this with existing appsi highs interface? @kutlay @FGlazov @mrmundt

[FGlazov]

@sk-surya In my case, I used the CBC solver and fed it with an initial solution. The APPSI CBC Solver is a bit buggy in my experience (it doesn't interrupt correctly when the solver times out), I used the command line CBC solver variant that you can get via pyo.SolverFactory("cbc"), and make sure you call solve with warmstart=True.

Specifically in my use case I still use HiGHS to do the rolling horizon heuristic, as it's quick enough even without initial solution, and then feed the result from HiGHS into CBC. Quite hacky and requires two full model passes to two diffrent solver backends, but it's OK for my use case since it barely makes a dent in my 6-hour time budget and gives me better solutions.

[sk-surya]

@FGlazov, Interesting. In my case, I'm doing multiple MIP resolves (~500) to get sensitivities. Each resolve is very cheap to do (~ms). I guess I could build the model in pyomo, generate LP file, read it in Highspy and do the resolving there. But would need to figure out mapping the pyomo block components and modify the corresponding value in the loaded model. It would be great to just hack and get it working in APPSI HIGHS for the interim before the new interface. I can contribute, if @mrmundt or someone could provide some initial directions.

[mrmundt]

Hi, everyone! Wow, this got a lot of traffic. We have the new interfaces high on our priority list, but we only have so much time. (Unfortunately, all of us only work on Pyomo part time.) If someone wants to work on getting this integrated into APPSI, we are a-okay with that and welcome PRs. I recommend reviewing the Contribution Guide. I will also work on getting some of our design discussions published somewhere accessible for you all to look over (we have a large set of running notes but they're internal to the team).