process_selection/
process_selection.rs

1//! Structural optimization of a process flowsheet
2//! This is the GAMS model library problem `PROCSEL`
3//! (SEQ=116), from:
4//!
5//!   Kocis & Grossmann (1987), "Relaxation Strategy for the Structural
6//!   Optimization of Process Flow Sheets", Ind. Eng. Chem. Res. 26(9),
7//!   1869-1880; also Morari & Grossmann (eds.), "Chemical Engineering
8//!   Optimization Models with GAMS" (1991).
9//!
10//! Chemical C is produced from B in unit 1. B is either purchased on the
11//! external market (`bp`) or produced from raw material A through one of two
12//! competing units (2 or 3). Binaries `y1,y2,y3` switch the three units on/off;
13//! the goal is to maximise annual profit.
14//!
15//!         A2    +-----+  B2      BP
16//!        +----->|  2  |----->+    |
17//!   A    |      +-----+      |    |  B1    +-----+    C1
18//!   ---->|                   +----+------->|  1  |-------->
19//!        |      +-----+      |             +-----+
20//!        +----->|  3  |----->+
21//!         A3    +-----+  B3
22//!
23//! The input-output relations of units 2 and 3 are the (convexified) nonlinear
24//! laws `exp(b2) - 1 = a2` and `exp(b3/1.2) - 1 = a3`, which make this an MINLP.
25//!
26//! Run with one nonlinear backend (enable exactly one, never both at once):
27//!   cargo run --example process_selection --features gams
28//!   cargo run --example process_selection --features gurobi
29
30#[cfg(any(feature = "gams", feature = "gurobi"))]
31mod model {
32    use oximo::prelude::*;
33
34    fn solve_and_report<S: Solver>(
35        label: &str,
36        mut solver: S,
37        opts: &S::Options,
38    ) -> Result<(), Box<dyn std::error::Error>> {
39        let m = Model::new("procsel");
40
41        // Positive variables (consumptions, capacities, purchases).
42        variable!(m, a2 >= 0.0);
43        variable!(m, a3 >= 0.0);
44        variable!(m, b2 >= 0.0);
45        variable!(m, b3 >= 0.0);
46        variable!(m, bp >= 0.0);
47        variable!(m, b1 >= 0.0);
48        variable!(m, 0.0 <= c1 <= 1.0);
49
50        // Binaries: existence of each process unit.
51        variable!(m, y1, Bin);
52        variable!(m, y2, Bin);
53        variable!(m, y3, Bin);
54
55        constraint!(m, inout1, c1 == 0.9 * b1);
56        constraint!(m, inout2, b2.exp() - 1.0 == a2);
57        constraint!(m, inout3, (b3 / 1.2).exp() - 1.0 == a3);
58        constraint!(m, mbalb, b1 == b2 + b3 + bp);
59        constraint!(m, log1, c1 <= 2.0 * y1);
60        constraint!(m, log2, b2 <= 4.0 * y2);
61        constraint!(m, log3, b3 <= 5.0 * y3);
62
63        // profit = sales - fixed investment - operating cost - purchases
64        objective!(
65            m,
66            Max,
67            11.0 * c1 - 3.5 * y1 - y2 - 1.5 * y3 - b2 - 1.2 * b3 - 1.8 * (a2 + a3) - 7.0 * bp
68        );
69
70        let result = solver.solve(&m, opts)?;
71
72        if let Some(obj) = result.objective() {
73            println!("--- {label} ---");
74            println!("status = {:?}", result.termination);
75            println!("profit = {obj:.4} M$/yr");
76
77            println!("units selected:");
78            for (name, y) in [("process 1", y1), ("process 2", y2), ("process 3", y3)] {
79                println!(
80                    "  {name}: {}",
81                    if (result.value_of(y).unwrap_or(0.0) - 1.0).abs() < f64::EPSILON {
82                        "ON"
83                    } else {
84                        "OFF"
85                    }
86                );
87            }
88
89            println!("flows:");
90            println!("  c1 (C produced)   = {:.4}", result.value_of(c1).unwrap_or(0.0));
91            println!("  b1 (B into 1)     = {:.4}", result.value_of(b1).unwrap_or(0.0));
92            println!(
93                "  b2, b3 (B from A) = {:.4}, {:.4}",
94                result.value_of(b2).unwrap_or(0.0),
95                result.value_of(b3).unwrap_or(0.0)
96            );
97            println!("  bp (B purchased)  = {:.4}", result.value_of(bp).unwrap_or(0.0));
98            println!(
99                "  a2, a3 (A used)   = {:.4}, {:.4}",
100                result.value_of(a2).unwrap_or(0.0),
101                result.value_of(a3).unwrap_or(0.0)
102            );
103        } else {
104            println!("--- {label} ---");
105            println!("status = {:?}", result.termination);
106            println!("no objective value");
107        }
108
109        Ok(())
110    }
111
112    #[cfg(feature = "gams")]
113    pub fn run_gams() -> Result<(), Box<dyn std::error::Error>> {
114        use oximo::gams::{GamsBaronOptions, GamsSolverConfig};
115        use oximo::solvers::Gams;
116        use std::time::Duration;
117
118        let opts = GamsOptions::default()
119            .time_limit(Duration::from_secs(120))
120            .solver(GamsSolverConfig::Baron(GamsBaronOptions::default().eps_r(1e-4)))
121            .verbose(true);
122        solve_and_report("GAMS + BARON", Gams::new(), &opts)
123    }
124
125    #[cfg(feature = "gurobi")]
126    pub fn run_gurobi() -> Result<(), Box<dyn std::error::Error>> {
127        use oximo::solvers::Gurobi;
128        use std::time::Duration;
129
130        let opts = GurobiOptions::default().time_limit(Duration::from_secs(120));
131        solve_and_report("Gurobi", Gurobi, &opts)
132    }
133
134    pub fn run() -> Result<(), Box<dyn std::error::Error>> {
135        #[cfg(feature = "gams")]
136        run_gams()?;
137        #[cfg(feature = "gurobi")]
138        run_gurobi()?;
139        Ok(())
140    }
141}
142
143#[cfg(any(feature = "gams", feature = "gurobi"))]
144fn main() -> Result<(), Box<dyn std::error::Error>> {
145    model::run()
146}
147
148#[cfg(not(any(feature = "gams", feature = "gurobi")))]
149fn main() {
150    println!("Enable a nonlinear-capable backend feature:");
151    println!("  cargo run --example process_selection --features gams");
152    println!("  cargo run --example process_selection --features gurobi");
153}
process_selection/process_selection.rs

process_selection/
process_selection.rs
