A Survey of Recent Trends in Multiobjective Optimal Control—Surrogate Models, Feedback Control and Objective Reduction

Abstract

Multiobjective optimization plays an increasingly important role in modern applications, where several criteria are often of equal importance. The task in multiobjective optimization and multiobjective optimal control is therefore to compute the set of optimal compromises (the Pareto set) between the conflicting objectives. The advances in algorithms and the increasing interest in Pareto-optimal solutions have led to a wide range of new applications related to optimal and feedback control, which results in new challenges such as expensive models or real-time applicability. Since the Pareto set generally consists of an infinite number of solutions, the computational effort can quickly become challenging, which is particularly problematic when the objectives are costly to evaluate or when a solution has to be presented very quickly. This article gives an overview of recent developments in accelerating multiobjective optimal control for complex problems where either PDE constraints are present or where a feedback behavior has to be achieved. In the first case, surrogate models yield significant speed-ups. Besides classical meta-modeling techniques for multiobjective optimization, a promising alternative for control problems is to introduce a surrogate model for the system dynamics. In the case of real-time requirements, various promising model predictive control approaches have been proposed, using either fast online solvers or offline-online decomposition. We also briefly comment on dimension reduction in many-objective optimization problems as another technique for reducing the numerical effort.

1. Introduction

There is hardly ever a situation where only one goal is of interest at the same time. When performing a purchase for example, we want to pay a low price while getting a high quality product. In the same manner, multiple goals are present in most technical applications, maximizing quality versus minimizing the cost being only one of many examples. This dilemma leads to the field of multiobjective optimization, where we want to optimize all relevant objectives simultaneously. However, this is obviously impossible as the above example illustrates. Generically, the different objectives contradict each other such that we are forced to choose a compromise. While we are usually satisfied with one optimal solution in the scalar-valued setting, there exists in general an infinite number of optimal compromises in the situation where multiple objectives are present. The set of these compromise solutions is called the Pareto set, the corresponding points in the objective space form the Pareto front.

Since the solution to a Multiobjective Optimization Problem (MOP) is a set, it is significantly more expensive to compute than the optimum of a single objective problem, and many researchers devote their work to the development of algorithms for the efficient numerical approximation of Pareto sets. These advances have opened up new challenging application areas for multiobjective optimization. In optimal control, the optimization variable is not finite-dimensional, but rather a function, typically depending on time. The goal is to steer a dynamical system in such a way that one (or multiple) objective is minimized. Two particularly challenging control problems are feedback control and control problems constrained by Partial Differential Equations (PDEs). In the first case, the time for computing the Pareto set is strictly limited, often to a small fraction of a second. In the latter case, even the solution of single objective problems is often extremely time consuming so that the development of new algorithmic ideas is necessary to make these problems computationally feasible.

In fact, in situations like these, surrogate models or dimension reduction techniques are a promising approach for significantly reducing the computational effort and thereby enabling real-time applicability. This article gives an overview of recent advances in surrogate modeling for multiobjective optimal control problems, where the approach is to replace the underlying system dynamics by a reduced order model, which can be solved much faster. On the other hand, it introduces an approximation error, which has to be taken into account when analyzing convergence properties. The article is structured as follows. In Section 2, we are going to review some basics about multiobjective optimization, including the most popular solution methods. The two main challenges are addressed individually in the next sections, starting with expensive models in Section 3. There are also surrogate modeling approaches for MOPs that directly provide a mapping from the control variable to the objective function values. Since there already exist extensive surveys for this case (cf. [1,2,3], for instance), these are summarized only very briefly. In Section 4, real-time feedback control is discussed. Finally, we briefly discuss the question of dimension reduction in the number of objectives in Section 5 before concluding with a summary of further research directions in Section 6.

2. Multiobjective Optimization

In this section, the concepts of multiobjective optimization and Pareto optimality will be introduced and some widely-used solution methods will be summarized. More detailed introductions to multiobjective optimization can be found in, e.g., [4,5].

2.1. Theory

In multiobjective optimization, we want to minimize multiple objectives at the same time. Consequently, the fundamental difference from scalar optimization is that the objective function $J:\mathcal{U}\to {\mathbb{R}}^k$ is vector-valued. Hence, the general problem is of the form:

$$\begin{array}{cc}\hfill \underset{u\in \mathcal{U}}{min}J\left(u\right)=& \underset{u\in \mathcal{U}}{min}\left(\begin{array}{c}{J}_1\left(u\right)\\ ⋮\\ J_k\left(u\right)\end{array}\right)\hfill \\ \hfill \mathrm{s}.\mathrm{t}.g_i\left(u\right)& \le 0,i=1,\dots ,l,\hfill \\ \hfill h_j\left(u\right)& =0,j=1,\dots ,m,\hfill \end{array}$$

(MOP)

where $u\in \mathcal{U}$ is the control variable and $g:\mathcal{U}\to {\mathbb{R}}^l$, $g\left(u\right)={(g_1\left(u\right),\dots ,g_l\left(u\right))}^{\top }$ and $h:\mathcal{U}\to {\mathbb{R}}^m$, $h\left(u\right)={(h_1\left(u\right),\dots ,h_m\left(u\right))}^{\top }$, are inequality and equality constraints, respectively. The space of the control variables $\mathcal{U}$ is also called the decision space (according to the term decision variable for u in classical multiobjective optimization), and the objective function maps u to the objective space. Depending on the problem setup, $\mathcal{U}$ can either be finite-dimensional, i.e., $\mathcal{U}={\mathbb{R}}^n$, or some appropriate function space.

Remark 1.

It is common in finite-dimensional optimization to use the notation x for the control or optimization variable and F for the objective function. In contrast to that, u and J are more common for control problems. In order to unify the notation, the latter will be used throughout this article for all optimization and optimal control problems.

In contrast to classical optimization, in optimal control, we have to compute an input in such a way that a dynamical system behaves optimally with respect to some specified cost functional. Hence, we have the system dynamics as an additional constraint, very often in the form of ordinary (ODEs) or partial differential equations (PDEs):

$$\begin{array}{cccc}\hfill \dot{y}(x,t)& =G\left(y\right(x,t),u(t\left)\right),\hfill & \hfill (x,t)& \in \mathsf{\Omega }\times (t_0,t_e],\hfill \\ \hfill a(x,t)\frac{\partial y}{\partial n}(x,t)+b(x,t)y(x,t)& =c(x,t),\hfill & \hfill (x,t)& \in \mathsf{\Gamma }\times (t_0,t_e],\hfill \\ \hfill y(x,t_0)& =y_0\left(x\right),\hfill & \hfill x& \in \mathsf{\Omega },\hfill \end{array}$$

(PDE)

where the domain of interest $\mathsf{\Omega }\subset {\mathbb{R}}^{n_x}$ is a connected open set with spatial dimension $n_x$ and the boundary is denoted by $\mathsf{\Gamma }=\partial \mathsf{\Omega }$ with outward normal vector n. The coefficients $a(x,t)$, $b(x,t)$ and $c(x,t)$ in the boundary condition are given by the problem definition. The operator G is a partial differential operator describing the evolution of the system. The cost functional $\widehat{J}:\mathcal{U}\times \mathcal{Y}\to {\mathbb{R}}^k$ of an optimal control problem consequently depends on the control u, as well as the system state y, which results in a multiobjective optimal control problem:

$$\begin{array}{cc}\hfill \underset{u\in \mathcal{U},y\in \mathcal{Y}}{min}\widehat{J}(u,y)=& \underset{u\in \mathcal{U},y\in \mathcal{Y}}{min}\left(\begin{array}{c}{\int }_{t_0}^{t_e}{C}_1(y(x,t),u\left(t\right))dt+{\mathsf{\Phi }}_1\left(y(x,t_e)\right)\\ ⋮\\ {\int }_{t_0}^{t_e}{C}_k(y(x,t),u\left(t\right))dt+{\mathsf{\Phi }}_k\left(y(x,t_e)\right)\end{array}\right)\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& \left(\mathrm{PDE}\right)\hfill \\ & g_i(y,u)\le 0,i=1,\dots ,l,\hfill \\ & h_j(y,u)=0,j=1,\dots ,m.\hfill \end{array}$$

