Multiobjective Model Predictive Control of a Parabolic Advection-Diffusion-Reaction Equation

: In the present paper, a multiobjective optimal control problem governed by a linear parabolic advection-diffusion-reaction equation is considered. The optimal controls are computed by applying model predictive control (MPC), which is a method for controlling dynamical systems over long or inﬁnite time horizons by successively computing optimal controls over a moving ﬁnite time horizon. Numerical experiments illustrate that the proposed solution approach can be successfully applied although some of the assumptions which are necessary to conduct the theoretical analysis cannot be guaranteed for the studied tests.


Introduction
Model predictive control (MPC) is a method for controlling dynamical systems over long or infinite time horizons by successively computing optimal controls over a moving finite time horizon; cf., for example, Reference [1]. The successive re-optimization on the one hand introduces a feedback mechanism, which makes the method attractive as a real-time control scheme and as such it is widely used in industry [2]. On the other hand, MPC provides a method to reduce the complexity of optimal control problems by splitting up problems on long or infinite horizons into smaller subproblems over shorter horizons. As such, MPC can be seen as a model reduction technique in time. Clearly, this approach only makes sense if the solution generated by MPC is in some sense close to the true optimal solution. Fortunately, this approximate optimality property can be rigorously proved for many classes of optimal control problems, see, for example, References [3,4] and the references therein.
In practice, it is often the case that a single optimization criterion is not sufficient for modelling the demands in a given application. This leads to the concept of multiobjective optimal control, in which typically there does not exist a single optimal value but rather a whole set of optimal values, the so called Pareto front. Due to its usefulness in practice it is no surprise that this concept has been used and investigated also in the MPC context, see, for example, References [5][6][7][8][9]. Particularly, in References [10,11] MPC algorithms were presented that allow for rigorous suboptimality estimates also in the multiobjective case. The main feature of this class of algorithms is a particular constraint that depends on the chosen point on the Pareto front when solving the first optimal control problem in the MPC scheme and that is applied in all subsequent optimal control problems. While the examples in References [10,11] are limited to finite dimensional systems, in this paper we show that this idea can be successfully applied also to optimal control problems governed by partial differential equations (PDEs) and that the necessary constraint can be efficiently implemented using a gradient descent scheme. While we will see that some of the assumptions made in References [10,11] are difficult to a.e. on Σ, a.e. in Ω. (1) The functions b i ∈ L ∞ (Γ), 1 ≤ i ≤ m, denote given control shape functions, for example, characteristic functions: By setting H = L 2 (Ω) and V = H 1 (Ω), one can show that (1) has a unique weak solution y = y(u, y • ) ∈ W(0, T) = H 1 (0, T; V ) ∩ L 2 (0, T; V) for all controls u ∈ L 2 (0, T; R m ) =: U and all initial conditions y • ∈ H. Recall that W(0, T) is a Hilbert space endowed with its common inner product. Moreover, W(0, T) is continuously embedded into C([0, T]; H); cf. Reference [24]. The solution operator S : U × H → W(0, T) → L 2 (0, T; H), which maps any control u ∈ U and initial value y • ∈ H to the solution y = S(u, y • ) of (1), is affine linear and continuous in both components.
is called the Pareto front of (MOCP).

