Asymptotic stability of POD based model predictive control for a semilinear parabolic PDE

Asymptotic Stability of POD based Model Predictive Control for a semilinear parabolic PDE

Alessandro Alla^∗

Dipartimento di Matematica, Sapienza Universit`a di Roma, P.le Aldo Moro 2, 00185 Rome, Italy

Stefan Volkwein^∗∗

Department of Mathematics and Statistics, University of Konstanz, 78457 Konstanz, Germany

Abstract

In this article a stabilizing feedback control is computed for a semilinear parabolic partial differential equation utilizing a nonlinear model predictive (NMPC) method. In each level of the NMPC algorithm the finite time horizon open loop problem is solved by a reduced-order strategy based on proper orthogonal decomposition (POD). A stability analysis is derived for the combined POD- NMPC algorithm so that the lengths of the finite time horizons are chosen in order to ensure the asymptotic stability of the computed feedback controls. The proposed method is successfully tested by numerical examples.

Keywords: Dynamic programming, nonlinear model predictive control, asymptotic stability, suboptimal control, proper orthogonal decomposition.

2000 MSC:35K58, 49L20, 65K10, 90C30.

1. Introduction

In many control problems it is desired to design a stabilizing feedback con- trol, but often the closed-loop solution can not be found analytically, even the unconstrained case since it involves the solution of the corresponding Hamilton- Jacobi-Bellman equations. One approach to circumvent this problem is the repeated solution of open-loop optimal control problems. The first part of the resulting open-loop input signal is implemented and the whole process is re- peated. Control approaches using this strategy are referred to as model pre-

∗This author wishes to acknowledge the support obtained by the ESF Grant no 4160.

∗∗This author gratefully acknowledges support by the DFG grant VO no 1658/2-1. S.

Volkwein is the corresponding author.

Email addresses: alla@mat.uniroma1.it(Alessandro Alla), stefan.volkwein@uni-konstanz.de(Stefan Volkwein)

Preprint submitted to Advances in Computational Mathematics December 4, 2013 Konstanzer Online-Publikations-System (KOPS)

dictive control(MPC), moving horizon controlor receding horizon control. In general one distinguishes between linear andnonlinearMPC (NMPC). In linear MPC, linear models are used to predict the system dynamics and considers linear constraints on the states and inputs. Note that even if the system is linear, the closed loop dynamics are nonlinear due to the presence of constraints. NMPC refers to MPC schemes that are based on nonlinear models and/or consider a nonquadratic cost functional and general nonlinear constraints. Although linear MPC has become an increasingly popular control technique used in industry, in many applications linear models are not sufficient to describe the process dynamics adequately and nonlinear models must be applied. This inadequacy of linear models is one of the motivations for the increasing interest in nonlinear MPC; see. e.g., [2, 10, 13, 21]. The prediction horizon plays a crucial role in MPC algorithms. For instance, the quasi infinite horizon NMPC allows an effi- cient formulation of NMPC while guaranteeing stability and the performances of the closed-loop as shown in [3, 11] under appropriate assumptions.

Since the computational complexity of NMPC schemes grows rapidly with the length of the optimization horizon, estimates for minimal stabilizing hori- zons are of particular interest to ensure stability while being computationally fast. Stability and suboptimality analysis for NMPC schemes without stabiliz- ing constraints are studied in [13, Chapter 6], where the authors give sufficient conditions ensuring asymptotic stability with minimal finite prediction horizon.

Note that the stabilization of the problem and the computation of the mini- mal horizon involve the (relaxed)dynamic programming principle (DPP); see [14, 20]. This approach allows estimates of the finite prediction horizon based on controllability properties of the dynamical system.

Since several optimization problems have to be solved in the NMPC method, it is reasonable to apply reduced-order methods to accelerate the NMPC al- gorithm. Here, we utilize proper orthogonal decomposition (POD) to derive reduced-order models for the nonlinear dynamical systems; see, e.g., [16, 25]

and [15]. The application of POD is justified by an a priori error analysis for the considered nonlinear dynamical system, where we combine techniques fro [17, 18] and [24]. Let us refer to [12], where the authors also combine successfully an NMPC scheme with a POD reduced-order approach. However, no analysis is carried out ensuring the asymptotic stability of the proposed NMPC-POD scheme. Our contribution focusses on the stability analysis of the POD-NMPC algorithm without terminal constraints, where the dynamical system is a semi- linear parabolic partial differential equation with an advection term. A minimal finite horizon is determined to guarantee stabilization of the system. Our ap- proach is motivated by the work [4]. The main difference here is that we have added an advection term in the dynamical system and utilize a POD suboptimal strategy to solve the open-loop problems. Since the minimal prediction horizon can be large, the numerical solution of the open-loop problems is very expen- sive within the NMPC algorithm. The application of the POD model reduction reduces efficiently the computational cost by computing suboptimal solutions.

But we involve this suboptimality in our stability analysis in order to ensure the asymptotic stability of our NMPC scheme.

The paper is organized in the following manner: In Section 2 we formulate our infinite horizon optimal control problem governed by a semilinear parabolic equation and bilateral control constraints. The NMPC algorithm is introduced in Section 3. For the readers convenience, we recall the known results of the stability analysis. Further, the stability theory is applied to our underlying nonlinear semilinear equations and bilateral control constraints. In Section 4 we investigate the finite horizon open loop problem which has to be solved at each level of the NMPC algorithm. Moreover, we introduce the POD reduced- order approach and prove an a-priori error estimate for the semilinear parabolic equation. Finally, numerical examples are presented in Section 5.

2. Formulation of the control system

Let Ω = (0,1)⊂R be the spatial domain. For the initial time t_◦ ∈ R⁺0 = {s ∈ R|s ≥0} we define the space-time cylinder Q = Ω×(t_◦,∞). By H = L²(Ω) we denote the Lebesgue space of (equivalence classes of) functions which are (Lebesgue) measurable and square integrable. We endowH by the standard inner product – denoted byh·,·iH – and the associated induced normkϕkH = hϕ, ϕi^1/2_H . Furthermore,V =H₀¹(Ω)⊂H stands for the Sobolev space

V =

ϕ∈H

Z

Ω

ϕ⁰(x)

2dx <∞andϕ(0) =ϕ(1) = 0

.

Recall that bothH andV are Hilbert spaces. InV we use the inner product hϕ, φi_V =

Z

Ω

ϕ⁰(x)φ⁰(x) dx forϕ, φ∈V

and setkϕkV =hϕ, ϕi^1/2_V forϕ∈V. For more details on Lebesgue and Sobolev spaces we refer the reader to [9], for instance. When the time t is fixed for a given function ϕ: Q → R, the expression ϕ(t) stands for a function ϕ(·, t) considered as a function in Ω only. Recall that the Hilbert spaceL²(Q) can be identified with the Bochner spaceL²(t◦,∞;H).

We consider the following control system governed by the following semi- linear parabolic partial differential equation: y =y(x, t) solves the semilinear initial boundary value problem

y_t−θy_xx+y_x+ρ(y³−y) =u in Q, (2.1a) y(0,·) =y(1,·) = 0 in (t◦,∞), (2.1b)

y(t_◦) =y_◦ in Ω. (2.1c)

In (2.1a) it is assumed that the controlu=u(x, t) belongs to the set of admis- sible control inputs

Uad(t_◦) =

u∈U(t_◦)

u(x, t)∈Uadfor almost all (f.a.a.) (x, t)∈Q , (2.2)

where U(t_◦) = L²(t_◦,∞;H) and Uad = {u ∈ R|ua ≤ u ≤ ub} with given ua ≤0≤ub . The parametersθ andρsatisfy

(θ, ρ)∈Dad=

(˜θ,ρ)˜ ∈R²

θa≤θ˜andρa≤ρ˜

with positive θa and ρa. Further, in (2.1c) the initial condition y_◦ =y_◦(x) is supposed to belong toH.

A solution to (2.1) is interpreted in the weak sense as follows: for given (t_◦, y_◦)∈ R⁺0 ×H and u∈Uad(t_◦) we call y a weak solution to (2.1) for fixed (θ, ρ)∈Dadify(t)∈V,yt(t)∈V⁰ hold f.a.a. t≥t◦ and y satisfiesy(t◦) =y◦

inH as well as d

dthy(t), ϕi_H+ Z

Ω

θyx(t)ϕ⁰+ yx(t) +ρ(y³(t)−y(t)) ϕdx=

Z

Ω

u(t)ϕdx (2.3) for allϕ∈V and f.a.a. t > t_◦. The following result is proved in [6], for instance.

Proposition 2.1. For given (t_◦, y_◦)∈R⁺0 ×H andu∈Uad(t_◦) there exists a unique weak solutiony=y_{[u,t◦,y◦]} to (2.1)for every (θ, ρ)∈Dad.

Let (t◦, y◦) ∈ R⁺0 ×H be given. Due to Proposition 2.1 we define the quadratic cost functional:

J(u;ˆ t_◦, y_◦) :=1 2

Z ∞ t◦

ky[u,t◦,y◦](t)−y_dk²_Hdt+λ 2

Z ∞ t◦

ku(t)k²_Hdt (2.4) for allu∈U(t_◦)⊃Uad(t_◦), wherey_{[u,t◦,y◦]} denotes the unique weak solution to (2.1). We suppose thaty_d=y_d(x) is a given desired stationary state inH (e.g., the equilibrium y_d = 0) and that λ > 0 denotes a fixed weighting parameter.

Then we consider the nonlinear infinite horizon optimal control problem min ˆJ(u;t_◦, y_◦) subject to (s.t.) u∈Uad(t_◦). (2.5) Suppose that the trajectoryyis measured at discrete time instances

tn=t◦+n∆t, n∈N,

where the time step ∆t >0 stands for the time step between two measurements.

Thus, we want to select a controlu∈Uad(t) such that the associated trajectory y[u,t◦,y◦] follows a given desired stateyd as good as possible. This problem is called atracking problem, and, ify_d= 0 holds, astabilization problem.

Since our goal is to be able to react to the current deviation of the stateyat timet=t_n from the given reference valuey_d, we would like to have the control in feedback form, i.e., we want to determine a mapping µ:H →Uad(t_◦) with u(t) =µ(y(t)) for t∈[t_n, t_n+1].

3. Nonlinear model predictive control

We present an NMPC approach to compute a mapping µ which allows a representation of the control in feedback form. For more details we refer the reader to the monographs [13, 21], for instance.

3.1. The NMPC method

To introduce the NMPC algorithm we write the weak form of our control system (2.1) as a parametrized nonlinear dynamical system. Let us introduce theθ-dependent linear operator Awhich maps the spaceV into its dual space V⁰ as follows:

Aϕ=−θϕxx+ϕx∈V⁰ forϕ∈V andθ≥θa. Moreover, letf be a mapping fromV intoV⁰ given by

f(ϕ) =ρ(ϕ³−ϕ)∈V⁰ forϕ∈V andρ≥ρa.

Setting F(ϕ, v) = Aϕ+f(ϕ)−v forϕ ∈ V, v ∈ H and (θ, ρ) ∈D_ad we can express (2.3) as the nonlinear dynamical system

y⁰(t) =F(y(t), u(t))∈V⁰ for allt > t_◦, y(t_◦) =y_◦in H (3.1) for given (t_◦, y_◦)∈R⁺0 ×H. The cost functional has been already introduced in (2.4). Summarizing, we want to solve the following infinite horizon minimization problem

min ˆJ(u;t_◦, y_◦) = Z ∞

t◦

` y_{[u,t◦,y◦]}(t), u(t)