(MOCP)

There are articles on multiobjective optimal control that specifically address the implications of multiple objectives for optimal control [6,7]; see also [8] for a short survey of methods. However, for many problems, there exists a unique solution y for every u such that (MOCP) can be simplified by introducing a so-called control-to-state operator $\mathcal{S}:\mathcal{U}\to \mathcal{Y}$; see [9] for details. By setting $\widehat{J}(u,y)=\widehat{J}(u,\mathcal{S}u)=J\left(u\right)$, the problem is transformed into (MOP). For this reason, we will from now on only consider (MOP).

In the situation where $\mathcal{U}$ is a function space (i.e., in the case of optimal control), the problem can be numerically transformed into a high-, yet finite-dimensional problem in a direct solution method via discretization, cf. [10,11]. This results in a large number of control variables, which can be very challenging on its own in multiobjective optimization. If the system dynamics are governed by a PDE, then the spatial discretization of the state y results in an even higher number of unknowns, which can easily reach several millions or more [12].

In contrast to single objective optimization problems, there exists no total order of the objective function values in ${\mathbb{R}}^k$ with $k\ge 2$ (unless they are not conflicting). Therefore, the comparison of values is defined in the following way [4]:

Definition 1.

Let $v,w\in {\mathbb{R}}^k$. The vector v is less than w (denoted by $v<w$), if $v_i<w_i$ for all $i\in \left\{1,\dots ,k\right\}$. The relation ≤ is defined in an analogous way.

A consequence of the lack of a total order is that we cannot expect to find isolated optimal points. Instead, the solution to (MOP) is the set of optimal compromises (also called the Pareto set or set of non-dominated points):

Definition 2.

Consider the multiobjective optimization problem (MOP). Then:

1. 

a point $u^{\ast }$ dominates a point u, if $J\left(u^{\ast }\right)\le J\left(u\right)$ and $J\left(u^{\ast }\right)\ne J\left(u\right)$.

2. 

a feasible point $u^{\ast }$ is called globally Pareto optimal if there exists no feasible point $u\in \mathcal{U}$ dominating $u^{\ast }$. The image $J\left(u^{\ast }\right)$ of a globally Pareto optimal point $u^{\ast }$ is called a globally Pareto optimal value. If this property holds in a neighborhood $U\left(u^{\ast }\right)\subset \mathcal{U}$, then $u^{\ast }$ is called locally Pareto optimal.

3. 

the set of non-dominated feasible points is called the Pareto set ${\mathcal{P}}_S$ and its image the Pareto front ${\mathcal{P}}_F$.

A consequence of Definition 2 is that for each point that is contained in the Pareto set (the red line in Figure 1a), one can only improve one objective by accepting a trade-off in at least one other objective. Figuratively speaking, in a two-dimensional problem, we are interested in finding the “lower left” boundary of the feasible set in objective space (cf. Figure 1b).

Similar to single objective optimization, a necessary condition for optimality is based on the gradients of the objective functions. The first order conditions were independently discovered by Karush in 1939 [13] and by Kuhn and Tucker in 1951 [14]. Due to this, they are widely known as the Karush–Kuhn–Tucker (KKT) conditions:

Theorem 1

([14]). Let $u^{\ast }$ be a Pareto-optimal point of Problem (MOP), and assume that $\nabla h_j\left(u^{\ast }\right)$ for $j=1,\dots ,m$ and $\nabla g_s\left(u^{\ast }\right)$ for $s=1,\dots ,l$ are linearly independent. Then, there exist non-negative scalars ${\alpha }_1,\dots ,{\alpha }_k\ge 0$ with ${\sum }_{i=1}^k{\alpha }_i=1$, $\gamma \in {\mathbb{R}}^m$ and $\mu \in {\mathbb{R}}^l$ such that:

$$\begin{array}{cc}\hfill \sum _{i=1}^k{\alpha }_i\nabla J_i\left(u^{\ast }\right)+\sum _{j=1}^m{\gamma }_j\nabla h_j\left(u^{\ast }\right)+\sum _{s=1}^l{\mu }_s\nabla g_s\left(u^{\ast }\right)& =0,\hfill \\ \hfill h_j\left(u^{\ast }\right)& =0,j=1,\dots ,m,\hfill \\ \hfill g_s\left(u^{\ast }\right)& \le 0,s=1,\dots ,l,\hfill \\ \hfill {\mu }_s{g}_s\left(u^{\ast }\right)& =0,s=1,\dots ,l,\hfill \\ \hfill {\mu }_s& \ge 0,s=1,\dots ,l.\hfill \end{array}$$

(KKT)

The set of points satisfying these conditions is called the set of substationary points ${\mathcal{P}}_{S,\mathsf{sub}}$ [4]. Obviously, ${\mathcal{P}}_{S,\mathsf{sub}}$ is a superset of the Pareto set ${\mathcal{P}}_S$. Many algorithms for MOPs compute the set of substationary points, in particular gradient-based methods, as we will see in the next section. This set can be reduced to the Pareto set in a consecutive step by performing a (comparatively inexpensive) non-dominance test.

2.2. Solution Methods

Many researchers in multiobjective optimization focus their attention on developing efficient algorithms for the computation of Pareto sets. Algorithms for solving MOPs can be compiled into several fundamentally different categories of approaches. The first category is based on scalarization techniques, where ideas from single objective optimization theory are extended to the multiobjective situation. All scalarization techniques have in common that the Pareto set is approximated by a finite set of Pareto-optimal points, which are computed by solving scalar subproblems. Consequently, the resulting solution methods involve solving multiple optimization problems consecutively. Scalarization can be achieved by various approaches such as the weighted-sum method, the $ϵ$-constraint method, normal boundary intersection or reference point methods [4,5,15].

Continuation methods make use of the fact that under certain conditions, the Pareto set is a smooth manifold of dimension $k-1$ [16]. This means that one can compute the tangent space in each point of the set, and a predictor step is performed in this space. The resulting point then has to be corrected to a Pareto-optimal solution using a descent method [17].

Another prominent approach is based on evolutionary algorithms [18,19], where the underlying idea is to evolve an entire population of solutions during the optimization process. Significant advances have been made concerning Multiobjective Evolutionary Algorithms (MOEAs) in recent years [20,21] (see also [22] for a survey) such that they are nowadays the most popular choice for solving MOPs due to the applicability to very complex problems and being easy to use in a black box fashion. Since convergence rates can be relatively slow for MOEAs, they can be coupled with locally fast methods close to the Pareto set. These approaches are known as memetic algorithms; see, e.g., [23,24,25,26].

Set-oriented methods provide an alternative deterministic approach to the solution of MOPs. Utilizing subdivision techniques, the desired Pareto set is approximated by a nested sequence of increasingly refined box coverings [27,28,29]. This way, a superset is computed, which converges to the desired solution, even in situations where the Pareto set is disconnected. However, their complexity depends on both the dimension of the Pareto set, as well as the decision space dimension. Due to this, one has to take additional steps to apply these algorithms for multiobjective optimal control problems.

Depending on the method of choice, gradient information can be used to accelerate convergence. While this is widely accepted in scalar-valued optimization, this is less the case when multiple objectives are present [30]. Nonetheless, many approaches exist where gradients are exploited, for example in order to create sequences converging to single points [31,32,33], to compute the entire set of valid decent directions [30,34], to obtain superlinear or quadratic convergence [35,36] or in combination with evolutionary approaches (memetic algorithms) [37,38,39,40,41]. In many of the gradient-based methods, the descent direction for all objectives is a convex combination of the individual gradients:

$$q\left(u\right)=\sum _{i=1}^k{\widehat{\alpha }}_i\nabla J_i\left(u\right).$$

(1)

