# Feedback control of isolation and contact for SIQR epidemic model via strict Lyapunov function

Hiroshi Ito, Michael Malisoff, Frederic Mazenc

► To cite this version:

Hiroshi Ito, Michael Malisoff, Frederic Mazenc. Feedback control of isolation and contact for SIQR epidemic model via strict Lyapunov function. Mathematical Control and Related Fields, 2023, 13 (4), pp.1438 - 1465. 10.3934/mcrf.2022043. hal-04337733

HAL Id: hal-04337733

https://inria.hal.science/hal-04337733v1

Submitted on 12 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
FEEDBACK CONTROL OF ISOLATION AND
CONTACT FOR SIQR EPIDEMIC MODEL VIA
STRICT LYAPUNOV FUNCTION

HIROSHI ITO<sup>1</sup>, MICHAEL MALISOFF<sup>2</sup> AND FRÉDÉRIC MAZENC<sup>3</sup>

<sup>1</sup>Department of Intelligent and Control Systems
Kyushu Institute of Technology
680-4 Kawazu, Iizuka 820-8502, Japan

<sup>2</sup>Louisiana State University, Department of Mathematics
Baton Rouge, LA 70803-4918, USA

<sup>3</sup>Inria EPI DISCO, L2S-CNRS-CentraleSupélec
3 rue Joliot Curie 91192, Gif-sur-Yvette, France

(Communicated by Jean-Michel Coron)

ABSTRACT. We derive feedback control laws for isolation, contact regulation, and vaccination for infectious diseases, using a strict Lyapunov function. We use an SIQR epidemic model describing transmission, isolation via quarantine, and vaccination for diseases to which immunity is long-lasting. Assuming that mass vaccination is not available to completely eliminate the disease in a time horizon of interest, we provide feedback control laws that drive the disease to a small endemic equilibrium. We prove the input-to-state stability (or ISS) robustness property on the entire state space, when the immigration perturbation is viewed as the uncertainty. We use an ISS Lyapunov function to derive the feedback control laws. A key ingredient in our analysis is that all compartment variables are present not only in the Lyapunov function, but also in a negative definite upper bound on its time derivative. We illustrate the efficacy of our method through simulations, and we discuss the usefulness of parameters in the controls. Since the control laws are feedback, their values are updated based on data acquired in real time. We also discuss the degradation caused by the delayed data acquisition occurring in practical implementations, and we derive bounds on the delays under which the ISS property is maintained when delays are present.

1. Introduction. We establish a basic framework for global feedback control con-
struction for epidemics exhibiting long-lasting immunity. We use the SIQR model
consisting of four compartments, namely, susceptible individuals, infected individu-
als, quarantined individuals, and recovered individuals [8, 14]. Our work is amenable

<footnote n="2020">2020 Mathematics Subject Classification. Primary: 93D30, 93C10, 93D09; Secondary: 92D25, 34D23.</footnote>

<footnote>*Key words and phrases.* Epidemic models, quarantine, contact regulations, vaccination, Lyapunov functions, robustness, input delays.</footnote>

<footnote>The work of H. Ito was supported by JSPS KAKENHI Grant Number JP20K04536. The work of M. Malisoff was supported by NSF Grant 1711299.</footnote>

<footnote>*Corresponding author: Hiroshi Ito.</footnote>
to cases where the disease cannot be eliminated due to insufficient availability of vac-
cines counteracting strong infectability. Popular approaches to tackling disease con-
trol problems entail Jacobian linearization. However, linearization is only valid in a
sufficiently small region. To prove global results, a common technique is to instead
use Lyapunov functions in conjunction with LaSalle’s invariance principle, but this
technique is not amenable to proving robustness properties [9, 10, 11]. As for sta-
bility analysis without uncertainty, this difficulty has sometimes been bypassed by
removing variables and introducing state space restrictions [8, 17, 18, 19, 37, 38, 40],
which leads to analyzes that assume a constant population or that use Lyapunov
functions in a restricted domain of a certain phase. While previous approaches
have been accepted in the mathematical biology community, they have hampered
the development of global robust feedback control synthesis, which call for strict
Lyapunov functions (as defined, e.g., in [23]).

Recently, a significant advance was presented to prove input-to-state stability (or ISS) for the SIQR model on the entire state space [12]. This entailed constructing a kind of strict Lyapunov function, called an ISS Lyapunov function. The construction in [12] used a nontrivial analog of the strictification technique from [23]. Unlike Lyapunov functions that are popular in epidemiology, all of the compartment variables appeared in not only the proposed Lyapunov function, but also in a negative definite upper bound on its time derivative in the Lyapunov function decay estimate. In the epidemic case, the Lyapunov function was used to propose a control law for mass vaccination. However, [12] treated the rate $\beta$ of transmission and contact [15], and the rate $\nu$ at which infected individuals are isolated via quarantine, as two constants. Hence, only the vaccination rate $\rho$ was used as a feedback control in [12]. On the other hand, endemic equilibria can naturally arise only when vaccine supply and logistics are insufficient. Therefore, here we treat all three of these rates $\beta$, $\nu$, and $\rho$ as feedback controls. Then, using a nontrivial variant of the global strict Lyapunov function analysis from [12], we design these feedbacks to globally asymptotically stabilize a class of endemic equilibria when no immigration perturbations are present, and to ensure ISS with respect to suitably bounded perturbations.

There are several ways to reduce the spread of infection [14]. In this paper, we consider isolation via quarantine and contact reduction, using the newly discovered ISS Lyapunov function to design global feedback controls that ensure robustness to a perturbation in the number of immigrants. In practice, implementation of countermeasures to infectious diseases is subject to time delay since feedback control updates input quantities based on the latest information; see e.g., [28, 29] for background on delay systems. Hence, we also demonstrate that the proposed feedback control laws are robust to delays in updating the control input quantities.

We use the following definitions and notation. We use $|f|_J$ to denote the usual sup norm of a bounded function $f$ over a subset $J$ of its domain, $|f|_\infty$ is the sup norm over its entire domain, and $|\cdot|$ is the usual Euclidean norm. Let $\mathcal{K}$ denote the set of all strictly increasing continuous functions $\alpha : [0, \infty) \to [0, \infty)$ such that $\alpha(0) = 0$; if, in addition, $\alpha$ is unbounded, then we say that $\alpha$ is of class $\mathcal{K}_\infty$. We say that a continuous function $\Phi : [0, \infty) \times [0, \infty) \to [0, \infty)$ is of class $\mathcal{KL}$ provided for each fixed $s > 0$, the function $\Phi(\cdot, s)$ belongs to class $\mathcal{K}$, and for each fixed $r \ge 0$, the function $\Phi(r, \cdot)$ is non-increasing and $\Phi(r, s) \to 0$ as $s \to \infty$. A system of the form $\dot{x}(t) = F(x(t), \varepsilon(t))$ with a state space (i.e., the set of all feasible states) $\mathcal{X} \subseteq \mathbb{R}^n$ satisfying $F(0, 0) = 0$ is called input-to-state stable (which is also abbreviated as
ISS) [16] on $\mathcal{X}$ with respect to a disturbance set $\mathcal{E}$ provided: There are $\Phi \in \mathcal{KL}$ and $\Gamma \in \mathcal{K}_{\infty}$ such that for each initial state $x(0) \in \mathcal{X}$ and each locally bounded piecewise continuous function $\varepsilon$ that is valued in $\mathcal{E}$, the corresponding solution $x: [0, \infty) \to \mathcal{X}$ of the dynamics is uniquely determined and satisfies $|x(t)| \le \Phi(|x(0)|, t) + \Gamma(|\varepsilon|_{[0,t]})$ for all $t \ge 0$.

## 2. SIQR model with vaccination and time-varying rates.
This paper considers

$$ \dot{S}(t) = B + \epsilon(t) - \rho(t)S(t) - \mu S(t) - \beta(t)I(t)S(t), \quad (1a) $$

$$ \dot{I}(t) = \beta(t)S(t)I(t) - (\gamma + \nu(t) + \mu)I(t), \quad (1b) $$

$$ \dot{Q}(t) = \nu(t)I(t) - (\tau + \mu)Q(t), \quad (1c) $$

$$ \dot{R}(t) = \gamma I(t) + \tau Q(t) - \mu R(t) + \rho(t)S(t), \quad (1d) $$

whose positive real valued state variables $S$, $I$, $Q$, and $R$ are numbers of susceptible, infected, quarantined, and recovered individuals, respectively [8, 14, 26]. The known constant $B > 0$ is the immigration rate including newborn. The unknown piecewise continuous locally bounded function $\epsilon$ represents the immigration perturbation. We assume that $\epsilon$ is valued in the perturbation set

$$ \mathcal{P} = (-B, \infty). \quad (2) $$

The non-negative-valued function $\beta(t)$ combines the rate of disease transmission and contact [15]. The positive parameters $\mu$ and $\gamma$ are the non-associated mortality rate and the rate of recovering from the disease, respectively. Equation (1d) together with (1a) imply that the SIQR model (1) assumes the acquired immunity to be long-lasting. The non-negative-valued functions $\rho(t)$ and $\nu(t)$ are the vaccination rate and the rate at which infected individuals are isolated via quarantine, respectively. The positive parameter $\tau$ is the reciprocal of the average time spent in isolation. The parameter $\rho(t)$ not only describes the amount of vaccine administration in the society, but also the reciprocal of the average time to acquire immunity. Due to the non-negativity of the parameters, the non-negative orthant $[0, \infty)^4$ is a forwardly invariant set for (1). In other words, each state component of (1) stays nonnegative for all $t \ge 0$ if the initial state for (1) is in $[0, \infty)^4$.

The goal of our control design is to derive a class of control functions that achieve ISS of an error system associated with (1) for the disturbance set $\mathcal{P}$, where the error system will be specified in the next section. The states of the error system are the differences between the states of (1) and their desired equilibrium values, after a change of variables, so this will imply convergence of the states to their desired values with an overshoot term depending on the immigration perturbation in a sense that is made precise in our ISS definition from Section 1. To specify the desired equilibrium values, our error system, and the class of controls, we will use the control inputs $u_V(t)$, $u_I(t)$, and $u_C(t)$ in

$$ \rho(t) = \hat{\rho} + u_V(t), \quad (3a) $$

$$ \nu(t) = \hat{\nu} + u_I(t), \quad \text{and} \quad (3b) $$

$$ \beta(t) = \hat{\beta} + u_C(t), \quad (3c) $$

for positive constants $\hat{\rho}$, $\hat{\nu}$, and $\hat{\beta}$ representing nominal rates, where $u_V$, $u_I$, and $u_C$ will be specified in our theorem in such a way that they will satisfy

$$ u_V(t) \in [-\hat{\rho}, \infty), \quad u_I(t) \in [-\hat{\nu}, \infty), \quad \text{and} \quad u_C(t) \in [-\hat{\beta}, 0] \quad (4) $$
for all $t \ge 0$. The nominal rate $\hat{\beta}$ describes the natural transmission rate of the disease without any contact regulation. Let

$$ \lambda = \gamma + \hat{\nu} + \mu, \quad \chi = \hat{\rho} + \mu, \quad \text{and} \quad \hat{R}_0 = \frac{\hat{\beta}B}{(\hat{\rho} + \mu)(\gamma + \hat{\nu} + \mu)}. \qquad (5) $$

The constant $\hat{R}_0$ is called the basic reproduction number. When $\epsilon = 0$ and $u_V = u_I = u_C = 0$, we can easily check that the dynamics (1) admits the equilibrium

$$ (S_{\star}, I_{\star}, Q_{\star}, R_{\star}) = \left( \frac{\lambda}{\hat{\beta}}, \frac{B}{\lambda} - \frac{\chi}{\hat{\beta}}, \frac{\hat{\nu}}{\tau + \mu} \left( \frac{B}{\lambda} - \frac{\chi}{\hat{\beta}} \right), \frac{1}{\mu} \left[ \left( \gamma + \frac{\tau \hat{\nu}}{\tau + \mu} \right) \left( \frac{B}{\lambda} - \frac{\chi}{\hat{\beta}} \right) + \frac{\hat{\rho} \lambda}{\hat{\beta}} \right] \right) \quad (6) $$

if and only if $B/\lambda \ge \chi/\hat{\beta}$, i.e., if and only if

$$ \hat{R}_0 \ge 1. \qquad (7) $$

The point (6) is called the endemic equilibrium [14]. This paper uses the notation $X_{\star} = (S_{\star}, I_{\star}, Q_{\star}, R_{\star})$. Notice that the system (1) also admits the equilibrium point

$$ \left( \frac{\lambda}{\hat{\beta}}, 0, 0, \frac{\hat{\rho}\lambda}{\mu\hat{\beta}} \right), \qquad (8) $$

which is called the disease-free equilibrium in the literature. The zeros in the second and the third components indicate that the disease is eliminated at the disease-free equilibrium. The two equilibria coincide with each other if $\hat{R}_0 = 1$. In this paper, we consider a fixed choice of the equilibrium (6), and we assume that

$$ \hat{R}_0 > 1 \qquad (9) $$

