MPC Controllers in SI I R Epidemic Models

: Infectious diseases are one of the most important problems of the modern world, for example, the periodic outbreaks of coronavirus infections caused by COVID-19, inﬂuenza, and many other respiratory diseases have signiﬁcantly affected the economics of many countries. Hence, it is therefore important to minimize the economic damage, which includes both loss of work and treatment costs, quarantine costs, etc. Recent studies have presented many different models describing the dynamics of virus spread, which help to analyze the epidemic outbreaks. In the current work we focus on ﬁnding solutions that are robust to noise and take into account the dynamics of future changes in the process. We extend previous results by using a nonlinear model-predictive-control (MPC) controller to ﬁnd effective controls. MPC is a computational mathematical method used in dynamically controlled systems with observations to ﬁnd effective controls.


Introduction
Infectious diseases are some of the modern world's most important problems, e.g., influenza, the periodic outbreaks of coronavirus infections caused by COVID-19, and many other respiratory diseases have significantly affected the economies of many countries (lockdowns, masks, and reduced production efficiency) [1,2]. It is therefore important to minimize the economic damage, which includes both loss of work and treatment costs. Recent studies have presented many different models describing the dynamics of virus spread, which help to analyze and explain epidemic outbreaks [3,4]. Previously, it was shown that several infections can circulate in the population at the same time. In this case, the controlled Susceptible-Infected-Recovered model with multiple infected subgroups can be used [5]. The solution of this model is usually found by using classical optimality criteria, e.g., the Pontryagin's maximum principle [6]. In this study, we focus on finding solutions that are robust to noise and take into account the dynamics of future changes in the process. We extend previous results by using a nonlinear model-predictive-control (MPC) controller to find effective controls [7]. MPC is a computational mathematical method used in dynamically controlled systems with observations to find effective controls.
It has been shown that several virus species, or several different virus strains of the same species, can circulate simultaneously in a population. In this case, the subpopulation of infected individuals can be divided into several subgroups. In our study, two types of virus were considered to simplify the calculations, but the results obtained can be generalized to a larger number of circulating strains.
Experience with influenza epidemics and the SARS-CoV-2 pandemic has shown that there can be several waves of illness in a population during a single epidemic season. To reduce the morbidity of an epidemic wave, the government and health authorities apply various measures to prevent and contain the epidemic. The use of different preventive and curative measures during a single wave reduces the number of cases and the rate at which the virus spreads in the population. Thus, the rise of the disease in the next wave of the epidemic begins with new initial states of the system. A review of the literature led to the idea of using predictive control models to develop effective control that takes into account changes in the environment. MPC controllers allow for the resulting control to be adapted by dividing the entire time interval into smaller sub-intervals. At each interval, the control system is checked to see if the constraints are satisfied, and if not, the control action is adjusted.
The main applications of this method are in the petrochemical [8], woodworking [9], and energy [10,11] industries, where, among other things, the stability of the obtained controllers and the possibility of their automatic adjustment are important. The MPC controller finds a solution in the form of piecemeal functions on the control interval, but the optimization problem is solved over a larger time interval-the forecast horizon. The MPC controller finds the solution in the form of piecewise constant functions on the control interval, but the optimization problem is solved on the forecast horizon.
As the main contribution of this work, we formulated a control problem and solved it using an MPC controller, a series of numerical experiments are carried out to confirm the effectiveness of the solutions obtained by the MPC controller for the SIIR model.
The paper is organized as follows. Section 2 presents the epidemic models and modified infection rate parameter. Section 3 shows the application of MPC-controllers to a SIIR epidemic model. In Section 4, we use numerical simulation to illustrate our results. The paper is concluded in Section 5.

