# Optimization Over the Pareto Front of Nonconvex Multi-objective Optimal
Control Problems

C. Yalçın Kaya111Mathematics, UniSA STEM, University of South Australia,
Mawson Lakes, S.A. 5095, Australia. E-mail<EMAIL_ADDRESS>. Helmut
Maurer222Institut für Numerische und Angewandte Mathematik, Westfälische
Wilhelms-Universität Münster, Münster, Germany. E-mail: helmut.maurer@uni-
muenster.de .

###### Abstract

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.

> Key words: Multi-objective optimization, Optimal control, Optimization over
> Pareto front, Optimization over efficient set, Numerical methods, Rayleigh
> problem, Tuberculosis, Time-delay problems.

##  1 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:

$\mbox{(P)}\hskip 28.45274pt\min_{(x,u,t_{f})\in X}\
(\varphi_{1}(x(t_{f}),t_{f}),\ldots,\varphi_{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(\cdot)$ and
$u(\cdot)$, 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(\cdot)$ and $u(\cdot)$. It is worth noting that although each of the
objective functionals $\varphi_{i}(x(t_{f}),t_{f})$, $i=1,\ldots,r$, in (P)
above constitutes the so-called Mayer form, other forms (Bolza and Lagrange)
can easily be converted into this form conveniently. Therefore, the general
model in (P) caters for a wide range of conflicting objectives; for instance,
minimization of the energy, the terminal time, the deviations from a reference
state trajectory, or the uncertainty in measurements, to name just a few.

Broadly speaking, the simultaneous or Pareto minimization in Problem (P) is
the process of finding a compromise solution, referred to as a Pareto minimum,
where the value of some cost cannot be improved (i.e., reduced) further,
without making the value of some other cost worse (i.e., higher). One typical
example is the case when one wants to minimize simultaneously the fuel
expenditure of an airplane travelling from one given city to another and the
time the airplane takes for this travel: A shorter travel time often requires
a higher fuel consumption. The set of all such compromise or trade-off
solutions form the Pareto set in the optimization space, or the Pareto front
in the value space. Pareto set and Pareto front are also commonly referred to
as the efficient set and the efficient front, respectively333These and other
definitions will be given in more precise terms in Section 2..

The authors of this paper have studied in [27] the problem of constructing the
Pareto front of Problem (P) involving ODEs and constraints of general form.
They discussed and demonstrated that for the nonconvex optimal control
problems like the one in Problem (P), it is better to use the so-called
weighted Chebyshev-norm scalarization (or just Chebyshev scalarization) to
guarantee that the whole Pareto front can be constructed, instead of using the
traditional weighted-sum scalarization, i.e., a convex combination of the
objective functionals. They discretized the scalarized problem directly and
utilized large-scale optimization software (the AMPL–Ipopt suite [23, 46]) to
find the Pareto fronts of two constrained optimal control problems as
examples, one involving tumour anti-angiogenesis and the other a fed-batch
bioreactor, by means of what they called a scalarize–discretize–then–optimize
approach. This approach is in contrast with the other existing
discretize–scalarize–then–optimize approach (see e.g. [28, 29, 30, 39]) which
scalarizes the discretized problem rather than the original (continuous-time)
problem.

An additional benefit of the Chebyshev scalarization is also reported and
illustrated in [27]: One can compute the whole Pareto front by using only
those weights of the objective functionals within what they name as the
essential subinterval of weights, instead of the whole interval. Having to
compute fewer Pareto solutions over a smaller number of grid points in a
subinterval is obviously a computational advantage. For further details and an
extensive list of references on multi-objective optimal control the reader is
referred to [27]. Other relevant studies on the topic in more recent years
have appeared in [16, 13].

Apart from certain trivial or special cases, the Pareto front consists of
infinitely many solutions to choose from. When a discrete approximation of the
front is found the number of solutions to choose from is still relatively
large since the approximate front is required to be accurate enough. Making a
decision as to which Pareto solution in the front is the most suitable (to the
needs of a practitioner) is often very hard for the following reasons.

  * •

In the case of three or more objectives, the Pareto front might be difficult
(if not impossible) to view and to carry out a visual inspection (or
“eyeballing”) for a decision.

  * •

Even with two objectives, a visual inspection alone may not be enough to
choose a desirable solution.

  * •

Constructing the whole Pareto front might just be too costly a thing to do
numerically.

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.

##  2 Problem Statement and Preliminaries

###  2.1 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.

$\mbox{(OCP) }\left\\{\begin{array}[]{rl}\displaystyle\min&\ \
\displaystyle(\varphi_{1}(x(t_{f}),t_{f}),\ldots,\varphi_{r}(x(t_{f}),t_{f}))\\\\[8.53581pt]
\mbox{subject to}&\ \ \dot{x}(t)=f(x(t),u(t),t)\,,\mbox{\ \ for a.e.
}t\in[0,t_{f}]\,,\\\\[5.69054pt] &\ \
\theta(x(0),x(t_{f}),t_{f})=0\,,\\\\[5.69054pt] &\ \
\widetilde{\theta}(x(0),x(t_{f}),t_{f})\leq 0\,,\\\\[5.69054pt] &\ \
C(x(t),u(t),t)\leq 0\,,\mbox{\ \ for a.e. }t\in[0,t_{f}]\,,\\\\[5.69054pt] &\
\ S(x(t),t)\leq 0\,,\mbox{\ \ for all }t\in[0,t_{f}]\,,\end{array}\right.$

where $r\in\\{2,3,4,\ldots\\}$ is fixed, the state variable $x\in
W^{1,\infty}(0,t_{f};{\rm{I\ \kern-5.39993ptR}}^{n})$, $\dot{x}:=dx/dt$, and
the control variable $u\in L^{\infty}(0,t_{f};{\rm{I\
\kern-5.39993ptR}}^{m})$, with $x(t):=(x_{1}(t),\ldots,x_{n}(t))\in{\rm{I\
\kern-5.39993ptR}}^{n}$ and $u(t):=(u_{1}(t),\ldots,u_{m}(t))\in{\rm{I\
\kern-5.39993ptR}}^{m}$. The functions $\varphi_{i}:{\rm{I\
\kern-5.39993ptR}}^{n}\times{\rm{I\ \kern-5.39993ptR}}_{+}\rightarrow{\rm{I\
\kern-5.39993ptR}}$, $f:{\rm{I\ \kern-5.39993ptR}}^{n}\times{\rm{I\
\kern-5.39993ptR}}^{m}\times{\rm{I\ \kern-5.39993ptR}}_{+}\rightarrow{\rm{I\
\kern-5.39993ptR}}^{n}$, $\theta:{\rm{I\ \kern-5.39993ptR}}^{n}\times{\rm{I\
\kern-5.39993ptR}}^{n}\times{\rm{I\ \kern-5.39993ptR}}_{+}\rightarrow{\rm{I\
\kern-5.39993ptR}}^{p_{1}}$, $\widetilde{\theta}:{\rm{I\
\kern-5.39993ptR}}^{n}\times{\rm{I\ \kern-5.39993ptR}}^{n}\times{\rm{I\
\kern-5.39993ptR}}_{+}\rightarrow{\rm{I\ \kern-5.39993ptR}}^{p_{2}}$,
$C:{\rm{I\ \kern-5.39993ptR}}^{n}\times{\rm{I\
\kern-5.39993ptR}}^{m}\times{\rm{I\ \kern-5.39993ptR}}_{+}\rightarrow{\rm{I\
\kern-5.39993ptR}}^{p_{3}}$, and $S:{\rm{I\
\kern-5.39993ptR}}^{n}\times{\rm{I\ \kern-5.39993ptR}}_{+}\rightarrow{\rm{I\
\kern-5.39993ptR}}^{p_{4}}$, are continuous in their arguments. In this
problem, $t_{f}$ is either fixed or free. Here, $L^{\infty}(0,t_{f};{\rm{I\
\kern-5.39993ptR}}^{m})$ corresponds to the space of essentially bounded,
measurable functions equipped with the essential supremum norm. Furthermore,
$W^{1,\infty}(0,t_{f};\mathbb{R}^{n})$ is the Sobolev space consisting of
functions $x:[0,t_{f}]\rightarrow\mathbb{R}^{n}$ whose first derivatives lie
in $L^{\infty}$.

Assume that $\varphi_{i}(x(t_{f}),t_{f})\geq 0$, for all $i=1,\ldots,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

$\displaystyle\dot{x}(t)=f(x(t),x(t-d_{x}),u(t),u(t-d_{u}),t)\,,\mbox{\ \ for
a.e. }t\in[0,t_{f}]\,,$ (1a) $\displaystyle x(t)=x_{0}(t)\,,\mbox{\ \ for all
}t\in[-d_{x},0)\,,$ (1b) $\displaystyle u(t)=u_{0}(t)\,,\mbox{\ \ for all
}t\in[-d_{u},0)\,,$ (1c)

where $d_{x},d_{u}>0$ are the time delays in the state and control variables,
respectively.

For technical convenience, let $t_{f}\leq t_{f}^{\max}$, where
$t_{f}^{\max}>0$ is some constant. Next, we define the feasible set, $X\subset
W^{1,\infty}(0,t_{f};{\rm{I\ \kern-5.39993ptR}}^{n})\times
L^{\infty}(0,t_{f};{\rm{I\ \kern-5.39993ptR}}^{m})\times{\rm{I\
\kern-5.39993ptR}}_{+}$, such that

$\displaystyle X$ $\displaystyle:=$
$\displaystyle\\{(x,u,t_{f})\,:\,\dot{x}(t)=f(x(t),x(t-d_{x}),u(t),u(t-d_{u}),t)\,,\mbox{\
\ for a.e. }t\in[0,t_{f}]\,;$ $\displaystyle\quad x(t)=x_{0}(t)\,,\mbox{\ \
for all }t\in[-d_{x},0];\ \ u(t)=u_{0}(t)\,,\mbox{\ \ for all
}t\in[-d_{u},0)\,;$ $\displaystyle\quad\theta(x(0),x(t_{f}),t_{f})=0\,;\
\widetilde{\theta}(x(0),x(t_{f}),t_{f})\leq 0\,;$ $\displaystyle\quad
C(x(t),u(t),t)\leq 0\,,\mbox{ for a.e. }t\in[0,t_{f}];\ S(x(t),t)\leq 0,\mbox{
for all }t\in[0,t_{f}]\\}\,.$

Note that, for the case of time delays in the state and control variables, we
have included Equations (1a)–(1c) instead of the ODEs
$\dot{x}(t)=f(x(t),u(t),t)$ in the set $X$.

Define the vector of objective functionals,
$\varphi(x(t_{f}),t_{f}):=(\varphi_{1}(x(t_{f}),t_{f}),\ldots,\varphi_{r}(x(t_{f}),t_{f}))$.
The triplet $(x^{*},u^{*},t_{f}^{*})\in X$ is said to be a Pareto minimum if
there exists no $(x,u,t_{f})\in X$ such that
$\varphi(x(t_{f}),t_{f})\neq\varphi(x^{*}(t_{f}^{*}),t_{f}^{*})$ and

$\varphi_{i}(x(t_{f}),t_{f}))\leq\varphi_{i}(x^{*}(t_{f}^{*}),t_{f}^{*})\,,\quad\mbox{for
all }i=1,\ldots,r\,.$

On the other hand, $(x^{*},u^{*},t_{f}^{*})\in X$ is said to be a weak Pareto
minimum if there exists no $(x,u,t_{f})\in X$ such that

$\varphi_{i}(x(t_{f}),t_{f}))<\varphi_{i}(x^{*}(t_{f}^{*}),t_{f}^{*})\,,\quad\mbox{for
all }i=1,\ldots,r\,.$

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$-dimensionalobjective value, or
outcome, space. Note that the coordinates of a point in the Pareto front are
simply $\varphi_{i}(x^{*}(t_{f}^{*}),t_{f}^{*})$, $i=1,\ldots,r$. Obviously,
when $r=2$ the Pareto front is in general a curve; and when $r=3$ the Pareto
front is in general a surface.

###  2.2 Scalarization

In [27], to compute a solution of Problem (OCP), the following single-
objective problem (Pw), i.e., scalarization, was employed.

$\mbox{(P${}_{w}$)}\qquad\min_{(x,u,t_{f})\in X}\
\max\\{w_{1}\,\varphi_{1}(x(t_{f}),t_{f}),\ldots,w_{r}\,\varphi_{r}(x(t_{f}),t_{f})\\}\,,$

where $w_{i}$, $i=1,\ldots,r$, are referred to as weights, with the vector of
weights $w$ defined as $w:=(w_{1},\ldots,w_{r})\in{\rm{I\
\kern-5.39993ptR}}^{r}$, such that $\sum_{i=1}^{r}w_{i}=1$. Problem (Pw) is
referred to as the weighted Chebyshev problem (or Chebyshev scalarization)
because of the weighted Chebyshev norm,
$\max_{i}|w_{i}\,\varphi_{i}(x(t_{f}),t_{f})|=\max_{i}w_{i}\,\varphi_{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

$Y:=\left\\{w\in{\rm{I\ \kern-5.39993ptR}}^{r}\ |\
\sum_{i=1}^{r}w_{i}=1\right\\}\,.$

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 (Pw) for some
$w_{1},\ldots,w_{r}>0$.

###### Remark 1

Suppose that $Z\subset 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
(Pw) for all $w\in 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. $\Box$

An ideal cost $\varphi^{*}_{i}$, $i=1,\ldots,r$, associated with Problem (Pw)
is the optimal value of the optimal control problem,

$\min_{(x,u,t_{f})\in X}\ \varphi_{i}(x(t_{f}),t_{f})\,.$ (2)

Let $(\overline{x},\overline{u},\overline{t}_{f})$ be a minimizer of the
single-objective problem in (2). Then
$\varphi^{*}_{i}:=\varphi_{i}(\overline{x}(\overline{t}_{f}),\overline{t}_{f})$
and we also define
$\overline{\varphi}_{j}:=\varphi_{j}(\overline{x}(\overline{t}_{f}),\overline{t}_{f})$,
for $j\neq i$ and $j=1,\ldots,r$.

In the case when $\varphi^{*}_{i}$ is negative, one can simply add a large
enough positive number to the $i$th 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
$\beta^{*}:=(\beta_{1}^{*},\ldots,\beta_{r}^{*})$, with
$\beta^{*}_{i}:=\varphi_{i}^{*}-\eta_{i}$ and $\eta_{i}>0$ for all
$i=1,\ldots,r$. Problem (Pw) can then be equivalently written as

$\min_{(x,u,t_{f})\in X}\
\max\\{w_{1}\,(\varphi_{1}(x(t_{f}),t_{f})-\beta_{1}^{*}),\ldots,w_{r}\,(\varphi_{r}(x(t_{f}),t_{f})-\beta_{r}^{*})\\}\,.$

In the case when the objective functionals and the constraints in Problem
(OCP) are differentiable in their arguments, it is worth reformulating Problem
(Pw) using a standard technique from mathematical programming in the following
(smooth) form.

$\mbox{(OCP${}_{w}$)
}\left\\{\begin{array}[]{rl}\displaystyle\min_{{\alpha\geq
0}\atop{(x,u,t_{f})\in X}}&\ \ \alpha\\\\[11.38109pt] \mbox{subject to}&\ \
w_{1}\,(\varphi_{1}(x(t_{f}),t_{f})-\beta_{1}^{*})\leq\alpha\,,\\\ &\hskip
68.2866pt\vdots\\\ &\ \
w_{r}\,(\varphi_{r}(x(t_{f}),t_{f})-\beta_{r}^{*})\leq\alpha\,.\end{array}\right.$

Problem (OCPw) is referred to as goal attainment method [38], as well as
Pascoletti-Serafini scalarization [22]. We will solve Problem (OCPw) 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.

$\mbox{(P${}_{ws}$)}\qquad\min_{(x,u,t_{f})\in X}\
\sum_{i=1}^{r}w_{i}\,\varphi_{i}(x(t_{f}),t_{f})\,.$

This deficiency is illustrated with a multi-objective optimal control problem,
for example, in the fed-batch bioreactor problem in [27].

###  2.3 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 (Pw), or equivalently
(OCPw), 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_{f}^{w}:=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 $\varphi_{0}(x^{w},u^{w},t_{f}^{w})$. 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\geq 2$ objectives is nothing but a bilevel programming problem
and can be written as

$\mbox{(OPF) }\left\\{\begin{array}[]{rll}\displaystyle\min_{w\in Y}&\ \
\varphi_{0}(x^{w},u^{w},t_{f}^{w})&\\\\[5.69054pt] \ \ \mbox{subject to}&\ \
\displaystyle\min_{{\alpha\geq 0}\atop{(x,u,t_{f})\in
X}}&\alpha\\\\[11.38109pt] &\ \ \mbox{subject
to}&w_{1}\,(\varphi_{1}(x(t_{f},w),t_{f})-\beta_{1}^{*})\leq\alpha\,,\\\
&&\hskip 42.67912pt\vdots\\\
&&w_{r}\,(\varphi_{r}(x(t_{f},w),t_{f})-\beta_{r}^{*})\leq\alpha\,.\end{array}\right.$

###### Remark 2

The lower-level problem in (OPF) for some given $w$ is simply Problem (OCPw).
A solution of (OCPw) 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=\cup_{w\in Y}Z_{w}$. Now
Problem (OPF) can equivalently be written as

$\left\\{\begin{array}[]{rll}\displaystyle\min_{w\in Y}&\ \
\varphi_{0}(x^{w},u^{w},t_{f}^{w})&\\\\[5.69054pt] \ \ \mbox{subject to}&\ \
(x^{w},u^{w},t_{f}^{w})\in Z_{w}\,.\end{array}\right.$

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 (OCPw) 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]. $\Box$

###  2.4 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\in P$, where $P$ is a Banach space:

$\mbox{(OCP(p))}\left\\{\begin{array}[]{rl}\min_{\,x,u,p}&g(x(t_{f}),t_{f},p)\\\\[5.69054pt]
\mbox{subject to}&\dot{x}(t)=\widetilde{f}(x(t),u(t),p)\,,\mbox{\ \ for a.e.
}t\in[0,t_{f}]\,,\\\\[2.84526pt]
&\psi(x(0),x(t_{f}),t_{f},p)=0\,,\quad\tilde{\psi}(x(0),x(t_{f}),t_{f},p)\leq
0\,,\\\\[2.84526pt] &\tilde{C}(x(t),u(t),p)\leq 0\,,\mbox{\ \ for a.e.
}\;t\in[0,t_{f}]\,,\\\\[2.84526pt] &\tilde{S}(x(t),p)\leq 0\,,\mbox{\ \ for
a.e. }\;t\in[0,t_{f}]\,.\end{array}\right.$

We note that problem (OCPw) 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\mapsto(x(\cdot,p),u(\cdot,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(\cdot,p),u(\cdot,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\in P$, the differentials

$z_{d}(t,p_{0})=\frac{\partial x}{\partial p}(t,p_{0})d,\quad
v_{d}(t,p_{0})=\frac{\partial u}{\partial p}(t,p_{0})d,$

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 bang-bang 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.

##  3 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 (OCPw) 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
(OCPw), the working of the method will also be illustrated on problems of more
general class as in Problem (OCPw).

In the scalarized problem (OCPw) with two objectives ($r=2$), by choosing
$w_{1}=w$, and $w_{2}=1-w$, where $w\in[0,1]$, one can simply consider the
single parameter $w$.

###  3.1 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}]\subset[0,1]$,
with $w_{0}>0$ and $w_{f}<1$, for the generation of the whole front. Figure 1
illustrates the geometry to compute the subinterval end-points, $w_{0}$ and
$w_{f}$. In the illustration, the points
$(\varphi_{1}^{*},\overline{\varphi}_{2})$ and
$(\overline{\varphi}_{1},\varphi_{2}^{*})$ represent the boundary of the
Pareto front. The equations of the “rays” which emanate from the utopia point
$(\beta_{1}^{*},\beta_{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

$w_{0}=\frac{(\varphi_{2}^{*}-\beta_{2}^{*})}{(\overline{\varphi}_{1}-\beta_{1}^{*})+(\varphi_{2}^{*}-\beta_{2}^{*})}\qquad\mbox{and}\qquad
w_{f}=\frac{(\overline{\varphi}_{2}-\beta_{2}^{*})}{(\varphi_{1}^{*}-\beta_{1}^{*})+(\overline{\varphi}_{2}-\beta_{2}^{*})}\,.$
(3)

Figure 1: Determination of the essential subinterval of weights
$[w_{0},w_{f}]$ [27].

From the geometry depicted in Figure 1, as also discussed in [27], one can
deduce that with every $w\in[0,w_{0}]$ the solution of (OCPw) will yield the
same boundary point $(\overline{\varphi}_{1},\varphi_{2}^{*})$ on the Pareto
front. Likewise with every $w\in[w_{f},1]$ the same boundary point
$(\varphi_{1}^{*},\overline{\varphi}_{2})$ is generated. This observation
justifies the avoidance of the weights $w\in[0,w_{0})\cup(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.

###  3.2 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\in[w_{0},w_{f}]$ such that the master objective function
$\varphi_{0}(x^{w},u^{w},t_{f}^{w})$ is minimized, where
$(x^{w},u^{w},t_{f}^{w})$ is found by solving (OCPw) for that $w$. For a
simpler setting, it is helpful to define a function $F:[0,1]\to{\rm{I\
\kern-5.39993ptR}}_{+}$ representing the function we want to minimize over the
Pareto front:

$F(w):=\varphi_{0}(x^{w},u^{w},t_{f}^{w})\,,$ (4)

such that $(x^{w},u^{w},t_{f}^{w})$ solves (OCPw). In other words, an
evaluation of the function $F(\cdot)$ at $w$ requires the solution of Problem
(OCPw) with that $w$.

Problem (OPF) can now be re-written in an even more concise form as

$\min_{w\in[w_{0},w_{f}]}\ F(w)\,,$ (5)

where $F(\cdot)$ 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.

###### Fact 1

Consider the minimization problem in (5) with $F(\cdot)$ of class $C^{1}$.

  * (a)

The interior point $w^{*}\in(w_{0},w_{f})$ is a strict local minimizer of
$F(\cdot)$ if, and only if,

$F^{\prime}(w^{*})=0\,,$ (6)

and, for arbitrarily small $\varepsilon>0$,

$F^{\prime}(w^{*}-\varepsilon)<0\quad\mbox{and}\quad
F^{\prime}(w^{*}+\varepsilon)>0\,.$ (7)

  * (b)

The end point $w_{0}$ (resp. $w_{f}$) is a strict local minimizer of
$F(\cdot)$ if, and only if, either

    * (i)

$F^{\prime}(w_{0})>0$ (resp. $F^{\prime}(w_{f})<0$) or

    * (ii)

$F^{\prime}(w_{0})=0$ (resp. $F^{\prime}(w_{f})=0$) and, for arbitrarily small
$\varepsilon>0$, $F^{\prime}(w_{0}+\varepsilon)>0$ (resp.
$F^{\prime}(w_{f}-\varepsilon)<0$).

######  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^{\prime}(w_{0})F^{\prime}(w_{f})<0$ : Since $F^{\prime}(\cdot)$ 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^{\prime}(w_{0})F^{\prime}(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^{\prime}(w_{0})F^{\prime}(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. $\Box$

###### 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.
$\Box$

The derivative of $F(\cdot)$ is defined at the end points of the interval
$[w_{0},w_{f}]$ as one-sided limits,

$F^{\prime}(w_{0}):=\lim_{\delta\to
0^{+}}\frac{F(w_{0}+\delta)-F(w_{0})}{\delta}\quad\mbox{and}\quad
F^{\prime}(w_{f}):=\lim_{\delta\to
0^{-}}\frac{F(w_{f}+\delta)-F(w_{f})}{\delta}\,,$

and in the interior, i.e., for $w\in(w_{0},w_{f})$, as

$F^{\prime}(w):=\lim_{\delta\to 0}\frac{F(w+\delta)-F(w)}{\delta}\,,$

where $F(\cdot)$ is evaluated as in (4). In computations, we will use the
forward, and backward, finite difference approximations of
$F^{\prime}(\cdot)$. Namely, for some small $\delta>0$, we will set

$F^{\prime}(w)\approx\left\\{\begin{array}[]{rl}\dfrac{F(w+\delta)-F(w)}{\delta}\,,&\
\mbox{ if }w\in[w_{0},w_{f}-\delta)\,,\\\\[8.53581pt]
\dfrac{F(w)-F(w-\delta)}{\delta}\,,&\ \mbox{ if
}w\in[w_{f}-\delta,w_{f}]\,.\end{array}\right.$ (8)

The step $\delta$ 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^{\prime}(w_{0})$ and $F^{\prime}(w_{f})$ and checks the cases I–III in
Remark 3, and then if $F^{\prime}(w_{0})F^{\prime}(w_{f})<0$ it uses the
bisection method, to find a numerical solution to Problem (OPF).

###### Algorithm 1

Step $0.0$

(Initialization) Choose utopia parameters, $\eta_{1},\eta_{2}>0$, a small
numerical differentiation step $\delta>0$, a stopping tolerance $\epsilon>0$,
and a maximum number of iterations $k_{\max}$ . Set $k:=1$.

Step $0.1$

(Boundary points of the front) Solve (2) to get
$(\overline{x}^{i},\overline{u}^{i},\overline{t}_{f}^{i})$, $i=1,2$. Set
$(\varphi_{1}^{*},\overline{\varphi}_{2}):=(\varphi_{1}(\overline{x}^{1}(\overline{t}_{f}^{1}),\overline{t}_{f}^{1}),\varphi_{2}(\overline{x}^{1}(\overline{t}_{f}^{1}),\overline{t}_{f}^{1}))$,
$(\overline{\varphi}_{1},\varphi_{2}^{*}):=(\varphi_{1}(\overline{x}^{2}(\overline{t}_{f}^{2}),\overline{t}_{f}^{2}),\varphi_{2}(\overline{x}^{2}(\overline{t}_{f}^{2}),\overline{t}_{f}^{2}))$
.

Step $0.2$

(Utopia point) Set $\beta^{*}:=(\beta_{1}^{*},\beta_{2}^{*})$ with
$\beta^{*}_{i}:=\varphi_{i}^{*}-\eta_{i}$, $i=1,2$.

Step $0.3$

(Essential interval) Determine the subinterval $[w_{0},w_{f}]\subset[0,1]$
using (3).

Step $0.4$

(Signs at end points) Compute $F^{\prime}(w_{0})$ and $F^{\prime}(w_{f})$
using (8), with $F(\cdot)$ evaluated as in (4).

  * •

If Fact 1(b)(i) or (ii) is satisfied then $w^{*}=w_{0}$ or $w^{*}=w_{f}$
appropriately; STOP.

  * •

If $F^{\prime}(w_{0})F^{\prime}(w_{f})=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]$.

Step $k.2$

(Stopping criterion) Compute $F^{\prime}(c)$ using (8), with $F(\cdot)$
evaluated as in (4).

  * •

If $F^{\prime}(c)=0$ or $(b-a)/2<\epsilon$ then set $w^{*}=c$ and STOP.

  * •

If $k=k_{\max}$ then declare “Maximum number of iterations exceeded.” and
STOP.

Step $k.3$

(New subinterval) Set $k:=k+1$ . If $F^{\prime}(a)F^{\prime}(c)>0$ then update
the subinterval as $[a,b]:=[c,b]$; otherwise, set $[a,b]:=[a,c]$. GO TO Step
$k.1$.

##  4 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 (OCPw), 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 (OCPw) 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.

###  4.1 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.

$\displaystyle\dot{x}_{1}(t)=x_{2}(t)\,,$
$\displaystyle\dot{x}_{2}(t)=-x_{1}(t)+x_{2}(t)\,(1.4-0.14\,x_{2}^{2}(t))+4\,u(t)\,,\mbox{\
\ for a.e. }t\in[0,t_{f}]\,,$

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\in[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

$(x_{1}(0),x_{2}(0))=(-5,-5)\quad\mbox{and}\quad(x_{1}(t_{f}),x_{2}(t_{f}))=(0,0)\,,$

and the dynamics are subject to constraints on the control variable such that

$-1\leq u(t)\leq 1\,,\mbox{\ \ for a.e. }t\in[0,t_{f}]\,.$

The optimal control problem is posed as a bi-objective problem with

$\min\ \left[t_{f}\,,\ \
\int_{0}^{t_{f}}\left(x_{1}^{2}(t)+u^{2}(t)\right)dt\right]\,,$

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

$\dot{x}_{3}(t)=x_{1}^{2}(t)+u^{2}(t)\,,\mbox{\ \ for a.e. }t\in[0,t_{f}]\,,\
\ x_{3}(0)=0\,.$

Then the two objective functionals as in Problem (OCP), or Problem (OCPw), can
be expressed as

$\varphi_{1}(x(t_{f}),t_{f})=t_{f}\quad\mbox{and}\quad\varphi_{2}(x(t_{f}),t_{f})=x_{3}(t_{f})\,.$

As we have stated above, the bi-objective Rayleigh problem is in the same form
as Problem (OCP) and, in particular, Problem (OCPw). The decision maker’s
objective for this problem will be to minimize a weighted distance to the
origin of the value space. We choose

$\varphi_{0}(x^{w},u^{w},t_{f}^{w}):=100\,\varphi_{1}^{2}(x^{w}(t_{f}),t_{f}^{w})+\varphi_{2}^{2}(x^{w}(t_{f}),t_{f}^{w})\,,$

where the scaling multiplier 100 is used to make the orders of magnitudes of
$\varphi_{1}$ and $\varphi_{2}$ the same. We aim to solve Problem (OPF), to
determine a scalar $w\in(0,1)$ with $w_{1}:=w$ and $w_{2}:=1-w$ that results
in the best Pareto solution in the sense that $\varphi_{0}(\cdot,\cdot,\cdot)$
is minimized, subject to the solution of Problem (OCPw).

In [35], Maurer and Oberle numerically illustrate that an optimal solution
does not exist for the single objective problem minimizing the quadratic
functional $\varphi_{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}\leq 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 (OCPw) can now explicitly be written for the Rayleigh problem as

$\left\\{\begin{array}[]{rll}\displaystyle\min_{{\alpha\geq
0}\atop{x(\cdot),u(\cdot),t_{f}}}&\ \ \alpha&\\\\[11.38109pt] \mbox{subject
to}&\ \ \dot{x}_{1}(t)=x_{2}(t)\,,&x_{1}(0)=-5\,,\
x_{1}(t_{f})=0\,,\\\\[5.69054pt] &\ \
\dot{x}_{2}(t)=-x_{1}(t)+x_{2}(t)\,(1.4-0.14\,x_{2}^{2}(t))+4\,u(t)\,,&x_{2}(0)=-5\,,\
x_{2}(t_{f})=0\,,\\\\[5.69054pt] &\ \
\dot{x}_{3}(t)=x_{1}^{2}(t)+u^{2}(t)\,,&x_{3}(0)=0\,,\\\\[5.69054pt] &\ \
-1\leq u(t)\leq 1\,,\mbox{\ \ for a.e. }t\in[0,t_{f}]\,,\ \ t_{f}\leq
5\,,&\\\\[5.69054pt] &\ \
w\,(t_{f}-\beta_{1}^{*})\leq\alpha\,,&\\\\[5.69054pt] &\ \
(1-w)\,(x_{3}(t_{f})-\beta_{2}^{*})\leq\alpha\,.&\end{array}\right.$

The Hamiltonian $H:{\rm{I\ \kern-5.39993ptR}}^{3}\times{\rm{I\
\kern-5.39993ptR}}\times{\rm{I\ \kern-5.39993ptR}}^{3}\to{\rm{I\
\kern-5.39993ptR}}$ for this problem simply is

$H(x,u,\lambda):=\lambda_{1}x_{2}+\lambda_{2}\left[(-x_{1}+x_{2}\,(1.4-0.14\,x_{2}^{2})+4u\right]+\lambda_{3}(x_{1}^{2}+u^{2})\,,$

where $\lambda(t):=(\lambda_{1}(t),\lambda_{2}(t),\lambda_{3}(t))\in{\rm{I\
\kern-5.39993ptR}}^{3}$ is referred to as the adjoint variable vector. Using
the convenient notation $H[t]:=H(x(t),u(t),\lambda(t))$, suppose that

$\displaystyle\dot{\lambda}_{1}(t):=-H_{x_{1}}[t]=\lambda_{2}(t)-2\lambda_{3}(t)x_{1}(t)\,,$
(9a)
$\displaystyle\dot{\lambda}_{2}(t):=-H_{x_{2}}[t]=-\lambda_{1}(t)-\lambda_{2}(t)(1.4-0.42x_{2}^{2}(t))\,,$
(9b) $\displaystyle\dot{\lambda}_{3}(t):=-H_{x_{3}}[t]=0\,,$ (9c)

for all $t\in[0,t_{f}]$, with certain transversality conditions as required by
the maximum principle. In (9a)–(9c), $H_{x_{i}}:=\partial H/\partial x_{i}$,
$i=1,2,3$. We will not go into the details of these (boundary) conditions
here. However we note that $\lambda_{3}(t)=\overline{\lambda}_{3}$, a
constant, for all $t\in[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 $\lambda(\cdot)$ satisfying (9a)–(9c), along with certain
transversality conditions, such that $\lambda(t)\neq 0$, for all
$t\in[0,t_{f}]$, and

$u(t)=\operatorname*{argmin}_{v\in[-1,1]}H(x(t),v,\lambda(t))=\operatorname*{argmin}_{v\in[-1,1]}\left(4\lambda_{2}(t)v+\lambda_{3}(t)v^{2}\right).$
(10)

for a.e. $t\in[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

$u(t)=\operatorname*{argmin}_{v\in[-1,1]}\lambda_{2}(t)v\,,$

resulting in

$u^{w}(t)=\left\\{\begin{array}[]{rl}1\,,&\ \ \mbox{if\ \
}\lambda^{w}_{2}(t)<0\,,\\\\[2.84526pt] -1\,,&\ \ \mbox{if\ \
}\lambda^{w}_{2}(t)>0\,,\\\\[2.84526pt] \mbox{undetermined}\,,&\ \ \mbox{if\ \
}\lambda^{w}_{2}(t)=0\,,\end{array}\right.$ (11)

for a.e. $t\in[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\in[w_{f},1]$. Recall
that if one does not have $\lambda^{w}_{2}(t)=0$ for all
$[t^{\prime},t^{\prime\prime}]\subset[0,t_{f}]$, where
$t^{\prime}<t^{\prime\prime}$, 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\in[w_{0},w_{f})$,

$u^{w}(t)=\left\\{\begin{array}[]{rl}1\,,&\ \ \mbox{if\ \
}2\lambda^{w}_{2}(t)<-\overline{\lambda}^{w}_{3}\,,\\\\[2.84526pt]
-2\lambda^{w}_{2}(t)/\overline{\lambda}^{w}_{3}\,,&\ \ \mbox{if\ \
}-\overline{\lambda}^{w}_{3}\leq
2\lambda^{w}_{2}(t)\leq\overline{\lambda}^{w}_{3}\,,\\\\[2.84526pt] -1\,,&\ \
\mbox{if\ \ }2\lambda^{w}_{2}(t)>\overline{\lambda}^{w}_{3}\,,\\\\[2.84526pt]
\end{array}\right.$ (12)

for all $t\in[0,t_{f}]$, provided $\overline{\lambda}^{w}_{3}\neq 0$. Again by
virtue of the discussion in Section 3.1, $u^{w}(t)$ in (12) is the same for
all $w\in[0,w_{0}]$. We define the switching function as

$\sigma^{w}(t):=\left\\{\begin{array}[]{rl}2\,\lambda^{w}_{2}(t)/\overline{\lambda}^{w}_{3}\,,&\
\ \mbox{if\ \ }0\leq w<w_{f}\,,\\\\[2.84526pt] 16\,\lambda^{w}_{2}(t)\,,&\ \
\mbox{if\ \ }w_{f}\leq w\leq 1\,.\end{array}\right.$ (13)

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.

$u^{w}(t)=\left\\{\begin{array}[]{rl}\left\\{\begin{array}[]{rl}1\,,&\ \
\mbox{if\ \ }\sigma^{w}(t)<-1\\\\[2.84526pt]
-2\lambda^{w}_{2}(t)/\overline{\lambda}^{w}_{3}\,,&\ \ \mbox{if\ \
}-1\leq\sigma^{w}(t)\leq 1\,,\\\\[2.84526pt] -1\,,&\ \ \mbox{if\ \
}\sigma^{w}(t)>1\,.\\\\[2.84526pt] \end{array}\right\\},&\ \ \mbox{if\ \
}0\leq w<w_{f}\,,\\\\[19.91692pt] \left\\{\begin{array}[]{rl}1\,,&\ \
\mbox{if\ \ }\sigma^{w}(t)<0\,,\\\\[2.84526pt] -1\,,&\ \ \mbox{if\ \
}\sigma^{w}(t)>0\,.\\\\[2.84526pt] \end{array}\right\\},&\ \ \mbox{if\ \
}w_{f}\leq w\leq 1\,.\end{array}\right.$ (14)

As to why $\sigma^{w}(\cdot)$ is referred to as the switching function should
now be more clear from (14): the value of $\sigma^{w}(\cdot)$ determines when
to switch from one case of the control function $u^{w}(\cdot)$ to another.

For Problem (OCPw) written for the Rayleigh problem above, we have chosen the
utopia vector as $(\beta_{1}^{*},\beta_{2}^{*})=(0,0)$, since
$\varphi_{i}(x(t_{f}),t_{f})>0$, for $i=1,2$. Figure 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).

The essential interval is found to be $[w_{0},w_{f}]=[0.8994,0.9269]$, with

$(\varphi_{1}^{w_{0}},\varphi_{2}^{w_{0}})=(5.000,44.71)\quad\mbox{and}\quad(\varphi_{1}^{w_{f}},\varphi_{2}^{w_{f}})=(3.668,46.50)\,,$

correct to four significant figures, where
$\varphi_{i}^{w}:=\varphi_{i}(x^{w}(t_{f}),t_{f}^{w})$, $i=1,2$, with
$w=w_{0}$ or $w_{f}$, or as will be the case below, $w=w^{*}$. Optimization
over the Pareto front results in $w^{*}=0.9247$, after 14 iterations of
Algorithm 1, yielding

$\varphi_{0}^{w^{*}}=58.71\quad\mbox{and}\quad(\varphi_{1}^{w^{*}},\varphi_{2}^{w^{*}})=(3.709,45.51)\,.$

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
$(\varphi_{1}^{w^{*}},\varphi_{2}^{w^{*}})=(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.

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

(b) Switching function as defined in (13).

(c) Phase plane trajectories.

(d) Control variable.

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.

The numerical Pareto-optimal state and control variable solutions are
presented in Figures 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\in[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\in[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 $\sigma^{w}(\cdot)$ plotted in Figure 2(b) by using
(13) (recall that discrete approximations of $\lambda^{w}_{2}(t)$ and
$\lambda^{w}_{3}(t)$ can readily be obtained from AMPL) furnishes the means to
verify the optimality condition for $u^{w}(\cdot)$ expressed in (14). It is
evident from the dashed (green) plot of the switching function that, for
$w\in[w_{f},1]$, when $\sigma^{w}(\cdot)$ crosses the time axis there is a
jump (from $1$ to $-1$ or vice versa) in the value of the corresponding
$u^{w}(\cdot)$ plot. Likewise, for $w\in[0,w_{0}]$ and for
$w=w^{*}\in[w_{0},w_{f})$, whenever $\sigma^{w}(\cdot)$ crosses one of the
lines $\sigma^{w}(t)=1$ and $\sigma^{w}(t)=-1$ (shown by two black lines in
Figure 2(b) for convenience) the expression for the control function
$u^{w}(\cdot)$ switches from one case in (14) to another, as required.

###  4.2 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):

$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.

The dynamical system is given by

$\left\\{\begin{aligned} \dot{S}(t)&=\mu N-\frac{\beta}{N}I(t)S(t)-\mu
S(t),\\\\[2.84526pt] \dot{L}_{1}(t)&=\frac{\beta}{N}I(t)\left(S(t)+\sigma
L_{2}(t)+\sigma_{R}R(t)\right)-\left(\delta+\tau_{1}+\epsilon_{1}u_{1}(t-d_{u_{1}})+\mu\right)L_{1}(t),\\\\[2.84526pt]
\dot{I}(t)&=\phi\,\delta\,L_{1}(t)+\omega
L_{2}(t)+\omega_{R}R(t)-\tau_{0}I(t-d_{I})\ +\mu I(t),\\\\[2.84526pt]
\dot{L}_{2}(t)&=(1-\phi)\delta
L_{1}(t)-\sigma\frac{\beta}{N}I(t)L_{2}(t)-(\omega+\epsilon_{2}u_{2}(t-d_{u_{2}})+\tau_{2}+\mu)L_{2}(t).\end{aligned}\right.$
(15)

The recovered population is defined by

$R(t):=N-S(t)-L_{1}(t)-I(t)-L_{2}(t)\,,$ (16)

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:

$\begin{array}[]{l}S(0)=76\,N/120,\ L_{1}(0)=36\,N/120,\ L_{2}(0)=2\,N/120,\
R(0)=N/120,\\\\[5.69054pt] I(t)=5\,N/120\quad\mbox{for}\ -\\!d_{I}\leq t\leq
0,\quad u_{k}(t)=0\quad\mbox{for}\ -\\!d_{u_{k}}\leq
t<0,\quad(k=1,2).\end{array}$ (17)

Symbol | Description | Value
---|---|---
$\beta$ | Transmission coefficient | variable
$\mu$ | Death and birth rate | $1/70\,yr^{-1}$
$\delta$ | Rate at which individuals leave $L_{1}$ | $12\,yr^{-1}$
$\phi$ | Proportion of individuals going to $I$ | $0.05$
$\omega$ | Endogenous reactivation rate for persistent latent infections | $0.0002\,yr^{-1}$
$\omega_{R}$ | Endogenous reactivation rate for treated individuals | $0.00002\,yr^{-1}$
$\sigma$ | Factor reducing the risk of infection as a result of acquired |
| immunity to a previous infection for $L_{2}$ | $0.25$
$\sigma_{R}$ | Rate of exogenous reinfection of treated patients | 0.25
$\tau_{0}$ | Rate of recovery under treatment of active TB | $2\,yr^{-1}$
$\tau_{1}$ | Rate of recovery under treatment of early latent individuals $L_{1}$ | $2\,yr^{-1}$
$\tau_{2}$ | Rate of recovery under treatment of persistent latent individuals $L_{2}$ | $1\,yr^{-1}$
$N$ | Total population | $30,000$
$\epsilon_{1}$ | Efficacy of treatment of early latent $L_{1}$ | $0.5$
$\epsilon_{2}$ | Efficacy of treatment of persistent latent TB $L_{2}$ | $0.5$
$t_{f}$ | Total simulation duration |  $5$ $years$
$d_{I}$ | delay in the diagnosis of $I$ |  $0.1$ $years$
$d_{u_{1}}$ | delay in the diagnosis of early latent individuals $L_{1}$ |  $0.2$ $years$
$d_{u_{2}}$ | delay in the prophylactic treatment of persistent latent individuals $L_{2}$ |  $0.2$ $years$
Table 1: Parameter values for the TB control model.

The control constraints are given by

$0\leq u_{k}(t)\leq 1\,,\quad\forall t\in[0,t_{f}]\,,\quad(k=1,2).$ (18)

We consider the following parametric objective functional with control weights
$a_{1},a_{2}\geq 0$:

$\int\limits_{0}^{t_{f}}(I(t)+L_{2}(t)+a_{1}u_{1}(t)+a_{2}u_{2}(t))\,dt\,.$
(19)

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,

$\begin{array}[]{l}x_{5}(t_{f}):=\displaystyle\int\limits_{0}^{t_{f}}\Big{(}I(t)+L_{2}(t)+a_{11}\,u_{1}(t)+a_{12}\,u_{2}(t)\Big{)}\,dt\,,\\\\[8.53581pt]
x_{6}(t_{f}):=\displaystyle\int\limits_{0}^{t_{f}}\Big{(}I(t)+L_{2}(t)+a_{21}\,u_{1}(t)+a_{22}\,u_{2}(t)\Big{)}\,dt\,.\end{array}$
(20)

with control weights $a_{11},a_{12},a_{21},a_{22}\geq 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

$\begin{array}[]{l}\dot{x}_{5}(t)=I(t)+L_{2}(t)+a_{11}\,u_{1}(t)+a_{12}\,u_{2}(t)\,,\
\ x_{5}(0)=0\,,\\\\[5.69054pt]
\dot{x}_{6}(t)=I(t)+L_{2}(t)+a_{21}\,u_{1}(t)+a_{22}\,u_{2}(t)\,,\ \
x_{6}(0)=0\,.\end{array}$ (21)

Denoting the (augmented) state vector by
$x(t)=(S(t),L_{1}(t),I(t),L_{2}(t),x_{5}(t),x_{6}(t))\in\mathbb{R}^{6}$ and
the control vector $u(t):=(u_{1}(t),u_{2}(t))\in\mathbb{R}^{2}$, the two
competing objectives in the general problem (P) are given by

$\varphi_{1}(x(t_{f}),t_{f})=x_{5}(t_{f})=:F_{1}(x,u)\quad\mbox{and}\quad\varphi_{2}(x(t_{f}),t_{f})=x_{6}(t_{f})=:F_{2}(x,u)\,,$

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

$\varphi_{0}(x^{w},u^{w},t_{f}^{w}):=\varphi_{1}^{2}(x^{w}(t_{f}),t_{f}^{w})+\varphi_{2}^{2}(x^{w}(t_{f}),t_{f}^{w})\,,$

Our aim is to solve Problem (OPF), to determine a scalar $w\in(0,1)$ with
$w_{1}:=w$ and $w_{2}:=1-w$ that results in the best Pareto solution in the
sense that $\varphi_{0}(\cdot,\cdot,\cdot)$ is minimized, subject to the
solution of Problem (OCPw).

Next we focus on the solution of Problem (OCPw): We aim to find a pair of
functions $(x,u)\in W^{1,\infty}([0,t_{f}],\mathbb{R}^{6})\times
L^{\infty}([0,t_{f}],\mathbb{R}^{2})$ that minimizes the parameter $\alpha$
subject to the time-delayed dynamics (15) and the auxiliary dynamics (21),
initial conditions (17), control constraints (18) and auxiliary weighted
inequalities involving $\varphi_{1}$ and $\varphi_{2}$.

We consider the necessary optimality conditions for the time-delayed optimal
control problem (OCPw); see Göllmann and Maurer [24], Vinter [45]. For this
purpose we introduce the delayed state variable
$y_{3}(t)=x_{3}(t-d_{I})=I(t-d_{I})$ and delayed control variables
$v_{k}(t)=u_{k}(t-d_{u_{k}})$, $k=1,2$. Denoting the adjoint variable vector
by
$\lambda(t):=(\lambda_{S}(t),\lambda_{L_{1}}(t),\lambda_{I}(t),\lambda_{L_{2}}(t),\lambda_{5}(t),\lambda_{6}(t))\in\mathbb{R}^{6}$
the Hamiltonian or Pontryagin function is given by

$\begin{array}[]{rl}H(x,y_{3},\lambda,u_{1},v_{1},u_{2},v_{2})=&\lambda_{s}\,(\mu
N-\frac{\beta}{N}IS-\mu S)\\\\[5.69054pt] &\ \ \,+\
\lambda_{L_{1}}\,(\frac{\beta}{N}I\left(S+\sigma
L_{2}+\sigma_{R}R\right)-\left(\delta+\tau_{1}+\epsilon_{1}v_{1}+\mu\right)L_{1})\\\\[5.69054pt]
&\ \ \,+\ \lambda_{I}\,(\,\phi\,\delta L_{1}+\omega
L_{2}+\omega_{R}\,R-\tau_{0}y_{3}+\mu I)\\\\[5.69054pt] &\ \ \,+\
\lambda_{L_{2}}\,((1-\phi)\delta
L_{1}-\sigma\frac{\beta}{N}IL_{2}-(\omega+\epsilon_{2}v_{2}+\tau_{2}+\mu)L_{2})\\\\[5.69054pt]
&\ \ \,+\ \lambda_{5}\,(I+L_{2}+a_{11}u_{1}+a_{12}u_{2})\\\\[5.69054pt] &\ \
\,+\ \lambda_{6}\,(I+L_{2}+a_{21}u_{1}+a_{22}u_{2})\,,\end{array}$ (22)

where $R$ is given as in (16). The Minimum Principle [24, 45] yields the
adjoint equations

$\dot{\lambda}_{S}(t)=-\frac{\partial H}{\partial
S}[t],\quad\dot{\lambda}_{L_{1}}(t)=-\frac{\partial H}{\partial
L_{1}}[t],\quad\dot{\lambda}_{L_{2}}(t)=-\frac{\partial H}{\partial
L_{2}}[t],$ $\dot{\lambda}_{x_{5}}(t)=-\frac{\partial H}{\partial
x_{5}}[t]=0\,,\quad\dot{\lambda}_{x_{6}}(t)=-\frac{\partial H}{\partial
x_{6}}[t]=0\,,$

and the advanced adjoint equation

$\dot{\lambda}_{I}(t)=-\frac{\partial H}{\partial
I}[t]-\chi_{[0,t_{f}-d_{I}]}(t)\frac{\partial H}{\partial I}[t+d_{I}]\,,$

where the argument $[t]$ stands for evaluating all arguments at time $t$. We
note that $\lambda_{5}^{w}(t)=\overline{\lambda}_{5}^{w}$ and
$\lambda_{6}^{w}(t)=\overline{\lambda}_{5}^{w}$, constants, for any fixed
$w\in[0,1]$. In the last equation, the term $\chi_{[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 $\sigma_{k}(t)$ for $k=1,2$:

$\begin{array}[]{rcl}\sigma_{k}^{w}(t)&=&\displaystyle\frac{\partial
H}{\partial u_{k}}[t]+\chi_{[0,t_{f}-d_{u_{k}}]}(t)\frac{\partial H}{\partial
v_{k}}[t+d_{u_{k}}]\\\\[11.38109pt]
&=&\left\\{\begin{array}[]{lcl}a_{1k}\overline{\lambda}_{5}^{w}+a_{2k}\overline{\lambda}_{6}^{w}-\
\epsilon_{k}\lambda_{L_{k}}^{w}(t+d_{u_{k}})L_{k}^{w}(t+d_{u_{k}})\,,&\mbox{if}&0\leq
t\leq t_{f}-d_{u_{k}}\,,\\\\[2.84526pt]
a_{1k}\overline{\lambda}_{5}^{w}+a_{2k}\overline{\lambda}_{6}^{w}\,,&\mbox{if}&t_{f}-d_{u_{k}}\leq
t\leq t_{f}\,.\end{array}\right.\end{array}$ (23)

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)

$u_{k}^{w}(t)=\left\\{\begin{array}[]{rcl}0\,,&\mbox{
if}&\sigma_{k}^{w}(t)>0\,,\\\ 1\,,&\mbox{
if}&\sigma_{k}^{w}(t)<0\,,\end{array}\right.\qquad k=1,2.$ (24)

for all $w\in[0,1]$. In particular, for positive weights $a_{1}>0$, $a_{2}>0$,
the switching functions (23) and the control law (24) imply

$u_{k}^{w}(t)=0\quad\forall\;t_{f}-d_{u_{k}}\leq t\leq t_{f}\,,$

for all $w\in[0,1]$.

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 $\varphi_{1}$
and $\varphi_{2}$.

For Problem (OCPw) written for the TB problem, we have chosen the utopia
vector as $(\beta_{1}^{*},\beta_{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 essential interval in this case is found to be
$[w_{0},w_{f}]=[0.5251,0.5709]$, with

$(\varphi_{1}^{w_{0}},\varphi_{2}^{w_{0}})=(28155,31133)\quad\mbox{and}\quad(\varphi_{1}^{w_{f}},\varphi_{2}^{w_{f}})=(26459,35205)\,,$

where $\varphi_{i}^{w}:=\varphi_{i}(x^{w}(t_{f}),t_{f}^{w})$, $i=1,2$, with
$w=w_{0}$ or $w_{f}$, or as will be the case below, $w=w^{*}$. Optimization
over the Pareto front results in $w^{*}=0.5358$, after 10 iterations of
Algorithm 1, yielding

$\varphi_{0}^{w^{*}}=41621\quad\mbox{and}\quad(\varphi_{1}^{w^{*}},\varphi_{2}^{w^{*}})=(27255,31455)\,.$

In Figure 3 only five iterations are displayed (labelled 1–5) for clarity in
viewing. The Pareto (master) solution with $w=w^{*}$ is represented by a
square.

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..

(a) Control variable $u_{1}^{w}$ (24) and scaled
switching function $\sigma_{1}^{w}$ (23) superposed.

(b) Control variable $u_{2}^{w}$ (24) and scaled
switching function $\sigma_{2}^{w}$ (23) superposed.

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.

The numerical Pareto-optimal control variable solutions $u_{1}^{w}(\cdot)$ and
$u_{2}^{w}(\cdot)$ are presented in Figures 4(a)–(b) for
$w=w_{0},w^{*},w_{f}$. As with Rayleigh, one of the boundary Pareto-optimal
solutions is shown using solid (blue) curves for $w=w_{0}$, the same solution
for all $w\in[0,w_{0}]$. The other boundary Pareto-optimal solution for
$w=w_{f}$, which holds for all $w\in[w_{f},1]$, is shown using dashed (green)
curves. Both of the control solutions are of bang–bang 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, $\sigma_{k}^{w}(\cdot)$,
$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 $\lambda^{w}_{L_{k}}(t)$,
$k=1,2$, $\lambda^{w}_{5}(t)$ and $\lambda^{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_{s_{k}}^{w}$, $k=1,2$, at which the respective $u_{k}^{w}(\cdot)$ 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_{5}^{w_{f}}(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_{6}^{w_{0}}(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.

Scalarization | Functional values | Switching times | Terminal state values
---|---|---|---
weight $w$ | $x_{5}^{w}(t_{f})$ | $x_{6}^{w}(t_{f})$ | $t_{s_{1}}^{w}$ | $t_{s_{2}}^{w}$ | $S^{w}(t_{f})$ | $L_{1}^{w}(t_{f})$ | $I^{w}(t_{f})$ | $L_{2}^{w}(t_{f})$ | $R^{w}(t_{f})$
$w_{0}=0.5251$ : | 28155 | 31133 | 0.145 | 2.864 | 1193.1 | 28.2 | 13.3 | 864.0 | 27901.4
$w^{*}=0.5358$ : | 27255 | 31455 | 0.809 | 3.439 | 1205.8 | 27.5 | 13.0 | 747.6 | 28006.1
$w_{f}=0.5709$ : | 26459 | 35205 | 4.083 | 4.752 | 1238.2 | 23.8 | 11.2 | 419.3 | 28307.5
Table 2: TB problem.

##  5 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 (OCPw) 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 (OCPw) 20 to 30 times. On the other hand, without the algorithm
we propose, it is necessary to construct the Pareto front by solving (OCPw)
thousands of times in order to obtain the same solution with the same
computational accuracy.

The proposed algorithm can be improved/modified in various ways. For example,
scalarization techniques other than Chebyshev might be employed; see for
example [8, 9] and the references therein. Bisection method might be replaced
by methods with higher convergence rates, for example regula falsi and secant
methods (see [10]), at the expense of approximating higher order derivatives
of course, although the latter would make the algorithm applicable to problems
with more than just two objective functionals.

## References

  * [1] Alt, W., Baier, R. , Lempio, F., Gerdts, M.: Approximations of linear control problems with bang–bang solutions. Optimization, 62, 9–32 (2013)
  * [2] Benson, H.P.: Optimization over the efficient set, J. Math. Anal. Appl.. 98, 562–580 (1984)
  * [3] Benson, H.P.: A finite, non-adjacent extreme point search algorithm for optimization over the efficient set. J. Optim. Theory Appl., 73, 47–64 (1992)
  * [4] Betts, J.T.: Practical Methods for Optimal Control Using Nonlinear Programming, Third Edition. Advances in Design and Control, SIAM Publications, Philadelphia (2020)
  * [5] Bolintinéanu, S.: Optimality conditions for minimization over the (weakly or properly) efficient set. J. Math. Anal. Appl., 173(2), 523–541 (1993)
  * [6] Bolintinéanu, S.: Minimization of a quasi-concave function over an efficient set. Math. Prog., 61, 89–110 (1993)
  * [7] Bonnel, H., Kaya, C.Y.: Optimization over the efficient set of multi-objective convex optimal control problems. J. Optim. Theory Appl., 147, 93–11 (2010)
  * [8] Burachik, R.S., Kaya, C.Y., Rizvi, M.M.: A new scalarization technique to approximate Pareto fronts of problems with disconnected feasible sets. J. Optim. Theory Appl., 162, 428–446 (2014)
  * [9] Burachik, R.S., Kaya, C.Y., Rizvi, M.M.: A new scalarization technique and new algorithms to generate Pareto fronts. SIAM J. Optim., 27, 1010–1034 (2017)
  * [10] Burden, R.L., Faires, J.D.: Numerical Analysis, 9th edition. Thompson Brooks/Cole, Belmont, CA, USA (2011)
  * [11] Büskens, C.: Optimierungsmethoden und Sensitivitätsanalyse für optimale Steuerprozesse mit Steuer– und Zustands–Beschränkungen. PhD Thesis, Institut für Numerische Mathematik, Universität Münster, Germany (1998)
  * [12] Büskens, C., Maurer, H.: SQP-methods for solving optimal control problems with control and state constraints: Adjoint variables, sensitivity analysis and real-time control. J. Comp. Appl. Math., 120, 85–108 (2000)
  * [13] Chorobura, A.P.: Multi-objective infinite horizon optimal control problems: characterization of the Pareto fronts and Pareto solutions. Comp. Appl. Math., 40, 258 (2021) URL: https://doi.org/10.1007/s40314-021-01633-0.
  * [14] Dauer J.P.: Optimization over the efficient set using an active constraint approach. Z. Oper. Res., 35, 185–195 (1991).
  * [15] Dauer, J.P., Fosnaugh, T.A.: Optimization over the efficient set, J. Global Optim., 7, 261–277 (1995).
  * [16] Désilles, A., Zidani, H.: Pareto front characterization for multiobjective optimal control problems using Hamilton–Jacobi approach. SIAM J. Control Optim., 57, 3884–3910 (2019)
  * [17] Dontchev, A.L., and Hager, W.W.: Lipschitz stability in in nonlinear control and optimization, SIAM J. Control Optim., 31, 569–603 (1993).
  * [18] Dontchev, A.L., Hager, W.W.: The Euler approximation in state constrained optimal control problems. Math. Comput., 70, 173–203 (2001)
  * [19] Dontchev, A.L., Hager, W.W., Malanowski, K.: Error bound for Euler approximation of a state and control constrained optimal control problem. Numer. Funct. Anal. Optim., 21(6), 653–682 (2000)
  * [20] Dontchev, A.L., Hager, W.W., Veliov, V.M.: Second-order Runge-Kutta approximations in control constrained optimal control, SIAM J. Num. Anal., 38(1), 202–-226 (2000)
  * [21] Dutta, J., Kaya, C.Y.: A new scalarization and numerical method for constructing the weak Pareto front of multi-objective optimization problems. Optimization, 60, 1091–1104 (2011)
  * [22] Eichfelder, G.: Adaptive Scalarization Methods in Multiobjective Optimization. Springer, Berlin, Heidelberg (2008)
  * [23] Fourer, R., Gay, D.M., Kernighan, B.W.: AMPL: A Modeling Language for Mathematical Programming, Second Edition. Brooks/Cole Publishing Company / Cengage Learning (2003)
  * [24] Göllmann, L., Maurer, H.: Theory and applications of optimal control problems with multiple time-delays, J. of Industrial and Management Optimization, 10, 413–441 (2014)
  * [25] Horst, R., Thoai, N.V.: Maximizing a concave function over the efficient or weakly-efficient set. European J. Oper. Res., 117, 239–252 (1999)
  * [26] Horst, R., Thoai, N.V., Yamamoto, Y., Zenke, D.: On optimization over the efficient set in linear multicriteria programming. J. Optim. Theory Appl., 134, 433–443 (2007)
  * [27] Kaya, C.Y., Maurer, H.: A numerical method for nonconvex multi-objective optimal control problems. Comp. Optim. Appl., 57, 685–702 (2014)
  * [28] Logist, F., van Erdeghem, P.M.M., van Impe, J.F.: Efficient deterministic multiple objective optimal control of (bio)chemical processes. Chem. Eng. Sci., 64, 2527–2538 (2009)
  * [29] Logist, F., Houska, B., Diehl, M., van Impe, J.: Fast Pareto set generation for nonlinear optimal control problems with multiple objectives. Struct. Multidisc. Optim., 42, 591-603 (2010)
  * [30] Logist, F., Vallerio, M., Houska, B., Diehl, M., van Impe, J.: Multi-objective optimal control of chemical processes using ACADO toolkit. Comp. Chem. Eng., 37, 191–199 (2012)
  * [31] Liu, Z., Ehrgott, M.: Primal and dual algorithms for optimization over the efficient set. Optimization, 67,1661–1686 (2018)
  * [32] Malanowski, K., Maurer, H.: Sensitivity analysis for parametric control problems with control-state constraints. Computational Optimization and Applications, 5, 253-283 (1996)
  * [33] Malanowski, K., Maurer, H.: Sensitivity analysis for state constrained optimal control problems. Discrete and Continuous Dynamical Systems, 4, 241-272 (1998)
  * [34] Maurer, H., Büskens, C., Kim, J.-H.R., Kaya, C.Y.: Optimization methods for the verification of second-order sufficient conditions for bang–bang controls. Optim. Contr. Appl. Meth., 26, 129–156 (2005)
  * [35] Maurer, H., Oberle, H.J.: Second order sufficient conditions for optimal control problems with free final time: the Riccati Approach. SIAM J. Control Optim., 41(2), 380–403 (2002)
  * [36] Maurer, H., Pesch, H.J.: Solution differentiability for nonlinear parametric control problems. SIAM J. Control Optim., 32, 1542–1554 (1994)
  * [37] Maurer, H., Pesch, H.J.: Solution differentiability for parametric nonlinear control problems with control-state constraints. J. Optim. Theory Appl., 86, 285–309 (1995)
  * [38] Miettinen, K.M.: Nonlinear Multiobjective Optimization, Kluwer (1999)
  * [39] Ober-Blöbaum, S., Ringkamp, M., zum Felde, G.: Solving multiobjective optimal control problems in space mission design using discrete mechanics and reference point techniques. Proceedings of the 51st IEEE Conference on Decision and Control, Dec. 10-13, Maui, Hawaii, USA, pp. 5711–5716 (2012)
  * [40] Osmolovskii, N.P., H. Maurer, H.: Applications to Regular and Bang-Bang Control: Second-Order Necessary and Sufficient Optimality Conditions in Calculus of Variations and Optimal Control. SIAM Advances in Design and Control, Vol. DC 24, SIAM Publications, Philadelphia, 2012.
  * [41] Philip, J.: Algorithms for the vector maximization problem. Math. Prog., 207–229 (1972).
  * [42] Pietrus, A., Scarinci, T., Veliov, V.: High order discrete approximations to Mayer’s problems for linear systems SIAM J. Control Optim., 56(1), 102–119 (2018)
  * [43] Silva, C. J., Maurer, H., Torres, D.F.M.: Optimal control of a Tuberculosis model with state and control delays. Math. Biosci. Eng., 14, 321–337 (2017)
  * [44] Tuberculosis. World Health Organization, 14 October 2021 and 27 October 2022. URL: https://www.who.int/news-room/fact-sheets/detail/tuberculosis. Accessed: 15 August 2022 and 14 January 2023.
  * [45] Vinter, R. B.: State constrained optimal control problems with time delays. J. Math. Anal. Appl., 457(2), 1696–1712 (2018)
  * [46] Wächter, A., Biegler, L.T.: On the Implementation of a primal-dual interior point filter line search algorithm for large-scale nonlinear programming. Math. Progr., 106, 25–57 (2006)
  * [47] Yamamoto, Y.: Optimization over the efficient set : overview. J. Global Optim., 22(1-4), 285–317 (2002).