key: cord-0781595-yyejvqa4
authors: Essoufi, E.-H.; Zafrar, A.
title: Boundary optimal control of time–space SIR model with nonlinear Robin boundary condition
date: 2021-11-06
journal: Int J Dyn Control
DOI: 10.1007/s40435-021-00886-1
sha: c0e0d6de56982068c639853df0a913a37c9e68f6
doc_id: 781595
cord_uid: yyejvqa4

A boundary optimal control problem arising in time-space SIR epidemic models is treated. In this work we aim with the control of the flux of infected individuals crossing part of boundary. On the other side of the domain, we suppose a nonlinear boundary condition of third kind: nonlinear Robin boundary condition, this condition models immersing individual crossing this part of the boundary of the domain of study. We give the existence and uniqueness of the solution of both state and optimal control problem ending some numerical tests throughout a simple example.

Mathematical modelling have been extensively employed for obtaining interesting quantitative and qualitative information in different fields especially in the healthcare sector. These informations are the basic criterion to provide useful guidelines to outbreak management and policy development. Epidemiology is one of essential population biology disciplines that is concerned with public health. Epidemiological models provide important powerful means for predicting the spread of a disease and then controlling and understanding the dynamics which effect its transmission and diffusion through a given population. A famous example of such model based on the compartmental analysis is the Susceptible-Infected-Recovered (SIR) model and its variants, for epidemic [1] [2] [3] , the Lotka-Volterra model for predator-prey dynamics [1, 4] and demographic and migration models established in sociology and demography [5, 6] .

The majority of models confronted in such literature are based on compartmental analysis and consisting of systems B A. Zafrar zafrar.abd@gmail.com E.-H. Essoufi e.h.essoufi@gmail.com 1 Laboratory MISI, Univ. Hassan 1, 26000 Settat, Morocco of ordinary differential equations (ODEs). In spite that these models are simple in their formulation, analysis, and simple to solve mathematically and numerically, they are inadequate to describe diffusion and dynamic in both space and time. A common way to deal with introduce spatial variation into such ODE models is by considering meta-population models (models describing a population consisting of n subpopulations which are connected by immigration or emigration) [7, 8] , in other words by defining regional compartments corresponding to different areas in physical space, with coupling terms added to the model equations to account for the movement of species among the different regions [8] [9] [10] [11] . This approach was lately developed in [12, 13] to model and describe the spread of the recent pandemic produced by COVID-19 virus among the different regions in Italy. Despite that this approach may be adequate for some application and description of complex spatial dynamics within compartments is difficult and possibly even non-feasible in this framework. Moreover, Some existing work consider an optimal control of SIR model by introducing density of vaccination in the model which is treated as control function [14] [15] [16] [17] and then minimize a given cost function. This procedure can be considered as introducing a new compartment of vaccination and look for a density that decrease the number of infected and susceptible individuals in considered population.

In this paper we consider a compartmental SIR-model based on partial differential equations (PDEs) which include spatial information more naturally [18] [19] [20] [21] [22] [23] . Generally, PDE models take into consideration a space continuous description of the suitable dynamics, enabling the description of dynamics in both time and space beyond all scales. This furnishes a significant interest over ODE models, whose capacity to describe spatial information is limited by the number of spatial compartments that the model includes. Our contribution consists in a boundary optimal control problem [24] [25] [26] [27] subject to an SIR-diffusion model describing a time-space spread of a disease. The constraint is a system of three PDEs with Dirichlet and Neumann boundary condition concerning the S and R unknown. Concerning the third unknown, we adopt a non linear Robin boundary condition (boundary condition of third kind), this is motivated by the fact that the Robin boundary conditions describe when part of the infected individuals I is "immersed" at the prescribed part of the boundary. This is the case where infected individuals can travel between different region and we interested in the influence of this hypothesis on the population of the specific region subject to our study. On the second part of the boundary, we control the flux of individual constituting the population.

The present paper is Briefly outlined as follows. In the next section we introduce the time splitting method in order to handle the non linearity and we state the existence and uniqueness of the solution. In Sect. 3 we study the optimal control problem and provide a necessary optimality condition. Section 4 is devoted to the numerical illustrations by mean of simple example.

We start this subsection with some useful notations. Let be an open bounded polyhedral domain subset of R d (d = 1, 2, 3) with a boundary = ∂ . We define the norm · L 2 ( ) for the Sobolev space L 2 ( ). We shall denote by Let the boundary consist of three disjoint parts = 1 ∪ 2 ∪ 3 , with 1 and 2 being the union of some (d − 1)-dimensional polyhedral domains. To consider a spatiallytemporally dependent problem, we define the spatial-temporal domains Q and i , for i = 1, 2, 3, respectively, by