Here, $\widehat{\alpha }$ is a fixed weight vector, which is determined in such a way that:

$$\langle q\left(u\right),\nabla J_i\left(u\right)\rangle <0\forall i=1,\dots ,k,$$

(2)

see, e.g., [31,32]. In the unconstrained case, there exists no u satisfying (2) only if $\parallel q\left(u\right)\parallel =0$, which implies that u is substationary, cf. (KKT).

3. Surrogate Models

The ever-increasing computational capabilities allow us to analyze more and more complicated systems with a very large number of degrees of freedom, also in the context of optimal control, where practical problems range from process control [42,43] and energy management [15] over space mission design [44] to mobility and autonomous driving [45,46,47,48].

The above-mentioned examples can all be described by ordinary differential equations with a finite-dimensional state space $\mathcal{Y}$. In contrast to these problems, many phenomena in physics such as mechanical strain, heat flow, electromagnetism, fluid flow or even multi-physics simulations are governed by partial differential equations. Using numerical discretization schemes for the approximation of the spatial domain (such as finite elements or finite volumes) results in a very large number of degrees of freedom and a heavy computational burden. For more complex systems (such as turbulent flows [49]), simulating the dynamics is already very costly. Consequently, optimal control of these systems is all the more challenging, and considering multiple objectives further increases the cost. Due to this reason, only a few problems have been addressed directly; see, e.g., [50,51,52,53]. A method exploiting the special structure in the system dynamics has been proposed in [54], and a special case of Pareto-optimal solutions, namely Nash equilibria, has been computed in [55,56] (When not using a priori selection methods such as scalarization or Nash equilibria, a decision maker selects the appropriate solution. This is called Multi-Criteria Decision Making (MCDM) [57] and is an entire area of research on its own. Thus, we will not go into further details about the decision making process here.).

A very popular approach to circumvent the problem of prohibitively large computational cost is the use of surrogate models. Here, the exact objective function $J\left(u\right)$ is replaced by a surrogate $J^r\left(u\right)$, where the superscript r stands for reduced. In many situations, this surrogate function can be evaluated faster by several orders of magnitude. The challenge here is to find a good trade-off between acceleration and model accuracy, and many approaches have been proposed over the past two decades.

In optimal control, there are two fundamentally different ways for model reduction. The first case, which is equivalently applicable to multiobjective optimization problems, is to directly derive a surrogate for the objective function, i.e., $J^r:\mathcal{U}\to {\mathbb{R}}^k$ is constructed by polynomials, radial basis functions or other means. In optimal control, an alternative way of reducing the computational effort is by introducing a reduced model for the system dynamics:

$$J^r\left(u\right)={\widehat{J}}^r(u,y)=\widehat{J}(u,{\mathcal{S}}^{r}u),$$

where the reduced control-to-state operator ${\mathcal{S}}^r$ indicates that the model reduction is due to a surrogate model for the system dynamics.

In both situations, we cannot expect that $J^r\left(u\right)=J\left(u\right)$ holds for all $u\in \mathcal{U}$. Instead, we introduce an error, which has to be taken into account. This is closely related to questions concerning uncertainty and noise. In this context, many researchers have addressed inaccuracies. In [58], the notion of $ϵ$-efficiency (cf. Definition 3) was first introduced in order to handle uncertain objective values. Uncertainty has also been considered in the context of multiobjective evolutionary computation (see, e.g., [59,60,61]. Alternative). Methods such as probabilistic [62,63], deterministic [64] or set-oriented approaches [65,66] have also been proposed. The special case of many-objective optimization is covered in [67], and applications are addressed in [68,69]. A different approach to uncertainties is via robust algorithms. Several examples from multiobjective optimization, as well as optimal control can be found in [43,70,71,72].

In the following, we will first introduce some results concerning inaccuracies in multiobjective optimization in Section 3.1 and then give an overview of existing methods for both of the above-mentioned approaches, i.e., surrogate models for the objective function (Section 3.2) or for the system dynamics (Section 3.3). Since the first approach has already been covered extensively in several surveys [1,2,3], we only give a brief overview of the existing methods and the corresponding references.

3.1. Inaccuracies and $ϵ$-Dominance

When using surrogate models in order to accelerate the solution process, we have to accept an error both in the objective function, as well as the respective gradients. Furthermore, inaccuracies may occur due to stochastic processes or due to unknown model parameters. In these situations, the objective function and the corresponding gradients are only known approximately, which has to be taken into account.

Suppose now that we only have approximations $J^r\left(u\right),\nabla J^r\left(u\right)$ of the objectives $J_i\left(u\right)$ and the gradients $\nabla J_i\left(u\right)$, $i=1,\dots ,k$, respectively. Furthermore, let us assume that upper bounds $ϵ,\kappa \in {\mathbb{R}}^k$ for these errors are known:

$$\begin{array}{cc}\hfill \parallel J_i^r\left(u\right)-J_i{\left(u\right)\parallel }_2& \le {ϵ}_i\forall u\in \mathcal{U},\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill \parallel \nabla J_i^r\left(u\right)-\nabla J_i{\left(u\right)\parallel }_2& \le {\kappa }_i\forall u\in \mathcal{U}.\hfill \end{array}$$

(4)

In this situation, we need to replace the dominance property 2 by an inexact version, also known as $ϵ$-dominance (see also [58,69]):

Definition 3

([66]). Consider the multiobjective optimization problem (MOP) where the objective function $J\left(u\right)$ is only known approximately according to (3). Then:

1. 

a point $u^{\ast }$ confidently dominates a point u, if $J^r\left(u^{\ast }\right)+ϵ\le J^r\left(u\right)-ϵ$ and $J_i^r\left(u^{\ast }\right)+{ϵ}_i<J_i^r\left(u\right)-{ϵ}_i$ for at least one $i\in 1,\dots ,k$.

2. 

The set of almost non-dominated points, which is a superset of the Pareto set ${\mathcal{P}}_S$, is defined as:

$$\begin{array}{c}\hfill {\mathcal{P}}_{S,ϵ}=\left\{u^{\ast }\in \mathcal{U}|\nexists u\in \mathcal{U}\mathit{which}\mathit{confidently}\mathit{dominates}{u}^{\ast }\right\}.\end{array}$$

(5)

The concept of $ϵ$-dominance is visualized in Figure 2. Theoretically, the true point could be anywhere inside the box defined by $ϵ$ such that in the cases (a)–(c), the lower left point does not confidently dominate the other point. The necessary condition is violated for one component in (a) and (b), respectively, and for both components in (c). The gray points in Figure 2a show a possible realization of the true points in which no point is dominated by the other. In (d), the orange point confidently dominates the black one, and in (e), we see the implications for the computation of Pareto fronts. Due to the inexactness, the number of points that are not confidently dominated is larger than in the exact case. This is also evident in Figure 3, where the exact and the inexact solution of an example problem from production [32] have been computed with an extension of the subdivision technique presented in [27], cf. [66] for details. Here, inexactness is introduced due to uncertainties in pricing. $ϵ$-dominance can be used for the development of algorithms for MOPs with uncertainties [59,61,62,63,65,67,68], for accelerating expensive MOPs [60,66,73], as well as for increasing the number of compromise solutions for the decision maker [63,64,65,69].

When considering gradient-based methods, inaccuracies in the gradients (Equation (4)) have to be taken into account. Since only approximations of the true gradient are known, these result in an inexact descent direction $q^r\left(u\right)$. The inaccuracy in the gradients introduces an upper bound for the angle between the individual gradients $\nabla J_i^r\left(u\right)$ and the descent direction $q^r\left(u\right)$. This is equivalent to a lower bound ${\widehat{\alpha }}_{\mathsf{min}}$ for the weight vector, i.e., ${\widehat{\alpha }}_i\in [0,1]$ in Equation (5) has to be replaced by ${\widehat{\alpha }}_i^r\in [{\widehat{\alpha }}_{\mathsf{min},i},1]$ for $i=1,\dots ,k$:

$$q^r\left(u\right)=\sum _{i=1}^k{\widehat{\alpha }}_i^r\nabla J_i^r\left(u\right);$$

(6)

see [66] for a detailed discussion. The additional constraint ensures that if $q^r\left(u\right)$ is a descent direction for the inexact problem, it is also a descent direction for the true problem. Furthermore, we obtain a criterion for the accuracy up to which we can compute ${\mathcal{P}}_{S,\mathsf{sub}}$ based on inexact gradient information.

Theorem 2

([66]). Consider the multiobjective optimization problem (MOP) without constraints. Only approximate gradients according to (4) are available, and consequently, the descent direction is also only known approximately according to Equation (6). Assume that $\parallel q^r{\left(u\right)\parallel }_2\ne 0$ and $\parallel \nabla J_i^r{\left(u\right)\parallel }_2\ne 0$ for $i=1,\dots ,k$. Let:

$$\begin{array}{c}\hfill {\widehat{\alpha }}_{\mathsf{min},i}=\frac{1}{\parallel \nabla J_i^r{\left(u\right)\parallel }_2^2}\left(\parallel q^r{\left(u\right)\parallel }_2{\kappa }_i-\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^k{\widehat{\alpha }}_j\left(\nabla J_i^r\left(u\right)·\nabla J_i^r\left(u\right)\right)\right),i=1,\dots ,k.\end{array}$$

Then, the following statements are valid:

(a) 

If ${\sum }_{i=1}^k{\widehat{\alpha }}_{\mathsf{min},i}>1$ then there exists no direction $q\left(u\right)$ with:

$$\langle q\left(u\right),\nabla J_i\left(u\right)\rangle <0\forall i=1,\dots ,k,$$

i.e., no descent direction for the exact problem.

(b) 

All points u with ${\sum }_{i=1}^k{\widehat{\alpha }}_{\mathsf{min},i}=1$ are contained in the set:

$$\begin{array}{c}\hfill {\mathcal{P}}_{S,\kappa }=\left\{u\in {\mathbb{R}}^n|||\sum _{i=1}^k{\widehat{\alpha }}_i\nabla J_i\left(u\right){||}_2\le 2{\parallel \kappa \parallel }_{\infty }\right\}.\end{array}$$

(7)

A combination of Theorem 2 with the subdivision algorithm from [27] is shown in Figure 4. The algorithm constructs a nested sequence of increasingly refined box coverings, which converges to the set of substationary points where in the unconstrained case, $\parallel q\left(u^{\ast }\right)\parallel =0$ holds for all $u^{\ast }\in {\mathcal{P}}_{S,\mathsf{sub}}$. The set ${\mathcal{P}}_{S,\mathsf{sub}}$ is shown in red in Figure 4a. Due to the inexactness, we can no longer guarantee $\parallel q\left(u^{\ast }\right)\parallel =0$. Instead, we obtain the set ${\mathcal{P}}_{S,\kappa }$, which is shown in green. The background is colored according to the optimality condition $\parallel q\left(u\right)\parallel$ of the exact problem, and the dashed white line indicates the error bound (7) from the above theorem.

3.2. Surrogate Models for the Objective Function

The most straight-forward approach for introducing a surrogate model is to directly construct a map from the control to the decision space. This means that only the essential input-output behavior is covered, whereas internal states, as well as the system dynamics are neglected. Using such a meta model, one can very quickly obtain the objective function value for every u. This approach is equally applicable in finite-dimensional multiobjective optimization and has been used extensively in this context. To this end, we will only briefly cover the main questions that have to be addressed when using such models.

In special cases, one can exploit some structure in the problem formulation that yields a simplified analytic expression. However, in most cases, even if the equations can be written down in closed form, this approach requires a deep understanding of the underlying system, which is often hard to obtain. Moreover, small changes in the problem setup may result in having to repeat this tedious process all over. Due to these reasons, data-based approaches are used much more frequently. They are often easy to apply and much more general. In these approaches, the original objective function (or even a real-world experiment) is evaluated for a small number of inputs that are contained in the set ${\mathcal{U}}_{\mathsf{ref}}=\{u^1,\dots ,u^p\}$, where p is the number of function evaluations (or experiments). The data points $J\left(u^1\right)$–$J\left(u^p\right)$ are then used to fit the meta model such as, e.g., the coefficients of a polynomial basis. Obviously, the choice of suitable ansatz functions is essential for the success of a meta modeling strategy. Popular choices are:

• Response Surface Models (RSM);
• Radial Basis Functions (RBF);
• statistical models such as Kriging or Gauss regression;
• machine learning methods such as artificial neural networks or support vector machines;

See [1] for an extensive survey in the context of multiobjective optimization. Additional surveys can be found in [74], where different statistical methods are compared, in [3,75] in relation to MOEAs or in [76], where RSM and RBF are compared for crashworthiness problems.

Besides selecting the correct meta model, questions concerning the training dataset have to be answered:

• How large does the set ${\mathcal{U}}_{\mathsf{ref}}$ have to be?
• How can we pick the correct elements for ${\mathcal{U}}_{\mathsf{ref}}$?
• Do we define ${\mathcal{U}}_{\mathsf{ref}}$ in advance or online during the model building process?

In addition to the meta model, Point 1 significantly depends on the problem under consideration. The more non-linear a problem is, the more data points are generally required to accurately construct a meta model. Obviously, the number also depends on Point 2, the locations $u^1$–$u^k$ of these evaluations. The question of choosing the correct location is closely related to the field of optimum experimental design or Design of Experiments (DoE); see, e.g., [77,78]. The relevant question there is how to optimally pick a set of experiments such that the overall approximation error becomes minimal. Such approaches have successfully been coupled with multiobjective optimization in [79,80]. Point 3 depends both on the meta modeling approach, as well as the problem under consideration. In some situations (e.g., for real-world experiments), it may not be possible to iteratively determine the experiments. Instead, a batch approach has to be used. In computer experiments, flexibility is often higher such that an interplay between model building and high-fidelity evaluations can help to further reduce the number of experiments. In the context of machine learning, this process is also known as active learning [81,82].

Many algorithms combining multiobjective optimization and meta modeling have been proposed, and there exists a vast literature concerning this topic.

Table 1. Overview of publications (in chronological order) where surrogate modeling and multiobjective optimization are combined. MOEA, Multiobjective Evolutionary Algorithm; RSM, Response Surface Model; POD, Proper Orthogonal Decomposition; TR, Trust Region.

Surveys
Tabatabei et al. [2], Chugh et al. [3]	Extensive surveys on meta modeling for MOEAs
Voutchkov and Keane [74], Knowles and Nakayama [1], Jin [75]	Surveys on meta modeling approaches from statistics (RSM, RBF) and machine learning in combination with MOEAs
Benner et al. [83], Taira et al. [85], Peherstorfer et al. [84]	Surveys on reduced order modeling of dynamical systems
Algorithms Using Meta Models for the Objective Function
Ong et al. [112], Ray et al. [113]	Combination of RBF and MOEA
Chung and Alonso [114], Keane [115]	Combination of kriging models and MOEA
Karakasis and Giannakoglou [116]	RBF as an inexpensive pre-processing step in a MOEA
Knowles [117]	Combination of DoE and an interactive method
Zhang et al. [118]	Combination of Gaussian process models and scalarization
Telen et al. [79]	Combination of DoE and scalarization and MOEA
Chugh et al. [119]	Kriging model in combination with reference vector approach for MOPs
Meta Models Specifically Tailored to Multiobjective Optimization
Shimoyama et al. [120]	Kriging surrogate for hypervolume approximation (MOEA)
Pan et al. [121]	Surrogate model for dominance relations with uncertainties
Algorithms Using Surrogate Models for the System Dynamics
Iapichino et al. [105]	Combination of POD and weighted sum
Banholzer et al. [108,109]	Combination of POD and reference point method
Iapichino et al. [106]	Combination of RB and weighted sum
Peitz [73]	Combination of TR-POD and reference point method
Beermann et al. [110,111]	Combination of POD and set-oriented method
Applications
Albunni et al. [53]	POD and MOEA applied to the Maxwell equation
Ma and Qu [80]	MO of a switched reluctance motor by coupling RSM and MOEA (particle swarm optimization)
Peitz et al. [107]	POD-based multiobjective optimal control of the Navier–Stokes equations via scalarization and set-oriented methods
Wang et al. [122]	MOEA with multi-fidelity surrogate-management and offline-online decomposition applied to a trauma system

