Prev Next

@(@\newcommand{\W}[1]{ \; #1 \; } \newcommand{\R}[1]{ {\rm #1} } \newcommand{\B}[1]{ {\bf #1} } \newcommand{\D}[2]{ \frac{\partial #1}{\partial #2} } \newcommand{\DD}[3]{ \frac{\partial^2 #1}{\partial #2 \partial #3} } \newcommand{\Dpow}[2]{ \frac{\partial^{#1}}{\partial {#2}^{#1}} } \newcommand{\dpow}[2]{ \frac{ {\rm d}^{#1}}{{\rm d}\, {#2}^{#1}} }@)@
ODE Fitting Using Fast Representation

Purpose
In this section we represent a more complex representation of the simultaneous forward and reverse ODE fitting problem (described above). The representation defines the problem using simpler functions that are faster to differentiate (either by hand coding or by using AD).

Objective Function
We use the following representation for the objective function : For @(@ k = 0 , \ldots , Nz - 1 @)@, we define the function @(@ r^k : \B{R}^{Ny+Na} \rightarrow \B{R} @)@ by @[@ \begin{array}{rcl} fg_0 (x) & = & \sum_{i=1}^{Nz} H_i ( y^{S(i)} , a ) \\ fg_0 (x) & = & \sum_{k=0}^{Nz-1} r^k ( u^{k,0} ) \end{array} @]@ where for @(@ k = 0 , \ldots , Nz-1 @)@, @(@ u^{k,0} \in \B{R}^{Ny + Na} @)@ is defined by @(@ u^{k,0} = ( y^{S(k+1)} , a ) @)@

Range Indices I(k,0)
For @(@ k = 0 , \ldots , Nz - 1 @)@, the range index in the vector @(@ fg (x) @)@ corresponding to @(@ r^k ( u^{k,0} ) @)@ is 0. Thus, the range indices are given by @(@ I(k,0) = \{ 0 \} @)@ for @(@ k = 0 , \ldots , Nz-1 @)@.

Domain Indices J(k,0)
For @(@ k = 0 , \ldots , Nz - 1 @)@, the components of the vector @(@ x @)@ corresponding to the vector @(@ u^{k,0} @)@ are @[@ \begin{array}{rcl} u^{k,0} & = & ( y^{S(k+1} , a ) \\ & = & ( x_{Ny * S(k+1)} \; , \; \ldots \; , \; x_{Ny * S(k+1) + Ny - 1} \; , \; x_{Ny * S(Nz) + Ny } \; , \; \ldots \; , \; x_{Ny * S(Nz) + Ny + Na - 1} ) \end{array} @]@ Thus, the domain indices are given by @[@ J(k,0) = \{ Ny * S(k+1) \; , \; \ldots \; , \; Ny * S(k+1) + Ny - 1 \; , \; Ny * S(Nz) + Ny \; , \; \ldots \; , \; Ny * S(Nz) + Ny + Na - 1 \} @]@

Initial Condition
We use the following representation for the initial condition constraint : For @(@ k = Nz @)@ we define the function @(@ r^k : \B{R}^{Ny} \times \B{R}^{Na + Ny} @)@ by @[@ \begin{array}{rcl} 0 & = & fg_i ( x ) = y_i^0 - F_i (a) \\ 0 & = & r_{i-1}^k ( u^{k,0} ) = y_i^0 - F_i(a) \end{array} @]@ where @(@ i = 1 , \ldots , Ny @)@ and where @(@ u^{k,0} \in \B{R}^{Ny + Na} @)@ is defined by @(@ u^{k,0} = ( y^0 , a) @)@.

Range Indices I(k,0)
For @(@ k = Nz @)@, the range index in the vector @(@ fg (x) @)@ corresponding to @(@ r^k ( u^{k,0} ) @)@ are @(@ I(k,0) = \{ 1 , \ldots , Ny \} @)@.

Domain Indices J(k,0)
For @(@ k = Nz @)@, the components of the vector @(@ x @)@ corresponding to the vector @(@ u^{k,0} @)@ are @[@ \begin{array}{rcl} u^{k,0} & = & ( y^0 , a) \\ & = & ( x_0 \; , \; \ldots \; , \; x_{Ny-1} \; , \; x_{Ny * S(Nz) + Ny } \; , \; \ldots \; , \; x_{Ny * S(Nz) + Ny + Na - 1} ) \end{array} @]@ Thus, the domain indices are given by @[@ J(k,0) = \{ 0 \; , \; \ldots \; , \; Ny - 1 \; , \; Ny * S(Nz) + Ny \; , \; \ldots \; , \; Ny * S(Nz) + Ny + Na - 1 \} @]@

Trapezoidal Approximation
We use the following representation for the trapezoidal approximation constraint : For @(@ k = 1 , \ldots , Nz @)@, we define the function @(@ r^{Nz+k} : \B{R}^{2*Ny+Na} \rightarrow \B{R}^{Ny} @)@ by @[@ r^{Nz+k} ( y , w , a ) = y - w - [ G( y , a ) + G( w , a ) ] * \frac{ \Delta t_k }{ 2 } @]@ For @(@ \ell = 0 , \ldots , N(k)-1 @)@, using the notation @(@ i = Ny * S(k-1) + \ell + 1 @)@, the corresponding trapezoidal approximation is represented by @[@ \begin{array}{rcl} 0 & = & fg_{Ny+i} (x) \\ & = & y^i - y^{i-1} - \left[ G( y^i , a ) + G( y^{i-1} , a ) \right] * \frac{\Delta t_k }{ 2 } \\ & = & r^{Nz+k} ( u^{Nz+k , \ell} ) \end{array} @]@ where @(@ u^{Nz+k,\ell} \in \B{R}^{2*Ny + Na} @)@ is defined by @(@ u^{Nz+k,\ell} = ( y^{i-1} , y^i , a) @)@.

Range Indices I(k,0)
For @(@ k = Nz + 1, \ldots , 2*Nz @)@, and @(@ \ell = 0 , \ldots , N(k)-1 @)@, the range index in the vector @(@ fg (x) @)@ corresponding to @(@ r^k ( u^{k,\ell} ) @)@ are @(@ I(k,\ell) = \{ Ny + i , \ldots , 2*Ny + i - 1 \} @)@ where @(@ i = Ny * S(k-1) + \ell + 1 @)@.

Domain Indices J(k,0)
For @(@ k = Nz + 1, \ldots , 2*Nz @)@, and @(@ \ell = 0 , \ldots , N(k)-1 @)@, define @(@ i = Ny * S(k-1) + \ell + 1 @)@. The components of the vector @(@ x @)@ corresponding to the vector @(@ u^{k,\ell} @)@ are (and the function @(@ fg (x) @)@ in cppad_ipopt_nlp ) @[@ \begin{array}{rcl} u^{k, \ell} & = & ( y^{i-1} , y^i , a ) \\ & = & ( x_{Ny * (i-1)} \; , \; \ldots \; , \; x_{Ny * (i+1) - 1} \; , \; x_{Ny * S(Nz) + Ny } \; , \; \ldots \; , \; x_{Ny * S(Nz) + Ny + Na - 1} ) \end{array} @]@ Thus, the domain indices are given by @[@ J(k,\ell) = \{ Ny * (i-1) \; , \; \ldots \; , \; Ny * (i+1) - 1 \; , \; Ny * S(Nz) + Ny \; , \; \ldots \; , \; Ny * S(Nz) + Ny + Na - 1 \} @]@

Source
The file ipopt_nlp_ode_fast.hpp contains source code for this representation of the objective and constraints.
Input File: cppad_ipopt/example/ode2.omh