The thermal conductivity q(x) ∈ C(¯ ) satisfies 0 < q 0 ≤ q ≤ q 1 < ∞ for some positive constants q 0 and q 1 . In the following sections, we shall frequently use the notation C to denote a generic constant, which may differ at every situation. When no confusion can occur, we shall ignore the symbols dx, ds and dt in the integrals for notational simplicity. Now we are in position to state our reaction-diffusion SIR model, which will be the direct problem with nonlinear Robin problem as follows:

In order to afford the existence and uniqueness, we suppose that the map g 2 (·, ·) satisfies the following properties

After all g(x, ·) is strictly increasing for a.e. x ∈ 1 , then it has an inverse which we express as g(x, ·). Let G, G :

It follows that for a.e. x ∈ 1 , G(x, ·) and G(x, ·) are continuous, strictly increasing, convex and G(x, 0) = G(x, 0) = 0. Definition 1.1 [28] The functions G(x, ·) and G(x, ·) are called complementary N -functions.

for all ξ ∈ R and a.e. x ∈ 1 . 2. The condition (G): We say that g satisfies the condition

for all ξ, η ∈ R and a.e. x ∈ 1 .

Let γ u be the trace of u on 1 and let g :

where we denote by M( 1 ) is the space of all measurable functions on 1 . Since G and G satisfy the ( 2 )-condition, it follows from [28] that L G ( 1 ) endowed with the Luxemburg norm

is a reflexive Banach space. In addition, by [28] one has the following generalized version of Hölder's inequality:

Now we introduce the state vector space V given by

The system of Eq. (1.1) are split into two systems of sub equations as follows, the nonlinear reaction equations

which are used for the first half of the time step, and the linear diffusion equations

which are used for the second half of the time step. The numerical method used to solve the equations is the forward Euler scheme.

The above equations transform to

⎦ , S n , I n and R n indicate the approximate values of S, I and R at position x and time nδt, n = 0, 1, . . . and where S n+ 1 2 , I n+ 1 2 , and R n+ 1 2 indicate the representative values at the half-time step.

Similarly, for the second half time step

Below our analysis will be reserved to the nonlinear elliptic boundary value problem given by (2.4) . For the sake of simplicity we shall omit the exponent n + 1 since the interval time is discretised we omit the dependence in time for the maps g, g, G and G.

The weak formulation of (2.4) is give as follows

Now, let us define the following nonlinear form

and we have the following result Lemma 2.1 Assume that g satisfies all assertions in (1.2) and that G and G satisfy the ( 2 )-condition. Then the bilinear form a(·, ·) is monotone and hemicontinuous and a(u, ·) ∈ V for all u ∈ V . Furthermore a(·, ·) is coercive,

is monotone. by the continuity argument of g(·) for a.e. x ∈ 1 , we obtain that for all u, v, w ∈ V ,

therefore a(·, ·) is hemicontinuous.

Using the fact that 0 ≤ G(x, g(ξ )) ≤ (1/c 1 ) G(x, ξ) , we get that a(u, ·) ∈ V for every u ∈ V . Finally, because of both G and G satisfy the ( 2 )-condition, it follows that

as consequence a(·, ·) is coercive.

The existence and uniqueness of the solution of the problem (2.4) is stated in the following Theorem 2.2 Assume that g satisfies (1.2) and that G and G satisfy the ( 2 ) -condition. Let f ∈ L p ( ) then the problem (2.4) has a unique weak solution. Moreover we have the following majoration

In the remainder part of this paper we restrict ourself to the N -functions where g(u) given by

Then the function g(·) obtained from G(·) have some interesting properties:

(1) For a nonnegative integer k, g(·) is k-differentiable;

(2) For u ∈ H 1 ( ), g(u) ∈ L 2 ( );

(3) For u ∈ H 1 ( ), g(|u|) ∈ L 1 ( ) and there exists a positive constant κ such that g(|u|) ≤ 1 + meas( )

where meas( ) is the measure of .

In fact (see [29] ), let u be given in H 1 ( ). Then u ∈ H 1/2 ( ).

Using embedding results for Orlicz-Sobolev spaces (see [28] ) since ⊂ R 2 , we have H 1/2 ( ) ⊂⊂ L G where the Nfunction G(u) = e u 2 − 1. Thus there exists a constant κ > 0 such that

Hence for a small (enough) > 0 and for r = u L G + we obtain 

making → 0 achieves the argument. Next our contribution is concerned in the formulation of optimal control problem. The classic way to do this is as follows: Given a measurement or observation data u ob of u, find an optimal control σ such that the desired profile is reached.

The regularized problems of optimal control problem then given as the following minimization problems:

subject to the parabolic equation (2.4), where C is a subset of L 2 ( ) defined by

ε > 0 is a regularisation parameter.