Multiobjective Model Predictive Control (Mompc)
In this section we explain how MPC is combined with multiobjective optimal control. Given a final (large) time horizon T > 0 we use an equidistant time grid The open-loop problems in the MPC algorithm are then solved on smaller time horizons ∆t with 1 < < L. We set By the MPC method we will finally compute an optimal control on the time horizon [0, t end ]. For n = 0, . . . , L − the optimal control to this problem on the interval (t n , t n + ∆t) = (t n , t n+1 ) is stored and used to compute the next part of the state trajectory on the interval (t n , t n+1 ). Then the time horizon is shifted by ∆t and the procedure is repeated until the final time horizon is reached. Notice that the end time for the computed MPC control is t end .
To deal with this framework we have to consider (1) on subintervals of [0, T]. To this end we introduce further notations: For an arbitrary initial time t a ∈ [0, T) and end time t b ∈ (t a , T], we study the state equation on the time horizon I = (t a , t b ): for some initial value y a ∈ H. Again it is possible to show that for each u ∈ U I = L 2 (I; R m ) and y a ∈ H there exists a unique weak solution y = y(u, y a ) ∈ W(I) to (3) and that the associated solution operator S I : U I × H → W(I) → L 2 (I; H) is affine linear and continuous in both components. The two objectivesĴ I 1 ,Ĵ I 2 : U I × H → R are given bŷ for i = 1, 2. Note that we would actually have to write C 1 I i ∈ L 2 (I; H) for i = 1, 2 to display the dependence of the operators on the time interval I. Since we define these operators in a canonical way by where E : L 2 (I; H) → L 2 (0, T; H) extends any function f ∈ L 2 (I; H) by 0 to a function in L 2 (0, T; H), we omit this dependence on I for the purpose of readability. Throughout we denote by U I opt (y a ) and J I opt (y a ) the Pareto set and Pareto front, respectively, of the MOCP min u∈U IĴ I (u, y a ) = min (MOCP I ) Especially for the open-loop problems in the MPC algorithm we make the following notations: Given n ∈ {0, . . . , L − } we set t n = t n + ∆t ≤ T for the final time instance and U n = L 2 (t n , t n ; R m ).
Furthermore, for brevity we define S n = S (t n ,t n ) ,Ĵ n i =Ĵ (t n ,t n ) i for i = 1, 2, and U n opt (y a ) = U (t n ,t n ) opt (y a ) as well as J n opt (y a ) = J (t n ,t n ) opt (y a ). In the context of MOMPC the notion of external stability turns out to be important, see Reference ( [10], Definition 6) and Theorem 2 below.
be a time interval and the initial condition y a ∈ H be arbitrary. The set J I (U I , y a ) is called externally stable, if for every y ∈Ĵ I (U I , y a ) there isȳ ∈ J I opt (y a ) withȳ ≤ y. This is equivalent toĴ I (U I , y a ) ⊂ J I opt (y a ) + R 2 ≥0 .
For the cost functions defined in (4) it is possible to show that the setĴ I (U I , y a ) is externally stable, as the following theorem shows.
be a time interval and the initial condition y a ∈ H be arbitrarily given. We assume that σ 3 i > 0 for i = 1, 2. Then the setĴ I (U I , y a ) is externally stable.
Proof. It is possible to show that the Pareto optimal points ofĴ I (U I , y a ) andĴ I (U I , y a ) + R 2 ≥0 are the same, that is, it holds J I opt (y a ) = (Ĵ I (U I , y a ) + R 2 ≥0 ) opt . So if we show that the setĴ I (U I , y a ) + R 2 ≥0 is externally stable, that is, that it holdsĴ I (U I , y a ) + R 2 ≥0 ⊂ (Ĵ I (U I , y a ) + R 2 ≥0 ) opt + R 2 ≥0 , we are done, since then we havê According to Reference ( [12], Theorem 2.21) we only need to show thatĴ I (U I , y a ) + R 2 ≥0 = ∅ and that for all y ∈Ĵ I (U I , y a ) + R 2 ≥0 the set (Ĵ I (U I , y a ) + R 2 ≥0 ) ∩ (y − R 2 ≥0 ) is compact (this is the notion of R 2 ≥0 -compactness). It is clear thatĴ I (U I , y a ) = ∅ holds. Since the setĴ I (U I , y a ) + R 2 ≥0 is bounded from below and y − R 2 ≥0 is bounded from above, the set (Ĵ I (U I , y a ) + R 2 ≥0 ) ∩ (y − R 2 ≥0 ) is bounded. The closedness ofĴ I (U I , y a ) + R 2 ≥0 follows from the convexity, the continuity and the property lim which follows from σ 3 i > 0 for i = 1, 2, of the cost functionsĴ I 1 (·, y a ),Ĵ I 2 (·, y a ), and was for example shown in the proof of Theorem 3.35 in Reference [22].
The pseudocode for a general MOMPC algorithm is shown in Algorithm 1. Here we use a fixed time horizon ∆t on which the computed control is applied before a re-optimization is performed. We note that it may be possible to use varying time horizons here, similar as discussed for the single-objective case, for example, in Reference ( [4], Section 10.4).

5:
Compute a Pareto optimal controlū n ∈ U n opt (y a ).

