#### Abstract

The large-scale structure systems in engineering are complex, high dimensional, and variety of physical mechanism couplings; it will be difficult to analyze the dynamic behaviors of complex systems quickly and optimize system parameters. Model order reduction (MOR) is an efficient way to address those problems and widely applied in the engineering areas. This paper focuses on the model order reduction of high-dimensional complex systems and reviews basic theories, well-posedness, and limitations of common methods of the model order reduction using the following methods: center manifold, Lyapunov–Schmidt (L-S), Galerkin, modal synthesis, and proper orthogonal decomposition (POD) methods. The POD is a powerful and effective model order reduction method, which aims at obtaining the most important components of a high-dimensional complex system by using a few proper orthogonal modes, and it is widely studied and applied by a large number of researchers in the past few decades. In this paper, the POD method is introduced in detail and the main characteristics and the existing problems of this method are also discussed. POD is classified into two categories in terms of the sampling and the parameter robustness, and the research progresses in the recent years are presented to the domestic researchers for the study and application. Finally, the outlooks of model order reduction of high-dimensional complex systems are provided for future work.

#### 1. Introduction

The large and complex structures exist widely in engineering field of aviation, aerospace, shipping, and so on which are complex, high degrees-of-freedom (DOFs), and coupled with a variety of physical mechanisms. Dynamical systems are the basic framework for modeling and control of these enormous varieties of complex structure systems [1]. Examples include fluid dynamics, design optimization, control, chemically reacting flows, data-driven systems, and vibration suppression in large structure systems and other complex underlying physical process. The mechanism model of any complex structure system can be established by classic mechanics in theory. However, the model is usually a large-scale partial differential system, an approximate simplified high-dimensional ordinary system, a coupling system with partial system, and ordinary system, which cannot be solved by theory directly. The common method for dealing with the abovementioned large complex system in engineering is to use the finite element, finite difference, finite volume, and other methods for numerical simulation analysis [2]. However, the number of DOFs of the complex system obtained by finite element methods may be tens of thousands. If the system has strong fluid-solid coupling effects [3], the number of DOFs may reach millions, even billions. Numerical simulation of large-scale dynamical systems plays a fundamental role in studying a wide range of complex physical phenomena. Although the computer is more advanced in software and hardware, however, the inherent large-scale nature of the models often leads to unmanageable demands on computational resources and a large amount of calculation time to obtain the accurate solution, such as several hours or even longer [3–6].

The MOR methods were developed in the area of systems and control theory. However, there are several definitions of MOR, and it depends on the context in which one is preferred. A flow chart (Figure 1) for modeling complex physical systems can be elaborated based on time domain, frequency domain, time-frequency domain and optimization techniques, and artificial intelligence as follows [7–12].

In time domain, Chebyshev-based approach for model order reduction of linear systems is presented based on Chebyshev rational functions [7]. Algorithms for the estimation of the moments matching of linear and nonlinear systems are proposed for model order reduction, and the estimates are exploited to construct families of reduced-order models [8]. Wavelet-based approach is proposed for model order reduction of linear circuits in the time-frequency domain [9]. Approaches for model order reduction based on artificial neural network aims at obtaining a reduced-order model out of a relatively complex model, generally obtained in a reasonable time and has accepted error [10]. A low-cost fuzzy rule-based implementation of Sammon’s method for structure preserving model order reduction is presented [11]. Particle swarm optimization is usually used to solve optimization problems when the number of parameters is low, and also to find a good solution typically involves multiple evaluations of the objective function [12].

A series of model order reduction methods were proposed to reduce the number of DOFs of the system to improve the efficiency of calculation in the field of science and engineering, for example, center manifold method, Lyapunov–Schmidt (L-S) method, Galerkin method, nonlinear Galerkin method, mode synthesize method, Krylov approximation method, balanced truncation method, and POD method [3–24]. Thus, the area of model reduction contains a broad set of mathematical methods to yield and evaluate the reduced models. These model order reduction methods have been applied in various engineering fields and become more mature, but each method has its adaptability and limitations. Many order reduction methods may fail when the complexity and DOF of the system increase, state parameters vary, and nonlinear factors couple with each other. Then, the order reduction model cannot reflect the real dynamic characteristics of original systems. In this paper, the mature model reduction method in high-dimensional systems will be reviewed, and basic principles, adaptability, and limitations will be expounded.

The POD method is an order reduction method suitable to process large and complex structures quickly and efficiently [2, 23, 24]. The POD method has been widely studied and applied by many scholars in recent years. The method can obtain the main structural components of complex systems with few POD reduced-order modes. The method can reduce the DOFs number of complex systems greatly and improve the computational efficiency significantly while ensuring the accuracy of reduction. The method is important for dynamic analysis and optimized design of parameters in large and complex systems.

The motivation of this paper is to summarize the review of the order reduction methods in large-scale structure systems and provides the classification of POD. In Section 2, the center manifold is introduced. The L-S reduction and Galerkin methods are discussed in Sections 3 and 4, respectively. The mode synthesis method is introduced in Section 5. Several issues of the POD method which are worthy to study in the future will be pointed out based on the characteristics of each order reduction method in Section 6. Finally, the conclusions and outlooks are drawn in Section 7.

#### 2. Center Manifold Order Reduction Method

The center manifold order reduction method is a local order reduction method based on the center manifold theory [13, 23–27]. The method projects the high-dimensional system onto low-dimensional central subspace, according to the differential homeomorphic mapping between central subspace and stable subspace near the equilibrium point, and reserves all the topological properties of the high-dimensional system near the equilibrium point. So, the asymptotic behavior in neighborhood of local bifurcation point and central subspace of the original dynamic system is equal in the high-dimensional nonlinear system [25]. The basic principle of the center manifold order reduction method will be described briefly as follows.

The dynamic equation of the high-dimensional complex nonlinear system in state space is expressed as follows:where ** a** is a system parameter,

*m*and

*n*are the number of DOFs in the parameter space and state space, and , where is

*r*order continuous differentiable function space. Formula (1) is expanded linearly at the equilibrium [23, 24].where is the Jacobian matrix of and is the nonlinear part. Assuming that the equilibrium point is nonhyperbolic,

**has eigenvalues with the real part of 0. For ease of understanding, assuming that the real part of the remaining eigenvalues is less than 0, formula (2) can be expressed via the coordinate transformation of eigen-space as follows [25]:where ,**

*A***is a matrix of , and**

*B***is a matrix of . Nonlinear function is at least second-order differentiable, so according to center manifold theory [13, 23–25], there is a differential homeomorphism mapping in partial neighborhood of so that formulas (3) and (4) can be expressed as follows:**

*C*Substituting formula (5) into (3) and (4), we obtain the following:

can be calculated through formula (6). However, it is difficult to obtain the analytical expression directly through formula (6) in the complex system [23, 24].