dt s.t. u∈Uad(t_◦), (P(t_◦)) where we have defined the running quadratic cost as

`(ϕ, v) = 1 2

kϕ−ydk²_H+λkvk²_H

forϕ, v∈H. (3.2) If we have determined a state feedbackµfor (P(t_◦)), the controlu(t) =µ(y(t)) allows a closed loop representation for t ∈ [t_◦,∞). Then, for a given initial conditiony₀ ∈ H we sett_◦ = 0, y_◦ = y₀ in (3.1) and insert µ to obtain the closed-loop form

y⁰(t) =F(y(t), µ(y(t))) inV⁰ fort∈(t_◦,∞),

y(t_◦) =y_◦ inH. (3.3)

Although an infinite horizon problem may be very hard to solve due to the dimensionality of the problem, it guarantees the stabilization of the problem.

This is a very important issue for optimal control problems. In an NMPC algorithm a state feedback law is computed for (P(t_◦)) by solving a sequence of finite time horizon problems. Let us mention that another important tool to compute a feedback law is given by the solution of the Hamilton-Jacobi-Bellman equation; see, e.g., [5, 9] and [19].

To formulate the NMPC algorithm we introduce the finite horizon quadratic cost functional as follows: for (t◦, y◦)∈R⁺0 ×H andu∈U^Nad(t◦) we set

Jˆ^N(u;t_◦, y_◦) = Z t^N_◦

t◦

` y_{[u,t◦,y◦]}(t), u(t) dt,

