FIELD OF THE INVENTION

[0001]
This invention relates generally to methods and systems for efficient model reduction and system identification of dynamic systems, and more specifically, for linear dynamic systems having multiple inputs using a new Eigensystem Realization Algorithm.
BACKGROUND OF THE INVENTION

[0002]
Most modern dynamic models are constructed based on finite spatial discretizations of continuous systems, often resulting in a considerable number of degrees of freedom in the model. Consequently, for fast and efficient estimation of dynamic behavior, as well as optimizations and closedloop control designs, a model reduction is typically performed. Desirable characteristics of a reducedorder dynamic model include that the size of the system must not be too large, the model must be robust and have a good fidelity, it must be in the statespace, time domain formulation for implementation of active control systems and nonlinear time analysis, and finally, the reduction process itself must not be too expensive.

[0003]
There have been many model reduction methods available, but most of them require modifying the main frame of the computational model, and are prone to a long model construction time if the model is subjected to many driving inputs. The latter in particular is true in unsteady Computational Fluid Dynamic (CFD) applications where the moving solid boundary is often described by many structural mode inputs. For example, a typical commercial airplane simulation may involve as many as 200 structural modes.

[0004]
Recently, Eigensystem Realization Algorithm (ERA) was used successfully in aeroelastic flutter predictions using a computational fluid dynamics code (CFD) developed by the NASA Langley Research Center known as CFL3D. (See Juang, J. N., Applied System Identification, Prentice Hall Englewood Cliffs, N.J., 1994; Silva, W. A., and Bartels, R. E., Development of ReducedOrder Models for Aeroelastic Analysis and Flutter Prediction Using CFL3Dv6.0 Code, AIAA20021596; and Hong, M. S., Bhatia, K. G., SenGupta, G., Kim, T., Kuruvila, G., Silva, W. A., Bartels, and R., Biedron, R., Simulations of a TwinEngine Transport Flutter Model In the Transonic Dynamics Tunnel, IFASD Paper 2003US44). The ERA method, which is usually used as a system identification technique, has a very attractive feature in that unlike model reduction methods based on Galerkin scheme (e.g., Dowell, E. H., Hall, K. C., and Romanowski, M. C., Eigenmode Analysis in Unsteady Aerodynamics: ReducedOrder Models, Applied Mechanics Review, Vol. 50, No. 6, 1997, pp. 371386), there is no need for online implementation of the algorithm. That is, the ERA method is a postprocessing tool that identifies and generates system matrices based on the input and output data alone.

[0005]
FIG. 1 is a schematic view of the ERA (a.k.a. Pulse/ERA) method 100 of modeling a dynamic system in accordance with the prior art. For simplicity, only its fundamental statespace realization theory which is attributed to Ho and Kalman is discussed. (see Ho, B. L. and Kalman, R. E., Effective Construction of Linear StateVariable Models from Input/Output Data, Proceedings of the 3rd Annual Allerton Conference on Circuit and System Theory, 1965, pp. 449459). For a general description of the Pulse/ERA method, see Juang, J.N. and Pappa, R. S., An Eigensystem Realization Algorithm for Modal Parameter Identification and Model Reduction, Journal of Guidance, Control, and Dynamics, Vol. 8, 1985, pp. 620627.

[0006]
In the Pulse/ERA method 100, it is assumed that the dynamic system under consideration can be described in the following finitedimensional, discretetime form:
x ^{n+1} =Ax ^{n} Bu ^{n} (1)
y ^{n} =Cx ^{n} +Du ^{n} (2)
where
x≡(L×1) state vector (3)
u≡(N_{i}×1) input vector (4)
y≡(N_{0}×1) output vector (5)
The system matrices (A, B), (A, C) are assumed controllable and observable. First, given M+1 equally distributed time steps,
$\begin{array}{cc}{t}^{n}\equiv n\text{\hspace{1em}}\Delta \text{\hspace{1em}}t& \left(n=0,1,2,\dots \text{\hspace{1em}},M\right),\end{array}$
for a single ith input vector the system output is sampled subjected to the unit pulse,
$\begin{array}{cc}{u}_{i}^{n}={\delta}^{\mathrm{ni}}\equiv \left\{\begin{array}{cc}1& \left(n=0\right)\\ 0& \left(\mathrm{ni}=1,2,\dots \text{\hspace{1em}},M\right)\end{array}\right\}& \left(6\right)\end{array}$
Assuming zero initial condition, x^{0}=x(0)=0, one obtains
$\begin{array}{cc}{y}^{n}\text{?}={d}_{i}\text{}{y}^{1}\text{?}={\mathrm{Cb}}_{i}\text{}{y}^{1}\text{?}={\mathrm{CA}}_{i}\text{}{y}^{3}\text{?}={\mathrm{CA}}^{2}{b}_{i}\text{}\text{\hspace{1em}}\vdots \text{\hspace{1em}}\vdots \text{}{y}^{M}\text{?}={\mathrm{CA}}^{M1}{b}_{i}\text{}\mathrm{where}\text{}& \left(7\right)\\ {b}_{i}\equiv \text{?}\text{}\mathrm{nth}\text{\hspace{1em}}\mathrm{column}\text{\hspace{1em}}\mathrm{of}\text{\hspace{1em}}B& \left(8\right)\\ {d}_{i}\equiv \text{?}\text{}\mathrm{nth}\text{\hspace{1em}}\mathrm{column}\text{\hspace{1em}}\mathrm{of}\text{\hspace{1em}}D\text{}\text{?}\text{indicates text missing or illegible when filed}& \left(9\right)\end{array}$
The constant matrices in the above sequence are known as the Markov parameters. As shown in FIG. 1, for all inputs 102 into the dynamic system, the Markov parameters are computed (block 104) creating the sequence of pulseresponse matrices:
$\begin{array}{cc}\begin{array}{c}\begin{array}{cc}{Y}^{n}\equiv \uf603\begin{array}{cccc}{y}_{1}^{n}& {y}_{2}^{n}& \dots & {y}^{n}\text{?}\end{array}\uf604& \left(n=0,1,2,\dots \text{\hspace{1em}}M\right)\end{array}\\ \text{?}\text{indicates text missing or illegible when filed}\end{array}& \left(10\right)\end{array}$
Next, based on the system Markov parameters (block 104), the Pulse/ERA method 100 defines N_{0}×(Ni×(M−1)) Hankel matrices:
$\begin{array}{cc}{H}_{0}\equiv \uf603\begin{array}{cccc}{Y}^{1}& {Y}^{2}& \dots & {Y}^{M1}\end{array}\uf604\text{}\text{\hspace{1em}}\equiv C\uf603\begin{array}{ccccc}I& A& {A}^{2}& \dots & {A}^{M2}\uf604B\end{array}& \left(11\right)\\ {H}_{1}\equiv \uf603{Y}^{2}\text{\hspace{1em}}{Y}^{3}\text{\hspace{1em}}\dots \text{\hspace{1em}}{Y}^{M}\uf604\text{}\text{\hspace{1em}}\equiv C\uf603\begin{array}{cccc}A& {A}^{2}& \dots & {A}^{M1}\uf604B\end{array}& \left(12\right)\end{array}$

[0007]
Singularvaluedecomposition (SVD) of H_{0 }yields
$\begin{array}{cc}{H}_{0}\equiv U\text{\hspace{1em}}\sum \text{\hspace{1em}}{V}^{T}\text{}\simeq \uf603\begin{array}{cc}{U}_{R}& {U}_{D}\end{array}\uf604\left[\begin{array}{cc}\sum _{R}& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}{V}_{R}^{T}\\ {V}_{D}^{T}\end{array}\right]\text{}={U}_{R}\underset{R}{\sum ^{1/2}}\underset{R}{\sum ^{1/2}}{V}_{R}^{T}& \left(13\right)\end{array}$
where R≡rank(H_{0}). Finally, a balanced realization of the system under question is obtained by pseudoinverting various submatrix components (block 108 of FIG. 1) as follows:
$\begin{array}{cc}D={Y}^{0}\left({N}_{0}\times 1\right)& \left(14\right)\\ C\simeq {U}_{R}\sum _{R}^{1/2}\left({N}_{0}\times R\right)& \left(15\right)\\ B\simeq \mathrm{the}\text{\hspace{1em}}\mathrm{first}\text{\hspace{1em}}{N}_{i}\text{\hspace{1em}}\mathrm{columns}\text{\hspace{1em}}\mathrm{of}\text{\hspace{1em}}\underset{R}{\sum ^{1/2}}{V}_{R}^{T}\text{\hspace{1em}}\left(\mathrm{Ri}\times {N}_{i}\right)& \left(16\right)\\ A\simeq \underset{R}{\sum ^{1/2}}{U}_{R}^{T}{H}_{1}{V}_{R}\underset{R}{\sum ^{1/2}}\text{\hspace{1em}}\left(R\times R\right)& \left(17\right)\end{array}$

[0008]
Since R<<L, the above model represents a reducedorder realization of the original system. Note that the realization is optimal in that it is balanced between inputs and outputs. However, the total number of samples taken is N_{i}×(M+1), which increases proportional to the number of inputs. Also, for an accurate system identification H_{0 }must have sufficient columns and rows, i.e., N_{i}×(M−1), ≧R and N_{0}≧R.

[0009]
For a very large data set Y^{n }with many time steps and a large number of inputs, the Eigensystem Realization Algorithm/Data Correlations (ERA/DC) method may be preferred. As described, for example, in Juang, J. N., Applied System Identification, Prentice Hall Englewood Cliffs, N.J., 1994, the ERA/DC method may be used to compress the amount of data and reduce the computation time required for the SVD of the Hankel matrix.

