A Quasiphysics Intelligent Model for a Long Range Fast Tool Servo

Abstract

Accurately modeling the dynamic behaviors of fast tool servo (FTS) is one of the key issues in the ultraprecision positioning of the cutting tool. Herein, a quasiphysics intelligent model (QPIM) integrating a linear physics model (LPM) and a radial basis function (RBF) based neural model (NM) is developed to accurately describe the dynamic behaviors of a voice coil motor (VCM) actuated long range fast tool servo (LFTS). To identify the parameters of the LPM, a novel Opposition-based Self-adaptive Replacement Differential Evolution (OSaRDE) algorithm is proposed which has been proved to have a faster convergence mechanism without compromising with the quality of solution and outperform than similar evolution algorithms taken for consideration. The modeling errors of the LPM and the QPIM are investigated by experiments. The modeling error of the LPM presents an obvious trend component which is about ±1.15% of the full span range verifying the efficiency of the proposed OSaRDE algorithm for system identification. As for the QPIM, the trend component in the residual error of LPM can be well suppressed, and the error of the QPIM maintains noise level. All the results verify the efficiency and superiority of the proposed modeling and identification approaches.

1. Introduction

In recent years, the fast tool servo (FTS) based single point diamond turning (SPDT) has been gaining increasing applications in extensive areas such as dynamic compensation of machining errors, generation of nonaxisymmetric optical and components, fabrication of microstructured functional surfaces [1–3]. The determination of trajectory tracking strategy for the cutting tool is an ongoing topic for the FTS system to achieve excellent control performance with high accuracy and high robustness, accordingly to achieve fine profile accuracy and high quality of finished surface.

Up to now, various trajectory tracking strategies have been developed. They could be mainly categorized into model-based control approaches relying on the models of the controlled systems and adaptive control approaches with no dependences on the models. As for the adaptive control approaches, the slide mode controller (SMC) has been mainly implemented for the tool positioning of FTS [4, 5]. The SMC does not need accurate model of controlled system and has relatively high robustness to external disturbances and model uncertainties or variations, but only at the expense of positioning accuracy. The model-based controllers which have been applied to FTS generally comprise of proportional integral derivate (PID) controller [6–8], adaptive feedforward cancellation (AFC) controller [3, 9, 10], repetitive controller (RC) [11, 12], and active disturbance rejection controller (ADRC) [13, 14]. To further enhance the robustness and the stability of FTS system, extended state observer (ESO) based compensation approach serving as one of the key components of the ADRC has been extensively embedded in both model-based controllers and adaptive controllers [10, 11, 13–15]. As for the optimal design process of the PID controller, the feedforward component and phase advance parameters of AFC controller, and the extensively employed ESO compensator, accurate model for describing the physical process of FTS should be critically established. Unfortunately, almost all of the FTSs are approximately equivalent to linear second-order dynamic systems [8–10, 12–14], even for the piezoelectric actuated FTS with inherent hysteresis and creep nonlinearity (HCN) effects [4, 15]. It is evident that there should exist inherent nonlinearities in any practical FTS system, the linear approximations of these servo systems would lead to significant modeling errors and consequently deteriorate control performances of FTS systems. Hence, considerable work should be further done to improve the modeling accuracy.

Generally, the system modeling approaches can be mainly divided into two categories, which are physical approach and phenomenological approach, respectively. On one hand, the description capability of physical model is limited by its assumptions. On the other hand, phenomenological model has its own drawbacks, such as outlier and overfitting. Such problems probably lead to conclusions against physical laws, especially when phenomenological model is used for extrapolation [16–19]. So it is believed that integrating these models has the potential to overcome the drawbacks of both physical and phenomenological models [16–19].

Motivated by this, a quasiphysics intelligent model (QPIM) integrating a linear physics model (LPM) and a radial basis function (RBF) based neural model (NM) is proposed to accurately describe the dynamic behaviors of a long range FTS (LFTS) driven by a voice coil motor (VCM), and a novel Opposition-based Self-adaptive Replacement Differential Evolution (OSaRDE) algorithm is developed for the parameters identification of the LPM governed by physics laws. The remainder of this paper is organised as follows: the mechanical structure and the modeling procedure of the LFTS is introduced in Section 2; the preliminary of the differential evolution (DE) algorithm and the proposed OSaRDE algorithm is introduced in Section 3; the experiment results and corresponding discussions are carried out in Section 4; the main conclusions of this paper are drawn in Section 5.

2. Modeling and Identification Procedure of the LFTS

2.1. Mechanical Structure of the LFTS

The authors of this paper have developed a LFTS for ultraprecision diamond turning; the mechanical structure of the LFTS is illustrated in Figure 1(a). The LFTS mainly contains two parts: the VCM actuator and the flexure-based guide mechanism. A BEI Kimco VCM is employed as the driving element. It can provide maximum continuous stall forces of 82.74 N (18.6 lbs) with a stroke of ±5.72 mm. Being benefited from the long stroke of the actuator, the stroke of the LFTS could reach up to 1 mm. Two groups of cross shape flexure hinges are specially designed for the guide of the moving component of the LFTS. The two hinges are connected by a linking member crossing the VCM stator. The mechanical design process and performance evaluations of this LFTS are detailed in [20].

(a) LFTS mechanical structure

(b) Dynamic model of the LFTS mechanism

Figure 1 

Schematic of the LFTS mechanism.

2.2. The Quasiphysics Intelligent Model of the LFTS

From an electromagnetic point of view, the VCM can be considered as an inductance component with a resistance [21, 22]. Meanwhile, from the dynamic point of view, the flexure-based mechanism can be equivalent as a damped mass-spring system [8, 23]. With consideration of the mass of the moving part including the flexure hinges and the bobbin of the VCM, the dynamic model of the LFTS mechanism can be obtained, as illustrated in Figure 1(b). The electromagnetic dynamic balance of the VCM in the LFTS, as equivalently shown in Figure 1(b), can be derived based on the Kirchhoff’s law, yielding [21, 22] 𝑉_𝑚=𝑅_𝑚𝑖_𝑚+𝐿_𝑚𝑑𝑖_𝑚𝑑𝑡+𝑉_𝑚𝑏=𝑅_𝑚𝑖_𝑚+𝐿_𝑚𝑑𝑖_𝑚𝑑𝑡+𝐾_𝑚𝑣𝑠𝑑𝑥,𝑑𝑡 (1) where 𝑉_𝑚 is the voltage applied to the VCM, 𝑉_𝑚𝑏 is the back electromotive force (EMF), and 𝐾_𝑚𝑣𝑠 is the equivalent back EMF constant. 𝑅_𝑚, 𝐿_𝑚, and 𝑖_𝑚 represent the equivalent resistance of voice coils, inductance and applied current of coil, respectively.

With the electrical dynamics in (1) derived, the Flemining’s left-hand rule is then employed to derive the electromagnetic forces for actuation, which can be expressed as [21, 22] 𝐹=𝑛_𝑚𝐵_𝑚𝑙_𝑚𝑖_𝑚=𝑛_𝑚𝐵_𝑚𝑙_𝑚𝑅_𝑚𝑉_𝑚−𝐿_𝑚𝑑𝑖_𝑚𝑑𝑡−𝐾_𝑚𝑣𝑠𝑑𝑥,𝑑𝑡 (2) where 𝑛_𝑚 is the number of coil loops, 𝐵_𝑚 is the magnetic flux densities within the air gap between bobbin and magnets, and 𝑙_𝑚 is the effective coil lengths considering all loops. Generally, the term 𝐿_𝑚𝑑𝑖_𝑚/𝑑𝑡 can be considered negligible due to the fact that the electrical dynamics of the conducted current is much faster than mechanical behaviors; thus 𝐹=𝑛_𝑚𝐵_𝑚𝑙_𝑚𝑖_𝑚=𝑛_𝑚𝐵_𝑚𝑙_𝑚𝑅_𝑚𝑉_𝑚−𝐾_𝑚𝑣𝑠𝑑𝑥.𝑑𝑡 (3)

Based on Newton’s second law of motion, the differential equation of dynamic motion for the LFTS mechanism can be given as follows: 𝑚_𝑓ℎ+𝑚_VCM𝑑²𝑥𝑑𝑡²+𝑐_𝑓ℎ𝑑𝑥𝑑𝑡+𝑘_𝑓ℎ𝑥=𝐹, (4) where 𝑚_𝑓ℎ and 𝑚_VCM represent the equivalent moving mass of the flexure-based mechanism and the bobbin of the VCM, respectively. 𝑐_𝑓ℎ and 𝑘_𝑓ℎ represent the equivalent damping coefficient and the equivalent stiffness of the flexure-based mechanism, respectively.

Submitting (3) into (4) yields 𝑚_𝑓ℎ+𝑚_VCM𝑑²𝑥𝑑𝑡²+𝑐_𝑓ℎ+𝐾_𝑚𝑣𝑠𝑛_𝑚𝐵_𝑚𝑙_𝑚𝑅_𝑚𝑑𝑥𝑑𝑡+𝑘_𝑓ℎ𝑥=𝑛_𝑚𝐵_𝑚𝑙_𝑚𝑅_𝑚𝑉_𝑚. (5)

If this mechanical system is a linear time-invariant (LTI) system, the transfer function of this LFTS system can be obtained as follows by applying the Laplace transform: 𝑃(𝑠)=𝑋(𝑠)𝑉(𝑠)=𝐾𝑠²+2𝜁𝜔_𝑛𝑠+𝜔²_𝑛, (6) where 𝐾=𝑛_𝑚𝐵_𝑚𝑙_𝑚/(𝑚_𝑓ℎ+𝑚_VCM)𝑅_𝑚, 2𝜁𝜔_𝑛=(𝑐_𝑓ℎ𝑅_𝑚+𝑛_𝑚𝐵_𝑚𝑙_𝑚𝐾_𝑚𝑣𝑠)/(𝑚_𝑓ℎ+𝑚_VCM)𝑅_𝑚, 𝜔²_𝑛=𝑘_𝑓ℎ/(𝑚_𝑓ℎ+𝑚_VCM).

Unfortunately, due to manufacturing tolerances on various crucial dimensions of both bobbin and magnets and possible misalignments of magnetic components, there often exist an uneven magnetic filed within the air gap, as a result, it would lead to a more complex and position-varying magnetic flux density 𝐵_𝑚 [21, 22]. On the other hand, the temperature in the VCM would increase with the motion ongoing, consequently leads to the variations of the equivalent resistance and inductance of the VCM. Meanwhile, the equivalent damping coefficient and stiffness of the flexure-based mechanism could not be permanent constants due to the inherent properties of the employed material for the mechanical structure. Let 𝑥^∗(𝑡) be the real output of the LFTS; it would be difficult for the description of 𝑥^∗(𝑡) due to the presence of noise (uncontrollable) factors, measurement error, and inherent nonlinearities of system. Hence, all the factors making the system to be nonlinear and random would significantly deteriorate the accuracy of the LTI model.

For nonlinear dynamic modelling, neural model has been an extremely popular choice due to its high adaptability to various nonlinear systems [24–27]. In the feedforward neural networks, RBF neural networks are local neural networks with series advantages of the faster training speed, small computing, and simple structure [27]. To enhance the accuracy of the pure physics model as obtained in (6), an RBF neural network (RBFNN) based NM is further employed to model the components with nonlinearities and uncertainties. By introducing this NM, the governing equation of the QPIM for the LFTS mechanism in the time domain can be modified as 𝑥^∗(𝑡)=𝐿⁻¹(𝑋(𝑠))+𝑁𝑁𝑉_𝑚,𝑉̇_𝑚, (7) where 𝐿⁻¹(⋅) represent the inverse Laplace transform operation and 𝑁𝑁(⋅) represent the RBF based NM. The command signal and its differentiation are employed as the input of the NM, which are three-layer feedforward networks consisting of an input layer, a hidden layer, and an output layer. Generally, the Gaussian function is employed as the basis function of the hidden layer, given by [27] 𝐅_𝑖−𝐱=exp^𝑖−𝐜_𝑖²2𝜎²_𝑖, (8) where 𝐱_𝑖, 𝐜_𝑖, and 𝜎_𝑖 are the input vector, centre vector, and width of the Gaussian, respectively. ‖⋅‖ denotes the Euclidean norm operation.

The output of the NM is a linear combiner, which is defined by [27] 𝑁𝑁𝑉_𝑚,𝑉̇_𝑚=^𝑛𝑐_𝑖=1𝑤_𝑖⋅Φ_𝑖𝑉_𝑚,𝑉̇_𝑚, (9) where 𝑤_𝑖 is the weight between the hidden layer and the output layer; 𝑛_𝑐 is the number of nodes in the hidden layer.

2.3. Parameters Identification Procedure

The basic idea in the LFTS system identification is to compare the input dependent response of the system and a corresponding mathematic model by certain performance criterions giving a measure to how well the model response fits the system response. The principle of the identification procedure for the parameter identification of the LFTS system is illustrated in Figure 2, where 𝑥(𝑡) denotes the response of the LFTS system; 𝑥₁(𝑡) and 𝑥₂(𝑡) denote the response of the LPM and the NM, respectively. As discussed above, the LPM depends on a set of three parameters, that is, Γ=[𝜁,𝜔_𝑛,𝐾]∈ℜ³. Hence, the objective of the LPM identification is to find a set of parameters Γ^∗=[𝜁^∗,𝜔^∗_𝑛,𝐾^∗]∈ℜ³ that minimize the difference between the response of the real LFTS system and the output of the LPM at each discrete sampling point. Therefore, the identification process lies on minimizing the following criterion: 𝐸(Γ)=^𝑁_𝑖=1𝑥₁𝑡_𝑖−𝑥𝑡_𝑖, (10) where 𝑖 represents the index of the 𝑖th discrete sampling point.

Figure 2 

The principle of system identification.

The criterion 𝐸(⋅) is commonly called a fitness function or objective function reflecting the goodness of solution in the parameter tuning process. The identification problem thus can be treated as a linearly constrained multidimensional optimization problem and formalized as minimize𝐸(Γ)=^𝑁_𝑖=1𝑥₁Γ,𝑡_𝑖−𝑥𝑡_𝑖,Γ=𝜁,𝜔_𝑛,𝐾s.t.Γ∈𝐁,𝐁=Γ,Γ_min,𝑗<Γ_𝑗<Γ_max,𝑗,∀𝑗=1,2,3, (11) where 𝐁 is the 3-dimensional feasible search space; Γ_min,𝑗 and Γ_max,𝑗 denote the upper bounds and the lower bounds of the 𝑗th parameter Γ_𝑗, respectively.

𝐸(Γ) serving as objective function maps decision variables into the objective space. Obviously, the fitness landscape of this identification task may have many local optima and a highly complex topology. Hence, an improved OSaRDE algorithm, which is a simple yet efficient heuristic approach, is employed as the searching tool and will be detailed in the following section.

As for the NM, the similar minimization criterion is defined as 𝑁_𝐸(𝐰)=1𝑁^𝑁_𝑖=1𝐸Γ^∗,𝑡_𝑖−𝑥₂𝐰,𝑡_𝑖. (12)

For the training process, the command signal and its differentiation are employed as the input signal, while the residual error of the LPM is employed as the target signal. To obtain the optimal weights of the RBFNN, the developed OSaRDE algorithm could also be implemented as shown in [26]. However, the training of the NM to obtain the optimal weights and structure is executed through the NEURB function in Matlab just for simplicity. Following the scheme provided by Matlab, new neuron is created at each iteration, and the criterion of the new network is checked. If it is low enough, the training is finished; otherwise, the next neuron is added. This procedure is repeated until the criterion goal is met or the maximum number of iteration is reached [27].

3. The Opposition-Based Self-Adaptive Replacement Differential Evolution Algorithm

The DE algorithm proposed by Storn and Price (1997) is a heuristic approach for minimizing or maximizing possibly nonlinear and nondifferentiable continuous space functions [28]. Due to some attractive characteristics such as simple differential operator, one-to-one competition scheme, and constructive cooperation between current individuals and memory of the best individuals, DE has been extensively implemented for the identification of both linear and nonlinear systems [29–31]. The obtained results demonstrate that DE is a very promising tool for optimal parameters determination and outperforms other evolutionary technologies [30–32].

3.1. The Classical DE Algorithm

The DE algorithm is a population-based algorithm like general evolution algorithms using the similar operators: crossover, mutation, and selection. In DE, a population of potential solution vectors is initialized covering the entire parameter space at the start, which is then evolved to find optimal solutions through the mutation, crossover, and selecting operation procedures [28, 31]. More specifically DE’s basic strategy can be formalized as follows.

(S1)  Initialization. Assuming that an optimization task consists of 𝐷 parameters, let 𝐁∈ℜ^𝐷 be the search space of the problem under consideration; the initial vector population randomly covering the entire parameter space can be expressed as 𝐗⁰_𝑖=𝑥⁰_𝑖,1,𝑥⁰_𝑖,2,𝑥⁰_𝑖,3,…,𝑥⁰_𝑖,𝐷∈𝐁𝑖=1,2,…,NP, (13) where NP denotes the population size;𝐗⁰_𝑖 denotes a potential solution for the task at start.

(S2)  Mutation. For each target vector 𝐗^𝑘_𝑖 at the 𝑘th generation, an associated mutant vector 𝐗^𝑘_𝑖=(𝑥^𝑘_𝑖,1,𝑥^𝑘_𝑖,2,𝑥^𝑘_𝑖,3,…,𝑥^𝑘_𝑖,𝐷) should be generated via a certain mutation operator. The objective of mutation is to enable search diversity in the parameter space as well as to direct the existing object vectors to better results at a suitable time. Motivated by this, various mutation strategies have been proposed, and certain typical strategies can be summarized as [33, 34] DE/rand/1∶𝐗^𝑘_𝑖=𝐗^𝑘_𝑟1+𝐹𝐗^𝑘_𝑟2−𝐗^𝑘_𝑟3,DE/best/1∶𝐗^𝑘_𝑖=𝐩^𝑘_best+𝐹𝐗^𝑘_𝑟1−𝐗^𝑘_𝑟2,DE/cur-to-best/1∶𝐗^𝑘_𝑖=𝐗^𝑘_𝑖+𝐹𝐩^𝑘_best−𝐗^𝑘_𝑖+𝐹𝐗^𝑘_𝑟1−𝐗^𝑘_𝑟2, (14)(15)(16) where 𝑟1, 𝑟2, and 𝑟3 represent the random and mutually different integers generated within the range [1,NP], and also different from index 𝐼; 𝐹 is a mutation scale factor within the range [0,2], usually less than 1; 𝐩^𝑘_best denotes the best individual in generation 𝑘.

(S3)  Crossover. After the mutation phase, the crossover operation is applied to each pair of the generated mutant vector 𝐗^𝑘_𝑖 and its corresponding target vector 𝐗^𝑘_𝑖 to generate a trial vector 𝐘^𝑘_𝑖 as follows: 𝐘^𝑘_𝑖=𝑦^𝑘_𝑖,1,𝑦^𝑘_𝑖,2,𝑦^𝑘_𝑖,3,…,𝑦^𝑘_𝑖,𝐷,𝑦^𝑘_𝑖,𝑗=𝑥^𝑘_𝑖,𝑗𝑗(𝑖),𝑥,ifrand≤CRor𝑗=rand𝑛_{𝑘𝑖,𝑗},ifrand𝑗>CRor𝑗≠rand𝑛(𝑖), (17) where 𝑗=1,2,…,𝐷; 𝑖=1,2,…,NP; rand(𝑗) is the 𝑗th independent random number uniformly distributed in the range of [0,1]. rand𝑛(𝑖) is a randomly chosen index from the set {1,2,…,𝐷}; CR is a user-specified crossover factor within [0,1] that controls the diversity of the population [28, 31].

(S4)  Selection. DE employs a greedy selection process that the better one of new offspring and its parent win the competition. If 𝐹_𝐶(⋅) denotes the objective function under consideration, the next generation can be determined by 𝐗_𝑖^𝑘+1=𝐘^𝑘_𝑖,if𝐹_𝐶𝐗^𝑘_𝑖<𝐹_𝐶𝐘^𝑘_𝑖,𝐗^𝑘_𝑖,otherwise. (18) For the minimum problem, the current best individual 𝐩^𝑘_best at the 𝑘th generation is defined as [35] 𝐩^𝑘_best∶=arg_𝐗𝑘min𝐹_𝐶𝐗^𝑘_𝑖,∀𝑖; (19)𝐠^𝑘_best which denotes the best individual in the previous entire generations is defined as [35] 𝐠^𝑘_best=𝐩^𝑘_best,if𝐠^𝑘−1_best>𝐩^𝑘_best,𝐠^𝑘−1_best,otherwise. (20)

3.2. The Improved OSaRDE Algorithm

However, the optimum performances of different problems are not only highly dependent on the control parameters and learning strategies involved in DE, but also dependent on the characteristics of the region of the search landscape being explored [36, 37]. For a given task, we should spend a huge amount of time to try through various strategies and fine tune the corresponding parameters to achieve the best optimum performance. Motivated by this dilemma, large number of researchers have devoted themselves to develope adaptive DE algorithm to solve general problems more efficiently. Up to now, a classification, into three macrogroups of the main adaptive strategies employed in the DE scheme can be proposed here. (1) Self-adapting parameter setting of 𝐹 and CR, including the pseudorandom based updating approach [34, 36, 38] and the local optimum search based updating approach [39, 40]. (2) Adaptive selection of mutation strategies according to their capacities of generating promising solutions [36, 37, 41]. (3) Dynamic population sizing strategy based on self-adaptation [42] or progressive reduction [43, 44].

The classical DE algorithm typically converges at a rapid pace in the initial stages of the search procedure and then gradually slows down or keeps unchanged as it approaches global optimum. In case there is no improvement in the fitness in successive generations, it would be an indication that the individuals are clustered together in a local region and the operations are very insignificant to allow any improvement of the evolution process [45]. As is evident from this phenomenon, the information of the evolutionary process contained in these individuals would overlap and become extensively redundant. However, we do not suggest that we should discard the redundant individuals by the dynamic population sizing strategy. If the redundant individuals were simply discarded, the information contained in these individuals would be totally lost. Moreover, the population sizing strategy also has no assistance to the performance of escaping from local attractors.

In order to avoid the redundancy and enhance the population diversity, Ali and Pant (2011) take advantage of Cauchy mutation to perturb the clustered individuals and force them to move to someother location, thereby providing them an opportunity to improve their performance [45]. Being totally different from Cauchy mutation based and other adaptive strategies based DE algorithm mentioned above, a novel Opposition-based Self-adaptive Replacement Differential Evolution (OSaRDE) algorithm is proposed in this research to enhance the performance of DE. In the OSaRDE scheme, the redundant inferior individuals are self-adaptively selected according to a criterion defined by Euclidean distance after a specified number of generations, and then each of them is replaced by its opposite based on the opposition-based learning (OBL) strategy. After a period of regular evolutions, the process will be repeated. More specifically OSaRDE’s basic strategy considering the constraint task will be formalized as follows.

3.2.1. The OSaRDE’s Basic Strategy

For the replacement operation in the 𝑘^𝑅_𝑛 generation, the current worst individual 𝐩^𝑘𝑅𝑛_worst can be defined as 𝐩^𝑘𝑅𝑛_worst∶=arg_{𝐗𝑘𝑅𝑛}max𝐹_𝐶𝐗^𝑘𝑅𝑛_𝑖,∀𝑖. (21)

The redundant inferior individuals are defined as 𝐔^𝑘𝑅𝑛_𝑚=𝐗^𝑘𝑅𝑛_𝑖if𝐗^𝑘𝑅𝑛_𝑖−𝐩^𝑘𝑅𝑛_worst,𝐩≤Θ⋅DisDis=^𝑘𝑅𝑛_best−𝐩^𝑘𝑅𝑛_worst,𝑖=1,2,…,NP,𝑚=1,2,…,𝑀,𝑛=1,2,…,𝑁, (22) where 𝑀 denotes the number of inferior individuals, which is self-determined according to (22); Θ representing the replacement ratio is a user-defined constant; ‖⋅‖ denotes the operator for Euclidean distance; 𝑁 denotes the total number of the replacement operations during the evolution process. Let 𝑇, 𝑇_𝐶, and 𝑇_𝑃 represent the number of total generations, the specified number of generations from which the replacement operation begins, and the evolution period during two adjacent replacement operations. Considering that the replacement operation should not be executed at the last generation, 𝑁 can be obtained as 𝑁=ceil𝑇−𝑇_𝐶𝑇_𝑃−1. (23)

The opposite of each selected inferior individual can be defined as [26, 46] 𝐎𝐔^𝑘𝑅𝑛_𝑚=𝐌𝐈𝐍^𝐏+𝐌𝐀𝐗^𝐏−𝐔^𝑘𝑅𝑛_𝑚, (24) where 𝐌𝐈𝐍^𝐏=[min^𝑝₁,min^𝑝₂,…,min^𝑝_𝐷] and 𝐌𝐀𝐗^𝐏=[max^𝑝₁,max^𝑝₂,…,max^𝑝_𝐷] represent the predefined boundaries of the optimum task.

By applying the replacement operation, the population in the 𝑘th generation can be obtained as 𝐗_𝑖=𝐎𝐔^𝑘𝑅𝑛_𝑚if𝐗^𝑘𝑅𝑛_𝑖−𝐩^𝑘𝑅𝑛_worst.𝐗≤Θ⋅Dis_𝑖otherwise. (25)

By applying the replacement approach to the self-adaptively selected inferior individuals, the evolutionary process dynamically calculates the corresponding opposite according to the current location of each selected individual, thereby forcing them to jump to new solution candidates. The evolution process would possess much higher diversity and accordingly escape from the problem of premature convergence. Then the evolution process could have more opportunities to find better individuals during a period of regular evolutions after a replacement operation.

3.2.2. The Constraint Handling Technique

Generally, most of the optimum problems in engineering application are constrained. Unfortunately, the constraint of the searching space is not considered in the conventional DE scheme. A simple and popular repair operation of the individuals beyond the boundary works as follows [47]: if the 𝑗th element 𝑥^𝑘_𝑖,𝑗 of the mutant vector 𝐗^𝑘_𝑖 is out of the search region [min^𝑝_𝑗,max^𝑝_𝑗], then 𝑥^𝑘_𝑖,𝑗 is reset as 𝑥^𝑘_𝑖,𝑗=minmax^𝑝_𝑗,2min^𝑝_𝑗−𝑥^𝑘_𝑖,𝑗if𝑥^𝑘_𝑖,𝑗<min^𝑝_𝑗,maxmin^𝑝_𝑗,2max^𝑝_𝑗−𝑥^𝑘_𝑖,𝑗if𝑥^𝑘_𝑖,𝑗>max^𝑝_𝑗. (26)

For sake of clarity, the framework of the proposed OSaRDE is emphasized in Algorithm 1.

Parameters:   (i) NP: Population size;   (ii) F: Scaling factor;   (iii) CR: Crossover control parameter;   (iv) 𝑇𝐶: The specified number of generations from when the replacement operation begins;   (v) 𝑇𝑃: The period during two replacement operations;   (vi) G: The number of maximum generations. Step 1. Initialization:   S 1.1. Set the low and up boundaries of the searching space 𝐁∈[𝐌𝐈𝐍𝐏, 𝐌𝐀𝐗𝐏];   S 1.2. Randomly sample NP vectors 𝐗0𝑖=[𝑥0𝑖,1,𝑥0𝑖,2,𝑥0𝑖,3…,𝑥0𝑖,𝐷], 𝑖=1,2,…NP from B;   S 1.3. Evaluate the fitness values of these individuals. Step 2.  for 𝑘=1,2,…,𝐺, do:   S 2.1. Mutation operation:     Generate a mutant vector 𝐗𝑘𝑖=(𝑥𝑘𝑖,1,𝑥𝑘𝑖,2,𝑥𝑘𝑖,3…,𝑥𝑘𝑖,𝐷) via (16).   S 2.2. Crossover operation:     Execute the crossover operation to generate a trail vector 𝐘𝑘𝑖=(𝑦𝑘𝑖,1,𝑦𝑘𝑖,2,𝑦𝑘𝑖,3…,𝑦𝑘𝑖,𝐷) via (16).   S 2.3. Constraint repair:     if 𝐘𝑘𝑖 is not feasible (i.e. it is not in B), make it feasible via (26).   S 2.4. Selection operation:     S 2.4.1 if 𝐹𝐶(𝐘𝑘𝑖)<𝐹𝐶(𝐗𝑘𝑖), set 𝐗𝑘𝑖=𝐘𝑘𝑖, 𝐇𝐈(𝑖)=𝐹𝐶(𝐘𝑘𝑖);         else, set 𝐗𝑘𝑖=𝐗𝑘𝑖, 𝐇𝐈(𝑖)=𝐹𝐶(𝐗𝑘𝑖).     S 2.4.2 set 𝐩𝑘best∶=arg𝐗𝑘min{𝐇𝐈,∀𝑖}, 𝐩𝑘worst∶=arg𝐗𝑘max{𝐇𝐈,∀𝑖}, 𝐉best(𝑘)=min(𝐇𝐈).   S 2.5. Replacement operation:      if 𝑘≥𝑇𝐶 and 𝑘≠𝑇𝐶 and (𝐺−𝑘)/𝑇𝑃 == ceil((𝐺−𝑘)/𝑇𝑃), do       if 𝐗𝑘𝑖−𝐩𝑘worst≤Θ⋅Dis, Dis=𝐩𝑘best−𝐩𝑘worst;         𝐗𝑘𝑖=𝐌𝐈𝐍𝐏+𝐌𝐀𝐗𝐏−𝐗𝑘𝑖; else 𝐗𝑘𝑖=𝐗𝑘𝑖.   end for Step 3. Output 𝐩𝐺best and 𝐉best.

Algorithm 1 

The framework of the OSaRDE.

3.3. A Brief Testing on Benchmark Functions

A comprehensive set of benchmark functions, including eight typical global optimization problems, has been used only for performance verification and comparison of the proposed OSaRDE. The definition of the selected benchmark functions and their global optimum(s) are listed in the appendix [45, 48].

The eight algorithms described in the appendix have been tested and compared with standard DE and Cauchy mutation based DE (MDE) proposed by Ali and Pant (2011) [45]. Experiments related to DE and MDE in [45] have been employed as a reference. In order to have a fair comparison, the control parameters of OSaRDE are set as constant, where 𝐹=0.7 and CR=0.3, which are the same as employed in [45]. The number of function evaluations (NFEs) of OSaRDE are kept as 1 × 10⁵, which is half of that employed in [45]. The DE/cur-to-best/1 is chosen as the mutation strategy. In order to have statistically significant results, each test has been run 30 times. The average absolute error and standard deviation are recorded for comparisons. The corresponding results are given in Table 1, and the best results are highlighted in bold face. As evidence from the results in Table 1, it can be clearly observed that both DE and MDE algorithms give the same results for lower dimension problem, while almost all of the absolute errors obtained by OSaRDE are significantly smaller than that obtained by DE and MDE except for 𝑓_SWF(𝐗). As for 𝑓_SWF(𝐗), the average absolute error obtained by OSaRDE is larger than that obtained by MDE, while it is quite smaller than those obtained by DE. The results demonstrate that the proposed OSaRDE has a faster convergence mechanism without compromising with the quality of solution, significantly outperforming DE and MDE for both low and high dimension optimum tasks. It could be a fine choice for the OSaRDE to be applied to the engineering practice.

Table 1 

Comparisons of DE, MDE, and OSaRDE in terms of average absolute error and standard deviation for 30 runs.

4. Experiment Results and Discussion

Figure 3 illustrates the testing condition of the experiment part, where the LFTS is installed on a precision machine tool. The command signal to be applied to the LFTS is generated by a computer and converted through a data-acquisition card from ADLINK. The displacement of the LFTS is measured by a Renishaw linear encoder with a resolution of 10 nm. The measured displacement which is a digital signal is then directly gathered and stored in the computer through the data-acquisition card with a sample time of 0.45 ms for further analysis.

Figure 3 

The testing and working condition of the LFTS. ((1) The developed LFTS; (2) the workpiece; (3) the spindle; (4) the lathe bed.)

4.1. Identification Results of the LPM

To catch some prior knowledge of the LFTS system and accordingly help achieve more accurate parameters of the system, the sweep check is carried out to estimate the damped natural frequency. Thus, a sweep harmonic signal with constant amplitude and varying frequency ranged from 0 Hz to 400 Hz is applied to the VCM actuator. Fourier transform is then applied to the corresponding response, and the result is illustrated in Figure 4. As shown in Figure 4, the first damped natural frequency is about 172.6 Hz. Generally, the damping factor is about 0.05 and less than 0.1. With the roughly estimated first damped natural frequency and damping factor, the proportionality factor 𝐾 can be investigated as 5 × 10⁶. Thus, the feasible search space of the physical model can be defined as Γ_min=0,2(1−𝛼)𝜋×172.6,1−𝛽×5×10⁵,Γ_max=0.1,2(1+𝛼)𝜋×172.6,1+𝛽×5×10⁵, (27) where 𝛼=0.1 and 𝛽=0.5 are the scale factors of the first damped natural frequency and damping factor, respectively.

Figure 4 

Fourier transform of the response with sweep excitation.

To identify the parameters of the LPM, the step command with 1 V voltage and the corresponding response of the LFTS are employed for the parameter identification due to the reason that the step response contains abundant information of the system dynamics, such as the natural frequency and the damping behavior. To guarantee that the identified results have reached up to the maximum description capacity of the LPM, the optimum procedure is carried out by applying OSaRDE with different population sizes and generations, and the obtained best fitness values are further compared with them obtained by standard DE. The best fitness values are summarized in Table 2, and the convergence processes are illustrated in Figure 5. As given in Table 2, the proposed OSaRDE outperforms DE for low NFEs and gradually tends to be the same for large NFEs from the perspective of best fitness values. It also should be noticed that the best fitness values obtained by OSaRDE with different NFEs remain the same. It means that the OSaRDE is of high robustness for optimum tasks and could afford a much faster convergence rate than DE, thereby effectively reducing the computational time.

Table 2 

Fitness values with different population sizes and generations.

(a) NP = 30, 𝐺 = 50

(b) NP = 50, 𝐺 = 50

(c) NP = 50, 𝐺 = 100

(d) NP = 100, 𝐺 = 100

Figure 5 

The convergence of the optimum process (the red line and the blue line represent the process of OSaRDE and DE, resp.).

The transfer function of the LFTS can be obtained by introducing the optimum parameters corresponding to the best fitness value, that is 𝐺(𝑠)=6.19×10⁶𝑠²+50.97𝑠+1.196×10⁶. (28)

4.2. Comprehensive Results of the QPIM

In this experiment, an input excitation scheme with varying frequency and amplitude as shown in Figure 6(a) is applied to the VCM, and the corresponding output of the LFTS is illustrated in Figure 6(b). The 13600 sampling data are divided into two parts in which the former 6000 groups of data are used for training the NM and the latter 7600 groups of data are used for prediction. The modeling and prediction performance of the trained RBF NM is illustrated in Figure 7. As is evident from Figure 7, the NM could well approximate to the residual error of the LPM and present a good prediction capacity. To give a comprehensive comparison of the performances between the LPM and the proposed QPIM with neural model calibration, the modeling errors of the two models are illustrated in Figure 8. It can be seen that the maximum modeling error of the LPM is about ±1.15% of the full span range; meanwhile, obvious trend component attributing to the unmodeled parts of the physical process which could not be described by the LPM can be observed. The relatively small error verifies the efficiency of the proposed OSaRDE algorithm for system identification. As for the QPIM, the trend component in the residual modeling error of the LPM is well suppressed and the corresponding error maintains noise level. The result demonstrates that the QPIM is able to coincide more effectively with the response of the LFTS system; namely, the proposed modeling and identification method is of excellent performance.

(a) The excitation signal

(b) The response of the LFTS system

Figure 6 

The excitation signal and the corresponding response of the LFTS.

Figure 7 

Modeling and prediction results of the NM (the blue line denotes the modeling error of the LPM and the red line denotes the output of the NM).

Figure 8 

Modeling errors of the LPM and the QPIM (the blue line and the red line denote the modeling error of the LPM and the QPIM, resp.).

5. Conclusion

A quasiphysics intelligent model (QPIM) integrating a linear physics model (LPM) and a radial basis function neural network (RBFNN) based neural model (NM) is developed in this paper to accurately describe the dynamic behaviors of a voice coil motor (VCM) actuated long range fast tool servo (LFTS). A linear second-order differential equation of the LFTS system is established according to physics laws, and a novel Opposition-based Self-adaptive Replacement Differential Evolution (OSaRDE) algorithm is further developed for the parameter identification of the LPM. The NM is developed for the calibration of the unmodeled parts of the physics process attributing to the nonlinearity effect of the system. Both numerical and experiment examinations are carried out to respectively assess the performance of the OSaRDE algorithm and the QPIM. Certain conclusions can be drawn as follows.(1)A comprehensive set of benchmark functions, including eight typical global optimization problems, has been used for performance verification and comparison of the proposed OSaRDE. In comparison with both standard DE and a modified DE (MDE), it is observed that the OSaRDE performs remarkably better for both low and high dimension optimum tasks. The observation demonstrates that the proposed OSaRDE has a faster convergence mechanism without compromising with the quality of solution and outperforms DE and MDE taken for consideration. We may suggest that it should be a fine choice for the OSaRDE to be implemented for engineering practice.(2)As is evident from the identification results of the LPM and the QPIM, the modeling error of the LPM presents an obvious trend component which is about ±1.15% of the full span range. The relatively small error verifies the efficiency of the proposed OSaRDE algorithm for system identification. The residual modeling error is further reduced by NM calibration; the trend component in the residual error of the LPM is well suppressed, and the error of the QPIM maintains noise level. The result demonstrates that the QPIM is able to coincide more effectively with the response of the LFTS system, namely, the proposed modeling and identification method is of excellent performance.

Appendix

Benchmark Problems

(1)  Eason Function (EP):𝑓_EP(𝐗)=−cos𝑥₁cos𝑥₂exp−𝑥₁−𝜋²−𝑥₂−𝜋², (A.1)

with −10≤𝑥_𝑖≤10 and min𝑓_EP(𝜋,𝜋)=−1.

It is a multimodal, nonseparable function having several optima.

(2)  Sphere Function (SP):𝑓_SP(𝐗)=^𝑛_𝑖=1𝑥²_𝑖, (A.2)

with −100≤𝑥_𝑖≤100 and min𝑓_SP(0,…,0)=0.

It is a simple, continuous unimodal, separable, and highly convex function. It serves as test case for validating the convergence speed of an algorithm.

(3)  Ackley’s Function (ACK):𝑓_ACK(𝐗)=−20exp−0.2^𝑛_𝑖=1𝑥²_𝑖𝑛1−exp𝑛^𝑛_𝑖=1cos2𝜋𝑥_𝑖+20+𝑒, (A.3)

with −30≤𝑥_𝑖≤30 and min𝑓_ACK(0,…,0)=0.

It has numerous local minima, but its complexity is moderated. Any search strategies that analyzes a wider region will be able to cross the valley among the optima and achieve the global optimum.

(4)  Griewank Function (GW):𝑓_GW(𝐗)=14000^𝑛_𝑖=1𝑥²_𝑖−^𝑛_𝑖=1𝑥cos^𝑖√𝑖+1, (A.4)

with −600≤𝑥_𝑖≤600 and min𝑓_GW(0,…,0)=0.

It is a highly multimodal, nonseparable function. It has many regularly distributed local minima. It tests both convergence speed and the stability to escape from a shallow local minimum.

(5)  Levy and Montalvo 2 Problem (LM2):𝑓_LM2(𝐗)=sin²3𝜋𝑥₁+^𝑛−1_𝑖=1𝑥_𝑖−1²1+sin²3𝜋𝑥_𝑖+1+𝑥_𝑛−11+sin²2𝜋𝑥_𝑛, (A.5)

with −50≤𝑥_𝑖≤50 and min𝑓_GW(1,…,1)=0.

It is a multimodal, nonseparable function having several local optima.

(6)  Schwefel’s Problem (SWF):𝑓_SWF(𝐗)=418.9829𝑛−^𝑛_𝑖=1𝑥_𝑖sin𝑥_𝑖, (A.6)

with −500≤𝑥_𝑖≤500 and min𝑓_GW(420.97,…,420.97)=0.

It is a multimodal function with very deep sinusoidal interactions. It is generally considered to be difficult for optimization algorithms.

(7)  Rosenbrock Problem (RB):𝑓_RB(𝐗)=^𝑛−1_𝑖=1100𝑥_𝑖+1−𝑥²_𝑖²+1−𝑥_𝑖², (A.7)

with −30≤𝑥_𝑖≤30 and min𝑓_RB(1,…,1)=0.

As the number of dimensions increase, it ceases to be unimodal. The optimum of RB function lies inside a long, narrow, parabolic-shaped float valley. It helps in testing the ability of an algorithm to prevent premature convergence.

(8)  Rastrigin’s Function (RG):𝑓_RG(𝐗)=10𝑛+^𝑛_𝑖=1𝑥²_𝑖−10cos2𝜋𝑥_𝑖, (A.8)

with −5.12≤𝑥_𝑖≤5.12 and min𝑓_RG(0,…,0)=0.

It is a multimodal, nonseparable function having large numbers of local optima.

Acknowledgments

The authors are grateful to the financial support from the NSF of China (51175221; 51075041; 50995077), the Ministry of Science and Technology (MoSt) of China (2008AA04Z125), the Ministry of Education (MoE) of China (20070183104), and the Department of Science and Technology (DoST) of Jilin Province, China (20090337; 20100359; 20130522155JH).