www.biomathforum.org/biomath/index.php/biomath ORIGINAL ARTICLE Enumerative numerical solution for optimal control using treatment and vaccination for an SIS epidemic model Vianney Mbazumutima∗, Christopher Thron†, Léonard Todjihounde‡ ∗Department of Mathematics Institute of Mathematics and Physical Sciences (IMSP), Porto-Novo, Bénin vianney.mbazumutima@imsp-uac.org † Department of Sciences and Mathematics Texas A & M University-Central Texas, Killeen TX 76549 USA. thron@tamuct.edu ‡ Department of Mathematics Institute of Mathematics and Physical Sciences (IMSP), Porto-Novo, Bénin leonardt@imsp-uac.org Received: 14 October 2019, accepted: 13 December 2019, published: 18 December 2019 Abstract—Optimal control problems in mathe- matical epidemiology are often solved by Hamil- tonian methods. However, these methods require conditions on the problem to guarantee that they give global solutions. Because of the improved com- putational power of modern computers, numerical approximate solutions that systematically try a large number of possibilities have become practical. In this paper we give an efficient implementation of an enumerative numerical solution method for an optimal control problem, which applies to cases where standard methods cannot guarantee global optimality. We demonstrate the method on a model where vaccination and treatment are used to control the level of prevalence of an infectious disease. We describe the solution algorithm in detail, and verify the method with simulations. We verify that the enumerative numerical method produces solutions that are locally optimal. Keywords-Epidemic model, vaccination, treat- ment, optimal control, numerical method, enumer- ative method, global optimum. I. INTRODUCTION Within the field of epidemiology, extensive re- search efforts have been devoted to establishing mathematical models that accurately character- ize disease dynamics, including the effects of disease controls. Numerous mathematical tech- niques have been developed since the ground- breaking 1926 paper by A. G. M’Kendrick [37]. The most basic and most widely used models in epidemiology are multi-compartment models such as SIS (susceptible-infected-susceptible), SIR Copyright: c©2019 Mbazumutima et al. This article is distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Citation: Vianney Mbazumutima, Christopher Thron, Léonard Todjihounde, Enumerative numerical solution for optimal control using treatment and vaccination for an SIS epidemic model, Biomath 8 (2019), 1912137, http://dx.doi.org/10.11145/j.biomath.2019.12.137 Page 1 of 22 http://www.biomathforum.org/biomath/index.php/biomath https://creativecommons.org/licenses/by/4.0/ https://creativecommons.org/licenses/by/4.0/ http://dx.doi.org/10.11145/j.biomath.2019.12.137 V. Mbazumutima, C. Thron, L. Todjihounde, Enumerative numerical solution for optimal control using ... (susceptible-infected-recovered) and similar mod- els. The SIS model is suitable to model dis- eases which do not confer immunity. Some of these diseases are discussed in [1] and [28]. In the basic SIS model, individuals move between two compartments, susceptible and infected [22]. Numerous elaborations of the basic model have been developed, such as age-structured epidemic models [15], stochastic models [39], and models with vaccination [23]. Important model properties include the basic reproduction number, equilibria and stability characteristics. In mathematical models with controls, finding optimal controls is an important problem with significant practical implications. Many general numerical techniques have been developed for this purpose [8]. In [10] it is suggested a new methodology to find solution for an optimal con- trol problems with delays by shooting method which is mixed with continuation on the delay. A formulation on binary indicator functions, direct and simultaneous adaptive collocation approach to optimal control are developed in [9] while in [13] proposes a numerical solution technique for constrained optimal control problems with pa- rameters where an extension penalty function is used to adjust the state, controls, and parameter inequality constraints. In [26] an approach with Haar wavelets method is applied for finding the piecewise constant feedback controls for a finite- time linear optimal control problem of a time- varying state-delayed system. A direct method based on hybrid of block-pulse functions and Legendre polynomials is discussed in[35] while a direct collocation method is used in [49] to solve numerically optimal control problems. In [50], a sequential quadratic programming is used to solve an optimal control problem which has convex control constraints. [43] studies dynamical tunneling versus fast diffusion for a non-convex Hamiltonian and find that dynamical tunneling results at an important quicker rate than classical transport while [14] develops for certain non- convex Hamilitonian-Jacobi Equations, their ho- mogenization and non-homogenization. An algo- rithm for solving a non-convex state-dependent Hamilton-Jacobi partial differential equations is established in [12]. Researchers have also been involved in finding solutions to the infectious diseases and techniques to control them. In [40], R. M. Neilan and S. Lenhart introduce the theory of optimal control applied to systems of ordinary differential equa- tions with an application on SEIR model while in [33] S. Lenhart and J. T. Workman give an interesting overview on optimal control applied to biological models. Many authors have focused on this topic of optimal control of different diseases. Treatment in the SIS model under learning have been considered in [34], while in [19], treatment is used in the controlled SIS model, but considers dif- ferent cost structures than the earlier literature. For controlling an epidemic spread, optimal quarantine programs are used in [47], while [16] and [1] consider non-vaccine prevention in the SI and SIS models respectively. The prevention by strategic individuals in linked sub-populations is analyzed in [45], while [46] considers prevention through social distancing. The prevention and treatment in an SIS framework are considered in [17]. In [51], vaccination and treatment are used in an SIR setting and simulate optimal paths while a similar approach is discussed in [5]for an HIV type disease. In [29], it was examined the fundamental role of three type of controls, personal protection, treatment and mosquito reduction strategies in controlling malaria. In addition to prevention and treatment strate- gies and the allocation of resources available to reduce these diseases, researchers have focused in the development of techniques to evaluate which are more effective for approximating the solutions of those epidemic models. In [52], A. Zeb et al. studied the SEIR and employed the multi- step generalized differential transform method and compared the results with those obtained by the fourth-order Runge-Kutta method and non- standard finite difference method in the case of integer. In [2] F.S. Akinboro et al. used differential transformation method and variational iteration Biomath 8 (2019), 1912137, http://dx.doi.org/10.11145/j.biomath.2019.12.137 Page 2 of 22 http://dx.doi.org/10.11145/j.biomath.2019.12.137 V. Mbazumutima, C. Thron, L. Todjihounde, Enumerative numerical solution for optimal control using ... method to obtain the approximate solution of SIR model with initial condition. O. J.Peter and al.[3] use the differential transform method to study the transmission dynamics of typhoid fever diseases in a population while [38] uses differential transfor- mation method and variational iteration method in finding the approximate solution of Ebola model. In [20], it is investigated the application of matrix nonstandard finite difference schemes to obtain numerical solutions of epidemic models while [4] applies Non Standard Finite Difference (NSFD) Scheme to a modified SIR epidemic model with the effect of time delay. Laplace-Adomian Decom- position Method (LADM) in [21] is used to study the fractional order childhood disease model and shows that this method provides excellent numeri- cal solutions for nonlinear fractional order models compared to homotopy analysis, homotopy pertur- bation method, and fourth order Runge-Kutta. In [7], the Homotopy Analysis Method (HAM ) is applied and finds that it is successfully for finding the approximate solution of fractional SIR model. In [6] A. J. Arenas and al. developed the non standard finite difference scheme with Conserva- tion Law (NSFDCL) for predictorcorrector type for epidemic models and compared the result with the RungeKutta type schemes. V. K. Srivastava and al. [48] compare the solution from Euler’s and RK4 methods with those obtained by the differential transform method when they studied HIV infection of CD4+ T cells. In this paper, we will investigate the effect of the vaccination and treatment on an SIS epidemic model. In this model, neither the Pontryagin theo- rem for local optimality nor the Arrow theorem for global optimality applies, so the usual analytical methods for finding locally or globally optimal solutions cannot be used. We develop instead a numerical algorithm that first finds the overall best control from a large class of controls, and then improve it successively until local optimal conditions are satisfied. The paper is organized as follows. In Section 2, we recall basic analytical results on optimal control problems that give necessary and suffi- cient conditions for a globally optimal solution. In Section 3, we establish an SIS epidemic model under treatment and vaccination controls, and in Section 4, we formulate objective functions for the model. In Section 5 we discuss the necessary conditions for locally optimal control vaccination and treatment. Section 6 describes an enumerative numerical method for finding near-optimal con- trols, while in Section 7, we present simulation results and conclude. II. SOME BASIC RESULTS ON OPTIMAL CONTROL PROBLEMS Two of the foundational results in optimal con- trol theory are Pontryagin’s Maximum Principle and the Arrow Sufficiency Theorem. Pontryagin’s theorem gives conditions that a locally optimal so- lution must satisfy; while Arrow’s theorem guaran- tees global optimality of a locally optimal solution. These theorems are stated below. Theorem 2.1: (Pontryagin’s maximum principle)[33] If u∗(t) and x∗(t) are optimal for the problem max u J[x(t),u(t)], where J[x(t),u(t)] = ∫ tf t0 f(t,x(t),u(t))dt, (1) subject to{ dx dt = g(t,x(t),u(t)), x(t0) = x0, . (2) where the functions f and g are continuously differentiable and x(t) is piecewise differentiable. Then there exists a piecewise differentiable adjoint function λ(t) such that H(t,x∗(t),u(t),λ(t))≤H(t,x∗(t),u∗(t),λ(t)), (3) for all controls u at each time t, where the Hamil- tonian H is given by H(t,x(t),u(t),λ(t)) = f(t,x(t),u(t)) + λ(t)g(t,x(t),u(t)), (4) and{ λ′(t) = −∂H(t,x ∗(t),u∗(t),λ(t)) ∂x , λ(tf ) = 0. (5) Biomath 8 (2019), 1912137, http://dx.doi.org/10.11145/j.biomath.2019.12.137 Page 3 of 22 http://dx.doi.org/10.11145/j.biomath.2019.12.137 V. Mbazumutima, C. Thron, L. Todjihounde, Enumerative numerical solution for optimal control using ... Theorem 2.1 is proved in [44]. The following theorem gives sufficient condi- tions under which local optima are also globally optimal. Theorem 2.2: (The Arrow theorem)[11] For the optimal control problem (2), the conditions of the maximum principle are sufficient for the global minimization of J[x(t),u(t)] if the maximized Hamiltonian function H, defined in(4), is convex in the variable x for all t in the time interval [t0, tf ] for the given λ. Theorem 2.2 is proved in [24] and [36]. When the conditions of Arrow’s theorem are not satisfied, it may be difficult to guarantee global optimality of a locally optimal solution. In this case, numerical algorithms that widely explore the space of possible controls may be used to avoid converging prematurely to a local optimum that is not globally optimal. The strategy used in this paper is to find the best control from a large class of controls and then improve it successively until local optimal conditions are satisfied. While this strategy does not necessarily give the global optimum, by first doing a preliminary search over a large class of controls it does help to prevent obtaining a suboptimal solution that is only locally optimal. III. SIS MODEL UNDER TREATMENT AND VACCINATION CONTROLS The SIS (Susceptible-Infected-Susceptible) model is the model where a susceptible individual is sick and when he recovers immediately becomes susceptible again. In this basic model, each individual belongs in one of the following two states: susceptible or infectious. In the literature, many studies have been made to analyze the importance of the use vaccination and treatment on the spread of infectious diseases by using the control theory ([25]; [27]; [31]; [30]; [32]). Those optimal controls techniques play the role of limiting the spread of the infectious disease from the concerned population. In this section, we will etablish the controlled SIS system . The classical SIS epidemic model under vac- cination and treatment has four groups or com- partments, whose populations are represented by four letters: S, the number of individuals who are healthy but susceptible to the infection; I, the number of individuals who have been contami- nated and can spread the infection to susceptibles; T , the number of individuals who have undergone treatment to cure an infection; and V the number of individuals who have been vaccinated when susceptible. We also denote the total population by N, so that N = S +I +V +T . Note S,I,V,T may all vary with time, so that all are represented as functions of time. In our model we suppose that individuals enter and leave the population (either by birth/death or immigration/emigration), but the total popu- lation remains constant. Individuals leave from each group in the same proportion, but all incom- ing individuals are susceptible. We suppose that susceptibles are infected by direct contact with infected individuals, so that susceptibles’ rate of infection is proportional to the number of infected individuals. We suppose that susceptible individu- als are vaccinated at a given rate (which may de- pend on time), and infected individuals are treated and become no longer infective at a different rate (also possibly time-dependent). Finally, we assume that the vaccination is not efficacious at hundred percent, the vaccinated individuals who contact infected individuals may become reinfected at the small rate. These assumptions are represented by the following system of ordinary differential equa- tions:   dS dt = µN −βSI + γI −µS −u1S, dI dt = βSI − (µ + γ + u2)I + β�V I, dV dt = u1S −µV −β�V I, dT dt = u2I −µT, N = S + I + T + V, (6) with initial conditions S(0) =S0,I(0) =I0,V (0) =V0,T(0) =T0, (7) Biomath 8 (2019), 1912137, http://dx.doi.org/10.11145/j.biomath.2019.12.137 Page 4 of 22 http://dx.doi.org/10.11145/j.biomath.2019.12.137 V. Mbazumutima, C. Thron, L. Todjihounde, Enumerative numerical solution for optimal control using ... and control conditions 0 ≤ ui(t) ≤ uimax ≤ 1, i = 1, 2. The parameters in (6) have the following sig- nificances: µ and γ are the population replace- ment rate and the recovery rate from infection respectively; β = λ N is the disease transmission coefficient, so that 1/β is the average number of infective contacts per unit time which result in the susceptible individual becoming infected; � is the fraction of vaccinated individuals for whom the vaccine is ineffective; and u1(t) and u2(t) are the proportionate vaccination and treatment control levels, so that for example u1(t)S(t) is the number of susceptible individuals vaccinated per unit time at time t. We may rewrite the system (6) in matrix format. Define ~x(t) ≡ (S(t),I(t),V (t),T(t))T , (8) and define the matrix A(~x(t)) as follows: A(~x(t)) = −βI(t)−u1(t)−µ γ 0 0βI(t) −(µ+γ+u2(t)) β�I(t) 0 u1(t) 0 −(µ+β�I(t)) 0 0 u2(t) 0 −µ  , (9) where N = S0 + I0 + V0 + T0. Then system (6) can be written as d dt ~x(t) = A(~x(t))~x(t) +   µN 0 0 0   , (10) with initial conditions ~x(0) = (S0,I0,V0,T0). Definition 3.1: [18]. Metzler matrices are square (real) matrices in which all the off-diagonal elements are non-negative: aij ≥ 0,∀i 6= j. The matrix A is a Metzler matrix according to Definition (3.1). Proposition 3.1: The set Γ = {(S,I,V,T) ∈ R4+} is positively invariant for the system (6). Proof: Whenever xj = 0 and xi ≥ 0 for i 6= j, since A is Metzler matrix it follows that dxj/dt ≥ 0. It follows that if xi(0) ≥ 0 ∀i, then none of the xi(t) will change sign for t ≥ 0. Therefore, we conclude that the system (6) determined by the matrix(9) preserves invariance of the non-negative cone R4+. IV. OBJECTIVE FUNCTIONS FORMULATION In this section, we establish objective functions in case of nonlinear cost function and the piece- wise linear cost function. A. Original objective function We suppose a nonlinear cost function in order to take into account various types of costs affecting the Susceptible and Infected populations. The cost function for system (6) takes the following form: J(u1,u2,S,I) =∫ tf 0 [ f1(u1(t),S(t))+f2(u2(t),I(t)) ] dt+zI(tf ), (11) where f1(u1,S) ={ c′0 if u1 = 0, c0 + c1u1S + c2u 2 1 if 0 < u1 ≤ u1max, (12) and f2(u2,I) ={ d′0 + d2I if u2 = 0, d0 + d1u2I + d2I if 0 < u2 ≤ u2max, (13) The different terms in (11)-(13) are motivated as follows. The functions f1 and f2 represent cost rates associated with vaccination and treatment re- spectively, while the final additive term in (11) rep- resents costs attributed to latent infections which remain after the treatment period is complete. As far as vaccination cost, we expect a fixed, low-level maintenance cost rate during time peri- ods when no active vaccination efforts are being made: this fixed cost rate is represented by the constant c′0 in (12). When vaccination efforts are being prosecuted, higher fixed costs are incurred, including salaries and facilities: this higher cost level is represented by c0 in (12), where c0 > c′0. We suppose that each vaccination has a fixed Biomath 8 (2019), 1912137, http://dx.doi.org/10.11145/j.biomath.2019.12.137 Page 5 of 22 http://dx.doi.org/10.11145/j.biomath.2019.12.137 V. Mbazumutima, C. Thron, L. Todjihounde, Enumerative numerical solution for optimal control using ... cost c1, which multiplies the total vaccination rate u1S in (12). Finally, the quadratic term c2u21 is included to model diminishing returns of higher levels of vaccination effort, including the cost of vaccinating harder-to-reach or less cooperative individuals. As far as treatment cost, as with vaccination we also expect a low-level maintenance cost rate even when no treatments are given, which is modeled by the constant d′0 term in (13). There is also a cost rate per infected individual (which may include both financial cost and quality-of-life cost), expressed by the coefficient d2 in (13). Lastly, each patient who receives treatment has a cost rate d2, which multiplies the treatment rate u2I to give the non-fixed cost rate associated with patient care. The values u1max and u2max represent maxi- mum population penetration rates for vaccination and treatment, respectively. For example, a value u1max = 0.05 indicates that at most 5% of the susceptible population can be vaccinated per basic time unit (which is typically taken as days); while a value u2max = 0.1 means that at most 10% of the infected population is treated per day. In summary, the optimization problem is to find the controls (u∗1,u ∗ 2) that minimize the objective function: (u∗1,u ∗ 2) = argmin u1,u2 J(u1,u2,S,I), (14) where (u1,u2) ∈ U such that 0 < u1 < u1max and 0 < u2 < u2max. (15) In this problem, the cost function is not contin- uously differentiable, and the Hamiltonian is not convex. So neither Theorem 2.1 nor Theorem 2.2 can be applied in this case. B. Piecewise linear objective function Since we are using an enumerative approach to numerical solution, it is preferred to have a finite number of possible optimal control values. In Section V, we will show that the optimal treatment level is either 0 or u2max, given the cost function (13). However, with the cost function (12) there is an infinite number of possible optimal values. For this reason, we consider a modification of the vaccination cost function f1 which closely approximates (12), but which leads to a finite number of optimal vaccination levels (as will be shown in the next section). The modified function is piecewise linear with the following mathematical form: f1(u1,S) = (16)  c′0 if u1 = 0, c0 + c1u1S + c ′ 2u1 + c ′ 3(u1 −u1mid) +, if 0 < u1 ≤ u1max, where the function x+ is the ramp function, x+ = max(x, 0). (17) The constants u1mid,c′2 and c ′ 3 may be chosen to approximate the quadratic term c2u21 in (12). Figure 1 shows various possibilities for piecewise linear approximations. Given the constant c2, the values u1mid = 0.5, c′2 = 0.5c2 and c ′ 3 = c2 produce a cost function that is an upper bound to the quadratic term c2u21; the values u1mid = 5/9, c′2 = 0.2c2, and c ′ 3 = 2c2 produce an effective lower bound; and the values u1mid = 0.5, c′2 = ( √ 2−1)c2 and c′3 = (3− √ 2)c2 gives an approx- imation that minimizes the maximum deviation (maximum deviation is equal to 0.043c2). Using these different functions, upper and lower bounds on the optimal cost for the original quadratic cost function (12) may be obtained. V. LOCAL OPTIMALITY CONDITIONS FOR THE NON-AUTONOMOUS VACCINATION AND TREATMENT CONTROL PROBLEM A solution ~x(t) with controls u1(t) and u2(t) is locally optimal if small perturbations of the controls u1 and u2 during small time intervals never decrease the cost. This means that there is no way to improve the solution by making slight adjustments to the controls. Local optimization is applicable to the global problem in that a globally optimal solution must also be locally optimal. Hence local optimality is a necessary condition for global optimality. In this section, we explicitly calculate the effect of local changes in the controls Biomath 8 (2019), 1912137, http://dx.doi.org/10.11145/j.biomath.2019.12.137 Page 6 of 22 http://dx.doi.org/10.11145/j.biomath.2019.12.137 V. Mbazumutima, C. Thron, L. Todjihounde, Enumerative numerical solution for optimal control using ... Figure 1. Approximation of the quadratic cost function with piecewise linear functions. The linear functions shown give upper and lower bounds to the quadratic cost function, as well as a best approximation that minimizes the maximum deviation between the quadratic cost function and piecewise linear approximation. u1(t) and u2(t) on the cost function. Changes in the two controls are considered separately, because in practice they can be varied independently. This gives us a way to identify local optima, which correspond to control levels for which differential changes yield no improvement. In previous sections, we have considered the au- tonomous problem in which the parameters β,�,γ, and λ and also the costs c0,c′0,d0,d ′ 0,c1,c2,d1,d2 are constants independent of time. In this section we consider the more general non-autonomous problem, in which all parameters can be continu- ous function of time. The reader should understand that parameters and costs in this section now rep- resent time dependent functions (e.g. β represents β(t)). A. Necessary conditions for the optimal control vaccination with simplified cost function First, we will perturb the control u1(s) and calculate the effect on the cost function. Given a control u1(t), the perturbed control u′1(s,t) differs from u1(t) by a small amount on an interval of length δ, as follows: u′1(s,t) ={ u1(t) + du1 for s < t < s + δ, u1(t) otherwise, (18) where s is a fixed value between 0 and tf . The system evolution x′(s,t) corresponding to controls u′1(s,t) and u2(t) may be written ~x ′(s,t) ≡ (S′(s,t),I′(s,t),V ′(s,t),T ′(s,t))T , (19) and satisfies the equations ~x′(s, 0) = ~x(0), ∂ ∂t ~x′(s,t) =   A(~x′(t))~x′(s,t) for t < s and t > s + δ, A(~x(t))~x′(s,t)−S′(s,t)du1~∆10 for s ≤ t ≤ s + δ, (20) Biomath 8 (2019), 1912137, http://dx.doi.org/10.11145/j.biomath.2019.12.137 Page 7 of 22 http://dx.doi.org/10.11145/j.biomath.2019.12.137 V. Mbazumutima, C. Thron, L. Todjihounde, Enumerative numerical solution for optimal control using ... where −→ ∆10 ≡ (1, 0,−1, 0)T . (21) System (20) corresponds to system (10) with per- turbed control u′1(s,t) on the interval s ≤ t ≤ s + δ. Note that ~x′(s,t) = ~x(t) for t ≤ s since ~x′(s, 0) = ~x(0) and ~x′(s,t) and ~x(t) satisfy the same differential equation on [0,s]. In order to compute the new cost associated with the per- turbed solution x′(s,t) with controls u′1(s,t) and u2(t), we need to solve (20) for x′(s,t). Since equations (20) have different forms on the intervals s < t < s + δ and t > s + δ, we will treat these two intervals separately. First we find the solution on s ≤ t ≤ s + δ. For t > s, we may define: (∆1S(s,t), ∆1I(s,t), ∆1V (s,t), ∆1T (s,t)) ≡ 1 a [ S(t) −S′(s,t),I(t) − I′(s,t), V (t) −V ′(s,t),T(t) −T ′(s,t) ] , (22) where a = S(s)du1δ. (23) We also use the notation: ~∆1(s,t)≡(∆1S(s,t),∆1I(s,t),∆1V(s,t), ∆1T(s,t))T. (24) Note that ~∆1(s,s) = ~x(s) −~x′(s,s) = ~0. For simplicity, in the remaining discussion we will write the functions S,S′,I,I′,V,V ′,T,T ′, ∆1, ∆1, ∆2 without arguments (e.g. we write I′(s,t) as I′). From (10) and (20) we have for interval s < t < s + δ,  dS dt =µN −βSI + γI −µS −u1S, dS′ dt =µN−βS′I′+γI′−µS′−u1S′−S′du1. (25) Using the notations (23) and (24) for the system (25), we obtain for s < t < s + δ dS dt − dS′ dt =−aβ(I∆1S +S∆1I)+aγ∆1I −aµ∆1S−au1∆1S +S′du1 +O(a2). (26) Using similar computations, we may obtain cor- responding equations for infected, vaccinated, and treated compartments on the interval s < t < s+δ: dI dt − dI′ dt =aβ(∆1SI+∆1IS)−(µ+γ+u2)a∆1I + β�(a∆1IV + a∆1V )I + O(a2). (27) dV dt − dV ′ dt =au1∆1S−aµ∆1V −aβ�(∆1V I+∆1IV )−S′du1 +O(a2). (28) dT dt − dT ′ dt = au2∆1I −aµ∆1T . (29) Equations (26)-(29) may be rewritten in matrix form: ∂ ∂t [ ~∆1(s,t) ] =(−βI−u1−µ −βS+γ 0 0 βI −(µ+γ+u2)+βS+β�V β�I 0 u1 −β�V −µ−β�I 0 0 u2 0 −µ ) ~∆1(s,t) + ~∆10 δ + O(a). (30) Using (30), we have on t ∈ [s,s + δ] ∂~∆1(s,t) ∂t =Ã(~x(t))~∆1(s,t)+ 1 δ −→ ∆10 +O(a), (31) where Ã(~x(t)) = A(~x(t)) + ( 0 −βS 0 0 0 β(S+�V ) 0 0 0 −β�V 0 0 0 0 0 0 ) . (32) We also note ~∆1(s,s) = ~0. To lowest order in δ, the solution of (31) gives ~∆1(s,s + δ) ≈ ~∆10. On the interval t ∈ [s+δ,tf ], ~∆1(s,t) is a solution to ∂~∆1(s,t) ∂t = Ã(~x(t))~∆1(s,t) (33) We are now ready to compute the difference between cost functions for x′ and x. We denote these cost functions by J′ and J respectively. For simplicity, we first consider the case where c′0 = c0 and d′0 = d0. From (11),(13) and (16), we obtain the sequence of equations (34)-(36), where H(x) in (36) denotes Biomath 8 (2019), 1912137, http://dx.doi.org/10.11145/j.biomath.2019.12.137 Page 8 of 22 http://dx.doi.org/10.11145/j.biomath.2019.12.137 V. Mbazumutima, C. Thron, L. Todjihounde, Enumerative numerical solution for optimal control using ... J′ = ∫ s 0 [ [c0 + c1u1(t)S(t) + c ′ 2u1 + c ′ 3(u1 −u1mid) +] + [d0 + d1u2(t)I(t) + d2I(t)] ] dt + ∫ s+δ s [ c0 + c1(u1(t) + du1)S ′(s,t) + c′2(u1(t) + du1) + c ′ 3((u1 + du1) −u1mid) + + d0 + (d1u2(t) + d2)I ′(s,t) ] dt + ∫ tf s+δ [ [c0 +c1u1(t)S ′(s,t)+c′2u1 +c ′ 3(u1−u1mid) +]+[d0 +(d1u2(t)+d2)I ′(s,t)] ] dt + I′(s,tf ); (34) J = ∫ s 0 [ [c0 + c1u1(t)S(t) + c ′ 2u1 + c ′ 3(u1 −u1mid) +] + [d0 + (d1u2(t) + d2)I(t)] ] dt + ∫ s+δ s [ [c0 + c1u1(t)S(t) + c ′ 2u1 + c ′ 3(u1 −u1mid) +] + [d0 + (d1u2(t) + d2)I(t)] ] dt + ∫ tf s+δ [ [c0 + c1u1(t)S(t) + c ′ 2u1 + c ′ 3(u1 −u1mid) +] + [d0 + (d1u2(t) + d2)I(t)] ] dt + I(tf ); (35) J′−J du1 = ∫ s+δ s [ −S(s)δ (c1u1(t)∆1S(s,t)+[d1u2(t)+d2]∆1I(s,t))+c1S′1(s,t)+c ′ 2 +c ′ 3H(u1−u1mid) ] dt −S(s)δ [∫ tf s [c1u1(t)∆1S(s,t) + (d1u2(t) + d2)∆1I(s,t)] ] dt−S(s)δ∆1I(s,tf ). (36) the Heaviside (step) function (note the Heaviside function is the derivative of the ramp function). Taking the limit as du1 → 0 for small δ we obtain ∂J ∂u1 =δ [ c1S+c ′ 2 +c ′ 3H (u1(s)−u1mid)−Ψ1 ] , (37) where Ψ1≡S(s) [∫ tf s [ c1u1(t)∆1S(s,t) + (d1u2(t)+d2)∆1I(s,t) ] dt+∆1I(s,tf ) ] . (38) In this case, the local optimum conditions on vaccination control u1(s) depend on 3 cases: (a) Ψ1 ≤ c1S + c′2: then ∂J ∂u1(s) ≥ 0 and J(u1(s)) is minimized when u1(s) = 0. (b) c1S + c′2 ≤ Ψ1 ≤ c1S + c ′ 2 + c ′ 3: in this case, it is necessary to check whether J(u1mid) < J(0). Using (37), we obtain J(u1mid) −J(0+) = δu1mid[c1S + c′2 − Ψ1] and J(0+) −J(0) = (c0 − c′0)δ. So, J(u1mid) < J(0) is equivalent to Ψ1 > (c0 − c′0)/u1mid + c ′ 2 + c1S. If J(u1mid) − J(0) < 0 then u1(s) = u1mid is locally opti- mal, otherwise u1(s) = 0 is locally optimal. (c) Ψ1 ≥ c1S+c′2+c ′ 3: in this case, it is necessary to check whether J(u1max) < J(0). From Biomath 8 (2019), 1912137, http://dx.doi.org/10.11145/j.biomath.2019.12.137 Page 9 of 22 http://dx.doi.org/10.11145/j.biomath.2019.12.137 V. Mbazumutima, C. Thron, L. Todjihounde, Enumerative numerical solution for optimal control using ... (37), we have J(u1max)−J(0+) = δ[u1max(c′2 − Ψ1) +c′3(u1max −u1mid)] and J(0+) −J(0) = δ(c0 − c′0). So J(umid) < J(0) is equivalent to Ψ1 > c1S + c ′ 2 + c ′ 3(1−u1mid/u1max) + δ(c0−c ′ 0). If J(u1max) −J(0) < 0 then u1(s) = u1max is optimal otherwise if J(u1max) −J(0) > 0 then u1(s) = 0 is optimal. B. Necessary conditions for the optimal control treatment To find the necessary conditions for the optimal control treatment function u2(t), we will perturb the control u2(s) and calculate the effect on the cost function. The perturbed control u′2 which differs from u2 by a small amount on an interval of length δ is u′2(s,t) = { u2(t) + du2 for s ≤ t ≤ s + δ, u2(t) otherwise . (39) For this purpose, we define ~x ′(s,t) ≡ (S′(s,t),I′(s,t),V ′(s,t),T ′(s,t))T , (40) which satisfies the equations: ~x ′(s, 0) = ~x(0), ∂ ∂t ~x ′(s,t) ={ A(~x(t))~x ′(s,t) for t < s or t > s + δ, A(~x(t))~x ′(s,t)−I′(s,t)du2~∆20 for s≤t≤s+δ, (41) where −→ ∆20 ≡ (0,−1, 0, 1)T . (42) System (41) corresponds to system (10) with perturbed control u′2(s,t) on the interval s ≤ t ≤ s + δ. Note that ~x′(s,t) = ~x(t) for t ≤ s since ~x′(s, 0) = ~x(0) and ~x′(s,t) and ~x(t) satisfy the same differ- ential equation on [0,s]. In order to compute the new cost associated with the perturbed solution x′(s,t) with controls u′2(s,t) and u1(t), we need to solve (20) for x′(s,t). Since equations (20) have different forms on the intervals s < t < s + δ and t > s + δ, we will solve for x′(s,t) on these two intervals separately. For s ≤ t ≤ s + δ, we define : (∆2S(s,t), ∆2I(s,t), ∆2V (s,t), ∆2T (s,t)) ≡ 1 b [S(t)−S′(s,t),I(t)−I′(s,t),V (t)−V ′(s,t), T(t) −T ′(s,t)], (43) where b ≡ I(s)du2δ. The following notation will be also used: ~∆2(s,t)≡(∆2S(s,t), ∆2I(s,t), ∆2V(s,t), ∆2T(s,t)) T . (44) Note that ~∆2(s,s) = ~0. Following the same arguments used from (25) to (30), to first order for s < t < s + δ, we obtain ∂~∆2(s,t) ∂t = Ã(~x(t))~∆2(s,t) + 1 δ −→ ∆20, (45) and for t > s + δ ∂~∆2(s,t) ∂t = Ã(~x(t))~∆2(s,t). (46) To lowest order in b, the solution of (45) gives ~∆2(s,s + δ) ≈ ~∆20, (47) and ~∆2(s,t) is solution to  ∂~∆2(s,t) ∂t = Ã(~x(t))~∆2(s,t) + 1 δ ~∆20 for t ∈ (s,s + δ), ∂~∆2(s,t) ∂t = Ã(~x(t))~∆2(s,t) for t ∈ [0,s] ∪ [s + δ,∞). By the same argument which has been used for (34) and (35), we have the series of equations (48)- (50), where b ≡ I(s)du2δ in (50). Biomath 8 (2019), 1912137, http://dx.doi.org/10.11145/j.biomath.2019.12.137 Page 10 of 22 http://dx.doi.org/10.11145/j.biomath.2019.12.137 V. Mbazumutima, C. Thron, L. Todjihounde, Enumerative numerical solution for optimal control using ... J = ∫ s 0 [ [c0 + c1u1(t)S(t) + c ′ 2u1(t) + c ′ 3(u1 −u1mid) +] + [d0 + (d1u2(t) + d2)I(t)] ] dt + ∫ s+δ s [ [c0 + c1u1(t)S(t) + c ′ 2u1(t) + c ′ 3(u1 −u1mid) +] + [d0 + (d1u2(t) + d2)I(t)] ] dt + ∫ tf s+δ [ [c0 + c1u1(t)S(t) + c ′ 2u1(t) + c ′ 3(u1 −u1mid) +] + [d0 + (d1u2(t) + d2)I(t)] ] dt + I(tf ); (48) J′ = ∫ s 0 [ [c0 + c1u1(t)S(t) + c ′ 2u1(t) + c ′ 3(u1 −u1mid) +] + [d0 + d1u2(t)I(t) + d2I(t)] ] dt + ∫ s+δ s [ [c0 + c1u1(t)S ′(s,t) + c′2u1(t) + c ′ 3(u1 −u1mid) +] + [d0 + d1(u2(t) + du2)I ′(s,t) + d2I ′(s,t)] ] dt + ∫ tf s+δ [ [c0 + c1u1(t)S ′(s,t) + c′2u1(t) + c ′ 3(u1 −u1mid) +] + [d0 + d1u2(t)I ′(s,t) + d2I ′(s,t)] ] dt + I′(s,tf ); (49) J′ −J du2 = ∫ s+δ s [ −S(s)δ (c1u1(t)∆2S(s,t) + [d1u2(t) + d2]∆2I(s,t)) + d1I′(s,t) ] dt − I(s)δ [∫ tf s+δ [(d1u2(t) + d2)∆2I(s,t) + c1u1(t)∆2S(s,t)]dt + ∆2I(s,tf ) ] . (50) Taking limits as before, we obtain 1 δ ∂J ∂u2(s) = I(s)(d1 − Ψ2), (51) where Ψ2≡ ∫ tf s [ (d1u2(t)+d2)∆2I(s,t) + c1u1(t)∆2S(s,t) ] dt −∆2I (s,tf ). (52) For this case, there is no effect from the control u2(s) on ∂J ∂u2(s) · The local optimum conditions on u2(s) depends 3 cases: (a) d1 > Ψ2 =⇒ ∂J∂u2(s) > 0, then J(u2(s)) is minimized when u2(s) = 0, (b) d1 = Ψ2 =⇒ ∂J∂u2(s) = 0, then u2(s) has no effect on J(u2(s)), (c) d1 < Ψ2 =⇒ ∂J∂u2(s) < 0, then J(u2(s)) is minimized at u2 when u2(s) = u2max. C. Summary of necessary conditions for optimal controls The following theorem summarizes the result of the previous discussion: Theorem 5.1: Suppose u∗1 and u ∗ 2 are locally optimal controls for the system (6) with objective function (11),where c0 > c′0 and d0 > d ′ 0. Let us consider Ψ1(s) and Ψ2(s) given by the expressions (38) and (52) respectively. Then (a) If Ψ1(s) < 0 or S(s) = 0 then u1(s) = 0 is locally optimal, Biomath 8 (2019), 1912137, http://dx.doi.org/10.11145/j.biomath.2019.12.137 Page 11 of 22 http://dx.doi.org/10.11145/j.biomath.2019.12.137 V. Mbazumutima, C. Thron, L. Todjihounde, Enumerative numerical solution for optimal control using ... (b) If c1S + c′2 ≤ Ψ1 ≤ c1S + c ′ 2 + c ′ 3 and Ψ1 > (c0−c′0)/u1mid + c1S + c ′ 2 and S(s) > 0 then u∗1(s) = u1mid is locally optimal, (c) If c1S + c′2 ≤ Ψ1 ≤ c1S + c ′ 2 + c ′ 3 and Ψ1 < (c0 −c′0)/u1mid + c1S + c ′ 2 then u ∗ 1(s) = 0 is locally optimal, (d) If Ψ1 ≥ c1S + c′2 + c ′ 3 and Ψ1 > c1S + c ′ 2 + c′3(1−u1mid/u1max)+δ(c0−c ′ 0) and S(s) > 0 then u∗1(s) = u ∗ 1max is locally optimal, (e) If Ψ1 ≥ c1S + c′2 + c ′ 3 and Ψ1 < c1S + c ′ 2 + c′3(1−u1mid/u1max)+δ(c0−c ′ 0) then u ∗ 1(s) = 0 is locally optimal, (f) If Ψ2 < d1 or I(s) = 0 then u∗2(s) = 0 is locally optimal, (g) If Ψ2 > d1 or I(s) > 0 and Ψ2 > (d0−d′0) I(s)u2max + d1 then u∗2(s) = u2max is locally optimal, (h) If Ψ2 > d1 or I(s) > 0 and Ψ2 = (d0−d′0) I(s)u2max + d1 then u∗2(s) = 0 or u ∗ 2(s) = u2max is locally optimal, (i) If Ψ2 > d1, I(s) > 0 and Ψ2 < (d0−d′0) I(s)u2max +d1 then u∗2(s) = 0 is locally opimal, (j) If I(s) = 0 then the value of u2(s) has no effect on the solution to system (6) or on the value of J. Corollary 5.1: For the system (6) with objec- tive function (11), given any locally optimal con- trol (u∗1,u ∗ 2), then for any 0 ≤ s ≤ tf , we have u∗1 ∈{0,u1mid,u1max} and u∗2 ∈{0,u2max} . Proof: Follows immediately from Theorem 5.1. Corollary 5.1 implies that at any given time instant, there are only 6 possible optimal control pairs. So, given that an optimal control is constant on a specified set N of intervals, then there are 6N of possible optimal controls. This fact is based to the discussion in the next section. VI. EFFICIENT NUMERICAL METHOD FOR FINDING NEAR-OPTIMAL CONTROLS Theorem 5.1 gives conditions for local optimal- ity, which does not necessary imply global opti- mality. Many algorithms employ a process which converges to a locally optimal solution given a starting point. If the starting point is close enough to a globally optimal solution, these algorithms will converge to a global optimum. Thus, it is important to identify near-optimal solutions. One way on doing this is to find optimal solutions from a large class of controls that are representative of the different possibilities. First, consider the class of control strategies that are constant in intervals of length T/N where T is the total time of the system and N is a positive in- teger. We also consider strategies that are restricted to the optimal values specified in Theorem 5.1. Then there are 6N strategies of which meet these conditions. If N is small enough, the best solution from this class can be found by simply evaluating the cost for all 6N strategies. This limits on the size of N for practical computation. In order to increase N, we may make further assumptions. We expect the vaccination level u1mid to occur as the system is transitioning from no control to full control. In order to reduce computation time, we consider only the two extreme vaccination strategies: 0,u1max. This leads to 4N strategies that are constant on the N intervals. On each interval 0 ≤ k ≤ N − 1, there are four strategy options: (u1k,u2k) = {(0, 0), (u1max, 0), (0,u2max), (u1max,u2max)}. These options may be indexed as follows: (u1k,u2k)∈(0, 0) ⇐⇒ aN−k−1 = 0, (u1k,u2k) = (u1max, 0) ⇐⇒ aN−k−1 = 1, (u1k,u2k) = (0,u2max) ⇐⇒ aN−k−1 = 2, (u1k,u2k) = (u1max,u2max) ⇐⇒ aN−k−1 = 3, where aN−k−1 is the index of the control on inter- val k with k = 0, · · · ,N−1. Then (a0, · · · ,aN−1) completely specifies the control. We associate this control with the index aN−k−14 N−k−1 + aN−k−24 N−k−2 + · · · + a0. Schematically, we have the Figure 2. For more explanation of the numerical method used, we consider the flowchart in Figure 3. Biomath 8 (2019), 1912137, http://dx.doi.org/10.11145/j.biomath.2019.12.137 Page 12 of 22 http://dx.doi.org/10.11145/j.biomath.2019.12.137 V. Mbazumutima, C. Thron, L. Todjihounde, Enumerative numerical solution for optimal control using ... Figure 2. Strategies computed on kth interval. Figure 3. Algorithm’s flowchart to find near-optimal strategy An enumerative algorithm for finding a near- optimal strategy from among these strategies pro- ceeds as follows: Algorithm 1 calculates the opti- Algorithm 1 Calculation of optimal piecewise constant control strategy part I. 1: N = number of intervals 2: Cbest = 1E100 3: for j in 0, · · · , 4N − 1 do 4: Generate strategy Sj = (a (j) 0 , · · · ,a (j) N−1) associated with index j 5: Evaluate Cj = cost of strategy Sj 6: if Cj < Cbest then 7: Cbest ← Cj 8: Sbest ← Sj 9: end if 10: end for 11: return Cbest,Sbest mal piece-wise constant control strategy by com- puting the costs of all 4N possible strategies. This is computationally expensive. It is possible to greatly reduce costs by reusing cost computations between strategies as follows. Suppose S1 and S2 are two strategies which agree on the first k intervals. This means that the cost contribu- tions from first k intervals are the same for both strategies. If these interval costs are known for S1, then it is not necessary to recompute them for S2. Thus in the calculation for S2, it is only necessary to compute the costs from the last N−k intervals. It may be shown that on average, only 2 intervals’ costs must be recomputed, instead of N intervals as in Algorithm 1. It follows that the total amount of computation required is reduced by a multiplicative factor of 2/N. A pseudo-code for the improved algorithm is shown in Algorithm 2. After part I, we consider on flowchart part II.a, where the intervals are divided in half to obtain 2N intervals. As shown in Figure 4, we keep the u20, · · · ,u2N−1 the same. We recompute u′10, · · · ,u ′ 12N−1 by evaluating the 2N strategies and choosing the optimum. The resulting strategy has vaccination constant on 22N intervals and treatment constant on N intervals. We denote the best vaccination strategy obtained from the previous procedure by (u∗10, . . . ,u ∗ 12N−1). After this, the next step shown Biomath 8 (2019), 1912137, http://dx.doi.org/10.11145/j.biomath.2019.12.137 Page 13 of 22 http://dx.doi.org/10.11145/j.biomath.2019.12.137 V. Mbazumutima, C. Thron, L. Todjihounde, Enumerative numerical solution for optimal control using ... Algorithm 2 Improved version of Algorithm 1. 1: N = number of intervals 2: j = 0 3: Generate strategy S0 = (a (0) 0 , · · · ,a (0) N−1) 4: Compute cost for each interval C0(0), · · · ,C0(N − 1) 5: Cbest = C0(0) + · · · + C0(N − 1) 6: Sbest = S0 7: for j = 1, · · · , 4N − 1 do 8: Compute largest ` such that a(j)` 6= a (j−1) ` 9: Cj(0) = Cj(k),k = 0 . . .`− 1 10: Recompute Cj(`), · · · ,Cj(N − 1) 11: if ∑ Cj(k) < Cbest then 12: Cbest ← ∑ Cj 13: Sbest ← (a (j) 0 , · · · ,a (j) N−1) 14: end if 15: end for 16: return Cbest,Sbest Figure 4. Schema of flowchart part II.a. in flowchart part II.b is to divide the intervals of constant treatment in half and recompute u′20, · · · ,u ′ 22N−1 keeping the vaccination control u∗10, · · · ,u ∗ 12N−1 fixed. The resulting strategy is represented in Figure 5. Figure 5. Schema of flowchart part II.b. At the end of part II, the resulting strategy is constant on 2N intervals of equal length. So in part III, the purpose is to improve the solution by adjusting the sizes of active treatment and vaccination intervals. Each iteration of this part III changes each active treatment or vaccination inter- val by at most 1. Also, this part III works by con- sidering all strategies that agree with the previous best strategy except at the endpoints of the active treatment or vaccination intervals. Figures 6 and 7 show how different vaccination and treatment controls are tried which differ from the previous best solution only where active treatment intervals begin or end. Note that the dashed line in Figures 6 and 7 represents the previous optimal vaccination and previous optimal treatment respectively. Part III of the flowchart Figure 3 has the fol- lowing pseudocode: Algorithm 3 Pseudocode for Part III of algorithm. 1: S0 = (a (0) 0 , · · · ,a (0) T/dt−1) is the previous best strategy 2: Identify change points in vaccination strategy. N1 = number of change points. 3: j = 0 4: time steps before and after changes =N1 5: Identify change points in treatment strategy. N2 = number of change points 6: compute the costs on all 2N intervals for pre- vious optimal strategy C0(0), · · · ,C0(2N−1) 7: Sbest = S0 8: for j = 1, · · · , 3N1 2N2 − 1 do do 9: Compute next candidate strategy Snew = a (j) 1 , . . .a (j) T/dt 10: Compute the first ` such that a(j)` 6= a (j−1) ` 11: Cj(0) = Cj(k),k = 0 . . .`− 1 12: Recompute Cj(`), · · · ,Cj(N − 1) 13: ∑ Cj(k) ← Cnew 14: if ∑ Cj(k) < Cbest then 15: Cbest ← Cnew 16: Sbest ← Snew 17: end if 18: end for 19: return Cbest,Sbest Biomath 8 (2019), 1912137, http://dx.doi.org/10.11145/j.biomath.2019.12.137 Page 14 of 22 http://dx.doi.org/10.11145/j.biomath.2019.12.137 V. Mbazumutima, C. Thron, L. Todjihounde, Enumerative numerical solution for optimal control using ... Figure 6. Vaccination strategy modification in stage III of algorithm. Figure 7. Treatment strategy modification in stage III of algorithm. VII. NUMERICAL SIMULATIONS AND DISCUSSION A. Numerical simulations In this section, we use this algorithm described in Section VI to solve the control problem given by Equations (11),(12) and (16). Table I, shows the baseline parameters for our simulations. In our analysis, we used the baseline configuration described in Table I and varied one cost parameter at a time. Figures 8 - 14 show the results for the parameters c0,c1,c2,d0,d1,d2 and z respec- tively. In each set of three sub-figures, the first two sub-figures give the optimal vaccination and treatment controls found by the algorithm. For all solutions found by the algorithm, local opti- mality was numerically verified. The third sub- figure shows the process of convergence during the algorithm. In each graph, the thin solid line shows the default configuration, and the other two Biomath 8 (2019), 1912137, http://dx.doi.org/10.11145/j.biomath.2019.12.137 Page 15 of 22 http://dx.doi.org/10.11145/j.biomath.2019.12.137 V. Mbazumutima, C. Thron, L. Todjihounde, Enumerative numerical solution for optimal control using ... Table I BASELINE PARAMETERS USED IN SIMULATION Parameters Interpretation Values units Source c1 Cost coefficient of u1S 10 dimensionless assumed c2 Cost coefficient of u21 10 5 dimensionless assumed d1 Cost coefficient of u2I 40 dimensionless assumed d2 Cost coefficient of I 5 dimensionless assumed z Final cost coefficient of I 50 dimensionless assumed c′2 Approximative constant of c2 when 0 < u1 ≤ u1mid 0 dimensionless assumed c′3 Approximative constant of c2 when u1 ≤ u1mid ≤ u1max 0 dimensionless assumed nCint Number of intervals on which control is constant 3 dimensionless assumed β Disease transmission coefficient 8×10−5 day−1 assumed γ Recovery rate 0.65 ∈ [0.25,1.5] day−1 [42] � Small rate infection of vaccinated individuals 0.0001 day−1 assumed µ rate of replacement including both birth/death and immigration/emigration 0.004 day−1 [41] u1max Maximum population rate for vaccination 0.05 day−1 assumed u2max Maximum population rate for treatment 0.1 day−1 assumed tf Final time 100 days assumed dt Time increment 0.1 day assumed N Population size 10000 humans assumed lines give configurations with different values of the chosen parameter. The labels I.a, II.a,b and III on each third subfigure correspond to different parts of the algorithm described in the flowchart in Figure 3. B. Discussion Figures 8, (9) and (10) show that increasing c0,c1 or c2 (which are fixed cost of vaccination and cost per vaccination respectively) reduces the active vaccination interval and increases the active treatment interval. For Figures (10) shows that no matter how much the quadratic vaccination cost c2 is increased, the vaccination time interval is reduced while the active treatment time interval is increasing until to a certain maximum level less than 20 . Figure 11 shows that increasing the value of d0, increases a little the vaccination interval, while the treatment interval decreases to zero. In Figure 12, when the value of d1 is increased, it follows that the vaccination interval is increased slightly while the treatment interval decreases. Figure 13 shows that no matter how much d2 is increased, both vaccination and treatment intervals don’t change much. Figure 14 shows that by increasing the value of z, the vaccination interval increases a little while the treatment interval is reduced. The rightmost sub-figure in each set of figures shows the rates of convergence, and the costs of the solutions found by the algorithm. The algo- rithm takes 3 to 20 iterations to converge, and increasing any cost parameter produces increased final cost. Typically initial large decreases in cost which represents the cost improvement from Al- gorithm I followed by Algorithm IIa. Usually little improvement iteration of part III modify each treatment or vaccination interval by stimulating variable reduction in the cost. The algorithm some times III brings rapid improvement followed by slower. During the rapid phase, both controls are being adjusted. The rapid improvement phase ends when one control has reached its optimal configu- ration and the control continues to adjust. In most cases, execution time to compute each optimal control was a minute or less on an Intel R©Core(TM) i3-2328M CPU at 2.20GHz, 2200 MHz, 2 cores, 4 logical processors with 4G RAM. CONCLUSIONS In this paper, we introduced an SIS epidemic model under vaccination and treatment controls. We formulated an objective function and simpli- fied it for numerical computations: the simplified objective function can configured to give an upper bound, lower bound, or best estimate for the cost. Biomath 8 (2019), 1912137, http://dx.doi.org/10.11145/j.biomath.2019.12.137 Page 16 of 22 http://dx.doi.org/10.11145/j.biomath.2019.12.137 V. Mbazumutima, C. Thron, L. Todjihounde, Enumerative numerical solution for optimal control using ... Figure 8. Optimal vaccination, treatment controls and strategy cost for different values of c0. Figure 9. Optimal vaccination, treatment controls and strategy cost for different values of c1. We determined necessary conditions for optimal control treatment and necessary conditions for optimal control vaccination with simplified cost function. We established an algorithm to optimize the strategy cost. This algorithm has been im- proved to reduce the execution time to find the best strategy cost. We also verified that the final strategy obtained by the algorithm in simulation Biomath 8 (2019), 1912137, http://dx.doi.org/10.11145/j.biomath.2019.12.137 Page 17 of 22 http://dx.doi.org/10.11145/j.biomath.2019.12.137 V. Mbazumutima, C. Thron, L. Todjihounde, Enumerative numerical solution for optimal control using ... Figure 10. Optimal vaccination, treatment controls and strategy cost for different values of c2. Figure 11. Optimal vaccination, treatment controls and strategy cost for different values of d0. satisfies the local optimum conditions given in Theorem 5.1 . Although, this does not guarantee global optimality, the fact that we have tried a large diversity of strategies in the algorithm makes it plausible that the final strategy is indeed the global optimum. Finally, some numerical simulations are presented to illustrate the performance of this algorithm. Biomath 8 (2019), 1912137, http://dx.doi.org/10.11145/j.biomath.2019.12.137 Page 18 of 22 http://dx.doi.org/10.11145/j.biomath.2019.12.137 V. Mbazumutima, C. Thron, L. Todjihounde, Enumerative numerical solution for optimal control using ... Figure 12. Optimal vaccination, treatment controls and strategy cost for different values of d1. Figure 13. Optimal vaccination, treatment controls and strategy cost for different values of d2. ACKNOWLEDGEMENTS The support of this research through the Ger- man Academic Exchange Service (DAAD) and the anonymous reviewers are hereby acknowledged. Biomath 8 (2019), 1912137, http://dx.doi.org/10.11145/j.biomath.2019.12.137 Page 19 of 22 http://dx.doi.org/10.11145/j.biomath.2019.12.137 V. Mbazumutima, C. Thron, L. Todjihounde, Enumerative numerical solution for optimal control using ... Figure 14. Optimal vaccination, treatment controls and strategy cost for different values of z. REFERENCES [1] David Aadland, David C. Finnoff, and Kevin X. D. Huang. Syphilis cycles. The B.E Journal of Economic Analysis & Policy, 13(1):297–348, 2013. [2] F. S. Akinboro, S. Alao, and F. O. Akinpelu. Numerical solution of SIR model using differential transformation method and variational iteration method. General Math- ematics Notes, 22(2):82–92, 2014. [3] Oluwaseun B. Akinduko, C. Y. Ishola, O. A. Afolabi, A. B. Ganiyu, et al. Series solution of typhoid fever model using differential transform method. Malaysian Journal of Computing, 3(1):67–80, 2018. [4] Muhammad Asghar Ali, Muhammad Rafiq, and Muhammad Ozair Ahmad. Numerical analysis of a modified SIR epidemic model with the effect of time delay. Punjab Univ. j. math, 51:79–90, 2019. [5] Christian Almeder, Gustav Feichtinger, Warren C. Sanderson, and Vladimir M. Veliov. Prevention and medication of HIV/AIDS: the case of botswana. Central European Journal of Operations Research, 15(1):47–61, 2007. [6] Abraham J. Arenas, Gilberto González-Parra, and Ben- ito M. Chen-Charpentier. A nonstandard numerical scheme of predictor–corrector type for epidemic mod- els. Computers & Mathematics with Applications, 59(12):3740–3749, 2010. [7] Omar Abu Arqub and Ahmad El-Ajou. Solution of the fractional epidemic model by homotopy analysis method. Journal of King Saud University-Science, 25(1):73–81, 2013. [8] John T. Betts. Practical methods for optimal control and estimation using nonlinear programming, volume 19. Siam, 2010. [9] Hans Georg Bock, Christian Kirches, Andreas Meyer, and Andreas Potschka. Numerical solution of optimal control problems with explicit and implicit switches. Optimization Methods and Software, 33(3):450–474, 2018. [10] Riccardo Bonalli, Bruno Hérissé, and Emmanuel Trélat. Solving nonlinear optimal control problems with state and control delays by shooting methods combined with numerical continuation on the delays. arXiv preprint arXiv:1709.04383, 2017. [11] Alpha C. Chiang. Elements of dynamic optimization. McGraw-Hill, New York, 1992. [12] Yat Tin Chow, Jérôme Darbon, Stanley Osher, and Wotao Yin. Algorithm for overcoming the curse of di- mensionality for state-dependent hamilton-jacobi equa- tions. Journal of Computational Physics, 387:376–409, 2019. [13] Brian C. Fabien. Numerical solution of constrained optimal control problems with parameters. Applied Mathematics and Computation, 80(1):43–62, 1996. [14] William M Feldman and Panagiotis E Souganidis. Homogenization and non-homogenization of certain non-convex hamilton–jacobi equations. Journal de Mathématiques Pures et Appliquées, 108(5):751–782, 2017. [15] Zhilan Feng, Wenzhang Huang, and Carlos Castillo- Chavez. Global behavior of a multi-group SIS epidemic model with age structure. Journal of Differential Equa- Biomath 8 (2019), 1912137, http://dx.doi.org/10.11145/j.biomath.2019.12.137 Page 20 of 22 http://dx.doi.org/10.11145/j.biomath.2019.12.137 V. Mbazumutima, C. Thron, L. Todjihounde, Enumerative numerical solution for optimal control using ... tions, 218(2):292–324, 2005. [16] Pierre-Yves Geoffard and Tomas Philipson. Rational epidemics and their public control. International eco- nomic review, pages 603–624, 1996. [17] Mark Gersovitz and Jeffrey S. Hammer. The economical control of infectious diseases. The Economic Journal, 114(492):1–27, 2003. [18] Giorgi Giorgio and Cesare Zuccotti. Metzlerian and generalized metzlerian matrices: Some properties and economic applications. Journal of Mathematics Re- search, 7(2):42, 2015. [19] Steven M. Goldman and James Lightwood. Cost op- timization in the SIS model of infectious disease with treatment. 1996. [20] Gilberto Carlos González Parra, Rafael Jacinto Vil- lanueva Micó, and Abraham José Arenas Tawil. Matrix nonstandard numerical schemes for epidemic models. WSEAS Transactions on Mathematics, 9(11):840–850, 2010. [21] Fazal Haq, Muhammad Shahzad, Shakoor Muhammad, Hafiz Abdul Wahab, et al. Numerical analysis of fractional order epidemic model of childhood diseases. Discrete Dynamics in Nature and Society, 2017, 2017. [22] Herbert W. Hethcote. Qualitative analyses of communi- cable disease models. Mathematical Biosciences, 28(3- 4):335–356, 1976. [23] Li Jianquan and Ma Zhien. Global analysis of SIS epidemic models with variable total population size. Mathematical and computer modelling, 39(11- 12):1231–1242, 2004. [24] Morton I. Kamien and Nancy L. Schwartz. Sufficient conditions in optimal control theory. Journal of Eco- nomic Theory, 3(2):207–214, 1971. [25] T. K. Kar and Ashim Batabyal. Stability analysis and optimal control of an SIR epidemic model with vaccination. Biosystems, 104(2-3):127–135, 2011. [26] Hamid Reza Karimi. A computational method for optimal control problem of time-varying state-delayed systems by haar wavelets. International Journal of Computer Mathematics, 83(02):235–246, 2006. [27] Hassan Laarabi, El Houssine Labriji, Mostafa Rachik, and Abdelilah Kaddar. Optimal control of an epidemic model with a saturated incidence rate. Nonlinear anal- ysis: modelling and control, 17(4):448–459, 2012. [28] Karen E. Lamb, David Greenhalgh, and Chris Robert- son. A simple mathematical model for genetic effects in pneumococcal carriage and transmission. Journal of Computational and Applied Mathematics, 235(7):1812– 1818, 2011. [29] Abid Ali Lashari, Shaban Aly, Khalid Hattaf, Gul Zaman, Il Hyo Jung, and Xue-Zhi Li. Presentation of malaria epidemics using multiple optimal controls. Journal of Applied mathematics, 2012, 2012. [30] Abid Ali Lashari, Khalid Hattaf, Gul Zaman, and Xue- Zhi Li. Backward bifurcation and optimal control of a vector borne disease. Applied Mathematics & Information Sciences, 7(1):301–309, 2013. [31] Abid Ali Lashari and Gul Zaman. Optimal control of a vector borne disease with horizontal transmission. Non- linear Analysis: Real World Applications, 13(1):203– 212, 2012. [32] Kwang Sung Lee and Abid Ali Lashari. Stability analysis and optimal control of pine wilt disease with horizontal transmission in vector population. Applied Mathematics and Computation, 226:793–804, 2014. [33] Suzanne Lenhart and John T. Workman. Optimal control applied to biological models. Crc Press, 2007. [34] James Lightwood and Steven M. Goldman. The SIS model of infectious disease with treatment. Unpublished Manuscript, 1995. [35] Khosrow Maleknejad and E Saeedipoor. An efficient method based on hybrid functions for fredholm integral equation of the first kind with convergence analysis. Applied Mathematics and Computation, 304:93–102, 2017. [36] Olvi L. Mangasarian. Sufficient conditions for the optimal control of nonlinear systems. SIAM Journal on control, 4(1):139–152, 1966. [37] AG M’Kendrick. Applications of mathematics to medi- cal problems. Proceedings of the Edinburgh Mathemat- ical Society, 44:98–130, 1925. [38] Ghazala Nazir, Shaista Gul, et al. Comparative study of mathematical model of ebola virus disease via using differential transform method and variation of iteration method. Matrix Science Mathematic (MSMK), 3(1):17– 19, 2019. [39] Peter Neal. Stochastic and deterministic analysis of SIS household epidemics. Advances in applied probability, 38(4):943–968, 2006. [40] Rachael Miller Neilan and Suzanne Lenhart. An in- troduction to optimal control with an application in disease modeling. In Modeling Paradigms and Analysis of Disease Trasmission Models, pages 67–82, 2010. [41] Jiafeng Pan, Alison Gray, David Greenhalgh, and Xuerong Mao. Parameter estimation for the stochastic SIS epidemic model. Statistical Inference for Stochastic Processes, 17(1):75–98, 2014. [42] Sarada S. Panchanathan, Diana B. Petitti, and Dou- glas B. Fridsma. The development and validation of a simulation tool for health policy decision making. Jour- nal of biomedical informatics, 43(4):602–607, 2010. [43] SM Pittman, E Tannenbaum, and EJ Heller. Dynamical tunneling versus fast diffusion for a non-convex hamil- tonian. The Journal of chemical physics, 145(5):054303, 2016. [44] Lev Semenovich Pontryagin. Mathematical theory of optimal processes. Routledge, 2018. [45] Timothy C. Reluga. An SIS epidemiology game with two subpopulations. Journal of Biological Dynamics, 3(5):515–531, 2009. [46] Timothy C. Reluga. Game theory of social distancing in response to an epidemic. PLoS computational biology, 6(5):e1000793, 2010. [47] Suresh P. Sethi. Optimal quarantine programmes for Biomath 8 (2019), 1912137, http://dx.doi.org/10.11145/j.biomath.2019.12.137 Page 21 of 22 http://dx.doi.org/10.11145/j.biomath.2019.12.137 V. Mbazumutima, C. Thron, L. Todjihounde, Enumerative numerical solution for optimal control using ... controlling an epidemic spread. Journal of the Opera- tional Research Society, 29(3):265–268, 1978. [48] Vineet K. Srivastava, Mukesh K. Awasthi, and Sunil Kumar. Numerical approximation forHIV infection of CD4+T cells mathematical model. Ain Shams Engineering Journal, 5(2):625–629, 2014. [49] Oskar Von Stryk. Numerical solution of optimal control problems by direct collocation. In Optimal Control, pages 129–143. Springer, 1993. [50] Daniel Wachsmuth. Numerical solution of optimal control problems with convex control constraints. In IFIP Conference on System Modeling and Optimization, pages 319–327. Springer, 2005. [51] Gul Zaman, Yong Han Kang, and Il Hyo Jung. Optimal vaccination and treatment in the SIR epidemic model. Proc. KSIAM, 3:31–33, 2007. [52] Anwar Zeb, Madad Khan, Gul Zaman, Shaher Momani, and Vedat Suat Ertürk. Comparison of numerical meth- ods of the SEIR epidemic model of fractional order. Zeitschrift für Naturforschung A, 69(1-2):81–89, 2014. Biomath 8 (2019), 1912137, http://dx.doi.org/10.11145/j.biomath.2019.12.137 Page 22 of 22 http://dx.doi.org/10.11145/j.biomath.2019.12.137 Introduction Some basic results on optimal control problems SIS model under treatment and vaccination controls Objective functions formulation Original objective function Piecewise linear objective function Local optimality conditions for the non-autonomous vaccination and treatment control problem Necessary conditions for the optimal control vaccination with simplified cost function Necessary conditions for the optimal control treatment Summary of necessary conditions for optimal controls Efficient numerical method for finding near-optimal controls Numerical simulations and discussion Numerical simulations Discussion References References