POD-Based Bicriterial Optimal Control by the Reference Point Method

2nd IFAC Workshop on Control of Systems Governed 2nd IFAC Control 2ndPartial IFAC Workshop Workshop on Control of of Systems Systems Governed Governed by Differentialon Equations 2nd IFAC Workshop on Control of Systems Governed by Partial Differential Equations Available online at www.sciencedirect.com by Partial Differential Equations June 13-15, 2016. Bertinoro, Italy by Partial Differential Equations June 13-15, 2016. Bertinoro, Italy June June 13-15, 13-15, 2016. 2016. Bertinoro, Bertinoro, Italy Italy

ScienceDirect

IFAC-PapersOnLine 49-8 (2016) 210–215

POD-Based Bicriterial Optimal Control POD-Based POD-Based Bicriterial Bicriterial Optimal Optimal Control Control the Reference Point Method the Reference Point Method the Reference Point Method

by by by

Stefan Banholzer ∗∗ Dennis Beermann ∗∗ Stefan Volkwein ∗∗ ∗∗ Stefan Banholzer ∗∗ Dennis Beermann ∗∗ Stefan Volkwein ∗∗ Stefan Stefan Banholzer Banholzer Dennis Dennis Beermann Beermann Stefan Stefan Volkwein Volkwein ∗∗ ∗ ∗ Department of Mathematics and Statistics, University of Konstanz, ∗ of Mathematics and Statistics, University of Konstanz, ∗ Department Department and University 78457 Konstanz, Germany Department of of Mathematics Mathematics and Statistics, Statistics, University of of Konstanz, Konstanz, 78457 Konstanz, Germany ∗∗ 78457 Konstanz, Germany Department of Mathematics and Statistics, University of Konstanz, 78457 Konstanz, Germany ∗∗ ∗∗ Department of Mathematics and Statistics, University of Konstanz, ∗∗ Department of Germany Mathematics and Statistics, Statistics, University of of Konstanz, Konstanz, 78457 Konstanz, (e-mail: [email protected]) Department of Mathematics and University 78457 Konstanz, Germany (e-mail: [email protected]) 78457 78457 Konstanz, Konstanz, Germany Germany (e-mail: (e-mail: [email protected]) [email protected]) Abstract: In the present paper a bicriterial optimal control problem governed by a parabolic Abstract: In the present paper a bicriterial optimal control problem governed by a parabolic Abstract: In paper a optimal control problem by partial differential equation (PDE) and bilateral control is considered. For the Abstract: In the the present present paper a bicriterial bicriterial optimal controlconstraints problem governed governed by aa parabolic parabolic partial differential equation (PDE) and bilateral control constraints is considered. For the partial differential equation (PDE) and bilateral control constraints is considered. For numerical optimization the reference point method is utilized. The PDE is discretized by partial differential equation (PDE) and bilateral control constraints is considered. For the the numerical optimization the reference point method is utilized. The PDE is discretized by numerical optimization theutilizing reference point method is utilized. utilized. The PDE PDE is discretized discretized by anumerical Galerkin optimization approximation thepoint method of proper orthogonal decomposition (POD). the reference method is The is by a Galerkin approximation utilizing the method of proper orthogonal decomposition (POD). a Galerkin Galerkin approximation utilizing the reduced-order method of of proper proper orthogonal for decomposition (POD). POD is a powerful approach to derive approximations evolution problems. a approximation utilizing the method orthogonal decomposition (POD). POD is a powerful approach to derive reduced-order approximations for evolution problems. POD is reduced-order approximations Numerical examples approach illustrate to thederive efficiency of the proposed strategy. for POD is a a powerful powerful approach to derive reduced-order approximations for evolution evolution problems. problems. Numerical examples illustrate the efficiency of the proposed strategy. Numerical Numerical examples examples illustrate illustrate the the efficiency efficiency of of the the proposed proposed strategy. strategy. © 2016, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Keywords: Optimal control, multiobjective optimization, reference point method, proper Keywords: Optimal control, multiobjective optimization, optimization, reference reference point point method, proper Keywords: Optimal control,a-posteriori multiobjective proper orthogonal decomposition, erroroptimization, analysis. Keywords: Optimal control, multiobjective reference point method, method, proper orthogonal decomposition, a-posteriori error analysis. orthogonal decomposition, a-posteriori error analysis. orthogonal decomposition, a-posteriori error analysis. 1. INTRODUCTION problems. Numerical examples are presented in Section 5. 1. INTRODUCTION INTRODUCTION problems. Numerical examples are presented in Section 5. 1. problems. Numerical examples are presented Finally, a conclusion drawn in 6. in 1. INTRODUCTION problems. Numerical is examples areSection presented in Section Section 5. 5. Finally, a conclusion is drawn in Section 6. aa conclusion is drawn in Section 6.2 Finally, conclusion is drawn in Section In real applications, optimization problems are often de- Finally, 1 6. n Throughout this paper, if x1 , x2 ∈ Rn are two In real real applications, applications, optimization optimization problems problems are are often often dede- Notation: In 2 ∈ Rn Throughout this paper, if x11 ,, x are two scribed by introducing several objective functions conIn real applications, optimization problemsfunctions are oftenconde- Notation: Notation: Throughout if ∈ R are two vectors, we write x11 ≤ this x22 ifpaper, x1i1 ≤ x n, and scribed by introducing several objective Notation: Throughout this paper, if2i22 x xfor ,x xii2 = ∈ 1, Rn..., are two scribed by introducing objective functions con1 2 1 vectors, we write x ≤ x if x ≤ x for = 1, ..., n, and flicting with each other.several This leads to multiobjective or similarily scribed by introducing several objective functions con1 2 i ≤ x2 i for i = 1, ..., n, and 1 2 1 vectors, we write x ≤ x if x for x < x . flicting with with each each other. other. This This leads leads to to multiobjective multiobjective or or similarily vectors, weforwrite 1 x 2≤ x if xii ≤ xii for i = 1, ..., n, and flicting multicriterial optimization problems; Ehrgott (2005), Miflicting with each other. This leads Ehrgott to multiobjective or similarily similarily for for x x11 < 0 the state equation is given by maximal comfort; Fong et al. (2006) and Kusiak et al. maximal comfort; et (2006) Kusiak al. (2011). Finding theFong optimal control thatand represents good time T > 00 the state equation is given by maximal comfort; Fong et al. al. (2006) and Kusiakaa et et al. For For time T > the is m (2011). Finding the optimal control that represents good For time T > 0 the state state equation equation is given given by by  (2011). Finding themain optimal control thatproblems. representsFor aa good m compromise is the issue in these that (2011). Finding the optimal control that represents good  m m  (t, x) − ∆y(t, x) = u (t)χ (x) for (t, x) ∈ Q, y compromise is is the the main main issue issue in in these these problems. problems. For For that that t i i  compromise reason the concept of Pareto efficient points is yyt (t, compromise is the main issue optimal in theseor For that (t, x) x) − − ∆y(t, ∆y(t, x) x) = = i=1 u uiii (t)χ (t)χiii (x) (x) for for (t, (t, x) x) ∈ ∈ Q, Q, reason the the concept concept of Pareto Pareto optimal orproblems. efficient points points is x) − ∆y(t, x) = u (t)χ (x) for (t, x) ∈ Q, ytt (t, reason of optimal or efficient is developed. In contrast to scalar-valued optimization probi=1 reason the concept of Pareto optimal or efficient points is (1) i=1 developed. In contrast to scalar-valued optimization prob∂y i=1 developed. In scalar-valued optimization problems, the computation a set of Pareto optimal points ∂y developed. In contrast contrast to toof optimization prob- (1) (1) (t, x) = 0 for (t, x) ∈ Σ, (1) ∂y lems, the the computation computation ofscalar-valued set of of Pareto Pareto optimal points points ∂y (t, (t, x) x) = = 00 for (t, (t, x) ∈ ∈ Σ, ∂n lems, of aa set optimal is required. Consequently, many scalar-valued constrained for lems, the computation of a set of Pareto optimal points (t, x) = 0 for (t, x) x) ∈ Σ, Σ, ∂n is required. Consequently, many scalar-valued constrained ∂n y(0, x) = y (x) for x ∈ Ω, is required. Consequently, many scalar-valued constrained ∂n ◦ optimization problems have to bescalar-valued solved. is required. Consequently, many constrained y(0, x) = y (x) for x ∈ Ω, ◦ optimization problems have to be solved. y(0, x) = yy◦◦ (x) (x) for x x∈ ∈ Ω, Ω, optimization problems have to be = 3}, for optimization be solved. solved. Ω ⊂ Rdd ,y(0, d ∈x){2, is a bounded domain with In this paperproblems we applyhave the to reference point method Ro- where d where Ω ⊂ R ,, d ∈ {2, 3}, is aa bounded domain with d where Ω ⊂ R d ∈ {2, 3}, is bounded domain with In this paper we apply the reference point method RoLipschitz-continuous boundary Γ = ∂Ω and n stands for where Ω ⊂ R , d ∈ {2, 3}, is a bounded domain with In this paper we apply the reference point method Romaus etpaper al. (2009) in the orderreference to transform a bicriterial In this we apply point method RoLipschitz-continuous boundary Γ = ∂Ω and n stands for Lipschitz-continuous boundary Γ = ∂Ω and n stands for maus et al. (2009) in order to transform a bicriterial the outward normal vector. We set Q = (0, T )×Ω and Σ = Lipschitz-continuous boundary Γ = ∂Ω and n stands for maus et al. (2009) in order to transform a bicriterial optimal control problem into ato sequence of scalar-valued maus et al. (2009) in order transform a bicriterial the outward normal vector. We set Q = (0, T )×Ω and Σ = 2 1 the outward normal vector. We set Q = (0, T )×Ω and Σ = optimal control problem into a sequence of scalar-valued (0, T ) × Γ. Let H = L (Ω) and V = H (Ω) be endowed the outward normal vector. We set Q = (0, T )×Ω and Σ = 2 1 optimal control problem into a sequence of scalar-valued problems and to solve them using welloptimal control problem into a sequence of scalar-valued 2 (Ω) and V = H 1 1 (Ω) be endowed (0, T ) × Γ. Let H = L 2 (0, T ) × Γ. Let H = L (Ω) and V = H (Ω) be endowed optimal control problems and to solve them using wellby the canonical inner products; see Dautray and Lions (0, T ) × Γ. Let H = L (Ω) and V = H (Ω) be endowed optimal control problems and to solve them using well- by the canonical inner products; see Dautray and Lions known optimal control techniques; see Tr¨ oltzsch (2010). optimal control problems and to solve them using wellthe inner products; known optimal optimal control control techniques; techniques; see see Tr¨ Tr¨ ltzsch (2010). (2010). by (2000). The variable u = (u1 , . . . ,see um )Dautray ∈ U = and L22 (0,Lions T )m the canonical canonical inner products; Dautray and Lions known oo Preliminary results combining reduced-order modeling and by known optimal control techniques; see Tr¨ oltzsch ltzsch (2010). 2 m (2000). The variable u = (u .. ,,see u ∈ U = L (0, T ))m 1 ,, .. .. ∞ m) 2 (2000). The variable u = (u u ) ∈ U = L (0, T Preliminary results combining reduced-order modeling and 1 m denotes the control and χ ∈ L (Ω), 1 = 1, . . . , m, are (2000). The variable u = (u , . . . , u ) ∈ U = L (0, T )m i ∞ Preliminary results combining reduced-order modeling and 1 m multiobjective PDE-constrained optimization are recently Preliminary results combining reduced-order modeling and ∞ denotes the control and χ ∈ L (Ω), 1 = 1, . . . , m, are i ∞ ∞ (Ω), 1 = 1, denotes the χ m, are multiobjective PDE-constrained PDE-constrained optimization optimization are are recently recently given control shape and functions. Furthermore, y◦.. ..∈.. ,, L denotes the control control and χii ∈ ∈ L L (Ω), 1 = 1, m, are ∞ (Ω) multiobjective derived Iapichino et al. (2013),optimization Iapichino et are al. recently (2015) given multiobjective PDE-constrained ∞ (Ω) L control shape functions. Furthermore, yy◦ ∈ ∞ ∈ L (Ω) given control shape functions. Furthermore, derived Iapichino et al. (2013), Iapichino et al. (2015) ◦ denotes a given initial heat distribution. We write y(t) ∈ L (Ω) given control shape functions. Furthermore, y derived Iapichino et al. (2013), Iapichino et al. (2015) ◦ write y(t) and Peitz et al. (2015). derived Iapichino et al. (2013), Iapichino et al. (2015) denotes a given initial heat distribution. We denotes initial heat distribution. y(t) and Peitz Peitz et et al. al. (2015). (2015). when y aais given considered a function in x We onlywrite for fixed given initialas heat distribution. We write y(t) and and Peitz (2015). as follows: In Section 2 we intro- denotes when yy is considered as aa function in x only for fixed when is considered as function in x only for fixed The paperetisal. organized t ∈ [0, T ]. Recall that when y is considered as a function in x only for fixed The paper paper is is organized organized as as follows: follows: In In Section Section 22 we we introintro- t ∈ that   The ∈ [0, [0, T T ]. ]. Recall Recall that 2 duce the bicriterial optimization problem under2consideraThe is organized as follows: In Section we intro- tt ∈ [0, ]. Recall that  WT(0, T) =  ϕ ∈ L2 (0, T ; V )  ϕt ∈ L22 (0, T ; V  ) duce paper the bicriterial bicriterial optimization problem under considera  duce the optimization problem under considera2 2  W (0, T )) = ϕ ∈ L T ;; V ))  ϕ L T ;; V tion. Thebicriterial reference point methodproblem is explained in consideraSection 3. duce the optimization under t ∈ 2 (0, 2 (0, ) W (0, T = ϕ ∈ L (0, T V ϕ ∈ L (0, T V ) t tion. The reference point method is explained in Section 3. W (0, T ) = ϕ ∈ L (0, T ; V ) ϕ ∈ L (0, T ; V ) t common inner prodis a Hilbert space endowed with the tion. The reference point method is explained in Section 3. Here we also derive the a-posteriori error estimator which tion. reference method is explained in Section 3. is aa Hilbert space endowed with the common inner prodHere The we also also derivepoint the a-posteriori a-posteriori error estimator estimator which is endowed with common inner prodsee, e.g.,space Dautray and Lions (2000). A weak solution is a Hilbert Hilbert space endowed with the the common inner prodHere we the error is essential inderive our reduced-order approach. Sectionwhich 4 is uct; Here we also derive the a-posteriori error estimator which uct; see, e.g., Dautray and Lions (2000). A weak solution is essential essential in in our our reduced-order reduced-order approach. approach. Section Section 44 is is uct; uct; see, e.g., Dautray and Lions a (2000). (2000). A weak weak solution y ∈ Y = W (0, T ) to (1) is called state and has to satisfy see, e.g., Dautray and Lions A solution is devoted to recall the POD method for optimal control is essential in our reduced-order approach. Section 4 is y ∈ Y = W (0, T ) to (1) is called a state and has to satisfy devoted to to recall recall the the POD POD method method for for optimal optimal control control for yy ∈ Y W T called aa state functions ϕ ∈is : ∈ all Y= =test W (0, (0, T )) to to (1) (1) isV called state and and has has to to satisfy satisfy devoted devoted to author recall gratefully the POD method financial for optimal control for all test functions ϕ ∈ V ::  for all test functions ϕ ∈ V The second acknowledge support by the m   for all test functions ϕ ∈ V :  d  ∇y(t) · ∇ϕ dx =  The second author gratefully acknowledge financial support by the m  d The author gratefully financial by the  m ui (t)χi , ϕH , y(t), ϕH +  project Hybrides Planungsverfahren zur energieeffizienten W¨ arme The second second author gratefully acknowledge acknowledge financial support support by the d m  d y(t), ϕ project Hybrides Planungsverfahren zur energieeffizienten W¨ a rme+ ∇y(t) ·· ∇ϕ dx = ui (t)χ , ϕH ,, dt (2) H i=1 Ω project Hybrides Planungsverfahren zur energieeffizienten W¨ a rmey(t), ϕ + ∇y(t) ∇ϕ dx = und Stromversorgung von st¨ a dtischen Verteilnetzen funded by the y(t), ϕ project Hybrides Planungsverfahren zur energieeffizienten W¨ armeuii (t)χ (t)χiii ,, ϕ ϕH H + Ω ∇y(t) · ∇ϕ dx = i=1u H, dt (2) H und Stromversorgung von st¨ a dtischen Verteilnetzen funded by the dt (2) i=1 Ω dt und Stromversorgung von st¨ a dtischen Verteilnetzen funded by the (2) German Ministry for Economic Affairs Verteilnetzen and Energy. funded by the y◦ , ϕH . y(0), ϕH = i=1 Ω und Stromversorgung von st¨ adtischen German Ministry for Economic Affairs and Energy. = y , ϕ . y(0), ϕ German y(0), German Ministry Ministry for for Economic Economic Affairs Affairs and and Energy. Energy. = y y◦◦◦ ,, ϕ ϕH y(0), ϕ ϕH H = H .. H H

Copyright © 2016, 2016 International Federation of 212Hosting by Elsevier Ltd. All rights reserved. 2405-8963 © IFAC (International Federation of Automatic Control) Copyright 2016 International Federation of 212 Copyright © 2016 International Federation of 212 Peer review© under responsibility of International Federation of Automatic Automatic Control Copyright © 2016 International Federation of 212Control. Automatic Control 10.1016/j.ifacol.2016.07.443 Automatic Control Automatic Control

212 212 212 212

IFAC CPDE 2016 June 13-15, 2016. Bertinoro, Italy

Stefan Banholzer et al. / IFAC-PapersOnLine 49-8 (2016) 210–215

It is shown in Dautray and Lions (2000) that (2) admits a unique solution y and   (3) yY ≤ C y◦ H + uU for a contant C ≥ 0. We introduce the linear operator S : U → Y, where y = Su is the solution to (2) for given u ∈ U with y◦ = 0. From (3) it follows that S is bounded. Moreover, let yˆ ∈ Y be the solution to (2) for u = 0. Then, the affine linear mapping U  u → y(u) = yˆ + Su ∈ Y is affine linear, and y(u) is the weak solution to (1). 2.2 The multiobjective optimal control problem For given ua , ub ∈ U with ua ≤ ub in U, the set of admissible controls is given as    Uad = u ∈ U  ua (t) ≤ u(t) ≤ ub (t) in [0, T ] . Introducing the bicriterial cost functional   2  y(T ) − y 1 Ω H J : Y × U → R2 , J(y, u) = 2 u2U

the multiobjective optimal control problem (MOCP) reads (P) min J(y, u) subject to (s.t.) (y, u) ∈ F(P) with the feasible set    F(P) = (y, u) ∈ Y × Uad  y solves (2) . Next we define the reduced cost function Jˆ = (Jˆ1 , Jˆ2 ) : ˆ U → R2 by J(u) = J(ˆ y + Su, u) for u ∈ U. Then, (P) can be equivalently formulated as ˆ ˆ (P) min J(u) s.t. u ∈ Uad .

ˆ involves the minimization of a vector-valued Problem (P) objective. This is done by using the concepts of order relation and Pareto optimality; see, e.g., Ehrgott (2005). In R2 we make use of the following order relation: For all z 1 , z 2 ∈ R2 we have    z 1 ≤ z 2 ⇔ z 2 − z 1 ∈ R2+ = z ∈ R2  zi ≥ 0 for i = 1, 2 . Definition 1. The point u ¯ ∈ Uad is called Pareto optimal ˆ if there is no other control u ∈ Uad \ {¯ u} with for (P) ˆ ˆ ˆ ˆ Ji (u) ≤ Ji (¯ u), i = 1, 2, and Jj (u) < Jj (¯ u) for at least one j ∈ {1, 2}. 3. THE REFERENCE POINT METHOD 3.1 The reference point problem

The theoretical and numerical challenge is to present the decision maker with an approximation of the Pareto front    ˆ  u ∈ Uad is Pareto optimal ⊂ R2 P = J(u) In order to do so, we follow the ideas laid out in Peitz et al. (2015) and make use of the reference point method : Given a reference point z = (z1 , z2 ) ∈ R2 that satisfies ˆ (4) z < J(u) for all u ∈ Uad we introduce the distance function Fz : U → R by     ˆ − z|2 = 1 Jˆ1 (u) − z1 2 + 1 Jˆ2 (u) − z2 2 . Fz (u) = 1 |J(u) 2

2

2

The mapping Fz measures the geometrical distance beˆ tween J(u) and z. Lemma 2. The mapping Fz is strictly convex. 213

211

2 Proof. The mapping Fz is of the form Fz = i=1 gi ◦ Jˆi where, because of (4), we have gi : (zi , ∞) → R+ 0 with gi (ξ) = (ξ − zi )2 /2. Because of the affine linearity of u → y(u), Jˆ1 is convex and Jˆ2 strictly convex. Further, gi is strictly convex and monotone increasing for i = 1, 2.  Altogether, Fz itself is strictly convex. Suppose that z is componentwise strictly smaller than every objective value which we can achieve within Uad . The goal is that – by approximating z as best as possible ˆ Therefore, we – we get a Pareto optimal point for (P). have to solve the reference point problem ˆ z) min Fz (u) s.t. u ∈ Uad (P

which is a scalar-valued minimization problem.

Theorem 3. For any z ∈ R2 the reference point problem admits a unique solution u ¯z ∈ Uad . Proof. By Lemma 2 the mapping Fz is convex. Now, the proof is identical to the proof of Theorem 2.14 in Tr¨ oltzsch (2010) and uses the strict convexity of Fz along with the fact that Uad is bounded and closed in U.  Theorem 4. Let (4) hold and u ¯z ∈ Uad be an optimal ˆ z ) for a given z ∈ R2 . Then u solution to (P ¯z is Pareto ˆ optimal for (P). Proof. We follow along the lines of Theorem 4.20 in Ehrgott (2005): Assume that u ¯z ∈ Uad is not Pareto ˆ optimal, then there exists a point u ∈ Uad with J(u) ≤ ˆ ˆ ˆ J(¯ uz ) and Jj (u) < Jj (¯ uz ) for j ∈ {1, 2}. Using (4) we get uz ) − zi for i = 1, 2 0 < Jˆi (u) − zi ≤ Jˆi (¯ and strictly smaller for i = j. Together, this yields Fz (u) < Fz (¯ uz ) which is a contradiction to the assumption that u ¯z ˆ z ). is optimal for (P 

ˆ z ) consecutively with an adaptive variation By solving (P of z, we are able to move along the Pareto front in a uniform manner. This way, we get a sequence {z k }k∈N ⊂ R2 of reference points along with optimal controls {uk }k∈N ⊂ ˆ z ) with z = z k as well as {Jˆk }k∈N ⊂ R2 Uad that solve (P k k ˆ ). To be more precise, the next reference with Jˆ = J(u point z k+1 is chosen as Jˆk − Jˆk−1 Jˆk − z k + hz for k ≥ 2, (5) z k+1 = Jˆk + hJ |Jˆk − Jˆk−1 | |Jˆk − z k | where hJ , hz ≥ 0 are chosen to control the coarseness of the approximation to the Pareto front. The algorithm is ˆ initialized by applying the weighted sum method to (P); 1 ˆ2 ˆ Zadeh (1963). This yields the first iterates J , J ∈ P. We therefore do not require z 1 , z 2 and compute z 3 by setting hz = 0 in (5). Note that the algorithm only moves in one direction: If Jˆ11 > Jˆ12 , then it turns to the upper left in the R2 -plane. Therefore, we perform the algorithm twice, the second time with switched roles of Jˆ1 , Jˆ2 to cover the other direction as well. 3.2 Optimality conditions Applying the chain rule, we get for any u ∈ U (6) ∇Fz (u) = (Jˆ1 (u) − z1 )∇Jˆ1 (u) + (Jˆ2 (u) − z2 )∇Jˆ2 (u).

IFAC CPDE 2016 212 June 13-15, 2016. Bertinoro, Italy

Stefan Banholzer et al. / IFAC-PapersOnLine 49-8 (2016) 210–215

It is well known (Hinze et al. (2009)) that the gradients of Jˆ1 and Jˆ2 take the form    (7) ∇Jˆ1 (u) i = χi p(·) dx, 1 ≤ i ≤ m, ∇Jˆ2 (u) = u, Ω

where p ∈ Y is the weak solution to the adjoint equation for (t, x) ∈ Q, −pt (t, x) = ∆p(t, x) ∂p (8) (t, x) = 0 for (t, x) ∈ Σ, ∂n p(T, x) = y(T, x) − yΩ (x) for x ∈ Ω. In particular, we have (9) pY ≤ Cy(T ) − yΩ H for a constant C ≥ 0 which does not depend on u.

3.3 A-posteriori error estimation

We infer from (6)    ∇Fz (u) i = (Jˆ1 (u) − z1 ) χi p(·) dx Ω

+ (Jˆ2 (u) − z2 )ui ∈ U for i = 1, . . . , m. The first-order necessary optimality condition for an optimal u ¯z ∈ Uad now reads as the variational inequality (10) 0 ≤ ∇Fz (¯ uz ), u − u ¯z U for all u ∈ Uad . Next, we investigate second-order derivatives: we find 2   ∇2 Fz (u)v = (Jˆi (u) − zi )∇2 Jˆi (u)v i=1

where



∇2 Jˆ1 (u)v



i

 + ∇Jˆi (u), vU ∇Jˆi (u) for v ∈ U,

=



χi q(·) dx, Ω

∇2 Jˆ2 (u)v = v

and q ∈ Y solves the second adjoint equation for (t, x) ∈ Q, −qt (t, x) = ∆q(t, x) ∂q (t, x) = 0 for (t, x) ∈ Σ, ∂n q(T, x) = y˜(T, x) for (t, x) ∈ Ω and y˜ = Sv ∈ Y is the solution to (2) for u = v and y◦ = 0. We are interested in whether the second derivative of Fz is coercive. Let u ∈ Uad and v ∈ U: ∇2 Fz (u)v, vU 2   2 = (Jˆi (u) − zi )∇2 Jˆi (u)v, vU + ∇Jˆi (u), vU  i=1

≥ (Jˆ1 (u) − z1 )



T

0

  m Ω

2 z2 )vU

k=1

Remark 6. It is a major advantage for the upcoming error estimation that a coercivity constant for ∇2 Fz can be chosen as Jˆ2 (u) − z2 . Usually, the smallest eigenvalue of ∇2 Fz (u) has to be computed to gain information about the coercivity constant; Trenz (2016). From a numerical point of view, this is both very expensive and unstable, whereas the value κz (u) is easily available. 

 χk vk q dxdt

+ (Jˆ2 (u) − 2 2 = (Jˆ1 (u) − z1 )˜ y (T )H + (Jˆ2 (u) − z2 )vU with y˜ = Sv. This yields the following result.

Theorem 5. Let (4) hold. For any u ∈ Uad the hessian ∇2 Fz (u) satisfies 2

∇2 Fz (u)v, vU ≥ κz (u) vU

for κz (u) = Jˆ2 (u) − z2 > 0. Moreover, if    (11) κ ¯ z = min Jˆ2 (u) − z2  u ∈ Uad > 0

˜U , where u ¯z ∈ Uad We want to estimate the error ¯ uz − u ˆ z ) and u is the (unknown) optimal control of (P ˜ ∈ Uad is an given admisible control. We follow along the lines of Kammann et al. (2005) and consider the perturbed problem (12) min Fz (u) + ζ, uU s.t. u ∈ Uad for a fixed function ζ ∈ U. Notice that the functional Fz + ζ, ·U is also strictly convex. Hence, (12) has a unique solution u ¯ζ for any ζ ∈ U. The first-order sufficient optimality condition for (12) reads uζ ) + ζ, u − u ¯ζ U ≥ 0 for all u ∈ Uad . (13) ∇Fz (¯ We will later show that a perturbation ζ = ζ(˜ u) can be computed numerically such that for a known admissible control u ˜ ∈ Uad the variational inequality (13) holds. This implies that u ˜ solves (12). Suppose that (11) holds. Inserting u ˜ in (10) as u and u ¯z in (13) as u yields uz ), u ˜−u ¯z U + ∇Fz (˜ u) + ζ, u ¯z − u ˜ U 0 ≤ ∇Fz (¯ ≤ −∇Fz (¯ uz ) − ∇Fz (˜ u), u ¯z − u ˜U + ζU ¯ uz − u ˜ U . Using the mean value theorem for the Fr´echet-differentiable ˆ ∈ U with function ∇Fz yields the existence of a u u)(¯ uz − u ˜), u ¯z − u ˜U + ζU ¯ uz − u ˜ U 0 ≤ −∇2 Fz (ˆ 2

≤ −¯ κz ¯ uz − u ˜U + ζU ¯ uz − u ˜ U . Summarizing, we deduce the rigorous estimate 1 ˜ U ≤ ζU . ¯ uz − u κ ¯z We still have to identify the function ζ ∈ U. Exactly as in Kammann et al. (2005), we use  u))i (t)]− if u ˜i (t) = (ua )i (t), [(∇Fz (˜ u))i (t)]+ if u ˜i (t) = (ub )i (t), (14) ζi (t) = −[(∇Fz (˜ −(∇Fz (˜ u))i (t) otherwise, where we have used the decomposition for a real number ξ ∈ R as ξ = [ξ]+ − [ξ]− with [ξ]+ = max(0, ξ) and [ξ]− = − min(0, ξ). Thus, we have proved the next theorem.

Theorem 7. Assume that (11) holds and u ¯z ∈ Uad is the ˆ z ). For an arbitrary control u unique solution to (P ˜ ∈ Uad let the function ζ ∈ U be defined as in (14). Then we can estimate the error by 1 ζU . (15) ¯ uz − u ˜U ≤ ∆(˜ u) with ∆(˜ u) = κ ¯z 3.4 Continuity of the mapping z → u ¯z

then ∇2 Fz is uniformly positive definite with coercivity constant κ ¯z . 214

An important aspect is how the optimal control u ¯z adapts to the choice of the reference point.

IFAC CPDE 2016 June 13-15, 2016. Bertinoro, Italy

Stefan Banholzer et al. / IFAC-PapersOnLine 49-8 (2016) 210–215

Theorem 8. Suppose that (11) holds. Let u ¯z be the optiˆ z ) for a given feasible reference point mal solution to (P z ∈ Z with    ¯ for all u ∈ Uad . Z = z˜ ∈ R2  z˜1 < Jˆ1 (u), z˜2 < Jˆ2 (u) − κ Then the mapping z → u ¯z is continuous from Z to Uad .

Proof. Let z1 = (z11 , z12 ), z2 = (z21 , z22 ) ∈ Z be chosen arbitrarily with according optimal controls u ¯1 = u ¯ z1 , u ¯2 = u ¯z2 . By using (10) we get u1 ), u ¯2 − u ¯1 U + ∇Fz2 (¯ u2 ), u ¯1 − u ¯ 2 U 0 ≤ ∇Fz1 (¯ u1 ) − ∇Fz1 (¯ u2 ), u ¯1 − u ¯ 2 U = −∇Fz1 (¯ (16) + ∇Fz2 (¯ u2 ) − ∇Fz1 (¯ u2 ), u ¯1 − u ¯ 2 U . Arguing as in the proof of Theorem 7 the first term is ¯2 2U . Hence, (16) implies that bounded by −¯ κ¯ u1 − u 2

(17) κ ¯ ¯ u1 − u u2 ) − ∇Fz1 (¯ u2 ), u ¯1 − u ¯ 2 U . ¯2 U ≤ ∇Fz2 (¯ From (6) we find that ∇Fz (u) depends linearly on z and ∇Fz2 (¯ u2 ) − ∇Fz1 (¯ u2 ), u ¯1 − u ¯ 2 U 2  = (z1i − z2i ) ∇Jˆi (¯ u2 ), u ¯1 − u ¯ 2 U

Now suppose that we have computed a POD basis {ψi }i=1 ⊂ X of rank ≤ n. The POD Galerkin solution (19)

y  (t) =

  i=1

ai (t)ψi for t ∈ [0, T ], a : [0, T ] → R

satisfies the following POD Galerkin scheme  d  y (t), ψj H + ∇y  (t) · ∇ψj dx dt Ω m  (20) u (t) χ , ψ  , t ∈ (0, T ), = i

i

j H

y  (0), ψj H = y◦ , ψj H for j = 1, ..., and t ∈ (0, T ). It is known that (20) admits a unique solution y  = y  (u) ∈ H 1 (0, T ; V  ) ⊂ Y; see, e.g., Gubisch and Volkwein (2016).

i=1

Inserting (19) into (20) we get the following system in R :

U

i=1

  ≤ C(z2 )z1 − z2 ¯ ¯ 2 U , u1 − u 2 where C(z2 ) = ( i=1 ∇Jˆi (¯ u2 )2U )1/2 is independent of u ¯2 and therefore independent of z2 . We derive from (17):  C(z2 )  ¯2 U ≤ ¯ u1 − u z1 − z2 . κ ¯ Fixing z2 and letting z1 → z2 in R2 now yields u ¯1 → u ¯2 in U which proves the claim.  Remark 9. The set Uad is bounded. Combining (3) and (9) it follows that the weak solution p to (8) is bounded in Y by a constant independent of u ∈ Uad . Consequently, the norms ∇Jˆ1 (u)U and ∇Jˆ2 (u)U are bounded by a constant independent of u ∈ Uad . Hence the constant C(z2 ) in the proof of Theorem 8 is bounded from above by a constant independent of u ∈ Uad . This implies that the mapping z → u ¯z is even Lipschitz-continuous.  4. REDUCED-ORDER MODELING (ROM) BY POD The previously described reference point algorithm makes ˆ z ) which ultimately goes it necessary to repeatedly solve (P back to computing the state equation (1) and its adjoint equation (8) for many different instances. This multiquery context makes model-order reduction techniques conceivable; we focus in particular on proper orthogonal decomposition (POD); see, e.g., Holmes et al. (2012): After having computed the first control u1 by means of a weighted sum problem, a finite-dimensional subspace V  = span {ψ1 , . . . , ψ } ⊂ V is created such that the projected trajectory of y(u1 ) has a least-squares deviation from its full-dimensional version at given time instances 0 = t1 < t2 < . . . < tn = T . Let X denote either the space H or V . Then we consider the POD cost function I : V  → R which is given by n   2     In (ψ1 , ..., ψ ) = αj y(tj ) − y(tj ), ψi X ψi  i=1

with positive (trapezoidal) weights αj . The POD optimization problem then reads (18) min In (ψ1 , ...ψ ) s.t. ψi , ψj X = δij (i, j = 1, ..., ). The solution to (18) is well-known by the singular value decomposition of a certain operator, compare e.g. Kammann et al. (2005).

i=1

2     2 1/2 ∇Jˆi (¯ u2 ), u ¯1 − u ¯2  ≤ z1 − z2 

j=1

213

X

215

M a˙  (t) + Sa (t) = Bu(t), a (0) = a◦ with the mass matrix M ∈ R× , the stiffness matrix S ∈ R× and the control operator B ∈ R×m which are given by for i, j = 1, ..., , Mij = ψi , ψj H  Sij = ∇ψi · ∇ψj dx for i, j = 1, ..., , Ω

Bij = ψj , χi H

for i = 1, ..., m, j = 1, ..., .

The POD Galerkin approximation to (P) is given by the minimization problem ˆ ) (P min Jˆ (u) s.t. u ∈ Uad ,

where we set Jˆ (u) = J(y  (u), u) and y  (u) denotes the solution to (20) for u ∈ Uad .

ˆ  ) we apply the reference point method utilizing To solve (P the corresponding distance function 1 Fz (u) = |Jˆ (u) − z|2 . 2 with a reference point z ∈ R2 . Thereby, we obtain a POD suboptimal control u ¯z ∈ Uad for any z. The resulting error is then estimated using (15) with u ˜=u ¯z . This indicates the quality of the current POD basis, which can then be recomputed if necessary. 5. NUMERICAL EXPERIMENTS 5.1 Setting All computations were carried out on a standard PC with Ubuntu 14.04 LTS, Intel(R) Core i7-4600U CPU @ 2.10GHz x4, 11.4 GiB RAM. For the state equation, we choose a two-dimensional domain Ω = (0, 1)2 and m = 9 shape functions which are indicator functions representing a uniform partition of the domain. The initial condition is y◦ ≡ 0 and we chose

IFAC CPDE 2016 214 June 13-15, 2016. Bertinoro, Italy

Stefan Banholzer et al. / IFAC-PapersOnLine 49-8 (2016) 210–215

1.2

4.5

Pareto points Reference points

1

4

0.8

3.5 0.6

J2

3 0.4

2.5

J2

0.2

2

0

1.5

-0.2

-0.4 -0.05

1

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

J

1

0.5

Fig. 2. Part of the Pareto front with according reference points, chosen close to the front itself.

0 0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

J1

ˆ z ) was solved using a projected Newton-CG Problem (P method; see Kelley (1999). For the coefficients in (5), two different settings were used as will be explained later in detail. The discretization was done using linear finite elements (FE) on a grid with Nx = 729 nodes. The time interval was discretized uniformly using Nt = 100 time instances. The time discretization of the state and adjoint equation were done by an implicit Euler scheme. For the POD basis computation we utilize X = V . 5.2 Results First results can be seen in Figure 1, where the Pareto front was computed using the POD-ROM ansatz described in Section 4 with  = 2 basis functions.

One can observe that the front becomes coarser towards the lower right. This is due to the fact that in this area, Jˆ2 is relatively small which corresponds to little to no control effort. For these objective values, a very small increment in the control variable can yield a large improvement towards the desired state. However, this is due to the nature of the problem and is also observed using different multiobjective optimization techniques, like, e.g., the weighted sum method. Apart from this, the results are quite satisfactory: We have achieved a rather uniform discretization of the Pareto front which can be presented to the decision maker. In Figures 2 and 3, the corresponding reference points can be seen. Note that the axis are not scaled equally so the pictures look distorted. We have utilized two different strategies for choosing the parameters in (5), one resulting in a very tight neighbourhood of the front, the other one in a curve further apart. This has consequences which will become apparent when considering the error estimation. Recall that in the estimation for the real error (15), the term κ ¯ z appears which is given by the condition (11). As 216

1

0.8

0.6

J2

T = 1. The bilateral constraints were set to ua ≡ 0 and ub ≡ 1. For the cost functional J, the terminal condition is yΩ ≡ 1.

Pareto points Reference points

1.2

Fig. 1. Discrete Pareto front consisting of 64 points computed using POD-ROM with  = 2.

0.4

0.2

0

-0.2

-0.4 -0.05

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

J1

Fig. 3. Part of the Pareto front with according reference points, chosen far from the front itself. it turns out, this condition is too strong for our numerical purposes. Instead, we only require that the inequality holds for all u in a neighbourhood of the optimal control u ¯. As an approximation, we therefore compute κz (¯ u ) = Jˆ (¯ u  ) − z2 , z

2

z

ˆ  ), and use is the suboptimal control solving (P where this value instead of κ ¯ z . Furthermore, an additional termination check is implemented to terminate the multiobjective optimization if κz (¯ uz ) falls below a certain tolerance 0 < εκ  1. u ¯z

By (11), if the reference points are further apart from the Pareto front, κ ¯ z should be larger which would correspond to a tighter estimator ∆(˜ u) in (15) for the choice u ˜=u ¯z . We have analyzed the behavior of the error estimator along the multiobjective optimization iteration: For each reference point z ∈ R2 , we obtain a POD-based suboptimal control u ¯z which is an approximation to the (unknown) optimal control u ¯z that is approximated by the FE control denoted by u ¯hz . In Figures 4 and 5, we observe the behavior of the true error ¯ uhz − u ¯z U as well as the error estimator ∆(¯ uz ) in the course of the multiobjective optimization, both for close and far reference points. It can be seen that indeed, we are able to achieve a much tighter estimate for the choice of distant reference points. Nevertheless, the estimator can overshoot the real error

IFAC CPDE 2016 June 13-15, 2016. Bertinoro, Italy

Stefan Banholzer et al. / IFAC-PapersOnLine 49-8 (2016) 210–215

215

REFERENCES True error Error estimator 10

0

10 -2

10

-4

10 -6

10

-8

10

20

30

40 50 Reference point iteration

60

70

80

uz ) Fig. 4. True error ¯ uhz − u ¯z U and error estimator ∆(¯ for  = 2, using reference points z close to the Pareto front.

True error Error estimator 10 -4

10 -5

10

-6

10

20

30 Reference point iteration

40

50

60

Fig. 5. True error ¯ uhz − u ¯z U and error estimator ∆(¯ uz ) for  = 2, using reference points z far from the Pareto front for factors of up to 103 . This can be explained by the fact that we have used some rather rough estimates to get to the bound κ ¯ z in (15). The neglected terms then explain the observed discrepancies. The zigzag behavior of the error estimates can be easily explained since in the implementation, we do not cover both directions in the objective space consecutively. Instead, we alternate between the two directions. We therefore alternate between different regions of the Pareto front where different properties of the reference point problem prevail.

6. CONCLUSION In the present paper it is illustrated that POD reducedorder strategies can be efficiently combined with the reference point method in order to solve bicriterial optimization problems. A rigorous a-posteriori error estimator allows us to control the POD error in order to ensure a desired tolerance. As a next step we plan to involve also convective terms in the state equation in order to model a heat convection in a building. 217

R. Dautray and J.-L. Lions: Mathematical Analysis and Numerical Methods for Science and Technology. Volume 5: Evolution Problems I. Springer-Verlag, Berlin, 2000. M. Ehrgott: Multicriteria Optimization. Springer, Berlin, 2005. K. F. Fong, V. I. Hanby, and T.-T. Chow: HVAC system optimization for energy management by evolutionary programming. Energy and Buildings, 38:220-231, 2006. M. Gubisch and S. Volkwein: Proper orthogonal decomposition for linear-quadratic optimal control. To appear in P. Benner, A. Cohen, M. Ohlberger, and K. Willcox (eds.), Model Reduction and Approximation: Theory and Algorithms. SIAM, Philadelphia, PA, 2016. M. Hinze, R. Pinnau, M. Ulbrich and S. Ulbrich: Optimization with PDE Constraints. Springer, 2009. P. Holmes, J.L. Lumley, G. Berkooz, and C.W. Rowley: Turbulence, Coherent Structures, Dynamical Systems and Symmetry. Cambridge Monographs on Mechanics. Cambridge University Press, Cambridge, 2nd ed. edition, 2012. L. Iapichino, S. Trenz and S. Volkwein: Reduced-order multiobjective optimal control of semilinear parabolic problems. http://nbn-resolving.de/urn:nbn:de:bsz:3520-313874, 2015. L. Iapichino, S. Ulbrich and S. Volkwein: Multiobjective PDE-constrained optimization using the reduced-basis method. http://nbn-resolving.de/urn:nbn:de:bsz:352250190, 2013. E. Kammann, F. Tr¨oltzsch and S. Volkwein: A posteriori error estimation for semilinear parabolic optimal control problems with application to model reduction by POD. ESAIM: Mathematical Modelling and Numerical Analysis, 47:555-581, 2013. C.T. Kelley: Iterative Methods for Optimization, SIAM Frontiers in Applied Mathematics, no 18, 1999. A. Kusiak, F. Tang and G. Xu: Multi-objective optimization of HVAC system with an evolutionary computation algorithm. Energy, 36:2440-2449, 2011. K. Miettinen: Nonlinear Multiobjective Optimization. International Series in Operations Research & Management Science, Springer, 1998. S. Peitz, S. Oder-Bl¨obaum and M. Dellnitz: Multiobjective optimal control methods for fluid flow using reduced order modeling. http://arxiv.org/pdf/1510.05819v2.pdf, 2015. C. Romaus, J. B¨ocker, K. Witting, A. Seifried, and O. Znamenshchykov: Optimal energy management for a hybrid energy storage system combining batteries and double layer capacitors. IEEE, pages 1640-1647, San Jose, CA, USA, 2009. W. Stadler: Multicriteria Optimization in Engineering and in the Sciences. Plenum Press, New York, 1988. S. Trenz: A-posteriori Error Estimates for Reduced Order Models in Nonlinear Optimization Problems. Ph.D. thesis, University of Konstanz, 2016. F. Tr¨oltzsch: Optimal Control of Partial Differential Equations. Theory, Methods and Applications. American Math. Society, Providence, volume 112, 2010. L. Zadeh. Optimality and non-scalar-valued performance criteria. IEEE Transactions on Automatic Control, 8, 1963.