Physically meaningful Lyapunov modal contributions in linear systems

This paper introduces a finite-time generalization of Lyapunov modal contributions (LMC) – a technique for energy-based selective modal analysis. The main focus is a case study involving a system with physically meaningful squared output and considering its unforced movement as transferring excessive energy towards its environment. This problem statement leads to LMCs estimating energy transfers measured in joules. The results facilitate solutions to several issues that occur in applications of methods based on Gramian spectral decomposition.


Introduction
Modal analysis and other spectral methods in control is an actively developing field of research. The classic modal analysis involves studying a system's eigenmodes, i.e. eigenvalues of its dynamics matrix. The introduction of participation factors (PF) allowed associating eigenvalues with sets of state variables resulting in the development of the selective modal analysis (Pérez-Arriaga et al., 1982). Such approach found a wide use primarily but not exclusively in the power industry for stability analysis (Song et al., 2019), optimal sensor and actuator placement (Singh et al., 2010), and many other applications. However, original PF definitions may lead to counter-intuitive or misleading results mainly caused by inconsistent dependencies on an initial condition and its perturbations (Hashlamoun et al., 2009;Konoval & Prytula, 2017). Therefore, a number of authors proposed alternative definitions (Hashlamoun et al., 2009;Iskakov, 2020;Konoval & Prytula, 2017), including energybased ones (Iskakov & Yadykin, 2021;Vassilyev et al., 2017) and ones generalized for nonlinear systems (Hamzi & Abed, 2020). Data-based approaches such as Prony analysis (Hauer et al., 1990) and others allow solving similar problems using real-world measurements instead of a pre-existing model. In addition, several techniques based on Koopman Mode Decomposition (KMD) allow connecting measured data with nonlinear models (Rowley et al., 2009). Furthermore, a vector version of Prony analysis is able to provide a finite approximation of KMD (Susuki & Mezić, 2015).
CONTACT D. E. Kataev dekataev@gmail.com Linear modal analysis and adjacent techniques found their applications and are mostly well-studied. Therefore, modern research focuses on fixing known flaws and generalizing existing approaches for nonlinear cases.
Sub-Gramian based techniques and methods, i.e. based on a spectral decomposition of a Gramian, a transfer function squared H 2 -norm, or a system output L 2norm are helpful for both lines of research. The base sub-Gramian method grants modal analysis with an ability to estimate modal interaction and its impact on system output (Yadykin, Kataev, et al., 2015). Its specific implementation also allows defining energy-based participation factors (Iskakov & Yadykin, 2021;Vassilyev et al., 2017). This approach provides additional means to ensure PF comparability and alleviate inconsistencies caused by initial conditions. The additive property of finite-time sub-Gramians (Kataev & Yadykin, 2016) hints at possible application to time-variant systems. A bilinear version of the sub-Gramian method (Iskakov & Yadykin, 2020;Yadykin & Iskakov, 2019) is a more developed way to generalize it for a subset of mildly nonlinear systems.
The main goal of this paper is to fill in the blanks left by existing sub-Gramian related research. A general understanding of sub-Gramians corresponding to system modes is obtainable by analysing the PFs of these modes. However, pairwise sub-Gramians and derived energy functionals correspond to mode pairs that lack defined PFs. Therefore, the published research treated pairwise sub-Gramians as estimations of modal interaction or a computational necessity that does not require detailed consideration. Industry representatives showed interest in a possible physical interpretation of sub-Gramian based functionals on several occasions, yet information on this matter is scarce. Scalar sub-Gramian based functionals share a problem of negative and complex values interpretation with PFs. It is especially significant in consequence of these functionals representing energy or energy-related quantities. Present lack of methods to thoroughly validate a software implementation of sub-Gramian based techniques is an obstacle for algorithm optimization and potential industrial deployment. All these shortcomings negatively impact the research process and limit currently feasible applications.