whereN is a natural number,t^N_◦ =t_◦+N∆tis the final time andN∆tdenotes the length of the time horizon for the chosen time step ∆t > 0. Further, we introduce the Hilbert space U^N(t◦) = L²(t◦, t^N_◦;H) and the set of admissible controls

U^Nad(t_◦) =

u∈U^N(t_◦)

u(x, t)∈U_ad f.a.a. (x, t)∈Q^N

with Q^N = Ω×(t_◦, t^N_◦) ⊂ Q; compare (2.2). In Algorithm 1 the method is presented. In each iteration over n we store the optimal control on the first Algorithm 1(NMPC algorithm)

Require: time step ∆t >0, finite horizonN ∈N, weighting parameterλ >0.

1: forn= 0,1,2, . . .do

2: Measure the statey(tn)∈V of the system attn=n∆t.

3: Sett_◦=tn =n∆t,y_◦=y(tn) and compute a global solution to

min ˆJ^N(u;t◦, y◦) s.t. u∈U^Nad(t◦). (P^N(t◦)) We denote the obtained optimal control by ¯u^N.

4: Define the NMPC feedback value µ^N(y(t)) = ¯u^N(t), t ∈ (t_◦, t_◦ + ∆t]

and use this control to compute the associated statey =y_[µN(·),t◦,y_◦] by solving (3.1) on [t_◦, t_◦+ ∆t].

5: end for

time interval [tn, tn+1] and the associated optimal trajectory of the sampling time. Then, we initialize a new finite horizon optimal control problem whose initial condition is given by the optimal trajectory ¯y(t) = y_[µN(·),t◦,y◦](t) at t=t_◦+ ∆tusing the optimal controlµ^N(y(t)) = ¯u(t) fort∈(t_◦, t_◦+ ∆t] . We iterate this process. Of course, the larger horizon the better approximation one can have, but we would like to have the minimal horizon which can guarantee stability [3]. Notice that (P^N(t_◦)) is an open loop problem on a finite time horizon [t_◦, t_◦+N∆t] which will be studied in Section 4.

3.2. Dynamic programming principle (DPP) and asymptotic stability

For the readers convenience we now recall the essential theoretical results from DPP and stability analysis. Let us first introduce the so called value functionv defined as follows for an infinite horizon optimal control problem:

v(t_◦, y_◦) := inf

u∈Uad(t◦)

Jˆ(u;t_◦, y_◦) for (t_◦, y_◦)∈R⁺0 ×H.

LetN ∈N be chosen. Due to the DPP the value function v satisfies for any k∈ {1, . . . , N}witht^k_◦=t_k+k∆t:

v(t_◦, y_◦)

= inf

u∈U^kad(t_◦)

(Z t^k_◦

t◦

` y_{[u,t◦,y◦]}(t), u(t)

dt+v t_◦+k∆t, y_{[u,t◦,y◦]}(t_◦+k∆t) )

which holds under very general conditions on the data; see, e.g., [5] for more details. The value function for the finite horizon problem (P^N(t◦)) is of the following form:

v^N(t_◦, y_◦) = inf

u∈U^Nad(t◦)

Jˆ^N(u;t_◦, y_◦) for (t_◦, y_◦)∈R⁺0 ×H.

The value functionv^N satisfies the DPP for the finite horizon problem fort◦+ k∆t, 0< k < N:

v^N(t_◦, y_◦)

= inf

u∈U^k_ad(t◦)

(Z t◦+k∆t t_◦

` y_{[u,t◦,y◦]}(t), u(t)

dt+v^N y_{[u,t◦,y◦]}(t_◦+k∆t) )

.

Nonlinear stability properties can be expressed by comparison functions which we recall here for the readers convenience [13, Definition 2.13].

Definition 3.1. We define the following classes of comparison functions:

K=

β:R⁺0 →R⁺0

β is continuous, strictly increasing andβ(0) = 0 , K_∞=

β:R⁺0 →R⁺0

β ∈ K, β is unbounded , L=n

