# A Quasiphysics Intelligent Model for a Long Range Fast Tool Servo

Authors: Qiang Liu; Xiaoqin Zhou; Jieqiong Lin; Pengzi Xu; Zhiwei Zhu
Journal: The Scientific World Journal
 (2013)

---

## 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 novelOpposition-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.


---

## Body

## 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 novelOpposition-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 Figure1(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].Schematic of the LFTS mechanism.
(a)
LFTS mechanical structure
(b)
Dynamic model 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)Vm=Rmim+Lmdimdt+Vmb=Rmim+Lmdimdt+Kmvsdxdt,
where Vm is the voltage applied to the VCM, Vmb is the back electromotive force (EMF), and Kmvs is the equivalent back EMF constant. Rm, Lm, and im 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)F=nmBmlmim=nmBmlmRm(Vm-Lmdimdt-Kmvsdxdt),
where nm is the number of coil loops, Bm is the magnetic flux densities within the air gap between bobbin and magnets, and lm is the effective coil lengths considering all loops. Generally, the term Lmdim/dt can be considered negligible due to the fact that the electrical dynamics of the conducted current is much faster than mechanical behaviors; thus
(3)F=nmBmlmim=nmBmlmRm(Vm-Kmvsdxdt).Based on Newton’s second law of motion, the differential equation of dynamic motion for the LFTS mechanism can be given as follows:(4)(mfh+mVCM)d2xdt2+cfhdxdt+kfhx=F,
where mfh and mVCM represent the equivalent moving mass of the flexure-based mechanism and the bobbin of the VCM, respectively. cfh and kfh represent the equivalent damping coefficient and the equivalent stiffness of the flexure-based mechanism, respectively.Submitting (3) into (4) yields
(5)(mfh+mVCM)d2xdt2+(cfh+KmvsnmBmlmRm)dxdt+kfhx=nmBmlmRmVm.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:(6)P(s)=X(s)V(s)=Ks2+2ζωns+ωn2,
where K=nmBmlm/(mfh+mVCM)Rm, 2ζωn=(cfhRm+nmBmlmKmvs)/(mfh+mVCM)Rm, ωn2=kfh/(mfh+mVCM).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 densityBm [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 x*(t) be the real output of the LFTS; it would be difficult for the description of x*(t) 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)x*(t)=L-1(X(s))+NN(Vm,V˙m),
where L-1(·) represent the inverse Laplace transform operation and NN(·) 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]
(8)Fi=exp(-∥xi-ci∥22σi2),
where xi, ci, and σi 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]
(9)NN(Vm,V˙m)=∑i=1ncwi·Φi(Vm,V˙m),
where wi is the weight between the hidden layer and the output layer; nc 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 Figure2, where x(t) denotes the response of the LFTS system; x1(t) and x2(t) denote the response of the LPM and the NM, respectively. As discussed above, the LPM depends on a set of three parameters, that is, Γ=[ζ,ωn,K]∈ℜ3. Hence, the objective of the LPM identification is to find a set of parameters Γ*=[ζ*,ωn*,K*]∈ℜ3 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:
(10)E(Γ)=∑i=1N∥x1(ti)-x(ti)∥,
where i represents the index of the ith discrete sampling point.Figure 2
The principle of system identification.The criterionE(·) 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
(11)minimizeE(Γ)=∑i=1N∥x1(Γ,ti)-x(ti)∥,Γ=[ζ,ωn,K]s.t.Γ∈B,B={Γ,Γmin,j<Γj<Γmax,j,∀j=1,2,3Γ,Γmin,j<Γj<Γmax,j},
where B is the 3-dimensional feasible search space; Γmin,j and Γmax,j denote the upper bounds and the lower bounds of the jth parameter Γj, respectively.E
(
Γ
) 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(12)NE(w)=1N∑i=1N∥E(Γ*,ti)-x2(w,ti)∥.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].


