Chattering and its approximation in control of psoriasis treatment

We consider a nonlinear system of differential equations describing a process of the psoriasis treatment. Its phase variables are the concentrations of T-lymphocytes, keratinocytes and dendritic cells. A scalar bounded control is introduced into this system to reflect the medication dosage. For such a control system, on a given time interval the minimization problem of the Bolza type functional is stated. Its terminal term is the concentration of keratinocytes at the terminal time, and its integral term is the product of the non-negative weighting coefficient with the total cost of the psoriasis treatment. This cost is linear in the control and proportional to the concentration of keratinocytes. For the analysis of such a problem, the Pontryagin maximum principle is used. As a result of this analysis, it is shown that if the weighting coefficient is zero, then the corresponding optimal control can contain a singular arc. We establish that it is a chattering control, and therefore does not make much sense as a type of a medical treatment. If the weighting coefficient is positive, then the corresponding optimal control is bang-bang, and it can be presented as a type of psoriasis treatment. In addition, when this coefficient tends to zero, such optimal controls can be considered as chattering approximations. Therefore, the convergence of these optimal controls, the corresponding optimal solutions of the original system, and the minimum values of the functional are studied. The obtained theoretical results are illustrated by numerical calculations and the corresponding conclusions are made.

Abstract. We consider a nonlinear system of differential equations describing a process of the psoriasis treatment. Its phase variables are the concentrations of T-lymphocytes, keratinocytes and dendritic cells. A scalar bounded control is introduced into this system to reflect the medication dosage. For such a control system, on a given time interval the minimization problem of the Bolza type functional is stated. Its terminal term is the concentration of keratinocytes at the terminal time, and its integral term is the product of the non-negative weighting coefficient with the total cost of the psoriasis treatment. This cost is linear in the control and proportional to the concentration of keratinocytes. For the analysis of such a problem, the Pontryagin maximum principle is used. As a result of this analysis, it is shown that if the weighting coefficient is zero, then the corresponding optimal control can contain a singular arc. We establish that it is a chattering control, and therefore does not make much sense as a type of a medical treatment. If the weighting coefficient is positive, then the corresponding optimal control is bang-bang, and it can be presented as a type of psoriasis treatment. In addition, when this coefficient tends to zero, such optimal controls can be considered as chattering approximations. Therefore, the convergence of these optimal controls, the corresponding optimal solutions of the original system, and the minimum values of the functional are studied. The obtained theoretical results are illustrated by numerical calculations and the corresponding conclusions are made.
1. Introduction. The study of psoriasis is increasingly important as the disease leads to the growth in morbidity, disability, and other debilitating circumstances. Through the efforts of scientists a certain success was achieved in revealing the mechanisms of the development of this disease. Psoriasis is one of the most common dermatoses. Its frequency varies from 0.1 to 3% in different countries. The incidence of psoriasis in the USA is from 0.5 to 1.5%. On a share of psoriasis among persons addressed with a dermal pathology it is necessary from 2 to 8%. The causes of psoriasis remain unclear, although in recent years significant progress in studying the mechanisms of pathogenesis of the disease has been achieved. Despite the fact determined the ways of possible concatenations of a singular arc with non-singular portions, we can find the specific optimal solution of the optimal control problem only numerically ( [38]).
The results obtained in [4,7,33] are very interesting, the presence of a square of the control under the integral of the objective functional, as it was noted above, makes the optimal control problem numerically simplistic but prohibits the attainment of the optimal analytical solution. Before solving an optimal control problem numerically, it would be beneficial to investigate the problem analytically in order to reveal some features and properties of its optimal solution. Therefore, in this paper we will consider the control model proposed in [33] with a different objective functional. Namely, one of two problems that we will study is the problem of minimizing the concentration of keratinocytes at the end of the time interval. It will be shown that the corresponding optimal control may have a singular arc of the order of two, that is, such a control is chattering. It is important to note that if as the objective functional, we consider the cumulative concentration of keratinocytes (remove the total cost of treatment from the functional considered in [33]), then a preliminary analysis shows that the corresponding optimal control has the same feature, only established as a result of more complicated and cumbersome calculations.
Let us now consider for the control model of the psoriasis treatment the problem of minimizing the objective functional in the form of the weighted sum of the concentration of keratinocytes at the end of the time interval and the total cost of treatment that is the integral of control. A preliminary analysis shows as well that the corresponding optimal control may have a singular arc of the order of one. At the same time, it is very difficult to verify for it the fulfillment of the necessary optimality condition ( [37]). Therefore, the second problem that we will study is the problem of minimizing the similar functional, only with the total cost of treatment that is linear in the control and proportional to the concentration of keratinocytes.
This paper is organized as follows. In Section 2, within the framework of the mathematical description of the process of the psoriasis treatment, on a given time interval a nonlinear control system of three differential equations describing relationships between the concentrations of T-lymphocytes, keratinocytes and dendritic cells is considered. The scalar bounded control is also included in this system to reflect medication intake. Such properties of its solutions as positiveness, boundedness, and continuation on a given time interval, are discussed. Then, the minimization problem of the Bolza type functional is stated. Its terminal term is the concentration of keratinocytes at the terminal time, and its integral term is the product of the non-negative weighting coefficient with the total cost of the psoriasis treatment. This cost is linear in the control and proportional to the concentration of keratinocytes. The existence in such a problem of an optimal solution consisting of the optimal control and corresponding optimal solutions of the system of differential equations, is established as well. For its analysis, in Section 3, the Pontryagin maximum principle is used. In Section 4, using the Lie brackets from the geometric control theory, it is shown that if the weighting coefficient is zero, then the corresponding optimal control can contain a singular arc. The order of such a singular arc is determined, the fulfillment of the necessary optimality condition for it is discussed, and the forms of a concatenation of a singular arc and non-singular bang-bang portions are found. As a result, it is established that the optimal control is chattering. Therefore, it does not make much sense as a type of a medical treatment for psoriasis. The corresponding results of numerical calculations and related conclusions are presented. If the weighting coefficient is positive, then, as reflected in Section 5, the corresponding optimal control is bang-bang, and it can be presented as a type of psoriasis treatment. The corresponding numerical calculations and related conclusions are made as well. In addition, when this coefficient tends to zero, such optimal controls can be considered as chattering approximations. Therefore, the convergence of these optimal controls, the corresponding optimal solutions of the original system, and the minimum values of the functional are studied. Finally, Section 6 contains our final conclusions. 2. Mathematical model and its optimal control problem. Let us consider on a given time interval [0, T ] a mathematical model describing the relationships between the concentrations l(t), k(t), and m(t) of T-lymphocytes, keratinocytes and dendritic cells, respectively, at a specific time t. The interaction between these types of cells causes the occurrence of psoriasis. We suppose that σ is the constant rate of accumulation for T-lymphocytes and ρ is the constant accumulation rate of dendritic cells; the rate of activation of T-lymphocytes by dendritic cells is δ and also β is the activation rate of dendritic cells by T-lymphocytes as well. Moreover, the per capita removal rates of T-lymphocytes, epidermal keratinocytes and dendritic cells are denoted by µ, λ, and ν, respectively. We consider that γ 1 is the rate of activation of keratinocytes due to T-lymphocytes and growth of keratinocytes by T-lymphocytes occurs at the rate γ 2 . Finally, we introduce here the medication dosage in term of a control function u(t) at the places of the interaction between T-lymphocytes and keratinocytes to restrict the excess growth of recent ones. Assembling all the above assumptions, we formulate the mathematical control model of psoriasis treatment as follows: with given initial conditions: In system (1) the control u(t) subjects to the following restrictions: The set of admissible controls U (T ) is formed by all Lebesgue measurable functions u(t), which for almost all t ∈ [0, T ] satisfy inequalities (3). Moreover, we identify functions u(·) ∈ U (T ) that coincide on the interval [0, T ] almost everywhere. Now, let us define the following positive constants: and introduce a region: It is easy to see that the inclusion holds. Here R 3 is the Euclidean space consisting of all column vectors and the sign means transposition.
The boundedness, positiveness and continuation of the solutions for system (1) is stated by the following lemma.
Lemma 2.1. For any admissible control u(t) the corresponding absolutely continuous solutions l(t), k(t), m(t) for system (1) are defined on the entire interval [0, T ] and satisfy the inclusion: Proof. Let us consider an arbitrary admissible control u(t) and the corresponding to it absolutely continuous solutions l(t), k(t), m(t) for system (1). We suppose that these solutions are defined on the interval [0, t 0 ) that is the maximum possible interval for the existence of these solutions. Writing the first and third equations of system (1) in the form of linear nonautonomous differential equations with positive inhomogeneities: we easily see that their integration with the corresponding positive initial conditions from (2) implies the inequalities: The second equation of system (1) written in the form: is also a linear non-autonomous differential equation with positive inhomogeneity due to inequalities (6). Its integration with the corresponding positive initial condition from (2) leads to the inequality: Hence, the positivity of the solutions l(t), k(t), m(t) on the interval [0, t 0 ) is proved. Now, we show the boundedness of these solutions. For this, the function W (l, k, m) = γ 2 l + γ 1 k + 1 + β −1 δ γ 1 m is introduced, and then for the solutions l(t), k(t), m(t) its derivative by virtue of system (1) is calculated: By inequalities (6), this formula implies the inequality: or, after standard arguments, the inequality: Integrating it on the interval [0, t], we find a chain of the inequalities: The definition of the function W (l, k, m) and inclusion (4) imply the inequality: Thus, we have established the boundedness of the solutions l(t), k(t), m(t) on the interval [0, t 0 ). Finally, when T < t 0 , the fact justified in this lemma is proved. For T ≥ t 0 , the existence of these solutions on the whole interval [0, T ] follows from their continuation due to inequalities (6)-(8) and Theorem 3.1 ( [16], Chapter 2). Hence, inclusion (5) is established and the proof of the lemma is complete. Now, let us consider for system (1) on the set of admissible controls U (T ) the following problem of minimization: where ξ is a non-negative weighting coefficient. When ξ = 0, problem (9) is to minimize the concentration of keratinocytes at the terminal time T : For ξ > 0, this problem implies minimizing the weighted sum of the already mentioned concentration of keratinocytes at the end of time interval [0, T ] and the total cost of psoriasis treatment. This cost is linear in the control u(t) and proportional to the concentration of keratinocytes k(t). The bigger psoriatic lesions of human skin correspond to higher concentration of keratinocytes, and hence to the higher cost of treatment. We note that the total cost of disease control constraints as the additional term of the type: were used previously in [5,11,12,13] for SIR and SEIR control models. Also, we see that the optimal control problem (9) differs from problems that are typically considered in the literature on the control of psoriasis models ( [4,7,33]) in that the functional of (9) does not include an integral of the square of the control u(t), which is responsible for the total cost of psoriasis treatment. The existence in problem (9) of the optimal control u * ξ (t) and the corresponding optimal solutions l * ξ (t), k * ξ (t), m * ξ (t) for system (1) follows from Lemma 2.1 and Theorem 4 ([20], Chapter 4).
Corollary 3.2. Formula (14) and Lemma 3.1 yield the following relationship for the optimal control u * 0 (t) :  (14) shows the possible types of the optimal control u * ξ (t). It can have a bang-bang type, and switches only between the values u min and 1. Such situation occurs, when at all time moments where the switching function L ξ (t) is zero, its derivativeL ξ (t) does not vanish. Or, in addition to the portions of bang-bang type, the control u * ξ (t) can also contain a singular arc ( [37,43]). This situation arises, when the switching function L ξ (t) vanishes identically over some open subinterval of the interval [0, T ]. The next two sections are devoted to finding the relationships between the parameters of system (1) and values of the weighting coefficient ξ under which the optimal control u * ξ (t) is of one type or another.
4. Investigation of a singular arc for ξ = 0. First, let us consider the case ξ = 0, and hence study the existence of a singular arc for optimal control u * 0 (t). According to [37,43], this means that can the corresponding switching function L 0 (t) become zero identically on some subinterval of the interval [0, T ]? To answer this question, we use the concepts and notations of the geometric control theory from [37,38].
We consider the Euclidean space R 3 consisting of all possible column vectors, in which r, s is the scalar product of its elements r, s ∈ R 3 . Let z = (l, k, m) ∈ R 3 . We rewrite system (1) in the form: where with f (z) the drift and g(z) the control vector fields of this system. Here z(t) is the column vector consisting of the solutions l(t), k(t), m(t) that correspond to the admissible control u(t), that is z(t) = (l(t), k(t), m(t)) ∈ R 3 . Let Df (z) and Dg(z) be the Jacobian matrices of the vector functions f (z) and g(z), respectively. By formulas (17), we find the relationships: Then, using the introduced concepts and notations, it is easy to see that z * is the optimal trajectory for system (16) corresponding to the optimal control u * 0 (t). Also, let ψ * 0 (t) = (ψ * 1,0 (t), ψ * 2,0 (t), ψ * 3,0 (t)) be the appropriate non-trivial solution for the adjoint system (13), or in the new notations: . The Hamiltonian (15) takes the form: . The switching function L 0 (t) becomes the scalar product of the adjoint solution ψ * 0 (t) with the control field g(z * 0 (t)): Next, according to [37,38], we introduce for the drift and control vector fields f (z) and g(z) the corresponding Lie brackets: where D[f, g](z), by analogy with formulas (18), is the Jacobian matrix of the vector function [f, g](z) defined by equality (20).
For this, first, we introduce a region Λ 0 in which, as follows from Lemma 2.1, the dynamics of system (16) changes: Next, in addition to the vector fields f (z) and g(z), we define the following two auxiliary vector fields: Calculating the determinant of the vectors g(z), h(z), and q(z), it is easy to find that it is not equal to zero for all z ∈ Λ 0 . Hence, the following lemma is true.
Lemma 4.1. In the region Λ 0 , the vector fields g(z), h(z), and q(z) are linearly independent.
Moreover, we introduce the auxiliary vector functions of variable z: Now, using formulas (17) and (18), by direct computations we find the Lie bracket (20): ELLINA GRIGORIEVA AND EVGENII KHAILOV which then can be rewritten in the more convenient form: where p f 11 (z) is the scalar function of the vector argument z: (25), we find the Jacobian matrix of the vector function [f, g](z) as where ∇p f 11 (z) is the column gradient of the function p f 11 (z). Substituting formulas (25) and (26) into relationships (21) and (22), after the necessary transformations, the following formulas are obtained: For this, we note that direct calculations allow us to find the convenient representations for the Jacobian matrices of the vector functions g(z), h(z), and f (z): Substituting formulas (29) into the expressions of the Lie brackets [g, h](z) and [f, h](z), similar to (20), after the necessary transformations, the following relationships can be obtain: Let us substitute formula (30) into relationship (27), and then perform the necessary transformations. As a result, we find the following representation of the Lie bracket [g, [f, g]](z): where p g 21 (z) and p g 22 (z) are the scalar functions of the vector argument z: . Now, we substitute the formulas of the vector functions c(z), h(z), f (z), and χ(z) into the corresponding terms of the last brackets of formula (31), and then transform the resulting expression. After that, we substitute such a modified formula (31) into relationship (28), and then perform the necessary transformations as well. As a result, we obtain the following representation of the Lie bracket [f, [f, g]](z): where p f 21 (z) and p f 22 (z) are the scalar functions of the vector argument z: and d(m) is the scalar function of the argument m: Finally, let us substitute relationships (25), (32), and (33) in the corresponding formulas (23) and (24) Now, let us return to the question posed at the beginning of this section. Let the switching function L 0 (t) vanish identically on a some interval ∆ ⊂ [0, T ]. This implies, on the one hand, by formula (19), the fulfillment of the equality: and, on the other hand, a similar vanishing on the interval ∆ of the first and second derivatives L 0 (t) of the switching function L 0 (t). As a consequence of this fact and formulas (34) and (35), the following equalities are true: By Lemma 4.1 and the non-triviality of the adjoint solution ψ * 0 (t), relationships (36) and (37) imply the equality: The formula of the function d(m) shows that the analysis of equality (38) is related to the behavior of the quadratic function: the discriminant of which has the form: Using inequalities (12), we see that depending on the values α and D w the following cases are possible.
Case A. If α < 0 and D w < 0, then the function w(m) takes negative values for all m > 0, and hence equality (38) is impossible. Thus, the switching function L 0 (t) cannot be zero on any subinterval of the interval [0, T ]. Therefore, the optimal control u * 0 (t) does not have a singular arc. By formula (14), we see that it is bangbang control taking the values {u min ; 1} on the interval [0, T ]. Here, we can consider the question of estimating the number of zeros of the function L 0 (t), and hence the estimate of the number of switchings of the control u * 0 (t). This question was studied in more detail in [14].
For simplification of consequent arguments let the following condition be valid.

