Identification Approach for Nonlinear MIMO Dynamics of Closed-Loop Active Magnetic Bearing System

Abstract

A systematic identification approach for the rotor/radial active magnetic bearing (rotor/RAMB) system is presented in this study. First, the system identification of the controller of commercial TMP is undertaken, and the corresponding linear dynamic models are constructed. To perfectly excite the nonlinearities of the rotor/RAMB system, a parallel amplitude-modulated pseudo-random binary sequence (PAPRBS) generator, which possesses the merits of no correlation among the perturbation signals, is employed. The dynamics of the rotor/RAMB system is identified with a Hammerstein–Wiener model. To reduce the difficulty of the identified two nonlinear blocks, the output nonlinear characteristics are estimated prior to the recursive process. Two conventional nonlinear model structures, i.e., NARX and NARMAX, are employed for comparison to verify the effectiveness of the identified Hammerstein–Wiener model. The averaged fit values of the Hammerstein–Wiener model, NARX model, and NARMAX model are 93.25%, 88.36%, and 76.91%, respectively.

1. Introduction

Compared to ball bearings, active magnetic bearings (AMBs) possess the merits of zero ware, energy conservation, and being contactless and lubrication free. Most of all, the stiffness and damping induced by the active magnetic bearings (AMBs) to support the rotor can be adjusted with the aid of controllers for AMBs. The adjustment of the stiffness and damping of the rotor/AMB system enables the rotor/AMB system to pass the critical/resonant speeds of system [1,2,3,4,5] and suppression of vibrations induced by the unbalance mass of rotor [6,7,8,9,10,11,12,13,14,15]. The performance of the rotor/AMB is dominated by the controller of AMB. One of the reliable approaches to developing controllers is to design the controller based on the identified mathematical model of a plant. Hence, the accuracy of dynamic models for the rotor/AMB system is important.

At present, many scholars have proposed numerous methods to identify the parameters or dynamic models of magnetic-bearing systems. Regarding the estimation of the damping and stiffness of the rotor/AMB system, Molina et al. proposed a gray-box modeling procedure to estimate the displacement stiffness and current stiffness of single-axis magnetic bearings [16]. Hua and Ming proposed an online identification method based on the unbalance vibration suppression algorithm to estimate the parameters of the active magnetic bearing [17]. Lauridsen et al. identified the parameters, i.e., current/force factors, displacement/force factors, and time constant of the first order approximation model through the closed-loop identification method and static-loading method [18]. Zhou et al. presented an unbalance response-based method for identifying closed-loop AMB stiffness and damping coefficients [19]. Xu et al. proposed an identification procedure for rotor/active magnetic-bearing dynamics. The model is established with the finite-element method (FEM) and is verified by various experimental rotor unbalance responses [20]. Prasad and Tiwari proposed a gyroscopic dynamic reduction method to identify the speed-dependent parameters of AMBs equipped with a flexible rotor [21].

On the other hand, scholars have proposed various identification procedures to identify the overall rotor/AMB dynamics. Martynenko investigated the nonlinear dynamics of rotor levitated by the passive–active magnetic bearing [22]. Khader et al. established a 2 × 2 MIMO system identification procedure. The resulted model is the combination of both rigid-mode and flexible-mode models [23]. M.H.R.A. Aziz and R. Mohd-Mokhtar proposed a two-stage identification approach with the continuous subspace method to model the MIMO magnetic-bearing system [24]. R. Mohd-Mokhtar and L. Wang proposed an identification approach to model the linear MIMO dynamics of the radial-magnetic-bearing system [25]. To identify more appropriate behaviors of a real rotor/AMB system, the nonlinear auto-regressive methods are applied. Liu et al. identified the model of zero-bias axial magnetic bearing with the NARX (nonlinear autoregressive with external input) neural network [26]. Miranda and Manzano estimated the coefficients of SISO (single-input, single-output) NARMAX (nonlinear auto-regressive moving average with exogenous variables) polynomial model of radial magnetic bearing [27]. Noshadi et al. proposed a genetic algorithms (GA) identification for AMB which can get the global optimum estimation [28].

Since the dynamics of an open-loop active magnetic bearing (AMB) system is inherently unstable, a controller is essential for an AMB system. However, the performance and the dynamics of the AMB system are dominated by the controller. As a result, the accuracy of the identified model employed to design the controller is important. That is, the more informative the system model is, the more suitable the model-based controller can be designed. However, few papers presented identification procedures for nonlinear MIMO dynamics of the closed-loop rotor/AMB system. Moreover, many nonlinear models of AMBs are identified with black-box identification approaches. It is hard to interpret the information from the black-box models. As most of the controllers for AMBs are model-based controllers, this study aims to propose an identification approach to estimate a high-accuracy dynamic model from a MIMO closed-loop system for controller design.

2. System Identification on Controllers for Radial Magnetic Bearings

The system to be identified in this work is the rotor/radial active magnetic bearing (rotor/RAMB) module of the commercial turbo-molecular pump (TMP), OSAKA K701, produced by Osaka Vacuum, Ltd. (Chuo-ku, Osaka). This section is aimed at identifying the controller embedded in the closed-loop TMP system. The control structure of the TMP is depicted in Figure 1. Since the control loop of TMP has integrated: (i) controller, (ii) power amplifier, (iii) shaft/rotor, and (iv) gap sensors, the following steps are to be undertaken so that the transfer function of the controller, $G_{\mathrm{C}}(s)$, can be obtained. The identification procedure of the controller is depicted in Figure 2.

Step #1.

Measure the upper and lower bounds of the control circuit’s input and output signals.

Step #2.

Disconnect the control circuit from the control loop of the TMP, but the power is still supplied to the control circuit.