Epidemic Process in Total Urban Population
Our model is based on the extended Susceptible-Infected-Recovered model with two types of viruses, which circulate simultaneously in an urban population [6]. A size of population is equal to N. According to these assumptions, we can consider four groups: the Susceptible, the Infected by virus V 1 , the Infected by virus V 2 , and the Recovered group. Susceptible is a group of people who are not infected by viruses, but can be infected by one or both types of viruses and who have no immunity to the viruses. Depending on the viral strength, we observe that the number of people infected by virus 1 or by virus 2 can be different, and people infected by virus V 1 or V 2 belong to the subgroup Infected. The Recovered subgroup consists of people who have recovered from the infection. The mixing of populations allows for viruses to spread quickly, and each person in the population is assumed to have an equal probability of coming into contact with others. Therefore, when an infected person interacts with a susceptible person, the virus can spread. A virus with higher virulence, we assume, will have a higher probability of success in spreading when an encounter occurs between the infected and the susceptible. The scheme of process is represented in the Figure 1.
We model the spread of the virus in the population based on the classical SIR epidemiological model, where a system of differential equations is used to describe the proportion of the population as a function of time. Then, at time t, n S (t), n I 1 (t), n I 2 (t), and n R (t) correspond to the fractions of the population that are Susceptible, Infected by virus v 1 , Infected by virus v 2 , and Recovered, respectively, and for all t the condition N as proportions of Susceptible, Infected, and Recovered. At the beginning of the epidemic t = 0, most people in the population belong to the susceptible subpopulation, a small group in the total population is infected and other people are in the recovered subpopulation. Therefore, initial states are 0 < S(0) = S 0 < 1, 0 < I 1 (0) = I 0 1 < 1, We have extended the simple SIR model introduced by [12] to describe the situation with two virus types:Ṡ where δ i are infection rates for virus V i , i = 1, 2, σ i are recovery rates (see Figure 1). Infection rate is defined as a product of the contact rate l and transmissibility of infection, i.e., probability of transmission infection during the contact, δ i 0 Infection rate is integrated into the evolution of mutation process to epidemics in the urban population. Medical treatment or quarantine reduces the number of the infected individuals in the urban population. These prevention measures can be interpreted as control parameters in the system defined as u = (u 1 , u 2 ); here, u i is a fraction of the infected which are quarantined or under intensive medical treatment. Recovered rates are inversely proportional to The objective function: We will minimize the overall cost in time interval [0, T] [13][14][15][16]. At any given t following costs exist in the system: These costs can be interpreted as the price the government or medical organization pays for the health system to treat the infected and increase the probability of recovery for the treated individuals. Infected individuals who are not detected because they are asymptomatic are not treated, and therefore there are no direct costs for them. • The treatment or medical cure of those infected produces costs for society. The functions h i (u i (t)) represents the costs of treating the infected to keep them out of the labor market until recovery. The function h i (u i (t)) is increasing on the argument u i (t), is twice differentiable, and is such that h i (u i (0)) = 0. • The function g(R(t)) is a benefit rate. These are the bonuses that society receives from Recovered individuals who have returned to work and started contributing to society. The function g(t) is non-decreasing and differentiable and g(0) = 0.
• Constants e I 1 , e I 2 , e R represent the costs and benefit for Infected and Recovered in the end of the epidemic.
Aggregated system costs is given by and the optimal control problem is to minimize the cost, i.e., min To simplify the analysis, we consider the case where e I 1 = e I 1 = e R = 0.
Here, we assume that as a control parameters which help to reduce numbers of infected in human population, medical treatment or isolation can be considered, then define variable u = (u 1 , u 2 ) as control 0 ≤ u 1 (t) ≤ 1, 0 ≤ u 2 (t) ≤ 1, for all t.

SIIRS Models
The SIIR models discussed in the previous section are best suited to describe epidemic processes in which viruses mutate slowly. Because of the small number of mutations, individuals can develop lifelong immunity to the virus. However, there are viruses that mutate rapidly and for these, SIIRS (Susceptible-Infected-Recovered-Susceptible) models have been developed (see Figure 2) to account for the decline of immunity in the population over time.
The system of differential equations describing the dynamics of state change is equal to Here, as in the SIIR model, the parameters δ i , σ i , i = 1, 2 are responsible for the intensity of infection and cure, u i -the proportion of infected people on medication, and γ -the intensity of immune decline in the population. The objective function for such tasks is no different from (2). Note that due to the gradual decline of immunity, the epidemic process in the SIIRS problem can continue indefinitely. For this reason, when solving this problems we should choose large values for T. However, in such a case there is a risk (due to the high mutation rate of viruses) that the model parameters will no longer be relevant. Therefore, methods that can adapt to changes in the process by predicting its behavior can be effective in solving such problems.