Preliminaries
Consider a linear time-invariant system with one outpuṫ where x ∈ R n is a system state vector, A ∈ R n×n is a system matrix, x 0 ∈ R n is an initial condition vector, n is a system order, y ∈ R is a system output, and C ∈ R n is an output matrix. Assume matrix A to be diagonalizable and to have n distinct eigenvalues {λ 1 , λ 2 , . . . , λ n }: where U and V are matrices composed of right and left eigenvectors, (·) T is the transpose operation. Eigenvectors u i and v i can be real-valued if a corresponding eigenvalue is real, but we operate under the assumption that they are complex-valued. We also impose the following conditions: where I is an identity matrix of order n. Let us consider a positive-definite solution P = P * > 0 of the following Lyapunov equation: where (·) * denotes conjugate transposition. A special case of P when Q = C T C is called an observability Gramian of the system (1). It allows expressing a squared L 2 -norm of the system output caused by the initial condition: Assuming A has a simple spectrum as specified earlier: where R i = u i v T i and R j are residue matrices. This equation defines a general spectral decomposition of a Gramian and allows us to introduce the following more specific definitions (Iskakov et al., 2018;Iskakov & Yadykin, 2021).
Definition 2.1: Observability sub-Gramians P i and pairwise sub-Gramians P ij are elements of an observability Gramian P spectral decomposition (Iskakov et al., 2018) Each sub-Gramian P i corresponds to an eigenvalue λ i of A and each pairwise sub-Gramian P ij to an eigenvalue pair λ i + λ j . Definition 2.2: Lyapunov modal contribution (LMC) E i characterizes the contribution of the ith mode to the squared L 2 -norm of the system output and is defined as Iskakov and Yadykin (2021).
Consequentially we also define pairwise LMC E ij : A finite observability Gramian P(t) is a solution of the following differential Lyapunov equation: According to Yadykin, Grobovoy, et al. (2015) and Kataev and Yadykin (2016), this solution has its own spectral decomposition, mostly similar to (3) but includes timevariant exponential parts: Definition 2.3: Finite observability sub-Gramians P i (t) and pairwise sub-Gramians P ij (t) are elements of a finite observability Gramian P(t) spectral decomposition (based on Yadykin, Grobovoy, et al., 2015) Naturally, we would suggest the existence of finitetime Lyapunov modal contributions E i (t) characterizing the contribution of the ith mode to the integral of the squared system output. We would presumably define such LMC as well as pairwise LMC E ij (t) in the following way:

Numerical validation and imaginary parts of sub-Gramians
While one can safely assume that analytically derived equations like (3) are correct, the problem of validating their software implementations posed some challenges. One way of numerical validation involves substituting a sum of all sub-Gramians into a corresponding Lyapunov equation. Another is to compare a sum of LMCs against a numerically calculated squared L 2 -norm of a system output. Such validation procedures do not guarantee that individual decomposition elements have intended values. Using finite sub-Gramians allows us to propose a method capable of validating each decomposition element.
Consider the system output as a sum of exponentials: The output y(t) is real-valued, but the exponentials y i (t) can be complex-valued, therefore: Integrating (10) leads to expressions identical to (7)- (8): It allows us to state the sub-Gramian validation problem as a curve fitting problem where m is an iterator, M is a number of elements in calculated time series, t is a time step for the calculated output and finite sub-Gramians.
Problem statement (13) works if all eigenvalues λ i are real; otherwise, fitting a real-valued signal with complexconjugate curves yields an infinite number of solutions. The most reasonable workaround lies in expanding (13) and eliminating its imaginary part because it would be mutually destroyed by its complex-conjugate pair during the summation. This leads us to the following problem statement: Corollary 3.1: Real and imaginary parts of a pairwise LMC are amplitudes of cosine and sine components of a corresponding finite pairwise LMC: Re(P q e qt ) = Re(P q ) e Re(qt) cos(Im(qt)) − Im(P q ) e Re(qt) sin(Im(qt)), Corollary 3.2: If y i,j (t) is oscillatory, then its phase shift for given x 0 is where zero shift corresponds to a cosine.