[0010]
Although desirable results have been achieved using the prior art computational methods, there is room for improvement. For example, if the unsteady CFD model is driven by multiple structural inputs, as described above, the computation time required to obtain all the pulse responses increases proportional to the number of the inputs, making the ERA method undesirably slow and inefficient. Therefore, a need exists for improved methods for model reduction and system identification of largescaled linear dynamic systems having multiple inputs
SUMMARY OF THE INVENTION

[0011]
The present invention is directed to methods and systems for model reduction and system identification of dynamic systems having multiple inputs. Embodiments of the methods in accordance with the present invention may advantageously reduce the model construction and system identification time, especially for a largescaled system, and compress the amount of output data. Furthermore, they do not require modifying the original code, and take only input and output data for the model construction. That is, it is 100% a postprocessing tool.

[0012]
In one embodiment, a method of model reduction and system identification includes generating a plurality of statistically independent random numbers for use as input signals, and performing a singularvaluedecomposition directly on the system response due to the simultaneous excitation of all the inputs, also known as the SingleCompositeInput (SCI). In alternate embodiments, the method further includes sampling individual pulse responses for the first time and second time steps. Alternately, the input signals are filtered through a lowpass filter. The plurality of input signals may further include applying multiple step inputs in a sequential manner, and applying multiple pulse inputs in a sequential manner.

[0013]
In an alternate embodiment, a method of model reduction and system identification includes generating a plurality of statistically independent random numbers for use as input signals, performing a singularvaluedecomposition directly on the system response due to the SingleCompositeInput (SCI), sampling individual pulse responses for the first time and second time steps. The method further includes defining the Hankellike matrices H_{c0 }and H_{c1 }as follows:
$\begin{array}{cc}{H}_{\mathrm{c0}}\equiv \uf603\begin{array}{cccc}{y}_{\mathrm{c0}}^{1}& {y}_{\mathrm{c0}}^{2}& \dots & {y}_{\mathrm{c0}}^{M1}\end{array}\uf604\text{}\text{\hspace{1em}}=C\uf603\begin{array}{cccc}{x}^{1}& {x}^{2}& \dots & {x}^{M1}\end{array}\uf604& \left(25\right)\\ {H}_{\mathrm{c1}}\equiv \uf603\begin{array}{cccc}{y}_{\mathrm{c1}}^{1}& {y}_{\mathrm{c1}}^{2}& \dots & {y}_{\mathrm{c1}}^{M1}\end{array}\uf604\text{}\text{\hspace{1em}}=\mathrm{CA}\uf603\begin{array}{cccc}{x}^{1}& {x}^{2}& \dots & {x}^{M1}\end{array}\uf604\text{}\mathrm{SVD}\text{\hspace{1em}}\mathrm{of}\text{\hspace{1em}}{H}_{\mathrm{c0}}\text{\hspace{1em}}\mathrm{yields}& \left(26\right)\\ {H}_{\mathrm{c0}}\equiv U\sum {V}^{T}\text{}\text{\hspace{1em}}\simeq \uf603\begin{array}{cc}{U}_{R}& {U}_{D}\end{array}\uf604\left[\begin{array}{cc}\sum _{R}& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}{V}_{R}^{T}\\ {V}_{D}^{T}\end{array}\right]\text{}\text{\hspace{1em}}={U}_{R}\underset{R}{\sum ^{1/2}}\underset{R}{\sum ^{1/2}}{V}_{R}^{T}& \left(27\right)\end{array}$
Finally, the method identifies the system matrices (A, B, C, D) by a leastsquare approximation as follows:
$\begin{array}{cc}D={Y}^{0}& \left(28\right)\\ C\simeq {U}_{R}\underset{R}{\sum ^{1/2}}& \left(29\right)\\ B\simeq \underset{R}{\sum ^{1/2}}{U}_{R}^{T}{Y}^{1}& \left(30\right)\\ A\simeq \underset{R}{\sum ^{1/2}}{U}_{R}^{T}{H}_{\mathrm{c1}}{V}_{R}\underset{R}{\sum ^{1/2}}& \left(31\right)\end{array}$
BRIEF DESCRIPTION OF THE DRAWINGS

[0014]
The preferred and alternative embodiments of the present invention are described in detail below with reference to the following drawings.

[0015]
FIG. 1 is a schematic view of the Pulse/ERA method in accordance with the prior art;

[0016]
FIG. 2 is a schematic view of the SCI/ERA method in accordance with an embodiment of the present invention;

[0017]
FIG. 3 is a schematic view of a vortex lattice model for computing an unsteady flow field around a flat, rectangular wing in accordance with an embodiment of the present invention;

[0018]
FIG. 4 shows three sets of random numbers generated as inputs to the vortex lattice formulation of FIG. 3;

[0019]
FIG. 5 shows sets of generalized aerodynamic forces computed using the vortex lattice formulation of FIG. 3;

[0020]
FIG. 6 shows a set of aeroelastic eigenvalues computed using the vortex lattice formulation of FIG. 3 at a first scale for V=80 m/sec;

[0021]
FIG. 7 shows a set of aeroelastic eigenvalues computed using the vortex lattice formulation of FIG. 3 at a second scale for V=80 m/sec;

[0022]
FIG. 8 shows an aeroelastic time response of the first structural mode computed using the vortex lattice formulation of FIG. 3 for V=80 m/sec;

[0023]
FIG. 9 shows an aeroelastic time response of the second structural mode computed using the vortex lattice formulation of FIG. 3 for V=80 m/sec;

[0024]
FIG. 10 shows a set of aeroelastic eigenvalues computed using the vortex lattice formulation of FIG. 3 at a first scale for V=121.2 m/sec;

[0025]
FIG. 11 shows a set of aeroelastic eigenvalues computed using the vortex lattice formulation of FIG. 3 at a second scale for V=121.2 m/sec;

[0026]
FIG. 12 shows an aeroelastic time response of the first structural mode computed using the vortex lattice formulation of FIG. 3 for V=121.2 m/sec;

[0027]
FIG. 13 shows an aeroelastic time response of the second structural mode computed using the vortex lattice formulation of FIG. 3 for V=121.2 m/sec;

[0028]
FIG. 14 shows a graph of model construction time of different ERA methods versus number of inputs for the vortex lattice formulation of FIG. 3;

[0029]
FIG. 15 shows a graph of CPU seconds spent constructing various ERA models in accordance with alternate embodiments of the present invention (on which FIG. 14 is based);

[0030]
FIG. 16 is a representative CFL3D simulation model in accordance with an embodiment of the present invention;

[0031]
FIG. 17 shows four structural mode shapes used for the CFL3D simulation model of FIG. 16; and

[0032]
FIG. 18 shows a Vg plot of two different reducedorder aeroelastic models based on the CFL3D simulation in accordance with an embodiment of the present invention.
DETAILED DESCRIPTION

[0033]
The present invention relates to methods and systems for model reduction and system identification of dynamic systems having multiple inputs. Many specific details of certain embodiments of the invention are set forth in the following description and in FIGS. 218 to provide a thorough understanding of such embodiments. One skilled in the art, however, will understand that the present invention may have additional embodiments, or that the present invention may be practiced without several of the details described in the following description.

[0034]
In general, embodiments of methods in accordance with the present invention include a new, efficient discretetime domain system identification and model reduction method for largescaled linear dynamic systems with multiple inputs. In one embodiment, a method in accordance with the present invention is based on a modification of the classical Eigensystem Realization Algorithm (ERA) and a simultaneous injection of multiple inputs, so called the SingleCompositeInput (SCD). Since the system response is sampled almost exclusively for the single representative input instead of the multiple individual inputs, embodiments of the present invention can significantly reduce the model construction time as well as the amount of the sampled data. Derivation of the new algorithm in accordance with the present disclosure may include performing a singular value decomposition using a single set of output measurements that are not necessarily attributed to pulse inputs. Representative simulations performed using embodiments of the present invention exhibit reduced computation times and excellent results obtained from the reducedorder models, thereby showing great potential of the present invention as a linear system identification and model reduction tool for largescaled systems subjected to multiple inputs.

[0035]
Unless otherwise stated, the following nomenclature is used throughout the following detailed description:
 A B C D System matrices
 A_{d1 }A_{d2 }Aeroelastic system matrices
 A_{s }B_{s }C_{s }Structural system matrices
 b Reference length
 C_{ij }Crosscorrelation coefficient
 E Expected value
 K Covariance matrix defined in (48)
 L Dimension of original system
 M Number of time or frequency samples or Mach number
 m c k Mass, damping, stiffness matrices
 p (R_{1}×1) generalized coordinate vector
 q Dynamic pressure (=½p V^{2})
 R Number of chosen singular modes or the dimension of realized model
 R_{1 }Number of chosen KL modes
 s Laplace variable
 t Real time
 u (N_{i}×1) input or generalized structural coordinate vector
 u_{i }ith input or structural coordinate
 V Air speed
 x (L×1) state or aerodynamic state vector
 X Fourier transform of x
 y (N_{0}×1) output vector or (N_{i}×1) generalzied aerodynamic force vector
 y^{n} _{i }Pulse response due to ith input
 φ_{i }KL mode
 Φ KL modal matrix
 ρ Air density
 τ Reduced time
$\left(\equiv \text{?}\right)$
$\text{?}\text{indicates text missing or illegible when filed}$
 ω Frequency (rad/sec)
 ω_{c }Maximum cutoff frequency

[0065]
Subscripts
 i Input
 o Output
 R Reduced
 ref Reference
 s Structure