SWIIRS Model
In 2020, the SWIIRS (Susceptible-Warned-Infected-Recovered-Susceptible) model was proposed, which takes into account citizens' awareness of the epidemic [17]. This model divides the population into five groups instead of four. An additional group of individuals is the group of warned W individuals who are aware of the virus (see Figure 3). This approach allows a more accurate description of epidemic processes in large countries, as the speed of spread of the virus often depends not only on the characteristics of the virus itself, but also on the behavior of members of the population. The changes in population groups for the controlled SWIIRS model are described by a system of five differential equations: In the system (4), the constants δ S i and δ W i describe the infection rates for groups S and W, respectively, and σ i -the cure rate i = 1, 2. The value σ 3 is responsible for the probability that a knowledgeable individual will acquire immunity to viruses (e.g., go to a health facility to get a vaccine). The intensity of the decline in immunity is described by γ. Two parameters are responsible for informing (warning) members of the population:

•
Parameter u 3 -percentage of susceptible individuals informed by the media (become Warned); • Parameter η-probability of informing the susceptible individual when communicating with the warned individual.

Controllable Infection Rate Parameter
Consider a model (1) with one virus I i and let p i be the probability of virus transmission V i from a diseased individual to a susceptible one, and l be the average number of contacts of an individual (with other members of the population) per unit of time. Consider a random contact of a susceptible individual. Infection with the virus will only occur if two events occur simultaneously: the contact was with the infected individual, and the virus was transmitted on contact.
Define two incompatible events: • Event Aa susceptible individual becomes infected through one of an l contacts; • Event B-with a similar number of contacts, infection will not occur.
Obviously that events A and B form a complete group of events (see Figure 4), then The event B will only occur if there is no infection at any of the l contacts. In other words From Equations (5) and (6) we see that the probability of transmission to a susceptible individual by l contact with infected individuals can be defined as The same formula can be obtained in another way. The probability of the event A is given by This probability p(A) is called the modified infection rate parameter and it is denoted by δ i (I i (t), l).
Let us redefine u i (·) (intensive virus treatment) as u i 1 (·). Let us introduce two types of control into the resulting system: isolation of infected individuals and reduction in the number of contacts in the whole system. Let u i 2 (·) be the fraction of infected individu-als isolated, and u l (·) be the fraction of contacts prevented by reducing the mobility of individuals in the whole population. In this case we obtain where 0 ≤ u i 2 (t), u l (t) ≤ 1. Then the system of differential equations of the corresponding single-variable SIR problem is If we perform similar transformations to (8) for the SIIR dual virus model, the calculated formulae for the controllable infection parameters are as follows where i = 1, 2. If we substitute these parameters into the system (1), the differential equations take the following form: For the introduced controls we also need to redefine the cost functions. Let h 1 2 (u 1 2 , I 1 ), h 2 2 (u 2 2 , I 2 ), h l (u l ) be non-reducing functions corresponding to the costs of the controls u 1 2 , u 2 2 , and u l . Then the function (2) takes the form

MPC Controllers
Model predictive control (MPC) is a set of methods for finding effective control based on model prediction. The basic idea of MPC is to estimate the future behavior of the controlled system at a finite time interval and to compute an effective control signal that minimizes a given objective function, given the constraints on the system (see Figure 5). The next algorithm of the model-predictive-control method includes the following steps: 1.
Modeling of an object or process using mathematical and computer modeling techniques; 2.
Defining the main objectives to be achieved by management, selecting the forecasting horizon; 3.
Predicting the behavior of the process under the influence of the control signal; 4.
Find an effective control signal, taking into account the full range of constraints imposed on the control and manipulated variables; 5.
Select the control horizon in which the effective control found is applied; 6.
Measuring the actual process states at the end of the control horizon, which are taken as the new initial conditions; 7.
Shift the forecast horizon by the size of the management horizon and repeat steps 3-6.

