A Comparative Stability Analysis of Neumann and Dirichlet Boundary MPC for the Heat Equation

Abstract We continue our study on stability properties of model predictive control without terminal constraints applied to the heat equation. Here, we focus on boundary control and in particular on the differences between Neumann and Dirichlet boundary conditions. We first illustrate the differences by means of numerical simulations and then explain our findings by a theoretical analysis.


INTRODUCTION
In this paper we continue our study of stability properties for model predictive control (MPC) applied to the one dimensional heat equation. We particularly investigate the case of Neumann boundary control and compare this case to our earlier findings for Dirichlet boundary control and distributed control.
MPC, also called receding horizon control (RHC), is a well known method in which a feedback controller is synthesized from the iterative solution of finite horizon optimal control problems. The application of MPC to parabolic PDEs for boundary control can be found, e.g., in Dubljevic and Christofides (2006) and Ito and Kunisch (2002). In Dubljevic and Christofides (2006), the authors use a modal decomposition technique in order to obtain low dimensional problems and focus on state and control constraints. In Ito and Kunisch (2002), a receding horizon method for infinite dimensional systems is presented. These references have in common that the stability of the closed loop solution is guaranteed by adding terminal constraints and control Lyapunov functionals as terminal costs to the finite horizon optimal control problems. However, as the construction of terminal constraints and costs is a challenging task, here we consider MPC without such constraints and costs, which is also often preferred in practical applications, cf. Qin and Badgwell (2003).
One important parameter in unconstrained MPC schemes is the length of the optimization horizon N . It is well known that under suitable conditions stability holds for sufficiently large N , see Grüne and Pannek (2011, Chapter 6 and the references therein). In this paper we investigate the dependence of the minimal stabilizing N on the system structure and parameters. Since the computational effort in each MPC step grows with N , this information is of particular interest. We investigate this problem for the linear heat equation, extending Altmüller and Grüne (2012) by including the case of Neumann boundary control and comparing this case with our earlier results on Dirichlet boundary control and distributed control.
Our paper is organized as follows. In Section 2 we recall stability results for MPC schemes without terminal constraints. The PDE control system under consideration is introduced in Section 3 in which we also present numerical simulation results motivating our study. The parameters needed for the stability analysis are derived in Section 4 in which we first summarize useful earlier results from Altmüller and Grüne (2012) and then present our new results on Neumann boundary control. Finally, in Section 5 we use these results in order to compare Neumann with Dirichlet boundary and distributed control and analytically explain the numerical findings. Section 6 concludes the paper.

