# Intelligent Approach to Robust Design Optimization of a Rotor System due to Its Support Stiffness Uncertainty

Authors: Bensheng Xu; Chaoping Zang; Genbei Zhang
Journal: Shock and Vibration
 (2020)

---

## Abstract

In this paper, an intelligent robust design approach combined with different techniques such as polynomial chaos expansion (PCE), radial basis function (RBF) neural network, and evolutionary algorithms is presented with a focus on the optimization of the dynamic response of a rotor system considering support stiffness uncertainty. In the proposed method, the PCE method instead of the traditional Monte Carlo uncertainty analysis is applied to analyze the uncertain propagation of system performance. The RBF network is introduced to establish the approximate models of the objective and constraint functions. Taking the low-pressure rotor of a gas turbine with support stiffness uncertainty as an example, the optimization model is established with the mean and variance of unbalanced response of the rotor system at different operating speeds as the objective function, and the maximum unbalance response is less than the upper limit as the constraint function. The polynomial chaos expansion is generated to facilitate a rapid analysis of robustness in the presence of support stiffness uncertainties that is defined in terms of tolerance with good accuracy. The optimal Hypercubus are used as experimental plans for building RBF approximation models of the objective and constraint functions. Finally, the robust solutions are obtained with the multiobject optimization algorithm NSGA-II. Monte Caro simulation analysis demonstrates that the qualified rate of maximum vibration responses of the low-pressure rotor system can be increased from 83.6% to over 99%. This approach to robust design optimization is shown to lead to designs that significantly decrease vibration responses of the rotor system and improved system performance with reduced sensitivity to support stiffness uncertainty.


---

## Body

## 1. Introduction

The presence of some degrees of uncertainty in characterizing any real engineering system is inevitable. The traditional design optimization to seek only optimality may not achieve its expected robust performance. The robust design optimization (RDO), which is inextricably linked with the name of Taguchi who initiated a highly influential design philosophy, has been developed to improve the quality of a product by minimizing the effects of variation without eliminating the causes in the past twenty years [1].The RDO methods have been applied to problems with static and dynamic systems. Dynamic systems are those whose target value depends on the input signal set by the system user or operator. The concept of dynamic robustness in structural design was proposed by Oka and Yamakawa [2]. Tui [3] investigated and compared different modelling approaches for dynamic robust design problems. Zang et al. [4] studied the robust design of a vibration absorber with mass and stiffness uncertainty based on optimization approach, which shows that robust design optimization has great potential for application in structural dynamics to deal with the uncertainty. At present, the ways of RDO of structural dynamics can be considered as three types: Taguchi method based [5, 6], response surface method (RSM) based [7], and multiobjective optimization based [8, 9]. Most of the Taguchi-based methods use direct experimentation, and the objective functions for the optimization are expressed as the signal-to-noise ratio (SNR). In addition, an ANOVA may be performed on the SNR to look for variables that affect both the product variance and the mean so that they may be used to minimize the variation. Taguchi’s techniques have been commonly used in industrial applications because of its simplicity and intuitiveness. However, there are certain limitations such that the SNR may lead to nonoptimal solutions, efficiency loss, and information loss [10]. RSM is a collection of statistical and mathematical techniques useful for developing, improving, and optimizing processes. The basic idea of RSM is to approximately replace the complex performance analysis model with a simple function (polynomial function) to improve the computational efficiency of the optimization process [11]. The multiobjective optimization-based approaches, which simultaneously optimize both the mean performance and the variance in performance for robust design, are based on nonlinear programming methods. A trade-off decision must be made, to choose the best design with the maximum robustness.Rotor system is one of the most critical subsystems of rotating machinery and influences significantly the performance, life, and reliability of the machine. In the process of machining and assembly, due to the changes in machining accuracy, assembly environment, and assembly process, there are random uncertainties, which are inevitable in mass production. To take into account such uncertainties, the RDO of the rotor system has been developed by some researchers and practitioners. Ritto et al. [12] proposed an RDO strategy to seek the optimal parameters for which the natural frequencies of the system are as far away as possible from the rotational speeds of the machine. Ghisu et al. [13] studied the robust design optimization of gas turbine compression systems. Stocki et al. [14] developed an RDO method that combines an acceptable structural weight with the robustness with respect to uncertainties of residual unbalances for the optimization of the typical single-span rotor shaft of the 8-stage centrifugal compressor. Lopez et al. [15] presented a new method that a penalization was modelled as a random variable. The robust optimization was performed by minimizing the expected value and variance of the penalization, resulting in a multiobjective optimization problem. However, a multiobjective optimization for robust design usually requires a significantly high number of simulations for evaluating the response statistics of the performance functions and a selected optimization method. Therefore, it is costly and time-consuming. In this paper, an intelligent robust design approach combined with different techniques such as polynomial chaos expansion (PCE), radial basis function (RBF) neural network, and evolutionary algorithms is presented with a focus on the optimization of the dynamic response of a rotor system considering support stiffness uncertainty. In the proposed method, the PCE method instead of traditional Monte Carlo uncertainty analysis is applied to analyze the uncertain propagation of system performance. The RBF network is introduced to establish the approximate models of the objective and constraint functions.The rest of this paper has been organized in the following sequence: Section2 presents a comprehensive illustration of the developed intelligent approach to robust design optimization. In Section 3, an example of the robust parameter design of the low-pressure rotor of a gas turbine is used to demonstrate the application of these techniques. The numerical results verify the feasibility of the proposed method. Section 4 provides a summary and conclusions as a whole.


## 2. Intelligent Approach to Robust Design Optimization

The framework of intelligent approach to robust design optimization, in which design of experiment (DOE), PCE, and RBF neural network are combined to obtain reliable approximate models (AMs) and reliable optimization results with a short number of simulations and a high reduction of computational cost, is depicted in Figure1.Figure 1
Framework of robust design optimization.In Figure1, the robust design optimization procedure from the framework is considered as four steps. The first step is the definition and execution of the DOE. It is well known that DOE is a systematic approach to get the maximum amount of information out of various types of experiments while minimizing the number of experiments. A set of samples containing the input variables and the associated outputs is obtained. In the second step, the PCE method is adopted to calculate the variance (the robustness measure of system performance) of the concerned response for each sample point of input variables performed during the DOE phase. The third step is to construct AMs of objective functions and objective robustness functions, constraint functions, and constraint robustness functions. Once the AMs have been assessed and found to be sufficient, the robust solutions can be quickly achieved on such AMs. The use of AM for optimization has to be done cautiously, ensuring sufficient quality of the AM. The quality indexes on the approximation capability of the AM have been determined by using multidetermination coefficient R2. The value of coefficient R2 close to 1 (0) indicate that the approximate analytical model will (not) perform well for points that have not been simulated. In the last step, the RDO model is established in the fourth step. RDO is inherently a multiobjective and nondeterministic optimization problem in which both the expected value and the standard deviation of the performance function are to be minimized [16]. In contrast to single-objective optimization, there is no single global solution for multiobjective optimization problem and it is often necessary to determine a set of points that all fit a predetermined definition. The set of points is called as the Pareto set. The proposed RDO method has been fully developed within the Isight Optimus software. To exploit the full potential of the software in automating the simulations and in the parallel computing, a simple workflow shown in Figure 2 has been created.Figure 2
Isight workflow.The Isight workflow consists of two application components (MATLAB and Approximation) and one process component (Optimizaiton 1). The application component, MATLAB, contains all the design input variables with their variability ranges and is used to perform dynamic analysis and uncertainty quantification. The “Approximation” component executes DOE and constructs approximate models with RBF neural network. In the Optimization 1 component, design variables, constraints, and object function are established. The Pareto set can be obtained by using the multiobject optimization algorithm integrated in the Isight Optimus software. At the end, the results are saved in an output file and are extracted in the Output Variables block. This workflow has allowed exploiting at the best the two processors’ workstation used for the calculations, carrying out up to 12 simulations in parallel with a large saving of computational time.From Figure1, it is clearly shown that the RDO is carried out with combined different techniques including DOE, PCE, and RBF neural network. In the proposed RDO framework, the PCE method is utilized to propagate uncertainty, RBF neural network is used to estimate both the mean and standard deviation functions, and the multiple-objective optimization is to ensure the integrity of the content through the mathematical form of robust. The more details can be illustrated in the following.


### 2.1. Polynomial Chaos Expansion

As a newly developed method of uncertainty analysis, PCE has attractive attributes for its strong mathematical rigor and ability to produce functional representations of any stochastic quantities. One important advantage of PCE is that it can conveniently generate the complete probabilistic distribution (e.g., probability density function) of the output response [17]. The roots of polynomial chaos (PC) are found in the homogeneous chaos expansion first proposed by Wiener [18], who employed Hermite polynomials in terms of Gaussian random variables to express stochastic processes with finite variance. An extension was proposed by Xiu and Karniadakis [19] to deal with more general random inputs more efficiently. Suppose that Y (τ) is a random process in the probability space (Ω, F, P), the generalized polynomial chaos expansion of Y (τ) can be expressed as(1)Yτ=α0I0+∑i1=1∞αi1I1ξi1τ+∑i1=1∞∑i2=1i1αi1i2I2ξi1τ,ξi2τ+∑i1=1∞∑i2=1i1∑i3=1i2αi1i2i3I3ξi1τ,ξi2τ,ξi3τ+⋯,where α0,αi1,αi1i2,αi1i2i3… are the undetermined coefficients. The symbol τ indicates that the quantity involved is a random variable (for simplicity, τ will be omitted later). Inξ (n ∈ N, N represents the collection of natural numbers) is an nth-order Wiener–Askey orthogonal polynomial. Inξ is a function of multi-dimensional standard random variables ξ=ξi1,ξi2,…,ξin. The dimension of the random variable ξ is usually consistent with the dimension of the original random input variable X=xi1,xi2,…,xin in the probability space, and there is a one-to-one correspondence between ξ and X. Considering that only limited items can be calculated in practical applications, equation (1) is rewritten as(2)Y=∑i=0Npca^iI^iξ,where α^i is the undetermined coefficient and I^iξ is the basis function of the orthogonal polynomial. In the Wiener–Askey scheme introduced by Xiu and Karniadakis, different polynomial types associate with different random variables. The Hermite polynomials are applicable for Gaussian distributions which exhibits a bell shape and symmetric probability density function. For non-Gaussian cases, a possible choice is the transformation technique with the procedure that the nonnormal distributions are transformed into standard normally distributed variables [20]. In equation (2), Npc (the number of terms) is calculated by(3)Npc=p+d!p!d!−1,where p is the order of the polynomial and d is the dimension of the random variables. In complex application, the PCE (equation (2)) can be used to represent the response of uncertain system in nonintrusive formulation, and the coefficients can be determined by the stochastic response surface method (SRSM) [20]. In SRSM, an efficient sampling scheme and regression analysis are employed to evaluate the PCE coefficients. The process of establishing PCE with SRSM is shown in Figure 3.Figure 3
The process of establishing PCE with SRSM.


### 2.2. RBF Neural Network

RBF neural network is a forward network and has a feed-forward architecture consisting of a single hidden layer of locally tuned units that are fully interconnected to an output layer of linear units [21]. In the proposed RDO framework, the AMs based on the RBF is used to estimate the mean and standard deviation of the characteristics of interest. For a given characteristic (denoted as f), both the mean and standard deviation functions are obtained with the same network structure, as shown in Figure 4.Figure 4
RBF network structure.In the RBF network structure, the outputs (μ˜f and σ˜f) are linear combination of Euclidean distances in which the weights are determined by the least squares method. In Figure 4, the terms p and q denote the number of sampling points and the number of hidden neurons, respectively. φi is the i-th basis function determined at the sampling point xi. In general, the following Gaussian function is used as the basis function:(4)ϕiX=exp−X−xiTX−xiri,where ri is the width of the i-th basis function φi.


### 2.3. Multiobjective Optimization

The RDO can be described as the following multiobjective optimization problem [22]:(5)minμfs,x,σfs,xμfs,x=Ezfs,x,Δxσf2s,x=Ezfs,x,Δx−μfs,x,Δx2,s.t.gjs,x,Δx≤0,j=1,2,…,m,where s denotes signal factor. x and Δx denote the control factor (also called as design variable) and its tolerance, respectively. The system performances are random variables due to the uncertainty in the design variable x. For given signal and control factors, the quantities of mean and the standard deviation of system performance may be calculated if the joint probability density function (PDF) of the noise factors is known. In practice, it is usually assumed that all variables have independent normal distributions. In equation (5), the constraints are approximated as the constraints must be satisfied as for all values of control for a worst case. Because of the absolute values, this approximation is likely to be very conservative, and statistical analysis is often used to handle constraints. For a statistical analysis, the constraint is not always satisfied, and the probability that the constraints are satisfied must be chosen a priori. The constraints become(6)μgjx,Δx+k⋅σgjx,Δx≤0,where μgj is the mean of the constraint function gj, σgj is an approximation to the standard deviation of the j-th constraint, and k is a constant that reflects the probability that the constraint will be satisfied. For example, k = 3 means that for a large number of fluctuations, the constraint will be satisfied 99.865 percent of the time, if the total constraint variation is normally distributed.


## 2.1. Polynomial Chaos Expansion

As a newly developed method of uncertainty analysis, PCE has attractive attributes for its strong mathematical rigor and ability to produce functional representations of any stochastic quantities. One important advantage of PCE is that it can conveniently generate the complete probabilistic distribution (e.g., probability density function) of the output response [17]. The roots of polynomial chaos (PC) are found in the homogeneous chaos expansion first proposed by Wiener [18], who employed Hermite polynomials in terms of Gaussian random variables to express stochastic processes with finite variance. An extension was proposed by Xiu and Karniadakis [19] to deal with more general random inputs more efficiently. Suppose that Y (τ) is a random process in the probability space (Ω, F, P), the generalized polynomial chaos expansion of Y (τ) can be expressed as(1)Yτ=α0I0+∑i1=1∞αi1I1ξi1τ+∑i1=1∞∑i2=1i1αi1i2I2ξi1τ,ξi2τ+∑i1=1∞∑i2=1i1∑i3=1i2αi1i2i3I3ξi1τ,ξi2τ,ξi3τ+⋯,where α0,αi1,αi1i2,αi1i2i3… are the undetermined coefficients. The symbol τ indicates that the quantity involved is a random variable (for simplicity, τ will be omitted later). Inξ (n ∈ N, N represents the collection of natural numbers) is an nth-order Wiener–Askey orthogonal polynomial. Inξ is a function of multi-dimensional standard random variables ξ=ξi1,ξi2,…,ξin. The dimension of the random variable ξ is usually consistent with the dimension of the original random input variable X=xi1,xi2,…,xin in the probability space, and there is a one-to-one correspondence between ξ and X. Considering that only limited items can be calculated in practical applications, equation (1) is rewritten as(2)Y=∑i=0Npca^iI^iξ,where α^i is the undetermined coefficient and I^iξ is the basis function of the orthogonal polynomial. In the Wiener–Askey scheme introduced by Xiu and Karniadakis, different polynomial types associate with different random variables. The Hermite polynomials are applicable for Gaussian distributions which exhibits a bell shape and symmetric probability density function. For non-Gaussian cases, a possible choice is the transformation technique with the procedure that the nonnormal distributions are transformed into standard normally distributed variables [20]. In equation (2), Npc (the number of terms) is calculated by(3)Npc=p+d!p!d!−1,where p is the order of the polynomial and d is the dimension of the random variables. In complex application, the PCE (equation (2)) can be used to represent the response of uncertain system in nonintrusive formulation, and the coefficients can be determined by the stochastic response surface method (SRSM) [20]. In SRSM, an efficient sampling scheme and regression analysis are employed to evaluate the PCE coefficients. The process of establishing PCE with SRSM is shown in Figure 3.Figure 3
The process of establishing PCE with SRSM.


## 2.2. RBF Neural Network

RBF neural network is a forward network and has a feed-forward architecture consisting of a single hidden layer of locally tuned units that are fully interconnected to an output layer of linear units [21]. In the proposed RDO framework, the AMs based on the RBF is used to estimate the mean and standard deviation of the characteristics of interest. For a given characteristic (denoted as f), both the mean and standard deviation functions are obtained with the same network structure, as shown in Figure 4.Figure 4
RBF network structure.In the RBF network structure, the outputs (μ˜f and σ˜f) are linear combination of Euclidean distances in which the weights are determined by the least squares method. In Figure 4, the terms p and q denote the number of sampling points and the number of hidden neurons, respectively. φi is the i-th basis function determined at the sampling point xi. In general, the following Gaussian function is used as the basis function:(4)ϕiX=exp−X−xiTX−xiri,where ri is the width of the i-th basis function φi.


## 2.3. Multiobjective Optimization

The RDO can be described as the following multiobjective optimization problem [22]:(5)minμfs,x,σfs,xμfs,x=Ezfs,x,Δxσf2s,x=Ezfs,x,Δx−μfs,x,Δx2,s.t.gjs,x,Δx≤0,j=1,2,…,m,where s denotes signal factor. x and Δx denote the control factor (also called as design variable) and its tolerance, respectively. The system performances are random variables due to the uncertainty in the design variable x. For given signal and control factors, the quantities of mean and the standard deviation of system performance may be calculated if the joint probability density function (PDF) of the noise factors is known. In practice, it is usually assumed that all variables have independent normal distributions. In equation (5), the constraints are approximated as the constraints must be satisfied as for all values of control for a worst case. Because of the absolute values, this approximation is likely to be very conservative, and statistical analysis is often used to handle constraints. For a statistical analysis, the constraint is not always satisfied, and the probability that the constraints are satisfied must be chosen a priori. The constraints become(6)μgjx,Δx+k⋅σgjx,Δx≤0,where μgj is the mean of the constraint function gj, σgj is an approximation to the standard deviation of the j-th constraint, and k is a constant that reflects the probability that the constraint will be satisfied. For example, k = 3 means that for a large number of fluctuations, the constraint will be satisfied 99.865 percent of the time, if the total constraint variation is normally distributed.


## 3. Robust Design Optimization of a Rotor System

The low-pressure rotor of a gas turbine is taken here as an example. The gas turbine low-pressure rotor system has nine-stage compressors and a single-stage turbine with three ball bearings supporting the rotor system. The rotor shaft has an outer diameter of 0.17 m and inner diameter of 0.16 m. The length of the rotor shaft is 2.4 m. The properties of the shaft areE = 81.2 GPa and ρ = 7810 kg/m3. Figure 5 shows a schematic of the low-pressure rotor system.Figure 5
Finite element model of the rotor system.In Figure5, the shaft is discretized into 20 Timoshenko beam elements with four degrees of freedom vwθφT at each node and a constant circular section with 0.12 m length (the translational degrees of freedom v,w in x and y directions, and the rotational degrees of freedom θ, ψ around x-axis and y-axis.). With 20 distinct shaft sections, there are 21 shaft nodes and hence 84 degrees of freedom. The compressors and turbine are treated as concentrated masses. Support_1, support_2, and support_3 indicate the compressor front support, the compressor rear support, and the turbine rear support, respectively. It is assumed that the support stiffness is identical in the vertical and horizontal directions (the effect of the support anisotropy is small because the bearings are quite stiff compared to the rotor). The nominal stiffness at support_1, support_2, and support_3 are 2 × 107 N/m, 1 × 109 N/m, and 1 × 109 N/m, respectively. The damping values in support_1, support_2 and support_3 are 6 × 103 N·s/m, 1 × 103 N·s/m, and 1 × 103 N s/m, respectively. The amount of residual unbalance is 3 × 10−4 kg·m at node 2, node 10, and node 20 shown in Figure 5. Figure 6 is the Campbell diagram of the rotor system.Figure 6
Campbell diagram for the rotor system.It is clear from the Campbell diagram that the rotor has five critical speeds in the range of 0∼13000 rev/min. The values of the two first forward and three first backward critical speeds of the rotor system are summarized in Table1.Table 1
Critical speeds of the rotor system.
Order
Speed (rev/min)
1st backward
3145.2
1st forward
3454.2
2nd backward
6855.9
2nd forward
9981.0
3rd backward
9663.3The mode shapes of the 1st and 2nd forward critical speeds are plotted in Figure7.Figure 7
The mode shapes of the (a) first forward critical speed and (b) second forward critical speed.
(a)
(b)The unbalance responses at support_1, support_2, and support_3 are calculated, and the frequency response functions are given in Figure8. To facilitate comprehension, the first, second, and third backward critical speeds do not appear on the unbalance responses due to the fact that the bearing stiffness is identical in the vertical and horizontal directions.Figure 8
The frequency response functions at support_1, support_2, and support_3.Compared with support_2 and support_3, it is clear from Figure8 that support_1 has the maximum unbalance response within the speed range from 0 to 1300 rev/min. In addition, the unbalanced response at support_1 is quite large when the speed is over 9000 rev/min. In practice, the range of speed of interest is set to 1000∼8500 rev/min, and the maximum unbalance response at support_1 should not exceed 8 mm/s. In the range of 1000∼8500 rev/min, there are nine operating conditions as shown in Table 2.Table 2
The operating speeds.
Operating speeds
Value (rev/min)
ω
1
3000
ω
2
3500
ω
3
4000
ω
4
4500
ω
5
5000
ω
6
5500
ω
7
6000
ω
8
6500
ω
9
7000For the RDO of the present rotor system, the objective function is to minimize the unbalance response and its variance at each operating condition. This can be represented by the P-diagram shown in Figure9.Figure 9
P-diagram for the rotor system.As shown in Figure9, the system responses are random variables due to the uncertainty (Δx) in the design variable x. The mathematical model of robust design optimization is formulated as(7)minμfs,x,σfs,xμfs,x=∑i=19μviωi,x,Δxσfs,x=∑i=19σviωi,x,Δx,s.t.gx,Δx+3⋅σgx,Δx≤8mm/skj∈kjL,kjU,j=1,2,3,Δk1≤0.1k1,Δk2≤0.3k2,Δk3≤0.3k3,where x = [k1, k2, k3] and Δx = [Δk1, Δk2, Δk3]. μvi(·) is the mean of unbalance response at support_1 for operating speed ωi, and g· is the maximum unbalance response at support_1 over the range 1000∼8500 rev/min. The standard deviations σvi(·) and σg(·) are used as the robustness measure of the objective and constraint. The initial ranges of each support stiffness parameter are listed in Table 3.Table 3
Initial ranges of the design variables (unit: N/m).
Design variableskj
Low boundarykjL
High boundarykjU
k
1
1 × 107
3 × 107
k
2
7.5 × 108
1.25 × 109
k
3
7.5 × 108
1.25 × 109


### 3.1. Optimal Latin Hypercube DOE

Latin hypercube design is a way to generate design points that can spread observations evenly over the range of each input variable. For a Latin hypercube design of sizen, the domain of each input variable is divided into n intervals and a set of n design points is chosen in such a way that the projections of design points onto each factor consist of exactly one observation for each interval. The optimal Latin hypercube design (LHD) DOE is selected due to the specific suitability to build approximation models and for the uniform coverage of the input domain. Figure 10 shows the optimal Latin hypercube sample in the initial design space (sample size n = 1000).Figure 10
Optimal Latin hypercube sampling fork1, k2, and k3 (sample size: 1000).According to the optimal Latin hypercube sample, 1000 simulation experiments have been performed based on the numerical model (FEM) in Figure5. Figure 11 shows the main effect plot (the main effect means the effect of a single factor on the system response) for k1.Figure 11
Main effect plot fork1.The ordinateg is the unbalanced response at support_1. It is clear from Figure 11 that the unbalanced response gdecreases first and then increases as k1 increases in the range [1 × 107 N/m, 3 × 107 N/m]. The unbalanced response g varies from 6.71 mm/s to 10.93 mm/s. As the unbalance response at support_1 should not exceed 8 mm/s, the qualified rate is about 33.7%. According to the LHD DOE result, the main effect plot which is used to analyze the influence of each design variable on the maximum unbalance response at support_1 within the initial design range of each input variable is shown in Figure 12.Figure 12
Main effect plots fork1, k2, and k3.It can be seen from the main effect plot that, compared withk2 and k3, the parameter k1 has a greater impact on the maximum unbalanced response at support_1. Therefore, the uncertainty (defined in terms of tolerance) of parameter k1 may be limited to a smaller range for uncertainty analysis of the unbalance response at support_1.


### 3.2. Uncertainty Analysis with PCE Method

The deterministic elementary stiffness matrix of the support_i (i = 1, 2, 3) is defined as(8)Kie=ki0000ki0000000000,where ki is the stiffness of support_i. It is assumed that the uncertain parameters k1, k2, and k3 are independent of each other and obey normal distribution. The randomness of ki is modelled by the truncated Gaussian random variable ξi, which yields(9)ki=k¯i+σki⋅ξi,i=1,2,3,where k¯i is the mean value (nominal value) and σki is the standard deviation. The uncertainty of design variable ki is defined in terms of tolerance (denoted as Δki), and the stiffness parameters k1, k2, and k3 are allowed to undergo 10%, 30%, and 30% variations respectively. Using worst-case formulation, the standard deviations of the stiffness k1, k2, and k3 are σk1=1/3⋅Δk1, σk2=1/3⋅Δk2, and σk3=1/3⋅Δk3. The expression of the assembled stiffness matrix [Ki] is(10)Ki=Ki0+ξi⋅Ki1,where [Ki]0 is the mean of [Ki] and [Ki]1 is its first expansion term. For random variables subject to a normal distribution, the Hermite orthogonal polynomial is chosen as the optimal polynomial basis function, and the optimal polynomial chaos expansion is written as(11)Yτ=∑i=0Npcα^iHiξ,where Y (τ) is the unbalance response. ξ = [ξ1, ξ2, ξ3], and Hi (ξ) is the Hermite orthogonal polynomial basis function. For common situations, the PCE model with order p = 2 or p = 3 has been reported to be accurate enough to represent the stochastic quantity [23]. Therefore, only order p = 2 is considered to trade off between accuracy and efficiency in this paper. From equation (3), Npc = 9 (p = 2, d = 3). The Hermite orthogonal polynomial basis functions are H0 (ξ) = 1, H1 (ξ) = ξ1, H2 (ξ) = ξ2, H3 (ξ) = ξ3, H4 (ξ) = ξ1ξ2, H5 (ξ) = ξ1ξ3, H6 (ξ) = ξ2ξ3, H7 (ξ) = ξ12−1, H8 (ξ) = ξ22−1, and H9 (ξ) = ξ32−1. In the stochastic space (Δk1, Δk2, and Δk3), the SRSM is used to determine the PCE coefficient α^i. It should be noted that the polynomial chaos expansion is associated with the speed of the rotor system. For the initial design (k1 = 2 × 107 N/m, k2 = 1 × 109 N/m, and k3 = 1 × 109 N/m), the obtained PCE associated with ω = 3000 rev/min is(12)Yτω=3000rev/min=1.2556−0.1523ξ1−4.7715×10−3ξ2+6.3096×10−5ξ3+8.9611×10−4ξ1ξ2+4.7545×10−6ξ1ξ3+2.4299×10−7ξ2ξ3+8.9611×10−4ξ12−1+2.5662×10−4ξ22−1−3.1730×10−6ξ32−1,where Y (τ) is the velocity response function (unit: mm/s). Figure 13 shows the uncertainty analysis result of frequency response function at support_1 obtained by conducting Monte Carlo simulation on the PCEs associated with different speed ω (the increment is Δω = 25 rev/min, and the number of simulation is 100).Figure 13
Uncertainty analysis of frequency response function at support_1.For comparison with the Monte Carlo method, Figure14 shows the overlay of Monte Carlo simulation results with 10000 samples (grey solid line) and those from the PCE (red solid line).Figure 14
The overlay of comparison of the variation results obtained with Monte Carlo (grey solid line) and the PCE (red solid line) methods.It can be seen that the results obtained with PCE closely matched with those of traditional Monte Carlo. The time consumption of the PCE method is about 5.8 seconds, while the time consumption of the Monte Carlo method is about 1038.5 seconds (a workstation with 32 GB of RAM and two 6-core processors (Intel Xeon (R) X5690 @ 3.46 GHz)). The integration of the PCE method into a procedure of subsequent optimization with uncertainties will lead to a drastic reduction of the calculation time while preserving an adequate level of precision.


### 3.3. Construction of AMs

The AMs of concerned responses, including the mean and standard deviation of unbalance response at support_1 are constructed with the Isight Optimus software. The applied basis functions in the hidden layer of the RBF neural network are considered to be the Gaussian function. The comparison between AM and simulation results of the maximum unbalance responseμ˜f at support_1 is shown in Figure 15.Figure 15
Comparison between AM (RBF) and simulation results forμ˜f.In Figure15, all sampling points of optimal LHD DOE are taken as inputs and another 100 random sampling points as error analysis points. The statistical coefficient R2 of μ˜f equals to 0.994. For σ˜f and G˜G˜=μ˜+3⋅σ˜g, the statistical coefficients are equal to 0.946 and 0.989, respectively. The error analysis plots for σ˜f and G˜ are omitted here.


### 3.4. Optimization Results and Discussion

The optimization is conducted based on the AMs obtained with RBF neural network. The multiobjective evolutionary algorithm NSGA-II is utilized to obtain the optimization results. In the NSGA-II, the population size is 40, the number of generations is 200, and the crossover probability is 0.9. The feasible solutions obtained are shown in Figure16.Figure 16
Feasible solutions.In Figure16, the red curve is the Pareto set. It is clear that the mean and variance of objective function display opposite changing trends. The Pareto set is displayed in Figure 17.Figure 17
The Pareto set.In the obtained Pareto set, the intervals of stiffness parameterk1, k2, and k3 are [1.538 × 107 N/m, 1.590 × 107 N/m], [9.287 × 108 N/m, 1.134 × 109 N/m], and [1.118 × 109 N/m, 1.229 × 109 N/m], respectively. From the Pareto set, four optimal solutions, RD1, RD2, RD3, and RD4 (shown in Figure 17) are selected for the further discussion. The values of the objective function μf are 10.630 mm/s, 10.650 mm/s, 10.668 mm/s, and 10.697 mm/s, and the corresponding standard deviations σf are 0.926 mm/s, 0.910 mm/s, 0.900 mm/s, and 0.889 mm/s, respectively.The deterministic optimization problems can be established by ignoring the uncertainty of parameters and the robustness of targets and constraints. That is,(13)min∑i=19μviωi,x,s.t.gx≤8mm/skj∈kjL,kjU,j=1,2,3.It can be readily obtained from equation (13) that the deterministic solution (DS) is k1 = 1.618 × 107 N/m, k2 = 7.972 × 108 N/m, and k3 = 8.063 × 108 N/m. For DS and the chosen RDs, the qualified rate is analyzed with Monte Carlo simulation. To obtain convergence of the Monte Carlo method, 10000 simulations were performed, and the results are shown in Table 4.Table 4
The results obtained by the proposed method and Monte Carlo method.
The proposed method
Monte Carlo method
Qualified rate (%)
μ
f (mm/s)
σ
f (mm/s)
μ
f (mm/s)
σ
f (mm/s)
DS
—
—
10.518
0.950
83.6
RD1
10.630
0.926
10.630
0.926
100
RD2
10.650
0.910
10.651
0.912
99.89
RD3
10.668
0.900
10.669
0.903
99.87
RD4
10.697
0.889
10.696
0.889
99.91From Table4, both the means and standard deviations of objective f obtained with the proposed method are consistent with those simulated by the Monte Carlo method. The qualified rates (shown in the fourth column of Table 4) of RD1, RD2, RD3, and RD4 are over 99%. Although the DS has the optimal mean, the standard deviation is larger than the RDs, and the qualified rate is only 83.6%. Figure 18 shows the frequency response functions corresponding to the initial design (k1 = 2 × 107 N/m, k2 = 1 × 109 N/m, and k3 = 1 × 109 N/m) and RD solutions.Figure 18
Initial design and the robust design solutions.In Figure18, the dotted-dashed line is the frequency response function corresponding to the initial design. The critical speed in initial design is 3452.2 rev/min, while the interval of critical speeds obtained by the proposed RDO method is [3225 rev/min, 3250 rev/min], which is more away from the working speeds (ω1 = 3000 rev/min, ω2 = 3500 rev/min). This is critical for safe and stable operation of the rotor system.


## 3.1. Optimal Latin Hypercube DOE

Latin hypercube design is a way to generate design points that can spread observations evenly over the range of each input variable. For a Latin hypercube design of sizen, the domain of each input variable is divided into n intervals and a set of n design points is chosen in such a way that the projections of design points onto each factor consist of exactly one observation for each interval. The optimal Latin hypercube design (LHD) DOE is selected due to the specific suitability to build approximation models and for the uniform coverage of the input domain. Figure 10 shows the optimal Latin hypercube sample in the initial design space (sample size n = 1000).Figure 10
Optimal Latin hypercube sampling fork1, k2, and k3 (sample size: 1000).According to the optimal Latin hypercube sample, 1000 simulation experiments have been performed based on the numerical model (FEM) in Figure5. Figure 11 shows the main effect plot (the main effect means the effect of a single factor on the system response) for k1.Figure 11
Main effect plot fork1.The ordinateg is the unbalanced response at support_1. It is clear from Figure 11 that the unbalanced response gdecreases first and then increases as k1 increases in the range [1 × 107 N/m, 3 × 107 N/m]. The unbalanced response g varies from 6.71 mm/s to 10.93 mm/s. As the unbalance response at support_1 should not exceed 8 mm/s, the qualified rate is about 33.7%. According to the LHD DOE result, the main effect plot which is used to analyze the influence of each design variable on the maximum unbalance response at support_1 within the initial design range of each input variable is shown in Figure 12.Figure 12
Main effect plots fork1, k2, and k3.It can be seen from the main effect plot that, compared withk2 and k3, the parameter k1 has a greater impact on the maximum unbalanced response at support_1. Therefore, the uncertainty (defined in terms of tolerance) of parameter k1 may be limited to a smaller range for uncertainty analysis of the unbalance response at support_1.


## 3.2. Uncertainty Analysis with PCE Method

The deterministic elementary stiffness matrix of the support_i (i = 1, 2, 3) is defined as(8)Kie=ki0000ki0000000000,where ki is the stiffness of support_i. It is assumed that the uncertain parameters k1, k2, and k3 are independent of each other and obey normal distribution. The randomness of ki is modelled by the truncated Gaussian random variable ξi, which yields(9)ki=k¯i+σki⋅ξi,i=1,2,3,where k¯i is the mean value (nominal value) and σki is the standard deviation. The uncertainty of design variable ki is defined in terms of tolerance (denoted as Δki), and the stiffness parameters k1, k2, and k3 are allowed to undergo 10%, 30%, and 30% variations respectively. Using worst-case formulation, the standard deviations of the stiffness k1, k2, and k3 are σk1=1/3⋅Δk1, σk2=1/3⋅Δk2, and σk3=1/3⋅Δk3. The expression of the assembled stiffness matrix [Ki] is(10)Ki=Ki0+ξi⋅Ki1,where [Ki]0 is the mean of [Ki] and [Ki]1 is its first expansion term. For random variables subject to a normal distribution, the Hermite orthogonal polynomial is chosen as the optimal polynomial basis function, and the optimal polynomial chaos expansion is written as(11)Yτ=∑i=0Npcα^iHiξ,where Y (τ) is the unbalance response. ξ = [ξ1, ξ2, ξ3], and Hi (ξ) is the Hermite orthogonal polynomial basis function. For common situations, the PCE model with order p = 2 or p = 3 has been reported to be accurate enough to represent the stochastic quantity [23]. Therefore, only order p = 2 is considered to trade off between accuracy and efficiency in this paper. From equation (3), Npc = 9 (p = 2, d = 3). The Hermite orthogonal polynomial basis functions are H0 (ξ) = 1, H1 (ξ) = ξ1, H2 (ξ) = ξ2, H3 (ξ) = ξ3, H4 (ξ) = ξ1ξ2, H5 (ξ) = ξ1ξ3, H6 (ξ) = ξ2ξ3, H7 (ξ) = ξ12−1, H8 (ξ) = ξ22−1, and H9 (ξ) = ξ32−1. In the stochastic space (Δk1, Δk2, and Δk3), the SRSM is used to determine the PCE coefficient α^i. It should be noted that the polynomial chaos expansion is associated with the speed of the rotor system. For the initial design (k1 = 2 × 107 N/m, k2 = 1 × 109 N/m, and k3 = 1 × 109 N/m), the obtained PCE associated with ω = 3000 rev/min is(12)Yτω=3000rev/min=1.2556−0.1523ξ1−4.7715×10−3ξ2+6.3096×10−5ξ3+8.9611×10−4ξ1ξ2+4.7545×10−6ξ1ξ3+2.4299×10−7ξ2ξ3+8.9611×10−4ξ12−1+2.5662×10−4ξ22−1−3.1730×10−6ξ32−1,where Y (τ) is the velocity response function (unit: mm/s). Figure 13 shows the uncertainty analysis result of frequency response function at support_1 obtained by conducting Monte Carlo simulation on the PCEs associated with different speed ω (the increment is Δω = 25 rev/min, and the number of simulation is 100).Figure 13
Uncertainty analysis of frequency response function at support_1.For comparison with the Monte Carlo method, Figure14 shows the overlay of Monte Carlo simulation results with 10000 samples (grey solid line) and those from the PCE (red solid line).Figure 14
The overlay of comparison of the variation results obtained with Monte Carlo (grey solid line) and the PCE (red solid line) methods.It can be seen that the results obtained with PCE closely matched with those of traditional Monte Carlo. The time consumption of the PCE method is about 5.8 seconds, while the time consumption of the Monte Carlo method is about 1038.5 seconds (a workstation with 32 GB of RAM and two 6-core processors (Intel Xeon (R) X5690 @ 3.46 GHz)). The integration of the PCE method into a procedure of subsequent optimization with uncertainties will lead to a drastic reduction of the calculation time while preserving an adequate level of precision.


## 3.3. Construction of AMs

The AMs of concerned responses, including the mean and standard deviation of unbalance response at support_1 are constructed with the Isight Optimus software. The applied basis functions in the hidden layer of the RBF neural network are considered to be the Gaussian function. The comparison between AM and simulation results of the maximum unbalance responseμ˜f at support_1 is shown in Figure 15.Figure 15
Comparison between AM (RBF) and simulation results forμ˜f.In Figure15, all sampling points of optimal LHD DOE are taken as inputs and another 100 random sampling points as error analysis points. The statistical coefficient R2 of μ˜f equals to 0.994. For σ˜f and G˜G˜=μ˜+3⋅σ˜g, the statistical coefficients are equal to 0.946 and 0.989, respectively. The error analysis plots for σ˜f and G˜ are omitted here.


## 3.4. Optimization Results and Discussion

The optimization is conducted based on the AMs obtained with RBF neural network. The multiobjective evolutionary algorithm NSGA-II is utilized to obtain the optimization results. In the NSGA-II, the population size is 40, the number of generations is 200, and the crossover probability is 0.9. The feasible solutions obtained are shown in Figure16.Figure 16
Feasible solutions.In Figure16, the red curve is the Pareto set. It is clear that the mean and variance of objective function display opposite changing trends. The Pareto set is displayed in Figure 17.Figure 17
The Pareto set.In the obtained Pareto set, the intervals of stiffness parameterk1, k2, and k3 are [1.538 × 107 N/m, 1.590 × 107 N/m], [9.287 × 108 N/m, 1.134 × 109 N/m], and [1.118 × 109 N/m, 1.229 × 109 N/m], respectively. From the Pareto set, four optimal solutions, RD1, RD2, RD3, and RD4 (shown in Figure 17) are selected for the further discussion. The values of the objective function μf are 10.630 mm/s, 10.650 mm/s, 10.668 mm/s, and 10.697 mm/s, and the corresponding standard deviations σf are 0.926 mm/s, 0.910 mm/s, 0.900 mm/s, and 0.889 mm/s, respectively.The deterministic optimization problems can be established by ignoring the uncertainty of parameters and the robustness of targets and constraints. That is,(13)min∑i=19μviωi,x,s.t.gx≤8mm/skj∈kjL,kjU,j=1,2,3.It can be readily obtained from equation (13) that the deterministic solution (DS) is k1 = 1.618 × 107 N/m, k2 = 7.972 × 108 N/m, and k3 = 8.063 × 108 N/m. For DS and the chosen RDs, the qualified rate is analyzed with Monte Carlo simulation. To obtain convergence of the Monte Carlo method, 10000 simulations were performed, and the results are shown in Table 4.Table 4
The results obtained by the proposed method and Monte Carlo method.
The proposed method
Monte Carlo method
Qualified rate (%)
μ
f (mm/s)
σ
f (mm/s)
μ
f (mm/s)
σ
f (mm/s)
DS
—
—
10.518
0.950
83.6
RD1
10.630
0.926
10.630
0.926
100
RD2
10.650
0.910
10.651
0.912
99.89
RD3
10.668
0.900
10.669
0.903
99.87
RD4
10.697
0.889
10.696
0.889
99.91From Table4, both the means and standard deviations of objective f obtained with the proposed method are consistent with those simulated by the Monte Carlo method. The qualified rates (shown in the fourth column of Table 4) of RD1, RD2, RD3, and RD4 are over 99%. Although the DS has the optimal mean, the standard deviation is larger than the RDs, and the qualified rate is only 83.6%. Figure 18 shows the frequency response functions corresponding to the initial design (k1 = 2 × 107 N/m, k2 = 1 × 109 N/m, and k3 = 1 × 109 N/m) and RD solutions.Figure 18
Initial design and the robust design solutions.In Figure18, the dotted-dashed line is the frequency response function corresponding to the initial design. The critical speed in initial design is 3452.2 rev/min, while the interval of critical speeds obtained by the proposed RDO method is [3225 rev/min, 3250 rev/min], which is more away from the working speeds (ω1 = 3000 rev/min, ω2 = 3500 rev/min). This is critical for safe and stable operation of the rotor system.


## 4. Conclusions

Robust design due to parameter uncertainty is a multiobjective and nondeterministic approach. This work concerns the development of a novel methodology for robust design optimization of rotor systems due to its support stiffness uncertainty. The optimization methodology presented in this paper is based on the combined use of simulations, DOE, PCE, RBF, and NSGA-II algorithm. Using the tolerance to define the stiffness uncertainty, the variance of objective and constraint function can be calculated by means of PCE in a reasonable computational time. Under the condition that the maximum unbalance response does not exceed the specified value, the robust optimization results of support stiffness of the rotor system are obtained with the multiobjective optimization model based on the building RBF approximation models of the objective and constraint functions. In practical applications, it is necessary to make a good trade-off among the nominal objective function, uncertainty of the design variables, and reliability of the constraints and thereby provide engineers with a reasonable robust design optimization result.



---

*Source: 2564679-2020-02-11.xml*