We have the following result due to [29] :

is a sequence such that u k → u a.e. on 1 , g satisfies (1.2) and

where C > 0 is a constant independent of k. Then

Proof By (1.2), g is continuous and u k → u a.e. on 1 , we deduce that g (u k ) → g(u) a.e. on 1 . Since g is strictly increasing map we have g(u)u ≥ 0 on 1 so that we may use Fatou's result to get

we conclude that |g(x)| ≤ g(x)x + S ∀x ∈ R Thus |g(u)| ≤ |g(u)||u| + S on 1 i.e., g(u) ∈ L 1 ( 1 ). Employing again the identity (2.12) we obtain, for each δ > 0 and for a.e. x ∈ 1 , either

where C δ = sup |x|≤δ −1 |g(x)|. For every measurable subset

for k greater than some given integer N > 0. Thus

where meas( s ) is the measure of s . Hence, the sequence of functions {g (u k )} has equi-absolutely continuous integrals. By Vitali's Convergence Theorem,

We are now prepared to prove the existence of an optimal solution.

Then there exists an element (û,σ ) ∈ C × K that minimizes (2.10) subject to (2.4) .

and

We deduce that {σ k } is bounded in L 2 ( 2 ) and u k H 1 is bounded. Hence we can extract a subsequence {(u k , σ k )} such that

Furthermore, trace theorems and compact embedding results imply u k → u in L 2 ( 1 ) ; this in turn implies u k → u pointwise a.e. on 1 after extracting subsequence. By setting v = u k in (2.13) we obtain

where M is a constant independent of k. By Lemma 2.3 ,

For each v ∈ C ∞ ( ), (2.14) allows us to pass to the limit in (2.13) to obtain

Then using the denseness of C ∞ ( ) in H 1 ( ) and the fact that g(û) ∈ L 2 ( 1 ), we obtain

Thus (û,σ ) ∈ C × K . Finally using the weak lower semicontinuity of J (·, ·), we conclude that (û,σ ) is indeed an optimal solution, i.e.

In this section we will derive a necessary condition that an optimal solution must satisfy. The existence of an optimal solution has been established in the above section. In Sect. 2 we have shown that for each σ ∈ L 2 ( 2 ) , there exists a unique u satisfying (2.5). Thus the state u is a well-defined function of σ and will be denoted by u = u(σ ).

Our aim is to state the necessary condition rigorously and express it in a more practical form. To this end we shall use the following theorem (see [30, Lemma 3.11, p. 1127 ] and [24] ) Theorem 3.1 Let B be a Banach space X and Y two reflexive Banach spaces. Let also be two C 1 functions

We suppose that, for all σ ∈ B, 1. There exists a unique u(σ ) such that F(σ, u(σ )) = 0, 2. ∂ 2 F(σ, u(σ )) is an isomorphism from X onto Y , Then, J (σ ) = L(σ, u(σ )) is differentiable and, for every

where p(σ ) ∈ Y is the adjoint state, unique solution to

In order to apply this result we shall need the following preparations. By Riesz's representation theorem we define

and

Theorem 3.2 The first order necessary optimality condition reads: any optimal control σ satisfies

where p is the unique solution of

Proof Let us put B = L 2 ( 2 ), X = V and Y = X . The first assertion is straightforward. The second assertion of Theorem 3.1 is fulfilled since the operator A in ∂ 2 F is coercive. Then, J (σ ) = L(σ, u(σ )) is differentiable and, for every ζ ∈ B,

where p(σ ) ∈ Y is the adjoint state, unique solution to

which is equivalent to

thus p, ∂ 2 F(σ, u(σ ))w = ∂ 2 L(σ, u(σ )), w .

and this is exactly (3.7). Now the necessary condition is given by

is then equivalent to β σ, ξ L 2 ( 2 ) − p, y(ξ ) L 2 ( 2 ) = 0, for all ξ which gives the condition (3.6).

We start by introducing the optimality system formed by state variational equation, adjoint variational equation and optimality condition (3.6) as follows:

A finite element discretization of the optimality system (4.1) is defined in the usual approach. We assume is a convex polygonal domain and we choose families of finite dimensional subspaces V h ⊂ H 1 ( ) satisfying the approximation property: there exists a constant C and an integer k such that

where u h is the finite approximation of u and h is the size of the approximation. In the rest of this paper we consider V h like a finite dimensional space with piecewise continuous functions v h ∈ C(¯ ) such that for every K ∈ T h , v h /K ∈ P 1 (K ), where T h is a family of element K forming and P 1 (K ) is the vectorial space of polynomial of degree less or equal 1. The family {V h } is a finite approximation of V defined above since if v h /K ∈ P 1 (K ) the trace T r(v h ) of v h is a scalar function which belongs to L G ( ). Now we provide a finite dimensional formulation of the optimality system (4.1) as follows: Find u h and p h such that