Formula (7) in neighborhood of equilibrium point is obtained by differential homeomorphism mapping. Formula (7) and the original system have completely topological equal dynamic behaviors.

The theory of the center manifold order reduction method is rigorous and complete, and the method has simple steps to solve the order reduction model. So, the method is applied widely in high-dimensional complex nonlinear dynamic systems. Cao et al. [25] applied the center manifold theory to study the reduction of Fold bifurcation in three points electrical system voltage collapse. Yoursef and Yoursef [26] and Zhang et al. [27] proposed general program to calculate arbitrary high-dimensional systems by Mathematica, and center manifold dimension is less than 6.

Boyaci et al. [28] established a perfectly balanced symmetric rotor which is supported by two identical floating ring bearings, and the analytical results obtained from the center manifold reduction are compared with numerical results by a continuation method. The random center manifold theory was established to lay the foundation for the theory analysis of the high-dimensional nonlinear random dynamic system in Ref. [29]. Liu et al. [30] used center manifold theory to reduce the nonlinear aeroelastic system with 9 state variables to 2-DOF at the equilibrium. Rendall [31] analyzed some dynamic systems of the spatial isotropic universe model by center manifold theory and obtained progressive behavior of all solution of Einstein equation with Bianchi-III formation.

Sinou et al. [32] applied center manifold theory to analyze the effects of friction of the high-dimensional nonlinear braking model to system stability. Center manifold theory was applied by Hupkees and Lunel [33] to analyze the solution behavior of nonlinear autonomous mixed functional differential equation near the equilibrium. In recent years, Kano et al. [34] investigated the order reduction and bifurcation analysis of a flexible rotor system supported by a full circular journal bearing, and bifurcation phenomena at around the instability point are investigated by applying the center manifold theory and using the normal form theory. Liu et al. [35] investigated the stability and bifurcation behavior of a kind of active magnetic bearing rotor. It is found that a Hopf bifurcation occurs in the system by using center manifold and normal form. Liu et al. [36] established a Jeffcott rotor model for the rotor system of the permanent magnet synchronous motors in electric vehicles, and center manifold theorem and Lyapunov method are used to determine the stabilities of multiple equilibrium points.

The center manifold order reduction method is applied in many fields, but the theory belongs to local redaction method. It is only established near the equilibrium and very difficult to discuss the global behavior of high-dimensional systems. At the same time, nonconvergence problems also occur when using series solutions [23, 24]. In addition, the DOF of the reduced-order system is determined by center manifold. For a nonlinear system with higher DOF, the center manifold dimensions may be still high and it is difficult to handle the method. So, the center manifold order reduction method has significant theoretical values for the complex nonlinear system below a dozen DOF.

#### 3. L-S Method

The Lyapunov–Schmidt (L-S) method [3] is similar to the center method, and it is also a local reduction method.

The difference can be expressed as follows. Center manifold theory reduces the high-dimensional complex nonlinear system to a low-dimensional dynamic system which can remain equivalent to the dynamics topology of the original system. However, the L-S method can process static bifurcation of stationary solution directly in the high-dimensional complex nonlinear system. From the implicit function theorem, one set of equations has unique solution, the equations are solved and the solution is substituted to the other set of equations, and then the solution of the high-dimensional static bifurcation issue can be reduced to a solution of low-dimensional equations.

The basic principle of the L-S method is to project high-dimensional nonlinear algebra equations to the two mutually orthogonal subspaces of its value space, and two sets of equations can be obtained [13, 24]. The basic principle of the L-S method will be briefly described as follows.

Consider the static bifurcation of the stationary solution of the high-dimensional nonlinear system; from the mapping of equation (1), the high-dimensional nonlinear algebraic equation can be expressed as follows:

*X* and *Y* are Banach space. The static bifurcation point of the original system can be changed to the coordinate origin by coordinate translation transformation, that is, . Assuming that is Fredholm operator, the zero space is and the value space is . And the dimension of the zero space is finite and greater than zero. Then, through the spatial straight sum decomposition near the coordinate origin, the following is obtained [37]:

Defining the mapping operator and the complement operator , I is identical operator. The state variable can be expressed as, based on equation (9). Equation (8) is equivalent to equations (10) and (11) from action of mapping operator and complement operator [23]:

Because , if ** A** is limited to

*M*, then

*A*is reversible. It can be known from the implicit function theory that the unique solution of equation (10) in the neighborhood of makes and . Submit the solution of the above equation to equation (11), there is a mapping , so that

The solution of equation (12) corresponds to the solution of equation (8), thereby the static bifurcation issue near the equilibrium in the high-dimensional complex system can be equivalent to the solution of the low-dimensional nonlinear algebraic equation for the reduction.

The L-S method has strict theoretical basis, and it is used widely to study the bifurcation theory in the high-dimensional complex nonlinear system. Chicone [37] used Melnikov integral and the L-S method to analyze the periodical solution bifurcation issue of the autonomic differential system with small parameters. The L-S numerical iterative method was proposed to solve the reduced-order bifurcation equation and Hopf bifurcation issue in Ref. [38]. Zhang and Stewart [39] applied the L-S method to discuss the existing of the bounded solution of the nonautonomous parabolic equation. Lukas [40] used the Wolfram symbol software package of Mathematica to develop the L-S reduction algorithm, and the algorithm can be applied to calculate the three-dimensional vortex structure of the dispersed nonlinear Schrodinger equation. Sandfry andHall [41] applied the L-S reduction to determine an analytic relationship between parameters to recognize conditions for which a jump phenomenon occurs. In recent year, Buica et al. [42] applied the L-S method to analyze the periodical solution bifurcation issue of the regular nondegenerate cluster in the Lipschitz system. The L-S method was applied by Chen and Zheng [43] to study the solution aggregation behavior of the fractional order nonlinear Schrodinger equation. Pogan et al. [44] studied the O(2) Hopf bifurcation of the viscous shock in tube, and they used L-S method to handle the quasi-linear hyperbolic-parabolic system with hot viscoelastic issue. Cao et al. [45] studied the dynamic property of the damped harmonic oscillator with delayed feedback, and they applied the L-S method to obtain the judging condition for bifurcation of periodical solution and number of branches. Li and Ma [46] promoted the L-S method to the fractional differential system. Guo and Ma [47] analyzed the stability and bifurcation behaviors of the reaction-diffusion equations with the Dirichlet boundary condition and applied the L-S method to prove the existence of the partial nonuniform solution.

