Optimal Control for Economic Development During the Pandemic

This research introduces a compartmental model designed to address the critical challenge of cost-effective pandemic modeling and analysis. The proposed framework extends the Susceptible, Exposed, Infectious, and Recovered (SEIR) model by incorporating environmental pathogen dynamics and a mortality factor within the compartmental model. To identify an optimal strategy for disease control, we formulate an optimization problem. To expedite the solution of this nonlinear optimization problem, we leverage Geometric Programming (GP), which is well-suited for handling convex optimization problems related to parameter control within the compartmental model. Additionally, we employ Particle Swarm Optimization (PSO) to explore potential solutions. Our simulation results underscore a key finding: the optimal disease control strategy is a dynamic function of time. This insight highlights the need to go beyond conventional tactics like managed isolation and quarantine, thereby improving our approach to pandemic management.


I. INTRODUCTION
COVID-19 resulted in severe economic hardship, unemployment, and a decline in Gross Domestic Product (GDP) [1], [2].Depending on how severe the social distancing measures are, how long they are implemented, and how much compliance there is, the negative economic effects may vary.As a result of the pandemic and its subsequent interventions, certain socio-demographic groups may suffer from mental health distress, economic inequality, and special difficulties [3].The solution for controlling an airborne disease must include multiple criteria [4] since the problem is deeply connected to the economic state of the system.For instance, although non-pharmaceutical interventions (NPI) may minimize the strain placed on the health care system, The associate editor coordinating the review of this manuscript and approving it for publication was Mauro Gaggero .these policies will have a direct effect on the economic status [3].
In general, finding a solution for controlling a pandemic, can be considered an optimization problem, with multiple contracting criteria, and finding the weight for different conditions or setting a restriction on the possible pressure on the health system needs a high-level knowledge of the system states.These difficulties can be abstracted in a mathematical presentation with a cost function in the optimization problem, and the weights of this function can represent the policy or the possible level of bending for the society governors in a different part of the system [5].
Providing accurate forecasts and simplicity are the most important attributes of mathematical models.It is worth noting that a complex system may consider different aspects of the disease in society, however, finding the accurate constant parameters and running this model for an appropriate amount of time needs a lot of computation power [6].Therefore, a useful solution to simulate virus behavior in the system is the utilization of the compartmental models.These models are well studied, e.g., simple SEIR (Susceptible, Exposed, Infectious, Recovered), for simulating Coronavirus Disease (COVID-19) [7].In these models, different states of the disease are modeled as compartments where each compartment represents the density of the population in that state.Flows in this model are simulated by simple differential equations [8].One of the main functions of these models is to help policymakers find and assess the best strategies and policies by providing useful information about the future state of the disease.In other words, these models can be used to evaluate the performance of a policy [9].
Many studies used different techniques including mathematical models, compartmental models, Simulation models, and Multi-criteria decision-making to analyze the dynamics of the disease and the effects of different policies.Koo et al. [10] studied the effect of isolation policies on the dynamics of COVID-19 by using a simulation model in Singapore.Gatyeni et al. [11] developed a dynamic COVID-19 model for South Africa that evaluates the existing infection prevention measures such as screening, testing, mask usage, and physical distancing by utilizing optimal control theory.Libotte et al. [12] formulated and solved an inverse problem, using the Differential Evolution algorithm to determine the parameters that define the mathematical model based on the compartmental SIR model and found the best strategy for distributing COVID-19 vaccines during the pandemic in China.Reference [13] presents a novel approach to modeling COVID-19 data using the chicken swarm optimization algorithm and a stochastic process model with the objective of extending the model's applicability by exploring different modeling methodologies, estimation methods, and fields of application.
Despite the efforts made to classify different policies, most studies do not consider the economic effects of the disease on the behavior of individuals and how these controlling policies might impact the overall economy of society [14], [15].Glover et al. [16] addressed this problem by constructing a model with a focus on the effects of the aged-based policy on the economy.To show the simultaneous effects of different isolation policies on both economic and health systems, Goldsztejn et al. [17] developed a predictive SEIR model.Rowthorn and Maciejowski [18] suggest a straightforward cost-benefit analysis that draws inspiration from optimal control theory and incorporates the SIR model of disease propagation.
The important aspect of finding controlling parameters is to search every policy and evaluate it with respect to criteria like total deaths, the pressure on the health system, or the total cost of the system.In response to this problem, Navascués et al. [19] developed a method to connect the minimization problem to the optimization over policies space, i.e. the result of the optimization represents a distribution over policy space such that minimizes the cost function.Moreover, by designing an optimization problem that can be converted into a Geometric Programming (GP) problem, which is easier to solve, Piéni and Szederkényi [5] addressed the problem of finding the optimal controlling parameter for COVID-19 hospitalizations while considering both economics and capacity of the hospital at the same time.
While many studies used the SEIR model for studying the dynamics of the disease and analyzing the economic effects of the disease, in this study a novel compartmental model, named the SEIR-PHD, is developed to include the economic aspects of the disease within a novel compartmental model.SEIR-PHD is the extended version of the SEIR model.i.e.Pathogen, Hospitalized, and Deceased states are added to enrich the simple SEIR model.The pathogen block in the SEIR-PHD is considered to clarify the dynamics of the virus by separating two different aspects (contact of individuals with infected ones or touching a contaminated surface) for keeping the virus in society.Moreover, the infectious group is divided into two subgroups to show a distinction between the symptomatic and asymptomatic cases in society.Also, the hospitalized block is added to model the additional cost of hospitalizing infected individuals with critical conditions.Additionally, by considering an upper bound on hospital capacity, defining this state helps with realizing that resources in hospitals are limited [20].
The SEIR-PHD has a simple mathematical equation governing the system and is highly nonlinear.Therefore, optimizing the control parameters over a large period with heuristic algorithms like PSO (Particle Swarm Optimization), genetic algorithm, or in general black-box optimization algorithms may result in an approximately accurate result.Hezam [21] developed a dynamic mathematical model integrating COVID-19 and Chikungunya outbreaks and solved the ordinary differential equations (ODEs) using PSO in Yemen.In another study by Zreiq et al. [22], the PSO technique was applied alongside phenomenological and compartmental models to model reported COVID-19 cases in Saudi Arabia However, the drawback of these algorithms is that they need a high amount of computation power [23].To overcome this difficulty and better analyze the economic factors in our model, the current study utilized a convex transformation.To be specific, the problem is transformed to the GP form by using Taylor approximation and refreezing each equality constraint to inequality form.Generally, the GP form can be transformed into a convex optimization problem by transforming the variables, the objective function, and the constraints.It is worth noting that an optimum point of a convex problem can be found with a fast algorithm [24] i.e., polynomial time with respect to the size of the problem.Moreover, a theorem (1) is proved to illustrate that this approximation is accurate enough and the result of the GP form can be used for the base problem.With this transformation, it is possible to find the optimized parameters in the space of all variables.An example of using GP can be found in the study by Hayhoe et al. [25] which utilizes multitask learning to determine the parameters of a compartmental discrete-time epidemic model using data from various sources.This method enables the design of control strategies for human-mobility restrictions that effectively contain the epidemic while also minimizing the economic costs of implementing non-pharmaceutical interventions.
Another difficulty in the study of these models is the calculation of constant parameters that are used for computing the state of the disease.Considering the COVID-19 virus and its responses to the policies of society, it can be inferred that the virus responds to the policy by a genetic mutation.Additionally, some factors in the model show seasonal behavior, and some depend on the availability of the healthcare capacity.Notwithstanding Piéni and Szederkényi [5] reduced the complexity of the model by considering all factors to be constant.Despite the difficulties that result from inspecting some factors to be a function of the system state within the optimization problem, in the current study, hospital admission rate and symptomatic death rate are calculated using the state of the system.
The remainder of this study is constructed as follows.Section II presents an overview of the proposed methodology including the SEIR-PHD model and its differential equations, discretization procedure, and aspects of convex optimization specifically geometric programming along with the PSO.The results of the study were presented in Section III.Section IV provides some discussions, limitations, and implications of the study.Ultimately, Section V presents conclusions and guidelines for future research.