on page 17 mentions both survey articles and a (non-exhaustive) selection of popular methods.

3.3. Surrogate Models for the Dynamical System

The above-mentioned meta modeling methods are widely studied and have been successfully applied in a large variety of multiobjective optimization problems. In the context of control, there exists an alternative option, which is to derive a surrogate model for the system dynamics, i.e., to replace the high-fidelity control-to-state operator $\mathcal{S}$ by a less expensive surrogate model ${\mathcal{S}}^r$. Since the largest part of the computational effort is due to solving the dynamical system, by this, the reduced cost function $J^r\left(u\right)={\widehat{J}}^r(u,y)=\widehat{J}(u,{\mathcal{S}}^{r}u)$ is also much less expensive to evaluate.

The use of Reduced Order Models (ROMs) is not limited to control, but has been used successfully in a large variety of multi-query problems [12], and several extensive surveys have been written on reduced order modeling and using ROMs in prediction, uncertainty quantification or optimization; see [83,84,85]. For nonlinear problems, the two most widely-used approaches are the reduced basis method and Proper Orthogonal Decomposition (POD) (also known as Principal Component Analysis (PCA) or Karhunen–Loève decomposition).

3.3.1. ROMS via Proper Orthogonal Decomposition or the Reduced Basis Method

Numerically solving a PDE is generally realized by discretizing the spatial domain with a numerical mesh using finite differences, finite elements or finite volumes. By this, the infinite-dimensional state space $\mathcal{Y}$ is transformed into a finite-dimensional space ${\mathcal{Y}}^N$ via Galerkin projection:

$$y(x,t)\approx y^N(x,t)=\sum _{i=1}^N{z}_i^N\left(t\right){\varphi }_i\left(x\right).$$

(8)

Here, N denotes the number of degrees of freedom, and ${\varphi }_i\left(x\right)$ are basis functions with local support such as indicator functions or hat functions. This transforms the PDE into an N-dimensional ordinary differential equation for the coefficients z. For complex domains, as well as complex dynamics, the dimension can easily reach the order of millions such that solving the problem in ${\mathcal{Y}}^N$ can quickly become very expensive, which is particularly challenging in the multi-query case.

The general concept in projection-based model reduction is therefore to find an appropriate space ${\mathcal{Y}}^r$ with dimension $\ell \ll N$ in which the system dynamics can nevertheless be approximated with sufficient accuracy. The two most common approaches to do this are the Reduced Basis (RB) method and Proper Orthogonal Decomposition (POD). In both cases, we compute s so-called snapshots of the high-dimensional system and then use the dataset $\{y_1^N,\dots ,y_s^N\}$ to construct a reduced basis $\psi =\{{\psi }_1,\dots ,{\psi }_{\ell }\}$ such that:

$$y(x,t)\approx y^N(x,t)\approx y^r(x,t)=\sum _{i=1}^{\ell }{z}_i^{\ell }\left(t\right){\psi }_i\left(x\right).$$

The following example describes this approach in more detail. For an extensive introduction, the reader is referred to [86,87].

Example 1 (Heat equation).

Suppose we want to solve the time-dependent heat equation on a domain Ω with homogeneous Neumann conditions on the boundary Σ:

$$\begin{array}{cccc}\hfill y_t(x,t)-\lambda \Delta y(x,t)& =0,\hfill & \hfill (x,t)& \in \mathsf{\Omega }\times (t_0,t_e],\hfill \\ \hfill y(x,0)& =y_0\left(x\right),\hfill & \hfill x& \in \mathsf{\Omega },\hfill \\ \hfill y_n(x,t)& =0,\hfill & \hfill (x,t)& \in \mathsf{\Sigma }\times (t_0,t_e].\hfill \end{array}$$

(9)

Here, y is the temperature, and λ is the heat conductivity. The subscripts t and n indicate the derivatives with respect to time and the outward normal vector of the boundary, respectively. We now derive the weak form of (9) by multiplying with a test function φ and integrating over the domain Ω. Using Gauss’s theorem, we obtain the following equation:

$$\begin{array}{ccc}\hfill {\int }_{\mathsf{\Omega }}{y}_t(·,t)\phi -\lambda \nabla y(·,t)·\nabla \phi dx& =0,\hfill & \hfill (x,t)\in \mathsf{\Omega }\times (t_0,t_e].\end{array}$$

(10)

If we want to solve (10) using the finite element method, we have to insert the Galerkin ansatz (8) into (10) and individually take each of the basis functions as a test function. By this, we obtain the following system of equations:

$${\int }_{\mathsf{\Omega }}\sum _{i=1}^N{z}_{i,t}^N\left(t\right){\varphi }_i{\varphi }_j-\lambda z_i^N\left(t\right)\nabla {\varphi }_i·\nabla {\varphi }_{j}dx=0,j=1,\dots ,N.$$

(11)

Introducing the mass matrix $M\in {\mathbb{R}}^{N\times N}$ and the stiffness matrix $K\in {\mathbb{R}}^{N\times N}$ with:

$$M_{i,j}={\int }_{\mathsf{\Omega }}{\varphi }_i{\varphi }_{j}dx,K_{i,j}={\int }_{\mathsf{\Omega }}\nabla {\varphi }_i·\nabla {\varphi }_{j}dx,$$

this yields the following N-dimensional linear system:

$$M{z}_t^N\left(t\right)-\lambda K{z}^N\left(t\right)=0.$$

If we now want to compute a reduced order model instead of a high-dimensional finite element approximation, we can apply the same procedure, except that now, we have to use the reduced basis in the Galerkin ansatz, as well as for test functions.

The most important difference between RB and POD is the area of application, although this is not a strict separation. RB is mostly applied to parameter-dependent, yet time-independent (i.e., elliptic) problems, whereas POD (introduced in [88]) is applied to time-dependent problems described by parabolic or hyperbolic PDEs. Consequently, in RB, the snapshots $\{y_1^N,\dots ,y_s^N\}$ are solutions corresponding to parameters $\{u_1,\dots ,u_s\}$, and in POD, they are snapshots in time, collected at the time instants $\{t_0,\dots ,t_{s-1}\}$. Using an equidistant time grid h, the snapshots are taken at $\{t_0,t_0+h\dots ,t_0+(s-1)h\}$. In RB, the snapshots $\{y_1^N,\dots ,y_s^N\}$ often directly serve as the basis $\psi$. For time-dependent problems, this can cause numerical difficulties since some snapshots might be very similar (e.g., for very slow systems or periodic dynamics) such that the snapshot matrix $S=(y_1^N,\dots ,y_s^N)$ is ill-conditioned. Due to this, a singular value decomposition is performed on S, and the ℓ leading left singular vectors are taken as the basis $\psi$. This results in an orthonormal basis, which can be shown to be optimal with respect to the $L^2$ projection error [87,88]. Furthermore, the truncation error is given by the sum over the neglected singular values:

$${ϵ}_{\ell }=\frac{{\sum }_{i=\ell +1}^s{\sigma }_i}{{\sum }_{j=1}^s{\sigma }_j}.$$

Whereas the error between the infinite-dimensional solution and the solution via a standard discretization approach can be neglected in many situations, the error of the ROM depends on several factors such as the reference data, the basis size and the parameter or control for which the ROM is evaluated. Consequently, this error can be significant such that proper care has to be taken. The most common approach is to derive bounds either for the error of the reduced state $\parallel y^r-y^N\parallel$ or, in the case of optimal control, of the optimal solution $\parallel u^r-u^N\parallel$ obtained using the ROM; see, e.g., [89,90,91,92,93] for POD or [94,95,96,97] for RB methods. In addition, there are other measures that can be taken such as deriving balanced input-output behavior [98,99] or introducing additional terms [100] or modifications [101] in the POD-based ROM. For more detailed introductions to RB and POD, the reader is referred to [97] and [87,88], respectively.