8: end for
The most important question is how to perform step 5 of Algorithm 1. Clearly, one possibility is to compute the entire Pareto set U n opt (y a ), which is nonempty by Theorem 1, and then to choose a control u n ∈ U n opt (y a ) according to some given preferences. However, this procedure is computationally costly, and thus, infeasible for many applications. On the other hand, one could compute only one Pareto optimal controlū n ∈ U n opt (y a ). Then the question is under which criteria this control should be computed. In Reference [10], the authors allowed an arbitrary choice forū 0 ∈ U 0 opt (y • ) and used the additional constraintĴ n i (ū n , y a ) ≤Ĵ n i (ũ, y a ), for i = 1, 2 and y a = y n (t n ), for all n ∈ {1, . . . , L − }, whereũ is given bỹ and the feedback κ ∈ U (t n−1 ,t n ) is chosen such that holds for i = 1, 2, with y κ = S n−1 (ū n−1 , y n−1 (t n−1 ))(t n−1 ). The external stability of the set J n (U n , y n (t n )) ensures that such a controlū n exists.
Proof. This follows directly from Theorem 1, where it was shown that the setĴ n (U n , y a ) is externally stable.
However, the existence of a feedback κ ∈ U (t n−1 ,t n ) fulfilling (7) cannot be shown in general for our problem, but has to be assumed. Such a feedback exists if and only if the minimal value of the minimization problem is smaller or equal than zero. Thus, by solving the minimization problem (8) we can on the one hand check, if such a feedback exists, and on the other hand compute it explicitly. If the conditions (5) and (7) are fulfilled throughout Algorithm 1 we can show the following performance theorem, compare Reference ( [10], Theorem 11), where also the proof was taken from and adapted to our situation. Theorem 3. Let y • ∈ H be arbitrary and assume that σ 3 i > 0 holds for i = 1, 2. Denote byū 0 ∈ U 0 opt (y • ) the initially chosen Pareto optimal control in Algorithm 1. Furthermore, assume that it is possible to choose a feedback κ ∈ U (t n−1 ,t n ) fulfilling (7) for every n ∈ {1, . . . , L − }. Let µ ∈ U (0,t end ) be the resulting MOMPC feedback control resulting from Algorithm 1. Then it holdŝ Proof. Note that it holds µ(t) =ū n (t) for almost all t ∈ (t n , t n+1 ) and all n ∈ {0, . . . , L − } and that Since we assume the existence of a feedback κ fulfilling (7), we can conclude for all n ∈ {0, . . . , L − } withũ n+1 ∈ U (t n+1 ,t n+1 ) defined as in (6). Due to (6) we find that Moreover, it follows that Plugging (10) into (9) and utilizing (11) and (12) we deduce that Since the controlū n+1 fulfills (5) for n ∈ {0, . . . , L − − 1}, we can further estimate The first two sums on the right-hand side are a telescopic sum. Together with the nonnegativity of the norms, this finally implies which is what we wanted to show.

Remark 1.
(1) The statement from Theorem 3 is important, since it gives us performance bounds on the cost functions values of the MOMPC feedback control µ already after choosing the initial controlū 0 ∈ U 0 opt (y • ), but before performing the MOMPC Algorithm 1. Thus, one strategy is to compute the entire Pareto set U 0 opt (y • ), which is computationally cheap due to the small time horizon (t 0 , t 0 ), and then to choose the initial control u 0 ∈ U 0 opt (y • ) according to the desired upper bounds on the cost functions.
Theorem 3 holds for arbitrary T > 0. In particular, by taking the limit T → ∞, the result can also be shown for the infinite-horizon case.
Now we want to present a gradient descent scheme for computing a controlū n ∈ U n opt (y a ) fulfilling (5). To this end we use the following result which is taken from for example, Reference ( [25], Theorem 2.1), where it is proved for the case of finite-dimensional controls. However, the extension to our infinite-dimensional setting U n = L 2 (t n , t n ; R m ) is straight-forward, since the proof of Theorem 2.1 in Reference [25] does not explicitly use the finite-dimensionality of the space U n , but only its property of being a Hilbert space.
Suppose that y a ∈ H, u ∈ U n are arbitrarily given and Then either is a descent direction for the cost functionsĴ n 1 (· , y a ) andĴ n 2 (· , y a ), or q n (u) = 0 holds. Since the cost functionŝ J n 1 (· , y a ),Ĵ n 2 (· , y a ) are strictly convex, the latter implies that u is a Pareto optimal point, c.f. Reference [12]. (13) is a quadratic optimization problem in the variable α ∈ R 2 with the constraints α ∈ [0, 1] 2 and α 1 + α 2 = 1. By substituting α 2 = 1 − α 1 , this problem can be reformulated as a box-constrained quadratic optimization problem in only one variable, which is easy to solve. In our numerical experiments we utilize the MATLAB routine quadprog.

