key: cord-0969313-axgmeqzn
authors: Marinoschi, Gabriela
title: Parameter Estimation of an Epidemic Model with State Constraints
date: 2021-09-08
journal: Appl Math Optim
DOI: 10.1007/s00245-021-09815-2
sha: f62e7aeaefb0ae3fbd41e295a91cab7be9686222
doc_id: 969313
cord_uid: axgmeqzn

We consider a mathematical model with five compartments relevant to depict the feature of a certain type of epidemic transmission. We aim to identify some system parameters by means of a minimization problem for a functional involving available measurements for observable compartments, which we treat by an optimal control technique with a state constraint imposed by realistic considerations. The proof of the maximum principle is done by passing to the limit in the conditions of optimality for an appropriate approximating problem. The proof of the estimates for the dual approximating system requires a more challenging treatment since the trajectories are not absolutely continuous, but only with bounded variation. These allow to pass to the limit to obtain the conditions of optimality for the primal problem and lead to a singular dual backward system with a generalized solution in the sense of measure. As far as we know, this approach developed here for the identification of parameters in an epidemic model considering a state constraint related to the actions undertaken for the disease containment was not addressed in the literature and represents a novel issue in this paper.

We study an identification problem in relation with a nonlinear system of five differential equations describing an epidemic evolution, focusing on the estimation of B Gabriela Marinoschi gabriela.marinoschi@acad.ro some parameters of the system connected with the actions undertaken for the epidemic containment.

We consider the following nonlinear system of differential equations

(1.1)

3) This system describes a mathematical model of an epidemic with five compartments, represented at time t by the susceptible S(t) (healthy individuals which can acquire the disease), infected asymptomatic A(t) (individuals who have acquired the disease, have no symptoms or only mild ones and have not yet been confirmed), infected with symptoms but not confirmed yet I (t), detected by medical tests and isolated (e.g., by hospitalization) L(t), and recovered R(t). In this model we assume that there are two classes of infected populations A and I which are both infectious, that is individuals from each these classes can spread the virus to the susceptible individuals with the rates β I (t) and β A (t), respectively. The number of individuals in these classes is not known since they are not still checked and confirmed. Since the identification of the parameters we propose here is considered to be performed on short intervals of time, we skip in this model the extinct (from any reason) and the newborns, because their number is negligible compared with the large number of people in the relevant compartments for a disease with large transmissibility and low mortality. Thus, in the perspective of the identification problem approached in this work, we assume that the model contains relevant classes of populations evolving during this epidemic. The study presented in this paper can be applied, of course, to a model completed with equations for other compartments of interest. This model is suitable to describe the infection with SARS-CoV-2, since it is known that the asymptomatic individuals infected with this virus (A in this model) can transmit the disease (see [1] ). The model we consider here may be seen as a reduced version of more general descriptions of a SARS-CoV-2 epidemic, see e.g., those presented and analyzed e.g., in [2] [3] [4] which include many compartments, in particular of hospitalized individuals exhibiting more or less serious symptoms. These classes are included in the compartment L of our model.

The model is given here in a normalized form, that is the real number of individuals at time t in each class is divided by the total population at time t, given by the sum of all components N (t) := S(t) + A(t) + I (t) + L(t) + R(t). By (1.1)-(1.5) we note that

so that we deduce that the total population N (t) remains unchanged in time, that is

Since each term in these sums is nonnegative (as representing a fraction of population), then it is bounded by N . Moreover, since we work with fractions of population, N is actually equal to 1, but we shall keep it written as N to precisely indicate where it occurs.

We are going to explain the meaning of the parameters in order to sustain the mathematical assumptions we shall set for them. The functions β I and β A represent the average number of individuals infected in unit time by an infective symptomatic and asymptomatic. The values σ, l A and μ A are, respectively, the rates at which an asymptomatic can become symptomatic, is detected and isolated, and is recovered. Similarly, l I and μ I are the rates at which a symptomatic is detected and recovered, respectively. The isolated individuals recover with the rate μ L . As in many other transmissible diseases, it was evidenced (see e.g., [5] ) that in the case of infection with SARS-CoV-2 the immunity is gained for a short time only, so that a recovered individual can be infected again with the virus after some time and return into the susceptible compartment with the rate ξ. The removal rates μ I , μ A , μ L and ξ are constant, being represented by the inverse of the average healing time for a symptomatic, asymptomatic, infectious treated in the hospital and the inverse of the loss immunity time, respectively.

