Feedback control problem of an SIR epidemic model based on the Hamilton-Jacobi-Bellman equation

1.   Introduction

Mathematical models are very useful for investigating the epidemiological characteristics of infectious diseases and for evaluating the effectiveness of control strategies. A susceptible-infective-recovered (SIR) model is the simplest compartment model of an epidemic proposed by O. Kermack and A. G. McKendrick in 1927 ^[1]. The model consists of three compartments: Susceptible (S), Infective (I), and Recovered (R). Previous studies have conducted several outstanding surveys of basic compartment models and have explored several features of modified models ^[2,3,4]. A qualitative analysis of a time-varying SIS epidemic model with incidence rate depending on the susceptible and infective populations is conducted in ^[5]. H. M. Yang and A. R. R. Freitas explore a mathematical model, taking into account the updated biological aspects of rubella vaccine and application to the recently introduced dengue vaccine ^[6].

In general, vaccination is an important control measure to reduce the size of infected populations. Many researchers have studied how to predict and evaluate the effectiveness of various vaccination strategies with a mathematical model. M. De la Sen et al. study a simple continuous-time linear vaccination-based control strategy for an SEIR epidemic model ^[7]. Using the method of Lyapunov functions, global stability of the disease-free equilibria of a two-group SI model with feedback control variables is investigated in ^[8]. In addition, optimal control techniques have been used in epidemic models to derive an efficient vaccination policy for influenza outbreaks in specific settings ^{[9,10,11,12,13,14]}. An optimal control problem has been studied to minimize the costs associated with treatment and preventive campaigns to avoid relapse in ^[9]. J. Lee et al. consider optimal control theory to derive optimal intervention strategies for infectious diseases, using vaccinations, antiviral therapy, and social distance as control functions ^[10]. Some authors discuss the optimal quarantine control issues with an SIQS model describing the epidemic dynamical process on complex networks ^[13]. An optimal vaccine distribution strategy under limited vaccine resources using an epidemic model with group mixing has been explored in ^[14]. Moreover, many researchers have also used optimal control theory to develope optimal treatment strategies for other diseases, such as HIV, tuberculosis, and vector-borne diseases ^{[15,16,17,18]}.

Two main control techniques lead to optimal vaccination strategies: open-loop control and closed-loop control, also called feedback control. Most of the above-mentioned optimal control studies are the result of open-loop control. Open-loop control can be found by solving the optimality system which is obtained from Pontryagin's principle ^[19]. However, open-loop control suffers the disadvantage of only depending on the initial value of the state and not being able to reflect changes in the state that may appear over time. Consequently, the results are generally not useful for practical applications. In the present paper, we study an optimal feedback control problem of an SIR epidemic model with the goal of inducing optimal vaccination strategies that are dependent on the states, i.e., that can immediately reflect any perturbation of the states.

Various methodologies have been proposed to solve the feedback control problem. When the model is a linear system and the objective functional is quadratic in the state and control, feedback control can easily be obtained by solving the Riccati equation ^[20,21]. However, most mathematical models for biological systems including the SIR model as studied in this paper are nonlinear systems. The state-dependent Riccati equation (SDRE) is a powerful approach to solve the feedback control problems with nonlinear dynamical systems ^[21,22]. The main drawback of the SDRE approach is that the proposed feedback controller is only suboptimal ^[23]. Another extension is proposed in ^[24] to design an LQR controller such that the properties guaranteed to the norm-bounded linear differential inclusions (NLDI) representation are also valid for the underlying nonlinear system. A state feedback control scheme is applied to a nonlinear HIV model to reduce the spread of infection ^[25]. The control design is based on the linearization in the neighborhood of a suitable equilibrium point, aiming at applying the linear quadratic (LQ) regulator theory.

In this paper, we consider a special nonlinear partial differential equation called the Hamilton-Jacobi-Bellman (HJB) equation. The HJB equation originates from the idea of dynamic programming theory proposed by Richard Bellman ^[26]. The HJB approach can be generalized to the case where the constraint equation is a stochastic system. Knowing the value function, which is the solution of the HJB equation, we can obtain the optimal feedback control function. However, there is generally no guarantee of the existence and uniqueness of the solution regardless of the smoothness of the coefficients of the HJB equation. To overcome these shortcomings, M. G. Crandall and P. L. Lions developed the concept of viscosity solutions in 1983 ^[27]. The basic idea of a viscosity solution theory is to add a small artificial viscosity term to the HJB equation. The existence and uniqueness of the viscosity solutions have been established and the convergence of the viscosity solutions to the solution of the HJB equation has been proven in ^[28,29]. Since it is almost impossible to find a viscosity solution analytically, a numerical algorithm is needed to approximate the viscosity solution corresponding to the feedback controller ^{[30,31,32,33,34]}. The HJB equation is a system combined with nonlinear PDE and optimization problems. We employ a successive approximation method to decouple the system and use the upwind finite difference method to solve a convection-diffusion equation that is a part of the decoupled system.

