Optimization over the Pareto front of nonconvex multi-objective optimal control problems

Simultaneous optimization of multiple objective functions results in a set of trade-off, or Pareto, solutions. Choosing a, in some sense, best solution in this set is in general a challenging task: In the case of three or more objectives the Pareto front is usually difficult to view, if not impossible, and even in the case of just two objectives constructing the whole Pareto front so as to visually inspect it might be very costly. Therefore, optimization over the Pareto (or efficient) set has been an active area of research. Although there is a wealth of literature involving finite dimensional optimization problems in this area, there is a lack of problem formulation and numerical methods for optimal control problems, except for the convex case. In this paper, we formulate the problem of optimizing over the Pareto front of nonconvex constrained and time-delayed optimal control problems as a bi-level optimization problem. Motivated by existing solution differentiability results, we propose an algorithm incorporating (i) the Chebyshev scalarization, (ii) a concept of the essential interval of weights, and (iii) the simple but effective bisection method, for optimal control problems with two objectives. We illustrate the working of the algorithm on two example problems involving an electric circuit and treatment of tuberculosis and discuss future lines of research for new computational methods.


Introduction
We continue our study of optimal control problems where one wishes to minimize simultaneously a number of conflicting objective functionals.These problems are referred to as multi-objective optimal control problems and can be expressed in the following concise form: (P) min (x,u,t f )∈X (ϕ 1 (x(t f ), t f ), . . ., ϕ r (x(t f ), t f )) .
The constraint or the feasible set X in Problem (P) involves a system of differential equations (DEs) in the state and control variables x(•) and u(•), respectively, over a time horizon [0, t f ].
The set X also typically involves point and path equality and inequality constraints.The DEs and constraints in X might even include time delays in the variables x(•) and u(•).It is Motivated by these drawbacks, minimization of an additional (single) objective function over the Pareto front has been of great interest to many researchers over the past decades-see, for example, [2,3,5,6,14,15,25,26,31,41,47].Despite this rich collection of works, to the knowledge of the authors, it was not before the reference [7] that optimization over the Pareto front was studied and a numerical method proposed for convex multi-objective optimal control problems.In the current paper, we extend the works in [7,27] to nonconvex multi-objective optimal control problems and propose a numerical method for carrying out optimization over the Pareto front.
We set the optimal control problem as a bi-level optimization problem as in [7]: One has to minimize a master objective functional subject to the minimization of a scalarization of Problem (P).The lower level problem uses the Chebyshev scalarization as in [27], as opposed to the weighted-sum scalarization in [7].The problems we consider is in much more general form in this paper: We consider nonconvex instead of convex problems compared to [7] and we consider problems with time-delay instead of those without time delay compared to [27].Just to re-iterate, [27] only proposes a technique to construct the Pareto front, otherwise it does not carry out optimization over the Pareto front.
As the optimization technique over the Pareto front, we propose the simplest possible technique, namely the bisection method, over the set of weights for the bi-objective problem, which are the parameters of the lower level optimal control problem.Even in this simplest case, it is necessary to obtain derivatives with respect to the weight, for which we employ difference approximations.However, is it guaranteed that these derivatives exist?This question is answered by [32,33,36,37] which studied the differentiability of a solution of a parametric optimal control problem with respect to the parameters.We add a discussion concerning these studies in the paper.
The main algorithm first finds the essential interval of weights over which the first step of the bisection method is taken to find a new subinterval.Then the subsequent steps of the bisection method are carried out until the stopping criterion is met.
The algorithm is illustrated on two challenging numerical examples: the Rayleigh problem, which comes from an electric circuit, and a compartmental optimal control model for tuberculosis.In the first problem there are constraints on the control variables, and the second problem not only has constraints on the two control variables but also time delays on both the control and state variables.
The paper is organized as follows.In Section 2, we introduce the multi-objective optimal control problem, discuss scalarization, introduce the problem of optimization over the Pareto front, and elaborate on solution differentiability.In Section 3, we first define and explain the essential interval of weights, and then introduce the bisection method for our problem and provide the detailed algorithm.In Section 4, we illustrate the algorithm on two example optimal control problems.Finally, in Section 5, we provide concluding remarks.