II. METHODOLOGY A. MATHEMATICAL MODELING
As a general modeling technique, compartmental models are capable tools for the prediction of the equilibrium point in complex systems over the course of time [26].With compartmental models, the studied population is divided into compartments, and the order of the labels shows how the population moves and flows from one compartment to another [27].For predicting the number of infected, recovered, and removed individuals due to infectious diseases, a number of mathematical models have been introduced.Kermack and McKendrick [28] developed a model, SIR, which stands for Susceptible, Infected, and Recovered, which defined the dynamics of the disease and predicted its progression.Since the development of the SIR, many variations of this basic model have been introduced like SIS [29] and SIRV(Susceptible, Infectious, Recovered, Vaccinated) [30].SEIR model is the basis of the model used in this study, including an exposed period, named incubation period which is the time interval when someone is infected and becomes infected [26].Despite the high forecasting accuracy of the SEIR model [31] this model does not consider the role of pathogens in infecting individuals.Moreover, there are no blocks for analyzing hospitalized individuals in addition to  those who have died from the disease.For this reason, the compartmental model used in this study is an extension of the SEIR model named SEIR-PHD where P, H, and D represent the Pathogen, Hospitalized, and Deceased states respectively.Figure 1, which is adapted from the model in Mwalili et al. [31], represents the SEIR-PHD model and the related parameters.
Table 1 and Table 2 respectively present the Parameters and the states of the proposed model.
It is worth noting that the block P indicates the pathogen block, which supposes that the disease can be transmitted from surfaces and the virus can still be dangerous outside of the living even though this type of transmission is not very common [32].The idea of considering a separate block for the pathogen is the same as the study conducted by Mwalili et al. [31] for simulating that a susceptible individual can become infected by touching, contacting the surface, or breathing the air containing the virus.In contrast with what Péni and Szederkényi [5] proposed, in this study an individual can get infected in two scenarios.In the first scenario, an individual is exposed to the virus by having close contact with the infectious group, while in the second scenario, pathogens are responsible for exposing a susceptible individual.In the second scenario, contact with pathogens, touching the same place, or breathing air contaminated with pathogens can move the individual to the exposed state.As a result of this separation, the virus dynamics become more aligned with current intuition [33].Additionally, this separation provides a way to set control parameters based on the status of the virus and infectious population in the system.
The reason for choosing this model is that this model extends the SEIR model to include the addition of H , I S ,I A , and P blocks, which can improve the accuracy of the model.As a result of the mathematical model in this study, the ways in which the disease spreads are greatly modeled.When the infectious group is divided into two distinct groups, further restrictions can be decided in society before the virus spreads, and by adding a pathogen block, policymakers can determine the extent of cleaning the environment of the virus based on the dynamics of the virus in the environment.Moreover, the hospitalized block can help to implement the need for care in cases where patients require expensive or limited resources.Eq. ( 1) to Eq. ( 8) represent the mathematical equations of the proposed novel compartmental model.Where Eq. ( 1) shows the transition from the susceptible group to the exposed group based on two factors.In the first part, the effect of being exposed to the pathogen environment is demonstrated with the rate proportional to βν 1 P, while in the second part, the effect of infected individuals on the spread of the virus is shown with the rate βν 2 (I A +I S ).It is worth noting that ν 1 and ν 2 are controlling parameters related to the virus and social behavior.
Eq.( 2) shows the changing rate of the exposed block where the first and second factors show the incoming susceptible individuals and the last factor represents the outgoing susceptible individuals, of which there are two types, asymptomatic and symptomatic individuals.The changes in the Asymptomatic and symptomatic groups are presented by Eq. ( 3) and Eq. ( 4) respectively.There is the possibility that people in these groups can swap, which means that their health status can change from moderate to critical and this explains the flows between the symptomatic and asymptomatic individuals presented in Figure 1.Eq. ( 5) shows the group of individuals that are hospitalized due to their critical conditions.The first part of the equation shows the proportion of symptomatic individuals reaching the hospital, whereas the second and third parts represent the proportion of recovered and deceased patients leaving the hospital.Eq. ( 7) represents individuals that have recovered from the disease and come from either the hospital or the asymptomatic groups.It is assumed in Eq. ( 6) that the infectious groups, I S and I A produce pathogens with the rates η S and η A respectively.A death rate was considered for the virus since it has a limited lifetime when existing outside the host, and this rate ranges between 4 hours and 6 days for surfaces like copper [34] and stainless steel respectively [35].Moreover, Eq. ( 8) represents the group of individuals who did not survive the disease, with the first part of the equation showing hospitalized patients that die from the disease and the second part demonstrating the mortality in the symptomatic individuals.If Eq. ( 6) is not considered and the remaining equations are summed then Eq. ( 9) can be obtained, which indicates the equilibrium state, where Ṅ represents the changing rate of the individual in the system.The main hypothesis behind the Eq. ( 9) is the equality of human birth and death in the model.This assumption is reasonable since only a limited time bound can be considered in the compartmental model to smooth out the chaotic behavior of a nonlinear system.In the short period, these two rates can be approximated to zero compared to other factors included in the system.
It is worth noting that in Eq. ( 5), natural death rate and the natural birth rate are not considered.It can be inferred from this system that limited resources in the hospital, both in terms of special equipment and medical, significantly affect the system dynamics.In other words, considering constant factors for the transition of individuals between blocks is far from reality.The hospitalized block has a specific role on the individuals that die from the disease, thus it is expected that in the case of sufficient resources (doctors and beds), the dynamics of the system would be in a more controlled state.
It is expected that even if the number of patients surpasses a specific threshold, the hospital should continue servicing patients with its maximum service rate.In this case, the death rate will rise among patients who require hospital treatment but cannot get it.The hospital service rate can be modeled with a simple function as demonstrated in Figure 2.
As it is clear from Figure 2 when the number of patients is lower than a predefined maximum(H 0 ), it serves all the incoming symptomatic infectious, but when the number of patients exceeds this value, the facility only serves a fixed fraction of them.It is expected that the service level of the hospital increases when the number of patients rises but the upper bound limits the service level.However, individuals in the I S group will have a high probability of being hospitalized if the hospital can accept new patients because their probability of death decreases.Therefore it is likely that before the hospital goes to full capacity, they enter the hospitalized block.It is reasonable to assume that γ H ≫ d S and if γ H increases, d S decreases and vice-versa because it shows that hospital services to patients with critical conditions actually work and their status gradually improves.It is assumed that d S increases gradually until reaching a constant value, however, the hospitalization rate, on the contrary, starts from a constant value and declines gradually until it becomes zero.The expected plots for these two rates are illustrated in Figure 3. Furthermore, the visualization of these relaxations can be found in Figure 4.
For simplicity, an approximation of these two functions is used (inverse of a linear function).Note that the fractures in the plots can result in a singular behavior of the dynamics of the system.The basic idea is to use a smooth function for relaxing these plots.It is worth noting that functions like the Sigmoid function and inverse tangent have higher accuracy, and thus can better fit the desired function.However, considering that the parameters of the function are in the range of [0, 1], and a smooth function has a high order derivative, it can be approximated by using their Taylor approximation around zero.Additionally, it is considered that terms with order more than the product of five state blocks be zero.This results in Taylor approximation of up to three term will be sufficient for our application which can be properly modeled with the inverse of the linear function.An approximation of these rates can be obtained in Eq. ( 10) and Eq.(11).

