feat: add support for initial solutions with HiGHS (#83)

* build: switch over to highs bindings from git

* feat: support initial solutions for HiGHS

* test: add first tests for HiGHS solver

* build: update highs to 1.7.0

* Revert "build: update highs to 1.7.0"

This reverts commit 8e0cb62.

* build: update highs and russcip

Cargo.toml

@@ -24,8 +24,8 @@ minilp = ["microlp"] # minilp is not maintained anymore, we use the microlp fork
 coin_cbc = { version = "0.1", optional = true, default-features = false }
 microlp = { version = "0.2.6", optional = true }
 lpsolve = { version = "0.1", optional = true }
-highs = { version = "1.5.0", optional = true }
-russcip = { version = "0.4.1", optional = true }
+highs = { version = "1.7.0", optional = true }
+russcip = { version = "0.5.1", optional = true }
 lp-solvers = { version = "1.0.0", features = ["cplex"], optional = true }
 cplex-rs = { version = "0.1", optional = true }
 clarabel = { version = "0.9.0", optional = true, features = [] }

src/solvers/highs.rs

@@ -9,9 +9,10 @@ use crate::{
     solvers::DualValues,
     variable::{UnsolvedProblem, VariableDefinition},
 };
-use crate::{Constraint, IntoAffineExpression, Variable};
+use crate::{Constraint, IntoAffineExpression, Variable, WithInitialSolution};
 use highs::HighsModelStatus;
 use std::collections::HashMap;
+use std::iter::FromIterator;
 
 /// The [highs](https://docs.rs/highs) solver,
 /// to be used with [UnsolvedProblem::using].
@@ -48,6 +49,7 @@ pub fn highs(to_solve: UnsolvedProblem) -> HighsProblem {
         sense,
         highs_problem,
         columns,
+        initial_solution: None,
         verbose: false,
         options: Default::default(),
     }
@@ -170,6 +172,7 @@ pub struct HighsProblem {
     sense: highs::Sense,
     highs_problem: highs::RowProblem,
     columns: Vec<highs::Col>,
+    initial_solution: Option<Vec<(Variable, f64)>>,
     verbose: bool,
     options: HashMap<String, HighsOptionValue>,
 }
@@ -250,6 +253,15 @@ impl SolverModel for HighsProblem {
     fn solve(mut self) -> Result<Self::Solution, Self::Error> {
         let verbose = self.verbose;
         let options = std::mem::take(&mut self.options);
+        let initial_solution = self.initial_solution.as_ref().map(|pairs| {
+            pairs
+                .iter()
+                .fold(vec![0.0; self.columns.len()], |mut sol, (var, val)| {
+                    sol[var.index()] = *val;
+                    sol
+                })
+        });
+
         let mut model = self.into_inner();
         if verbose {
             model.set_option(&b"output_flag"[..], true);
@@ -265,6 +277,10 @@ impl SolverModel for HighsProblem {
             }
         }
 
+        if initial_solution.is_some() {
+            model.set_solution(initial_solution.as_deref(), None, None, None);
+        }
+
         let solved = model.solve();
         match solved.status() {
             HighsModelStatus::NotSet => Err(ResolutionError::Other("NotSet")),
@@ -305,6 +321,16 @@ impl SolverModel for HighsProblem {
     }
 }
 
+impl WithInitialSolution for HighsProblem {
+    fn with_initial_solution(
+        mut self,
+        solution: impl IntoIterator<Item = (Variable, f64)>,
+    ) -> Self {
+        self.initial_solution = Some(Vec::from_iter(solution));
+        self
+    }
+}
+
 /// The solution to a highs problem
 #[derive(Debug)]
 pub struct HighsSolution {
@@ -350,3 +376,67 @@ impl WithMipGap for HighsProblem {
         self.set_mip_rel_gap(mip_gap)
     }
 }
+
+#[cfg(test)]
+mod tests {
+    use crate::{constraint, variable, variables, Solution, SolverModel, WithInitialSolution};
+
+    use super::highs;
+    #[test]
+    fn can_solve_with_inequality() {
+        let mut vars = variables!();
+        let x = vars.add(variable().clamp(0, 2));
+        let y = vars.add(variable().clamp(1, 3));
+        let solution = vars
+            .maximise(x + y)
+            .using(highs)
+            .with((2 * x + y) << 4)
+            .solve()
+            .unwrap();
+        assert_eq!((solution.value(x), solution.value(y)), (0.5, 3.))
+    }
+
+    #[test]
+    fn can_solve_with_initial_solution() {
+        // Solve problem initially
+        let mut vars = variables!();
+        let x = vars.add(variable().clamp(0, 2));
+        let y = vars.add(variable().clamp(1, 3));
+        let solution = vars
+            .maximise(x + y)
+            .using(highs)
+            .with((2 * x + y) << 4)
+            .solve()
+            .unwrap();
+        // Recreate same problem with initial values slightly off
+        let initial_x = solution.value(x) - 0.1;
+        let initial_y = solution.value(x) - 1.0;
+        let mut vars = variables!();
+        let x = vars.add(variable().clamp(0, 2));
+        let y = vars.add(variable().clamp(1, 3));
+        let solution = vars
+            .maximise(x + y)
+            .using(highs)
+            .with((2 * x + y) << 4)
+            .with_initial_solution([(x, initial_x), (y, initial_y)])
+            .solve()
+            .unwrap();
+
+        assert_eq!((solution.value(x), solution.value(y)), (0.5, 3.))
+    }
+
+    #[test]
+    fn can_solve_with_equality() {
+        let mut vars = variables!();
+        let x = vars.add(variable().clamp(0, 2).integer());
+        let y = vars.add(variable().clamp(1, 3).integer());
+        let solution = vars
+            .maximise(x + y)
+            .using(highs)
+            .with(constraint!(2 * x + y == 4))
+            .with(constraint!(x + 2 * y <= 5))
+            .solve()
+            .unwrap();
+        assert_eq!((solution.value(x), solution.value(y)), (1., 2.));
+    }
+}