so that the nominal vaccination rate $\hat{\rho}$ is not large enough with respect to the natural disease transmission rate $\hat{\beta}$ to eliminate the disease. However, we do not assume any proportionality conditions (related, e.g., to the number $S(t)$ of susceptible individuals at each time $t$) on $\hat{\rho}$. Note that for $\epsilon = 0$, $\rho = \hat{\rho}$, $\nu = \hat{\nu}$, and $\beta = \hat{\beta}$, under the assumption $\hat{R}_0 > 1$, the disease-free equilibrium cannot be contained in the domain of attraction of the endemic equilibrium. Therefore, the domain of the initial state $(S(0), I(0), Q(0), R(0))$ is set to

$$ \mathcal{D} = (0, \infty)^4 \qquad (10) $$

for (1). Irrespective of the control inputs $u_V$, $u_I$, and $u_C$, the set $\mathcal{D}$ is forwardly invariant for (1). ISS will be established for the state space $\mathcal{X} = \mathcal{D}$. We summarize the notation in Table I below.

### 3. Control laws and guarantees.
The work [12] provided a vaccination controller assuming vaccination was the only available control, and the vaccination control in [12] was chosen to provide the required decay condition for the novel ISS Lyapunov function construction from [12]. By contrast, here we address the more complex problem that we described in the introduction above based on a novel use of a balancing parameter in the Lyapunov function construction from [12] (which we describe in further detail in Remark 3.1 below after we introduce the required notation) and three controls (instead of only one) for vaccination, isolation, and transmission and contact, which will also be chosen to provide the required Lyapunov function decay conditions along solutions of the feedback controlled dynamics.
<table><caption>TABLE 1. Parameters, Functions, and Sets from Sections 2-3</caption><thead><tr><th>Symbols</th><th>Meanings</th></tr></thead><tbody><tr><td>S(t)</td><td>number of susceptible individuals in (1)</td></tr><tr><td>I(t)</td><td>number of infected individuals in (1)</td></tr><tr><td>Q(t)</td><td>number of individuals isolated after infection in (1)</td></tr><tr><td>R(t)</td><td>number of recovered individuals in (1)</td></tr><tr><td>B</td><td>immigration rate in (1)</td></tr><tr><td>ε(t)</td><td>immigration perturbation in (1)</td></tr><tr><td>P</td><td>perturbation set (-B, ∞) from (2)</td></tr><tr><td>β(t)</td><td>transmission and contact rate in (1)</td></tr><tr><td>μ</td><td>nonassociative mortality in (1)</td></tr><tr><td>γ</td><td>recovery rate in (1)</td></tr><tr><td>ρ(t)</td><td>vaccination rate in (1)</td></tr><tr><td>ν(t)</td><td>rate of isolation in (1)</td></tr><tr><td>τ</td><td>reciprocal of average isolation time in (1)</td></tr><tr><td>ρ̂</td><td>nominal vaccination rate in (3)</td></tr><tr><td>u<sub>V</sub>(t)</td><td>vaccination rate control from (11)</td></tr><tr><td>ν̂</td><td>nominal isolation rate in (3)</td></tr><tr><td>u<sub>I</sub>(t)</td><td>isolation rate control from (11)</td></tr><tr><td>β̂</td><td>nominal transmission and contact in (3)</td></tr><tr><td>u<sub>C</sub>(t)</td><td>transmission and contact control from (11)</td></tr><tr><td>β̿</td><td>tuning constant in u<sub>C</sub>(t) from (11c)</td></tr><tr><td>X<sub>*</sub></td><td>equilibrium (S<sub>*</sub>, I<sub>*</sub>, Q<sub>*</sub>, R<sub>*</sub>) from (6)</td></tr><tr><td>λ and χ</td><td>λ = γ + ν̂ + μ and χ = ρ̂ + μ used in equilibrium (6)</td></tr><tr><td>D</td><td>set of feasible states (0, ∞)<sup>4</sup> for system (1)</td></tr><tr><td>H<sub>i</sub>, i = 1, 2, 3</td><td>functions (13) used in controls (11)</td></tr><tr><td>c<sub>◇</sub> and c</td><td>constants c<sub>◇</sub> ∈ (0, 2c<sub>◇</sub>) and c > 0 from (13)</td></tr><tr><td>c̄<sub>◇</sub></td><td>bound (12) related to tuning constant c<sub>◇</sub></td></tr><tr><td>(S̃, ξ̃, Q̃, R̃)</td><td>error states (S − S<sub>*</sub>, ln(I) − ln(I<sub>*</sub>), Q − Q<sub>*</sub>, R − R<sub>*</sub>) of (15)</td></tr><tr><td>ψ<sub>*</sub></td><td>constant λI<sub>*</sub> in (15)</td></tr><tr><td>D̃</td><td>feasible states (−S<sub>*</sub>, ∞) × ℝ × (−Q<sub>*</sub>, ∞) × (−R<sub>*</sub>, ∞) of (15)</td></tr><tr><td>P̃</td><td>perturbation set [−ψ<sub>*</sub>/4, ψ<sub>*</sub>/4] ∩ (−B, ∞) from Theorem 3.1</td></tr></tbody></table>

In addition to its mathematical novelty, the new results of this section are significant from the practical point of view, owing to the usefulness of manipulating not only vaccination rates, but also isolation and transmission and contact when
reducing infections. We use the following feedback control laws for the three control inputs:

$$ 
\begin{aligned}
u_V(t) &= f_V(S(t), I(t)) \\
&= \max \{-\hat{\rho}, \omega_V H_1(S(t), I(t))\},
\end{aligned} 
\quad (11a) $$

$$ 
\begin{aligned}
u_I(t) &= f_I(S(t), I(t), Q(t)) \\
&= \max \{-\hat{\nu}, \omega_I(H_2(S(t), I(t)) - I(t)H_3(Q(t)))\},
\end{aligned} 
\quad (11b) $$

$$ 
\begin{aligned}
u_C(t) &= f_C(S(t), I(t)) \\
&= \max \{\underline{\beta} - \hat{\beta}, \min \{0, \omega_C S(t)(I(t)H_1(S(t), I(t)) - H_2(S(t), I(t)))\}\}.
\end{aligned} 
\quad (11c) $$

The constant $\underline{\beta} \in [0, \hat{\beta}]$ and the nonnegative constants $\omega_V, \omega_I$, and $\omega_C$ are parameters with which one can tune the three control inputs (as we illustrate in Section 4 below), where in terms of the components (6) of our equilibrium point, we use the constant

$$ \bar{c}_{\diamond} = \frac{(\tau + \mu)c\lambda}{\hat{\nu}^2}, \quad (12) $$

and the functions

$$ H_1(S, I) = (1+c)\tilde{S} + c \left( \frac{\hat{\rho} + \mu}{\hat{\beta}} \ln \frac{I}{I_*} + \tilde{I} \right), \quad (13a) $$

$$ H_2(S, I) = c \left( \tilde{S} + \frac{\hat{\rho} + \mu}{\hat{\beta}} \ln \frac{I}{I_*} + \tilde{I} \right) \left( \frac{\hat{\rho} + \mu}{\hat{\beta}} + I \right) + \frac{(c+1)(\gamma + \hat{\nu} + \mu)}{\hat{\beta}} \tilde{I}, \quad (13b) $$

$$ \text{and } H_3(Q) = c_{\diamond} \tilde{Q}, \text{ where} \quad (13c) $$

$$ \tilde{S} = S - S_{\star}, \quad \tilde{I} = I - I_{\star}, \quad \text{and} \quad \tilde{Q} = Q - Q_{\star}, \quad (13d) $$

and where the constants $c_{\diamond} \in (0, 2\bar{c}_{\diamond})$ and $c > 0$ can also be used to tune the controls. The motivation for the controls (11) in combination with (13) is that they provide a negative definite decay condition on our Lyapunov function when the immigration perturbation $\epsilon(t)$ is the zero function, and an ISS Lyapunov function decay condition for nonzero perturbations, which signifies that they force convergence toward the desired equilibrium with an overshoot depending on the magnitude of the immigration perturbation in the sense of the ISS definition in Section 1; see Corollary 3.1. See Section 6 for more motivation for our controls, and our illustrations below, where the controls are shown to provide desirable population profiles using realistic parameter values from disease transmission data, which also motivates our controls (11). Then $H_1(S_\star, I_\star) = 0$, $H_2(S_\star, I_\star) = 0$, and $H_3(Q_\star) = 0$ hold. Also, in terms of the definitions

$$ \tilde{\xi} = \ln \frac{I}{I_*} \quad \text{and} \quad \psi_* = \lambda I_* \quad (14) $$

and the variables $\xi_* = \ln(I_\star)$ and $\xi = \ln(I)$, we have $\tilde{\xi} = \xi - \xi_*$, and the equations (13d) and (14) transform the model (1) into the dynamics

$$ \dot{\tilde{S}}(t) = \epsilon(t) - (\chi + \hat{\beta}e^{\tilde{\xi}(t)+\xi_{\star}}) \tilde{S}(t) + \psi_{\star}(1 - e^{\tilde{\xi}(t)}) - u_V(t)S(t) + u_C(t)e^{\xi(t)}S(t), \quad (15a) $$

$$ \dot{\tilde{\xi}}(t) = \hat{\beta}\tilde{S}(t) - u_I(t) + u_C(t)S(t), \quad (15b) $$

$$ \dot{\tilde{Q}}(t) = \hat{\nu}e^{\xi_{\star}}(e^{\tilde{\xi}(t)} - 1) - (\tau + \mu)\tilde{Q}(t) + u_I(t)e^{\xi(t)}, \quad (15c) $$

$$ \dot{\tilde{R}}(t) = \gamma e^{\xi_{\star}}(e^{\tilde{\xi}(t)} - 1) + \tau \tilde{Q}(t) - \mu \tilde{R}(t) + \hat{\rho} \tilde{S}(t) + u_V(t)S(t) \quad (15d) $$
which are defined on the state space

$$ \tilde{D} = (-S_{\star}, \infty) \times \mathbb{R} \times (-Q_{\star}, \infty) \times (-R_{\star}, \infty), \quad (16) $$

and where the formula (15) follows from our choices (5)-(6), i.e., the argument that produced [12, Equation (9)] in the $u_I = u_C = 0$ case, except with the additional control terms in the $\nu$ and $\beta$ formulas. We also use the set

$$ \tilde{P} = [-\psi_{\star}/4, \psi_{\star}/4] \cap (-B, \infty), \quad (17) $$

which will serve a disturbance set containing all admissible values of the immigration perturbation $\epsilon$ in our first theorem.

Our goal in this section is to specify conditions on the controls $u_C$, $u_I$, and $u_V$ such that (15) is ISS on its state space $\tilde{D}$ with respect to the disturbance set $\tilde{P}$ as defined in (17), using the ISS definition from Section 1 (but see Theorem 3.2 for an extension that allows the perturbation set (2)). Although our proof of the theorem will show that we can achieve the preceding ISS objectives with $u_V$, $u_I$, and $u_C$ all chosen to be the zero function, our proof and illustrations below will demonstrate how nonzero choices of $u_V$, $u_I$, and $u_C$ can lead to different solutions of the dynamics having faster convergence of the states towards their desired equilibrium values. This motivates our analysis of more general nonzero choices of $u_V$, $u_I$, and $u_C$ that can improve the convergence performance as compared to the results that would be obtained had we instead chosen these functions to be the zero function.

Using the function

$$ U(\tilde{S}, \tilde{\xi}, \tilde{Q}) = \frac{1}{2}\tilde{S}^2 + \frac{c}{2}\left[\tilde{S} + \frac{\chi}{\beta}\tilde{\xi} + I_{\star}\left(e^{\tilde{\xi}} - 1\right)\right]^2 + \frac{(c+1)\psi_{\star}}{\beta}\left(e^{\tilde{\xi}} - 1 - \tilde{\xi}\right) + \frac{c_{\diamond}}{2}\tilde{Q}^2 \quad (18) $$

of the difference variables from (13d) and (14), we get

$$ \frac{\partial U}{\partial \tilde{S}} = H_1(\tilde{S} + S_{\star}, I_{\star}e^{\tilde{\xi}}), \quad \frac{\partial U}{\partial \tilde{\xi}} = H_2(\tilde{S} + S_{\star}, I_{\star}e^{\tilde{\xi}}), \quad \text{and} \quad \frac{\partial U}{\partial \tilde{Q}} = H_3(\tilde{Q} + Q_{\star}), \quad (19) $$

where we used the fact that $\lambda = \psi_{\star}/I_{\star}$. The preceding choice of $U$ is inspired by, but different from, the Lyapunov function for a smaller SIQ error dynamics from [12], which in turn was obtained from a Matrosov type argument that transforms a Lyapunov function for a smaller SI error dynamics by adding an additional non-negative valued term. Although it has no specific physical interpretation, it is an essential ingredient for our mathematical analysis that builds our ISS Lyapunov function $V$ below. See also Remark 3.1 below for more comparisons with [12]. In terms of (18) and the preceding difference variables, we define the functions

$$ V(\tilde{S}, \tilde{\xi}, \tilde{Q}, \tilde{R}) = U(\tilde{S}, \tilde{\xi}, \tilde{Q}) + W(\tilde{S}, \tilde{\xi}, \tilde{Q}, \tilde{R}), \quad \text{where} \quad (20a) $$