3.3.2. Optimal Control Using Surrogate Models

There is a rich literature on optimal control of PDEs using surrogate models. The approaches can be summarized into three main categories:

• build a model once,
• construction of regular updates in a trust region approach,
• construction of regular updates using error estimators.

Whereas the first category is the most efficient one (see, e.g., [102] for optimal control of the Navier–Stokes equations), it is in general not possible to prove the convergence of the resulting algorithm.

In the second approach (which was developed by Fahl [103] for POD-based ROMs and one objective), one defines a trust region within which the current surrogate model is considered as trustworthy; see Figure 5 for an illustration. The ROM-based optimal control problem is then solved with the additional constraint that the solution has to remain within the trust region, i.e., $\parallel u^i-u^{\mathsf{ref}}\parallel \le {\delta }^i$, where i is the current step of the iterative optimization scheme and ${\delta }^i$ is the current trust region radius. After having obtained $u^i$, the high-dimensional system is evaluated, and the improvement of the full system is determined:

$$\rho =\frac{|J^N\left(u^i\right)-J^N\left(u^{i-1}\right)|}{|J^r\left(u^i\right)-J^r\left(u^{i-1}\right)|}.$$

If $\rho$ is close to one, then the ROM is sufficiently accurate, and the iterate $u^i$ is accepted. We then use the high-dimensional solution to construct the next ROM at $u^{\mathsf{ref}}=u^i$. If, on the other hand, $\rho$ is close to zero, then the ROM accuracy was bad, and the iterate $u^i$ is rejected. Instead, the trust region radius ${\delta }^i$ is reduced, and the optimal control problem is solved again. Using the Trust Region (TR-POD) approach, one can ensure convergence to the optimal solution of the high-dimensional problem. In the case of the Navier–Stokes equations, this has been shown for different problem setups [103,104].

In the third approach, error estimators ${\Delta }^J$ for the current iterate $u^i$ are required. By evaluating ${\Delta }^J\left(u^i\right)$, it is possible to efficiently estimate the error between the high- and the low-dimensional solution, since:

$$|J^r\left(u^i\right)-J^N\left(u^i\right)|<{\Delta }^J.$$

If this error estimate is larger than some prescribed upper bound $ϵ$, then the ROM has to be updated using data from a high-dimensional solution. Detailed information on error estimates can be found in, e.g., [91,92,93,94,95,96,97].

3.4. ROM-Based Multiobjective Optimal Control of PDEs

All three approaches for using ROMs in optimal control have recently been extended to multiobjective optimal control problems. Besides different ROM techniques, different algorithms for MOPs have been used, as well, such that a variety of methods has evolved, each of which is well suited for certain situations.

3.4.1. Scalarization

A natural and widely-used approach to MOPs is via scalarization. By this, the vector of objectives is synthesized into a scalar objective function, and the MOP is transformed into a sequence of scalar optimization problems for different scalarization parameters. In terms of ROM-based optimal control, this is advantageous because many techniques from scalar-valued optimal control can be extended. The main difference is now that the objective function may have a more complicated structure. In [105,106], the weighted sum method has been used in combination with RB in order to solve MOPs constrained by elliptic PDEs. In the weighted sum method, scalarization is achieved via convex combination of the individual objectives using the weight vector $\alpha$:

$$\underset{u\in \mathcal{U}}{min}\overline{J}\left(u\right)=\underset{u\in \mathcal{U}}{min}\sum _{i=1}^k{\alpha }_i{J}_i\left(u\right).$$

(12)

The weighted sum method is probably the most straight-forward approach for including ROMs in MOPs. However, the method has strong limitations in the situation of non-convex problems, where it is impossible to compute the entire Pareto set [5].

A more advanced approach is the so-called reference point method [5] (cf. Figure 6 for an illustration), where the distance to an infeasible target point T with $T<J\left(u\right)$ has to be minimized:

$$\underset{u\in \mathcal{U}}{min}\overline{J}\left(u\right)=\underset{u\in \mathcal{U}}{min}\parallel T-J_i\left(u\right)\parallel .$$

(13)

By adjusting the target, we can move along the Pareto front and hence obtain an approximately equidistant covering of the front. The reference point method has been coupled with all three of the above-mentioned ROM approaches. In [107], it was used for multiobjective optimal control of the Navier–Stokes equations using one reduced model (cf. Figure 7a). Here, the objectives are to stabilize a periodic solution (the well-known von Kármán vortex street) and to minimize the control cost at the same time. In [73], the trust region framework by Fahl (TR-POD, [103]) was extended (cf. Figure 7b–c for a heat flow problem with a tracking and a cost minimization objective). The third ROM approach was used in [108,109]. The difficulty here is that the minimization of the distance to the target point results in a more complicated objective function, which has to be treated carefully.

Scalarization techniques generally have the same limitations on the decision space dimension as scalar-valued optimal control problems. This means that very efficient techniques (both direct and indirect) exist for high-dimensional controls. However, the number of objectives is limited because the parametrization in the scalarization step becomes extremely tedious, and it is almost impossible to obtain a good approximation of the entire Pareto set for more than three objectives.

3.4.2. Set-Oriented Approaches with $ϵ$-Dominance

In contrast to scalarization, in set-oriented techniques, the Pareto set is approximated by a box covering [27,28,29]. Here, the limitations are contrary, meaning that the dimension of the decision space is rather limited, while the number of objectives does not pose any problems for the algorithms. In practice, however, the computational cost increases exponentially with the number of objectives (i.e., the dimension of the Pareto set) such that we are still limited to a moderate number of objectives.

First results coupling the subdivision algorithm developed in [27] with error estimates for POD-based ROMs have recently appeared [110,111]. In the subdivision algorithm, the decision space is divided into boxes, which are alternatingly subdivided and selected. In the subdivision step, each existing box is subdivided into two smaller boxes. In the selection step, all boxes are eliminated that are dominated, i.e., they do not cover any part of the Pareto set. Numerically, this is realized by representing a box by a finite number of sample points and then marking a box as dominated if all sample points are dominated by samples from another box in the covering; see [27] for details.

The subdivision algorithm can be extended to inaccuracies by replacing the strict dominance test by an $ϵ$-dominance test as presented in Section 3.1 (see also Figure 3). After fixing an upper bound $ϵ$, we have to ensure that the surrogate models we use do not violate this bound anywhere in the control domain. Since this cannot be achieved with a single ROM, one has to use multiple, locally valid ROMs instead (cf. Figure 8c). The covering by local ROMs is managed in such a way that all points in the neighborhood around the reference $u^{\mathsf{ref}}$ at which the data for the ROM were collected satisfy the prescribed error bound $ϵ$. This way, the number of solutions of the high-dimensional system can be reduced significantly. A comparison between the exact and the ROM-based solution is shown in Figure 8 for a semilinear heat flow MOCP with two tracking type objectives and a cost minimization objective. Due to the ROM approach, the number of evaluations of the FEM model could be reduced by a factor of ≈1000.

3.5. Summary

Before moving on to feedback control, a summary of the relevant publications where multiobjective optimization and meta modeling interact is given in Table 1. The references are categorized into surveys, algorithms using surrogate models for the objective function, the system dynamics and specific reduction approaches for MOPs. Furthermore, some applications are referenced.

4. Feedback Control

Even when the objective function is not very expensive to evaluate, MOCPs often have a large computational cost; see, e.g., [123] for various examples. This becomes a limiting factor in situations where the solution time is critical as is the case in real-time applications. Due to the increasing computational power, as well as the advances in algorithms, Model Predictive Control (MPC) (see [124,125] for extensive introductions) has become a very powerful and widely-used method for realizing model-based feedback control of complex systems.

