International Journal of Rotating Machinery

Volume 2017 (2017), Article ID 7529716, 14 pages

https://doi.org/10.1155/2017/7529716

## Application of System Identification for Modeling the Dynamic Behavior of Axial Flow Compressor Dynamics

^{1}Department of Mechanical Engineering, Idaho State University, Mail Stop 8060, Pocatello, ID 83209, USA^{2}Institute of Engineering Thermophysics, CAS, 11 Beisihuanxi Road, Beijing 100190, China

Correspondence should be addressed to Marco P. Schoen

Received 1 February 2017; Accepted 11 April 2017; Published 7 May 2017

Academic Editor: Pietro Zunino

Copyright © 2017 Marco P. Schoen and Ji-Chao Lee. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

Identification of a one-stage axial compressor system is addressed. In particular, we investigate the underlying dynamics of tip air injection and throttle activation to the overall compressor dynamics and the dynamics around the tip of the compressor blades. A proposed subspace system identification algorithm is used to extract three mathematical models: relating the tip air injection to the overall dynamics of the compressor and to the flow dynamics at the tip of the compressor blade and relating the movement of the throttle to the overall compressor dynamics. As the system identification relays on experimental data, concerns about the noise level and unmodeled system dynamics are addressed by experimenting with two model structures. The identification algorithm entails a heuristic optimization that allows for inspection of the results with respect to unmodeled system dynamics. The results of the proposed system identification algorithm show that the assumed model structure for the system identification algorithm takes on an important role in defining the coupling characteristics. A new measure for the flow state in the blade passage is proposed and used in characterizing the dynamics at the tip of the compressor blade, which allows for the inspection of the limits for the utilized actuation.

#### 1. Introduction

Operation of axial compressor systems, as found in a number of aerospace applications, faces difficulties due to instabilities at the peak of its performance curves. To mitigate these instabilities, a number of different control mechanisms have been introduced. Among them are air injection into the tip of the blade passage [1–3] and the use of grooved casings [4–6]. Modeling of the overall compressor dynamics has been a challenging undertaking, resulting in the well-known Moore-Greitzer model. However, there does not exist a model relating the influence of air injection to the overall system dynamics. One of the aims of this work is to present a system identification based approach to characterize the influence of the air injection—commonly used to control the compressor—to the overall compressor dynamics.

As the system identification relays on experimental data, concerns about the noise level need to be addressed as well as unmodeled system dynamics. The unmodeled system dynamics can become an important issue when the experimental data do not include all system modes during the system identification experiment or the assumed model structure does not allow for accommodating one or more particular dynamic characteristics. We propose to utilize two different model structures to investigate the unmodeled system dynamics and the effects of the proposed optimization. These structures are the autoregressive moving average with exogenous input (ARMAX) model and the autoregressive with exogenous input (ARX) model. An innovation sequence is extracted by using an estimated parameter sequence. The ARMAX model and ARX model parameters are then utilized to construct a Hankel matrix, representing the sampled impulse response of the system relating tip air injection to pressure rise coefficient. A singular value decomposition of the Hankel matrix is performed for an eigensystem realization (ERA, [7, 8]) resulting in a state-space description of the coupling dynamics. However, prior to performing the state-space realization, a heuristic optimization algorithm is introduced to optimize the singular values with respect to the measured output data. This allows for minimization of the noise influence to the extracted state-space model as well as the effect the assumed model structure has on fitting to the recorded characteristics from the system identification experiment. The proposed system identification is used for data collected on an isolated rotor axial compressor system whose blade geometry allows for spike stall inception. The tip air injection is applied in a random fashion in order to allow for sufficient system mode excitation. The resulting pressure rise and associated flow coefficient are computed using a set of pressure sensors. The optimization of the Hankel matrix is accomplished by using a Tabu Search (TS) algorithm [9]. The TS searches the vicinity around each singular value defined by a percentage of the largest singular value. The eigensystem realization is done using a balanced realization, an input-normal form, and an output normal form realization in order to assess the impact of the optimized singular values. Another aim of this work is to characterize the dynamics within the blade passage at the tip of the rotor blade due to air injection at the leading edge of the rotor blade. Finally, the proposed hybrid system identification algorithm using TS is used to model the compressor dynamics excited by movements of the throttle and the resulting pressure rise changes.

In the following, the details of the system identification scheme as well as the proposed realization and optimization are introduced.

#### 2. System Identification Algorithm

Consider a linear, time-invariant, discrete time system:where , , , and are the system matrices and , , are the output and input data sequences recorded during a system experiment with some sample frequency , resulting in discrete data points for each of the two sequences. is the number of outputs, and is the number of inputs; is the process noise and is the measurement noise, while is the state vector.

