feos_campd/solver/
nlp.rs

1use super::dual_vec_multiple::DualSVecMult;
2use super::{Gradient, MixedIntegerNonLinearProgram, OptimizationResult, OuterApproximation};
3#[cfg(feature = "ipopt")]
4use ipopt::{Ipopt, IpoptOption, SolveStatus};
5use ipopt_ad::{ADProblem, BasicADProblem, CachedADProblem};
6use nalgebra::{Const, DVector, SMatrix, SVector, U1};
7use num_dual::{Derivative, DualNum};
8#[cfg(feature = "ripopt")]
9use ripopt::{solve, SolverOptions};
10use std::fmt::Debug;
11
12struct Nlp<'a, P, const N_Y1: usize, const N_Y2: usize> {
13    y: SMatrix<f64, N_Y1, N_Y2>,
14    minlp: &'a P,
15}
16
17impl<'a, P, const N_Y1: usize, const N_Y2: usize> Nlp<'a, P, N_Y1, N_Y2> {
18    pub fn new(y: SMatrix<f64, N_Y1, N_Y2>, minlp: &'a P) -> Self {
19        Self { y, minlp }
20    }
21}
22
23impl<
24        P: MixedIntegerNonLinearProgram<N_X, N_Y1, N_Y2>,
25        const N_X: usize,
26        const N_Y1: usize,
27        const N_Y2: usize,
28    > BasicADProblem<N_X> for Nlp<'_, P, N_Y1, N_Y2>
29{
30    fn bounds(&self) -> ([f64; N_X], [f64; N_X]) {
31        let vars = self.minlp.x_variables().data.0[0];
32        (vars.map(|(l, _, _)| l), vars.map(|(_, u, _)| u))
33    }
34
35    fn initial_point(&self) -> [f64; N_X] {
36        self.minlp.x_variables().data.0[0].map(|(_, _, i)| i)
37    }
38
39    fn constraint_bounds(&self) -> (Vec<f64>, Vec<f64>) {
40        self.minlp
41            .constraints()
42            .into_iter()
43            .map(|c| (c.lower_bound(), c.upper_bound()))
44            .unzip()
45    }
46}
47
48impl<
49        P: MixedIntegerNonLinearProgram<N_X, N_Y1, N_Y2>,
50        const N_X: usize,
51        const N_Y1: usize,
52        const N_Y2: usize,
53    > CachedADProblem<N_X> for Nlp<'_, P, N_Y1, N_Y2>
54{
55    type Error = P::Error;
56
57    fn evaluate<D: DualNum<f64> + Copy>(&self, x: [D; N_X]) -> Result<(D, Vec<D>), P::Error> {
58        self.minlp.evaluate(SVector::from(x), self.y.map(D::from))
59    }
60}
61
62impl<
63        M: MixedIntegerNonLinearProgram<N_X, N_Y1, N_Y2>,
64        const N_X: usize,
65        const N_Y1: usize,
66        const N_Y2: usize,
67    > OuterApproximation<'_, M, N_X, N_Y1, N_Y2>
68where
69    M::Error: Debug,
70{
71    pub fn solve_nlp(
72        &mut self,
73        y: SMatrix<f64, N_Y1, N_Y2>,
74        s: Vec<f64>,
75    ) -> Option<&OptimizationResult<N_X, N_Y1, N_Y2>> {
76        #[cfg(feature = "ipopt")]
77        return self.solve_nlp_with_options(y, s, &[]);
78        #[cfg(feature = "ripopt")]
79        self.solve_nlp_with_options(y, s, &Default::default())
80    }
81
82    pub fn solve_nlp_with_options(
83        &mut self,
84        y: SMatrix<f64, N_Y1, N_Y2>,
85        s: Vec<f64>,
86        #[cfg(feature = "ipopt")] options: &[(&str, IpoptOption)],
87        #[cfg(feature = "ripopt")] options: &SolverOptions,
88    ) -> Option<&OptimizationResult<N_X, N_Y1, N_Y2>> {
89        let key = self.minlp.y_to_string(&y);
90        if self.known_solutions.contains_key(&key) {
91            return self.known_solutions.get(&key);
92        }
93
94        let optim = Nlp::new(y, self.minlp);
95        let Ok(problem) = ADProblem::new_cached(optim) else {
96            return None;
97        };
98        #[cfg(feature = "ipopt")]
99        let mut ipopt = Ipopt::new(problem).unwrap();
100        #[cfg(feature = "ipopt")]
101        let (status, x, lambda) = {
102            for &(s, o) in options {
103                ipopt.set_option(s, o);
104            }
105            let res = ipopt.solve();
106            (
107                matches!(res.status, SolveStatus::SolveSucceeded)
108                    || matches!(res.status, SolveStatus::SolvedToAcceptableLevel),
109                res.solver_data.solution.primal_variables,
110                res.solver_data.solution.constraint_multipliers,
111            )
112        };
113        #[cfg(feature = "ripopt")]
114        let res = solve(&problem, options);
115        #[cfg(feature = "ripopt")]
116        let (status, x, lambda) = {
117            (
118                matches!(res.status, ripopt::SolveStatus::Optimal)
119                    || matches!(res.status, ripopt::SolveStatus::Acceptable),
120                &res.x,
121                &res.constraint_multipliers,
122            )
123        };
124        if status {
125            let x = SVector::from_column_slice(x);
126            let (objective, constraints) = Self::gradients(self.minlp, y, x);
127            let result = OptimizationResult::new(
128                key.clone(),
129                y,
130                s,
131                objective,
132                constraints,
133                x,
134                DVector::from_column_slice(lambda),
135            );
136            self.known_solutions.insert(key.clone(), result);
137            self.known_solutions.get(&key)
138        } else {
139            // println!("{}", res.status);
140            None
141        }
142    }
143
144    fn gradients(
145        minlp: &M,
146        y: SMatrix<f64, N_Y1, N_Y2>,
147        x: SVector<f64, N_X>,
148    ) -> (Gradient<N_X, N_Y1, N_Y2>, Vec<Gradient<N_X, N_Y1, N_Y2>>) {
149        let mut y_dual = y.map(DualSVecMult::from_re);
150        let mut x_dual = x.map(DualSVecMult::from_re);
151        for (i, y) in y_dual.iter_mut().enumerate() {
152            y.eps1 = Derivative::derivative_generic(Const::<N_Y1>, Const::<N_Y2>, i);
153        }
154        for (i, x) in x_dual.iter_mut().enumerate() {
155            x.eps2 = Derivative::derivative_generic(Const::<N_X>, U1, i);
156        }
157        let (f, con) = minlp
158            .evaluate(x_dual, y_dual)
159            .expect("Unexpected error ocurred in the calculation of gradients!");
160        let f = (
161            f.re,
162            f.eps1.unwrap_generic(Const::<N_Y1>, Const::<N_Y2>),
163            f.eps2.unwrap_generic(Const::<N_X>, U1),
164        );
165        let con = con
166            .into_iter()
167            .map(|con| {
168                (
169                    con.re,
170                    con.eps1.unwrap_generic(Const::<N_Y1>, Const::<N_Y2>),
171                    con.eps2.unwrap_generic(Const::<N_X>, U1),
172                )
173            })
174            .collect();
175        (f, con)
176    }
177}
feos_campd/solver/nlp.rs

feos_campd/solver/
nlp.rs
