Optimal Control and Computational Method for the Resolution of Isoperimetric Problem in a Discrete-Time SIRS System

We consider a discrete-time susceptible-infected-removed-susceptible “again” (SIRS) epidemic model, and we introduce an optimal control function to seek the best control policy for preventing the spread of an infection to the susceptible population. In addition, we define a new compartment, which models the dynamics of the number of controlled individuals and who are supposed not to be able to reach a long-term immunity due to the limited effect of control. Furthermore, we treat the resolution of this optimal control problem when there is a restriction on the number of susceptible people who have been controlled along the time of the control strategy. Further, we provide sufficient and necessary conditions for the existence of the sought optimal control, whose characterization is also given in accordance with an isoperimetric constraint. Finally, we present the numerical results obtained, using a computational method, which combines the secant method with discrete progressive-regressive schemes for the resolution of the discrete two-point boundary value problem.


Introduction
Many mathematical models in epidemiology are used to assist in finding the most appropriate control strategies for a given group of individuals who belong to different classes.These classes are often represented in epidemic systems, using compartments that are usually named susceptible (S), exposed (E), infectious or infected (I) and removed or recovered (R) [1].In this paper, we are interested in the study of a population infected by an epidemic and whose dynamics are described using a discrete-time SIRS system.The SIRS models in the continuous-time case have been widely studied by many researchers as in [2], where Acedo et al. proposed an analytical approach to find the exact global solution of the classical SIRS epidemic system.Furthermore, there are Alexander and Moghadas in [3] and Hu et al. in [4], who all provided bifurcation analysis of the SIRS model with different incidence rates.The authors who contributed with Teng in [5] and in [6] found significant results from the study of the persistence and extinction of disease using SIRS models.As for Jin et al. in [7], Liu and Zhou in [8] and Chen in [9], they obtained stability conditions for other SIRS systems.A stability analysis of the SIRS model in the discrete-time case is not often available, but there exist interesting analyses done for some classes of this type of model; see for example, the work of Hu et al. in [10].As an application of such models in a particular case of disease, Mukhopadhyay and Tapaswi published their paper about Japanese encephalitis in [11].Other authors studied SIRS dynamics when the model framework was in the form of a discrete metapopulation-like system [12].
In parallel, there are many researchers who have benefited from modeling approaches in epidemics, in order to determine the best prevention strategies against the spread of infection to susceptibles, using different optimization techniques such as optimal control methods; see examples in [13][14][15][16][17][18][19][20].On the other hand, some models as in [21][22][23] discussed the impact of limited public health resources in the propagation of infectious diseases, but there are very few optimization problems that have been adapted to such subjects.Here, we try to resolve this issue by exploiting studies published in [24][25][26], where medical constraints have been modeled differently with a constraint called "isoperimetric".More precisely, we propose an anti-epidemic control strategy that targets susceptible people, under the isoperimetric condition that we could not control all individuals of this category due to restricted health resources.
We consider a simple discrete-time epidemic compartmental model devised in the form of difference equations, which describe the dynamics of a discrete-time SIRS model with a temporary controlled class, meaning that the controlled people cannot acquire long-lived immunity to move towards the removed compartment due to the temporary effect of the control parameter.Thereafter, we characterize the sought optimal control, and we show the effectiveness of this limited control policy.This optimal control problem leads to the execution of two numerical methods all combined together at the same time, namely the forward-backward sweep method to generate the optimal state and control functions and the secant method adapted to the isoperimetric restriction.

Materials
Let us define a discrete-time model with the four following main compartments: • S: the number of susceptible people to infection or who are not yet infected, • C S : the number of susceptible people who are temporarily controlled, so they cannot move to the removed class due to the limited effect of control.It can represent the compartment of vaccinated people in case a vaccination is not 100% effective due to the difficulty of producing a perfect vaccine, the heterogeneity of the population or a vaccine not conferring a lifelong immunity [17,27], • I: the number of infected people who are capable of spreading the epidemic to those in the susceptible and temporarily controlled categories, • R: the number of removed people from the epidemic, but can return to the susceptible class because of the short-term removal individuals' immunity.
In our modeling approach, we aim to describe the dynamics of variables S, C S , I and R at time i based on the following difference equations: with initial conditions S 0 > 0, C S 0 ≥ 0, I 0 ≥ 0 and R 0 ≥ 0 and where Π i = µN i with N i = S i + C S i + I i + R i , gives the newborn people, aθ (0 ≤ a ≤ 1) is the recruitment rate of susceptibles to the controlled class with θ defining the control parameter as a constant between 0 and 1 (see such consideration in the case of vaccination in [27]) and "a" modeling the reduced chances of a susceptible individual to be controlled, β = δ N i with δ the infection transmission rate, µ the natural death rate, bθ (0 ≤ b ≤ 1) the recruitment rate of controlled people to the infected class even in the presence of θ with "b" modeling the reduced chances of a temporarily controlled individual to be infected, γ the recovery rate and σ the losing removal individuals' immunity rate.We note that the population size N i is constant at any time i because

Methods
Now, we consider the mathematical model (1) with θ as a discrete control function.Motivated by the desire to reduce the number of infected people as much as possible while minimizing the value of the control θ over N times, our objective is to seek an optimal control θ * such that: where J is the functional defined by: and where the control space Θ is defined by the set: A and B represent constant severity weights associated with functions I and θ, respectively.Managers of the anti-epidemic resources cannot well predict whether their control strategy will reach all the susceptible population over N times.To model the situation in which a restricted resource of control is available, we consider that the number of susceptible people we can control is equal to a constant C > 0 for N days.Hence, we try to find θ * under the definition of the following isoperimetric restriction: In [25,26], the authors defined an isoperimetric constraint on the control variable only, to model the total tolerable dosage amount of a therapy along the treatment period.In their conferences talks [28,29], Kornienko et al. and de Pinho et al. introduced state constraints in an optimal control problem that is subject to an S-exposed-I-Rdifferential system to model the situation of the limited supply of vaccine based on the work in [24] and where the isoperimetric constraint is defined on the product of the control and state variables.Our study aims to highlight more the importance of such optimal control approaches by considering a discrete model rather than a continuous one.This would be interesting since data are often collected at discrete times, as noted in [30].
In our case, to take into account the constraint (4) for the resolution of the optimal control problem (2), we consider a new variable Z defined as: with The discrete-time system of (1) becomes: In the following, we announce two theorems for proving the existence and the characterization of the sought optimal control θ * .Theorem 1. (Sufficient conditions) For the isoperimetric optimal control problem given by (2) along with the discrete state equations in (6), there exists a control θ * ∈ Θ such that J(θ * ) = min θ∈Θ J(θ).
Proof.In order to prove the existence of a solution θ * in Θ, we try to prove that min θ∈Θ J(θ) exists.
We can deduce then that in f θ∈Θ J(θ) is finite since J(θ) is bounded, and there exists a finite number j of uniformly bounded sequences θ j ∈ Θ such that lim j→∞ J(θ j ) = in f θ∈Θ J(θ) and corresponding sequences of states S j , C j S , I j and R j .Thus, there exists θ * ∈ Θ and S * , C * S , I * , R * ∈ R N such that on a subsequence, Finally, due to the finite dimensional structure of the system (6) and the objective function J(θ), θ * is an optimal control with corresponding states S * , C * S , I * and R * [26].Therefore, taking into account the structure of J being a convex function, in f θ∈Θ J(θ) is achieved.
In order to derive the necessary conditions of optimality, we employ the discrete version of Pontryagin's maximum principle stated in Theorem A1 in Appendix A.
In addition, the optimal control θ * is characterized at each iteration i by: Proof.With the application of a discrete version of Pontryagin's maximum principle in Appendix A and as done in [26,31,32], we can determine the discrete optimal control θ * for the problem (6) and its associated trajectories S * , C * S , I * and R * .
We define a discrete Hamiltonian H i as a brief notation of the function H defined for i = 0, ..., N − 1. as follows: The discrete-time adjoint system is resolved using the following formulations: that we associate with the following transversality conditions: with φ N representing the payoff term function in (3), namely AI N .Then, we obtain the following discrete-time adjoint system: with the transversality conditions λ 1,N = 0, λ 2,N = 0, λ 4,N = 0, λ 3,N = A and λ 5,N is unknown.
In order to find the transversality condition λ 5,N = constant, we use the secant-method as the appropriate numerical technique for finding the zero of the function λ 5,N → V(λ 5,N ) = ZN − Z N where ZN is the value of Z at final iteration N for various values of λ 5,N and Z N is the value fixed by C [25,33].
Since θ i is a bounded control, we can then define a Lagrangian L as follows: where ω 1,i , ω 2,i ≥ 0 ∀i verifying at Let L i be the brief notation of L and L * i be the brief notation of L at S * , C * S , I * , R * and θ * .The condition of minimization is defined as: In order to find the solution θ * i of ( * * ), we differentiate the Lagrangian L i with respect to θ i on the set Θ to obtain the optimality equation: Furthermore, we find If: If: ω 1,i = 0, therefore: implying that: Knowing that ω 2,i ≥ 0 and B > 0, we obtain θ * i ≤ If: ω 2,i = 0, thus: Using these standard optimality arguments, we characterize the control u * k by: or by a more reduced form, we can write

Numerical Results and Discussion
In this section, we resolve the discrete two-point value problem defined by System ( 6) with initial conditions along with Equations ( 7)- (11) with final conditions, using a discrete version of the forward-backward sweep method (FBSM) [26,33] with the incorporation of a discrete progressive iterative scheme to stock at each iteration i, the values of the state variables corresponding to the forward discrete-time system (6), to use them in a second discrete regressive iterative scheme incorporated for stocking at each time i, the values of the adjoint state variables corresponding to the backward discrete-time adjoint system ( 7)- (11).In fact, at each time i, the values stocked of both state and adjoint state variables were utilized in the characterization of the optimal control θ * .In brief, our algorithm is defined by the following four steps of numerical calculus (Algorithm 1).
Step 0: Guess an initial estimation of θ.
Find the optimal states S * , C * S , I * , R * and Z * , which iterate forward in the discrete two-point boundary value problem (6).
Step 2: Use the stocked values by θ and the transversality conditions λ l,N+1 for l = 1, 2, 3, 4 while searching the constant λ 5,N+1 using the secant-method.More precisely, the secant method is used to obtain the zero of the function λ 5,N → V(λ 5,N ) = ZN − Z N where ZN is the value of Z at final iteration N for various values of λ 5,N , and Z N is the value fixed by C. In addition, due its structure in (4), we choose the constant C in a way that it cannot exceed an upper bound N 0 × N where N 0 is the initial population size and N is the number of iterations.

Step 3:
Update the control utilizing new S, C S , I, R, Z and λ l for l = 1, 2, 3, 4, 5 in the characterization of θ * as presented in (12).

Step 4:
Test the convergence.If the values of the sought variables in this iteration and the final iteration are sufficiently small, check out the recent values as solutions.If the values are not small, go back to Step 1.
Figure 1 depicts the behavior of the number of susceptible people in the absence and presence of the control, and we can see that the number of susceptible people had decreased from its initial condition once the control had been introduced, while there was no significant decrease of the S function compared to the case when there was yet no control.With these parameters used, it reached only three people because of the maximal value of one taken by the optimal control θ * in almost alltimes of the control strategy, as seen in the last figure.
In Figure 2, we can well understand the increase of the number of the removed people because of a natural recovery, but it cannot represent a significant recovery because it has not reached even 14 people, and this means that only a very small number of people have been removed based on the initial condition I(0) considered.In the presence of the control, the R function increased towards a much higher number of removed people and showed it can even reach more than 31 individuals recovered from the disease in the first 17 days; this number decreased thereafter because of the results in the next figure, which will show that infection will disappear as we move forward in time.In Figure 3, the simulation shows that the number of infected people could decrease only because of a natural recovery or death, while the infection was still serious and remained present in more than 28 individuals.After the introduction of the control, the I function started to decrease once the anti-epidemic was followed, and it tended to zero values after 37 days; this means most people would recover from the disease at the end of the control strategy.As regards Figure 4, we sought to verify the condition C = 110, which meets the value 2 × S 0 , and we presented associated simulations of the number of controlled people, which increased to 53 individuals once the optimal control had been introduced and increased thereafter, again showing that it could even exceed that number, approximately towards 80 people when we go forward towards the end.In the same figure, we show the values of the optimal control θ * , which take the value of one as the maximal peak for almost alldays, and we can also see that the imposed isoperimetric constant has been verified in the final instant with some error = 3.6589.In fact, it is not evident that it reached any imposed value while verifying convergence tests of both methods used.Sometimes, the program did not stop iterating or could not show the plot because of a NANvalue, and then, the only solution was to fix the number of iterations of the secant method in which the imposed initial guess of C was approximately reached.
In Figure 5, we exhibit the value of the sought constant missing transversality condition λ 5,N , which will be essential to verify the necessary conditions announced in Theorem 2. As we can observe from this figure, the value obtained equals −1.3830 × 10 8 .Figure 6 presents a numerical simulation of the Z function when we did not seek the verification of the condition Z N = C, and we let Z N free, so we could prove that our algorithm in the case of the isoperimetric constraint helped to approximate Z N to C or even verify the equality between them, far from the value that could reach Z N when it was free; as we can see from the mentioned figure, Z N = 105.5729,which led to an important error of about 4.4271 from C that was sought in Figure 4.

Conclusions
In this paper, an optimal control approach with an isoperimetric constraint has been applied to a discrete-time SIRS model, which was in the form of a four-compartmental epidemic model where it was supposed that the controlled population did not reach the removed class due to the temporary effect of the control.The isoperimetric restriction, which has been proposed to define the number of susceptible people who receive the control along the anti-epidemic measures period, allowed us to find the optimal control needed to fight against a disease when there were limited resources.
Let us consider a discrete-time optimal control problem over times 0, ..., N, defined by: subject to the discrete-time system: X 0 given (A3) We now define the Hamiltonian function H i to be: and in the optimal control and state by H * i = H(X * i , θ * i , i).Then, based on results of the discrete version of the maximum principle discussed in [34], we can derive the following necessary conditions for our problem (A1) based on the following theorem.Theorem A1. (A discrete version of the maximum principle) Given a discrete optimal control θ * i in the sense of sufficient conditions and given solutions X * i , then the necessary conditions for θ * i to be optimal for Problem (A1)-(A4) are: X * i+1 = X * i + f (X * i , θ * i , i), X 0 given

Figure 3 .
Figure 3. Number of infected people in the absence and presence of the control in the two cases θ = 0 and θ = 0, with the same parameter values, initial conditions and severity weights as in Figure 1.

Figure 4 .
Figure 4. Number of controlled people, the optimal control θ * and the variable Z with an imposed constant C = 110, with the same parameter values, initial conditions and severity weights as in Figure 1.