β :R⁺0 →R⁺0

β is continuous, strictly decreasing, lim

t→∞β(t) = 0o , KL=

β:R⁺0 ×R⁺0 →R⁺0

β is continuous,β(·, t)∈ K, β(r,·)∈ L . Utilizing a comparison functionβ∈ KLwe introduce the concept of asymp- totic stability; see, e.g. [13, Definition 2.14].

Definition 3.2. Let y_{[µ(·),t◦,y◦]} be the solution to (3.3)andy_∗∈H an equilib- rium for (3.3), i.e., we have F(y_∗, µ(y_∗)) = 0. Then, y_∗ is said to be locally asymptotically stable if there exist a constant η > 0 and a function β ∈ KL such that the estimate

ky_[µ(·),t◦,y_◦](t)−y∗k_H ≤β ky◦−y∗k_H, t) holds for ally_◦∈H satisfyingky◦−y_∗kH< η and all t≥t_◦.

Let us recall the main result about asymptotic stability via DPP; see [14].

Proposition 3.3. Let N ∈Nbe chosen and the feedback mappingµ^N be com- puted by Algorithm1. Assume that there exists anα^N ∈(0,1]such that for all (t_◦, y_◦)∈R⁺0 ×H therelaxed DPP

v^N(t◦, y◦)≥v^N t◦+ ∆t, y_[µN(·),t◦,y◦](t◦+ ∆t)

+α^N` y◦, µ^N(y◦))

(3.4) holds. Furthermore, we have for all(t_◦, y_◦)∈R⁺0 ×H:

α^Nv(t_◦, y_◦)≤α^NJˆ(µ^N(y_[µN(·),t◦,y◦]);t_◦, y_◦)≤v^N(t_◦, y_◦)≤v(t_◦, y_◦), (3.5)

where y_[µN(·),t◦,y◦] solves the closed-loop dynamics (3.3) with µ = µ^N. If, in addition, there exists an equilibriumy∗∈H andα1, α2∈ K∞ satisfying

`∗(y◦) = min

u∈Uad

`(y◦, u)≥α1 ky◦−y∗k_H

, (3.6a)

α2 ky◦−y∗k_H

≥v^N(t◦, y◦) (3.6b)

hold for all(t_◦, y_◦)∈R⁺0 ×H, theny_∗ is a globallyasymptotically stable equi- libriumfor (3.3)with the feedback mapµ=µ^N and value functionv^N. Remark 3.4. 1) Our running cost`defined in (3.2) satisfies condition (3.6a)

for the choice yd = y_∗. Further, (3.6b) follows from the finite horizon quadratic cost functional ˆJ^N, the definition of the value functionv^N and our a-priori analysis presented in Lemma 3.6 below. Therefore, we only have to check the relaxed DPP (3.4).

2) It is proved in [14] that lim

N→∞α^N = 1. Hence, we would like to findα^N close to one to have the best approximation ofv in terms ofv^N. On the other hand, a large N implies that the numerical solution of (P^N(t◦)) is

much more involved. ♦

In order to estimate α^N in the relaxed DPP we require the exponential controllability property for the system.

Definition 3.5. System (3.1) is called exponentially controllable with respect to the running cost`if for each(t_◦, y_◦)∈R⁺0 ×H there exist two real constants C >0,σ∈[0,1) and an admissible control u∈Uad(t_◦)such that:

`(y<sub>[u,t◦,y◦]</sub>(t), u(t))≤Cσ<sup>t−t◦</sup>`_∗(y_◦) f.a.a. t≥t_◦. (3.7) We have the next a-priori estimate for the uncontrolled solution to (3.1), i.e., the solution foru= 0.

Lemma 3.6. Let (t_◦, y_◦)∈R⁺0 ×H and u= 0∈Uad(t_◦). Then, the solution y=y_{[0,t◦,y◦]} to (3.1)satisfies the a-priori estimate

ky(t)k_H≤e^{−γ(t−t◦)}ky◦k_H f.a.a. t≥t_◦ (3.8) withγ=γ(θ, ρ) =θ/C_V −ρ.

Proof. Recall thatV is continuously (even compactly) embedded intoH. Due to the Poincar´e inequality [9] there exists a constantC_V >0 such that

kϕk_H≤CVkϕk_V for allϕ∈V. (3.9) Using (3.9), choosingu(t) = 0 andϕ=y(t) in (2.3) we obtain

d

dtky(t)k²_H+ 2θ CV

ky(t)k²_H≤2ρky(t)k²_H f.a.a. t≥t◦

which implies d

dtky(t)k²_H ≤2

ρ− θ CV

ky(t)k²_H f.a.a. t≥t_◦. Thus, by Gronwall’s inequality we derive (3.8) withγ=θ/CV −ρ.

Remark 3.7. Forθ > ρCV we haveγ >0. Then, (3.8) implies thatky(t)kH <

ky◦kH for any t > t_◦. Moreover, it is easy to check that the origin y_◦ = 0 is

unstable forγ <0. ♦

Let us choose y_d = 0. Suppose that we have a particular class of state feedback controls of the formu(x, t) =−Ky(x, t) with a positive constantK;

see [4]. This assumption helps us to derive the exponential controllability in terms of the running cost ` and to compute a minimal finite time prediction horizonN∆tensuring asymptotic stability. In this case, (3.8) has to be modified because we do not setu= 0, but u=−Ky. Utilizing similar arguments as in the proof of Lemma 3.6 we find for a givenK >0 that the statey=y_{[−Ky,t◦,y◦]} satisfies

ky(t)k_H≤e^{−γ(K)(t−t◦)}ky_◦k_H f.a.a. t≥t_◦ (3.10) withγ(K) =θ/CV +K−ρ. Thus, if K > ρ−θ/CV holds, ky(t)kH tends to zero fort→ ∞. Combining (3.10) with the desired exponential controllability (3.7) and usingy_d= 0 we obtain for allt≥t_◦ (see [4]):

