Minimal curvature trajectories: Riemannian geometry concepts for slow manifold computation in chemical kinetics
Research Highlights
► Extremum principle for trajectories near slow attracting manifolds. ► Allows formulation of optimization problem for computing slow manifolds. ► Efficient numerical solution of optimization problem possible. ► Curvature-based geometric criteria serve as optimization criteria. ► Geometric point of view shows analogy to fundamental curvature-force.
Introduction
The need for model reduction in chemical kinetics is mainly motivated by the fact that the computational effort for a full simulation of reactive flows, e.g. of fluid transport involving multiple time scale chemical reaction processes, is extremely high. For detailed chemical reaction mechanisms involving a large number of chemical species and reactions, a simulation in reasonable computing time requires reduced models of chemical kinetics [1].
However, model reduction is often also of general interest for theoretical purposes in mathematical modeling. Reduced models are intended to describe some essential characteristics of dynamical system behavior while reducing the state space dimension. Therefore, they often allow a better insight into complicated reaction pathways, e.g. in biochemical systems [2], and their nonlinear dynamics.
In dissipative ordinary differential equation systems modeling chemical reaction kinetics different time scales cause anisotropic phase volume contraction along solution trajectories. This leads to a bundling of trajectories near “manifolds of slow motion” of successively lower dimension as time progresses, illustrated in Fig. 1. Many model reduction methods exploit this for simplifying chemical kinetics via a time scale separation into fast and slow modes. The aim is to approximate the system dynamics with a dimension-reduced model after eliminating the fast modes by enslaving them to the slow ones via computation of slow attracting manifolds.
Very early model reduction approaches like the quasi steady-state (QSSA) and partial equilibrium assumption (PEA) [1] performed “by hand”, have set the course for modern numerical model reduction methods that automatically compute a reduced model without need for detailed expert knowledge of chemical kinetics by the user. Many of these modern techniques are explicitly or implicitly based on a time-scale analysis of the underlying ordinary differential equation (ODE) system with the purpose to identify a slow attracting manifold in phase space where – after a short initial relaxation period – the system dynamics evolve. For a comprehensive overview see e.g. [3], and references therein.
In 1992 Maas and Pope introduced the ILDM-method [4] which became very popular and widely used in the reactive flows community, in particular in combustion applications. Based on a singular perturbation approach, a local time-scale analysis is performed on the Jacobian of the system of ordinary differential equations modeling chemical kinetics. Fast time scales are assumed to be fully relaxed and after a suitable coordinate transformation fast variables are computed as a function of the slow ones by solving a nonlinear algebraic equation system. For recent developments and extensions of the ILDM method see e.g. [5], and references therein.
Another popular technique is computational singular perturbation (CSP) method proposed by Lam [6], [7]. The basic concept of this method is a representation of the dynamical system in a set of “ideal” basis vectors such that fast and slow modes are decoupled. Some iterative methods came into application that are based on an evaluation of functional equations suitably describing the central characteristics of a slow attracting manifold, for example invariance and stability. Examples are Fraser’s algorithm[8], [9], [10] and the method of invariant grids [11], [12], [3]. Other widely known and applied methods are e.g. the constrained runs algorithm [13], [14], the rate-controlled constrained equilibrium (RCCE) method [15], the invariant constrained equilibrium edge preimage curve (ICE-PIC) method [16], [17], and flamelet-generated manifolds [18], [19]. In [20] finite time Lyapunov exponents and vectors are analyzed for evaluation of timescale information. Mitsos et al. formulated an integer linear programming problem explicitly minimizing the number of species in the reduced model subject to a given error constraint [21]. In [22] the authors connect fixed points via heteroclinic orbits in a two-dimensional system for identifying the one-dimensional slow manifold. In particular, there seems to be a connection between the methods of Gear et al. [13], [23], Adrover et al. [24] and our method which might be analyzed in future work.
It is obvious that non-local information on phase space dynamics should be taken into account to get accurate approximations of slow attracting manifolds in the general case. Reaction trajectories in phase space that are solutions of the ODE system describing chemical kinetics and uniquely determined by their initial values bear such information. Based on Lebiedz’ idea to search for an extremum (variational) principle that distinguishes trajectories on or near slow attracting manifolds, we apply an optimization approach for computing such trajectories [25], [26]. Various optimization criteria have been suggested [27], systematically investigated and the trajectory-based approach has been extended to the computation of manifolds of arbitrary dimension via parameterized families of trajectories.
In this context it is important to remark that equations of motion in classical mechanics can also be derived from a variational principle, Hamilton’s principle of least action. In Lagrange–Hamilton mechanics, see e.g. [28], the trajectory of a mechanical system is determined in such a way that the action (which is defined as the integral of the Lagrangian over time) is minimal where the Lagrangian is the difference of kinetic energy and potential energy. Whereas the Hamiltonian formalism leads to a provable exact description of classical mechanics, the variational principle we propose here is empirical and leads to approximations of slow invariant manifolds in dissipative dynamical systems.
This paper derives and comprehensively discusses various geometrically motivated objective criteria for computing trajectories approximating slow attracting manifolds in chemical kinetics as a solution of an optimization problem. The corresponding objective functionals are supposed to implicitly incorporate essential characteristics of slow attracting manifolds related to a minimal remaining relaxation of chemical forces along trajectories on these manifolds. We consider the picture of abstract “dissipative chemical forces” imagined to drive the single elementary reaction steps [29]. Due to the second law of thermodynamics these forces successively relax while the chemical system is approaching equilibrium. The successive relaxation of these forces causes curvature in the reaction trajectories (in the sense of velocity change along the trajectory). A slow 1-D manifold in this picture would correspond to a minimally curved reaction trajectory along which the remaining relaxation of chemical forces is minimal while approaching chemical equilibrium.
In particular, we propose and motivate an optimization criterion suitably measuring curvature which is rooted in a thermodynamically motivated Riemannian geometry specifically defined for chemical reaction kinetics and based on the second law of thermodynamics. This metric provides the phase space of chemical reaction kinetics with a geometry specifically capturing the structure of chemical kinetic systems [30].
The proposed method for the computation of a slow manifold is automatic. The user has to provide only the desired dimension of the reduced model and the range of concentrations of the reaction progress variables supposed to parameterize the reduced model. The choice of number and kind of species for the reaction progress variables parametrizing the manifold should not be made arbitrarily, see e.g. [31]. In the present manuscript we choose products and reactants of the overall reactions, which monotonously increase/decrease when the reaction progresses, as it is often done when using model reduction in computational fluid dynamics of combustion processes. From a theoretical point of view it would be desirable to identify the number of reaction progress variables in a way guaranteeing a large time gap and the kind of chemical species such that they parameterize the slow manifold in a locally unique sense. This can be done using e.g. CSP-pointers and other methods mentioned in [31] or methods for complexity reduction and time scale analysis [2], [32], [26]. A bad choice of the number and kind of the reaction progress variables (where a good choice can vary in phase space) may lead to problems with the identification of the slow manifold. Our method aims at computing the manifold as good as possible under the given constraints, but it cannot overcome a potential constraint concerning dimensionality and choice of progress variables provided by the user that stands in opposition to the intrinsic properties of the dynamics of the model. In any case the choice of manifold dimensionality and choice of parameterization by the reaction progress variables has to be made by the user before starting our computational algorithm.
For the application examples presented in this work, the numerical optimization algorithm shows fast convergence independent of the initial values chosen for numerical initialization. The optimization problems seem to be convex for the example systems presented in this paper (see “optimization landscapes” in Sections 3.1.2 Optimization landscapes, 3.4.3 Optimization landscapes) and would then have a unique solution corresponding to a global minimum of the objective functional. An advantage of the presented trajectory optimization approach over local time scale separation methods like QSSA and ILDM is the fact that it produces smooth manifolds and whole 1-D manifolds (trajectories) as a solution of a single run of the optimization algorithm. Methods based on explicit local time scale separation might yield non-smooth manifolds when the fast–slow spectral decomposition changes its structure. Furthermore, the formulation as an optimization problem assures results even under conditions where the time scale separation is small and many common slow manifold computation methods fail or numerical solutions are difficult to obtain.
Section snippets
Trajectory-based optimization approach
As described in the introduction, the key of the method presented here is the exploitation of global phase space information contained in the behavior of trajectories on their way towards chemical equilibrium. This information can be used within an optimization framework for identifying suitable reaction trajectories approximating slow invariant and attracting manifolds (SIM). A suitable formulation as the numerical solution of an optimization problem assures the existence of a reduced model
Results
In this section, results for our slow manifold computation method based on trajectory optimization are presented. We choose four different chemical reaction mechanisms to demonstrate its application: the Davis–Skodje model system, a simple nonlinear three-component reaction mechanism, a simplified reaction mechanism for the combustion of H2, and a realistic temperature dependent mechanism for ozone decomposition. For these four mechanisms all previously discussed objective functionals (5), (7),
Summary, discussion, and outlook
We present various geometrically motivated criteria for the numerical computation of trajectories approximating slow (attracting) invariant manifolds (SIM) in chemical reaction kinetics. The key idea of our approach is to approximately span the SIM by trajectories being solutions of an optimization problem for initial values of these trajectories. The objective functional is supposed to characterize the extent of relaxation of chemical forces (being minimal in the optimal solution) along a
Acknowledgments
This work was supported by the German Research Foundation (DFG) through the Collaborative Research Center (SFB) 568. The authors also thank Prof. Dr. Dr. h.c. Hans Georg Bock (IWR, Heidelberg) for providing the MUSCOD-II-package along with DAESOL and Prof. Dr. Moritz Diehl (K.U. Leuven/Belgium) for support with the initial value embedding add-on for MUSCOD-II. We especially thank Andreas Potschka (IWR, Heidelberg) for continuous support concerning MUSCOD-II and many helpful hints related to
References (55)
- et al.
Simplifying chemical kinetics: intrinsic low-dimensional manifolds in composition space
Combust. Flame
(1992) - et al.
Rate-controlled partial-equilibrium method for treating reacting gas mixtures
Combust. Flame
(1971) - et al.
Optimal automatic reaction and species elimination in kinetic mechanisms
Combust. Flame
(2008) - et al.
Stretching-based diagnostics and reduction of chemical kinetic models with diffusion
J. Comput. Phys.
(2007) - et al.
Constructive methods of invariant manifolds for kinetic problems
Phys. Rep.
(2004) - et al.
CSP analysis of a transient flame-vortex interaction: time scales and manifolds
Combust. Flame
(2003) - et al.
An efficient multiple shooting based reduced SQP strategy for large-scale dynamic process optimization. Part I: theoretical aspects
Comput. Chem. Eng.
(2003) - et al.
An efficient multiple shooting based reduced SQP strategy for large-scale dynamic process optimization. Part II: Software aspects and applications
Comput. Chem. Eng.
(2003) - et al.
Real-time optimization and nonlinear model predictive control of processes governed by differential-algebraic equations
J. Process. Contr.
(2002) - et al.
Quasi-equilibrium grid algorithm: geometric construction for model reduction
J. Comput. Phys.
(2008)
The application of chemical reduction methods to a combustion system exhibiting complex dynamics
Proc. Combust. Inst.
Combustion: Physical and Chemical Fundamentals, Modeling and Simulation, Experiments, Pollutant Formation
Automatic network coupling analysis for dynamical systems based on detailed kinetic models
Phys. Rev. E
Invariant Manifolds for Physical and Chemical Kinetics
Simple global reduction technique based on decomposition approach
Combust. Theor. Model.
Singular perturbation for stiff equations using numerical methods
The CSP method for simplifying kinetics
Int. J. Chem. Kinet.
Geometric investigation of low-dimensional manifolds in systems approaching equilibrium
J. Chem. Phys.
The steady state and equilibrium approximations: a geometrical picture
J. Chem. Phys.
Geometrical picture of reaction in enzyme kinetics
J. Chem. Phys.
Comparison of invariant manifolds for model reduction in chemical kinetics
Commun. Comput. Phys.
Invariant grids: method of complexity reduction in reaction networks
Complexus
Projecting to a slow manifold: singularly perturbed systems and legacy codes
SIAM J. Appl. Dyn. Syst.
Analysis of the accuracy and convergence of equation-free projection to a slow manifold
ESAIM: Math. Model. Numer. Anal.
The invariant constrained equilibrium edge preimage curve method for the dimension reduction of chemical kinetics
J. Chem. Phys.
Cited by (24)
Flow curvature manifold and energy of generalized Liénard systems
2022, Chaos, Solitons and FractalsCitation Excerpt :This method gives an implicit non intrinsic equation, because it depends on the euclidean metric. A ‘kinetic energy metric’ has been introduced in [19] for chemical kinetic systems and an extremum principle for computing slow invariant manifolds has been formulated [20,21] which can be viewed as minimum curvature geodesics. In [14] a curvature-based differential geometry formulation for the slow manifold problem has been used for the purpose of a coordinate-independent formulation of the invariance equation.
Investigation on Ginzburg-Landau equation via a tested approach to benchmark stochastic Davis-Skodje system
2021, Alexandria Engineering JournalCitation Excerpt :These methods are based on the exponential Milstein scheme. Davis and Skodje introduced a model includes of a two-dimensional system to compare various reduction methods [48,49]. The Davis-Skodje mechanism is a measurement criterion for model reduction methods because of the adjustable time scale separation.
Approximation of slow and fast dynamics in multiscale dynamical systems by the linearized Relaxation Redistribution Method
2012, Journal of Computational PhysicsDevelopment of reduced and optimized reaction mechanisms based on genetic algorithms and element flux analysis
2012, Combustion and FlameCitation Excerpt :Besides the classical quasi steady-state (QSS) and partial equilibrium approximations (PE) [37–39], and computational singular perturbation (CSP) [40,41], which, allowing fast and slow time scales analysis, has extensively been used also for generation of reduced and skeletal mechanisms [42–45], other methodologies have been recently developed, such as intrinsic, trajectory-generated, constraint-defined low-dimensional manifolds (ILDM, TGLDM, CDLDM) [46–48], rate-controlled constrained equilibrium (RCCE) [49,50], invariant manifold methods (MIM) [51], and invariant constrained equilibrium edge preimage curve method (ICE-PIC) [52]. In the minimal-curvature-trajectory based approach, an optimization procedure is adopted for the generation of optimal low-dimensional manifolds or 1D trajectories [53]. The reduction and optimization method proposed in the present paper aims at the generation of reduced reaction mechanisms for hydrocarbon combustion which are able to capture ignition delay times over broad validity ranges, defined by engine-relevant conditions.
Analysis of fast and slow dynamics of chemical kinetics using singular value decomposition
2023, Reaction Kinetics, Mechanisms and CatalysisVARIATIONAL PROBLEMS FOR COMBUSTION THEORY EQUATIONS
2022, Journal of Applied Mechanics and Technical Physics