B. DISCRETIZATION
In this section, the objective function and the methods for solving are explained in detail.To solve the proposed model, discretization is used.This is the first step towards solving the set of differential equations, in a way that is suitable for numerical evaluation on the computer.The equation ḟ = df dt = g(t, x) can be written as Eq. ( 12): Suppose that the time interval is 1 day, then Eq. ( 13) can be obtained: Eq. ( 14) to Eq. ( 20) are obtained by assuming that the number of individuals in group M in time t is denoted by X M ,t and the controlling parameter is time-dependent.

C. CONVEX OPTIMIZATION
Optimization is one of the essential parts for characterizing multi-criteria problems like COVID-19 [36].In situations where different goals are connected via a problem, optimization techniques can be used to transform the problem into an optimization problem with the cost function that is parameterized with concerning goals.solving the transformed problem results in finding weights that minimize the cost function, each weight can be interpreted as a measure of its goal importance [5].Convex optimization is a special class of mathematical optimization problems.Besides, its generality for converting the optimization problem in this format, it benefits by having fast algorithms like the Simplex algorithm for solving this class of problems in polynomial time [24].One instance of an optimization problem that can be easily converted to a convex problem is Geometric Programming (GP).This transformation is done by using a change of variables and a transformation of the objective and constraint functions [24].For defining this class properly, two types of functions need to be defined first which are demonstrated in Eq. ( 21) and Eq. ( 22): such that c > 0 and a i ∈ R is called a monomial function And the following sum i.e. the sum of the monomial function is called a polynomial function.
First, the primitive version of the optimization problem will be illustrated then it will be modified to fit the Geometric Programming conditions.The objective function will minimize 1 v 1 , 1 v 2 , the total number of deaths and patients with critical conditions i.e.I S + H + D. With the proper determination of the weights, the objective function can be described in Eq. (24).
The control parameters can be changed when the duration of L is passed with respect to the last time it is modified.This is Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
a reasonable assumption because to change these parameters, the policymaker should decide by testing the state of the society with current restriction parameters.this switch from one policy to another is assumed to take at least L units of time.
After that, the optimization problem is reformulated into a convex problem by applying Geometric Programming (GP).After reformulation, Eq. ( 25) to Eq. ( 31) can be obtained as follows.
The linear coefficients of the equations are positive except for the right side of Eq. ( 25) therefore that the equations can be re-written in the GP form.The Eq. ( 25) can be left out since X S is considered to be constant in every time interval of length L that the policy changes.Moreover, since X S is linear, if X S,0 is set equal to X S,K of every dynamic with the optimum parameters in every time interval then a strict bound would be obtained for the controlling parameters.If N is the length of the whole period, then changing the controlling parameters can occur every month (L = 1 month).For a higher accuracy, consider L to be 2 weeks, since there is a correlation between L and X S .It is worth noting that only monomial equality and posynomial inequity are accepted in geometric programming.By converting equality conditions from Eq. ( 26) to Eq. ( 31), to inequity, and the hypothesis that all the coefficients in the reminded equation are positive, the desired GP form is achieved.It needs to be proved that these changes will not affect the final result i.e. each equation takes the equality in the optimal value or if one equation is holding in inequity form, it is possible to reduce it further without affecting the cost function.As it appears some coefficients in these equations have to be approximated for fitting in the proper GP form.Taylor approximation is used, because the normalized value of each block is used, this approximation gives good confidence about the actual function value.First, the coefficient of X I S ,i in the Eq. ( 28) is approximated.
now using Taylor approximation around zero gives Eq. ( 33). 1 based on the expectation d * S ≪ γ * H ≪ 1, the coefficient related to the X H ,i is positive, rewriting equation Eq. ( 28) by using this approximation, results in Eq.34.
Eq. ( 34) is in the proper form of Geometric Programming.using same idea, for rewriting coefficient of the term X I S ,i in the Eq. ( 31) results in Eq. ( 35). 1 using this procedure for reforming Eq. ( 29) does not work out the problem, because the Taylor approximation of the fraction 1 X H ,i +α results in a negative factor of X H ,i .for resolving this issue, the constant approximation of state X H ,i is used i.e. for each L units of time, the last value of X H ,i is used for next period.by considering this, the equation can be written as illustrated in Eq. (36).
by considering the above modification, the final optimization problem is represented in Eq. (37).minimize Additional constraints are added to control the population in the hospital block.worth noting, the constant factor ϵ XH represents an upper bound on the hospital capacity.The next theorem shows that modifying the equation by converting the equality to inequity constrain, doesn't affect the final result, and in the optimum value, all inequity holds as equality.

) doesn't affect the final result of the dynamic i.e. the optimum value holds when the equality version of all inequity hold.
Proof: Consider ⃗ X * as a solution to the optimization problem 37 which some of the inequalities strictly hold.Additionally, this solution in the set of optimum solutions has the maximum number of conditions that hold inequality form.the theorem can be proved by contradiction.Worth noting that if all the inequality holds in equality form in the ⃗ X * solution then the theorem holds and there is nothing to prove, hence it is considered that at least one of the conditions in the Eq.( 20) is held strictly.The ⃗ X * can be converted to a better or the same solution(measuring the cost function) which has more inequality conditions in the exact form which contradict our hypotheses.
Note that all the coefficients in the constraints 37 are positive.now assume, for sake of simplicity first equation(X E,t * ) holds strictly for a specific time.If X E,t * is reduced to its minimum value (condition holds inequality form) without changing anything else, every other condition is proposed to be held.By reducing the X E,t * only conditions with a timeline greater than t * have a chance to not hold anymore.But, this value is reduced and the fact that the conditions are increment functions of the state parameters implies that all the altered conditions still hold.Therefore, changing the form of the equation will not affect the final state of system dynamics.The new solution ⃗ X ′ which differ with ⃗ X * in only one position i.e.X E,t * has more conditions that hold in equality form than ⃗ X * , this contradict to the hypotheses.Hence, all the inequality conditions must hold inequality form in the solution ⃗ X * .□

