Reduced order methods for parametrized non-linear and time dependent optimal flow control problems, towards applications in biomedical and environmental sciences
Maria Strazzullo, Zakia Zainib, Francesco Ballarin, Gianluigi Rozza
RReduced order methods for parametrized non-linear and timedependent optimal flow control problems, towards applications inbiomedical and environmental sciences
Maria Strazzullo (cid:93) , Zakia Zainib (cid:93) , Francesco Ballarin (cid:93) , and Gianluigi Rozza (cid:93) (cid:93) mathlab, Mathematics Area, International School for Advanced Studies (SISSA), Via Bonomea 265, I-34136Trieste, Italy
Abstract
We introduce reduced order methods as an efficient strategy to solve parametrized non-linear andtime dependent optimal flow control problems governed by partial differential equations. Indeed,the optimal control problems require a huge computational effort in order to be solved, most of allin physical and/or geometrical parametrized settings. Reduced order methods are a reliable andsuitable approach, increasingly gaining popularity, to achieve rapid and accurate optimal solutionsin several fields, such as in biomedical and environmental sciences. In this work, we employ aPOD-Galerkin reduction approach over a parametrized optimality system, derived from the Karush-Kuhn-Tucker conditions. The methodology presented is tested on two boundary control problems,governed respectively by ( i ) time dependent Stokes equations and ( ii ) steady non-linear Navier-Stokesequations. Parametrized optimal flow control problems (OFCP( µ )s) constrained to parametrized partial differen-tial equations (PDE( µ )s) are a very versatile mathematical model which arises in several applications,see e.g. [7, 9, 4]. These problems are computationally expensive and challenging even in a sim-pler non-parametrized context. The computational cost becomes unfeasible when these problemsinvolve time dependency [1, 14] or non-linearity [5, 6, 4], in addition to physical and/or geometricalparametrized settings that describe several configurations and phenomena. A suitable strategy tolower this expensive computational effort is to employ reduced order methods (ROMs) in the contextof OCP( µ )s, which recast them in a cheap, yet reliable, low dimensional framework [8, 12]. We ex-ploit these techniques in order to solve boundary OCP( µ )s on a bifurcation geometry [13] which canbe considered as ( i ) a riverbed in environmental sciences and as ( ii ) a bypass graft for cardiovascularapplications. In the first research field, reduced parametrized optimal control framework (see e.g.[10, 11]) can be of utmost importance. It perfectly fits in forecasting and data assimilated models andit could be exploited in order to prevent possibly dangerous natural situations [15]. The presented testcase is governed by time dependent Stokes equations, which are an essential tool in marine sciencesin order to reliably simulate evolving natural phenomena.Furthermore, discrepancies between computational modelling in cardiovascular mechanics and realityusually ought to high computational cost and lack of optimal quantification of boundary conditions,especially the outflow boundary conditions. In this work, we present application of the aforementionednumerical framework combining OFCP( µ ) and reduced order methods in the bifurcation geometry.The aim is to quantify the outflow conditions automatically while matching known physiological datafor different parameter-dependent scenarios [2, 17]. In this test case, Navier-Stokes equations willmodel the fluid flow.The work is outlined as follows: in section 2, the problem formulation and the methodology aresummarized. Section 3 shows the numerical results for the two boundary OCP( µ )s, based on [9, 13].Conclusions follow in section 4. a r X i v : . [ m a t h . NA ] M a y Proper Orthogonal Decomposition for OCP( µ )s In this section, we briefly describe the problem and the adopted solution strategy for time dependentnon-linear boundary OCP( µ )s: in the cases mentioned in section 1, the reader shall take the non-linearterm and time-dependent terms to be zero accordingly [16, 17]. The goal of OCP( µ )s is to find aminimizing solution for a quadratic cost functional J under a PDE( µ ) constraint thanks to an externalvariable denoted as control . In the next section, we will show numerical results over a bifurcationgeometry Ω with physical and/or geometrical parametrization represented by µµµ ∈ D ⊂ R d , d ∈ N .Thus, considering the space-time domain Q = Ω × [ , T ] with a sufficiently regular spatial boundary ∂ Ω ∗ , let us define the Hilbert spaces S = V × P , Z = Z V × Z P and U for state and adjoint velocityand pressure, and control variables denoted by sss = ( vvv , p ) ∈ S , uuu ∈ U and zzz = ( , q ) ∈ Z , respectively.The stability and uniqueness of the optimal solution will be guaranteed if S ≡ Z , which will be ourassumption in this work. We introduce X = S × U such that x = ( sss , uuu ) ∈ X . Then, the problem reads: given µµµ ∈ D , find ( x , zzz ; µµµ ) ∈ X × S s.t.: A ( x , y ; µµµ ) + B ( zzz , y ; µµµ ) + E ( vvv , , y vvv ; µµµ ) + E ( , vvv , y vvv ; µµµ ) = ∫ T (cid:104)H ( µµµ ) , y (cid:105) dt , ∀ y ∈ X , B ( x , κκκ ; µµµ ) + E ( vvv , vvv , κκκ ; µµµ ) = ∫ T (cid:104)G ( µµµ ) , κκκ (cid:105) dt , ∀ κκκ ∈ S , (1)where y = (cid:0) y vvv , y p , y uuu (cid:1) and κκκ = ( κκκ , κκκ q ) . The bilinear forms A : X × X → R and H ( µµµ ) are relatedto the minimization cost functional J , while the bilinear form B : X × Z → R represents the linearpart of state-constraints and E is the non-linear convection term, which will be zero for the timedependent linear case.In order to solve the optimality system (1) we exploit Galerkin Finite Elements (FE) snapshots-based Proper orthogonal decomposition (POD)–Galerkin (see [8]), summarized in table 1, where thenumber of time steps are denoted by N t . Offline phase: Input: µ for lifting, N , n .Output: Reduced order solution spaces.1 . Compute snapshots δ N δ ×N t ( µ n ) for 1 ≤ n ≤ | Λ | , δ = v , p , u , w , q and state and adjointsupremizers. The global dimension of FE space discretization is N = N v + N p + N u .2 . Solve eigenvalue problems A δ ρ δ n = λ δ n ρ δ n , n = , · · · , | Λ | , where A δ is correlationmatrix of snapshots.3 . If relative energy of eigenvalues is greater than 1 − (cid:15) tol , < (cid:15) tol (cid:28)
1, keep correspondingeigenvalue-eigenvector pairs ( λ n , ρ n ) .4 . Construct orthonormal POD basis from the retained N eigenvectors and add the POD modesof the supremizers to state and adjoint velocities.Online phase: Input: Online parameter µ ∈ D .Output: Reduced order solution.1 . Perform Galerkin projection to calculate reduced order coefficients such that δδδ ≈ X δδδ δδδ N where, X δδδ denotes reduced bases matrices containing all the time instances.2 . Solve the reduced order version of the optimality system (1).
Table 1:
Algorithm: POD–Galerkin for OCP( µ )In order to guarantee the efficiency of the POD–Galerkin approach, we rely on the affine assumptionover the forms, i.e. every form can be written as a linear combination of µ − dependent functions and µ − independent quantities. In this way, the system resolution is divided into parameter independent( offline ) and dependent ( online ) phases (see table 1 for details) such that the expensive calculationsare absorbed in the former stage and only online stage is repeated every time the parameter µµµ changes.From the perspective of the problem stability, to ensure uniqueness of pressure at the reduced orderlevel, we enrich the state and adjoint velocity space with supremizers and, to guarantee the fulfillmentof Brezzi’s inf-sup condition [3] at the reduced level, we use aggregated equivalent state and adjointspaces. Thus, dimension of the reduced problem reduces from N × N t to 13 N . ∗ For the steady case, T = and Q ≡ Ω Results µ ) governed by Stokes Equations In this section, inspired by [9, 13], we propose an OFCP( µµµ ) governed by a time dependent Stokesequation. First of all, let us introduce the smooth domain Ω ( µ ) . The parameter stretches the lengthof the reference domain shown in figure 2, which will be indicated with Ω from now on. We wantto recover a measurement vvv d ( µ ) ∈ L ( , T ; [ L ( Ω )] ) over the one dimensional observation domain Γ OBS controlling the Neumann flux over Γ C , with the inflow v in ( µ ) = ( ( x − ) ( − x ) , ) . Thesetting is suited for environmental applications: we control the flow in order to avoid potentiallydangerous situations in an hypothetical real time monitoring plan on the domain, which can representa riverbed. The space-time domain is Q = Ω × [ , ] . Let us consider the following function spaces: V = L ( , T ; [ H Γ D ( Ω )] ) ∩ H ( , T ; [ H Γ D ( Ω ) ∗ ] ) , P = L ( , T ; L ( Ω )) and U = L ( , T ; [ L ( Ω )] ) for state and adjoint velocity, state and adjoint pressure and for control, respectively. Then, we define X = ( V × P ) × U . For a given µ ∈ D = [ . , ] × [ , ] × [ . , ] , we want to find the solution oftime dependent Stokes equations which minimizes: J : = ∫ T ∫ Γ OBS ( vvv ( µ ) − vvv d ( µ )) dsdt + α ∫ T ∫ Γ C uuu ( µ ) dsdt + α ∫ T ∫ Γ C |∇ uuu ( µ ) ttt | dsdt , (2)where α = − , α = − and ttt is the unit tangent vector to Γ C and vvv d ( µ ) = [ µ ( ( x − x − x + ) + (− x − x + x + )) , ] . The cost functional penalizes not only the magnitude of the control, butalso its rapid variations over the boundary. The constrained minimization problem (2) is equivalentto the resolution of problem (1) where the considered forms are defined by: A( x , y ) = ∫ T ∫ Γ OBS vvv · y vvv dsdt + α ∫ T ∫ Γ C uuu · y uuu dsdt + α ∫ T ∫ Γ C ∇ uuuttt · ∇ y uuu ttt dsdt , B( x , z ; µ ) = ∫ Q ∂ vvv ∂ t · w dxdt + µ ∫ Q ∇ vvv · ∇ w dxdt − ∫ Q p (∇ · w ( µ )) dxdt − ∫ Q q (∇ · ( v ( µ )) dxdt − ∫ T ∫ Γ C uuu · w dsdt , (cid:104)H ( µ ) , y (cid:105) = ∫ Γ OBS vvv d ( µ ) · y vvv ds , (cid:104)G( µ ) , q (cid:105) = , ∀ y ∈ X , for every x , y ∈ X and κ ∈ S . We built the reduced space with N =
35 over a training set of 70snapshots of global dimension 131400, for N t =
20. In time dependent applications, ROMs are ofgreat advantage: in table 2 the speedup index is shown with respect to N . The speedup representshow many ROM systems one can solve in the time of a FE simulation. Nevertheless, we do notpay in accuracy as figure 1 and figure 3 show: it represents the relative error between FE and ROMvariables. The relative error between FE and ROM J is presented in table 2 Table 2:
Speedup analysis and relative error J . N Speedup Relative error J
15 66338 10 −
20 47579 10 −
25 34335 10 −
30 22477 10 −
35 17420 10 − igure 1: FE vs ROM mean relative error over 50 param-eters.
Figure 2:
Physical domain.
Figure 3:
FE (top) vs ROM (bottom) comparison of state velocity and state pressure, for t=0.05,0.5,1 and µ = ( . , . , ) . µ ) governed by Navier-Stokes equations In this section, we will demonstrate the numerical results for second test case with optimal boundarycontrol problem governed by non-linear incompressible steady Navier-Stokes equations. We considera bifurcation domain Ω as employed in the previous example (see figure 4), which can be consideredas an idealized model of arterial bifurcation in cardiovascular problems [9, 13, 17]. Fluid shall enterthe domain from Γ in and shall leave through the outlets Γ c . In this example, physical parameterizationis considered for the inflow velocity given by v in ( µ ) = µ ( x ( − x ) , ) and the desired velocity,denoted by v d ∈ L ( Ω ) and prescribed at the 1-D observation boundary Γ obs through the followingexpression: v d ( µ ) = (cid:32) µ (cid:16) . (cid:16) ( x − ) − ( x − ) − ( x − ) + (cid:17) + . (cid:16) − ( x − ) − ( x − ) + x (cid:17)(cid:17) (cid:33) . The cost-functional J is defined as: J ( v , u ; µ ) = (cid:107) v ( µ ) − v d ( µ ) (cid:107) L ( Γ obs ) + α (cid:107) u ( µ ) (cid:107) L ( Γ c ) + . α (cid:107)∇ u ( µ ) t (cid:107) L ( Γ c ) , (3)where t is the tangential vector to Γ c . The mathematical problem reads: Given µ ∈ D = [ . , . ] , find ( v ( µ ) , p ( µ ) , u ( µ )) that minimize J and satisfy the Navier-Stokes equations with v in ( µ ) prescribedat the inlet Γ in , no-slip conditions at the walls Γ w and u ( µ ) implemented at Γ c through Neumannconditions. At the continuous level, we consider X ( Ω ) = H Γ in ∪ Γ w ( Ω ) × L ( Ω ) × L ( Γ c ) , where H Γ in ∪ Γ w ( Ω ) = (cid:104) v ∈ (cid:2) H ( Ω ) (cid:3) : v | Γ in = v in and v | Γ w = (cid:105) . igure 6: Eigenvalues of N =
10 POD modes.
Figure 7:
Relative error for solution variables and J .Thus, A ( x , y ) = ∫ Ω v ( µ ) · y v d Ω + α ∫ Γ c u ( µ ) · y u d Γ c + α ∫ Γ c (∇ u ( µ )) t · ∇ ( y u ) t d Γ c , B ( x , z ) = η ∫ Ω ∇ v ( µ ) · ∇ w d Ω − ∫ Ω p ( µ ) (∇ · w ) d Ω − ∫ Ω q (∇ · v ( µ )) d Ω − ∫ Γ c u ( µ ) · w d Γ c , E ( v , v , w ) = ∫ Ω ( v ( µ ) · ∇) v ( µ ) · w d Ω and (cid:104)H ( µ ) , y (cid:105) = ∫ Ω v d ( µ ) · y v d Ω . Figure 4:
Domain ( Ω ) . Mesh size
No. of reduced order bases N D [ . , . ]| Λ | offline phase . × secondsonline phase × seconds Figure 5:
Computational details of POD–Galerkinfor Navier-Stokes constrained OCP( µ ). | v || u | Figure 8:
State velocity and control for µ = . , . , . ( Λ ) of 100 parameter valuesand solving the problem (1) through Galerkin Finite Element method, we construct the snapshot atrices for the solution variables v , p , u , w , q . For N =
10, eigenvalues energy of the state, controland adjoint variables is demonstrated in figure 6. Evidently, N eigenvalues capture 99 .
9% of theGalerkin FE discretized solution spaces and the reduced order spaces are thus built with dimensions13N + = µ = . , . , . J infigure 7. The former decreases upto 10 − along with the latter decreasing upto 10 − . In this work, we propose ROMs as a suitable tool to solve a parametrized boundary OCP( µ )s for timedependent Stokes equations and steady Navier-Stokes equations. The framework proposed is suitedfor several many query and real time applications both in environmental marine sciences and bio-engineering. The reduction of the KKT system is performed through a POD-Galerkin approach, whichleads to accurate surrogate solutions in a low dimensional space. This work aims at showing howROMs can have an effective impact in the management of parametrized simulations for social life andactivities, such as coastal engineering and cardiovascular problems. Indeed, the proposed frameworkdeals with faster solving of parametrized optimal solutions which can find several applications inmonitoring planning both in marine ecosystems and patient specific geometries. Acknowledgements
We acknowledge the support by European Union Funding for Research and Innovation – Horizon2020 Program – in the framework of European Research Council Executive Agency: ConsolidatorGrant H2020 ERC CoG 2015 AROMA-CFD project 681447 “Advanced Reduced Order Methodswith Applications in Computational Fluid Dynamics” (P.I. Prof. G. Rozza). We also acknowledge theINDAM-GNCS project “Advanced intrusive and non-intrusive model order reduction techniques andapplications”. The computations in this work have been performed with RBniCS library, developedat SISSA mathLab, which is an implementation in FEniCS of several reduced order modellingtechniques; we acknowledge developers and contributors to both libraries.
References [1] V. Agoshkov, A. Quarteroni, and G. Rozza. A mathematical approach in the de-sign of arterialbypass using unsteady stokes equations.
Journal of Scientific Computing , 2006.[2] F. Ballarin, E. Faggiano, S. Ippolito, A. Manzoni, A. Quarteroni, G. Rozza, and R. Scrofani.Fast simulations of patient-specific haemodynamics of coronary artery bypass grafts based on apod–galerkin method and a vascular shape parametrization.
Journal of Computational Physics ,315:609–628, 2016.[3] F. Brezzi. On the existence, uniqueness and approximation of saddle-point problems arising fromLagrangian multipliers.
Revue française d’automatique, informatique, recherche opérationnelle.Analyse numérique , 8(2):129–151, 1974.[4] J. C. De los Reyes and F. Tröltzsch. Optimal control of the stationary Navier-Stokes equationswith mixed control-state constraints.
SIAM Journal on Control and Optimization , 46(2):604–629, 2007.[5] A. V. Fursikov, M. D. Gunzburger, and L. Hou. Boundary value problems and optimal boundarycontrol for the Navier–Stokes system: the two-dimensional case.
SIAM Journal on Control andOptimization , 36(3):852–894, 1998.[6] M. D. Gunzburger, L. Hou, and Th. P. Svobodny. Analysis and finite element approximationof optimal control problems for the stationary Navier-Stokes equations with distributed andNeumann controls.
Mathematics of Computation , 57(195):123–151, 1991.[7] J. Haslinger and R. A. E. Mäkinen.
Introduction to shape optimization: theory, approximation,and computation . SIAM, Philadelphia, 2003.
8] J. S. Hesthaven, G. Rozza, and B. Stamm. Certified reduced basis methods for parametrizedpartial differential equations.
SpringerBriefs in Mathematics , 2015, Springer, Milano.[9] F. Negri, A. Manzoni, and G. Rozza. Reduced basis approximation of parametrized optimalflow control problems for the Stokes equations.
Computers & Mathematics with Applications ,69(4):319–336, 2015.[10] A. Quarteroni, G. Rozza, L. Dedè, and A. Quaini. Numerical approximation of a control problemfor advection-diffusion processes. In
IFIP Conference on System Modeling and Optimization ,pages 261–273, Ceragioli F., Dontchev A., Futura H., Marti K., Pandolfi L. (eds) SystemModeling and Optimization. CSMO 2005. vol 199. Springer, Boston, 2005.[11] A. Quarteroni, G. Rozza, and A. Quaini. Reduced basis methods for optimal control of advection-diffusion problems. In
Advances in Numerical Mathematics , pages 193–216. RAS and Universityof Houston, 2007.[12] G. Rozza, D.B.P. Huynh, and A. Manzoni. Reduced basis approximation and a posteriori errorestimation for Stokes flows in parametrized geometries: Roles of the inf-sup stability constants.
Numerische Mathematik , 125(1):115–152, 2013.[13] G. Rozza, A. Manzoni, and F. Negri. Reduction strategies for PDE-constrained oprimizationproblems in Haemodynamics. pages 1749–1768, ECCOMAS, Congress Proceedings, Vienna,Austria, September 2012.[14] M. Stoll and A. Wathen. All-at-once solution of time-dependent Stokes control.
J. Comput.Phys. , 232(1):498–515, January 2013.[15] M. Strazzullo, F. Ballarin, R. Mosetti, and G. Rozza. Model reduction for parametrized optimalcontrol problems in environmental marine sciences and engineering.
SIAM Journal on ScientificComputing , 40(4):B1055–B1079, 2018.[16] M. Strazzullo, F. Ballarin, and G. Rozza. POD-galerkin model order reduction for parametrizedtime dependent linear quadratic optimal control problems in saddle point formulation. submitted.https://arxiv.org/abs/1909.09631., 2019.[17] Z. Zainib, F. Ballarin, G. Rozza, P. Triverio, L. Jiménez-Juan, and Fremes. S. Reduced ordermethods for parametric optimal flow control in coronary bypass grafts, towards patient-specificdata assimilation. submitted, https://arxiv.org/abs/1911.01409., 2019., 40(4):B1055–B1079, 2018.[16] M. Strazzullo, F. Ballarin, and G. Rozza. POD-galerkin model order reduction for parametrizedtime dependent linear quadratic optimal control problems in saddle point formulation. submitted.https://arxiv.org/abs/1909.09631., 2019.[17] Z. Zainib, F. Ballarin, G. Rozza, P. Triverio, L. Jiménez-Juan, and Fremes. S. Reduced ordermethods for parametric optimal flow control in coronary bypass grafts, towards patient-specificdata assimilation. submitted, https://arxiv.org/abs/1911.01409., 2019.