The COVID-19 epidemic started in Wuhan, China, at the end of 2019, was declared to be a pandemic by the World Health Organization (WHO) on March 11, 2020. The most common symptoms of the disease include fever, cough, fatigue, shortness of breath, and loss of the sense of smell, where complications may include pneumonia and respiratory distress known as a severe acute respiratory syndrome (SARS). At the time of writing this paper, the primary mode of treatment has been symptomatic and supportive therapy [1, 2] , and there has been no approved vaccine or specific antiviral treatment.
The pandemic shook the economy of the whole world with a significant reduction of exports, a decline in tourism, mass unemployment, and business closures [3] . Governments and health authorities worldwide have responded by implementing non-pharmaceutical intervention (NPI) policies that include travel restrictions, lockdown strategies, social distancing measures, workplace hazard controls, closing down of schools and workplaces, curfew strategies, and cancellation of public events. Many countries have also upgraded existing infrastructure and personnel to increase testing capabilities and facilities for focused isolation. The public has been instructed to wash hands several times a day, cover mouth and nose when coughing or sneezing, maintain a certain physical distance from other people, wear a face mask in public places/gatherings, and to monitor and self-isolate if one suspects herself/himself to be infected. The extent to which such policies and measures have been implemented in a certain country is called the stringency index of its government's response [4, 5] . Each government responded in its capacity to find a suitable balance between saving lives and saving livelihoods, which [6] termed as a problem of health versus wealth. Livelihoods can be saved through the implementation of suitable relief and recovery measures for people and small businesses. On the other hand, lives can be saved through the implementation of good testing and NPI policies. In other words, there is a direct relation between the stringency index of the government and 'saving lives'.
All the above strategies and policies implemented by the governments are considered to be the control mechanisms for the epidemic. Such strategies fall under two categories: mitigation and suppression [7, 8] . The mitigation strategies slow down the rate of transmission of disease or, in other words, 'flatten the curve.' However, they do not necessarily stop the spreading of the disease, which is the goal of suppression strategies. Given the required objectives (e.g., minimizing the number of deaths caused by the epidemic) and constraints (e.g., socio-economic costs), the problem of finding optimal strategies has been recently studied through the framework of optimal control theory.
To understand, predict, and control the evolution of the COVID-19 epidemic, a huge effort has been devoted by the researchers to design models as accurate and as effective as possible. Each model, by and large, is a variant and/or an extension of SIR (susceptible, infected, recovered) and SEIR (susceptible, exposed, infected, recovered) models, which describe the flow of population through three or four mutually exclusive stages of infection, respectively (see [9] and [10] for a comprehensive review). These basic models have few parameters that are easy to identify [11] , and are considered as population models that view the epidemic from the macroscopic perspective. This is in contrast with the approaches that capture heterogeneity of population structure such as network epidemic models [12, 13, 14] or metapopulation epidemic models [15, 16, 17] , that view the epidemic from the microscopic perspective. In what follows, however, we study the epidemic through the macroscopic perspective of population models.
Following the outbreak of COVID- 19 , there has been an effort to produce comprehensive population models with a focus on different facets of the epidemic. Such models are more complex than simple SIR and SEIR, and include several intermediate stages that accurately portray the dynamics of the epidemic. For instance, [18] develops an extension of the SEIR model that incorporates the governmental actions (e.g., preventive measures and restrictions) and the individual behavioral reactions, whereas [19] develops an extension of the SIR model that incorporates the number of deaths due to the epidemic. Another quite interesting model is the one developed in [20] that considers an eightcompartment model called SIDARTHE, which includes eight stages of infection: susceptible (S), infected (I), diagnosed (D), ailing (A), recognized (R), threatened (T), healed (H), and extinct (E). A distinguishing feature of this model is that it differentiates between the infected individuals based on the severity of their symptoms and whether they are diagnosed by a health authority. It is crucial, as also emphasized in [21, 22] , to differentiate between diagnosed and undiagnosed individuals because the former are typically isolated and are less likely to spread the infection. Similar models have been adopted and extended to study optimal control policies for the epidemic such as the implementation of social distancing measures [23, 24, 25] , lockdown strategies [26, 27, 28] , and heterogeneous policy responses based on age-groups [29, 30] .
In addition to the above NPI strategies, testing and isolating the infected population from the susceptible population is one of the most important strategies to control the epidemic spread. For instance, it has been reported that COVID-19 was eliminated from the Italian village Vo'Euganeo through testing both symptomatic and asymptomatic cases [31, 32] . Moreover, on his media briefing 1 of March 16, 2020, Dr. Tedros Adhanom Ghebreyesus, the Director-General of WHO, urged the following:
"Social distancing measures can help to reduce transmission and enable health systems to cope. Hand-washing and coughing into your elbow can reduce the risk for yourself and others. But on their own, they are not enough to extinguish this pandemic. It's the combination that makes the difference. As I keep saying, all countries must take a comprehensive approach. But the most effective way to pre-vent infections and save lives is breaking the chains of transmission. And to do that, you must test and isolate. You cannot fight a fire blindfolded. And we cannot stop this pandemic if we don't know who is infected. We have a simple message for all countries: TEST, TEST, TEST." COVID-19 can be detected through two types of tests known as type-1 (RT-PCR) and type-2 (serology). In the type-1 test, a swab is inserted into the subject's nose to qualitatively detect nucleic acid from SARS-CoV-2 in the upper and lower respiratory specimens [33] , which enables one to detect whether the subject is currently infected with COVID-19. Type-2 test, on the other hand, is a serum test in order to detect relevant antibodies, which enables one to know whether the subject was infected in the past with COVID-19 and now he/she has recovered. Both types of tests are important in the control of an epidemic. Type-1 tests help to limit the disease spread by the identification of infected individuals and their contact tracing [34] . Type-2 tests, on the other hand, are useful in reducing the size of the testable population for type-1 tests [35] that helps to increase the testing specificity. However, the type-1 test, up to now, is considered to be the only recommended method for the identification and laboratory confirmation of COVID-19 cases according to the WHO [36] . Moreover, only type-1 tests can provide information in real-time related to describe the outburst of the epidemic, which is the reason that the datasets related to testing only include type-1 tests 2 .
Following the recommendation of the WHO director, with different levels of setups, many governments increased their testing capacities, while others feared the economic burden of intensive testing policy. However, [37, 38] show that such a burden is only short-term and, on the contrary, intensive testing reduces the 'overall' cost of the epidemic in the long run because it enables the government to gain rapid control of the epidemic and revive the economy of a country. Testing enables the health authority to identify and isolate the infected people from the susceptible population, which limits the transmission of the disease. Therefore, testing is considered to be a crucial control mechanism for the epidemic [39] . However, few attempts have been dedicated to study the testing policy for an epidemic from a control-theoretic perspective.
In somewhat similar to a resource allocation problem [40, 41] in epidemic control, [42] poses the optimal test allocation as a well-known sensor selection problem in control theory, whereas [43] poses it as a welfare maximization problem by considering specificity and sensitivity of tests. The main assumption in these papers, however, is the availability of information portfolios of all individuals in a society, which enables the decision makers to compute the infection probability of individuals and utility loss for each individual in case of decision errors. On the other hand, [44, 45] study the problem of testing policy from an economic perspective, where the goal is to find an optimal testing pol-icy that minimizes the total number of quarantined people to incur minimal cost on the economic activity of a country while also mitigating the epidemic spread. Without such a testing policy, governments usually resort to indiscriminate quarantining of people that burdens the economy of a country without any reason. Therefore, testing allows to identify and isolate the positive cases to allow for case-dependent quarantining. Another aspect of testing policy is studied in [46] , which numerically computes an optimal trade-off between testing effort and lockdown intervention under the constraint of limited Intensive Care Units (ICU).
We introduce a new model named SIDUR model -susceptible (S), undiagnosed infected (I), diagnosed infected (D), unidentified recovered (U), and identified removed (R) -to study the control of an epidemic through testing. Similar to [20, 21, 22] , we differentiate between the undiagnosed and diagnosed infected population. We assume that the diagnosed infected population are either quarantined and/or hospitalised and only the undiagnosed infected population is responsible for the disease transmission to the susceptible population. The identified removed population consists of people who recover or die after being diagnosed and the unidentified recovered population consists of people who recover without getting diagnosed.
The control input in the SIDUR model is defined as the number of tests performed per day, where the influence of the control is directly linked with the testing specificity. The testing specificity determines the probability of detecting an undiagnosed infected person through a test, which, for instance, can be increased through efficient contact tracing. As stated previously, only type-1 tests can provide information in real-time related to describe the outburst of the epidemic. Therefore, the control input in the SIDUR model only accounts for the type-1 (RT-PCR) tests.
In this paper, we consider the COVID-19 case of France as a benchmark example. That is, we estimate and validate the model on French COVID-19 data. Then, we propose two static testing policies in order to control the epidemic: 1) Best effort strategy for testing (BEST) and 2) Constant optimal strategy for testing (COST).
The best-effort strategy for testing (BEST) is a suppression strategy for an epidemic that provides the minimum number of tests to be performed per day in order to stop the epidemic spread. Thus, BEST is meaningful only during the spreading phase of the disease. We provide an algorithm to compute the number of tests required by BEST policy. Since BEST is a suppression strategy that stops the epidemic growth immediately, it usually requires a lot of tests to be performed per day. However, it requires less number of tests if implemented sooner, which is illustrated for the case of France by plotting the number of tests required by BEST with respect to time.
The constant optimal strategy for testing (COST) is a mitigation strategy when the total stockpile of tests is limited. It provides the optimal number of tests per day that must be allocated in a certain time interval from the onset of the epidemic in order to minimize the peak of infected population. When the stockpile finishes at the terminal time, then the tests are not performed anymore. The main idea of COST is the consideration of two peaks of the epidemic; one that occurs before the stockpile of tests finishes and one that occurs after the stockpile finishes. If the stockpile of tests is allocated constantly such that it reduces the first peak, it will however result in the increase of second peak, and vice versa. Thus, the optimal allocation of tests per day is the one that minimizes both peaks of the infected population, which occurs when both peaks are equal. Both BEST and COST policies are compared with the actual COVID-19 testing scenario of France through the prediction of the number of active ICU cases and the cumulative number of deaths.
In Section 2, we describe the model, inputs, and outputs, and compute the basic and effective reproduction numbers of the model. Section 3 illustrates the French COVID-19 data and provides the data imputations to infer the missing data from the raw data. Then, Section 4 provides the estimation and validation of model for the French COVID-19 case. Finally, in Section 5, we propose two testing policies BEST and COST, and evaluate them by comparing the predicted number of active ICU cases and the predicted cumulative number of deaths with the actual data.
We devise a five-compartment model with the purpose of evaluating and devising testing policies. We assume that testing allows for diagnosing and isolating the infected people from the population to prevent the transmission of the disease to the susceptible population. The proposed model is named SIDUR, which corresponds to the five compartments: susceptible (S), undiagnosed infected (I), diagnosed infected (D), unidentified recovered (U), and identified removed (R). The model is characterized by four parameters and one input, which is the testing rate corresponding to RT-PCR tests. • x R (t): Number of removed people who either recover or die after being diagnosed.
The development of the model is based on the following three assumptions: Assumption 1. The population remains constant during the evolution of the epidemics:
where N stands for the total population.
Only the undiagnosed infected population x I (t) is responsible for the disease transmission to the susceptible population x S (t).
Assumption 3. All the deaths from COVID-19 are identified and reported; they are included in the removed population x R (t) along with the people who recover after being diagnosed.
Assumption 4. The efficiency of the acquired immunity is sustainable enough. That is, the unidentified recovered population x U (t) and the removed population x R (t) cannot get infected again.
Based on the above assumptions, the model is given aṡ
where β, γ, and ρ are the infection, recovery, and removal rates, respectively, u(t) is the testing rate,
is the testable population, and θ is the testing specificity parameter which takes values in the interval [0, 1]. The recovery rate γ is the inverse of the average recovery time 1/γ after which an undiagnosed infected person recovers, and the removal rate ρ is the inverse of the average removal time 1/ρ after which a diagnosed infected person recovers or dies. The average recovery time is expected to be shorter than the average removal time, i.e., γ ≥ ρ, because the undiagnosed infected population that comprises the undetected asymptomatic cases and cases with mild symptoms recover faster than the diagnosed population that comprises mostly the cases with severe symptoms.
The infection rate β is the product of the frequency of contacts among the susceptible and infected populations and the probability of disease transmission after a contact has been made. Thus, the parameters β, γ, and ρ are related to the disease biology. However, the value of β can also be partially impacted by a country's government through non-pharmaceutical interventions (NPI) such as social distancing, lockdown, confinement, travel restrictions, and preventive policies (i.e., to maintain a certain distance from other people, to wear a face mask in public spheres, to wash/sanitize hands more often, etc.). The value of β is expected to be smaller when NPI's are implemented than the value of β when no NPI is implemented. Depending on the time periods during which different NPI policies are implemented, we assume the infection rate β to be piecewise constant.
The testing specificity parameter θ, on the other hand, is solely dependent on the testing policy implemented by the public health authority. Given that the testing rate is constant, the value of θ will be larger when the tests are allocated efficiently through contact tracing than the value of θ when the tests are performed randomly. However, there are other factors that can also influence θ, for example, if only the people with severe symptoms are tested, then the probability x I /x T of detecting an infected person from the testable population is equal to one, i.e., the testing specificity parameter θ = 1. This is to indicate that the larger value of θ doesn't necessarily imply the efficiency of testing policy, rather it only signifies the specificity of tests. Depending on the time periods during which different testing policies are implemented, we assume the testing specificity parameter θ to be piecewise constant.
We consider the testing rate u(t) to be the model input which depends on three factors: the daily testing capacity c(t), the remaining stockpile of tests r(t), and the testable population x T (t). We consider c(t) to be time-varying in order to take into account the fact that the capacity of testing can change on a daily basis. Moreover, on a given day t, one cannot do more tests than the daily testing capacity c(t), the remaining stockpile of tests r(t), or the testable population x T (t). Therefore, we have
In case the total stockpile of tests r max is limited, the remaining stockpile of tests at time t is given by
In case new tests can be produced and supplied easily, (3) is simply given by u(t) := min (c(t), x T (t)) .
Usually, the testable population x T (t) is much larger than the daily testing capacity, thus
In order to diagnose the infected people at time t, the tests are allocated to a proportion of the testable population x T (t), which is a sample from the total population N . From (2) , it is obvious that the infected population x I (t) ≤ x T (t) at any given time t. Thus, given the testing specificity parameter θ ∈ [0, 1], the probability of detecting an infected person per test in a homogeneous population structure is given by
The testing specificity parameter θ allows for the adjustment of the testable population to accommodate for the detection rate of tests. In most countries, at the beginning of an epidemic outbreak, the number of available tests are limited. Thus, the available tests are usually utilized to confirm the symptomatic infected cases or to diagnose certain people such as medical care agents, politicians, athletes, etc. In such a case, the testable population is close to the infected population and the value of θ increases to approximately one. Once the capacity of testing is increased, the size of the testable population is also increased that can include, for example, contacts of diagnosed people, the whole population of a city where a cluster is identified, travellers, etc. As a consequence, the value of θ decreases.
The SIDUR model is described by the one-way transfer of population between compartments, where an outflow from one compartment is the inflow to the other compartments. Thus, it suffices to describe only the outflows from the compartments to describe the dynamics of the model.
In the beginning of the epidemic, most of the population is in the susceptible compartment (S) with the exclusion of those who are initially infected and/or diagnosed. Some of the susceptible people in S may get infected and leave this compartment when they come in contact with an infected person. The rate of the outflow from this compartment is according to the infection transmission rate, which depends on the product of the number of susceptible and infected populations, and is given as
where β is the infection rate. The term x I (t)/N is the proportion of undetected infected population at any time t in a homogeneous population structure. Note that in light of Assumption 2, diagnosed population x D (t) does not participate in the infection transmission because they are either quarantined and/or hospitalized, i.e., they are temporarily removed from the population. Finally, by Assumption 4, there is no inflow to the susceptible compartment.
Detection. The outflow from the infected compartment (I) is either due to detection (i.e., transfer to the diagnosed compartment (D)) or recovery without detection (i.e., transfer to the unidentified recovered compartment (U)). The first outflow is due to the testing rate u(t), i.e., the number of tests performed per day. Since the probability of detecting an infected person from a testable population by a single test is x I (t)/x T (t), therefore we have
the rate of diagnosing the infected population in I compartment.
Recovery. The second outflow from the I compartment consists of those people who are not diagnosed and recover naturally with an average recovery period of 1/γ. The unidentified recovered compartment (U) accumulates the infected people who recover naturally without being detected with a recovery rate γ.
Removal. The diagnosed compartment (D) admits u(t)x I (t)/x T (t) as an inflow, whereas the outflow is ρx D (t) with ρ being the removal rate. That is, 1/ρ is the average time period after which a typical diagnosed person either recovers or dies. The removed compartment (R) accumulates the diagnosed people who die or recover with a removal rate ρ.
The outputs signals y i (t), t ∈ R ≥0 , i = 1, 2, . . . , m, from the model (or model outputs) correspond to the sampled output measurements y i (k), k ∈ Z ≥0 , from the data (or data outputs), which can be approximated by a continuous signals y i (t), respectively. The data outputs and model outputs are related as follows:
where w i (t) is the measurement noise/error.
The model outputs are functions of the states of SIDUR model, which are of two types: (i) outputs whose functions are known and (ii) outputs whose functions are not known. First, we define three model outputs whose function is known.
• Cumulative number of diagnosed people
• Cumulative number of removed people
• Number of positively tested people (or positive test results) per day
These model outputs are fitted with the data outputs in order to estimate the model parameters β, θ, γ, ρ in Section 4. Note that these model outputs are related to each other. Since the number of diagnosed infected people at any time t can be obtained as x D (t) = y 1 (t) − y 2 (t), which is also known as the number of active diagnosed cases, we obtain the following relation between y 1 (t) and y 2 (t) from (1e)ẏ
On the other hand, the number of positive test results per day y 3 (t) is related to the cumulative number of diagnosed cases y 1 (t) by the following relation
The cumulative number of diagnosed people y 1 (t) can be obtained by integrating the daily number of positive test results as
These output relations (8), (9), (10) are used to infer the missing data from the available data in Section 3.
Second, on can note that the number of active ICU patients, denoted B(t), (or ICU beds occupied) and the cumulative number of deaths, denoted E(t), (or extinct cases) are positively correlated with the total number of active infected cases A(t) = x I (t) + x D (t) and the cumulative number of infected cases I(t) = N − x S (t), respectively. Therefore, we define two additional model outputs:
• Number of active Intensive Care Unit (ICU) cases (or the number of ICU beds currently occupied by COVID-19 patients):
where
is the number of active infected cases, ψ −1 is the average time period a typical COVID-19 ICU case takes from getting infected to being admitted to ICU, and the function g is to be defined.
• Cumulative number of deaths due to COVID-19 (or extinct cases):
where
is the cumulative number of infected cases, φ −1 is the average time period a typical COVID-19 extinct case takes from getting infected to death, and the function h is to be defined.
The functions g and h will be estimated from available data on the number of active ICU cases and the cumulative number of deaths, respectively, in Section 4. Then, the model outputs B(t) and E(t) will be used as performance outputs to evaluate the testing strategies proposed in Section 5.
An important quantity to assess the epidemic potential of a disease is the basic reproduction number R 0 , which is defined as the expected number of secondary infected cases produced by a single infected person in a completely susceptible population [10] . If R 0 > 1, then each generation of infected cases produces more secondary cases in the next generation and the disease has a potential of becoming an epidemic. If R 0 < 1, then each generation of infected cases produces less secondary cases in the next generation and the disease will eventually die out. It is worth noticing, however, that the definition of R 0 assumes that the people around a primary infected case are all susceptible. This suggests that determining R 0 is important only at the onset of an epidemic. However, in the later stages, more people get infected and not all people around an infected person are necessarily susceptible. As more people get infected, the conditions favoring the disease to propagate change and the number of susceptible people that an infected person infects is actually less than that what R 0 predicts. Thus, a more suitable quantity during the later stages of the epidemic is the effective reproduction number R t , which takes into account the proportion of susceptible people in the total population [47] .
For the SIDUR model, we define the effective reproduction number R t to be the ratio of the inflow and the outflow of the undiagnosed infected compartment (I). That is, R t is the ratio of the number of newly infected people and the number of newly diagnosed and recovered people at time t. If R t < 1, this means that more people are being diagnosed and recovered than the people being infected at time t, which implies that x I (t) will decrease. If R t > 1, this means that more people are being infected than the people being diagnosed and recovered at time t, which implies that x I (t) will increase. Notice that the definition of the basic reproduction number R 0 is same as R t for t = 0.
To derive the expression of R t , we consider the model equation (1b), where the undiagnosed infected population x I (t) satisfieṡ
The positive rate (inflow) βx S (t)x I (t)/N tells how many new infections will be generated in the next moment, and the negative rates (outflows) u(t)x I (t)/x T (t) and γx I (t) tell how many infected people will be diagnosed or recovered in the next moment, respectively. Therefore, the effective reproduction number is the following ratio
To derive the expression of R 0 , we consider the expression of R t at t = 0, which corresponds to the onset of the epidemic. For t = 0, we can assume few infected cases, which implies that x S (0) ≈ N and x T (0) ≈ (1−θ)N . Under these approximations, we have
This expression can also be obtained by following the methodology of [48] . Notice that the reproduction number R 0 depends on the initial testing policy u(0). This indicates that it is possible to suppress the epidemic in the beginning by having an intensive testing policy, which can be seen, for example, in the case of South Korea [49] . In general, however, we have u(0) ≈ 0, which gives R 0 ≈ β/γ.
The data related to COVID-19 in France is collected from the French government's platform for publicly available data 3 for the time period of January 24 to July 01, 2020. In particular, we use datasets provided by the French Ministry of Social Affairs and Health (Ministère des Solidarités et de la Santé (MSS)) and the French Public Health Agency (Santé Publique France (SPF)). From MSS, we obtain the data about different categories of people affected by COVID-19, i.e., diagnosed, hospitalized, recovered from hospitals, and dead. From SPF, we obtain the data for the number of PCR tests performed and positive test results obtained per day.
The data obtained from both sources is incomplete in several aspects. For instance, the data for the number of recovered people does not record those who recover from their homes after being diagnosed. These people do not show severe symptoms of the disease and, therefore, are not hospitalized, but are quarantined in their homes for some days. Only those who are hospitalized after being diagnosed are recorded as recovered when they are discharged from the hospital. On the other hand, the data for COVID-19 PCR tests is also incomplete. To illustrate this, we consider three intervals of time: (1) January 24 to March 09, 2020, (2) March 10 to May 12, 2020, and (3) May 13 to July 01, 2020. There is no data available for the tests during the first interval. During the second interval, the testing data is collected only from the medical laboratories and not from the hospitals. However, we have reliable data only during the third interval which is collected both from the medical laboratories and the hospitals. Therefore, the data obtained from the above sources can be considered as a raw data which needs to be imputed.
This subsection illustrates the data obtained from MSS and SPF without any modification.
Cumulative number of diagnosed cases. We denote the data for the cumulative number of diagnosed cases by y 1 , which is illustrated in Figure 2 and corresponds to the model output y 1 (t) in (5) . It is also known as the total "confirmed" cases. This is a cumulative data for all the cases diagnosed with the disease through RT-PCR tests 4 . Thus, it includes both the active cases (those who are either admitted to the hospitals and/or quarantined) and the inactive cases (those who either recovered or died after being diagnosed). That is, y 1 (k) corresponds to the sum of people in the diagnosed (D) and removed (R) compartments of the SIDUR model ( Figure 1 ) on a given day k, as given in (5).
There is also an additional data for the diagnosed cases from French retirement homes (EHPAD). However, the French government database 5 and several other international databases 6 7 do not add the diagnosed cases from EHPAD to the cumulative number of diagnosed (confirmed) cases. That is, the data for the cumulative number of diagnosed cases is considered to be inclusive of the diagnosed cases from EHPAD. However, in all the above databases, the data on cumulative number of deaths is collected separately from both the hospitals and EHPAD. Number of active hospitalized and ICU cases. The data on the number of active hospitalized cases is denoted as H and is illustrated in Figure 3 along with the number of active ICU cases B. This data corresponds to the number of people who are admitted to the hospitals and/or ICU on a given day. That is, it is not a cumulative data. Moreover, it doesn't include those who were diagnosed but not hospitalized. That is, this data corresponds to a certain proportion of people in the diagnosed compartment (D) of the SIDUR model. This data is available from March 17, 2020, onward. Cumulative number of recovered cases from hospitals. This data is illustrated in Figure 3 . It corresponds to people who, after recovering from the disease, were discharged from the hospitals. Obviously, prior to recovering, they were diagnosed with the disease and hospitalized due to their severe symptoms.
Cumulative number of deaths. The data on the cumulative number of diagnosed cases is considered to be inclusive of the diagnosed cases from the French retirement homes (EHPAD). However, the case for data on the cumulative number of deaths is different. Those who died at the hospitals and those who died in the retirement homes (EHPAD) are considered to be distinct. Thus, the cumulative number of deaths is the sum of both data, which are illustrated in Figure 4 .
Number of tests and positive tests per day. We have two types of data related to COVID-19 PCR tests. The first type of data is collected by SPF on the number of tests performed and positive test results per day from March 10 to May 26, 2020. However, this data is collected only from the central sampling laboratories: Eurofins Biomnis and Cerba. Figure 5 (a) (blue) illustrates the number of tests performed per day and Figure 5 (b) (blue) illustrates the number of positive test results per day. The second type of data was made available after the deployment of a new information screening system (SI-DEP) by the SPF. This data is available from May 13, 2020, onward. It is collected from both the laboratories and the hospitals. However, the data reported by SI-DEP is the number of 'tested people' per day instead of the number of 'tests performed' per day. SI-DEP guarantees that only one test is counted per person. In the case of, for instance, multiple negative test results for a certain person, SI-DEP considers only the first date on which the PCR test was performed. Later, if that person gets a positive test result, then only this new result is reported in the data and the previous data is erased. Figure 5 
In the raw data, we only have the data for those who recover or die in the hospitals after being diagnosed with the disease. However, the removed compartment of the SIDUR model also comprises the diagnosed cases who were not hospitalized but were quarantined in their homes. There is no data that records the recovery of these people. Moreover, the data on PCR tests is also incomplete; there is no data on PCR tests from January 24 to March 09, 2020, and the data from March 10 to May 12, 2020, doesn't include the tests performed in the hospitals. Therefore, in order to infer the missing data, we impute the raw data by making reasonable assumptions.
From the data on the total number of recovered people from hospitals shown in Figure 3 , we see that 76,540 people have recovered from the hospitals as of July 01, 2020. If we subtract this number and the total number of deaths (Figure 4) , i.e., 29,860, from the total number of diagnosed cases (Figure 2) , i.e., 165,700, we obtain 165, 700 − 76, 540 − 29, 860 = 59, 300 people. Further subtracting the currently hospitalized cases (Figure 3) , i.e., 8336 as of July 01, we obtain 59, 300 − 8336 = 50, 964 people, who might still be infected or have recovered. These people were diagnosed but were not hospitalized; they were quarantined in their homes. However, there is no data that provides a correct answer for how many people have recovered and how many of them are still infected. Therefore, using the relevant raw data, we infer the cumulative number of removed cases y 2 by estimating the number of diagnosed cases who recovered from home.
We use the following notations for simplicity and brevity: y 1 Total diagnosed and hospitalized y 1 Total diagnosed but not hospitalized y 2 Total recovered/died in a hospital y 2 Total recovered from home after diagnosis By definition, we have y 1 (k) = y 1 (k) + y 1 (k) y 2 (k) = y 2 (k) + y 2 (k) (15) where k is from January 24 to July 01, 2020. Note that y 1 (k) is illustrated in Figure 2 and
where y 2 (k) is the sum of the total number of recovered cases from hospitals ( Figure 3 ) and the total number of deaths (Figure 4) , and H is the number of active hospitalized cases ( Figure 3 ). Thus, we can compute the total diagnosed cases who were not hospitalized as
Since there is no data for the diagnosed people who recovered from home, therefore y 2 (k) is unknown. Thus, we assume the following:
That is, the ratio of the diagnosed cases who recovered in homes to the total diagnosed cases who were quarantined at homes is equal to the ratio of the diagnosed cases who recovered or died in hospitals to the total diagnosed cases who were hospitalized. In other words, we assume that the removal rate of people who were not hospitalized is equal to the removal rate of people who were hospitalized. Thus, from (15) and (16), we obtain
which corresponds to the model output y 2 (t) in (6).
From Figure 5 , we see that the first type of data, which is available from March 10 to May 26, 2020, considers the number of tests performed and positive test results per day. On the other hand, the second type of data, which is available from May 13, 2020, onward, considers the number of tested people and positively tested people per day. However, no person is usually tested more than once per day. Therefore, we assume that the number of tested people per day is same as the number of tests performed per day. Similarly, the number of positively tested people per day is same as the number of positive test results per day. Note that if a person is tested more than once but on different days, then this assumption is not violated.
We consider three time intervals: (i) January 24-March 09, when there is no data on PCR tests; (ii) March 10-May 12, when there is incomplete data; (iii) May 13-July 01, when there is complete data. Let u and y 3 denote the number of tests performed and the number of positive test results per day, respectively, for the entire interval January 24 to July 01, 2020. Let u , u , u and y 3 , y 3 , y 3 be the number of tests performed and the number of positive test results obtained for the first, second, and third time intervals, respectively. Since the data in the third interval is reliable, we do not make any imputations for u and y 3 . For the other two intervals, we make reasonable assumptions to complete the data.
(i) January 24-March 09: This interval corresponds to the beginning of the epidemic in France and the data for tests performed and positive test results per day for this interval is u and y 3 , respectively. During this interval, only those people were tested who showed symptoms. Moreover, recall the output relation (9) . Then, we compute the data as follows: u (k) ≈ y 3 (k) and y 3 (k) = y 1 (k + 1) − y 1 (k), where k is from January 24 to March 09. In other words, during the first interval, the number of tests performed per day is assumed to be approximately equal to the number of positive test results obtained per day. Moreover, the number of positive test results obtained per day is equal to the number of diagnosed cases that day.
(ii) March 10-May 12: In the second interval, we have the data on PCR tests that is reported only by the laboratories and not by the hospitals. During this interval, we compute the data as follows: u is same as the data ( Figure 5 ) and y 3 (k) = y 1 (k + 1) − y 1 (k), where k is from March 10 to May 12.
Based on the above data imputations, we obtain the number of tests performed per day u = u u u , which corresponds to the control input u(t), and the number of positive tests obtained per day y 3 = y 3 y 3 y 3 , which corresponds to the model output y 3 (t) in (7), for the complete time interval January 24 to July 01, 2020. 
The input u(t) to the model corresponds to the number of tests performed per day u. Let {1, 2, . . . , τ } be the index set of the 160 days from January 24 to July 01, 2020, where τ = 160. Then, for k = 1, 2, . . . , τ , we define the input signal as u(t) = u(k), for t ≤ k < t which is illustrated in Figure 6 . We denote by y i (t), for i = 1, 2, 3, the outputs obtained from the data. That is, the output signals from the data and the model are related as follows
where w i (k), for i = 1, 2, 3, represents the measurement noise. The outputs y 1 (t), y 2 (t), and y 3 (t) correspond to the cumulative number of diagnosed cases, the cumulative number of removed cases, and the number of positive test results per day, respectively. Similar to the input u(t), we define the output signals from the data as    y 1 (t) = y 1 (k), for t ≤ k < t y 2 (t) = y 2 (k), for t ≤ k < t y 3 (t) = y 3 (k), for t ≤ k < t which are illustrated in Figure 7 .
In this section, we validate the SIDUR model by estimating the model parameters ρ, β, θ, and γ for the case of COVID-19 in France.
The removal rate ρ can be directly estimated from the data outputs y 1 and y 2 . Consider a daily sampling of the model equation (1e), which leads to
where ∆ stands for the forward difference operator, i.e., ∆x R (k) = x R (k + 1) − x R (k) for some nonnegative integer k. Therefore, from the relation between y 1 and y 2 in (8), we obtain ∆y 2 (k) = ρ y 12 (k) + e(k) (17) where y 12 (k) = y 1 (k) − y 2 (k) and e(k) is the error term due to measurement noise. Then, the problem of estimating ρ can be formulated as follows: Find ρ * such that
Notice that the solution of this problem can be obtained through least-square estimation [50, Chapter 7].
We formulate a problem of fitting the model outputs y 1 (t), y 2 (t), y 3 (t) to the data outputs y 1 (k), y 2 (k), y 3 (k), where k = 1, 2, . . . , τ with τ being the final time. The model fitting is done by optimizing the parameters β, θ, γ for the time interval [0, τ ] under the assumption that γ is constant whereas β and θ are piecewise constants.
To limit the rate of spread of COVID-19, the French government announced to place a lockdown all over France from March 17 to May 10, 2020, which included restricted human mobility, strict social distancing measures, and closure of schools, offices, and marketplaces. However, the essential services and public establishments were authorized to remain open under strict preventive measures. People were allowed to leave their homes with face masks only for necessary groceries, brief exercise within a certain radius of their homes, or for urgent medical reasons. Such an intervention from the public authority is necessary to mitigate the rate of spread of the disease and to reduce the value of infection rate β. Therefore, in relation to the case of France, we divide the time into three intervals: (i) Before lockdown (January 24 to March 16), (ii) During lockdown (March 17 to May 10), and (iii) After lockdown (May 11 to July 01). We consider a different value of the infection rate β during each of these intervals, i.e., Let p = [ β 1 β 2 β 3 θ 1 θ 2 γ ] T be the parameter vector. Then, the goal is to find p * such that
where the cost function is given by Note that we consider the data from January 24 to July 01, 2020, therefore we have τ = 160 days.
To solve this problem, one can also pose it as a least-square estimation, as we did for the removal rate ρ, by defining relations between the data outputs y 1 , y 2 , y 3 . However, such relations include the difference operator ∆ applied twice to the data outputs, which is usually not recommended when the data is noisy because it amplifies the measurement noise. Moreover, the gradient-based estimation algorithms [51, Chapter 4] are also not suitable due to the difficulty of computing the gradient of the cost function J online with respect to the parameter vector p. This is because the model outputs y 1 , y 2 , y 3 do not depend directly on the parameters but through the solution trajectories of the SIDUR model. For simplicity, therefore, we choose the particle swarm optimization (PSO) [52] described in Appendix A, which is a 'derivative-free' algorithm, to estimate the parameter vector p.
Infection rate β 1 = 0.3708 β 2 = 0.0707 β 3 = 0.3717 Testing specificity θ 1 = 0.9948 θ 2 = 0.9967 Recovery rate γ = 0.1589 Removal rate ρ = 0.0499 The estimated parameter values are provided in Table 1 . The estimated recovery rate γ and removal rate ρ show that an undiagnosed person recovers in an average period of about 6.3 days and a diagnosed person recovers or dies in an average period of about 20 days. The testing specificity parameter changes slightly from θ 1 = 0.9948 to θ 2 = 0.9967, which can have significant impact on the positive test results because it multiplies with the sum of the susceptible and unidentified recovered population in (2) that is in the order of 10 7 in the case of France.
The infection rate β changes its value twice. First, it drops from β 1 = 0.3708 to β 2 = 0.0707 when the lockdown is implemented in France on March 17, which significantly decreased the rate of the epidemic spread. Then, it rises from β 2 = 0.0707 to β 3 = 0.3717 when the lockdown is lifted on May 10. Many restrictions like social distancing and wearing of face masks were still in place after May 10 in order to prevent the spread of COVID-19 in France. However, the increase in the value of β can be explained by the summer vacations when people were allowed to travel everywhere across France and Europe 9 . This made the places with tourist attractions very crowded and resulted in a higher infection rate.
Using the estimated values of the model parameters in Table 1 , we run the model from January 24 to July 01, 2020. The model fits the output signals data as shown in Figure 8 .
The basic reproduction number R 0 at the outbreak of the COVID-19 epidemic in France is computed using (14) and the average value of effective reproduction number R t during three phases (before, during, and after lockdown) are computed using (13) . For the 'after lockdown' phase, we chose July 01, 2020, to compute R t because it is the date up to which our data is considered. These computed values are shown in Table 2 with a comparison to the ones reported by the government 10 11 . The change in the value of R t also influences the evolution of diagnosed population x D (t). This is because larger value of R t results in a larger infected population x I (t) and smaller value of R t results in a smaller infected population x I (t), which respectively increases and decreases the probability of detecting an infected person x I (t)/x T (t) by a single test. Keeping the number of tests performed per day same, the larger probability of detection x I (t)/x T (t) results in a larger diagnosed population x D (t).
In Table 2 , we see that the placement and lifting of lockdown on March 17 and May 11, respectively, had a significant impact on the value of R t . Such an effect on R t impacted the evolution of the diagnosed population x D (t) = y 1 (t) − y 2 (t), which can be interpreted as the number of active confirmed cases and is illustrated in Figure 9 . The placement of lockdown reduced the value of R t and resulted in less number of active confirmed cases as compared to the scenario in Figure 9 where the lockdown was not placed on March 17. In this scenario, as shown in Figure 9 , the number of active confirmed cases would have increased to a point that could have challenged the available medical facilities such as hospital beds, ventilators, and ICUs. On the other hand, the lifting of lockdown increased the value of R t and resulted in more number of active confirmed cases as compared to the scenario where the lockdown was not lifted on May 11. 
The number of active ICU patients B(t) is a function of the number of active infected people A(t). Since an infected person starts to show symptoms after the average incubation period of approximately 5 days, [53] , and a person takes on average 12 days from being diagnosed to being addmitted to ICU, [54] , we assume ψ −1 = 5 + 12 = 17 days to be the average time delay from getting infected to being admitted to ICU for a typical COVID-19, critically ill case. Thus, we model the number of active ICU patients B(t) as a function of A(t − ψ −1 ), which is approximated by:
where b 1 and b 2 are the parameters (Appendix B, Table B .4) that are determined via the least-square solution to fit (20) to the data on the number of ICU patients. This is illustrated in Figure 10 . Similar to the case of the number of active ICU patients, a typical nonsurviving case has an average incubation period of 5 days and, in addition to that, an average removal period of ρ −1 ≈ 20 days, where ρ is the removal rate, Table 1 . Thus, assuming φ −1 = 5 + ρ −1 = 25 days to be the average time delay from getting infected to death of a typical non-surviving COVID-19 case, we model the number of deaths E(t) as a function of I(t − φ −1 ), which is approximated by the following polynomial:
where e i , for i = 1, . . . , 10, are the parameters (Appendix B, Table B .4) that are determined via the least-square solution to fit (21) to the data on the number of deaths. This is illustrated in Figure 12 . Using the relations (20) and (21), we illustrate the model fit of the number of active ICU cases and the cumulative number of deaths with the data in Figure 11 and 13.
In this section, we use the model validated with the COVID-19 data of France to design two types of testing policies. The first one assumes that the total stockpile of tests is unlimited, but the testing capacity per day is limited. This results in the so called Best Effort Strategy for Testing (BEST), which gives the minimum number of tests needed to be performed per day in order to stop the epidemic from growing. In other words, if the BEST is applied, then the number of new infections stop to grow with respect to time. The second strategy, named COST (Constant Optimal Strategy for Testing) assumes that the total stockpile of tests in a country is limited. In such a case, we investigate the optimal testing policy per day that results in a minimum epidemic peak. In contrast with BEST, there is an optimal value of tests to be performed per day which is smaller than the maximum testing capacity per day.
Assume that the total stockpile of tests is unlimited during the whole epidemic period. This assumption is valid for a country that can manufacture or buy tests continuously during the time of the epidemic. Based on this, we provide a testing policy recommendation on the daily testing capacity starting from a certain time t * in order to change the course of the epidemic in a sense that is defined below. For simplicity, we further assume that the number of tests performed per day is considered to be the daily testing capacity. In other words, the daily testing capacity is utilized completely each day, i.e., u(t) = c(t).
We say that, at any time t, an epidemic is spreading if the number of undiagnosed infected population x I (t) is increasing, i.e., the effective reproduction number R t > 1. On the other hand, an epidemic is non-spreading if x I (t) is not increasing, i.e., the effective reproduction number R t ≤ 1.
Definition 1 (BEST). The best effort strategy for testing (BEST) at a given time t * is the minimum number of tests to be performed per day from time t * onward such that the epidemic switches from spreading to non-spreading at t * .
In other words, BEST provides the smallest lower bound on the number of tests performed per day, sufficient to change at given time t * the course of the epidemic from spreading to non-spreading. In order to compute the BEST at a given time t * , we first define the following function:
where, by definition, for any scalar z, |z| + = z if z > 0 and |z| + = 0 otherwise.
Assume that the infection rate β is non-increasing while the testing specificity parameter θ is non-decreasing on a time interval [t * , t 1 ), for some t * < t 1 . Then, the best effort strategy for testing (BEST) at time t * is given by
Proof. In order to prove that c * (t * ) for t ∈ [t * , t 1 ) is the BEST at t * , we show the following:
(i) If u(t) > c * (t) (resp., u(t) ≥ c * (t)) for any t ∈ [t * , t 1 ), then x I is decreasing (resp., non-increasing) on [t * , t 1 ).
(ii) If u(t) > c * (t * ) (resp., u(t) ≥ c * (t * )) for any t ∈ [t * , t 1 ), then x I is decreasing (resp., non-increasing) on [t * , t 1 ).
Assume that u(t) > c * (t) on [t * , t 1 ). Then, Φ(t) := β(t) xS(t) N − u(t) xT(t) − γ < 0 which implies that x I is decreasing sinceẋ I (t) = Φ(t)x I (t) almost everywhere. If only the weaker assumption u(t) ≥ c * (t) on [t * , t 1 ) is fulfilled, then, by using the continuity of the solutions of ODE with respect to perturbations of the right-hand side, one gets that x I is non-increasing.
Assume now that u(t) > c * (t * ) on [t * , t 1 ), where c * (t * ) is constant. Then, by continuity, u(t) > c * (t) on a certain interval [t * , t 2 ), for some t 2 ∈ (t * , t 1 ]. As a consequence of the result (i) shown previously, x I decreases on [t * , t 2 ). Moreover, assume that t 2 is the maximal point in (t * , t 1 ] having this property. In order to show that t 2 = t 1 , it is sufficient to show that u(t 2 ) > c * (t 2 ), otherwise one may consider a larger value for t 2 which will lead to a contradiction with the fact that it is maximal. Since x I decreases on [t * , t 2 ) and θ is non decreasing, from Lemma 6 (see Appendix C), we can conclude that x T also decreases on this interval. On the other hand, since x S is always decreasing and β is non increasing, one can conclude that c * (t) also decreases on [t * , t 2 ). This is obtained by upper bounding c * (t). Thus, one has c * (t * ) > c * (t), which implies that that u(t 2 ) > c * (t 2 ). Therefore, as t 2 = t 1 , we have established that x I decreases on the whole interval [t * , t 1 ). For the case where u(t) ≥ c * (t * ), we can use the same argument of continuity of the trajectories.
From the previous results, one deduces that the BEST is given by c * (t * ) and the testing rate u(t) ≥ c * (t * ) for t ∈ [t * , t 1 ). If u(t) < c * (t * ), for t ∈ [t * , t 1 ), then one can show easily that the epidemic goes on spreading in the interval [t * , t 1 ). Hence, u(t) = c * (t * ) is the BEST policy at t * .
Algorithm 1: Computation of the BEST policy at time t * . 1. Inputs: N , β, γ, θ, t * , x S (t * ) and x T (t * ).
2. Compute the BEST policy c * (t * ) using (22) . 4. Return to step 2 if β increases or θ decreases. Proposition 1 states that the peak of x I (t) is uniquely determined by the BEST policy c * (t * ), where the peak is achieved at time t * . Therefore, Algorithm 1 can be used to set the peak time t * once parameters β, γ, and θ are learned from the data. Remark 1. Requiring that β must not increase and θ must not decrease in the interval (t * , t 1 ) for some t 1 > t * is necessary for the BEST policy. It is thus important to keep the external conditions that determine the values of β and θ either constant or such that β decreases (e.g., through the implementation of lockdown) and/or θ increases (e.g., through efficient contact tracing).
Remark 2. The case where β decreases and/or θ increases at some time t 1 > t * has the effect of speeding up the suppression of the epidemic under BEST policy.
Remark 3. From (1b), we can note that if x S (t)/N < γ/β, the the epidemic naturally decreases. In this case, doing no testing u(t) = 0 is the BEST policy, which, by definition, gives a minimum number of tests to be performed in order to stop the growth of the infected population x I . However, if testing is resumed in this case, i.e., u(t) > 0, it will further speed up the decrease of the infected population.
Giving data from France, we first compute c * (t * ) for different values of t * from January 24 to March 13. Figure 14 shows the number of tests per day required by the BEST policy if it is implemented on day k and the corresponding value of peak of infected cases x I (k * ). One can note that the later BEST is applied the higher is the required number of tests. An exponential increase can even be observed from February 28 which corresponds to an acceleration of the infection. Now, we consider a scenario where BEST is implemented on March 01. Figure 15 depict the number of active cases when u(t) is the actual testing scenario (see Figure 6 ) and when u(t) is given by BEST. To evaluate BEST, we use u(t) as given by recorded data before March 01 then use u(t) given by BEST from March 01. In the first case, the peak of the infected population x I (t), which are the active undiagnosed cases, is about 6 million. In the second case, the peak of infected population in this case is 363, 169. The required number of tests per day to be performed for the implementation of BEST on March 01 is c * ≈ 147, 000. The impact in terms of ICU occupation and number of deaths is now evaluated using the equations (20) and (21) respectively. The results are illustrated in Figure 16 and 17.
We observe that the peak of the number of active ICU patients could have been reduced by 34.71% and the number of deaths could have been reduced by 74.45% if the BEST policy was applied from March 01, 2020.
Consider now that the total stockpile of tests is limited and given by r max . In this case, the testing rate should be chosen carefully so as to not finish the stockpile of tests too early by performing too many tests per day, or to be unable to control the epidemic spread by performing too few tests per day. In the former case, when the stockpile of tests finishes too early due to intensive testing, the infections, even if significantly reduced in the beginning, will start to spread again and result in a second wave of the epidemic with an infection peak much higher than before. On the other hand, if the number of tests performed per day is too little, this would result in the infected population reaching a very high peak during the first wave of the epidemic, which could challenge the available medical facilities of a country. Thus, in this section, we determine a constant optimal allocation of a limited stockpile of tests r max such that both peaks of the infected population are minimized.
Given the total stockpile of tests r max , we assume that the number of tests performed per day is given by where the time period T := r max /C represent the duration of the testing policy once C is determined.
Definition 2 (COST). The constant optimal strategy for testing (COST) with the total stockpile of tests r max is the policy of class (23) that minimizes the peak of the infected population x I (t).
In other words, the COST policy allocates the limited stockpile of tests r max as C tests per day for the time interval [0, T ], where T is measured in days, such that the maximum value of the infected population x I (t) is minimized.
In the following, we first study the SIDUR model in a new coordinate of 'infection' time ξ with the aim of finding analytic solutions. Then, we use those solutions to compute the two peak values of x I (t), where the first peak arrives at t ≤ T and the second at t > T . Finally, we compute the optimal value of C that minimizes those peak values by equating them, which is depicted in Figure 18 .
The SIDUR model has an interesting property: most of the terms on the right-hand side of the model (1) depend linearly on x I (t). This gives an intuition that the epidemic process goes faster when the number of infected people x I is high. Indeed, this idea was used very early to analyze the basic SIR model, see [9] . It is possible to define a new time variable ξ, which we call infection time, as an integral of undetected infected people up to time t, i.e., Figure 18 : Under the limited stockpile of tests rmax, the number of infected cases xI(t) will have two peaks. The first peak happens before the stockpile of tests is finished and the second peak happens after the stockpile finishes. The second peak is the result of the second wave of the epidemic during which no testing is done. If the testing policy is chosen such that the first (respectively, second) peak is reduced, then the second (respectively, first) peak will increase. The minimum of both peaks is achieved when the peak values of both waves of the epidemic are equal. This is ensured by COST policy C.
time by:
Note that the infection timescale preserves the peak value of epidemic, which is of interest to us. That is,
and dx I /dt = 0 implies dx I /dξ = 0 and vice versa (as long as x I = 0). Therefore, one can perform a model analysis based on the infection timescale ξ. Moreover, by writing the equations in ξ, the partial solution of SIDUR model can be obtained analytically.
Proposition 2. Evolution of the unidentified recovered population x U is affine with respect to the infection time ξ, whereas the susceptible population x S is exponential with respect to ξ. That is,
Proof. It is straightforward to see that the solution of (1d) is given by (25) . On the other hand, (1a) yields
whose solution is given by (26) .
It is not so trivial to find analytic solutions for x I , x D , and x R . However, under certain hypothesis on the testable population, it is possible to get a reasonable approximation of the solutions.
Assumption 5. The testable population x T is approximated as
This makes it possible to find an analytic solution for the infected population x I (t).
Proposition 3. Evolution of the infected population x I under Assumption 5 is given by
Proof. Equation (28) can be directly obtained by integrating (1b), if one first transforms it using the infection time and then substitutes the definition for x T from Assumption 5 and the solution for x S from (26).
Using (23), we can now write a complete evolution of the infected population:
where ξ * denotes the moment of stopping the testing in the infection time coordinates, which by definition of the infection time is given by an implicit formula
derived from the relation
where we substitute the solution for x I (ξ) given by (28) , and the fact that u(t) = C during the testing interval [0, T ].
From the definition of the effective reproductive number (13) , for convenience, we can introduce two reproduction numbers at the origin: R C with the COST policy C and R W without testing:
In general the evolution of the infected population given by (29) can have at most two peaks, one during testing and one after stopping. To have guarantees that the peaks exist, we further state the following assumption:
The COST policy C does not suppress the epidemic at the origin, i.e., R C > 1.
To explain this assumption, suppose in contrary that the COST policy C suppresses the epidemic at the origin. However, such a policy will not develop the herd immunity, and, in the case of the limited stockpile of tests r max , the second peak will arrive once the testing is stopped. Therefore, it is reasonable to develop at least a partial herd immunity using the controlled number of tests C in the first wave of epidemic.
When the two peaks of the infected population exist, their values can be obtained by setting the derivatives of (29) with respect to ξ to zero:
where the peaks positions in the infection time ξ are given as:
Writing them in terms of reproduction numbers R C and R W , we obtain
Since R W > R C > 1, by Assumption 6, we get ξ peak2 > ξ peak1 > 0. Further one should notice that the peak ξ peak1 occurs only if ξ peak1 < ξ * , otherwise it is ill-defined. Similarly, ξ peak2 > ξ * .
The peaks values themselves are given by
It is obvious that x I,peak1 decreases as C increases, because x I,peak1 represents a peak under still performed testing.
Further, it is possible to prove that Cξ * decreases as C increases for all sufficiently large C. Indeed, (30) can be rewritten as
Taking the derivative with respect to ξ * and substituting again (30) in place of the integral, we see that
which is the same as
Therefore,
where ξ * is rewritten by the definition of the infection time (24) , while r max /C is represented by the integral of 1 over the time T = r max /C. The result is negative unless x I (T ) is smaller than an average number of infected people all the time before, which means that T is long after the peak has come. Thus for all reasonable sufficiently large C such that T is not too large C, x I,peak1 decreases, while x I,peak2 increases. This property can be used to optimize a maximum between two peaks. Indeed, a maximum between a decreasing and an increasing function is minimized when they are equal each other. Thus, we equate two peak values in (35) and solve it with respect to ξ * , obtaining an optimality condition as
Equation (36) provides an optimality condition which, being combined with the connection between ξ * and C given by (30) , constitutes a complete system determining the COST C.
It is possible to show that the optimal strategy provides switching in between of two peaks ξ peak1 and ξ peak2 , thus both peaks exist. Indeed, first we prove a simple lemma:
and the equality is possible only if x = 1.
Proof. First of all, it is obvious that for x = 1 all three parts of the inequality are equal to zero. Further, let us take the derivative of this inequality:
which is always true for all x ≥ 1, and it holds strictly for all x > 1. This concludes the proof. Now, consider the difference ξ * − ξ peak1 :
where we denote x = R W /R C , and the last inequality comes from Lemma 4. Thus the first peak is well-defined. Further, writing the difference ξ peak2 − ξ * in the same manner and using the same definition for x, we see that
Thus we have proven that 0 < ξ peak1 ≤ ξ * ≤ ξ peak2 , which means the stopping of testing happens in between of the two peaks.
Finally, the COST policy can be obtained by solving (30) with respect to C, where (36) are used as an upper limit of the integral. This can be formally stated as follows:
Proposition 5. Consider the testing rate defined by (23) . Then, under Assumptions 5 and 6, the COST policy C is given by solving (30)-(36), i.e.,
The time duration for positive testing rate (23) is given by T = r max /C. Note that, by construction, a solution to the system of equations (37) and (38) exists and is unique as long as Assumption 6 holds. This system of equations can be solved numerically for C by using the Newton's method (if the initial estimate C 0 is chosen sufficiently close to the real optimal value C). The process is sketched in Algorithm 2.
Algorithm 2: Computation of the COST policy.
1. Take the initial estimate for COST C 0 .
2. On n-th iteration step, given C n , find ξ * n by (38) .
3. Find the value of the right-hand side of (37). This can be done by computing the integral numerically. Denote the result as f (C n ).
Find the derivative f (C n ) of the right-hand side of (37). This can be done analytically: (29) and, from (38),
5. Once f (C n ) and f (C n ) are computed, update the COST by the Newton's method:
6. Repeat, from step 2, until the desired accuracy is achieved.
Infection rate β 1 = 0.2643 β 2 = 0.0006 β 3 = 0.0642 Testing specificity θ 1 = 0.9415 θ 2 = 0.7993 Recovery rate γ = 0.0542 Removal rate ρ = 0.0499 
We first validate Assumption 5. This assumption leads to a new set of model parameters given by Table 3 . With these parameters, the model validation is illustrated in Figure 19 . Table 3 .
Similar to BEST, we evaluate the COST through the prediction of the number of active ICU cases B(t) and the cumulative number of deaths E(t) using equations (20) and (21), respectively. However, we fit these equations to the data B(k) and E(k) by using the model parameters of Table 3 that are estimated under the Assumption 5. Assuming the total stockpile of tests r max = 2, 038, 037, which is the total number of tests performed from January 24 to July 01, 2020, in France, we obtain the optimal number of tests to be performed per day as C = 17, 144 from Algorithm 2, where the terminal time T = 118 (i.e., May 20) . This is illustrated in Figure 20 along with the actual number of tests per day from the data.
Using the actual number of tests per day, we fit the model outputs (20) and Table B .5. Then, we predict the number of active ICU cases B(t) and the cumulative number of deaths when the COST is applied instead of actual number of tests, as illustrated in Figure 21 and Figure 22 , respectively. We observe that COST reduces the peak of B(t) by 11.12% and the total number of cumulative deaths E(t) by 37.52%. Notice that the total stockpile of tests in the actual testing and COST is the same. However, most of the tests in the actual testing are consumed after May 13, whereas all the tests in COST are allocated equally in the time interval January 24 to May 20. The number of tests performed per day in the actual testing exceeds 50, 000 per day after May 13, which is about three times more than what is required by COST, i.e., 17, 144. Thus, COST is optimal in a sense that it is practical, it decreases the burden on medical facilities, and it reduces the number of deaths significantly.
We proposed a SIDUR model for the control of the COVID-19 epidemic through the number of RT-PCT tests performed per day. Such tests enable the government to diagnose and isolate the infected people from the susceptible population. We estimated and validated the model on the French COVID-19 data, and proposed two testing policies to control the epidemic: 1) best effort strategy for testing (BEST) and 2) constant optimal strategy for testing (COST). BEST provides the minimum number of tests to be performed from a certain day onward in order to make the increasing infected population non-increasing immediately. That is, it changes the course of epidemic from spreading to nonspreading. On the other hand, COST considers a limited stockpile of tests that are optimally allocated in a time interval starting from the beginning of the epidemic in order to minimize the peak of infected population.
The control input in SIDUR model corresponds to the number of RT-PCR tests performed per day. However, another type of test, a serology test, which is not considered in the current model because of the unavailability of its data, is also very important. A serology test determines the relevant antibodies in a subject's serum in order to detect whether he/she was infected in the past. By performing serology tests on the testable population, one can detect the unidentified recovered population and transfer them in the identified removed compartment of the model. This reduces the size of the testable population, which in turn increases the testing specificity of RT-PCR tests. In other words, the serology tests complement the RT-PCR tests [34, 35] . Therefore, as a future prospect, it will be interesting to consider two control inputs corresponding to both types of test in the SIDUR model.
The model is estimated and validated by fitting the model outputs with the available data of France. This allows us to predict the unmeasured states of the model. However, there is no certainty whether the predicted states correspond to the reality. Another prospect is to design an observer for the SIDUR model in order to estimate the true states of the model.
Both BEST and COST policies are easy to compute and implement, however they are static. Thus, their influence on the control of epidemic is limited. In future, it will be interesting to solve a finite/infinite-horizon optimal control problem to minimize the peak and/or cumulative number of the infected population by a dynamic control input. and by Inria, France, in the framework of the COVID-19 fast-track research program.
The raw data used in this paper is available at https://www.data.gouv.fr/ fr/datasets/ with the following pointers: Chiffres-clès concernant l'épidémie de COVID-19 en France, Données relatives aux tests de dépistage de COVID-19 réalisés en laboratoire de ville, and Données relatives aux résultats des tests virologiques COVID-19 (SI-DEP). The imputed data and codes are available at http://scale-freeback.eu/.
In this appendix, we briefly describe the particle swarm optimization (PSO) algorithm, [52] , which is considered to be one of the most powerful algorithms. It considers a foraging swarm of n particles who collectively search for an optimal solution of (18) in the parameter space. At time step h = 0, 1, 2, . . . , each particle i visits a positionp i h by moving with velocity v i h . Initially, when h = 0, the positionsp i 0 , for all i ∈ {1, . . . , n}, are chosen randomly in the parameter space and the velocities v i 0 = 0. Each particle i stores its personal best pair (p i * h , J i * h ) and the social best pair (s * h , J s * h ) in memory, where J i * h = J (p i * h ) and J s * h = J (s * h ) are the costs (19) of personal best positionp i * h and social best position s * h = arg minpi * h ,i∈{1,...,n} J (p i * h ), respectively. Notice that J s * h ≤ J i * h for all i ∈ {1, . . . , n}. The personal best pair of a particle corresponds to the best position in the parameter space it has visited so far. The social best pair, on the other hand, corresponds to the best position in the parameter space that anyone in the swarm has visited so far. At every time step, each particle updates its velocity, position, its personal best pair, and the social best pair. The velocity and position are updated as follows:
where w is the inertia weight, c 1 , c 2 are the acceleration coefficients, and r h,1 , r h,2 are uniformly distributed random numbers in [0, 1] generated at each time step h. There are many ways of choosing these parameters [55, 56, 57] . Each particle i computes the cost J i h+1 = J (p i h+1 ) at its current position and updates its personal best pair as Appendix C. Decreasing property of the testable population Lemma 6. The testable population x T decreases on any interval on which x I is decreasing and θ is non-decreasing.
Proof. Let us consider an interval (t, t ), t < t on which x I is decreasing while θ is non-decreasing. First, one can note that
= θ(t )x I (t) + (1 − θ(t ))(N − x D (t) − x R (t)), (C.1) because x I (t) + x D (t) + x R (t) − N < 0 and θ is non-negative and non-decreasing by assumption. Since x I is supposed to be decreasing, then θ(t )x I (t) ≥ θ(t )x I (t ).
On the other handṄ
meaning that N − x D − x R is a decreasing function. As 1 − θ(t ) ≥ 0, one gets
Adding the two inequalities, one deduces from (C.1) that x T (t) ≥ x T (t ) whenever x I is decreasing on (t, t ). A tighter examination shows that, as both expressions θ(t ) and 1−θ(t ) cannot be zero together, at least one of the two terms of the sum indeed decreases between t and t . Therefore, x T (t) > x T (t ).