pounce_algorithm/sqp/
warm_start.rs

1//! Warm-start helpers — building a [`pounce_qp::WorkingSet`] from
2//! a converged IPM (or any) iterate so the next SQP solve can pick
3//! up where the IPM left off (Phase 5c §7.5 + sensitivity
4//! corrector handoff).
5//!
6//! The classifier is the **multiplier-sign + primal-distance**
7//! heuristic standard in mixed IPM/SQP warm-start pipelines
8//! (Wächter-Biegler 2006 §6; Forsgren-Gill-Wright 2002 §5). It is
9//! intentionally lossy at degenerate active sets — the QP solver
10//! will detect and correct any misclassification in the first
11//! step of the next QP, so correctness is preserved.
12
13use pounce_common::types::{NLP_LOWER_BOUND_INF, NLP_UPPER_BOUND_INF};
14use pounce_common::Number;
15use pounce_qp::{BoundStatus, ConsStatus, WorkingSet};
16
17/// Classify the active set at iterate `(x, λ_x, λ_g)` against the
18/// supplied bounds and constraint-bound vectors.
19///
20/// Inputs:
21/// - `lambda_x`: packed signed bound multipliers (`z_l − z_u`) of
22///   length `n`. Positive ⇒ lower bound active; negative ⇒ upper.
23/// - `lambda_g`: stacked constraint multipliers `[y_c ; y_d]` of
24///   length `m = m_eq + m_ineq`. Same sign convention.
25/// - `m_eq`: number of equality rows at the start of `lambda_g`.
26///   Used to flag rows as [`ConsStatus::Equality`] without
27///   consulting `g_l`/`g_u`.
28/// - `x`, `x_l`, `x_u`: primal iterate and variable bounds, length
29///   `n`. The bound-classifier double-checks the primal is close
30///   to the bound (within `primal_tol`) — guards against the case
31///   where a multiplier is large but the primal hasn't actually
32///   reached the bound (e.g. near-degenerate KKT or a bad
33///   multiplier estimate).
34/// - `g`, `g_l`, `g_u`: constraint values and bounds, length `m`.
35///   Used identically for constraint rows.
36/// - `mult_tol`: multiplier-magnitude threshold; a row whose
37///   `|λ|` falls below this is classified as `Inactive`
38///   regardless of primal distance.
39/// - `primal_tol`: distance threshold between `x[i]` and `x_l[i]`
40///   / `x_u[i]` (resp. `g[i]` vs `g_l[i]` / `g_u[i]`) below which
41///   a row is treated as "at the bound".
42///
43/// Variable bounds with `x_l[i] == x_u[i]` are classified
44/// [`BoundStatus::Fixed`]; constraint rows in the first `m_eq`
45/// slots are [`ConsStatus::Equality`]. Both are unconditionally
46/// active.
47#[allow(clippy::too_many_arguments)]
48pub fn classify_working_set(
49    lambda_x: &[Number],
50    lambda_g: &[Number],
51    m_eq: usize,
52    x: &[Number],
53    x_l: &[Number],
54    x_u: &[Number],
55    g: &[Number],
56    g_l: &[Number],
57    g_u: &[Number],
58    mult_tol: Number,
59    primal_tol: Number,
60) -> WorkingSet {
61    let n = lambda_x.len();
62    let m = lambda_g.len();
63    debug_assert_eq!(x.len(), n);
64    debug_assert_eq!(x_l.len(), n);
65    debug_assert_eq!(x_u.len(), n);
66    debug_assert_eq!(g.len(), m);
67    debug_assert_eq!(g_l.len(), m);
68    debug_assert_eq!(g_u.len(), m);
69    debug_assert!(m_eq <= m);
70
71    // Bound-finiteness uses the same `NLP_*_BOUND_INF` sentinels
72    // pounce uses everywhere else (default ±1e19). Naive
73    // `.is_finite()` would falsely include `−1e19` as a real lower
74    // bound and tag any unbounded variable at that value as
75    // `AtLower` (PR #50 review A4).
76    let mut bounds = Vec::with_capacity(n);
77    for i in 0..n {
78        let lo_fin = x_l[i] > NLP_LOWER_BOUND_INF;
79        let up_fin = x_u[i] < NLP_UPPER_BOUND_INF;
80        if lo_fin && up_fin && (x_u[i] - x_l[i]).abs() < primal_tol {
81            bounds.push(BoundStatus::Fixed);
82            continue;
83        }
84        let mu = lambda_x[i];
85        let at_lo = lo_fin && (x[i] - x_l[i]).abs() < primal_tol;
86        let at_up = up_fin && (x_u[i] - x[i]).abs() < primal_tol;
87        let status = if mu > mult_tol && at_lo {
88            BoundStatus::AtLower
89        } else if mu < -mult_tol && at_up {
90            BoundStatus::AtUpper
91        } else if at_lo && mu >= 0.0 {
92            BoundStatus::AtLower
93        } else if at_up && mu <= 0.0 {
94            BoundStatus::AtUpper
95        } else {
96            BoundStatus::Inactive
97        };
98        bounds.push(status);
99    }
100
101    let mut constraints = Vec::with_capacity(m);
102    for i in 0..m {
103        if i < m_eq {
104            constraints.push(ConsStatus::Equality);
105            continue;
106        }
107        let lo_fin = g_l[i] > NLP_LOWER_BOUND_INF;
108        let up_fin = g_u[i] < NLP_UPPER_BOUND_INF;
109        if lo_fin && up_fin && (g_u[i] - g_l[i]).abs() < primal_tol {
110            constraints.push(ConsStatus::Equality);
111            continue;
112        }
113        let mu = lambda_g[i];
114        let at_lo = lo_fin && (g[i] - g_l[i]).abs() < primal_tol;
115        let at_up = up_fin && (g_u[i] - g[i]).abs() < primal_tol;
116        let status = if mu > mult_tol && at_lo {
117            ConsStatus::AtLower
118        } else if mu < -mult_tol && at_up {
119            ConsStatus::AtUpper
120        } else if at_lo && mu >= 0.0 {
121            ConsStatus::AtLower
122        } else if at_up && mu <= 0.0 {
123            ConsStatus::AtUpper
124        } else {
125            ConsStatus::Inactive
126        };
127        constraints.push(status);
128    }
129
130    WorkingSet {
131        bounds,
132        constraints,
133    }
134}
135
136#[cfg(test)]
137mod tests {
138    use super::*;
139
140    #[test]
141    fn classify_treats_nlp_bound_inf_sentinel_as_unbounded() {
142        // PR #50 review A4 regression. Variables with `x_l =
143        // NLP_LOWER_BOUND_INF` (the −1e19 sentinel) are unbounded
144        // below; even a primal value at exactly that sentinel must
145        // be tagged `Inactive`, not `AtLower`. Prior to the fix
146        // `is_finite()` would treat `−1e19` as a real bound.
147        let ws = classify_working_set(
148            &[0.0],
149            &[],
150            0,
151            &[-1.0e19],
152            &[NLP_LOWER_BOUND_INF],
153            &[NLP_UPPER_BOUND_INF],
154            &[],
155            &[],
156            &[],
157            1e-8,
158            1e-6,
159        );
160        assert_eq!(ws.bounds[0], BoundStatus::Inactive);
161    }
162
163    #[test]
164    fn classify_all_inactive_when_strictly_interior() {
165        // 1-D unconstrained, x* in the interior, no multipliers.
166        let ws = classify_working_set(
167            &[0.0],
168            &[],
169            0,
170            &[0.5],
171            &[-1.0],
172            &[1.0],
173            &[],
174            &[],
175            &[],
176            1e-8,
177            1e-8,
178        );
179        assert_eq!(ws.bounds[0], BoundStatus::Inactive);
180        assert!(ws.constraints.is_empty());
181    }
182
183    #[test]
184    fn classify_lower_bound_active_when_primal_at_bound_and_mult_positive() {
185        let ws = classify_working_set(
186            &[2.0],
187            &[],
188            0,
189            &[0.0],
190            &[0.0],
191            &[1.0],
192            &[],
193            &[],
194            &[],
195            1e-8,
196            1e-8,
197        );
198        assert_eq!(ws.bounds[0], BoundStatus::AtLower);
199    }
200
201    #[test]
202    fn classify_upper_bound_active_when_primal_at_bound_and_mult_negative() {
203        let ws = classify_working_set(
204            &[-2.0],
205            &[],
206            0,
207            &[1.0],
208            &[0.0],
209            &[1.0],
210            &[],
211            &[],
212            &[],
213            1e-8,
214            1e-8,
215        );
216        assert_eq!(ws.bounds[0], BoundStatus::AtUpper);
217    }
218
219    #[test]
220    fn classify_fixed_when_bounds_equal() {
221        let ws = classify_working_set(
222            &[0.0],
223            &[],
224            0,
225            &[2.0],
226            &[2.0],
227            &[2.0],
228            &[],
229            &[],
230            &[],
231            1e-8,
232            1e-8,
233        );
234        assert_eq!(ws.bounds[0], BoundStatus::Fixed);
235    }
236
237    #[test]
238    fn classify_equality_constraint_always_active() {
239        // 1 eq constraint at row 0, no ineqs.
240        let ws = classify_working_set(
241            &[],
242            &[1.0],
243            1,
244            &[],
245            &[],
246            &[],
247            &[5.0],
248            &[5.0],
249            &[5.0],
250            1e-8,
251            1e-8,
252        );
253        assert_eq!(ws.constraints[0], ConsStatus::Equality);
254    }
255
256    #[test]
257    fn classify_inequality_at_lower_bound() {
258        let ws = classify_working_set(
259            &[],
260            &[3.0],
261            0,
262            &[],
263            &[],
264            &[],
265            &[1.0],
266            &[1.0],
267            &[10.0],
268            1e-8,
269            1e-8,
270        );
271        assert_eq!(ws.constraints[0], ConsStatus::AtLower);
272    }
273
274    #[test]
275    fn classify_inequality_at_upper_bound() {
276        let ws = classify_working_set(
277            &[],
278            &[-3.0],
279            0,
280            &[],
281            &[],
282            &[],
283            &[10.0],
284            &[0.0],
285            &[10.0],
286            1e-8,
287            1e-8,
288        );
289        assert_eq!(ws.constraints[0], ConsStatus::AtUpper);
290    }
291
292    #[test]
293    fn classify_inactive_when_primal_off_bound_despite_large_multiplier() {
294        // Bound multiplier is large but primal is mid-range —
295        // tag as Inactive, not AtLower. This guards against
296        // stale-multiplier carry from a slightly mis-aligned
297        // perturbation.
298        let ws = classify_working_set(
299            &[2.0],
300            &[],
301            0,
302            &[0.5],
303            &[0.0],
304            &[1.0],
305            &[],
306            &[],
307            &[],
308            1e-8,
309            1e-8,
310        );
311        assert_eq!(ws.bounds[0], BoundStatus::Inactive);
312    }
313}
pounce_algorithm/sqp/warm_start.rs

pounce_algorithm/sqp/
warm_start.rs