PRELIMINARIES
In this section we summarize results on MPC and its stability. We consider nonlinear discrete time control systems given by z(n + 1) = f (z(n), u(n)) (1) with state space Z and set of control values U , which are both possibly infinite dimensional normed spaces. Continuous time PDE models can be converted into this discrete form by sampling, for details see Section 3. The solution trajectory for given initial state z 0 and control sequence (u(n)) n∈N0 ⊆ U is denoted by z u (·) = z u (·; z 0 ) and the space of control sequences u : N 0 → U by U. Our goal is to compute an optimal asymptotically stabilizing feedback law for (1) by MPC. To this end we define the finite horizon cost functional with optimization horizon N ∈ N ≥2 and the optimal value function Given an initial value z 0 , the MPC feedback µ N : Z → U and the corresponding solution z µ N of the closed loop system z µ N (n + 1; z 0 ) = f (z µ N (n; z 0 ), µ(z µ N (n; z 0 ))) (4) are then obtained in the following iterative way: 1. Fix some N ∈ N, set z µ N (0; z 0 ) := z 0 and n := 0 2. Minimize (2) with initial value z µ N (n; z 0 ) in order to obtain the optimal control sequence u ⋆ (0), u ⋆ (1), . . . , u ⋆ (N − 1) and set µ N (z µ N (n; z 0 )) := u ⋆ (0) 3. Compute z µ N (n + 1; z 0 ) by (4), set n := n + 1 and go to 2.
In order to achieve stability of the closed loop system additional terminal constraints or terminal costs are often introduced in (2), cf. Rawlings and Mayne (2009) or Grüne and Pannek (2011, Chapter 5) and the references therein. In contrast to this approach, here we use MPC without additional constraints or costs. We remark that state and control constraints coming from modelling considerations can be included by choosing appropriate subsets Z ⊆ Z and U ⊆ U but since this topic is not in the focus of this paper we will not go into details concerning feasibility issues.
Next we state a stability result from Grüne et al. (2010) which requires the following assumption. Assumption 1. The system (1) is called exponentially controllable with respect to the stage costs ℓ if there exist an overshoot bound C ≥ 1 and a decay rate σ ∈ (0, 1) such that for each z ∈ Z there exists u z ∈ U satisfying ℓ(z uz (n; z), u z (n)) ≤ Cσ n min u∈U ℓ(z, u) =: Cσ n ℓ ⋆ (z) (5) for all n ∈ N 0 .
Based on this assumption, the following stability theorem is proven in Grüne et al. (2010). Theorem 2. (Stability Theorem). Let the Assumption 1 hold with overshoot constant C ≥ 1 and decay rate σ ∈ (0, 1). Let the prediction horizon N be such that the stability condition holds with γ i := C i−1 n=0 σ n . Assume, moreover, that K ∞functions η, η exist satisfying η( z − z ⋆ ) ≤ ℓ ⋆ (z) and η( z − z ⋆ ) ≥ V N (z) for all z ∈ Z. Then the closed loop system (4) is asymptotically stable and thus in particular each trajectory converges to the equilibrium z ⋆ as n → ∞.
Our goal is now to use this theorem in order to determine the minimal horizon N which is needed in order to ensure (6). While quantitative estimates of horizons N for which stability holds derived from Theorem 2 tend to be conservative (see, e.g., Altmüller et al. (2010)), it was observed in Altmüller and Grüne (2012) that Assumption 1 and Condition (6) allow for a very precise explanation of the qualitative behavior of such minimal horizons N in dependence on system parameters and system structure. This leads to the following road map for our analysis: • Find (not necessarily optimal) controls u z such that the exponential controllability condition is fulfilled. • Calculate C and σ for this particular control.
• Use this information to explain the qualitative behavior of the minimal horizon N guaranteeing (6).
In order to give an intuition for Condition (6), in Figure 1 we visualize the minimal N satisfying (6) depending on C and σ: the areas in Figure 1 visualize the pairs (C, σ) corresponding to the minimal stabilizing horizon N = . . . indicated within the respective area. The key information from this visualization for our subsequent analysis is the different impact of the constants: For fixed σ it is always possible to obtain stability with the shortest possible horizon N = 2 by reducing C. In contrast to this, for fixed C it is in general impossible to arbitrarily reduce N by reducing σ. Therefore, the overshoot constant C plays an important role in the stability analysis.

HEAT EQUATION
In this section we introduce the equation which we want to stabilize using MPC. We consider a one dimensional linear diffusion-reaction equation, henceforth briefly denoted as heat equation, with Dirichlet and Neumann boundary control. We assume that the control can only act on the right boundary, i.e., we consider (9) On the left side we impose a homogenous Dirichlet condition, the domain is given by Ω = (0, 1) and y 0 ∈ H 1 0 (Ω) denotes the initial condition. The uncontrolled equation (v(t) ≡ 0) is known to be unstable if the reaction parameter µ becomes larger than the least eigenvalue of (−∂ xx ), cf. Krstic and Smyshlyaev (2008). Thus, the equation is unstable for µ ≥ π 2 in the Dirichlet controlled case and for µ ≥ π 2 /4 for the Neumann controlled PDE. Figure  results of the optimal control problem for both Dirichlet and Neumann boundary control can be found in Lasiecka and Triggiani (2000).
Although we focus on boundary control in this paper, for our analysis it will turn out to be useful to also consider the heat equation with distributed control, i.e., the equation in which the control v(x, t) acts on the whole domain Ω. Imposing homogeneous Dirichlet boundary conditions this equation reads The uncontrolled PDE is again unstable for µ ≥ π 2 . The existence and regularity theory of the distributed optimal control problem can be found in Lions (1971).
In order to use the theory presented in Section 2 we have to rewrite the continuous time PDE as a discrete time system (1). Let ϕ(x, t; v, y 0 ) be the solution of (7, 8), (7,9) or (10) with control v and initial function y 0 . Then we define the sampled data system with sampling period T > 0 as z(n + 1) := ϕ(x, T ; u(n), z(n)) (11) with z(0) = y 0 . In the boundary control case we define the discrete time control as u(n) In the sequel, this relation between v and u will be implicitly used when verifying (5). The discrete time n corresponds to the continuous time nT which implies z(n) = y(·, nT ) ∈ H 1 (Ω) =: Z.
It is known that the choice of the running cost plays an important role for stability of the MPC closed loop. Results of the influence of different stage cost on the horizon can be found in Altmüller and Grüne (2012). Since in this paper we do not focus on the stage cost but rather on the different types of boundary control, we exclusively use the stage cost i.e., we penalize the state in the L 2 (Ω)-norm and the squared absolute value of the control where λ > 0 denotes a so called regularization or Tikhonov parameter. The analogous choice for the distributed problem reads (13) In order to motivate our subsequent study we first numerically evaluate the stability properties for Dirichlet and Neumann boundary control. The following simulations illustrate certain phenomena which we will later explain in our theoretical analysis. To this end, using the sampling time T = 0.01 and the initial function y 0 (x) = 1 5 sin(πx), by performing numerical simulations of the MPC closed loop we have determined the minimal horizon N N for which asymptotic stability could be observed. In these simulations, the optimization problem within each MPC step was solved by a Newton-CG method, cf. Borzì and Schulz (2012). The arising PDEs were discretized by a finite difference scheme. The results are shown in Table 1. Note that with our numbering N N = 2 is the smallest possible horizon. For comparison, Table 2 shows the corresponding values N D for Dirichlet boundary control and N dis for distributed control, see also Altmüller and Grüne (2012). Note that in both tables we have only considered those values of µ for which the respective uncontrolled equation is unstable. 4  2  2  2  5  3  2  2  7  4  2  2  9  6  2  2  11  7  3  3  13  8  5  4  15  9  6  5   Table 1. Minimal stabilizing horizon N N for Neumann boundary control with varying λ and µ determined from numerical simulations of the MPC closed loop  One observes that in all cases a smaller value of the regularization parameter λ leads to a smaller horizon. However, for Neumann boundary control only for values of µ below µ ≈ 10 it is possible to obtain stability for the shortest possible horizon N N = 2 by reducing λ, while for larger values of µ the decrease stops at a value N N ≥ 3 even for very small values of λ.
Comparing this behavior with that of the Dirichlet boundary control and the distributed control from Table 2, one sees that for values of µ below µ ≈ 10 the dependence of N N on λ is qualitatively similar to distributed control while for larger values of µ it is similar to Dirichlet boundary control. It is the purpose of the remainder of the paper to theoretically explain this numerically observed behavior.

