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