Remark 2. Problem
With Theorem 4 we can set up Algorithm 2, which takes an arbitrary point u ∈ U n as an input and results in a Pareto optimal pointū, compare Reference ( [26], Algorithm 1).
From Reference ( [26] Lemma 4), it follows that a step length t k > 0 can always be found in step 5 of Algorithm 2. Moreover, in Reference ( [26], Theorem 1) the authors proved the following convergence result. Again, they assumed a finite-dimensional space U n , but the arguments of the proof transfer directly to our setting U n = L 2 (t n , t n ; R m ). Theorem 5. Every accumulation point of the sequence (u k ) k∈N produced by Algorithm 2 is a Pareto optimal point.

Remark 3.
(1) Note that we cannot prove that the sequence (u k ) k∈N has an accumulation point in the infinite-dimensional case. However, we will not encounter this problem in our numerical implementation, since the space U n will be discretized. Therefore, Algorithm 2 will in practice terminate in a finite number of steps.
By construction of the algorithm it holdsĴ n i (ū n , y a ) ≤Ĵ n i (u, y a ), i = 1, 2, for any initial control u ∈ U n , so that (5) is fulfilled, if we choose u =ũ.
Finally, the MOMPC algorithm looks as follows.

Algorithm 3 Multiobjective Model Predictive Control
Given: Initial state y • ∈ H. 1: for n = 0 to L do 2: Set y a = y n (t n ).

Numerical Tests
Throughout the numerical experiments let the domain Ω be given by Ω = (0, 1) 2 with points x = (x 1 , x 2 ). This domain is discretized by using linear finite elements with 494 degrees of freedom. We choose t end = 1 and ∆t = 1/198. The number of total time steps L is set to L = 197 + in dependence of the smaller time horizons of the MPC open-loop problems. Accordingly, the final time T is given by T = L∆t. For the time integration of the PDE we use the Crank-Nicolson method. In Algorithm 2 we set the tolerance ε = 10 −4 . In Algorithm 3 the initial controlū 0 is chosen as an element of the Pareto set U 0 opt (y • ). Therefore, in the following we will always compute a finite approximation U 0 opt,appr (y • ) of this set, which depends on the MPC horizon length , by the Euclidean reference point method.
We note that in all our examples we use fixed temporal and spatial discretizations. Clearly, this is not necessarily efficient and one may rather prefer to use adaptive discretizations that only use fine grids in space and/or time when this is really relevant for obtaining a sufficiently accurate numerical approximation for the control function. For single-objective PDE governed MPC problems, such methods were developed in, for example, References [27,28]. We conjecture that they could be adapted to the multiobjective setting.
Given these data we can apply Theorem 1 and show the external stability of the setsĴ n (U n , y n (t n )), which implies the feasibility of the steps 4 and 6 in Algorithm 3.
However, we can in general not expect that there is a feedback κ fulfilling (7). The reason for this is that such a κ would have to fulfill the inequality due to the choices σ 1 2 = σ 2 2 = 0 in the second cost function, which implies κ = 0. Plugging this into (7) for i = 1 we see that the inequality is only fulfilled, if the uncontrolled system would move towards the desired temperatures. This is unlikely to happen in our setting, since the desired temperatures are larger than the initial condition and increasing in time. Nevertheless we can still compute a minimizer of (8) and use it to defineũ in (6). Note that in this case, the assumptions of Theorem 3 are not fulfilled, so that this choice of κ andũ is of heuristic nature.
Another heuristic approach for determiningũ that we want to test in the following is motivated by the criteria used in step 2 of Algorithm 4 in Reference [11], which was designed for problems without terminal condition and translates to our problem as follows: During the n-th loop iteration of Algorithm 3 (n ∈ {0, . . . , L − }) compute κ ∈ U (t n−1 ,t n+1 ) opt (y n (t n−1 )) in step 9 such that κ 2 L 2 (t n−1 ,t n ;R 4 ) ≤ ū n 2 L 2 (t n−1 ,t n ;R 4 ) holds. Then we setũ Here we do not impose the inequality (15) on κ explicitly. The reason for this is that demanding (15) would not guarantee us any performance results for our framework, but only increase the computational time.
Therefore, we use again the gradient descent method presented in Algorithm 2 for computing a control κ ∈ U (t n−1 ,t n+1 ) opt . As initial control u we choose Again, it can be shown that for any accumulation pointū of the sequence (u k ) k∈N produced by Algorithm 2, it holdsū ∈ U (t n−1 ,t n+1 ) opt (y n (t n−1 )). Moreover, although the inequality (15) is not guaranteed directly by this method, we still getĴ (t n−1 ,t n+1 ) (κ, y n (t n−1 )) ≤Ĵ (t n−1 ,t n+1 ) (u, y n (t n−1 )).