D. PSO
The PSO algorithm, draws inspiration from natural phenomena such as the flocking of birds and schooling of fish.This bio-inspired optimization technique is employed to discover the optimal solution for a given optimization problem [38].
In PSO, a population-based approach is adopted, wherein a collective of particles inhabits a multi-dimensional search space.Each particle within the swarm symbolizes a potential solution to the optimization problem at hand.The algorithm continuously updates the positions and velocities of these particles in a coordinated manner, with the ultimate aim of guiding them toward the most promising solution identified thus far.In this section, we outline the key components and steps of the PSO algorithm.

1) MATHEMATICAL FORMULATION
Let D be the dimensionality of the problem, and N be the number of particles in the swarm.Each particle i can be represented by its position vector X i , and velocity vector V i which both have the size of D. The objective function f (X i ) evaluates the quality of a particle's position in the search space.PSO aims to find the minimum (or maximum) of this function.The position and velocity of each particle are updated using Eq. ( 38): where: : velocity of particle i in dimension j at iteration t + 1, ω: the inertia weight, c 1 and c 2 : acceleration coefficients, r 1 and r 2 : random values between 0 and 1, Pbest ij : personal best position of particle i in dimension j, Gbest j : global best position in dimension j.
For updating the position of a particle, Eq. ( 39) is used: The termination of the PSO algorithm is determined by predefined stopping criteria, which may include reaching a maximum number of iterations, attaining a specified target fitness value, or satisfying a convergence threshold.Implementation of the PSO algorithm allows for flexibility in adjusting parameters, adopting strategies to update the global best position, managing constraints, and adapting the inertia weight.The selection of these parameters and strategies is contingent upon the nature of the problem under consideration.The pseudo-code for the PSO algorithm is provided below: Algorithm 1 Particle Swarm Optimization (PSO)

