Optimal Control for an Epidemic Model of COVID-19 with Time-Varying Parameters

Abstract

The coronavirus disease 2019 (COVID-19) pandemic disrupted public health and economies worldwide. In this paper, we investigate an optimal control problem to simultaneously minimize the epidemic size and control costs associated with intervention strategies based on official data. Considering people with undetected infections, we establish a control system of COVID-19 with time-varying parameters. To estimate these parameters, a parameter identification scheme is adopted and a mixed algorithm is constructed. Moreover, we present an optimal control problem with two objectives that involve the newly increased number of infected individuals and the control costs. A numerical scheme is conducted, simulating the epidemic data pertaining to Shanghai during the period of 2022, caused by the Omicron variant. Coefficient combinations of the objectives are obtained, and the optimal control measures for different infection peaks are indicated. The numerical results suggest that the identification variables obtained by using the constructed mixed algorithm to solve the parameter identification problem are feasible. Optimal control measures for different epidemic peaks can serve as references for decision-makers.

1. Introduction

The coronavirus disease 2019 (COVID-19) pandemic disrupted health systems and economies throughout the world [1,2]. During 2020, about one-third of the US population was infected, and the United States experienced the highest numbers of reported cases and deaths [1]. COVID-19 has led to delays in essential medical services in Europe, including treatment for cardiovascular diseases [3]. Since its arrival, disruptions to childhood vaccine delivery have further jeopardised childhood vaccination efforts [4]. Long COVID has resulted in individuals losing nearly half of their income [5]. Overall, the COVID-19 pandemic had an unprecedented effect on lives, livelihoods, and economies around the world [6]. Therefore, it is crucial to further investigate modeling and control measures of COVID-19 on total health expenditure and development assistance for health, which can contribute to limiting the transmission of future high-consequence pathogens and future pandemic preparedness [7,8].

Since the outbreak of the COVID-19 epidemic, a large number of models have been proposed to describe, simulate, and forecast its transmission rule. Some models are individual-based [9] or multi-scale [10], some models are based on cells and viruses [11] or stochastic processes [12,13], and others are developments of classical SIR and SEIR models [14,15,16,17,18,19,20]. On the other hand, there are a lot of studies about control measures. Vaccination [21,22] and testing [23,24] are pharmaceutical interventions. In the absence of effective vaccines and treatment, non-pharmaceutical intervention measures are essential [25]. For example, using face masks [26], quarantining [27], and maintaining social distancing [28,29,30] are all effective methods. Albi et al. [31] put forward an optimal control problem of a socially structured epidemic model in the presence of uncertain data to reduce the spread of epidemics by applying non-pharmaceutical intervention measures. Sereno et al. [32] investigate how to minimize an epidemic’s final scale while keeping the infected peak prevalence controlled at any time. Kuddus et al. [33] perform a mathematical analysis of COVID-19 and find out the most effective control measure with the cost-benefit of the health economy. Zhang et al. [34] establish a stochastic SAIR epidemic model to investigate the dynamics and control strategies of COVID-19. Reza et al. [35] delve into the numerical solutions of a six-compartment fractional model with a Caputo derivative. They report on the sensitivity of the most critical parameter and its influence on COVID-19 dynamics, along with its impact on the basic reproduction number.

In this paper, we aim to simultaneously minimize the newly increased number of infected individuals and the control costs associated with the intervention strategies based on official data. Different from the works mentioned above, the original contributions of the content studied in this paper are as follows:

(i)

Undetected infections are considered.

(ii)

The parameter identification method to estimate the time-varying transmission rate functions in the epidemic model is employed.

(iii)

The optimal control measures for different infection peaks are shown.

We establish a control epidemic system of COVID-19 that involves individuals with undetected infections and time-varying parameters. To determine these time-varying parameters, we adopt the parameter identification scheme and construct a mixed algorithm. Moreover, we put forward an optimal control problem that includes two objectives: one is the newly increased number of infected people, and the other is the control cost associated with intervention strategies. Coefficient combinations of these objectives are indicated to attain different infection peaks. A numerical scheme is carried out using the official data in Shanghai 2022.

The structure of this paper is outlined as follows. In Section 2, the control epidemic system of COVID-19 with time-varying parameters is established, the control reproduction number is derived, and the parameter identification method and the mixed algorithm are constructed. Section 3 puts forward a bi-objective optimal control model and the characteristics of optimal controls. Numerical simulations are performed with the official data of Shanghai 2022 in Section 4. Section 5 concludes the paper.

2. COVID-19 Control System

Let ${\mathbb{Z}}^{+}$ be the set of positive integers. For any $n\in {\mathbb{Z}}^{+},$ define ${\mathbb{Z}}_n=\{1,2,\cdots ,n\}$ as the set of natural numbers from 1 to $n.$

People are divided into eight groups: susceptible individuals ($S\left(t\right)$ for short), infected individuals with infectiousness in the incubation period ($P\left(t\right)$ for short), asymptomatic individuals ($A\left(t\right)$ for short), symptomatic individuals ($I\left(t\right)$ for short), tested asymptomatic individuals ($T_A\left(t\right)$ for short), tested symptomatic individuals ($T_I\left(t\right)$ for short), undetected removed individuals ($R_U\left(t\right)$ for short), and reported removed individuals ($R_R\left(t\right)$ for short). Asymptomatic individuals have no symptoms and do not seek medical attention, so some of them may not be detected and are self-healing. Let $S\left(t\right)+P\left(t\right)+A\left(t\right)+I\left(t\right)+T_A\left(t\right)+T_I\left(t\right)+R_U\left(t\right)+R_R\left(t\right)\triangleq N\left(t\right)$ be the total number of people, t be time, and T be the terminal time.

Since the tested asymptomatic and tested symptomatic individuals are isolated, they are not infectious. Here, we only consider that the infected individuals with infectiousness in the incubation period, the asymptomatic individuals, and the symptomatic individuals are infectious; let ${\beta }_P\left(t\right),$ ${\beta }_A\left(t\right)$, and ${\beta }_I\left(t\right)$ be the infection rates of the infectious people $P\left(t\right),$ $A\left(t\right)$, and $I\left(t\right)$ at time t, respectively. $\frac{1}{\alpha }$ is the period from the beginning of being infectious to the time at which symptoms first appear. Although asymptomatic individuals never have symptoms, we treat asymptomatic individuals as a special case of symptomatic individuals for simplicity. The ratios of the symptomatic individuals and the asymptomatic individuals at time t are $\sigma \left(t\right)$ and $1-\sigma \left(t\right).$ $\lambda \left(t\right)$ is the tested rate; $\mu \left(t\right)$ and ${\mu }_T\left(t\right)$ are the transmission rates from the asymptomatic to the symptomatic people and from the tested asymptomatic to the tested symptomatic people at time $t,$ respectively. ${\gamma }_U\left(t\right),$ ${\gamma }_A\left(t\right)$, and ${\gamma }_I\left(t\right)$ are the removed rates of the asymptomatic, the tested asymptomatic, and the tested symptomatic individuals at time $t.$

