key: cord-0162699-7euag84a
authors: DIM, Emilio Molina; CMM,; LJLL,; MISTEA,; Rapaport, Alain; Ramirez, Hector
title: Equivalent formulations of optimal control problems with maximum cost and applications
date: 2022-02-25
journal: nan
DOI: nan
sha: 4584d37493500f4549b9d5f19cdb12f766e5a0e8
doc_id: 162699
cord_uid: 7euag84a

We revisit the optimal control problem with maximum cost with the objective to provide different equivalent reformulations suitable to numerical methods. We propose two reformulations in terms of extended Mayer problems with constraint, and another one in terms of a differential inclusion with upper-semi continuous right member but without constraint. For this last one we also propose an approximation scheme of the optimal value from below. These approaches are illustrated and discussed on several examples.

We consider the optimal control problem which consists in minimizing the maximum of a scalar function over a where y(t) = θ(t, ξ(t)) with ξ(·) solution of a controlled dynamicsξ = φ(ξ, u), ξ(t 0 ) = ξ 0 . This problem is not in the usual Mayer, Lagrange or Bolza forms of the optimal control theory, and therefore is not suitable to use the classical necessary optimality conditions of Pontryagin Maximum Principle or existing solving algorithms (based on direct method, shooting or Hamilton-Bellman Jacobi equation). However, this problem falls into the class of optimal control with L ∞ criterion, for which several characterizations of the value function have been proposed in the literature [3, 4, 9] . Typically, the value function is solution, in a general sense, of a variational inequality of the form

without boundary condition. Nevertheless, although necessary optimality conditions and numerical procedures have been formulated [2, 6, 7, 8] , there is no practical numerical tool to solve such problems as it exists for Mayer problems, to the best of our knowledge. The aim of the present work is to study different reformulations of this problem into Mayer form in higher dimension with possibly state or mixed constraint, for which existing numerical methods can be used. Indeed, it has already been underlined in the literature that discrete-time optimal control problems with maximum cost do not satisfy the Principle of Optimality but can be transformed into problems of higher dimension with additively separable objective functions [10, 11] . We pursue here this idea but in the continuous time framework, which faces the lack of differentiability of the max function. The paper is organized as follows. In Section 2, we give the setup, hypotheses and define the problem. In Section 3, we give equivalences with two Mayer problems with fixed initial condition, under state or mixed constraint. In Section 4, we propose another formulation without constraint in terms of differential inclusion, and then show how the optimal value can be approximated from below by a sequence of more regular Mayer problems. Section 5 is devoted to numerical illustrations. We first consider a very particular class of problems for which we are able to give explicitly the optimal solution, which allows to compare the numerical performances of the different formulations. We then consider a more sophisticated problem from epidemiology, and discuss the various issues in numerical implementations of the different formulations. We also compare numerically with L p approximations. Finally, we discuss in Section 6 about the potential merits of the different formulations as practical methods to compute optimal solution of L ∞ control problems.

We shall consider autonomous dynamical systems defined on a invariant domain D of R n+1 of the form

(where g is a scalar function) with u ∈ U ⊂ R p . Throughout the paper, we shall assume that the following properties are fulfilled.

i. U is a compact set.

ii. The maps f and g are C 1 on D × U .

iii. The maps f and g have linear growth, that is there exists a number C > 0 such that

For instance, y(·) can be a smooth output of a dynamicṡ

which can be rewritten as

Let U be the set of measurable functions u(·) : [0, T ] → U and consider (x 0 , y 0 ) ∈ D, T > 0. Under the usual arguments of the theory of ordinary differential equations, Assumption 1 ensures that for any u(·) ∈ U there exists an unique absolutely continuous solution (x(·), y(·)) of (1) on [0, T ] for the initial condition (x(0), y(0)) = (x 0 , y 0 ). Define then the solutions set S := {(x(·), y(·)) ∈ AC([0, T ], R n+1 ), sol. of (1) for u(·) ∈ U with (x(0), y(0)) = (x 0 , y 0 )} We consider then the optimal control problem which consists in minimizing the "peak" of the function y(·):

