Jump to content

Van der Pol OED: Difference between revisions

From mintOC
 
(26 intermediate revisions by the same user not shown)
Line 5: Line 5:
}}
}}


The '''Van der Pol problem''' is a variation of the [[:Van der Pol]] problem. It looks for optimal time intervals to measure the two states in order to minimize the uncertainty of a follow-up parameter estimation problem for the two unknown parameters.
The '''Van der Pol OED problem''' is a variation of the [[:Van der Pol Oscillator]] problem. It looks for optimal time intervals to measure the two states in order to minimize the uncertainty of a follow-up parameter estimation problem for the two unknown parameters.


The mathematical equations form a small-scale [[:Category:ODE model|ODE model]]. It also includes state sensitivities, the Fisher information matrix entries and integrated sampling states.
The mathematical equations form a small-scale [[:Category:ODE model|ODE model]]. It also includes state sensitivities, the Fisher information matrix entries and integrated sampling states.
Line 15: Line 15:
<math>
<math>
\begin{array}{rcl}
\begin{array}{rcl}
\dot{x_1}(t) &=& p_1 \cdot (1 - x_2(t)^2) \cdot x_1(t) - x_2(t) + u(t), && t \in [0,t_f], \quad x_1(0) = 0, \\
\dot{x_1}(t) &=& (p_1 - x_2(t)^2) \cdot x_1(t) - x_2(t) + u(t), && t \in [0,t_f], \quad x_1(0) = 0, \\
\dot{x_2}(t) &=& p_2 \cdot x_1(t), && t \in [0,t_f], \quad x_2(0) = 1.
\dot{x_2}(t) &=& p_2 + x_1(t), && t \in [0,t_f], \quad x_2(0) = 1.
\end{array}  
\end{array}  
</math>
</math>
Line 22: Line 22:
Additionally, we add the constraint
Additionally, we add the constraint
<math>
<math>
  x_1(t) \geq -0.25 \; t\in [0,t_f]
  x_1(t) \geq -0.5, \; t\in [0,t_f].
</math>
</math>


The initial values and <math>t_f = 10</math> are fixed. We are interested in how to fish and when to measure, with an upper bound <math>M</math> on the measuring time. We can measure the states directly, <math>h^1(x(t)) = x_1(t)</math> and <math>h^2(x(t)) = x_2(t)</math>. We use two different sampling functions, <math>w^1(\cdot)</math> and <math>w^2(\cdot)</math> in the same experimental setting. This can be seen either as a two-dimensional measurement function <math>h(x(t))</math>, or as a special case of a multiple experiment, in which <math>u(\cdot), x(\cdot)</math>, and <math>G(\cdot)</math> are identical. The experimental design problem then reads
The initial values and <math>t_f = 10</math> are fixed. We are interested in how to choose the control <math>u</math> and when to measure, with an upper bound <math>M</math> on the measuring time. We can measure the states directly, <math>h^1(x(t)) = x_1(t)</math> and <math>h^2(x(t)) = x_2(t)</math>. We use two different sampling functions, <math>w^1(\cdot)</math> and <math>w^2(\cdot)</math> in the same experimental setting. This can be seen either as a two-dimensional measurement function <math>h(x(t))</math>, or as a special case of a multiple experiment, in which <math>u(\cdot), x(\cdot)</math>, and <math>G(\cdot)</math> are identical.


Now we formulate the OED problem with <math>\theta := (p_1, p_2)</math>:
<p>
<p>
<math>
<math>
  \begin{array}{lll}
  \begin{array}{lll}
  \displaystyle \min_{x,G,F,z^1,z^2,u,w^1,w^2} && \text{trace} \; \left( F^{-1}(t_f) \right) \\
  \displaystyle \min_{x,G,F,z,w,u} && \text{trace} \; \left( F^{-1}(t_f) \right) \\
  \text{subject to} \\
  \text{subject to} \\
