Research Article | Open Access

# Optimal Controls of the Highly Active Antiretroviral Therapy

**Academic Editor:**Valery Y. Glizer

#### Abstract

In this paper, we study generic properties of the optimal, in a certain sense, highly active antiretroviral therapy (or HAART). To address this problem, we consider a control model based on the -dimensional Nowak–May within-host HIV dynamics model. Taking into consideration that precise forms of functional responses are usually unknown, we introduce into this model a nonlinear incidence rate of a rather general form given by an unspecified function of the susceptible cells and free virus particles. We also add a term responsible to the loss of free virions due to infection of the target cells. To mirror the idea of highly active anti-HIV therapy, in this model we assume six controls that can act simultaneously. These six controls affecting different stage of virus life cycle comprise all controls possible for this model and account for all feasible actions of the existing anti-HIV drugs. With this control model, we consider an optimal control problem of minimizing the infection level at the end of a given time interval. Using an analytical mathematical technique, we prove that the optimal controls are bang-bang, find accurate estimates for the maximal possible number of switchings of these controls and establish qualitative types of the optimal controls as well as mutual relationships between them. Having the estimate for the number of switchings found, we can reduce the two-point boundary value problem for Pontryagin Maximum Principle to a considerably simpler problem of the finite-dimensional optimization, which can be solved numerically. Despite this possibility, the obtained theoretical results are illustrated by numerical calculations using BOCOP–2.0.5 software package, and the corresponding conclusions are made.

#### 1. Introduction

##### 1.1. Antiretroviral Drugs and Guidelines

The World Health Organization and the United Nations Programme on HIV/AIDS (UNAIDS) estimated that over millions of people were living with HIV/AIDS in 2017; other million were newly infected with HIV in the same year [1, 2]. Of these millions of the infected, millions were living on antiretroviral therapy [1, 2]. These figures mean that the therapy is currently covers less than percent of the infected. One (and probably the major) of the reasons for the insufficient coverage is the high cost of the antiretroviral therapy, which ranges from almost US $ per patient per year in the United States (where costs of the therapy are the highest) [3] to US $– per a person per year in other countries [4] (p. 90). Such costs make the treatment inaccessible for the majority of patients in developing countries, where more than percent of infections and HIV-associated deaths occurred.

Starting from 1987, when the first antiretroviral drugs, Nucleoside and Nucleotide Reverse Transcriptase Inhibitors, were approved for HIV treatment, over of antiretroviral drugs were introduced in medical practice. A typical action of these drugs is inhibition of one of the several phases of HIV reproduction cycle. Depending on the stage that is inhibited, the drugs can be broadly classified into one of the following six groups [1, 2, 5]:(1)Entry or Fusion Inhibitors (Enfuvirtide, Maraviroc) Interfere with Binding, Fusion, and Entry of Virus to the Host Cell by Blocking One of Several Targets.(2)Nucleoside reverse transcriptase inhibitors (NRTI) and nucleotide reverse transciptase inhibitors (NtRTI) is the oldest and the largest class of antiretroviral drugs that includes drugs such as Zidovudine, Didanosine, Stavudine, Lamivudine, Abacavir, Tenofovir, and combinations of these. The NRTI and NtRTI inhibit reverse transcription by being incorporated into the newly synthesized viral DNA strand as faulty nucleotides.(3)Nonnucleoside reverse transcriptase inhibitors, or NNRTI, (Nevirapine, Delavirdine, Efavirenz, Etravirine, Rilpivirine) also inhibit reverse transcriptase, but in a different way, by binding to an allosteric site of the enzyme.(4)Integrase inhibitors (Raltegravir, Elvitegravir, Dolutegravir) inhibit the enzyme integrase, which is responsible for integration of the viral DNA into the DNA of the infected cell.(5)Protease inhibitors (PI) inhibit the activity of protease, an enzyme that cut nascent proteins into specific pieces for the final assembly of new functioning virions. There are ten approved protease inhibitors, namely: Saquinavir, Indinavir, Ritonavir, Nelfinavir, Amprenavir, Lopinavir/Ritonavir, Atazanavir, Fosamprenavir, Tipranavir, and Darunavir.(6)Maturation inhibitors inhibit the last step in gag processing where the viral capsid polyprotein is cleaved, blocking the conversion of the polyprotein into the mature capsid protein. As a result, the virions released have a defective core and are mainly noninfectious.

