interiors 0.1.1

Primal-Dual Interior Point Method for Nonlinear Programming
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119

use float_cmp::assert_approx_eq;
use std::iter::zip;

use full::Arr;
use sparsetools::coo::Coo;
use sparsetools::csc::CSC;
use sparsetools::csr::CSR;
use spsolve::rlu::RLU;

use crate::{nlp, Lambda, NonlinearConstraint, ObjectiveFunction, Options};

/// Hock & Schittkowski test problem #71
struct Constrained4DNonlinear {}

impl ObjectiveFunction for Constrained4DNonlinear {
    fn f(&self, x: &[f64], hessian: bool) -> (f64, Vec<f64>, Option<CSR<usize, f64>>) {
        let f = x[0] * x[3] * x[..3].iter().sum::<f64>() + x[2];
        let df = vec![
            x[0] * x[3] + x[3] * x[..3].iter().sum::<f64>(),
            x[0] * x[3],
            x[0] * x[3] + 1.0,
            x[0] * x[..3].iter().sum::<f64>(),
        ];
        if !hessian {
            (f, df, None)
        } else {
            let d2f = CSR::from_dense(&[
                vec![2.0 * x[3], x[3], x[3], 2.0 * x[0] + x[1] + x[2]],
                vec![x[3], 0.0, 0.0, x[0]],
                vec![x[3], 0.0, 0.0, x[0]],
                vec![2.0 * x[0] + x[1] + x[2], x[0], x[0], 0.0],
            ]);
            (f, df, Some(d2f))
        }
    }
}

impl NonlinearConstraint for Constrained4DNonlinear {
    fn gh(
        &self,
        x: &[f64],
        gradients: bool,
    ) -> (
        Vec<f64>,
        Vec<f64>,
        Option<CSR<usize, f64>>,
        Option<CSR<usize, f64>>,
    ) {
        let x = Arr::with_vec(x.to_vec());

        // let g = vec![x.iter().map(|&x| x.powi(2)).sum() - 40.0];
        let g = vec![x.pow(2).sum() - 40.0];
        let h = vec![-x.prod() + 25.0];
        if !gradients {
            (g, h, None, None)
        } else {
            let dg = CSC::from_dense(&[(&x * 2.0).vec()]).t();
            let dh = CSC::from_dense(&[(1.0 / &x * -x.prod()).vec()]).t();
            (h, g, Some(dh), Some(dg))
        }
    }

    fn hess(&self, x: &[f64], lam: &Lambda, sigma: f64) -> CSR<usize, f64> {
        let lambda = lam.eq_non_lin[0];
        let mu = lam.ineq_non_lin[0];

        let (_, _, d2f) = self.f(x, true);

        let l_xx: CSR<usize, f64> = d2f.unwrap() * sigma + Coo::identity(4).to_csr() * lambda * 2.0
            - CSR::from_dense(&vec![
                vec![0.0, x[2] * x[3], x[1] * x[3], x[1] * x[2]],
                vec![x[2] * x[3], 0.0, x[0] * x[3], x[0] * x[2]],
                vec![x[1] * x[3], x[0] * x[3], 0.0, x[0] * x[1]],
                vec![x[1] * x[2], x[0] * x[2], x[0] * x[1], 0.0],
            ]) * mu;

        l_xx
    }
}

#[test]
fn constrained_4d_nonlinear() {
    let f7 = Constrained4DNonlinear {};

    let x0 = vec![1.0, 5.0, 5.0, 1.0];
    let size = x0.len();
    let xmin = vec![1.0; size];
    let xmax = vec![5.0; size];

    let solver = RLU::default();
    let opt = Options::default();
    let (x, f, converged, _iterations, lambda) = nlp(
        &f7,
        &x0,
        &CSR::with_size(0, size),
        &vec![],
        &vec![],
        &xmin,
        &xmax,
        Some(&f7),
        &solver,
        &opt,
        None,
    )
    .unwrap();

    assert!(converged);
    assert_approx_eq!(f64, f, 17.0140173, epsilon = 1e-6);
    zip(x, vec![1.0, 4.7429994, 3.8211503, 1.3794082])
        .for_each(|x| assert_approx_eq!(f64, x.0, x.1, epsilon = 1e7));
    zip(lambda.eq_non_lin, vec![0.1614686])
        .for_each(|x| assert_approx_eq!(f64, x.0, x.1, epsilon = 1e7));
    zip(lambda.ineq_non_lin, vec![0.55229366])
        .for_each(|x| assert_approx_eq!(f64, x.0, x.1, epsilon = 1e8));
    zip(lambda.lower, vec![1.08787121024, 0.0, 0.0, 0.0])
        .for_each(|x| assert_approx_eq!(f64, x.0, x.1, epsilon = 1e8));
    zip(lambda.upper, vec![0.0; size])
        .for_each(|x| assert_approx_eq!(f64, x.0, x.1, epsilon = 1e8));
}