Step #3.

Excite the control circuit with the ‘chirp’ perturbation signals and record the corresponding responses.

Step #4.

Construct the dynamic model of the controller via the aid of the commercial software, MATLAB.

The internal view of the electronic control unit (ECU) for TMP, Osaka TD700/1100, is shown in Figure 3. The detailed structures of the control circuit of AMBs are depicted in Figure 4. The notations, e.g., S1, S3, …, S9, CP2, CP4, …, CP8 and CP10 are referred to as test points. Embedded in the control circuits, we employ the notch filters to lower down the control force under critical speeds and make the shaft of TMP to spin around its inertia axis instead of its centroidal axis. By applying the notch filters, vibrations induced by the unbalance mass of the rotor can be, to some extent, suppressed.

The individual control circuit for AMBs, composed of resistors, capacitors, and operational amplifiers (Op-Amps), is of the analog type. As a result, the controller can be identified by a linear system approach. In addition, because these five sets of control circuits are not connected to each other, the scheme by linear single-input-single-output (SISO) system identification is employed.

The block diagram of system identification on the controller is depicted in Figure 5. As to the rotor/AMB unit, the load to the power amplifiers is nothing but the AMB coils. Moreover, the AMB coil is of inductive load whose dynamic behavior can be expressed by a first-order equation. To identify the dynamics of the controller, the AMB coils are replaced by resistors whose resistance values are the same as those of the AMB coils. Furthermore, an anti-aliasing filter, i.e., a second-order Butterworth filter, is cascaded ahead of the ADC (analog-to-digital converter) to restrict the bandwidth of the signal to satisfy the Nyquist–Shannon sampling theorem over the band of interest [29]. The test rig for system identification on the controller is shown in Figure 6. To identify the controller, the input signal is the type of chirp in a sinusoidal waveform whose frequency is increased exponentially. The initial frequency of the chirp signal is 0 Hz, and the final frequency is 5000 Hz. The time duration of the chirp signal is set as 10 s, and the sampling period of DS1104, $T_{\mathrm{S}}$, is set as 0.5 ms.

2.1. System Identification on Controllers for URAMB

There are two controllers applied to regulate the position deviations of a rotor by upper radial AMB (URAMB) along the X-axis and Y-axis, respectively. The controller for URAMB along the X-axis, $G_{\mathrm{UX}\_\mathrm{C}}$, is composed of a notch filter and a regulator, i.e., $G_{\mathrm{UX}\_\mathrm{Notch}}$ and $G_{\mathrm{UX}\_\mathrm{Reg}}$, respectively, shown in Figure 4. That is, $G_{\mathrm{UX}\_\mathrm{C}}=G_{\mathrm{UX}\_\mathrm{Reg}}{G}_{\mathrm{UX}\_\mathrm{Notch}}$. By direction of the signal flow, the notch filter is cascaded ahead of the regulator. To identify the controller for URAMB along the X-axis, the chirp signal is exerted to Point S3, shown in Figure 4, i.e., the output of the upper gap sensor along the X-axis and the corresponding response is recorded at Point CP4, i.e., the output of the regulator for URAMB along the X-axis. The resulting Bode plot of the controller for URAMB along the X-axis is depicted in Figure 7. It is observed that the frequency of the notch filter is around 550 Hz.

The order of the dynamic model of controller for URAMB along the X-axis is 4. The corresponding state space representation and the discrete-time transfer function are summarized as:

$${\mathrm{A}}_{\mathrm{UX}\_\mathrm{C}}=\left[\begin{array}{cccc}-3072.915& -1593.095& 0& 0\\ 1593.095& -3072.915& 0& 0\\ 0& 0& -635.947& -3357.384\\ 0& 0& 3357.384& -635.947\end{array}\right]$$

(1)

$${\mathrm{B}}_{\mathrm{UX}\_\mathrm{C}}=\left[\begin{array}{c}-147.249\\ 182.816\\ 13.242\\ -19.935\end{array}\right]$$

(2)

$${\mathrm{C}}_{\mathrm{UX}\_\mathrm{C}}=\left[\begin{array}{cccc}-14.930& -43.428& 24.018& -37.681\end{array}\right]$$

(3)

$${\mathrm{D}}_{\mathrm{UX}\_\mathrm{C}}=0.06795$$

(4)

$$\begin{array}{l}{G}_{\mathrm{UX}\_\mathrm{C}}=(0{.06795 \mathrm{z}}^4-0{.4431 \mathrm{z}}^3+0{.9146 \mathrm{z}}^2-0.7769 \mathrm{z}+0.2373)/\\ { (\mathrm{z}}^4-3{.62 \mathrm{z}}^3+4{.94 \mathrm{z}}^2-3.009 \mathrm{z}+0.6901), \\ \mathrm{where} T_{\mathrm{S}}=0.05 \mathrm{ms}\end{array}$$

(5)

Since the purpose of the controller for URAMB along the X-axis is the same as the controller for URAMB along the Y-axis, i.e., to levitate the rotor stably and to suppress the position deviation of the rotor at the position of URAMB at a certain frequency, the Bode plot of the controller for URAMB along the Y-axis is almost the same as Figure 7. On the other hand, the state space representation and the discrete-time transfer function of the controller for URAMB along the Y-axis are listed as:

$${\mathrm{A}}_{\mathrm{UY}\_\mathrm{C}}=\left[\begin{array}{cccc}-3117.165& -1570.445& 0& 0\\ 1570.445& -3117.165& 0& 0\\ 0& 0& -615.177& -3332.377\\ 0& 0& 3332.377& -615.177\end{array}\right]$$

(6)