$$ W(\tilde{S}, \tilde{\xi}, \tilde{Q}, \tilde{R}) = \frac{g}{2} \left[ \tilde{S} + I_{\star} \left( e^{\tilde{\xi}} - 1 \right) + \tilde{Q} + \tilde{R} \right]^2 \quad \text{and} \quad \tilde{R} = R - R_{\star}, \quad (20b) $$

for another tuning parameter $g > 0$. The following theorem establishes a property of the time derivative of the function $V$ along all solutions of the controlled SIQR model (1), which we use later to prove our ISS result:

**Theorem 3.1.** Let $c > 0$, $g > 0$, $c_{\diamond} \in (0, 2\bar{c}_{\diamond})$, $\bar{\beta} \in [0, \bar{\beta}]$, and the nonnegative values $\omega_V$, $\omega_I$, and $\omega_C$ be given constants. Then, the time derivative of the function $V$ in (20) satisfies

$$ \dot{V}(t) \leq -\alpha(V(\tilde{S}(t), \tilde{\xi}(t), \tilde{Q}(t), \tilde{R}(t))) + \sigma(|\epsilon(t)|) \quad (21) $$
along all trajectories of (1) in closed loop with the feedback controls that are given by (3), (11), and (13) for all $t \ge 0$ and all piecewise continuous functions

$$ \epsilon : [0, \infty) \to \tilde{P}, \qquad (22) $$

where

$$ \alpha(r) = \min \left\{ \frac{1}{2}, 1 - \sqrt{\frac{c_{\diamond}}{2\bar{c}_{\diamond}}} \right\} \left[ \sqrt{k_4 + \min \left\{ \frac{1}{k_3}, 4\sqrt{k_4\mu} \right\}} r - \sqrt{k_4} \right], \quad (23a) $$

$$ \sigma(r) = \left[ \frac{\ln(2)c\chi}{\hat{\beta}} + cI_{\star} \right] r + \left[ \frac{(1+c)^2}{2\chi} + \frac{g}{2\mu} \right] r^2, \qquad (23b) $$

$$ k_1 = \max \left\{ \frac{1+2c}{\chi}, \frac{2}{c\psi_{\star}I_{\star}} \left[ 2cI_{\star}^2 + \frac{(c+1)\psi_{\star}}{2\hat{\beta}} \right] \right\}, \qquad (23c) $$

$$ k_2 = k_1 + \left( \frac{4c\chi^2}{\hat{\beta}} + (c+1)\psi_{\star} \right) \frac{4}{c\psi_{\star}\left(2\chi + \hat{\beta}I_{\star}\right)}, \qquad (23d) $$

$$ k_3 = \left( \frac{2c\chi^2}{\hat{\beta}} + \frac{(c+1)\psi_{\star}}{2} \right) \frac{16\hat{\beta}}{c^2\psi_{\star}^2\chi^2}, \quad \text{and} \qquad (23e) $$

$$ k_4 = \frac{k_2^2}{4k_3^2}. \qquad (23f) $$

Also, for each value of $t \ge 0$ and for the preceding parameters, $\dot{V}$ is nonincreasing in each of the parameters $\omega_V, \omega_I$, and $\omega_C$. $\square$

*Proof.* We indicate the changes that are needed in the proof of [12, Theorem 3.1] to prove this theorem. The function $V$ in (20) has the same structure as the function $V_c$ that was used in [12, Section 4.3]. Therefore, its time derivative satisfies

$$ \dot{V}(t) \le -\frac{1}{2} W_c(\tilde{\xi}(t), \tilde{S}(t)) - \mu W(\tilde{S}(t), \tilde{\xi}(t), \tilde{Q}(t), \tilde{R}(t)) \\ + c_{\diamond} \left( \dot{\tilde{v}}e^{\xi_{\star}}(e^{\tilde{\xi}} - 1)\tilde{Q} - (\tau + \mu)\tilde{Q}^2 \right) + Z(\tilde{S}(t), \tilde{\xi}(t), \tilde{Q}(t)) + \sigma(|\epsilon(t)|) \quad (24) $$

along all trajectories of (15) in closed loop with the controls (3) for all $t \ge 0$ and for all piecewise continuous functions $\epsilon$ satisfying (22), where $W_c$ is the function

$$ W_c(a_1, a_2) = \left( \chi + \hat{\beta}e^{a_1+\xi_{\star}} \right) a_2^2 + c\psi_{\star} \left[ \frac{\chi}{\hat{\beta}}a_1 + e^{\xi_{\star}}(e^{a_1}-1) \right] (e^{a_1}-1) $$

in the proof of [12, Theorem 3.1], and where the added function $Z$ is defined by

$$ \begin{aligned} Z(\tilde{S}, \tilde{\xi}, \tilde{Q}) &= \frac{\partial U}{\partial \tilde{S}}(-f_V(S, I)S - f_C(S, I)IS) \\ &\quad + \frac{\partial U}{\partial \tilde{\xi}}(-f_I(S, I, Q) + f_C(S, I)S) + \frac{\partial U}{\partial \tilde{Q}}f_I(S, I, Q)I \\ &= -SH_1(S, I)f_V(S, I) \\ &\quad + (SH_2(S, I) - ISH_1(S, I))f_C(S, I) \\ &\quad + (IH_3(Q) - H_2(S, I))f_I(S, I, Q). \end{aligned} $$

In fact, the inequality (24) follows from (19) and the first inequality in [12, Equation (71)] except with the expression

$$ \mathcal{H}(\tilde{Q}, \tilde{\xi}) = c_{\diamond} \frac{\tau+\mu}{2} \tilde{Q}^2 + c_{\diamond} \frac{\nu^2 e^{2\xi_{\star}}}{2(\tau+\mu)} [e^{\tilde{\xi}} - 1]^2 $$
that appears in [12, Equation (71)] replaced by the formula

$$ \dot{F}_3 = -(\tau + \mu)\tilde{Q}^2 + \nu e^{\xi_*} (e^{\tilde{\xi}} - 1) \tilde{Q} $$

for the time derivative of $F_3(\tilde{Q}) = \frac{1}{2}\tilde{Q}^2$ from [12, Equation (69)].

Young's inequality $rs \le nr^2/2 + s^2/(2n)$ with $r = \tilde{Q}$, $s = \hat{\nu}e^{\xi_*}(e^{\tilde{\xi}} - 1)$, and $n = b(\tau+\mu)$ gives

$$ c_{\diamond} (\hat{\nu}e^{\xi_*} (e^{\tilde{\xi}} - 1)\tilde{Q} - (\tau + \mu)\tilde{Q}^2) \le -(2-b) \frac{c_{\diamond}(\tau+\mu)}{2} \tilde{Q}^2 + \frac{c_{\diamond}\hat{\nu}^2 e^{2\xi_*}}{2b(\tau+\mu)} (e^{\tilde{\xi}} - 1)^2 \quad (25) $$

for any constant $b > 0$. Choose the constants

$$ b = \sqrt{\frac{2c_{\diamond}}{\bar{c}_{\diamond}}} $$

and $\delta = 2 - b$. Then $\delta = 2(1 - \sqrt{c_{\diamond}/(2\bar{c}_{\diamond})}) > 0$ holds, and our definition of $\bar{c}_{\diamond}$ from (12) and our choice $\psi_{\star} = \lambda I_{\star}$ from (14) give

$$ \frac{c_{\diamond} \hat{\nu}^2 e^{2\xi_{\star}}}{2b(\tau + \mu)} = \frac{\ell c \psi_{\star} e^{\xi_{\star}}}{2}, \quad (26) $$

where $\ell = \sqrt{c_{\diamond}/2\bar{c}_{\diamond}}$. Hence, with the choice $\delta = 2 - b$, we conclude from (25) that

$$ c_{\diamond} \left( \hat{\nu} e^{\xi_{\star}} (e^{\tilde{\xi}} - 1) \tilde{Q} - (\tau + \mu) \tilde{Q}^2 \right) \le -\delta \frac{c_{\diamond}(\tau + \mu)}{2} \tilde{Q}^2 + \frac{\ell c \psi_{\star} e^{\xi_{\star}}}{2} (e^{\tilde{\xi}} - 1)^2. \quad (27) $$

Since $\tilde{\xi}(e^{\tilde{\xi}} - 1) \ge 0$, we get $W_c(\tilde{\xi}(t), \tilde{S}(t)) \ge c\psi_{\star}e^{\xi_{\star}}(e^{\tilde{\xi}(t)} - 1)^2$, so from (24) we obtain

$$ \begin{aligned} \dot{V}(t) &\le -\frac{1-\ell}{2} W_c(\tilde{\xi}(t), \tilde{S}(t)) - \mu W(\tilde{S}(t), \tilde{\xi}(t), \tilde{Q}(t), \tilde{R}(t)) \\ &\quad - \frac{\delta c_{\diamond}(\tau+\mu)}{2} \tilde{Q}^2(t) + Z(\tilde{S}(t), \tilde{\xi}(t), \tilde{Q}(t)) + \sigma(|\epsilon(t)|) \\ &\le -\min\{2(1-\ell), 1, \delta\} \left[ \frac{1}{4} W_c(\tilde{\xi}(t), \tilde{S}(t)) + \mu W(\tilde{S}(t), \tilde{\xi}(t), \tilde{Q}(t), \tilde{R}(t)) \right. \\ &\qquad \left. + \frac{c_{\diamond}(\tau+\mu)}{2} \tilde{Q}^2(t) \right] + Z(\tilde{S}(t), \tilde{\xi}(t), \tilde{Q}(t)) + \sigma(|\epsilon(t)|) \\ &\le -\alpha(V(\tilde{S}(t), \tilde{\xi}(t), \tilde{Q}(t), \tilde{R}(t))) + Z(\tilde{S}(t), \tilde{\xi}(t), \tilde{Q}(t)) + \sigma(|\epsilon(t)|), \end{aligned} \quad (28) $$

where the last inequality follows from $\delta = 2(1 - \ell)$ and the same reasoning that gave [12, Equation (72)]. Next note that our formulas (11) for $f_V$, $f_I$, and $f_C$ give

$$ -SH_1(S, I)f_V(S, I) \le 0, \quad (29a) $$

$$ (IH_3(Q) - H_2(S, I))f_I(S, I, Q) \le 0, \quad (29b) $$

$$ (SH_2(S, I) - ISH_1(S, I))f_C(S, I) \le 0. \quad (29c) $$

Therefore, (21) follows from (28). Finally, the nonincreasing property from the statement of the theorem also follows from (29). □

Theorem 3.1 holds for any value of the tuning constant $g > 0$ that we used in V. The magnitude and the convergence speed of V estimated by (21) depend on $g$ as $g$ appears in $\sigma$. However, the control laws defined by (11) and (13) are independent of $g$. By contrast, the parameters $c$ and $c_{\diamond}$ used in V appear in the control laws. For the transformed system (15), the facts that $\alpha \in K_{\infty}$, combined with the proper and positive definiteness of V and the inequality (21) in Theorem 3.1, guarantee that V is an ISS Lyapunov function on $\tilde{D}$ for the perturbation set $\tilde{P}$ [35]. Using
ISS theory [34, 35] that ensures that the existence of an ISS Lyapunov function is sufficient for ISS, the next corollary on the error variable

$$ (\tilde{S}, \tilde{\xi}, \tilde{Q}, \tilde{R}) = (S - S_{\star}, \ln I - \ln I_{\star}, Q - Q_{\star}, R - R_{\star}) \quad (30) $$

is a direct consequence of Theorem 3.1, using $V$ as the ISS Lyapunov function:

**Corollary 3.1.** For any constants $c > 0$, $c_{\diamond} \in (0, 2\bar{c}_{\diamond})$, and $\beta \in [0, \hat{\beta}]$, and for any nonnegative real values $\omega_V$, $\omega_I$, and $\omega_C$, the dynamics for the error variable (30) for the SIQR model (1) with the controller consisting of (3), (11) and (13) are ISS on the entire state space $\mathcal{D}$ for the perturbation $\epsilon$ satisfying (22).

Corollary 3.1 restricts the perturbation $\epsilon$ to (22), i.e., the disturbance set $\tilde{\mathcal{P}}$ in (17). However, transforming the result of Theorem 3.1 back into the original error variable

$$ (\tilde{S}, \tilde{I}, \tilde{Q}, \tilde{R}) = (S - S_{\star}, I - I_{\star}, Q - Q_{\star}, R - R_{\star}), \quad (31) $$

we can also prove the following, which relaxes the constraint on the perturbation values, and (by the definition of ISS that we provided in Section 1) provides an ‘unrestrictively’ local ISS result that applies for initial states that are contained in the set $\mathcal{D}_a$ (where the set $\mathcal{D}_a$ is bounded but arbitrarily large, and the disturbance set $\mathcal{P}$ is unrestricted):

**Theorem 3.2.** Let $\mathcal{D}_a$ be any compact subset of $\mathcal{D}$ that contains $X_{\star}$. For any constants $c > 0$, $c_{\diamond} \in (0, 2\bar{c}_{\diamond})$, and $\beta \in [0, \hat{\beta}]$, and any nonnegative values $\omega_V$, $\omega_I$, and $\omega_C$, the dynamics for the error variable (31) for the SIQR model (1) with the controller consisting of (3), (11) and (13) are ISS on $\mathcal{D}_a$ for the perturbation $\epsilon$ valued in $\mathcal{P} = (-B, \infty)$. $\square$