[0071]
FIG. 2 is a schematic view of the SingleCompositeInput (SCI)/ERA method 200 in accordance with an embodiment of the present invention. As described more fully below, the SCI/ERA method 200 is based on the premise that for a linear system, one can apply the multiple inputs simultaneously and get all the system responses that are necessary for the model reduction. Since the computational model needs to be executed almost exclusively for the representative input, the model construction time can be significantly reduced.

[0072]
Embodiments of the present invention include implementing the SCI/ERA method for fast and efficient model reduction of linear, finitedimensional, discretetime systems with multiple inputs. To accommodate the SCI within the framework of the Pulse/ERA, it is necessary to modify the original algorithm. In particular, as shown below, the new formulation does not rely on the system Markov parameters explicitly. Instead, it performs the singularvaluedecomposition (SVD) directly on the output measurements that are in general not attributed to pulse inputs. Statistically independent random numbers are used in lieu of the pulses for the multiple input signals. Naturally, embodiments of the present invention can also be used towards system identification provided that all the time measurements are available from experiments. Application of the SCI/ERA method 200 to computational fluid dynamic systems and formulation of reducedorder aeroelastic models are presented below, where it is shown that depending on how the displacement and velocity inputs on the moving boundary are treated, two different kinds of reducedorder aerodynamic and aeroelastic models may be generated.

[0073]
For a demonstration of the proposed method, two CFD models are considered in the examples below. They are the Vortex Lattice Model (VLM) for inviscid, subsonic, incompressible flow, and the CFL3D for viscous, transonic, compressible flow. Reducedorder aeroelastic models are also constructed by combining the reduced aerodynamic models and the structural system. It is shown that not only does the new method shorten the model construction time substantially, but the accuracy of the resulting reducedorder models is excellent. The proposed scheme has a great potential as a linear system identification and model reduction tool for largescaled systems subjected to multiple inputs.

[0074]
With reference to FIG. 2, the SCI/ERA method 200 proceeds as follows. First, at a block 202, individual pulse responses are sampled for the first two time steps:
$\begin{array}{cc}{Y}^{0}=\text{?}\begin{array}{cccc}{y}_{1}^{0}& {y}_{2}^{0}& \dots & {y}_{{N}_{i}}^{0}\end{array}\text{?}& \left(18\right)\\ {Y}^{1}=\text{?}\begin{array}{cccc}{y}_{1}^{1}& {y}_{2}^{1}& \dots & {y}_{{N}_{i}}^{1}\end{array}\text{?}\text{}\text{?}\text{indicates text missing or illegible when filed}& \left(19\right)\end{array}$

[0075]
Next, at a block 204, an SCI is constructed as follows:
$\begin{array}{cc}{b}_{\mathrm{SCI}}^{n}\equiv \sum _{\text{?}=1}^{{N}_{i}}{b}_{i}{r}_{i}^{n\text{\hspace{1em}}}\text{\hspace{1em}}\left(\mathrm{for}\text{\hspace{1em}}\mathrm{states}\right)& \left(20\right)\\ {d}_{\mathrm{SCI}}^{n}\equiv \sum _{\text{?}=1}^{{N}_{i}}{d}_{i}{r}_{i}^{n\text{\hspace{1em}}}\text{\hspace{1em}}\left(\mathrm{for}\text{\hspace{1em}}\mathrm{outputs}\right)\text{}\mathrm{where}& \left(21\right)\\ {r}^{n}\text{?}\equiv a\text{\hspace{1em}}\mathrm{sequence}\text{\hspace{1em}}\mathrm{of}\text{\hspace{1em}}\mathrm{arbitrary}\text{\hspace{1em}}\mathrm{numbers}\text{}\text{?}\text{indicates text missing or illegible when filed}& \left(22\right)\end{array}$

[0076]
To ensure independency of the inputs, it is desirable that the signals be as uncorrelated as possible. In an ideal case they would be statistically uncorrelated random signals, i.e., C_{ij}(m)=E[r^{n} _{i }r^{n−m} _{i}]0 for i≠j, but they are hard to construct for numerical analysis.

[0077]
Subject to the SCI, at a block 206, the SCI/ERA method 200 samples the system response y^{n }for n=0, 1, 2, . . . , M, to obtain:
$\begin{array}{cc}{y}_{\mathrm{c0}}^{n}\equiv {\mathrm{Cx}}^{n}\text{}\text{\hspace{1em}}={y}^{n}\sum _{\text{?}=1}^{{N}_{i}}{y}^{n}\text{?}{r}^{n}\text{?}& \left(23\right)\\ {y}_{\mathrm{c1}}^{n}\equiv {\mathrm{CAx}}^{n}\text{}\text{\hspace{1em}}={y}^{n+1}\sum _{\text{?}=1}^{{N}_{i}}{y}^{0}\text{?}{r}^{n+1}\text{?}\sum _{\text{?}=1}^{{N}_{i}}{y}^{1}\text{?}{r}^{n}\text{?}\text{}\text{?}\text{indicates text missing or illegible when filed}& \left(24\right)\end{array}$

[0078]
Note that
${y}_{\mathrm{c0}}^{n},{y}_{\mathrm{c1}}^{n}$
are measurements of the states in the reduced dimension of C and C A.

[0079]
Next, at a block 208, similar to the Hankel matrices, the SCI/ERA method 200 defines the following Hankellike matrices:
$\begin{array}{cc}{H}_{\mathrm{c0}}\equiv \uf603\begin{array}{cccc}{y}_{\mathrm{c0}}^{1}& {y}_{\mathrm{c0}}^{2}& \dots & {y}_{\mathrm{c0}}^{M1}\end{array}\uf604\text{}\text{\hspace{1em}}=C\uf603\begin{array}{cccc}{x}^{1}& {x}^{2}& \dots & {x}^{M1}\end{array}\uf604& \left(25\right)\\ {H}_{\mathrm{c1}}\equiv \uf603\begin{array}{cccc}{y}_{\mathrm{c1}}^{1}& {y}_{\mathrm{c1}}^{2}& \dots & {y}_{\mathrm{c1}}^{M1}\end{array}\uf604\text{}\text{\hspace{1em}}=\mathrm{CA}\uf603\begin{array}{cccc}{x}^{1}& {x}^{2}& \dots & {x}^{M1}\end{array}\uf604\text{}\mathrm{SVD}\text{\hspace{1em}}\mathrm{of}\text{\hspace{1em}}{H}_{\mathrm{c0}}\text{\hspace{1em}}\mathrm{yields}& \left(26\right)\\ {H}_{\mathrm{c0}}\equiv U\sum {V}^{T}\text{}\text{\hspace{1em}}\simeq \uf603\begin{array}{cc}{U}_{R}& {U}_{D}\end{array}\uf604\left[\begin{array}{cc}\sum _{R}& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}{V}_{R}^{T}\\ {V}_{D}^{T}\end{array}\right]\text{}\text{\hspace{1em}}={U}_{R}\underset{R}{\sum ^{1/2}}\underset{R}{\sum ^{1/2}}{V}_{R}^{T}& \left(27\right)\end{array}$
where R≡rank(H_{co}). The size of the above matrices is N_{0}×(M−1), N_{i }times smaller than the previous H_{0}, H_{1 }defined in the Pulse/ERA.

[0080]
Thus, at a block 210, the new realization is then
$\begin{array}{cc}D={Y}^{0}& \left(28\right)\\ C\simeq {U}_{R}\sum _{R}^{1/2}\text{\hspace{1em}}& \left(29\right)\\ B\simeq \sum _{R}^{1/2}\text{\hspace{1em}}{U}_{R}^{T}{Y}^{1}& \left(30\right)\\ A\simeq \sum _{R}^{1/2}\text{\hspace{1em}}{U}_{R}^{T}{H}_{\mathrm{c1}}{V}_{R}\sum _{R}^{1/2}\text{\hspace{1em}}& \left(31\right)\end{array}$

[0081]
Unlike the Pulse/ERA method 100 described above with reference to FIG. 1, the SCI/ERA method 200 may be optimal in that it may be balanced between states and outputs. As in the previous method 100, for an accurate realization H_{co }must have (M−1)≧R and N_{0}≧R. However, using the SCI/ERA method 200 in accordance with the present invention, the total number of samples taken may be only M+1+2×N_{i }which is much less than the previous N_{i}×(M+1) when M samples of pulse response are taken for each input.

[0082]
It will be appreciated that compared to the prior art Pulse/ERA method 100, embodiments of the present invention may advantageously require a much smaller set of time measurements, thereby reducing significantly both the computation time and the bulk of the output data. Furthermore, the new H_{co}, H_{ci }are constructed based on the sampled system states subjected to combined random inputs, and as such they are not directly related to the Markov parameters. However, at block 202 embodiments of the present invention (e.g. the SCI/ERA method 200 shown in FIG. 2) do require the first two pulse responses
${y}_{i}^{0}$
and
${y}_{i}^{1}$
for each input in order to estimate the state measurements according to Equations (23) and (24).

[0083]
Although the use of random signals is described, other types of signals can also be used for the SCI provided that they are statistically independent. For a linear system, any arbitrary response contains the fundamental characteristics under the assumption that the system is controllable and observable. Embodiments of the present invention may advantageously employ this principle, along with the principle of superposition.