Linear Prediction Model
We consider a linear controllable and observable finite-difference system: where A(t), B(t), C(t) are some matrix functions of dimensions n × n, n × m, and n × k, respectively, U is set of valid controls (is a subset of piecemeal functions), and [0, T] is a given time interval.
Denote by x r , u r the reference values of state and control variables. The MPC is a computational tool that finds the minimum either of its built-in target functionality (15) or of a custom function that the user must describe. Efficient control is determined using numerical Knows-What-It-Knows algorithms (KWIK) [18].
where w x and w u weight vectors for variables x and u. However, the system (1) should be transformed to have a possibility to apply the MPC controller. In our case we find a control u(t) that provides the minimum to the objective function (2). The system (1) must be linearized before the MPC can be used. The predicted values are determined by the formula: • x(t), u(t)-the known current states and controls of the system, The solution of the obtained problem is carried out using standard methods for solving quadratic programming problems. Let solution of the problem (15), then v * (t + 1) will be called the effective control for the original system (14) on the interval [t, t + 1] obtained using MPC method. After applying the control, we obtain the new initial data x(t + 1) and u(t + 1) to make the following forecast and search for an effective control on the interval [t + 1, t + 2].

Linearization of SI 1 I 2 R Problem
In the current section, we discuss the application of the MPC controller to the extended SI 1 I 2 R system [19].
Suppose, that the initial states for all states S(0), I 1 (0), I 2 (0), R(0) are known. Otherwise, we can always measure it. Define SI IR(0) as a vector of initial states.
Another assumption is that the initial control u(0) = u 0 is defined. We divide the initial period of time [0, T] into k equal parts, then the objective function (2) by the integral property can be written as: Here, the total length of time interval is divided into several segments The initial conditions for the first iteration are defined in (18), for each subsequent problem they can be found from the previous task as: (1) is also valid for initial values, i.e., Consider some disturbance of this system: where ∆u(t)-control variable variation vector.
Then,ẋ i (t) + ∆ẋ(t) = f ( x, u, t). The right hand sides of system (22) is decomposed into a Taylor row:ẋ Here, O((∆x) 2 , (∆u) 2 ) is the second-order residual term of smallness. Subtract from (22) and obtain∆ As a result, the linearized controlled system is obtained. We substitute the solution (x, u) of the linearized system into the original model (1), (20) and obtain the following value of the objective function J i (x, u) = J i . Let (x * , u * ) be a solution that provides a minimum to the functional (20), then we can denote by an error in the solution the value of the ε i , which is determined by the formula: Going to the limit k → ∞, we have: This means that the greater the number of split points, the more accurate the solution will be.
We have conducted a series of numerical simulations to check whether the MPC controller solution converges to the optimal solution as the parameter k increases (the other parameters remain the same). It can be seen that all the costs (standard, efficient, and optimal) increase as k increases. This is due to the fact that the accuracy of the integral calculation in the function increases. In test no. 2, the value of the effective costs is extremely high, due to the fact that it is recommended to use the value k ≥ T. Starting from test no. 5, the values of effective costs and optimal costs differ insignificantly, so in practice we use the values k = 2T or k = 4T. However, the program execution time in the second case (test no. 9) is 10 times longer than in the first case (test no. 5). Figure 6 shows the work of the nonlinear MPC controller [20]. The MPC controller construct the efficient control (u * ) i and predictive values of the objective function using a vector of outputs y i of the system. The main difference between nonlinear and linear MPC controller is the linearization of the optimization problem before the application of the solver [21,22]. The received control is inputted to the dynamic system, in the next iteration MPCcontroller obtains another output. The process of constructing a solution is complete when the control horizon is reached. The final effective control strategy is chosen to minimize the objective function on the predictive horizon. Since numerical algorithms are used to find the solution, the resulting effective value can be suboptimal.