A first approach considers the family of constrained sets of solutions

and to look for the optimization problem inf{z; S z = ∅} This problem can be reformulated as a Mayer problem

where z(0) is free. Direct methods can be used for such a problem. However, as z(0) is free, solutions are not sought among solutions of a Cauchy problem, which prevents using other methods based on dynamic programming such as the Hamilton-Jacobi-Bellman equation.

We propose another extended dynamics in D × R with an additional control v(·) ∈ [0, 1]

Let V be the set of measurable functions v : [0, T ] → [0, 1]. Note that under Assumption 1, for any (x 0 , y 0 , z 0 ) ∈ D × R and (u, v) ∈ U × V, there exists an unique absolutely solution (x(·), y(·), z(·)) of (2) on [0, T ] for the initial condition (x(0), y(0), z(0)) = (x 0 , y 0 , z 0 ). Here, we fix the initial condition with z 0 = y 0 and consider the Mayer problem

under the constraint C and shows its equivalence with problem P. We first consider fixed controls u(·). 

and is reached for a control v(·) that takes values in {0, 1}.

Proof. From equations (2), one get that any solution z(·) is non decreasing, and as z satisfies the constraint z ≥ y, we deduce that one has z(T ) ≥ max

for any solution of (2), and thus max t∈[0,T ]

Let x(·), y(·) be the solution of (1) for the control u(·) and let I be the set of invisible points from the left of y, that is I := {t ∈ (0, T ); y(t ) > y(t) for some t < t} .

Consider then the control

When I is empty, y(·) is a non decreasing function, and with the control v = 1, one has z(t) = y(t) for any t ∈ [0, T ]. Therefore one has z(T ) = y(T ) = max

When I is non empty, there exists, from the sun rising Lemma [13] , a countable set of disjoint non-empty intervals I n = (a n , b n ) of [0, T ] such that -the interior of I is the union of the intervals I n , -one has y(a n ) = y(b n ) if b n = T , -if b n = T , then y(a n ) ≥ y(b n ).

Note that when t / ∈ int I, one has y(t) ≥ y(t ) for any t ≤ t. Therefore, the solution z with control (6) verifies z(t) = y(t), t / ∈ int I y(a n ), t ∈ I n for some n (see Figure 1 as an illustration). Lett ∈ [0, T ] be such that Remark 3.1. The proof of Proposition 3.2 gives an optimal construction of z(·) which is the lower envelope of non decreasing continuous functions above the function y(·), as depicted on Figure 1 . However, there is no uniqueness of the optimal control v(·). Any admissible solution z(·) that is above y(·) and such that z(t) =ŷ for t ≥t = min{t ∈ (0, T ], y(t) =ŷ}, whereŷ := max s∈[0,T ] y(s), is also optimal.

We then obtain the equivalence between problems P 1 and P in the following sense.

, v (·)) is optimal for Problem P 1 , then u (·) is optimal for Problem P. Conversely, if u (·) is optimal for Problem P, then (u (·), v (·)) is optimal for Problem P 1 where v (·) is optimal for the problem (3) for the fixed control u (.).

Let us give another equivalent Mayer problem but with a mixed constraint (this will be useful in the next section). We consider again the extended dynamics (2) with control v ∈ [0, 1] and the initial condition (x(0), y(0), z(0)) = (x 0 , y 0 , y 0 ), and define the mixed constraint

with the optimal control problem

Problems P 1 and P 2 are equivalent.

Proof. One can immediately see that for any admissible solution that satisfies constraint C, the constraint C m is necessarily fulfilled as max(y − z, 0) is identically null. Conversely, fix an admissible control u(·) and consider a control v(·) that satisfies C m . We show that this implies that the solution (y(·), z(·)) verifies necessarily z(t) ≥ y(t) for any t ∈ [0, T ]. If not, consider the non-empty set

which is open as z − y is continuous. Note that one hasż(t) −ẏ(t) ≥ 0 for a.e. t ∈ E. Therefore z − y is non decreasing in E and we deduce that for any t ∈ E, the interval [0, t] is necessarily included in E, which then contradicts the initial condition z(0) = y(0).