Condition 4.2.
We exclude from consideration the case when α < 0 the equality D w = 0 holds.
Next, as we see from formulas (35), (36), and also the equality: following from the first relationship of (37), the second derivative L 0 (t) of the switching function L 0 (t) does not contain the control u * 0 (t). This means that in Cases B and C the order q sing of the singular arc is greater than one ( [37,43]). Therefore, we take the third derivative L It is easy to see that this derivative is defined by the Lie bracket [f, [f, [f, g]]](z).
As in the case of finding the first two derivatives L 0 (t), we find the representation of this bracket in terms of vector fields g(z), h(z), and q(z). To do this, by analogy with formulas (20)-(22), we write the corresponding relationship: By formula (33), the Jacobian matrix of the vector function [f, [f, g]](z) has the form: Substituting relationships (33) and (43) into formula (42) and performing the necessary transformations, we find the formula: In this expression, the Lie brackets [f, g](z) and [f, h](z) are already found and given by formulas (25) and (31). Let us compute the Lie bracket [f, q](z). Substituting the third formula of (29) and the easily verifiable formula: into the formula of the Lie bracket [f, q](z), similar to (20), after the necessary transformations, the following relationship is obtained: Next, we substitute formulas (25), (31) and (46) where p f 31 (z), p f 32 (z), and p f 33 (z) are the scalar functions of the vector argument z: Now, let us substitute the relationship (47) into formula (41) of the third derivative L that is valid everywhere on the interval ∆. Then, this equality to zero of L 0 (t) leads to the expression: Condition 4.2 and equality (38) imply the inequality: Finally, by the formula of p f 33 (z), the first formula of (17), and also formula (38), equality (49) is equivalent to the relationship: , which allows us to find the value l sing for the variable l on the singular arc: Here, m sing is the value for the variable m on the singular arc as well. For Case B, this is m 0 sing , and for Case C, one of the values m 1 sing , m 2 sing . In addition, the positivity of the value l sing leads to the restriction on the value m sing :

