pounce_algorithm/sqp/
qp_assembly.rs

1//! Build a [`pounce_qp::QpProblem`] from the NLP linearization at
2//! the current SQP iterate `(x, λ_g)`.
3//!
4//! Standard SQP QP subproblem (Nocedal-Wright §18.1):
5//!
6//! ```text
7//!     min  ½ pᵀ ∇²L(x, λ) p + ∇f(x)ᵀ p
8//!     s.t.   bl_c ≤ c(x) + ∇c(x) p ≤ bu_c
9//!            xl − x ≤ p ≤ xu − x
10//! ```
11//!
12//! The QP's general bounds are shifted RHSs: `bl_qp = bl_c − c(x)`
13//! and `bu_qp = bu_c − c(x)` (treating equalities as
14//! `bl_c = bu_c = 0`). The QP's variable bounds are `xl − x` and
15//! `xu − x`, so the QP primal `p` directly equals the SQP step.
16//!
17//! `SqpQpData` owns the sparse storage and exposes a borrowed
18//! `QpProblem` view; this is the analog of
19//! `pounce_qp::ElasticReformulation::as_qp`.
20
21use pounce_common::types::{Index, NLP_LOWER_BOUND_INF, NLP_UPPER_BOUND_INF, Number};
22use pounce_linalg::triplet::{GenTMatrix, GenTMatrixSpace, SymTMatrix, SymTMatrixSpace};
23use pounce_qp::{HessianInertia, QpProblem};
24use std::rc::Rc;
25
26/// Owned linearization data for a single SQP iteration.
27pub struct SqpQpData {
28    pub n: usize,
29    pub m: usize,
30
31    pub h: SymTMatrix,
32    pub g: Vec<Number>,
33    pub a: GenTMatrix,
34    pub bl: Vec<Number>,
35    pub bu: Vec<Number>,
36    pub xl: Vec<Number>,
37    pub xu: Vec<Number>,
38    pub hessian_inertia: HessianInertia,
39}
40
41/// Sparse-triplet view of a derivative matrix. Indices are
42/// 1-based per the pounce-linalg convention; values are owned.
43#[derive(Clone)]
44pub struct Triplet {
45    pub n_rows: usize,
46    pub n_cols: usize,
47    pub irow: Vec<Index>,
48    pub jcol: Vec<Index>,
49    pub vals: Vec<Number>,
50}
51
52impl SqpQpData {
53    /// Assemble from concrete linearization arrays at iterate `x`.
54    ///
55    /// * `grad_f` — `∇f(x)`, length `n`.
56    /// * `c_vals` — `c(x)`, length `m` (may contain inequality
57    ///   slack values too; convention is `bl_c[i] ≤ c[i] ≤ bu_c[i]`).
58    /// * `bl_c`, `bu_c` — original NLP constraint bounds.
59    /// * `xl_orig`, `xu_orig` — original NLP variable bounds.
60    /// * `jac_c` — `∇c(x)` triplet, `m × n`.
61    /// * `hess_lag` — `∇²L(x, λ_g)` triplet, `n × n` symmetric.
62    pub fn build(
63        x: &[Number],
64        grad_f: &[Number],
65        c_vals: &[Number],
66        bl_c: &[Number],
67        bu_c: &[Number],
68        xl_orig: &[Number],
69        xu_orig: &[Number],
70        jac_c: Triplet,
71        hess_lag: Triplet,
72        hessian_inertia: HessianInertia,
73    ) -> Self {
74        let n = grad_f.len();
75        let m = c_vals.len();
76        assert_eq!(x.len(), n);
77        assert_eq!(bl_c.len(), m);
78        assert_eq!(bu_c.len(), m);
79        assert_eq!(xl_orig.len(), n);
80        assert_eq!(xu_orig.len(), n);
81        assert_eq!(jac_c.n_rows, m);
82        assert_eq!(jac_c.n_cols, n);
83        assert_eq!(hess_lag.n_rows, n);
84        assert_eq!(hess_lag.n_cols, n);
85
86        let h_space = SymTMatrixSpace::new(n as Index, hess_lag.irow, hess_lag.jcol);
87        let mut h = SymTMatrix::new(Rc::clone(&h_space));
88        h.set_values(&hess_lag.vals);
89
90        let a_space = GenTMatrixSpace::new(m as Index, n as Index, jac_c.irow, jac_c.jcol);
91        let mut a = GenTMatrix::new(Rc::clone(&a_space));
92        a.set_values(&jac_c.vals);
93
94        // QP general bounds: bl_c − c(x), bu_c − c(x) — but
95        // preserve ±∞ markers so the QP solver's one-sided
96        // ratio test still treats them as unbounded.
97        let mut bl = Vec::with_capacity(m);
98        let mut bu = Vec::with_capacity(m);
99        for i in 0..m {
100            bl.push(shift_bound(bl_c[i], c_vals[i], true));
101            bu.push(shift_bound(bu_c[i], c_vals[i], false));
102        }
103
104        // QP step bounds: xl_orig − x, xu_orig − x.
105        let mut xl = Vec::with_capacity(n);
106        let mut xu = Vec::with_capacity(n);
107        for i in 0..n {
108            xl.push(shift_bound(xl_orig[i], x[i], true));
109            xu.push(shift_bound(xu_orig[i], x[i], false));
110        }
111
112        Self {
113            n,
114            m,
115            h,
116            g: grad_f.to_vec(),
117            a,
118            bl,
119            bu,
120            xl,
121            xu,
122            hessian_inertia,
123        }
124    }
125
126    /// Borrowed `QpProblem` view ready for
127    /// `pounce_qp::QpSolver::solve`.
128    pub fn as_qp(&self) -> QpProblem<'_> {
129        QpProblem {
130            n: self.n,
131            m: self.m,
132            h: &self.h,
133            g: &self.g,
134            a: &self.a,
135            bl: &self.bl,
136            bu: &self.bu,
137            xl: &self.xl,
138            xu: &self.xu,
139            hessian_inertia: self.hessian_inertia,
140        }
141    }
142}
143
144/// Shift a bound by subtracting the current value. Preserves
145/// `NLP_*_BOUND_INF` sentinels so the QP solver's `is_finite`
146/// checks still trigger correctly.
147fn shift_bound(bound: Number, current: Number, is_lower: bool) -> Number {
148    if is_lower {
149        if bound <= NLP_LOWER_BOUND_INF {
150            NLP_LOWER_BOUND_INF
151        } else {
152            bound - current
153        }
154    } else if bound >= NLP_UPPER_BOUND_INF {
155        NLP_UPPER_BOUND_INF
156    } else {
157        bound - current
158    }
159}
pounce_algorithm/sqp/qp_assembly.rs

pounce_algorithm/sqp/
qp_assembly.rs
