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
|
# ifndef CPPAD_CPPAD_IPOPT_EXAMPLE_ODE_SIMPLE_HPP
# define CPPAD_CPPAD_IPOPT_EXAMPLE_ODE_SIMPLE_HPP
// 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 ipopt_nlp_ode_simple.hpp dev}
ODE Fitting Using Simple Representation
#######################################
{xrst_literal
// BEGIN C++
// END C++
}
{xrst_end ipopt_nlp_ode_simple.hpp}
*/
// BEGIN C++
# include "ode_problem.hpp"
// define in the empty namespace
namespace {
using namespace cppad_ipopt;
class FG_simple : public cppad_ipopt_fg_info
{
private:
bool retape_;
SizeVector N_;
SizeVector S_;
public:
// derived class part of constructor
FG_simple(bool retape_in, const SizeVector& N)
: retape_ (retape_in), N_(N)
{ assert( N_[0] == 0 );
S_.resize( N.size() );
S_[0] = 0;
for(size_t i = 1; i < N_.size(); i++)
S_[i] = S_[i-1] + N_[i];
}
// Evaluation of the objective f(x), and constraints g(x)
// using an Algorithmic Differentiation (AD) class.
ADVector eval_r(size_t not_used, const ADVector& x)
{ count_eval_r();
// temporary indices
size_t i, j, k;
// # of components of x corresponding to values for y
size_t ny_inx = (S_[Nz] + 1) * Ny;
// # of constraints (range dimension of g)
size_t m = ny_inx;
// # of components in x (domain dimension for f and g)
assert ( x.size() == ny_inx + Na );
// vector for return value
ADVector fg(m + 1);
// vector of parameters
ADVector a(Na);
for(j = 0; j < Na; j++)
a[j] = x[ny_inx + j];
// vector for value of y(t)
ADVector y(Ny);
// objective function -------------------------------
fg[0] = 0.;
for(k = 0; k < Nz; k++)
{ for(j = 0; j < Ny; j++)
y[j] = x[Ny*S_[k+1] + j];
fg[0] += eval_H<ADNumber>(k+1, y, a);
}
// initial condition ---------------------------------
ADVector F = eval_F(a);
for(j = 0; j < Ny; j++)
{ y[j] = x[j];
fg[1+j] = y[j] - F[j];
}
// trapezoidal approximation --------------------------
ADVector ym(Ny), G(Ny), Gm(Ny);
G = eval_G(y, a);
ADNumber dy;
for(k = 0; k < Nz; k++)
{ // interval between data points
Number T = s[k+1] - s[k];
// integration step size
Number dt = T / Number( N_[k+1] );
for(j = 0; j < N_[k+1]; j++)
{ size_t Index = (j + S_[k]) * Ny;
// y(t) at end of last step
ym = y;
// G(y, a) at end of last step
Gm = G;
// value of y(t) at end of this step
for(i = 0; i < Ny; i++)
y[i] = x[Ny + Index + i];
// G(y, a) at end of this step
G = eval_G(y, a);
// trapezoidal approximation residual
for(i = 0; i < Ny; i++)
{ dy = (G[i] + Gm[i]) * dt / 2;
fg[1+Ny+Index+i] =
y[i] - ym[i] - dy;
}
}
}
return fg;
}
// The operations sequence for r_eval does not depend on u,
// hence retape = false should work and be faster.
bool retape(size_t k)
{ return retape_; }
};
}
// END C++
# endif
|