In order to estimate the state variable , a Kalman filter is used, provided the states are observable [10]. Hence, the system given by (1) and (2) can be reformulated. Defining the Kalman filter gain as where is the solution of the steady-state algebraic Riccati equation. Hence,whereThe system in (1) is in the process form, while the system in (4) is in the innovation form. Using (6) in (4) results and defining and , we arrive at the predictor form of the given system:The predictor form of (8) can be used to derive an AutoRegressive with eXogenous input (ARX) model: Note, is an asymptotically stable square matrix; hence if is sufficiently large, we haveEquation (9) becomesNote that are the system Markov parameters and are the Kalman filter Markov parameters. Also, the existence of is guaranteed if the system is detectable and is stabilizable [7].

Equation (11) is of the formwhere and . Assuming the system is observable and controllable, the following output vector can be created (for simplicity, we assume* D* = 0 or the first input = 0):Henceis the observability matrix. The system is observable if is of . To test for controllability, one can create a state vector:Hereis a square matrix. The discrete time system is controllable if and only if . Suppose the system of (4) and (5) is controllable and observable with rank* n*; then One notices that the Hankel matrix is composed of the system Markov parameters. Using the ARX model parameters , can be given asDefining the system Markov parameters asthen can be written asTo realize a state-space system of the form a minimum realization is utilized (see [7, 10]). The minimum realization allows for investigating the primary or dominant modes of the system to be identified.

Perform a SVD on , where and the singular values are ordered asTruncate and to have columns, which results into and . Hence the following equality exists:Hence, the Hankel singular values are the square roots of the eigenvalues of :Now there are different ways to interpret the realization. They are organized by what combination and association is being taken:(i) input-normal form:(ii) balanced form:(iii) output normal form: Choosing the balanced form, we haveHaving realized and , we can utilize their original definitionHence, is first columns of or of ; is first rows of or of , . To compute , one forms Expressing in terms of :Using SVD on the last equation yieldsFrom this, one can compute :

#### 3. Estimation of ARMAX Model

Equation (12) is an ARX model. To incorporate a moving average term into (12), a third convolution term is added,where , , and are the model parameter matrices. The idea of using an ARMAX model is to accommodate the innovation and noise influence. The parameters of the ARX model can be estimated using standard least-squares techniques. For the ARMAX model, the following approach is adapted.

Definefor , where the subscript is used for denoting the time index. The innovation sequence can be estimated by using whereUsing and with , and definingwhere the subscript is used for the time index, the estimated ARMAX parameters are given byTo obtain a realization of the form of (39), the estimated model parameters of the ARX or ARMAX model are used to compute the system Markov parameters (see [11]) and construct a block Hankel matrix of the form given by (20) or (30). The realization then can be established by any of the three options given by (25), (26), or (27).

#### 4. Tabu Search Optimization

The Hankel matrix in (20) is composed of the sampled impulse response of the system. In addition to the system’s response, there is noise, unmodeled system dynamics, and responses from auxiliary systems that are coupled with the system to be identified. As there is no method to separate the given impulse response from coupled noise, unmodeled dynamics, and auxiliary system responses, we propose to imbed an optimization into the eigensystem realization (ERA) algorithm given by [8]. The optimization is done using an Enhanced Tabu Search (ETS) algorithm. The ETS employs a two-stage process, where during the first stage a list of promising areas in the search field is established. During the second stage, the promising areas are further searched by a regular TS algorithm. Details on the employed ETS and TS are found in [9]. The ERA is updated such that after selecting the number of states of the system—based on the magnitude of the singular values in —the entries of the selected singular values are searched within a close vicinity of their nominal value established by the SVD. The search areas are defined as a percentage of the largest singular value. The cost function for the ETS and TS is given bywhere is the estimated output from the resulting realized system. An easy extension to the given cost function is to incorporate a priori information of the system. For example, if the fundamental natural frequency of the system to be identified is known, (40) can be extended to where are weighting factors and is the resulting estimated natural frequency.

#### 5. Experimental Setup

In this work, we utilize a one-stage compressor system with a blade geometry that allows for spike inception. The proposed identification scheme is used to identify the dynamics at different flow conditions, including near stall conditions. The dynamics of the compressor is given by the changes within the blade passage area as well as the changes in the pressure rise coefficient computed by the input and output pressures. Hence, we extract three separate models, one for the dynamics within the blade passage due to air injection, one for the input/output relationship of the compressor due to throttle movements, and one that characterizes the dynamics due to air injection pressure as the input of the system and the pressure rise coefficient as the output of the system. The compressor employed for this research is a low-speed rotor with characteristic parameters as given in Table 1. The compressor exhibits stall at 40–50% of its rotating frequency. The operating range of this compressor is given by a flow coefficient between 0.58 and 0.49, where the latter bound represents the vicinity of the stall point. The experiments were carried out using a smooth casing; that is, there are no grooves within the tip clearance casings. The average tip clearance for the setup is 1 mm. A number of pressure sensors are used, as shown in Figure 1.