*Proof.* Define $\tilde{x}(t) = (\tilde{S}(t), \tilde{\xi}(t), \tilde{Q}(t), \tilde{R}(t))$, and let $\tilde{\mathcal{D}}$ be defined as (16). We also use the variables $\tilde{X}(t) = (\tilde{S}(t), \tilde{I}(t), \tilde{Q}(t), \tilde{R}(t))$, $X(t) = \tilde{X}(t) + X_{\star}$, and the set $\tilde{\mathcal{D}}_A = \{\tilde{X} \in (-\infty, \infty)^4 : \tilde{X} + X_{\star} \in \mathcal{D}_a\}$. As Corollary 3.1 states the ISS property, inequality (21) ensures the existence of $\Phi \in \mathcal{KL}$ and $\Gamma \in \mathcal{K}_{\infty}$ such that

$$ |V(\tilde{x}(t))| \le \Phi(|V(\tilde{x}(0)|, t) + \Gamma(|\epsilon|_{[0,t]}), \quad (32) $$

for all $t \ge 0$, $\tilde{x}(0) \in \tilde{\mathcal{D}}$, and $\epsilon$'s that are valued in $\tilde{\mathcal{P}}$ [35]. Define

$$ \bar{\theta}(s) = \max_{\tilde{X} \in \tilde{\mathcal{D}}_A : |\tilde{X}| \le s} V \left( (\tilde{X}_1, \ln(\tilde{X}_2 + I_{\star}) - \ln I_{\star}, \tilde{X}_3, \tilde{X}_4) \right), $$

$$ \underline{\theta}(s) = \min_{\tilde{X} \in \tilde{\mathcal{D}}_A : |\tilde{X}| \ge s} V \left( (\tilde{X}_1, \ln(\tilde{X}_2 + I_{\star}) - \ln I_{\star}, \tilde{X}_3, \tilde{X}_4) \right) $$

for $s \in [0, \infty)$. By the construction of $V$ in (20), we have $\underline{\theta} \in \mathcal{K}_{\infty}$. Note that for large $s$ for which $\{\tilde{X} \in \tilde{\mathcal{D}}_A : |\tilde{X}| \ge s\}$ is empty, the function $\underline{\theta}(s)$ is extended freely to form a class $\mathcal{K}_{\infty}$ function. There also exists $\hat{\theta} \in \mathcal{K}_{\infty}$ such that $\underline{\theta}(s) \le \hat{\theta}(s)$ for all $s \in [0, \infty)$. Hence, for $\Phi(s, t) = \underline{\theta}^{-1}(2\Phi(\hat{\theta}(s), t))$ and $\Gamma(s) = \underline{\theta}^{-1}(2\Gamma(s))$, which are of class $\mathcal{KL}$ and $\mathcal{K}_{\infty}$, respectively, from (32) we obtain

$$ |\tilde{X}(t)| \le \Phi(|\tilde{X}(0)|, t) + \Gamma(|\epsilon|_{[0,t]}), \quad (33) $$

for all $t \ge 0$, $\tilde{X}(0) \in \tilde{\mathcal{D}}_A$, and $\epsilon$'s that are valued in $\tilde{\mathcal{P}}$ by the subadditivity property $\underline{\theta}^{-1}(a + b) \le \underline{\theta}^{-1}(2a) + \underline{\theta}^{-1}(2b)$ for all $a \ge 0$ and $b \ge 0$. Since the SIQR model (1) satisfies

$$ \frac{d}{dt}(S(t) + I(t) + Q(t) + R(t)) = B + \epsilon(t) - \mu(S(t) + I(t) + Q(t) + R(t)), $$
we have

$$S(t) + I(t) + Q(t) + R(t) \leq e^{-\mu t}(S(0) + I(0) + Q(0) + R(0)) \\ + \sup_{\tau \in [0,t]} \frac{B + \epsilon(\tau)}{\mu}$$

and so also

$$|X(t)| \leq 4e^{-\mu t}|X(0)| + \sup_{\tau \in [0,t]} \frac{B+\epsilon(\tau)}{\mu}$$

by the subadditivity of the square root, for all $X(0) \in [0, \infty)^4$ and all choices of the function $\epsilon$ that are valued in $\mathcal{P}$. Hence, since $|\tilde{X}(t)| \leq |X(t)| + |X_\star|$ for all $t \geq 0$ and $|X(0)| \leq |\tilde{X}(0)| + |X_\star|$, we get

$$|\tilde{X}(t)| \leq 4e^{-\mu t}|\tilde{X}(0)| + 5|X_\star| + \sup_{\tau \in [0,t]} \frac{B+\epsilon(\tau)}{\mu} \quad (34)$$

for all $t \geq 0$ and all choices of the function $\epsilon$ that are valued $\mathcal{P}$. Pick any $\bar{\Phi} \in \mathcal{KL}$ and $\bar{\Gamma} \in \mathcal{K}_\infty$ satisfying

$$\bar{\Phi}(s,t) \geq \max\{\Phi(s,t), 4se^{-\mu t}\} \text{ for all } t \geq 0 \text{ and } s \geq 0, \\ \bar{\Gamma}(s) \geq \Gamma(s) \text{ for all } s \geq 0, \text{ and}$$

$$\bar{\Gamma}(s) \geq \frac{s+B}{\mu} + 5|x_{\star}| \text{ for all } s \in [\min\{B, \psi_{\star}/4\}, \infty).$$

By separately considering the cases where $\epsilon$ satisfies $|\epsilon|_\infty \in [0, \min\{B, \psi_\star/4\})$ and where the preceding inclusion is violated, we can combine (33) with (34) to obtain

$$|\tilde{X}(t)| \leq \bar{\Phi}(|\tilde{X}(0)|, t) + \bar{\Gamma}(|\epsilon|_\infty) \quad (36)$$

for all $t \geq 0$, $\tilde{X}(0) \in \tilde{\mathcal{D}}_A$, and $\epsilon$'s that are valued in $\mathcal{P}$. Since the solution $\tilde{X}(t)$ at each time $t$ does not depend on $\epsilon(\ell)$ values with $\ell > t$, we can replace $|\epsilon|_\infty$ by $|\epsilon|_{[0,t]}$ in (36) to get the desired ISS estimate. $\square$

**Theorem 3.2** by itself cannot meet the goal of this paper, since in particular, it does not give us any way of constructing controllers. Instead, Theorem 3.2 gives us a useful guarantee for a given controller. The controller consisting of (3), (11) and (13) is derived from the ISS Lyapunov function developed in Theorem 3.1. The Lyapunov function judiciously determines the control action on the entire state space to achieve the guarantee stated in Theorem 3.2.

**Remark 3.1.** The construction of the Lyapunov function $V$ as in (20) employs the basic idea proposed in [12]. In [12], the value $c_\diamond = \bar{c}_\diamond/2$ was used. This paper relaxes the choice to $c_\diamond \in (0, 2\bar{c}_\diamond)$. This relaxation made it possible to introduce control inputs to isolation and contact regulation. Having the flexible parameter to pursue an effective balance between the three control inputs can be useful for better disease control. The parameter $c$ in $V$ can also be used for this purpose.

**Remark 3.2.** If the transformed SIQR model (15) is written as $\dot{x} = F(x) + G_V(x)u_V + G_I(x)u_I + G_C(x)u_C$ with $x = (\tilde{S}, \tilde{\xi}, \tilde{Q}, \tilde{R})$, then, in terms of the $V$ given by (20), the control laws (11b) and (11c) are

$$-\omega_I \frac{\partial V}{\partial \tilde{x}} G_I, \text{ and } -\omega_C \frac{\partial V}{\partial \tilde{x}} G_C, \quad (37)$$
respectively, with maximization and minimization meeting the input constraints. For the control law (11a), this paper employs the weighting 1/S as

$$ -\omega_V \frac{1}{S} \frac{\partial V}{\partial \tilde{x}} G_V, \qquad (38) $$

since 1/S is canceled out by $G_V$. The weighting avoids the vaccination rate in [12] that becomes very small too quickly.

**Remark 3.3.** The function $V$ in this paper captures the behavior of all the variables $S$, $I$, $Q$, and $R$ as in (20). By the strict Lyapunov function decay condition (21), the function also characterizes the attractiveness of the endemic equilibrium without using LaSalle's invariance principle. The construction and the use of such a comprehensive Lyapunov function $V$ allows us to derive control laws from its gradient (37). By contrast, in [8], the Lyapunov-like function

$$ V(\tilde{S}, \tilde{I}) = S_{\star} \int_{S_{\star}}^{\tilde{S}+S_{\star}} \frac{\hat{\beta}r - (\gamma + \hat{\nu} + \mu)}{r} dr + \int_{I_{\star}}^{\tilde{I}+I_{\star}} \frac{\hat{\beta}S_{\star}r + \mu S_{\star} - B}{r} dr \quad (39) $$

is used for stability analysis assuming constant immigration and no manipulation, i.e., $\epsilon(t) = 0$ and $\rho(t) = u_V(t) = u_I(t) = u_C(t) = 0$ for all $t \ge 0$. Our model (1) agrees with the SIQR model in [8, Section 5], except with the disease related death in [8] combined in $R$ as a standard practice. The function $V$ in (39) does not contain the variables $Q$ and $R$. As noted in the proof of [8, Theorem 5], its time derivative

$$ \dot{V}(t) = - \frac{\hat{\beta} B \tilde{S}^2(t)}{\tilde{S}(t) + S_{\star}} \qquad (40) $$

involves only $\tilde{S}$. Hence, the Lyapunov-like function does not characterize the behavior of $Q$ and $R$. Moreover, concluding the analysis of $S$ and $I$ needs LaSalle's invariance principle. Even if the immigration perturbation $\epsilon(t) \neq 0$ and the manipulation inputs $\rho(t) \neq 0$ and $u_I(t) \neq 0$ are added to the SIQR model (to obtain the dynamics (1) except with $u_C = 0$), using

$$ V(\tilde{S}, \tilde{I}) = S_{\star} \int_{S_{\star}}^{\tilde{S}+S_{\star}} \frac{\hat{\beta}r - (\gamma + \hat{\nu} + \mu)}{r} dr + \int_{I_{\star}}^{\tilde{I}+I_{\star}} \frac{\hat{\beta}S_{\star}r + (\hat{\rho} + \mu)S_{\star} - B}{r} dr \quad (41) $$

only yields

$$ \begin{aligned} \dot{V}(t) &= \frac{S_{\star}}{S(t)} (\hat{\beta}S(t) - \hat{\beta}S_{\star}) (B + \epsilon(t) - (\hat{\rho} + \mu + u_V(t))S(t) - \hat{\beta}I(t)S(t)) \\ &\quad + (\hat{\beta}S_{\star}I(t) + (\hat{\rho} + \mu)S_{\star} - B) (\hat{\beta}S(t) - \hat{\beta}S_{\star} - u_I(t)) \\ &= - \frac{\hat{\beta}B\tilde{S}(t)^2 - \hat{\beta}S_{\star}\tilde{S}(t)(\epsilon(t) - u_V(t)S(t))}{\tilde{S}(t) + S_{\star}} \\ &\quad - (\hat{\beta}S_{\star}I(t) + (\hat{\rho} + \mu)S_{\star} - B)u_I(t). \end{aligned} \qquad (42) $$

Here, $\hat{\beta}S_{\star} = \gamma + \hat{\nu} + \mu$ is used. Hence, letting $u_V$ and $u_I$ be of the gradient form (37) with (41) does not produce the required strict Lyapunov function decay condition. The same is true if we allow $u_C \neq 0$. Hence, the gradient of (41) does not lend itself to deriving control laws.

**Remark 3.4.** It is easy to verify that Theorem 3.1 (hence, Corollary 3.1 and Theorem 3.2) remain true even if (11a) and (11b) are replaced by

$$ u_V(t) = f_V(S(t), I(t)) $$
$$ = \max \{-\hat{\rho}, \min \{\bar{\rho} - \hat{\rho}, \omega_V H_1(S(t), I(t))\}\}, \quad (43a) $$

$$ \begin{aligned} u_I(t) &= f_I(S(t), I(t), Q(t)) \\ &= \max \{-\hat{\nu}, \min \{\bar{\nu} - \hat{\nu}, \omega_I(H_2(S(t), I(t)) - I(t)H_3(Q(t)))\}\} \quad (43b) \end{aligned} $$

where the saturation levels $\bar{\rho} > \hat{\rho}$ and $\bar{\nu} > \hat{\nu}$ can be chosen arbitrarily. This applies throughout this paper. These saturated inputs reasonably represent limitations of resources, but the saturation is not necessary from the mathematical point of view.

4. **Simulation and discussion.** In this section and Section 6, we illustrate the effectiveness of the proposed control design through numerical simulations. We use

$$ \hat{\beta} = 0.126/N, \quad \gamma = 0.03, \quad \text{and} \quad \tau = 0.03, \quad (44) $$

for which the SIQR model (1) is used to analyze the COVID-19 outbreak in Japan in [26]. The unit of population is in millions, and the time $t$ is in days. The number $N$ denotes the total population. The demographic parameters

$$ N = 126, \quad \mu = 0.0000307, \quad \text{and} \quad B = 3110 \times 10^{-6} \quad (45) $$

are borrowed from [13, 36] since the inflow and the outflow are not of interest in [26]. We chose the immigration perturbation

$$ \epsilon(t) = -311 \times 10^{-6} \cos(\pi t/150) \quad (46) $$

to describe the 20% perturbation of immigrants and newborns. Unless otherwise stated, for the control laws (11) in (3), we use

$$ \hat{\rho} = 0.00005, \quad \omega_V = 0.000015, \quad c = 0.02, \quad (47a) $$

$$ \hat{\nu} = 0.005, \quad \omega_I = 0.00006, \quad c_{\diamond} = 1.8\bar{c}_{\diamond}, \quad (47b) $$

$$ \underline{\beta} = \hat{\beta}/4, \quad \text{and} \quad \omega_C = 0.0000001. \quad (47c) $$

The basic reproduction number $\hat{R}_0$ becomes 1.1001.

4.1. **Efficacy of the proposed control laws.** Figure 1a shows the four populations of (1) without vaccination, isolation, and contact regulation, i.e, (1) with the preceding parameter values but with $\rho = 0$, $\nu = 0$, and $\beta = \hat{\beta}$, while Figure 1b shows the populations controlled by (11) in (3). By the control, the infected population $I$ decreases fast. The peaks of $I$ and $Q$ are also reduced significantly. The remaining susceptible population $S$ in Fig. 1b is much larger than that in Fig. 1a. Note that $R$ in Fig. 1b includes the population of vaccinated individuals who are immune. In Fig. 1b, the isolated population $I$ is too small to be seen clearly since the no-isolation $\nu = 0$ implies the convergence of $I$ to zero.

4.2. **Vaccination of susceptible individuals.** The simulation is performed with the constant vaccination $\omega_V = 0$ in Fig. 2. Compared to the case of $\omega_V \neq 0$ in Fig. 1b, the population increase of recovered individuals is delayed. However, the amount of $I(t) + Q(t)$ of Fig. 3b is still substantially smaller than that of Fig. 1a. Nevertheless, the benefit of updating the vaccine rate is clearly seen in the significant reduction of the sum of $I(t) + Q(t)$ by comparing Fig. 3a with Fig. 3b. The accumulated total $\int_0^t \rho(\ell)S(\ell)d\ell$ of vaccinated individuals is also shown in Figs. 3a and 3b. It is observed that the control input $u_V$ slows the vaccination rapidly right before the infection peak. In other words, $\nu$ and $\beta$ influence (1) after the surge of $I$.
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Comparison of uncontrolled and controlled populations of (1).</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIGURE 2. Populations of (1) with constant vaccination $\omega_V = 0$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIGURE 3. Infected plus isolated population and accumulated total of vaccinated individuals.</figcaption></figure>

**4.3. Isolation of infected individuals.** Unlike vaccination, the effect of isolation of infected individuals has an immediate effect to reduce the disease spread [14]. Isolation is limited by the number of hospital beds and facilities including human
resources. Under the assumption of availability of isolation facilities including quarantine at home for humans, updating the isolation rate can be easier and quicker than updating vaccination speed.

The Lyapunov function $V$ defined in (20) has the $Q$-term

$$V_Q = \frac{c_{\diamond}}{2} \tilde{Q}^2. \qquad (48)$$

The proposed control laws in (11) navigate the system trajectory to decrease the $V$ values. Hence, the proposed control decreases $\tilde{Q}$ to zero. The weighting coefficient $c_{\diamond}$ is free in the interval $(0, 2\bar{c}_{\diamond})$. The larger $c_{\diamond}$ is, the more the control law in (11b) puts weight on the reduction of $\tilde{Q}$ than the other variables $\tilde{S}, \tilde{I}$, and $\tilde{R}$. Among the three control inputs (11), only (11b) involves $\tilde{Q}$ and $c_{\diamond}$.

The positive value $\bar{c}_{\diamond}$ given by (12) is a decreasing function of $\hat{v}$. If the resources for isolation are limited, the function $\bar{c}_{\diamond}$ gives a guideline telling how small the nominal isolation rate $\hat{v}$ should be. This mechanism is illustrated by the difference between the simulations in Figs. 1b and 4.

<figure><img src="image_5.png" /><figcaption>FIGURE 4. Populations of (1) with isolation coefficient $c_{\diamond} = 0.18\bar{c}_{\diamond}$.</figcaption></figure>

### 4.4. Manipulation of contact rate.
The parameter $\beta$ describes not only the disease transmission rate, but also the contact rate. For human diseases, a government can help control the contact rate by issuing a lock-down or a similar regulation and advisory. With a positive constant $\omega_C$ in the controls, the proposed control law (11c) updates $\beta$ in (3) so that the system trajectory leads to smaller $V$ values. Importantly, for human diseases, no government wants to spread the disease faster than a natural infection speed even temporarily. Therefore, the control input $\beta$ is one-sided, i.e., $\beta(t) = \hat{\beta}(t) + u_C(t)$, where $u_C(t)$ is non-positive in (11c). The target equilibrium $S_*$ of the susceptible population $S$ corresponding to the target number $I_*$ of the infected population is normally chosen to be much smaller than the initial population $S$. In fact, if the society can eliminate the disease completely, the equilibrium $S_*$ is generated only by the inflow, and it is almost zero, compared with the initial $S$, which must be the total population of the society. Rapidly increasing $\beta$ can be effective for large $S$ [39]. Since it is a prohibited strategy for humans, the contact rate in (11c) is nonzero only when $\beta$ can be decreased. As seen in Figs. 5 and 6, the total number of infected individuals is reduced about 27% at 250 [days] by regulating society or environment, i.e., lowering $u_C$ at appropriate timing. Here, the total number of infected individuals plotted on Fig. 5 is the integral of
$\beta(t)S(t)I(t)$ from $t = 0$. The contact regulation curbs the infection as seen in Fig. 1b compared to Fig. 7.

<figure><img src="image_6.png" /><figcaption>FIGURE 5. Total of infected individuals of (1).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 6. Control input of (1) for contact rate.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 7. Populations of (1) with constant contact $\omega_C = 0$.</figcaption></figure>
4.5. **Delays in the feedback control inputs.** Plots in Fig. 8 are computed for 2 days delays in the control inputs $u_V$, $u_I$, and $u_C$, meaning $t$ is replaced by $t-2$ in the controls for $t \ge 2$. The delays result in only a small change from the no delay case (Fig. 1b). As we saw in the proof of Theorem 3.2, system (1) is population behavior that inherits a conservation mechanism irrespective of control laws, the infection remains bounded for a bounded immigration perturbation. However, the delay can hamper convergence of the trajectories to an arbitrarily small neighborhood of the target equilibrium even for a sufficiently small immigration perturbation. Therefore, the remainder of this paper shows the existence of a non-zero upper bound on the input delays for which ISS is guaranteed with respect to inflow perturbation.

<figure><img src="image_9.png" /><figcaption>FIGURE 8. Populations of (1) controlled by (11) in (3) with 2 day delay.</figcaption></figure>

5. **Delay robustness.** We next consider the effect of delays in the control inputs, which we represent by

$$ \rho(t) = \hat{\rho} + u_V(t - L_V), \quad (49a) $$

$$ \nu(t) = \hat{\nu} + u_I(t - L_I), \quad (49b) $$

$$ \beta(t) = \hat{\beta} + u_C(t - L_C) \quad (49c) $$

for nonnegative constant delay lengths $L_V$, $L_I$, and $L_C$, where the control law functions $u_V$, $u_I$, and $u_C$ are given by (11). Define $L = \max\{L_V, L_I, L_C\}$ and $\tilde{x} = (\tilde{S}, \tilde{\xi}, \tilde{Q}, \tilde{R})$. Let $\tilde{x}_0(\tau) = \tilde{x}(\tau)$ for $\tau \in [-L, 0]$ be an arbitrary continuous function that is valued in the state space $\tilde{D} = (-S_*, \infty) \times \mathbb{R} \times (-Q_*, \infty) \times (-R_*, \infty)$ that we defined in (16). Using the definition from [29], we say that the SIQR model with (49) is ISS on a set $D_a$ containing $X_*$ for a perturbation set $P_a$ provided there exist $\Phi \in \mathcal{KL}$ and $\Gamma \in \mathcal{K}_\infty$ satisfying $|\tilde{x}(t)| \le \Phi(|\tilde{x}_0|_{[-L,0]}, t) + \Gamma(|\epsilon|_{[0,t]})$ for all $t \ge 0$, for all initial functions $\tilde{x}_0$ valued in $D_a$ and for all locally bounded piecewise continuous functions $\epsilon$ that are valued in $P_a$. Using the notation we introduced above, we can then prove the following (but see Remark 5.1 below for comments on applications that allow longer delays for which our condition (50) is not satisfied):

**Theorem 5.1.** Let $D_a$ be any compact subset of $D = (0, \infty)^4$ that contains $X_* = (S_*, I_*, Q_*, R_*)$. Then for any constants $c > 0$, $c_\diamond \in (0, 2\bar{c}_\diamond)$, $\beta \in [0, \hat{\beta}]$, and $\bar{L} > 0$, we can construct positive numbers $A$ and $\bar{\epsilon}$ such that if nonnegative constants $\omega_V$, $\omega_I$, $\omega_C$, $L_V$, $L_I$, and $L_C$ satisfy $\max\{L_V, L_I, L_C\} \le \bar{L}$ and

$$ (1 + \max\{\omega_V, \omega_I, \omega_C\})(\omega_V + \omega_I + \omega_C) \max\{L_V, L_I, L_C\} < A, \quad (50) $$
then the dynamics for the error variable (30) for the SIQR model (1) with the controller consisting of (11), (13), and (49) are ISS on $\mathcal{D}_a$ for the perturbation set $(-\min\{B, \bar{\epsilon}\}, \bar{\epsilon})$, where the positive number $\bar{\epsilon}$ is a non-increasing function of $\omega_V$, $\omega_I$, $\omega_C$, $L_V$, $L_I$, and $L_C$.

Proof. First, recall that the SIQR model (1) satisfies $\dot{S}+\dot{I}+\dot{Q}+\dot{R}=B+\epsilon-\mu(S+I+Q+R)$ independently of the control inputs and the initial functions. Therefore, for each choice of the set $\mathcal{D}_a$, we obtain a bounded set $\mathcal{Z} \subset \mathcal{D}$ that contains the solutions $X(t) = (S(t), I(t), Q(t), R(t))$ of (1) for all $t \ge 0$, all perturbations that are valued in $\tilde{\mathcal{P}}$, and all initial functions $X(\tau)$ for $\tau \in [-L, 0]$ that are valued in $\mathcal{D}_a$. The set $\mathcal{Z}$ is independent of $X(\tau)$ for $\tau \in [-L, 0)$. Define

$$ \begin{aligned} \tilde{\mathcal{D}}_a &= \{ (\tilde{S}, \tilde{\xi}, \tilde{Q}, \tilde{R}) \in \tilde{\mathcal{D}} : (\tilde{S} + S_\star, I_\star e^{\tilde{\xi}}, \tilde{Q} + Q_\star, \tilde{R} + R_\star) \in \mathcal{D}_a \}, \\ \tilde{\mathcal{Z}} &= \{ (\tilde{S}, \tilde{\xi}, \tilde{Q}, \tilde{R}) \in \tilde{\mathcal{D}} : (\tilde{S} + S_\star, I_\star e^{\tilde{\xi}}, \tilde{Q} + Q_\star, \tilde{R} + R_\star) \in \mathcal{Z} \}. \end{aligned} $$

The set $\tilde{\mathcal{Z}}$ can be unbounded, while $\tilde{\mathcal{D}}_a$ is bounded. The possible unboundedness of $\tilde{\mathcal{Z}}$ follows because $\tilde{\xi}(t) = \ln I(t) - \ln I_\star$ converges to $-\infty$ as $I(t) \to 0^+$, and because $I(t)$ is not guaranteed to have a positive lower bound. We use the variable $\tilde{x}(t) = (\tilde{S}(t), \tilde{\xi}(t), \tilde{Q}(t), \tilde{R}(t))$. For $\tilde{x}(t)$ valued in $\tilde{\mathcal{Z}}$, define

$$ \begin{aligned} \mathcal{G}_V(t) &= -S(t)H_1(S(t), I(t))(u_V(t - L_V) - u_V(t)), \\ \mathcal{G}_I(t) &= (I(t)H_3(Q(t)) - H_2(S(t), I(t))(u_I(t - L_I) - u_I(t)), \\ \mathcal{G}_C(t) &= (S(t)H_2(S(t), I(t)) - I(t)S(t)H_1(S(t), I(t))(u_C(t - L_C) - u_C(t)) \end{aligned} $$

with (11) and (13). These functions are the zero functions when the corresponding delays $L_V$, $L_I$, and $L_C$ are zero. From (13) one can verify the existence of non-negative constants $k_{a,i}$, $k_{b,i}$, and $k_{c,i}$ such that

$$ |H_1(S, I)| \le k_{a,1}|\tilde{S}| + k_{b,1}|\tilde{\xi}| + k_{c,1}|\tilde{Q}|, \quad (51a) $$

$$ |H_2(S, I)| \le k_{a,2}|\tilde{S}| + k_{b,2}|\tilde{\xi}| + k_{c,2}|\tilde{Q}|, \quad (51b) $$

$$ |H_3(Q)| \le k_{a,3}|\tilde{S}| + k_{b,3}|\tilde{\xi}| + k_{c,3}|\tilde{Q}| \quad (51c) $$

for all $\tilde{x} \in \tilde{\mathcal{Z}}$. The existence of the constants in (51) follows from the boundedness of the partial derivatives of $H_1$, $H_2$, and $H_3$ on $\tilde{\mathcal{Z}}$, which in turn follows from

$$ \frac{d}{d\tilde{\xi}} \tilde{I} = \frac{d}{d\xi} e^{\tilde{\xi}} = e^{\tilde{\xi}} = \tilde{I} $$

and its boundedness on the set $\{\tilde{\xi} : \tilde{x} \in \tilde{\mathcal{Z}}\}$. Let the feedback laws of the control inputs $u_V$, $u_I$, and $u_C$ in (11) be denoted by $\mathcal{F}_V(\tilde{x}) = f_V(S, I)$, $\mathcal{F}_I(\tilde{x}) = f_I(S, I, Q)$, and $\mathcal{F}_C(\tilde{x}) = f_C(S, I)$, respectively. Then the functions have non-negative (Lipschitz) constants $k_{d,V}$, $k_{d,I}$, and $k_{d,C}$ such that

$$ |\mathcal{F}_V(\tilde{x}) - \mathcal{F}_V(\tilde{y})| \le \omega_V k_{d,V} |\tilde{x} - \tilde{y}|, \quad (52a) $$

$$ |\mathcal{F}_I(\tilde{x}) - \mathcal{F}_I(\tilde{y})| \le \omega_I k_{d,I} |\tilde{x} - \tilde{y}|, \quad (52b) $$

$$ |\mathcal{F}_C(\tilde{x}) - \mathcal{F}_C(\tilde{y})| \le \omega_C k_{d,C} |\tilde{x} - \tilde{y}| \quad (52c) $$

for all $\tilde{x}$ and $\tilde{y}$ in $\tilde{\mathcal{Z}}$. Again $k_{d,V}$, $k_{d,I}$, and $k_{d,C}$ are chosen uniformly in the possibly unbounded set $\tilde{\mathcal{Z}}$. To verify the existence of $k_{d,V}$, $k_{d,I}$, and $k_{d,C}$, in addition to the bounded partial derivatives of $H_1$, $H_2$, and $IH_3$ with respect to $\tilde{x}$ on $\tilde{\mathcal{Z}}$, we can use the following for obtaining $k_{d,C}$. We can find a constant $I_a > 0$ such that $S(IH_1(S, I) - H_2(S, I)) > 0$ (and so also $f_C(S, I) = 0$) for all $(S, I)$
occurring as the first two components of tuples in $\mathcal{Z}$ and for which $I < I_a$. This
follows because $I \ln I \to 0$ as $I \to 0^+$. This provides a global Lipschitz constant for
$\min\{0, \omega_C S(IH_1(S, I) - H_2(S, I))\}$ on $\mathcal{Z}$ hence also the required constant $k_{d,C}$.

Using (51), (52), and the fact that $\tilde{x}(t)$ remains in the set $\tilde{\mathcal{Z}}$ for all $t \ge 0$, we can
find positive constants $q_{I1}$ and $q_{I2}$ such that

$$
\begin{align}
|\mathcal{G}_I(t)| & \le q_{I1}\omega_I(|H_3(\mathcal{Q}(t))| + |H_2(S(t), I(t))|) \nonumber \\
& \phantom{\le} \times \left\{ \int_{t-L_I}^t |\dot{\tilde{S}}(\ell)|d\ell + \int_{t-L_I}^t |\tilde{\xi}(\ell)|d\ell + \int_{t-L_I}^t |\dot{\tilde{Q}}(\ell)|d\ell \right\} \tag{53} \\
& \le q_{I2}\omega_I(|\tilde{S}(t)| + |\tilde{\xi}(t)| + |\tilde{Q}(t)|) \nonumber \\
& \phantom{\le} \times \left\{ \int_{t-L_I}^t |\dot{\tilde{S}}(\ell)|d\ell + \int_{t-L_I}^t |\tilde{\xi}(\ell)|d\ell + \int_{t-L_I}^t |\dot{\tilde{Q}}(\ell)|d\ell \right\}. \nonumber
\end{align}
$$

Here and in the sequel, all equalities and inequalities are along all solutions of the
SIQR model (1) with the controller consisting of (11), (13), and (49) with the initial
function $\tilde{x}_0$ valued in $\mathcal{D}_a$ and the perturbations $\epsilon$ valued in $\tilde{\mathcal{P}}$. Using the equation
(15) and the global Lipschitzness of $e^{\tilde{\epsilon}}$ on $\tilde{\mathcal{Z}}$, uniformly in $\tilde{\mathcal{Z}}$, we can construct
constants $q_3 > 0$, $q_4 > 0$ and $q_5 > 0$ such that

$$
\begin{align*}
|\dot{\tilde{S}}(t)| &\le q_3(1+\omega) (|\tilde{S}(t)| + |\tilde{\xi}(t)| + |\tilde{S}(t-L_V)| + |\tilde{\xi}(t-L_V)| \\
&\phantom{\le} + |\tilde{S}(t-L_C)| + |\tilde{\xi}(t-L_C)|) + |\epsilon(t)|, \\
|\dot{\tilde{\xi}}(t)| &\le q_4(1+\omega) (|\tilde{S}(t)| + |\tilde{\xi}(t)| + |\tilde{S}(t-L_I)| + |\tilde{\xi}(t-L_I)| + |\tilde{Q}(t-L_I)|) \\
&\phantom{\le} + |\tilde{S}(t-L_C)| + |\tilde{\xi}(t-L_C)|), \text{ and} \\
|\dot{\tilde{Q}}(t)| &\le q_5(1+\omega) (|\tilde{\xi}(t)| + |\tilde{Q}(t)| + |\tilde{S}(t-L_I)| + |\tilde{\xi}(t-L_I)| + |\tilde{Q}(t-L_I)|),
\end{align*}
$$

where $\omega = \max\{\omega_V, \omega_I, \omega_C\}$. Hence, by extending the domains of the initial func-
tions $\tilde{x}$ to $[-2L, 0]$ by stipulating that they are constant on $[-2L, -L]$ (which does
not affect the solutions of the dynamics, because the solutions only depend on values
of the initial functions on $[-L, 0]$), we can obtain

$$
\begin{align*}
\int_{t-L_I}^{t} |\dot{\tilde{S}}(\ell)| d\ell &\le 6q_3(1+\omega) \int_{t-2L}^{t} |(\tilde{S}(\ell), \tilde{\xi}(\ell))| d\ell + L |\epsilon|_{[t-L_I, t]}, \\
\int_{t-L_I}^{t} |\dot{\tilde{\xi}}(\ell)| d\ell &\le 6q_4(1+\omega) \int_{t-2L}^{t} |(\tilde{S}(\ell), \tilde{\xi}(\ell), \tilde{Q}(\ell))| d\ell, \text{ and} \\
\int_{t-L_I}^{t} |\dot{\tilde{Q}}(\ell)| d\ell &\le 5q_5(1+\omega) \int_{t-2L}^{t} |(\tilde{S}(\ell), \tilde{\xi}(\ell), \tilde{Q}(\ell))| d\ell.
\end{align*}
$$

Hence,

$$
\begin{align}
|\mathcal{G}_I(t)| &\le 3q_{I2}\omega_I|(\tilde{S}(t), \tilde{\xi}(t), \tilde{Q}(t))| \nonumber \\
&\quad \times (1+\omega)(6q_3+6q_4+5q_5) \int_{t-2L}^{t} |(\tilde{S}(\ell), \tilde{\xi}(\ell), \tilde{Q}(\ell))| d\ell \nonumber \\
&\quad + 3q_{I2}\omega_I|(\tilde{S}(t), \tilde{\xi}(t), \tilde{Q}(t))| L |\epsilon|_{[t-L_I,t]} \nonumber \\
&\le 18\omega_I(1+\omega)q_{I2}(q_3+q_4+q_5) \left\{ \sqrt{q_7}|(\tilde{S}(t), \tilde{\xi}(t), \tilde{Q}(t))| \right\} \nonumber \\
&\quad \times \left\{ \frac{1}{\sqrt{q_7}} \int_{t-2L}^{t} |(\tilde{S}(\ell), \tilde{\xi}(\ell), \tilde{Q}(\ell))| d\ell \right\} \nonumber \\
&\quad + \left\{ \sqrt{\frac{2\zeta}{3}} |(\tilde{S}(t), \tilde{\xi}(t), \tilde{Q}(t))| \right\} \left\{ \sqrt{\frac{3}{2\zeta}} 3L\omega_I q_{I2} |\epsilon|_{[t-L_I,t]} \right\} \tag{54}
\end{align}
$$
for any positive constants $q_7$ and $\zeta$. Applying Young's inequality to the pairs of terms in braces in (54) and then Jensen's inequality, we obtain

$$ egin{aligned} |\mathcal{G}_I(t)| & \le 9\omega_I(1+\omega)q_{I2}(q_3+q_4+q_5) \{ q_7 |(\tilde{S}(t), \tilde{\xi}(t), \tilde{Q}(t))|^2 \\ & \quad + \frac{1}{q_7} \left( \int_{t-2L}^t |(\tilde{S}(\ell), \tilde{\xi}(\ell), \tilde{Q}(\ell))| d\ell \right)^2 \} \\ & \quad + \frac{\zeta}{3} |(\tilde{S}(t), \tilde{\xi}(t), \tilde{Q}(t))|^2 + \frac{27}{4\zeta} L^2 \omega_I^2 q_{I2}^2 |\epsilon|_{[t-L_I,t]}^2 \\ & \le 9\omega_I(1+\omega)q_{I2}(q_3+q_4+q_5) \{ q_7 |(\tilde{S}(t), \tilde{\xi}(t), \tilde{Q}(t))|^2 \\ & \quad + \frac{2L}{q_7} \int_{t-2L}^t |(\tilde{S}(\ell), \tilde{\xi}(\ell), \tilde{Q}(\ell))|^2 d\ell \} \\ & \quad + \frac{\zeta}{3} |(\tilde{S}(t), \tilde{\xi}(t), \tilde{Q}(t))|^2 + \frac{27}{4\zeta} L^2 \omega_I^2 q_{I2}^2 |\epsilon|_{[t-L_I,t]}^2. \end{aligned} \quad (55) $$

In the same way, we obtain

$$ egin{aligned} |\mathcal{G}_V(t)| & \le 9\omega_V(1+\omega)q_{V2}(q_3+q_4+q_5) \{ q_7 |(\tilde{S}(t), \tilde{\xi}(t), \tilde{Q}(t))|^2 \\ & \quad + \frac{2L}{q_7} \int_{t-2L}^t |(\tilde{S}(\ell), \tilde{\xi}(\ell), \tilde{Q}(\ell))|^2 d\ell \} \\ & \quad + \frac{\zeta}{3} |(\tilde{S}(t), \tilde{\xi}(t), \tilde{Q}(t))|^2 + \frac{27}{4\zeta} L^2 \omega_V^2 q_{V2}^2 |\epsilon|_{[t-L_I,t]}^2, \text{ and} \end{aligned} \quad (56) $$

$$ egin{aligned} |\mathcal{G}_C(t)| & \le 9\omega_C(1+\omega)q_{C2}(q_3+q_4+q_5) \{ q_7 |(\tilde{S}(t), \tilde{\xi}(t), \tilde{Q}(t))|^2 \\ & \quad + \frac{2L}{q_7} \int_{t-2L}^t |(\tilde{S}(\ell), \tilde{\xi}(\ell), \tilde{Q}(\ell))|^2 d\ell \} \\ & \quad + \frac{\zeta}{3} |(\tilde{S}(t), \tilde{\xi}(t), \tilde{Q}(t))|^2 + \frac{27}{4\zeta} L^2 \omega_C^2 q_{C2}^2 |\epsilon|_{[t-L_I,t]}^2. \end{aligned} \quad (57) $$

Also, for an arbitrary constant $g > 0$, the function $V$ in (20) satisfies

$$ \dot{V}(t) \le -\alpha(V(\tilde{S}(t), \tilde{\xi}(t), \tilde{Q}(t), \tilde{R}(t))) + \mathcal{G}_V(t) + \mathcal{G}_I(t) + \mathcal{G}_C(t) + \sigma(|\epsilon(t)|) \quad (58) $$

along the trajectories of (1) with (11), (13), and (49) for all $t \ge 0$ and all piecewise continuous functions (22). Let

$$ \bar{q} = 9(\omega_V q_{V2} + \omega_I q_{I2} + \omega_C q_{C2})(1+\omega)(q_3+q_4+q_5). $$

Define

$$ V_{\#}(\tilde{S}_t, \tilde{\xi}_t, \tilde{Q}_t, \tilde{R}_t) = V(\tilde{S}(t), \tilde{\xi}(t), \tilde{Q}(t), \tilde{R}(t)) \\ + \frac{2L\bar{q}(1+p)}{q_7} \int_{t-2L}^{t} \int_{\ell}^{t} |(\tilde{S}(r), \tilde{\xi}(r), \tilde{Q}(r))|^2 dr d\ell \quad (59) $$

for $p > 0$. Then the derivative of this functional $V_{\#}$ with respect to $t$ is

$$ \dot{V}_{\#}(\tilde{S}_t, \tilde{\xi}_t, \tilde{Q}_t, \tilde{R}_t) = \dot{V}(\tilde{S}(t), \tilde{\xi}(t), \tilde{Q}(t), \tilde{R}(t)) - \frac{2L\bar{q}(1+p)}{q_7} \int_{t-2L}^{t} |(\tilde{S}(\ell), \tilde{\xi}(\ell), \tilde{Q}(\ell))|^2 d\ell \\ + \frac{4L^2\bar{q}(1+p)}{q_7} |(\tilde{S}(t), \tilde{\xi}(t), \tilde{Q}(t))|^2. $$

Thus, along the trajectories of (1), we have

$$ \begin{aligned} \dot{V}_{\#} \le & -\alpha(V(\tilde{S}(t), \tilde{\xi}(t), \tilde{Q}(t), \tilde{R}(t))) + \bar{q}q_7 |(\tilde{S}(t), \tilde{\xi}(t), \tilde{Q}(t))|^2 \\ & - \frac{2L\bar{q}p}{q_7} \int_{t-2L}^{t} |(\tilde{S}(\ell), \tilde{\xi}(\ell), \tilde{Q}(\ell))|^2 d\ell + \frac{4L^2\bar{q}(1+p)}{q_7} |(\tilde{S}(t), \tilde{\xi}(t), \tilde{Q}(t))|^2 \end{aligned} $$
$$ 
\begin{aligned}
& + \zeta |(\tilde{S}(t), \tilde{\xi}(t), \tilde{Q}(t))|^2 \\
& + \sigma(|\epsilon(t)|) + \frac{27}{4\zeta} L^2 (\omega_V^2 q_{V2}^2 + \omega_I^2 q_{I2}^2 + \omega_C^2 q_{C2}^2) |\epsilon|_{[t-L_I,t]}^2.
\end{aligned} 
\quad (60) $$

Assume that $N > 0$ satisfies

$$ \bar{L}N \geq \frac{2L\bar{q}(1+p)}{q_7} \int_{t-2L}^{t} \int_{\ell}^{t} M^2 dr d\ell \quad (61) $$

with $M = \max_{\tilde{x} \in \tilde{D}_a} |\tilde{x}|$. Let

$$ z = \bar{L}N + \max_{\tilde{x} \in \tilde{D}_a} V(\tilde{x}). \quad (62) $$

Due to (23a), there exists $m_0 > 0$ such that $\alpha(r)/r \ge m_0$ holds for all $r \in [0, z]$. There also exists $m_1 > 0$ such that $V(\tilde{S}, \tilde{\xi}, \tilde{Q}, \tilde{R}) \ge m_1|(\tilde{S}, \tilde{\xi}, \tilde{Q})|^2$ for all $(\tilde{S}, \tilde{\xi}, \tilde{Q}, \tilde{R}) \in \tilde{D}_a$. Indeed, the existence follows by finding a positive definite quadratic lower bound for $e^{\tilde{\xi}} - 1 - \tilde{\xi}$ in $\tilde{\xi}$ (which exists because $\tilde{D}_a$ is bounded), which we use to find a quadratic lower bound for $U$ as a function of $\tilde{S}, e^{\tilde{\xi}} - 1, \tilde{\xi}, \tilde{Q}, \tilde{R}$ with positive coefficients. Hence, if

$$ V(\tilde{x}(t)) \le z, \quad (63) $$

then we obtain

$$ 
\begin{aligned}
\dot{V}_{\#} & \le -b_1 V(\tilde{S}_t, \tilde{\xi}_t, \tilde{Q}_t, \tilde{R}_t) \\
& - b_2 \frac{2L\bar{q}(1+p)}{q_7} \int_{t-2L}^{t} \int_{\ell}^{t} |(\tilde{S}(r), \tilde{\xi}(r), \tilde{Q}(r))|^2 dr d\ell \\
& + \sigma(|\epsilon(t)|) + \frac{27}{4\zeta} L^2 (\omega_V^2 q_{V2}^2 + \omega_I^2 q_{I2}^2 + \omega_C^2 q_{C2}^2) |\epsilon|_{[t-L_I,t]}^2 \\
& \le -b_3 V_{\#}(\tilde{S}_t, \tilde{\xi}_t, \tilde{Q}_t, \tilde{R}_t) \\
& + \sigma(|\epsilon(t)|) + \frac{27}{4\zeta} L^2 \omega^2 (q_{V2}^2 + q_{I2}^2 + q_{C2}^2) |\epsilon|_{[t-L_I,t]}^2,
\end{aligned} 
\quad (64) $$

where

$$ b_1 = m_0 - \frac{\zeta}{m_1} - \frac{\bar{q}}{m_1} \left( q_7 + \frac{4L^2(1+p)}{q_7} \right), \quad b_2 = \frac{p}{2L(1+p)}, \\ \text{and } b_3 = \min\{b_1, b_2\}. $$

Let $\hat{\epsilon} > 0$ be such that

$$ zb_3 = \sigma(\hat{\epsilon}) + \frac{27}{4\zeta} L^2 \omega^2 (q_{V2}^2 + q_{I2}^2 + q_{C2}^2) \hat{\epsilon}^2. $$

Such a $\hat{\epsilon} > 0$ exists if $b_1 > 0$. Suppose that

$$ b_1 > 0 \quad (65) $$

holds. Define $\bar{\epsilon} = \min\{\psi_*/4, \hat{\epsilon}\}$. Then the definitions (59) and (62) with the assumption (61) yield $V_{\#} \le z$ at $t = 0$, and by virtue of (64), $V_{\#} \le z$ continues to hold for all $t > 0$ as long as $\epsilon$ belongs to $(-\min\{B, \bar{\epsilon}\}, \bar{\epsilon})$, and $\tilde{x}_0$ is valued in $\tilde{D}_a$. Since (59) implies $V \le V_{\#}$, the assumption (63) is met for all $t \ge 0$. Furthermore, for the perturbation set $(-\min\{B, \bar{\epsilon}\}, \bar{\epsilon})$ and for all initial functions $\tilde{x}_0$ valued in $\tilde{D}_a$,
$V_\#$ is an ISS Lyapunov-Krasovskii functional for the controlled SIQR model controlled with (49) [29]. We now specialize the preceding analysis to the case where $q_7 = m_0m_1/(2\bar{q})$. Then there exists $\zeta > 0$ such that (65) holds if and only if

$$m_0^2 m_1^2 > 16 L^2 \bar{q}^2 (1 + p). \qquad (66)$$

The inequality (61) holds if and only if

$$\bar{L}N \geq \frac{4L\bar{q}^2(1+p)}{m_0m_1} 2L^2M^2. \qquad (67)$$

Therefore, due to $L \le \bar{L}$, the pair of (66) and (67) provides the constant $A > 0$ that is required for (50) to hold. □

For each choice of $D_a$, condition (50) is always met if $\max\{\omega_V, \omega_I, \omega_C\} = 0$ or $\max\{L_V, L_I, L_C\} = 0$. Thus, Theorem 5.1 recovers Corollary 3.1 as a special case. Importantly, the delay robustness established by Theorem 5.1 is not local since $D_a$ can be arbitrarily large in the state space $D$. The larger the domain $D_a$ of interest is, the smaller $A$ can become. Nevertheless, there always exists $A > 0$ that is independent of $\omega_V, \omega_I, \omega_C, L_V, L_I$, and $L_C$. We can also remove the restriction $[-\bar{\epsilon}, \bar{\epsilon}]$ on the perturbation $\epsilon$, as follows:

**Theorem 5.2.** Let $D_a$ be any compact subset of $D$ containing $X_\star$. Then for any constants $c > 0$, $c_\diamond \in (0, 2\bar{c}_\diamond)$, and $\beta \in [0, \hat{\beta}]$, we can construct a positive constant $A$ such that if nonnegative constants $\omega_V, \omega_I, \omega_C, L_V, L_I$, and $L_C$ satisfy (50), then the dynamics for the error variable (31) for the SIQR model (1) with the controller consisting of (11), (13), and (49) are ISS on $D_a$ for the perturbation set $P$.

*Proof.* Theorem 5.1 ensures the existence of $A > 0$ in (50) implying the ISS property for all initial functions $\tilde{x}_0$ that are valued in $\tilde{D}_a$ and the perturbation set $(-\min\{B, \bar{\epsilon}\}, \bar{\epsilon})$. On the other hand, since $\dot{S} + \dot{I} + \dot{Q} + \dot{R} = B + \epsilon - \mu(S + I + Q + R)$ is satisfied by (1) independently of the control inputs, we have (34) for all $X(0) = (S(0), I(0), Q(0), R(0)) \in [0, \infty)^4$ and all $\epsilon$'s that are valued in $P$. Therefore, noting that $|\tilde{x}_0(0)| \le |\tilde{x}_0|_{[-L,0]}$, we can apply the argument used in the proof of Theorem 3.2 to the SIQR model (1) with the delayed inputs for $D_a$ and $\epsilon \in P$. □

**Remark 5.1.** We found that our controls continue to achieve our ISS objectives even when the condition (50) is not satisfied. For instance, using the same data that we used in our delay simulations from Section 4.5, we found the same qualitative convergence properties when the delays were 7 days instead of 2 days; see Fig. 9, which should be compared with Fig. 8. This indicates a desirable robustness property of our work with respect to violations of our condition (50), and it indicates that our controls (11) are useful even if the delays are larger than those that are allowed by the preceding two theorems.

**6. Practical interpretations and remarks.** Control of an infectious disease is the management of disease spread in a way a society can sustain, as discussed in [14, Chapter 8]. Roles and choices of models for investigating control strategies are also clarified there. For vaccination, this paper uses a standard model from [14, Chapter 8]. The vaccination in (1) is the rate $\rho$ with respect to the population $S$ of susceptible individuals since vaccination removes susceptible individuals and moves them to the group $R$ of individuals who gained immunity. The rate can be considered to be proportional to the amount of medical resources prepared for
<figure><img src="image_10.png" /><figcaption>FIGURE 9. Populations of (1) controlled by (11) in (3) with 7 day delay.</figcaption></figure>

vaccination taking into account the efficacy of vaccines. The resources by themselves define neither the nominal nor the maximum numbers of vaccinated individuals or vaccine shots per day since the numbers must depend on the population of individuals who are neither vaccinated nor infected, i.e., vaccination resources need to encounter susceptible individuals. The number of vaccine shots must decrease to zero as $S$ tends to zero. Mass vaccination ordered by a government is unnecessary, i.e., $u_V = 0$ if susceptible and infected populations are at their desired equilibrium, i.e., $\tilde{S} = \tilde{I} = 0$. When $S$ and $I$ are larger than the equilibrium, the vaccination needs to be increased, and the decision puts more weight on $\tilde{S}$ than $\tilde{I}$. This weighting is consistent of the idea of herd immunity [5]. The control law given by (3a) and (11a) represents such a policy.

To let a population become immune, vaccination is far slower than disease transmission due to the essential difference in their growth rates. In disease transmission, transmitters are doubling, while in vaccination, only professionals administer vaccines. For a human disease, no society would encourage individuals to contract the disease. If eliminating a disease from a region eventually is the only objective, increasing vaccination as much as possible is the best [14]. The basic production number shows that a disease is eliminated if the vaccination rate is larger than a threshold [5]. However, even when vaccination is well above the threshold, it always takes time to eliminate an infectious disease [14]. The elimination is a steady-state property. Intensive vaccination does not mean the removal of infection peaks unless vaccines are administrated well in advance. The effects of isolation and contact reduction are quicker and alter the transient response, although they cannot be used for infection reduction at a steady state. This motivates the idea of managing disease spread without solely relying on vaccination [14, Chapter 8]. This paper employs ISS since it is not a steady-state property. Utilizing ISS, this paper proposes a method to design the simultaneous use of three controls to reduce the infection peaks in a coordinated way to reduce the burden of the medical resources for the vaccination and the treatment of patients in the presence of immigration uncertainty. To illustrate this point, the simulations in Section 4 use a small gain parameter $\omega_V$ for vaccination. The parameters $\omega_V$, $\omega_I$, and $\omega_C$ in (11) are tuning parameters that can be selected taking into account the cost of the three strategies. Figure 10 plots a simulation result where $\omega_V$ is increased by 3 times, while the other parameters are the same as those for Fig. 1b. The infection peak is pressed down
further, although the number of vaccinated people per day increases significantly as shown in Fig. 11.

Since the proposed controller consisting of the three inputs (11) is constructed from a Lyapunov function, the controller enjoys an optimality property, which is known as inverse optimality offered by a control Lyapunov function [6, 32]. The optimal action occurs whenever the control inputs do not hit their range limits. The inverse optimality gives a cost function a posteriori. The ISS property achieved by the proposed controller provides such optimality. We do not have to select a cost function a priori, which is not at all an obvious task. Alternative control approaches for COVID-19 models include optimal control, using predetermined cost functions. For instance, the work [3] uses model predictive control to minimize the costs of mitigation strategies while ensuring that the capacity of a regional healthcare system network is not exceeded, [4] establishes the existence of optimal controls under transmission and treatment uncertainty, [20] derives optimal controls using social distancing as a control policy, the data-driven optimal control approach in [21] uses learning methods to estimate model parameters to forecast the evolution of an outbreak over short time periods and to provide scheduled controls, [22] uses an age-structured population compartmental finite-dimensional optimal control model to optimize vaccination policy to minimize deaths, the work [24] finds optimal strategies as combinations of implementing multiple non-pharmaceutical interventions, [25] derives optimal social distancing strategies using on-off social isolation strategies, [27] uses optimal control to optimize timings for two-dose vaccine roll outs, [30] uses optimal control methods to study trade-offs between lives saved versus reduced time under control, and [33] uses optimal control to minimize the epidemic final size while keeping the infected peak prevalence controlled at each time. The preceding works are notable, because policy makers generally have certain control objectives in mind that are not necessarily expressed in a mathematically rigorous way. However, our analysis is a departure from, and adds complementary value relative to, the preceding optimal control works. This is because our novel global strict Lyapunov function enables us to arrive at feedback controls quantifying the effects of input delays and immigration uncertainty, using ISS and a new family of feedback controls for three controlled quantities (namely, isolation, contact regulation, and vaccination), while also enabling us to compare the effects of different feedback control parameters using actual data from the COVID-19 pandemic. These important features are outside the scope of the preceding optimal control approaches that do not use feedback control. Furthermore, the proposed feedback control laws in the closed form (11) are free from the cost of real-time updates relying on numerical computations and storage of data.

The proposed controller (11) generates control inputs that decrease the deviation from the desired equilibrium all the time. The control law (11b) of isolation involves the population Q of isolated individuals since a large isolation rate directly implies an increase of isolated individuals that cannot exceed its resources. By contrast, the vaccination (11a) and the contact regulation (11c) do not directly increase the number of isolated individuals. Therefore, the control laws of the vaccination (11a) and the contact regulation (11c) use only the populations of the susceptible and infected individuals, which are directly relevant to the disease transmission. In this way, the proposed controller (11) mitigates the burden on societies by reducing contact and isolating mild cases which do not spend medical resources.
Duration of immunity varies with types of diseases [14]. This paper deals with the management of diseases with long-lasting immunity. Variant strains of a disease are sometimes treated better as a different disease. Sending recovered individuals back to the susceptible group in a disease model is insufficient in situations where the transmission rate, the efficacy of vaccination, the duration of immune protection and other rates alter. This paper uses data from the COVID-19 pandemic for the simulations since the usefulness of the SIQR model with long-lasting immunity is demonstrated in [26]. Nevertheless, extending the proposed mathematical idea to cover waning immunity is an important direction for future research. For mathematical simplicity, this paper allows control inputs to be updated all the time. The delay robustness is established in this paper for practical control implementation which does not update the control at the ideal timing. Figure 12 plots the simulation result in which the values of the control inputs (11) are updated only at every two weeks. The result is similar to the one with the continuously updated control inputs shown in Figs. 1b and 6b using the same controller parameters. It indicates practical usefulness of the proposed controller. Providing theoretical guidelines and guarantees for (approximately) piecewise constant implementation of the control inputs is also a practically important topic of future study.

<figure><img src="image_11.png" /><figcaption>FIGURE 10. Populations of (1) controlled by (11) in (3) with $\omega_V = 0.000045$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIGURE 11. Semi-log plot with number of vaccinated individuals per day for two different $\omega_V$'s.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>FIGURE 12. Effects of using piecewise constant control which updates values of $u_V$, $u_I$, and $u_C$ in (11) used in the controls (3) every 14 days and keeps the control values constant between updates.</figcaption></figure>

**7. Conclusion.** By adding isolation and contact regulation to insufficient mass vaccination, this paper demonstrated the usefulness of an ISS Lyapunov function for feedback design for global disease control. The proposed feedback laws of the three control inputs can reduce the peak infection levels and expedite the convergence to the equilibrium. To remove the inconvenience of Lyapunov functions that have been used in the literature, this paper focused on a new Lyapunov function and its time-derivative that contain all the population variables on the entire state space. Such construction of the Lyapunov function allows its gradient to yield the three control laws directly in a uniform manner. The Lyapunov function constructed in this paper has free parameters. The coefficient weighting the isolated population is useful for achieving a better balance between the three control inputs to reduce the burden on the resources of isolation facilities and medical treatments. Since the constructed Lyapunov function is also an ISS Lyapunov function, the control guarantees robustness for immigration perturbation. Furthermore, this paper has demonstrated robustness of the closed-loop for delays in control inputs. In future work, we will study control laws providing guarantees for longer time delays by delay compensation based on exact predictors, chain predictors [1, 2], or other dynamic extensions [31]. Investigating the possibility of incorporating new free parameters into the ISS Lyapunov function would also be important to pursue control laws that can reduce the infection levels and the burden on resources further.

## REFERENCES

[1] I. Bhogaraju, M. Farasat, M. Malisoff and M. Krstic, Sequential predictors for delay-compensating feedback stabilization of bilinear systems with uncertainties, Systems Control Lett., **152** (2021), 104933, 9 pp.

[2] F. Cacace and A. Germani, Output feedback control of linear systems with input, state and output delays by chains of predictors, Automatica J. IFAC, **85** (2017), 455-461.

[3] R. Carli, G. Cavone, N. Epicoco, P. Scarabaggio and M. Dotoli, Model predictive control to mitigate the COVID-19 outbreak in a multi-region scenario, Annual Reviews in Control, **50** (2020), 373-393.

[4] Y. Ding and H. Schellhorn, Optimal control of the SIR model with constrained policy, with an application to COVID-19, Mathematical Biosciences, **344** (2022), 108758, 15 pp.

[5] P. E. M. Fine, Herd immunity: History, theory, practice, Epidemiologic Reviews, **15** (1993), 265-302.
[6] R. A. Freeman and P. V. Kokotović, *Robust Nonlinear Control Design: State-Space and Lyapunov Techniques*, Birkhäuser, Boston, 1996.

[7] A. Gumel, S. Ruan, T. Day, J. Watmough, F. Brauer, P. van den Driessche, D. Gabrielson, C. Bowman, M. Alexander, S. Ardal, J. Wu and B. Sahai, *Modelling strategies for controlling SARS outbreaks*, *Proceedings of the Royal Society of London Series B*, **271** (2004), 2223-2232.

[8] H. Hethcote, M. Zhien and L. Shengbing, *Effects of quarantine in six endemic models for infectious diseases*, *Mathematical Biosciences*, **180** (2002), 141-160.

[9] H. Ito, *Input-to-state stability and Lyapunov functions with explicit domains for SIR model of infectious diseases*, *Discrete Contin. Dyn. Syst. Ser. B*, **26** (2021), 5171-5196.

[10] H. Ito, *A construction of strict Lyapunov functions for a bilinear balancing model*, *IFAC-PapersOnLine*, **54** (2021), 161-166.

[11] H. Ito, *Vaccination with input-to-state stability for SIR model of epidemics*, *Proceedings of the 60th IEEE Conference on Decision and Control*, (2021), 2812-2817.

[12] H. Ito, M. Malisoff, and F. Mazenc, *Strict Lyapunov functions and feedback controls for SIR models with quarantine and vaccination*, *Discrete and Continuous Dynamical Systems-B*, **27** (2022), 6969-6988.

[13] Japanese Ministry of Health, Labour and Welfare, https://www.mhlw.go.jp/english/, (2021).

[14] M. J. Keeling and P. Rohani, *Modeling Infectious Diseases in Humans and Animals*, Princeton University Press, Princeton, 2008.

[15] W. Kermack and A. McKendrick, *A contribution to the mathematical theory of epidemics*, *Proceedings of the Royal Academy of London Series A*, **115** (1927), 700-721.

[16] H. Khalil, *Nonlinear Systems, Third Edition*, Prentice-Hall, Englewood Cliffs, 2002.

[17] A. Korobeinikov, *Lyapunov functions and global properties for SEIR and SEIS epidemic models*, *Mathematical Medicine and Biology*, **21** (2004), 75-83.

[18] A. Korobeinikov, *Lyapunov functions and global stability for SIR and SIRS epidemiological models with non-linear transmission*, *Bulletin of Mathematical Biology*, **30** (2006), 615-626.

[19] A. Korobeinikov and G. Wake, *Lyapunov functions and global stability for SIR, SIRS, and SIS epidemiological models*, *Applied Mathematics Letters*, **15** (2002), 955-960.

[20] R. Kovacevic, N. Stilianakis, and V. Veliov, *A distributed optimal control model applied to COVID-19 pandemic*, *SIAM Journal on Control and Optimization*, **60** (2022), 221-245.

[21] H. Liu and X. Tian, *Data-driven optimal control of a SEIR model for COVID-19*, *Communications on Pure and Applied Analysis*, (2022), to appear.

[22] Q. Luo, R. Weightman, S. McQuade, M. Diaz, E. Trelat, W. Barbour, D. Work, S. Samaranayake and B. Piccoli, *Optimization of vaccination for COVID-19 in the midst of a pandemic*, *Networks and Heterogeneous Media*, **17** (2022), 443-466.

[23] M. Malisoff and F. Mazenc, *Constructions of Strict Lyapunov Functions*, Springer-Verlag, London, Ltd., London, 2009.

[24] S. T. McQuade, R. Weightman, N. Merrill, A. Yadav, E. Trelat, S. Allred and B. Piccoli, *Control of COVID-19 outbreak using an extended SEIR model*, *Mathematical Models and Methods in Applied Sciences*, **31** (2021), 2399-2424.

[25] M. M. Morato, S. B. Bastos, D. O. Cajueiro and J. E. Normey-Rico, *An optimal predictive control strategy for COVID-19 (SARS-CoV-2) social distancing policies in Brazil*, *Annual Reviews in Control*, **50** (2020), 417-431.

[26] T. Odagaki, *Analysis of the outbreak of COVID-19 in Japan by SIQR model*, *Infectious Disease Modelling*, **5** (2020), 691-698.

[27] R. Parino, L. Zino, G. Calafiore and A. Rizzo, *A model predictive control approach to optimally devise a two-dose vaccination rollout: A case study on COVID-19 in Italy*, *International Journal of Robust and Nonlinear Control*, (2021), to appear.

[28] P. Pepe, *A nonlinear version of Halanay's inequality for the uniform convergence to the origin*, *Math. Control Relat. Fields*, **12** (2022), 789-811.

[29] P. Pepe and Z.-P. Jiang, *A Lyapunov-Krasovskii methodology for ISS and iISS of time-delay systems*, *Systems Control Lett.*, **55** (2006), 1006-1014.

[30] T. Perkins and G. España, *Optimal control of the COVID-19 pandemic with non-pharmaceutical interventions*, *Bull. Math. Biol.*, **82** (2020), Paper No. 118, 24 pp.

[31] A. Selivanov and E. Fridman, *Predictor-based networked control under uncertain transmission delays*, *Automatica J. IFAC*, **70** (2016), 101-108.

[32] R. Sepulchre, M. Janković and P. Kokotovic, *Constructive Nonlinear Control*, Communications and Control Engineering Series, Springer-Verlag, Berlin, 1997.
[33] J. Sereno, A. Anderson, A. Ferramosca, E. Hernandez-Vargas and A. Gonzalez, *Minimizing the epidemic final size while containing the infected peak prevalence in SIR systems*, *Automatica J. IFAC*, **144** (2022), Paper No. 110496.

[34] E. D. Sontag, *Smooth stabilization implies coprime factorization*, *IEEE Trans. Automat. Control*, **34** (1989), 435-443.

[35] E. D. Sontag and Y. Wang, *On characterizations of input-to-state stability property*, *Systems Control Lett.*, **24** (1995), 351-359.

[36] Statistics Bureau of Japan, https://www.stat.go.jp/english/data/jinsui/tsuki/index.ht, (2021).

[37] S. Tang, Y. Xiao and D. Clancy, *New modelling approach concerning integrated disease control and cost-effectivity*, *Nonlinear Analysis*, **63** (2005), 439-471.

[38] C. Tian, Q. Zhang and L. Zhang, *Global stability in a networked SIR epidemic model*, *Applied Mathematics Letters*, **107** (2020), 106444, 6 pp.

[39] W. W. C. Topley and G. S. Wilson, *The spread of bacterial infection. The problem of herd-immunity*, *Journal of Hygiene*, **21** (1923), 243-249.

[40] G. G. Walter and M. Contreras, *Compartmental Modeling with Networks*, Modeling and Simulation in Science, Engineering and Technology. Birkhäuser Boston, Inc., Boston, MA, 1999.

Received January 2022; 1st revision June 2022; 2nd revision September 2022; early access September 2022.