[0084]
An alternate embodiment to the SCI/ERA method discussed above will now be described. As stated above, an important requirement in any ERA method is that to improve the accuracy of the model construction, one must have a sufficient number of output measurements, more than the number of singular modes that are extractable from the Hankel or Hankellike matrices. Failure to satisfy this requirement implies that we don't have enough modes to approximate the system response. When the requirement is not met, assuming again that (A, C) is observable we can expand the data matrices by sampling additional responses as follows:
$\begin{array}{cc}\begin{array}{c}{H}_{c\text{\hspace{1em}}\text{?}}\equiv \left[\begin{array}{c}C\\ \mathrm{CA}\\ {\mathrm{CA}}^{2}\\ \vdots \\ {\mathrm{CA}}^{K}\end{array}\right]\left[{x}^{1}{x}^{2}\dots \text{\hspace{1em}}{x}^{M1}\right]\\ =\left[\begin{array}{cccc}{y}_{c0}^{1}& {y}_{c0}^{2}& \dots & {y}_{c0}^{M1}\\ {y}_{\mathrm{c1}}^{1}& {y}_{\mathrm{c1}}^{2}& \dots & {y}_{\mathrm{c1}}^{M1}\\ \dots & \dots & \dots & \dots \\ {y}_{\mathrm{cK}}^{1}& {y}_{\mathrm{cK}}^{2}& \dots & {y}_{\mathrm{cK}}^{M1}\end{array}\right]\end{array}& \left(32\right)\\ \begin{array}{c}{H}_{c\text{\hspace{1em}}\text{?}}\equiv \left[\begin{array}{c}C\\ \mathrm{CA}\\ {\mathrm{CA}}^{2}\\ \vdots \\ {\mathrm{CA}}^{K}\end{array}\right]A\left[{x}^{1}{x}^{2}\dots \text{\hspace{1em}}{x}^{M1}\right]\\ =\left[\begin{array}{cccc}{y}_{\mathrm{c1}}^{1}& {y}_{\mathrm{c1}}^{2}& \dots & {y}_{\mathrm{c1}}^{M1}\\ {y}_{\mathrm{c2}}^{1}& {y}_{\mathrm{c2}}^{2}& \dots & {y}_{\mathrm{c2}}^{M1}\\ \dots & \dots & \dots & \dots \\ {y}_{\mathrm{cK}+1}^{1}& {y}_{\mathrm{cK}+1}^{2}& \dots & {y}_{\mathrm{cK}+1}^{M1}\end{array}\right]\end{array}\text{}\mathrm{where}& \left(33\right)\\ \begin{array}{c}{y}_{\mathrm{ck}}^{n}={\mathrm{CA}}^{k}x\text{?}\\ ={y}^{n+k}\sum _{\text{?}=\text{?}}^{\text{?}}\text{\hspace{1em}}y\text{?}\sum _{\text{?}=\text{?}}^{\text{?}}\text{\hspace{1em}}y\text{?}\dots \text{\hspace{1em}}\\ \sum _{\text{?}=\text{?}}^{\text{?}}\text{\hspace{1em}}y\text{?}\\ {y}^{n+k}\sum _{\text{?}=\text{?}}^{\text{?}}\text{\hspace{1em}}\sum _{\text{?}=\text{?}}^{\text{?}}\text{\hspace{1em}}y\text{?}\end{array}\text{}\text{?}\text{indicates text missing or illegible when filed}& \left(34\right)\end{array}$

 SVD of the new H_{c0 }leads to
$\begin{array}{cc}\begin{array}{c}{H}_{\mathrm{c01}}\equiv {U}_{1}\sum _{i}\text{\hspace{1em}}{V}_{1}^{T}\\ \simeq \left[{U}_{1R}{U}_{1D}\right]\left[\begin{array}{cc}\sum _{1R}\text{\hspace{1em}}& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}{V}_{1R}^{T}\\ {V}_{1D}^{T}\end{array}\right]\\ ={U}_{1R}\sum _{1R}^{1/2}\text{\hspace{1em}}\sum _{1R}^{1/2}{V}_{1R}^{T}\end{array}& \left(35\right)\end{array}$
 from which we obtain
$\begin{array}{cc}D=\mathrm{the}\text{\hspace{1em}}\mathrm{first}\text{\hspace{1em}}N\text{?}\text{\hspace{1em}}\mathrm{rows}\text{\hspace{1em}}\mathrm{of}\text{\hspace{1em}}{Y}^{0}& \left(36\right)\\ C\simeq \mathrm{the}\text{\hspace{1em}}\mathrm{first}\text{\hspace{1em}}N\text{?}\text{\hspace{1em}}\mathrm{rows}\text{\hspace{1em}}\mathrm{of}\text{\hspace{1em}}{U}_{1R}\sum _{1R}^{1/2}\text{\hspace{1em}}& \left(37\right)\\ B\simeq \sum _{1R}^{1/2}\text{\hspace{1em}}{U}_{1R}^{T}{Y}_{c}^{1}& \left(38\right)\\ A\simeq \sum _{1R}^{1/2}\text{\hspace{1em}}{U}_{1R}^{T}{H}_{\mathrm{c11}}{V}_{1R}\sum _{1R}^{1/2}\text{\hspace{1em}}\text{}\mathrm{where}& \left(39\right)\\ {Y}_{c}^{n}\equiv \left[\begin{array}{cccc}{y}_{1}^{n}& {y}_{2}^{n}& \dots & {y}^{n}\text{?}\\ {y}_{1}^{n+1}& {y}_{2}^{n+1}& \dots & y\text{?}\\ \dots & \dots & \dots & \dots \\ {y}_{1}^{n+K+1}& {y}_{2}^{n+K+1}& \dots & {y}^{n+K+1}\text{?}\end{array}\right]\left(n=0.1\right)\text{}\text{?}\text{indicates text missing or illegible when filed}& \left(40\right)\end{array}$

[0087]
Using this alternate approach, the total number of measurements available is now (K+1)×N_{0}. It is noted that the additional time samples are required for the pulse response as well as for the response due to the SCI. More specifically, if K blocks of outputs are to be added pulse responses due to each input must be sampled at K additional time steps in addition the first two time steps,
$\text{?}=0\text{\hspace{1em}}\mathrm{and}\text{\hspace{1em}}\text{?}=\Delta \text{\hspace{1em}}\text{?}.\text{}\text{?}\text{indicates text missing or illegible when filed}$
Also, the response due to the SCI must be sampled at K additional steps beyond the Mth step. The total number of samples to be taken is thus M+1+K+(2+K)×N_{i}. K should desirably be sufficiently large enough to satisfy the measurement requirement, (K+1)×N_{0}≧R. Needless to say, this alternate scheme requires extra computation time because of the additional time samples required in the data set.

[0088]
A few candidate signals for the SCI will now be described. For an ideal linear model any of these SCIs can be used and all of them should yield ROMs of essentially the same quality.

[0089]
As noted above, random signals may be used for construction of SCI. This approach may be referred to as random inputs SCI (RSCI). That is, use
$\begin{array}{cc}{r}_{i}^{n}\equiv a\text{\hspace{1em}}\mathrm{sequence}\text{\hspace{1em}}\mathrm{of}\text{\hspace{1em}}\mathrm{random}\text{\hspace{1em}}\mathrm{numbers}& \left(41\right)\end{array}$
in Equations (20) and (21) above.

[0090]
Alternately, a filtered inputs based SCI (FSCI) may be employed. More specifically, to inject smooth inputs, one can filter the random signals through a lowpass filter. That is,
$\begin{array}{cc}\begin{array}{c}\text{?}=\text{?}\\ \equiv a\text{\hspace{1em}}\mathrm{sequence}\text{\hspace{1em}}\mathrm{of}\text{\hspace{1em}}\mathrm{filtered}\text{\hspace{1em}}\mathrm{random}\text{\hspace{1em}}\mathrm{numbers}\end{array}\text{}\text{?}\text{indicates text missing or illegible when filed}& \left(42\right)\end{array}$

[0091]
Using such low frequency signals may allow better convergence when applying the SCI to CFD models. Furthermore, since the frequency content is limited, it is possible to generate a smaller ROM directly from the SCI/ERA without any further reduction. A potential drawback is that if the filtered signals become too narrowly banded, they may not be as uncorrelated as desired. However, based on the theory of ergodicity [9], the statistical independence could be fortified by using longer signals and sampling the response for a longer period of time. (see, e.g., Papoulis, A., Probability, Random Variables, and Stochastic Processes, McGrawHill Book Company, New York, N.Y., 1982, incorporated herein by reference).

[0092]
In another approach, a step input based SCI (SCCI) may be used. In analogy to a single step input, one can apply multiple step inputs in a sequential manner:
$\begin{array}{cc}\begin{array}{c}\text{?}=\text{?}\\ \equiv a\text{\hspace{1em}}\mathrm{step}\text{\hspace{1em}}\mathrm{input}\text{\hspace{1em}}\mathrm{applied}\text{\hspace{1em}}\mathrm{at}\text{\hspace{1em}}\text{?}\text{?}\text{\hspace{1em}}\mathrm{step}\end{array}\text{}\mathrm{where}& \left(43\right)\\ \text{?}\equiv \left\{\begin{array}{cc}0& \left(n=1,2,\dots \text{\hspace{1em}},\text{?}1\right)\\ 1& \left(n\ge \text{?}\right)\end{array}\right\}\text{}\text{?}\text{indicates text missing or illegible when filed}& \left(44\right)\end{array}$

[0093]
To assure independence of the inputs, the onsets of the signals should be apart from each other by a sufficient number of steps, i.e., k_{i }must be large enough.