ELLINA GRIGORIEVA AND EVGENII KHAILOV
Now, let us evaluate on the interval ∆ the fourth derivative L (4) 0 (t) of the switching function L 0 (t). Using [37,43], we see that it is written as follows: . We will argue, as in the case of finding the first three derivatives of the function L 0 (t). Namely, we obtain a representation of these brackets through vector fields g(z), h(z), and q(z). To do this, by analogy with formulas (20)- (22) and (42), we write the corresponding relationships: By formula (47) Substituting relationships (47) and (56) into formulas (54) and (55), after performing the necessary transformations, we find the formulas: In formulas (57) (46), respectively. Let us compute the Lie bracket [g, q](z). Substituting the first formula of (29) and formula (45) into the expression for this Lie bracket, which is similar to (20), after the necessary transformations the following relationship can be obtained: [g, q](z) = − a(z), q(z) g(z) + a(z), g(z) q(z).
We substitute formulas (25) where p f 4i (z) and p g 4i (z), i = 1, 2, 3, are the scalar functions of the vector argument z: . Now, we substitute formulas (60) and (61) into formula (53) of the fourth derivative L (4) 0 (t) of the switching function L 0 (t), and then equate it to zero. Using relationships (36), (38), (40), and (49), as a result, we obtain the equality: which is valid everywhere on the interval ∆. Let us write the equality (62) in a more convenient form for subsequent analysis. For this, first, by the formulas of the functions p f 11 (z) and p f 22 (z), we transform the formula of the function p f 33 (z) to the form: Next, we calculate the column gradient of this function. After that, together with formulas (17) of the vector fields f (z) and g(z), we substitute it into equality (62), and use as well equalities (38), (49), and also the formula: Here, the function w(m), which plays a key role in the analysis of the singular arc, is defined by formula (39). As a result, we obtain the equality: where u sing (t) is control and k sing (t) is the value of variable k on the singular arc. From analysis of this formula we can make the following conclusions. Firstly, the order q sing of the singular arc equals two. Secondly, the necessary condition of the optimality of the singular arc (Kelly-Cope-Moyer condition [43]) has the type: w (m sing ) ψ * 0 (t), q(z * 0 (t)) ≤ 0, t ∈ ∆. By inequalities (48), (50), and formula (63), the Kelly-Cope-Moyer condition is either valid in a strengthened form: or it is not valid at all, i.e.
Formula (51) allows us to rewrite relationship (67) as follows: Here r(m) is a quadratic function given by the formula: By inclusion (52), we consider this function on the interval [0, ν −1 ρ]. The relationships: r(0) = −µρ < 0, r ν −1 ρ = ν −1 ρσβ > 0 lead to the conclusion that the function r(m) has exactly one root m 0 on the interval 0, ν −1 ρ . In turn, this fact implies the validity of the relationship: which we use in analysis of formula (68). Indeed, the positivity of the left hand side of this formula and inclusion (52) imply the inclusion: which is the necessary condition for the existence of the singular arc. Let us analyze inclusion (69) for Cases B and C considered earlier.
Case B. If α > 0, then one of the following situations is possible: • the value m 0 sing / ∈ m 0 , ν −1 ρ . Therefore, the singular arc is not possible. Optimal w (m 1 sing ) > 0, w (m 2 sing ) < 0, and the function ψ * 0 (t), q(z * 0 (t)) , as it follows from (48), is sign-definite on the interval ∆, then for one of them inequality (65) will be satisfied, and for the other inequality (66). It means that only one of them satisfies the necessary condition of the optimality of the singular arc (Kelly-Cope-Moyer condition). Namely, this value becomes the value m sing of variable m on the interval ∆.
Finally, we discuss another important question related to the behavior of optimal control u * 0 (t) on the interval [0, T ]. When u sing (t) ∈ [u min , 1], the function u * 0 (t) = u sing (t) is the admissible control. Corollary 3.2 shows that the singular arc of the optimal control u * 0 (t) is concatenated with non-singular portion, where this control is bang-bang. Let τ ∈ (0, T ) be the time moment, where such portions are concatenated. Then, as it follows from [37,43], when u sing (t) ∈ (u min , 1) for all t ∈ ∆, the non-singular portion contains at least the countable number of switchings of the control u * 0 (t), accumulating to the time moment τ . This behavior of the optimal control u * 0 (t) on non-singular portions is called a chattering ( [37,43]), and will be observed on both sides of the interval ∆.
Thus, the following lemma is valid.
For case ξ = 0, the optimal control u * 0 (t) can contain a singular arc of order two, which concatenates with bang-bang portions of this control using chattering.

