A survey of recent trends in multiobjective optimal control – surrogate models, feedback control and objective reduction

: 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 conﬂicting 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 inﬁnite 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 over 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 ﬁrst case, surrogate models yield signiﬁcant 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 ofﬂine-online decomposition. We also brieﬂy comment on dimension reduction in many-objective optimization problems as another technique for reducing the numerical effort.


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 a rather function, typically depending on time.The goal is to steer a dynamical system in such a way that one (ore multiple) objectives are 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 over 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 which 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.

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

Theory
In multiobjective optimization we want to minimize multiple objectives at the same time.Consequently, the fundamental difference to scalar optimization is that the objective function J : U → R k is vector-valued.Hence, the general problem is of the form min u∈U J(u) = min u∈U    J 1 (u) . . .
where u ∈ U is the control variable and g : U → R l , g(u) = (g 1 (u), . . ., g l (u)) and h : U → R m , h(u) = (h 1 (u), . . ., h m (u)) , are inequality and equality constraints, respectively.The space of the control variables 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, U can either be finite-dimensional, i.e.U = 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): ẏ(x, t) = G(y(x, t), u(t)), (x, t) ∈ Ω × (t 0 , t e ], a(x, t) ∂y ∂n (x, t) + b(x, t)y(x, t) = c(x, t), (x, t) ∈ Γ × (t 0 , t e ], y(x, t 0 ) = y 0 (x), x ∈ Ω, where the domain of interest Ω ⊂ R n x is a connected open set with spatial dimension n x and the boundary is denoted by Γ = ∂Ω 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 Ĵ : U × Y → 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: . . .

(MOCP)
There are articles on multiobjective optimal control which 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 S : U → Y, see [9] for details.By setting Ĵ(u, y) = Ĵ(u, Su) = J(u), the problem is transformed into (MOP).For this reason, we will from now on only consider (MOP).
In the situation where 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 R k with k ≥ 2 (unless they are not conflicting).Therefore, the comparison of values is defined in the following way [4]: 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 * dominates a point u, if J(u * ) ≤ J(u) and J(u * ) = J(u).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 * be a Pareto optimal point of problem (MOP) and assume that ∇h j (u * ) for j = 1, . . ., m, and ∇g s (u * ) for s = 1, . . ., l, are linearly independent.Then there exist non-negative scalars α 1 , . . . , The set of points satisfying these conditions is called the set of substationary points P S,sub [4].Obviously, P S,sub is a superset of the Pareto set 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 cheap) non-dominance test.

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 descend 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: Here, α is a fixed weight vector which is determined in such a way that see, e.g., [31,32].In the unconstrained case, there exists no u satisfying (2) only if q(u) = 0 which implies that u is substationary, cf.(KKT).

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 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 few problems have been addressed directly, see, e.g., [50][51][52][53].A method exploiting special structure in the system dynamics has been proposed in [54], and a special case of Pareto optimal solutions -namely Nash equilibria -have been computed in [55,56]  1 .
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(u) is replaced by a surrogate J r (u), 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 : U → 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: where the reduced control-to-state operator 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 (u) = J(u) holds for all u ∈ 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], 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 over the existing methods and the corresponding references.

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 (u), ∇J r (u) of the objectives J i (u) and the gradients ∇J i (u), i = 1, . . ., k, respectively.Furthermore, let us assume that upper bounds , κ ∈ R k for these errors are known: 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(u) is only known approximately according to (3).Then When not using a priori selection methods such as scalarization or Nash equilibria, a decision maker to select the appropriate solution.This is called Multicriteria Decision Making (MCDM) [57] and is an entire area of research of its own.Thus, we will not go into further details about the decision making process here.
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) to (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 2 (a) 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 (u).The inaccuracy in the gradients introduces an upper bound for the angle between the individual gradients ∇J r i (u) and the descent direction q r (u).This is equivalent to a lower bound αmin for the weight vector, i.e., αi ∈ [0, 1] in Equation ( 1) has to be replaced by αr see [66] for a detailed discussion.The additional constraint ensures that if q r (u) 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 P S,sub based on inexact gradient information:  ).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 q r (u) 2 = 0 and ∇J r i (u) 2 = 0 for i = 1, . . ., k.Let Then the following statements are valid: (a) If ∑ k i=1 αmin,i > 1 then there exists no direction q(u) with i.e. no descent direction for the exact problem.(b) All points u with ∑ k i=1 αmin,i = 1 are contained in the set 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, q(u * ) = 0 holds for all u * ∈ P S,sub .The set P S,sub is shown in red in Figure 4 (a).Due to the inexactness, we can no longer guarantee q(u * ) = 0. Instead, we obtain the set P S,κ which is shown in green.The background is colored according to the optimality condition q(u) of the exact problem and the dashed white line indicates the error bound (7) from the above theorem.

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 which 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 get by.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 which are contained in the set U ref = {u 1 , . . ., u p }, where p is the number of function evaluations (or experiments).The data points J(u 1 ) to J(u p ) 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 data set have to be answered: 1. how large does the set U ref have to be? 2. how can we pick the correct elements for U ref ? 3. do we define U 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 to 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 on page 17 mentions both survey articles and a (non-exhaustive) selection of popular methods.

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 S by a cheaper surrogate model 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 (u) = Ĵr (u, y) = Ĵ(u, S r u) is also much cheaper 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).

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 Y is transformed into a finite-dimensional space Y N via Galerkin projection: Here, N denotes the number of degrees of freedom and φ i (x) 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 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 Y r with dimension 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 data set {y N 1 , . . ., y N s } to construct a reduced basis ψ = {ψ 1 , . . ., ψ } such that y(x, t) ≈ y N (x, t) ≈ y r (x, t) = ∑ i=1 z i (t)ψ i (x).
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 Σ: 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: 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: Introducing the mass matrix M ∈ R N×N and the stiffness matrix this yields the following N-dimensional linear system: 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 N 1 , . . ., y N s } are solutions corresponding to parameters {u 1 , . . ., u s } and in POD they are snapshots in time, collected at the time instants {t 0 , . . ., t s−1 }.Using an equidistant time grid h, the snapshots are taken at {t 0 , t 0 + h . . ., t 0 + (s − 1)h}.In RB, the snapshots {y N 1 , . . ., y N s } often directly serve as the basis ψ.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 N 1 , . . ., y N s ) 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 ψ.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: 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 y r − y N or -in the case of optimal control -of the optimal solution u r − u N 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.