III. FINDINGS AND RESULTS
In order to consider the different values of each parameter, a specific range is used.For example, for the parameter, ν, a mean(µ) and standard derivation(σ ) are computed, and for each simulation, a uniform random ν ← [µ − σ, µ + σ ] is computed.A large number of parameters are being considered in the results within strange behaviors of the system or a variety of bifurcations, however, investigating on system response to each subspace of the parameters space is out of the context of this work.Table 3 shows the mean and standard derivation of each constant parameter in the system dynamics.
The initial conditions, which represent the state of the system at the beginning of the pandemic, are set to: {S, E, I A , I S , H , R, D, P} = {0.98,1.e −5 , 5.e −6 , 2.e −6 , 1.e −6 , 1.e −8 , 1.e −8 , 0.2} (40) Displaying the final state of each block or the maximum population of a state during the dynamic process can give valuable insights into the dynamics of the system.For example, information about the maximum number of people in the hospital gives a meaningful idea of the possible death count or at least the approximate curve of the death rate.One of the best methods for observing the final state of the system is to project the parameter space to a line, this can be done by considering all except two parameters to be constant and relating these two parameters with a linear function.Figure 5 is obtained by calculating the final state of the system when all the parameters except β 1 , β 2 are extracted from Table 3, and these two are related with the linear equation β 2 = 0.01β 1 .Besides the information about the final state of the blocks, Figure 5 shows that around β 2 ≤ 0.1 the disease-free state is the stable equilibrium.Additionally, around β 2 ≈ 0.1 a bifurcation occurred, the stable point became unstable and the system goes to another equilibrium.The discovery of these bifurcations provides insight into dynamics behavior in different parameter regimes.Figure 5 shows that the disease-free equilibrium is stable when β 2 ≤ 0.1.Figure 5b represents the same results, the disease is completely terminated from the system when β 2 ≤ 0.1, in other words, the hospital is never used when the β 2 is less than 0.1.Around β 2 ≈ 0.12 another bifurcation occurred.The straight line at the end of Figure 5b can be justified with the hypothesis that the maximum amount of the hospital capacity is set to 0.01 of all the population, and this maximum is used when the β 2 passed 0.12.Considering the case where no control parameters are used for restricting the outbreak of the disease, Figure 6 shows the progress of the dynamic.
As it appears in Figure 6a, the disease spreads very fast in the system, and all the susceptible individuals become infected within 20 days.In symptomatic, asymptomatic, and hospitalized individuals, the curve peaks around 20 days, then the hospitalized continues with maximum constant usage for about 100 days, and the symptomatic and asymptomatic states, with a high and moderate slope, respectively, both reach zero within 120 days.This is reasonable, because if the population goes up in the I S state, with a high probability they go to the hospital, therefore with a very short latency they reach their peaks.With the same argument about the connection of I A and I S , it is possible to justify the same behavior of these two curves.The slope of the dead block increased around day 20.This change is related to the percentage of the population in critical condition.It shows that sudden usage of all resources in the hospital can influence the total fatality of the system because the hospital can only serve at a bounded rate and if the symptomatic individual does not get medical service, the probability of going to the dead block is high.For simulating the system dynamic when the control parameters are considered, a lower bound is set to restrict continuous change of these parameters, to be specific, control parameters are optimized every 2 weeks based on the dynamics of the system.In time t * the parameters are optimized with respect to bound [t * , t * + 4.L].In the convex optimization method, the value of L is set to one week because the Taylor approximation begins to deviate from the actual function when considering longer timeframes.Conversely, when employing the PSO method, all the control parameters are predicted for the entire period of consideration, which spans two years, in advance.
Additionally, the cost for controlling the β 1 is set to 0.05 the cost for controlling the β 2 , because the restriction on the people's movement or quarantine overpasses the cost of synthesizing a place.Figure 7 shows the dynamics of the system when considering the control parameters for both algorithms.
When applying the PSO method, it's important to consider the randomness involved in selecting the initial values for each swarm, determining the maximum number of iterations, and setting the number of particles.In this study, 100 swarms are utilized, and the maximum number of iterations is capped at 500.It's worth noting that since this algorithm is heuristic, there's potential for achieving better results in terms of controlling the spread by increasing the number of particles or iterations.However, with the current hyperparameters, the execution time for this algorithm is approximately 2 hours on a device equipped with an Intel i7-6500U CPU running at 2.50GHz.In contrast, when the convex method is applied, the execution time is reduced to around 5 minutes.
Given that these algorithms need to be applicable in non-homogeneous systems like graph structures, computational time becomes a crucial resource.Therefore, increasing the hyperparameters for the PSO algorithm may not always be feasible.
In the convex optimization method, the disease is moderately under control.To be more specific, the susceptible population decreases gradually with a low slope and stabilizes at approximately 0.1.Conversely, when the PSO method is applied, all susceptible individuals eventually get infected, and it takes roughly 100 days for the disease to spread completely and reach a stable point, whereas the convex optimization method takes approximately 200 days to reach a stable condition of around 0.1.This significant difference in settling time has a notable impact on reducing the strain on hospitals.
In both algorithms, attempts were made to control the pathogen by setting a lower value for v 1 .However, despite these efforts, the pathogen reached its highest value in both algorithms.The crucial difference lies in the time it took for the pathogen to reach this peak value.This discrepancy suggests that the convex optimization method was more successful in regulating the rate of infection, achieving a slower increase in pathogen levels compared to the PSO algorithm.
In Figure 7b, the number of symptomatic infections increased in both cases, reaching a peak at approximately 0.2.In contrast to the PSO method, in the convex optimization method, the same peak value is attained with a gentle slope, allowing the system to gradually adapt to changes.Worth noting, that in the PSO algorithm, the discrete lines in the Symptomatic infected category indicate that the algorithm didn't quite reach its optimal values.Specifically, there is a sudden decline around the 90-day mark.Figure .7crepresents the change of the control parameters over time.In the first stage of the outbreak, the lowest value (which means a high restriction on both the pathogen and the people) was considered for v 1 in both algorithms.It is rational because it is the only control tool that can affect the susceptible block behavior, and if pathogens exist in the environment, there is a possibility that the susceptible block will make contact with the pathogen block and cause an outbreak in the system(in the first stage the probability of infection is very low, therefore it is not cost-efficient to quarantine hole system i.e. using v 2 ).As time passes, the importance of the v 2 becomes clear.This parameter directly affects the propagation of the virus, unlike the pathogen which needs an infected individual to start its spread.Hence, using a higher restriction on this parameter is reasonable.In the convex optimization method after the initial value, it decreases and reaches its minimum.On the contrary, in the PSO algorithm, these control parameters tend to fluctuate showing the randomness in this algorithm.