Multi-objective optimal control problem
We consider the following general multi-objective optimal control problem (similar to that in [27] but made look slightly more general here) to underlie our study on minimization over its Pareto front.The ensuing notation and definitions can also be found in [27] but given here for completeness as well as convenience. (OCP) subject to ẋ(t) = f (x(t), u(t), t) , for a.e.t ∈ [0, t f ] , θ(x(0), x(t f ), t f ) = 0 , θ(x(0), x(t f ), t f ) ≤ 0 , C(x(t), u(t), t) ≤ 0 , for a.e.t ∈ [0, t f ] , S(x(t), t) ≤ 0 , for all t ∈ [0, t f ] , where r ∈ {2, 3, 4, . ..} is fixed, the state variable x ∈ W 1,∞ (0, t f ; IR n ), ẋ := dx/dt, and the control variable u ∈ L ∞ (0, t f ; IR m ), with x(t) := (x 1 (t), . . ., x n (t)) ∈ IR n and u(t and S : IR n × IR + → IR p 4 , are continuous in their arguments.In this problem, t f is either fixed or free.Here, L ∞ (0, t f ; IR m ) corresponds to the space of essentially bounded, measurable functions equipped with the essential supremum norm.Furthermore, W 1,∞ (0, Assume that ϕ i (x(t f ), t f ) ≥ 0, for all i = 1, . . ., r.Note that this assumption can easily be met by adding a large enough positive number to each objective functional.
Note that Problem (OCP) is in general a nonsmooth problem, because it does not require differentiability of the objective functionals or the constraints.Moreover, although we have stated Problem (OCP) in very broad terms, it can further be generalized, for example by adding multi-point constraints, partial differential equations, time delays, etc.In other words, although Problem (OCP) is already in a more general form than what one usually encounters in applications, it can be further made look more general.
Of the possible extensions mentioned above, time delays in the state and control variables, for instance, can be incorporated into Problem (OCP) by replacing the ODEs in Problem (OCP) with where d x , d u > 0 are the time delays in the state and control variables, respectively.
For technical convenience, let Note that, for the case of time delays in the state and control variables, we have included Equations (1a)-(1c) instead of the ODEs ẋ(t) = f (x(t), u(t), t) in the set X.

Define the vector of objective functionals
, for all i = 1, . . ., r .
On the other hand, (x * , u * , t * f ) ∈ X is said to be a weak Pareto minimum if there exists no (x, u, t f ) ∈ X such that The set of all the Pareto and weak Pareto minima is said to be the Pareto set.On the other hand, the set of all vectors of objective functional values at the Pareto and weak Pareto minima is said to be the Pareto front (or the efficient set) of Problem (OCP) in the r-dimensional objective value, or outcome, space.Note that the coordinates of a point in the Pareto front are simply Obviously, when r = 2 the Pareto front is in general a curve; and when r = 3 the Pareto front is in general a surface.
where w i , i = 1, . . ., r, are referred to as weights, with the vector of weights w defined as w := (w 1 , . . ., w r ) ∈ IR r , such that r i=1 w i = 1.Problem (P w ) is referred to as the weighted Chebyshev problem (or Chebyshev scalarization) because of the weighted Chebyshev norm, max i |w i ϕ i (x(t f ), t f )| = max i w i ϕ i (x(t f ), t f ), appearing in the objective.This type of scalarization is typically used for nonconvex multi-objective finite-dimensional optimization problems, as opposed to the weighted sum scalarization which is effective for convex problems but not the nonconvex ones-see, for example, [38].

Define the set of weights
The following theorem was originally presented in [27,Theorem 1] for the case when there was no delay in the state and control variables.It still holds with the set X modified with the delayed state equations.
Theorem 1 (Bijection between sets of weights and Pareto minima [27]) The triplet (x * , u * , t * f ) is a weak Pareto minimum of (OCP) if, and only if, (x * , u * , t * f ) is a solution of (P w ) for some w 1 , . . ., w r > 0.
Remark 1 Suppose that Z ⊂ X denotes the Pareto set, namely the set of all Pareto minima of (OCP).Then Theorem 1 establishes that there is a bijection between the set of weights Y and the Pareto set Z.This implies that by solving (P w ) for all w ∈ Y , one can obtain the whole Pareto set Z and in turn get the Pareto front.With numerical computations on the other hand, one would of course carry out some discretization of the weight space Y and typically get a discrete approximation of the Pareto front.The bijection between Y and Z will also help us devise our algorithm for optimization over the Pareto front.✷ An ideal cost ϕ * i , i = 1, . . ., r, associated with Problem (P w ) is the optimal value of the optimal control problem, min Let (x, u, t f ) be a minimizer of the single-objective problem in (2).Then ϕ * i := ϕ i (x(t f ), t f ) and we also define ϕ j := ϕ j (x(t f ), t f ), for j = i and j = 1, . . ., r.
In the case when ϕ * i is negative, one can simply add a large enough positive number to the ith objective, to make the objective positive.In general, it is useful to add a positive number to each objective in order to obtain an even spread of the Pareto points approximating the Pareto front -see for example [21] for further discussion and geometric illustration.To serve this purpose, it is common practice to define the so-called utopian objective values.
A utopian objective vector associated with Problem (OCP) is given as β * := (β * 1 , . . ., β * r ), with β * i := ϕ * i − η i and η i > 0 for all i = 1, . . ., r. Problem (P w ) can then be equivalently written as min In the case when the objective functionals and the constraints in Problem (OCP) are differentiable in their arguments, it is worth reformulating Problem (P w ) using a standard technique from mathematical programming in the following (smooth) form.
Problem (OCP w ) is referred to as goal attainment method [38], as well as Pascoletti-Serafini scalarization [22].We will solve Problem (OCP w ) in an algorithm we present in the next section, for the two examples we want to study.
We re-iterate that the "popular" weighted-sum scalarization, given below, fails to generate the "nonconvex parts" of a Pareto front.
(P ws ) min This deficiency is illustrated with a multi-objective optimal control problem, for example, in the fed-batch bioreactor problem in [27].

Optimization over the Pareto front
The main task in this paper is to devise a numerical algorithm for solving the problem of decision making as to which Pareto point should be chosen.This obviously depends on the criterion a decision maker uses in making his/her choice.As pointed in Remark 1, the whole Pareto front can be parameterized in terms of the vector of weights w.Therefore, Problem (P w ), or equivalently (OCP w ), can be regarded as a parametric optimal control problem, and it also makes sense to express the decision maker's objective as the minimization of a function of w.
Before going ahead with the statement of this problem, we re-write the variables of the optimal control problem, with a slight abuse of notation, as x w (t) := x(t, w), u w (t) := u(t, w), and t w f := t f (w) to emphasize their dependence on the vector of weights w.We call the decision maker's objective function the master objective function, expressed by ϕ 0 (x w , u w , t w f ).With the weight vector w of the scalarization treated now as a variable, the problem of optimization over the Pareto front reduces to the problem of finding an optimal weight w * .Then the corresponding Pareto minimum is a solution of Problem (OCP w * ).
The problem of optimizing a master objective function over the Pareto front of (OCP) with r ≥ 2 objectives is nothing but a bilevel programming problem and can be written as

Remark 2
The lower-level problem in (OPF) for some given w is simply Problem (OCP w ).
A solution of (OCP w ) is nothing but a point in the Pareto set Z of (OCP) and is described by the triplet Z w := (x * (t, w), u * (t, w), t * f (w)).Then the (whole) Pareto set can be expressed as Z = ∪ w∈Y Z w .Now Problem (OPF) can equivalently be written as We note that the optimization variable of the upper-level problem is the "unknown" parameter w.If the solution (x * (t, w), u * (t, w), t * f (w)) of Problem (OCP w ) is differentiable in the parameter w, then powerful differentiable optimization techniques can be employed in solving Problem (OPF) (or in a more concise form the above problem).This is what was done in [7] for convex multi-objective optimal control problems.In this paper, we are extending the work in [7] to the nonconvex setting by also incorporating the Chebyshev scalarization and the concept of essential interval of weights given in [27].✷

Solution differentiability
We briefly review results on solution differentiability or C 1 -sensitivity of solutions to the following parametric optimal control problems depending on a parameter p ∈ P , where P is a Banach space: We note that problem (OCP w ) is a special case of the parametric problem (OCP(p)) by simply taking the parameter as the weight, p = w, which then appears only in the terminal inequality constraints.The problem (OCP(p 0 )) corresponding to a reference parameter p 0 is considered as the nominal or unperturbed problem.It is assumed that a local solution (x 0 , u 0 ) of the reference solution exists.Let p be a parameter in a neighbourhood of the nominal parameter p 0 and denote the solution to (OCP(p)) by (x(t, p), u(t, p)).Dontchev and Hager [17] gave conditions under which the mapping p → (x(•, p), u(•, p)) is Lipschitz.
Malanowski and Maurer [32,33] and Maurer and Pesch [36,37] investigated the solution differentiability or C 1 -sensitivity of the optimal solution.The authors derived conditions such that an optimal solution (x(•, p), u(•, p)) of the perturbed control problem OCP(p) exists for all parameters p in a neighborhood of p 0 and, moreover, the solution (x(t, p), u(t, p)) is a C 1 function with respect to both arguments (t, p).In broad descriptions, these conditions include certain smoothness of the functions in Problem (OCP1), satisfaction of the strict Legendre-Clebsch condition, uniqueness of the optimal control minimizing the Hamiltonian, nonsingularity of the Jacobian of an associated boundary-value problem, and boundedness of the symmetric solution of an associated Riccati ODE.
Fixing an increment d ∈ P , the differentials satisfy a linear boundary value problem that contains only information obtained in the process of computing the unperturbed solution.The computations of these sensitivity differentials can also be performed by discretization methods applied to the parametric optimal control problem; see Büskens [11] and Büskens and Maurer [12].The sensitivity differentials can be conveniently used in the minimization of a master function defined on the Pareto front; see Section 2.3.
The above mentioned conditions for showing solution differentiability exclude optimal control problems with control appearing linearly, since for this class of problems the strict Legendre-Clebsch condition does not hold.Here, optimal controls are combinations of bangbang and singular arcs.In case of finitely many switching times and junction times with the boundary of a mixed control-state constraint or a pure state constraint, one can set up a finite-dimensional optimization problem, the Induced Optimization Problem, where the switching and junction times are optimized directly; see Maurer et al. [34] and Osmolovskii and Maurer [40].If second-order sufficient conditions hold for the Induced Optimization Problem (see [40]), one immediately obtains the result that the switching and junction times locally are differentiable functions of the parameter p.
To our knowledge extensions of these results on solution differentiability to optimal control problems with control and state delays can not be found in the literature.

An Algorithm For Optimization Over the Pareto Front
As discussed in Section 2.4, the results [36, Theorem 3.1] and [37, Theorem 5.1] lay the ground for devising and implementing numerical methods for solving Problem (OPF).Bonnel and Kaya propose in [7] a barrier method for convex bi-objective optimal control problems with pure control constraints.Their method relies on twice continuous differentiability of the solution (class C 2 ) in the weight w, using the result in [36,Theorem 3.1].
In this paper, we propose a bisection method also for the case of two objectives, which relies on the solution of Problem (OCP w ) being of class C 1 w.r.t. the weight w, and thus taking the result in [37, Theorem 5.1] as a basis.Although a mathematical justification of the applicability of our proposed method, i.e., solution differentiability, is given only for Problem (OCP w ), the working of the method will also be illustrated on problems of more general class as in Problem (OCP w ).

PSfrag replacements
Pareto front Figure 1: Determination of the essential subinterval of weights In the scalarized problem (OCP w ) with two objectives (r = 2), by choosing w 1 = w, and w 2 = 1 − w, where w ∈ [0, 1], one can simply consider the single parameter w.

Essential interval of weights
With the Chebyshev scalarization, it would usually be enough for the weight w to take values over a (smaller) subinterval [w 0 , w f ] ⊂ [0, 1], with w 0 > 0 and w f < 1, for the generation of the whole front.Figure 3.1 illustrates the geometry to compute the subinterval end-points, w 0 and w f .In the illustration, the points (ϕ * 1 , ϕ 2 ) and (ϕ 1 , ϕ * 2 ) represent the boundary of the Pareto front.The equations of the "rays" which emanate from the utopia point (β * 1 , β * 2 ) and pass through the boundary points are also shown.By substituting the boundary values of the Pareto curve into the respective equations, and solving each equation for w 0 and w f one simply gets and From the geometry depicted in Figure 3.1, as also discussed in [27], one can deduce that with every w ∈ [0, w 0 ] the solution of (OCP w ) will yield the same boundary point (ϕ 1 , ϕ * 2 ) on the Pareto front.Likewise with every w ∈ [w f , 1] the same boundary point (ϕ * 1 , ϕ 2 ) is generated.This observation justifies the avoidance of the weights w ∈ [0, w 0 ) ∪ (w f , 1] in order not to keep getting the boundary points of the Pareto front, as otherwise one would end up wasting valuable computational effort and time. As a result of the above argument, the bisection method, implemented in the algorithm described in the next section, starts with the essential interval [w 0 , w f ] rather than [0, 1].It is worth re-iterating that our main concern here, unlike in [27], is not really to construct the Pareto front, but rather do a search (in this case using the bisection method) over the Pareto front, at the same time avoiding the task of constructing the front, so as to find in some sense the best solution point in the Pareto front.

Bisection method for solving Problem (OPF)
The problem of finding a best point in the Pareto front/set has now been transformed into a problem of finding best w, by virtue of the surjection from the set of weights to the set of Pareto minima furnished by Theorem 1.This has resulted in Problem (OPF) and its concise form: Find some weight w ∈ [w 0 , w f ] such that the master objective function ϕ 0 (x w , u w , t w f ) is minimized, where (x w , u w , t w f ) is found by solving (OCP w ) for that w.For a simpler setting, it is helpful to define a function F : [0, 1] → IR + representing the function we want to minimize over the Pareto front: such that (x w , u w , t w f ) solves (OCP w ).In other words, an evaluation of the function F (•) at w requires the solution of Problem (OCP w ) with that w.
Problem (OPF) can now be re-written in an even more concise form as min where F (•) is evaluated as in (4).In [7], a log-barrier method is proposed and implemented to solve (5), with an underlying convex and smooth optimal control problem with no state constraints for which the solution can be assumed to be of class C 2 , and so Newton-like methods are used with heuristic barrier parameter updates.For the general form we have in Problem (OCP), which is nonconvex and has state constraints, we assume that the solution is of class C 1 .As elaborated in Section 2.4, under certain regularity conditions which can in many cases be checked, this assumption is guaranteed to hold.Therefore we apply the bisection method [10] as an effective and simple approach to solving (5) in the case of this paper.
Albeit elementary and standard, a statement of the optimality conditions in the fact below will be useful in formulating a computational algorithm later in this section.
and, for arbitrarily small ε > 0, (b) The end point w 0 (resp.w f ) is a strict local minimizer of F (•) if, and only if, either Remark 3 (Three Cases for the End Points of [w 0 , w f ]) We will apply the bisection method starting with the essential interval [w 0 , w f ].Before introducing the pertaining algorithm, we consider below the cases for the end points of this interval.
Case I. F ′ (w 0 )F ′ (w f ) < 0 : Since F ′ (•) is assumed to be continuous the bisection method is guaranteed to find a numerical solution to (6) by the intermediate value theorem.Condition (7) needs to be check to see if w * is a strict local minimizer.
Case II.F ′ (w 0 )F ′ (w f ) > 0 : By the conditions in Fact 1(b)(i), at least one of w 0 and w f is a strict local minimizer.
Case III.F ′ (w 0 )F ′ (w f ) = 0 : If one of the inequalities in Fact 1(b)(ii) is satisfied, then w 0 or w f is a strict local minimizer.It is possible that both w 0 and w f are, or only one or neither is, a local minimizer.
In Case I, the bisection method starts with the interval [w 0 , w f ] and terminates with an approximate solution in the interior of the interval.In Case II, a local minimum is found immediately, and so in principle there is no need to do a further search.In Case III, however, the conclusion might be that neither w 0 nor w f is a strict local minimizer, in which case it would be necessary to start the bisection method with a subinterval of [w 0 , w f ], and consider Cases I-III again.✷ Remark 4 In any of the scenarios elaborated in Remark 3, consideration of another subinterval of [w 0 , w f ] might as well yield a better (lower-value) solution, since the problem is nonconvex and we can only hope to get a locally optimal solution.In our approach here, however, we do not endeavour to obtain a global minimum.As a result of our discussion in Remark 3, we will consider only Case I, which clearly prompts us to use the bisection method directly.As suggested above, in the event of Case III not yielding a solution, the new subinterval could be chosen in such a way that one would fall into Case I. ✷ The derivative of F (•) is defined at the end points of the interval [w 0 , w f ] as one-sided limits, and in the interior, i.e., for w ∈ (w 0 , w f ), as where F (•) is evaluated as in (4).In computations, we will use the forward, and backward, finite difference approximations of F ′ (•).Namely, for some small δ > 0, we will set The step δ in the difference approximation formula ( 8) is small for an accurate estimation of the derivative but not too small in order not to divide one very small number by another and cause numerical instabilities.
In what follows we provide an algorithm to solve Problem (OPF).The algorithm first finds the essential interval [w 0 , w f ], computes the signs of F ′ (w 0 ) and F ′ (w f ) and checks the cases I-III in Remark 3, and then if F ′ (w 0 )F ′ (w f ) < 0 it uses the bisection method, to find a numerical solution to Problem (OPF).
• If k = k max then declare "Maximum number of iterations exceeded."and STOP.
Step k.

Numerical Examples
In this section, we illustrate the working of Algorithm 1 on two optimal control problems, one involving an electric circuit in Section 4.1 and the other a tuberculosis (TB) epidemic in Section 4.2.
In computations, we use direct discretization of optimal control problems for which convergence theory has been an active topic of research in the literature (see for example [1,4,[18][19][20]42], and see [27] for additional references and discussion).
We employ the scalarize-discretize-then-optimize approach that was previously used in [27].Under this approach, one first scalarizes the multi-objective problem in the infinite-dimensional space, and then discretizes the scalarized problem directly and applies a usually large-scale finite-dimensional optimization method to find a discrete approximate solution of the scalarized problem.By the existing theory of discretization mentioned above, under certain assumptions, the discrete approximate solution converges to a solution of the continuous-time scalarization of the original problem, yielding a Pareto minimum of the original problem.When possible, we will also check a posteriori to see if the necessary optimality conditions are satisfied by an accurate-enough numerical solution.
In Step 0.4 of Algorithm 1, a direct discretization of Problem (OCP w ), for example employing a Runge-Kutta scheme, such as Euler's method or the Trapezoidal rule, is solved by using Ipopt, version 3.12.13,four times.In Step k.2, Problem (OCP w ) is solved in a similar way two times.Ipopt is a popular optimization software based on an interior point method; see [46].We use AMPL [23] as an optimization modelling language, which employs Ipopt as a solver.

Example: Tunnel-diode oscillator (Rayleigh problem)
The tunnel-diode oscillator problem, also referred to as the Rayleigh problem in the literature, involves dynamics represented by the following differential equations.
where the state variable x 1 (t) denotes electric current, and the control variable u(t) stands for a suitable transformation of the voltage at a generator, both at time t ∈ [0, t f ]-see [35] for a detailed exposition of the problem.In this particular instance of the problem, the initial and terminal values of the state variables are specified as and the dynamics are subject to constraints on the control variable such that The optimal control problem is posed as a bi-objective problem with min t f , where the competing objectives are the minimization of the final time t f and the minimization of the sum of the square L 2 -norms, or in some sense the magnitudes, of the current and the generator voltage.Define a new state variable x 3 such that ẋ3 (t) = x 2 1 (t) + u 2 (t) , for a.e.t ∈ [0, t f ] , x 3 (0) = 0 .
Then the two objective functionals as in Problem (OCP), or Problem (OCP w ), can be expressed as As we have stated above, the bi-objective Rayleigh problem is in the same form as Problem (OCP) and, in particular, Problem (OCP w ).The decision maker's objective for this problem will be to minimize a weighted distance to the origin of the value space.We choose where the scaling multiplier 100 is used to make the orders of magnitudes of ϕ 1 and ϕ 2 the same.We aim to solve Problem (OPF), to determine a scalar w ∈ (0, 1) with w 1 := w and w 2 := 1 − w that results in the best Pareto solution in the sense that ϕ 0 (•, •, •) is minimized, subject to the solution of Problem (OCP w ).
In [35], Maurer and Oberle numerically illustrate that an optimal solution does not exist for the single objective problem minimizing the quadratic functional ϕ 2 (x(t f ), t f ), in that t f tends to infinity.They carry out a numerical test for checking the second-order sufficient conditions (SSC) of optimality and show that the test fails to confirm the SSC.Therefore, we will impose a bound on the terminal time, namely set t f ≤ 5. On the other hand, they illustrate also in [35] that for certain instances of the weighted-sum problem, the SSC of optimality are satisfied.

Problem (OCP w ) can now explicitly be written for the Rayleigh problem as
The Hamiltonian H : IR 3 × IR × IR 3 → IR for this problem simply is , where λ(t) := (λ 1 (t), λ 2 (t), λ 3 (t)) ∈ IR 3 is referred to as the adjoint variable vector.Using the convenient notation H[t] := H(x(t), u(t), λ(t)), suppose that for all t ∈ [0, t f ], with certain transversality conditions as required by the maximum principle.In (9a)-(9c), H x i := ∂H/∂x i , i = 1, 2, 3. We will not go into the details of these (boundary) conditions here.However we note that λ 3 (t) = λ 3 , a constant, for all t ∈ [0, t f ].Then the maximum principle states that if (x, u, t f ) is an optimal solution triplet then there exists a continuous function λ(•) satisfying (9a)-(9c), along with certain transversality conditions, such that λ(t) = 0, for all t ∈ [0, t f ], and for a.e.t ∈ [0, t f ].If w = 1, then the problem is a single-objective one, referred to as a time-optimal control problem, and the condition (10) reduces to for a.e.t ∈ [0, t f ].By the discussion given in Section 3.1 (also see [27]), u w (t) given in (11) is the same for all w ∈ [w f , 1].Recall that if one does not have λ w 2 (t) = 0 for all [t ′ , t ′′ ] ⊂ [0, t f ], where t ′ < t ′′ , then u w (t) in ( 11) is referred to as optimal control of bang-bang type.We assume (and therefore will numerically double-check) that the optimal control for the particular instance of the problem is of bang-bang type.
The optimality condition (10) can be shown to yield, for any given w ∈ [w 0 , w f ), for all t ∈ [0, t f ], provided λ w 3 = 0. Again by virtue of the discussion in Section 3.1, u w (t) in ( 12) is the same for all w ∈ [0, w 0 ].We define the switching function as The constant coefficients 2 and 16 above are used for scaling purposes, so that the graphs in Figure 2(b) can be viewed more easily.Now, using (13), we can summarize and combine the expressions for the optimal control in ( 11) and ( 12) as follows.
As to why σ w (•) is referred to as the switching function should now be more clear from ( 14): the value of σ w (•) determines when to switch from one case of the control function u w (•) to another.
For Problem (OCP w ) written for the Rayleigh problem above, we have chosen the utopia vector as 2(a) depicts the Pareto front for the instance of the multi-objective Rayleigh problem we consider here.It also displays the iterations of Algorithm 1.The Rayleigh problem is discretized using the trapezoidal rule, the number of grid points is set to be N = 5000, and the Ipopt's tolerance to 10 −10 , so as to get solutions for w accurate at least up to four decimal places (dp).
If there is a need to save the computational resources further, the algorithm can be asked to yield a less accurate result, say correct to three dp, which then yields w * = 0.925 in eight iterations with (ϕ w * 1 , ϕ w * 2 ) = (3.71,45.5).In Figure 2(a) only five iterations are displayed (labels 1-5 appearing to the right of each iteration) for clarity in viewing.The Pareto (master) solution with w = w * is represented by a square.
The numerical Pareto-optimal state and control variable solutions are presented in Fig- ures 2(c)-(d) for w = w 0 , w * , w f .One of the boundary Pareto-optimal solutions is shown using solid (blue) curves for w = w 0 , which is the same solution for all w ∈ [0, w 0 ], as previously discussed in Section 3.1.On the other hand, the other boundary Pareto-optimal solution for w = w f , which holds for all w ∈ [w f , 1], is shown using dashed (green) curves.The latter is nothing but a time-optimal control solution for the Rayleigh problem (a solution with the smallest t f ), resulting in a bang-bang type function with the sequence of values {1, −1, 1}, namely with two switchings.The master Pareto solution is given for w = w * using dashed-and-dotted (red) curves.
The switching function σ w (•) plotted in Figure 2(b) by using (13) (recall that discrete approximations of λ w 2 (t) and λ w 3 (t) can readily be obtained from AMPL) furnishes the means    to verify the optimality condition for u w (•) expressed in (14).It is evident from the dashed (green) plot of the switching function that, for w ∈ [w f , 1], when σ w (•) crosses the time axis there is a jump (from 1 to −1 or vice versa) in the value of the corresponding u w (•) plot.Likewise, for w ∈ [0, w 0 ] and for w = w * ∈ [w 0 , w f ), whenever σ w (•) crosses one of the lines σ w (t) = 1 and σ w (t) = −1 (shown by two black lines in Figure 2(b) for convenience) the expression for the control function u w (•) switches from one case in (14) to another, as required.

Example: Compartmental model for tuberculosis
In 2020 and 2021, tuberculosis (TB) was the second leading cause of death from an infectious disease worldwide after COVID-19 [44].Active TB refers to disease that occurs in someone infected with Mycobacterium tuberculosis.It is characterized by signs or symptoms of active disease, or both, and is distinct from latent tuberculosis infection, which occurs without signs or symptoms of active disease.Only individuals with active TB can transmit the infection.Many people with active TB do not experience typical TB symptoms in the early stages of the disease.These individuals are unlikely to seek care early, and may not be properly diagnosed when seeking care.Delays to diagnosis of active TB present a major obstacle to the control of a TB epidemic, it may worsen the disease, increase the risk of death and enhance tuberculosis transmission to the community.Both patient and the health system may be responsible for the diagnosis delay.
We study the control model with control and state delays presented in Silva et al. [43].In this model, reinfection and post-exposure interventions for tuberculosis are considered.The population is divided into five categories (compartments) (i.e., the control system has five state variables): The dynamical system is given by The recovered population is defined by with N = 30000.The system and delay parameters in the model (15) along with their values are listed in Table 1.In view of the delays the initial conditions and functions are: The control constraints are given by We consider the following parametric objective functional with control weights a 1 , a 2 ≥ 0: Denoting the (augmented) state vector by x(t) = (S(t), L 1 (t), I(t), L 2 (t), x 5 (t), x 6 (t)) ∈ R 6 and the control vector u(t) := (u 1 (t), u 2 (t)) ∈ R 2 , the two competing objectives in the general problem (P) are given by where F 1 (x, u) and F 2 (x, u) denote the two functionals in Lagrange form.
The bi-objective TB problem is now in the same form as Problem (OCP) and, in particular, Problem (OCPsd).The decision maker's objective for this problem will be to minimize the distance to the origin of the value space.We therefore choose Our aim is to solve Problem (OPF), to determine a scalar w ∈ (0, 1) with w 1 := w and w 2 := 1 − w that results in the best Pareto solution in the sense that ϕ 0 (•, •, •) is minimized, subject to the solution of Problem (OCP w ).
Next we focus on the solution of Problem (OCP w ): We aim to find a pair of functions ( ) that minimizes the parameter α subject to the time-delayed dynamics (15) and the auxiliary dynamics (21), initial conditions (17), control constraints (18) and auxiliary weighted inequalities involving ϕ 1 and ϕ 2 .
We consider the necessary optimality conditions for the time-delayed optimal control problem (OCP w ); see Göllmann and Maurer [24], Vinter [45].For this purpose we introduce the delayed state variable y 3 Denoting the adjoint variable vector by λ(t) := (λ S (t), λ L 1 (t), λ I (t), λ L 2 (t), λ 5 (t), λ 6 (t)) ∈ R 6 the Hamiltonian or Pontryagin function is given by where R is given as in (16).The Minimum Principle [24,45] where the argument [t] stands for evaluating all arguments at time t.We note that λ w 5 (t) = λ As in the Rayleigh problem, the superscript "w" above denotes dependence on the scalarization parameter/weight w.Then the controls minimizing the Hamiltonian are characterized by the switching conditions (control law) In what follows we choose the control weights as a 11 = a 12 = 10 (small) and a 21 = a 22 = 1000 (large) in the objective functionals ϕ 1 and ϕ 2 .
For Problem (OCP w ) written for the TB problem, we have chosen the utopia vector as (β * 1 , β * 2 ) = (0, 0). Figure 3 depicts the Pareto front for the TB problem we consider here.The plot also displays the iterations of Algorithm 1.The TB problem is discretized using the trapezoidal rule, the number of grid points is set to be N = 5000, and the Ipopt's tolerance to 10 −10 , so as to get solutions for w accurate at least up to four decimal places (dp).
The numerical Pareto-optimal control variable solutions u w 1 (•) and u w 2 (•) are presented in Figures 4(a)-(b) for w = w 0 , w * , w f .As with Rayleigh, one of the boundary Paretooptimal solutions is shown using solid (blue) curves for w = w 0 , the same solution for all w ∈ [0, w 0 ].The other boundary Pareto-optimal solution for w = w f , which holds for all w ∈ [w f , 1], is shown using dashed (green) curves.Both of the control solutions are of bangbang type (as required by ( 24)), with one switching (the number of switchings not dictated   by (24) alone).The master Pareto solution is given for w = w * using dashed-and-dotted (red) curves, in which the controls are also of bang-bang type with one switching.
The switching functions for each control and case, σ w k (•), k = 1, 2, scaled as indicated, are plotted with (black) dotted curves and superposed with the control plots in Figures 4(a)-(b).We remind that, by using (23) (recall that discrete approximations of λ w L k (t), k = 1, 2, λ w 5 (t) and λ w 6 (t) can readily be obtained as constraint multipliers from AMPL), one verifies the optimality condition in (24).
In each strategy, the two control efforts are "on" until the times t w s k , k = 1, 2, at which the respective u w k (•) is switched "off" (down to zero).These types of bang-bang controls are also referred to as on-off controls.In Table 2 the switching times for the boundary as well as the optimal weights are listed.Under these controls, the resulting terminal values of the state variables are also listed in Table 2.The plots of these variables are not provided as they are difficult to distinguish at earlier times (as expected) and that they become distinguishable/comparable only near the terminal time.
Under the controls minimizing x 5 (t f ) (with w = w f = 0.5709 and minimum x w f 5 (t f ) = 26459) the number of persistent latent individuals L 2 (t f ) turns out to be about 419 (in a population of 30000).This number is more than doubled to 864 if x 6 (t f ) is minimized (with w = w 0 = 0.5709 and minimum x w 0 6 (t f ) = 31133).The optimal Pareto solution minimizing the distance in value space to the origin yields with w = w * = 0.5358 the optimal L 2 (t f ) as 748.

Conclusion
We have proposed an algorithm to solve the problem of optimization over the Pareto front.The algorithm employs bisection method which starts with an essential interval of weights of the Chebyshev scalarization.It is applicable to a wide range of optimal control problems, including state-and control-constrained problems with time delay.Numerical solution of two challenging optimal control problems has demonstrated the effectiveness of the algorithm.
The main motive behind the algorithm we have proposed is that one can find the optimal solution minimizing a master objective functional without having to construct the Pareto front.The algorithm solves the challenging optimal control problem (OCP w ) a relatively smaller number of times than the case of constructing the Pareto front.In the examples we have studied the algorithm had to solve (OCP w ) 20 to 30 times.On the other hand, without the algorithm we propose, it is necessary to construct the Pareto front by solving (OCP w ) thousands of times in order to obtain the same solution with the same computational accuracy.

Fact 1
Consider the minimization problem in (5) with F (•) of class C 1 .(a) The interior point w * ∈ (w 0 , w f ) is a strict local minimizer of F (•) if, and only if, 0 and neither of the inequalities in Fact 1(b)(ii) is satisfied then declare "Algorithm failed.Change the interval [w 0 , w f ]." and STOP.Let a := w 0 and b := w f .Step k.1 (Bisection) Find the midpoint c := a + (b − a)/2 of the interval [a, b].
Pareto front, and iterations of Algorithm 1: Master solution is depicted by a (red) square and iterates by (light blue) circles.Switching function as defined in(13).

Figure 2 :
Figure 2: Rayleigh problem-Boundary Pareto solutions, corresponding to w 0 = 0.8994 and w f = 0.9269, are shown with (blue) solid curves and (green) dashed curves, respectively.Master Pareto solution, corresponding to w * = 0.9247, is shown with dashed-and-dotted (red) curves.

S : susceptible individuals, L 1 :
early latent individuals, recently infected (less than two years), I : infectious individuals, who have active TB, L 2 : persistent latent individuals, R : recovered individuals, N : total population N = S + L 1 + I + L 2 + R , assumed constant.The model has two control variables and three delays: u 1 : effort on early detection and treatment of recently infected individuals L 1 , d u 1 : delay on the diagnosis of latent TB, and commencement of latent TB treatment, u 2 : chemotherapy or post-exposure vaccine to persistent latent individuals L 2 , d u 2 : delay in the prophylactic treatment of persistent latent L 2 , d I : delay in I, i.e., delay in diagnosis.

w 5 and λ w 6 (t) = λ w 5 ,
constants, for any fixed w ∈ [0, 1].In the last equation, the term χ [0,t f −d I ] (t) denotes the characteristic function of the interval [0, t f − d I ] at time t.The minimization of the Hamiltonian with respect to the controls u 1 , u 2 and delayed controls v 1 , v 2 involves the switching functions σ k (t) for k = 1, 2:

Figure 3 :
Figure 3: TB problem-Pareto front, and iterations of Algorithm 1: Master solution is depicted by a (red) square and iterates by (light blue) circles..

Figure 4 :
Figure 4: TB problem-Boundary Pareto solutions, corresponding to w 0 = 0.5251 and w f = 0.5709, are shown with (blue) solid curves and (green) dashed curves, respectively.Master Pareto solution, corresponding to w * = 0.5358, is shown with dashed-and-dotted (red) curves.

Table 1 :
Parameter values for the TB control model.Depending on the priorities, the weights a 1 , a 2 can be chosen in different ways (for example, both can be chosen to be very small or very large) giving rise to competing objectives.Namely, 11 , a 12 , a 21 , a 22 ≥ 0, constitute two competing objective functionals.Both functionals are given in Lagrange form.The standard method to obtain an optimal control problem of Bolza type is to introduce additional state variables x 5 and x 6 defined by ẋ5 (t) = I(t) + L 2 (t) + a 11 u 1 (t) + a 12 u 2 (t) , x 5 (0) = 0 , ẋ6 (t) = I(t) + L 2 (t) + a 21 u 1 (t) + a 22 u 2 (t) , x 6 (0) = 0 .
w (t f ) L w 1 (t f ) I w (t f ) L w 2 (t f ) R w (t f )