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
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
|
// SPDX-License-Identifier: EPL-2.0 OR GPL-2.0-or-later
// SPDX-FileCopyrightText: Bradley M. Bell <bradbell@seanet.com>
// SPDX-FileContributor: 2003-22 Bradley M. Bell
// ----------------------------------------------------------------------------
/*
{xrst_begin abs_min_quad.cpp}
abs_min_quad: Example and Test
##############################
Purpose
*******
The function
:math:`f : \B{R}^3 \rightarrow \B{R}` defined by
.. math::
f( x_0, x_1 )
=
( x_0^2 + x_1^2 ) / 2 + | x_0 - 5 | + | x_1 + 5 |
For this case, the :ref:`abs_min_quad-name` object should be equal
to the function itself.
In addition, the function is convex and
:ref:`abs_min_quad-name` should find its global minimizer.
The minimizer of this function is
:math:`x_0 = 1`, :math:`x_1 = -1`.
Source
******
{xrst_literal
// BEGIN C++
// END C++
}
{xrst_end abs_min_quad.cpp}
-------------------------------------------------------------------------------
*/
// BEGIN C++
# include <cppad/cppad.hpp>
# include "abs_min_quad.hpp"
namespace {
CPPAD_TESTVECTOR(double) join(
const CPPAD_TESTVECTOR(double)& x ,
const CPPAD_TESTVECTOR(double)& u )
{ size_t n = x.size();
size_t s = u.size();
CPPAD_TESTVECTOR(double) xu(n + s);
for(size_t j = 0; j < n; j++)
xu[j] = x[j];
for(size_t j = 0; j < s; j++)
xu[n + j] = u[j];
return xu;
}
}
bool abs_min_quad(void)
{ bool ok = true;
//
using CppAD::AD;
using CppAD::ADFun;
//
typedef CPPAD_TESTVECTOR(size_t) s_vector;
typedef CPPAD_TESTVECTOR(double) d_vector;
typedef CPPAD_TESTVECTOR( AD<double> ) ad_vector;
//
size_t level = 0; // level of tracing
size_t n = 2; // size of x
size_t m = 1; // size of y
size_t s = 2 ; // number of data points and absolute values
//
// record the function f(x)
ad_vector ad_x(n), ad_y(m);
for(size_t j = 0; j < n; j++)
ad_x[j] = double(j + 1);
Independent( ad_x );
AD<double> sum = 0.0;
sum += ad_x[0] * ad_x[0] / 2.0 + abs( ad_x[0] - 5 );
sum += ad_x[1] * ad_x[1] / 2.0 + abs( ad_x[1] + 5 );
ad_y[0] = sum;
ADFun<double> f(ad_x, ad_y);
// create its abs_normal representation in g, a
ADFun<double> g, a;
f.abs_normal_fun(g, a);
// check dimension of domain and range space for g
ok &= g.Domain() == n + s;
ok &= g.Range() == m + s;
// check dimension of domain and range space for a
ok &= a.Domain() == n;
ok &= a.Range() == s;
// --------------------------------------------------------------------
// Choose the point x_hat = 0
d_vector x_hat(n);
for(size_t j = 0; j < n; j++)
x_hat[j] = 0.0;
// value of a_hat = a(x_hat)
d_vector a_hat = a.Forward(0, x_hat);
// (x_hat, a_hat)
d_vector xu_hat = join(x_hat, a_hat);
// value of g[ x_hat, a_hat ]
d_vector g_hat = g.Forward(0, xu_hat);
// Jacobian of g[ x_hat, a_hat ]
d_vector g_jac = g.Jacobian(xu_hat);
// trust region bound
d_vector bound(n);
for(size_t j = 0; j < n; j++)
bound[j] = 10.0;
// convergence criteria
d_vector epsilon(2);
double eps99 = 99.0 * std::numeric_limits<double>::epsilon();
epsilon[0] = eps99;
epsilon[1] = eps99;
// maximum number of iterations
s_vector maxitr(2);
maxitr[0] = 10; // maximum number of abs_min_quad iterations
maxitr[1] = 35; // maximum number of qp_interior iterations
// set Hessian equal to identity matrix I
d_vector hessian(n * n);
for(size_t i = 0; i < n; i++)
{ for(size_t j = 0; j < n; j++)
hessian[i * n + j] = 0.0;
hessian[i * n + i] = 1.0;
}
// minimize the approxiamtion for f (which is equal to f for this case)
d_vector delta_x(n);
ok &= CppAD::abs_min_quad(
level, n, m, s,
g_hat, g_jac, hessian, bound, epsilon, maxitr, delta_x
);
// check that the solution
ok &= CppAD::NearEqual( delta_x[0], +1.0, eps99, eps99 );
ok &= CppAD::NearEqual( delta_x[1], -1.0, eps99, eps99 );
return ok;
}
// END C++
|
