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, Number, NLP_LOWER_BOUND_INF, NLP_UPPER_BOUND_INF};
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.
43pub struct Triplet {
44    pub n_rows: usize,
45    pub n_cols: usize,
46    pub irow: Vec<Index>,
47    pub jcol: Vec<Index>,
48    pub vals: Vec<Number>,
49}
50
51impl SqpQpData {
52    /// Assemble from concrete linearization arrays at iterate `x`.
53    ///
54    /// * `grad_f` — `∇f(x)`, length `n`.
55    /// * `c_vals` — `c(x)`, length `m` (may contain inequality
56    ///   slack values too; convention is `bl_c[i] ≤ c[i] ≤ bu_c[i]`).
57    /// * `bl_c`, `bu_c` — original NLP constraint bounds.
58    /// * `xl_orig`, `xu_orig` — original NLP variable bounds.
59    /// * `jac_c` — `∇c(x)` triplet, `m × n`.
60    /// * `hess_lag` — `∇²L(x, λ_g)` triplet, `n × n` symmetric.
61    pub fn build(
62        x: &[Number],
63        grad_f: &[Number],
64        c_vals: &[Number],
65        bl_c: &[Number],
66        bu_c: &[Number],
67        xl_orig: &[Number],
68        xu_orig: &[Number],
69        jac_c: Triplet,
70        hess_lag: Triplet,
71        hessian_inertia: HessianInertia,
72    ) -> Self {
73        let n = grad_f.len();
74        let m = c_vals.len();
75        assert_eq!(x.len(), n);
76        assert_eq!(bl_c.len(), m);
77        assert_eq!(bu_c.len(), m);
78        assert_eq!(xl_orig.len(), n);
79        assert_eq!(xu_orig.len(), n);
80        assert_eq!(jac_c.n_rows, m);
81        assert_eq!(jac_c.n_cols, n);
82        assert_eq!(hess_lag.n_rows, n);
83        assert_eq!(hess_lag.n_cols, n);
84
85        let h_space = SymTMatrixSpace::new(n as Index, hess_lag.irow, hess_lag.jcol);
86        let mut h = SymTMatrix::new(Rc::clone(&h_space));
87        h.set_values(&hess_lag.vals);
88
89        let a_space = GenTMatrixSpace::new(m as Index, n as Index, jac_c.irow, jac_c.jcol);
90        let mut a = GenTMatrix::new(Rc::clone(&a_space));
91        a.set_values(&jac_c.vals);
92
93        // QP general bounds: bl_c − c(x), bu_c − c(x) — but
94        // preserve ±∞ markers so the QP solver's one-sided
95        // ratio test still treats them as unbounded.
96        let mut bl = Vec::with_capacity(m);
97        let mut bu = Vec::with_capacity(m);
98        for i in 0..m {
99            bl.push(shift_bound(bl_c[i], c_vals[i], true));
100            bu.push(shift_bound(bu_c[i], c_vals[i], false));
101        }
102
103        // QP step bounds: xl_orig − x, xu_orig − x.
104        let mut xl = Vec::with_capacity(n);
105        let mut xu = Vec::with_capacity(n);
106        for i in 0..n {
107            xl.push(shift_bound(xl_orig[i], x[i], true));
108            xu.push(shift_bound(xu_orig[i], x[i], false));
109        }
110
111        Self {
112            n,
113            m,
114            h,
115            g: grad_f.to_vec(),
116            a,
117            bl,
118            bu,
119            xl,
120            xu,
121            hessian_inertia,
122        }
123    }
124
125    /// Borrowed `QpProblem` view ready for
126    /// `pounce_qp::QpSolver::solve`.
127    pub fn as_qp(&self) -> QpProblem<'_> {
128        QpProblem {
129            n: self.n,
130            m: self.m,
131            h: &self.h,
132            g: &self.g,
133            a: &self.a,
134            bl: &self.bl,
135            bu: &self.bu,
136            xl: &self.xl,
137            xu: &self.xu,
138            hessian_inertia: self.hessian_inertia,
139        }
140    }
141}
142
143/// Shift a bound by subtracting the current value. Preserves
144/// `NLP_*_BOUND_INF` sentinels so the QP solver's `is_finite`
145/// checks still trigger correctly.
146fn shift_bound(bound: Number, current: Number, is_lower: bool) -> Number {
147    if is_lower {
148        if bound <= NLP_LOWER_BOUND_INF {
149            NLP_LOWER_BOUND_INF
150        } else {
151            bound - current
152        }
153    } else if bound >= NLP_UPPER_BOUND_INF {
154        NLP_UPPER_BOUND_INF
155    } else {
156        bound - current
157    }
158}
pounce_algorithm/sqp/qp_assembly.rs

pounce_algorithm/sqp/
qp_assembly.rs