DERIVATION OF C AND σ
Following our road map from Section 2 we will now derive estimates for the constants C and σ in Assumption 1. We start by recalling results from Altmüller and Grüne (2012) for the distributed control case, as these will turn out to be useful for the analysis of the boundary control case we discuss afterwards.

Distributed Control
In order to perform the first step of the road map presented in Section 2 we have to find a control u z for any initial value z such that the exponential controllability condition (5) is satisfied. A key issue in this method is that we do not need the optimality of this control. It is known that the linear feedback law v(x, t) := −Ky(x, t) (14) stabilizes (10) for a sufficient large constant K ∈ R. This property is used in the following theorem in order to establish (5). Theorem 3. The heat equation (10) with control (14), K > µ − π 2 , and stage costs (13) fulfills the exponential controllability condition (5). The corresponding constants are given by σ = e −2T (π 2 −µ+K) ∈ (0, 1) and C = 1 + λK 2 ∈ R.

Boundary Control
In this subsection we derive the exponential controllability constants for the boundary controlled PDE. Like in the distributed case we intend to find an exponentially stabilizing feedback law for (7). One approach to obtain such a feedback is the method of backstepping, cf. Krstic and Smyshlyaev (2008) for details. The idea behind this technique is to transform (7) into a so called "target system" which is exponentially stable. The transformation is done with a Volterra kernel A suitable target system is given by w t (x, t) = w xx (x, t) + (µ − K)w(x, t) (17) w(x, 0) = w 0 , w(0, t) = 0 and w(1, t) = 0 in the case of Dirichlet control or w x (1, t) = 0 for the Neumann control. One advantage of this backstepping approach is that the resulting Volterra kernel is given as an explicit formula. For our system (7) and the target system (17) the transformation and the inverse transformation are given by where I 1 (x) denotes the first modified Bessel function and J 1 (x) the first classical Bessel function. The feedback control law coming from the Volterra transformation is given by for the Dirichlet control and by (21) for the Neumann controlled system. A derivation that this transformation converts the system (7) into (17) can be found in Krstic and Smyshlyaev (2008). For a proof concerning the kernel properties see Liu (2003). In the subsequent Theorem 5 we use this control in order to establish (5). For its proof we need the following lemma whose proof is straightforward and thus omitted. Lemma 4. Let u(x) ∈ L 2 (0, 1) and for i = 1, 2 define u( holds with the constant L i : Dirichlet Boundary Control Now we specialize our results to the case of Dirichlet boundary control, for which we can prove the following theorem, cf. also Altmüller and Grüne (2012). Theorem 5. The Dirichlet boundary controlled heat equation (7) with control (20), K > µ − π 2 , and stage costs (12) satisfies the exponential controllability condition (5). The corresponding constants are given by σ = e −2T (π 2 −µ+K) ∈ (0, 1) and C = (1 + λK 2 η(K))ξ(K) ∈ R with η(K) := Proof. Again we explain the main steps of the proof of Altmüller and Grüne (2012, Theorem 5.2). First we note that the target system is exactly the reduced system (15) we used in the proof of Theorem 3. Thus, we can use the information about the decay rate. By combining Theorem 3 and Lemma 4 we obtain y(·, t) 2 L 2 (Ω) ≤ (1 + L 2 ) 2 w(·, t) 2 and thus with ξ(K) := (1 + L 1 (K)) 2 (1 + L 2 (K)) 2 and σ := e −2T (π 2 −µ+K) . Furthermore, we get ℓ(z(n), u(n)) = 1 2 y(·, nT ) 2 withC := (1 + λK 2 η(K)) and Combining (23) and (24) then yields the assertion.
Neumann Boundary Control Finally we now consider the heat equation with Neumann boundary control. In contrast to before, here we perform a different analysis for small and large values of the reaction parameter µ. First we will derive the exponential constants for sufficiently small values of µ. Afterwards, we consider the case of arbitrary values for µ.
In order to construct a control we use again the method of backstepping. The stage cost is given by (12).
It is known that the feedback y x (1, t) = −Ky(1, t) (25) stabilizes the problem for small values of the reaction parameter µ. In order to derive the corresponding constants C and σ, we need the following generalization of the Friedrichs inequality. Theorem 6. For u ∈ H 1 (Ω) the following inequality holds where C D is the least eigenvalue of the Dirichlet eigenproblem −∆u = µu in Ω, u = 0 in ∂Ω.
Note that Theorem 8 is valid for larger range of reaction parameters µ than Theorem 7, however, at the expense that the derived constants attain larger values.

EXPLANATION OF NUMERICAL OBSERVATIONS
We can now use Theorem 7 and 8 to explain the observations from Table 1. From the spectral properties we know that the uncontrolled equation is unstable for µ ≥ π 2 /4. Since C D = π 2 holds in Theorem 7, we can apply this theorem for π 2 /4 ≤ µ < π 2 . Consequently, for these values of µ we obtain the overshoot constant C = 1 + λK 2 M which we can bring arbitrarily close to one by reducing the value of the regularization parameter λ. Therefore, Figure 1 reveals that for µ < π 2 we can always achieve stability with the smallest possible horizon N = 2. This explains our numerical findings from Table 1 for µ < π 2 . The similarity of the derived overshoot constant with that from the distributed control case derived in Theorem 3, i.e., C = 1 + λK 2 , moreover explains the qualitatively similarity of the two cases for µ < π 2 .
For µ ≥ π 2 we have to resort to Theorem 8 which yields the larger overshoot constant (1 + λK 2 ( √ M /2 + η 2 (K)) 2 )ξ(K). This constant is similar to that derived in Theorem 5 for Dirichlet boundary control (C = (1 + λK 2 η(K))ξ(K)), which explains why the dependence of N N on λ in Table 1 is qualitatively similar to N D in Table 2. Particularly, we obtain C → ξ(K) for λ → 0. Since ξ(K) is monotonically increasing with ξ(K) > 1 for K ∈ R + we cannot obtain C → 1 as λ → 0 and therefore we cannot expect stability for N = 2 even for small values of λ, cf. Figure 1.
An interesting question is whether the bound µ = π 2 which limits the applicability of the "better" Theorem 7 in our theoretical stability analysis is actually a tight bound. In order to get evidence for this fact, we have performed further numerical simulations of the MPC closed loop. Since our analysis depends on both C and σ but we only investigated the behavior of C, we need to keep the influence of σ small in our numerical simulations. One way to achieve this is be reducing the sampling time T , since for T → 0 we obtain σ → 1 and thus we can observe the pure influence of C. Figure 3 displays the largest value of µ for which stability is obtained with N = 2 depending on the sampling time T (solid blue line) while the theoretical bound π 2 is indicated by the dashed line. It is clearly visible that the curve tends to the theoretical bound for small values of T . For even smaller values of T , the simulations did not yield reasonable information because the numerical errors became predominant. Still, the simulations strongly indicate that the bound µ = π 2 separating the two different stability behaviors is indeed tight.

CONCLUSIONS
We have shown that the control structure and the stage cost significantly influence the length of the minimal stabilizing horizon in MPC without stabilizing terminal constraints. Our controllability based analysis is able to explain several -at the first glance puzzling -numerically observed phenomena for the heat equation. While it is not clear whether a similar analysis is feasible for more complex PDE models, we nevertheless think that our analysis provides insight for the choice of stage costs guaranteeing good performance of MPC with short optimization horizons, which can also be valuable in practical applications.