5.
Numerical results for ξ = 0. To illustrate a chattering at the optimal control u * 0 (t), we demonstrate the results of a numerical solution of the minimization problem (10) corresponding to the considered case ξ = 0 for Case B. As it was shown earlier, a chattering is possible for Case C as well. The corresponding numerical calculations were also made, but there were no qualitative differences from results obtained for Case B. Therefore, here we do not give them.
In Case B, when α > 0, for numerical calculations the following values of the parameters of system (1), initial conditions (2), and control constraints (3) taken from [33] were used: These numerical calculations were conducted using "BOCOP-2.0.5" ( [2]). It is an optimal control interface, implemented in MATLAB, for solving optimal control problems with general path and boundary constraints and free or fixed final time. By a time discretization, such problems are approximated by finite-dimensional optimization problems, which are then solved by well-known software IPOPT, using sparse exact derivatives computed by ADOL-C. IPOPT is the open-source software package for large-scale nonlinear optimization.
Considering the time interval of 100 days (T = 100.0), a time grid with 5000 nodes was created, i.e. for t ∈ [0, 100.0] we get t = 0.02. Since our problem is solved by a direct method and, consequently, using an iterative approach, we impose at each step the acceptable convergence tolerance of ε rel = 10 −14 . Moreover, we use the sixth-order Lobatto III C discretization rule. In this respect, for more details we refer to [2].
The corresponding results of the numerical calculations are presented in Figure 1. We give the graphics of the optimal control u * 0 (t) and the corresponding optimal solutions l * 0 (t), k * 0 (t), and m * 0 (t). J * 0 = 2.860143 is the minimum value of the functional J 0 (u) of (10). IPOPT used 50 iterations in order to solve the minimization problem (10).
It is important to note that the control u(t) is auxiliary. It is introduced into system (1) to simplify analytical analysis. The corresponding physical control v(t) in the same system is related to the control u(t) by the formula v(t) = 1 − u(t). Therefore, where the auxiliary optimal control u * 0 (t) has a maximum value of 1, the 2268 ELLINA GRIGORIEVA AND EVGENII KHAILOV appropriate physical optimal control v * 0 (t) takes a minimum value of 0, and vice versa.
The conducted numerical calculations shown the performance of the software "BOCOP-2.0.5" in the study of such a complex from the computational point of view of the phenomenon as chattering. Namely, Figure 1 shows that the physical optimal control v * 0 (t) has a period of psoriasis treatment with a smooth increase in the dose of used medication. To the beginning and the end of this period of psoriasis treatment there are the periods with increasing number of switchings with lower intensity to greatest intensity and vice versa, which does not make much sense as a type of medical treatment. At the same time, the whole process of this treatment ends with the period with the greatest intensity. Moreover, Figure 1 shows that the optimal concentration of keratinocytes k * 0 (t) decreases to an acceptable level to the end T of psoriasis treatment.  6. Investigation of a singular arc for ξ > 0. Now, let us consider the case ξ > 0. Again we study the existence of a singular arc for optimal control u * ξ (t), which implies the study of the possibility of the identity vanishing of the switching function L ξ (t) on some subinterval of the interval [0, T ]. Here we also use the concepts and notations of geometric control theory ( [37,38]), widely used in Section 4. Namely, we fix system (1) written in form (16) with the drift and control vector fields f (z) and g(z) given by formulas (17). Their Jacobian matrices are expressed by formulas (18). Again we consider that z(t) is the column vector composed of the solutions l(t), k(t), m(t) corresponding to the admissible control u(t), that is z(t) = (l(t), k(t), m(t)) . Then, z * ξ (t) = (l * ξ (t), k * ξ (t), m * ξ (t)) is the optimal trajectory for system (16) corresponding to the optimal control u * ξ (t); function ψ * ξ (t) = (ψ * 1,ξ (t), ψ * 2,ξ (t), ψ * 3,ξ (t)) is the appropriate non-trivial solution of the the adjoint system (13), which in the new notations has the form: where e = (0, 1, 0) . The Hamiltonian (15) takes the form: . The switching function L ξ (t) is rewritten as follows: The study that we carry out in this section is related to the investigation of the first and second derivatives L  ξ (t) of the switching function L ξ (t). Therefore, using system (16) written for the trajectory z * ξ (t), system (71), as well as the Lie brackets given by relationships (20)-(22), we find the corresponding formulas for these derivatives: Here ψ = (ψ 1 , ψ 2 , ψ 3 ) is column vector of adjoint variables. Let us transform the function Φ ξ (z, ψ) as follows. We substitute into it the value Df (z)g(z) from (20), and also formulas (25) and (32), and take into account the equality e, h(z) = 0. Then, for Φ ξ (z, ψ) the following formula can be obtained: where p f 11 (z), p g 21 (z), and p g 22 (z) are the scalar functions of the vector argument z, defined in Section 4. Now, we back to the purpose of our study. Let the switching function L ξ (t) vanish identically on a some interval ∆ ⊂ [0, T ]. This means, on the one hand, by formula (72), that the following equality holds: and, on the other hand, the analogous vanishing on the interval ∆ of the first and second derivatives L ξ (t) of the switching function L ξ (t). As a consequence of this fact and formulas (73) and (74), the following equalities are true: ψ * ξ (t), [f, g](z * ξ (t)) + ξ e, f (z * ξ (t)) + g(z * ξ (t)) = 0, We note that in the last relationship the function Φ ξ (z, ψ) is defined by formula (75). Now, it is important to establish whether this function is zero on the interval ∆, or not. Therefore, we will continue studying it. Let us transform equality (77) as follows. We substitute into it formula (25). After that, the value ψ * ξ (t), h(z * ξ (t)) can be expressed from the resulting expression. By equality (76), we find the relationship: Finally, let us substitute formulas (76) and (79) into formula (75) of Φ ξ (z * ξ (t), ψ * ξ (t)). This means the consideration of such a function directly on the interval ∆. After the necessary transformations using the first formula of (29), as well as the formula of the function p g 22 (z), we obtain the relationship: Direct calculations in the right hand side of this expression with the substitution of the formulas of the functions p f 11 (z), p g 21 (z), p g 22 (z), formulas (17) and again the first formula of (29), lead to the final relationship: Analysis of expression (78) with the use of (80) shows that on the interval ∆ the necessary optimality condition for a singular arc (Kelly condition) is not valid (see [43]). This means that the optimal control u * ξ (t) is bang-bang taking only values {u min ; 1}.
This means that the following lemma is true.
Lemma 6.1. For case ξ > 0, the optimal control u * ξ (t) does not contain a singular arc. It is a bang-bang function with values {u min ; 1}.
If we consider the problem of maximization of the functional J ξ (u) of (9), then repeating similar arguments and calculations we can see that the transformations of the function Φ ξ (z * ξ (t), ψ * ξ (t)) would change but not dramatically. As a result of such transformations, in the formula of the function Φ ξ (z * ξ (t), ψ * ξ (t)), similar to (80), disappears minus sign. Hence, in the case of possible existence of a singular arc the corresponding necessary optimality condition will be valid in the strengthened form: which means the concatenation of a singular arc with non-singular bang-bang portions of the corresponding optimal control u * ξ (t). The existing situation in problem (9) with the behavior of the optimal control u * ξ (t) is called elliptic ( [3,38]), and it does not guarantee any estimate of the number of switchings of the control u * ξ (t). Although, some useful considerations on this subject are present in [38]. In addition, analyzing the formula for the switching function L ξ (t), we can not determine the interval with which a value of the control u * ξ (t) (with u min or 1), will be adjacent to the end T (see for comparing Corollary from Lemma 3.1 for case ξ = 0).