The guidelines for use of the antiretroviral drugs are a source of controversy within the medical community. Thus, there is diversity of opinion regarding the stage of infection and the CD4+ T cells count (which are the primary target cells for HIV, and which depletion leads to the immune system collapse and development of AIDS) when the therapy should be started. It is generally accepted that the objectives of the treatment are to maximally suppress the viral load and to maintain CD4+ T cells at an acceptable level. Other factors that have to be also taken into consideration are: (i) a possibility of individual intolerance and adverse side-effects; (ii) a possibility of developing drug resistance for a patient, and (iii) a possibility of emergence of drug-resistant strains, which thereby can be passed onto others. Moreover, the cost of the therapy, which is high, should be also taken into consideration. To reduce the possibility of developing drug resistance and of the emergence of drug-resistant strains and to minimize the side effects, a combination of several antiretroviral drugs (typically, of three or four) is usually used. This approach is known as highly active antiretroviral therapy, or HAART. Also, it is generally accepted that, in order to reduce a possibility of emergence of drug resistant strains, antiretroviral therapy, once started, should never be stopped: as current US Guidelines states, “Patients initiating antiretroviral therapy should be willing and able to commit to lifelong treatment” [6].

The above-mentioned controversy in formulating the guidelines implies that some sort of optimization can be applied, and that the problem should be quantified and explored mathematically. Of course, at the present stage, it would be hardly possible to take into consideration details such as individual sensitivity of a specific patient to a specific drug. However, some common principles that can assist in developing better guidelines can be established. A natural approach for this sort of problems is to use the tools and methods of the classic optimal control theory. The problem of HIV control and the applicability of the optimal control theory to this problem attracted attention of mathematicians fairly early: seminal works by Butler et al. [7] and Kirschner et al. [8] initiated extensive literature in this direction. A mathematical techniques suggested in these publications proved to be successful and was subsequently applied to a variety of mathematical model of HIV (see, e.g., [9–22] and bibliography therein). These results brought important insight into the problem. However, since this advance utilized essentially the same mathematical techniques, all these results also share some drawbacks, that we have to discuss before we proceed to our study.

##### 1.2. Form of the Objective Functional

Starting from [7, 8], all the above-mentioned publication (and, in fact, the vast majority of the literature on the optimal control in mathematical biology and medicine) deal with optimal control problems of maximization of the cumulative (integral) concentration of healthy CD4+ T cells and simultaneous minimization of the cost of the treatment (represented by the weighted squares of the considered controls). Control problems with such functionals are mathematically convenient because they never deal with singular control phenomenon and can be solved numerically as a two-point boundary value problem. A setback of such objective functionals is that the optimal controls, obtained for such objective functionals, are very sensitive to a specific form of the functional and, in particular, to the value of the exponent. That is, the optimal controls obtained for such objective functionals are not robust, and a small variation of the exponent or even the weights leads to significant changes of the optimal controls. Furthermore, the above-mentioned approach does not guarantee the optimality of the obtained solutions. Instead, it ensures only necessary conditions of the optimality and a local extremum of the functional.

At the same time, while the objective functionals that are independent of the controls are free of the mentioned drawbacks, their use involves dealing with considerably more difficult problems, as the control-independent functionals usually lead to bang-bang controls and can have singular arcs. There are also other issues of mathematical nature that are related to a form of the objective functional. A comprehensive discussion of mathematical aspects of this problem can be found in recently published book by Schättler and Ledzewicz [23] (p. 48–50). Nevertheless, despite this complexity, there is a certain progress in application of the control-independent functionals in biomedical control problems. Results that affirmatively demonstrate that the use of such functionals is possible and potentially highly beneficial, include publications of Ledzewicz, Schättler and their coauthors on the optimal control in anticancer therapy (see, e.g., [23, 24] and bibliography therein), and of the authors on the optimal waste water biotreatment [25, 26–28] and on the optimal control of epidemics [29–34].

