good_lp/solvers/
clarabel.rs

1//! A solver that uses [clarabel](https://oxfordcontrol.github.io/ClarabelDocs/stable/), a pure rust solver.
2
3use crate::affine_expression_trait::IntoAffineExpression;
4use crate::expression::LinearExpression;
5use crate::variable::UnsolvedProblem;
6use crate::{
7    constraint::ConstraintReference,
8    solvers::{ObjectiveDirection, ResolutionError, Solution, SolverModel},
9    SolutionStatus,
10};
11use crate::{Constraint, DualValues, SolutionWithDual, Variable};
12
13use clarabel::algebra::CscMatrix;
14use clarabel::solver::implementations::default::DefaultSettingsBuilder;
15use clarabel::solver::SupportedConeT::{self, *};
16use clarabel::solver::{DefaultSolution, SolverStatus};
17use clarabel::solver::{DefaultSolver, IPSolver};
18
19/// The [clarabel](https://oxfordcontrol.github.io/ClarabelDocs/stable/) solver,
20/// to be used with [UnsolvedProblem::using].
21pub fn clarabel(to_solve: UnsolvedProblem) -> ClarabelProblem {
22    let UnsolvedProblem {
23        objective,
24        direction,
25        variables,
26    } = to_solve;
27    let coef = if direction == ObjectiveDirection::Maximisation {
28        -1.
29    } else {
30        1.
31    };
32    let mut objective_vector = vec![0.; variables.len()];
33    for (var, obj) in objective.linear_coefficients() {
34        objective_vector[var.index()] = obj * coef;
35    }
36    let constraints_matrix_builder = CscMatrixBuilder::new(variables.len());
37    let mut settings = DefaultSettingsBuilder::default();
38    settings.verbose(false).tol_feas(1e-9);
39    let mut p = ClarabelProblem {
40        objective: objective_vector,
41        constraints_matrix_builder,
42        constraint_values: Vec::new(),
43        variables: variables.len(),
44        settings,
45        cones: Vec::new(),
46    };
47    // add trivial constraints embedded in the variable definitions
48    for (var, def) in variables.iter_variables_with_def() {
49        if def.is_integer {
50            panic!("Clarabel doesn't support integer variables")
51        }
52        if def.min != f64::NEG_INFINITY {
53            p.add_constraint(var >> def.min);
54        }
55        if def.max != f64::INFINITY {
56            p.add_constraint(var << def.max);
57        }
58    }
59    p
60}
61
62/// A clarabel model
63pub struct ClarabelProblem {
64    constraints_matrix_builder: CscMatrixBuilder,
65    constraint_values: Vec<f64>,
66    objective: Vec<f64>,
67    variables: usize,
68    settings: DefaultSettingsBuilder<f64>,
69    cones: Vec<SupportedConeT<f64>>,
70}
71
72impl ClarabelProblem {
73    /// Access the problem settings
74    pub fn settings(&mut self) -> &mut DefaultSettingsBuilder<f64> {
75        &mut self.settings
76    }
77
78    /// Convert the problem into a clarabel solver
79    /// panics if the problem is not valid
80    pub fn into_solver(self) -> DefaultSolver<f64> {
81        let settings = self.settings.build().expect("Invalid clarabel settings");
82        let quadratic_objective = &CscMatrix::zeros((self.variables, self.variables));
83        let objective = &self.objective;
84        let constraints = &self.constraints_matrix_builder.build();
85        let constraint_values = &self.constraint_values;
86        let cones = &self.cones;
87        DefaultSolver::new(
88            quadratic_objective,
89            objective,
90            constraints,
91            constraint_values,
92            cones,
93            settings,
94        ).expect("Invalid clarabel problem. This is likely a bug in good_lp. Problems should always have coherent dimensions.")
95    }
96}
97
98impl SolverModel for ClarabelProblem {
99    type Solution = ClarabelSolution;
100    type Error = ResolutionError;
101
102    fn solve(self) -> Result<Self::Solution, Self::Error> {
103        let mut solver = self.into_solver();
104        solver.solve();
105        match solver.solution.status {
106            SolverStatus::PrimalInfeasible | SolverStatus::AlmostPrimalInfeasible => {
107                Err(ResolutionError::Infeasible)
108            }
109            SolverStatus::Solved
110            | SolverStatus::AlmostSolved
111            | SolverStatus::AlmostDualInfeasible
112            | SolverStatus::DualInfeasible => Ok(ClarabelSolution {
113                solution: solver.solution,
114            }),
115            SolverStatus::Unsolved => Err(ResolutionError::Other("Unsolved")),
116            SolverStatus::MaxIterations => Err(ResolutionError::Other("Max iterations reached")),
117            SolverStatus::MaxTime => Err(ResolutionError::Other("Time limit reached")),
118            SolverStatus::NumericalError => Err(ResolutionError::Other("Numerical error")),
119            SolverStatus::InsufficientProgress => Err(ResolutionError::Other("No progress")),
120            SolverStatus::CallbackTerminated => Err(ResolutionError::Other("Callback terminated")),
121        }
122    }
123
124    fn add_constraint(&mut self, constraint: Constraint) -> ConstraintReference {
125        self.constraints_matrix_builder
126            .add_row(constraint.expression.linear);
127        let index = self.constraint_values.len();
128        self.constraint_values.push(-constraint.expression.constant);
129        // Cones indicate the type of constraint. We only support nonnegative and equality constraints.
130        // To avoid creating a new cone for each constraint, we merge them.
131        let next_cone = if constraint.is_equality {
132            ZeroConeT(1)
133        } else {
134            NonnegativeConeT(1)
135        };
136        let prev_cone = self.cones.last_mut();
137        match (prev_cone, next_cone) {
138            (Some(ZeroConeT(a)), ZeroConeT(b)) => *a += b,
139            (Some(NonnegativeConeT(a)), NonnegativeConeT(b)) => *a += b,
140            (_, next_cone) => self.cones.push(next_cone),
141        };
142        ConstraintReference { index }
143    }
144
145    fn name() -> &'static str {
146        "Clarabel"
147    }
148}
149
150/// The solution to a clarabel problem
151pub struct ClarabelSolution {
152    solution: DefaultSolution<f64>,
153}
154
155impl ClarabelSolution {
156    /// Returns the clarabel solution object. You can use it to dynamically add new constraints
157    pub fn into_inner(self) -> DefaultSolution<f64> {
158        self.solution
159    }
160
161    /// Borrow the clarabel solution object
162    pub fn inner(&self) -> &DefaultSolution<f64> {
163        &self.solution
164    }
165}
166
167impl Solution for ClarabelSolution {
168    fn status(&self) -> SolutionStatus {
169        SolutionStatus::Optimal
170    }
171    fn value(&self, variable: Variable) -> f64 {
172        self.solution.x[variable.index()]
173    }
174}
175
176impl<'a> SolutionWithDual<'a> for ClarabelSolution {
177    type Dual = &'a ClarabelSolution;
178
179    fn compute_dual(&'a mut self) -> Self::Dual {
180        self
181    }
182}
183
184impl DualValues for &ClarabelSolution {
185    fn dual(&self, constraint: ConstraintReference) -> f64 {
186        self.solution.z[constraint.index]
187    }
188}
189
190struct CscMatrixBuilder {
191    /// Indicates the row index of the corresponding element in `nzval`
192    rowval: Vec<Vec<usize>>,
193    /// All non-zero values in the matrix, in column-major order
194    nzval: Vec<Vec<f64>>,
195    n_rows: usize,
196    n_cols: usize,
197}
198
199impl CscMatrixBuilder {
200    fn new(n_cols: usize) -> Self {
201        Self {
202            rowval: vec![Vec::new(); n_cols],
203            nzval: vec![Vec::new(); n_cols],
204            n_rows: 0,
205            n_cols,
206        }
207    }
208    fn add_row(&mut self, row: LinearExpression) {
209        for (var, value) in row.linear_coefficients() {
210            self.rowval[var.index()].push(self.n_rows);
211            self.nzval[var.index()].push(value);
212        }
213        self.n_rows += 1;
214    }
215    fn build(self) -> clarabel::algebra::CscMatrix {
216        let mut colptr = Vec::with_capacity(self.n_cols + 1);
217        colptr.push(0);
218        for col in &self.rowval {
219            colptr.push(colptr.last().unwrap() + col.len());
220        }
221        clarabel::algebra::CscMatrix::new(
222            self.n_rows,
223            self.n_cols,
224            colptr,
225            fast_flatten_vecs(self.rowval),
226            fast_flatten_vecs(self.nzval),
227        )
228    }
229}
230
231fn fast_flatten_vecs<T: Copy>(vecs: Vec<Vec<T>>) -> Vec<T> {
232    // This is faster than vecs.into_iter().flatten().collect()
233    // because it doesn't need to allocate a new Vec
234    // (we take ownership of the first Vec and add the rest to it)
235    let size: usize = vecs.iter().map(|v| v.len()).sum();
236    let mut iter = vecs.into_iter();
237    let mut result = if let Some(v) = iter.next() {
238        v
239    } else {
240        return Vec::new();
241    };
242    result.reserve_exact(size - result.len());
243    for v in iter {
244        result.extend_from_slice(&v);
245    }
246    result
247}
248
249#[cfg(test)]
250mod tests {
251
252    use super::*;
253    use crate::variables;
254
255    #[test]
256    fn test_csc_matrix_builder() {
257        variables! {vars:
258            x;
259            y;
260            z;
261        }
262        let mut builder = CscMatrixBuilder::new(3);
263        builder.add_row((y + 2 * z).linear);
264        builder.add_row((3 * x + 4 * y + 5 * z).linear);
265        let matrix = builder.build();
266        /* The matrix is:
267        [ 0 1 2 ]
268        [ 3 4 5 ]
269        */
270        assert_eq!(matrix.m, 2); // 2 rows
271        assert_eq!(matrix.n, 3); // 3 columns
272        assert_eq!(matrix.get_entry((0, 0)), None); // get_entry((row, col))
273        assert_eq!(matrix.get_entry((0, 1)), Some(1.));
274        assert_eq!(matrix.get_entry((0, 2)), Some(2.));
275        assert_eq!(matrix.get_entry((1, 0)), Some(3.));
276        assert_eq!(matrix.get_entry((1, 1)), Some(4.));
277        assert_eq!(matrix.get_entry((1, 2)), Some(5.));
278    }
279}
good_lp/solvers/clarabel.rs

good_lp/solvers/
clarabel.rs