[0094]
Similarly, in another alternate approach, a pulse input based SCI (PSCI) may be used. In this approach, one can also apply multiple pulse inputs in a sequential manner:
$\begin{array}{cc}\begin{array}{c}\text{?}=\text{?}\\ \equiv a\text{\hspace{1em}}\mathrm{step}\text{\hspace{1em}}\mathrm{input}\text{\hspace{1em}}\mathrm{applied}\text{\hspace{1em}}\mathrm{at}\text{\hspace{1em}}\text{?}\text{?}\text{\hspace{1em}}\mathrm{step}\end{array}\text{}\mathrm{where}& \left(45\right)\\ \text{?}\equiv \left\{\begin{array}{cc}1& \left(n=\text{?}\right)\\ 0& \left(n=\mathrm{all}\text{\hspace{1em}}\mathrm{other}\text{\hspace{1em}}\mathrm{points}\right)\end{array}\right\}\text{}\text{?}\text{indicates text missing or illegible when filed}& \left(46\right)\end{array}$

[0095]
Again, k_{i }should be large enough to ensure independency of the applied signals.

[0096]
Yet another alternate embodiment to the SCI/ERA methods discussed above will now be described. In applications of discretetime computational models, there exist two conflicting requirements for the incremental time step Δt. For numerical convergence, one should adopt a sufficiently small Δt_{1}. On the other hand, given the highest frequency of interest, ω_{c, }the Nyquist criterion requires
$\Delta \text{\hspace{1em}}{t}_{2}\le \frac{\text{?}}{{\omega}_{c}}.\text{}\text{?}\text{indicates text missing or illegible when filed}$
Usually, for practical purposes
${\mathrm{\Delta z}}_{2}\u2aa2{\mathrm{\Delta z}}_{1}.$

[0097]
For instance, in a structural model that is coupled with a CFD model for aeroelastic applications, the highest mode usually has a much lower natural frequency than the highest frequency content in the aerodynamic model. If the signals used in the SCI/ERA methods are sharp as in the random, step, or pulse inputs, the SCI will excite all the system dynamics and hence this characteristic will be carried over to the ERA based reducedorder model. As a result, the ERA reducedorder model (ROM) may be prone to have a large dimension to span the same high frequency range as the original fullorder model (FOM), which suggests that there may be room for further order reduction in the ROM.

[0098]
To perform a second order reduction, one can apply the FrequencyDomain KarhunenLoeve (FDKL) method to the ERA ROM defined by matrices, (28)(31), wherein frequency samples of the system within the given frequency range, (−ω
_{c}, ω
_{c}) are used to extract optimal modes, and a new reducedorder model is constructed via the Galerkin's approximation. (see, e.g., Kim, T.,
Discrete
Time Eigen Analysis and Optimal Model Reduction for Flutter &
Aeroelastic Damping/Frequency Prediction Based On CFL3
D
ERA, BADVTECHLLLM02013, BCAG, February 2003, incorporated herein by reference). In this embodiment, the optimal or so called KL modes Ø
_{i}, are the eigenmodes of the covariance matrix K:
$\begin{array}{cc}K\text{\hspace{1em}}\varphi \text{?}={\lambda}_{i}\varphi \text{?}\text{}\mathrm{where}& \left(47\right)\\ K\text{?}\equiv X\left({\omega}_{i}\right)X\left({\omega}_{j}\right)\text{?}& \left(48\right)\\ \begin{array}{c}{\omega}_{\mathrm{ij}}\equiv \mathrm{sampling}\text{\hspace{1em}}\mathrm{frequencies}\\ \left[{\omega}_{1}\text{\hspace{1em}}{\omega}_{2}\text{\hspace{1em}}\dots \text{\hspace{1em}}\omega \text{?}\right]\end{array}\text{}\text{?}\text{indicates text missing or illegible when filed}& \left(49\right)\end{array}$

 where ω_{1}=−ω_{c }and ω_{M}=ω_{c}. X(ω_{i}) are the frequency solutions of the ERA ROM subjected to the singlecompositeinput described by (20) and (21) except that it is prescribed in the frequency domain. Once the optimal modes are obtained, assume
$\begin{array}{cc}x\simeq \Phi \text{\hspace{1em}}p\text{\hspace{1em}}\mathrm{where}& \left(50\right)\\ \Phi \equiv \uf603{\varphi}_{1}{\varphi}_{2}\dots \text{\hspace{1em}}\varphi \text{?}\uf604& \left(51\right)\\ p\equiv \left\{\begin{array}{c}{p}_{1}\\ {p}_{2}\\ p\text{?}\end{array}\right\}\text{}\text{?}\text{indicates text missing or illegible when filed}& \left(52\right)\end{array}$

[0100]
R_{1 }is set to be equal to the rank of the covariance matrix which is usually smaller than R.

[0101]
After inserting (50) into (1) and (2) with the ERA ROM matrices, premultiplying both sides by Φ^{T }yields a new reducedorder model as
p ^{n+1} =A _{1} p ^{n} +B _{i} u ^{n} (53)
Y ^{n} =C _{1} p ^{n} +D _{u} ^{n} (54)
where
A_{1≡Φ} ^{T}AΦ (55)
B_{1}≡Φ^{T}B (56)
C_{1}≡CΦ (57)

[0102]
The dimension of the new model is (R_{1}×R_{1}).

[0103]
Application of various embodiments of methods and systems to representative, largescaled CFD models will now be described. Unlike the general system described by Equations (1) and (2), an unsteady fluid dynamic system is driven by displacement and velocity of its moving boundary surface simultaneously as they both contribute to the normal downwash on the surface. If one considers a statically nonlinear, dynamically linearized flow field, the unsteady fluid motion can be described as
x ^{n+1} =Ax ^{n} +B _{0} u ^{n} +B _{1} {dot over (u)} ^{n} (58)
y ^{n} =q(Cx ^{n} +D _{0} u ^{n} +D _{1} {dot over (u)} ^{n}) (59)
where
x≡(L×1) fluid states (60)
u≡(N_{i}×1) generalized displacements (61)
{dot over (u)}≡(N_{i}×1) generalized velocities (62)
y≡(N_{i}×1) generalized aerodynamic forces (63)
q≡dynamic pressure (64)

[0104]
It is noted that the above equations progress in nondimensional time,
$\tau \equiv \frac{V\text{?}}{\text{?}},\text{}\text{?}\text{indicates text missing or illegible when filed}$
rather than in the real time t and ({dot over ( )}) is the first derivative with respect to T. In fact, the dependency of the normal downwash on air speed disappears when the equations are discretized in T, as in Equations (58) and (59). The structural degrees of freedom, u_{i }are the generalized coordinates associated with structural modes. These modes are used to describe the motion of the lifting surface. Two different types of reducedorder fluid dynamic models can be obtained depending on how the inputs are treated. If necessary, the FDKL/SCI can be performed for additional reduction.

[0105]
In one particular embodiment, a method in accordance with the invention may be applied to an aerodynamic ROM with displacement and velocity inputs. In this embodiment, one can treat u
^{n }and {dot over (u)}
^{n }separately as independent inputs. Thus, for the pulse inputs
$\begin{array}{cc}u\text{?}=\stackrel{.}{u}\text{?}=\delta \text{?}\equiv \left\{\begin{array}{cc}1& \left(n=0\right)\\ 0& \left(n=1,2,\dots \text{\hspace{1em}},M\right)\end{array}\right\}\text{}\text{?}\text{indicates text missing or illegible when filed}& \left(65\right)\end{array}$

 for i=1, 2, . . . , N_{i}, we obtain the corresponding responses
${y}_{i}^{0},{y}_{i}^{1}$
at the first two time steps. Let us define
$\begin{array}{cc}{Y}^{0}=\left[{y}_{1}^{2}{y}_{2}^{0}\dots \text{\hspace{1em}}{y}_{{N}_{i}}^{0}{y}_{{N}_{i}+1}^{0}{y}_{{N}_{i}+2}^{0}\dots \text{\hspace{1em}}{y}_{2{N}_{i}}^{0}\right]& \left(66\right)\\ {Y}^{1}=\left[{y}_{1}^{1}{y}_{2}^{1}\dots \text{\hspace{1em}}{y}_{{N}_{i}}^{1}{y}_{{N}_{i}+1}^{1}{y}_{{N}_{i}+2}^{1}\dots \text{\hspace{1em}}{y}_{2{N}_{i}}^{1}\right]& \left(67\right)\end{array}$
 where the first N_{i }samples are due to the pulses in u^{n }and the next N_{i }ones are due to the pulses in {dot over (u)}^{n}. Next, we prepare the following inputs,
$\begin{array}{cc}{b}_{\mathrm{SCI}}^{n}\equiv \sum _{\text{?}}^{{N}_{i}}{b}_{0\stackrel{.}{u}}\text{?}+\sum _{\stackrel{.}{u}=1}^{{N}_{i}}{b}_{1\stackrel{.}{u}}\text{?}& \left(68\right)\\ {d}_{\mathrm{SCI}}^{n}\equiv \sum _{\text{?}}^{{N}_{i}}{d}_{0\stackrel{.}{u}}\text{?}+\sum _{\stackrel{.}{u}=1}^{{N}_{i}}{d}_{1\stackrel{.}{u}}\text{?}\text{}\text{?}\text{indicates text missing or illegible when filed}& \left(69\right)\end{array}$