The objective functionals that include squares of the controls naturally occur in problems originated in mechanics and engineering, where square of a control is usually proportional to the energy required for the control action. A typical motivation for the use of the quadratic objective in HIV treatment (and in biomedical applications in general) is minimizing the “costs” in terms of both financial costs and the adverse side effects. That is, this specific form of objective functional implies that the cost and side effects increase as squares of the drug concentrations in patient’s bloodstream. However, while there is little doubt that the severity of side effects grows faster than linear with the growth of drug concentration, and, hence, an exponential functional response may be a good representation for this, the actual value of the exponent is uncertain. Moreover, it depends on patient’s individual response and varies for different patient and different drugs. This uncertainty combined with the outcome sensitivity to the exponent implies that for clinical practice the value of results, obtained for this particular type of objective functional, is questionable. Furthermore, the quadratic “cost” is hardly able to adequately represent the financial cost of the therapy either.

Taking into consideration the uncertainty of the form of the real-life “cost” function and the sensitivity of the results to this, we consider an optimal control problem for HIV treatment disregarding the “costs”. Instead, we assume that the “costs”, in terms of both the side effects and the financial cost of the therapy, is limited by the maximal acceptable levels of medication. An argument in support of such a formulation is that the actual objectives of the therapy are: (i) prevention of the development of AIDS by a patient, and (ii) control of the epidemic in a society. For an individual patient, these objectives imply reducing virus load to an acceptable safe level. Considering the severity of the possible consequences of the control failure, individual “costs” of the therapy can be disregarded. We also assume that the both, financial cost of the therapy and the side effects are known and accepted prior the start of the therapy.

An initial advance in this direction was made in [35], where the authors considered an optimal control problem for a reasonably simple -dimensional antiretroviral therapy model. The model in this paper was based on the -dimensional Wodarz model of HIV dynamics [36], where an equation representing the concentration of a drug inhibiting the incidence rate was added. The drug intake was considered as the single control. The optimal control problem for the model with a control-independent objective functional was state and solved, and the principal properties of the optimal controls were established.

##### 1.3. Interaction of Several Simultaneously Acting Controls

The majority of above-mentioned papers on the optimal anti-HIV therapy consider models with a single control, which correspond to the therapies where a single drug (a single action) is used. Very few of existing publications consider a possibility of simultaneous use of two drugs [15, 20, 37–40], whereas works exploring interaction of three or more controls are rather exceptional [41, 42]. On the other hand, the current anti-HIV therapy guidelines necessitate the use of the highly active antiretroviral therapy (HAART) that implies administering a combination of several (typically, three or four) antiretrovirals. Moreover, over thirty anti-HIV drugs are currently available, which, by their action, can be classified into six groups that affect six different stages of viral life cycle. (NRTI and NNRTI are in different groups but affect the same stage of viral replication). Such a surplus of drugs’ names and actions implies a large number of their possible combinations. Currently, the drug combinations used in HAART are arranged with the aims of (i) minimizing the adverse side effects and (ii) preventing the development of drug resistance. However, the existing surplus of the drugs and their actions presents an opportunity for designing combinations, as well as schedules of their administering, aimed to optimize the HAART efficacy as well. Apparently, this implies the need to explore interaction and interference of the six (that is, equal to the number of the HIV replication cycle stages that the existing antiretroviral drugs can inhibit) simultaneously applied antiretroviral actions.

##### 1.4. Nonlinearity of Functional Responses

The use of a numerical technique in papers [7–10, 12–22] means that in each case results are obtained for a model with a specific parametrization. However, due to intrinsic complexity of biological systems, the actual forms of functional responses and the parametrization are usually unknown in detail, and one hardly can expect that the information available on a bio-system is, or will be, sufficient. This implies that outcome of a study can depend on factors, which are neither studied, nor understood with a sufficient degree of certainty. This, in turn, leads to a question about the reliability and relevance of the results obtained for models with specific forms of functional responses and parametrization. In the framework of the traditional formulation that was used in these papers, this problem cannot be avoided.