IV. DISCUSSION, LIMITATION, AND IMPLICATIONS
According to the outcomes of the current study, control parameters have a significant role in controlling the disease.In other words, control parameters need to be determined at the right time in order to avoid the overuse of hospital resources and prevent death due to disease.Two parameters, β 1 and β 2 , that are particularly important to spreading diseases, are controlled in this study with control parameters ν 1 and ν 2 .The results in both algorithms demonstrate that these two control parameters are a function of time, which is not surprising.Once the number of cases reaches a high point, hospitals will be under pressure and there will be an increase in the number of people who die in hospitals due to lack of service.Since these parameters are time-dependent, measures like distancing, quarantining, and cleaning the environment can be implemented more strictly when cases increase in numbers.Policymakers need time to decide which policy needs to be implemented, therefore a certain policy will be implemented for a specific period of time before it is changed and this is close to reality.Another worthwhile point is transforming the optimization problem into the GP form.In this way, the optimization problem is turned into a convex optimization problem that can be solved in a polynomial time, unlike most optimization problems which are solved in an exponential time and the final solution is not necessarily optimum.Therefore by transforming the problem into the convex optimization form, the time period for the optimization problem can be extended.This can help decision-makers reach an optimum solution in a short amount of time.Also, by using discretization, the time period can be reduced to the desired amount which helps decision-makers by providing more accurate predictions.Furthermore, the PSO algorithm as a reference for heuristic methods, highlights the superiority of the convex optimization approach in situations where time is a critical resource.
The mathematical model in this study simulates the spread of the disease by considering two critical points.First dividing the infectious group into two distinctive groups is helpful in deciding on further restrictions in society before the virus is spread.Second by adding the pathogen block which considers the dynamics of the virus in the environment and helps policymakers in determining the extent of cleaning the environment from the virus.Moreover, there is a constraint on the maximum number of patients in the hospital.This restriction can be seen in most countries when the disease reaches its peak and the number of deaths rises since there are not enough resources in the hospital to support the patients.For this reason, two rates, γ H and d s , which represent the limitations in the hospital are defined as a function of hospital capacity.Despite the presence of these limitations in society, as the results show, the disease can be controlled.The extensions of the SEIR model bring the model closer to reality and make it more practical to use.In the real world, people get infected in two ways, pathogen and infectious individuals, which is not considered in the SEIR model.Moreover, infected individuals get infected with different levels of severity, some have symptoms and require more attention while others do not show any symptoms.Some proportion of the infected individuals are hospitalized and they usually put pressure on the resources in the hospital.These blocks represent the real world and are not included in the basic SEIR model, hence by adding them decision-makers can have an accurate presentation of reality.
Considering an upper bound limitation on hospital services, in general, is an important part of analyzing the spread of disease in society.On the other hand, considering a constraint that will not allow the control parameters to take their extreme values is due to the economic nature of the problem.In this study, the main goal, controlling the spread of the disease, is defined as an economic optimization problem.According to the results, control parameters can be considered as a function of time, and the virus is still being controlled in society.This is different from the extreme restriction methods implemented in reality.Apart from the theoretical aspect, results show that the decision-makers should be aware of the different parts of the society (different blocks of the disease including the hospital, recovered, etc.) dynamically and be able to change the restrictions in a short amount of time therefore that the economical damages would be minimized while the disease is being controlled.
Although this study met its objectives, there are some limitations.In this study, some rates are considered to be a function of the state of the system.However, with a deeper look at the system, it could be seen that these rates are a function of time and the blocks of the system.Moreover, there is no bound for the Taylor approximation.To be specific, by using Taylor's approximation we obtained a convex optimization form of the problem, however, the cost of this approximation to the optimum solution didn't measure.Also, there was no exact value for some rates, hence the rates η A , and β 1 , were obtained by assuming a value for them.With regards to the weights of the cost function, these weights were considered to be constant while it is expected that the cost function should be a function of time.

