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.