The L-S method is complete in theory, and the method can process the bifurcation issue of the high-dimensional nonlinear system effectively, but it is difficult for the method to process the static bifurcation issue of the large and complex nonlinear system. The main reason is as follows. It is difficult to solve the low-dimensional nonlinear algebraic equations in the processing of the high-dimensional space decomposing. The dimension of the obtained low-dimensional nonlinear algebraic equations is still very high after the high-dimensional complex system is decomposed by the L-S method. All the literature about the L-S method known by the authors of this paper can process the DOF of the high-dimensional system which is less than ten, so the L-S method is suitable for studying the bifurcation of the high-dimensional complex system with DOF less than twenty.

#### 4. Galerkin Method

The first two-order reduction methods are mainly used to analyze few high-dimensional nonlinear system models with several DOFs theoretically, but the complex dynamic system in engineering is usually continuous partial differential system, simplified approximate high-dimensional ordinary differential system, or coupling system with both partial and ordinary differential. Galerkin mapping is a bridge from partial differential to ordinary differential, and it is an effective order reduction method to handle this system [2, 15–17]. The final purpose of the Galerkin mapping is to obtain a low-dimensional dynamic system, which can reflect the dynamic characteristics of the original system. The basic principle of the standard Galerkin method is to map the original system to the mode space and make up the system by intercepting low-order modes and ignoring the effects of the high-order mode to achieve the purpose of reducing the order [15, 17]. The basic principle of the Galerkin method will be described briefly as follows.

The dispersed dynamic equations of the high-dimensional complex system obtained via the finite element or other method are

**M, C,** and **K** are total mass matrix, total damping matrix, and total stiffness matrix, respectively, and is the nonlinear force. Through the mode coordinate translation, , where is the system mode matrix ( and are the matrices constituted by the first *m* and the remaining *l* eigenvectors, respectively) and represents coordinates. Equations (14) and (15) can be obtained by projecting equation (13) into mode space as follows:

The standard Galerkin method directly ignores the effects of the high-order mode, so the original system can be reduced to the following [15]:

When the truncation mode *m* is taken to a certain extent, the original system can be expressed approximately as follows: .

Because of neglecting the effects of the high-order mode, the reduced model of complex nonlinear systems has large errors. So, some scholars proposed the nonlinear Galerkin method, which finds the approximate relationship between the high and the low mode coordinate by constructing inertial manifold [15–17, 24, 49] based on the inertial manifold theory [48]. The nonlinear Galerkin method can improve the approximate accuracy, but the method needs more calculation time [48] because each integral includes the nonlinear factor of the approximate inertial manifold method. Garcia-Archila [50, 51] proposed the posterior Galerkin method to solve the problem, which can save the calculation time [17]. For this method, it is not necessary to calculate the approximate inertia manifold in each step of the integration step and only need to consider the effects of higher order modes in the output of each step.

After decades of the development and application, the Galerkin method has become an important way to obtain the numerical solution of large and complex systems rapidly, and its related application literatures are endless. For example, Wang and Cao [15, 16] proposed the predictive correction Galerkin method, and they used the method to reduce the order of the large rotor-bearing system. Ding and Zhang [52] analyzed an isotropic flexible shaft acted by nonlinear fluid-induced forces and reduced dimensions of the rotor system by both the standard Galerkin method and the nonlinear Galerkin method. An adaptive discontinuous Galerkin method was proposed to obtain a high-resolution numerical solution efficiently in Ref. [53]. Sembera and Benes [54] applied the nonlinear Galerkin method to analyze the reaction-diffusion system in bounded invariant domain and proved the convergence of the method. Chatzisavvas et al. [55] investigated the effect of hydrodynamic thrust bearings on the nonlinear vibrations and the bifurcations occurring in rotor-bearing systems by using a global Galerkin approach. The postprocessing Galerkin method was used to reduce the mode order for nonlinear wind turbines [56]. Boelens et al. [57] used the second-order accurate discontinuous Galerkin finite element method to discretize the governing equations on a hexahedral mesh. Amabili et al. and Sarkar et al. [58, 59] applied Donnell nonlinear cylindrical shell theory to deeply research the great vibration for liquid filled cylindrical shell, and they obtained the 16-DOF accurate response of the system on simple hormonic excitation by using the Galerkin mode truncation method. The Galerkin and the partial POD methods were combined to study two one-dimensional parabolic equations [60]. Xie et al. and Xie et al. [61, 62] applied the Von Karman large deflection plate theory, Galerkin method, and POD method to study the wing aerodynamic flutter problem with supersonic flow condition.

The Galerkin method has been widely applied in engineering, but the method is insufficient, for example, the existence of inertial manifold in high-dimensional complex system, structure issue of inertial manifold, the truncation principle of modal order, and the convergence of the algorithm [23, 24, 63].

#### 5. Modal Synthesis Method

The modal synthesis method is another effective order reduction method; it is essentially a classic Galerkin method, so the method is applied widely to reduce the order of the large and complex system [15–17, 64–78]. The principle of the modal synthesis method is dividing large and complex system into substructures, ignoring the high-mode of sub-structures, and then synthesizing the system composed of low-order modalities of each substructure, thereby achieving the purpose of reducing system DOF [17, 19, 70, 73]. According to the different connections of substructures, the method can be divided into fixed-interface modal synthesis method, free interface modal synthesis method, and mixed interface modal synthesis method [73, 75]. The basic principle of the method will be introduced briefly as follows.

The discrete dynamic equation of the large and complex system is as follows:

In equation (17), , , , , and are the mass, damping, stiffness matrix, point displacement, and nonlinear force, respectively, which are composed of *n* substructure systems. are the mass, damping, and stiffness matrix of the *i*-th substructure, respectively. of substructure are divided into internal coordinate and boundary coordinate, and the dynamic equation of the substructure can be written as follows:where , , , , and .

To calculate constrained main mode and constrained mode with boundary or interface constraint, are the number of coordinates of the substructure internal nodes, the number of boundary coordinates, and the number of retained constraint main modes [1, 72]. can be calculated by the constrained boundary, and then the translation relationship between the physical coordinate and the modal coordinate can be obtained [1, 75].

The structural dynamic reduction modal can be obtained via substituting equation (19) to equation (18):where , , , and . The total reduction modal of the high-dimensional complex structural system can be obtained by integrating the reduced modal of each substructure.where , , , , and .

