Nonlinear Programming
Mathwrist Code Example Series

Copyright ©Mathwrist LLC 2023
(October 5, 2023)
Abstract

This document gives some code examples on how to use the nonlinear programming (NLP) features in Mathwrist’s C++ Numerical Programming Library (NPL).

Disk-constrained Rosenbrock

This example illustrates minimizing the Rosenbrock’s function f⁢(x) subject to unit disk constraint x12+x2≤1 using the sequential quadratic programming (SQP) solver.

f⁢(𝐱)=100⁢(x2-x12)2+(1-x1)2
1 // Objective function.
2 struct Rosenbrock : public FunctionND
3 {
4   FunctionND* clone() const { return 0; }
6   double eval(const Matrix& x, Matrix* _g=0, Matrix* _H=0) const
7   {
8     double x1 = x(0, 0);
9     double x2 = x(1, 0);
10     double temp1 = x2 - x1 * x1;
11     const double retval = 100 * temp1 * temp1 +
12       (1 - x1) * (1 - x1);
14     if (_g != 0)
15     {
16       (*_g)(0, 0) = 400 * (x1 * x1 * x1 - x1 * x2) + 2 * x1 - 2;
17       (*_g)(1, 0) = 200 * (x2 - x1 * x1);
18     }
20     if (_H != 0)
21     {
22       (*_H)(0, 0) = 1200 * x1 * x1 - 400 * x2 + 2;
23       (*_H)(0, 1) = -400 * x1;
24       (*_H)(1, 0) = -400 * x1;
25       (*_H)(1, 1) = 200;
26     }
28     return retval;
29   }
31   size_t dimension_size() const { return 2; }
32 };
34 // unit disk constraint x1^2 + x2^2 <= 1.
35 struct FDisk : public FunctionND
36 {
37   FunctionND* clone() const { return 0; }
39   double eval(const Matrix& x, Matrix* _g=0, Matrix* _H=0) const
40   {
41     double x1 = x(0, 0);
42     double x2 = x(1, 0);
43     const double retval = x1 * x1 + x2 * x2;
45     if (_g != 0)
46     {
47       (*_g)(0, 0) = 2 * x1;
48       (*_g)(1, 0) = 2 * x2;
49     }
51     if (_H != 0)
52     {
53       (*_H)(0, 0) = 2;
54       (*_H)(0, 1) = 0;
55       (*_H)(1, 0) = 0;
56       (*_H)(1, 1) = 2;
57     }
59     return retval;
60   }
62   size_t dimension_size() const { return 2; }
63 };
65 Rosenbrock rosenbrock;
67 // Nonlinear constraint functions.
68 VtrValueFunctionND_M cx(1);
69 cx[0] = FunctionND::Handle(new FDisk);
71 NonlinearProg::Problem nlp(rosenbrock, cx);
73 // Unit disk bounded below 1.
74 nlp.bounds_nonlinear[0] = Bound::upper(1);
76 // Initial guess x = 0
77 Matrix x(2, 1);
79 SQP_ActiveSet sqp(nlp);
80 sqp.set_max_iter(100);
82 sqp(x);
84 // Compare with known solution.
85 assert(std::abs(x(0, 0) - 0.786415) < 1.e-6);
86 assert(std::abs(x(1, 0) - 0.617698) < 1.e-6);