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
|
// 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_check.cpp dev}
Correctness Check for Both Simple and Fast Representations
##########################################################
{xrst_literal
// BEGIN C++
// END C++
}
{xrst_end ipopt_nlp_ode_check.cpp}
*/
// BEGIN C++
# include "ode_run.hpp"
bool ode_check(const SizeVector& N, const NumberVector& x)
{ bool ok = true;
size_t i, j;
// number of components of x corresponding to values for y
size_t ny_inx = x.size() - Na;
// compute the partial sums of the number of grid points
// and the maximum step size for the trapezoidal approximation
SizeVector S(Nz+1);
S[0] = 0;
Number max_step = 0.;
for(i = 1; i <= Nz; i++)
{ S[i] = S[i-1] + N[i];
max_step = std::max(max_step, Number(s[i] - s[i-1]) / Number(N[i]) );
}
// split out return values
NumberVector a(Na), y_0(Ny), y_1(Ny), y_2(Ny);
for(j = 0; j < Na; j++)
a[j] = x[ny_inx+j];
for(j = 0; j < Ny; j++)
{ y_0[j] = x[j];
y_1[j] = x[Ny + j];
y_2[j] = x[2 * Ny + j];
}
// Check some of the optimal a value
Number rel_tol = max_step * max_step;
Number abs_tol = rel_tol;
Number check_a[] = {a0, a1, a2}; // see the y_one function
for(j = 0; j < Na; j++)
{
ok &= CppAD::NearEqual(
check_a[j], a[j], rel_tol, abs_tol
);
}
// check accuracy of constraint equations
rel_tol = 1e-9;
abs_tol = 1e-9;
// check the initial value constraint
NumberVector F = eval_F(a);
for(j = 0; j < Ny; j++)
ok &= CppAD::NearEqual(F[j], y_0[j], rel_tol, abs_tol);
// check the first trapezoidal equation
NumberVector G_0 = eval_G(y_0, a);
NumberVector G_1 = eval_G(y_1, a);
Number dt = (s[1] - s[0]) / Number(N[1]);
Number check;
for(j = 0; j < Ny; j++)
{ check = y_1[j] - y_0[j] - (G_1[j]+G_0[j])*dt/2;
ok &= CppAD::NearEqual( check, 0., rel_tol, abs_tol);
}
//
// check the second trapezoidal equation
NumberVector G_2 = eval_G(y_2, a);
if( N[1] == 1 )
dt = (s[2] - s[1]) / Number(N[2]);
for(j = 0; j < Ny; j++)
{ check = y_2[j] - y_1[j] - (G_2[j]+G_1[j])*dt/2;
ok &= CppAD::NearEqual( check, 0., rel_tol, abs_tol);
}
//
// check the objective function (specialized to this case)
check = 0.;
NumberVector y_i(Ny);
for(size_t k = 0; k < Nz; k++)
{ for(j = 0; j < Ny; j++)
y_i[j] = x[S[k+1] * Ny + j];
check += eval_H<Number>(k + 1, y_i, a);
}
Number obj_value = 0.; // optimal object (no noise in simulation)
ok &= CppAD::NearEqual(check, obj_value, rel_tol, abs_tol);
// Use this empty namespace function to avoid warning that it is not used
static size_t ode_check_count = 0;
ode_check_count++;
ok &= count_eval_r() == ode_check_count;
return ok;
}
// END C++
|