2.
Compute κ by using Algorithm 2 as described above and setũ according to (16).
The obtained results when using strategy 1 for computing κ are shown in the top row of Figure 1a-c. It is clearly visible that an increasing MPC horizon length has the positive effect that the MOMPC points produced by Algorithm 3 are located closer to the Pareto front J (0,1) opt (y • ). This can be explained by the fact that a larger MPC horizon allows for a better prediction of the future behavior of the system dynamics. However, we see a clustering of MOMPC points in the middle of the Pareto front, which only slightly improves, if the MPC horizon length is increased. So it is not possible to obtain the whole extent of the Pareto front by varying the initial controlū 0 ∈ U 0 opt,appr (y • ). On the other hand, we can see the results of strategy 2 in the bottom row of Figure 1a-c. Let us first note that the inequality (15) is fulfilled in most steps of Algorithm 3 for all initial controls u 0 ∈ U 0 opt,appr (y • ). Only for some initial controlsū 0 , for whichĴ 0 (ū 0 , y • ) is located on the top left part of the Pareto front J 0 opt (y • ), that is, for initial controls corresponding to the part of the Pareto front, where the main goal is to minimizeĴ 0 1 (·, y • ) almost regardless of the function values ofĴ 0 2 (·, y • ), the condition (15) is not fulfilled in some steps of the MOMPC algorithm. The reason for this is the following: When using Algorithm 2 for computing a feedback κ ∈ U (t n−1 ,t n+1 ) opt , the resulting minimizer α of (13) displays the weight of the cost functionsĴ (t n−1 ,t n+1 ) (·, y n (t n−1 )) at the Pareto optimal point κ ∈ U (t n−1 ,t n+1 ) opt (y n (t n−1 )). If we start with a pointū 0 , for whichĴ 0 (ū 0 , y • ) is located on the top left part of the Pareto front J 0 opt (y • ), the weighting of the cost functions will beα ≈ (1, 0) . Therefore, the descent direction q n (u) from Algorithm 2 will mostly point in the direction of the negative gradient of the first cost function, which might lead to (15) not being satisfied.
Looking at the results, we observe again that a larger MPC horizon length leads to a better result in the sense that the MOMPC points are closer to the Pareto front. In contrast to strategy 1, we do not see a clustering of MOMPC points in the middle of the Pareto front, but rather in the lower right part of the Pareto front. This clustering improves, if the MPC horizon length is increased, so that the entire scale of the Pareto front is obtained for a horizon length of = 99. The difference in the clustering behavior of strategy 1 compared to strategy 2 can be seen by looking at the following relation: Every initial control vectorū 0 ∈ U 0 opt,appr (y • ) is the solution to a weighted-sum problem for some weight vector α init ∈ R 2 ≥0 with α init 1 + α init 2 = 1; see, for example, Reference [12]. This weight vector can be determined easily when using the Euclidean reference point method to compute U 0 opt,appr (y • ); cf. Reference ( [19], Lemma 5). After executing Algorithm 3 we can check to which Pareto optimal pointȳ =Ĵ (0,1) (ū, y • ) ∈ J (0,1) opt (y • ) the MOMPC pointĴ (0,1) (µ, y • ) has the smallest distance. Again,ū is the solution to a weighted-sum problem min u∈U (0,1) for a weight vector α end ∈ R 2 ≥0 with α end 1 + α end 2 = 1, which can again be determined. The mapping can be seen in Figure 2 for both strategies. Since the mapping displays the first weight vector, a small value in the plot corresponds to a small weighting of the first cost function (the tracking term) and a large weighting of the second cost function (the control costs), and vice versa.
Note that the 'ideal' result of this mapping would be the identity, since this would imply that the weight α init of the initial controlū 0 remains constant throughout the MOMPC algorithm.
From this plot the clustering of the MOMPC points for both strategies can be deduced. For strategy 1 it can be seen that the clustering of the points happen at a weight of α end 1 ∈ [0.7, 0.8] for all MPC horizon lengths. Up to an initial weight of α init 1 ≈ 0.8 all the initial controls lead more or less to the same result of the MOMPC algorithm. Only for initial controls with an initial weight of α init 1 > 0.8, we observe that the upper part of the Pareto front can be reached. (·, y • ). Moreover, for these two MPC horizons, we can see a clear cut-off behavior at a value of 0.8. If α init 1 is larger than 0.8, the value of α end 1 is larger than 0.9. For an MPC horizon length of = 99, we do not observe such a clear cut-off behavior. Although the plot is still far from being the identity, it is visible that one can control the outcome of the MOMPC algorithm more precisely by varying the initial control.
In conclusion we can say that for strategy 2 a larger MPC horizon leads to a better distribution of the MOMPC points on the Pareto front, and that it is possible to approximate the full scale of the Pareto front by choosing different initial controls u 0 ∈ U 0 opt (y • ). Strategy 1, however, seems not to be well-suited for this problem framework, since it is not possible to reach to whole extent of the Pareto front by choosing different initial controls. Quite on the contrary, the MOMPC algorithm steers almost all initial controls to the same region of the Pareto front, independently of the MPC horizon length. Test 2. Now we want to show why the use of MPC in the multiobjective context is needed. For this we consider the following setup: Imagine that we want to compute the Pareto front J (0,1) opt (y • ). However, there is only a prediction of the advection field available, which is given by In the end it turns out that the prediction of the advection field is not very accurate, and the true advection is given by the function c in (14), that is, it deviates from the prediction c pred after t = 0.5.
Denote by U  opt,pred (y • ), y • ). One can clearly see that the controls in U (0,1) opt,pred (y • ) are far from being Pareto optimal for the problem with advection c, especially in the upper part of the Pareto front. This can be explained by the fact that the changing advection totally changes the strategy that has to be used to come close to the desired temperature, which is not captured by the open-loop problem using the prediction c pred . We compare this to the results obtained by using the MOMPC Algorithm 3 together with strategy 2 from the first test, where we assumed that the true advection is known to the open-loop problems in the MPC algorithm. This is reasonable since the MPC horizons are reasonably small to allow a precise prediction of the future behavior of the advection field. Just as an example we look at an MPC horizon length of = 50 in Figure 3b and see that in contrast to the predicted Pareto optimal controls U (0,1) opt,pred (y • ) the MOMPC feedback control µ comes quite close to the true Pareto front. This underlines that the use of MOMPC is necessary in situations, in which there are only predictions of data available.