Optimal control using surrogate models
There is a rich literature on optimal control of PDEs using surrogate models.The approaches can be summarized in three main categories: 1. build a model once, 2. construction of regular updates in a trust region approach, 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 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. u i − u ref ≤ δ i , where i is the current step of the iterative optimization scheme and δ 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: If ρ is close to 1, 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 ref = u i .If, on the other hand, ρ is close to 0, then the ROM accuracy was bad and the iterate u i is rejected.Instead, the trust region radius δ 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 case of the Navier-Stokes equations, this has been shown for different problem setups [103,104].In the third approach, error estimators ∆ J for the current iterate u i are required.By evaluating ∆ J (u i ) it is possible to efficiently estimate the error between the high-and the low-dimensional solution, since 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].

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.

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 α: 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(u) has to be minimized: 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 7 (a)).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 7 (b)-(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.

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 which 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 is cannot be achieved with a single ROM, one has to use multiple, locally valid ROMS instead (cf. Figure 8 (c)).The covering by local ROMs is managed in such a way that all points in the neighborhood around the reference u ref at which the data for the ROM was 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.

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.

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.

Surveys
Tabatabei et al. [2], Chugh et al. [3] Extensive surveys on meta modeling for MOEAs Voutchkov & Keane [74], Surveys on meta modeling approaches from statistics (RSM, Knowles & Nakayama [1], Jin [75] RBF) and machine learning in combination with MOEAs Benner et al. [83], Taira et al. [85], Surveys on reduced order modeling of dynamical systems Peherstorfer et al. [84] Algorithms using meta models for the objective function Ong et al. [112], Ray et al. [113] Combination of RBF and MOEA Chung & Alonso [114], Keane [115] Combination of Kriging models and MOEA Karakasis & Giannakoglou [116] RBF as a cheap 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 In MPC, an optimal control problem is solved on 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.

Preprints
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 2 .
When implementing a 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 which is determined in the first MPC step.Due to this, upper bounds for the objective function values can be guaranteed.

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: 1. compute a single Pareto optimal solution according to some predefined preference, 2. compute only a rough approximation of the Pareto set, 3. compute an arbitrary Pareto optimal control which 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 which one can build on.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 2 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.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 proved for the resulting MPC algorithm.
In [134], this approach is extended by providing gradient information of the objectives with respect to the weight vector α.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 * , 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 lexocographic 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 a MOEA and stop the computation after 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] and [143], where a game theoretic approach is used.

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 which fall in 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 to solve 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 which 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 which 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 min 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 α.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.

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 10 (a) 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 to 100 km/h, cf.scenario (I I) in Figure 10 (a).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(t 0 ) = . . ., 59.9, 60.0, 60.1, . ..), 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 10 (b) for an example track, where the black lines correspond to constant weighting of the two criteria and the

Algorithms with offline phase
Fonseca [145], Offline computation of Pareto optimal controller parameters Herreros et al. [146] 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] Offline computation of Pareto optimal controllers parameters via Xiong et al. [150] simple cell mapping As has been mentioned before, an alternative to interactively choosing a weight is to implement some heuristic which automatically chooses a weight based on the current situation.Such an approach is visualized in Figure 10 (c), 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.