Pairwise frequencies and sub-Gramians
The squared output representation (10-11) shows that elements y i,j (t) are mutually amplitude-modulated oscillations. These are known to produce difference frequencies like (Im(λ i ) ± Im(λ j )) in: However, we propose a more general and systemagnostic term for such frequencies -a pairwise frequency. It prevents potentially misleading associations from appearing; simultaneously, it aligns with concepts of pairwise sub-Gramians and eigenvalue pairs as elements of a pairwise spectrum.
In order to show a physical side of the phenomenon, let us consider a linear quarter-car suspension model defined as: where z u is an unsprung mass (UM) coordinate, v u -UM velocity, M u -UM, k b -suspension stiffness, z b -a vehicle body coordinate, k u -tyre stiffness, c b -suspension damping, v b -body velocity, c u -tyre damping and M b -body mass. Providing realistic parameter values would result in system (16a) having a matrix A with spectrum (λ u , λ * u , λ b , λ * b ). λ u and λ b represent eigenmodes of the unsprung mass and body respectively. Let's assume they have frequencies of 14 and 1 Hz.
According to (9a), any observed movement y(t) caused by x 0 would be defined by 1 and 14 Hz oscillations. However, decomposing y 2 (t) with sub-Gramians means considering the system inside an 'energy space' instead. More specifically, instead of movement-produced signals, we observe their signal energies. Picking the proper output matrix or matrices C may provide observations of physically meaningful energy (J) or energy-based values, e.g. power (W) or action (J·s). System (16a) physically contains two kinetic energies of a body E K b and an UM E K u as well as two potential ones of a tyre E P u and a spring E P s . These signal energies also have a new spectrum of pairwise frequencies: 2, 13, 15, and 28 Hz. Some energies are straightforwardly dependent on a single state variable: It is relatively safe to assume these energies inherit double frequencies of respective state variables as predominant ones: 2 Hz for E K b and 28 Hz for E K u and E P u . However, the spring compression does not explicitly exist as a state variable, and its potential energy is expressed as follows: Considering the summand 2z b z u as a function of time and expanding it as a product of two sums of exponentials (9a) results in dominance of mutually amplitude-modulated 1 and 14 Hz oscillations or an equivalent sum of 14 ± 1 Hz ones. As a result, it allows us to state that E P s is the only physical energy in the system tangibly but not exclusively defined by pairwise 14 ± 1 Hz frequencies.
The preceding research papers on sub-Gramians present vaguely defined pairwise sub-Gramians and Gramian-based energy functionals as numerical estimations of modal interaction. This example shows that it is possible and presumably productive to give them detailed application-specific definitions. The particular case above suggests modal interaction in mechanical systems to reflect the process of multiple bodies interchanging energy using a common intermediate energy container.

Case study
Suppose a stable system (1) has several outputs. For all of them, y 2 (t) quantifies mechanical power this system applies to its environment. Then x T 0 P(t)x 0 is a time integral of power and would quantify work, i.e. energy leaving the system. Providing that these outputs cover all such energy transfers, the energy conservation law must be in effect. Suppose this system has initial condition x 0 and contains an excess of energy E sys 0 (x 0 ). Then the system does work x T 0 P(t)x 0 in the process of unforced movement: We expect the system studied this way to provide detailed information on energy propagation in it.

Linear road vehicle suspension model
We use a linear vehicle suspension model with four wheels, a rigid body, tyre deformation, and no boundary conditions. Figure 1 provides a body diagram of the model. It is four quarter-car models (16a) linked together except the vehicle body now has three degrees of freedom: where F fr , F fl , Frr, Frl are resulting forces of front right (FR), front left (FL), rear right (RR), rear left (RL) suspensions (i.e. springs and dampers but not tyres), I y -body moment of inertia about the y-axis, θ p -pitch angle, L f -distance from the centre of gravity (COG) to the front axle, L r -COG to rear axle distance, I x -body moment of inertia about the x-axis, θ r -roll angle, L fr , L rr , L fl , L rl -lateral distances from COG to respective suspension attachment points. Equations (16a) need adjustments to account for the new degrees of freedom. Table 1 contains all referred parameter values. A linear state space representation of this model is a 14th-order system. The chosen initial condition x 0 means a 0.01 m FL unsprung mass (UM) displacement. Such x 0 value provides the system with 15.4 J of excess energy E sys 0 . The tyre stores 14 J and the spring another 1.4 J. This energy can leave the system only through its eight damping elements: four tyres and four dampers. For each, we introduce a separate output. The values of respective matrices C provide y 2 (t) quantifying corresponding damping force power measured in watts. It is possible because for all damping in this system: where N(t) is power of a damping force applied to a body, F(t) -this damping force, c -damping of an element, v(t) -relative velocity of a body, C -a corresponding output matrix. Appendix contains a state space system and clarifies its state vector structure. The authors confirm that the data supporting the findings of this study are available within the article and its supplementary materials. Table 2 presents the spectrum of the system involved.