The rest of this paper is organized as follows. In section 2, we formulate an optimal feedback control problem for a nonlinear SIR model to derive an efficient vaccination strategy for influenza outbreaks. To obtain a feedback control function, we first define a value function, and then derive the HJB equation with the value function, according to Bellman's principle. A viscosity solution that is an approximate solution to the solution of the HJB equation is also discussed. In section 3, we introduce a successive approximation method combined with the upwind finite difference method to find the viscosity solution. In section 4, we present numerical results of the optimal vaccination strategy which can immediately reflect changes in the number of susceptible and infective individuals in various situations. We provide concluding remarks in section 5.

2.   Feedback control problem and HJB

In this section we formulate a feedback control problem of an SIR model and derive the HJB equation to find the optimal feedback control. The SIR model is described by the following ordinary differential equation (ODE) system:

Here $S(t)$, $I(t)$, and $R(t)$ denote the density of susceptible, infective, and recovered individuals at time $t$, respectively. The control function $u(t)$ indicates the rate at which susceptible individuals are vaccinated. The constants $\beta$ and $\mu$ are the rate of disease transmission and vaccination efficacy, respectively. Infective individuals leave the compartment at rate $\gamma$. We may consider only equations of the states $S$ and $I$, since the state $R$ is decoupled with respect to the system of the states $S$ and $I$. The control variable $u(t)$ is defined as the rate of susceptible individuals receiving the vaccination, so the value will vary within the interval $[0,1]$. Let $\mathcal{U}$ be the set of admissible control functions,

where $s$ is an initial time and $T$ is a terminal time.

When developing a response plan for disease outbreaks, policymakers look for strategies to minimize the disease incidence, while minimizing the cost of intervention strategies. So our goal in the optimal control problem is to minimize the number of infected people while minimizing vaccine immunization costs. We now formulate the control problem in the form

subject to the model

where $y = (y_1, y_2)$ is an initial condition. The functional $\mathcal{J}(s, y, u)$ is called the objective or cost functional. The second and third terms in this functional represent the number of administered vaccines and the vaccination cost, respectively. The last term in the objective functional is also added to facilitate the mathematical analysis of control theory. $w_1, w_2,$ and $w_3$ are weight constants for balancing the relative costs in the objective functional. When the initial time $s = 0$ and the initial state $y$ are fixed, the above control problem is called an open-loop optimal control problem. The solution in the open-loop control problem is optimal only along the optimal trajectory starting from the initial state $y$. In this paper, we assume that $s$ and $y$ in (2.2) and (2.3) are not fixed, so that the solution of our control problem provides a global optimal control defined over an entire time-space region. That is, even when the state deviates from the optimal trajectory due to disturbance, an optimal control function corresponding to the system can be found immediately.

In general, finding the optimal feedback control is difficult, but Bellman's dynamic programming approach can provide a useful starting point for finding the control function ^[26]. The optimal feedback control can be given by solving the HJB equation which consists of the value function of the optimal control. To achieve this, we consider the abstract form of the optimal control problem and introduce the value function:

subject to

where $x(t) = (S(t), I(t))^T$, the integrand $f(t, x(t), u(t)) = w_{1}I^2(t)+w_{2} u(t) S(t) + w_{3}u^2(t)$, and the right-hand side function

Define the value function

The optimal control problem (2.4) and (2.5) can be reformulated into the HJB equation by Bellman's principle ^[35]. By the dynamic programming approach, we get

Taking the limit as $h$ goes to $0$, we arrive at

Finally, the HJB equation is derived as follows:

where $\nabla$ is the gradient with respect to $x$. In our case, the dimension of space is two, i.e., $n = 2$. For more details about deriving the HJB equation, see the references ^[29,35].

The relationship between the value function and the optimal feedback control is well explained in the following two propositions and theorem ^[29,30,31]. In particular, Theorem 2.3 provides a method for constructing the optimal feedback control function numerically.