`(y(t), u(t)) = 1

2 ky(t)k²_H+λku(t)k²_H

= 1

2(1 +λK²)ky(t)k²_H

≤ 1

2C(K)e^{−2γ(K)(t−t◦)}ky◦k²_H =C(K)σ(K)^t−t◦`_∗(y_◦)

(3.11)

f.a.a. t≥t_◦ and for every (t_◦, y_◦)∈R⁺0 ×H, where

C(K) = (1 +λK²), σ(K) =e^−2γ(K), γ(K) =θ/C_V +K−ρ. (3.12) In the following theorem we provide an explicit formula for the scalarα^N in (3.4). A complete discussion is given in [14].

Theorem 3.8. Assume that the system (3.1) and` statisfy the controllability condition (3.7). Let the finite prediction horizonN∆tbe given withN ∈Nand

∆t >0. Then the parameterα^N depends onK and is given by:

α^N(K) = 1− ηN(K)−1 QN

i=2 ηi(K)−1 QN

i=2ηi(K)−QN

i=2 ηi(K)−1 (3.13) whereηi(K) =C(1−σⁱ)/(1−σ)and the constants C=C(K),σ=σ(K)are given by Definition3.5.

K y_◦a <0 y_◦a>0 y_◦b<0 K < ub/|y_◦b| no constraints y◦b>0 K <min

ua/y◦b, ub/|y◦a| K <|ua|/y◦b

Table 3.1: Constraints for the feedback factor K in u(x, t) = −Ky(x, t) considering the bilateral control constraints (3.14) and the initial condition (3.15).

Remark 3.9. Theorem 3.8 suggests how we can compute the minimal horizon N ensuring asympotic stability. Due to (3.12) we maximize

1− ηN(K)−1 QN

i=2 ηi(K)−1 QN

i=2ηi(K)−QN

i=2 ηi(K)−1, ηi(K) = (1 +λK²)1−e^{−2i(θ/CV+K−ρ)} 1−e^{−2(θ/CV+K−ρ)} with respect to K > max(0, ρ−θ/C_V) and N ∈ N in order to get α^N > 0.

Further, we suppose thatu∈ U^N(t_◦) holds. Hence, we have to guarantee the bilateral control constraints

ua≤ −Ky(x, t)≤ub f.a.a. (x, t)∈Q^N (3.14) withua ≤0≤ub. Under these assumptions, the computation of K andN has to take into account the influence of the control constraints. Since we determine K in such a way that γ(K) =θ/CV +K−ρ >0 is satisfied, we derive from (3.10) that

ky(t)k_H≤ ky◦k_H f.a.a. t≥t◦.

Let us suppose that we have y◦ 6= 0 andky(t)k_C(Ω) ≤ ky◦k_C(Ω) f.a.a. t ≥t◦. Then, we define

y_◦a= min

x∈Ω

y_◦(x), y_◦b= max

x∈Ω

y_◦(x). (3.15)

Then, K has to satisfy K > max(0, ρ−θ/CV) and the restrictions shown in Table 3.1. Summarizing,K has always an upper bound due to the constraints u_a, u_b and a lower bound due to the stabilization related toγ(K)>0. ♦ 4. The finite horizon problem (P^N(t_◦))

In this section we discuss (P^N(t◦)), which has to be solved at each level of Algorithm 1.

4.1. The open loop problem

Recall that we have introduced the final timet^N_◦ =t_◦+N∆tand the control spaceU^N(t_◦) =L²(t_◦, t^N_◦;H). The spaceY^N(t_◦) =W(t_◦, t^N_◦) is given by

W(t◦, t^N_◦) =

ϕ∈L²(t◦, t^N_◦;V)

ϕt∈L²(t◦, t^N_◦;V⁰) ,

which is a Hilbert space endowed with the common inner product [8, pp. 472- 479]. We define the Hilbert spaceX^N(t◦) =Y^N(t◦)×U^N(t◦) endowed with the standard product topology. Moreover, we introduce the Hilbert spaceZ^N(t◦) = Z^N1(t◦)×HwithZ^N1(t◦) =L²(t◦, t^N_◦;V) and the nonlinear operatore= (e1, e2) : X^N(t◦)→Z^N(t◦)⁰ by

he1(x), ϕi_ZN

1(t◦)⁰,Z^N₁(t◦)= Z t^N_◦

t_◦

hyt(t), ϕ(t)i_V0,Vdt +

Z t^N_◦ t◦

Z

Ω

θyx(t)ϕ(x) +

yx(t) +ρ y(t)³−y(t)

−u(t)

ϕ(t) dxdt, he2(x), φi_H=hy(t◦)−y◦, φi_H

forx= (y, u)∈X^N(t_◦), (ϕ, φ)∈Z^N(t_◦), where we identify the dualZ^N(t_◦)⁰ of Z^N(t_◦) with L²(t_◦, t^N_◦;V⁰)×H andh·,·i_ZN

1(t_◦)⁰,Z^N1(t_◦) denotes the dual pairing betweenZ^N1(t◦)⁰ andZ^N1(t◦). Then, for givenu∈U^N(t◦) the weak formulation for (2.3) can be expressed as the operator equation e(x) = 0 inZ^N(t◦)⁰. Fur- ther, we can write (P^N(t◦)) as a constrained infinite dimensional minimization problem

minJ(x) = Z t^N_◦

t◦

`(y(t), u(t)) dt s.t. x∈F^Nad(t_◦) (4.1) with the feasible set

F^Nad(t_◦) =

x= (y, u)∈X^N(t_◦)