7.
Numerical results for ξ > 0. The considered case (ξ > 0) we support by the conducted numerical calculations using as well "BOCOP-2.0.5" for the values of the parameters of system (1), initial conditions (2), and control constraints (3) of (70). The weighting coefficient ξ in functional J ξ (u) of (9) takes values {0.001; 0.0001; 0.000001}. The numerical results are demonstrated in Figures 2-4. Here J * ξ is the minimum value of the functional J ξ (u). Figures 2-4 show that the physical optimal control v * ξ (t) has a large initial period of psoriasis treatment is performed with the lower intensity. Then, the switching occurs, and psoriasis treatment is provided with the greatest intensity, but in a smaller period. Then again, after the switching, this treatment is performed with the lower intensity, but in an even shorter period. Finally, the switching occurs, and psoriasis treatment is returned to the greatest intensity. In addition, these figures show that the optimal concentration of keratinocytes k * ξ (t) decreases to an acceptable level to the end T of psoriasis treatment.   M ξ specifies the number of iterations required by IPOPT in order to solve the minimization problem (10). Table 1. Minimum values J * ξ of the functional J ξ (u) and the corresponding numbers of iterations M ξ .
The analysis of the numerical calculations for all values of ξ leads to the following conclusions. First, the minimum value J * ξ of the functional J ξ (u) decreases monotonically from 3.030960 at ξ = 10.0 up to 2.860158 at ξ = 0.000001. Let us note that for ξ = 0 the corresponding minimum value J * 0 is 2.860143. Secondly, optimal solutions l * ξ (t), k * ξ (t), m * ξ (t) do not change qualitatively with decreasing the value ξ. Thirdly, the optimal control u * ξ (t) with decreasing the value ξ varies  from a constant function that takes the value of 1 for ξ ∈ {10.0; 1.0} to a piecewise constant function that takes the values {u min ; 1} for ξ ∈ {0.001; 0.0001; 0.000001} (see Figures 2-4); switching appears only at ξ = 0.1. Corresponding conclusion can be made for the physical optimal control v * ξ (t). The appearance of switchings with decreasing the value ξ is due to the fact that, as it shown in Section 4, in the minimization problem (10) the optimal control u * 0 (t) contains the singular arc to which non-singular bang-bang portions with increasing number of switchings are concatenated. As it was already noted in this section, this optimal control does not make much sense as a type of medical treatment. Therefore, the optimal controls u * ξ (t) for ξ ∈ {0.001; 0.0001; 0.000001} can be considered as approximating controls, and the minimization problem (9) for ξ > 0 as an approximating problem for ξ = 0 (see [41]). With regard to this, it remains only to investigate when the weighting coefficient ξ tends to zero, the convergence of the optimal controls u * ξ (t) to the optimal control u * 0 (t), the optimal trajectories z * ξ (t) to the optimal trajectory z * 0 (t), and, finally, the minimum values J * ξ of the functional J ξ (u) to the minimum value J * 0 of the functional J 0 (u).