Let $\mathbf{u}\left(t\right)$ be the control measures, such as the use of face masks, proper hand-washing, maintaining physical distancing at time t, and so on.

Based on the actual spreading process of COVID-19, we establish a control epidemic model:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{S}^{\prime }\left(t\right)=-S\left(t\right)[1-\mathbf{u}\left(t\right)]\frac{{\beta }_P\left(t\right)P\left(t\right)+{\beta }_A\left(t\right)A\left(t\right)+{\beta }_I\left(t\right)I\left(t\right)}{N\left(t\right)},\hfill \\ P^{\prime }\left(t\right)=S\left(t\right)[1-\mathbf{u}\left(t\right)]\frac{{\beta }_P\left(t\right)P\left(t\right)+{\beta }_A\left(t\right)A\left(t\right)+{\beta }_I\left(t\right)I\left(t\right)}{N\left(t\right)}-\alpha P\left(t\right),\hfill \\ A^{\prime }\left(t\right)=\alpha \left[1-\sigma \left(t\right)\right]P\left(t\right)-\left[\mu \left(t\right)+{\gamma }_U\left(t\right)+\lambda \left(t\right)\right]A\left(t\right),\hfill \\ I^{\prime }\left(t\right)=\alpha \sigma \left(t\right)P\left(t\right)+\mu \left(t\right)A\left(t\right)-\lambda \left(t\right)I\left(t\right),\hfill \\ T_A^{\prime }\left(t\right)=\lambda \left(t\right)A\left(t\right)-\left[{\mu }_T\left(t\right)+{\gamma }_A\left(t\right)\right]T_A\left(t\right),\hfill \\ T_I^{\prime }\left(t\right)=\lambda \left(t\right)I\left(t\right)-{\gamma }_I\left(t\right)T_I\left(t\right)+{\mu }_T\left(t\right)T_A\left(t\right),\hfill \\ R_U^{\prime }\left(t\right)={\gamma }_U\left(t\right)A\left(t\right),\hfill \\ R_R^{\prime }\left(t\right)={\gamma }_A\left(t\right)T_A\left(t\right)+{\gamma }_I\left(t\right)T_I\left(t\right).\hfill \end{array}\right.\end{array}$$

(1)

Considering the biological sense, all variables are non-negative and the initial conditions are

$$\begin{array}{c}\hfill S\left(0\right)\ge 0,P\left(0\right)\ge 0,A\left(0\right)\ge 0,I\left(0\right)\ge 0,T_A\left(0\right)\ge 0,T_I\left(0\right)\ge 0,R_U\left(0\right)\ge 0,R_R\left(0\right)\ge 0.\end{array}$$

(2)

The model (1) and the initial conditions (2) are called the control system of COVID-19, the diagram of which is indicated in Figure 1.

Let

$$\begin{array}{c}\hfill \mathbf{x}\left(t\right)={(S\left(t\right),P\left(t\right),A\left(t\right),I\left(t\right),T_A\left(t\right),T_I\left(t\right),R_U\left(t\right),R_R\left(t\right))}^T\triangleq {({\mathbf{x}}_1\left(t\right),\cdots ,{\mathbf{x}}_8\left(t\right))}^T.\end{array}$$

According to refs. [36,37], we can obtain the existence, uniqueness, and non-negativeness of solution $\mathbf{x}\left(t\right)$ of the COVID-19 system, whose set is denoted as $\mathbb{X}.$

2.1. Control Reproduction Number

Let $P\left(t\right)=0,A\left(t\right)=0,I\left(t\right)=0,T_A\left(t\right)=0,T_I\left(t\right)=0$; we can obtain the disease-free equilibria $\left(N\right(0),0,0,0,0,0,0,0).$ And we calculate the control reproduction number according to ref. [38]:

$$\begin{array}{c}\hfill R_c\left(\mathbf{u}\right)=\frac{{\beta }_P\left(t\right)[1-\mathbf{u}\left(t\right)]}{\alpha }+\frac{{\beta }_A\left(t\right)[1-\mathbf{u}\left(t\right)][1-\sigma \left(t\right)]}{\mu \left(t\right)+{\gamma }_U\left(t\right)+\lambda \left(t\right)}\\ \hfill +\frac{{\beta }_I\left(t\right)[1-\mathbf{u}\left(t\right)]\{\mu \left(t\right)+\sigma \left(t\right)[{\gamma }_U\left(t\right)+\lambda \left(t\right)]\}}{\lambda \left(t\right)\left[\mu \right(t)+\eta (t)+\lambda (t\left)\right]}.\end{array}$$

For $\mathbf{u}=0,$ the control reproduction number $R_c\left(\mathbf{u}\right)$ transforms into the basic reproduction number:

$$\begin{array}{c}\hfill R_0\left(t\right)=\frac{{\beta }_P\left(t\right)}{\alpha }+\frac{{\beta }_A\left(t\right)[1-\sigma \left(t\right)]}{\mu \left(t\right)+{\gamma }_U\left(t\right)+\lambda \left(t\right)}+\frac{{\beta }_I\left(t\right)\{\mu \left(t\right)+\sigma \left(t\right)[{\gamma }_U\left(t\right)+\lambda \left(t\right)]\}}{\lambda \left(t\right)[\mu \left(t\right)+{\gamma }_U\left(t\right)+\lambda \left(t\right)]}.\end{array}$$

where $R_c\left(\mathbf{u}\right)$ and $R_0\left(t\right)$ are the functions of time. We will show their change trend on the numerical aspect later.

2.2. Parameter Identification

In the control system of COVID-19, the variables ${\beta }_P\left(t\right),$ ${\beta }_A\left(t\right),$ ${\beta }_I\left(t\right),$ ${\gamma }_U\left(t\right),$ ${\gamma }_A\left(t\right),$ ${\gamma }_I\left(t\right),$ $\sigma \left(t\right),$ $\mu \left(t\right)$ and $\lambda \left(t\right)$ are unknown. To estimate them, the control variable $\mathbf{u}\left(t\right)$ is set to zero in the control system. Let

$$\begin{array}{ccc}& & \mathbf{p}\left(t\right)={\left({\beta }_P\left(t\right),{\beta }_A\left(t\right),{\beta }_I\left(t\right),{\gamma }_U\left(t\right),{\gamma }_A\left(t\right),{\gamma }_I\left(t\right),\sigma \left(t\right),\mu \left(t\right),\lambda \left(t\right)\right)}^T.\hfill \end{array}$$

Next, we use the parameter identification method to determine $\mathbf{p}\left(t\right).$

We aim to make the COVID-19 system without control describe the spreading process better. For this, we put forward the following objective functional:

$$\begin{array}{c}\hfill {\mathcal{L}}_1\left(\mathbf{p}\right)={\int }_0^T\left[\frac{\left|T_A(\mathbf{p};t)-A_{data}\left(t\right)\right|}{A_{data}\left(t\right)}+\frac{\left|T_I(\mathbf{p};t)-I_{data}\left(t\right)\right|}{I_{data}\left(t\right)}+\frac{\left|R_R(\mathbf{p};t)-R_{data}\left(t\right)\right|}{R_{data}\left(t\right)}\right]dt,\end{array}$$