## 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 Figure1(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].Schematic of the LFTS mechanism.
(a)
LFTS mechanical structure
(b)
Dynamic model 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)Vm=Rmim+Lmdimdt+Vmb=Rmim+Lmdimdt+Kmvsdxdt,
where Vm is the voltage applied to the VCM, Vmb is the back electromotive force (EMF), and Kmvs is the equivalent back EMF constant. Rm, Lm, and im 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)F=nmBmlmim=nmBmlmRm(Vm-Lmdimdt-Kmvsdxdt),
where nm is the number of coil loops, Bm is the magnetic flux densities within the air gap between bobbin and magnets, and lm is the effective coil lengths considering all loops. Generally, the term Lmdim/dt can be considered negligible due to the fact that the electrical dynamics of the conducted current is much faster than mechanical behaviors; thus
(3)F=nmBmlmim=nmBmlmRm(Vm-Kmvsdxdt).Based on Newton’s second law of motion, the differential equation of dynamic motion for the LFTS mechanism can be given as follows:(4)(mfh+mVCM)d2xdt2+cfhdxdt+kfhx=F,
where mfh and mVCM represent the equivalent moving mass of the flexure-based mechanism and the bobbin of the VCM, respectively. cfh and kfh represent the equivalent damping coefficient and the equivalent stiffness of the flexure-based mechanism, respectively.Submitting (3) into (4) yields
(5)(mfh+mVCM)d2xdt2+(cfh+KmvsnmBmlmRm)dxdt+kfhx=nmBmlmRmVm.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:(6)P(s)=X(s)V(s)=Ks2+2ζωns+ωn2,
where K=nmBmlm/(mfh+mVCM)Rm, 2ζωn=(cfhRm+nmBmlmKmvs)/(mfh+mVCM)Rm, ωn2=kfh/(mfh+mVCM).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 densityBm [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 x*(t) be the real output of the LFTS; it would be difficult for the description of x*(t) 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)x*(t)=L-1(X(s))+NN(Vm,V˙m),
where L-1(·) represent the inverse Laplace transform operation and NN(·) 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]
(8)Fi=exp(-∥xi-ci∥22σi2),
where xi, ci, and σi 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]
(9)NN(Vm,V˙m)=∑i=1ncwi·Φi(Vm,V˙m),
where wi is the weight between the hidden layer and the output layer; nc 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 Figure2, where x(t) denotes the response of the LFTS system; x1(t) and x2(t) denote the response of the LPM and the NM, respectively. As discussed above, the LPM depends on a set of three parameters, that is, Γ=[ζ,ωn,K]∈ℜ3. Hence, the objective of the LPM identification is to find a set of parameters Γ*=[ζ*,ωn*,K*]∈ℜ3 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:
(10)E(Γ)=∑i=1N∥x1(ti)-x(ti)∥,
where i represents the index of the ith discrete sampling point.Figure 2
The principle of system identification.The criterionE(·) 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
(11)minimizeE(Γ)=∑i=1N∥x1(Γ,ti)-x(ti)∥,Γ=[ζ,ωn,K]s.t.Γ∈B,B={Γ,Γmin,j<Γj<Γmax,j,∀j=1,2,3Γ,Γmin,j<Γj<Γmax,j},
where B is the 3-dimensional feasible search space; Γmin,j and Γmax,j denote the upper bounds and the lower bounds of the jth parameter Γj, respectively.E
(
Γ
) 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(12)NE(w)=1N∑i=1N∥E(Γ*,ti)-x2(w,ti)∥.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 D parameters, let B∈ℜD be the search space of the problem under consideration; the initial vector population randomly covering the entire parameter space can be expressed as
(13)Xi0=[xi,10,xi,20,xi,30,…,xi,D0]∈Bi=1,2,…,NP,
where NP denotes the population size;Xi0 denotes a potential solution for the task at start.(S2)  Mutation. For each target vector Xik at the kth generation, an associated mutant vector X^ik=(x^i,1k,x^i,2k,x^i,3k,…,x^i,Dk) 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]
(14)DE/rand/1:X^ik=Xr1k+F(Xr2k-Xr3k),(15)DE/best/1:X^ik=pbestk+F(Xr1k-Xr2k),(16)DE/cur-to-best/1:X^ik=Xik+F(pbestk-Xik)+F(Xr1k-Xr2k),
where r1, r2, and r3 represent the random and mutually different integers generated within the range [1,NP], and also different from index I; F is a mutation scale factor within the range [0,2], usually less than 1; pbestk denotes the best individual in generation k.(S3)  Crossover. After the mutation phase, the crossover operation is applied to each pair of the generated mutant vector X^ik and its corresponding target vector Xik to generate a trial vector Yik as follows:
(17)Yik=(yi,1k,yi,2k,yi,3k,…,yi,Dk),yi,jk={x^i,jk,if(rand(j)≤CRorj=randn(i)),xi,jk,if(rand(j)>CRorj≠randn(i)),
where j=1,2,…,D; i=1,2,…,NP; rand(j) is the jth independent random number uniformly distributed in the range of [0,1]. randn(i) is a randomly chosen index from the set {1,2,…,D}; 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 FC(·) denotes the objective function under consideration, the next generation can be determined by
(18)Xik+1={Yik,ifFC(Xik)<FC(Yik),Xik,otherwise.
For the minimum problem, the current best individual pbestk at the kth generation is defined as [35]
(19)pbestk:=argXkmin{FC(Xik),∀i};gbestk which denotes the best individual in the previous entire generations is defined as [35]
(20)gbestk={pbestk,ifgbestk-1>pbestk,gbestk-1,otherwise.


### 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 F 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 theknR generation, the current worst individual pworstknR can be defined as
(21)pworstknR:=argXknRmax{FC(XiknR),∀i}.The redundant inferior individuals are defined as(22)00000000UmknR=XiknRif(∥XiknR-pworstknR∥≤Θ·Dis),00000000000000000000000Dis=∥pbestknR-pworstknR∥,i=1,2,…,NP,m=1,2,…,M,n=1,2,…,N,
where M 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; N denotes the total number of the replacement operations during the evolution process. Let T, TC, and TP 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, N can be obtained as
(23)N=ceil(T-TCTP)-1.The opposite of each selected inferior individual can be defined as [26, 46]
(24)OUmknR=MINP+MAXP-UmknR,
where MINP=[min1p,min2p,…,minDp] and MAXP=[max1p,max2p,…,maxDp] represent the predefined boundaries of the optimum task.By applying the replacement operation, the population in thekth generation can be obtained as
(25)Xi={OUmknRif(∥XiknR-pworstknR∥≤Θ·Dis).Xiotherwise.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 jth element x^i,jk of the mutant vector X^ik is out of the search region [minjp,maxjp], then x^i,jk is reset as
(26)x^i,jk={min{maxjp,2minjp-x^i,jk}if(x^i,jk<minjp),max{minjp,2maxjp-x^i,jk}if(x^i,jk>maxjp).For sake of clarity, the framework of the proposed OSaRDE is emphasized in Algorithm1.Algorithm 1:The framework of the OSaRDE.
Parameters:
(i) NP: Population size;
(ii)F:  Scaling factor;
(iii) CR: Crossover control parameter;
(iv)TC: The specified number of generations from when the replacement operation begins;
(v)TP: 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 B∈[MINP, MAXP];
S 1.2. Randomly sample NP vectors Xi0=[xi,10,xi,20,xi,30…,xi,D0], i=1,2,…NP from B;
S 1.3. Evaluate the fitness values of these individuals.
Step
2.  fork=1,2,…,G, do:
S 2.1.Mutation operation:
Generate a mutant vectorX^ik=(x^i,1k,x^i,2k,x^i,3k…,x^i,Dk) via (16).
S 2.2.Crossover operation:
Execute the crossover operation to generate a trail vectorYik=(yi,1k,yi,2k,yi,3k…,yi,Dk) via (16).
S 2.3.Constraint repair:
ifYik is not feasible (i.e. it is not in B), make it feasible via (26).
S 2.4.Selection operation:
S 2.4.1ifFC(Yik)<FC(Xik), set Xik=Yik, HI(i)=FC(Yik);
else, set Xik=Xik, HI(i)=FC(Xik).
S 2.4.2set pbestk:=argXkmin{HI,∀i}, pworstk:=argXkmax{HI,∀i}, Jbest(k)=min(HI).
S 2.5.Replacement operation:
ifk≥TC and k≠TC and (G-k)/TP == ceil((G-k)/TP), do
if∥Xik-pworstk∥≤Θ·Dis, Dis=∥pbestk-pworstk∥;
Xik=MINP+MAXP-Xik; elseXik=Xik.
end for
Step
3.
Output pbestG and Jbest.


### 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 F=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 × 105, 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 fSWF(X). As for fSWF(X), 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.
Fun.
D
Absolute error
Standard deviation
DE
MDE
OSaRDE
DE
MDE
OSaRDE
f
EP
(
X
)
2
1.3424
e
-
06
1.3424
e
-
06
0.0000
e
-
00
6.9820
e
-
07
1.2595
e
-
06
0.0000
e
-
00
f
SP
(
X
)
30
8.2416
e
-
05
9.2528
e
-
07
6.0296
e
-
18
9.7617
e
-
06
1.9492
e
-
07
2.9392
e
-
17
f
ACK
(
X
)
30
2.5799
e
-
03
5.5200
e
-
04
2.5364
e
-
10
1.7069
e
-
4
8.7203
e
-
05
9.3680
e
-
10
f
GW
(
X
)
30
4.8967
e
-
04
1.8229
e
-
07
3.2746
e
-
09
9.3176
e
-
05
6.7985
e
-
08
1.1802
e
-
08
f
LM
2
(
X
)
30
9.3176
e
-
05
1.5778
e
-
07
4.8193
e
-
13
1.5190
e
-
06
204184
e
-
08
1.8087
e
-
12
f
SWF
(
X
)
30
5.9698
e
+
02
6.8199
e
-
04
2.9260
e
+
02
5.1068
e
+
02
6.8199
e
-
04
1.0148
e
+
01
f
RB
(
X
)
30
3.0242
e
+
01
2.4416
e
+
01
2.0256
e
+
01
9.7011
e
-
01
1.5638
e
-
01
8.1206
e
+
00
f
RG
(
X
)
30
1.3773
e
+
02
8.6468
e
+
01
8.2609
e
+
01
7.4039
e
-
00
1.2761
e
+
01
7.2384
e
+
00


## 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 D parameters, let B∈ℜD be the search space of the problem under consideration; the initial vector population randomly covering the entire parameter space can be expressed as
(13)Xi0=[xi,10,xi,20,xi,30,…,xi,D0]∈Bi=1,2,…,NP,
where NP denotes the population size;Xi0 denotes a potential solution for the task at start.(S2)  Mutation. For each target vector Xik at the kth generation, an associated mutant vector X^ik=(x^i,1k,x^i,2k,x^i,3k,…,x^i,Dk) 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]
(14)DE/rand/1:X^ik=Xr1k+F(Xr2k-Xr3k),(15)DE/best/1:X^ik=pbestk+F(Xr1k-Xr2k),(16)DE/cur-to-best/1:X^ik=Xik+F(pbestk-Xik)+F(Xr1k-Xr2k),
where r1, r2, and r3 represent the random and mutually different integers generated within the range [1,NP], and also different from index I; F is a mutation scale factor within the range [0,2], usually less than 1; pbestk denotes the best individual in generation k.(S3)  Crossover. After the mutation phase, the crossover operation is applied to each pair of the generated mutant vector X^ik and its corresponding target vector Xik to generate a trial vector Yik as follows:
(17)Yik=(yi,1k,yi,2k,yi,3k,…,yi,Dk),yi,jk={x^i,jk,if(rand(j)≤CRorj=randn(i)),xi,jk,if(rand(j)>CRorj≠randn(i)),
where j=1,2,…,D; i=1,2,…,NP; rand(j) is the jth independent random number uniformly distributed in the range of [0,1]. randn(i) is a randomly chosen index from the set {1,2,…,D}; 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 FC(·) denotes the objective function under consideration, the next generation can be determined by
(18)Xik+1={Yik,ifFC(Xik)<FC(Yik),Xik,otherwise.
For the minimum problem, the current best individual pbestk at the kth generation is defined as [35]
(19)pbestk:=argXkmin{FC(Xik),∀i};gbestk which denotes the best individual in the previous entire generations is defined as [35]
(20)gbestk={pbestk,ifgbestk-1>pbestk,gbestk-1,otherwise.


## 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 F 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 theknR generation, the current worst individual pworstknR can be defined as
(21)pworstknR:=argXknRmax{FC(XiknR),∀i}.The redundant inferior individuals are defined as(22)00000000UmknR=XiknRif(∥XiknR-pworstknR∥≤Θ·Dis),00000000000000000000000Dis=∥pbestknR-pworstknR∥,i=1,2,…,NP,m=1,2,…,M,n=1,2,…,N,
where M 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; N denotes the total number of the replacement operations during the evolution process. Let T, TC, and TP 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, N can be obtained as
(23)N=ceil(T-TCTP)-1.The opposite of each selected inferior individual can be defined as [26, 46]
(24)OUmknR=MINP+MAXP-UmknR,
where MINP=[min1p,min2p,…,minDp] and MAXP=[max1p,max2p,…,maxDp] represent the predefined boundaries of the optimum task.By applying the replacement operation, the population in thekth generation can be obtained as
(25)Xi={OUmknRif(∥XiknR-pworstknR∥≤Θ·Dis).Xiotherwise.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 jth element x^i,jk of the mutant vector X^ik is out of the search region [minjp,maxjp], then x^i,jk is reset as
(26)x^i,jk={min{maxjp,2minjp-x^i,jk}if(x^i,jk<minjp),max{minjp,2maxjp-x^i,jk}if(x^i,jk>maxjp).For sake of clarity, the framework of the proposed OSaRDE is emphasized in Algorithm1.Algorithm 1:The framework of the OSaRDE.
Parameters:
(i) NP: Population size;
(ii)F:  Scaling factor;
(iii) CR: Crossover control parameter;
(iv)TC: The specified number of generations from when the replacement operation begins;
(v)TP: 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 B∈[MINP, MAXP];
S 1.2. Randomly sample NP vectors Xi0=[xi,10,xi,20,xi,30…,xi,D0], i=1,2,…NP from B;
S 1.3. Evaluate the fitness values of these individuals.
Step
2.  fork=1,2,…,G, do:
S 2.1.Mutation operation:
Generate a mutant vectorX^ik=(x^i,1k,x^i,2k,x^i,3k…,x^i,Dk) via (16).
S 2.2.Crossover operation:
Execute the crossover operation to generate a trail vectorYik=(yi,1k,yi,2k,yi,3k…,yi,Dk) via (16).
S 2.3.Constraint repair:
ifYik is not feasible (i.e. it is not in B), make it feasible via (26).
S 2.4.Selection operation:
S 2.4.1ifFC(Yik)<FC(Xik), set Xik=Yik, HI(i)=FC(Yik);
else, set Xik=Xik, HI(i)=FC(Xik).
S 2.4.2set pbestk:=argXkmin{HI,∀i}, pworstk:=argXkmax{HI,∀i}, Jbest(k)=min(HI).
S 2.5.Replacement operation:
ifk≥TC and k≠TC and (G-k)/TP == ceil((G-k)/TP), do
if∥Xik-pworstk∥≤Θ·Dis, Dis=∥pbestk-pworstk∥;
Xik=MINP+MAXP-Xik; elseXik=Xik.
end for
Step
3.
Output pbestG and Jbest.


## 3.2.1. The OSaRDE’s Basic Strategy

For the replacement operation in theknR generation, the current worst individual pworstknR can be defined as
(21)pworstknR:=argXknRmax{FC(XiknR),∀i}.The redundant inferior individuals are defined as(22)00000000UmknR=XiknRif(∥XiknR-pworstknR∥≤Θ·Dis),00000000000000000000000Dis=∥pbestknR-pworstknR∥,i=1,2,…,NP,m=1,2,…,M,n=1,2,…,N,
where M 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; N denotes the total number of the replacement operations during the evolution process. Let T, TC, and TP 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, N can be obtained as
(23)N=ceil(T-TCTP)-1.The opposite of each selected inferior individual can be defined as [26, 46]
(24)OUmknR=MINP+MAXP-UmknR,
where MINP=[min1p,min2p,…,minDp] and MAXP=[max1p,max2p,…,maxDp] represent the predefined boundaries of the optimum task.By applying the replacement operation, the population in thekth generation can be obtained as
(25)Xi={OUmknRif(∥XiknR-pworstknR∥≤Θ·Dis).Xiotherwise.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 jth element x^i,jk of the mutant vector X^ik is out of the search region [minjp,maxjp], then x^i,jk is reset as
(26)x^i,jk={min{maxjp,2minjp-x^i,jk}if(x^i,jk<minjp),max{minjp,2maxjp-x^i,jk}if(x^i,jk>maxjp).For sake of clarity, the framework of the proposed OSaRDE is emphasized in Algorithm1.Algorithm 1:The framework of the OSaRDE.
Parameters:
(i) NP: Population size;
(ii)F:  Scaling factor;
(iii) CR: Crossover control parameter;
(iv)TC: The specified number of generations from when the replacement operation begins;
(v)TP: 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 B∈[MINP, MAXP];
S 1.2. Randomly sample NP vectors Xi0=[xi,10,xi,20,xi,30…,xi,D0], i=1,2,…NP from B;
S 1.3. Evaluate the fitness values of these individuals.
Step
2.  fork=1,2,…,G, do:
S 2.1.Mutation operation:
Generate a mutant vectorX^ik=(x^i,1k,x^i,2k,x^i,3k…,x^i,Dk) via (16).
S 2.2.Crossover operation:
Execute the crossover operation to generate a trail vectorYik=(yi,1k,yi,2k,yi,3k…,yi,Dk) via (16).
S 2.3.Constraint repair:
ifYik is not feasible (i.e. it is not in B), make it feasible via (26).
S 2.4.Selection operation:
S 2.4.1ifFC(Yik)<FC(Xik), set Xik=Yik, HI(i)=FC(Yik);
else, set Xik=Xik, HI(i)=FC(Xik).
S 2.4.2set pbestk:=argXkmin{HI,∀i}, pworstk:=argXkmax{HI,∀i}, Jbest(k)=min(HI).
S 2.5.Replacement operation:
ifk≥TC and k≠TC and (G-k)/TP == ceil((G-k)/TP), do
if∥Xik-pworstk∥≤Θ·Dis, Dis=∥pbestk-pworstk∥;
Xik=MINP+MAXP-Xik; elseXik=Xik.
end for
Step
3.
Output pbestG and Jbest.


## 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 F=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 × 105, 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 fSWF(X). As for fSWF(X), 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.
Fun.
D
Absolute error
Standard deviation
DE
MDE
OSaRDE
DE
MDE
OSaRDE
f
EP
(
X
)
2
1.3424
e
-
06
1.3424
e
-
06
0.0000
e
-
00
6.9820
e
-
07
1.2595
e
-
06
0.0000
e
-
00
f
SP
(
X
)
30
8.2416
e
-
05
9.2528
e
-
07
6.0296
e
-
18
9.7617
e
-
06
1.9492
e
-
07
2.9392
e
-
17
f
ACK
(
X
)
30
2.5799
e
-
03
5.5200
e
-
04
2.5364
e
-
10
1.7069
e
-
4
8.7203
e
-
05
9.3680
e
-
10
f
GW
(
X
)
30
4.8967
e
-
04
1.8229
e
-
07
3.2746
e
-
09
9.3176
e
-
05
6.7985
e
-
08
1.1802
e
-
08
f
LM
2
(
X
)
30
9.3176
e
-
05
1.5778
e
-
07
4.8193
e
-
13
1.5190
e
-
06
204184
e
-
08
1.8087
e
-
12
f
SWF
(
X
)
30
5.9698
e
+
02
6.8199
e
-
04
2.9260
e
+
02
5.1068
e
+
02
6.8199
e
-
04
1.0148
e
+
01
f
RB
(
X
)
30
3.0242
e
+
01
2.4416
e
+
01
2.0256
e
+
01
9.7011
e
-
01
1.5638
e
-
01
8.1206
e
+
00
f
RG
(
X
)
30
1.3773
e
+
02
8.6468
e
+
01
8.2609
e
+
01
7.4039
e
-
00
1.2761
e
+
01
7.2384
e
+
00


## 4. Experiment Results and Discussion

Figure3 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 Figure4. 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 K can be investigated as 5 × 106. Thus, the feasible search space of the physical model can be defined as
(27)Γmin=[0,2(1-α)π×172.6,(1-β)×5×105],Γmax=[0.1,2(1+α)π×172.6,(1+β)×5×105],
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 Table2, 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.
G
=
50
G
=
100
P
=
30
P
=
50
P
=
50
P
=
100
OSaRDE
1.517
e
+
04
1.517
e
+
04
1.517
e
+
04
1.517
e
+
04
DE
1.592
e
+
04
1.519
e
+
04
1.517
e
+
04
1.517
e
+
04The convergence of the optimum process (the red line and the blue line represent the process of OSaRDE and DE, resp.).
(a)
NP = 30,G = 50
(b)
NP = 50,G = 50
(c)
NP = 50,G = 100
(d)
NP = 100,G = 100The transfer function of the LFTS can be obtained by introducing the optimum parameters corresponding to the best fitness value, that is(28)G(s)=6.19×106s2+50.97s+1.196×106.


### 4.2. Comprehensive Results of the QPIM

In this experiment, an input excitation scheme with varying frequency and amplitude as shown in Figure6(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.The excitation signal and the corresponding response of the LFTS.
(a)
The excitation signal
(b)
The response of the LFTS systemFigure 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.).


## 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 Figure4. 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 K can be investigated as 5 × 106. Thus, the feasible search space of the physical model can be defined as
(27)Γmin=[0,2(1-α)π×172.6,(1-β)×5×105],Γmax=[0.1,2(1+α)π×172.6,(1+β)×5×105],
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 Table2, 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.
G
=
50
G
=
100
P
=
30
P
=
50
P
=
50
P
=
100
OSaRDE
1.517
e
+
04
1.517
e
+
04
1.517
e
+
04
1.517
e
+
04
DE
1.592
e
+
04
1.519
e
+
04
1.517
e
+
04
1.517
e
+
04The convergence of the optimum process (the red line and the blue line represent the process of OSaRDE and DE, resp.).
(a)
NP = 30,G = 50
(b)
NP = 50,G = 50
(c)
NP = 50,G = 100
(d)
NP = 100,G = 100The transfer function of the LFTS can be obtained by introducing the optimum parameters corresponding to the best fitness value, that is(28)G(s)=6.19×106s2+50.97s+1.196×106.


## 4.2. Comprehensive Results of the QPIM

In this experiment, an input excitation scheme with varying frequency and amplitude as shown in Figure6(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.The excitation signal and the corresponding response of the LFTS.
(a)
The excitation signal
(b)
The response of the LFTS systemFigure 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.



---

*Source: 641269-2013-09-18.xml*