Summary
We again conclude the section by giving a summary of publications where multiobjective optimization is applied in a real-time context.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.

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 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 are 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 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 which are worse in all remaining objectives by a constant δ or more.This approach is also exploited in [180].In [181,182], POD (cf.Section 3.3) is used to identify such a subset, 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 on 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 in 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.

Future directions
This survey has given an overview over 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 is to merge the two surrogate modeling categories from Sections 3.2 and 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.

Figure 1 .
Figure 1.The red lines are the Pareto set (a) and Pareto front (b) of an exemplary multiobjective optimization problem (two paraboloids) of the form min u∈R J(u), J : R 2 → R 2 .The point J * = (0, 0) is called the utopian point.

Figure 2 .
Figure 2. Example for the -dominance property.A point-wise comparison is illustrated in (a) -(d).The uncertainties are marked by the dashed boxes.Only in case (d) the lower left point confidently dominates the other point.In (e) the entire Pareto fronts for the exact problem (P F ) and the inexact problem (P F, ) are shown in red and orange, respectively.

Figure 3 .
Figure 3. Example for an MOP from production [32] where inexactness is introduced due to uncertainties in pricing.The Pareto set P S for the exact problem is shown in (a), the inexact set P S, is shown in (b).

Figure 4 .
Figure 4. Exact and inexact solution (P S and P S,κ ) for a simple example with J : R 2 → R 2 , cf. [66] for details.(a) The sets P S and P S,κ (for a random error with κ = (0.01, 0.01) ) are shown in red and green, respectively.The background is colored according to the optimality condition q(u) which has to be 0 for all substationary points.The dashed white line shows the error bound as derived in Theorem 2. (b) The corresponding Pareto fronts.

Figure 5 .
Figure 5. Trust region method.(a) The ROM-based optimal control problem is solved within the trust region δ 0 .(b) If the improvement is poor for the full system (i.e., ρ is small), then the trust region radius is reduced and we repeat the computation with the same problem.(c) If the improvement is acceptable (intermediate values of ρ), then we compute a new model and proceed with a smaller trust region δ 1 < δ 0 .(d) If the improvement is good (i.e., ρ ≈ 1), then the trust region radius is increased.

PreprintsFigure 6 .
Figure 6.Reference point method.(a) Determination of a Pareto optimal solution by solving (13).(b) Determination of the consecutive point on the Pareto front by adjusting the target and solving the next scalar problem.

Figure 7 .
Figure 7. (a) Pareto front for an MOCP involving the Navier-Stokes equations (flow stabilization vs. cost), solved by coupling of a ROM (created once in advance) with the reference point method.Although we observe acceptable agreement, convergence cannot be guaranteed.(b) Pareto front for a heat flow MOCP (reference tracking vs. cost), solved by the TR-POD approach coupled with the reference point method.Convergence is achieved while reducing the number of expensive finite element (FEM) evaluations by a factor of ≈ 22, cf.(c).

Figure 8 .
Figure 8.(a) Pareto set of a semilinear heat flow MOP with four controls (coloring according to u 4 ), solved directly including an FEM model in the subdivision algorithm.(b) Pareto set of the same problem, solved with localized ROMs.(c) The reference controls for which the local ROMs have been computed are shown in black, the colored dots are sample points at which the objective function was evaluated.The colorings denote assignments to a specific ROM.(d) The corresponding Pareto fronts, where the FEM solution is shown in green and the ROM solution in red.

Figure 9 .
Figure 9. Sketch of the MPC method.Due to the real-time constraints, the optimization problem has to be solved faster than the sample time h.

Figure 10 .
Figure10.Results for the offline-online MPC approach from[48].(a) Different constraint scenarios ((I) to (V I)), i.e., constant velocity, acceleration, deceleration, stopping.(b) Example track driven with the MPC algorithm.The red lines define the velocity bounds, the black dashed lines are trajectories corresponding to a constant weight α, and the green line is a trajectory where the weight is changed from 0 (energy efficient) over 0.5 (average) to 1 (fast).(c) Comparison between the MPC algorithm (coupled with a simple heuristic for the weighting) and the global optimum obtained via dynamic programming.

PreprintsFigure 11 .
Figure 11.Visualization of the hierarchical structure of Pareto sets.(a) Pareto set of an example problem with J : R 3 → R 4 .(b) The four Pareto sets taking only three objectives into account form the boundary of the original Pareto set.(c) The Pareto sets in (b) are again bounded by the respective bi-objective subproblems.
Acknowledgement: This work is supported by the Priority Programme SPP 1962 "Non-smooth and Complementarity-based Distributed Parameter Systems" of the German Research Foundation (DFG).

Table 1 .
Overview of publications (in chronological order) where surrogate modeling and multiobjective optimization are combined.

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