This problem was recognized at a very early stage of existence of mathematical biology. One way of dealing with this problem was suggested in 1930s by Kolmogorov [43], who suggested to “reverse” the problem: instead of studying properties of a biological system with specifically defined functional responses and parametrization, Kolmogorov suggested to formulate a model where the functional responses are given by unspecified functions and then establish the properties (such as positivity, monotonicity or concavity/convexity, etc.) that these functions should have in order to ensure the possession of certain properties by the considered system. This Kolmogorov’s concept was successfully used in the global analysis, particularly in combination with the direct Lyapunov method [44–48]. However, it is obvious that this concept hardly can be employed in a combination with numerical analysis.

In biomedical applications of the optimal control theory this problem was recognized as well. To the best authors’ knowledge, Horst Behncke was the first who employed the Kolmogorov’s concept to an optimal control problem in mathematical biology [49]. He considered the classic -dimensional SIR Kermack–McKendrick epidemic model of an isolated epidemic assuming that the infection rate (which is the most important of all functional responses in epidemic models) is given by an unspecified function constrained by a few conditions. Behncke considered this model under action of one of three controls; namely, vaccination of the susceptible individuals, isolation of the infectious individuals, and education that reduces a probability of infection: these three controls comprise all controls possible for such a simple model. Considering arguments as above, that is taking into consideration the disparity between the costs of a control policy and the costs and losses inflicted by a full-scale epidemic, Behncke employed objective functionals independent of the controls and considered minimization of the cumulative number of the infected over a given time interval. Assuming that in each case only one control is employed and that the controls are bounded, he found qualitative forms of these optimal controls.

A similar approach was later applied to more complicated problems. Thus, the optimal controls for a SIR model of endemically persistent diseases (i.e., a SIR model with demographic processes) where the incidence rate was given by an unspecified function were considered in [33]. In this control model, four simultaneously acting controls were considered and their principal qualitative types established. It is noteworthy that there are more variants of the optimal controls for this problem compared with the single epidemic models considered by Behncke. Moreover, comparison of these results with the controls found for the same model with the standard bilinear transmission rate in [30, 32] shows that the nonlinearity of incidence rate does not affect the principal properties of the control. A more complicated -dimensional SEIR model with an unspecified nonlinear incidence rate and five simultaneously acting controls (which comprise all controls possible for this model) was considered in [34].

##### 1.5. Objectives and Mathematical Technique of this Paper

The objective of this paper is to study and establish properties of the optimal, highly active antiretroviral therapy, at the same time avoiding the above-mentioned drawbacks of the preceding studies. Having this goal in mind, we consider a model of HIV within-host dynamics based on the classical Nowak–May model, which describe interaction of the susceptible and infected target cells and free virus particles [50, 51]. The Nowak–May model postulates that the susceptible target cells can be infected by free virus particles and after the instance of infection move into the infected population; in their turn, the infected cells produce free virus particles that can infect other susceptible cells. Using the Nowak–May model as basis, we assume that the rate of infection is given by a unspecified function of the concentrations of susceptible cells and free virus particles, constrained by a few conditions. Into the model, we also add a term responsible for the loss of free virus due to infection (see [52] for discussion of this term).

In order to mirror the HAART, we introduce six bounded controls that can act simultaneously into this model. These six controls comprise all controls that are possible for this model and correspond to the six stages of the viral replication cycle that the current antiretrovirals can affect. For this model, we consider a control problem of minimizing the infection level and the viral load at the end point of a finite time interval. The objective functional that we consider is independent of the controls.

While for our type of objective functional, the optimal controls are likely to be of the bang-bang type, theoretically, they can contain singular arcs. Thus in this paper, we prove that our optimal controls are only bang-bang and that no singular arcs are possible. Next, the maximal number of switchings for each control can be analytically estimated and then the specific number of the switching moments for each control can be found numerically as a problem of the finite-dimensional optimization.