If one considers the compartment L as containing all confirmed isolated individuals, the time dependent parameters l A and l I are the rates of detection, relative to asymptomatic and symptomatic cases, respectively (see e.g., [3, 6] ), that is the percentage of infections that are discovered. They may also reflect the number of tests performed over the population and can be modified by enforcing sustained actions, as for example a massive testing campaign (see [7] ). The rate l I may be larger than l A , because a symptomatic individual is more likely to be tested. Here, we assume that not all confirmed individuals remain isolated, so that L contains only a part of the confirmed individuals and note that L does not returns individuals in the infected classes in our model. Thus, we interpret the rates l A and l I as detection-isolation rates.

The possibility of a better prediction of the epidemic evolution can be enhanced by a more accurate knowledge of the transmission parameters and by a correct calibration of the model. Some of the parameters cannot be obtained by direct measurements. For instance, the detection rates are notoriously hard to be measured during an epidemic (see [6] ). This leads to the necessity of calibrating the model by identification techniques using available measured data for the state. There are many recent papers dealing with the parameter estimation by various techniques using a best-fit approach, namely by finding the parameters that locally minimize the sum of the squares of the errors, for models with many compartments well evidenced during the SARS-CoV-2 outbreak, see e.g., [2] [3] [4] . We also refer to the numerical simulations presented in [6] emphasizing the sensitivity of the predictions of the detection rates under various scenarios corresponding to the intervention policies.

We see that if l A = l I = μ L = 0, the class of isolated individuals disappear and so all asymptomatic and symptomatic remain in circulation. Thus, the class L is relevant as a control class for the disease containment and a better isolation action can be done via the coefficients l A and l I , within a time interval [0, T ]. These parameters standing for the rates of isolation, as explained before, cannot be precisely identified by direct observations, since there is not a reliable control on the confirmed but unhospitalized people. Thus, we focus here only on the determination of the parameters l A and l I using the measured confirmed and recovered cases.

Let us assume a certain epidemiological evolution within the interval [0, T ], expressed by known initial data, the observed values L obs (t) and R obs (t), for the classes L and R, and known system parameters, except for l A and l I . In view of a correct analysis of the epidemic evolution a basic requirement is to find the parameters l A and l I which determined the components L and R of the system to evolve according the observed behaviors L obs (t) and R obs (t) within the interval [0, T ].

From the mathematical point of view this means to minimize the functional As a component of the state system, L(t) ∈ [0, N ]. However, as is it unlikely to test or isolate all the population at the same time, it is reasonable to assume that L(t) does not reach the value N , but is bounded by a certain threshold less than N . In particular, the maximum value of L cannot exceed the number of tests done at time t, and we denote this time dependent threshold the testing capacity L test (t). Thus, we minimize the functional along a constrained state

. We assume that L test ∈ C([0, T ]) and assume

We note that the lower bound L ≥ 0 is not actually a constraint because this follows from a property which will be proved for the solution to the state system.

In the sequel, we assume the following conditions for the coefficients of the system:

and denote

We approach problem (P) by an optimal control technique. In Sect. 2, after proving the existence and uniqueness of the solution to the state system (1.1)-(1.6), we show that there exists at least a solution (l * A , l * I ) to problem (P) in Proposition 2.2. For this problem the optimality conditions cannot be directly calculated since the cost functional is not regular because of the constraint. So, one introduces, in Sect. 3.1, an approximating problem (P ε ) indexed along a positive parameter ε, which contains an appropriate penalized term replacing the state constraint. The convergence of a sequence of solutions to (P ε ) precisely to a certain chosen solution to (P) is proved in Proposition 3.2 and the approximating optimality conditions are provided in Proposition 3.3. Certain particular arguments and calculations are led to find appropriate estimates for the solution to the approximating dual system in Proposition 4.1. Actually, they ensure the fact that one component is a bounded in the class of BV (bounded variations) functions. Relying on these, local conditions of optimality for problem (P) are obtained by passing to the limit in the approximating optimality condition system, in Theorem 4.3.

