reaction_path/
reaction_path.rs

1//! Logical inference for reaction path synthesis (MILP).
2//!
3//! Given 22 possible chemical reactions over 34 chemicals, determine whether
4//! acetone (y06, ch3coch3) can be synthesized from a fixed set of raw materials
5//! and catalysts.
6//!
7//! Binary variable `y[v] = 1` if chemical `v` is synthesizable. Available raw
8//! materials are fixed to 1, unavailable chemicals to 0. For each reaction that
9//! can produce `v` from reactants `vv`:
10//!
11//! ```text
12//! sum_vv (1 - y[vv]) >= 1 - y[v]
13//! ```
14//!
15//! This forces: if all reactants are present, the product must be present.
16//! Minimizing `y[y06]` reveals whether acetone is synthesizable (optimal = 1).
17//!
18//! This example is translated from GAMS code in the REACTION model (SEQ=121).
19//!
20//! Reference:
21//! Raman, R, and Grossmann, I E, "Relation between MINLP Modeling
22//! and Logical Inference for Chemical Process Synthesis", Computers and
23//! Chemical Engineering 15, 2 (1991), 73–84.
24//!
25//! Run with HiGHS (default):
26//! ```text
27//! cargo run --example reaction_path
28//! ```
29//!
30//! Run with GAMS / CPLEX:
31//! ```text
32//! cargo run --example reaction_path --features gams
33//! ```
34//! Requires a licensed GAMS installation with CPLEX on PATH.
35
36#[cfg(any(feature = "gams", feature = "highs"))]
37use oximo::prelude::*;
38
39#[cfg(feature = "gams")]
40use oximo::gams::{GamsCplexOptions, GamsSolverConfig};
41#[cfg(feature = "gams")]
42use oximo::solvers::Gams;
43
44#[cfg(all(feature = "highs", not(feature = "gams")))]
45use oximo::solvers::Highs;
46
47#[cfg(any(feature = "gams", feature = "highs"))]
48#[expect(clippy::cast_precision_loss)]
49fn main() -> Result<(), Box<dyn std::error::Error>> {
50    // 34 chemicals: index = yXX - 1 (y01 -> 0, ..., y34 -> 33).
51    const CHEMICALS: [&str; 34] = [
52        "y01", "y02", "y03", "y04", "y05", "y06", "y07", "y08", "y09", "y10", "y11", "y12", "y13",
53        "y14", "y15", "y16", "y17", "y18", "y19", "y20", "y21", "y22", "y23", "y24", "y25", "y26",
54        "y27", "y28", "y29", "y30", "y31", "y32", "y33", "y34",
55    ];
56
57    // (reaction label, product index, reactant indices) - all 0-based.
58    // Transcribed from logicc set in GAMS REACTION model (SEQ=121).
59    let logicc: &[(&str, usize, &[usize])] = &[
60        ("rxn01", 3, &[0, 1, 2]),      // y04 <- y01 + y02 + y03
61        ("rxn02", 5, &[3, 4]),         // y06 <- y04 + y05
62        ("rxn03", 6, &[3, 4]),         // y07 <- y04 + y05
63        ("rxn04", 2, &[3, 4]),         // y03 <- y04 + y05
64        ("rxn05", 10, &[7, 8, 9]),     // y11 <- y08 + y09 + y10
65        ("rxn06", 5, &[10, 11, 12]),   // y06 <- y11 + y12 + y13
66        ("rxn07", 14, &[13, 8, 9, 4]), // y15 <- y14 + y09 + y10 + y05
67        ("rxn08", 5, &[14, 15, 16]),   // y06 <- y15 + y16 + y17
68        ("rxn09", 5, &[17, 18, 11]),   // y06 <- y18 + y19 + y12
69        ("rxn10", 19, &[17, 18, 11]),  // y20 <- y18 + y19 + y12
70        ("rxn11", 8, &[20, 21]),       // y09 <- y21 + y22
71        ("rxn12", 23, &[8, 22]),       // y24 <- y09 + y23
72        ("rxn13", 17, &[23, 16]),      // y18 <- y24 + y17
73        ("rxn14", 20, &[24, 25]),      // y21 <- y25 + y26
74        ("rxn15", 26, &[24, 25]),      // y27 <- y25 + y26
75        ("rxn16", 13, &[2, 27, 28]),   // y14 <- y03 + y28 + y29
76        ("rxn17", 31, &[29, 30, 11]),  // y32 <- y30 + y31 + y12
77        ("rxn18", 7, &[29, 30, 11]),   // y08 <- y30 + y31 + y12
78        ("rxn19", 29, &[24, 32]),      // y30 <- y25 + y33
79        ("rxn20", 12, &[24, 32]),      // y13 <- y25 + y33
80        ("rxn21", 0, &[33, 2]),        // y01 <- y34 + y03
81        ("rxn22", 33, &[13, 27]),      // y34 <- y14 + y28
82    ];
83
84    // y02, y03, y05, y10, y12, y13, y17, y22, y25, y26, y28, y31, y33: fixed to 1.
85    let available: &[usize] = &[1, 2, 4, 9, 11, 12, 16, 21, 24, 25, 27, 30, 32];
86    // y16, y19: fixed to 0.
87    let unavailable: &[usize] = &[15, 18];
88
89    let m = Model::new("reaction_path");
90    let chemicals = Set::strings(CHEMICALS);
91
92    variable!(m, y[v in chemicals], Bin);
93    // Fix availability: raw materials/catalysts to 1, unavailable chemicals to 0.
94    for &i in available {
95        m.fix(y[CHEMICALS[i]], 1.0);
96    }
97    for &i in unavailable {
98        m.fix(y[CHEMICALS[i]], 0.0);
99    }
100
101    // sum_vv (1 - y[vv]) >= 1 - y[v]
102    //    <=>  y[v] - sum_vv y[vv] >= 1 - |reactants|
103    for &(_rx, prod, reactants) in logicc {
104        let n = reactants.len() as f64;
105        constraint!(m, y[CHEMICALS[prod]] - sum!(y[CHEMICALS[vv]] for vv in reactants) >= 1.0 - n);
106    }
107
108    objective!(m, Min, y["y06"]); // acetone
109
110    #[cfg(feature = "gams")]
111    let result = {
112        let opts = GamsOptions::default()
113            .time_limit(std::time::Duration::from_secs(60))
114            .solver(GamsSolverConfig::Cplex(GamsCplexOptions::default()))
115            .verbose(true);
116        let mut solver = Gams::new();
117        solver.solve(&m, &opts)?
118    };
119
120    #[cfg(all(feature = "highs", not(feature = "gams")))]
121    let result = Highs.solve(&m, &HighsOptions::default().verbose(true))?;
122
123    println!("Status : {:?}", result.termination);
124    if let Some(obj) = result.objective() {
125        println!(
126            "Acetone (y06, ch3coch3): {}",
127            if (obj - 1.0).abs() < 1e-6 { "Synthesizable" } else { "Not synthesizable" }
128        );
129    }
130
131    let synthesizable: Vec<&str> = CHEMICALS
132        .iter()
133        .copied()
134        .filter(|name| (result.value_of(y[*name]).unwrap_or(0.0) - 1.0).abs() < 1e-6)
135        .collect();
136    if !synthesizable.is_empty() {
137        println!("Synthesizable chemicals: {}", synthesizable.join(", "));
138    }
139
140    Ok(())
141}
142
143#[cfg(not(any(feature = "gams", feature = "highs")))]
144fn main() {
145    println!("Enable at least one solver feature:");
146    println!("  cargo run --example reaction_path                  # HiGHS (default)");
147    println!("  cargo run --example reaction_path --features gams  # GAMS/CPLEX");
148}
reaction_path/reaction_path.rs

reaction_path/
reaction_path.rs