In MPC, an optimal control problem is solved in a short time horizon (the prediction horizon) while the real system (the plant) is running. Then, the first entry of the optimal control is applied to the plant, and the process is repeated with the time frame moving forward by one sample time h; see Figure 9 for an illustration. This way, a closed-loop behavior is achieved. On the downside, we have to solve the optimal control problem within the sample time h. This can be in the order of seconds or minutes (in the case of chemical processes) down to a few microseconds, for example in power electronics applications.

In many MPC problems, stabilization of the system with respect to some reference state is the most important aspect. Nevertheless, there exist variations where stability is not an issue such that other (more economic) objectives can be pursued. These methods are known as economic MPC [126,127]. Another variation is the so-called explicit MPC [128], where the optimal control is computed in advance for a large number of different states and stored in a library. This way, the computational effort is shifted to an offline phase, and during operation, we only have to select the optimal control from the library.

As in multiobjective optimal control, there are numerous applications where multiple objectives are interesting in feedback control. This means that we would have to solve the problem (MOCP) with $t_0=t_s$ and $t_e=t_{s+p}$ repeatedly within the sample time h. It is immediately clear that even the simplest MOCPs cannot be solved fast enough to allow for real-time applicability. Consequently, efficient algorithms have to be developed, which can be divided into approaches where one Pareto optimal solution is computed online and approaches with an offline phase during which the MOCP is solved (in this article, algorithms where Artificial Neural Networks (ANN) have to be trained beforehand are nonetheless assigned to the first category if the optimization is performed entirely online).

When implementing a multiobjective MPC (MOMPC) algorithm, one should keep in mind that, regardless of the algorithm used, the resulting trajectory need not be Pareto optimal, even if each single step is, cf. [129]. A remedy to this issue is presented in [130], where the selection of compromise solutions is restricted to a part of the Pareto front that is determined in the first MPC step. Due to this, upper bounds for the objective function values can be guaranteed.

4.1. Online Multiobjective Optimization

In the classical MPC framework, the optimal control problem is solved online within the sample time h. Since it is in general impossible to approximate the entire Pareto set sufficiently accurately within this time frame, there are three alternatives:

• compute a single Pareto-optimal solution according to some predefined preference,
• compute only a rough approximation of the Pareto set,
• compute an arbitrary Pareto-optimal control that satisfies additional constraints (e.g., stability).

In the first approach, the objective function is scalarized using, for instance, the weighted sum method (12) or the reference point method (13). In this situation, well-established approaches from scalar-valued MPC exist on which one can build. First results using the weighted sum method have appeared in [131]. In [132,133], the authors use the same scalarization for an MPC problem with convex objective functions. In this situation, it is guaranteed that any Pareto-optimal solution can be computed using weighted sums, and the weights can be adapted online according to a decision maker’s preference. Due to the convexity, stability can be proven for the resulting MPC algorithm. In [134], this approach is extended by providing gradient information of the objectives with respect to the weight vector $\alpha$. This way, the weights can be adapted in such a way that a desired change in the objective space is realized. For non-convex problems, the weighted sum method is incapable of computing the entire Pareto set. Therefore, in [135], a variation of the reference point method is applied, where the target T is the utopian point $J^{\ast }$, i.e., the vector of the individual minima. This way, also non-convex problems can be treated. In fact, due to the reference point method, the objective function is always convex [109], which can be exploited during the optimization. Alternative scalarization methods are the $ϵ$-constraint method [136] or lexicographic ordering [137].

A disadvantage of a priori scalarization is that it is often difficult to select the scalarization parameter in such a way that a desired trade-off solution is obtained, and the remedy proposed in [134] is only applicable to a specific class of problems. Therefore, an alternative approach is to quickly compute a rough approximation of the entire Pareto set and then select the desired control online. Such methods have been proposed by many authors. The general approach is to us an MOEA and stop the computation after a few iterations. In the next step, one of these suboptimal solutions is selected. This selection is realized by specifying a weight vector for the objectives in [138,139,140,141] and by the satisficing trade-off method in [142].

As a third option, we can compute a single Pareto-optimal point without specifying which one we are specifically interested in as long as it satisfies additional constraints such as the stability of the system. Approaches of this type have been developed in [136,143], where a game theoretic approach is used.

4.2. Offline-Online Decomposition

A well-known trick to avoid heavy online computations is to introduce an offline-online decomposition (very similar to meta modeling approaches where surrogate models are constructed before solving the MOCP). This means that the Pareto set is computed beforehand, and in the online phase, an optimal compromise is selected according to a decision maker’s preference or some heuristic based on the system state or the environment.

Many of the approaches that fall into this category use a standard feedback controller instead of MPC; see [144] for a short review concerning methods using scalarization and offline PID controller optimization. In the offline phase, a Pareto set is computed for the controller parameters. Possible objectives are, among many others, overshooting behavior, energy efficiency or robustness. Algorithms of this type have been proposed in [145,146,147] using MOEAs and in [148,149,150] using set-oriented methods.

An alternative approach is motivated by explicit MPC, i.e., the idea of solving many MOPs offline such that the correct solution can be extracted from a library in the online phase. Such a method has been proposed in [48]. In the offline phase, one has to identify all possible scenarios that can occur in the online phase. Such a scenario consists of both system states, as well as constraints. This results in a large number of MOPs that have to be solved. In order to reduce this number, symmetries in the problem are exploited. To this end, a concept known as motion primitives [151,152] is extended. In short, this means that if:

$$arg\underset{u\in \mathcal{U}}{min}{\mathrm{MOP}}_1=arg\underset{u\in \mathcal{U}}{min}{\mathrm{MOP}}_2,$$

where MOP${}_1$ and MOP${}_2$ are two problem instances from the offline library, then we only have to solve one of the problems in order to have a Pareto-optimal solution for both. Moreover, if two problem instances vary only slightly, one can use a previously-computed solution as a good initial guess for the next MOCP to further decrease the computational effort [73]. In the online phase, the correct Pareto set is selected from the library (according to the system state and the constraints), and an optimal compromise is selected according so a decision maker’s preference $\alpha$. In contrast to the affine linear solutions, which can be computed in explicit MPC for linear-quadratic problems, one has to rely on interpolation between solutions in the nonlinear setting.

4.2.1. Example: Autonomous Driving

We here want to demonstrate the superiority of multiobjective approaches over scalar-valued MPC using the example of autonomously-driving electric vehicles [47,48,153]. The problem there is to find the set of optimal engine torque profiles such that the velocity is maximized while the energy consumption is minimized. Additional constraints have to be taken into account such as speed limits or stop signs. The system dynamics are described by a four-dimensional, highly nonlinear ODE for the vehicle velocity, the battery state of charge and two battery voltage drops, cf. [153] for details. Numerical investigations reveal several symmetries in the system such that the only relevant state for a scenario is the current velocity, whereas all other states only have a minor influence on the solution of the MOP. Consequently, the velocity, as well as the constraints form the above-mentioned scenarios; see Figure 10a for an illustration. For example, a scenario could be that the current velocity is 60 km/h and that the speed limit is currently increasing from 50–100 km/h, cf. Scenario (II) in Figure 10a. We then solve the MOP for this scenario and store the Pareto set in a library. By discretizing the velocity into steps of $0.1$ (i.e., $v\left(t_0\right)=\dots ,59.9,60.0,60.1,\dots$), we have to solve 1727 MOPs in total in the offline phase.

In the online phase, we now select the relevant Pareto set from the library and, according to a decision maker’s preference, apply one of the Pareto-optimal controls to the electric vehicle. This is done repeatedly such that a feedback loop is realized. The result is illustrated in Figure 10b for an example track, where the black lines correspond to constant weighting of the two criteria and the green line corresponds to a varying weight. This way, a flexible cruise control is established where the driver can quickly adapt, for instance, to changing energy requirements.