To deal with this optimal control problem, we employ a following analytical technique:(i)Using Pontryagin Maximum Principle, we found differential equations for the adjoint variables and for the switching functions that entirely determine types of the corresponding optimal controls.(ii)Then, analyzing the system of differential equations for switching functions, we prove that the optimal controls are bang-bang.(iii)Our next task is to find accurate estimates for the number of zeros of the switching functions which correspond to switching points of controls. In order to do this, we apply the generalized Rolle’s theorem [53] to the equations of the system for switching functions. However, a direct application of the theorem to these equations is impossible, and we have to convert the matrix of this system into the upper triangular form. In order to do this, we have to find a suitable transformation, which is usually defined by a system of nonautonomous quadratic differential equations, and to prove that solutions of these equations are defined on the entire time interval. One of the results of this paper is that we managed to reduce a -dimensional system of nonautonomous quadratic differential equations to two -dimensional systems: the analysis of these two -dimensional systems is simpler than that of a single -dimensional system.(iv)Having the maximal possible number of switchings found, we can immediately find the principal type of the optimal controls.(v)Finally, to solve the problem completely, one can find the specific number of switchings and their precise positions, which can be done numerically. Thus, our analytical technique allows to reduce the optimal control problem to a considerably simpler problem of finite-dimensional optimization. Despite this possibility, the obtained theoretical results are illustrated by numerical calculations using BOCOP–2.0.5 software package, and the corresponding conclusions are made.

A crucial for our objectives feature of this approach is that the estimations for the number of switchings and the values of the controls on the subintervals can be found analytically. This implies that there is a significant flexibility with respect to particular forms of the functional responses in the model and, hence, a possibility to use this approach in combination with the Kolmogorov’s “reversion of a problem” concept; we exploit this possibility in the paper.

An analytical methods used in this paper can be employed to analyze optimal controls for other possible objectives, as well as for other models of virus dynamics. These findings have the immediate practical relevance providing answers to long-existing questions about the interference of drugs with different actions combinations of which are used in HAART and about the conditions and objectives for which the current therapy guidelines are optimal.

The paper is organized as follows: In Section 2, we formulate the model and state the optimal control problem. In Section 3, we apply Pontryagin Maximum Principle to find the adjoint system and the system of nonautonomous differential equations for the switching functions. In Sections 4, we prove that the optimal controls are bang-bang. In Section 5, we find the estimates for the number of zeros of the switching functions, and thereby, in Section 6, we describe the principal types of the optimal controls that are possible for this problem. In Section 7, we provide results of computations that illustrate our analytical results. Section 8 contains conclusion and discussion of the results. Appendix A describes details of the technique of the reduction of a -dimensional matrix of a system of linear nonautonomous homogeneous differential equations for the switching functions to the upper triangular form. Appendix B is devoted to the verification of one of our auxiliary assumptions.

#### 2. Model and the Optimal Control Problem

We consider the following model of the within-host HIV dynamics:

Here , and are the appropriate concentrations of uninfected (susceptible) target cells, infected target cells and free virus particles at time . The model postulates that the susceptible cells (for HIV these are CD4+ T helper cells) enter into the system (into the patient blood from the thymus, where CD4+ cells mature) at a constant rate , that the susceptible cells are infected by the free virus particles at a rate and the infected cells immediately move to the infected population, that the infected cells produce free virus particles at the rate , and that average life-spans of the susceptible cells, infected cells and free virus particles are, respectively, , and (naturally, for HIV ). The third equation of the model also includes a removal term representing the loss of virus due to infection of cells. (This term is usually omitted in mathematical models of HIV dynamics as very small compared to the virus death rate . However, for low viral loads typical for patients under antiviral therapy, the impact of the removal of free virus due to infection can be of importance).

The model (2.1) is nearly classical. In the literature it usually appears with a bilinear incidence rate and without the term representing the loss of free virus due to infection and is referred to as the Nowak–May model [51, 54]. To the best of the authors’ knowledge, a version of the model with the unspecified general nonlinear incidence rate firstly appeared in [46]. The dynamics of this model with both bilinear and unspecified nonlinear incidence rates is studied in detail both analytically and numerically in [50, 51, 54–56, 44–47]. In particular, and what is of importance for further analysis, it is known that the nonnegative octant of the -dimensional space is a positive invariant set of system (2.1) and that all solutions of the system initiated in this set are bounded.

In this paper, we assume that the incidence rate satisfies the following conditions.

*Assumption 1. *Let the function satisfy:(i) is a continuous differentiable function for all ;(ii) for all , and for all ;(iii)Partial derivatives and are positive for all and nonnegative for all and .

Please note that this set of assumptions is less restrictive than these that are usually used in the literature (see, e.g., [44–47, 56]).

