Prev Next Index-> contents reference index search external Up-> CppAD Appendix deprecated cppad_ipopt_nlp ipopt_nlp_ode ipopt_nlp_ode_problem deprecated-> include_deprecated FunDeprecated CompareChange omp_max_thread TrackNewDel omp_alloc memory_leak epsilon test_vector cppad_ipopt_nlp old_atomic zdouble autotools cppad_ipopt_nlp-> ipopt_nlp_get_started.cpp ipopt_nlp_ode ipopt_ode_speed.cpp ipopt_nlp_ode-> ipopt_nlp_ode_problem ipopt_nlp_ode_simple ipopt_nlp_ode_fast ipopt_nlp_ode_run.hpp ipopt_nlp_ode_check.cpp ipopt_nlp_ode_problem-> ipopt_nlp_ode_problem.hpp Headings-> Notation Forward Problem Measurements ---..Simulation Analytic Solution ---..Simulation Parameter Values ---..Simulated Measurement Values Inverse Problem Trapezoidal Approximation ---..Trapezoidal Time Grid Black Box Method ---..Two levels of Iteration ---..Derivatives Simultaneous Method Source

An ODE Inverse Problem Example

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}$$

Simulation Parameter Values
 $\bar{a}_0 = 1$   initial value of $y_0 (t, a)$ $\bar{a}_1 = 2$   transfer rate from compartment zero to compartment one $\bar{a}_2 = 1$   transfer rate from compartment one to outside world $\sigma = 0$   standard deviation of measurement noise $e_i = 0$   simulated measurement noise, $i = 1 , \ldots , Nz$ $s_i = i * .5$   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.