parametric_mpc/
parametric_mpc.rs

1//! Sliding-horizon parametric step — what MPC does between full
2//! re-solves.
3//!
4//! MPC typically alternates: every K steps, do a full IPM solve of the
5//! current horizon; between those, use a first-order parametric step
6//! `Δx ≈ ∂x*/∂p · Δp` to update the primal as the parameter slides.
7//! The parametric step is the cheap part because it's a back-solve
8//! against the cached factor; the IPM solve is the expensive part.
9//!
10//! This example demonstrates the **parametric-step half** of that
11//! loop: one IPM solve, followed by 10 parametric-step queries against
12//! the held factor as the parameter sweeps. Each query is a single
13//! back-solve, so the cost is essentially flat after the initial
14//! solve.
15//!
16//! What's deliberately *not* shown here: the periodic re-solve with
17//! symbolic-factor reuse across IPM runs. That's tracked by Phase 3b's
18//! `BackendPool` + `resolve()` work in
19//! `dev-notes/backend-pool-resolve.md`.
20//!
21//! Run with:
22//!   `cargo run --release --example parametric_mpc -p pounce-sensitivity`
23
24use std::cell::RefCell;
25use std::rc::Rc;
26use std::time::Instant;
27
28use pounce_algorithm::application::IpoptApplication;
29use pounce_common::types::{Index, Number};
30use pounce_nlp::tnlp::{
31    BoundsInfo, IndexStyle, IpoptCq, IpoptData, NlpInfo, Solution, SparsityRequest, StartingPoint,
32    TNLP,
33};
34use pounce_sensitivity::Solver;
35
36/// Same parametric NLP as `examples/sensitivity_session.rs`: a small
37/// 5-variable, 4-constraint problem with two pin-constraint parameters
38/// (`eta1`, `eta2`) we can sweep.
39struct ParametricTNLP {
40    eta1: Number,
41    eta2: Number,
42}
43
44impl TNLP for ParametricTNLP {
45    fn get_nlp_info(&mut self) -> Option<NlpInfo> {
46        Some(NlpInfo {
47            n: 5,
48            m: 4,
49            nnz_jac_g: 10,
50            nnz_h_lag: 5,
51            index_style: IndexStyle::C,
52        })
53    }
54    fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
55        for k in 0..3 {
56            b.x_l[k] = 0.0;
57            b.x_u[k] = 1.0e19;
58        }
59        b.x_l[3] = -1.0e19;
60        b.x_u[3] = 1.0e19;
61        b.x_l[4] = -1.0e19;
62        b.x_u[4] = 1.0e19;
63        b.g_l[0] = 0.0;
64        b.g_u[0] = 0.0;
65        b.g_l[1] = 0.0;
66        b.g_u[1] = 0.0;
67        b.g_l[2] = self.eta1;
68        b.g_u[2] = self.eta1;
69        b.g_l[3] = self.eta2;
70        b.g_u[3] = self.eta2;
71        true
72    }
73    fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
74        sp.x[0] = 0.15;
75        sp.x[1] = 0.15;
76        sp.x[2] = 0.0;
77        sp.x[3] = 0.0;
78        sp.x[4] = 0.0;
79        true
80    }
81    fn eval_f(&mut self, x: &[Number], _new_x: bool) -> Option<Number> {
82        Some(x[0] * x[0] + x[1] * x[1] + x[2] * x[2])
83    }
84    fn eval_grad_f(&mut self, x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
85        g[0] = 2.0 * x[0];
86        g[1] = 2.0 * x[1];
87        g[2] = 2.0 * x[2];
88        g[3] = 0.0;
89        g[4] = 0.0;
90        true
91    }
92    fn eval_g(&mut self, x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
93        let (x1, x2, x3, e1, e2) = (x[0], x[1], x[2], x[3], x[4]);
94        g[0] = 6.0 * x1 + 3.0 * x2 + 2.0 * x3 - e1;
95        g[1] = e2 * x1 + x2 - x3 - 1.0;
96        g[2] = e1;
97        g[3] = e2;
98        true
99    }
100    fn eval_jac_g(
101        &mut self,
102        x: Option<&[Number]>,
103        _new_x: bool,
104        mode: SparsityRequest<'_>,
105    ) -> bool {
106        match mode {
107            SparsityRequest::Structure { irow, jcol } => {
108                let rs: [Index; 10] = [0, 0, 0, 0, 1, 1, 1, 1, 2, 3];
109                let cs: [Index; 10] = [0, 1, 2, 3, 0, 1, 2, 4, 3, 4];
110                irow.copy_from_slice(&rs);
111                jcol.copy_from_slice(&cs);
112            }
113            SparsityRequest::Values { values } => {
114                let x = x.expect("Values without x");
115                values[0] = 6.0;
116                values[1] = 3.0;
117                values[2] = 2.0;
118                values[3] = -1.0;
119                values[4] = x[4];
120                values[5] = 1.0;
121                values[6] = -1.0;
122                values[7] = x[0];
123                values[8] = 1.0;
124                values[9] = 1.0;
125            }
126        }
127        true
128    }
129    fn eval_h(
130        &mut self,
131        _x: Option<&[Number]>,
132        _new_x: bool,
133        obj_factor: Number,
134        lambda: Option<&[Number]>,
135        _new_lambda: bool,
136        mode: SparsityRequest<'_>,
137    ) -> bool {
138        match mode {
139            SparsityRequest::Structure { irow, jcol } => {
140                let rs: [Index; 5] = [0, 1, 2, 4, 0];
141                let cs: [Index; 5] = [0, 1, 2, 0, 0];
142                irow.copy_from_slice(&rs);
143                jcol.copy_from_slice(&cs);
144            }
145            SparsityRequest::Values { values } => {
146                let lam = lambda.expect("Values without lambda");
147                values[0] = 2.0 * obj_factor;
148                values[1] = 2.0 * obj_factor;
149                values[2] = 2.0 * obj_factor;
150                values[3] = lam[1];
151                values[4] = 0.0;
152            }
153        }
154        true
155    }
156    fn finalize_solution(&mut self, _sol: Solution<'_>, _d: &IpoptData, _q: &IpoptCq) {}
157}
158
159fn make_app() -> IpoptApplication {
160    let mut app = IpoptApplication::new();
161    app.options_mut()
162        .set_integer_value("print_level", 0, true, false)
163        .unwrap();
164    app.options_mut()
165        .set_string_value("sb", "yes", true, false)
166        .unwrap();
167    app.initialize().unwrap();
168    app
169}
170
171fn main() {
172    let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
173        eta1: 5.0,
174        eta2: 1.0,
175    }));
176
177    // One full IPM solve at the nominal parameter.
178    let mut solver = Solver::new(make_app(), tnlp);
179    let t0 = Instant::now();
180    let status = solver.solve();
181    let solve_dt = t0.elapsed();
182    println!(
183        "IPM solve at nominal (eta1=5.0, eta2=1.0): status={status:?}, {:.3} ms",
184        solve_dt.as_secs_f64() * 1e3
185    );
186    assert!(solver.converged().is_some());
187
188    // 10 parametric steps as eta2 sweeps from 1.0 to 1.5. Each one is
189    // a back-solve against the held factor — no IPM iteration.
190    let pins = vec![2 as Index, 3];
191    let mut total_step_dt = std::time::Duration::ZERO;
192    let mut steps = Vec::new();
193    for k in 1..=10 {
194        let d_eta2 = 0.05 * k as f64;
195        let t = Instant::now();
196        let dx = solver
197            .parametric_step(&pins, &[0.0, d_eta2])
198            .expect("parametric_step ok");
199        total_step_dt += t.elapsed();
200        steps.push((d_eta2, dx));
201    }
202
203    println!("\nparametric steps against the held factor:");
204    for (d_eta2, dx) in &steps {
205        println!(
206            "  Δeta2 = {d_eta2:+.2}  -> Δx_primal = [{:+.5}, {:+.5}, {:+.5}]",
207            dx[0], dx[1], dx[2]
208        );
209    }
210    let avg_step_us = total_step_dt.as_secs_f64() * 1e6 / steps.len() as f64;
211    println!(
212        "\n10 parametric steps: total {:.3} ms, mean {avg_step_us:.1} µs/step",
213        total_step_dt.as_secs_f64() * 1e3
214    );
215    println!(
216        "\nRatio: each parametric step is roughly {:.0}x cheaper than the IPM solve.",
217        solve_dt.as_secs_f64() / (total_step_dt.as_secs_f64() / steps.len() as f64),
218    );
219}
parametric_mpc/parametric_mpc.rs

parametric_mpc/
parametric_mpc.rs