To make system (2.1) controllable, we assume that on a pre-defined time interval the system is an object of six independent controls , , , , and , which have the following biological meanings:(i) is a fraction of immune cells among the new target cells entering the system;(ii) is a per capita rate of “immunization” of the susceptible cells that makes the “immunized” cells entirely unsusceptible to the virus;(iii) is a magnitude of indirect measures that inhibit the incidence rate;(iv) is a magnitude of indirect measures aimed that inhibit the free virus particles production by the infected cells;(v) is a per capita rate of immobilization (or removal) of free virus particles from the system;(vi) is a per capita rate of killing (or of removal by other means) of the infected cells.

We would like to stress that these six controls comprise all interventions that are feasible or potentially possible for model (2.1).

We assume that the controls satisfy the following constraints:

Here the maximal values of control functions , , , , and are pre-defined positive constants.

Assuming that these six controls are applied to model (2.1) on the time interval , we obtain the following control model:

For convenience of notation, we introduce modified controls:

These modified controls satisfy the following constraints:

where

In addition to Assumption 1, we make the following assumption.

*Assumption 2. *Maximal and minimal levels , , , , , , , , of the controls satisfy the conditions:

As a set of admissible controls, , we consider the set of all Lebesque measurable functions , , , , and satisfying inequalities (2.5) for almost all .

Using modified controls (2.4), we rewrite control system (2.3) as following:

The following lemma ensures the boundedness, positiveness, and continuity of solutions for system (2.11).

Lemma 3. *For any admissible controls , , ,, , and , the corresponding absolutely continuous solution , , to system (2.11) are defined on the entire interval , and for its components , , the inequalities:**where**hold for all .*

Proof of this lemma is fairly straightforward and hence we omit it. Proofs of similar statements are given, for example, in [23, 57–59].

On the set of admissible controls for system (2.11), we consider the problem of minimizing a weighted sum of concentrations of the infected cells and free virus particles at the end of time interval :

Such objective functional usually corresponds to the task of complete elimination of the infection. In (2.15), and are positive weights, such that . Please note that the objective functional is independent of the controls.

For optimal control problem (2.11), (2.15), inequalities (2.12) in Lemma 3 and Theorem 4 of Chapter 4 in [60] ensure the existence of the optimal controls , , , , and and the corresponding optimal solution of system (2.11). In the following sections we establish the principal properties of these controls.

#### 3. Pontryagin Maximum Principle

To analyze the optimal control problem (2.11), (2.15), we use Pontryagin Maximum Principle [61]. The Hamiltonian of this problem is

where , and are adjoint variables. The Hamiltonian satisfies:

Here and are the partial derivatives of function with respect to variables and , respectively.

According to Pontryagin Maximum Principle, for the optimal controls , , , , , and the corresponding optimal solution of system (2.11), there exists a vector-function such that:(i)The vector-function is a nontrivial solution of the adjoint system

In this system, it is easy to note that

and(ii)The controls , , , , , maximize the Hamiltonian

with respect to variables:

for almost all , and, therefore, the following relationships hold:

In these relationships, by Lemma 3, functions

are the switching functions, which determine a type of the optimal controls , , , , and .

Please note that in (3.14)

Hence, it suffices to study the properties of only two of these four functions, say and . Therefore, further we focus on the analysis of the switching functions , , and . By (3.3) and (3.14), we obtain the following Cauchy problem for the switching functions , , and :

We actively use this Cauchy problem in subsequent arguments.

#### 4. Properties of the Switching Functions

We start the analysis by establishing some properties of switching functions , , , .

Lemma 4. *There exist values such that**hold.*

*Proof. *The switching functions and are the solutions of the third and fourth equations of system (3.16), respectively, and, hence, these functions are continuous functions. Since , and , by Theorem on the Stability of the Sign of a Continuous Function (see, e.g., [62]) there are intervals and where inequalities (4.1) necessary hold.

Lemma 5. *Depending on values of the weighting coefficients α and β, for switching functions and , the following statements hold:*(i)*if , then there exist , , such that*(ii)*if , then there exist , , such that*

*Proof. *One of the inequalities for function immediately follows from the arguments similar to those presented in the proof of Lemma 4. Indeed, switching function satisfies the first equation of system (3.16), and, hence, is absolutely continuous and, therefore, is also a continuous function. Then, for , by Theorem on the Stability of the Sign of a Continuous Function, there is such that, depending on the sign of , either the first inequality in (4.2), or the first inequality in (4.3), holds.