The modal synthesis method has simple principles and good robustness, and the method is easy to operate. So, the method is widely used to reduce the order of the engineering complex structure system. Yang [1] used finite element software to establish discrete dynamic modal of five points reverse rotation dual-rotor system, and they divided dual-rotor into two substructures of inner and outer rotors. They reduced the order of the system by the fixed-interface mode synthesis method and researched the moment response in processing of accelerating of five points reverse rotation dual-rotor system, and then the calculation scale is reduced. Glasgow [64] calculated the critical speed and the modal shape of the dual-rotor-bearing system by using the substructural modal synthesis method, and the calculation scale is reduced greatly, and then they discussed the calculation accuracy and the error of the method in detail. The mode synthesis method was used to reduce the linear part of the rotor-bearing system modally, and as the boundary coordinate, the nonlinear part together with the reduced modal constitutes a reduced-order system [66]. Wang and Kirkhope [67] improved the free interface mode synthesis method and applied it to reduce the order of multirotor-bearing system. Sundararajan and Noah [68] used the modal synthesis method to reduce the order of the nonlinear rotor system and analyzed the stability and the bifurcation of the reduced modal. Iwatsubo et al. [69] proposed an effective procedure using the component mode synthesis and the method of multiple scales or the harmonic balance method for the nonlinear vibration analysis of rotor systems. The mixed interface modal synthesis method was used to reduce the order for the rotary symmetry impeller structure [71]. Shanmugam and Padmanabhan [72] proposed the mixed interface modal synthesis method, which is more suitable to the dynamic analysis of the two rotors system. In recent years, Beck et al. [73] applied the modal synthesis method to reduce the order of integrally bladed rotors and analyzed the dynamic characteristics of the reduced modal. Zheng et al. [74] proposed a generalized and efficient method for parametric response analysis of large-scale asymmetric rotor in order to avoid costly approach, and the fixed-interface component mode synthesis is employed to form a reduced-order model. Zhang [75] combined the Ansys software and the mixed constraint modal synthesis method to study the large civil structure with nonlinear part. Luo et al. [76] used finite element software and the free interface modal synthesis method to establish the dynamic modal of the high-dimensional two rotors system with collision friction failure. The nonlinear factors of intermediate bearing and the squeeze film damper are taken into account in the modal. They used unit impulse response and the Duhamel integral method to obtain the numerical solution of the equation, and they studied the collision friction response characteristics of the reverse rotary two rotors system. The assessment method was studied by Kim et al. [77] to select the modal in the synthesis method. Joannin et al. [78] combined the complex mode and the modal synthesis method to solve the steady state response in the unconservative system.

The modal synthesis method is suitable for dealing with the large and complex structure, but the method mainly performs modal reduction for linear substructure and ignores the effect of high-order mode and system partial mode, so the reduction modal may have great error. The method is not accurate [23, 24], and the effect of reducing the order is not obvious for the large deformation and strong coupling nonlinear structural system.

#### 6. Proper Orthogonal Decomposition (POD)

In this section, the basic principles will be introduced in Section 6.1. The classification of POD is provided in Sections 6.2 and 6.3 based on sampling and parameter adaption. In Section 6.4, some improved POD methods and related issues are discussed. The POD method is the specific recommendation method which will be widely applied in actual engineering. The POD method contains sample, parameter adaptation problems, and other improved methods, and the classifications are listed in Table 1. The POD methods based on sampling and parameter adaptation problems are introduced briefly as follows.

The construction of reduced-order models (ROMs) plays a key issue in POD methods. Different sampling signals can yield different ROMs via POD, so many researchers have focused on determining the sampling methods that will yield the optimal ROMs. ROMs can be obtained from the following types of response signals: chaotic response signals, random response signals, transient response signals, and signals obtained through combined sampling methods.

The responses of a complex dynamic system are closely related to the system parameters and initial conditions. The ROMs obtained via POD usually lack robustness against variations in the system parameters. In principle, new ROMs should be constructed for new parameters; however, the number of calculations required to calculate each possible parameter response of a complex system to construct the corresponding ROMs would be very large. To resolve the problem of model order reduction for a complex system over a broad parametric domain and ensure the robustness of the ROMs, many modified methods have been proposed by researchers, including global POD methods, local POD methods, and adaptive POD methods.

##### 6.1. Basic Principles, Advantages, and Disadvantages

The proper orthogonal decomposition (POD) method is also a powerful order reduction method to deal with the large and complex system. Other names of the method are Karhunen–Loève decomposition, principle component analysis (PCA), singular value decomposition (SVD), and Hotelling translation [79, 80]. In the middle of the twentieth century, these methods were proposed, respectively, by Kosambi [81], Karhunen [82], Lorenz [83], Pougachev [84], and Hotelling [85].

The POD method is a projection-based order reduction method, which is similar as the projection/Galerkin method which maps the high-dimensional system to a low-dimensional subspace. The difference between the two methods is as follows. POD is a statistics method, and the POD reduction modal or the reduction basic function is obtained by solving eigenvector of the autocorrelation matrix. The autocorrelation matrix is statistically significant 2nd center moment [79]. The basic principle of the POD method is that the autocorrelation matrix is constructed from numerical simulation snapshot signal or experimental data snapshot signal in the original system. The POD reduced modal or reduced basic function is obtained by solving the eigenvector of autocorrelation matrix and then the high-dimensional system is projected onto the subspace of the reduced-order modal corresponding to first maximum eigenvalues, thereby achieving the purpose of reducing the order. The basic principle will be briefly described as follows.

For the state variables of the high-dimensional complex system (*n*th dimension) , the discrete time series can be obtained from numerical simulation signal or experimental signal, and the time series is projected onto the space opened by complete orthogonal specification base , that is:

We hope to find a group of orthogonal specification base that satisfies the minimum value constraint under the square norm [23, 24]:

In equation (23), indicates the average operator during the sampling period, denotes mapping operator, represents the inner product on Hilbert space, and is the norm.

The orthogonal specification base satisfying the above conditions can be obtained by the Lagrange multiplier method. The objective function is defined as follows [3, 23, 24]:

In equation (24), is the Lagrange operator. The extreme value of formula (24) is found, and let , we can get the following [79, 80, 86]:

Formula (25) is second class Fredholm equation, and its kernel function is autocorrelation: . We can prove that the orthonormal basis which is also called POD order reduction modes which are the eigenvectors of the autocorrelation matrix based on equations (22) and (25) and normal orthogonal condition.

The high-dimensional state variables are projected onto the subspace of first few POD modes (eigenvectors are arranged in descending order of corresponding eigenvalues), thereby the best square approximation of original system state variables during the sampling period can be obtained with the least reduced-order mode number *m*, that is:

On the other hand, the signal matrix composed of discrete time series can be obtained by the numerical simulation snapshot signal and the experiment signal. The signal matrix can be decomposed by singular value decomposition (SVD) theory [87].

In equation (27), and are orthogonal matrices, which satisfy , diagonal matrix is a singular matrix of the sampling signal matrix, and . Comparing equations (22), (26), and (27), the matrix translation relationship can be obtained between the eigenvalue of autocorrelation matrix and singular value decomposition of signal matrix.

The coordinate translation relationship can be defined by formula (28). The high-dimensional complex system is projected onto subspaces of a few POD reduced-order modes; and , respectively, are reduced-order modal translation matrix and the POD modal coordinate. Substitute equation (28) into dynamic equation of the large and complex structure system (17), and then simplifying the formula, we can obtain the following [87]:

In equation (29), , , , and . A low-dimensional dynamic modal can be obtained. This is the traditional POD reduction method.

In many literatures [79, 87–89], the dimensional number of the reduced-order *m* is determined by defining the singular values or energy percentage of autocorrelation matrix eigenvalues.where is the eigenvalue of the autocorrelation matrix, and the eigenvalue is determined , sometimes [87, 89], for ensuring the reduced-order signal and the original sampled signal matrixnorm approximation.

The DOF of the reduction modal obtained by the POD method is very low, and DOF is the best square approximation of the original system, so the method has obvious advantages in reducing DOF of the high-dimensional system and improving computation efficiency. However, from the basic theory of the method, the disadvantages of the method can be obtained. The advantages and disadvantages of the method will be described as follows.

According to the basic theory of the POD method, the method is easy to operate. The most significant advantage is as follows. The most important components of the infinite dimensional or high-dimensional complex system can be obtained with very few POD reduced-order modes, and the method has the optimal approximation under square norm. The DOF of the original system is greatly reduced [79, 80, 87] by using the POD method. So, the POD method is widely used in area of engineering, for example, fluid dynamics [86, 90–95], signal and image processing [96–99], optimal design [100–103], ocean engineering [104], and structural dynamic mechanics [58, 59, 61–63, 105–109].

However, according to the basic theory of the POD method, we can easily find that the POD modes come from the eigenvectors of the autocorrelation matrix, which is constructed by numerical simulation signal or experiment signal. As a result, for different sampling data, the POD modes are different with different sampling parameters, method, or length, which significantly influence the reduced modal. And the actual constructed modes just are the optimal approximation [110] for the original system under recent sampling length, which is not the optimal approximation for the original system of all states, so the constructing of the POD reduced modes is key for the POD method. In response to these problems, scholars have proposed a number of improved POD methods. Through the author’s literature research, it is found that these methods can be divided into two categories: one to solve the sampling problem of the POD method and the other to solve the parameter autocorrelation problems of the POD method. The two types of POD methods will be, respectively, reviewed as follows.

##### 6.2. Sampling Study of POD Method

According to previous description of the POD method, the most important part is the construction of reduced-order mode. In view of the sampling point, the different sampling signal can obtain different POD reduced modes; therefore, what kind of sampling signal and how to sample can obtain the optimal reduced-order mode have become a research issue for scholars. By reviewing the relative literature about POD sampling, it is found that some scholars obtained POD reduced modes from chaos response signal [58, 59, 79, 87, 88, 111, 112], random response signal [87, 113–115], moment response signal [116–125], different sampling method [126–133], and constructing different autocorrelation matrix [134].

Kerschen et al. [88] used the finite element method to discretize fixed beam; they used two permanent magnets to act a nonlinear constraint force on free ends of the beam, and then they used the POD method to analyze the reduced order of the discrete high-dimensional nonlinear system. When the outer excitation frequency of this nonlinear system is fixed, they change the excitation amplitude, and then the periodic, almost periodic, and chaos movements have appeared. They found that the POD reduced-order mode obtained from the chaotic response signal is more effective than other nonchaotic response signals in reduction.

Amabili et al. [58, 59] used Donnell nonlinear cylindrical shell theory to go further into the large vibration issue of the simple supporting fluid filled cylindrical shell, and they used the Galerkin modal truncation method to obtain the accurate response of the original system with 16-DOF under the periodic excitation condition. However, they combined the POD method and Galerkin method and then found that only 3-DOF can be used to approximate the periodic, almost periodic response of a 16-DOF system, but extracting the POD modal response must be noticed. They also found that extracting POD reduced-order modal from chaos movement is more effective than from periodic and almost periodic response. Meanwhile, they used the method to analyze the bifurcation, set the external excitation frequency near the inherent fundamental frequency of the system, and selected the excitation amplitude as the bifurcation parameter. Then, it is found that when the external excitation amplitude varies within a large range, even if the reduced-order mode is extracted from the chaotic response, it cannot be guaranteed that the reduced-order system is similar to the bifurcation structure of the original system. The results indicate that the POD method is not robust in the wide-range parameter domain.

Bizon et al. [111] applied the POD method to research the reduced-order reaction dynamic modal of the one-dimensional tubular reactor with additional thermal cycling; according to calculation maximum orbital entropy index, they proved that the reduced modal fundamental function constructed from chaos sampling is optimal. Xie et al. [61, 62] applied Von Karman large deflection plate theory and combined the POD method and Galerkin method to study the wing aerodynamic vibration under the supersonic flow. They found that POD reduced-order modes obtained from chaotic response signal is more accurate than from period response to approximate the original system. Moreover, the reduced-order mode obtained from the chaotic motion under a certain parameter can also be applied to reduce other parameters in a certain parameter range. The research by the above scholars shows that the chaotic response signal in the complex system includes abundant physical process and state information of the original system, so we extract POD reduced-order modes from chaotic response signal in recent [62, 88].

The complex nonlinear system of actual engineering may appear the chaos motion in some parameter domain, but the complex movements are often avoided in design, processing, and maintenance in these complex structural systems. So, the chaos motion of the system is difficultly obtained in actual engineering. When constructing the signal, the POD reduced-order mode is sampled from the periodic signal of the normal operation of the complex system. The reduced-order modes obtained by this way may include less structural information of the original system, so it is difficult that a few reduced-order modes can obtain the optimal approximation of the original system. Under the random excitation, the system power spectrum is continuous, and the random response is as same as chaos response which has abundant physical processing and state information of the original system, and it is easy to obtain the random excitation, so some scholars [87, 113–115] researched that the POD reduced-order mode is constructed from the random response signal.

Kumar [87, 113] analyzed two high-dimensional nonlinear systems: one is nonlinear fixed beam system [88] and another is multiple spring-mass-damping chain system with Duffing nonlinear factors. They found that under a random excitation of certain bandwidth, the reduced-order modes extracted from random response can obtain better order reduction effect, but the frequency bandwidth of random excitation must cover the main frequency domain, and the extracted response contains many high-frequency components and requires more POD modes to obtain accurate reductions. The random signals were used by Yu and Chakravorty [114] to obtain a global optimal POD reduction mode. Segala [115] used the POD method to process the parameter modal reduction of the nonlinearly supported beams; they also found that the reduced-order mode constructed by random excitation signal can contain the main structural components of the original system with few modal numbers, thus obtaining a lower-dimensional reduced-order model. Therefore, when we do not know that chaotic motion occurs under what parameters in the high-dimensional system, the POD reduction mode can be obtained from the random response.

