DETERMINATION OF THE OPTIMAL CONTROLS FOR AN EBOLA EPIDEMIC MODEL

. A control SEIR type model describing the spread of an Ebola epidemic in a population of a constant size is considered on the given time interval. This model contains four bounded control functions, three of which are distancing controls in the community, at the hospital, and during burial; the fourth is burial control. We consider the optimal control problem of minimizing the fraction of infectious individuals in the population at the given terminal time and analyze the corresponding optimal controls with the Pon-tryagin maximum principle. We use values of the model parameters and con-trol constraints for which the optimal controls are bang-bang. To estimate the number of zeros of the switching functions that determine the behavior of these controls, a linear non-autonomous homogenous system of diﬀerential equations for these switching functions and corresponding to them auxiliary functions are obtained. Subsequent study of the properties of solutions of this system allows us to ﬁnd analytically the estimates of the number of switchings and the type of the optimal controls for the model parameters and control constraints re-lated to all Ebola epidemics from 1995 until 2014. Corresponding numerical calculations conﬁrming the results are presented.


1.
Introduction. Ebola haemorrhagic fever is a rare and dangerous disease caused by a strain of Ebola virus. Nowadays, the epidemic outbreak in 2014 in West Africa is the largest and most complex from the moment when Ebola virus was discovered in 1976 in Zaire. The virus causes an acute heavy disease, which often ends in death. According to the World Health Organization (WHO), the mortality rate can reach 90%. By the report of the USA Center for Disease Control and Prevention (CDC), promulgated in September 2014, in the worst case the number of cases of Ebola disease exceeded one million for four months. Drugs or vaccines against the virus that are approved by FDA for Food and Drug Administration does not yet exist (symptomatic therapy is used for treatment) that leads to death in 60 − 90% of cases and the spread of Ebola disease through contact with various bodily fluids of infected people.
The mathematical modeling of Ebola is a powerful tool for studying the mechanisms of its spread. Appropriate epidemiological models are the basis for the prediction and evaluation of the dynamics of this disease. In order to curb, control, and eventually stop an Ebola epidemic, it is important to consider qualitatively adequate mathematical models. Due to advancements in computer technology and mathematical modeling it is now a realistic task.
There are various mathematical models describing diseases such as Ebola. Their classification and the corresponding references are given in [13]. Since the object of our investigations is to examine a deterministic model describing the spread of an Ebola epidemic by ordinary differential equations, we will consider compartmental models which include SIR and SEIR models as well as their modifications. Among these models, control SIR and SEIR models, and also control SEIR type models are allocated. They involve time-dependent control functions, which reflect a variety of preventive measures (monitoring, public education, contact tracing, social distancing, isolation, quarantine, screening, burial efforts) aimed at reducing contact between susceptible individuals and other community groups already infected by the Ebola virus in varying degrees (infectious, hospitalized, infectious deceased and others). For these models various optimal control problems are studied in order to find such optimal solutions, which provide a rational basis for preventive policies and design to control the spread of Ebola virus. A detailed survey of such problems and relevant references are also given in [13]. Here we allocate only papers published at the end of 2016 ( [1,8,18,25]). Note that in all these papers, the considered functionals include the total costs of the control constraints on the given time interval that are defined by integrals of the squares of the controls. This paper is a continuation of our studies presented in [11,14,20], in which the minimization of the linear terminal functional for a non-linear control system of six differential equations describing the spread of an Ebola epidemic in a population of constant size over a given time interval is considered. This system is a control SEIR type model; its equations describe the interaction between different compartments of the population such as susceptible, exposed, infectious, hospitalized, infectious deceased, and removed individuals. Each variable is a fraction of the corresponding compartment in the total population. This model also contains four bounded control functions which describe the efforts to reduce contact between the susceptible and infectious, between the susceptible and hospitalized, and finally, between the susceptible and infectious deceased individuals. The last control represents burial efforts. Our objective function is the fraction of the infected individuals at the end of the given time interval that is to be minimized.
For the analysis of optimal controls in this minimization problem, the Pontryagin maximum principle is applied. As a result, we find the linear non-autonomous homogeneous system of differential equations for the switching functions and corresponding auxiliary functions. The analysis of this system shows that the optimal controls corresponding to these switching functions are bang-bang. Estimation of the number of zeros of the switching functions is the next stage of our investigations. Here, the key point is the reduction of the matrix of the linear non-autonomous system to an upper triangular form on the entire given time interval, because the subsequent application of the generalized Rolle's theorem to each equation of the transformed linear system allows us to estimate the possible number of zeros of the switching functions. The result of this reduction procedure is a non-autonomous non-homogeneous quadratic system of ten differential equations.
Due to the large order of the quadratic system, consideration of the coefficients of its equations independently from each other and difficulty in calculations to allocate any regularities, in [11] we failed to apply effectively the technique presented in [19] in order to justify the existence of the solutions of this system defined on the entire given time interval. The existence of such solutions was established only for specific values of the parameters of the considered control SEIR type model, which correspond to the 2014 Ebola epidemic in Liberia and Sierra Leone. Successful application of the aforementioned technique presented in [12].
As a result of careful study of the considered quadratic system, we found in [20] that it can be split into several closed subsystems of differential equations. For each of these subsystems, sometimes after adding additional equations, the corresponding quadratic subsystems of the lower order can be obtained. More specifically, the ten differential equations of the original quadratic system split into four subsystems consisting of four, three, two and one differential equations, respectively. The first subsystem is initially closed, and therefore it remains unchanged; in the second and third subsystems, by adding three and four corresponding additional equations, it is possible to obtain two closed quadratic subsystems consisting of six differential equations. Finally, in the last subsystem, by adding three corresponding additional equations, the closed quadratic subsystem of four differential equations is found. After that, we apply the technique described in [12] to each of the obtained quadratic subsystems. It allows us, using less complex calculations, to find the estimates of the number of zeros of the switching functions, and therefore the corresponding estimates of the number of switchings of the appropriate optimal controls in the original problem of minimizing for a wider range of values of the parameters of the control SEIR type model. The resulting estimates enable us to reduce the considered optimal control problem to a considerably simpler problem of the finite-dimensional constrained minimization and, as a result of the corresponding numerical calculations, obtain the types of the optimal intervention strategies for the Ebola epidemics from 1995 until 2014. This paper is organized as follows. In Section 2, we introduce a control SEIR type model, describe its variables, control functions, parameters, and then study the properties of these variables. Also, we state the optimal control problem and discuss the existence of the optimal solutions in it. In Section 3, for the analysis of these solutions the Pontryagin maximum principle is used. By this principle, we see that the types of the optimal controls are determined by the corresponding switching functions. Based on the adjoint system of the optimal control problem, we find the linear non-autonomous homogeneous system of differential equations for these switching functions and corresponding to them auxiliary functions. Using this linear system, we study some properties of the switching functions and corresponding optimal controls. In particular, we find and then use only such values of the model parameters and control constraints for which the corresponding optimal controls are bang-bang. In Section 4, we investigate the switching functions in order to obtain the estimates of the number of zeros of these functions. The basis of such study consists of reducing the matrix of the linear system to an upper triangular form on the entire given time interval. To perform such a transformation, the arising non-autonomous non-homogeneous quadratic system of differential equations is investigated in order to find restrictions on the coefficients of these equations, for which there exist the corresponding solutions of this quadratic system defined on the given time interval. In order to facilitate such studies, splitting the considered quadratic system into several smaller order quadratic subsystems is used. This splitting leads to the fact that we are able to find the required restrictions and, after all, to obtain analytically the estimates of the number of switchings and the type of the optimal controls for the model parameters and control constraints relating to all Ebola epidemics from 1995 until 2014. Obtained results in previous section are the basis for the method of the solution to the optimal control problem, a description of which is given in Section 5. Additionally, the results of numerical calculations and their discussion are presented in Section 6. Finally, Section 7 contains our conclusions.
2. Control SEIR type model and its optimal control problem. Let us consider on a given time interval [0, T ] a nonlinear SEIR type model describing the spread of Ebola epidemic in a population of a constant size N , defined by the following system of differential equations: where β I , β H , β F , α, γ, σ, δ, ρ, χ, µ, T are the given positive parameters of the model and S 0 , E 0 , I 0 , H 0 , F 0 , R 0 are its initial conditions, which characterize the Ebola epidemic ( [23,29]). According with model (1), a population of size N is divided into six compartments S(t), E(t), I(t), H(t), F (t), R(t) such that Here the quantity S(t) is the number of susceptible individuals who can be infected by the Ebola virus after having contact with infectious individuals. The quantity E(t) denotes exposed individuals who have been infected by the Ebola virus but are not yet infectious. Individuals that are infected with the disease and are suffering from symptoms of Ebola are classified as infectious, the number of them is denoted by I(t). The quantity H(t) is the number of hospitalized individuals who are infectious. The quantity F (t) is the number of deceased individuals who may transmit the disease during funerals. Finally, the number of deceased and recovered individuals is denoted by the same variable R(t) because there has not been a case in which a person who survived Ebola contracts the disease again. This is the first factor that implies that the population size is constant.
In the system (1), the value β I N −1 S(t)I(t) is the number of individuals infected due to direct contact with an infected individual, β H N −1 S(t)H(t) is the number of individuals infected due to direct contact with a hospitalized individual, and, finally, the value β F N −1 S(t)F (t) is the number of individuals infected due to direct contact with a deceased individual during the funeral process. Here β I , β H , β F are the transmission rate in the community, the transmission rate at the hospital, and the transmission rate during traditional funerals, respectively. The value αE(t) is the number of individuals in the exposed stage, which show the symptoms of the disease and pass on to the infectious stage; α is the infectious rate. The values γI(t), δI(t), σI(t) are the number of individuals in the infectious stage, which are hospitalized, die and are buried, die and are not yet buried, respectively. Here γ, δ, σ are the corresponding rates of hospitalization, mortality with burial, and mortality without burial. The values ρH(t), χH(t) are the number of individuals in the hospitalization stage, which die and are not yet buried, and die and are buried, respectively. Here ρ, χ are the corresponding rates of mortality of hospitalized without burial and mortality of hospitalized with burial. Finally, the value µF (t) is the number of individuals in the funeral process which are buried; µ is the traditional burial rate. The equations of system (1) have no terms related to natural mortality or fertility because of the short time-span of an Ebola epidemic. This is the second factor which explains the validity of the assumption of constancy of the population size.
The control interventions of public health measures are crucial factors for stopping Ebola transmission. Therefore, in order to control the spread of infection and prevent the Ebola epidemic, we add to model (1) four bounded control functions u(t), v(t), w(t), η(t). The controls u(t), v(t), w(t) represent the efforts that reduce the contact between the susceptible and infectious individuals, between the susceptible and hospitalized, and, lastly, between susceptible and infectious deceased, respectively. The control η(t) represents the burial efforts. These controls satisfy the following constraints: Then, we have the following SEIR type control model: For the phase variables S(t), E(t), I(t), H(t), F (t), R(t) of this system equality (2) is still valid. Next, in system (4) it is convenient to replace the variables S(t), E(t), I(t), H(t), F (t), R(t) to the new variables s(t), e(t), i(t), h(t), f (t), r(t), defined as with the corresponding initial values: Variables (5) are the fractions of the compartments S(t), E(t), I(t), H(t), F (t), R(t) in a population of size N , and therefore, by equality (2), we have the formula: Also, it is convenient to consider a system for the variables s(t), e(t), i(t), h(t), f (t); then, by formula (6), the variable r(t) can be found from the equality: Then, for the variables s(t), e(t), i(t), h(t), f (t) we obtain the following nonlinear control system: For this system the set of all admissible controls is formed by all possible Lebesgue measurable functions u(t), v(t), w(t), η(t), which for almost all t ∈ [0, T ] satisfy constraints (3). Now, we introduce the region: The following lemma, which establishes the boundedness, positiveness and continuation of solutions for system (7) is valid.
Lemma 2.1. Let the inclusion (s 0 , e 0 , i 0 , h 0 , f 0 ) ∈ Λ hold. Then for any admissible controls u(t), v(t), w(t), η(t) the corresponding solutions s(t), e(t), i(t), h(t), f (t) for system (7) are defined on the entire interval [0, T ] and satisfy the inclusion : Relationship (8) implies that the region Λ is a positive invariant set for system (7). The proof of Lemma 2.1 is fairly straightforward and we omit it. Proofs of such statements are given, for example, in [9,21].
For system (7), on the set of admissible controls we consider the problem of minimization of the fraction of the infected individuals at the terminal time T : We note that the optimal control problem (7),(9) differs from the problems that are typically considered in the literature on control of SEIR models in that the functional in (9) does not include integrals of controls u(t), v(t), w(t), η(t) or their squares, which are responsible for the total cost of all epidemiological measures. This is equivalent to an assumption that when we use controls u(t), v(t), w(t), η(t) within given bounds (3), the costs of epidemic control policies are acceptable and do not influence the functional. The existence in problem (7),(9) of the optimal controls u * (t), v * (t), w * (t), η * (t) and the corresponding optimal solutions s * (t), e * (t), i * (t), h * (t), f * (t) follows from Lemma 2.1 and Theorem 4 ( [22], Chapter 4).
Let us introduce the auxiliary functions:
Assumption 1. Let the parameters of system (7) and control constraints (3) satisfy the inequalities: Then, the following lemma can be stated. The proof of this lemma is fairly straightforward and based on the arguments from the opposite, which contradicts to the Assumption 1. A detailed proof of Lemma 3.1 is given in [14]. Corollary 1. Formulas (11)- (14) and Lemma 3.1 imply that the optimal controls u * (t), v * (t), w * (t), η * (t) are bang-bang functions taking values u min ; u max , v min ; v max , w min ; w max , η min ; η max , respectively, for almost all t ∈ [0, T ]. These controls have no singular portions (see [30] for details) and they are entirely determined by the corresponding switching functions L 0 (t) and L η (t). Furthermore, formulas (11)- (13) show that the controls u * (t), v * (t), w * (t) switch from maximum to minimum values and vice versa at the same moments of switching.
Next, the following lemma can be stated.
There exists a value t 0 ∈ [0, T ) such that the inequalities : hold for all t ∈ (t 0 , T ).
Proof. By the second initial condition of system (15) and the continuity of the auxiliary function G(t), there exists such a value t ∈ [0, T ) that the inequality: is true for all t ∈ (t , T ]. Let us integrate the first differential equation of this system together with the corresponding initial condition. We obtain the formula: which, by inequality (17), implies the existence of value t 0 ∈ [0, t ) that for all t ∈ (t 0 , T ] the first inequality in (16) holds. Now, we integrate the fifth differential equation of system (15) together with the corresponding initial condition. As a result, we find the formula: from which, by the first inequality in (16), the validity of the second inequality for the specified values of t follows. This completes the proof.
Corollary 2. Formulas (11)- (14) and Lemma 3.2 yield the following relationships for the optimal controls u * (t), v * (t), w * (t), η * (t) : Let us study in detail the properties of zeros of the switching functions L 0 (t) and L η (t). The following lemma holds. Proof. Let τ ∈ (0, T ) be a zero of the function L 0 (t), that is Substituting this equality into the first equation of system (15), we have the formula: In this expression, if G(τ ) = 0, then the required fact is established. Indeed, the inequality L 0 (τ ) = 0, following from (20), implies the existence of a neighborhood of τ at which except the point τ itself the function L 0 (t) takes only nonzero values. Now, let us consider that G(τ ) = 0. Using this equality and (19) in the second equation of system (15), we find that Again, we have two possible situations. If Q(τ ) = 0, then applying arguments, which are similar to those presented above, we find that τ is an isolated zero of the function G(t), and hence the function L 0 (t). The required fact is established. Now, let us consider that Q(τ ) = 0. (21) Then, using equalities (19) and (21) in the third equation of system (15), we find the formula: . Again, there are two possible situations. If γP (τ ) − ση * (τ )L η (τ ) = 0, then, as before, by Q (τ ) = 0, we conclude that τ is the isolated zero of the function Q(t). Therefore, it is the isolated zero of G(t) and hence the function L 0 (t). This is what we wanted to establish. Now, let us consider that If P (τ ) = 0, then from (22) we find that L η (τ ) = 0. This leads to the fact that the linear homogeneous system (15) has on the interval [0, T ] only a trivial solution. It contradicts the corresponding initial conditions of this system. Similar arguments are true when L η (τ ) = 0.
Thus, further we consider that simultaneously the following inequalities hold: Function L η (t) is absolutely continuous, and therefore there exists a neighborhood ∆ τ of τ at which this function retains its sign. Hence, on this interval the control η * (t) takes only one of the possible values {η min ; η max }. By the corresponding equations of system (15), we obtain on the interval ∆ τ the differential equation for the function In addition, by (22), it is easy to see that for the function Y (t) the following equality holds: Let us consider the equation (23) at the value t = τ . Using (19), (24) and the definition of the constant λ, we obtain: By the inequalities of Assumption 1, we find that Y (τ ) = 0. This means that τ is the isolated zero of the function Y (t). Then, using on the interval ∆ τ the corresponding equations of system (15) in turn, we find that τ is the isolated zero of the function L 0 (t). This completes the proof.   Proof. The function L 0 (t) is absolutely continuous, and therefore the set Z 0 of its zeros is closed. Hence, the set [0, T ] \ Z 0 is a bounded open set that can be represented as a union of not more than a countable number of disjoint intervals (see [2]): At each interval (a k , b k ) the function L 0 (t) takes values of the same sign, either positive or negative; the endpoints of the intervals (a k , b k ), k ≥ 1 are zeros of this function. If the number of such intervals is finite, then the number of zeros of the function L 0 (t) is finite as well, that simultaneously means the finiteness of the number of elements of the set Z 0 . If the number of the intervals (a k , b k ), k ≥ 1 is countable infinite, then both the number of zeros of the function L 0 (t) and the number of terms in the series will be countable infinite as well. Series (25) defines the Lebesgue measure of the set [0, T ] \ Z 0 . Since it is finite, this series converges. Hence, its common term (b k − a k ) tends to zero when k → +∞, that is Let us note that the set of limit points of the sequence {a k } k≥1 is nonempty and consists of either one element, {0} or {T }, or two elements {0; T }. Indeed, we will show that there is not a point of the interval (0, T ) that would be a limit point of this sequence. Let us assume the contrary, that is, there exists a subsequence {a km } of the sequence {a k } k≥1 converging to a point t 0 ∈ (0, T ). Then, by (26), the corresponding subsequence {b km } converges to t 0 as well. Since, as it is noted above, the endpoints of the intervals (a km , b km ) are zeros of the function L 0 (t), then the sequence {t km } of these zeros also converges to the point t 0 . Closedness of the set Z 0 implies the validity of the inclusion t 0 ∈ Z 0 , which contradicts Lemma 3.3.
The assumption was wrong, and only one of the values of {0; T }, or simultaneously both can be the limit points of the sequence {a k } k≥1 . Now, let us show that: Case (i). If 0 is the unique limit point of the sequence {a k } k≥1 , then we can arrange the intervals (a k , b k ), k ≥ 1 in descending order, that is a 1 > a 2 > . . . ; Case (ii). If T is the unique limit point of the sequence {a k } k≥1 , then we can arrange the intervals (a k , b k ), k ≥ 1 in ascending order, that is a 1 < a 2 < . . . ; Case (iii). If the points 0 and T are at the same time the limit points of the sequence {a k } k≥1 , then we can arrange the intervals (a k , b k ), k ≥ 1 in both descending and ascending order, dividing them into two groups. The intervals of the first group, tending to 0, will be rearranged in descending order, that is a 1 > a 2 > . . . , and the intervals of the second group, tending to T , will be arranged in ascending order, that is a 1 < a 2 < . . . . Next, let us consider each case separately.
Case (i). Let us fix an arbitrary number ∈ (0, 0.5T ). For this there is a number K such that for all k > K the inequality a k < is valid. On the interval [0, ) we introduce a set of disjoint intervals [(n + 1) −1 , n −1 ), n ≥ 1. It is obvious that Let us note that only a finite number of elements of the sequence {a k } k≥K will fall into each interval. Otherwise, considering the interval [(n 0 + 1) −1 , n −1 0 ) for some n 0 as an example, we would again have a subsequence {a km } converging to a certain value t 0 ∈ [(n 0 + 1) −1 , n −1 0 ] ⊂ (0, T ). Since each a km is both the left endpoint of the interval (a km , b km ) and a zero of the function L 0 (t), we again have a sequence {t km } of zeros of this function that converges to its zero t 0 ∈ (0, T ). This fact contradicts Lemma 3.3. As a result, we arrange in decreasing order the endpoints a k of the intervals (a k , b k ) that fall into each interval [(n + 1) −1 , n −1 ). Since a number of such intervals is countable infinite, in this way we arrange the intervals (a k , b k ), k ≥ K in descending order, that is a K+1 > a K+2 > . . . . Also, in descending order we arrange the endpoints a k of the intervals (a k , b k ), k = 1, K. Finally, we obtain the intervals (a k , b k ), k ≥ 1, arranged in descending order, that is a 1 > a 2 > . . . . Each zero τ ∈ (0, T ) of the function L 0 (t) is simultaneously the right endpoint of the interval (a m+1 , b m+1 ) and the left endpoint of the interval (a m , b m ) for some m. Since the left endpoints of the intervals (a k , b k ), k ≥ 1 are arranged in descending order, by numbering the zeros of the function L 0 (t) in the appropriate way, we obtain a decreasing sequence of elements of the set Z 0 that accumulate at the point 0.
Case (ii). The arranging of the intervals (a k , b k ), k ≥ 1 in ascending order is carried out as in Case (i). It is easy to see that each zero τ ∈ (0, T ) of the function L 0 (t) is simultaneously the right endpoint of the interval (a m , b m ) and the left endpoint of the interval (a m+1 , b m+1 ) for some m. Since the left endpoints of the intervals (a k , b k ), k ≥ 1 are arranged in ascending order, by numbering the zeros of the function L 0 (t) in the appropriate way, we obtain an increasing sequence of elements of the set Z 0 that accumulate at the point T . Case (iii). Again, we fix an arbitrary number ∈ (0, 0.5T ). For this there is such a number K that for all k n , k m > K the inequalities: a kn < , a km > T − are simultaneously valid. Further, according to the algorithm described in Case (i), we first arrange in descending order the intervals of the first group, and then in ascending order the intervals of the second group. Finally, we arrange the intervals (a k , b k ) whose left endpoints a k fall into the interval [ , T − ], relating each such interval to the first group, if the inequality a k ≤ 0.5T holds, and to the second group, if the opposite inequality is true. Further, repeating the arguments of Cases (i) and (ii), we find two sequences of elements of the set Z 0 : the elements of one decreasing accumulate at 0, and the elements of the other increasing accumulate at T . This completes the proof.    [11,13,14,15,16]. Now, we assume that the following condition holds. Assumption 2. Let at a moment of switching θ * ∈ (0, T ) the values of the controls u * (t), v * (t), w * (t), and at a moment of switching τ * ∈ (0, T ) the value of the control η * (t) be equal to the corresponding limits from the left, that is Besides, let these controls be continuous at the ends of the interval [0, T ].
Analyzing the equations of system (15), we see that its third equation is different from the other equations. The difference is that the righthand side of the equation contains the function Q(t), for which this differential equation is written, and the other two functions P (t), L η (t) of this system. Whereas, the right-sides of other equations contain only the function, for which the differential equation is written, and at least one of the following functions of system (15). If this drawback is not eliminated, the subsequent application of the Generalized Rolle's Theorem ( [7]) to each equation of the system for estimating the number of zeros of the switching function L 0 (t) will be problematic. In order to eliminate this drawback, we introduce into system (15) a new auxiliary function: Then, the third equation of the system can be written as In order to obtain the differential equation for the function P (t), we have to be sure that P (t) is differentiable function almost everywhere on the interval [0, T ]. Functions P (t), L η (t) are absolutely continuous functions as the solutions of the corresponding differential equations of system (15). Moreover, by Assumption 2 and Theorem 3 ( [22], Chapter 4), the function η * (t)L η (t) is the absolutely continuous function as well. Finally, Corollary 5 guarantees the almost everywhere differentiability of the function η * (t). Thus, the function P (t) is almost everywhere differentiable on the interval [0, T ], and hence it satisfies the differential equation: Next, we omit the tilde at the function P (t) in equations (27) and (28), and then replace by them the corresponding equations of system (15) together with the appropriate initial conditions. As a result, we obtain the system: which further will be analyzed. (29) is a linear non-autonomous homogeneous system of differential equations, the matrix of which has the lower Hessenberg form:

Remark 2. System
Thus, we transformed the matrix of the linear system (15) to the such type. We note that linear non-autonomous systems with matrices in this form are used, for example, in the observation theory. Another way of reducing the matrix of a linear non-autonomous system to the Hessenberg form is presented in [3,4]. 4. Estimation of the number of zeros. Our goal now is to estimate the number of zeros of the switching functions L 0 (t). For this, we apply to the linear nonautonomous system of differential equations (29) the technique developed earlier in [12,19].
Let us reduce the matrix Θ(t) of system (29) to an upper triangular form: This procedure is called a triangulation and consists of the consecutive execution for such a system the following replacements of the variables: Here, the functions q lj (t), l = 5, 2, j = 1, (l − 1) satisfy a non-autonomous quadratic system of differential equations: The first replacement of the variables in (30) allows us to make zero elements of the fifth row of the matrix Θ(t) from the first and to prior entry located on the main diagonal. In this case, the fourth row also changes, but not essentially. The remaining rows of this matrix do not change at all. In the transformed matrix, the second replacement of the variables in (30) makes zero elements of the fourth row from the first to the prior of the entry located on the main diagonal. Again, the third row also changes, but not essentially. The remaining rows of this matrix do not change. This process will continue as long as we do not change by the aforementioned way the second row of the matrix. The first row we do not change because from the very beginning it had the required form.
Here the sign −1 means inversion.
Due to the large order of the quadratic system (31), consideration of the coefficients of its equations independently from each other and difficulty in calculations to allocate any regularities, in [11], we failed to apply effectively the technique presented in [12,19] in order to justify the following important lemma. The existence of such solutions was established only for the specific numerical values of the parameters for system (7) and the control constraints (3). We used values summarized in Table 1 below, which were used for modeling of the 2014 Ebola epidemics in Sierra Leone and Liberia ( [10,29]). Successful application of the mentioned technique presented in [12].
In order to overcome this drawback, we focused on the quadratic system (31). As a result of its careful study, we found that it can be split into four closed subsystems of differential equations. Further, for each of these subsystems, sometimes after adding additional equations, the corresponding quadratic subsystems of the lower order can be obtained. After that, we were able to establish the validity of Lemma 4.1 based on these subsystems for a wider range of values of the parameters for system (7) and the control constraints (3).
Remark 3. Of course, the equations of the non-autonomous quadratic system (31) can be rewritten aṡ where r(t) = (r 1 (t), r 2 (t), . . . , r 10 (t)) ; A j (t) is a square symmetric matrix, b j (t) is a vector, and c j (t) is a scalar, j = 1, 10. The results, presented in [11], showed that such representation of the equations of this quadratic system has led to the "impoverishment" of the results (they were obtained only for Liberia and Sierra Leone). So, further, we will study the quadratic systems (47)-(50). Now, let us study the quadratic subsystem (47). We assume the contrary, i.e. whatever a solution q 5,0 (t) for system (40) we do not consider, it is defined on the interval [0, t 5 1 ), where t 5 1 ≤ T , and this interval is the maximum possible interval of the existence of such a solution ( [17,22]). Then, from Lemma ( [5], § 14, Chapter 4) for solution q 5,0 (t) it follows the relationship: where ||·|| is the Euclidean norm of a vector. This relationship implies the existence of the number φ > 0 and the value t 5 ∈ [0, t 5 1 ) for which the inequality ||q 5,0 (t)|| ≥ φ holds for all t ∈ [t 5 , t 5 1 ). Note that the values φ and t 5 will be defined below. Now, on the interval [t 5 , t 5 1 ) we evaluate the derivative of the function ||q 5,0 (t)|| with respect to system (47). We have the equality: where A 1 , B 1 , C 1 are positive constants depending on the parameters of system (7) and the control constraints (3) by the following formulas: which are obtained using Lemma 2.1 and the control constraints (3). Now, we consider the quadratic equation: and then require the positivity of its discriminant, i.e.
Thus, for the considered numerical values of the parameters of system (7) and the control constraints (3)  Now, we consider for system (33) the different cases related to the inequalities from Assumption 1.
Case 1. Let one of the following two inequalities: hold. It means that one of the following two relationships: holds as well. Therefore, the function g(t) takes only either positive values or negative values on the interval [0, T ]. Then, the following lemma can be stated. Proof. We assume the contradiction. Let the function have on the interval [0, T ] at least five distinct zeros θ 1 j , j = 1, 5, where 0 ≤ θ 1 1 < θ 1 2 < θ 1 3 < θ 1 4 < θ 1 5 ≤ T . Applying the Generalized Rolle's Theorem ( [7]) to the first differential equation of system (33), we conclude that the functionG(t) has on the interval (0, T ) at least four distinct zeros θ 2 j , j = 1, 4, where 0 < θ 2 1 < θ 2 2 < θ 2 3 < θ 2 4 < T . Applying now this theorem to the second differential equation from (33), we find that the functionQ(t) has on the interval (0, T ) at least three zeros θ 3 j , j = 1, 3, where 0 < θ 3 1 < θ 3 2 < θ 3 3 < T . Applying again the Generalized Rolle's Theorem ( [7]) to the third differential equation of system (33), we obtain that the functionP (t) has on the interval (0, T ) at least two distinct zeros θ 4 j , j = 1, 2, where 0 < θ 4 1 < θ 4 2 < T . Finally, applying this theorem to the fourth differential equation from (33), we find that the functionL η (t) has on the interval (0, T ) at least one zero θ 5 1 . However, this function satisfies the linear homogeneous differential equation, and thereforeL η (t) = 0 everywhere on the interval [0, T ]. Considering then the fourth equation of system (33), we conclude that the functionP (t) is also zero everywhere on the interval [0, T ]. Similarly, considering initially the third equation, and then the second equation of this system, we find thatQ(t) = 0 andG(t) = 0 everywhere on the interval [0, T ], respectively. Therefore, using the first equation of system (33), the similar conclusion is obtained for the functionL 0 (t), i.e.L 0 (t) = 0 for all t ∈ [0, T ]. By equality (62), we find that L 0 (t) = 0 everywhere on the interval [0, T ]. This fact contradicts Lemma 3.2. Hence, our assumption was wrong, and the switching function L 0 (t) has at most four distinct zeros on the interval [0, T ]. The proof is complete. Now, using the first initial condition of system (29) and Lemma 4.2, we find that the switching function L 0 (t) has on the interval [0, T ) at most three distinct zeros. By the corresponding relationships from (18) and formulas (11)-(13), we obtain the type of the optimal controls u * (t), v * (t), w * (t) as where θ * j ∈ [0, T ), j = 1, 3 are the corresponding moments of switching. Next, let us apply the Generalized Rolle's Theorem ( [7]) to the last differential equation of system (29). Using the estimation of the number of zeros of the function L 0 (t) from Lemma 4.2 and the last initial condition of system (29), we conclude that the switching function L η (t) has on the interval [0, T ) at most three distinct zeros. By the corresponding relationship from (18) and formula (14), we find the type of the optimal control η * (t) as follows where τ * j ∈ [0, T ), j = 1, 3 are the corresponding moments of switching. Case 2. Let the inequalities: hold. It means that the inclusion: holds as well. Therefore, the function g(t) takes positive values at η * (t) = η max and negative values at η * (t) = η min . Then, the product g(t)L η (t) is the continuous function and takes positive values, when η * (t) = η min , or η * (t) = η max . Moreover, at those points, which are the zeros of the function L η (t), the product g(t)L η (t) is zero as well. Then, we apply arguments, which are similar to those presented above, to the first four equations of the system (29). As a result, we obtain the triangular system similar to (33), but consisting of four differential equations and having the product g(t)L η (t) in the last equation. Now, applying to this system arguments, which are similar to the arguments from Lemma 4.2, and then the arguments that follow this lemma, we find that the switching function L 0 (t) has on the interval [0, T ] at most four distinct zeros. Therefore, the following arguments and conclusions are similar to the previous case. Thus, in Case 2 for the optimal controls u * (t), v * (t), w * (t), and η * (t) we have formulas (63) and (64), respectively. Remark 4. Using the perturbation method, it can be shown that it is possible to abandon Assumption 1 of our previous arguments. Indeed, let us suppose that this condition is violated. Then we introduce a small positive perturbation parameter in system (7). Thus, we have a perturbed optimal control problem in which Assumption 1 is satisfied and the appropriate types of the corresponding optimal controls are known from the previous arguments. Let the perturbation parameter approach zero. Then the corresponding sequence of the piecewise constant optimal controls of the perturbed problem converges to the piecewise constant controls of the unperturbed problem. To show that the constructed controls also are optimal in the unperturbed problem, we investigate the properties of the upper and lower semi-continuity of the optimal value of the functional in the perturbed problem. A detailed justification of the discussed fact presented in [14].
Thus, the optimal control problem (7),(9) is now reduced to a problem of constrained minimization: min This problem is considerably simpler than the optimal control problem (7),(9) and it is solved numerically. In order to find the global minimum of the function F (θ, τ ), the overlapping method is used ( [32]) because the Lipschitz constant of this function is easily found (see for example [12]). The corresponding numerical algorithm is written in MAPLE.
Professor H. Maurer from University of Muenster (Muenster, Germany) kindly conducted numerical calculations for the same initial conditions and parameters using the optimal control package NUDOCCCS due to C. Büskens (Bremen) (formulation with AMPL and interior-point solver IPOPT) (see [24,26]), and got following more precise results. The optimal controls u * (t), v * (t), w * (t) and η * (t) are the piecewise constant functions as well. But the first three controls, u * (t), v * (t), w * (t) switch at the time θ * 3 = 28.4797 from the maximum values u max , v max , w max to the minimum values u min , v min , w min (θ * j = 0, j = 1, 2). The last control, η * (t), switches at the time τ * 3 = 27.3521 from the minimum value η min to the maximum value η max (τ * j = 0, j = 1, 2). The corresponding smallest value of the functional is J(u * , v * , w * , η * ) = 18.6985.
However, regarding Liberia we also calculated the dynamics of system (4) when all controls are constant: This means that applying all precocious measures at maximum does not improve the outcome. The value of the functional in (9) was 19.4621 that was greater than that obtained under the optimal controls.
Here we will not discuss the epidemiological sense of the type of the obtained optimal controls u * (t), v * (t), w * (t), and η * (t) because they do not represent actual intervention strategies and are only controls suitable for analytical studies. The corresponding conclusions that have epidemiological explanations are presented in [14]. 7. Conclusions. On a given time interval we considered a control SEIR type model, which describes the spread of Ebola epidemic in a population of a constant size. This model contained four bounded controls: one burial control, and three distancing controls in the community, at the hospital and during burial. These controls were introduced into the model in order to investigate optimal intervention strategies for minimizing the number of the infectious individuals at the end of the pre-defined period of time. We used such values of the model parameters and control constraints, for which the corresponding optimal controls were bang-bang. An approach for estimating the number of zeros of the corresponding switching functions, different from the one that was used in our previous papers, was applied. Having the maximal possible number of switchings, we were finally able to reformulate the original optimal control problem as a problem of constrained minimization, which solved numerically by standard software. The aforementioned approach allowed us to expand the range of the used values of the parameters of the control SEIR type model and its control constraints and to find the type of the optimal intervention strategies for all Ebola epidemics from 1995 until 2014. Corresponding results of numerical calculations were presented. Figure 2. Graphs of the optimal solutions for the Ebola epidemic in Sierra Leone: top row: S * (t), E * (t); middle row: I * (t), H * (t); bottom row: F * (t), R * (t).