Numerical Simulations
This section presents the results of a numerical simulation that corroborates our theoretical results. The numerical experiments are aimed at obtaining effective treatment strategies using nonlinear MPC controllers. We selected MATLAB as our platform for programming and numerical calculations for a series of experiments. MATLAB offers an extensive library of the model-predictive-control method, including various examples. Nonlinear-MPC controllers were employed in our study. These controllers utilize the linearization algorithm described in the article to find solutions. Nonlinear MPC is adequate to address these issues. However, larger systems may demand more precise predictive models based on machine learning and neural networks. For first three experiments, two variants of the SARS-CoV-2 virus were selected for study: Alpha [23] and Omicron [24]. The epidemic situation was considered for a period of 180 days (6 months). During this period, six MPC controllers were used sequentially with a forecast horizon of 30 days (1 month) and a control horizon of 21 days (3 weeks). Time dependence of MPC controller accuracy and performance on k parameter are shown in Table 1. The values of the main system parameters are shown in Table 2.
The serial application of multiple nonlinear MPC controllers allows the constraints, system parameters and objective function coefficients to be adjusted step by step. To simplify the calculations, we will use the same parameters in each iteration of the process.

Experiment 1.
In the first experiment, we considered the epidemic process with a small number of contacts, which is possible in cases where most of the infected comply with the quarantine. It can be seen that applying control is having a positive impact on stopping the virus spreading (Figure 7). On the effective control graph (Figure 8), you can notice "jumps" due to the fact that a new MPC controller is connected each new month (30th day). However, the resulting control greatly reduces the damage from the epidemic compared to a situation where we do not control the epidemic process. So the cost in the uncontrolled case is J 0 = 3159 monetary units (m.u.) and with the effective strategy J e f f = 1016 m.u. The benefit of using control is therefore 2143 m.u.

Experiment 2.
A second experiment (Figures 9 and 10) was run with the average number of contacts per day k = 3. This experiment is the most different from the others. The main wave of the epidemic does not occur in the first few days, but only at the beginning of the third month. The highest costs are in the same period. After t = 60, the effective control is similar to that of experiments 1 and 3. With this strategy, it is possible to reduce the cost of the epidemic from J 0 = 3613 m.u. to J e f f = 1045 m.u., which saves 2568 monetary units. It can also be seen that as the average number of contacts increased, so did the values of the objective functions. Just as importantly, the benefits of using MPC controllers have also increased.

