Download Code

In this tutorial, we propose the Hum method to approximate numerically the control in a null controllability problem for a non linear population dynamics model structuring in age and spatial diffusion.

For given the positive function $F$, we consider the following Population Dynamics Model Structuring by age and spatial diffusion:

$$(1):\left\{ \begin{array}{ccc} \partial_{t} y+\partial_{a}y-Δ y+\mu y&=\chi_{\Theta}V&\text{ in }\Omega\times(0,A)\times (0,T) ,\\ \dfrac{\partial y}{\partial \nu}&=0& \partial \Omega \times (0,A)\times (0,T)\\ y\left( x,0,t\right) &=F\left(\int_{0}^{A}\beta(a)yda\right) &\text{ in } \Omega\times(0,T) \\ y\left(x,a,0\right)&=y_{0}& \text{ in } \Omega\times(0,A). \end{array} \right.$$

where, $\Omega$ is a boundary subset of $\mathbb{R}^{N}$, $A$ is the maximal age and $T>0$ the final time. Also, $\Theta=\omega× (a_1,a_2)× (0,T),$ is the control support. Here, $\beta$ is the fertility rate depending also the age a, $\mu$ the mortality rate depending of the age a, $y_0$ the initial condition in $L^2(\Omega× (0,A)),$ $V$ is the control term and $y$ represent the density of the population of age $a$ at time $t$ and location $x∈ \Omega.$

We assume that the fertility rate $\beta$ and the mortality rate $\mu$ satisfy the demographic property:

$$(H_1)= \left\{ \begin{array}{c} \mu(a)\geq 0 \text{ for every } a\in (0,A)\\ \mu\in L^{1}_{loc}(0,A)\quad \int_{0}^{A}\mu(a)da=+\infty \end{array} \right.,$$

and

$$(H_2)=\left\{ \begin{array}{ccc} \beta(a)\geq 0\quad a.e \quad (0,A)\\ \beta(a)\in C^{1}([0,A])\\ \beta(a)=0 \text{ } \forall a\in (0,\hat{a}) \text{ where } \hat{a}\in (0,A), \\ \end{array} \right.$$

where

$$0≤ a_1<a_2≤ A \text{ and } a_1≤ \hat{a}.$$

The function $F,$ verify:

$$(H_3)= \left\{ \begin{array}{c} F \quad \text{ positive, globally Lipchitz }\\ F(0)=0,\quad F \text{ continuous and differentiable in } 0. \end{array} \right..$$

Numerical Simulations

In this part, the objective is to illustrate numerically the approximate null controllability of the nonlinear problem.

Indeed, we consider the following system:

$$(4):\left\{ \begin{array}{ccc} \partial_{t} y+\partial_{a}y-\partial^{2}_{xx}y+\mu y&=\chi_{\Theta}V&\text{ in }(0,1)\times(0,A)\times (0,T) ,\\ y'(0)=y'(1)&=0& \text{ in } (0,A)\times (0,T)\\ y\left( x,0,t\right) &=F\left(\int_{0}^{A}\beta(a)yda\right) &\text{ in } (0,1)\times(0,T) \\ y\left(x,a,0\right)&=y_{0}& \text{ in } (0,1)\times(0,A). \end{array} \right.$$

Discretization and simulation of uncontrolled system

The idea in this part is to highlight the numerical simimulation of the nonlinear problem. The first parts is to reduce the PDE to the finite dimensional system of the form

$$\dot{U_l}=A_lU_l+vectorF+BV,$$

where $A_l$ and $B$ are matrices, and

$$\begin{pmatrix} U_l(t) \\ V(t) \\ \end{pmatrix}$$

is the finite dimensional state and control vector. Here, VectorF is the contribution of the nonlinear part, which comes from births. So let’s consider the following system

$$(5):\left\{ \begin{array}{ccc} \dfrac{\partial y}{\partial t}+\textbf{A}y+\mu y&=\chi_{\Theta}V&\text{ in }\Omega\times(0,A)\times (0,T) ,\\ \dfrac{\partial y}{\partial \nu}&=0& \partial \Omega \times (0,A)\times (0,T)\\ y\left( x,0,t\right) &=F\left(\int_{0}^{A}\beta(a)yda\right) &\text{ in } \Omega\times(0,T) \\ y\left(x,a,0\right)&=y_{0}& \text{ in } \Omega\times(0,A). \end{array} \right.$$

where

$$\textbf{A}y=\partial_{a}y-\partial^{2}_{xx}y.$$

To solve $(4)$, the space (dimension 1) and age discretization is performed with finite difference method on rectangular grid on $[0,L]\times [0,A]$. For a given rectangular grid $\mathcal{T}$ with vertex $(x_i,a_j)\quad 1\leq i\leq N,\quad 1\leq j\leq M$ and uniform step size (without loss of generality) $\Delta x$ in $x-$direction and $\Delta a$ in $a-$direction, we denote by the diameter of the grid The finite difference approximation of the diffusion term of the operator $\textbf{A}$ is given by

$$\dfrac{\partial^2 y}{\partial x^2}(x_i,a_j,t)=\dfrac{y(x_{i+1},a_j,t)-2y(x_i,a_j,t)+y(x_{i-1},a_j,t)}{(\Delta x)^2}.$$

The finite difference approximation of the aging term of the operator A is given by

$$\dfrac{\partial y}{\partial a}(x_i,a_j,t)=\dfrac{y(x_{i},a_{j},t)-y(x_i,a_{j-1},t)}{\Delta a}.$$

Let $y_{i,j}(t)$ be the approximation of $y(x_i,a_j)$ and

$$U_{l}(t) = (y_{i,j}(t))_{1\leq i\leq N,\quad 1\leq j\leq M}$$

where $y_{i,j}$ is at the position $i+j*(N-1)$ and $A_l$ the matrix of the mortality approximation, the diffusion approximation and the aging approximation.

Take into account the newborns

We denote by $\beta(a),$ the fertility rate, the newborn is given by:

$$y(x,0,t)=F\left(\int_{0}^{A}\beta(a)y(x,a,t)da\right).$$

We approximate $\int_{0}^{A}\beta(a)y(x,a,t)da$ by

$$\int_{0}^{A}\beta(a)y(x,a,t)da=(A/n)\sum_{j=1}^{n-1}\beta(a_j)y(x,a_j,t)+\dfrac{\beta(0)y(x,0,t)+\beta(A)y(x,A,t)}{2}.$$

Then

$$y_{i,0}=A/n\sum_{k=1}^{n-1}\beta(a_k)y(x_i,a_k,t)+\dfrac{\beta(0)y(x_i,0,t)+\beta(A)y(x_i,A,t)}{2}\text{ for } i=1:n.$$

But as

$$\dfrac{\partial y_{i,1}(t)}{\partial a}=\dfrac{y_{i,1}(t)-y_{i,0}(t)}{\Delta a}\text{ } i=1:n,$$

then

$$\dfrac{\partial y_{i,1}(t)}{\partial a}=\dfrac{-A/n\sum_{k=1}^{n-1}\beta(a_k)y_{i,k}(t)+y_{i,1}(t)+\dfrac{\beta(0)y(x_i,0,t)+\beta(A)y(x_i,A,t)}{2}}{\Delta a}\quad i=1:n,$$

We create also the nonlinear vectorF from the births. Indeed, if we denote by

$$(N_l)_{1\leq i\leq n*n;\leq j\leq n*n}$$

the matrix of the births,

$$P_l=N_lU_l\quad where\quad \left(\dfrac{A\sum_{k=1}^{n-1}\beta(a_k)y_{i,k}(t)+\dfrac{\beta(0)y(x_i,0,t)+\beta(A)y(x_i,A,t)}{2}}{n\Delta a}\right)_{1\leq i\leq n}=\left(P_l\right)_{1\leq i\leq n}\text{ } and \text{ } \left(P_l\right)_{n+1\leq i\leq n*n}=0$$

and we define the vectorF by:

$$vectorF=F(N_lU_l).$$

The corresponding discrete system is given for a continuous initial solution $y_0$ by

$$\dfrac{dU_l}{dt}=A_lU_l+vectorF\quad y(0)=y_0(x,a)$$

Example 1

Simulation of the uncontrolled system:

For the simulation, we take $L=1,\text{ } A=10\text{, }\Delta x=1/54\quad \Delta a= 5/54$ and $T=40$. Moreover $F(t)=t\Phi(t)$ a Globally Lipchitz function where

$$\Phi(t)=0.75\exp(-t)$$

The fertility $\beta$ is given by:

$$\beta(a) = \begin{cases} 0 & \text { if } a<A/18 \\ \frac{\alpha(a-5/18)^{\gamma-1}e^{-(a-5/18)/v}}{v^{\gamma}\Gamma(\gamma)} \text{ if } A/18\leq a\leq 17A/18\\ 0 &\text{ if } a>17A/18. \end{cases}$$

here $v=1$, $\gamma=5$, $\alpha=7$ and the mortality rate $\mu$ by

$$\mu(a)=\dfrac{1}{10(A-a)}$$

---

Video 1. Evolution of the state of the uncontrolled system between t=0 and t=30

We observe in this video the growth (aging) of news borns over time. We also notice the high birth rate of this population.

Construction of the control and numerical simulation

We construct the control problem, which consists in minimizing the functional and we choose the classical Hum functional and the control matrix $B=\chi_{\Theta}$ where $\Theta=(x_1,x_2)\times (a_1,a_2)\times (0,T),$ and we suppose that $U_T$ is the desired state:

$$J_{\epsilon}(U_l,V_l)=\dfrac{1}{\epsilon}\int_{0}^{A}\int_{0}^{L}\left(U_l(x,a,T)-U_T\right)^2dxda+\int_{0}^{T}\int_{a_1}^{a_2}\int_{x_1}^{x_2}V_{l}^{2}dxdadt.$$

Here $U_T=0.$\ The approximate null controllability become the minimization of the functional $J,$ where $(U_l,V_l)$ verify the following system

$$\dfrac{dU_l}{dt}=A_lU_l+vectorF+BV_l\quad U_l(0)=(y_0(x_i,a_j))_{1\leq i\leq N,1\leq j\leq M}$$

In this part the fertility $\beta$ is given by:

$$\beta(a) = \begin{cases} 0 & \text { if } a<A/5 \\ \frac{\alpha(a-/8)^{\gamma-1}e^{-(a-/8)/v}}{v^{\gamma}\Gamma(\gamma)} \text{ if } A/5\leq a\leq 4A/5\\ 0 &\text{ if } a> 4A/5, \end{cases}$$

here $v=1$, $\gamma=5$, $\alpha=7$ and

$$\mu(a)=\dfrac{1}{15(A-a)}$$

and $VectorF=(P_lU_l).^2$.This means that we consider that $F(t)=t^2$ which is globally lipchizt on a compact which is our case.

The system becomes:

$$\dfrac{dU_l}{dt}=A_lU_l+(P_lU_l).^2+BV_l\quad U_l(0)=(y_0(x_i,a_j))_{1\leq i\leq N,1\leq j\leq M}.$$

Example 2

In this example, we take $ε=0.05$, $\Delta x=0.05$ and $\Delta a= 0.5$, with the initial condition

$$y_0(x,a)=e^{-(10(a-3/10)^2+70(x-4/10))}$$

the following numerical results were obtained for $T=12$ and $\Delta t=120$.

Evolution of the control V between t=0 and t=12.

Video 2. Evolution of the state of controlled system between t=0 and t=12

We notice that, with a control domain that takes into account the most populated part of the domain, the Hum method applied to the Casadi algorithm gives us an interesting approximative controllability result while maintaining the positivity of the density. Evolution of the controlled state (density of the population) between t=0 and t=12.

Video 3. Evolution of the control function between t=0 and t=12

Example 3

Here, we keep the same parameters of the Example 1 but we change the final time (here T=20) , and we take $\Delta x=1/18$, $Δ a=5/18$, $\epsilon=0.05$ the following numerical results were obtained:

Video 4. Evolution of the state of the controlled system between t=0 and t=20 Video 4. Evolution of the state of the controlled system between t=0 and t=20

Example 4

In this part we consider a smaller control domain. Indeed we take $\Theta =(0,2/5)\times(0,4)\times(0,T)$. In addition we also remove $\epsilon=0.001$ and $T=20.$ The rest is identical to the data of example 2.

We obtain the following numerical results:

Evolution of the control V between t=0 and t=20.

Video 4. Evolution of the state of the controlled system between t=0 and t=20

Here although the result is acceptable, we notice that the state of the controlled system is negative for most of the time.\ Evolution of the controlled state (the density of the population) between t=0 and t=20.

Video 5. Evolution of the control function between t=0 and t=20

We notice here that the solution to the final time tends to 0, but during the time of control much of the state of the system was negative. We notice that numerically we have approximation and the positivity of the state, if the support of the control covers the most populated part of the domain (space, age).

However it should be noted that we could not take a positivity constraint in our simulations. It should be noted that we have better approximations when we do not consider the cost of control.

References

[1] Dycon Toolbox

[2] CasaDi toolbox

[3] Luca Gerardo-Giorda: NUMERICAL APPROXIMATION OF DENSITY DEPENDENT DIFFUSION IN AGE-STRUCTURED POPULATION DYNAMICS

---