[0108]
Subject to the SCI, we sample the system response y
^{n }and get
$\begin{array}{cc}{y}_{\mathrm{c0}}^{n}\equiv {y}^{n}\sum _{\stackrel{.}{u}=1}^{2{N}_{i}}{y}_{\stackrel{.}{u}}^{0}r\text{?}& \left(70\right)\\ {y}_{\mathrm{c1}}^{n}\equiv {y}^{n+1}\sum _{\stackrel{.}{u}=1}^{2{N}_{i}}{y}_{\stackrel{.}{u}}^{0}r\text{?}\sum _{\stackrel{.}{u}=1}^{2{N}_{i}}{y}_{\stackrel{.}{u}}^{0}r\text{?}\text{}\mathrm{Define}& \left(71\right)\\ {H}_{\mathrm{c0}}\equiv \uf603\begin{array}{cccc}{y}_{\mathrm{c0}}^{1}& {y}_{\mathrm{c0}}^{2}& \dots & {y}_{\mathrm{c0}}^{M1}\end{array}\uf604& \left(72\right)\\ {H}_{\mathrm{c1}}\equiv \uf603\begin{array}{cccc}{y}_{\mathrm{c1}}^{1}& {y}_{\mathrm{c1}}^{2}& \dots & {y}_{\mathrm{c1}}^{M1}\end{array}\uf604\text{}\text{?}\text{indicates text missing or illegible when filed}& \left(73\right)\end{array}$

 where SVD of H_{c0 }yields
$\begin{array}{cc}{H}_{\mathrm{c0}}\equiv U\sum {V}^{T}\text{}\text{\hspace{1em}}\simeq \uf603\begin{array}{cc}{U}_{R}& {U}_{D}\end{array}\uf604\left[\begin{array}{cc}\sum _{R}& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}{V}_{R}^{T}\\ {V}_{D}^{T}\end{array}\right]\text{}\text{\hspace{1em}}={U}_{R}\underset{R}{\sum ^{1/2}}\underset{R}{\sum ^{1/2}}{V}_{R}^{T}& \left(74\right)\end{array}$
 with R≡rank(H_{c0}). Hence, the reducedorder model is given by
$\begin{array}{cc}{D}_{0}=\mathrm{the}\text{\hspace{1em}}\mathrm{first}\text{\hspace{1em}}{N}_{i}\text{\hspace{1em}}\mathrm{columns}\text{\hspace{1em}}\mathrm{of}\text{\hspace{1em}}{Y}^{0}& \left(75\right)\\ {D}_{1}=\mathrm{the}\text{\hspace{1em}}\mathrm{second}\text{\hspace{1em}}{N}_{i}\text{\hspace{1em}}\mathrm{columns}\text{\hspace{1em}}\mathrm{of}\text{\hspace{1em}}{Y}^{0}& \left(76\right)\\ C\simeq {U}_{R}\underset{R}{\sum ^{1/2}}& \left(77\right)\\ {B}_{0}\simeq \mathrm{the}\text{\hspace{1em}}\mathrm{first}\text{\hspace{1em}}{N}_{i}\text{\hspace{1em}}\mathrm{columns}\text{\hspace{1em}}\mathrm{of}\text{\hspace{1em}}\underset{R}{\sum ^{1/2}}{U}_{R}^{T}{Y}^{1}& \left(78\right)\\ {B}_{1}\simeq \mathrm{the}\text{\hspace{1em}}\mathrm{second}\text{\hspace{1em}}{N}_{i}\text{\hspace{1em}}\mathrm{columns}\text{\hspace{1em}}\mathrm{of}\text{\hspace{1em}}\underset{R}{\sum ^{1/2}}{U}_{R}^{T}{Y}^{1}& \left(79\right)\\ A\simeq \underset{R}{\sum ^{1/2}}{U}_{R}^{T}{H}_{\mathrm{c1}}{V}_{R}\underset{R}{\sum ^{1/2}}& \left(80\right)\end{array}$

[0111]
In yet another particular embodiment, a method in accordance with the invention may be applied to an aerodynamic ROM with only the displacements as the system inputs. This is possible by applying simultaneously the pulse and double pulse inputs,
$\begin{array}{cc}{u}_{i}^{n}={\delta}^{n}\equiv \left\{\begin{array}{cc}1& \left(n=0\right)\\ 0& \left(n=1,2,\dots \text{\hspace{1em}},M\right)\end{array}\right\}& \left(81\right)\\ {\stackrel{.}{u}}_{i}^{n}={\delta}^{n}\equiv \left\{\begin{array}{cc}\frac{1}{\Delta \text{\hspace{1em}}z}& \left(n=0\right)\\ \frac{1}{\mathrm{\Delta \tau}}& \left(n=1\right)\\ 0& \left(n=2,3,\dots \text{\hspace{1em}},M\right)\end{array}\right\}& \left(82\right)\end{array}$

 and get the corresponding responses
${y}_{d\text{\hspace{1em}}i}^{0}{y}_{d\text{\hspace{1em}}i}^{1}$
at the first two time steps:
$\begin{array}{cc}{Y}_{d}^{0}=\left[{y}_{d\text{\hspace{1em}}1}^{0}\text{\hspace{1em}}{y}_{d\text{\hspace{1em}}2}^{0}\text{\hspace{1em}}\dots \text{\hspace{1em}}{y}_{{\mathrm{dN}}_{i}}^{0}\right]& \left(83\right)\\ {Y}_{d}^{1}=\left[{y}_{d\text{\hspace{1em}}1}^{1}\text{\hspace{1em}}{y}_{d\text{\hspace{1em}}2}^{1}\text{\hspace{1em}}\dots \text{\hspace{1em}}{y}_{{\mathrm{dN}}_{i}}^{1}\right]& \left(84\right)\end{array}$

[0113]
For a new SCI, we use
$\begin{array}{cc}{b}_{\mathrm{SCI}}\text{?}\equiv \sum _{\stackrel{.}{u}=1}^{{N}_{i}}{b}_{0\stackrel{.}{u}}\text{?}+\sum _{\stackrel{.}{u}=1}^{{N}_{i}}{b}_{1\stackrel{.}{u}}\text{?}& \left(85\right)\\ {d}_{\mathrm{SCI}}\text{?}\equiv \sum _{\stackrel{.}{u}=1}^{{N}_{i}}{d}_{0\stackrel{.}{u}}\text{?}+\sum _{\stackrel{.}{u}=1}^{{N}_{i}}{d}_{1\stackrel{.}{u}}\text{?}\text{}\text{?}\text{indicates text missing or illegible when filed}& \left(86\right)\end{array}$

 where
${i}_{i}^{n}$
are discretetime derivative of
${i}_{i}^{n}.$
To be consistent with the use of the double pulse defined in (82),
${i}_{i}^{n}$
are obtained by filtering
${i}_{i}^{n}\text{\hspace{1em}}\mathrm{via}\text{\hspace{1em}}{\delta}_{i}^{n},i.e.,$
$\begin{array}{cc}{i}_{i}^{n}\equiv \mathrm{conv}\left({i}_{i}^{k},{\delta}_{i}^{k}\right)& \left(87\right)\end{array}$
 which is equivalent to the backward difference equation,
$\begin{array}{cc}{r}_{\stackrel{.}{u}}^{n}\equiv \frac{r\text{?}r\text{?}}{\Delta \text{\hspace{1em}}\tau}\text{}\text{?}\text{indicates text missing or illegible when filed}& \left(88\right)\end{array}$

[0116]
Subject to the new SCI we sample the system response y
^{n }and get
$\begin{array}{cc}{y}_{\mathrm{dc0}}^{n}\equiv {y}^{n}\sum _{\stackrel{.}{u}=1}^{{N}_{i}}{y}_{d\stackrel{.}{u}}^{0}{r}_{\stackrel{.}{u}}^{n}& \left(89\right)\\ {y}_{\mathrm{dc1}}^{n}\equiv {y}^{n+1}\sum _{\stackrel{.}{u}=1}^{{N}_{i}}{y}_{d\stackrel{.}{u}}^{0}{r}_{\stackrel{.}{u}}^{n+1}\sum _{\stackrel{.}{u}=1}^{{N}_{i}}{y}_{d\stackrel{.}{u}}^{1}{r}_{\stackrel{.}{u}}^{n+1}\text{}\mathrm{Defining}& \left(90\right)\\ {H}_{\mathrm{dc0}}\equiv \uf603\begin{array}{cccc}{y}_{\mathrm{dc0}}^{1}& {y}_{\mathrm{dc0}}^{2}& \dots & {y}_{\mathrm{dc0}}^{M1}\end{array}\uf604& \left(91\right)\\ {H}_{\mathrm{dc1}}\equiv \uf603\begin{array}{cccc}{y}_{\mathrm{dc1}}^{1}& {y}_{\mathrm{dc1}}^{2}& \dots & {y}_{\mathrm{dc1}}^{M1}\end{array}\uf604& \left(92\right)\end{array}$

 the SVD of H_{dc0 }yields
$\begin{array}{cc}{H}_{\mathrm{dc0}}\equiv U\sum {V}^{T}\text{}\text{\hspace{1em}}\simeq \uf603\begin{array}{cc}{U}_{R}& {U}_{D}\end{array}\uf604\left[\begin{array}{cc}\sum _{R}& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}{V}_{R}^{T}\\ {V}_{D}^{T}\end{array}\right]\text{}\text{\hspace{1em}}={U}_{R}\underset{R}{\sum ^{1/2}}\underset{R}{\sum ^{1/2}}{V}_{R}^{T}& \left(93\right)\end{array}$
with R≡rank(H_{dc0}). The new reducedorder model has only N_{i }input channels and is in the form
$\begin{array}{cc}x\text{?}={\mathrm{Ax}}^{n}+{\mathrm{Bu}}^{n}& \left(94\right)\\ {y}^{n}=q\left({\mathrm{Cx}}^{n}+{\mathrm{Du}}^{n}\right)\text{\hspace{1em}}\mathrm{where}& \left(95\right)\\ D={Y}^{0}\text{?}& \left(96\right)\\ C\simeq {U}_{R}{\Sigma}_{R}^{1/2}& \left(97\right)\\ B\simeq {\Sigma}_{R}^{1/2}{U}_{R}^{T}{Y}^{1}\text{?}& \left(98\right)\\ A\simeq {\Sigma}_{R}^{1/2}{U}_{R}^{T}H\text{?}{V}_{R}{\Sigma}_{R}^{1/2}\text{}\text{?}\text{indicates text missing or illegible when filed}& \left(99\right)\end{array}$