e(x) = 0 inZ^N(t_◦)⁰ andu∈U^Nad(t_◦) . For given fixed controlu∈U^Nad(t◦) we consider the state equatione(y, u) = 0∈ Z^N(t◦)⁰, i.e.,y satisfies

d

dthy(t), ϕi_H+ Z

Ω

θyx(t)ϕ⁰+ yx(t) +ρ(y(t)³−y(t)) ϕdx

= Z

Ω

u(t)ϕdx f.f.a. t∈(t_◦, t^N_◦ ], hy(t_◦), ϕi_H =hy_◦, ϕi_H

(4.2)

for allϕ∈V. The following result is proved in [26, Theorem 5.5].

Proposition 4.1. For given (t_◦, y_◦)∈R⁺0 ×H andu∈U^Nad(t_◦) there exists a unique weak solutiony∈Y^N(t_◦)to (4.2)for every(θ, ρ)∈Dad. If, in addition, y◦ is essentially bounded in Ω, i.e., y◦ ∈L^∞(Ω) holds, we have y ∈ L^∞(Q^N) satisfying

kykY^N(t◦)+kyk_L∞(Q^N)≤C kuk

U^N(t◦)+ky_◦k_L∞(Ω)

(4.3)

for aC >0, which is independent ofuandy_◦.

Utilizing (4.3) it can be shown that (4.1) possesses at least one (local) optimal solution which we denote by ¯x^N = (¯y^N,u¯^N)∈F^Nad(t_◦); see [26, Chapter 5]. For

the numerical computation of ¯x^N we turn to first-order necessary optimality conditions for (4.1). To ensure the existence of a unique Lagrange multiplier we investigate the surjectivity of the linearization e⁰(¯x^N) :X^N(t◦)→ Z^N(t◦)⁰ of the operator e at a given point ¯x^N = (¯y^N,u¯^N) ∈X^N(t◦). Notice that the Fr´echet derivativee⁰(¯x^N) = (e⁰₁(¯x^N), e⁰₂(¯x^N)) ofeat ¯x^N is given by

he⁰₁(¯x^N)x, ϕi

Z^N1(t◦)⁰,Z^N1(t◦)= Z t^N_◦

t◦

hy_t(t), ϕ(t)i_V0,V dt +

Z t^N_◦ t_◦

Z

Ω

θyx(t)ϕ(x) +

yx(t) +ρ 3¯y^N(t)²−1

y(t)−u(t)

ϕ(t) dxdt, he⁰₂(¯x^N)x, φi_H =hy(t_◦), φi_H

forx= (y, u)∈X^N(t_◦), (ϕ, φ)∈Z^N(t_◦). Now, the operatore⁰(¯x^N) is surjective if and only if for an arbitraryF = (F1, F2) ∈Z^N(t◦)⁰ there exists a pair x= (y, u)∈X^N(t◦) satisfying e⁰(¯x^N) =F in Z^N(t◦)⁰ which is equivalent with the fact that there exists an u ∈ U^N(t◦) and an y ∈ Y^N(t◦) solving the linear parabolic problem

yt−θyxx+yx+ρ(3¯y²−1)y=F1 inZ^N1 (t◦)⁰, y(t◦) =F2in H. (4.4) Utilizing standard arguments [8] it follows that there exists for anyu∈U^N(t_◦) a unique y ∈ Y^N(t_◦) solving (4.4). Thus, e⁰(¯x^N) is a surjective operator and the local solution ¯x^N to (4.1) can be characterized by first-order optimality conditions. We introduce the Lagrangian by

L(x, p, p_◦) =J(x) +he(x),(p, p_◦)i

Z^N(t_◦)⁰,Z^N(t_◦)

for x∈ X^N(t_◦) and (p, p_◦) ∈ Z^N(t_◦). Then, there exists a unique associated Lagrange multiplier pair (¯p^N,p¯_◦) to (4.1) satisfying the optimality system

∇yL(¯x^N,p¯^N,p¯^N_◦)y= 0 ∀y∈Y^N(t_◦) (adjoint equation)

∇_uL(¯x^N,p¯^N,p¯^N_◦)(u−¯u^N)≥0 ∀u∈U^Nad(t_◦) (variational inequality), he(¯x^N),(p, p_◦)i

Z^N(t_◦)⁰,Z^N(t_◦)= 0 ∀(¯p,p¯0)∈Z^N(t_◦)(state equation).

It follows from variational arguments that the strong formulation for the adjoint equation is of the form

−¯p^N_t −θp¯^N_xx−p¯^N_x −ρ 1−3(¯y^N)²

¯

p^N =yd−y¯^N in Q^N,

¯

p^N(0,·) = ¯p^N(1,·) = 0 in (t◦, t^N_◦),

¯

p^N(t^N_◦) = 0 in Ω.

(4.5)

Moreover, we have ¯p^N_◦ = ¯p^N(t_◦). The variational inequality base the form Z t^N_◦

t◦

Z

Ω

(λ¯u^N−p¯^N)(u−u¯^N) dxdt≥0 for allu∈U^Nad(t_◦). (4.6) Using the techniques as in [27, Proposition 2.12] one can proof that second- order sufficient optimality conditions can be ensured provided the residuum k¯y^N −ydk_L2(t◦,t^N_◦;H) is sufficiently small.

4.2. POD reduced order model for open-loop problem

To solve (4.1) we apply a reduced-order discretization based on proper or- thogonal decomposition (POD); see [15]. In this subsection we briefly introduce the POD method, present an a-priori error estimate for the POD solution to the state equation e(x) = 0 ∈Z^N(t_◦)⁰ and formulate the POD Galerkin approach for (4.1).

4.2.1. The POD method for dynamical systems