We note that (P), with the state constraint L(t) ∈ [0, L test ], requires a more challenging treatment involving a generalization of the basic duality theory to problems in which the trajectories are not absolutely continuous, but only with bounded variation. As a matter of fact, this is a typical phenomenon for state constraint optimal control problems (see [8] ). Actually, the maximum principle (the first order conditions of optimality) leads to a singular dual backward system with a generalized solution in the sense of measure, which allows for jumps caused by the presence of the state constraint in the primal problem. As far as we know, this approach of identifying the parameters of a SARS-CoV-2 epidemic considering a state constraint imposed in connection with the real actions undertaken for the disease containment, is not found in the literature devoted to epidemics and represents a novelty in this paper.

We begin with the proof of the well-posedness of the state system. Its solution will be sometimes denoted by X := (S, A, I , L, R). 5 and each component belongs to [0, N ] . The solution is continuous with respect to the data (l A , l I ).

Proof In system (1.1)-(1.6) we apply the Banach fixed point theorem, using the set

Let us pick a ∈ M and fix it in the system

2)

with the initial conditions (1.6). We define corresponding to a and a respectively, with the same initial condition. By calculating

where C 1 (T ) is a polynomial in T with coefficients consisting in sums of the constant systems parameters and the L ∞ -norms of the time dependent system parameters. This shows that is a contraction for small T , that is we obtain a local solution. Since all solution components are bounded by N it follows that the solution is global (see e.g., [9] , p. 41, Theorem 2.15). By a similar calculation one can easily prove the continuity with respect to l A and l I .

Proof It is obvious that an admissible pair exists. For example, for l A = l I = 0, by (1.4) and (1.9) we get L(t) ≤ L 0 ≤ L test (t) and J (l A , l I ) < ∞. Let d := inf J (l A , l I ) ≥ 0 and consider a minimizing sequence (l n A , l n I ) n , satisfying the restrictions in (P). The minimizing sequence also satisfies

This implies that l n A → l * A and l n I → l * I weak* in L ∞ (0, T ), as n → ∞ and l * A , l * I ∈ [0, 1] a.e. We denote by X n the solution to the state system corresponding to the minimizing sequence. By Proposition 2.1 this solution exists and is bounded in

Then, on a subsequence, X n → X * , (X n ) → (X * ) weak* in (W 1,∞ (0, T )) 5 . By Arzelà theorem it follows that X n → X * uniformly in [0, T ], and so X n (0) → X * (0) = X 0 and X n (T ) → X * (T ). By passing to the limit in (1.1)-(1.6) we deduce that (X * ) (t) = lim n→∞ (X n ) (t), whence it follows that X * is the solution to (1.1)-

Finally, by using the weakly lower semicontinuity of the norms we get lim n→∞ J (l n A , l n I ) = J (l * A , l * I ) and so (l * A , l * I ) turns out to be optimal in (P).

Let ε > 0 and let (l * A , l * I ) be optimal in (P). We introduce the adapted approximating cost functional

and study the following approximating problem (P ε ),

subject to (1.1)-(1.7). We observe that the state constraint in (P) is replaced here by the penalization of the L 2 -norm of the positive part of (L(t) − L test (t)), where α 2 > 0. The last two penalization terms in (3.1) will ensure the convergence of the approximating solution to the chosen optimal controller (l * A , l * I ) in (P). Proposition 3.1 Problem (P ε ) has at least one solution, (l * A,ε , l * I ,ε ) with the corresponding state X * ε ∈ (W 1,∞ (0, T )) 5 , having the components in [0, N ].

After observing that there is at least an admissible triplet (let it be an optimal one in (P)), the proof is led then exactly as in Proposition 2.2, using in addition the continuity of the function L → L + .

A , l * I ), X * } and {(l * A,ε , l * I ,ε ), X * ε } ε be optimal in (P) and in (P ε ), respectively. Then,

for all (l A , l I ) satisfying the constraints in (P ε ). In particular, we can set (l A , l I ) = (l * A , l * I ) a.e. on (0, T ) which is the optimal triplet chosen in (P) and then the previous inequality becomes

because the optimal L * ≤ L test in (P). Thus, the left-hand side is bounded and we deduce that, on a subsequence, l * A,ε → l * * A , l * I ,ε → l * * I weakly in L 2 (0, T ), weak* in L ∞ (0, T ) as ε → 0, and l * * A , l * * I ∈ [0, 1]. As in Proposition 2.2, we infer that X ε → X * * weak* in (W 1,∞ (0, T )) 5 