We posit Π = (x, y, z) ∈ D × R and consider the differential inclusioṅ

Let Π 0 = (x 0 , y 0 , y 0 ) and denote by S the set of absolutely continuous solutions of (7) with Π(0) = Π 0 ∈ D × R. We consider the Mayer problem

Under Assumption 2, problem P 3 admits an optimal solution. Moreover, any optimal solution Π(·) = (x(·), y(·), z(·)) verifies z(T ) = max

with (x(·), y(·)) solution of (1) for some control u(·) ∈ U that is optimal for problem P.

Proof. We fix the initial condition Π(0) = Π 0 and consider the augmented dynamicṡ

Under Assumption 2, the values of F † are convex compact. One can straightforwardly check that the set-valued map F † is upper semi-continuous 1 with linear growth. Therefore, the reachable set S † (T ) (where S † denotes the set of absolutely continuous solutions of (8) with Π(0) = Π 0 ) is compact (see for instance [1, Proposition 3.5.5]). Then, there exists a solution Π (·) = (x (·), y (·), z (·)) of (8) which minimizes z(T ).

Note that any admissible solution (x(·), y(·), z(·)) of system (2) that satisfies the constraint C m belongs to S ⊂ S † . We then get the inequality

Let us show that any solution Π(·) = (x(·), y(·), z(·)) in S verifies

We show that one has z(t) ≥ y(t) for any t ∈ [0, T ]. We proceed by contradiction, as in the proof of Proposition

by continuity, that one has z(0) − y(0) < 0 which contradicts the initial condition z(0) = y(0). Moreover, as the map h is non-negative, z(·) is non decreasing and we conclude that (10) is verified. On another hand, thanks to Assumptions 1 and 2, we can apply Filippov's Lemma to the set-valued map G, which asserts that (x(·), y(·)) is solution of (1) for a certain u(·) ∈ U. With(10), we obtain

where (x (·), y (·)) is solution of (1) for a certain u (·) ∈ U.

Finally, inequalities (9) and (11) with Propositions 3.2 and 3.3 show that z (T ) is reached by a solution of (2) under the constraint C m , and that u (·) is optimal for problem P. We also conclude that the optimal value z (T ) is reached by a solution in S , which is thus optimal for problem P 3 .

Remark 4.1. Let us stress that the function h is not continuous, which does not allow to use Filippov's Lemma for the set valued map F . This means that one cannot guarantee a priori that an absolutely continuous solution Π(·) = (x(·), y(·), z(·)) can be synthesized by a measurable control (u(·), v(·)). Proposition 4.1 shows that (x(·), y(·)) is indeed a solution of system (1) for a measurable control u(·), but one cannot guarantee a priori that z(·) can be generated by a measurable control v(·), what does not matter for our purpose.

We propose now an approximation from below of the optimal cost with a continuous dynamics. In minimization problems, approximations from below of the optimal value are useful to frame the optimal value of the problem, upper bounds being given by any sub-optimal control of problem P 0 , P 1 , P 2 or P 3 (provided typically by a numerical scheme). This will be illustrated in Section 5. Let us consider the family of dynamics parameterized by

with h θ (x, y, z, u, v) = max(g(x, y, u), 0)(1 − v e −θ max(y−z,0) ) (where the expression e −θ max(y−z,0) plays the role of an approximation of 1 R + (z − y) when θ tends to +∞). We then define the family of Mayer problems P θ 3 : inf

where S θ denotes the set of absolutely continuous solutions Π(·) = (x(·), y(·), z(·)) of (12) for the initial condition Π(0) = Π 0 . Let us underline that for these problems with Lipschitz dynamics without constraint, necessary conditions based on Pontryagin Maximum Principle can be derived, leading to shooting methods that are known to very accurate and that could be initialized from numerical solutions of problems P 1 or P 2 obtained for instance with direct methods.