Some scholars used the POD method to reduce the order of high-dimensional complex systems and found that the better order reduction effect also can be obtained when sampling signal contains the transient response of the original system. In 1998, when Park and Lee [116] studied the flow optimization control issue for two-dimensional viscous flow, they found that the POD reduced-order mode sampled from the flow field displacement signal can only simulate effectively the flow of the original system if it contains the transient flow field signal of the original system. Terragni and Valero [117–121] applied the POD method to reduce the order of a class of continuous dissipative systems; they compared and analyzed the complex bifurcation characteristics of the original system in the parameter domain by the low-dimensional reduced-order system. They found that the reduced-order mode extracted from the snapshot containing the transient motion of the system has a larger parameter applicable range than the reduced-order mode extracted from the attractor.

Yang [121] applied the POD method to research the model reduction and parameter identification of high-dimensional chaotic/linear systems; it is found that the POD reduced-order mode constructed by unstable transient period response signal can approximate to the chaotic motion of the original system. Yu et al. and Lu et al. [122, 123] applied the POD method to reduce the order of high-dimensional nonlinear rotor systems with loose, crack, rubbing, and other faults. The system transient response, which has more motion information, contains free and forced vibrations under given initial conditions. Therefore, they also found that the POD reduced-order mode extracted from the transient motion signal can obtain better order reduction effect. Lu et al. [124, 125] also obtained the POD reduced-order sampled from the system transient response and reduced the order of the multidisk rotor-bearing system.

The above is the study about types of sampling signals of the POD method. However, the snapshot signals of the physical quantities of the original system obtained by numerical simulation or experiment involve how to sample, so some scholars have studied the sampling methods of POD methods [126–135]. Standard sampling methods include uniform sampling, random sampling, and hierarchical sampling [133]. Uniform sampling, which samples the motion state of each parameter in the system parameter space, requires a large computational cost. Especially, when the dimension of the parameter space is very high, uniform sampling will lose the application value. Random sampling may not be able to obtain information on some important areas of the parameter space [126]. It is difficult for the standard sampling method to obtain the construction of optimal reduced-order mode in high-dimensional parameter space, so scholars proposed some autocorrelation sampling methods: model constraint (MC-POD) method [126], greedy sampling (GS-POD) method [127–130], confidence domain (TR-POD) method [131], optimal component (OS-POD) method [132, 130], and some other sampling methods [132–135]. These methods mainly optimize snapshot signals by various sampling methods, which allow more state information of the original system to be included in different parameter domains or some important parameter domains, thereby obtaining an optimal POD reduced-order mode basis function.

For multifield, multiphysics coupled high-dimensional complex systems, we can extract multiple sets of state vector signals from the original system, such as velocity and displacement of the flow field, temperature field, pressure field, velocity, and displacement of the structure. It should be considered that which state vector should be chosen to construct the autocorrelation matrix and how to construct it. Kirby et al. [136] proposed putting together groups of state variables of the coupled system to construct a total state vector. By constructing the autocorrelation matrix by this state vector, the POD reduced mode of the coupled system is obtained and then this reduced-order mode was used to reduce the governing equations of the coupled system. This method has been successfully applied to model reduction of complex multifield coupled systems [92, 137–139]. However, the constructed total autocorrelation matrix contains the coupling of the state vectors of each parts of the system. This kind of coupling is generated by the total state vector in the process of numerically calculating the autocorrelation matrix, which is inconsistent with the coupling of the actual system, and the problem of nonconvergence of the reduced-order model often occurs. Brenner et al. [134] proposed that the state vectors of the components of the coupled system separately construct the autocorrelation matrix, and the mutual coupling term between the components is set to zero. The reduced-order model obtained by the total autocorrelation matrix constructed by this method is more accurate and convergent than the method proposed by Kirby et al.

In summary, the sampling issue of the POD method is very important for constructing reduced-order modes. The signal that constructs the reduced-order mode should contain more abundant and detailed physical processes and state information of the original system, such as chaotic signals, random signals, transient signals, and signals obtained by various sampling methods. The POD reduced-order modes sampled from these signals can obtain high-precision approximation to the original system with few DOFs. Now, scholars generally extract POD reduced modes from chaotic signals. When the system chaotic motion signal cannot be obtained in advance, a random response signal of a certain bandwidth or a transient motion signal can be used to extract the reduced mode. The above is a review of the POD method in the sampling issue, and then we will elaborate on another type of issue of the POD method.

##### 6.3. Parameter Adaptation of POD Method

The response of the complex dynamic system is related closely to system parameters, initial conditions, and so on. Therefore, the reduced-order mode obtained by the POD method usually lacks robustness [140, 141] when the system parameters change. In principle, the reduced-order mode under the new parameters should be reconstructed. However, constructing the POD reduced mode of the corresponding parameter by constructing the response of each parameter in the complex system is not allowed from the point of calculation cost. In order to reduce the order of the complex system in parameter range and ensure the parameter robustness of POD reduced-order mode, scholars have proposed many improved POD methods: global POD method [142–145], local POD method [60, 117, 146–149], adaptive POD method-POD modal interpolation method [141, 152–156], subspace angle interpolation method [157–162], Grassmann manifold tangent space interpolation method [163–171], and other adaptive POD methods.

The global POD method aims to construct a global POD reduced-order mode that covers the entire parameter domain by using a snapshot set composed of different parameter values in a certain parameter domain and then uses the reduced-order mode to obtain the reduced-order modal of system parameter domain. This method is very simple and easy to implement. However, this method extracts snapshot sets with different parameters and needs proper positioning in the parameter domain while how to reasonably determine the parameters of the snapshot signal is irregular. There may be multiple solutions in the parameter domain of complex strong nonlinear systems. This method requires a lot of POD reduced-order modes to obtain a more accurate reduced-order model of the original system. However, in practice, it has been proved that in many cases, this method is not reliable. Because there are numerous parameters affecting the dynamics of high-dimensional complex systems, the globally optimal POD reduction modes may not exist, which causes the method to lose the best approximation [161, 163–166].

The local POD method is a local reduction method. The whole parameter domain is divided into multiple subparameter domains. The POD method is used to reduce the system of each subparameter domain. The POD reduced-order mode is constructed by using the snapshot signals of some local parameter domains obtained in advance, and then the original system is projected onto the subspaces of these local reduced-order modes to obtain the price reduction model. Rapún and Vega [60] combined the local POD method with the Galerkin method to study the two one-dimensional parabolic equations (nonautonomous Fisher-like equations and complex Ginzburg–Landau equations). They proved the effectiveness of the method and its robustness in the numerical way in the range of local parameters.