The state problem (4.3) is nonlinear, Newton Raphson numerical method is employed the compute the numerical state solution to this problem which is assembled in GNU Octave-6.1.0 using the usual piecewise finite element approximation. To compute the adjoint state p h in the second equation in the (4.3) LU decomposition is employed to inverse the resulting matrix system. In Figs. 1, 2 , 3, 4, 5 and 6 we visualize the solutions of the direct problem (1.1) for 

On the boundary part 2 we define the control σ ob as follows:

In Figs. 1, 2, 3, 4, 5 and 6 we illustrate the solution (S, I , R) with the above choice of parameters γ, β, β 1 and γ 1 at each time-step t 2 and t 50 , respectively. At each Fig. 9 we illustrate the given parameter σ ob compared with the computed one. Figure 10 shows the observed I ob and the computed I in the zone of the control and we remark that the computed control (optimal control) realises perfectly the desired profile.

On the boundary part 2 we define the control σ ob as follows:

In Figs. 11 and 12 we illustrate the computed I and the observed I for this example. In Fig. 13 we visualize the opti-Boundary optimal control of time-space SIR model with nonlinear Robin boundary condition 

In this paper we study a boundary optimal control problem subject to SIR model with nonlinear Robin boundary condition. The existence of the solution of the constraint problem Fig. 7 The computed I at the final time T corresponding to the example 1 Fig. 8 The observed I at final time T corresponding to the example 1 is obtained in Orlicz space for a specified class of functions defined the part of boundary where the Robin boundary condition is supposed. We provide an optimality condition and a numerical example illustration de solutions. Further-coming works will handle the error estimates corresponding the problem and a possible generalisation to a problem with general nonlinear Robin boundary condition, e.g. nonlinear boundary condition with general function g on the corresponding boundary. Fig. 9 The exact parameter σ ob (in orange x) and the computed optimal control (in blue *). Since we impose the Dirichlet boundary condition in the model we observe that the optimal control vanishes at boundaries {0} and {1}   Fig. 10 The exact and the computed I on the boundary of control. We remark that the desired profile is perfectly reached Fig. 13 The exact parameter σ ob (in orange x) and the computed optimal control (in blue *) on the boundary 2 . Since we impose the Dirichlet boundary condition in the model we observe that the optimal control vanishes at boundaries {0} and {1} 

Mathematical models in population biology and epidemiology

The mathematics of infectious diseases

A contribution to the mathematical theory of epidemics

Mathematical biology: I-an introduction

Differential equations: a modeling approach

Applied mathematical demography

Optimal control of vaccine distribution in a rabies metapopulation model

Spatiotemporal dynamics of epidemics: synchrony in metapopulation models

Networks and epidemic models

Compartmental models of ecological and environmental systems

On the analysis of a multiregions discrete sir epidemic model: an optimal control approach

Spread and dynamics of the covid-19 epidemic in Italy: effects of emergency containment measures

Modelling the covid-19 epidemic and implementation of population-wide interventions in Italy

Optimal control of epidemics with limited resources

An optimal control problem for a spatiotemporal sir model

Optimal control of vaccination and treatment for an sir epidemiological model

Optimal control for a sir epidemiological model with time-varying populations

The effects of spatial heterogeneity in population dynamics

Spatial ecology via reactiondiffusion equations

Partial differential equations in ecology: spatial interactions and population dynamics

Numerical simulation of a susceptible-exposed-infectious space-continuous model for the spread of rabies in raccoons across a realistic landscape

Spatial models and biomedical applications

Qualitative analysis of a Lotka-Volterra predatorprey system with migration

Contrôle optimal de systemes gouvernés par des équations aux dérivées partielles

Mixed finite element method for Dirichlet boundary control problem governed by elliptic PDES

Finite element error estimates on the boundary with application to optimal control

Computational methods for boundary optimal control and identification problems

Sobolev spaces

Analysis and finite element approximation of an optimal control problem in electrochemistry with current density controls

Optimal control of an elastic contact problem involving tresca friction law

Acknowledgements This is supported by CNRST Morocco. The authors of this work would like to express their gratitude to Editors and anonymous Referees for their prolific remarks and valuable suggestions to improve the quality of the paper.

Funding There was no funding.

Code Availability Not applicable.

Conflict of interest I declare that there is no conflict of interest in the publication of this article, and that there is no conflict of interest with any other author or institution for the publication of this article.Ethical statements I hereby declare that this manuscript is the result of my independent creation under the reviewers' comments. Except for the quoted contents, this manuscript does not contain any research achievements that have been published or written by other individuals or groups. I am the only author of this manuscript. The legal responsibility of this statement shall be borne by me.