Proposition 4.2. Under Assumption 2, for any increasing sequence of numbers θ n (n ∈ N) that tends to +∞, the problem P θn 3 admits an optimal solution, and for any sequence of optimal solutions (x n (·), y n (·), z n )(·)) of P θn 3 , the sequence (x n (·), y n (·)) converges, up to sub-sequence, uniformly to an optimal solution (x (·), y (·)) of Problem P, and its derivatives weakly to (ẋ (·),ẏ (·)) in L 2 . Moreover, z n (T ) is an increasing sequence that converges to max t∈[0,T ] y (t).

Proof. As in the proof of Proposition 4.1, we consider for any θ > 0 the convexified dynamics

where α ∈ [0, 1]. Then, there exists an absolutely continuous solution (x θ (·), y θ (·), z θ (·)) with a measurable control (u θ (·), v θ (·), α θ (·)) which minimizes z(T ). For the control (u θ (·), v θ (·), 0), the solution is given by (x θ (·), y θ (·),z θ (·)) wherez θ (·) is solution of the Cauchy probleṁ

One can check that the inequalityl

, z ∈ R is fulfilled, which gives by comparison of solutions of scalar ordinary differential equations (see for instance [14] ) the inequalityz

We deduce that (x θ (·), y θ (·), z θ (·)) is necessarily a solution of (12) . (1) By Proposition 4.1, we know that there exists an optimal solution (x(·), y(·), z(·)) of problem P 3 such that z(T ) =ȳ. Clearly, this solution belongs to S θ for any θ, and we thus get

Let 

Consider an increasing sequence of numbers θ n (n ∈ N), and denote Π n (·) = (x n (·), y n (·), z n (·)) an optimal solution of problem P θn 3 . Note that one has S θn+1 ⊂ S θn · · · ⊂ S θ0

Therefore, the sequenceΠ n (·) is bounded, and Π n (·) as well. As F is upper semi-continuous, we obtain that Π n (·) converges uniformly on [0, T ], up to a sub-sequence, to a certain Π (·) = (x (·), y (·), z (·)) which belongs to S l (see for instance [5, Th. 3.1.7]). From property (15), we obtain that z n (T ) is a non decreasing sequence that converges to z (T ), and from (13), we get passing at the limit

On another hand, (x (·), y (·), z (·)) belongs to S l and we get from Proposition 4.1 the inequality z (T ) ≥ȳ Therefore, one has z (T ) =ȳ and (x (·), y (·), z (·)) is then an optimal solution of problem P 3 . From Proposition 4.1, we obtain that one has necessarily z (T ) = max

Finally, the sequence (ẋ n (·),ẏ n (·)) being bounded, it converges, up to a sub-sequence, weakly to (ẋ (·),ẏ (·)) in L 2 tanks to Alaoglu's Theorem.

We begin by illustrating the different formulations on a problem for which the optimal solution is known.

We consider dynamics of the form

is optimal for problem P.

Proof. For a given x 0 in R n , let x(·) be the solution ofẋ = f (x), x(0) = x 0 independently to the control u(·). Then, for any solution y(·), one has

Let y (·) be defined as

where y (·) is a solution of Σ for any measurable control u (·) such that g(x(t), u (t)) = min We conclude that y (·) is an optimal trajectory of problem P for the control generated by the feedback φ .

As a toy example, we have considered the system

is an optimal control which minimizes max t∈[0,T ] y(t). Remark that this problem can be equivalently written with a scalar non-autonomous dynamicsẏ

for which the open-loop control

is optimal.

For T = 5, we have first computed the exact optimal solution of problem P with the open-loop u (·), by integrating the dynamics with Scipy in Python software (see Figure 2 ). Impacts of perturbations on the switching times on the criterion are presented in Table 1 , which show a quite high sensitivity of the optimal control for this problem. Then, we have solved numerically problems P 0 to P 2 with a direct method (Bocop software using Gauss II integration scheme) for 500 time steps and an optimization relative tolerance equal to 10 −10 . For problem P 3 , as the dynamics is not continuous, direct methods do not work well and we have used instead a numerical scheme based on dynamic programming (BocopHJB software) with 500 time steps and a discretization of 200 × 200 points of the state space. For the additional control v, we have considered only two possible values 0 and 1 as we know Table 2 , while Figure 3 presents the corresponding trajectories. We note that the direct method give very accurate results, and the computation time for problem P 0 is the lowest because it has only one control. The computation time for problem P 2 is slightly higher than for P 1 because the mixed constraint C m is heavier to evaluate. The numerical method for problem P 3 is of completely different nature as it computes the optimal solution for all the initial conditions on the grid, which explains a much longer computation time. The accuracy of the results is also directly related to the size of the discretization grid and can be improved by increasing this size but at the price of a longer computation time.