$${\mathrm{B}}_{\mathrm{UY}\_\mathrm{C}}=\left[\begin{array}{c}152.878\\ -187.223\\ 13.747\\ -19.596\end{array}\right]$$

(7)

$${\mathrm{C}}_{\mathrm{UY}\_\mathrm{C}}=\left[\begin{array}{cccc}15.508& 43.751& 24.061& -37.467\end{array}\right]$$

(8)

$${\mathrm{D}}_{\mathrm{UY}\_\mathrm{C}}=0.06957$$

(9)

$$\begin{array}{l}{G}_{\mathrm{UY}\_\mathrm{C}}=(0{.06957 \mathrm{z}}^4-0{.452 \mathrm{z}}^3+0{.9319 \mathrm{z}}^2-0.7913 \mathrm{z}+0.2417)/\\ { (\mathrm{z}}^4-3{.619 \mathrm{z}}^3+4{.936 \mathrm{z}}^2-3.005 \mathrm{z}+0.6885), \\ \mathrm{where} T_{\mathrm{S}}=0.05 \mathrm{ms}\end{array}$$

(10)

2.2. System Identification on Controllers for LRAMB

Similarly, the controller for lower radial AMB (LRAMB) along the X-axis, $G_{\mathrm{LX}\_\mathrm{C}}$, comprises a notch filter and a regulator, i.e., $G_{\mathrm{LX}\_\mathrm{Notch}}$ and $G_{\mathrm{LX}\_\mathrm{Reg}}$, respectively, shown in Figure 4. That is, $G_{\mathrm{LX}\_\mathrm{C}}=G_{\mathrm{LX}\_\mathrm{Reg}}{G}_{\mathrm{LX}\_\mathrm{Notch}}$. The chirp signal is externally imported to Point S7 shown in Figure 4, i.e., the output of lower gap sensor along the X-axis, and the corresponding response is recorded at Point CP8, i.e., the output of the regulator for LRAMB along the X-axis. The Bode plot of the controller for LRAMB along the X-axis is depicted in Figure 8. Unlike the notch frequency of the controllers at URAMB, the notch frequency is designed as 446 Hz for the controllers at LRAMB.

Since the order of the dynamic model of controller for LRAMB along the X-axis is 4, the state space representation and its discrete-time transfer function are as:

$${\mathrm{A}}_{\mathrm{LX}\_\mathrm{C}}=\left[\begin{array}{cccc}-4559.940& 0& 0& 0\\ 0& -701.529& -2635.884& 0\\ 0& 2635.884& -701.529& 0\\ 0& 0& 0& -\mathrm{25,679}.137\end{array}\right]$$

(11)

$${\mathrm{B}}_{\mathrm{LX}\_\mathrm{C}}=\left[\begin{array}{c}-67.522\\ -12.302\\ 2.900\\ -96.335\end{array}\right]$$

(12)

$${\mathrm{C}}_{\mathrm{LX}\_\mathrm{CC}}=\left[\begin{array}{cccc}-39.352& -14.682& 29.212& 221.189\end{array}\right]$$

(13)

$${\mathrm{D}}_{\mathrm{LX}\_\mathrm{CC}}=0.1937$$

(14)

$$\begin{array}{l}{G}_{\mathrm{LX}\_\mathrm{CC}}=(0{.1937 \mathrm{z}}^4-1{.048 \mathrm{z}}^3+1{.961 \mathrm{z}}^2-1.558 \mathrm{z}+0.4515)/\\ { (\mathrm{z}}^4-2{.987 \mathrm{z}}^3+3{.207 \mathrm{z}}^2-1.422 \mathrm{z}+0.2055),\\ \mathrm{where} T_{\mathrm{S}}=0.05 \mathrm{ms}\end{array}$$

(15)

On the other hand, the controllers for LRAMB along the X-axis and along the Y-axis are to levitate the lower part of the rotor together; the characteristics of the controller for LRAMB along the Y-axis are similar to those of the controller for LRAMB along the X-axis. The Bode plot of the controller for LRAMB along the Y-axis is close to that of the controller for LRAMB along the X-axis. The identified models of the controller for LRAMB along the Y-axis are:

$${\mathrm{A}}_{\mathrm{LY}\_\mathrm{C}}=\left[\begin{array}{cccc}-4600.290& 0& 0& 0\\ 0& -674.501& -2639.121& 0\\ 0& 2639.121& -674.501& 0\\ 0& 0& 0& -\mathrm{24,903}.283\end{array}\right]$$

(16)

$${\mathrm{B}}_{\mathrm{LY}\_\mathrm{C}}=\left[\begin{array}{c}-67.886\\ -11.936\\ 2.836\\ -95.533\end{array}\right]$$

(17)

$${\mathrm{C}}_{\mathrm{LY}\_\mathrm{C}}=\left[\begin{array}{cccc}-39.548& -14.883& 28.915& 213.690\end{array}\right]$$

(18)

$${\mathrm{D}}_{\mathrm{LY}\_\mathrm{C}}=0.1811$$

(19)

$$\begin{array}{l}{G}_{\mathrm{LY}\_\mathrm{C}}=(0{.1811 \mathrm{z}}^4-0{.9952 \mathrm{z}}^3+1{.879 \mathrm{z}}^2-1.502 \mathrm{z}+0.437)/\\ { (\mathrm{z}}^4-2{.999 \mathrm{z}}^3+3{.238 \mathrm{z}}^2-1.45 \mathrm{z}+0.2138), \\ \mathrm{where} T_{\mathrm{S}}=0.05 \mathrm{ms}\end{array}$$

(20)

3. Experimental Setup for Identification of Rotor/RAMB Dynamic Model

