Notation
The table below contains
the name of a variable,
the meaning of the variable value,
and the value for this particular example.
If the value is not specified in the table below,
the corresponding value in ipopt_nlp_ode_problem.hpp
can be changed and the example should still run
(with no other changes).
Name
Meaning
Value
@(@
Na
@)@
number of parameters to fit
3
@(@
Ny
@)@
number components in ODE
2
@(@
Nz
@)@
number of measurements
4
@(@
N(i)
@)@
number of grid points between i-1-th and i-th measurement
@(@
S(i)
@)@
number of grid points up to an including the i-th measurement
Forward Problem
We consider the following ordinary differential equation:
@[@
\begin{array}{rcl}
\partial_t y_0 ( t , a ) & = & - a_1 * y_0 (t, a )
\\
\partial_t y_1 (t , a ) & = & + a_1 * y_0 (t, a ) - a_2 * y_1 (t, a )
\end{array}
@]@
with the initial conditions
@[@
y_0 (0 , a) = F(a) = \left( \begin{array}{c} a_0 \\ 0 \end{array} \right)
@]@
where @(@
Na
@)@ is the number of parameters,
@(@
a \in \B{R}^{Na}
@)@ is an unknown parameter vector.
The function and @(@
F : \B{R}^{Na} \rightarrow \B{R}^{Ny}
@)@
is defined by the equation above
where @(@
Ny
@)@ is the number of components in @(@
y(t, a)
@)@.
Our forward problem is stated as follows: Given @(@
a \in \B{R}^{Na}
@)@
determine the value of @(@
y ( t , a )
@)@, for @(@
t \in R
@)@, that
solves the initial value problem above.
Measurements
We use @(@
Nz
@)@ to denote the number of measurements.
Suppose we are also given a measurement vector @(@
z \in \B{R}^{Nz}
@)@
and for @(@
i = 1, \ldots, Nz
@)@,
we model @(@
z_i
@)@ by
@[@
z_i = y_1 ( s_i , a) + e_i
@]@
where @(@
s_i \in \B{R}
@)@ is the time for the i-th measurement,
@(@
e_i \sim {\bf N} (0 , \sigma^2 )
@)@ is the corresponding noise,
and @(@
\sigma \in \B{R}_+
@)@ is the corresponding standard deviation.
Simulation Analytic Solution
The following analytic solution to the forward problem is used
to simulate a data set:
@[@
\begin{array}{rcl}
y_0 (t , a) & = & a_0 * \exp( - a_1 * t )
\\
y_1 (t , a) & = &
a_0 * a_1 * \frac{\exp( - a_2 * t ) - \exp( -a_1 * t )}{ a_1 - a_2 }
\end{array}
@]@
time corresponding to the i-th measurement,
@(@
i = 1 , \ldots , Nz
@)@
Simulated Measurement Values
The simulated measurement values are given by the equation
@[@
\begin{array}{rcl}
z_i
& = & e_i + y_1 ( s_i , \bar{a} )
\\
& = &
e_i + \bar{a}_0 * \bar{a}_1 *
\frac{\exp( - \bar{a}_2 * s_i ) - \exp( -\bar{a}_1 * s_i )}
{ \bar{a}_1 - \bar{a}_2 }
\end{array}
@]@
for @(@
k = 1, \ldots , Nz
@)@.
Inverse Problem
The maximum likelihood estimate for @(@
a
@)@ given @(@
z
@)@
solves the following inverse problem
@[@
\begin{array}{rcl}
{\rm minimize} \;
& \sum_{i=1}^{Nz} H_i [ y( s_i , a ) , a ]
& \;{\rm w.r.t} \; a \in \B{R}^{Na}
\end{array}
@]@
where the functions
@(@
H_i : \B{R}^{Ny} \times \B{R}^{Na} \rightarrow \B{R}
@)@ is
defined by
@[@
H_i (y, a) = ( z_i - y_1 )^2
@]@
Trapezoidal Approximation
This example uses a trapezoidal approximation to solve the ODE.
This approximation procedures starts with
@[@
y^0 = y(0, a) = \left( \begin{array}{c} a_0 \\ 0 \end{array} \right)
@]@
Given a time grid @(@
\{ t_i \}
@)@ and
an approximate value @(@
y^{i-1}
@)@ for @(@
y ( t_{i-1} , a )
@)@,
the a trapezoidal method approximates @(@
y ( t_i , a )
@)@
(denoted by @(@
y^i
@)@ ) by solving the equation
@[@
y^i = y^{i-1} +
\left[ G( y^i , a ) + G( y^{i-1} , a ) \right] * \frac{t_i - t_{i-1} }{ 2 }
@]@
where @(@
G : \B{R}^{Ny} \times \B{R}^{Na} \rightarrow \B{R}^{Ny}
@)@ is the
function representing this ODE; i.e.
@[@
G(y, a) = \left( \begin{array}{c}
- a_1 * y_0
\\
+ a_1 * y_0 - a_2 * y_1
\end{array} \right)
@]@
This @(@
G(y, a)
@)@ is linear with respect to @(@
y
@)@, hence
the implicit equation defining @(@
y^i
@)@ can be solved
inverting the a set of linear equations.
In the general case,
where @(@
G(y, a)
@)@ is non-linear with respect to @(@
y
@)@,
an iterative procedure is used to calculate @(@
y^i
@)@
from @(@
y^{i-1}
@)@.
Trapezoidal Time Grid
The discrete time grid, used for the trapezoidal approximation, is
denoted by @(@
\{ t_i \}
@)@ which is defined by:
@(@
t_0 = 0
@)@ and
for @(@
i = 1 , \ldots , Nz
@)@ and for @(@
j = 1 , \ldots , N(i)
@)@,
@[@
\begin{array}{rcl}
\Delta t_i & = & ( s_i - s_{i-1} ) / N(i)
\\
t_{S(i-1)+j} & = & s_{i-1} + \Delta t_i * j
\end{array}
@]@
where @(@
s_0 = 0
@)@,
@(@
N(i)
@)@ is the number of time grid points between
@(@
s_{i-1}
@)@ and @(@
s_i
@)@,
@(@
S(0) = 0
@)@, and
@(@
S(i) = N(1) + \ldots + N(i)
@)@.
Note that for @(@
i = 0 , \ldots , S(Nz)
@)@,
@(@
y^i
@)@ denotes our approximation for @(@
y( t_i , a )
@)@
and @(@
t_{S(i)}
@)@ is equal to @(@
s_i
@)@.
Black Box Method
A common approach to an inverse problem is to treat the forward problem
as a black box (that we do not look inside of or try to understand).
In this approach, for each value of the parameter vector @(@
a
@)@
one uses the
trapezoidal approximation
(on a finer grid that @(@
\{ s_i \}
@)@)
to solve for @(@
y_1 ( s_i , a )
@)@ for @(@
i = 1 , \ldots , Nz
@)@.
Two levels of Iteration
As noted above, the trapezoidal approximation
often requires an iterative procedure.
Thus, in this approach, there are two levels of iterations,
one with respect to the parameter values during the minimization
and the other for solving the trapezoidal approximation equation.
Derivatives
In addition, in the black box approach, differentiating the ODE solution
often involves differentiating an iterative procedure.
Direct application of AD to compute these derivatives
requires a huge amount of memory and calculations to differentiate the
iterative procedure.
(There are special techniques for applying AD to the solutions of iterative
procedures, but that is outside the scope of this presentation).
Simultaneous Method
The simultaneous forward and inverse method
uses constraints to include the solution of
the forward problem in the inverse problem.
To be specific for our example,
@[@
\begin{array}{rcl}
{\rm minimize}
& \sum_{i=1}^{Nz} H_i ( y^{N(i)} , a )
& \; {\rm w.r.t} \; y^1 \in \B{R}^{Ny} , \ldots , y^{S(Nz)} \in \B{R}^{Ny} ,
\; a \in \B{R}^{Na}
\\
{\rm subject \; to}
& y^j = y^{j-1} +
\left[ G( y^{j-1} , a ) + G( y^j , a ) \right] * \frac{ t_j - t_{j-1} }{ 2 }
& \; {\rm for} \; j = 1 , \ldots , S(Nz)
\\
& y^0 = F(a)
\end{array}
@]@
where for @(@
i = 1, \ldots , Nz
@)@,
@(@
N(i)
@)@ is the number of time intervals between
@(@
s_{i-1}
@)@ and @(@
s_i
@)@ (with @(@
s_0 = 0
@)@)
and @(@
S(i) = N(1) + \ldots + N(i)
@)@.
Note that, in this form, the iterations of the optimization procedure
also solve the forward problem equations.
In addition, the functions that need to be differentiated
do not involve an iterative procedure.
Source
The file ipopt_nlp_ode_problem.hpp
contains
source code that defines the example values and functions defined above.
Input File: cppad_ipopt/example/ode2.omh