\quad \dot{x_1}(t) & = & p_1 \cdot (1 - x_2(t)^2) \cdot x_1(t) - x_2(t) + u(t),\\
\quad \dot{x}(t) & = & f(x(t),u(t),\theta) \\
\quad \dot{x_2}(t) & = &  p_2 \cdot x_1(t),\\
\quad \dot{G}(t) & = & f_x(x(t),u(t),\theta) G(t) + f_\theta(x(t),u(t),\theta) \\
\quad \dot{G_{11}}(t) & = & f_{x11}(\cdot) \; G_{11}(t) + f_{x12}(\cdot) \; G_{21}(t) + f_{p12}(\cdot), \\
\quad \dot{F}(t) & = & \sum_{i=1}^{n_o} w_i(t)(h^i_x(x(t))G(t))^T(h^i_x(x(t))G(t)) \\
\quad \dot{G_{12}}(t) & = & f_{x11}(\cdot) \; G_{12}(t) + f_{x12}(\cdot) \; G_{22}(t), \\
\quad \dot{z}(t) & = & w(t), \\
\quad \dot{G_{21}}(t) & = & f_{x21}(\cdot) \; G_{11}(t) + f_{x22}(\cdot) \; G_{21}(t), \\
\quad x(0) & = & x_0 \\
\quad \dot{G_{22}}(t) & = & f_{x21}(\cdot) \; G_{12}(t) + f_{x22}(\cdot) \; G_{22}(t) + f_{p22}(\cdot), \\
\quad G(0) & = & \frac{\partial x(0)}{\partial \theta} \\
\quad \dot{F_{11}}(t) & = & w^1(t) G_{11}(t)^2 + w^2(t) G_{21}(t)^2, \\
\quad F(0) & = & I \cdot \varepsilon_{\mathrm{reg}}, \\  
\quad \dot{F_{12}}(t) & = & w^1(t) G_{11}(t) G_{12}(t) + w^2(t) G_{21}(t) G_{22}(t), \\
\quad z(0) & = & 0 \\
\quad \dot{F_{22}}(t) & = & w^1(t) G_{12}(t)^2 + w^2(t) G_{22}(t)^2, \\
\quad x_1(t) & \in & \mathcal{X} \\
\quad \dot{z^1}(t) & = & w^1(t), \\
\quad u(t) & \in & \mathcal{U} \\
\quad \dot{z^2}(t) & = & w^2(t), \\[1.5ex]
\quad w(t) & \in & \mathcal{W} \\
\quad x(0) &=& (0, 1), \\
\quad z_i(t_f) & \leq & M_i
\quad G(0) &=& F(0) = 0, \\
\quad z^1(0) &=& z^2(0) = 0, \\[1.5ex]
\quad u(t) & \in & \mathcal{U}, \; w^1(t) \in \mathcal{W}, \; w^2(t) \in \mathcal{W}, \\
\quad 0    & \le & M - z(t_f)
   \end{array}
   \end{array}
</math>
</math>
</p>
</p>


with
The evolution of the symmetric matrix <math>F: \left[0,t_f \right] \rightarrow \mathbb{R}^{2 \times 2}</math> is given by the weighted sum of observability Gramians
<p>
<math>h^i_x (x(t)) G(t), \ i = 1,2</math> for each observed function of states.
  f_{x11}(\cdot) &= \partial f_1(\cdot) / \partial x_1 = p_1 \cdot (1-x_2(t)^2), \\
  f_{x12}(\cdot) &= -2 \cdot p_1 \cdot x_2(t) \cdot x_1(t) - 1, \\
  f_{x21}(\cdot) &= p_2, \\
  f_{x22}(\cdot) &= 0,\\
  f_{p12}(\cdot) &= \partial f_1(\cdot) / \partial p_1 = (1-x_2(t)^2) \cdot x_1(t)</math>, \\
  f_{p24}(\cdot) &= \partial f_2(\cdot) / \partial p_2 = x_1(t)</math>. 
</p>
Note that the state <math>F_{21}(\cdot) = F_{12}(\cdot)</math> has been left out for reasons of symmetry.


== Parameters ==
== Parameters ==
These fixed values are used within the model:


We use <math>t_f=12</math>, <math>p_1 = p_2 = p_3 = p_4 = 1</math>, and <math>p_5 = 0.4</math>, <math>p_6 = 0.2</math>. The upper bound on the measurement time intervals is chosen as <math>M=4</math>.
{| border="1" align="center" cellpadding="5" cellspacing="0"
|- bgcolor=#c7c7c7
! Symbol !! Value !! Description
|-
| align=center | <math>p_1</math> || align=right | 1  || Unknown parameter
|-
| align=center | <math>p_2</math> || align=right | 1 || Unknown parameter
|-
| align=center | <math>t_\mathrm{f}</math> || align=right | 10 || Horizon of the control problem
|-
| align=center | <math>\varepsilon_{\mathrm{reg}}</math> || align=right | 0.1 || Regularization of Fisher matrix
|-
| align=center | <math>\mathcal{X}</math> || align=right | [-0.5,<math>\infty</math>] || Bounds of <math>x_1</math>
|-
| align=center | <math>\mathcal{U}</math> || align=right | [-1,1] || Bounds of control function
|-
| align=center | <math>\mathcal{W}</math> || align=right | [0,1] || Bounds of measurement function
|-
| align=center | <math>M_1, M_2</math> || align=right | 2 || Maximum measurement time
|}