Example 2
The parameter choices of the previous experiment imply that the choice of a feedback κ fulfilling (7) is not possible. In this section we want to present numerical results for a setup that allows such a feedback in principle.
To this end we choose the following parameter values: The diffusion coefficient is set to d = 0.5, the reaction coefficient to r = 0.5, and the advection field is chosen as As initial condition we choose which is a smooth approximation of the discontinuous function It can be shown that y = 0 is a stable steady state of the PDE (1) to the control u = 0. Therefore, we choose the desired temperatures y 1 1 = y 1 2 = 0 ∈ L 2 (0, 1; H) and y 2 1 = y 2 2 = 0 ∈ H together with the parameter values σ 1 1 = σ 1 2 = 1, σ 2 1 = σ 2 2 = 0.1 and σ 3 1 = σ 3 2 = 0.1. The linear functionals in the cost functions are given by Thus, both cost functions measure the deviation of the state from the steady state y = 0 together with some control costs. To make the cost functions conflicting, the linear functionals are chosen such that the first cost function penalizes deviations from the steady state mostly in the left half of the domain, whereas the second cost function penalizes the deviation in the right half of the domain.
If we run Algorithm 3 with this parameter setting, we still observe that in many iterations there is no feedback κ fulfilling (7), especially for the larger MPC horizons of = 25 and = 50. The reason for this is that the initial Pareto front J 0 opt (y • ) is quite close or even dominates the Pareto front J  Since the MOMPC points cannot perform better than any point on the Pareto front J (0,1) opt (y • ), it is not possible to achieve the performance result from Theorem 3 in those parts. As all the other assumptions of Theorem 3 are satisfied, the reason for not achieving the performance result has to be that it is not possible to find a feedback κ fulfilling (7) in all steps of the algorithm. This is what we observe numerically.
However, if we look at Figure 5 we see that we still obtain the performance result from Theorem 3 for many points, although a feedback κ cannot be found for all these points. The black lines indicate which performance valueĴ corresponds to which initial cost. Therefore, if the black lines are pointing to the top right (starting from the red point), this means that the performance result is satisfied for this point. We observe that the performance result holds for all points for an MPC horizon length of = 13 and for most points for = 25. Even for = 50 there are some points in the middle of the Pareto front, for which the performance result holds, although there are only very few steps, where a feedback κ fulfilling (7) can be found. Figure 6 indicates that this performance behavior does not change drastically, if we increase the time horizon t end . The reason for this is that at time t end = 1 the steady state y = 0 is almost reached, so that almost no further costs are needed after this point.  Figure 7, we observe that increasing the MPC horizon has two positive effects on the results: Firstly, the MOMPC points are located closer to the Pareto front, and secondly, they are spread more evenly over the entire Pareto front. Already for an MPC horizon length of = 50 we can see that the Pareto front is almost perfectly approximated.