By X we denote either the function space H orV. Then, for ℘∈Nlet the so-calledsnapshots ortrajectoriesy^k(t)∈X are given f.a.a. t∈[t◦, t^N_◦] and for 1≤k≤℘. At least one of the trajectoriesy^k is assumed to be nonzero. Then we introduce the linear subspace

V= spann

y^k(t)|t∈[t_◦, t^N_◦] a.e. and 1≤k≤℘o

⊂X (4.7)

with dimension d ≥ 1. We call the set V snapshot subspace. The method of POD consists in choosing a complete orthonormal basis inX such that for every l≤dthe mean square error betweeny^k(t) and their corresponding l-th partial Fourier sum is minimized on average:









 min

℘

X

k=1

Z t^N_◦ t◦

y^k(t)−

l

X

i=1

hy^k(t), ψ_ii_Xψ_i

2 X

dt s.t. {ψi}^l_i=1⊂X andhψi, ψji_X=δij, 1≤i, j≤l,

(P^l)

where the symbolδ_ijdenotes the Kronecker symbol satisfyingδ_ii= 1 andδ_ij = 0 fori6=j. An optimal solution{ψ¯_i}^l_i=1 to (P^l) is called a POD basis of rankl.

The solution to (P^l) is given by the next theorem. For its proof we refer the reader to [15, Theorem 2.13].

Theorem 4.2. Let X be a separable real Hilbert space andy^k₁, . . . , y_n^k ∈X are given snapshots for 1 ≤ k ≤ ℘. Define the linear operator R : X → X as follows:

Rψ=

℘

X

k=1

Z t^N_◦ t◦

hψ, y^k(t)i_Xy^k(t) dt forψ∈X. (4.8) Then, R is a compact, nonnegative and symmetric operator. Suppose that {λ¯i}_i∈N and {ψ¯i}_i∈N denote the nonnegative eigenvalues and associated or- thonormal eigenfunctions ofRsatisfying

Rψ¯i= ¯λiψ¯i, ¯λ1≥. . .≥¯λd>λ¯d+1=. . .= 0, ¯λi→0asi→ ∞. (4.9) Then, for everyl≤d the first ` eigenfunctions {ψ¯_i}^l_i=1 solve (P^l). Moreover, the value of the cost evaluated at the optimal solution{ψ¯_i}^l_i=1 satisfies

E(l) =

℘

X

k=1

Z t^N_◦ t◦

y^k(t)−

l

X

i=1

hy^k(t),ψ¯ii_Xψ¯i

2 X

dt=

d

X

i=l+1

¯λi. (4.10)

In real computations, we do not have the whole trajectories y^k(t) at hand f.a.a. t∈[t◦, t^N_◦] and for 1≤k≤℘. Therefore, we suppose that we are given a time grid 0 =t1< . . . < tn =t^N_◦ forn∈N. To ease the presentation we suppose that the time grids for all snapshots are the same. This can be generalized in a straightforward way. Let y^k_j denote an approximation for y^k(tj) ∈ X for 1≤j ≤nand 1≤k≤℘. We assume that at least one of they^k_j’s is nonzero.

Let us introduce the linear space Vⁿ = span y_j^k

1 ≤j ≤nand 1 ≤k ≤℘ with dimension dⁿ = dimVⁿ ∈ {1, . . . , n℘}. Analogous to (P^l) the discrete variant of the POD method consists in choosing a complete orthonormal basis in X such that for every l ≤dⁿ the mean square error between the y_j^k’s and their correspondingl-th partial Fourier sum is minimized on average:









 min

℘

X

k=1 n_k

X

j=1

αⁿ_j y_j^k−

l

X

i=1

hy^k_j, ψ_ii

Xψ_i

2 X

s.t. {ψi}^l_i=1⊂X andhψi, ψji_X=δij, 1≤i, j≤l,

(P^l,n)

where theαⁿ_j’s stand for the trapezoidal weights αⁿ₁ =t₂−t₁

2 , αⁿ_j = t_j+1−t_j−1

2 for 1≤j≤n, αⁿ_n= t_n−t_n−1

2 .

The solution to (P^l,n) is given by the solution to the eigenvalue problem Rⁿψ_iⁿ=

℘

X

k=1 n_k

X

j=1

α^k_jhy^k_j, ψ_iⁿi

Xy^k_j =λⁿ_iψ_iⁿ, 1≤i≤l,

where Rⁿ : X → Vⁿ ⊂ X is a linear, compact, selfadjoint and nonnegative operator; see, e.g., [15, Theorem 2.7]. Thus, there exists an orthonormal set {ψ¯ⁿ_i}_i∈N of eigenfunctions and corresponding nonnegative eigenvalues {λ¯ⁿ_i}_i∈N such that

Rⁿψ¯_iⁿ= ¯λⁿ_iψ¯ⁿ_i, ¯λⁿ₁ ≥¯λⁿ₂ ≥. . .≥¯λⁿ_dn>¯λⁿ_dn+1=. . .= 0. (4.11) We refer to [15, Section 2.3], where the relationship between (4.9) and (4.11) is investigated. Further, in [15, Remark 2.1] the equivalence of (4.11) with the singular value decomposition is discussed forX =R^m,℘= 1 andα_jⁿ= 1.

4.2.2. The Galerkin POD scheme for the state equation

Suppose that (t_◦, y_◦)∈R⁺0 ×H andt^N_◦ =t_◦+N∆twith prediction horizon N∆t >0. For given fixed control u∈U^Nad(t_◦) we consider the state equation e(y, u) = 0∈Z^N(t_◦)⁰, i.e.,y satisfies (4.2). Let us turn to a POD discretization of (4.2). To keep the notation simple we apply only a spatial discretization with POD basis functions, but no time integration by, e.g., the implicit Euler method.

Therefore, we utilize the continuous version of the POD method introduced in Section 4.2.1. In this section we distinguish two choices for X: X = H and X =V. We choose the snapshots y¹ =y and y² =yt, i.e., we set ℘= 2. By

Proposition 4.1 the snapshotsy^k,k= 1, . . . , ℘, belong toL²(0, T;V). According to (4.9) let us introduce the following notations:

RVψ=

℘

X

k=1

Z t^N_◦ t◦

hψ, y^k(t)i_V y^k(t) dt forψ∈V,

RHψ=

℘

X

k=1

Z t^N_◦ t◦

hψ, y^k(t)i_Hy^k(t) dt forψ∈H.

To distinguish the two choices for the Hilbert spaceXwe denote by the sequence {(λ^V_i , ψ_i^V)}i∈N ⊂R⁺0 ×V the eigenvalue value decomposition for X =V, i.e., we have

RVψ_i^V =λ^V_i ψ^V_i for alli∈N. Furthermore, let{(λ^H_i , ψ^H_i )}_i∈N⊂R⁺0 ×H in satisfy

R_Hψ_i^H=λ^H_i ψ^H_i for alli∈N.

Then, d = dimRV(V) = dimRH(H) ≤ ∞; see [24]. The next result – also taken from [24] – ensures that the POD basis{ψ^H_i }^l_i=1 of ranklbuild a subset of the test spaceV.

Lemma 4.3. Suppose that the snapshotsy^k ∈L²(0, T;V),k= 1, . . . , ℘. Then, we haveψ^H_i ∈V fori= 1, . . . , d.

Let us define the two POD subspaces V^l= span

ψ₁^V, . . . , ψ_l^V ⊂V, H^l= span

ψ^H₁ , . . . , ψ^H_l ⊂V ⊂H, whereH^l⊂V follows from Lemma 4.3. Moreover, we introduce the orthogonal projection operatorsP_H^l :V →H^l⊂V andP^l:V →V^l⊂V as follows:

v^l=P_H^l ϕfor anyϕ∈V iffv^l solves min

w^l∈H^lkϕ−w^lk_V, v^l=P_V^lϕfor anyϕ∈V iffv^l solves min

w^l∈V^lkϕ−w^lk_V. (4.12) It follows from the first-order optimality conditions for (4.12) that v^l = P_H^l ϕ satisfies

hv^l, ψ_i^Hi_V =hϕ, ψ_i^Hi_V, 1≤i≤l. (4.13) Writing v^l ∈ H^l in the formv^l = Pl

j=1v_j^lψ^H_j we derive from (4.13) that the vector v^l= (v^l₁, . . . ,v^l_l)^>∈R^l satisfies the linear system

l

X

j=1

hψ^H_j , ψ^H_i i_Vv^l_j =hϕ, ψ_i^Hi_V, 1≤i≤l.

For the operatorP_V^l we have the explicit representation P_V^lϕ=

l

X

i=1

hϕ, ψ^V_i i_Vψ_i^V forϕ∈V.

Moreover, we introduce the orthogonal projection operatorP^l:V →V^`by P_V^lϕ=