== Reference Solutions ==
== Reference Solutions ==
Line 71: Line 80:
Here is one local solution to the above control problem.
Here is one local solution to the above control problem.


<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2">
<gallery caption="Reference solution plots" widths="500px" heights="300px" perrow="1">
  Image:vplotkaUStates.png| States and fishing control.
  Image:Van_der_Pol_OED.png| States, control, and sampling functions for a local optimum.
Image:vplotkaUSensitivities.png| Sensitivities G().
Image:vplotkaUMeas1.png| Sampling function for first state.
Image:vplotkaUMeas2.png| Sampling function for second state.
</gallery>
</gallery>
== Source Code ==
Model descriptions are available in
* [[:Category:AMPL | AMPL]] at [[Lotka Experimental Design (AMPL)]]
* [[:Category: VPLAN | VPLAN code]] at [[Lotka Experimental Design (VPLAN)]]
== Variants ==
There are several alternative formulations and variants of the above problem, in particular
* a prescribed time grid for the control function <bib id="Sager2006" />, see also [[Lotka Experimental Design (AMPL)]],
* no fishing, i.e., <math>u \equiv 0</math>,
* different fishing control functions for the two species,
* different parameters and start values.
== Miscellaneous and Further Reading ==
The Lotka Volterra fishing problem was introduced by Sebastian Sager in a proceedings paper <bib id="Sager2006" /> and revisited in his PhD thesis <bib id="Sager2005" />. These are also the references to look for more details. The experimental design problem has been described in the habilitation thesis of Sager, <bib id="Sager2011d" />.


== References ==
== References ==
Line 105: Line 92:
[[Category:ODE model]]
[[Category:ODE model]]
[[Category:Bang bang]]
[[Category:Bang bang]]
[[Category:Population dynamics]]

Latest revision as of 09:25, 26 March 2026

Van der Pol OED
State dimension: 1
Differential states: 11
Discrete control functions: 3


The Van der Pol OED problem is a variation of the Van der Pol Oscillator problem. It looks for optimal time intervals to measure the two states in order to minimize the uncertainty of a follow-up parameter estimation problem for the two unknown parameters.

The mathematical equations form a small-scale ODE model. It also includes state sensitivities, the Fisher information matrix entries and integrated sampling states.

Mathematical formulation

We are interested in estimating the parameters p1 and p2 of the initial value problem

x1˙(t)=(p1x2(t)2)x1(t)x2(t)+u(t),t[0,tf],x1(0)=0,x2˙(t)=p2+x1(t),t[0,tf],x2(0)=1.

Additionally, we add the constraint x1(t)0.5,t[0,tf].

The initial values and tf=10 are fixed. We are interested in how to choose the control u and when to measure, with an upper bound M on the measuring time. We can measure the states directly, h1(x(t))=x1(t) and h2(x(t))=x2(t). We use two different sampling functions, w1() and w2() in the same experimental setting. This can be seen either as a two-dimensional measurement function h(x(t)), or as a special case of a multiple experiment, in which u(),x(), and G() are identical.

Now we formulate the OED problem with θ:=(p1,p2):

minx,G,F,z,w,utrace(F1(tf))subject tox˙(t)=f(x(t),u(t),θ)G˙(t)=fx(x(t),u(t),θ)G(t)+fθ(x(t),u(t),θ)F˙(t)=i=1nowi(t)(hxi(x(t))G(t))T(hxi(x(t))G(t))z˙(t)=w(t),x(0)=x0G(0)=x(0)θF(0)=Iεreg,z(0)=0x1(t)𝒳u(t)𝒰w(t)𝒲zi(tf)Mi

The evolution of the symmetric matrix F:[0,tf]2×2 is given by the weighted sum of observability Gramians hxi(x(t))G(t), i=1,2 for each observed function of states.

Parameters

These fixed values are used within the model:

Symbol Value Description
p1 1 Unknown parameter
p2 1 Unknown parameter
tf 10 Horizon of the control problem
εreg 0.1 Regularization of Fisher matrix
𝒳 [-0.5,] Bounds of x1
𝒰 [-1,1] Bounds of control function
𝒲 [0,1] Bounds of measurement function
M1,M2 2 Maximum measurement time

Reference Solutions

Here is one local solution to the above control problem.

References

There were no citations found in the article.