(3)

which means the accumulative relative errors between the output of the control system and the official data. $A_{data}\left(t\right),$ $I_{data}\left(t\right),$ and $R_{data}\left(t\right)$ denote the official data of the asymptomatic, symptomatic, and removed individuals at time $t.$ We introduce

$$\begin{array}{c}\hfill {\mathcal{L}}_2\left(\mathbf{p}\right)=\underset{0}{\overset{T}{\bigvee }}{\beta }_P\left(t\right)+\underset{0}{\overset{T}{\bigvee }}{\beta }_A\left(t\right)+\underset{0}{\overset{T}{\bigvee }}{\beta }_I\left(t\right),\end{array}$$

(4)

to describe the total variation in the infection rates, the definition of which is shown in refs. [39,40]. Thus, the objective functional is

$$\begin{array}{c}\hfill \mathcal{L}\left(\mathbf{p}\right)={\mathcal{L}}_1\left(\mathbf{p}\right)+a{\mathcal{L}}_2\left(\mathbf{p}\right),\end{array}$$

(5)

where $a\ge 0$ is a given weighting factor.

Moreover, the infection rate of symptomatic individuals is the largest among the three kinds of infected individuals, that is,

$$\begin{array}{c}\hfill {\beta }_P\left(t\right)<{\beta }_I\left(t\right),{\beta }_A\left(t\right)<{\beta }_I\left(t\right).\end{array}$$

(6)

In addition, since the removed time for symptomatic individuals is longer than that for asymptomatic people, we have

$$\begin{array}{c}\hfill {\gamma }_I\left(t\right)<{\gamma }_U\left(t\right),{\gamma }_I\left(t\right)<{\gamma }_A\left(t\right).\end{array}$$

(7)

Therefore, based on (3)–(7), we propose an identification problem (IP for short):