To get the inequality for , we integrate the second equation of system (3.16) on the interval . The integration yieldsThen, taking into consideration the constraints imposed on control by (2.5), Assumption 1 and the already established inequality for , we immediately see that the corresponding inequality for function in either (4.2) or (4.3) also necessary hold on the interval . From the absolute continuity of switching function and the Theorem on the Stability of the Sign of a Continuous Function it is immediately follows that there is value , , such that, depending on the sign of , the second inequality either in (4.2) or in (4.3) holds. This completes the proof.

Corollary 6. *Lemmas 4 and 5 combined with formulas (3.9), (3.10), (3.12), (3.13) and (3.15), immediately lead to the following relationships for the optimal controls , , , , and :*(i)*if , then*(ii)*if , then*

The following Lemma is of principal importance for further analysis.

Lemma 7. *The switching functions , , and are not identically equal to zero on any subinterval of interval .*

*Proof. *Firstly, we consider switching function . Let us suppose the contrary, that is that there is an interval on which

Then

almost everywhere on this interval as well and for equations of system (3.16) are

Hence,

and

where , are some constants.

Now we have to consider the possible cases depending on the particular values of and .

*Case 1*. If and , then, by (4.10) and (4.11),

everywhere on interval . Substituting (4.7) and (4.12) into (3.14), we get

everywhere on interval . System (3.3) is linear and homogeneous, and hence, equalities (4.13) necessary hold on the entire interval . This contradicts to the the boundary conditions at and, hence, this case is impossible.

*Case 2*. If and , then, from (4.10) and (4.11),

everywhere on interval . Hence, by (3.12) and (3.15), on the controls and simultaneously take either the values and , or the values and . Moreover, by the inequality from (4.14) and the last equation in (4.9), on the interval the switching function has at most one zero. (For analysis, it is sufficient that it is not equal to zero identically.) Hence, by (3.13), on . Substituting (4.14) and the above-found values of and into the first equality in (4.9), we conclude that either

or

necessary holds. However, either of these equalities contradicts to the control restrictions (2.9) and (2.10) of Assumption 2. Therefore, this case is also impossible.

*Case 3*. If and , then, by (4.10) and (4.11),

everywhere on the interval . Hence, by (3.10), for all , the control takes either value , or value . Equality in (4.17) substituted into the last equation in (4.9) yields

where is a constant. Here, the equality is impossible, since, by (3.14), in the case we come to the contradictory equalities (4.13). Therefore, , and then we have to conclude that on the interval the control takes either the value or value . Substituting (4.17) and the above-found values of control into the first equality of (4.9), we immediately conclude that either

or

must hold. However, neither of these equalities is possible as contradicting to the control restrictions (2.7) and (2.8) of Assumption 2. Therefore, this case is also impossible.

*Case 4*. If and , then, by (4.10) and (4.11),

everywhere on interval . Hence, by (3.10) and (3.12), for all , control takes either value , or value , and controls and simultaneously take either values and , or values and . By the same arguments as in Case 2, we have to conclude that control takes either value or value on the entire interval or on its subinterval .

Differentiating the first equality of (4.9) on the interval or on its subinterval , and then using the corresponding differential equations, we get equality

where , and . The first equality in (4.9) and (4.22) form a linear homogeneous system with respect to unknown and . By (4.21), the solution of this system must be nontrivial, and, hence, the determinant of this system must be equal to zero. That is, the equality

must hold. However, by the control restrictions (2.7)–(2.10) of Assumption 2, this equality is impossible. Therefore, this case is also impossible.

As all possible cases lead to a contradiction, we have to conclude that the hypothesis was incorrect and that the switching function cannot be identically equal to zero on any subinterval of the interval .

Now we have to consider whether switching function can be equal to zero on a subinterval of the interval . Let us assume that there is an interval on which

Then, almost everywhere on this interval

holds as well. Substituting (4.24) and (4.25) into the second equation of system (3.16) and using the constraints on control from (2.5) and the properties of function from Assumption 1, we have to conclude that the contradictory equality holds for all . Hence, our assumption is incorrect and the switching function cannot be identically equal to zero on any subinterval of the interval .