[0118]
Embodiments of aeroelastic systems in accordance with the present invention can be constructed using the embodiments of reducedorder aerodynamic models described above. For example, in one embodiment, we first note that structural model is normally described in real, continuoustime:
mü+c{dot over (u)}+ku=y (100)
({dot over ( )}) and ({umlaut over ( )}) respectively, represent the first and the second derivatives with respect to t. Hence, to construct aeroelastic model the continuoustime equation (100) is discretized in time:
$\begin{array}{cc}{z}^{n+1}={A}_{s}{z}^{n}+{B}_{s}{y}^{n}& \left(101\right)\\ {u}^{n}=C\text{?}{z}^{n}\text{\hspace{1em}}\mathrm{where}& \left(102\right)\\ z\equiv \left\{\begin{array}{c}u\\ \stackrel{.}{u}\end{array}\right\}& \left(103\right)\\ C\text{?}\equiv \left[\begin{array}{cc}I& 0\end{array}\right]\text{}\text{?}\text{indicates text missing or illegible when filed}& \left(104\right)\end{array}$
Note that given ΔT and V the consistent incremental time step
$\Delta \text{?}=\frac{\Delta \text{?}}{V},\text{}\text{?}\text{indicates text missing or illegible when filed}$
must be used in the conversion to the discretetime.
Aeroelastic Model I

[0119]
In this approach, we treat the displacement and velocity, u^{n}, {dot over (u)}^{n }as independent inputs and apply the b_{SCi} ^{n}, d_{SCI} ^{n }given by Eqs. (68) and (69) at a reference dynamic pressure,
${\text{?}}_{\mathrm{ref}}\equiv \frac{1}{2}{\text{?}}_{\mathrm{ref}}{V}_{\mathrm{ref}}^{2}.\text{}\text{?}\text{indicates text missing or illegible when filed}$
The corresponding samples are taken and scaled by
$\frac{1}{{\text{?}}_{\mathrm{ref}}}.\text{}\text{?}\text{indicates text missing or illegible when filed}$
The ROM of the first kind described above is then obtained by applying the SCI/ERA. An aeroelastic system that is valid at all V can be constructed by combining the resulting aerodynamic ROM with the structural equations:
$\begin{array}{cc}{X}^{n+1}=A\text{?}{X}^{n}\text{\hspace{1em}}\mathrm{where}& \left(105\right)\\ X\equiv \left\{\begin{array}{c}x\\ z\end{array}\right\}& \left(106\right)\\ A\text{?}\equiv \lceil \begin{array}{cc}A& \uf603B\text{?}{B}_{1}\uf604\\ {\mathrm{qB}}_{5}C& {A}_{s}+\mathrm{qB}\text{?}\uf603D\text{?}{D}_{1}\uf604\end{array}]\text{}\text{?}\text{indicates text missing or illegible when filed}& \left(107\right)\end{array}$

[0120]
Denoting the eigenvalues of system matrix (107) as
${\lambda}_{\mathrm{dl},\stackrel{.}{z}:},$
the aeroelastic eigenvalues in the Laplace domain are obtained as
$\begin{array}{cc}\lambda =\frac{\mathrm{log}\left(\lambda \text{?}\right)}{\Delta \text{?}}\text{}\text{?}\text{indicates text missing or illegible when filed}& \left(108\right)\end{array}$

[0121]
For flutter instability, one must have
$\mathrm{Real}\left({.}^{\mathrm{\lambda dl},\stackrel{.}{z}:\text{\hspace{1em}}}\right)>0\text{\hspace{1em}}\mathrm{and}\text{\hspace{1em}}I~~{\uf605}^{{\lambda}_{\mathrm{dl},\stackrel{.}{z}:}}\uf606>1$
for any i.

[0122]
Aeroelastic Model II.

[0123]
One can also apply the
${b}_{\mathrm{SCI}:}^{n}{d}_{\mathrm{SCI}}^{n}$
given by Equations (85) and (86) and get the ROM of the second kind described above. This will produce aerodynamic system matrices, A, B, C, and D where B, D each has only N_{i }columns. The new resulting reducedorder aeroelastic model can be obtained as
$\begin{array}{cc}X\text{?}=A\text{?}X\text{?}\text{}\mathrm{where}& \left(109\right)\\ A\text{?}\equiv \left[\begin{array}{cc}A& \mathrm{BC}\text{?}\\ \mathrm{qB}\text{?}C& A\text{?}+\mathrm{qB}\text{?}\mathrm{DC}\text{?}\end{array}\right]\text{}\text{?}\text{indicates text missing or illegible when filed}& \left(110\right)\end{array}$

[0124]
Although q can change in Equation (110), this model must be used only at the reference air speed V_{ref}. Given a local angle of attack u, plunging rate {dot over (u)}, and local air speed V the total aerodynamic downwash at the moving boundary is {dot over (u)}+V u and as such it is impossible to separate out and account for the effect of V without having both the displacement and velocity channels. However, this drawback is easily remedied by adjusting the incremental time step according to
$\Delta \text{?}=\frac{\Delta \text{?}}{V}$
$\text{?}\text{indicates text missing or illegible when filed}$
when the air speed changes from one value to another and discretizing the structural model based on the new Δt. That is, if one leaves the V dependency in the structure,
$\begin{array}{cc}A\text{?}\equiv \left[\begin{array}{cc}A& \mathrm{BC}\text{?}\\ \mathrm{qB}\text{?}\left(V\right)C& \mathrm{qB}\text{?}\left(V\right)\mathrm{DC}\text{?}\end{array}\right]\text{}\text{?}\text{indicates text missing or illegible when filed}& \left(111\right)\end{array}$
then the Aeroelastic Model II becomes valid for all air speeds.

[0125]
To demonstrate an embodiment of the present invention, an unsteady vortex lattice aerodynamic model subjected to several structural mode inputs may be considered. For example, FIG. 3 is a schematic view of a vortex lattice formulation 300 for computing an unsteady flow field around a flat, rectangular wing 302 in accordance with an embodiment of the present invention. In this embodiment, the unsteady flow field is modeled as an incompressible, subsonic air flow. In this example, the wing 302 is 3 inches wide and 12 inches long, has 10 and 20 vortex elements in the chordwise and spanwise directions, respectively. A free wake 304 has 40 and 20 vortex elements in the streamwise and spanwise directions, creating a total of 800 degrees of freedom. (see Kim, T., Nam, C., and Kim, Y., 1997, ReducedOrder Aeroservoelastic Model with an Unsteady Aerodynamic Eigen Formulation, AIAA Journal, Vol. 35, No. 6, pp. 10871088, incorporated herein by reference). The wing structure is modeled using 6 vibrational (3 bending and 3 torsional) modes (see Crawley, E. F., and Dugundji, J., 1980, Frequency Determination and NonDimensionalization for Composite Cantilever Plates, Journal Sound and Vibration, Vol. 72, No. 1, pp. 110, incorporated herein by reference).

[0126]
No structural damping was introduced in this example. Thus, the size of the fullorder aeroelastic model is (812×812).

[0127]
For the Aeroelastic Model I the reference air density and speed 306 were set at 1.23 kg/r^{3}, 80 m/sec, respectively. The incremental time at this reference speed is
$\Delta \text{?}=\frac{\mathrm{dz}}{\text{?}}$
$\text{?}\text{indicates text missing or illegible when filed}$
9.525×10^{−5 }sec. For the sampling of the vortex model, 480 extra outputs were extracted in addition to the 6 generalized aerodynamic forces at 481 time steps. Applying 12 sets of random signal inputs simultaneously, 6 for u^{n}, 6 for {dot over (u)}^{n}, yielded a single set of sampled data. 12 pulse inputs were also applied individually at the first two time steps to generate Y^{o }and Y^{1}.

[0128]
FIG. 4 shows three sets of random numbers 402, 404, 406 generated as inputs to the vortex lattice formulation of FIG. 3 using the wellknown MATLAB computer program. Out of 481 time samples, the SVD produced 413 linearly independent singular modes. This number was determined by the rank of H_{c0 }matrix. Thus, the size of the aerodynamic ROM became (413×413). The reducedorder aerodynamic model was then coupled with the structural model to create (425×425) aeroelastic model (ROM I.).

[0129]
For the Aeroelastic Model II. the reference air density and speed were again set at 1.23 kg/m^{3}, and 80 m/sec. 6 sets of random signals and 6 sets of discretetime derivatives of the random signals were applied for u^{n }and {dot over (u)}^{n }using 481 time steps and the 486 output measurements. This yielded (329×329) aerodynamic ROM which when combined with the structural system, produced (341×341) aeroelastic model (ROM II.). It is noted that ROM II. is approximately 20% smaller than ROM I. as a result of using only the half of the input channels.