V. CONCLUSION
In this study, by extending the basic SEIR model with the addition of P, H , I A , and I S blocks, the spread of the disease is modeled with higher accuracy.Moreover, some rates that are usually considered constant for convenience are assumed to be a function of other states of the system.The rates γ H and d s , which represent the death rate of hospitalized individuals and symptomatic individuals respectively, are considered as a function of hospital facilities, and the details for computing them with respect to hospital facilities function are also explained.At last, by using Taylor's approximation, an analytical function of these two rates is obtained.In this study, the optimization problem, which minimizes the cost of controlling the disease, helps in controlling the disease with minimum cost and control on the people.At last, with optimization techniques, the problem is transformed into a convex optimization problem.This transformation was beneficial in finding the optimum solution without computational limitations.Heuristic methods are time-consuming and do not ensure that the solution is optimum.Thus, the turning point in this model is the use of convex optimization and proving a theorem that shows the convex form does not change the optimum solution.
Based on the findings of this study, the following propositions may be made for future developments.The rates are a function of time and the state of the system, therefore future studies can find more accurate functions for finding the rates and also find ways to transform the problem into the convex optimization form.In this study, Taylor's approximation was used to transform the problem into the convex form, but this method didn't specify the gap between the solution and the optimum solution which was not found in this study.Finding the bifurcations can be useful in the decisions related to finding the control parameters.At some values of the rates, the disease cannot be controlled and by finding these bifurcations, the cost function used in the mathematical equations can be changed in a way to prevent entering these bifurcations that result in an overall high cost to society.In other words, the cost function can be considered a function of system rates.

