Fractional Modelling and Optimal Control of COVID-19 Transmission in Portugal

Abstract

A fractional-order compartmental model was recently used to describe real data of the first wave of the COVID-19 pandemic in Portugal [Chaos Solitons Fractals 144 (2021), Art. 110652]. Here, we modify that model in order to correct time dimensions and use it to investigate the third wave of COVID-19 that occurred in Portugal from December 2020 to February 2021, and that has surpassed all previous waves, both in number and consequences. A new fractional optimal control problem is then formulated and solved, with vaccination and preventive measures as controls. A cost-effectiveness analysis is carried out, and the obtained results are discussed.

1. Introduction

In January 2020, the World Health Organization (WHO) announced the existence of a significant number of pneumonia cases in Wuhan. Against all the predictions, COVID-19 (COrona VIrus Disease-19) spread quickly across the globe and, on 11th of March, was declared as a pandemic [1]. Caused by SARS-CoV-2 (Severe Acute Respiratory Syndrome Corona Virus 2), COVID-19 is the first pandemic in the digital era from which very few territories of the world are untouched. Many governments were forced to decree measures that seemed to be outdated, such as the isolation of individuals and the complete lockdown of regions and even countries, that compromise individual freedoms, damage business and economy, and threaten a significant number of jobs.

To fight COVID-19 and its harmful effects, a multidisciplinary approach is needed. In particular, mathematical modelling plays an important role in the prediction of possible scenarios and in its effective control [2,3]. Readers interested in fractional modelling are referred to [4,5] and references therein.

The pandemic numbers have put national health systems under pressure. Many reported cases were not reported on time, but with a delay of days. Hence, in this paper, we do not use the number of daily reported cases but the means of the previous five days of reported cases, as suggested in [6].

The mean of five days of daily reported cases induces memory into the model. Fractional derivatives have been intensively used to obtain models of infectious diseases that take into account the memory effects. Many researchers have focused particular attention in modelling real-world phenomena using non-integer order derivatives. Those dynamics have been modelled and studied by using the concept of fractional-order derivatives. These problems appear in a range of diversified fields of applied sciences, including biology, physics, ecology, and engineering [7,8].

A classical compartmental model with a super-spreader class was firstly applied to give an estimation of infected subjects and fatalities in Wuhan, China, in [9]. The first-order derivative was then substituted by a derivative of a fractional order, resulting in a model investigated with Caputo fractional derivatives [6]. Here, this fractional model is corrected and then used to model the third wave of COVID-19 in Portugal. We start by the estimation of parameters that best fit the real data. The sensitivity analysis to the fractional-order model is performed in order to identify which model parameters are most influential on the dynamics of the disease. Afterwards, fractional optimal control is applied, showing the effectiveness of our approach.

This paper is organized as follows. In Section 2, the fractional order model is formulated. Our main results are then given in Section 3: parameter estimation of the COVID-19 model with real data of Portugal (Section 3.1); sensitivity analysis of the parameters of the model, without forgetting the effect of the order of fractional differentiation (Section 3.2); fractional optimal control of the model (Section 3.3); and, finally, numerical simulations and cost-effectiveness of the fractional model (Section 3.4). We end with Section 4, which states the conclusions. Our main theoretical contributions consist of fractional-order model consistency and the mathematical problem rearrangement according to the Pontryagin theory. Moreover, we solve the fractional optimal control problem numerically by a method of Adams–Basforth–Moulton.

2. Fractional-Order COVID-19 Model

Since the time interval considered in this study is small, we assume that population is invariant. In addition, population is divided in eight classes: (i) susceptible individuals, S; (ii) exposed individuals, E; (iii) symptomatic infectious individuals, I; (iv) super-spreader individuals, P; (v) asymptomatic infectious individuals, A; (vi) hospitalized individuals, H; (vii) recovered individuals, R; and (viii) fatality class, F. Our fractional-order model is derived from the one presented in Ndaïrou et al. [6], which gives a generalization of an integer-order model that was used to study the start of the pandemic in Wuhan [9]. We use the fractional derivative in the sense of Caputo. Fractional differential equations are an active area of research and are adequate to incorporate the history of the processes. The model system of equations for COVID-19 proposed in [6] is given by