l

X

i=1

hϕ, ψiiV ψi forϕ∈V. (4.14) Further, we conclude from (4.10) that

℘

X

k=1

Z T 0

ky^k(t)− P_V^ly^k(t)k²_V dt=

d

X

i=l+1

λ^V_i . (4.15) Next we review an essential result from [24, Theorem 6.2], which is essential in our a-priori error analysis for the choice X = H. Recall thatH^l ⊂ V holds.

Consequently,kψ_i^H− P_H^l ψ_i^HkV is well-defined for 1≤i≤l.

Theorem 4.4. Suppose thaty^k ∈L²(0, T;V)for1≤k≤℘. Then,

℘

X

k=1

Z T 0

ky^k(t)− P_H^l y^k(t)k²_V dt=

d

X

i=l+1

λ^H_i kψ_i^H− P_H^l ψ_i^Hk²_V.

Moreover, P_H^l y^k converges to y^k in L²(0, T;V) as l tends to ∞ for each k ∈ {1, . . . , ℘}.

Let us define the linear spaceX^l⊂V as X^l= span

ψ1, . . . , ψl ,

where ψi = ψ_i^V in case of X = V and ψi = ψ^H_i in case of X = H. Hence, X^l = V^l and X^l = H^l for X = V and X = H, respectively. Now, a POD Galerkin scheme for (4.2) is given as follows: findy^l(t)∈X^l f.a.a. t∈[t_◦, t^N_◦] satisfying

d

dthy^l(t), ψi_H+ Z

Ω

θy^l_x(t)ψ⁰+ y_x^l(t) +ρ(y^l(t)³−y^l(t)) ψdx

= Z

Ω

u(t)ψdx f.f.a. t∈(t_◦, t^N_◦], hy^l(t_◦), ψi_H =hy_◦, ψi_H

(4.16)

for allψ∈X^l. It follows by similar arguments as in the proof of Proposition 4.1 that there exists a unique solution to (4.16). Ify_◦∈L^∞(Q^N) holds,y^lsatisfies the a-priori estimate

ky^lk

Y^N(t◦)+ky^lk_L∞(Q^N)≤C ky_◦k_L∞(Ω)+kuk

U^N(t◦)

. (4.17) where the constant C > 0 is independent of l and y_◦. Let P^l denote P_V^l in case of X = V and P_H^l in case of X = H. To derive an error estimate for ky−y^lk_YN(t◦)we make use of the decomposition

y(t)−y^l(t) =y(t)− P^ly(t) +P^ly(t)−y^l(t) =%^l(t) +ϑ^l(t) f.a.a. t∈[t_◦, t^N_◦ ]