In the inequality (3.4) without the last three terms on the left-hand side we pass to the limit using the weakly lower semicontinuity of the norms and deduce that J (l * * A , l * * I ) ≤ J (l * A , l * I ), hence (l * * A , l * * I ) is optimal in (P), too. Moreover, still by (3.4) we easily infer that l * A,ε → l * A and l * I ,ε → l * I strongly in L 2 (0, T ), hence (l * * A , l * * I ) = (l * A , l * I ) a.e. This ends the proof .

Let λ > 0 and set the variations ∈ (0, 1) .

Let us denote

is the solution to (1.1)-(1.7) corresponding to (l λ A,ε , l λ I ,ε ), and X * ε is the optimal state corresponding to (l * A,ε , l * I ,ε ). We introduce the linearized system for problem (P) where

First, by the known results for linear systems we infer that (3.5) has a unique global solution x = (s, a, i, l, r ) ∈ (W 1,∞ (0, T )) 5 . By a direct calculation, using the continuity with respect to the data of the solution to the state system it can be easily proved that x λ ε → x ε := (s, a, i, l, r ) strongly in C([0, T ]), as λ → 0, so that (3.5) stands for the system in variations.

We introduce the backward dual system for the variables ( p ε , q ε , d ε , e ε , f ε ) as 

We recall that the projection of a point z on a set K ⊂ R d , d ∈ N, is defined by

for all k > 0, I d being the identity operator.

where ( p ε , q ε , d ε , e ε , f ε ) is the solution to the backward dual system (3.8)- (3.13) .

Proof Let us multiply the equations for s, a, i, l, r in (3.5) by p ε , q ε , d ε , e ε , f ε , respectively and integrate over (0, T ). By integrating by parts and taking into account the dual system and the initial conditions in the system in variations, we obtain

On the other hand, for (l * A,ε , l * I ,ε ) optimal in (P ε ) we can write

By replacing the expression of the cost functional J ε , performing some algebra, dividing by λ and passing to the limit as λ → 0 we obtain

(3.16)

By comparison with (3.15) we deduce that

We recall the setting of ω A and ω I and choose, in particular, l A = l * A,ε , meaning that we keep l * A,ε fixed and give a variation only to l * I ,ε . Then, (3.17) yields

whence we have the second relation in (3.14). Then, let us set l I = l * I ,ε . By (3.17 ) we obtain that

for l A (t) ∈ [0, 1] a.e. t ∈ (0, T ), which yields the first relation in (3.14).

We begin by proving the boundedness of the solution to the dual system.