Because the open-loop dynamics of the rotor/RAMB unit is inherently unstable, the dynamic model of rotor/RAMB unit can be identified accurately as the rotor is fully levitated by the AMBs under the regulation of the controller. To identify the rotor/RAMB dynamic model, the closed-loop system must be excited with the perturbation signals and to record the corresponding responses. However, the circuits of the controller and the power amplifier are firmly connected to each other. As the perturbation signals are directly injected into the closed-loop system without isolation of the impedances at the input points of controller and at the DAC (digital-to-analog converter) of DS1104, the perturbation signals would be distorted. The corresponding terminology of the aforesaid phenomenon is called “loading effects”. To minimize the loading effects, a summer module, shown in Figure 9, is additionally designed by the author and inserted between the gap sensors and the controller in series.

3.1. Design of Summer Module

The layout of the summer module is depicted in Figure 10. The summer module can be divided into three sections: namely buffer section, summer section, and inverting section, respectively. In the buffer section, two voltage followers, namely voltage buffers, are applied to minimize the load effect. Afterwards, one summer is employed to import the perturbation signal into the closed-loop rotor/AMB unit for system identification of the rotor/AMB dynamic system. Finally, a unity-gain inverting amplifier is applied to restore the polarity of the signal processed by the summer module. The photograph of the summer module is shown in Figure 11. The operational amplifier IC, namely, TA74075P, possesses the merits of (i) low-input bias current, (ii) Low-input offset current, (iii) low noise, and (iv) wide bandwidth. Therefore, TA74075P is selected for the summer module.

3.2. Parallel Amplitude-Modulated Pseudo-Random Binary Sequence (PAPRBS) Generator

Numerous types of excitation signals can be applied for system identification, e.g., pseudo-random binary sequence (PRBS), chirp signal, multi-sine signal, and stepped sinusoidal wave [30]. It is noted that the rated rotational speed of the TMP rotor is pretty high, e.g., 30,000 RPM for the TMP studied by this work, and the rated speed is higher than the first bending frequency of the TMP rotor. Consequently, the chirp signal is not suitable for system identification upon the rotor/AMB dynamics since the levitated rotor could collide with the auxiliary bearings as the frequency of the perturbation signal is close to/or across the resonance frequencies of the rotor/AMB unit. Among these commonly employed perturbation signals, PRBS, namely maximum length sequence or M-sequence, possesses the merits of deterministic, lowest crest factor and white-noise-alike properties. Therefore, the accuracy of system identification by PRBS can be guaranteed to some extent. In addition, since the system to be identified, i.e., the rotor/RAMB (rotor/radial active magnetic bearing) dynamics, is a MIMO system, the correlation degree between any two perturbation signals has to be reduced as much as possible. Fortunately, once the parallel PRBS generator is well-designed, the aforesaid correlations can be almost eliminated. In this study, hence, the parallel PRBS is adopted as the perturbation signal to identify the MIMO rotor/RAMB dynamic system. Generally, the number of amplitude types of a PRBS is 2, i.e., +V and −V. However, to explore the nonlinear dynamic properties of the rotor/RAMB unit, the amplitude of the perturbation signal has to vary so that the nonlinearities of the studied system, e.g., saturation, hysteresis, or high-order modes, can be appropriately excited. Therefore, in this study, the conventional parallel PRBS generator is extended to a parallel “amplitude-modulated” PRBS (PAPRBS) generator by multiplying a random number sequence. The proposed PAPRBS is designed and realized via the commercial software, Simulink. The layout of the Simulink code for the PAPRBS generator is shown in Figure 12. The order of the PAPRBS is set as 15. As a result, the maximum length of the 15th-order PAPRBS is $2^{15}-1=32767$. Four separated sets of perturbation signals can be generated by the PAPRBS generator simultaneously to identify the rotor/RAMB dynamics. To make the four sets of the perturbation signals uncorrelated, the start points of the four sets of perturbation signals are set at: (i) the beginning, (ii) the one-fourth, (iii) the one-half, and (iv) the three-fourth of the 15th-order PRBS, respectively. The length of each perturbation signal sequence is made up by 8000 points such that the correlation between any two perturbation signals can be almost eliminated. The first 200 points of these four perturbation signal sequences are shown in Figure 13. It can be observed that the patterns of these four signals are quite different. To be more specific, the correlations of these four perturbation signals are evaluated and shown in Figure 14. From Figure 14, no cross-correlation between any two of these four perturbation signals is found. It is noted that the bit duration of PAPRBS has to be properly chosen so that the bandwidth of the perturbation signal can cover the operational range of the object/system to be identified. The relation between the maximum frequency of the identified model, i.e., $f_{\max }$, and the bit duration of PAPRBS, i.e., $\Delta$, has to satisfy the following condition [31]:

$$2\pi /3\Delta >f_{\max }$$

(21)

Since the maximum frequency of the identified model of rotor/AMB unit is set as 2000 Hz, the bit duration of PAPRBS signal has to meet the following inequality:

$$\Delta <2\pi /(3*2000)=1.047 \mathrm{ms}$$

(22)

Finally, the bit duration is selected as 0.5 ms, and the corresponding period of each PAPRBS as four seconds (i.e., 8000 × 0.5 ms).

3.3. Parallel Amplitude-Modulated Pseudo-Random Binary Sequence (PAPRBS) Generator

The test rig for the system identification on the multi-input-multi-output (MIMO) rotor/RAMB dynamic system is shown in Figure 15. The corresponding schematic diagram is depicted in Figure 9. In the experiments, the perturbation signals, i.e., PAPRBS signals, are imported to the input points of the controller, i.e., S3, S5, S7 and S9 shown in Figure 4, as the rotor is fully levitated by the commercial controller. Afterwards, the perturbation signals and their corresponding responses, i.e., position deviations of the rotor, are recorded and will be used later.