Results and discussion
Total work of all eight damping elements calculated as x T 0 P(∞)x 0 equals the initial energy value E sys 0 . Example (16a) showed that we could split the pairwise frequencies of the system into three groups. The first one contains 27.6-28.6 Hz UM oscillations. These are weakly interconnected and primarily represent UM kinetic energy variation. The second group contains 2.1-2.7 Hz body oscillations. Same-mode pairs represent variations of the body's three kinetic energies: one translational and two rotational. System equations suggest that intermodal pairs occur due to suspension attachment points states depending on all three degrees of freedom. Combining pairwise LMCs by energy frequencies instead of natural state space ones provides a different perspective while using the same method, e.g. explicitly accounting for inside processes that may be not fully observable in the state space or direct measurements. The third group of frequencies contains 12.4-15.7 Hz oscillations produced by the interaction of a body and UMs. The role of these modal interactions is to partially define energy variations in suspension springs. While ∼ 14 Hz components occur in every possible squared output of the system, they are significant enough to be visually distinguishable only in squared outputs explicitly containing spring energy (see Figure 2). Table 3 presents a partial spectral decomposition of system work for these three frequency groups. Stating that LMCs E ij are elementary energy transfers forming total system work on its environment does not align well with the existence of negative-valued LMCs. For this reason, we investigate interior energy transfers, i.e. works of elastic forces, both positive and negative. More specifically, we consider only LF tyre and LF suspension doing over 99% of total work. Figure 3 demonstrates that total elastic forces work converges to the sum of all LMCs E ij since both follow the energy conservation law. An absolute amount of internal energy transfers involving these elastic forces is defined for each force-body pair as: whereĒ is an absolute amount of transferred energy in a given force-body pair, v -absolute velocity of a relevant body, and k -stiffness of a relevant spring. Figure 3 shows that the sum of all positive LMCs (E ij > 0) converges to the total internal energy transfer value. It allows us to preliminarily conclude that positive-valued LMCs meeting conditions presented at the start of this section do reflect internal energy transfers. Therefore, the sum of all negative LMCs being a difference between absolute energy transfers and done work should account for all energy transferred more than once. That is a good enough explanation when dealing with infinite time intervals. However, this research considers evolutionary equations, so it is helpful to be more specific. According to Figure 3, absolute and physically correct works of elastic forces start to diverge only after 21 ms while negative LMCs are non-zero at any moment except zero. It concludes that finite-time negative LMCs account for transfers of energy resulting in that energy needing further transfers to leave the system.

Conclusion
This paper introduces a finite-time generalization of Lyapunov modal contribution (LMC) presented in (Iskakov & Yadykin, 2021). Such generalization allowed to create an independent validation method for software implementation of sub-Gramian based energy functionals. The method involves linear fitting of a squared system output with exponentials, demonstrating a certain relation to Prony analysis. Consequently, an adequately modified vector Prony analysis may be a way to define LMC for time-variant systems. We clarified the potential roles of pairwise sub-Gramians, LMCs, and corresponding frequencies in a linear system. These roles are often able to be more specific than generalized modal interaction. In certain systems, its cause is multiple forces interacting with a common intermediate energy container. To put it another way, state derivatives depending on values that do not explicitly exist in a state vector play a crucial role in defining modal interaction.
The main novelty of this paper is its case study involving a system with physically meaningful squared output and its unforced movement considered as transferring excessive energy towards its environment. Specifically, it results in LMCs having a meaning of elementary mechanical work or energy transfer measured in joules. Internal energy transfers in the studied system are reversible, while external ones are one-way. It provided some insight on negative sub-Gramian based energy estimations. The particular case study allows us to conclude the following: • The sum of all LMCs yields an estimation of physically correct external work, i.e. irreversible external energy transfer; • The sum of all non-negative LMCs estimates unsigned total internal work, i.e. amount of internally transferred energy; • Negative LMCs estimate an amount of energy transferred so that it would require further transfers, i.e. a fraction of internal work that the system would eventually undo.
This work may allow new energy-based problem statements for tuning and controlling of many kinds of mechanical suspension systems. It may help relevant applications to formally address certain poorly formalized issues related to human-machine interaction and mechanical reliability.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Appendix. Linear system
An approximate value of A is presented below. All unspecified elements are zeroes. The first matrix block with excluded first 4 rows is transposed to better fit the page. The state vector: FL UM z, pitch angle (rear to front), roll angle (right to left), body COG z, pitch velocity, roll velocity, body COG vertical velocity, FL UM velocity, FR UM z, FR UM velocity, RR UM z, RR UM velocity, RL UM z, RL UM velocity.   An approximate spectrum of original A for reference: 32.8 ± 86.8i; 32.7 ± 78.2i; 33.1 ± 89.8i; 33.4 ± 89.0i; 1.40 ± 6.66i; 2.42 ± 8.62i; 1.87 ± 7.65i.
Matrices C for y 2 (t) to quantify power (all unspecified elements are zero):