$$\begin{array}{ccc}\hfill \left(\mathrm{IP}\right)& & min\mathcal{L}\left(\mathbf{p}\right),\hfill \\ & & \mathrm{s}.\mathrm{t}.\left\{\begin{array}{c}\mathbf{x}(\mathbf{p};t)\in \mathbb{X},\hfill \\ {\beta }_P\left(t\right)\le {\beta }_I\left(t\right),{\beta }_A\left(t\right)\le {\beta }_I\left(t\right),\hfill \\ {\gamma }_I\left(t\right)\le {\gamma }_U\left(t\right),{\gamma }_I\left(t\right)\le {\gamma }_A\left(t\right),\hfill \\ {\sigma }_{min}\le \sigma \left(t\right)\le {\sigma }_{max},\hfill \\ {\lambda }_{min}\le \lambda \left(t\right)\le {\lambda }_{max},\hfill \\ \mu \left(t\right)\le {\mu }_{max},\hfill \\ {\gamma }_U^{min}\le {\gamma }_U\left(t\right)\le {\gamma }_U^{max},\hfill \\ {\gamma }_A^{min}\le {\gamma }_A\left(t\right)\le {\gamma }_A^{max},\hfill \\ {\gamma }_I^{min}\le {\gamma }_I\left(t\right)\le {\gamma }_I^{max},\hfill \\ 0\le \mathbf{p}\left(t\right)\le 1,t\in [0,T].\hfill \end{array}\right.\hfill \end{array}$$

where $0<{\sigma }_{min}<{\sigma }_{max}<1,$ $0<{\lambda }_{min}<{\lambda }_{max}\le 1,$ $0<{\mu }_{max}<1,$ $0<{\gamma }_U^{min}<{\gamma }_U^{max}<1,$ $0<{\gamma }_A^{min}<{\gamma }_A^{max}<1,$ $0<{\gamma }_I^{min}<{\gamma }_I^{max}<1.$

To solve the above problem, we construct a mixed algorithm; its steps are as follows:

Step 1°.

Guess $\mathbf{x}\left(0\right),$ $a\ge 0,$ ${\epsilon }_1\ge 0,{\epsilon }_2>0,$ $m_1\in {\mathbb{Z}}^{+},m_2\in {\mathbb{Z}}^{+},m\in {\mathbb{Z}}^{+},N_d\in {\mathbb{Z}}^{+},m_1>m_2.$

Step 2°.

Divide the interval $[0,T]$ into m segments: $0=t_0\le t_1\le \cdots \le t_m=T.$ Introduce the following objective functional:

$$\begin{array}{ccc}\hfill \tilde{\mathcal{L}}\left(\mathbf{p}\right)& =& \sum _{i=0}^m\left[\frac{\left|T_A(\mathbf{p};t_i)-A_{data}\left(t_i\right)\right|}{A_{data}\left(t_i\right)}+\frac{\left|T_I(\mathbf{p};t_i)-I_{data}\left(t_i\right)\right|}{I_{data}\left(t_i\right)}+\frac{\left|R_R(\mathbf{p};t_i)-R_{data}\left(t_i\right)\right|}{R_{data}\left(t_i\right)}\right]\hfill \\ & & +a\sum _{i=0}^{m-1}\left[\left|{\beta }_P\left(\mathbf{p};t_{i+1}\right)-{\beta }_P\left(\mathbf{p};t_i\right)\right|+\left|{\beta }_A\left(\mathbf{p};t_{i+1}\right)-{\beta }_A\left(\mathbf{p};t_i\right)\right|+\left|{\beta }_I\left(\mathbf{p};t_{i+1}\right)-{\beta }_I\left(\mathbf{p};t_i\right)\right|\right].\hfill \end{array}$$

Then, the problem (IP) can be transformed into the following problem (NIP for short):

$$\begin{array}{ccc}\hfill \left(\mathrm{NIP}\right)& & min\tilde{\mathcal{L}}\left(\mathbf{p}\right),\hfill \\ & & \mathrm{s}.\mathrm{t}.\left\{\begin{array}{c}\mathbf{x}(\mathbf{p};t_i)\in \mathbb{X},\hfill \\ {\beta }_P\left(t_i\right)\le {\beta }_I\left(t_i\right),{\beta }_A\left(t_i\right)\le {\beta }_I\left(t_i\right),\hfill \\ {\gamma }_I\left(t_i\right)\le {\gamma }_U\left(t_i\right),{\gamma }_I\left(t_i\right)\le {\gamma }_A\left(t_i\right),\hfill \\ {\sigma }_{min}\le \sigma \left(t_i\right)\le {\sigma }_{max},\hfill \\ {\lambda }_{min}\le \lambda \left(t_i\right)\le {\lambda }_{max},\hfill \\ \mu \left(t_i\right)\le {\mu }_{max},\hfill \\ {\gamma }_U^{min}\le {\gamma }_U\left(t_i\right)\le {\gamma }_U^{max},\hfill \\ {\gamma }_A^{min}\le {\gamma }_A\left(t_i\right)\le {\gamma }_A^{max},\hfill \\ {\gamma }_I^{min}\le {\gamma }_I\left(t_i\right)\le {\gamma }_I^{max},\hfill \\ 0\le \mathbf{p}\left(t_i\right)\le 1,i\in {\mathbb{Z}}_m.\hfill \end{array}\right.\hfill \end{array}$$

Step 3°.

Randomly generate $m_1$ groups of initial values

$$\begin{array}{c}\hfill {\tilde{\mathbf{p}}}_i^0\left(0\right){=}^i{\left({\beta }_P\left(0\right),{\beta }_A\left(0\right),{\beta }_I\left(0\right),{\gamma }_U\left(0\right),{\gamma }_A\left(0\right),{\gamma }_I\left(0\right),\sigma \left(0\right),\mu \left(0\right),\lambda \left(0\right)\right)}^T\triangleq {\tilde{\mathbf{p}}}_i^0,i\in {\mathbb{Z}}_{m_1}.\end{array}$$

Let ${\tilde{\mathbf{p}}}_i^k\left(t\right)\left(i\in {\mathbb{Z}}_{m_1}\right)$ denote the kth iteration identification variable by the ith initial value ${\tilde{\mathbf{p}}}_i^0.$ Set $k=0$.

Step 4°.

For every ${\tilde{\mathbf{p}}}_i^0\left(i\in {\mathbb{Z}}_{m_1}\right),$ use the interior algorithm to solve the problem (NIP), and obtain ${\tilde{\mathbf{p}}}_i^{k+1}.$

Solving constrained optimization problems using interior-point methods is equivalent to solving a series of approximate minimization problems. In the context of the (NIP) problem, let $\mathbf{x}(\mathbf{p};t)\in \mathbb{X}$ denote the equality constraint $h\left(\mathbf{p}\right)$, and let other inequality constraints be denoted as $g\left(\mathbf{p}\right)$. Construct the following series of approximation problems:

$$\begin{array}{c}\hfill \underset{\mathbf{p},s}{min}{\tilde{\mathcal{L}}}_{{\mu }_j}(\mathbf{p},s)=\underset{\mathbf{p},s}{min}\tilde{\mathcal{L}}\left(\mathbf{p}\right)-\mu \sum _{j}ln\left(s_j\right),s.t.s\le 0,h\left(\mathbf{p}\right)=0,g\left(\mathbf{p}\right)+s=0.\end{array}$$

During the process of solving the approximation problems, if the problem exhibits local convexity around the current iteration, Newton’s method is employed for resolution. However, if it is non-convex, a trust region algorithm is applied instead.

Judge whether the stop condition is met: if the optimality tolerance is smaller than ${\epsilon }_1$ or $k>N_d$ holds, output ${\tilde{\mathbf{p}}}_i^{★}={\tilde{\mathbf{p}}}_i^{k+1}$, and turn to Step $5^{\mathrm{o}};$ otherwise, set $k=k+1.$

Step 5°.

From the $m_1$ groups of ${\tilde{\mathbf{p}}}_i^{★}(i\in {\mathbb{Z}}_{m_1}),$ choose $m_2$ groups which satisfy the constraint conditions and take them as the initial population. Let

$$\begin{array}{c}\hfill {\widehat{\mathbf{p}}}^l\left(t\right)={\left({\beta }_P\left(t\right),{\beta }_A\left(t\right),{\beta }_I\left(t\right),{\gamma }_U\left(t\right),{\gamma }_A\left(t\right),{\gamma }_I\left(t\right),\sigma \left(t\right),\mu \left(t\right),\lambda \left(t\right)\right)}_l^T,\end{array}$$

be the lth iteration identification variable. Set $l=0.$

Step 6°.

Use the genetic algorithm to solve the problem (NIP), and find ${\widehat{\mathbf{p}}}^{l+1}.$ Judge whether the stop condition is satisfied: if $\left|\tilde{\mathcal{L}}\left({\widehat{\mathbf{p}}}^{l+1}\right)-\tilde{\mathcal{L}}\left({\widehat{\mathbf{p}}}^l\right)\right|\le {\epsilon }_2,$ output the optimal solution ${\mathbf{p}}^{*}={\widehat{\mathbf{p}}}^{l+1};$ otherwise, set $l=l+1.$

3. Optimal Control Problem

To control the disease with effect, we propose two objectives during the epidemic period:

$$\begin{array}{ccc}\hfill {\mathbf{J}}_1\left(\mathbf{u}\right)& =& {\int }_0^T\frac{S\left(t\right)\left[{\beta }_P\left(t\right)P\left(t\right)+{\beta }_A\left(t\right)A\left(t\right)+{\beta }_I\left(t\right)I\left(t\right)\right]}{N\left(t\right)}dt,\hfill \\ \hfill {\mathbf{J}}_2\left(\mathbf{u}\right)& =& {\int }_0^T{\mathbf{u}}^2\left(t\right)dt.\hfill \end{array}$$

where ${\mathbf{J}}_1\left(\mathbf{u}\right)$ denotes the newly increased number of infected individuals and ${\mathbf{J}}_2\left(\mathbf{u}\right)$ denotes the costs of control measures. We introduce the following the combination functional with two coefficients $C_1(\ge 0)$ and $C_2(\ge 0):$

$$\begin{array}{c}\hfill \mathbf{J}\left(\mathbf{u}\right)=C_1\frac{{\mathbf{J}}_1\left(\mathbf{u}\right)}{{\mathbf{J}}_1\left(0\right)}+C_2\frac{{\mathbf{J}}_2\left(\mathbf{u}\right)}{{\mathbf{J}}_2\left(1\right)}.\end{array}$$

where ${\mathbf{J}}_1\left(0\right)$ is the newly increased number of infected individuals without control and ${\mathbf{J}}_2\left(1\right)$ is the control cost with $\mathbf{u}\equiv 1.$

The goal is to minimize the epidemic size and control costs simultaneously. Therefore, we put forward the following optimal control problem (OCP for short):

$$\begin{array}{ccc}\hfill \left(\mathrm{OCP}\right)& & min\mathbf{J}\left(\mathbf{u}\right),\hfill \\ & & \mathrm{s}.\mathrm{t}.\left\{\begin{array}{c}\mathbf{x}(\mathbf{u};t)\in \mathbb{X},\hfill \\ 0\le \mathbf{u}\left(t\right)\le 1.\hfill \end{array}\right.\hfill \end{array}$$

Based on ref. [41], we can obtain that there exists an optimal control ${\mathbf{u}}^{*}\left(t\right)$ which is the optimal solution of the problem (OCP). Next, we will present the characterization of the optimal control.

Theorem 1.

Let ${\mathbf{u}}^{*}\left(t\right)$ be the optimal control of the problem (OCP) and ${\mathbb{X}}^{*}={(S^{*}\left(t\right),P^{*}\left(t\right),A^{*}\left(t\right),I^{*}\left(t\right),T_A^{*}\left(t\right),T_I^{*}\left(t\right),R_U^{*}\left(t\right),R_R^{*}\left(t\right))}^T$ denote the corresponding solution of the control system. Then, there exist continuous functions $y_k\left(t\right)(k\in \{S,P,A,I,TA,TI,RU,RR\})$ satisfying the following system:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{y}_S^{\prime }\left(t\right)=\left\{\left[1-{\mathbf{u}}^{*}\left(t\right)\right]\left[y_S\left(t\right)-y_P\left(t\right)\right]-\frac{C_1}{{\mathbf{J}}_1\left(0\right)}\right\}·\frac{\left[N^{*}\left(t\right)-S^{*}\left(t\right)\right]\left[{\beta }_P\left(t\right)P^{*}\left(t\right)+{\beta }_A\left(t\right)A^{*}\left(t\right)+{\beta }_I\left(t\right)I^{*}\left(t\right)\right]}{{\left(N^{*}\right)}^2\left(t\right)},\hfill \\ y_P^{\prime }\left(t\right)=\left\{\left[1-{\mathbf{u}}^{*}\left(t\right)\right]\left[y_S\left(t\right)-y_P\left(t\right)\right]-\frac{C_1}{{\mathbf{J}}_1\left(0\right)}\right\}·\frac{S^{*}\left(t\right)\left\{{\beta }_P\left(t\right)[N^{*}\left(t\right)-P^{*}\left(t\right)]+{\beta }_A\left(t\right)A^{*}\left(t\right)+{\beta }_I\left(t\right)I^{*}\left(t\right)\right\}}{{\left(N^{*}\right)}^2\left(t\right)},\hfill \\ y_A^{\prime }\left(t\right)=\left\{\left[1-{\mathbf{u}}^{*}\left(t\right)\right]\left[y_S\left(t\right)-y_P\left(t\right)\right]-\frac{C_1}{{\mathbf{J}}_1\left(0\right)}\right\}·\frac{S^{*}\left(t\right)\left\{{\beta }_A\left(t\right)[N^{*}\left(t\right)-A^{*}\left(t\right)]+{\beta }_P\left(t\right)P^{*}\left(t\right)+{\beta }_I\left(t\right)I^{*}\left(t\right)\right\}}{{\left(N^{*}\right)}^2\left(t\right)}\hfill \\ +\mu \left(t\right)\left[y_A\left(t\right)-y_I\left(t\right)\right]+{\gamma }_U\left(t\right)\left[y_A\left(t\right)-y_{RU}\left(t\right)\right]+\lambda \left(t\right)\left[y_A\left(t\right)-y_{TA}\left(t\right)\right],\hfill \\ y_I^{\prime }\left(t\right)=\left\{\left[1-{\mathbf{u}}^{*}\left(t\right)\right]\left[y_S\left(t\right)-y_P\left(t\right)\right]-\frac{C_1}{{\mathbf{J}}_1\left(0\right)}\right\}·\frac{S^{*}\left(t\right)\left\{{\beta }_I\left(t\right)[N^{*}\left(t\right)-I^{*}\left(t\right)]+{\beta }_P\left(t\right)P^{*}\left(t\right)+{\beta }_A\left(t\right)A^{*}\left(t\right)\right\}}{{\left(N^{*}\right)}^2\left(t\right)}\hfill \\ +\lambda \left(t\right)\left[y_I\left(t\right)-y_{TI}\left(t\right)\right],\hfill \\ y_{TA}^{\prime }\left(t\right)=\left\{\left[1-{\mathbf{u}}^{*}\left(t\right)\right]\left[y_P\left(t\right)-y_S\left(t\right)\right]+\frac{C_1}{{\mathbf{J}}_1\left(0\right)}\right\}·\frac{S^{*}\left(t\right)\left[{\beta }_P\left(t\right)P^{*}\left(t\right)+{\beta }_A\left(t\right)A^{*}\left(t\right)+{\beta }_I\left(t\right)I^{*}\left(t\right)\right]}{{\left(N^{*}\right)}^2\left(t\right)}\hfill \\ +{\mu }_T\left(t\right)\left[y_{TA}\left(t\right)-y_{TI}\left(t\right)\right]+{\gamma }_A\left(t\right)\left[y_{TA}\left(t\right)-y_{RR}\left(t\right)\right],\hfill \\ y_{TI}^{\prime }\left(t\right)=\left\{\left[1-{\mathbf{u}}^{*}\left(t\right)\right]\left[y_P\left(t\right)-y_S\left(t\right)\right]+\frac{C_1}{{\mathbf{J}}_1\left(0\right)}\right\}·\frac{S^{*}\left(t\right)\left[{\beta }_P\left(t\right)P^{*}\left(t\right)+{\beta }_A\left(t\right)A^{*}\left(t\right)+{\beta }_I\left(t\right)I^{*}\left(t\right)\right]}{{\left(N^{*}\right)}^2\left(t\right)}\hfill \\ +{\gamma }_I\left(t\right)\left[y_{TI}\left(t\right)-y_{RR}\left(t\right)\right],\hfill \\ y_{RU}^{\prime }\left(t\right)=\left\{\left[1-{\mathbf{u}}^{*}\left(t\right)\right]\left[y_P\left(t\right)-y_S\left(t\right)\right]+\frac{C_1}{{\mathbf{J}}_1\left(0\right)}\right\}·\frac{S^{*}\left(t\right)\left[{\beta }_P\left(t\right)P^{*}\left(t\right)+{\beta }_A\left(t\right)A^{*}\left(t\right)+{\beta }_I\left(t\right)I^{*}\left(t\right)\right]}{{\left(N^{*}\right)}^2\left(t\right)},\hfill \\ y_{RR}^{\prime }\left(t\right)=\left\{\left[1-{\mathbf{u}}^{*}\left(t\right)\right]\left[y_P\left(t\right)-y_S\left(t\right)\right]+\frac{C_1}{{\mathbf{J}}_1\left(0\right)}\right\}·\frac{S^{*}\left(t\right)\left[{\beta }_P\left(t\right)P^{*}\left(t\right)+{\beta }_A\left(t\right)A^{*}\left(t\right)+{\beta }_I\left(t\right)I^{*}\left(t\right)\right]}{{\left(N^{*}\right)}^2\left(t\right)},\hfill \end{array}\right.\end{array}$$

with transversality conditions

$$\begin{array}{c}\hfill y_S\left(T\right)=y_P\left(T\right)=y_A\left(T\right)=y_I\left(T\right)=y_{TA}\left(T\right)=y_{TI}\left(T\right)=y_{RU}\left(T\right)=y_{RR}\left(T\right)=0.\end{array}$$

And the optimal control ${\mathbf{u}}^{*}\left(t\right)$ is presented by

$$\begin{array}{ccc}\hfill & & {\mathbf{u}}^{*}\left(t\right)=min\left\{max\left\{0,\frac{S^{*}\left(t\right){\mathbf{J}}_2\left(1\right)\left[y_P\left(t\right)-y_S\left(t\right)\right]\left[{\beta }_P\left(t\right)P^{*}\left(t\right)+{\beta }_A\left(t\right)A^{*}\left(t\right)+{\beta }_I\left(t\right)I^{*}\left(t\right)\right]}{2C_2{N}^{*}\left(t\right)}\right\},1\right\}.\hfill \end{array}$$

(8)

Proof. 

We construct a Hamilton function:

$$\begin{array}{ccc}\hfill \mathcal{H}& =& \frac{S\left(t\right)\left[{\beta }_P\left(t\right)P\left(t\right)+{\beta }_A\left(t\right)A\left(t\right)+{\beta }_I\left(t\right)I\left(t\right)\right]}{N\left(t\right){\mathbf{J}}_1\left(0\right)}+\frac{{\mathbf{u}}^2\left(t\right)}{{\mathbf{J}}_2\left(1\right)}+y_S\left(t\right)S^{\prime }\left(t\right)+y_P\left(t\right)E^{\prime }\left(t\right)\hfill \\ & & +y_A\left(t\right)A^{\prime }\left(t\right)+y_I\left(t\right)I^{\prime }\left(t\right)+y_{TA}\left(t\right)T_A^{\prime }\left(t\right)+y_{TI}\left(t\right)T_I^{\prime }\left(t\right)+y_{RU}\left(t\right)R_U^{\prime }\left(t\right)+y_{RR}\left(t\right)R_R^{\prime }\left(t\right)\hfill \\ & =& \frac{S\left(t\right)\left[{\beta }_P\left(t\right)P\left(t\right)+{\beta }_A\left(t\right)A\left(t\right)+{\beta }_I\left(t\right)I\left(t\right)\right]}{N\left(t\right){\mathbf{J}}_1\left(0\right)}+\frac{{\mathbf{u}}^2\left(t\right)}{{\mathbf{J}}_2\left(1\right)}\hfill \\ & & +y_S\left(t\right)\left\{-S\left(t\right)[1-\mathbf{u}\left(t\right)]\frac{{\beta }_P\left(t\right)P\left(t\right)+{\beta }_A\left(t\right)A\left(t\right)+{\beta }_I\left(t\right)I\left(t\right)}{N\left(t\right)}\right\}\hfill \\ & & +y_P\left(t\right)\left\{S\left(t\right)[1-\mathbf{u}\left(t\right)]\frac{{\beta }_P\left(t\right)P\left(t\right)+{\beta }_A\left(t\right)A\left(t\right)+{\beta }_I\left(t\right)I\left(t\right)}{N\left(t\right)}-\alpha P\left(t\right)\right\}\hfill \\ & & +y_A\left(t\right)\left\{\alpha \left[1-\sigma \left(t\right)\right]P\left(t\right)-\left[\mu \left(t\right)+{\gamma }_U\left(t\right)+\lambda \left(t\right)\right]A\left(t\right)\right\}\hfill \\ & & +y_I\left(t\right)\left[\alpha \sigma \left(t\right)P\left(t\right)+\mu \left(t\right)A\left(t\right)-\lambda \left(t\right)I\left(t\right)\right]\hfill \\ & & +y_{TA}\left(t\right)\left\{\lambda \left(t\right)A\left(t\right)-\left[{\mu }_T\left(t\right)+{\gamma }_A\left(t\right)\right]T_A\left(t\right)\right\}\hfill \\ & & +y_{TI}\left(t\right)\left[\lambda \left(t\right)I\left(t\right)-{\gamma }_I\left(t\right)T_I\left(t\right)+{\mu }_T\left(t\right)T_A\left(t\right)\right]\hfill \\ & & +y_{RU}\left(t\right)\left[{\gamma }_U\left(t\right)A\left(t\right)\right]+y_{RR}\left(t\right)\left[{\gamma }_A\left(t\right)T_A\left(t\right)+{\gamma }_I\left(t\right)T_I\left(t\right)\right].\hfill \end{array}$$

Applying the Pontryagin’s minimum principle, we have

$$\begin{array}{ccc}& & y_S^{\prime }\left(t\right)=-\frac{\partial \mathcal{H}}{\partial S}=\left\{\left[1-\mathbf{u}\left(t\right)\right]\left[y_S\left(t\right)-y_P\left(t\right)\right]-\frac{C_1}{{\mathbf{J}}_1\left(0\right)}\right\}·\frac{\left[N\left(t\right)-S\left(t\right)\right]\left[{\beta }_P\left(t\right)P\left(t\right)+{\beta }_A\left(t\right)A\left(t\right)+{\beta }_I\left(t\right)I\left(t\right)\right]}{N^2\left(t\right)},\hfill \\ & & y_P^{\prime }\left(t\right)=-\frac{\partial \mathcal{H}}{\partial P}=\left\{\left[1-\mathbf{u}\left(t\right)\right]\left[y_S\left(t\right)-y_P\left(t\right)\right]-\frac{C_1}{{\mathbf{J}}_1\left(0\right)}\right\}·\frac{S\left(t\right)\left\{{\beta }_P\left(t\right)[N\left(t\right)-P\left(t\right)]+{\beta }_A\left(t\right)A\left(t\right)+{\beta }_I\left(t\right)I\left(t\right)\right\}}{N^2\left(t\right)},\hfill \\ & & y_A^{\prime }\left(t\right)=-\frac{\partial \mathcal{H}}{\partial A}=\left\{\left[1-\mathbf{u}\left(t\right)\right]\left[y_S\left(t\right)-y_P\left(t\right)\right]-\frac{C_1}{{\mathbf{J}}_1\left(0\right)}\right\}·\frac{S\left(t\right)\left\{{\beta }_A\left(t\right)[N\left(t\right)-A\left(t\right)]+{\beta }_P\left(t\right)P\left(t\right)+{\beta }_I\left(t\right)I\left(t\right)\right\}}{N^2\left(t\right)}\hfill \\ & & +\mu \left(t\right)\left[y_A\left(t\right)-y_I\left(t\right)\right]+{\gamma }_U\left(t\right)\left[y_A\left(t\right)-y_{RU}\left(t\right)\right]+\lambda \left(t\right)\left[y_A\left(t\right)-y_{TA}\left(t\right)\right],\hfill \\ & & y_I^{\prime }\left(t\right)=-\frac{\partial \mathcal{H}}{\partial I}=\left\{\left[1-\mathbf{u}\left(t\right)\right]\left[y_S\left(t\right)-y_P\left(t\right)\right]-\frac{C_1}{{\mathbf{J}}_1\left(0\right)}\right\}·\frac{S\left(t\right)\left\{{\beta }_I\left(t\right)[N\left(t\right)-I\left(t\right)]+{\beta }_P\left(t\right)P\left(t\right)+{\beta }_A\left(t\right)A\left(t\right)\right\}}{N^2\left(t\right)}\hfill \\ & & +\lambda \left(t\right)\left[y_I\left(t\right)-y_{TI}\left(t\right)\right],\hfill \\ & & y_{TA}^{\prime }\left(t\right)=-\frac{\partial \mathcal{H}}{\partial T_A}=\left\{\left[1-\mathbf{u}\left(t\right)\right]\left[y_P\left(t\right)-y_S\left(t\right)\right]+\frac{C_1}{{\mathbf{J}}_1\left(0\right)}\right\}·\frac{S\left(t\right)\left[{\beta }_P\left(t\right)P\left(t\right)+{\beta }_A\left(t\right)A\left(t\right)+{\beta }_I\left(t\right)I\left(t\right)\right]}{N^2\left(t\right)}\hfill \\ & & +{\mu }_T\left(t\right)\left[y_{TA}\left(t\right)-y_{TI}\left(t\right)\right]+{\gamma }_A\left(t\right)\left[y_{TA}\left(t\right)-y_{RR}\left(t\right)\right],\hfill \\ & & y_{TI}^{\prime }\left(t\right)=-\frac{\partial \mathcal{H}}{\partial T_I}=\left\{\left[1-\mathbf{u}\left(t\right)\right]\left[y_P\left(t\right)-y_S\left(t\right)\right]+\frac{C_1}{{\mathbf{J}}_1\left(0\right)}\right\}·\frac{S\left(t\right)\left[{\beta }_P\left(t\right)P\left(t\right)+{\beta }_A\left(t\right)A\left(t\right)+{\beta }_I\left(t\right)I\left(t\right)\right]}{N^2\left(t\right)}\hfill \\ & & +{\gamma }_I\left(t\right)\left[y_{TI}\left(t\right)-y_{RR}\left(t\right)\right],\hfill \\ & & y_{RU}^{\prime }\left(t\right)=-\frac{\partial \mathcal{H}}{\partial R_U}=\left\{\left[1-\mathbf{u}\left(t\right)\right]\left[y_P\left(t\right)-y_S\left(t\right)\right]+\frac{C_1}{{\mathbf{J}}_1\left(0\right)}\right\}·\frac{S\left(t\right)\left[{\beta }_P\left(t\right)P\left(t\right)+{\beta }_A\left(t\right)A\left(t\right)+{\beta }_I\left(t\right)I\left(t\right)\right]}{N^2\left(t\right)},\hfill \\ & & y_{RR}^{\prime }\left(t\right)=-\frac{\partial \mathcal{H}}{\partial R_R}=\left\{\left[1-\mathbf{u}\left(t\right)\right]\left[y_P\left(t\right)-y_S\left(t\right)\right]+\frac{C_1}{{\mathbf{J}}_1\left(0\right)}\right\}·\frac{S\left(t\right)\left[{\beta }_P\left(t\right)P\left(t\right)+{\beta }_A\left(t\right)A\left(t\right)+{\beta }_I\left(t\right)I\left(t\right)\right]}{N^2\left(t\right)}.\hfill \end{array}$$

And,

$$\begin{array}{c}\hfill \frac{\partial \mathcal{H}}{\partial \mathbf{u}}=\frac{2\mathbf{u}{C}_2}{{\mathbf{J}}_2\left(1\right)}+\frac{S\left(t\right)\left[y_S\left(t\right)-y_P\left(t\right)\right]\left[{\beta }_P\left(t\right)P\left(t\right)+{\beta }_A\left(t\right)A\left(t\right)+{\beta }_I\left(t\right)I\left(t\right)\right]}{N\left(t\right)}=0.\end{array}$$

Solving for $\mathbf{u}\left(t\right)$ yields

$$\begin{array}{c}\hfill \mathbf{u}\left(t\right)=\frac{S\left(t\right){\mathbf{J}}_2\left(1\right)\left[y_P\left(t\right)-y_S\left(t\right)\right]\left[{\beta }_P\left(t\right)P\left(t\right)+{\beta }_A\left(t\right)A\left(t\right)+{\beta }_I\left(t\right)I\left(t\right)\right]}{2C_{2}N\left(t\right)}.\end{array}$$

Using the standard argument for bounds $0\le \mathbf{u}\left(t\right)\le 1$, one obtains (8). □

4. Numerical Results

4.1. Data Description and Parameter Value

We collected the epidemic data of Shanghai from 1 March to 30 May 2022 from the official website of Shanghai Municipal Health Commission [42]. These data include the daily confirmed cases in the hospital, the daily asymptomatic carriers under medical observation, the daily cases from the asymptomatic to the confirmed, the cumulative asymptomatic cases discharged from medical observation, the cumulative cured cases discharged from the hospital, and the cumulative death cases. We take the time unit as one day.

According to the 2022 Shanghai Statistics Yearbook [43], the total population of Shanghai at the end of 2021 was 24,894,300. So, we set the initial value of the susceptible individuals as 24,894,300. At the beginning of the epidemic, there was one confirmed case in Putuo District and ten asymptomatic cases; one asymptomatic case was discharged from medical observation from 0:00 to 24:00 on 1 March 2022. So, $T_I\left(0\right)$ is set as 1, $T_A\left(0\right)$ is set as 10, and $R_R\left(0\right)$ is set as 1. Other initial values are indicated in

Table 1. Initial values of the COVID-19 system.

Parameter	Value	Source
$S\left(0\right)$	24,894,300	[43]
$P\left(0\right)$	0	Estimated
$A\left(0\right)$	14	Estimated
$I\left(0\right)$	2	Estimated
$T_A\left(0\right)$	10	[42]
$T_I\left(0\right)$	1	[42]
$R_R\left(0\right)$	1	[42]
$R_U\left(0\right)$	0	Estimated

.

Ref. [44] points out that the incubation period of the Omicron variant mostly ranges from 2 to 4 days. Thus, $\alpha$ is set to $\frac{1}{3}.$ The ratio $\sigma$ is calculated based on the data of the daily confirmed cases and the daily asymptomatic cases. Its minimum and maximum are 0.0192 and 0.4458, respectively. Ref. [45] introduces the fact that several rounds of nucleic acid screening were performed during the epidemic period in Shanghai; thus, we estimate the range of the tested rates to be from $\frac{2}{3}$ to 1. According to the data of the daily confirmed cases and the daily asymptomatic cases, ${\mu }_{max}$ is taken as 0.02693. For the asymptomatic individuals, we set ${\gamma }_U^{max}={\gamma }_A^{max}=\frac{1}{3}$ and ${\gamma }_U^{min}={\gamma }_A^{min}=\frac{1}{40}.$ For the symptomatic individuals, the length of stay in the hospital is 6 days [12], so we estimate ${\gamma }_A^{max}=\frac{1}{6};$ moreover, symptom disappearance takes 2 to 3 months [46], so we set ${\gamma }_I^{min}=\frac{1}{60}.$ In reference to the parameters associated with genetic algorithms, this paper assigns the population size as $p_s$, the crossover rate as $c_r$, and the mutation rate as $m_r$, and denotes the maximum number of iterations for algorithm termination as $N_m$. The parameter values in the control system, problem (IP), and genetic algorithm are shown in

Table 2. Parameter values throughout this article.

Parameter	Value	Source	Parameter	Value	Source
$\alpha$	$\frac{1}{3}$	[44]	${\gamma }_U^{max}$	$\frac{1}{3}$	Estimated
${\sigma }_{min}$	0.0192	[42]	${\gamma }_U^{min}$	$\frac{1}{40}$	Estimated
${\sigma }_{max}$	0.4458	[42]	a	1	Estimated
${\lambda }_{min}$	$\frac{2}{3}$	Estimated	${\epsilon }_1$	${10}^{-6}$	Estimated
${\lambda }_{max}$	1	Estimated	${\epsilon }_2$	${10}^{-6}$	Estimated
${\mu }_{max}$	0.02693	[42]	$p_s$	100	Estimated
${\gamma }_I^{min}$	$\frac{1}{60}$	[46]	$c_r$	0.8	Estimated
${\gamma }_I^{max}$	$\frac{1}{6}$	[12]	$m_r$	0.2	Estimated
${\gamma }_A^{min}$	$\frac{1}{40}$	Estimated	$N_m$	$5\times {10}^3$	Estimated
${\gamma }_A^{max}$	$\frac{1}{3}$	Estimated	$N_d$	$10\times {10}^3$	Estimated

.

The numerical simulations in this study were conducted on a personal computer equipped with an AMD Ryzen 5 5600 G processor, operating at a frequency of 3.9 GHz, alongside 32 GB of DDR4 memory, and a 500 GB SSD hard drive. The computer ran on the Windows 10 operating system (64-bit). MATLAB R2021a served as the primary numerical computing tool. All experiments were performed locally on the computer, without the use of remote computing resources.

4.2. Parameter Identification

Using the mixed algorithm and the parameter values in Table 1 and Table 2, we obtain the results of the identified parameters indicated in Figure 2. In Figure 2, we can observe that the infection rate of symptomatic individuals is the highest among the three infection rates.

To validate the rationality of the identified results, we simulate the real-time number of infected individuals with these values and compare them to the actual data shown in Figure 3. From Figure 3, we can see that the simulation results with the identification parameters are consistent with the actual process. In addition, we calculate the coefficients of determination as

$$\begin{array}{c}\hfill R_{TA}^2=0.9693,R_{TI}^2=0.9906,R_{RR}^2=0.9989.\end{array}$$

Thus, we can conclude that the obtained identification parameters are practical and applicable.

4.3. Optimal Control

Employing the Runge–Kutta scheme and the forward–backward sweep method, we present the optimal controls with different combinations of coefficients based on Theorem 1.

Our purpose is to show the control costs corresponding to the different infection peaks by adjusting $C_1$ and $C_2$. Five combinations of $C_1$ and $C_2$ are shown in

Table 3. Combinations of $C_1,C_2$, and optimal values.

	Case 1	Case 2	Case 3	Case 4	Case 5
$C_1$	1	28	451	1158	1258
$C_2$	1	1	1	1	1
$\rho$	92.7395%	50.1279%	10%	4.9980%	4.6941%
${\mathbf{J}}_1\left({\mathbf{u}}^{*}\right)$	655,816	359,886	73,119	36,700	34,466
${\mathbf{J}}_2\left({\mathbf{u}}^{*}\right)$	0.0002	0.0440	0.5146	0.8742	0.9116

. The optimal controls ${\mathbf{u}}^{*}\left(t\right)$ corresponding to the five combinations of $C_1$ and $C_2$ are indicated in Figure 4. We introduce the factor $\rho ,$ which is the ratio of the maximum number of tested symptomatic individuals with optimal controls to the official maximum. The maximum number of reported confirmed cases is 25,010 [42], so $\rho =\frac{max\left\{T_I\left({\mathbf{u}}^{*}\right)\right\}}{25,010}.$ We can assess the effects of optimal control measures on epidemic peaks using the factor $\rho .$ Case 2, Case 3, and Case 4 can control the epidemic peak to about $50\%,$$10\%$, and $5\%$ of the reported maximum, which correspond to 12,537; 2501; and 1250 people. The maximums of the tested symptomatic individuals for Case 1 and Case 5 are 23,194 and 1174, which correspond to two extremes. For $C_1<1,$ the control measures do not work; for $C_1>1258,$ the infection peak cannot decline. And the optimal newly increased number of infected individuals and the optimal costs for five cases are also shown in Table 3. We can see that the epidemic peak becomes smaller with the optimal epidemic size decreasing and the optimal control costs enlarging. From Figure 4, we can observe that the optimal controls decrease with time in general. With the increases in $C_1,$ the control costs increase. From Table 3 and Figure 4, we can conclude that greater control costs indicates better effects in the interval $[1,1258]$ of $C_1.$

Using the five coefficient combinations in Table 3 and the optimal controls in Figure 4, we show comparisons of the infected and removed individuals in Figure 5 and Figure 6. In Figure 5, the number of infected individuals with optimal controls is lower than the reported data during the epidemic period. We present the maximums of the number of infected individuals with optimal controls for five cases. For Case 1, the maximum numbers of pre-symptomatic, asymptomatic, tested asymptomatic, and symptomatic individuals are about 123,798; 36,025; 245,117; and 3067. For Case 2, the maximum numbers of pre-symptomatic, asymptomatic, tested asymptomatic, and symptomatic individuals are about 65,636; 19,113; 132,385; and 1626. For Case 3, the maximum numbers of pre-symptomatic, asymptomatic, tested asymptomatic, and symptomatic individuals are 12,371; 3609; 26,353; and 304. For Case 4, the maximum numbers of pre-symptomatic, asymptomatic, tested asymptomatic, and symptomatic individuals are about 6007; 1754; 13,127; and 147. For Case 5, the maximum numbers of pre-symptomatic, asymptomatic, tested asymptomatic, and symptomatic individuals are 5624; 1642; 12,319; and 138. We find that all maximums with optimal controls are less than those without control or the official data. In Figure 6, we show the removed individuals with different combinations. The maximum numbers of undetected and reported individuals are about 77,780 and 597,816 for Case 1; 41,943 and 323,661 for Case 2; 8062 and 63,496 for Case 3; 3941 and 31,379 for Case 4; and 3692 and 29,428 for Case 5. From Figure 5 and Figure 6, we can see that the optimal control measures with five cases are effective and can provide reference for decision-makers.

Figure 7 shows the time-variant control reproduction number with different combinations of $C_1$ and $C_2.$ The black curve denotes the basic reproduction number $R_0$. The time-variant control reproduction numbers are higher than one. The control reproduction numbers with different combinations of $C_1$ and $C_2$ are lower than the basic reproduction number during the epidemic period.

5. Conclusions

The purpose of this paper is to minimize the epidemic size and the control costs simultaneously. A control epidemic system of COVID-19 was established. We derived the control reproduction number. To estimate the time-varying parameter in the control system, we applied the parameter identification method and constructed a mixed algorithm to solve it. An optimal control problem with two objective functions was presented. We used the official data of Shanghai in 2022 to perform the numerical simulations. The optimal control measures with five coefficient combinations were presented.

The parameter identification method utilized in this study can be extrapolated to other infectious disease models, facilitating a deeper examination of disparities between parameters estimated via statistical methods. Drawing from official data from Shanghai in 2022, this paper deliberates optimal control strategies, elucidating the varying impacts of different control measures on epidemic peaks. These research findings furnish invaluable benchmarks for potential future outbreaks of similar diseases, equipping decision-makers with strategic insights. Table 3 showcases optimal control strategies across five distinct combinations, paving the way for the delineation of diverse optimal control objectives in the future and the exploration of their ramifications on epidemic peaks.