To construct the rotor/RAMB dynamic model by experiments, the input signals to the plant and the corresponding responses of the rotor/RAMB unit, i.e., Point B and Point C shown in Figure 9, respectively, are necessary to be recorded initially. The plant is referred to as the rotor/RAMB unit in this study. In order not to introduce the noises to the power amplifier directly, the input signals to the plant are not recorded at Point B but at Point A. The detailed steps to generate the corresponding responses at Point B, shown in Figure 9, by applying the perturbation signals recorded at Point A, as shown in Figure 9, to the dynamic model of controllers of AMBs are summarized as:

(i).

Four perturbation signals are imported at the input points of the controllers, as the TMP rotor is fully levitated.

(ii).

Meanwhile, record the four perturbation signals imported at the input points of the controller, i.e., at Point A, shown in Figure 9.

(iii).

Import the recorded four perturbation signals to the dynamic model of the controllers of AMBs, i.e., Equation (5), Equation (10), Equation (15) and Equation (20), respectively.

(iv).

Record the corresponding responses of the plant model with the controllers.

The flowchart of responses generation with perturbation signals is depicted in Figure 16. The sampling period, $T_{\mathrm{S}}$, is set as $0.05 \mathrm{ms}$. The input signals imported to the dynamic model of the controllers and the corresponding responses are shown in Figure 17. As the input signals to the rotor/RAMB dynamics and the corresponding responses, i.e., Point B and Point C, shown in Figure 9, respectively, are obtained, the rotor/RAMB dynamic model can therefore be identified.

4. System Identification of Rotor/RAMB Dynamics

The 5-DOF (degree of freedom) rotor/active magnetic bearing (rotor/AMB) unit at TMP can be divided into a 4-DOF rotor/radial active magnetic bearing (rotor/RAMB) module and a 1-DOF rotor/axial active magnetic bearing (rotor/AAMB) module. Except the controller for active magnetic bearings (AMBs), the hardware within a closed-loop AMB system comprises power amplifiers, rotor/AMB module, and gap sensors. Compared with the dynamic responses of the mechanical part, i.e., rotor/RAMB module, the dynamic responses of the power amplifier and the gap sensors are much fast. That is, it is reasonable to describe the input/output characteristics of the power amplifier and the gap sensors with static models in the rotor/AMB system. It is noted that the scheme of the Hammerstein–Wiener model [31] is close to the composition of the hardware within a closed-loop AMB system. Hence, the model of the rotor/RAMB unit is described by a MIMO Hammerstein–Wiener model representation in this study. The schematic diagram of the MIMO Hammerstein–Wiener model is depicted in Figure 18. The MIMO Hammerstein–Wiener model comprises a MIMO nonlinear input element, a MIMO linear element, and several SISO nonlinear output elements. The nonlinearities and the dynamic properties of the rotor/RAMB module can be well described by the combination of the Hammerstein model and the Wiener model.

As to the Hammerstein model, the nonlinear static element, $\mathrm{H}$, is cascaded ahead of the linear dynamic element, $\mathrm{L}$. In contrast, the linear dynamic element, $\mathrm{L}$, is cascaded ahead of the nonlinear static elements, $w_k, k=1,2,\mathrm{..},j$, in the Wiener model. The nonlinear element of the Hammerstein model, $\mathrm{H}$, is to represent the actuator nonlinearities or any other potential nonlinear effects embedded on the system input side. In comparison, the nonlinear elements of the Wiener model, $w_k, k=1,2,\mathrm{..},j$, are referred to as the output nonlinearities, which are the nonlinear characteristics of gap sensors in this study. The symbols, $\mathrm{U}$ and $y_k, k=1,2,\mathrm{..},j$, denote the inputs and the outputs of the Hammerstein–Wiener model. The symbols, ${\omega }_k, k=1,2,\mathrm{..},j$ and ${\upsilon }_k, k=1,2,\mathrm{..},j$, are the white noises and colored noises, respectively. The relations between the elements of the Hammerstein–Wiener model are summarized as:

$$\mathrm{X}=\mathrm{H}(U)=U+{\displaystyle \sum _{l=2}^{n_h}{a}_l}{h}_l(U), where U={[u_1,u_2,\cdots ,u_i]}^T$$

(23)

$$\begin{array}{l}\mathbf{V}=\mathbf{L}·\mathbf{X}+\mathbf{N}·\varpi \\ =\left({\displaystyle \sum _{k=1}^n{\mathbf{C}}_k}{z}^{-n}\right)\mathbf{X}/\left(1+{\displaystyle \sum _{l=1}^n{\mathbf{D}}_l}{z}^{-n}\right)+\left(1+{\displaystyle \sum _{k=1}^m{\mathbf{E}}_k}{z}^{-m}\right)\varpi /\left(1+{\displaystyle \sum _{l=1}^m{\mathbf{F}}_l}{z}^{-m}\right)\end{array}$$

(24)

$${\mathrm{v}}_k(t)={\mathrm{w}}_k^{-1}(y_k)={\displaystyle \sum _{l=1}^{n_w}{b}_{kl}}{y}_k^l(t),k=1,2,\cdots ,j$$

(25)

where $i$ is the number of inputs, $j$ is the number of outputs, $n_h$ is the degree of the nonlinear input functions, and $n_w$ is the degree of inverse nonlinear output functions. $n$ is the order of the linear dynamic model, and $m$ is the order of the colored-noise model. $a_l$ and $b_{kl}$ are the scale parameters, while ${\mathrm{C}}_k$ denotes the parameter matrices. The diagonal parameter matrices, ${\mathrm{D}}_l$, ${\mathrm{E}}_k$, and ${\mathrm{F}}_l$, are represented as:

$${\mathrm{D}}_l=\left[\begin{array}{cccc}{d}_{1l}& & & 0\\ & d_{2l}& & \\ & & \ddots & \\ 0& & & d_{jl}\end{array}\right], where l=1,2,\mathrm{..},n$$

(26)

$${\mathrm{E}}_k=\left[\begin{array}{cccc}{e}_{1k}& & & 0\\ & e_{2k}& & \\ & & \ddots & \\ 0& & & e_{jk}\end{array}\right], where k=1,2,\mathrm{..},m$$

(27)

$${\mathrm{F}}_l=\left[\begin{array}{cccc}{f}_{1l}& & & 0\\ & f_{2l}& & \\ & & \ddots & \\ 0& & & f_{jl}\end{array}\right], where l=1,2,\mathrm{..},m$$

(28)

The Hammerstein–Wiener model is a gray-box model, which processes the merits of the white-box model and the black-box model. To reduce the uncertainties of the Hammerstein–Wiener model but keep the nonlinearities of the dynamic model, we can measure the characteristics of the gap sensors prior to the process of parameter estimation. The parameters of the inverse equations of static nonlinear elements shown in Equation (25) are known. The remaining parameters to be estimated are the part of the Hammerstein model, i.e., Equations (23) and (24). Equations (23) and (24) are organized into one multivariable regression equation system:

$$\mathbf{V}\left(\mathbf{t}\right)={\displaystyle \sum _{k=1}^n{\mathsf{\phi }}_{d_k}(t){\mathsf{\theta }}_{d_k}}+{\displaystyle \sum _{l=1}^{n_h}{\mathsf{\phi }}_{a{c}_l}(t)}{\mathsf{\theta }}_{a{c}_l}+\upsilon (t)$$

(29)

where $\upsilon (t)$ are the colored noises. Normally, colored-noise terms can be regarded as the white noises filtered with the transfer functions. In this study, the colored-noise is assumed to be the ARMA-type. Therefore, Equation (29) can be represented as the auxiliary model:

$$\mathbf{V}\left(\mathbf{t}\right) ={\displaystyle \sum _{k=1}^n{\mathsf{\phi }}_{d_k}(t){\mathsf{\theta }}_{d_k}}+{\displaystyle \sum _{l=1}^{n_h}{\mathsf{\phi }}_{a{c}_l}(t)}{\mathsf{\theta }}_{a{c}_l}+{\displaystyle \sum _{l=1}^m{\mathsf{\phi }}_{e_l}(t)}{\mathsf{\theta }}_{e_l}+{\displaystyle \sum _{l=1}^m{\mathsf{\phi }}_{f_l}(t)}{\mathsf{\theta }}_{f_l}+\varpi (t)$$

(30)

where

$$\mathbf{V}(t)=\left[\begin{array}{c}{v}_1(t)\\ v_2(t)\\ ⋮\\ v_j(t)\end{array}\right],$$

(31)

$${\mathsf{\phi }}_{d_k}(t)=-\left[\begin{array}{cccc}{v}_1(t-k)& & & 0\\ & v_2(t-k)& & \\ & & \ddots & \\ 0& & & v_j(t-k)\end{array}\right], k=1,2,\mathrm{..},n$$

(32)

$${\mathsf{\theta }}_{d_k}={\left[d_{1k}\hspace{0.33em}{d}_{2k}\hspace{0.33em}\cdots d_{jk}\right]}^T, k=1,2,\mathrm{..},n$$

(33)

$${\mathsf{\phi }}_{a{c}_l}(t)=\left[\begin{array}{cccc}{h}_1{(U(t-j))}^T& & & 0\\ & h_2{(U(t-j))}^T& & \\ & & \ddots & \\ 0& & & h_i{(U(t-j))}^T\end{array}\right], l=1,2,\mathrm{..},n_h$$

(34)

$${\mathsf{\theta }}_{ac}{}_{{}_l}={\left[{\mathrm{G}}_l^T{(:,1)}^T\hspace{0.33em}{\mathrm{G}}_l^T{(:,2)}^T\hspace{0.33em}\cdots \hspace{0.33em}{\mathrm{G}}_l^T{(:,n)}^T\hspace{0.33em}\right]}^T, l=1,2,\mathrm{..},n\times n_h$$

(35)

$${\mathsf{\phi }}_{e_l}(t)=\left[\begin{array}{cccc}{\omega }_1(t-l)& & & 0\\ & {\omega }_2(t-l)& & \\ & & \ddots & \\ 0& & & {\omega }_j(t-l)\end{array}\right], l=1,2,\mathrm{..},m$$

(36)

$${\mathsf{\theta }}_e{}_{{}_l}={\left[e_{1l} e_{2l}\hspace{0.33em}\cdots \hspace{0.33em}{e}_{jl}\right]}^T, l=1,2,\mathrm{..},m$$

(37)

$${\mathsf{\phi }}_{f_l}(t)=-\left[\begin{array}{cccc}{\upsilon }_1(t-l)& & & 0\\ & {\upsilon }_2(t-l)& & \\ & & \ddots & \\ 0& & & {\upsilon }_j(t-l)\end{array}\right], l=1,2,\mathrm{..},m$$

(38)

$${\mathsf{\theta }}_f{}_{{}_l}={\left[f_{1l}{f}_{2l}\hspace{0.33em}\cdots \hspace{0.33em}{f}_{jl}\right]}^T, l=1,2,\mathrm{..},m$$

(39)

$$\varpi (t)=\left[\begin{array}{c}{\omega }_1(t)\\ {\omega }_2(t)\\ ⋮\\ {\omega }_j(t)\end{array}\right],$$

(40)

where the vectors ${\mathrm{G}}_l^T{(:,n)}^T$ shown in Equation (35) are extracted from ${\mathrm{G}}_{n\times n_h}=C_n{a}_{{\mathrm{n}}_h}$. Equation (29) is rewritten in the multivariable linear regression form as:

$$\begin{array}{l}\mathbf{V}(t)=\mathsf{\phi }(t)\mathsf{\theta }+\varpi (t)\\ =\left[\begin{array}{cc}{\mathsf{\phi }}_L(t)& {\mathsf{\phi }}_N(t)\end{array}\right]\left[\begin{array}{l}{\mathsf{\theta }}_L\\ {\mathsf{\theta }}_N\end{array}\right]+\varpi (t)\end{array}$$

(41)

where ${\mathsf{\phi }}_L(t)=\left[{\mathsf{\phi }}_{d_1}(t)\cdots \hspace{0.33em}{\mathsf{\phi }}_{d_n}(t) {\mathsf{\phi }}_{a{c}_1}\cdots \hspace{0.33em}{\mathsf{\phi }}_{a{c}_{n_h}}(t)\right]$, ${\mathsf{\phi }}_N(t)=\left[{\mathsf{\phi }}_{e_1}(t)\hspace{0.33em}\cdots \hspace{0.33em}{\mathsf{\phi }}_{e_{\mathrm{m}}}(t)\hspace{0.33em}{\mathsf{\phi }}_{f_1}(t)\hspace{0.33em}\cdots \hspace{0.33em}{\mathsf{\phi }}_{f_{\mathrm{m}}}(t)\hspace{0.33em}\right]$, ${\mathsf{\theta }}_L={\left[{\mathsf{\theta }}_{d_1}^T\cdots \hspace{0.33em}{\mathsf{\theta }}_{d_n}^T {\mathsf{\theta }}_{a{c}_1}^T\hspace{0.33em}\cdots \hspace{0.33em}{\mathsf{\theta }}_{a{c}_{n_h}}^T\right]}^T$, and ${\mathsf{\theta }}_N={\left[{\mathsf{\theta }}_{e_1}^T\cdots \hspace{0.33em}{\mathsf{\theta }}_{e_m}^T {\mathsf{\theta }}_{f_1}^T\cdots \hspace{0.33em}{\mathsf{\theta }}_{f_m}^T\right]}^T$. It is noted that the terms shown in Equation (36) and Equation (38) cannot be measured from the dynamic system directly. By applying the following equations, the estimated values of these terms can be derived by [33]:

$$\begin{array}{l}\widehat{\varpi }(\mathit{t})=\mathbf{V}(\mathit{t}) - \mathsf{\phi }(\mathit{t})\mathsf{\theta }\\ =\mathbf{V}(\mathit{t}) - \left[\begin{array}{cc}{\mathsf{\phi }}_L(t)& {\widehat{\mathsf{\phi }}}_N(t)\end{array}\right]\left[\begin{array}{l}{\widehat{\mathsf{\theta }}}_L\\ {\widehat{\mathsf{\theta }}}_N\end{array}\right]\end{array}$$

(42)

$$\widehat{\upsilon }(t)=\mathbf{V}(t) - {\mathsf{\phi }}_L(t){\widehat{\mathsf{\theta }}}_L$$

(43)

Therefore, the parameters of the MIMO dynamic system can be estimated by the unbiased identification algorithm:

$$\widehat{\mathsf{\theta }}(t)=\widehat{\mathsf{\theta }}(t-1)+K(t)\left[\mathrm{V}(t)-\mathsf{\phi }(t)\widehat{\mathsf{\theta }}(t-1)\right],$$

(44)

$$K(t)=P(t-1){\mathsf{\phi }}^T(t){\left[\mathsf{\phi }(t)P(t-1){\mathsf{\phi }}^T(t)+I\right]}^{-1},$$

(45)

$$P(t)=\left[I-K(t)\mathsf{\phi }(t)\right]P(t-1),$$

(46)

The products of the parameter matrices, $C_i (i=1,\cdots ,n)$ and ${\widehat{a}}_l (l=2,\cdots ,n_h)$, shown in the matrix, ${\mathsf{\theta }}_P$, are defined as:

$${\mathsf{\theta }}_P=\left[\begin{array}{cccc}{C}_1^T& {\mathrm{C}}_2^T& \cdots & {\mathrm{C}}_n^T\\ a_2^T{C}_1^T& a_2^T{C}_2^T& \cdots & a_2^T{\mathrm{C}}_n^T\\ ⋮& ⋮& \ddots & ⋮\\ a_{n_h}^T{C}_1^T& a_{n_h}^T{C}_2^T& \cdots & a_{n_h}^T{\mathrm{C}}_n^T\end{array}\right],$$

(47)

$$a={\left[1,a_2,\cdots ,a_{n_h}\right]}^T,$$

(48)

$$C={\left[{\mathrm{C}}_1^T,{\mathrm{C}}_2^T,\cdots ,{\mathrm{C}}_n^T\right]}^T,$$

(49)

The estimation of $C_l (l=1,\cdots ,n)$, i.e., ${\widehat{C}}_l (l=1,\cdots ,n)$, can be directly obtained from:

$${\widehat{C}}_l={\widehat{\mathrm{G}}}_{(l-1)n_h+1}, l=1,\cdots ,n$$

(50)