Conclusions and Outlook
In this paper we have adopted the framework of MOMPC presented in References [10,11] for ODEs to parabolic advection-diffusion-reaction equations. A key ingredient for the MOMPC algorithm is a gradient-based descent scheme for multiobjective optimization taken from Reference [25,26]. On the one hand, this scheme results in a Pareto optimal point, and on the other hand, it intrinsically implies that the crucial condition (5) is satisfied. Assuming the existence of a feedback κ satisfying (7) we could show that the performance result from ( [10], Theorem 11) also holds for our framework of the MOMPC algorithm. Numerically, it turns out that it is hard to guarantee that such a feedback exists. However, even in the case that such a feedback does not exist, the results of the MOMPC algorithm are comparable with the Pareto optimal points of the full-horizon open-loop problem. Especially when some parameters of the PDE can only be predicted, the MOMPC algorithm turns out to be a useful tool, if we still want to compute many compromises between the cost functions.
Both the analysis in Section 3 and the numerical experiments in Section 4 were only carried out for two cost functions, but it is possible to extend both to an arbitrary number of cost function as long as all cost functions are of the form (2). For the analysis this extension is straight-forward, since all the arguments can be transferred directly to this more general case. On the numerical side we can still use the Euclidean reference point method for computing the initial Pareto set U 0 opt (y • ), see for example, Reference [29]. In general, the effort for computing the Pareto set U 0 opt (y • ) scales exponentially with the number of cost functions. However, this computation has only to be performed once in step 4 of Algorithm 3, which is in many cases not underlying any time restrictions. While performing the rest of the MOMPC algorithm the number of cost functions only influences the complexity of steps 6 and 9. In step 6 we need to perform Algorithm 2, in which the dimension of the variable space of the minimization problem (13) is equal to the number of cost functions. Since we only have to compute one Pareto optimal controlū n ∈ U n opt (y a ) this step is computationally not expensive. In step 9 the number of cost functions in the min-max problem 8 corresponds to the total number of cost functions. Again the effort for solving this problem does not cause complex calculations.
During the MOMPC Algorithm 3 numerous open-loop problems on the MPC horizon have to be solved. Since this involves the solution of the PDE and its associated adjoint equation, this is numerically expensive. Therefore, in a future work one can apply model-order reduction techniques (e.g., the proper orthogonal decomposition method, cf. References [30,31]) for lowering the computational effort. Strategies presented in References [16,17] could be used to show convergence to the full solution.
Author Contributions: All authors have contributes equally. All authors have read and agreed to the published version of the manuscript.
Funding: S.B. and G.F. gratefully acknowledge partially support by the German DFG-Priority Program 1962 and German Science Fund DFG grant "Reduced-Order Methods for Nonlinear Model Predictive Control", respectively. Moreover, S.B. is partially funded by the German Landesgraduiertenförderung of Baden-Württemberg.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

MO
Multiobjective optimization MOP Multiobjective optimization problem MOCP Multiobjective optimal control problem MPC Model predictive control MOMPC Multiobjective model predictive control ODE Ordinary differential equation PDE Partial differential equation