Terragni and Valero [117] also combined the local POD method and the Galerkin method to study the two-dimensional roof-driven cavity flow issue. Compared with commercial software solution, this method can significantly reduce the calculation cost and be applied to reduce the order of different parameters within the range of system adjustable parameters, which is convenient for checking complex systems. Sahyoun and Djouadi [150] used the local POD method based on clustering vector space to reduce the order of high-dimensional nonlinear systems. Compared with the global POD method, the local POD method further reduces the system dimension and ensures the calculation accuracy. However, the method is divided into multiple subparameter domains, so the approximate solution is not smooth, and the robustness of the method in the parameter domain and the scope of its application are not yet rigorously proved [117, 148].

The adaptive POD method [150, 151] essentially achieves the best approximation for each parameter of the original system by updating the POD reduced-order mode within the parameter range. This method has better parameter robustness than the above two methods. In the literature about the POD method, there are many methods called adaptive POD methods, such as POD modal interpolation method [141, 152–156], subspace angular interpolation method [157–162], Grassmann manifold spatial interpolation method [163–171], and other adaptive POD methods [59, 118–120, 172].

The POD modal interpolation method is to use some interpolation methods (such as empirical discrete interpolation and Lagrange interpolation) to construct reduced-order modes of other parameters [141, 152–156] in the parameter domain by obtaining the reduced order of other parameters in the parameter domain. Xu and Lin [153] combined the POD method and empirical discrete interpolation method through greedy sampling to obtain interpolation points and reduced order of nonlinear parameter systems. Opmeer [154] combined the balanced-POD method and rational interpolation method to reduce the order of control systems.

Yao and Marques [155] used the POD method and the empirical discrete interpolation method to obtain the optimal POD reduction mode under the new parameters by training the radial basis function artificial neural network in the parameter space. They used the radial basis function artificial neural network of each interpolation point to reconstruct the flow field under the new parameters. The POD modal interpolation method is simple. However, the reduced-order mode is the orthogonal norm base vector, and the reduced-order modes, which are obtained by some direct interpolation methods under the new parameters, are no longer orthogonal normative. Therefore, in many cases, the accurate reduction model cannot be obtained by using the reduced-order mode obtained by direct interpolation. For example, when Lieu and Lesoinne [156] applied this method to analyze the aerodynamic problems of the F16 machine, accurate results were obtained at subsonic flight, but the method did not obtain correct results at transonic and supersonic speeds.

Generally, in the linear interpolation process of any two base vectors, their angles are not guaranteed to be linear interpolation. The subspace angle interpolation method is proposed based on the concepts of the main vector and the protagonist of the two subspaces [157]. The method can guarantee that the angles of the two interpolation basis vectors are also linear interpolation. By obtaining the POD reduced-order modal vector of some two parameters in the parameter domain in advance and then using the subspace angle interpolation method to obtain the POD reduced-order mode of other parameter values, the method has been successfully applied to the F16 machine in different freedoms by Lieu [158–161] for pneumatic analysis under the Mach number and angle of attack parameters. However, the subspace angle interpolation method is a low-order interpolation method, an accurate reduced-order model that cannot be obtained when two parameter values are far apart in the parameter space. The computational efficiency will be too low and lose the significance for order reduction when the distance is too close [160, 161, 165, 166].

In recent years, Amsallem proposed a more robust adaptive POD method, Grassmann manifold tangent spatial interpolation method [163–171]. The method is based on some concepts and mathematical conclusions in differential geometry, such as the Grassmann manifold, the calculation of the tangent space on the manifold, and the geodesic path. Amsallem [160, 161] used this method to study the aerodynamic problems of F16 and F18. Compared with POD modal interpolation and subspace angle interpolation, this method is not only suitable for subsonic and transonic aerodynamic analysis but also for supersonic aerodynamic analysis, which has good robustness. At the same time, Amsallem and Farhat [165, 166] also proved that the two-point interpolation of this method is equivalent to the subspace angle interpolation method. Amsallem and Cortial [167] also used this method to reduce the order of 24-DOF spring-mass-damping system and continuous wing structure system. Comparing the reduced-order model with the response of the original system, they proved its effectiveness and also showed that the method is universal and suitable for the reduction of other complex structural systems. Paquay [171] used this method to reduce the model of the nonlinear magnetic dynamic system and compared it with the direct POD model reduction and POD modal interpolation method, showing the superiority of the method.

The above are several major adaptive POD methods, and of course there are some other adaptive POD methods, such as Terragni and Valero [118–120] proposed an adaptive POD method in 2014 to analyze the distribution characteristics of complex systems in certain parameters. Because it is very difficult to analyze the bifurcation characteristics of complex systems in the parameter domain by direct numerical calculation, some scholars have considered the bifurcation characteristics of the low-dimensional reduced-order model to reflect the bifurcation characteristics of complex high-dimensional systems. Amabili et al. [59] used the POD method to analyze the bifurcation characteristics of complex structural dynamic systems in the parameter domain in 2006. However, this method has robustness problems in the parameter domain. They found that if the POD reduced-order mode is extracted from the chaotic motion signal, the method can be applied within a certain parameter range. However, for a larger parameter range, even if the method is sampled from chaotic signals, it is difficult to obtain an accurate reduced-order model, so they warned scholars not to use this method blindly.

In literature [172], Amabili et al. also compared the POD method with another reduced-order method of nonlinear structural dynamic systems, i.e., the nonlinear modal methods (NNMs). The results show that the POD method can obtain a more accurate reduced-order model than the nonlinear mode method (NNMs). However, Amabili and others used traditional POD methods, not adaptive POD methods. Terragni and Valero [118–120] determined the transient motion response of complex systems in advance with different parameter values through the Galerkin method, as a snapshot signal for updating the POD reduced-order mode in the parameter domain; they used the truncation error function to detect when the reduced-order modes should be updated and how many reduced-order modes need to be selected. The strategy for updating the POD reduced-order mode is to mix the modal vectors of the old and new parameters with appropriate weights. Then, they combined the POD method and the Galerkin method to obtain a bifurcation diagram of the complex Ginzburg–Landau equation with periodic parameters, almost periodic periods and chaotic motions on larger parameter range. Compared with the original system bifurcation diagram obtained by direct numerical calculation, it is proved that the proposed method has better robustness in a larger parameter range and can swiftly analyze the bifurcation characteristics of complex systems in the parameter range.

##### 6.4. Other Improved POD Methods and Related Issues