FIGURE 2 .
FIGURE 2. Hospital Service rate as a function of the number of patients in the hospital.

FIGURE 3 .
FIGURE 3. Functionality of a.rate γ H to number of patients in hospital b.rate d S to number of patients in hospital.

FIGURE 4 .
FIGURE 4. approximations of a. γ H (H) function b. shows an approximation of d S (H) function.

FIGURE 5 .
FIGURE 5. a.The maximum population in the states I S , H, and total value of death caused by the disease b.The final susceptible population .

FIGURE 6 .
FIGURE 6. Dynamics of the a. S, E , I A , D, R and P population over time b.I S and H blocks over time.

FIGURE 7 .
FIGURE 7. Dynamics of the a. E , I S , I A , H, D and R population over time b.S and P population over time c. two control parameters v 1 , v 2 over time.

TABLE 1 .
The states of the SEIR-PHD model.

TABLE 2 .
Parameters of the SEIR-PHD model.

10 end 11 end 12 return Gbest;
Update particle's velocity using Eq.(38); 1 Initialize particles' positions and velocities randomly; 2 Initialize Pbest and Gbest with the best initial positions; 3 while not termination condition met do 4 for each particle do 5 Evaluate fitness; 6 Update Pbest if current position is better; 7 Update Gbest if current position is better; 8 9Update particle's position using Eq.(39);

TABLE 3 .
Parameters used in the SEIR-PHD model.