Experiment 3.
In the third simulation the number of contacts is k = 10, which approximates the real situation. Figure 11 shows that with a large number of contacts, almost all agents are infected by t = 5. However, when efficient control is used, the infection peaks are significantly reduced (from 0.63 to 0.29 and from 0.26 to 0.14). As in the previous experiments, the use of efficient control is economically beneficial, reducing costs from J 0 = 4282 m.u.
to J e f f = 860 m.u. The benefit of using an effective strategy is 3422 m.u. (see Figure 12).  It should also be noted that despite the largest value of k the costs of using effective management in Experiment 3 turned out to be the smallest. This may be due to the fact that when solving the problem, they are linearized and the accuracy of the solution can vary from the parameters of the problem. At the same time, all three experiments showed that effective control can reduce the economic damage from the epidemic, which shows the practical importance of using the model predictive control in SIIR problem. Experiment 4. In the fourth experiment we examined the SIIR model with the data taken from the paper "Optimal control of heterogeneous mutating viruses" [5]. Table 3 shows the parameter values used in the simulation.
Cost of treatment u 2 Initial states  Figure 13 shows the population states in the uncontrolled and controlled cases. It can be seen that the use of the control found with MPC significantly reduces the number of infections. In the uncontrolled case, the maximum proportion of the number of infections is I 1 max = 0.594, I 2 max = 0.386, and in the controlled case, I 1 max = 0.32, I 2 max = 0.18, corresponding to the initial proportion of infected in the population. After the computer simulation, the solution of the formulated problem was obtained and the values of the functional were found (see Figure 14).
The strategies that increased to u 1 = 0.26 and u 2 = 0.24 during the first three days of the epidemic and then decreased to zero were found to be effective. The resulting control signals are similar in appearance to the continuous controls presented in the experiments of this paper. The aggregated costs are equal to 41.6 m.u. in the controlled case and 227.51 m.u. in the uncontrolled case. Thus, the benefit of using the control is 185.91 monetary units, which is slightly less than using the optimal control 189.17 m.u. In Figure 15, we can see the structure of the optimal control for the model and the data from this experiment. It can be observed that the control structure in the case of effective control is similar to the structure in the case of optimal control.  Experiment 4 showed that the developed software package finds near-optimal solutions to classical SIIR problems and can be used as an alternative to classical optimal control methods. Experiment 5. As the software implementation of the MPC method showed good accuracy for standard tasks, Experiment 5 was conducted using a modified SIIRS model (with controlled infection rate parameter). Data for the experiment were obtained from sources such as the World Health Organization, Our World in Data, and Worldometer. The parameters of the experiment are shown in Table 4.
Using effective strategies (see Figure 16), it was possible to reduce the maximum proportion infected with V 2 virus from 0.6 to 0.3, but the proportion infected with V 1 virus increased from 0.34 to 0.38, which is due to the fact that the proportion of people (30% of the population) who were not infected with V 2 virus (due to the measures introduced) managed to become infected with V 1 virus. It can also be seen that from t = 15, the proportion infected with both viruses does not exceed 1%.   By running the computer simulation program, effective controls for the formulated problem were obtained (see Figure 17). It can be observed that the treatment of individuals for both viruses takes a maximum value starting from t = 1 and decreases for virus V 1 after t = 21 and for virus V 1 after t = 18. For V 1 , a small fraction (u 1 2 = 0.1) of those infected during [5,8] are isolated, and after t = 15 the control u 1 2 turns on again, reaching a maximum value, and turns off after t = 24. For V 2 , maximum isolation of the infected occurs at [0, 4], [5,8], [12,21], and partial isolation at [9,11], [29,30]. Reducing the total number of contacts was most effective in periods [15,25]. For problems with a modified infection rate parameter, no optimality conditions have yet been found, so the only way to assess the quality of the solution is the cost savings provided by the constructed efficient control. For this reason, the values of the target functional and the total control for the two viruses were calculated for this experiment (see Figure 18). With the combined effect of all controls it was possible to significantly reduce the value of the target function from 33,286 to 11,180 m.u. The benefit of using effective controls was therefore 22,106 m.u. Experiment 6. The SWIIRS model was used for Experiment 6, and the main objective was to construct long term strategies (T ≥ 90) using multiple MPC controllers. The experimental parameter values are available in Table 5.  Figures 19 and 20 show the results of Experiment 6. The effective control plots show that when several MPC controllers are connected, there are pulses (switching points) in the control. They are related to the fact that the prediction model is recalculated to a new prediction interval. You may notice that from the third controller (t = 60) the constructed controllers have a similar structure. This is reflected both in the behavior of the process and in the behavior of the total cost function (they become periodic).  By applying efficient management, the value of the target function was significantly reduced from 2918.8 to 177.8 m.u. Thus, the benefit of applying management was 2748 monetary units. It can be seen that model-predictive-control solutions remain effective over long time intervals, but require the connection of several MPC controllers.

Conclusions
In this paper, we have investigated the application of MPC controllers to an extended SIR model with two different viruses. Model predictive control is one of the ways to achieve efficient control in dynamically controlled systems. We linearized the system of differential equations describing the dynamics of the viruses and showed that as the linearization step is reduced, the accuracy of the solution increases. We carried out a series of numerical experiments in which we obtained effective strategies whose form is close to the optimal strategies obtained by Pontryagin's maximum principle. The experiments also showed that the effective control has a similar form for different numbers of contacts. It can be seen that efficient control is much simpler and faster than optimal control and uses fewer resources. This method can also be used to assess the situation during the spread of epidemics, when time is of the essence.

Conflicts of Interest:
The authors declare no conflict of interest.