On Figure 3 , one may notice some difference between the obtained trajectories. Let us underline that after the peak of y(·), there is no longer uniqueness of the optimal control.

The SIR model is one of the most basic transmission model in epidemiology for a directly transmitted infectious disease (for a complete introduction, see for instance [15] ) and it retakes great importance nowadays due to covid-19 epidemic.

Consider on a time horizon [0, T ] variables S(t), I(t) and R(t) representing the fraction of susceptible, infected and recovery individuals at time t ∈ [0, T ], so that one has S(t) + I(t) + R(t) = 1 with S(t), I(t), R(t) ≥ 0. Let β > 0 be the rate of transmission and γ > 0 the recovery rate. Interventions as lock-downs and curfew are modeled as a factor in rate transmission that we denote u and which represents our control variable taking values in [0, u max ] with u max ∈ (0, 1), where u = 0 means no intervention and u = u max the most restrictive one which reduces as much as possible contacts among population. The SIR dynamics including the control is then given by the following equations:Ṡ = −(1 − u)βSI (16)

When the reproduction number R 0 = β/γ is above one and the initial proportion of susceptible is above the herd immunity threshold R −1 0 , it is well known that there is an epidemic outbreak. Then, the objective is to minimize the peak of the infected population max

with respect to control u(·) subject to a L 1 budget

on a given time interval [0, T ] where T is in general chosen large enough to ensure the herd immunity of the population is reached at date T . Note that one can drop the R dynamics to study this problem. If the constraint (19) were not imposed, then the optimal solution would be the trivial control u(t) = u max , t ∈ [0, T ], which is in general unrealistic from a operational point of view. A similar problem has been considered in [12] but under the constraint that intervention occurs only once on a time interval of given length, that we relax here. Note that the constraint (19) can be reformulated as a target condition, considering the augmented dynamicṡ

with initial condition C(0) = Q and target {C ≥ 0}. Extension of the results of Sections 3 and 4 to problems with target do not present any particular difficulty, and is left to the reader. The parameters considered for the numerical simulations are given in Table 3 . Adding the z-variable, we end β γ T Q S(0) I(0) 0.21 0.07 300 28 1 − 10 −6 10 −6 Table 3 : SIR parameters considered in numerical computations up with a dynamics in dimension four, which is numerically heavier than for the previous example. In particular, methods based on the value function are too time consuming to obtain accurate results for refined grids in a reasonable computation time. So we have considered direct methods only. We do not consider here problem P 3 , but instead its regular approximations P θ 3 suitable to direct methods. For direct methods that use algebraic differentiation of the dynamics, convergence and accuracy are much better if one provides differentiable dynamics. This is why we have approximated the max(·, 0) operator for problems P 1 and P 2 by the Laplace formula

with λ = 100 for the numerical experiments. For problem P θ 3 , one has to be careful about the interplay between the approximations of max(·, 0) and the sequence θ n → +∞, to provide approximations from below of the optimal value. The function h θ is thus approximated by the expression h θ (x, y, z, u, v) log e λ1g(x,y,u) + 1

log(e λ 2 (y−z) +1) which depends on three parameters λ 1 , λ 2 and θ. Posit for convenience α := θ λ 2 and consider the function ω α,λ2 (ξ) := e −α log(e −λ 2 ξ +1) , ξ ∈ R which approximates the indicator function 1 R + . One has the following properties. 2. For any ε ∈ (0, 1), one has ω α,λ2 −ε 2 = ε and ω α,λ2 (0) = 1 − ε exactly for