Next, we consider a possibility that the switching function is identically equal to zero on any subinterval of the interval . Applying the same arguments as those for the function , but using the third equation of system (3.16) instead of the second equation, we come to the conclusion that the switching function cannot be identically equal to zero on any subinterval as well.

Finally, we consider the switching function . Let us suppose the contrary and assume that there is the interval on which

Then, almost everywhere on this interval

holds. Again, we substitute (4.26) and (4.27) into the fourth equation of system (3.16). Using the constraints on control from (2.5) we find the contradictory equality for . Hence the above-made assumption was incorrect and the switching function is not identically equal to zero on any subinterval of the interval . This completes the proof.

Lemma 7 combined with equalities (3.8)–(3.13) and (3.15) immediately leads to the following Corollary.

Corollary 8. *The optimal controls , , , , and are the bang-bang functions taking values , , , , and , respectively. Moreover, each of the pairs composed of controls and , and and switches from the minimum values to the maximum values and vice versa at the same switching moments, which are specific for each of these pairs.*

In addition to this corollary, we note that the controls , , , , and have no singular arcs (see [23, 24] for details).

#### 5. Estimating the Number of Zeros of the Switching Functions

In previous section, we proved that the optimal controls , , , , and are the bang-bang functions with values , , , , and , respectively. In this section, we proceed to estimating the number of switchings of these controls. Equalities (3.8)–(3.13) and (3.15) show that for this task it suffices to estimate the number of zeros of the switching functions , , and defined by system (3.16). From the analysis of the second and third equations of this system, we conclude that if we have an estimate for the number of zeros of function , then, applying the Generalized Rolle’s Theorem (see, e.g., [53]) to each of these equations, we can find estimates for the number of zeros of functions and as well. The use of this theorem with the fourth equation of system (3.16) gives us an estimate of the number of zeros of function . Thus, the analysis of switching function is the basis for completing our task. The method that allows to estimate the number of zeros of the function consists in the subsequent application of the above-mentioned Generalized Rolle’s Theorem to the equations of system (3.16).

Let us consider a possibility of application of the Generalized Rolle’s Theorem to the first equation of system (3.16). The direct application of the theorem to this equation is impossible, because, apart of the term containing the function , the equation has two terms independent of . Therefore, we introduce an auxiliary function

which, by Lemma 7 and relationships (2.7)–(2.10) of Assumption 2, is not equal to zero on any subinterval of . Then, the first equation of system (3.16) can be rewritten as

In order to formulate the differential equation for function , we note that the switching functions and are absolutely continuous (as the solutions of system (3.16) and, therefore, are almost everywhere differentiable functions on the interval [63]. Theorem 3 of Chapter 4 in [60], gives us absolute continuity and hence, also almost everywhere differentiability of functions and . Finally, arguments similar to these used by the authors in [34] allow to conclude that the functions and are also differentiable almost everywhere on the interval . Then, using the second and third equations of system (3.16), we formulate the required differential equation:

Let us introduce another auxiliary function

Then, equation (5.3) can be rewritten as follows:

By the above-presented arguments, function is also almost everywhere differentiable and satisfies

The last two terms of this equation can be expressed in terms of the auxiliary functions and . Then,

Let us denote

Due to the constraints on the controls , , and in (2.5) and the properties of function in Assumption 1, we conclude that functions , and are positive almost everywhere on the interval . Moreover, functions and are piecewise-constant and functions , and are piecewise-continuous. Combining equations (5.2), (5.5) and (5.7), we write down the following system of differential equations for the switching function and auxiliary functions and , which is more convenient for the subsequent analysis than system (3.16):

Let us rewrite system (5.9) in the following matrix form:

Now, using a special change of variables, we reduce the matrix of system (5.10) to the upper-triangular form. This transformation is fairly standard, and is described in detail in Appendix A. Comparing systems (5.10) and (A.1) and applying arguments from Appendix A, we find the desired system in the form

where is the switching function, and , are the new auxiliary functions corresponding to . According to the arguments of Appendix A, the functions , , and perform such a transformation and, by (A.21) and (A.26), satisfy the following systems of differential equations: