Design and Identification Problems of Rotor Bearing Systems Using the Simulated Annealing Algorithm

1. Introduction

The study of rotating machinery appears in the context of machines and structures due to the significant number of phenomena typical to their operation that impact their dynamic behavior and maintenance. Consequently, rotor bearing systems face numerous problems that affect a wide variety of machines, e.g., compressors, pumps, motors, centrifuge machines, large and small turbines. This type of machine finds various applications in the industry, such as, automotive, aerospace and power generation. In most applications an unpredictable stoppage can lead to considerable financial losses and risks. Therefore, there is an evident need for the complete modelling of rotating systems, including the components of the interface between fixed and moveable parts, such as the hydrodynamic bearings. Bench-scale experimental analyses provide more complete models of the main components of the rotor, with strong emphasis on the modelling of the bearings of rotary machines, since they constitute the rotor-foundation structure connecting elements.

The machinery parameters are needed to study the dynamic behavior of the system, namely the Campbell diagram, stability analysis, critical speeds, excitation responses, control and health monitoring. The determination of unknown parameters in rotating machinery is a difficult task. To overcome this difficulty, the use of optimization techniques to solve the inverse problem represents an important alternative approach.

In the literature, various works have been proposed to determine unknown parameters of dynamic systems. Edwards et al. [1] presented a procedure to determine unbalance and support parameters simultaneously based on the least-squares method. Xu et al. [2] proposed a rotor balancing method by using optimization techniques, which does not need trial weights. Assis and Steffen [3] developed strategies in order to use optimization techniques for determining the parameters of gyroscopic systems and they commented about the difficulties that arise in using classical optimization algorithms due to their difficulty in avoiding local minima. The properties of the supports located at the ends of the rotor were considered as variables in the optimization procedure. An inverse problem was developed by using a hybrid cascade-type optimization scheme considering a single unbalance distribution. Castro et al. [4] proposed an optimization method based on genetic algorithms to tune displacements of the rotor supported by hydrodynamic bearings. Castro et al. [5] applied a hybrid algorithm based on genetic algorithm and simulated annealing to tune the orbits of the rotary system in the critical region. In this search algorithm, the genetic algorithm is applied in order to make an approximation of the optimal result, while the simulated annealing refines this result. Tiwari and Chakravarthy [6] presented an identification algorithm for simultaneous estimation of the residual unbalances and the bearing dynamic parameters by using the impulse response measurements for multi-degree-of-freedom flexible rotor-bearing systems. Kim et al. [7] presented a bearing parameter identification of rotor–bearing system using clustering-based hybrid evolutionary algorithm. Castro et al. [8] applied multi-objective genetic algorithm to identify unbalance parameters. Nauclér and Söderstöm [9] consider the problem of unbalance estimation of rotating machinery based on the development of a novel method which takes disturbances into account, leading to a nonlinear estimator. More recently, Saldarriaga et al. [10] proposed a methodology for the experimental determination of the unbalance distribution on highly flexible rotating machinery using Genetic Algorithms. Modal analysis techniques were previously performed to obtain an initial guess for the unknown parameters. A pseudo-random optimization-based approach was used first to identify the parameters of the system in such a way that a reliable rotor model was obtained. Satisfactory results encouraged the use of the proposed approach in the industrial context. Sudhakar and Sekhar [11] proposed a method dedicated to fault identification in a rotor bearing system by minimizing the difference between equivalent loads estimated in the system due to the fault and theoretical fault model loads. This method has a limitation since the error found in the identified fault parameters increases when decreasing the number of measured experimental data.

In this context, the present chapter discusses the possibility of using the Simulated Annealing algorithm (SA) for the identification of unknown parameters of a rotor model from the unbalanced response of the system. Basically, the SA algorithm exploits the analogy between the search for a minimum in the optimization problem and the process of gradual cooling of a metal in a crystalline structure of minimal energy. A desirable characteristic of a minimum search method is the ability to avoid the convergence to a local optimal point, e.g., in terms of the physical process of annealing a meta-stable structure is obtained in the end. Thus, the paradigm of SA is to offer means of escaping from local optima through the analysis of the neighbourhood of the current solution, which can assume, within a given probability, worse solutions, but makes the finding of a new path to the global optimum possible. Metropolis et al. [12] presented an algorithm that simulates the evolution of a crystalline structure in the liquid state up to its thermal equilibrium. Metropolis’ algorithm can be used to generate sequences of configurations in a combinatorial optimization problem. SA is seen as a sequence of Metropolis algorithms, executed with a decreasing sequence of the control parameter. The temperature (control parameter) is continually reduced after a certain number of neighbourhood searches in the current state.

It is worth mentioning that although the SA is a powerful and important optimization tool, often it is not applied according to strict adherence to sufficiency conditions, permitting the researcher to truly claim that the optimal solution has been (statistically) found. According to Ingber [13], the reason typically given is simply that many variants of this technique are considered to be too consuming of resources to be applied in such strict fashion. There exist faster variants of SA canonical, but these apparently are not as quite easily coded and so they are not widely used. Many modifications of SA are really quenching, and should aptly be called simulated quenching (SQ).

In the present contribution, the canonical SA, e.g., based on the algorithm proposed by Kirkpatrick et al. [14]to include a temperature schedule for efficient searching, is used for the design and identification of rotor bearing systems. The goal for the first problem presented is to increase the difference between two critical speeds of a rotor-bearing system that was previously modelled by using the finite element method. In this case, the design variables are the parameters of the rotor-bearing system. To solve this multi-criteria optimization problem a methodology based on a combination of SA, non-dominated sorting strategy and crowding distance operator for guaranteeing convergence and diversity of potential candidates in the population is proposed. The second problem studied is related to the identification of unknown parameters of flexible rotor-bearing systems, modelled mathematically by using the finite element method. The difference between the unbalance experimental responses of the rotor and the simulated unbalance responses (obtained by using the mathematical model) is used to write the objective function to be minimized, so that the damping and stiffness parameters are found. For illustration purposes, the experimental (synthetic) data used were generated by using the solution of the direct problem to which artificial noise was added.

This chapter is organized as follows. The rotor bearing formulation is revisited in Section 4. In Sections 5 and 6 the main characteristics of the SA and multi-objective optimization are briefly presented, respectively. The Multi-objective Optimization Simulated Annealing (MOSA) proposed in this work is described in Section 7. The results and discussion are presented in Section 8. Finally, the conclusions and suggestions for future work complete the chapter.

2. Rotor bearing modelling

The mathematical model used to calculate the unbalance forces, natural frequencies and vibration mode shapes is obtained by using the Finite Element Method. The discrete rotor model is composed of symmetric rigid disc elements, symmetric Timoshenko beam elements, nonsymmetric coupling elements, and nonsymmetric viscous damped bearings, as presented in Figure 1.

Figure 1.

Rotor references frames.

Two reference systems are considered, namely the inertial frame (X,Y,Z) and the frame (x,y,z) that is fixed to the disk [3].By using the Lagrange’s equations in steady-state conditions, the rotor model is represented by the following matrix differential equation [15]:

$$M\ddot{q}+C\dot{q}+Kq=F_1+F_2\sin (\Omega t)+F_3\cos (\Omega t)+F_4\sin (a\Omega t)+F_5\cos (a\Omega t)$$E1

(1)

where q is the N order generalized coordinate displacement vector; K is the stiffness matrix which takes into account the symmetric matrices of the beam and the nonsymmetric matrices of the bearings; C is the matrix containing the antisymmetric matrices due to gyroscopic effects and the nonsymmetric matrices due to bearings viscous damping; F₁ is the constant body force such as gravity; F₂ and F₃ are the forces due to unbalance; F₄ and F₅ are the forces due to the nonsynchronous effect; and a is a constant.

3. Simulated annealing

SA resembles the cooling process of molten metal through annealing (slow cooling process). At high temperature (T), the atoms in the molten metal can move freely with respect to each other, but as the temperature is reduced, the movement of the atoms gets restricted. The atoms start to get arranged and finally form crystals having the minimum possible energy which depends on the cooling rate. If the temperature is reduced at a very fast rate, the crystalline state may not be achieved at all and, instead, the system may end up in a polycrystalline state, which may have a higher energy state than the crystalline state. Therefore, in order to achieve the absolute minimum energy state, the temperature should be reduced at a slow rate [16].

From the optimization point of view, this physical process is analogous to the determination of near-global or global optimum solutions. The energy of the atoms represents the objective function and the final ground state corresponds to the global minimum of the objective function. The analogy between the physical system and the optimization problem is shown in Table 1 [17].

Physical system	Optimization problem
State	Feasible Solution
Energy	Cost Function
Ground state	Optimal solution
Rapid quenching	Local search
Careful annealing	Simulated annealing

Table 1.

Analogy between simulated annealing and optimization.

The basic steps of canonical SA are presented in Figure 2 and described in the following subsections [18].

Figure 2.

Simulated Annealing algorithm flowchart (NI is the number of iterations).

3.1. Initial population

In this iterative technique, an initial guess is randomly generated according to the design space. It should be emphasized that other forms of generating the initial population can be used to initialize the optimization process.

3.2. Initial temperature

The control of the ‘temperature’ parameter must be carefully defined since it controls the acceptance rule defined by the Boltzmann distribution. T has to be large enough to enable the algorithm to move off a local minimum but small enough not to move off a global minimum. According to Chibante et al. [18], the value of T should be defined in an application based approach (ad-hoc) since it is related with the magnitude of the objective function value.

3.3. Perturbation mechanism

This operator permits the creation of new solutions from the current one. In other words it deals with the exploration of the neighbourhood of the current solution by adding small changes to the current solution.

A solution s is defined as a vector s = (x₁,..., x_n) representing a point in the search space. A new solution is generated by using a vector σ = (σ₁,..., σ_n) of standard deviations to create a perturbation from the current solution. A neighbour solution is then produced from the present solution by:

$$x_{i+1}=x_i+N(0,{\sigma }_i)$$E2

(2)

where N(0,σ_i) is a random Gaussian number with zero mean and σ_i standard deviation.

3.4. Temperature update

The most common cooling schedule is the geometric rule for temperature variation:

$$T_{i+1}=\exp \left(\frac{\log \left(sto{p}_{temp}/star{t}_{temp}\right)}{n_{temp}-1}\right)T_i$$E3

(3)

where stop_tempand start_temp are the final temperature (standard deviation) and the initial temperature, respectively, and n_temp is the number of temperatures considered. However, other schedules have been proposed in the literature [19]. Another parameter is the number of iterations for each temperature, which is often related with the size of the search space or with the size of the neighbourhood. This number of iterations can even be constant or, alternatively, can be defined as a function of the temperature or based on a feedback from the process [18].

3.5. Termination criterion

Among the several strategies proposed for the termination of the algorithm, we can cite some very common approaches: the maximum number of iterations; the minimum temperature value; the minimum value of the objective function; the minimum value of the acceptance rate and the maximum computational time.

4. Multi-objective optimization

Real-world design problems involve the simultaneous optimization of two or more (often conflicting) objectives, known as multi-objective optimization problems (MOOP). The solution of such problems is different from the one of the single-objective optimization problems. The main difference is that MOOP normally have not one but a set of solutions, which should be equally satisfactory [20,21].

Traditionally, the treatment of such problems is done by transforming the original MOOP into a scalar single-objective problem. Several studies dealing with multi-objective optimization techniques have been reported over the past decades, based on the Kuhn-Tucker's criterion. These techniques follow the preference-based approach in which a relative preference vector is used to rank multiple objectives. Classical searching and optimization methods use a point-to-point approach, in which the solution is successively modified so that the outcome of the classical optimization method is a single optimized solution. However, Evolutionary Algorithms (EA) can find multiple optimal solutions in one single simulation run due to their population-based search approach. Thus, EA are ideally suited for multi-objective optimization problems.

When dealing with MOOP, the notion of optimality needs to be extended. The most common one in the current literature is that originally proposed by Edgeworth [22] and later generalized by Pareto [23]. This notion is called Edgeworth-Pareto optimality, or simply Pareto optimality, and refers to finding good tradeoffs among all the objectives. This definition leads to a set of solutions that is known as the Pareto optimal set, whose corresponding elements are called non-dominated or non-inferior. The concept of optimality in the single objective context is not directly applicable in MOOPs. For this reason a classification of the solutions is introduced in terms of Pareto optimality, according to the following definitions [20]:

• Definition 1 - The Multi-objective Optimization Problem (MOOP) can be defined as:

$$f\left(x\right)=\left(f_1\left(x\right),f_2\left(x\right),\mathrm{...},f_m\left(x\right)\right),m=\text{1},\mathrm{...},M$$E4

(4)

subject to

$$h\left(x\right)=\left(h_1\left(x\right),h_2\left(x\right),\mathrm{...},h_i\left(x\right)\right),i=\text{1},\mathrm{...},H$$E5

(5)

$$g\left(x\right)=\left(g_1\left(x\right),g_2\left(x\right),\mathrm{...},g_j\left(x\right)\right),j=\text{1},\mathrm{...},J$$E6

(6)

$$x=\left(x_1,x_2,\mathrm{...},x_n\right),n=\text{1},\mathrm{...},N,x\in X$$E7

(7)

where x is the vector of design (or decision) variables, f is the vector of objective functions and X is denoted as the design (or decision) space. The constraints h and g (≥ 0) determine the feasible region.

• Definition 2 - Pareto Dominance: for any two decision vectors u and v, u is said to dominate v, if u is not worse than v in all objectives and u is strictly better than v in at least one objective.

• Definition 3 - Pareto Optimality: when the set P is the entire search space, or P = S, the resulting non-dominated set P’ is called the Pareto-optimal set. Like global and local optimal solutions in the case of single-objective optimization, there could be global and local Pareto-optimal sets in multi-objective optimization.

In the multi-objective context, various Multiple-Objective Evolutionary Algorithms (MOEAs) can be found. This group of algorithms conjugates the basic concepts of dominance described above with the general characteristics of evolutionary algorithms. Basically, the main features of these MOEAs are [20,21]:

• Mechanism of adaptation assignment in terms of dominance - between a non-dominated solution and a dominated one, the algorithm will favour the non-dominated solution. Moreover, when both solutions are equivalent in dominance, the one located in a less crowded area will be favoured. Finally, the extreme points (e.g. the solutions that have the best value in one particular objective) of the non-dominated population are preserved and their adaptation is better than any other non-dominated point, to allow for maximum front expansion.

• Incorporation of elitism - the elitism is commonly implemented using a previously stored secondary population of non-dominated solutions. When performing recombination (selection-crossover-mutation), parents are taken from this file in order to produce the offspring.

In the literature, various multi-objective algorithms based on SA have been proposed. Basically, the first extensions were proposed by Serafini [24,25] and by Ululgu and Teghem [26], where various ways of defining the probability in the multi-objective framework and how they affect the performance of SA based multi-objective algorithms. Czyzak et al. [27] combined mono-criterion SA and genetic algorithm to provide efficient solutions for multi-criteria shortest path problem. Ulungu et al. [28] designed a MOSA (Multi-objective Optimization Simulated Annealing) algorithm and tested its performance using multi-objective combinatorial optimization problems. Suppapitnarm et al. [29] used the neighbourhood perturbation method to create a new point around an old point using MOSA. In this algorithm, the single objective SA is modified to give a set of non-dominated solutions by using archiving of solutions generated earlier, and using a sorting procedure (based on non-dominance and crowding). Kasat et al. [30] used the concept of jumping genes in natural genetics to modify the binary-coded non-dominated sorting genetic algorithm (NSGA-II) to give NSGA-II-JG. Smith et al. [31] compared the candidate to the current solution according to the cardinalities of their dominant subsets in the file. Marcoulaki and Papazoglou [32] proposed a new multiple objective optimization approach by using a Monte Carlo-based algorithm stemmed from SA. Since the expected result in a multiple objective optimization task is usually a set of Pareto-optimal solutions, the optimization problem states assumed here are themselves sets of solutions.

5. Multi-objective optimization simulated annealing – MOSA

Due to the success obtained by the SA in different science and engineering applications, their extension to the multi-objective context is desirable. In this work, the Multi-objective Optimization Simulated Annealing (MOSA) algorithm is proposed. This approach is based on the classical SA associated with the so-called Fast Non-Dominated Sorting operator and has the following structure:

• An initial population of size NP is randomly generated;

• All dominated solutions are removed from the population through the operator Fast Non-Dominated Sorting. In this way, the population is sorted into non-dominated fronts μ_j (sets of vectors that are non-dominated with respect to each other) [20,21];

• Following, SA is applied to generate the new population (potential candidates to solve the MOOP);

• If the number of individuals of the population is larger than a number defined by the user, it is truncated according to the Crowding Distance criterion [20,21].

The steps presented are repeated until a determined stopping criterion is reached. The operators used in the MOSA are described below.

5.1. Fast non-dominated sorting

The so-called Fast Non-Dominated Sorting operator was proposed by Deb et al. [21] in order to sort a population of size N according to the level of non-domination. Each solution must be compared with every other solution in the population to find if the solution is dominated. This requires O(MN) comparisons for each solution, where M is the number of objective functions. When this process is continued to find the members of the first non-dominated class for all population members, the total complexity is O(MN²). At this point, all individuals in the first non-dominated front are found. In order to obtain the individuals in the next front, the solutions of the first front are temporarily discarded and the above procedure is repeated. In the worst case, the task of obtaining the second front also requires O(MN²) computations. The procedure is repeated so that subsequent fronts are found.

5.2. Crowding distance operator

This operator describes the density of solutions surrounding a vector. To compute the Crowding Distance for a set of population members the vectors are sorted according to their objective function value for each objective function. To the vectors with the smallest or largest values, an infinite Crowding Distance (or an arbitrarily large number for practical purposes) is assigned. For all other vectors, the Crowding Distance (dist_xi) is calculated according to [20,21]:

$$dis{t}_{x_i}={\displaystyle \sum _{j=1}^m\frac{f_{j,i+1}-f_{j,i-1}}{\left|f_{j\max }-f_{j\min }\right|}}$$E8

(8)

where f_j corresponds to the j-th objective function and m equals the number of objective functions. This operator is important to avoid many points close together in the Pareto’s Front and to promote the diversity in terms of space objectives [21].

5.3.Consideration. of constraints

In this work, the treatment of constraints is made through the Static Penalization Method, proposed by Castro [33]. This approach consists in assigning limit values to each objective function to play the role of penalization parameters. According to Castro [33], it is guaranteed that any non-dominated solution dominates any solution that violates at least one of the constraints. In the same way, any solution that violates only one constraint will dominate any solution that presents two constraint violations, and so on. For a constrained problem the vector containing the objective functions to be accounted for, is given by:

$$f(x)\equiv f(x)+r_p{n}_{viol}$$E9

(9)

where f(x) it is the vector of objective functions, r_p it is the vector of penalty parameters that depends on the type of problem considered, and n_viol is the number of violated constraints.

6. Applications

6.1. Rotor-dynamics design

Modern design of rotor-bearing systems usually aims at increasing power output and improved overall efficiency. The demanding requirements placed on modern rotating machines, such as turbines, electric motors, electrical generators, compressors, turbo-pumps, have introduced a need for higher speeds and lower vibration levels [34]. This problem can be formulated as a multi-objective problem aiming at minimizing, for instance, the total weight of the shaft, the transmitted forces at the bearings and the positions of the critical speeds [35]. In this context, the present application considers the maximization of the difference between the 6^th and 5^th critical speeds for the system whose finite element model is composed of rigid disks with seventeen elements, two bearings and two additional masses, as shown in Figure 3.

Figure 3.

Finite element model of the rotor-bearing system.

The material used for the shaft and disks is the steel-1020 (density = 7800 Kg/m³, Elasticity modulus = 2.1E11 N/m² and Poisson coefficient = 0.3). The shaft geometry is so that the diameter and length are 10 mm and 552 mm, respectively. The geometric characteristics of the disks are presented in Table 2.

Disc	Mass (Kg)	Moment of Inertia (Kg m^2)	External Diameter (mm)	Thickness (mm)
1	0.818	0.0008	90	16.0
2	1.600	0.0045	150	11.2
3	0.981	0.0018	120	10.6

Table 2.

Geometric characteristics of the disks.

Mathematically, the optimization problem can be formulated as:

$$\min f_1\left(x\right)={\left(RO{T}_{\inf }-a_1{v}_c\left(5\right)\right)}^2$$E10

(10)

$$\min f_2\left(x\right)={\left(v_c\left(6\right)-a_{2}RO{T}_{\sup }\right)}^2$$E11

(11)

where v_c is the critical speeds vector, ROT_i are the permissible rotations (i=inf or sup), a₁=1.3 and a₂=1.3.

For evaluating the methodology proposed in this work, some practical points regarding the application of this procedure should be emphasized:

1. The design variables are the following: radius of bar elements (x_i), where the design space is given by: 0.4 mm ≤ x_i ≤ 0.8 mm.

2. ROT_inf = 1400 Hz and ROT_inf = 1900 Hz.

3. To solve the optimization problem the following heuristics are used:

4. Non-dominated Sorting Genetic Algorithm (NSGAII) parameters [20,36]: population size (50), crossover probability (0.8), mutation probability (0.01). For the considered parameters, the number of objective function evaluations is 12550.

5. Multi-objective Optimization Differential Evolution (MODE) parameters [37]: population size (50), perturbation rate (0.8), crossover probability (0.8), DE/rand/1/bin strategy for the generation of potential candidates, reduction rate (0.9) and number of pseudo-curves (10). For the considered parameters, the number of objective function evaluations is 15050.

6. Multi-objective Optimization Simulated Annealing (MOSA) parameters [14]: population size (50), initial temperature (5.0), cooling rate (0.75), number of temperatures (20), number of times the procedure is repeated before the temperature is reduced (25), and tolerance (10^-6). For the considered parameters, the number of objective function evaluations is 12550.

7. Stopping criterion: maximum number of generations (250).

8. Each algorithm was run 20 times by using 20 different seeds for the random generation of the initial population.

9. Objective Function (OF_A) is the best value of objective function considering the first objective proposed. Objective Function (OF_B) is the best value of objective function considering the second objective function. Objective Function (OF_C) is calculated using the origin of the coordinated axes as a reference, e.g., the point (0,0) is used to obtain the distance between this point and each one of the solution points along the Pareto's Front. Thus, the smallest distance obtained was defined as the choice criterion.

Figure 4 shows the Pareto’s Front obtained by NSGA II, MODE and MOSA algorithms.

Figure 4.

Pareto’s Front.

In this figure it is possible to observe that all evolutionary algorithms are able to obtain, satisfactorily, the Pareto’s Front for a similar number of objective function evaluations.

Table 3 present some points of Pareto’s Front obtained by the MOSA algorithm by considering the criteria specified earlier.

6.2. Identification

As mentioned earlier, the identification of unknown parameters in rotating machinery is a difficult task and optimization techniques represent an important alternative for this goal [3,38,39]. The machine parameters are needed to perform the dynamic analysis and prediction of rotor-bearing systems: Campbell diagram, stability, critical speeds, excitation responses [15]. Another important aspect is when one desires to tune a finite element model to match experimental data generated by tests of an actual rotor system [10].

	OFA	OFB	OFC
x1 (mm)	0.55186	0.40015	0.50692
x2 (mm)	0.49994	0.42153	0.52587
x3 (mm)	0.53139	0.40623	0.50733
x4 (mm)	0.79128	0.40000	0.50851
x5 (mm)	0.49810	0.40000	0.50726
x6 (mm)	0.55709	0.40000	0.50025
x7 (mm)	0.41059	0.40000	0.50811
x8 (mm)	0.71018	0.40000	0.50683
x9 (mm)	0.57777	0.40151	0.50713
x10 (mm)	0.44676	0.41028	0.50192
x11 (mm)	0.79883	0.40000	0.49997
x12 (mm)	0.67159	0.40000	0.50671
x13 (mm)	0.57684	0.40000	0.50592
x14 (mm)	0.46109	0.40000	0.51012
x15 (mm)	0.52452	0.40000	0.50257
x16 (mm)	0.61355	0.40399	0.50747
x17 (mm)	0.57261	0.40281	0.50533
f1 (Hz^2)	0.00364	1110741.5	479693.69
f2 (Hz^2)	1842938.0	225219.04	602024.24

Table 3.

Results obtained using MOSA (all the algorithms were executed 20 times so that average values were calculated).

Furthermore, identification procedures try to establish an unequivocal relation in between the damage and specific mechanical parameters, based on a suitable model and can be used to fault detection and machinery diagnosis as in Seibold and Fritzen [40]. On a simple manner parameter identification of rotor-bearing systems can be performed as follows: i) the frequency response function (or unbalance response) is measured for different operation speeds; ii) the design variables (unknown parameters) are initialized; iii) and an error function between experimental and simulated data is minimized.

Figure 5.

Rotor system finite element model.

In this application, a simple flexible rotor containing two disks and two bearings is studied. The Figure 5 shows the finite element model of the system with 10 nodes, 2 disks and two bearings.

The characteristics of bearing and disks are given in Table 4. It can be observed that damping parameters are also taken into account in this application.

Bearing	Stiffness(N/m)			Damping (Nm/s)
	$$k_{xx}$$	$$k_{zz}$$	$$k_{xz}$$	$$c_{xx}$$	$$c_{zz}$$	$$c_{xz}$$
A=B	1E06	2E06	1E04	1E03	2E03	1E02
Disk	External Diameter (m)			Thickness (m)
1	0.5			0.005
2	0.3			0.005

Table 4.

Parameters of bearings and disks.

The goal of this application is to identify the unknown parameters of the rotor-bearing system, e.g., the stiffness and damping parameters. For this purpose the following steps are established:

1. The objective function consists in the determination of the stiffness and damping values through the minimization of the difference between the experimental and calculated values given by the solution of the direct problem. To mimic real experimental data, sets of synthetic experimental data were generated from eq. (12):

$${\theta }_i^{\exp }={\theta }_i^{cal}\left(Z_{exact}\right)+\kappa \lambda$$E12

(12)

were θ represents the calculated values of the unbalance response by using known values of the physical properties Z_exact (k_xx, k_zz, k_xz, c_xx, c_zz, c_xz). In real applications these values are not available (they can be obtained through the solution of the corresponding inverse problem). κ simulate the standard deviation of the measurement errors, and λ is a pseudo-random number from the interval [-1, 1].

$${\theta }_i^{\exp }={\theta }_i^{cal}\left(Z_{exact}\right)+\kappa \lambda$$E13

(13)

In order to examine the accuracy of the inverse problem approach for the estimation of the physical parameters, the influence of noise (κ =0.02, e.g., corresponding to 5% error) was compared to the case without noise (κ =0).

The design variables considered to generate the synthetic experimental data are presented in Table 5. The following ranges for the design space are considered: 5.0E05 N/m ≤ k_xx ≤ 1.2E07 N/m, 1.0E06 N/m ≤ k_zz ≤ 2.4E07 N/m, 5.0E03 N/m ≤ k_xz ≤ 1.2E05 N/m, 5.0E02 Nm/s ≤ c_xx ≤ 1.2E04 Nm/s, 1E03 Nm/s ≤ c_zz ≤ 2.4E04 Nm/s, 1E01 Nm/s ≤ c_xz ≤ 1.2E03 Nm/s.

1. To solve the optimization problem the following heuristics are used:

2. Genetic Algorithm (GA) parameters [41]: population size (50), type of selection (normal geometric in the range [0 0.08]), type of crossover Arithmetic, 2), type of mutation (non-uniform [2 100 3]).

3. Differential Evolution (DE) parameters [42]: population size (25), perturbation rate and crossover probability both equal to 0.8 and DE/rand/1/bin strategy.

4. Particle Swarm Optimization (PS) parameters [43]: population size (25), maximum velocity (100), upper limits (2.0), and a linearly decreasing inertia weight starting at 0.7 and ending at 0.4 was used.

5. Simulated Annealing (SA) parameters [14]: initial design (generated randomly in the design space), initial temperature (5.0), cooling rate (0.75), number of temperatures (20), number of times the procedure is repeated before the temperature is reduced (25), and tolerance (10^-6).

6. Stopping criterion: maximum number of objective function evaluations equal to 5000.

7. Each algorithm was run 20 times by using 20 different seeds for the random generation of the initial population.

Table 5 presents the results obtained by the algorithms considered (pristine condition and noisy data).

Considering κ=0 (see Table 5), all the optimization strategies were able to estimate the parameters satisfactorily as shown by the values obtained for the objective function. However, the SA algorithm shows to be very competitive, in averege, with the smallest standard deviation of the objective function. When noise is taken into account (κ=0.002, e.g., error corresponding to 5%), all the algorithms were able to obtain good estimates, as presented in Figure 6.

	Error	kxx; (N/m)	kzz; (N/m)	kxz; (N/m)	cxx; (Nm/s)	czz; (Nm/s)	cxz (Nm/s)	OF
SA	0 %	984674.3	1802257.4	24197.8	1157.4	1953.4	168.4	4.38
5 %	962861.5	1817105.5	18378.3	1046.6	2843.4	35.8	7.01
GA	0 %	998204.7	2030908.5	3589.1	1078.8	2068.1	149.7	4.42
5 %	986445.7	2114381.8	8695.9	1010.5	2164.4	25.8	07.11
DE	0 %	1000380.9	2005635.1	17449.8	1102.4	2092.8	183.7	4.35
5 %	970336.4	1964371.6	46943.8	1037.2	2347.5	84.5	6.96
PS	0 %	987254.1	2098581.6	57858.5	1079.1	1425.6	119.8	4.40
5 %	990895.2	2042193.6	15630.8	1118.9	1937.8	189.7	7.09

Table 5.

Estimation Results.

Figure 6.

Boxplots showing the influence of different optimizaion strategies to solve the inverse problem.

7. Conclusions

In the present contribution, the mono and multi-objective algorithms based on Simulated Annealing were used in the design and identification of rotor bearing systems. For illustration purposes, two simple test-cases were studied by using the proposed methodology. The goal for the first application was to increase the difference between two critical speeds of the rotor-bearing system through the formulation of a multi-objective problem, where the radii of bar elements were taken as design variables. To solve this multi-objective problem the Multi-objective Optimization Simulated Annealing (MOSA) algorithm was proposed. This evolutionary strategy is based on the Simulated Annealing algorithm associated with the non-dominated sorting and crowding distance operators. The second application consists in the identification of unknown parameters of flexible rotor-bearing systems. The objective function was defined as the difference between the unbalance experimental responses of the rotor and the simulated unbalance responses so that the parameters of damping and stiffness are obtained by an inverse problem approach. The experimental (synthetic) data used were generated by using the solution of the direct problem and adding artificial noise. In all applications, the finite element method was used to obtain the mathematical model of the system.

It is important to emphasize that the results obtained in both test-cases are considered satisfactory as compared with those obtained by other evolutionary strategies. In addition, it is possible to conclude that the proposed methodology represents an interesting alternative for design and identification of mechanical systems.

Further research work will be focused on the influence of the optimization parameter values on the solution of the optimization problem. Also, strategies to dynamically update the SA parameters will be evaluated. Finally, the authors will study the performance of the Simulated Quenching algorithm aiming at proposing a hybrid approach involving the Simulated Annealing and Simulated Quenching algorithms.

Acknowledgement

The authors acknowledge the financial support provided by FAPEMIG and CNPq (INCT-EIE). The fourth author is grateful to the financial support provided by CNPq and FAPERJ.