As has been mentioned before, an alternative to interactively choosing a weight is to implement some heuristic that automatically chooses a weight based on the current situation. Such an approach is visualized in Figure 10c, where the weighting depends on the vehicle velocity, as well as on current and future velocity constraints. For a simpler track, it is possible to compute a globally optimal solution for a scalarized objective using dynamic programming. We see that with the heuristic, the MOMPC approach yields trajectories close to the global optimum while only having finite horizon information.

4.3. Summary

We again conclude the section by giving a summary of publications where multiobjective optimization is applied in a real-time context, cf.

Table 2. Overview of publications (in chronological order) for multiobjective feedback control.

Algorithms without Offline Phase: Computation of Single Points
Kerrigan et al. [131], Wojsznis et al. [154]	Scalarization via Weighted Sum (WS)
Kerrigan and Maciejowski [155], He et al. [137]	Scalarization via lexicographic ordering
Bemporad and Muñoz de la Peña [132,133]	Scalarization via WS for convex objectives, guaranteed stability for large gain vs. noise robust stabilizing objectives
Geisler and Trächtler [134]	WS, online adaptation of weights using gradient information
Maestre et al. [143]	Scalarization via game-theoretic approach
Zavala and Flores-Tlacuahuac [135]	Scalarization via reference point approach
Hackl et al. [129]	Scalarization via WS for Linear Time-Invariant (LTI) systems
Zavala [136]	Scalarization via $ϵ$-constraint: economic objective, stability as constraint
Grüne and Stieler [130]	Economic objectives, performance bounds via selection criterion
Algorithms without Offline Phase: Approximation of the Entire Pareto Set
Laabidi et al. [138,140], Garcìa et al. [141]	ANN for state prediction, optimization via MOEA, selection of Pareto point via WS
Bouani et al. [139]	ANN for state prediction, comparison of two MOEAs and WS for MOP
Nakayama et al. [142]	Few MOEA iterations online, selection via satisficing trade-off method
Algorithms with Offline Phase
Fonseca [145], Herreros et al. [146]	Offline computation of Pareto optimal controller parameters using MOEA
Scherer et al. [156]	Robust control using a common Lyapunov function for multiple stability criteria
Ben Aicha et al. [147]	Offline computation of Pareto optimal controllers parameters via EA and WS, online selection according to geometric mean of objectives
Krüger et al. [148]	Offline computation of Pareto optimal controllers parameters via Set oriented methods, parametric model reduction for increased efficiency
Hernández et al. [149], Xiong et al. [150]	Offline computation of Pareto-optimal controllers parameters via simple cell mapping
Peitz et al. [48]	Offline-online decomposition similar to explicit MPC
Applications
Zambrano and Camacho [157]	MOMPC of a solar refrigeration plant via scalarization
Porfírio et al. [158]	MOMPC of an industrial splitter using a min-max reformulation
Pedersen and Yang [159]	MO PID controller design for magnetic levitation systems via MOEA
Li et al. [160]	Multiobjective adaptive cruise control for vehicles
Hu et al. [161]	MOMPC of high-power converters via WS
Núñez et al. [162]	MOMPC of dynamic pickup and delivery problems using MOEA
Peitz et al. [163]	MOMPC of an industrial laundry, scalarization of a traveling salesman problem via WS

. The publications are divided into four categories. The first two contain algorithms where the MOP is solved online. The categories differ in whether a single point is computed or the entire set is approximated. Consequently an offline phase is not required except in the case where surrogate models are trained in order to accelerate the online computations. The third category then contains the methods with an offline optimization phase, and some applications are mentioned in the fourth category.

5. Reduction Techniques for Many-Objective Optimization Problems

Another important restricting factor in multiobjective optimization is the number of objectives [164]. For MOPs with four or more objectives, the term Many-Objective Optimization (MaOP) has been coined, and over the past few years, many researchers have dedicated their work to address MaOPs and the issues arising from the curse of dimensionality, cf. [165] for an overview, [166,167] for new concepts for identifying non-dominated solutions and [168,169,170,171,172] for evolutionary approaches.

A popular approach for MaOPs is comprised of interactive methods [173,174,175,176,177,178]. These methods do not compute the entire set of optimal compromises, but instead interactively explore the Pareto set. Starting at the current Pareto-optimal solution, a decision maker can choose in which direction to proceed, i.e., which objective to improve at the expense of some other, currently less important objective. The approach in [178], for example, allows for Pareto-optimal movements both in the decision and objective space. One of the main advantages of interactive methods is the reduced computational effort, especially in the presence of many criteria, since it is not affected significantly by the dimension of the Pareto set. Moreover, this way, decision making from a vast number of Pareto-optimal solutions is avoided, which can be overwhelming for a decision maker. Consequently, interpretability and usability are increased.

Besides interactive methods, several reduction techniques have been proposed in the context of many-objective optimization, and although it is not the main theme of this review article, we want to give a brief overview of these reduction approaches since they also aim at increasing the efficiency of solving MOPs. These reduction techniques can be divided into two main categories. The first one is objective reduction, where the aim is to reduce the number of objectives while (approximately) preserving the Pareto set. The observation behind this is that not all objectives are of equal importance to the structure of the Pareto set, which is measured by the degree of conflict [179]. Consequently, when it is possible to identify the main contributors to the Pareto set, then one can solve a reduced problem taking into account only these most important objectives. Different approaches have been proposed for this identification step, all of which use a set of sample points. In [179], both exact and inexact algorithms are proposed for selecting a subset of objectives such that only those points in the Pareto set are lost that are worse in all remaining objectives by a constant $\delta$ or more. This approach is also exploited in [180]. In [181,182], POD (cf. Section 3.3) is used to identify such a subset, and a related concept is implemented in [183] using hyperplanes. In [184], an entropy-based approach is presented, and in [185], the relevant subset is selected multiple times within an evolutionary procedure. A slightly different approach is pursued in [186], where the authors split large decision spaces into several smaller ones according to the relevance of decision variables for specific objectives.

The method proposed in [187] possesses characteristics of the first category, namely objective reduction, as well as of the second category, which is the exploitation of the structure of Pareto sets and front. Therein, first the corners of the Pareto front are identified, and this information is used to select the relevant objectives. Algorithms of the second type all exploit this hierarchical structure. This means that under certain assumptions, the Pareto front is bounded by the Pareto fronts of subproblems where one or more objectives have been neglected, cf. [188,189]. This way, the solution can be computed by a hierarchical approach where, starting with scalar problems, the boundary is computed before, finally, the interior is obtained. Very recently, results about the hierarchical structure of the Pareto set, i.e., in decision space, have appeared; see [73,190] for details. This approach is illustrated in Figure 11 where the solution to an MOP with four objectives is shown, as well as the Pareto sets of the subproblems with three and two objectives, respectively.

6. Future Directions

This survey has given an overview of recent advances in the context of accelerating multiobjective optimization. These are surrogate models, feedback control and objective reduction techniques. Similar to almost every other field of science, it can be expected that the immense developments in data-based methods will also have a major impact on research in multiobjective optimization, in particular in the context of surrogate modeling. A very large number of researchers from the dynamical systems community is working on data-based methods using the Koopman operator, which is an infinite-dimensional, but linear operator describing the dynamics of observables [191,192]. Significant effort has been put into the development of numerical methods for approximating this operator from data; see, e.g., [193,194]. This way, the dynamics of observations can be reconstructed entirely from data and without any knowledge of the underlying system dynamics. In a way, this allows us to merge the two surrogate modeling categories from Section 3.2 and Section 3.3 since we can approximate the dynamics not only of the state, but directly of the objectives. Several methods have recently been proposed to use the Koopman operator for data-based controller design, both in simulations [195,196,197,198], as well as experiments [199,200]. The results are very promising, such that it is just a matter of time until these methods are utilized for multiobjective optimal control.

In the same manner, machine learning techniques [201] will very likely gain more and more attention, both in multiobjective optimization, as well as optimal control. There are already many papers on this topic or related ones, and the number is growing quickly.