The pseudo inverse of the non-square matrix $\mathrm{C}$ shown in Equation (49), i.e., ${\mathrm{C}}^{†}$ in Equation (47), can be derived by the Moore–Penrose inverse [34]. Hence, ${\widehat{a}}_l ( l=2,\dots ,n_h)$ can be derived by the following equation:

$${\widehat{a}}_l={\left({\mathsf{\theta }}_P{\mathrm{C}}^{†}\right)}^T,\hspace{0.33em}l=2,\dots ,n_h$$

(51)

The flowchart of the identification procedure of the Wiener–Hammerstein model is depicted in Figure 19. The identification procedure is accomplished as the following condition is met:

$$\underset{\forall i}{\max }\left|\left[{\widehat{\theta }}_i(t)-{\widehat{\theta }}_i(t-1)\right]/{\widehat{\theta }}_i(t-1)\right|<\Delta$$

(52)

The parameter $\Delta$ shown in Equation (52) is set as 0.01.

As the perturbations fed in the closed-loop system, the perturbation signals and the responses of the plant, i.e., radial position deviations of the rotor, are recorded. The data length for validation is 4000, which is equivalent to the time length of 0.2 s. First, the estimated controller outputs are generated by feeding the recorded perturbation signals generated by the proposed parallel amplitude-modulated PRBS (PAPRBS) to the controller models presented in Section 2. The corresponding controller outputs are shown in Figure 20. The UX, UY, LX, and LY denote the upper RAMB along the X-axis, the upper RAMB along the Y-axis, the lower RAMB along the X-axis, and the lower RAMB along the Y-axis, respectively.

To verify the effectiveness of the identified Hammerstein–Wiener model, two conventional nonlinear model structures, i.e., NARX and NARMAX, are employed for comparison purposes. The NARX model is constructed by the MATLAB System Identification Toolbox. Both numbers of the delayed outputs, delayed inputs are selected as 4. The nonlinear estimator of NARX is selected as the Sigmoid network. The other nonlinear model for comparison is the NARMAX model with the FROLS (forward regression orthogonal least squares) algorithm [34]. The maximum numbers of the delayed outputs, delayed inputs are both selected as 4, and the maximum delayed noise is selected as 3. After several trials, the minimum number of NARMAX terms is selected as 14. By feeding the inputs to the identified models, shown in Figure 20, the responses of the proposed Hammerstein–Wiener model, NARX model, and NARMAX model are depicted in Figure 21. It is observed that the predicted responses of the Hammerstein–Wiener model show the best fit to the validation data, while the worst case is the NARMAX model. Although the predicted responses of the NARX model shows a good approximation result, the responses by the NARX model sometimes diverge with other sets of validation data.

To quantify the accuracy of the identified models, the following equation applies to evaluate the validation fit:

$$Fit=100\%\times (1-\frac{\sqrt{{\displaystyle \sum _{k=1}^N{\left[(y_k^{Valid}-{\overline{y}}^{Valid})-(y_k^{Est}-{\overline{y}}^{Est})\right]}^2}}}{‖y^{Valid}-{\overline{y}}^{Valid}‖})$$

(53)

where $N$ denotes data length, i.e., 4000. $y_k^{Valid}$ and $y_k^{Est}$ are the kth validation data and kth identified output, respectively. The upper bar represents the mean values. The fit values are summarized in

Table 1. Fit values between validation data and responses of Hammerstein–Wiener model, NARX model, and NARMAX model.

Fit (%)	Hammerstein–Wiener Model	NARX Model	NARMAX Model
UX	92.11%	87.63%	75.91%
UY	93.38%	88.88%	77.88%
LX	91.51%	85.57%	78.32%
LY	95.99%	91.34%	75.54%
Average	93.25%	88.36%	76.91%

. The averaged fit values of the Hammerstein–Wiener model, NARX model, and NARMAX model are 93.25%, 88.36%, and 76.91%, respectively.

The significant limitations of this study are addressed as follows:

• By applying the proposed modeling approach presented in Section 4, the mathematical model, i.e., the nonlinear element of the Wiener model, should be estimated prior to the recursive identification algorithm being undertaken. For instance, the static relations between the inputs and outputs of the gap sensor are measured with a high-precision positioning platform.
• As the perturbation signals are being injected into the rotor/RAMB system, the responses of the rotor/RAMB system will be distorted by the vibrations of the surrounding objects. How the engineering of vibration isolation performs will affect the accuracy of the identified model.

5. Conclusions

In this study, a systematic procedure to identify the nonlinear multi-input multi-output (MIMO) dynamic model from the closed-loop rotor/radial active magnetic bearing (rotor/RAMB) system is proposed. To reduce the potential correlation problem caused by the perturbation signals, a generator of perturbation signals for the MIMO system, named parallel amplitude-modulated pseudo-random binary sequence (PAPRBS), is proposed. By applying the proposed parallel amplitude-modulated pseudo-random binary sequence (PAPRBS) generator, the accuracy of the identified MIMO dynamics can be guaranteed to some extent. Since the composition of the hardware within a closed-loop AMB system is close to the scheme of the Hammerstein–Wiener model, the dynamic behaviors are described as a MIMO Hammerstein–Wiener model by the recursive least squares algorithm. To reduce the difficulty of the identified two nonlinear blocks of the Hammerstein–Wiener model, the output nonlinear characteristics are estimated prior to the recursive process. Two conventional nonlinear model structures, i.e., NARX and NARMAX, are employed for comparison to verify the effectiveness of the identified Hammerstein–Wiener model. The averaged fit values of the Hammerstein–Wiener model, NARX model, and NARMAX model are 93.25%, 88.36%, and 76.91%, respectively. Future work includes development of the online recursive identification algorithm for estimating the nonlinear MIMO dynamics of closed-loop systems.