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
|
# ifndef CPPAD_CPPAD_IPOPT_EXAMPLE_ODE_RUN_HPP
# define CPPAD_CPPAD_IPOPT_EXAMPLE_ODE_RUN_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_run.hpp dev}
Driver for Running the Ipopt ODE Example
########################################
{xrst_literal
// BEGIN C++
// END C++
}
{xrst_end ipopt_nlp_ode_run.hpp}
*/
// BEGIN C++
# include "ode_problem.hpp"
namespace { // BEGIN empty namespace -----------------------------------------
using namespace cppad_ipopt;
template <class FG_info>
void ipopt_ode_case(
bool retape ,
const SizeVector& N ,
NumberVector& x )
{
size_t i, j;
// compute the partial sums of the number of grid points
assert( N.size() == Nz + 1);
assert( N[0] == 0 );
SizeVector S(Nz+1);
S[0] = 0;
for(i = 1; i <= Nz; i++)
S[i] = S[i-1] + N[i];
// number of components of x corresponding to values for y
size_t ny_inx = (S[Nz] + 1) * Ny;
// number of constraints (range dimension of g)
size_t m = ny_inx;
// number of components in x (domain dimension for f and g)
size_t n = ny_inx + Na;
// the argument vector for the optimization is
// y(t) at t[0] , ... , t[S[Nz]] , followed by a
NumberVector x_i(n), x_l(n), x_u(n);
for(j = 0; j < ny_inx; j++)
{ x_i[j] = 0.; // initial y(t) for optimization
x_l[j] = -1.0e19; // no lower limit
x_u[j] = +1.0e19; // no upper limit
}
for(j = 0; j < Na; j++)
{ x_i[ny_inx + j ] = .5; // initiali a for optimization
x_l[ny_inx + j ] = -1.e19; // no lower limit
x_u[ny_inx + j ] = +1.e19; // no upper
}
// all of the difference equations are constrained to the value zero
NumberVector g_l(m), g_u(m);
for(i = 0; i < m; i++)
{ g_l[i] = 0.;
g_u[i] = 0.;
}
// object defining the objective f(x) and constraints g(x)
FG_info fg_info(retape, N);
// create the CppAD Ipopt interface
cppad_ipopt_solution solution;
Ipopt::SmartPtr<Ipopt::TNLP> cppad_nlp = new cppad_ipopt_nlp(
n, m, x_i, x_l, x_u, g_l, g_u, &fg_info, &solution
);
// Create an Ipopt application
using Ipopt::IpoptApplication;
Ipopt::SmartPtr<IpoptApplication> app = new IpoptApplication();
// turn off any printing
app->Options()->SetIntegerValue("print_level", 0);
app->Options()->SetStringValue("sb", "yes");
// maximum number of iterations
app->Options()->SetIntegerValue("max_iter", 30);
// approximate accuracy in first order necessary conditions;
// see Mathematical Programming, Volume 106, Number 1,
// Pages 25-57, Equation (6)
app->Options()->SetNumericValue("tol", 1e-9);
// Derivative testing is very slow for large problems
// so comment this out if you use a large value for N[].
app->Options()-> SetStringValue( "derivative_test", "second-order");
app->Options()-> SetNumericValue( "point_perturbation_radius", 0.);
# ifndef NDEBUG
bool ok = true;
//
// Initialize the application and process the options
Ipopt::ApplicationReturnStatus status = app->Initialize();
ok &= status == Ipopt::Solve_Succeeded;
//
// Run the application
status = app->OptimizeTNLP(cppad_nlp);
ok &= status == Ipopt::Solve_Succeeded;
//
assert(ok);
# else
app->Initialize();
app->OptimizeTNLP(cppad_nlp);
# endif
// return the solution
x.resize( solution.x.size() );
for(j = 0; j < x.size(); j++)
x[j] = solution.x[j];
return;
}
} // END empty namespace ----------------------------------------------------
// END C++
# endif
|