|e ε (t)| + T 0 e ε (t) dt ≤ C, f or all t ∈ [0, T ], (4.1) | p ε (t)| + |q ε (t)| + |d ε (t)| + | f ε (t)| ≤ C, f or all t ∈ [0, T ],(4.

with C independent of ε.

Proof By (1.9) we have L min = min t∈[0,T ] {L test (t)} and recall that L min ≥ L 0 . Let us choose y 1 ∈ (0, L min ), ρ > 0 and θ ∈ R, with |θ | = 1 such that y 1 + ρθ ∈ (0, L min ). This can be equivalently written 0 ≤ y 1 ± ρ < L min and holds if 0 < ρ < min{y 1 , L min − y 1 }.

The values y 1 and ρ are at our free choice and different situations for the positions of L 0 , L min and y 1 will be discussed later, at the appropriate place. For the moment, let us multiply (3.11) by (L * ε − (y 1 + ρθ)) and integrate over (0, T ). We have

Let us define ϕ(y) = 1 2 (y + ) 2 , ϕ : R → [0, ∞). Its subdifferential ∂ϕ(y) = y + and ϕ(y) − ϕ(z) ≤ ∂ϕ(y)(y − z) for all z ≥ 0. In the second term on the left-hand side of (4.4) we set

where 1 {t;e ε (t) =0} is the characteristic function of the set indicated as subscript. Using the relation for the subdifferential of ϕ and integrating by parts the first term on the left-hand side of (4.4) we have

Since y 1 + ρθ < L min ≤ L test (t), e ε (T ) = 0, using Eq. (1.4) and making some rearrangements we are led to the relation

We go back to (4.7) and using (4.8) we write

This implies

We have to check if the coefficient of E ε is positive. We note that the function t → 2(μ L + 1)(1 + 4 √ te Ct )N t is positive, strictly increasing for t > 0 and vanish at 0. Thus, if

then there exists T loc > 0 such that the inequality

takes place on (0, T loc ). Actually, T loc is the solution to the equation

It remains to check if one can choose y 1 and ρ to satisfy (4.10), taking into account the possible situations for L 0 and L min which are fixed, with L 0 ≤ L min .

I. Let L 0 ≤ L min 2 . We choose L 0 < y 1 < L min 2 and ρ ∈ (y 1 − L 0 , y 1 ). Then, min{y 1 , L min − y 1 } = y 1 and both (4.3) and (4.10) are satisfied.

II. Let L 0 ∈ L min 2 , L min . We choose y 1 ∈ L min 2 , L 0 ⊂ L min 2 , L min and ρ ∈ (L 0 − y 1 , L min − y 1 ) . Then, min{y 1 , L min − y 1 } = L min − y 1 and both (4.3) and (4.10) are satisfied.

Going back to (4.9) we can write that

for T ∈ (0, T loc ). (4.13) with C independent of ε. This and (4.6) imply (4.1 ). Recalling (4.8) we obtain (4.2) . This ends the proof. We also recall by the Lebesgue decomposition theorem that every μ ∈ (L ∞ (0, T )) * can be uniquely written as

where μ a ∈ L 1 (0, T ) and μ s is a singular measure (see e.g. [10] ). This means that for and v s is a singular part (for instance it can be a jump function with bounded variation or a function with bounded variation with a.e. zero derivative). We note that if v ∈ BV ([0, T ]), then its distributional derivative dv dt := μ belongs to (L ∞ (0, T )) * , and in virtue of the Lebesgue decomposition, it is represented by the sum of the absolutely continuous part and the singular part

where D (0, T ) is the space of Schwartz distributions on (0, T ).

In the next theorem we shall pass to the limit in the approximating optimality conditions. To this end, we introduce the system where ν := α 2 lim

We also define

We denote by N K (t) (z) the normal cone to K (t) at z,

and by (4.12) , the optimality conditions for problem (P) read: where ν a (t) ∈ N K (t) (L * (t)), a.e. t ∈ (0, T ), (4.31)

Since ν s is singular with respect to the Lebesgue measure, this inequality implies separately that each term in (4.39) is greater or equal to 0. If t 0 ∈ (0, T ) is such that L * (t 0 ) ≤ L test (t 0 ), then choosing in the first part

where z 0 ≤ L test (t 0 ), and ε > 0 is arbitrary small, we get 1 ε t 0 +ε t 0 −ε ν a (t)(L * (t) − z 0 )dt ≥ 0, for all ε > 0.

Letting ε → 0, we obtain ν a (t 0 )(L * (t 0 ) − z 0 ) ≥ 0 for almost all Lebesgue points t 0 and since z 0 is arbitrary, the relation ν a (t) ∈ N K (t) (L * (t)), a.e. t ∈ (0, T ) is inferred. This phenomenon is typical for state constraints optimal control problems (see [8] ). Therefore, we can interpret that the optimal rates l * A and l * I are absolutely continuous on the set {t; L * (t) < L test (t)} where the number of positive detected and isolated individuals is less than the number of tests, and become singular, having possible jumps at the moments when all tested individuals are found infected and are isolated.

We solved an inverse problem for the identification of the detection-isolation rates of an epidemic model, on the basis of the observations made for the confirmed isolated and recovered compartments. Using these determined rates one can estimate the number of the individuals in all compartments, especially in the infectious classes. This can offer information about the epidemiologic situation, which may help the policies for modifying these rates in view of a better disease containment. In this respect, an additional estimation of the rates l A and l I that may contribute to the minimization

Features, Evaluation and Treatment of Coronavirus (COVID-19)

Identifying the measurements required to estimate rates of COVID-19 transmission, infection, and detection, using variational data assimilation

Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy

Modeling the transmission of COVID-19 in the US: a case study

Incubation period and other epidemiological characteristics of 2019 novel coronavirus infections with right truncation: a statistical analysis of publicly available case data

The significance of case detection ratios for predictions on the outcome of an epidemic: a message from mathematical modelers

Covid-19 mass testing facilities could end the epidemic rapidly

State constraints in convex control problems of Bolza

Differential Equations

Integrals which are convex functionals

Convexity and Optimization in Banach Spaces

The stochastic reflection problem in Hilbert spaces

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations

because L * ε , L obs and |y 1 + ρθ| are smaller than N . We denote E ε := T 0 e ε (t) dt, F 0 := |L 0 − y 1 | , (4.5) and since e ε (T ) = 0, we note thatThen, we have (because l * A,ε , l * I ,ε ≤ 1)Let us denote S 2 ε (t) = | p ε (t)| 2 + |q ε (t)| 2 + |d ε (t)| 2 + | f ε (t)| 2 . We multiply Eqs. (3.8)-(3.10), (3.12) by p ε , q ε , d ε , f ε respectively, integrate over (t, T ) and sum up. After some calculations, and a variable change t → T − t, we obtainHere, ε consists in sums of the constant coefficients of the equations in the dual system plus the L ∞ -norms of the time dependent coefficients (or square coefficients) and G ε = 2( ε + l * 2 A,ε + l * 2 I ,ε + α 2 R ) ≤ C, since all terms in G ε are bounded. Further, C may denote several positive constants. By the Gronwall's lemma, the previous inequality yields the estimate S 2 Proof We shall establish some estimates in order to pass to the limit in the approximating optimality conditions determined in Proposition 3.3. Let T < T loc , let (l * A , l * I ) be optimal in (P) with the corresponding state X * and let us consider {(l * A,ε , l * I ,ε ), X * ε } optimal in (P ε ). Recalling Proposition 3.2 and the continuity property from Propo- 5 and uniformly on [0, T ].By (4.1) and (4.2) each component of the solution to the dual system follows to be bounded in L ∞ (0, T ). Moreover, by (3.8)-(3.10), (3.12) it follows thatThe component (e ε ) ε is bounded in L ∞ (0, T ) by (4.6) and its derivative is bounded in L 1 (0, T ), by (4.1). These imply that e ε ∈ BV ([0, T ]) and by Helly's theorem (see e.g., [11] Going back to (3.11) we deduce that 1. This is pointed out in the proof of Corollary 2B in [10] , but this assertion does not follow directly from Alaoglu theorem. An argument can be found in [12] , and we resume it below.Let us consider the linear operator : C([0, T ]) → L ∞ (0, T ), z = , which maps a continuous function into the corresponding class of equivalence (of all functions a.e. equal). Its adjoint * :and using the Alaoglu theorem it follows that ( * μ n ) n is weak* sequentially compact in M([0, T ]). Therefore, it follows that (μ n ) n is weak* sequentially compact in M([0, T ]). Passing to the limit in μ n ( z) = ( * μ n )(z) we get μ( z) = ( * μ)(z) for any ∈ L ∞ (0, T ) which is of the form z with z ∈ C([0, T ]). Then, due to the Hahn-Banach theorem, μ can be extended to all L ∞ (0, T ) and so we conclude that (μ n ) n is weak* sequentially compact in (L ∞ (0, T )) * .Therefore, one can extract a subsequence such thatThus, relying on (4.33)-(4.35), we can pass to the limit in (3.8 )-(3.13) and obtain (4.17)-(4.23). Now, we move to (3.14), or more exactly in (3.18) and (3.19 ) and pass to the limit. The left-hand side of (3.18) converges and since the normal cone is maximal monotone, hence strongly-strongly closed it follows thatwhich implies the second relation in (4.27) . Similarly, we proceed in (3.19 ) and obtainwhence we get the first relation in (4.27). Finally, we want to interpret Eq. (4.20). Since z + is the subdifferential of the function 1 2 z + 2 , we can writefor all z ∈ K. At limit we obtain ν(L * − z) ≥ 0, for all z ∈ K, (4.36) whence, ν ∈ N K (L * ).Since ν ∈ (L ∞ (0, T )) * ⊂ M([0, T ]), recalling (4.14) and (4.16), we can represent ν = ν a + ν s , where ν a is the absolutely continuous part (in the sense of measure) and ν s is the singular part of ν. Also, e ∈ M([0, T ]) and e = (e ) a + (e ) s where (e ) a = e a ∈ L 1 (0, T ). Then, (4.20) can be rewritten as in (4.29)-(4.30). Relation (4.36) implies that ν a (t) ∈ N K (t) (L * (t)), a.e. t ∈ (0, T ), that is(4.37)As regards ν s we haveHere is the argument. Relation (4.36) can be still written This problem can be similarly treated as problem (P) with slight differences in some calculations due to the new cost functional. The appropriate approximating problem isThe approximating optimality conditions are calculated as in Proposition 3.3, the approximating optimal controls are still given by (3.14), but the dual system is