[0130]
Next, the dimensions of the reducedorder aerodynamic models were further decreased using the FDKL/SCI method. As mentioned earlier, the incremental time step embedded in both the FOM and SCI/ERA ROM is too small to be effective for various aeroelastic simulations which usually involve a low frequency range. Considering that the highest free vibrational frequency of the structural modes is 4,160 rad/sec, the sampling range in the FDKL method was restricted to (−4,500, 4,500) rad/sec. For the ROM I., out of 174 frequency samples within the range 129 KL modes were selected based on the rank of the covariance matrix K Hence, the size of the new reducedorder aerodynamic and aeroelastic models (ROM I.FDKL) became 129 and 141, respectively. Likewise, for the ROM II. 97 KL modes out of 130 frequency samples in the same frequency range were selected yielding new (109×109) aeroelastic model (ROM II.FDKL). For computational efficiency, these reducedorder models are to be preferred over the ROM I. and ROM II.

[0131]
FIG. 5 shows three sets of generalized aerodynamic forces 500 computed using the vortex lattice formulation of FIG. 3. Specifically, FIG. 5 shows (6×6) generalized aerodynamic forces for V=80 m/sec obtained from the FOM, ROM I.FDKL, and ROM II.FDKL models (FOM (800), ROM I.FDKL (129), and ROM II.FDKL (97)), in the nondimensional time, T. It is seen that despite the cutoff frequency range present in the latter two models, they reproduce the pulse aerodynamic responses of the original model very well.

[0132]
FIGS. 6 and 7 show two different scales of the aeroelastic eigenvalues 600, 700 of the various models in the s domain at V=80 m/sec. It is seen that many eigenvalues of the reducedorder aeroelastic models match very well with those of the full model (FIG. 6). More specifically, the 12 complex eigenvalues associated with 6 structural modes agree very well between the FOM and ROMs, although the higher structural modes (5th and 6th) in the ROM II. and ROM II.FDKL are slightly mismatched (FIG. 7). Also noteworthy is that all the eigenvalues of the ROM I.FDKL and ROM II.FDKL are approximately within the specified bound, (−4,500, 4,500) rad/sec as the model was obtained based on frequency samples in the range.

[0133]
FIGS. 8 and 9 show time responses of the first two structural modes 800, 900 due to an initial condition in ü_{1}. It can be seen that the three sets of curves are practically indistinguishable from each other.

[0134]
Similarly, FIGS. 1013 show the aeroelastic results at V=121.2 m/sec. Specifically, FIG. 10 shows a set of aeroelastic eigenvalues 1000 computed using the vortex lattice formulation of FIG. 3 at a first scale for V=121.2 m/sec. FIG. 11 shows a set of aeroelastic eigenvalues 1100 computed using the vortex lattice formulation of FIG. 3 at a second scale for V=121.2 m/sec. FIG. 12 shows a time response 1200 of the first structural mode computed using the vortex lattice formulation of FIG. 3 for V=121.2 m/sec. And FIG. 13 shows a time response 1300 of the second structural mode computed using the vortex lattice formulation of FIG. 3 for V=121.2 m/sec. As can be seen from these figures, the wing is on the verge of flutter at this speed. It may be noted how accurately the ROM I.FDKL is able to reproduce neutrally stable, sinusoidal time responses (FIGS. 1213). However, the ROM II.FDKL exhibits a noticeable but minor error in producing the transient response.

[0135]
FIG. 14 shows a graph 1400 of model construction time of VLM ERA ROMs versus number of inputs for the vortex lattice formulation of FIG. 3. Using FIG. 14, the model construction time may be compared between the two ERA methods. In order to obtain accurate and consistent singular modes, H_{c0}, H_{c1 }matrices were kept as square as possible by keeping the number of time samples approximately equal to the total number of measurements which is the sum of the number of generalized aerodynamic forces and the number of auxiliary measurements. The same numbers of time steps and auxiliary outputs were used both in the Pulse/ERA and Sd/ERA. Thus, in the first case where only the first bending mode alone excites the flow field, M=131, N_{0}=131. In the second case where the first bending and first torsional modes were included, M=251, N_{0}=252, and in the third case where the first bending and torsional as well as the second bending modes excite the aerodynamic field, M=281, N_{0}=283. 4 and 5 inputs were also used with M=331, N_{0}=334 and M=411, N_{0}=415, respectively.

[0136]
FIG. 15 shows CPU seconds spent constructing various models in accordance with embodiments of the present invention (Note that FIG. 14 was obtained based on the data in FIG. 15). Specifically, Table 1 of FIG. 15 shows CPU seconds spent in constructing ROM I. on a SGI machine. Also presented in parenthesis are the dimensions of the corresponding reducedorder models. Table 2. shows CPU seconds consumed for ROM II. on the SGI machine. The numbers shown in FIGS. 14 and 15 represent total CPU seconds spent not only in sampling the response but also processing the data in the subsequent ERA schemes. As seen in FIGS. 14 and 15, the new method clearly has an advantage over the Pulse/ERA in reducing the model construction time yielding saving factors of multiple numbers. Needless to say, as the number of inputs increases so does the saving. It is interesting that for a given number of inputs both ERA methods generate ROMs of very similar sizes. As expected, ROM II of the SCI/ERA are as small as 80% of the corresponding ROM I. Despite the different input channels, both SCI/ERA ROM I. and ROM II. require approximately the same CPU time implying that in the SCI/ERA the overall computation time is not very sensitive to the number of inputs;

[0137]
In yet another embodiment of the present invention, the SCI/ERA method has also been applied to unsteady aerodynamic systems modeled by CFL3D code. CFL3D is a finite element program that based on the NavierStokes equations models nonlinear viscous, compressible fluid motion in subsonic as well as transonic flow fields. (see Krist, S. L., Biedron, R. T., and Rumsey, C. L., CFLSD User's Manual (Version 5.0), The NASA Langley Research Center, Hampton, Va.). Although CFL3D describes statically and dynamically nonlinear flow, when subjected to a small amplitude it predicts dynamically linearized behavior of the flow field around the nonlinear static position in the form of Equations (58) and (59). For most practical aeroelastic analyses such as flutter prediction and dynamic gust loads calculation, the small amplitude approximation yields sufficiently accurate results.

[0138]
The aeroelastic formulation for the CFD code is slightly different than for the vortex lattice case in that the second aerodynamic model discussed above is used but whenever the air speed V is changed the incremental time step Δt is adjusted accordingly when discretizing the structural model, Equation (100). As in the Aeroelastic Model I, the resulting model can account for the effect of changing the free stream speed.

[0139]
For example, FIG. 16 is a representative CFD simulation model 1600 in accordance with an embodiment of the present invention. In this embodiment, the airplane configuration under consideration is that of the TwinEngineTransportFlutterModel (TETFM). The aerodynamic grid is given by the so called the WingPencilNacelle (WPN) model with the strut between the wing and nacelle omitted. The structural motion is described by 10 structural modes, resulting in total of 10 generalized aerodynamic forces per a mode shape.

[0140]
The computational aerodynamic model consists of approximately 700,000 cells and 30 blocks. For detailed description of the modeling, see Hong, M. S., Bhatia, K. G., SenGupta, G., Kim, T., Kuruvila, G., Silva, W. A., Bartels, and R., Biedron, R., Simulations of a TwinEngine Transport Flutter Model In the Transonic Dynamics Tunnel, IFASD Paper 2003US44, incorporated herein by reference. Also, the details of the computational model and construction of the aerodynamic and aeroelastic ROMs based on the SCI/ERA using various types of input signals mentioned earlier is described more fully in Kim, T., Hong, M. S., Bhatia, K. G., SenGupta, G., Aeroelastic Model Reduction for an Affordable CFD Based Flutter Analysis, AIAA Paper 20042040, published subsequent to the filing of the present application and incorporated herein by reference.

[0141]
FIG. 17 shows four structural mode shapes 1700 used for the CFD simulation model 1600 of FIG. 16.

[0142]
FIG. 18 shows a Vg plot 1800 of the two different aeroelastic models in accordance with the previous and an embodiment of the present invention, respectively. The WPN model was examined at M=0.831. With Δt=3.34×10^{−4 }sec and 995 time steps, aerodynamic ROMs were created using both the Pulse/ERA and SSCI/ERA methods based on the displacement method described earlier (ROM II.). The size of the aerodynamic ROMs is (512×512). Given the 10 mode inputs, the total number of time steps consumed in the time marching required for the SCI/ERA process was 995+81×10=1805. On the other hand, the traditional ERA required total of 10×995=9,950. Using a single CPU of IBM/Regatta machine, the total CPU hours for the Pulse/ERA calculation was 366 while the SSCI/ERA required only 64 hours. Thus, the computational cost was reduced by a factor of 5.7. As shown in FIG. 18, the two ROMs match very well in flutter characteristics predicting the two flutter points within just 0.13% and 0.36% differences, respectively.

[0143]
In accordance with embodiments of the present invention, efficient timedomain model reduction/system identification techniques have been presented and demonstrated for linear dynamic systems that are subjected to multiple righthandside inputs. Methods and systems in accordance with the present invention do not require modifying the original code and take only input and output data for the model construction. Furthermore, methods in accordance with the present invention are based on a direct singular value decomposition of the output measurements that are not necessarily attributed to pulse inputs but due to multiple signal inputs applied simultaneously at the input channels. Compared to the Pulse/ERA, the SCI/ERA methods disclosed herein can significantly reduce the model construction time and compress the amount of output data. Therefore, such methods and systems are very attractive for largescaled dynamic systems with multiple driving inputs such as CFD models wherein the moving boundary input is often described by many structural modes.

[0144]
While preferred and alternate embodiments of the invention have been illustrated and described, as noted above, many changes can be made without departing from the spirit and scope of the invention. Accordingly, the scope of the invention is not limited by the disclosure of the preferred and alternate embodiments. Instead, the invention should be determined entirely by reference to the claims that follow.