$$\left\{\begin{array}{cc}{}_0^C{D}_t^{\alpha }S=\hfill & -\beta {\displaystyle \frac{IS}{N}}-l\beta {\displaystyle \frac{HS}{N}}-{\beta }^{\prime }{\displaystyle \frac{PS}{N}},\hfill \\ {}_0^C{D}_t^{\alpha }E=\hfill & \beta {\displaystyle \frac{IS}{N}}+l\beta {\displaystyle \frac{HS}{N}}+{\beta }^{\prime }{\displaystyle \frac{PS}{N}}-\kappa E,\hfill \\ {}_0^C{D}_t^{\alpha }I=\hfill & \kappa {\rho }_{1}E-({\gamma }_a+{\gamma }_i)I-{\delta }_{i}I,\hfill \\ {}_0^C{D}_t^{\alpha }P=\hfill & \kappa {\rho }_{2}E-({\gamma }_a+{\gamma }_i)P-{\delta }_{p}P,\hfill \\ {}_0^C{D}_t^{\alpha }A=\hfill & \kappa (1-{\rho }_1-{\rho }_2)E,\hfill \\ {}_0^C{D}_t^{\alpha }H=\hfill & {\gamma }_a(I+P)-{\gamma }_{r}H-{\delta }_{h}H,\hfill \\ {}_0^C{D}_t^{\alpha }R=\hfill & {\gamma }_i(I+P)+{\gamma }_{r}H,\hfill \\ {}_0^C{D}_t^{\alpha }F=\hfill & {\delta }_{i}I+{\delta }_{p}P+{\delta }_{h}H,\hfill \end{array}\right.$$

(1)

where ${}_0^C{D}_t^{\alpha }$ denotes the left Caputo fractional-order derivative with derivative order $\alpha$ ($0<\alpha ⩽1$). A description of the parameters of Model (1) can be found in

Table 1. Description and values of the parameters of Model (1) taken from [9,10].

Name	Description	Value
$\beta$	human-to-human transmission coefficient	2.55
l	transmissibility of hospitalized patients	1.56
${\beta }^{\prime }$	transmission coefficient of super-spreaders	7.65
$\kappa$	rate at which an individual leaves the exposed	0.25
	class to become infectious
${\rho }_1$	proportion of progression from class E	0.58
	to symptomatic infectious class I
${\rho }_2$	rate at which exposed ind. become super-spreaders	0.001
${\gamma }_a$	rate at which symptomatic and super-spreaders	0.94
	become hospitalized
${\gamma }_i$	recovery rate without being hospitalized	0.27
${\gamma }_r$	recovery rate of hospitalized patients	0.5
${\delta }_i$	disease induced death rate due to infected ind.	1/23
${\delta }_p$	disease induced death rate due to super-spreader ind.	1/23
${\delta }_h$	disease induced death rate due to hospitalized ind.	1/23

.

We note that the equations of Model (1) do not have appropriate time dimensions. Indeed, on the left-hand side, dimension is (time)${}^{-\alpha }$, while on the right-hand side, dimension is (time)${}^{-1}$. This means that Model (1) is only consistent when $\alpha =1$. For the importance to be consistent with dimensions, we refer the reader, e.g., to [11,12]. Thus, here, we correct System (1) as follows:

$$\left\{\begin{array}{cc}{}_0^C{D}_t^{\alpha }S=\hfill & -{\beta }^{\alpha }{\displaystyle \frac{IS}{N}}-l{\beta }^{\alpha }{\displaystyle \frac{HS}{N}}-{{\beta }^{\prime }}^{\alpha }{\displaystyle \frac{PS}{N}},\hfill \\ {}_0^C{D}_t^{\alpha }E=\hfill & {\beta }^{\alpha }{\displaystyle \frac{IS}{N}}+l{\beta }^{\alpha }{\displaystyle \frac{HS}{N}}+{{\beta }^{\prime }}^{\alpha }{\displaystyle \frac{PS}{N}}-{\kappa }^{\alpha }E,\hfill \\ {}_0^C{D}_t^{\alpha }I=\hfill & {\kappa }^{\alpha }{\rho }_{1}E-({\gamma }_a^{\alpha }+{\gamma }_i^{\alpha })I-{\delta }_i^{\alpha }I,\hfill \\ {}_0^C{D}_t^{\alpha }P=\hfill & {\kappa }^{\alpha }{\rho }_{2}E-({\gamma }_a^{\alpha }+{\gamma }_i^{\alpha })P-{\delta }_p^{\alpha }P,\hfill \\ {}_0^C{D}_t^{\alpha }A=\hfill & {\kappa }^{\alpha }(1-{\rho }_1-{\rho }_2)E,\hfill \\ {}_0^C{D}_t^{\alpha }H=\hfill & {\gamma }_a^{\alpha }(I+P)-{\gamma }_r^{\alpha }H-{\delta }_h^{\alpha }H,\hfill \\ {}_0^C{D}_t^{\alpha }R=\hfill & {\gamma }_i^{\alpha }(I+P)+{\gamma }_r^{\alpha }H,\hfill \\ {}_0^C{D}_t^{\alpha }F=\hfill & {\delta }_i^{\alpha }I+{\delta }_p^{\alpha }P+{\delta }_h^{\alpha }H.\hfill \end{array}\right.$$

(2)

A flow diagram of Model (2) is given in Figure 1. Note that, in the particular case $\alpha =1$, both Models (1) and (2) coincide with the classical COVID-19 model first introduced in [9]. System (2) will be investigated in Section 3.

3. Main Results

We begin by discussing the adherence of the corrected Model (2) introduced in Section 2 with respect to COVID-19 and real data from the third wave that occurred in Portugal. For that, we need an adequate estimation of the parameter values.

3.1. Parameter Estimation

The uncorrected model (1) was used to study the dynamics of COVID-19 at the early stage of the pandemic, between 3 March 2020 and 27 April 2020, in [6]. During that period, government decreed a general lockdown of the country that was well accepted by the population due to the impact of the disease in other countries where it started earlier. With limited knowledge, this revealed to be an effective measure to reduce contacts and control the disease in a relatively short period of time.

That initial stage ended in May 2020, with the gradual release of the country from COVID-19 container measures and the cancelling of the State of Emergency. In that date, preventive measures were adopted to control the disease. One of the most important was the use of masks in confined spaces, made mandatory by Portuguese Decreto-Lei n.º 20/2020.

In October 2020, due to autumn weather conditions, the number of infected individuals started to be worrisome. In order to control it and to protect the Portuguese National Health Service, in November 2020, a new series of States of Emergency began. In January 2021, the cold weather, some relaxation among the population concerning preventive measures (being in crowds or poorly ventilated spaces, the misuse of masks, …) combined with new variants of COVID-19—more contagious than ever— that started circulating, forced the government to declare a new lockdown with the closure of schools in 22 January 2021.

Due to the implementation of the preventive measures, during our model fitting, we assume that parameters have the same values of the first wave, with the exception of contact rates. The period we chose starts in 27 December 2020 and ends 16 February 2021, covering the third wave of the pandemic in Portugal. During that period, the closure of schools was declared, and that most certainly had an impact on contact rates. Hence, we consider that the transmission coefficient $\beta$ is replaced by $\beta (1-m(t\left)\right)$, and that ${\beta }^{\prime }$ is replaced by ${\beta }^{\prime }(1-m\left(t\right))$, where $m\left(t\right)$ is a continuous function that represents the rate of reduction of the contacts, and that varies with time according to Figure 2.

The values of two parameters were determined by the fitting of the model: (i) derivative order of the model, $\alpha$, and (ii) a scaling factor s. See

Table 2. Results of model fitting.

Derivative Order	s	Absolute Error	Relative Error (%)
1.0	21.08	8595	14.13
0.99	19.87	8135	13.37

for the resulting values and errors associated with the fitting.

The fitting curves are presented in Figure 3, where real data was obtained from [13,14].

It is known that the available data has some mistakes. Frequently, days reporting few cases are followed by days reporting many new cases, without correspondence with reality. Due to the stress that the pandemic imposed over doctors and health professionals, these data do not flow as quickly and effectively as expected. Therefore, the numbers of new cases are not correctly determined by reported daily numbers. Thus, the number of new cases considered in this manuscript is, as suggested in [6], the mean of the previous five days of reported cases.

Following [15], our data fitting consisted in minimization of the $l_2$ norm of the difference between the real values and predictive cases of COVID-19 infection given by Model (2). Due to the oscillation of the number of confirmed cases per day, the fitting curves observed in Figure 3 (left) have significant gaps in comparison with real data. Figure 3 (right) presents the difference between the number of confirmed cases per day and the number of estimated cases, showing that the proposed model approximates well the average number of cases.

Consequently, fitting errors are quite high, being 14.13% for the classical integer-order model ($\alpha =1$) and 13.37% for the fractional-order model with $\alpha =0.99$, according to Table 2. The difference between the two derivative orders, in terms of absolute errors, justifies the preference for the fractional-order model with respect to the classical one.

3.2. Sensitivity Analysis

An important threshold, while studying infectious disease models, is the basic reproduction number. Following [9,10], we conclude that the basic reproduction number of the COVID-19 model (2) is

$$R_0=\frac{{\beta }^{\alpha }{\rho }_1({\gamma }_a^{\alpha }l+a_h)}{a_i{a}_h}+\frac{({\beta }^{\alpha }{\gamma }_a^{\alpha }l+{{\beta }^{\prime }}^{\alpha }{a}_h){\rho }_2}{a_p{a}_h},$$

(3)

where $a_i={\gamma }_a^{\alpha }+{\gamma }_i^{\alpha }+{\delta }_i^{\alpha }$, $a_p={\gamma }_a^{\alpha }+{\gamma }_i^{\alpha }+{\delta }_p^{\alpha }$ and $a_h={\gamma }_r^{\alpha }+{\delta }_h^{\alpha }$.

The impact of the variation of the derivative order $\alpha$, in the evolution of $R_0$, is presented in Figure 4. We observe in Figure 4 that $R_0⩾1$, that is, we have an endemic scenario regardless of the value of $\alpha$.

Sensitivity analysis measures the importance of each parameter of the model in the disease transmission.

Definition 1

(The sensitivity index [16,17]). Let $R_0$ be differentiable with respect to a given parameter p. The sensitivity of the model with respect to that parameter p can be measured by the index ${{\rm Y}}_p^{R_0}$ given by

$${{\rm Y}}_p^{R_0}=\frac{\partial R_0}{\partial p}\frac{p}{R_0}.$$

The sensitivity analysis of the classical COVID-19 model was presented in [9], that is, in the particular situation $\alpha =1$ in (2). In it, the authors concluded that the most sensitive parameters to the basic reproduction number $R_0$ are $\beta$, ${\rho }_1$, and l. Therefore, special attention should be paid to the estimation of those parameters. In contrast, the estimation of ${\beta }^{\prime }$, ${\rho }_2$, ${\gamma }_a$, ${\delta }_i$, ${\delta }_p$, and ${\delta }_h$ does not require as much attention because of its low sensitivity.

In the general situation of our fractional-order model (2), it should be emphasized that the sensitivity depends on the derivative order $\alpha$ of the fractional operator. We can observe this in Figure 5, where one clearly sees that the variation of $\alpha$ influences the sensitivity index of parameters $\beta$, ${\rho }_1$, ${\gamma }_i$, ${\gamma }_r$, and l.

We see in Figure 5 that the sensitivity indexes of $\beta$ and ${\rho }_1$ exhibit a quite similar evolution with respect to $\alpha$, being very sensitive to the variation of $\alpha$. On the opposite side, the sensitivity indexes of ${\gamma }_i$ and l are much less sensitive to the variation of the fractional-order $\alpha$. The graphics with respect to the remaining parameters of the model are omitted here because, using the same scale, their curves do not go far from the x-axis.

The evolution of the sensitivity index for the basic reproduction number $R{}_0$ with the variation of the derivative order $\alpha$ is presented in Figure 6. We observe that the sensitivity index decreases with the decrease of $\alpha$.

Figure 5 and Figure 6 also show that the sensitivity of the model decreases in absolute value with the decrease in the derivative order. This means that the fractional-order model is less sensitive than the classical one, which is an advantage of our model.

3.3. Fractional Optimal Control of the Model

In this section, we aim to minimize the number of COVID-19-infected individuals and reduce the cost of control measures. This is carried out through (i) the use of vaccination as an effective measure time-dependent control $v\left(t\right)$ and (ii) the use of the preventive measure control $m\left(t\right)$ (representing the use of masks, limitations of the number of individuals in closed confined spaces, etc.), in order to control person-to-person contact. Therefore, the following fractional-order optimal control problem is considered:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill min\mathcal{J}(I,P,v,& m)={\int }_0^{t_f}\left(k_{1}I+k_{2}P+k_3{v}^2+k_4{m}^2\right)dt\hfill \end{array}\end{array}$$

(4)

subject to

$$\left\{\begin{array}{c}{}_0^C{D}_t^{\alpha }S=-{\beta }^{\alpha }(1-m){\displaystyle \frac{IS}{N}}-l{\beta }^{\alpha }(1-m){\displaystyle \frac{HS}{N}}-{{\beta }^{\prime }}^{\alpha }(1-m){\displaystyle \frac{PS}{N}}-vS,\hfill \\ {}_0^C{D}_t^{\alpha }E={\beta }^{\alpha }(1-m){\displaystyle \frac{IS}{N}}+l{\beta }^{\alpha }(1-m){\displaystyle \frac{HS}{N}}+{{\beta }^{\prime }}^{\alpha }(1-m){\displaystyle \frac{PS}{N}}-{\kappa }^{\alpha }E,\hfill \\ {}_0^C{D}_t^{\alpha }I={\kappa }^{\alpha }{\rho }_{1}E-({\gamma }_a^{\alpha }+{\gamma }_i^{\alpha })I-{\delta }_i^{\alpha }I,\hfill \\ {}_0^C{D}_t^{\alpha }P={\kappa }^{\alpha }{\rho }_{2}E-({\gamma }_a^{\alpha }+{\gamma }_i^{\alpha })P-{\delta }_p^{\alpha }P,\hfill \\ {}_0^C{D}_t^{\alpha }A={\kappa }^{\alpha }(1-{\rho }_1-{\rho }_2)E,\hfill \\ {}_0^C{D}_t^{\alpha }H={\gamma }_a^{\alpha }(I+P)-{\gamma }_r^{\alpha }H-{\delta }_h^{\alpha }H,\hfill \\ {}_0^C{D}_t^{\alpha }R={\gamma }_i^{\alpha }(I+P)+{\gamma }_r^{\alpha }H+vS,\hfill \\ {}_0^C{D}_t^{\alpha }F={\delta }_i^{\alpha }I+{\delta }_p^{\alpha }P+{\delta }_h^{\alpha }H,\hfill \end{array}\right.$$

(5)

with given initial conditions

$$\left(S\left(0\right),E\left(0\right),I\left(0\right),P\left(0\right),A\left(0\right),H\left(0\right),R\left(0\right),F\left(0\right)\right)=\left(S_0,E_0,I_0,P_0,A_0,H_0,R_0,F_0\right)⩾0.$$

(6)

Parameters $0<k_1,\dots ,k_4<+\infty$ are positive numbers that balance the size of the linear and quadratic terms in the cost, and $t_f$ is the duration of the control program. Moreover, $k_3$ and $k_4$ represent the costs of applying the control measures v and m, respectively. The set of admissible control functions is

$$\mathcal{U}=\left\{(v(·),m(·))\in L^{\infty }(0,t_f):0⩽v\left(t\right)⩽v_{max},0⩽m\left(t\right)⩽m_{max}\forall t\in (0,t_f)\right\}.$$

(7)

Pontryagin’s minimum principle (PMP) for fractional optimal control [18] is used to determine the solution of the problem. The Hamiltonian associated with our optimal control problem is given by

$$\begin{array}{cc}\hfill H& =k_{1}I+k_{2}P+k_3{v}^2+k_4{m}^2+{\xi }_1(-{\beta }^{\alpha }(1-m){\displaystyle \frac{IS}{N}}-l{\beta }^{\alpha }(1-m){\displaystyle \frac{HS}{N}}\hfill \\ \hfill & -{{\beta }^{\prime }}^{\alpha }(1-m){\displaystyle \frac{PS}{N}}-vS)+{\xi }_2({\beta }^{\alpha }(1-m){\displaystyle \frac{IS}{N}}+l{\beta }^{\alpha }(1-m){\displaystyle \frac{HS}{N}}\hfill \\ \hfill & +{{\beta }^{\prime }}^{\alpha }(1-m){\displaystyle \frac{PS}{N}}-{\kappa }^{\alpha }E)+{\xi }_3({\kappa }^{\alpha }{\rho }_{1}E-({\gamma }_a^{\alpha }+{\gamma }_i^{\alpha })I-{\delta }_i^{\alpha }I)\hfill \\ \hfill & +{\xi }_4({\kappa }^{\alpha }{\rho }_{2}E-({\gamma }_a^{\alpha }+{\gamma }_i^{\alpha })P-{\delta }_p^{\alpha }P)+{\xi }_5\left({\kappa }^{\alpha }(1-{\rho }_1-{\rho }_2)E\right)\hfill \\ \hfill & +{\xi }_6({\gamma }_a^{\alpha }(I+P)-{\gamma }_r^{\alpha }H-{\delta }_h^{\alpha }H)+{\xi }_7({\gamma }_i^{\alpha }(I+P)+{\gamma }_r^{\alpha }H+vS)\hfill \\ \hfill & +{\xi }_8({\delta }_i^{\alpha }I+{\delta }_p^{\alpha }P+{\delta }_h^{\alpha }H)\hfill \end{array}$$

(8)

and the adjoint system of PMP asserts that the co-state variables ${\xi }_i\left(t\right)$, $i=1,\dots ,8$, satisfy

$$\left\{\begin{array}{c}{}_0^C{D}_t^{\alpha }{\xi }_1\left(t^{\prime }\right)=\frac{(m-1)({\beta }^{\alpha }(I+lH)+{{\beta }^{\prime }}^{\alpha }P)({\xi }_1-{\xi }_2)}{N}+({\xi }_7-{\xi }_1)v,\hfill \\ {}_0^C{D}_t^{\alpha }{\xi }_2\left(t^{\prime }\right)={\kappa }^{\alpha }(-{\xi }_2+{\xi }_3{\rho }_1+{\xi }_4{\rho }_2-{\xi }_5({\rho }_1+{\rho }_2-1)),\hfill \\ {}_0^C{D}_t^{\alpha }{\xi }_3\left(t^{\prime }\right)=k_1-({\gamma }_a^{\alpha }+{\gamma }_i^{\alpha }){\xi }_3+{\gamma }_a^{\alpha }{\xi }_6+{\gamma }_i^{\alpha }{\xi }_7+{\delta }_i^{\alpha }({\xi }_8-{\xi }_3)+\frac{{\beta }^{\alpha }(m-1)({\xi }_1-{\xi }_2)S}{N},\hfill \\ {}_0^C{D}_t^{\alpha }{\xi }_4\left(t^{\prime }\right)=k_2-({\gamma }_a^{\alpha }+{\gamma }_i^{\alpha }){\xi }_4+{\gamma }_a^{\alpha }{\xi }_6+{\gamma }_i^{\alpha }{\xi }_7+{\delta }_p^{\alpha }({\xi }_8-{\xi }_4)+\frac{{{\beta }^{\prime }}^{\alpha }(m-1)({\xi }_1-{\xi }_2)S}{N},\hfill \\ {}_0^C{D}_t^{\alpha }{\xi }_5\left(t^{\prime }\right)={}_0^C{D}_t^{\alpha }{\xi }_7\left(t^{\prime }\right)={}_0^C{D}_t^{\alpha }{\xi }_8\left(t^{\prime }\right)=0,\hfill \\ {}_0^C{D}_t^{\alpha }{\xi }_6\left(t^{\prime }\right)={\gamma }_r^{\alpha }({\xi }_7-{\xi }_6)+{\delta }_h^{\alpha }({\xi }_8-{\xi }_6)-\frac{l{\beta }^{\alpha }(m-1)({\xi }_1-{\xi }_2)S}{N},\hfill \end{array}\right.$$

(9)

with $t^{\prime }=t_f-t$. In turn, the minimality conditions of PMP establish that the optimal controls $v^{*}$ and $m^{*}$ are given by

$$\begin{array}{cc}\hfill v^{*}\left(t\right)& =min\left\{max\left\{0,\frac{({\xi }_1\left(t\right)-{\xi }_7\left(t\right))S\left(t\right)}{2k_3}\right\},v_{max}\right\},\hfill \\ \hfill m^{*}\left(t\right)& =min\left\{max\left\{0,\frac{({\beta }^{\alpha }(I\left(t\right)+lH\left(t\right))+{{\beta }^{\prime }}^{\alpha }P\left(t\right))({\xi }_2\left(t\right)-{\xi }_1\left(t\right))S\left(t\right)}{2k_{4}N}\right\},m_{max}\right\}.\hfill \end{array}$$

(10)

Finally, according to PMP, the following transversality conditions also hold true:

$${\xi }_i\left(t_f\right)=0,i=1,\dots ,8.$$

(11)

In Section 3.4, we use the necessary optimality conditions (9)–(11) to numerically solve the optimal control problem (4)–(7), both in classical and fractional cases.

3.4. Numerical Results and Cost-Effectiveness of the Fractional Optimal Control Problem

Pontryagin Minimum Principle (PMP) is utilized to numerically solve the optimal control problem as discussed in Section 3.3. For that we use the predict-evaluate-correct-evaluate (PECE) method of Adams–Basforth–Moulton [19], implemented by us in MATLAB. Firstly, we solve System (5) by the PECE procedure, with the initial values for the state variables (6) given, and a guess for the controls over the time interval $[0,t_f]$, thus obtaining the values of the state variables. Analogously to [20], a change in variables is employed in the adjoint system and in the transversality conditions, obtaining a fractional initial value problem (IVP) from (9)–(11). The resulting IVP is also solved with the PECE algorithm, and the values of the co-state variables ${\xi }_i$, $i=1,\dots ,8$, are obtained. The controls are then updated by a convex combination of the controls of the previous iteration and the current values, computed according to (10). This procedure is iteratively repeated until the values of all the variables and the values of the controls are almost coincident with the ones of the previous iteration, that is, until convergence is achieved. The solutions of the classical model (i.e., the case $\alpha =1$) were successfully confirmed by a classical forward-backward scheme, also implemented by us in MATLAB.

As estimated in [21], we consider that 10% of infected individuals are super-spreaders.

According with [22], the initial number of asymptomatic persons is estimated to be

$$\frac{I_0+P_0}{0.15}.$$

In what follows, we assume that the total population is equal to 10,280,000 (N). Based on real data obtained from [13] and on the above assumptions, the initial conditions are $R_0=278776$, $E_0=92069$, $P_0=68208\times 0.1$, $I_0=68208\times 0.9$, $A_0=68208/0.15$, $F_0=34$, $H_0=2366$ and $S_0=N-R_0-E_0-P_0-I_0-A_0-H_0$.

The maximum number of effective daily vaccinated is estimated to be 30,000 (60,000 vaccine doses for two-dose vaccines), which corresponds to 0.003 in percentage, and this is the value of $v_{max}$.

In addition, we consider that the maximum rate of reduction in contacts is the maximum value of function $m\left(t\right)$, $m_{max}$, exhibited in Figure 2 and considered during the modelling phase.

Following [23], the relevance of the two control measures considered in the control of the disease is calculated using the sensitivity index, as presented in Definition 1. With this purpose, the basic reproduction number of the model, now incorporating the two controls, needs to be determined. For that, we look to the controls as parameters. Using the next-generation matrix method of [24], we obtain that the basic reproduction number is now given by

$$R_0=\frac{(a_h+{\gamma }_a^{\alpha }l){\beta }^{\alpha }(m-1){\rho }_1}{a_h{a}_{i}v}+\frac{(a_h{{\beta }^{\prime }}^{\alpha }+{\beta }^{\alpha }{\gamma }_a^{\alpha }l)(m-1){\rho }_2}{a_h{a}_{p}v},$$

(12)

where $a_i={\gamma }_a^{\alpha }+{\gamma }_i^{\alpha }+{\delta }_i^{\alpha }$, $a_p={\gamma }_a^{\alpha }+{\gamma }_i^{\alpha }+{\delta }_p^{\alpha }$, and $a_h={\gamma }_r^{\alpha }+{\delta }_h^{\alpha }$. The sensitivity indices are presented as functions of the control parameters in Figure 7, using the parametric values from Table 1 and the classical derivative.

The plots of Figure 7 show the following: (i) the curve of vaccination is constant and equal to $-1$, meaning that the vaccination program has a substantial impact, even for small rates of application; (ii) the curve of preventive measures rapidly moves away from zero, meaning that a preventive programme has a substantial impact only if its rate of application is high. Thus, the use of masks and the limitation of individuals in confined spaces (shops, schools, and other public spaces) have to be mandatory and should be monitored by the authorities, when it is possible.

The weights of the cost functional (4) balance the relative importance of quadratic control terms. Since super-spreaders have a greater impact in the dynamics of the disease, we consider that super-spreaders are more expensive to control than regular infected individuals. Because the preventive measures need a high rate of application to be effective, we consider that preventive measures are more expensive than vaccination. In our numerical experiments, weights are $k_1=1$, $k_2=5$, $k_3=1$, and $k_4=10$.

In Figure 8, the trajectories of the fractional optimal control problem (FOCP) for $\alpha =0.99$ are exhibited along with the solution of the classical optimal control problem (i.e., with $\alpha =1$) and the original (uncontrolled) model (2).

The corresponding Pontryagin controls are shown in Figure 9.

We can see in Figure 8 and Figure 9 that a change in the value of $\alpha$ corresponds to variations in the state and control variables. Moreover, comparing the solution of the original/uncontrolled model with the solution of the optimal control problem obtained from the application of the Pontryagin principle, we conclude that the considered control measures are effective in the management of COVID-19.

Figure 10 exhibits the efficacy function, defined in [25] by

$$E_f\left(t\right)=\frac{i\left(0\right)-i^{*}\left(t\right)}{i\left(0\right)}=1-\frac{i^{*}\left(t\right)}{i\left(0\right)},$$

(13)

where $i^{*}\left(t\right)=I^{*}\left(t\right)+P^{*}\left(t\right)$ is the optimal solution associated with the fractional optimal control, and $i\left(0\right)=I\left(0\right)+P\left(0\right)$ is the correspondent initial condition.

The efficacy function $E_f\left(t\right)$ measures the proportional variation in the number of infected individuals after the application of the control measures, $\{v^{*},m^{*}\}$, by comparing the number of infectious individuals at time t with its initial value $i\left(0\right)$.

To assess the cost and the effectiveness of the proposed fractional control measures during the intervention period, some summary measures are now presented.

The total cases averted by the intervention during the time period $t_f$ is defined by

$$AV=t_{f}i\left(0\right)-{\int }_0^{t_f}{i}^{*}\left(t\right)dt,$$

(14)

where $i^{*}\left(t\right)$ is the optimal solution associated with the fractional optimal controls, and $i\left(0\right)$ is the correspondent initial condition [25]. Note that this initial condition is obtained as the equilibrium proportion $\overline{i}$ of System (2), which is independent of time, so that $t_{f}i\left(0\right)={\int }_0^{t_f}\overline{i}dt$ represents the total infectious cases over a given period of $t_f$ days.

Effectiveness is defined as the proportion of cases averted to the total cases possible under no intervention [25]:

$$\overline{F}=\frac{AV}{i\left(0\right)t_f}=1-\frac{{\displaystyle {\int }_0^{t_f}{i}^{*}\left(t\right)dt}}{i\left(0\right)t_f}.$$

(15)

The total cost associated with the intervention is defined in [25] by

$$TC={\int }_0^{t_f}\left(C_1{v}^{*}\left(t\right)s^{*}\left(t\right)+C_2{m}^{*}\left(t\right)i^{*}\left(t\right)\right)dt,$$

(16)

where $s^{*}\left(t\right)=S^{*}\left(t\right)$, and $C_i$ corresponds to the per person unit cost of the two possible interventions: (i) vaccination at time t of susceptible individuals ($C_1$) and (ii) the implementation of preventive measures, such as the use of masks and the physical distancing of susceptible individuals ($C_2$).

Finally, the average cost-effectiveness ratio is given by

$$ACER=\frac{TC}{AV}$$

(17)

(see [25,26]).

Table 3. Summary of cost-effectiveness measures (13)–(17) for classical ($\alpha =1$) and fractional ($0<\alpha <1$) COVID-19 disease optimal control problems. Parameters according to Table 1 and $C_1=C_2=1$ in (16).

$\mathit{\alpha }$	A	$\mathit{TC}$	$\mathit{ACER}$	$\overline{\mathit{F}}$
0.99	1870.08	1116.43	0.596998	0.79967
1.0	1865.95	1114.53	0.597296	0.79791

summarizes the presented cost-effectiveness measures (13)–(17). The results clearly show the effectiveness of the controls in the reduction of COVID-19 infections and the advantage of using the fractional model.

To conclude, Figure 11 exhibits the dynamics of the infected population, $I+P+H$, and of fatalities, F, in three scenarios: (i) when the two controls are used; (ii) when only control m is used ($v=0$); and (iii) when only control v is used ($m=0$). Due to low vaccination rates, considering $v=0$, one obtains almost the same solution as the one obtained using both controls. On the other hand, erasing preventive measures leads to a serious healthcare problem. Preventive measures are in this case more effective than vaccination.

4. Conclusions

A classical compartmental model with super-spreaders was firstly proposed and applied to provide an estimation of infected individuals and deaths in Wuhan, China, in [9], and this model was later extended to the fractional-order case in order to include memory effects and better describe the realities of Spain and Portugal [6]. Unfortunately, the fractional-order model of [6] is inconsistent, in the sense that it does not satisfy appropriate time dimensions. Here, the fractional model of [6] was corrected and then used, for the first time in the literature, to model the third wave of the COVID-19 pandemic in Portugal, which occurred between 27 December 2020 and 16 February 2021. Our data fitting consisted in the minimization of the $l_2$ norm of the difference between the real values reported from the Health Authorities and the predictive cases of COVID-19 infection given by our model (2), showing that, in terms of absolute errors, the fractional-order model is better than the classical integer-order one. Another advantage of the fractional-order model was found from a sensitivity analysis, measuring the importance of each parameter in COVID-19 transmission, which allowed us to show that the sensitivity of the model decreases in absolute value with the decrease in the fractional-order $\alpha$ of differentiation, i.e., the fractional-order model is less sensitive to disturbances in the parameters than the classical one. Moreover, we introduced into the corrected model the use of vaccination and preventive control measures, investigating the use of fractional-order optimal control theory to minimize the number of COVID-19-infected individuals and reduce the associated costs. A post-optimal cost-effectiveness analysis has shown the effectiveness of controls to combat COVID-19 and the advantage of using the fractional-order model. Finally, it is shown that preventive measures are essential in the control of the pandemic.

The model investigated here is deterministic. For future work, it would be interesting to take into account the effect of noise. For integer-order models, one can refer to the works [27,28], where some COVID-19 stochastic differential equations are proposed that take into account noises derived from environmental fluctuations. Fractional stochastic models for COVID-19 are scarce and still need further investigations.