The above is the classification of various improved POD methods. Of course, there are some other improved POD methods, such as the POD-Galerkin method [58–62, 117–120, 173], which regards the POD reduced-order modal basis function as the orthogonal modal function of the Galerkin method. The POD method based on frequency domain obtains the POD reduced-order modal basis function in the frequency domain, balanced-POD method (BPOD) [3, 154, 176–178]. These improved POD methods are all proposed by combining traditional POD methods with other methods. In this paper, various improved POD methods are mainly divided into two categories related to sampling and parameters. On the one hand, the POD method is essentially a statistical analysis method, which involves how to sample, which method should be used to sample it, when to sample the signal, which part of the signal should be extracted, and how long the signal should be extracted. These are directly related to sampling. On the other hand, from the basic theory of the POD method, the original system is projected onto the subspace of the first few reduced modes by the POD reduced mode of the snapshot signal constructed by a certain parameter value of the original system. The resulting reduced-order model is only the best approximation of the original system for this particular parameter. Since the infinitesimal neighborhood of the parameter point can be linearized, the method can also obtain a better approximation for the system near the specific parameter value [79, 179]. We hope that the POD reduced-order mode obtained from a certain parameter can also be applied to the model reduction of other parameters in the system so that the reduced-order model under other parameters also has the optimal approximation characteristics. However, such a conclusion is not given in the basic theory of the POD method, so the POD method lacks robustness in the parameter domain.

Some POD methods for solving parameter robustness have been proposed, but each has its own advantages and disadvantages. At present, the adaptive POD method is more mature in dealing with parameter robustness, but it still needs to be continuously developed. For example, the adaptive POD method, Grassmann manifold tangent space interpolation method, proposed in recent years has achieved good results in dealing with parameter robustness problems. However, we know that manifold tangent spaces have local properties, and the method still has limitations in large parameter domain.

The POD method projects a high-dimensional system onto the space spanned by a few reduced-order modes, which is essentially a projection reduction method. The adaptive POD method can solve the problem of parameter robustness of the reduced-order system to certain extent. However, the reduced-order model on the parameter domain is obtained by continuously updating the POD reduced-order mode or adjusting the number of modes, resulting in the reduced-order model being some discrete numerical equations, and there is no invariant reduced-order mathematical model in the parameter domain. Even the reduced-order model dimensions of different subparameter domains are inconsistent. Therefore, a low-dimensional model that can approximately reflect the dynamic characteristics of complex high-dimensional systems cannot be obtained through the adaptive POD method in the parameter domain so that it is difficult to carry out in-depth theoretical analysis of high-dimensional systems. How to use the POD method to obtain the invariant reduced-order model of the high-dimensional system in the parameter domain? The POD method is a projection order reduction method, so it is impossible to update the POD reduced-order mode in the parameter domain, and only the reduced-order mode of the parameter domain can be obtained by constructing a reduced-order mode which is optimal in the entire parameter domain. What conditions can the POD reduced-order modes reduce the parameter domain or is there an optimal POD reduced-order mode for the entire parameter domain? How many DOF can the invariant reduced-order model reflect the dynamics of the original system over the entire parameter domain? In what manner does the method interact with uncertainty [180–184] for the high-dimensional mechanical systems?

Computers nowadays can handle real-time calculations using the finite element method. The finite element model and POD order reduction techniques will be applied to study real-time computer modeling [185], e.g., surgical simulations [186, 187] (for doctor training, sometimes in VR) and hybrid simulation [188, 189] (HIL of mechanical systems where mode superposition or other methods are applied to reduce the model order). In real-time simulations, high computational effort to reduce the model is worth it even if the real-time computations can be slightly minimized.

Machine learning can also be used for model order reduction [190, 191]; manifold learning and manifold processes for model order reduction should be discussed in detail in our future work. Chaos features can improve the efficiency of reduction models [192–194] and the combination method of uncertainty quantification [195, 196], and model order reduction will be efficient to study the complex dynamic systems. Another point is the systems with substructures. The way it works and method for choosing the approximation order are elaborated by researchers [197, 198]. As a commonly used method in model updating, the quadratic inverse eigenvalue problem (QIEP) aims to construct the mass, stiffness, and damping matrices and can be employed to assist model order reduction for large-scale engineering systems [199, 200].

The above questions are all about the POD order reduction methods that deserve further study in the future.

#### 7. Conclusions and Outlooks

The model order reduction for the high-dimensional complex nonlinear system is one of the important issues in the field of engineering research, and it is one of the advanced issues in the area of nonlinear science research. Many scholars from various countries have obtained some mature order reduction methods for a long research process. However, modern engineering structural systems are more and more complex, operating conditions are complex and variable, and various nonlinear factors are coupled to each other. As a result, many model reduction methods are no longer applicable. This is a challenge to the study of dynamic characteristics of high-dimensional complex systems and system parameter optimization design. According to the characteristics of each order reduction method, the worthy problems of further study on the order reduction of high-dimensional complex system models in the future are as follows:(1)Multiple order reduction methods are combined for second-level order reduction on high-dimensional complex systems. For example, the complex system is divided into several substructures, and the modal synthesis method is used to reduce the order, and then the center manifold method or the L-S method is used for further analysis. For another example, the complex system is reduced by the Galerkin method or the POD method, and then the center manifold or the L-S method is used for the order reduction study. The center manifold and the L-S method can reserve the topology properties of the original system. Therefore, the reduced-order model of the complex structural system can be obtained by the Galerkin method, the modal synthesis method, or the POD method, and then the theoretical research is carried out by using the center manifold and the L-S method.(2)The adaptive model order reduction method with system parameter variation and various nonlinear model reduction methods should be further studied. For example, the adaptive order reduction method based on Grassmann manifold tangent space interpolation has local properties in manifold tangent space, and it is worthy of further study to solve robustness problem in large parameter domain. Some new nonlinear model order reduction methods have been proposed based on neural networks, local reduced-order bases, and manifold learning, but each has its own advantages and disadvantages, which still need further study.(3)A low-dimensional model in the parameter domain can approximately reflect the dynamic characteristics of complex high-dimensional systems which cannot be obtained by the adaptive POD method, and it is difficult to carry out in-depth theoretical analysis of high-dimensional systems. How to use the POD method to obtain the invariant order reduction model of the high-dimensional system in the parameter domain? The POD method is a projection reduction method, and the POD reduced-order modes cannot be updated in the parameter domain. Therefore, order reduction of the parameter domain can be obtained only by constructing a reduced-order mode which is optimal in the entire parameter domain. Under what conditions can the POD reduced-order modes reduce the dimension of parameter domain or is there an optimal POD reduction mode for the entire parameter domain? How many DOF can the invariant order reduction model reflect the dynamics of the original system over the entire parameter domain? The POD order reduction method can be applied for further study for these problems.(4)Special attention should be taken to the model reduction of interconnected systems in order to preserve the integrity and interconnection structure among different subsystems. For example, the dual-rotor system contains high pressure and low pressure rotors, the system is complex and high dimensional, and the nonlinearity is strong between the joints of high and low pressure rotors. It will be a challenge to apply the POD method in this kind of system.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This study was funded by the National Natural Science Foundation of China (grant nos. 12072263, 11802235, and 11972295), the Natural Science Foundation of Shaanxi Province (grant no. 2020JQ-129), and the Aviation Engine Innovation Center of National Defense Science, Technology and Industry (grant no. CXZX-2019-001).