Support iterative solve of the same solver_model (when using highs)

[dannyopts]

I'm not sure what the right approach is here, so I'm going to write what I have done and found.

I want to be able to solve the same model repeatedly with an altered cost function to understand model sensitivity.

My specific use case is that I have a PyPSA model which takes about 50 seconds to solve using highs.

I want to solve this model many time with different costs to understand If I update the cost function.

One way is in a loop, just keep resolving the problem, but I want to take advantage of the warm start power of highs.

I have done this by accessing the solver_model and calling

m.solver_model.changeColsCost(len(coeffsi), np.arange(len(coeffsi), dtype=np.int32), coeffsi)

See below for 15 times taken for this kind of sequential solve, this gives approximately a 5x speed up vs the independent sequential solve (the first is kind of representative of the time for solving the problem without warm start)

[58.32, 23.88, 2.74 8.54, 14.34, 17.98, 2.55, 22.21, 20.98, 24.55, 17.32, 13.22, 5.00, 4.12, 3.07]

But then I cant apply this solution back to the linopy model without extracting a bunch of code from model.solve related to applying the solution from a solver_model to the model.

def apply_solution(model, h):
    CONDITION_MAP = {}
    condition = h.modelStatusToString(h.getModelStatus()).lower()
    termination_condition = CONDITION_MAP.get(condition, condition)
    status = Status.from_termination_condition(termination_condition)
    status.legacy_status = condition

    def get_solver_solution() -> Solution:
        objective = h.getObjectiveValue()
        solution = h.getSolution()

        sol = pd.Series(solution.col_value, h.getLp().col_names_, dtype=float).pipe(
            set_int_index
        )
        dual = pd.Series(solution.row_dual, h.getLp().row_names_, dtype=float).pipe(
            set_int_index
        )

        return Solution(sol, dual, objective)

    solution = safe_get_solution(status, get_solver_solution)

    result = Result(status, solution, h)
    
    model.objective._value = result.solution.objective
    model.status = result.status.status.value
    model.termination_condition = result.status.termination_condition.value
    model.solver_model = result.solver_model
    
    sol = result.solution.primal.copy()
    sol.loc[-1] = np.nan

    for name, var in model.variables.items():
        idx = np.ravel(var.labels)
        try:
            vals = sol[idx].values.reshape(var.labels.shape)
        except KeyError:
            vals = sol.reindex(idx).values.reshape(var.labels.shape)
        var.solution = xr.DataArray(vals, var.coords)

By pulling out the application of a solution from the solve method as a seperate public function, it would make this usage much easier.

I know that the run_{solver} methods take in a warmstart_fn, but this isnt supported for highs, maybe simply supporting this would be sufficent? But having experimented with this I think that the last release of highspy doesnt yet have the bindings for setSolution, but master does (presumably this is why it isnt supported?)

This would incur an additional overhead of writing and reading the basis file on each loop, but maybe that is fine for most usecases (I havent been able to test how this impacts performance, since I dont have v1.6 of highs).

Am I doing something stupid, or is this kind of iterative solve not supported with highs via linopy right now? Would it be something that you want to support?

If yes, I would be happy to work on a PR if you could help me scope out what the right interfaces are.

For a bit more context here is an example of how I am using this code with pypsa:

def solve_range(start_price, step_size, iterations):
    price = start_price
    n = create_network(price)

    # solve the initial model
    n.optimize(solver_name="highs", parallel="on")

    m = n.model
    del n.model
    yield n, time.time() - start

    for i in range(iterations):
        price += step_size
        # create the similar but different network (with n.model already created)
        n_i = create_model(price)
       
        mi = n_i._network.model
        dfi = mi.objective.flat
        coeffsi = dfi.coeffs.values
        # update the highs model to include the new cost col - the constraints dont vary between runs
        m.solver_model.changeColsCost(len(coeffsi), np.arange(len(coeffsi), dtype=np.int32), coeffsi)
        
        # solve the highs model
        m.solver_model.run()
        
        # write the solution back to the model
        apply_solution(m, m.solver_model)
        n_i.model = m
        # write the solution back to the network (functions imported from pypsa)
        assign_solution(n_i)
        assign_duals(n_i)
        post_processing(n_i)
        del n_i.model

        yield h2_grid_model_i, end - start

[FabianHofmann]

Hey @dannyopts, you are totally right. I find this is indeed a missing feature in linopy. It would not require too much work, but I am a little bit constrained atm. The solve function would be needed to be subdivided, ideally converted into classes "Solver" which has underlying functions like, "get_solution", "assign_parameters", etc. If you want to propose something, you are more than welcome, otherwise this will lay on my todo list for some time still.

[dannyopts]

I will take a look and see how I get on

[dannyopts]

@FabianHofmann I am thinking something like:

We replace the f"run_{solver_name}" function with an abstract Solver class which is implemented for each solver.

Roughly I see this looking something like this:

class Solver:
    """A persistent instance of a solver model which can be updated and resolved
    
    Concrete versions of this class replace f"run_{solver_name}"
    """
    def __init__(self, model, sanitize_zeros=True, problem_fn=None, log_fn=None, **solver_options):
        self.sanitize_zeros = sanitize_zeros
        self.problem_fn = problem_fn
        self.log_fn = log_fn
        self.solver_model = self.create_solver_model(model)
        self.solver_options = solver_options
        
    @staticmethod
    def get_solver(solver_name, **kwargs):
        if solver_name == "highs":
            return HighsSolver(**kwargs)
        if solver_name == "gurobi":
            return GurobiSolver(**kwargs)
        ...
    
    def run(self) -> SolverStatus:
        """Solver the solver_model"""
        raise NotImplementedError()
    
    def compute_IIS(self):
        raise NotImplementedError()
    
    def update_cost_coeffs(self, c) -> None:
        """Set the cost vector on the solver_model"""
        raise NotImplementedError()
    
    def update_constraint(self, constraint_name: str, values: xr.DataArray) -> None:
        """Set the cost vector on the solver_model"""
        pass
    
    def update_rhs(self, b) -> None:
        """Set the b vector"""
        pass
    
    def get_result(self) -> Result: 
        pass

    def close(self) -> None:
        """Free any resources (eg release the license if using gurobi)""""
        pass

(possibly making it a context manager would be a nice touch?)

In the result class we would then add the responsibility of writing the result back to the model:

@dataclass
class Result:
    ...
       
    def update_model(self, model):
        model.objective._value = self.solution.objective
        model.status = self.status.status.value
        model.termination_condition = self.status.termination_condition.value
        model.solver_model = self.solver_model

        if not self.status.is_ok:
            return
        self.assign_solution(self, model)
        self.assign_dual(self, model)

    def assign_solution(self, model):
        # map solution and dual to original shape which includes missing values
        sol = self.solution.primal.copy()
        sol.loc[-1] = nan

        for name, var in model.variables.items():
            idx = np.ravel(var.labels)
            try:
                vals = sol[idx].values.reshape(var.labels.shape)
            except KeyError:
                vals = sol.reindex(idx).values.reshape(var.labels.shape)
            var.solution = xr.DataArray(vals, var.coords)

    def assign_dual(self, model):
        # map solution and dual to original shape which includes missing values
        if not self.solution.dual.empty:
            dual = self.solution.dual.copy()
            dual.loc[-1] = nan

            for name, con in model.constraints.items():
                idx = np.ravel(con.labels)
                try:
                    vals = dual[idx].values.reshape(con.labels.shape)
                except KeyError:
                    vals = dual.reindex(idx).values.reshape(con.labels.shape)
                con.dual = xr.DataArray(vals, con.labels.coords)

And in the Model class solve would become something like:

def solve(
        self,
        ...
    ):
        if remote:
            ...
            return self.status, self.termination_condition

        solver = Solver.get_solver(
            solver_name=solver_name, 
            io_api=io_api, 
            env=env, 
            sanitize_zeros=sanitize_zeros, 
            problem_fn=problem_fn,
            solution_fn=solution_fn,
            log_fn=log_fn,
            basis_fn=basis_fn,
            warmstart_fn=warmstart_fn,
            **solver_options
        )
        with solver:
          result = solver.run()
          result.info()
          
          result.update_model(self)
          
        return result

But a user could also create a Solver instance outside of solve and have a persistent handle to the same underlying model instance and update at their leisure.

Note: this would also resolve the issue which was mentioned in #192 about using a remote gurobi server and wanting to compute an IIS.

Does this sound like a reasonable path forward @FabianHofmann? Any feedback?

[FabianHofmann]

Oh, I like this idea! Some thoughts (hopefully later more):

• I thought about having a Solver class as a skeleton which is then inherited by the individual solvers, like

class Solver:
    def __init__(self, **kwargs):
        ....

    def solve(self):
        raise NotImplementedError

class SolverA(Solver):
    ....
    def solve(self):
        # Implement the solve method for SolverA
        pass

class SolverB(Solver):
    ....
    def solve(self):
        # Implement the solve method for SolverB
        pass

• I would be in favor of moving the assign_solution to the Model class, following a composition approach (ie the model pull from the solution)

[Cellophil]

Are there news on this front in 2024? I very much love the idea! I dabbled a bit with iterative solving with highs ( adding new constraints), and my two cents are that highspy keeps a found solution so rerunning an "h" instance acts like a warmstart. And, that highs indexes are not necessarily identical to linopys, in case variables are deleted.
I would definitely make use of a cool solution as proposed, looking forward.