8.
Behavior of the optimal solutions as ξ → 0. Here we will study the convergence of the optimal controls u * ξ (t) to the optimal control u * 0 (t), the optimal trajectories z * ξ (t) to the optimal trajectory z * 0 (t), and the minimum values J * ξ of the functional J ξ (u) to the minimum value J * 0 of the functional J 0 (u), as the weighting coefficient ξ tends to zero.
First, we establish the validity of the following lemma.
Lemma 8.1. For the minimum values J * ξ and J * 0 of the corresponding functionals J ξ (u) and J 0 (u) the relationship: is true.
Remark 8.2. Formula (81) denotes the continuity on the right for ξ = 0 of the minimum value J * ξ of the functional J ξ (u). Proof. It is convenient to transform the considered problem (9) for system (1) to the optimal control problem with a terminal functional. For this, for an arbitrary admissible control u(t), to the corresponding solutions l(t), k(t), m(t) of system (1) with initial conditions (2) we add the auxiliary solution n(t), which satisfies the differential equationṅ(t) = (1−u(t))k(t) and the initial condition n(0) = 0. Next, in the Euclidean space R 4 we introduce a non-zero column vector ϑ ξ = (0, 1, 0, ξ) , and define the trajectory y(t) = (l(t), k(t), m(t), n(t)) satisfying the following system of differential equations:ẏ with given initial conditions: Here F (y) and G(y) are the corresponding drift and control vector fields, which are obtained from formulas (17) of the vector fields f (z) and g(z) by adding the fourth coordinates k and (−k), respectively, and y = (l, k, m, n) ∈ R 4 . Then, the functional J ξ (u) of (9) is rewritten in the form of the corresponding scalar product of the vector ϑ ξ with y(T ), that is: Thus, on the set of admissible controls U (T ) we consider for system (82) with initial conditions (83) the minimization of the functional J ξ (u) given by formula (84). Due to the clear and simple relationship of the formulated and original optimal control problems, it is easy to see the relationship between the corresponding optimal solutions: the optimal controls are identical and represent the control u * ξ (t), the minimum values of the functionals are also equal to each other and is J * ξ , the first three components of the optimal trajectories coincide and define the solutions l * ξ (t), k * ξ (t), m * ξ (t). Now, let us consider an arbitrary numerical sequence {ξ r } r≥1 → +0 and the corresponding sequence of optimal controls {u * ξr (·)}. Since {u * ξr (·)} ⊂ U (T ) in the space L 2 ([0, T ]) and the set U (T ) is weakly compact in this space (Theorem 13 ([30])), then the sequence {u * ξr (·)} has a subsequence {u * ξr s (·)}, which weakly converges in L 2 ([0, T ]) to some control u * 0 (·) ∈ U (T ). We will assume that the sequence {u * ξr (·)} itself weakly converges in L 2 ([0, T ]) to this control. Then, due to the well-known result of F. Riesz on the weak convergence of a sequence in space L p ([0, T ]) for p > 1 ( [23]), we find the relationship: lim ξr→+0 ||u * ξr (·) − u * 0 (·)|| ω = 0, where for any element u(·) ∈ U (T ) the norm || · || ω is defined as follows ||u(·)|| ω = max Let {y * ξr (·)} be the sequence of solutions of system (82) corresponding to the sequence of optimal controls {u * ξr (·)}, and y * 0 (·) be the solution of this system corresponding to the control u * 0 (·) as well. Moreover, let us consider a continuous function y(t) = (l(t), k(t), m(t), n(t)) on the interval [0, T ] that is a trajectory of system (82) and corresponds to a control u(·) ∈ U (T ) as an element of a Banach space C([0, T ], R 4 ) with the norm: where ||y|| R 4 is the Euclidean norm of a column vector y = (l, k, m, n) in R 4 .
Relationship (91) means the uniform convergence on the interval [0, T ] of the optimal solutions l * ξ (t), k * ξ (t), m * ξ (t) to the corresponding optimal solutions l * 0 (t), k * 0 (t), m * 0 (t) for ξ = 0. Relationship (92) implies the convergence of the optimal controls u * ξ (t) to the optimal control u * 0 (t) in the so-called metric of a sliding mode. Such a metric is convenient for analyzing optimal control problems ( [9]), especially in the study of problems with a singular arc.
This convergence can be strengthened by Theorem 1 from [42], and the convergence of the controls u * ξ (t) to the control u * 0 (t) is guaranteed as ξ → +0 in the space metric L 1 ([0, T ]), and hence in the space metric L 2 ([0, T ]), if the function u * 0 (t) is bang-bang. In our studies this is possible in Case A (see Section 4).
Thus, as it can be seen from the relationships (81), (91), and (92), as the weight coefficient ξ tends to zero, the optimal controls u * ξ (t) and the corresponding optimal solutions l * ξ (t), k * ξ (t), m * ξ (t) can be considered as the approximations of the appropriate optimal control u * 0 (t) and optimal solutions l * 0 (t), k * 0 (t), m * 0 (t). The confirmation of this fact easily follows from the analysis of the corresponding graphs of Figures 1-4.
Remark 8.5. Other ways of approximating a chattering by piecewise constant controls are presented in [38,44,45,47]. 9. Conclusions. We considered a nonlinear control system of three differential equations that reflects the process of psoriatic treatment. Its phase variables are concentrations of T-lymphocytes, keratinocytes and dendritic cells; the scalar bounded control determines the medication dosage. For this control system, on a given time interval the minimization problem with the Bolza type functional was stated. In this functional, the terminal term determines the concentration of keratinocytes at the terminal time, and the integral term (linear in control and proportional to the concentration of keratinocytes) represents the so-called weighted cost of psoriasis treatment. Further, the issues relating to the properties of the original system and the existence of an optimal solution in this minimization problem were discussed. Analysis of such an optimal solution, consisting of the optimal control and the corresponding to it optimal solutions of the system, was conducted using the Pontryagin maximum principle. Using the Lie brackets technique from the geometric control theory, further we investigated both cases, when the weighting coefficient before the integral term vanishes, and the case when it is positive. As a result, it was shown that in the first case, the corresponding optimal control can be chattering, and therefore does not make much sense as a type of medical treatment for psoriasis. In the second case, it was found that the corresponding optimal control is