Proposition 2.1. $\phi\in C^1([t_0,t_f] \times \mathbb{R}^n)$ be the value function. $u^*(\cdot)$ is an optimal control where $x^*$ is the state corresponding to $u^*$ if there exists a control $u^*(t)\in \mathcal{U}$ such that

The optimal control is denoted as

for almost all $t\in[t_0,t_f]$.

Proposition 2.2. Let $\phi(t,x)\in C^1([t_0,t_f]\times \mathbb{R}^n)$ be the value function. Then $(u^*(t), x^*(t))$ is an optimal pair in feedback form if and only if

Theorem 2.3. Let $\phi(t,x)\in C^1([t_0,t_f]\times \mathbb{R}^n)$ be the value function. Assume $u(t,x)$ satisfies

Then,

is the feedback table of the optimal control problem where $x_{y}(t)$ is the solution of

In general, the existence and uniqueness of the solution of the HJB Eq (2.6) are not guaranteed regardless of the smoothness of the functions $g$ and $f$. These drawbacks can be overcome by applying the idea of the viscosity solution, as introduced by M. G. Crandall and P. L. Lions in the early 1980s ^[27]. Consider the following convection-diffusion equation:

This equation differs from (2.6) by the addition of the diffusion term $\epsilon \nabla^2 \phi(t, x)$, where $\epsilon >0$ represents the diffusion parameter or viscosity. When $\epsilon$ is sufficiently small, the Eq (2.7) is a viscosity approximation to (2.6) and its solution is called a viscosity solution. The existence and uniqueness of the viscosity solution have been studied by many researchers ^[28,29]. It was proven that the solution of (2.7) converges to the solution of (2.6) as $\epsilon \rightarrow 0$ ^[29,34]. Nevertheless, finding an analytical solution of (2.7) remains impossible.

Although the Eq (2.7) is defined on the unbounded domain $\mathbb{R}^n$, we often limit our consideration to a bounded domain $\Omega \subset \mathbb{R}^n$ as follows:

The feedback control function can be derived by solving the Eq (2.8) with the given Dirichlet boundary condition in the bounded domain numerically. One way to get the boundary condition is to interpolate the solutions of a sequence of open-loop control problems at some points on the boundary of domain, $\partial\Omega$. However, this approach necessitates very high computational cost. Alternatively, we introduce the heuristic approach to impose the Dirichlet boundary condition. We extend $\Omega$ to a larger domain $\tilde{\Omega}$ satisfying $\Omega \subset \tilde{\Omega}$ and reformulate the Eq (2.8) with the larger domain $\tilde{\Omega}$.

where $\tilde{h}(t, x)$ is an artificial Dirichlet boundary condition. A natural choice is that $\tilde{h}(t, x) = \phi(T, x)$. For details about this heuristic approach to impose the boundary condition, see the references ^[33,34,36]. While the approach for finding an approximate optimal control is heuristic, comparing feedback control with the open-loop control in Figure 4 provides means to verify that the approach is effective. In the absence of disturbances in the state, the feedback control is almost identical to the open-loop control obtained by applying the Pontryagin’s Maximum Principle.

3.   Successive approximation method

We may determine an approximate solution to the optimal feedback control problem (2.4) and (2.5) by solving the HJB Eq (2.9). Since the HJB Eq (2.9) is the coupled system with $\phi(t,x)$ and $u(t,x)$, the equation can be decoupled by the following algorithm:

This approach is known as the successive approximation method ^[37]. After decoupling, the HJB Eq (2.9) consists of a convection-diffusion Eq (3.1) and a simple optimization problem (3.2). Among many approaches for solving the convection-diffusion Eq (3.1), we use the upwind finite difference method.

Let $\varDelta t = 1/N_t$, where $N_t$ is the number of subintervals in the time interval and $t_k = k \varDelta t$ for $k = 0,1,\cdots, N_t$. For spatial discretization, let $\varDelta x_1$ and $\varDelta x_2$ be step sizes of $\Omega \subset \mathbb{R}^2$ for each dimension. Define $x_1^{(n)} = n\varDelta x_1$ and $x_2^{(m)} = m\varDelta x_2$ for $n = 1,2,\cdots, N_{x_1}$, $m = 1,2,\cdots, N_{x_2}$, where $N_{x_2}$ is the number of discretization in each dimension. We also define

The upwind finite difference method for solving the convection-diffusion Eq (3.1) is as follows: for the given $u(t,x)$,

for $n = 1, 2, 3, \cdots, N_{x_1}-1$, $m = 1, 2, 3, \cdots, N_{x_2}-1$, and $k = 0, 1,2,\cdots, N_t$.

The successive approximation algorithm for solving the HJB equation is given in Algorithm 1:

4.   Numerical simulations

In the numerical simulations, we assume that the SIR model is defined over a time interval $[0, \ 30]$ and on a spatial domain $\Omega = \left\{ {x = (S,I)\in \mathbb{R}^2 : 0 \leq S, \, \ 0 \leq I, \, \ 0 \leq S + I \leq N} \right\}$. Here, $N$ is the total population, i.e., $N = S+I+R$. It is also assumed that the total population, $N$, is a constant set at 100. Consequently, the spacial domain $\Omega$ has a triangular shape. As we discussed in the previous section, we need to extend the domain $\Omega$ to $\tilde{\Omega}$ for the numerical simulation. We set $\tilde{\Omega} = [-2N,2N]\times[-2N,2N]$, and then $\Omega \subset \tilde{\Omega}$. In addition, we assign $\tilde{h}(t,x) = 0$ as a Dirichlet boundary condition on $\partial \tilde{\Omega}$. To balance the differences in the magnitude of the number of infected individuals, the number of vaccinated people, and control function in the objective functional (2.2), the weight constants $w_1$, $w_2$, and $w_3$ are selected to be $w_1 = 1$, $w_2 = 10$, and $w_3 = 100$, respectively. The parameters and the number subintervals used in solving the HJB equation are chosen as follows:

Figure 1 displays the value of the feedback control function $u(t,S,I)$ on the whole domain $\Omega$ at $t = 0.$ The feedback control function provides complete policy information for all combinations of the state values, $S$ and $I$, and may therefore be used as a complete map for planning because the value of the control function can be obtained immediately by knowing the values of $S$ and $I$. It is also helpful in decision making as it may compensate for the uncertainty in the states by referring to the results of the neighbors. The numerical result shows that vaccination should be done at the highest rate for effective control of the epidemic in most cases, for example, $(S, I) = (80, 10)$, $(60, 30)$, $(40, 40)$, $(20, 70)$, and so on. However, no vaccination should be done when the number of susceptible or infective individuals is low, for example, $(S, I) = (20, 5)$, $(10, 10)$, $(5, 15)$, and so on. To explain this in more detail, let's consider a case in which the value of $I$ is fixed at 10. In this case, complete vaccination is needed if S is greater than 40 but no vaccine is needed if S is less than 25 because many recovered individuals ($R$) already have immunity. For instance, $R = 70$ when $S = 20$ and $I = 10$ among the total population $N = 100$. Therefore, the number of individuals with immunity can be a good index to determine the level of vaccination needed.

To examine how the optimal feedback control changes over time, the functions for different initial time values $s = 0, 7, 14, 21, 28, 29$ are presented in Figure 2. We also give the cost functional values (2.2) for different initial time values in Table 1. The numerical results show that the shape of the feedback control functions is almost identical and the cost functional values are almost same for the first three weeks, after which the domain that requires vaccination shrinks gradually. For instance, vaccination should be done at the highest rate for the first three weeks, but vaccines are not needed at time $t = 29$ when $S = 80, I = 10$. In most domains, vaccination is no longer needed at time t = 29, but is still required when the number of infected individuals is very high.

In order to investigate the change in the shape of the control function $u(t,S,I)$ with respect to time $t$ more precisely, we provide $L_2$ norm of $u_t$ in Figure 3. The value of $\| u_t \|_2$ is almost zero until the 24th day, which means that there is no meaningful change in the shape of the control function until the 24th day. This implies that the optimal vaccination strategy is independent of the time variable except for the last few days and is determined only by the $S$ and $I$ values.

In Figure 4, we present the optimal vaccination strategies and the corresponding states with given initial time $s = 0$ and initial states $(S_0, I_0) = (95, 5)$ using both open-loop control and feedback control approaches for comparison. In the open-loop control approach, we apply the Pontryagin's Maximum Principle to induce an optimality system from which the optimal solution is obtained. A gradient-based algorithm is used to solve the optimality system numerically. For the details of this algorithm, we refer the interested reader to ^[38,39]. The first row is the result of using the open-loop control approach and the second row is the result of using the feedback control approach: the control functions and the corresponding states are almost the same. The graphs indicate that for the effective control of the endemic, it is essential to vaccinate at the highest rate for the first few days and then gradually lower the rate. The small difference in the two approaches is the cost functional value, which is 1510.4473 (open-loop) and 1508.5036 (feedback).

Now we consider the case where a small change in the value of the state variable occurs due to environmental factors or measurement errors. We add a randomly selected integer in $[-10, 10]$ to the $S(t)$ and a randomly selected integer in $[-5, 5]$ to the $I(t)$ at time $t = 1, 2, 3, 7,$ and $9$. The black square dots indicate the number of the perturbed susceptible and infective individuals in Figure 5. As shown in Figure 5, the control function obtained using the open-loop approach does not respond to changes in the states at all. On the other hand, the feedback control function properly reflects the change of the states. As a result, vaccination should be done at the highest rate for a longer period of time in feedback control. In the case of open-loop control, the numerical results also show that vaccination should be completely stopped after about three days because the peak of the epidemic has passed, whereas the feedback control function indicates that a certain level of vaccination should be used after three days due to a change in the state variables.

The impact of the weight constants in the cost functional (2.2) is explored by varying the weights $w_2$ and $w_3$. Figure 6 shows the feedback control functions $u(t,S,I)$ at $t = 0$ with the choice of $w_2 = 10, 30, 50,$ and $70$. As shown in Figure 6, the areas of $S$ and $I$ that need to use vaccination at the highest rate decrease as the weight constant $w_2$ increases. Note that the transition area, where the value of the control function changes from 0 to 1, is very narrow and almost takes the form of a straight line ($S+I = k$). This means that the sum of $S$ and $I$ act as a threshold quantity in the sense that $u = 1$ if it is greater than the specified value ($k$) and $u = 0$ otherwise. So policymakers can decide whether to give vaccines by taking into account the number of recovered individuals ($R$) who are immune since $R = N-(S+I)$. For example, in the case of $w_2 = 70$, vaccination is unnecessary if there are more than twenty people who already have immunity, but vaccination should be applied if the number is less than twenty.

Figure 7 shows the feedback control functions for different values of $w_3 = 100,300,500,$ and $700$. As expected, the areas in which the value of the control function is set to one are decreasing as the weight constant increases, but their patterns are different from those in the previous case. The transition area is no longer a straight line and the area widens. The value of the control function changes gradually, rather than suddenly switching from 1 to 0, depending on the values of $S$ and $I$. In practice, this suggests that if the cost of the vaccine is so high that its use should be reduced, then vaccination should be gradually reduced but not completely stopped. Even in these cases, the feedback control functions, $u(t, S, I)$, are not affected by time variables except for the last few days. Most of the time, the trend of the optimal vaccination strategies depending on the state values is unchanged.

Finally, we study the effects of the disease transmission rate, $\beta$, on optimal vaccination strategies with the choice of $\beta = 0.005, 0.01, 0.015,$ and $0.02$. Figure 8 demonstrates that the areas of $S$ and $I$ that do not require the use of vaccines widen as the rate decreases. It is noteworthy that the area is biased downwards. That is, when the disease transmission rate is significantly low, vaccination is unnecessary, no matter how many people are susceptible, unless there are many infective individuals.

5.   Conclusion

This paper formulates a feedback control problem of an SIR model to develop optimal vaccination strategies depending on the state values. The HJB equation is derived to solve the control problem using the idea of Bellman's dynamic programming approach. The relationship between the optimal feedback control and the HJB equation is discussed. We introduce a successive approximation method combined with the upwind finite difference method to find the feedback control functions.

The numerical results show that the feedback control functions provide complete policy information for all combinations of the state values, $S$ and $I$. In addition, the optimal vaccination strategies are independent of the time variable except for the last few days and are determined only by the state values. We present the optimal vaccination strategies and the corresponding states using both open-loop control and feedback control approaches to compare them. In both approaches, it is very effective to vaccinate at the highest rate for the first few days and then gradually lower the rate. However, the feedback control function properly and immediately reflects the change in state, but the open-loop control approach does not. We explore the impact of the weight constants in the cost functional on optimal vaccination strategies by varying the weights. The areas of the states that require vaccination at the highest rate decrease as the weight constants increase. Numerical simulations also indicate that it is not necessary to administer the vaccine, no matter how many people are susceptible, as long as the disease transmission rate is significantly low and the number of infective individuals is not high.

Acknowledgements

The work of Hee-Dae Kwon was supported by NRF-2016R1D1A1B04931897 and the NST grant by the Korea government (MSIP) (No. CRC-16-01-KRICT). The work of Jeehyun Lee was supported by NRF-2015R1A5A1009350 and NRF-2016R1A2B4014178.

Conflict of Interest

The authors declare that there is no conflict of interest.