Proof. One has first Finally, with simple algebraic manipulation of the conditions ω α,λ2 −ε 2 = ε and ω α,λ2 (0) = 1 − ε, one obtains straightforwardly the expressions (23).

We have taken λ 1 = 5000 and considered a sequence of approximations of the indicator function for the values given in Table 4 according to expressions (23) of Lemma 5.1 (see Figure 4 ). Computations have been performed with Bocop software on a standard laptop computer (with a Gauss II integration scheme, 600 time steps and relative tolerance 10 −10 ). As one can see in Figure 5 and Table 5 problems P 0 , P 1 , P 2 present similar performances for peak values and computation time. In Figure 6 and Table 6 , the numerical solutions of P θ 3 are illustrated for the values of α and λ 2 given in Table 4 . As expected, the numerical computation of the family of problems P θ 3 provides an increasing sequence of approximation from below of the optimal value and thus complements the computation of problems P 0 , P 1 or P 2 . From Figures of Tables 5 and 6 , one can safely guarantee that the optimal value belongs to the interval [0.1010, 0.1015]. However, the trajectories found for P θ 3 are not as closed as the ones of problems P 0 , P 1 or P 2 . This can be explained by the fact that problems P θ 3 are not subject to the constraint z(t) ≥ y(t) and thus provides trajectories for which z(T ) is indeed below max t y(t). Figure 6 : Comparison of the numerical results for problem P θ 3 p ∈ {2, 5, 10, 15} (see Figure 7) . Besides, to ensure convergence it was necessary take 1200 time step instead of 600 as in previous simulations. The total time of the process is 78s after summing computation times given in Table 7 : Comparison of the numerical results with the L p approximation other methods. Moreover, the same method for p = 15 but initialized from the solution found for p = 2 gives poor results for a computation time of 50s (see Figure 8 ). We conclude that the L p approximation is not practically reliable for this kind of problems.

In this work, we have presented different formulations of optimal control problems with maximum cost in terms of extended Mayer problems with fixed initial condition, and tested them numerically. We have proposed two classes of problems: one with state or mixed constraint suitable to direct methods, and another one without constraint but less regular and suitable to dynamical programming methods. Moreover, for this last class, we have proposed and approximation scheme with a sequence of regular Mayer problems without constraint, which turned out to give better results than approximations with the L p norms. Although this second approach requires larger computation time, it complements the first one providing approximations of the optimal value from above. Finally, we summarize some advantages and drawbacks of the different formulations for the use of numerical methods in Table 8 .

P 0 P 1 or P 2 P 3 P θ 3 suitable to direct methods yes yes no yes suitable to Hamilton-Jacobi-Bellman methods no yes yes yes suitable to shooting methods without constraint no no no yes provides approximations from below no no no yes This first work puts in perspective the study of necessary optimality conditions for the maximum cost problems with the help of these formulations, which will be the matter of a future work.

Viability Theory

The Pontryagin maximum principle for minimax problems of optimal control

The Bellman equation for minimizing the maximum cost

The L ∞ control problem with continuous control functions

Optimization and Nonsmooth Analysis

A numerical procedure for minimizing the maximum cost

Minimax optimal control problems. Numerical analysis of the finite horizon case

Solving minimax control problems via nonsmooth optimization

A Bellman's equation for minimizing the maximum cost

Extensions of the Dynamic Programming Framework: Battery Scheduling, Demand Charges, and Renewable Integration

A generalization of Bellman's equation with application to path planning, obstacle avoidance and invariant set estimation

Optimal, near-optimal, and robust epidemic control

An introduction to measure theory

Ordinary Differential Equations

The SIR model and the foundations of public health

The authors are grateful to Pierre Martinon for fruitful discussions and advices. This work was partially supported by ANID-PFCHA/Doctorado Nacional/2018-21180348, FONDECYT grant 1201982 and Centro de Modelamiento Matemático (CMM), ACE210010 and FB210005, BASAL funds for center of excellence, all of them from ANID (Chile)