feos_campd/
lib.rs

1#[cfg(all(feature = "ipopt", feature = "ripopt"))]
2compile_error!("Features 'ipopt' and 'ripopt' cannot be enabled at the same time.");
3
4#[cfg(not(any(feature = "ipopt", feature = "ripopt")))]
5compile_error!("Either feature 'ipopt' or 'ripopt' must be enabled.");
6
7use feos::core::FeosError;
8use good_lp::{Constraint, Variable};
9use nalgebra::{SMatrix, SVector};
10use num_dual::DualNum;
11use quantity::MolarWeight;
12use std::array;
13use std::collections::{HashMap, HashSet};
14
15mod molecule;
16mod process;
17mod property;
18mod solver;
19pub use molecule::{CoMTCAMD, Disjunction, MolecularRepresentation, SuperMolecule};
20pub use process::{ContinuousVariable, ProcessModel};
21pub use property::{GcPcSaftPropertyModel, PcSaftPropertyModel, PropertyModel};
22pub use solver::{
23    GeneralConstraint, MixedIntegerNonLinearProgram, OptimizationOptions, OptimizationResult,
24    OuterApproximation,
25};
26
27/// Input for group-contribution models that allows for derivatives.
28pub struct ChemicalRecord<D> {
29    pub groups: HashMap<&'static str, D>,
30    pub bonds: HashMap<[&'static str; 2], D>,
31    pub molar_weight: MolarWeight<D>,
32}
33
34impl<D> ChemicalRecord<D> {
35    pub fn new(
36        groups: HashMap<&'static str, D>,
37        bonds: HashMap<[&'static str; 2], D>,
38        molar_weight: MolarWeight<D>,
39    ) -> Self {
40        Self {
41            groups,
42            bonds,
43            molar_weight,
44        }
45    }
46}
47
48/// A full optimization problem consisting of a [MolecularRepresentation], a [PropertyModel], and a [ProcessModel].
49#[derive(Clone)]
50pub struct IntegratedDesign<M, R, P> {
51    molecule: M,
52    property: R,
53    process: P,
54}
55
56impl<M, R, P> IntegratedDesign<M, R, P> {
57    pub fn new(molecule: M, property: R, process: P) -> Self {
58        Self {
59            molecule,
60            property,
61            process,
62        }
63    }
64}
65
66impl<
67        M: MolecularRepresentation<N_Y>,
68        R: PropertyModel<N>,
69        P: ProcessModel<N_X, N>,
70        const N_X: usize,
71        const N_Y: usize,
72        const N: usize,
73    > MixedIntegerNonLinearProgram<N_X, N_Y, N> for IntegratedDesign<M, R, P>
74{
75    type Error = FeosError;
76
77    fn x_variables(&self) -> SVector<(f64, f64, f64), N_X> {
78        SVector::from(self.process.variables().map(|v| (v.lobnd, v.upbnd, v.init)))
79    }
80
81    fn y_variables(&self) -> SMatrix<(i32, i32), N_Y, N> {
82        SMatrix::from(array::from_fn(|_| self.molecule.structure_variables()))
83    }
84
85    fn linear_constraints(&self, y: SMatrix<Variable, N_Y, N>) -> Vec<Constraint> {
86        y.data
87            .0
88            .iter()
89            .flat_map(|&y| self.molecule.constraints(y))
90            .collect()
91    }
92
93    fn constraints(&self) -> Vec<GeneralConstraint> {
94        self.process.constraints()
95    }
96
97    fn evaluate<D: DualNum<f64> + Copy>(
98        &self,
99        x: SVector<D, N_X>,
100        y: SMatrix<D, N_Y, N>,
101    ) -> Result<(D, Vec<D>), FeosError> {
102        let y_set: HashSet<_> = y.data.0.iter().map(|y| y.map(|y| y.re() as i32)).collect();
103        if y_set.len() != N {
104            Err(FeosError::IncompatibleComponents(N, y_set.len()))
105        } else {
106            let cr = y.data.0.map(|y| self.molecule.build_molecule(y));
107            let eos = self.property.build_eos(cr.each_ref());
108            self.process.evaluate(&eos, cr.each_ref(), x.data.0[0])
109        }
110    }
111
112    fn y_to_string(&self, y: &SMatrix<f64, N_Y, N>) -> String {
113        y.data.0.map(|y| self.molecule.smiles(y)).join(";")
114    }
115
116    fn exclude_solutions(&self, s: &[f64]) -> Vec<Vec<f64>> {
117        let mut exclude = vec![s.to_vec()];
118        if N == 2 {
119            let (s1, s2) = s.split_at(s.len() / 2);
120            exclude.push([s2, s1].concat())
121        }
122        exclude
123    }
124}
125
126#[cfg(test)]
127mod test {
128    use super::*;
129    use approx::assert_relative_eq;
130    use good_lp::highs;
131    use process::OrganicRankineCycle;
132
133    #[test]
134    fn test_solve() {
135        // Heptane
136        let mut y = [0.0; 28];
137        y[0] = 1.0;
138        y[6] = 2.0;
139        y[7] = 5.0;
140
141        let property =
142            PcSaftPropertyModel::full("../feos/parameters/pcsaft/rehner2023_homo.json", None)
143                .unwrap();
144
145        let campd = IntegratedDesign::new(CoMTCAMD, property, OrganicRankineCycle::default());
146        let solver = OuterApproximation::new(&campd);
147        solver.solve_ranking(SVector::from(y), highs, 5, &Default::default());
148    }
149
150    #[test]
151    fn test_solve_process() {
152        // Isopentane
153        let mut y = [0.0; 28];
154        y[0] = 1.0;
155        y[6] = 3.0;
156        y[7] = 1.0;
157        y[8] = 1.0;
158
159        let property =
160            PcSaftPropertyModel::full("../feos/parameters/pcsaft/rehner2023_homo.json", None)
161                .unwrap();
162
163        let campd = IntegratedDesign::new(CoMTCAMD, property, OrganicRankineCycle::default());
164        let mut oa = OuterApproximation::new(&campd);
165        let result = oa.solve_nlp(SVector::from(y), vec![]).unwrap();
166        println!("{:8.5} {:.5?}", result.objective.0, result.x.data.0[0]);
167    }
168
169    #[test]
170    fn test_supermolecule_fixed() {
171        const MOLECULE: SuperMolecule = SuperMolecule::alkene(0, 5);
172        const N_Y: usize = MOLECULE.variables();
173
174        // Propene
175        let mut y = [0.0; N_Y];
176        y[0] = 1.0;
177        let campd = IntegratedDesign::new(
178            MOLECULE,
179            GcPcSaftPropertyModel,
180            OrganicRankineCycle::default(),
181        );
182        let mut oa = OuterApproximation::new(&campd);
183        let result = oa.solve_nlp(SVector::from(y), vec![]).unwrap();
184        println!(
185            "{:8.5} {:.5?} {}",
186            result.objective.0,
187            result.x.data.0[0],
188            MOLECULE.smiles(&result.y.data.0[0])
189        );
190        assert_relative_eq!(result.objective.0, -0.4378352970105434, max_relative = 1e-8);
191    }
192
193    #[test]
194    fn test_supermolecule() {
195        let process = OrganicRankineCycle::default();
196        let molecule = SuperMolecule::alkene(0, 5);
197        let campd = IntegratedDesign::new(molecule, GcPcSaftPropertyModel, process);
198        let solver = OuterApproximation::new(&campd);
199
200        // Propene
201        let y = [1.0, 1.0, 1.0, 0.0, 0.0, 0.0];
202        solver.solve_ranking(SVector::from(y), highs, 1, &Default::default());
203    }
204
205    #[test]
206    fn test_supermolecule_disjunction_fixed() {
207        let process = OrganicRankineCycle::default();
208        let molecule = SuperMolecule::non_associating(0, 5);
209        let campd = IntegratedDesign::new(molecule.clone(), GcPcSaftPropertyModel, process);
210        let mut oa = OuterApproximation::new(&campd);
211
212        // Propene
213        let y0 = [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0];
214        let s = SuperMolecule::get_initial_values(&molecule, "alkene", &y0);
215        let y: [f64; 12] = array::from_fn(|i| s[i]);
216        let result = oa.solve_nlp(SVector::from(y), vec![]).unwrap();
217        println!(
218            "{:8.5} {:.5?} {}",
219            result.objective.0,
220            result.x.data.0[0],
221            molecule.smiles(&result.y.data.0[0])
222        );
223        assert_relative_eq!(result.objective.0, -0.4378352970105434, max_relative = 1e-8);
224    }
225
226    #[test]
227    fn test_supermolecule_disjunction() {
228        let process = OrganicRankineCycle::default();
229        let molecule = SuperMolecule::non_associating(0, 5);
230        let campd = IntegratedDesign::new(molecule.clone(), GcPcSaftPropertyModel, process);
231        let solver = OuterApproximation::new(&campd);
232
233        // Propene
234        let y0 = [1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0];
235        let s = SuperMolecule::get_initial_values(&molecule, "alkene", &y0);
236        let y: [f64; 12] = array::from_fn(|i| s[i]);
237        solver.solve_ranking(SVector::from(y), highs, 3, &Default::default());
238    }
239}
feos_campd/lib.rs

feos_campd/
lib.rs
