# Improved Beetle Antennae Search Algorithm-Based Lévy Flight for Tuning of PID Controller in Force Control System

Authors: Yuqi Fan; Junpeng Shao; Guitao Sun; Xuan Shao
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

To enhance the anti-interference capability of an electrohydraulic force servo control system and increase the efficiency of the PID controller, this paper proposes a LBAS-PID controller. In LBAS, the random step created by the Lévy flight trajectory was used in the original algorithm to enhance the diversity of the population and convergence speed. In the force servo control system, LBAS-PID can enhance the performance of the system. First, the basic mathematical model of an electrohydraulic force servo control system was built based on theoretical analysis. The transfer function model was obtained by identifying the system parameters. Second, the introduced Lévy flight beetle antennae search algorithm was introduced and applied to ten benchmark functions, and the results were compared with those of other algorithms. Then, the proposed algorithm was applied in the PID controller to tune PID parameters in the force servo control system. To comprehensively evaluate performances of an electrohydraulic force servo control system that is controlled by the LBAS-PID controller, the frequency response analysis and temporal response analysis were obtained. The numerical analysis results indicate that an electrohydraulic force servo control system with an LBAS-PID controller could substantially increase the control characteristics of the system and restrain the external disturbances when different interference signals are examined.


---

## Body

## 1. Introduction

Electrohydraulic systems (EHSs) have been one of the most sought-after subjects in the industrial and academic fields because of high force/mass ratio, fast linear movement, and large torque/force output capability [1–3]. EHSs mainly consist of a servo valve, a hydraulic cylinder, and different kinds of detection sensors and other auxiliary parts [4, 5]. An electrohydraulic system (EHS), which includes a position servo control system and a force servo control system, is a typical closed-loop torque control system [6, 7], and EHSs also are multidomain energy system integrating characteristic of the mechanics, pattern recognition, digital technology, and control engineering [8, 9].In an electrohydraulic system, force generation and amplification are achieved through the force servo control system. Because of its superior load adaptability in the complex stress environments, the force servo control system plays a major role in engineering fields including electricity generator [10], flight simulator [11], gas turbine [12], robotic arm [13], and vehicles [14]. With the development of the electrohydraulic system, the force servo control system has been one of the hottest topics in engineering research. Due to changes in the loading environments, the interferences from unknown external loads, and the uncertain hydraulic parameters that are produced by variable load stiffness and oil temperature variations, the force servo control system is time-variant and nonlinear. To achieve the desired signal, not only synchronization precision of controllers but also tracking precision for the desired signal should be ensured. However, due to the eccentricity of mass load, cross-coupled characteristic, parameter uncertainty, and environmental interference, systems designed by simple controllers generally have the poor tracking and synchronization precision [15, 16].Among the common control methods in force, servo systems are the robust control method, adaptive control method, and PID control method. The robust control method is suitable for the application of stability as the primary objective and is most commonly used in overload environment [17]. Robust control systems do not have to adjust the parameters of the controller online, and robust control can guarantee good control performances when the dynamic characteristics of systems largely fluctuate [18]. The robust control method is based on some worst cases, so the controlled system does not work in the optimal state. When the prior knowledge of the controlled system is relatively small, the most commonly used control method is an adaptive control method. The adaptive control system can effectively find signal changes in real time by continuously measuring and monitoring controlled systems [19]. When the controlled system is a complex system, the adaptive control method has to identify numerous parameters by the long-running job, which can cause damage in transmission system parts and excessive bearing wear [20].The PID control system, which is also called the proportional-integral-derivative controller, is the favored control method [21]. Compared with two other control methods, the PID controller only has three parameters, which provides it with easy tuning space to get a better control point [22, 23]. It is shown that a rapid response speed and high control ability can be obtained by implementing appropriate PID parameters [24]. Therefore, the PID parameter tuning method has become one of the research hotspots in engineering fields. Currently, scientists have found that metaheuristic optimization algorithms are useful for solving professional engineering problems, and these algorithms have become very prevalent and have attracted the attention of many researchers. Several studies have been conducted that have enhanced the ability to obtain reasonable PID parameters, such as genetic algorithm (GA) [25, 26], particle swarm optimization algorithm (PSO) [27, 28], bat algorithm (BA) [29], water wave optimization (WWO) [30], artificial immune system algorithm (AIS) [31], grey wolf optimization algorithm (GWO) [32], cuckoo search optimization algorithm (CS) [33], and kidney-inspired algorithm (IKA) [34]. PID controller designed by other nature-inspired optimization algorithms is also applied in servo systems. Jaen-Cuellar et al. [35] introduced the GA-PID controller in servo systems of the industrial PUMA robot. Zhang et al. [36] used a differential evolution algorithm to design a PID controller in the pneumatic rotary actuator position servo system. Pršić et al. [37] applied the firefly algorithm to tune three PID parameters in the parallel robot platforms with six degrees of freedom gains. The literature [38] proposed a PID controller based on the self-growing Lévy flight salp swarm algorithm in hydraulic systems.The beetle antennae search algorithm (BAS) is a powerful and efficient metaheuristic optimization technique. This was proposed by Jiang and Li in 2017 [39]. The algorithm not only has individual and environment recognition abilities but also can ensure a compatible balance in the exploratory stages. Beetle antennae search algorithm has been used in many engineering fields, such as unmanned aerial vehicle avoidance [40] , path planning and control systems [41–43]. Although BAS has been widely used in engineering fields, it still has some shortcomings. Firstly, BAS has a uniformly reduced search step that makes the optimization solution jump from one region to another region, which can cause low search accuracy. Second, the initial fixed-step size search that is set according to artificial experiences can result in slow convergence speed and limited exploration efficiency in the later optimization period.To further increase the optimization capabilities of BAS while simultaneously strengthening the PID-working abilities of hydraulic force actuator systems, this paper proposes the Lévy flight beetle antennae search algorithm (LBAS) by Lévy random walking type whose step lengths are distributed according to heavy power-law tails, and then, this paper also proposes the LBAS-PID controller whose parameters are tuned by LBAS in force servo systems.The major contributions of this research work are summarized as follows:(1)
The Lévy flight beetle antennae search algorithm is proposed. In LBAS, the random step created by the Lévy flight trajectory will replace the fixed-step reduction trajectory of the original algorithm, and the random reduction factor will replace the reduction factor of the original algorithm.(2)
For the electrohydraulic system, a LBAS-PID controller is proposed and used in the force control system to enhance system performance. The parameter tuning problem in the PID controller will be converted into an algorithm optimization problem. Three parameters of the LBAS-PID controller will be tuned by LBAS.The remainder of the paper is arranged as follows.Section2 establishes the mathematical model of the electrohydraulic force servo control system. Section 3 introduces the original beetle antennae and the Lévy flight trajectory. A new algorithm which is called the Lévy flight beetle antennae search algorithm is presented. In Section 4, we not only present the LBAS-PID controller but also select the evaluation function ITAE. Section 5 tests ten benchmark functions by using different algorithms. Section 6 shows the results of ITAE value analysis and the temporal response analysis; the frequency response analysis of the electrohydraulic force servo control system added different PID controllers. Conclusions are given in the last section.


## 2. System Model

The basic structure of an electrohydraulic force servo control system includes a link mechanism, an electrohydraulic servo valve, a servo amplifier, a valve-controlled hydraulic cylinder, and a force sensor. The inner loop, which is driven by a hydraulic actuator, is a force closed-loop. The force sensor is modeled as an elastic internal link. The other parts are modeled as rigid links. The equivalent simplified physical model is presented in Figure1. The transfer function model of the electrohydraulic force servo system can be gained by using system parameter identification.Figure 1
The physical model of the electrohydraulic force servo control system.The working principle of the force servo control system is defined as follows. The servo amplifier works by transforming the voltage signal into an electric current signal which is input into the servo valve. The slide valve of the main valve will be shifted by depending on the electric current signal. The piston rod of the hydraulic cylinder will be drawn back or stretched by the movement of the slide valve. The force sensor which is linked with the piston rod works by measuring the load pressure. Finally, the signal will be fed back to the controller to drive the system because of the analog-digital conversion. Parameter identification is a technology that combines the mathematical model and experimental data for prediction. The technology can get system parameters according to the experimental data and the established model. When the error between results calculated by the model and the actual value is too large, the established model is amended, and the parameters are reselected. When the results are closer to the measurement results, the model is reliable. In hydraulic systems, some parameters, which cannot be got by numerical calculation, are time-varying. So the basic system model is applied in parameter identification.The schematic diagram of the valve-controlled hydraulic cylinder is shown in Figure2. The valve-controlled hydraulic cylinder can be simplified as a double-ended hydraulic actuator with the three-and-four-way servo valve.Figure 2
The schematic diagram of the valve-controlled hydraulic cylinder.Q
1 (m3/s) represents the oil ﬂow of the servo valve inlet. Q2 (m3/s) represents the oil ﬂow of the servo valve return. ps (Pa) and po (Pa) represent the inlet oil pressure and the return oil pressure. A (m2) is the piston area. p1 (Pa) and p2 (Pa) illustrate the working loads of the rodless cavity and the rod cavity, respectively. pL (Pa) means the pressure drop. xv (m) is the spool displacement of the main valve. y (m) is the displacement of the piston rod. m (kg) is the total equivalent mass referred to the piston. Bc (N/(m/s)) means the viscous damping coefficient. K (kN/m) denotes the spring stiffness. F (N) means an external load. y (m) is the displacement of the piston. The inlet oil pressure ps, which is the pressure threshold, is output by the pump. When Q1 and Q2 are influenced and throttled by the servo valve, the external load F is compensated by the load pressure pL of the hydraulic cylinder.To illustrate the model of the valve-controlled hydraulic cylinder, the miniaturized flow equation, hydraulic cylinder flow continuity equation, and the force-balance equation of the hydraulic cylinder need to be calculated. The linearize load flowQL for the servo valve can be expressed as(1)QL=Kqxv−KcpL,where Kq (m2/s) is the flow gain coefficient and Kc (m5/(N·s)) is the flow pressure coefficient, which can be written as(2)Kq=∂QL∂xv=Cdwps−pLρ,Kc=∂QL∂pl=Cdwxvps−pL/ρ2ps−pL,where Cd (m3/(s·pa)) is the discharge coefficient, wv (m) is the constant area gradient of valve orifices, and ρ (kg/m3) is the hydraulic oil density.The oil compressibility, volume changes in the roadless cavity and the rod cavity, flow equations of the oil ﬂow of servo valve inlet, and the oil ﬂow of servo valve return can be described as(3)Q1=Adydt+CippL+Cepp1+V1+Ayβedp1dt,Q2=Adydt+CippL−Cepp1−V2−Ayβedp1dt,where Cip (m3/(s·pa)) and Cep (m3/(s·pa)), respectively, represent the internal leakage coefficient and the external leakage coefficient, βe (N/(m2·pa)) is the oil effective bulk modulus, and V1 (m3) and V2 (m3) are the initial oil-in cavity volume and the initial oil-out cavity volume, respectively.As the piston rod of the actuator is located in the center of the hydraulic cylinder, linearize load flowQL is the average value of Q1 and Q2. The initial oil-in cavity and oil-out cavity have the same volume. The hydraulic cylinder flow continuity equation can be expressed as(4)QL=Adydt+CtppL+Vt4βedpLdt,where Ctp (m5/(N·s)) is the total leakage coefficient defined as Ctp = Cip + Cep/2 and Vt (m3) is the total volume of the cavity.The force-balance equation exerted on the payload of the hydraulic cylinder can be deduced by applying Newton’s second motion law in the noninertial reference frame and ignoring the friction of the cylinder and the mass of oil:(5)ApL=md2ydt2+Bcdydt+Ky.Equations (1), (4), and (5) are differential equations of dynamics in the valve-controlled hydraulic cylinder model, and their Laplace transform can be written as follows:(6)QL=KqXv−KcPL,QL=A1Ys+CtcPL+Vt4βePLs,APL=ms2Y+BcsY+KY.The electrohydraulic servo valve, which is used as an electrohydraulic transfer apparatus, works by transforming the lower input signal into the useful hydraulic pressure energy to accomplish targets of hydraulic systems. The transfer function of the servo valve can be expressed by(7)Gsvs=XvsIs=Ksv,where I (A) is the output current and Ksv (m/A) is the no-load flow gain.The servo amplifier works by transforming the input voltage and the input power, and it is composed of an electronic tube, the power transformer, and other electrical parts. The transfer function can be seen as the proportion of the current and the voltage. The transfer function can be expressed by(8)Ka=IsUs,where U (V) is the input voltage and Ka (A/V) is the amplification coefficient.The force sensor is a measuring device that converts the force signal into an electrical signal. It can measure mechanical characteristics, such as tension, pressure, torque, internal stress, and strain. The force sensor, which is used in the feedback stage, has the advantages of high receiving sensitivity, strong ability of anti-interference ability, and good performance in low-frequency detection. The response frequency of the force sensor is much larger than that of the system, which means that the transfer function can be seen as a pure proportional stage. The transfer function can be expressed by(9)Kf=UsFgs,where Fg (N) is the input force and Kf (V/N) is the force sensor coefﬁcient.The open-loop transfer function of the force servo control system from the servo valve spool displacement to the exerted force can be deduced from (1) to (9). The transfer function can be expressed by(10)Gforces=Kq/Ams2+Bcs+KmVt/4βeA2s3+mKce/A2+VtBc/4βeA2s2+1+BcKce/A2+KVt/4βeA2s+KKce/A2.


## 3. The Proposed Heuristic Optimization Algorithm

### 3.1. Beetle Antennae Search Algorithm

Beetles are the richest species of the order Coleoptera. The family contains approximately 35000 worldwide recorded species worldwide. Beetles have two long antennae that are often longer than the body. The two antennae can be used to enlarge the detecting area for detecting food resources and detecting the sex pheromone of a potentially suitable mate. In addition, antennae can also act as an exploration apparatus when exploring unknown areas. The exploration behaviors of beetles using two antennae can be modeled as the metaheuristic algorithm called the beetle antennae search algorithm (BAS). The position of each beetle, respectively, represents an achievable solution, so the optimum solution is regarded as the minimum distance from food. The beetle antennae search algorithm (BAS) can be optimized without the gradient information.The specific search process can be described as follows:Step 1. Define all beetle antennae search algorithm parameters. Randomly initialize N positions of beetles xn (n = 1, 2, …, N). Set the maximum number of iterations Kmax. Set k = 0.Step 2. In order to expand the initial exploration environment, initial antennae positions of beetles can be built to be normalized in randomly dimensional space. A random searching direction can be normalized as follows:(11)b⟶=rnddim,1rnddim,1,wherernd (.) denotes a random function whose dim represents dimensions of the solution.Step 3. Beetles use their antennae to determine the location of food when foraging. When the antenna on one side is closer to food, the food odor that is received by that antenna is stronger, and the individual will move to that same antenna side. The right antenna and the left antenna can be normalized as follows:(12)xrk=xk+dk⋅b⟶,(13)xlk=xk−dk⋅b⟶,wherek is the iteration number; xrk and xlk, respectively, denote positions of the right antenna and the left antenna; xk is the position of the beetle; and dk is the sensing length of the antennae.Step 4. The beetle determines the next position by detecting the odor; therefore, so we can explore the next location of one beetle according to the strength of the odor. They will move to the right if the right antenna receives a stronger scent than the left antenna; otherwise, it will move to the left. The beetle’s location is updated by(14)xk+1=xk+a⋅δk⋅b⟶⋅signfxrk−fxlk,whereδk is the step size of searching, f(.) is an evaluation function, and a is the movement direction of the beetle. Sign (.) represents a sign function:(15)signfxrk−fxlk=1,fxrk−fxlk>0,0,fxrk−fxlk=0,−1,fxrk−fxlk<0.Step 5. The sensing length of the antennae dk and the step size of searching δk can be updated as follows:(16)dk+1=w1∗dk+0.01,(17)δk+1=w2∗δk,wherew1 and w2 represent fixed reduction factors (typically between 0 and 1).Step 6. Calculate the evaluation function value of each individual using the evaluation function. Compare all feasible solutions to determine the optimal solution. Update the number of iterations k = k + 1 and return to Step 2. Repeat until k = Kmax.Step 7. Output the optimal solution.


### 3.2. Lévy Flight Trajectory

The Lévy flight trajectory, which was initially introduced by the French mathematician Lévy in 1925, is a particular type of stochastic random walk whose movements are based on a probability distribution. The Lévy flight trajectory defines a scale-invariance random walk model whose small steps connect longer relocations. The probability density distribution is as follows:(18)Ls∼s−β,where s is the Lévy flight length and β is the power-law exponent, 1 ≤ β ≤ 3.The Lévy flight walk, which is an optimal search theory, is a finite-velocity random walk whose step length is determined based on a fixed time, and it follows a dynamical motion process. Studies have proven that the flight and predation actions of insects obey the typical characteristics of the Lévy flight trajectory because insects walking paths in nature are considered to be as a random or quasi-random manner in nature. The Lévy flight trajectory is theorized to be the most efficient motion model for locating prey positions when the searching power-law exponentβ ∈ (1.5, 2]. The Lévy flight length s can be calculated by Mantegna’s algorithm, as follows:(19)s=uv1/β,where u and v obey normal distribution, u ∼ N (0, σu2), v ∼ N (0, σv2), σv = 1, and σu can be described as(20)σu=Γ1+β×sinπβ/2Γ1+β/2×β×2β−1/21/β.The two-dimensional Lévy flight trajectory whose steps were 1000 was implemented one independent run in MATLAB whenβ = 1.5. The result is shown in Figure 3. From Figure 3, we can see that there are numerous short-distance searching paths and occasional long-distance searching paths. This picture illustrates that the observed movement model across natural landscapes with different expected resource distributions conform to the algorithm optimization path. From the figure, we can see that the length and direction of the searching paths exhibit strong randomizations. The Lévy flight mechanism can enhance the local searching ability and avoid falling into a local minimum.Figure 3
The two-dimensional Lévy flight trajectory.


### 3.3. Lévy Flight Beetle Antennae Search Algorithm

The beetle search lengths and step sizes are fixed based on natural selection. Therefore, it is easy to move from one field to another, which favors the global search in the early optimization stage. However, the high jumping ability makes the local search not a sufficiently thorough optimization procedure, and information near the local optimum is insufficiently utilized. Therefore, the BAS has the drawbacks of a poor local search ability, premature convergence, and low optimization precision in later stages. To improve the exploitation ability and convergence velocity of the BAS, this paper proposes the Lévy flight beetle antennae search algorithm.Lévy flight can help the algorithm conduct efficient global exploration, avoid falling into the local minimum, and minimize the viciousness of the searching agents. To obtain a better optimization result and avoid the stochastic blindness of the heavy-tailed distribution, we introduce a self-learning searching method into the searching path. The self-learning searching method works by regulating the displacement difference between the optimal position and the individual position in each iteration. The searching scope is gradually changed from large to small by using the result of the last local optimization, and the repeating-cycle asymptotic searching is successively adopted. Therefore, the searching path in this paper consists of two parts: one part is the Lévy flight updating strategy, and the other part is the self-learning searching method. The Lévy flight updating strategy can increase the migration speed of other solutions approaching the optimum solution. The self-learning searching method can enhance the diversity of the predatory individual position, which makes every individual have the self-evolution capability.The new step size of searchingδ can be expressed as follows:(21)δnewk=s⋅xk−xbest,where s is the Lévy flight length and xbest is the current optimal solution.The new next exploration location of one beetle can be seen by(22)xnewk+1=xk+δnewk⋅a⋅b⟶⋅signfxrk−fxlk.The ideal optimization ability consists of a high searching ability in the early phase and a strong searching accuracy in the later phase. Thus, the new searching step size makes each individual have a high searching ability in the early phase and strong convergence accuracy in the later phase, which allows the individuals to be able to quickly obtain a better predation position and then conduct more careful and thorough searches around the better predation positions. Based on the above analysis, a random disturbance method that gradually changes the search scope from large to small is added to the sensing length of the antennae. The new sensing length of the antennae can be updated by(23)dnewk+1=dnewk∗w+0.01.The fix reduction factor is replaced by a random reduction factorw2 which will be chosen randomly in the range of [0, 1] in each iteration. To maximizing the searching capability, the initial value d0 is set to the upper bound of the search range. Random parameters can greatly improve search accuracy and efficiency. Initially, random parameters, which are created by the Lévy flight trajectory and random reduction factor, can increase the diversity of feasible solutions for continuous nonlinear optimization problems with higher dimensions within shorter computation times compared to the other evolutionary algorithms. In addition, the LBAS has fewer parameters than the BAS to adjust because the Lévy flight trajectory and random reduction factor will automatically select the random parameters in the LBAS. LBAS can be applied in various fields due to the high computational speed obtained by characteristics of random parameters.The LBAS main step can be summarized in the pseudocode shown as Algorithm1.Algorithm 1: The framework of LBAS.
Input: N positions of beetles xn (n = 1, 2, …, N). The search range [γmax, γmin]. Kmax·k = 0. Initial optimum solution xbest. Initial optimum value fbest. All initial parameters.
Output: xbest, fbest.
(1)
while (k < Kmax)
(2)
forn = 1 : N
(3)
Generate the searching direction by equation (11)
(4)
Calculate the right-left antenna position and by equations (12) and (13)
(5)
Calculate the Lévy flight lengths by equation (19)
(6)
Calculate the step size of searching by equation (21)
(7)
Calculate the next position of theith beetle by equation (22)
(8)
Calculate the function valuef(xnk)of ith beetle.
(9)
iff(xnk)is better than fbest
(10)
xbest=xnk
(11)
fbest = f(xnk)
(12)
end if
(13)
end for
(14)
Update the sensing length of the antennaedk by equation (23)
(15)
k = k + 1
(16)
end while


## 3.1. Beetle Antennae Search Algorithm

Beetles are the richest species of the order Coleoptera. The family contains approximately 35000 worldwide recorded species worldwide. Beetles have two long antennae that are often longer than the body. The two antennae can be used to enlarge the detecting area for detecting food resources and detecting the sex pheromone of a potentially suitable mate. In addition, antennae can also act as an exploration apparatus when exploring unknown areas. The exploration behaviors of beetles using two antennae can be modeled as the metaheuristic algorithm called the beetle antennae search algorithm (BAS). The position of each beetle, respectively, represents an achievable solution, so the optimum solution is regarded as the minimum distance from food. The beetle antennae search algorithm (BAS) can be optimized without the gradient information.The specific search process can be described as follows:Step 1. Define all beetle antennae search algorithm parameters. Randomly initialize N positions of beetles xn (n = 1, 2, …, N). Set the maximum number of iterations Kmax. Set k = 0.Step 2. In order to expand the initial exploration environment, initial antennae positions of beetles can be built to be normalized in randomly dimensional space. A random searching direction can be normalized as follows:(11)b⟶=rnddim,1rnddim,1,wherernd (.) denotes a random function whose dim represents dimensions of the solution.Step 3. Beetles use their antennae to determine the location of food when foraging. When the antenna on one side is closer to food, the food odor that is received by that antenna is stronger, and the individual will move to that same antenna side. The right antenna and the left antenna can be normalized as follows:(12)xrk=xk+dk⋅b⟶,(13)xlk=xk−dk⋅b⟶,wherek is the iteration number; xrk and xlk, respectively, denote positions of the right antenna and the left antenna; xk is the position of the beetle; and dk is the sensing length of the antennae.Step 4. The beetle determines the next position by detecting the odor; therefore, so we can explore the next location of one beetle according to the strength of the odor. They will move to the right if the right antenna receives a stronger scent than the left antenna; otherwise, it will move to the left. The beetle’s location is updated by(14)xk+1=xk+a⋅δk⋅b⟶⋅signfxrk−fxlk,whereδk is the step size of searching, f(.) is an evaluation function, and a is the movement direction of the beetle. Sign (.) represents a sign function:(15)signfxrk−fxlk=1,fxrk−fxlk>0,0,fxrk−fxlk=0,−1,fxrk−fxlk<0.Step 5. The sensing length of the antennae dk and the step size of searching δk can be updated as follows:(16)dk+1=w1∗dk+0.01,(17)δk+1=w2∗δk,wherew1 and w2 represent fixed reduction factors (typically between 0 and 1).Step 6. Calculate the evaluation function value of each individual using the evaluation function. Compare all feasible solutions to determine the optimal solution. Update the number of iterations k = k + 1 and return to Step 2. Repeat until k = Kmax.Step 7. Output the optimal solution.


## 3.2. Lévy Flight Trajectory

The Lévy flight trajectory, which was initially introduced by the French mathematician Lévy in 1925, is a particular type of stochastic random walk whose movements are based on a probability distribution. The Lévy flight trajectory defines a scale-invariance random walk model whose small steps connect longer relocations. The probability density distribution is as follows:(18)Ls∼s−β,where s is the Lévy flight length and β is the power-law exponent, 1 ≤ β ≤ 3.The Lévy flight walk, which is an optimal search theory, is a finite-velocity random walk whose step length is determined based on a fixed time, and it follows a dynamical motion process. Studies have proven that the flight and predation actions of insects obey the typical characteristics of the Lévy flight trajectory because insects walking paths in nature are considered to be as a random or quasi-random manner in nature. The Lévy flight trajectory is theorized to be the most efficient motion model for locating prey positions when the searching power-law exponentβ ∈ (1.5, 2]. The Lévy flight length s can be calculated by Mantegna’s algorithm, as follows:(19)s=uv1/β,where u and v obey normal distribution, u ∼ N (0, σu2), v ∼ N (0, σv2), σv = 1, and σu can be described as(20)σu=Γ1+β×sinπβ/2Γ1+β/2×β×2β−1/21/β.The two-dimensional Lévy flight trajectory whose steps were 1000 was implemented one independent run in MATLAB whenβ = 1.5. The result is shown in Figure 3. From Figure 3, we can see that there are numerous short-distance searching paths and occasional long-distance searching paths. This picture illustrates that the observed movement model across natural landscapes with different expected resource distributions conform to the algorithm optimization path. From the figure, we can see that the length and direction of the searching paths exhibit strong randomizations. The Lévy flight mechanism can enhance the local searching ability and avoid falling into a local minimum.Figure 3
The two-dimensional Lévy flight trajectory.


## 3.3. Lévy Flight Beetle Antennae Search Algorithm

The beetle search lengths and step sizes are fixed based on natural selection. Therefore, it is easy to move from one field to another, which favors the global search in the early optimization stage. However, the high jumping ability makes the local search not a sufficiently thorough optimization procedure, and information near the local optimum is insufficiently utilized. Therefore, the BAS has the drawbacks of a poor local search ability, premature convergence, and low optimization precision in later stages. To improve the exploitation ability and convergence velocity of the BAS, this paper proposes the Lévy flight beetle antennae search algorithm.Lévy flight can help the algorithm conduct efficient global exploration, avoid falling into the local minimum, and minimize the viciousness of the searching agents. To obtain a better optimization result and avoid the stochastic blindness of the heavy-tailed distribution, we introduce a self-learning searching method into the searching path. The self-learning searching method works by regulating the displacement difference between the optimal position and the individual position in each iteration. The searching scope is gradually changed from large to small by using the result of the last local optimization, and the repeating-cycle asymptotic searching is successively adopted. Therefore, the searching path in this paper consists of two parts: one part is the Lévy flight updating strategy, and the other part is the self-learning searching method. The Lévy flight updating strategy can increase the migration speed of other solutions approaching the optimum solution. The self-learning searching method can enhance the diversity of the predatory individual position, which makes every individual have the self-evolution capability.The new step size of searchingδ can be expressed as follows:(21)δnewk=s⋅xk−xbest,where s is the Lévy flight length and xbest is the current optimal solution.The new next exploration location of one beetle can be seen by(22)xnewk+1=xk+δnewk⋅a⋅b⟶⋅signfxrk−fxlk.The ideal optimization ability consists of a high searching ability in the early phase and a strong searching accuracy in the later phase. Thus, the new searching step size makes each individual have a high searching ability in the early phase and strong convergence accuracy in the later phase, which allows the individuals to be able to quickly obtain a better predation position and then conduct more careful and thorough searches around the better predation positions. Based on the above analysis, a random disturbance method that gradually changes the search scope from large to small is added to the sensing length of the antennae. The new sensing length of the antennae can be updated by(23)dnewk+1=dnewk∗w+0.01.The fix reduction factor is replaced by a random reduction factorw2 which will be chosen randomly in the range of [0, 1] in each iteration. To maximizing the searching capability, the initial value d0 is set to the upper bound of the search range. Random parameters can greatly improve search accuracy and efficiency. Initially, random parameters, which are created by the Lévy flight trajectory and random reduction factor, can increase the diversity of feasible solutions for continuous nonlinear optimization problems with higher dimensions within shorter computation times compared to the other evolutionary algorithms. In addition, the LBAS has fewer parameters than the BAS to adjust because the Lévy flight trajectory and random reduction factor will automatically select the random parameters in the LBAS. LBAS can be applied in various fields due to the high computational speed obtained by characteristics of random parameters.The LBAS main step can be summarized in the pseudocode shown as Algorithm1.Algorithm 1: The framework of LBAS.
Input: N positions of beetles xn (n = 1, 2, …, N). The search range [γmax, γmin]. Kmax·k = 0. Initial optimum solution xbest. Initial optimum value fbest. All initial parameters.
Output: xbest, fbest.
(1)
while (k < Kmax)
(2)
forn = 1 : N
(3)
Generate the searching direction by equation (11)
(4)
Calculate the right-left antenna position and by equations (12) and (13)
(5)
Calculate the Lévy flight lengths by equation (19)
(6)
Calculate the step size of searching by equation (21)
(7)
Calculate the next position of theith beetle by equation (22)
(8)
Calculate the function valuef(xnk)of ith beetle.
(9)
iff(xnk)is better than fbest
(10)
xbest=xnk
(11)
fbest = f(xnk)
(12)
end if
(13)
end for
(14)
Update the sensing length of the antennaedk by equation (23)
(15)
k = k + 1
(16)
end while


## 4. LBAS-PID Control System

### 4.1. PID Controller and Evaluation Function

PID controllers have a simple structure and high efficiency and thus are the most popular controller type in engineering. There are three parameters which include the proportional parameterKp, the differential parameter Kd, and the integral parameter Ki in the PID controller. Kp, Ki, and Kd, respectively, affect the response speed, the dynamic performance, and the steady-state error. It has been shown that a rapid response speed and high controllability can be obtained when three parameters are tuned fittingly. The current PID controller can be mainly divided into two modes: the continuous form and the discrete form. The continuous form can be expressed as follows:(24)ut=Kpet+Ki∫0tetdt+Kddetdt,where t is the time, Kp is the proportional constant, Ki is the integral constant, Kd is the derivative constant, u(t) is the signal of the PID controller, and e(t) is the error signal.The uninterrupted time is replaced by the sampling instant. The integral item and the differentiation item are discretized. The discrete form of the PID controller is written as follows:(25)uK=KpeK+KiT∑k=0meK+KdTeK−eK−1,where T, K, and m, respectively, represent the sampling period, the sampling number, and the total sampling numbers.The control signal is adjusted by using three parameters for the PID controller, and the signal is sent to the controller module. The PID parameter selection process can actually be converted into an optimization process. The PID parameters can be considered to be the position vector in the three-dimensional space. The result of the evaluation function is used to measure the performance of the control system. It is crucial to define the evaluation function before optimization. Evaluation functions of control systems can be mainly classified into two categories: the single objective function and the hybrid objective function. Single objective functions include the integration of the absolute value of error (IAE, integration absolute value error) and the integral of the squared value of error (ISE, integration square value error). Hybrid objective functions include the integral of time multiplied by the absolute value of error (ITAE, integration time absolute value error), and the mean of the square of the error (MSE, mean square error). The IAE, which uses importance to the absolute error, only takes the single factor into account. The single factor cannot thoroughly reflect the true state of the electrohydraulic system, and therefore, this mode cannot satisfy the multi-index analysis. The ISE, which is a single objective function, uses the squared of the error. Because the square of a large error is much larger than smaller errors, large errors are punished more. The accumulated small errors cause poor accuracy in the later stages of systems. The MSE reduces the weaknesses of the ISE by calculating the average of the ISE and time. However, to minimize the shortcomings of the ISE, the system has to run for a long time. The ITAE, which is a great measure of system performance, not only weighs long-term errors but also uses an additional time multiplication, which penalizes accumulated errors. Thus, ITAE is written as follows:(26)ITAE=∫0∞tetdt.


### 4.2. LBAS-PID Controller

The working principle of the PID controller that is tuned by LBAS-PID is that the PID controller drives the system by using three parameters that are found by the Lévy flight beetle antennae search algorithm. The ITAE value will be constantly calculated automatically by LBAS when the controlled system starts running. PID parameters, which are three quantitative values, includeki, kp, and kd. PID gain indicates the increased range of PID parameters when the system reaches critical stability. The minimum ITAE is used to update the optimum PID parameters in each iteration. To obtain optimal parameters when using the LBAS algorithm, the PID parameter optimization process can be modified into the three-dimensional optimization question in this paper. The three designed can be seen as coded solutions and beetle positions which are in a three-dimensional space, and each solution and position is represented by a real number. The ITAE is used as the evaluation function. The optimal system performance is regarded as the minimum of the objective function ITAE which is calculated by the LBAS. The beetle position, which minimizes the ITAE, is seen as the optimum PID parameters. The flow chart of LBAS-PID implementation for the force servo control system is shown in Figure 4. Tuning steps of the LBAS-PID controller are shown as follows:Step 1. Randomly generate N beetles Xnn¯=1,2,&,Nk=Kp,Ki,Kd, where each parameter uses the real number coding. Select discrete ITAE as the evaluation function. Determine the ranges of Kp, Ki, and Kd. Set the maximum iteration MAXiteration. Set the initial optimum solution xbest [0, 0, 0]. Set k = 0. Set the power-law exponent β. Set the search range [γmax,γmin]. Set d0 = γmax.Step 2. To make diversities of searching directions maximum, a searching direction in three-dimensional space can be calculated using equation (11).Step 3. Beetles use their right-left antennae to determine their next position when finding the minimum ITAE value. When the antenna on the one side is closer to the smaller ITAE, the beetle will move to this side. Positions of one beetle right-left antennae are calculated using equations (12) and (13).Step 4. Operate the control system. Calculate the displacement difference between the current optimal position which minimizes the ITAE and every beetle position. The step size of searching can be updated using equation (21).Step 5. Operate the control system. Right-left antennae positions of one beetle are transported into the control system to determine ITAE values. Calculate the ITAE values of right-left antennae positions. Update the next searching position using equation (22). Then, new beetle positions, which can be seen as new PID parameters, are carried over into the control system to calculate the ITAE value. Record all ITAE values.Step 6. By comparing all ITAE values, find the current minimum ITAE and the currently optimal beetle position. After comparing the current minimum ITAE with minimum ITAE of the last iteration, update the global minimum ITAE and global optimum beetle position. Then, the global optimum beetle position is seen as the global optimum PID parameters.Step 7. Randomly select the random reduction factor w in the range of [0, 1] in each iteration. Update the sensing length of the antennae using equation (23).Step 8. Calculate k = k + 1. Judge k = MAXiteration. If k is not equal to MAXiteration, return to Step 2. If k is equal to MAXiteration, stop.Figure 4
The flow chart of LBAS-PID implementation for the force servo control system.


## 4.1. PID Controller and Evaluation Function

PID controllers have a simple structure and high efficiency and thus are the most popular controller type in engineering. There are three parameters which include the proportional parameterKp, the differential parameter Kd, and the integral parameter Ki in the PID controller. Kp, Ki, and Kd, respectively, affect the response speed, the dynamic performance, and the steady-state error. It has been shown that a rapid response speed and high controllability can be obtained when three parameters are tuned fittingly. The current PID controller can be mainly divided into two modes: the continuous form and the discrete form. The continuous form can be expressed as follows:(24)ut=Kpet+Ki∫0tetdt+Kddetdt,where t is the time, Kp is the proportional constant, Ki is the integral constant, Kd is the derivative constant, u(t) is the signal of the PID controller, and e(t) is the error signal.The uninterrupted time is replaced by the sampling instant. The integral item and the differentiation item are discretized. The discrete form of the PID controller is written as follows:(25)uK=KpeK+KiT∑k=0meK+KdTeK−eK−1,where T, K, and m, respectively, represent the sampling period, the sampling number, and the total sampling numbers.The control signal is adjusted by using three parameters for the PID controller, and the signal is sent to the controller module. The PID parameter selection process can actually be converted into an optimization process. The PID parameters can be considered to be the position vector in the three-dimensional space. The result of the evaluation function is used to measure the performance of the control system. It is crucial to define the evaluation function before optimization. Evaluation functions of control systems can be mainly classified into two categories: the single objective function and the hybrid objective function. Single objective functions include the integration of the absolute value of error (IAE, integration absolute value error) and the integral of the squared value of error (ISE, integration square value error). Hybrid objective functions include the integral of time multiplied by the absolute value of error (ITAE, integration time absolute value error), and the mean of the square of the error (MSE, mean square error). The IAE, which uses importance to the absolute error, only takes the single factor into account. The single factor cannot thoroughly reflect the true state of the electrohydraulic system, and therefore, this mode cannot satisfy the multi-index analysis. The ISE, which is a single objective function, uses the squared of the error. Because the square of a large error is much larger than smaller errors, large errors are punished more. The accumulated small errors cause poor accuracy in the later stages of systems. The MSE reduces the weaknesses of the ISE by calculating the average of the ISE and time. However, to minimize the shortcomings of the ISE, the system has to run for a long time. The ITAE, which is a great measure of system performance, not only weighs long-term errors but also uses an additional time multiplication, which penalizes accumulated errors. Thus, ITAE is written as follows:(26)ITAE=∫0∞tetdt.


## 4.2. LBAS-PID Controller

The working principle of the PID controller that is tuned by LBAS-PID is that the PID controller drives the system by using three parameters that are found by the Lévy flight beetle antennae search algorithm. The ITAE value will be constantly calculated automatically by LBAS when the controlled system starts running. PID parameters, which are three quantitative values, includeki, kp, and kd. PID gain indicates the increased range of PID parameters when the system reaches critical stability. The minimum ITAE is used to update the optimum PID parameters in each iteration. To obtain optimal parameters when using the LBAS algorithm, the PID parameter optimization process can be modified into the three-dimensional optimization question in this paper. The three designed can be seen as coded solutions and beetle positions which are in a three-dimensional space, and each solution and position is represented by a real number. The ITAE is used as the evaluation function. The optimal system performance is regarded as the minimum of the objective function ITAE which is calculated by the LBAS. The beetle position, which minimizes the ITAE, is seen as the optimum PID parameters. The flow chart of LBAS-PID implementation for the force servo control system is shown in Figure 4. Tuning steps of the LBAS-PID controller are shown as follows:Step 1. Randomly generate N beetles Xnn¯=1,2,&,Nk=Kp,Ki,Kd, where each parameter uses the real number coding. Select discrete ITAE as the evaluation function. Determine the ranges of Kp, Ki, and Kd. Set the maximum iteration MAXiteration. Set the initial optimum solution xbest [0, 0, 0]. Set k = 0. Set the power-law exponent β. Set the search range [γmax,γmin]. Set d0 = γmax.Step 2. To make diversities of searching directions maximum, a searching direction in three-dimensional space can be calculated using equation (11).Step 3. Beetles use their right-left antennae to determine their next position when finding the minimum ITAE value. When the antenna on the one side is closer to the smaller ITAE, the beetle will move to this side. Positions of one beetle right-left antennae are calculated using equations (12) and (13).Step 4. Operate the control system. Calculate the displacement difference between the current optimal position which minimizes the ITAE and every beetle position. The step size of searching can be updated using equation (21).Step 5. Operate the control system. Right-left antennae positions of one beetle are transported into the control system to determine ITAE values. Calculate the ITAE values of right-left antennae positions. Update the next searching position using equation (22). Then, new beetle positions, which can be seen as new PID parameters, are carried over into the control system to calculate the ITAE value. Record all ITAE values.Step 6. By comparing all ITAE values, find the current minimum ITAE and the currently optimal beetle position. After comparing the current minimum ITAE with minimum ITAE of the last iteration, update the global minimum ITAE and global optimum beetle position. Then, the global optimum beetle position is seen as the global optimum PID parameters.Step 7. Randomly select the random reduction factor w in the range of [0, 1] in each iteration. Update the sensing length of the antennae using equation (23).Step 8. Calculate k = k + 1. Judge k = MAXiteration. If k is not equal to MAXiteration, return to Step 2. If k is equal to MAXiteration, stop.Figure 4
The flow chart of LBAS-PID implementation for the force servo control system.


## 5. Benchmark Function Analysis

### 5.1. Experimental Environment

To test the performance of the proposed LBAS algorithm from different perspectives, ten benchmark functions that have been widely used in engineering are discussed. The detailed information of these functions is summarized in Table1. In Table 1, dim represents the dimension of each function, the range represents the searching space of each function, and fmin represents the ideal minimum of each function.Table 1
Basic information on benchmark functions.
Name
Function
Dim
Range
f
min
Wayburn Seader 01
f
1
x
=
x
1
6
+
x
2
4
−
17
2
+
2
x
1
+
x
2
−
4
2
2
[−100, 100]
0
Three-hump camel
f
2
x
=
2
x
1
2
−
1.05
x
1
4
+
x
1
6
/
6
+
x
1
x
2
+
x
2
2
2
[−100, 100]
0
Egg cate
f
3
x
=
x
1
2
+
x
2
2
+
25
sin
2
x
1
+
sin
2
x
2
2
[−100, 100]
0
Booth
f
4
x
=
x
1
+
2
x
2
−
7
2
+
2
x
1
+
x
2
−
5
2
2
[−100, 100]
0
Beale
f
5
x
=
1.5
−
x
1
+
x
1
x
2
2
+
2.25
−
x
1
+
x
1
x
2
2
2
+
2.625
−
x
1
+
x
1
x
2
3
2
2
[−100, 100]
0
Sphere
f
6
x
=
∑
i
=
1
dim
x
i
2
20
[−2, 2]
0
Schwefel 2.23
f
7
x
=
∑
i
=
1
dim
x
i
10
20
[−2, 2]
0
G function
f
8
x
=
∏
i
=
1
dim
4
x
i
−
2
+
i
−
2
/
2
/
i
/
2
20
[−2, 2]
0
Dixon price
f
9
x
=
x
1
−
1
2
+
∑
i
=
2
dim
i
2
x
i
2
−
x
i
−
1
2
20
[−2, 2]
0
Csendes
f
10
x
=
∑
i
=
1
dim
x
i
6
2
+
sin
1
/
x
i
20
[−2, 2]
0These functions are divided into low-dimensional functions and high-dimensional functions. Low-dimensional functions (f1−f5) can investigate the convergence speed and abilities of different algorithms. The complexity of the solution spaces increases as the dimensionality increases because the high-dimensional function (f6−f10) has several local optima; therefore, high-dimensional functions solve nonlinear and complex problems.To verify the performance of LBAS, this paper also compares the proposed method with other popular artificial intelligence algorithms, including PSO [44], GA [45], BA [46], and beetle antennae search algorithm (BAS). PSO is one of the evolutionary algorithms which is inspired not by the evolutionary mechanism of natural selection, but rather by the social swarm behavior such as birds. Each particle in the particle swarm represents a possible solution. Through the simple behavior and the information interaction of one particle, the problem can be solved. GA is a computational model which imitated the natural selection mechanism and Darwinian biological evolution mechanism. Genetic algorithm starts from individuals which represent potential solutions, and a population is composed of a number of individuals encoded by the gene. The optimal individual in the last population can be decoded as the optimal solution. Bat algorithm is an iterative optimization technique, which imitates catching food behavior of bats in nature. A group of random solutions is initialized in BA, and then, the local new solution is generated by a random flight around the optimal solution, which strengthens the local searching ability. PSO, GA, and BA are classical optimization algorithms which are the research hotspots in engineering fields. They are all metaheuristics, and they all use an iterative searching mechanism to find the optimal solution. There are numerous articles and numerous studies about PID controllers based on GA [25, 26], PSO [27, 28], and BA [29]. So this paper chooses PSO, GA, and BA as the comparison group.The initial parameter of all algorithms selected optimal parameters by original algorithm literature. For the LBAS, parameters were fixed as the power-law exponentβ = 1.5. For the BAS, parameters were fixed as the initial sensing length of the antennae d0 = 4, the initial step size δ0 = 4, and w1 = w2 = 0.98. For the PSO, parameters were set as learning factors c1 = c2 = 1 and the inertia weight w = 1. For the BA, parameters were set as the loudness coefficient α = 0.9, the rate coefficient γ = 0.9, and the pulse frequency f select random number uniformly distributed over [0, 1]. For the GA, the selection probability was generated by the elitism roulette wheel selection mechanism; the crossover probability Pcross = 0.7, and the mutation probability Pm = 0.1. To make a fair comparison, each algorithm was independently run 10 times in MATLAB, the maximum number of iterations was 500, and the population size was 50.


### 5.2. Results of the Benchmark Function Experiments

Calculation results that contain best, worst, mean, and std are listed in Table2. In Table 2, best, worst, mean, and std represent the optimal function index, the worst function value, the average value obtained by 10 independent runs, and the standard deviation of values obtained by 10 independent runs. We can find that the values that are calculated by the LBAS are much closer to the ideal values. In addition, when the functions become more complex as the dimension increases, it is obvious to notice that the LBAS algorithm can obtain the best indices in Table 2. Testing can prove that the LBAS and BAS have the ability to provide a consistent and reliable optimal solution, but the solution quality of the LBAS is obviously better than that of the BAS in any environment. Therefore, it can be speculated that the LBAS algorithm has stronger computing power.Table 2
Comparison of results for various benchmark functions.
Function
Index
Algorithm
LBAS
BAS
PSO
BA
GA
f
1
Best
3.8639E − 16
1.0833E − 05
0.0063
6.3953E − 04
1.8311E − 04
Worst
3.5277E − 12
0.0830
0.8862
1.5578
0.3170
Mean
1.2562E − 12
0.0175
0.2484
0.2618
0.0616
Std
1.3435E − 12
0.0261
0.2769
0.4633
0.1026
f
2
Best
6.0983E − 41
2.0554E − 15
2.0579E − 04
2.2892E − 07
1.8008E − 06
Worst
7.6642E − 39
1.7528E − 14
0.0416
0.2987
3.4544E − 04
Mean
1.5018E − 39
5.6255E − 15
0.0139
0.0896
6.4699E − 05
Std
2.4509E − 39
4.5862E − 15
0.0141
0.1442
1.0509E − 04
f
3
Best
4.3498E − 39
1.2813E − 14
0.0112
2.4815E − 06
2.3647E − 06
Worst
2.2569E − 37
1.1819E − 12
0.7171
0.1037
0.0013
Mean
4.6175E − 38
4.8338E − 13
0.2507
0.0197
3.1942E − 04
Std
6.8052E − 38
5.9753E − 13
0.2136
0.0326
4.7989E − 04
f
4
Best
0
5.8376E − 29
6.9732E − 04
5.4196E − 06
5.0591E − 07
Worst
0
1.9822E − 26
0.3066
6.6743E − 04
8.6271E − 04
Mean
0
5.2070E − 27
0.0659
1.1809E − 04
2.0593E − 04
Std
0
6.6760E − 27
0.0921
2.0348E − 04
2.8323E − 04
f
5
Best
1.3012E − 19
8.6122E − 05
0.0058
0.0072
3.0196E − 04
Worst
1.2215E − 14
0.0052
1.0116
0.7818
0.1794
Mean
2.4014E − 15
0.0024
0.2478
0.3334
0.0328
Std
4.1688E − 15
0.0015
0.3026
0.2520
0.0537
f
6
Best
2.2719E − 10
0.0028
2.5178
0.6364
3.2472
Worst
5.5551E − 09
0.0445
9.1503
11.5886
5.2708
Mean
1.8718E − 09
0.0218
4.8465
3.4634
4.3997
Std
1.9900E − 09
0.0118
2.2263
3.5860
0.6088
f
7
Best
7.6355E − 11
0.0030
0.8001
0.0501
1.3554
Worst
4.0793E − 08
0.0621
23.4165
11.5856
5.6635
Mean
6.1733E − 09
0.0210
5.9295
3.9514
3.3672
Std
1.2642E − 08
0.0199
6.7150
4.3766
1.5632
f
8
Best
0
2.1247E − 07
0.0670
2.0867E − 04
0.1550
Worst
0
3.8616E − 05
2.2452
1.0323
0.5583
Mean
0
1.3238E − 05
0.7034
0.1666
0.3313
Std
0
1.3386E − 05
0.7085
0.3342
01341
f
9
Best
0.6674
0.8371
41.1291
3.9887
73.1997
Worst
0.7717
1.0727
209.7346
178.8130
104.0203
Mean
0.6830
0.9553
149.6278
58.2281
89.4864
Std
0.0320
0.0817
52.1859
57.3489
11.5272
f
10
Best
1.5734E − 11
1.6408E − 04
1.7869
0.4243
2.6959
Worst
5.0438E − 08
0.0179
12.7341
24.2207
10.2832
Mean
7.7986E − 09
0.0032
5.2805
9.2224
6.3464
Std
1.5499E − 08
0.0054
3.3581
8.7022
2.5984The average convergence curves of different algorithms when addressing all functions over 10 independent runs are shown in Figures5–14. It must be mentioned here that all the curves in Figures 5–14 are the average convergence curves. The LBAS has the fastest convergence speed in all algorithms, showing that it has the best performance when finding the global optimum. For the compared algorithms, the BA and PSO have poor performances for all functions, and they clearly cannot jump out of the local optimum in the last iteration stage. With respect to the five kinds of algorithm iteration curves, the LBAS achieves a distinguished convergence balance and the ability to jump out of the local optimum. It is certain that the LBAS has the fastest convergence rate in the early optimization stage and the highest searching accuracy in the later optimization stage. The underlying reason for the excellent performance of the proposed LBAS is that the step size and the orientation of the LBAS searching route are highly random based on the Lévy flight mechanism; therefore, the LBAS can easily to jump out of the local optimum when falling into local convergence, and it more quickly reaches the optimal value soon after the iterations start. Thus, we can conclude that the LBAS can absolutely obtain optimal solutions.Figure 5
The convergence curves forf1.Figure 6
The convergence curves forf2.Figure 7
The convergence curves forf3.Figure 8
The convergence curves forf4.Figure 9
The convergence curves forf5.Figure 10
The convergence curves forf6.Figure 11
The convergence curves forf7.Figure 12
The convergence curves forf8.Figure 13
The convergence curves forf9.Figure 14
The convergence curves forf10.The box plot chart, which is also known as a box-and-whiskers chart or a box-line chart, is a statistical chart that can represent a set of dispersed data. It can be used to compare several samples, and it can analyze the symmetry and the distribution of data. The box plot chart includes the maximum value, the minimum value, the median, the upper quartile, and the lower quantile. The stability and optimization ability of algorithms can be confirmed by evaluating the median, upper quartile, and lower quantile in the box plot chart. We can also deduce the performance of an algorithm using a box plot chart. The box plot charts of the values that were calculated by the 5 algorithms after 10 independent runs are shown in Figures15–24. It must be mentioned here that all values of the box plot charts are the values that are obtained from 10 independent runs. It is clear that the values that are calculated by the LBAS result in a narrowest box plot chart with the fewest outliers and the weakest dispersion, which proves that the LBAS has a stable convergence ability and the powerful optimization precision when dealing with different functions. The median, upper quartile, and lower quantile of the LBAS computing are lower than others, which also proves that the LBAS has an admirable optimization capability.Figure 15
The box plot charts forf1.Figure 16
The box plot charts forf2.Figure 17
The box plot charts forf3.Figure 18
The box plot charts forf4.Figure 19
The box plot charts forf5.Figure 20
The box plot charts forf6.Figure 21
The box plot charts forf7.Figure 22
The box plot charts forf8.Figure 23
The box plot charts forf9.Figure 24
The box plot charts forf10.


## 5.1. Experimental Environment

To test the performance of the proposed LBAS algorithm from different perspectives, ten benchmark functions that have been widely used in engineering are discussed. The detailed information of these functions is summarized in Table1. In Table 1, dim represents the dimension of each function, the range represents the searching space of each function, and fmin represents the ideal minimum of each function.Table 1
Basic information on benchmark functions.
Name
Function
Dim
Range
f
min
Wayburn Seader 01
f
1
x
=
x
1
6
+
x
2
4
−
17
2
+
2
x
1
+
x
2
−
4
2
2
[−100, 100]
0
Three-hump camel
f
2
x
=
2
x
1
2
−
1.05
x
1
4
+
x
1
6
/
6
+
x
1
x
2
+
x
2
2
2
[−100, 100]
0
Egg cate
f
3
x
=
x
1
2
+
x
2
2
+
25
sin
2
x
1
+
sin
2
x
2
2
[−100, 100]
0
Booth
f
4
x
=
x
1
+
2
x
2
−
7
2
+
2
x
1
+
x
2
−
5
2
2
[−100, 100]
0
Beale
f
5
x
=
1.5
−
x
1
+
x
1
x
2
2
+
2.25
−
x
1
+
x
1
x
2
2
2
+
2.625
−
x
1
+
x
1
x
2
3
2
2
[−100, 100]
0
Sphere
f
6
x
=
∑
i
=
1
dim
x
i
2
20
[−2, 2]
0
Schwefel 2.23
f
7
x
=
∑
i
=
1
dim
x
i
10
20
[−2, 2]
0
G function
f
8
x
=
∏
i
=
1
dim
4
x
i
−
2
+
i
−
2
/
2
/
i
/
2
20
[−2, 2]
0
Dixon price
f
9
x
=
x
1
−
1
2
+
∑
i
=
2
dim
i
2
x
i
2
−
x
i
−
1
2
20
[−2, 2]
0
Csendes
f
10
x
=
∑
i
=
1
dim
x
i
6
2
+
sin
1
/
x
i
20
[−2, 2]
0These functions are divided into low-dimensional functions and high-dimensional functions. Low-dimensional functions (f1−f5) can investigate the convergence speed and abilities of different algorithms. The complexity of the solution spaces increases as the dimensionality increases because the high-dimensional function (f6−f10) has several local optima; therefore, high-dimensional functions solve nonlinear and complex problems.To verify the performance of LBAS, this paper also compares the proposed method with other popular artificial intelligence algorithms, including PSO [44], GA [45], BA [46], and beetle antennae search algorithm (BAS). PSO is one of the evolutionary algorithms which is inspired not by the evolutionary mechanism of natural selection, but rather by the social swarm behavior such as birds. Each particle in the particle swarm represents a possible solution. Through the simple behavior and the information interaction of one particle, the problem can be solved. GA is a computational model which imitated the natural selection mechanism and Darwinian biological evolution mechanism. Genetic algorithm starts from individuals which represent potential solutions, and a population is composed of a number of individuals encoded by the gene. The optimal individual in the last population can be decoded as the optimal solution. Bat algorithm is an iterative optimization technique, which imitates catching food behavior of bats in nature. A group of random solutions is initialized in BA, and then, the local new solution is generated by a random flight around the optimal solution, which strengthens the local searching ability. PSO, GA, and BA are classical optimization algorithms which are the research hotspots in engineering fields. They are all metaheuristics, and they all use an iterative searching mechanism to find the optimal solution. There are numerous articles and numerous studies about PID controllers based on GA [25, 26], PSO [27, 28], and BA [29]. So this paper chooses PSO, GA, and BA as the comparison group.The initial parameter of all algorithms selected optimal parameters by original algorithm literature. For the LBAS, parameters were fixed as the power-law exponentβ = 1.5. For the BAS, parameters were fixed as the initial sensing length of the antennae d0 = 4, the initial step size δ0 = 4, and w1 = w2 = 0.98. For the PSO, parameters were set as learning factors c1 = c2 = 1 and the inertia weight w = 1. For the BA, parameters were set as the loudness coefficient α = 0.9, the rate coefficient γ = 0.9, and the pulse frequency f select random number uniformly distributed over [0, 1]. For the GA, the selection probability was generated by the elitism roulette wheel selection mechanism; the crossover probability Pcross = 0.7, and the mutation probability Pm = 0.1. To make a fair comparison, each algorithm was independently run 10 times in MATLAB, the maximum number of iterations was 500, and the population size was 50.


## 5.2. Results of the Benchmark Function Experiments

Calculation results that contain best, worst, mean, and std are listed in Table2. In Table 2, best, worst, mean, and std represent the optimal function index, the worst function value, the average value obtained by 10 independent runs, and the standard deviation of values obtained by 10 independent runs. We can find that the values that are calculated by the LBAS are much closer to the ideal values. In addition, when the functions become more complex as the dimension increases, it is obvious to notice that the LBAS algorithm can obtain the best indices in Table 2. Testing can prove that the LBAS and BAS have the ability to provide a consistent and reliable optimal solution, but the solution quality of the LBAS is obviously better than that of the BAS in any environment. Therefore, it can be speculated that the LBAS algorithm has stronger computing power.Table 2
Comparison of results for various benchmark functions.
Function
Index
Algorithm
LBAS
BAS
PSO
BA
GA
f
1
Best
3.8639E − 16
1.0833E − 05
0.0063
6.3953E − 04
1.8311E − 04
Worst
3.5277E − 12
0.0830
0.8862
1.5578
0.3170
Mean
1.2562E − 12
0.0175
0.2484
0.2618
0.0616
Std
1.3435E − 12
0.0261
0.2769
0.4633
0.1026
f
2
Best
6.0983E − 41
2.0554E − 15
2.0579E − 04
2.2892E − 07
1.8008E − 06
Worst
7.6642E − 39
1.7528E − 14
0.0416
0.2987
3.4544E − 04
Mean
1.5018E − 39
5.6255E − 15
0.0139
0.0896
6.4699E − 05
Std
2.4509E − 39
4.5862E − 15
0.0141
0.1442
1.0509E − 04
f
3
Best
4.3498E − 39
1.2813E − 14
0.0112
2.4815E − 06
2.3647E − 06
Worst
2.2569E − 37
1.1819E − 12
0.7171
0.1037
0.0013
Mean
4.6175E − 38
4.8338E − 13
0.2507
0.0197
3.1942E − 04
Std
6.8052E − 38
5.9753E − 13
0.2136
0.0326
4.7989E − 04
f
4
Best
0
5.8376E − 29
6.9732E − 04
5.4196E − 06
5.0591E − 07
Worst
0
1.9822E − 26
0.3066
6.6743E − 04
8.6271E − 04
Mean
0
5.2070E − 27
0.0659
1.1809E − 04
2.0593E − 04
Std
0
6.6760E − 27
0.0921
2.0348E − 04
2.8323E − 04
f
5
Best
1.3012E − 19
8.6122E − 05
0.0058
0.0072
3.0196E − 04
Worst
1.2215E − 14
0.0052
1.0116
0.7818
0.1794
Mean
2.4014E − 15
0.0024
0.2478
0.3334
0.0328
Std
4.1688E − 15
0.0015
0.3026
0.2520
0.0537
f
6
Best
2.2719E − 10
0.0028
2.5178
0.6364
3.2472
Worst
5.5551E − 09
0.0445
9.1503
11.5886
5.2708
Mean
1.8718E − 09
0.0218
4.8465
3.4634
4.3997
Std
1.9900E − 09
0.0118
2.2263
3.5860
0.6088
f
7
Best
7.6355E − 11
0.0030
0.8001
0.0501
1.3554
Worst
4.0793E − 08
0.0621
23.4165
11.5856
5.6635
Mean
6.1733E − 09
0.0210
5.9295
3.9514
3.3672
Std
1.2642E − 08
0.0199
6.7150
4.3766
1.5632
f
8
Best
0
2.1247E − 07
0.0670
2.0867E − 04
0.1550
Worst
0
3.8616E − 05
2.2452
1.0323
0.5583
Mean
0
1.3238E − 05
0.7034
0.1666
0.3313
Std
0
1.3386E − 05
0.7085
0.3342
01341
f
9
Best
0.6674
0.8371
41.1291
3.9887
73.1997
Worst
0.7717
1.0727
209.7346
178.8130
104.0203
Mean
0.6830
0.9553
149.6278
58.2281
89.4864
Std
0.0320
0.0817
52.1859
57.3489
11.5272
f
10
Best
1.5734E − 11
1.6408E − 04
1.7869
0.4243
2.6959
Worst
5.0438E − 08
0.0179
12.7341
24.2207
10.2832
Mean
7.7986E − 09
0.0032
5.2805
9.2224
6.3464
Std
1.5499E − 08
0.0054
3.3581
8.7022
2.5984The average convergence curves of different algorithms when addressing all functions over 10 independent runs are shown in Figures5–14. It must be mentioned here that all the curves in Figures 5–14 are the average convergence curves. The LBAS has the fastest convergence speed in all algorithms, showing that it has the best performance when finding the global optimum. For the compared algorithms, the BA and PSO have poor performances for all functions, and they clearly cannot jump out of the local optimum in the last iteration stage. With respect to the five kinds of algorithm iteration curves, the LBAS achieves a distinguished convergence balance and the ability to jump out of the local optimum. It is certain that the LBAS has the fastest convergence rate in the early optimization stage and the highest searching accuracy in the later optimization stage. The underlying reason for the excellent performance of the proposed LBAS is that the step size and the orientation of the LBAS searching route are highly random based on the Lévy flight mechanism; therefore, the LBAS can easily to jump out of the local optimum when falling into local convergence, and it more quickly reaches the optimal value soon after the iterations start. Thus, we can conclude that the LBAS can absolutely obtain optimal solutions.Figure 5
The convergence curves forf1.Figure 6
The convergence curves forf2.Figure 7
The convergence curves forf3.Figure 8
The convergence curves forf4.Figure 9
The convergence curves forf5.Figure 10
The convergence curves forf6.Figure 11
The convergence curves forf7.Figure 12
The convergence curves forf8.Figure 13
The convergence curves forf9.Figure 14
The convergence curves forf10.The box plot chart, which is also known as a box-and-whiskers chart or a box-line chart, is a statistical chart that can represent a set of dispersed data. It can be used to compare several samples, and it can analyze the symmetry and the distribution of data. The box plot chart includes the maximum value, the minimum value, the median, the upper quartile, and the lower quantile. The stability and optimization ability of algorithms can be confirmed by evaluating the median, upper quartile, and lower quantile in the box plot chart. We can also deduce the performance of an algorithm using a box plot chart. The box plot charts of the values that were calculated by the 5 algorithms after 10 independent runs are shown in Figures15–24. It must be mentioned here that all values of the box plot charts are the values that are obtained from 10 independent runs. It is clear that the values that are calculated by the LBAS result in a narrowest box plot chart with the fewest outliers and the weakest dispersion, which proves that the LBAS has a stable convergence ability and the powerful optimization precision when dealing with different functions. The median, upper quartile, and lower quantile of the LBAS computing are lower than others, which also proves that the LBAS has an admirable optimization capability.Figure 15
The box plot charts forf1.Figure 16
The box plot charts forf2.Figure 17
The box plot charts forf3.Figure 18
The box plot charts forf4.Figure 19
The box plot charts forf5.Figure 20
The box plot charts forf6.Figure 21
The box plot charts forf7.Figure 22
The box plot charts forf8.Figure 23
The box plot charts forf9.Figure 24
The box plot charts forf10.


## 6. Results and Discussion

### 6.1. Parameter Identification of the Force Servo Control System

The electrichydraulic servo valve is FF102-30 whose all parameters in the model are obtained from the servo valve product book, numerical derivation, and the experimental test. Inherent parameters were taken from the product book of the FF102-30, including effective the piston rod area, total stroke of the cylinder piston, and the total volume of the cavity. Working parameters were taken from the experimental and test numerical derivation, including no-load flow gain, the amplification coefficient, and the sensor coefﬁcient. Other parameters are determined by engineering practice experience. Piston diameterD = 50 mm. Piston rod diameter d = 35 mm. Inlet oil pressure Ps = 21 MPa. Rated flow Q = 30 (L/min). Amplification coefficient Ka = 0.016 (A/V). Force sensor coefﬁcient Kf = 0.0051 (V/N). No-load flow gain Ksv = 0.00579 (m/A). The platform of the force servo control system in Figure 25.Figure 25
The platform of the force servo control system.To illustrate the ability of the LBAS-PID controller, the PID controllers that use other algorithms were compared. PID controllers based on different tuning methods include GA-PID, PSO-PID, BA-PID, and BAS-PID. The parameters of all algorithms were the same as those in the test function environment. Three parameters of each PID controller were delivered by [0, 200]. In other words, initial population positions of all algorithms were generated in the range of [0, 200]. The total sampling number was 100, and the sampling period was 0.01s. Select discrete ITAE as the evaluation function. The disturbance signal was the step signal. The population size and total iterations were selected as 50 and 500 for all algorithms. Each optimization method was implemented for one independent run using the MATLAB software (MathWorks, Natick, MA, USA). The frequency response analysis and the temporal response analysis were all calculated by MATLAB.


### 6.2. Temporal Response Analysis

Optimal parameters of PID controllers and minimum ITAE values obtained with different algorithms can be seen in Table3. In Table 3, the LBAS-PID controller has the lowest ITAE. The simulation result evidence that the LBAS-PID controller has a more stable performance. When the input signal is the step signal, the time-varying process of the output signal is known as the temporal response of the system. The temporal response of the system is made up of the transient response and the steady-state response. The transient response, which reflects the system stability and speed, refers to the response process of the output signal from the initial state to the stable state. The transient response indices include overshoot Mp which reflects the balance capability and the settling time ts which reflects the ability to maintain stability ability. The overshoot indicated the difference value between the maximum value and the steady-state value. The settling time is the running time when the steady-state value of the response curve is within five percent of the ideal value. To determine the stability and the effectiveness of the temporal response of the proposed LBAS-PID controller, we compared the step response curves and the values of the temporal response indices for the different PID controllers when the input signal is the step signal. The results are also given in Table 3. As can be clearly seen from Table 3, temporal response indices using the LBAS-PID controller has the lowest overshoot Mp and the fastest settling time ts. The results can show that the system using the LBAS-PID controller can still achieve continue excellent dynamic characteristics and stability when subject to external interference signals.Table 3
PID parameters and performances.
Parameter
PID tuning methods
LBAS
BAS
PSO
BA
GA
K
p
1.7227
2.6289
1.6723
2.4663
0.9162
K
i
81.1255
90.0008
147.3600
59.9670
89.9934
K
d
0.0030
0.0051
0.0149
0.0001
0.0117
ITAE
7.498E − 04
9.077E − 04
0.001252
0.002353
0.001897
M
p
0
0.0604
0.1562
0
0.1122
t
s
0.06
0.11
0.15
0.21
0.16To determine the stability and the effectiveness of the temporal response of the proposed LBAS-PID controller, we compared the step response curves and the values of the temporal response indices and steady-state response indices for the different PID controllers when the input signal is the step signal. The comparison of the response curves that were derived from the step signal is shown in Figure26. Three-dimensional curves are shown in Figure 27. As seen from Figures 26 and 27, the LBAS-PID controller that is designed for an electrohydraulic servo system has the least overshoot, reflecting perfect robustness. The overshoots of the systems that were controlled by the GA-PID and the PSO-PID are too large. The response curve of BA-PID has zero overshoot, but it requires more time than the LBAS-PID to reach system stability. Though the response curves of the system that was controlled by the BAS-PID have less zero overshoots, there are still oscillations when the system starts to run.Figure 26
The step response curves.Figure 27
The three-dimensional step response curves.


### 6.3. Frequency Response Analysis

The frequency response is the response of the system under the sinusoidal signals. We use the amplitude-frequency characteristic which can reflect the resonance frequency, mechanical impedance, and dynamic stiffness in the system. The amplitude-frequency characteristic is defined as the amplitude ratio of the input signal to the ideal signal. To further prove the reliability of the LBAS-PID controller, the response results of the PID controllers based on the different algorithms are presented in Figures28–33 when the disturbance signals are six sinusoidal signals. For different sinusoidal signals, the angular velocity was, respectively, set as 2π, 4π, 5π, 6π, 8π, and 10π, the initial phase was zero, and the amplitude was set to 1, 10, and 20, respectively. Figures 28–33 show the response curves of different sinusoidal signals and the local amplification of those response curves. The actual amplitudes of the systems that are controlled by the GA-PID and PSO-PID are larger than the ideal amplitude in all figures, and the actual amplitudes of the systems that were controlled by other PID controllers are lower than the ideal amplitude. The system controlled by the GA-PID has the maximum amplitude overshoot. The system controlled by BA-PID has the maximum amplitude difference between the output amplitude and the ideal amplitude. Amplitude differences increase as angular velocity increase. Local enlarged drawings for the LBAS-PID controller are closest to the ideal amplitude, and the response curves for the LBAS-PID controller have the smallest distance between the ideal amplitude and the actual amplitude. Local amplification figures apparently display that the LBAS-PID controller has the performance of suppressing interference signals and reducing overshoot. For unusual signals combined with interference, the proposed PID controller can keep reliable inspection and eliminate the error. Therefore, we can deduce that the LBAS-PID controller has remarkable shock resistance and interference-resistant property under different signals.Figure 28
Response curves of the sinusoidal signal whose angular velocity is 2π: amplitude is (a) 1, (b) 10, and (c) 20.
(a)
(b)
(c)Figure 29
Response curves of the sinusoidal signal whose angular velocity is 4π: amplitude is (a) 1, (b) 10, and (c) 20.
(a)
(b)
(c)Figure 30
Response curves of the sinusoidal signal whose angular velocity is 5π: amplitude is (a) 1, (b) 10, and (c) 20.
(a)
(b)
(c)Figure 31
Response curves of the sinusoidal signal whose angular velocity is 6π: amplitude is (a) 1, (b) 10, and (c) 20.
(a)
(b)
(c)Figure 32
Response curves of the sinusoidal signal whose angular velocity is 8π: amplitude is (a) 1, (b) 10, and (c) 20.
(a)
(b)
(c)Figure 33
Response curves of the sinusoidal signal whose angular velocity is 10π: amplitude is (a) 1, (b) 10, and (c) 20.
(a)
(b)
(c)The Bode plot is composed of a magnitude characteristic plot and a phase characteristic plot. Bode plot is used to show the stability of a system in the frequency domain by observing its phase margin. The phase margin can be regarded as the phase changing before the system enters the unstable state. If the phase margin is greater than zero, the system is stable. The Bode plot for the proposed approach is shown in Figure34. In Figure 34, the blue dot represents minimum stability margins. In addition, the phase margin and delay margin obtained from Bode analysis are 85.9 and 0.0503, respectively. According to these results, the system is stable and has very good robustness.Figure 34
The Bode plot.


## 6.1. Parameter Identification of the Force Servo Control System

The electrichydraulic servo valve is FF102-30 whose all parameters in the model are obtained from the servo valve product book, numerical derivation, and the experimental test. Inherent parameters were taken from the product book of the FF102-30, including effective the piston rod area, total stroke of the cylinder piston, and the total volume of the cavity. Working parameters were taken from the experimental and test numerical derivation, including no-load flow gain, the amplification coefficient, and the sensor coefﬁcient. Other parameters are determined by engineering practice experience. Piston diameterD = 50 mm. Piston rod diameter d = 35 mm. Inlet oil pressure Ps = 21 MPa. Rated flow Q = 30 (L/min). Amplification coefficient Ka = 0.016 (A/V). Force sensor coefﬁcient Kf = 0.0051 (V/N). No-load flow gain Ksv = 0.00579 (m/A). The platform of the force servo control system in Figure 25.Figure 25
The platform of the force servo control system.To illustrate the ability of the LBAS-PID controller, the PID controllers that use other algorithms were compared. PID controllers based on different tuning methods include GA-PID, PSO-PID, BA-PID, and BAS-PID. The parameters of all algorithms were the same as those in the test function environment. Three parameters of each PID controller were delivered by [0, 200]. In other words, initial population positions of all algorithms were generated in the range of [0, 200]. The total sampling number was 100, and the sampling period was 0.01s. Select discrete ITAE as the evaluation function. The disturbance signal was the step signal. The population size and total iterations were selected as 50 and 500 for all algorithms. Each optimization method was implemented for one independent run using the MATLAB software (MathWorks, Natick, MA, USA). The frequency response analysis and the temporal response analysis were all calculated by MATLAB.


## 6.2. Temporal Response Analysis

Optimal parameters of PID controllers and minimum ITAE values obtained with different algorithms can be seen in Table3. In Table 3, the LBAS-PID controller has the lowest ITAE. The simulation result evidence that the LBAS-PID controller has a more stable performance. When the input signal is the step signal, the time-varying process of the output signal is known as the temporal response of the system. The temporal response of the system is made up of the transient response and the steady-state response. The transient response, which reflects the system stability and speed, refers to the response process of the output signal from the initial state to the stable state. The transient response indices include overshoot Mp which reflects the balance capability and the settling time ts which reflects the ability to maintain stability ability. The overshoot indicated the difference value between the maximum value and the steady-state value. The settling time is the running time when the steady-state value of the response curve is within five percent of the ideal value. To determine the stability and the effectiveness of the temporal response of the proposed LBAS-PID controller, we compared the step response curves and the values of the temporal response indices for the different PID controllers when the input signal is the step signal. The results are also given in Table 3. As can be clearly seen from Table 3, temporal response indices using the LBAS-PID controller has the lowest overshoot Mp and the fastest settling time ts. The results can show that the system using the LBAS-PID controller can still achieve continue excellent dynamic characteristics and stability when subject to external interference signals.Table 3
PID parameters and performances.
Parameter
PID tuning methods
LBAS
BAS
PSO
BA
GA
K
p
1.7227
2.6289
1.6723
2.4663
0.9162
K
i
81.1255
90.0008
147.3600
59.9670
89.9934
K
d
0.0030
0.0051
0.0149
0.0001
0.0117
ITAE
7.498E − 04
9.077E − 04
0.001252
0.002353
0.001897
M
p
0
0.0604
0.1562
0
0.1122
t
s
0.06
0.11
0.15
0.21
0.16To determine the stability and the effectiveness of the temporal response of the proposed LBAS-PID controller, we compared the step response curves and the values of the temporal response indices and steady-state response indices for the different PID controllers when the input signal is the step signal. The comparison of the response curves that were derived from the step signal is shown in Figure26. Three-dimensional curves are shown in Figure 27. As seen from Figures 26 and 27, the LBAS-PID controller that is designed for an electrohydraulic servo system has the least overshoot, reflecting perfect robustness. The overshoots of the systems that were controlled by the GA-PID and the PSO-PID are too large. The response curve of BA-PID has zero overshoot, but it requires more time than the LBAS-PID to reach system stability. Though the response curves of the system that was controlled by the BAS-PID have less zero overshoots, there are still oscillations when the system starts to run.Figure 26
The step response curves.Figure 27
The three-dimensional step response curves.


## 6.3. Frequency Response Analysis

The frequency response is the response of the system under the sinusoidal signals. We use the amplitude-frequency characteristic which can reflect the resonance frequency, mechanical impedance, and dynamic stiffness in the system. The amplitude-frequency characteristic is defined as the amplitude ratio of the input signal to the ideal signal. To further prove the reliability of the LBAS-PID controller, the response results of the PID controllers based on the different algorithms are presented in Figures28–33 when the disturbance signals are six sinusoidal signals. For different sinusoidal signals, the angular velocity was, respectively, set as 2π, 4π, 5π, 6π, 8π, and 10π, the initial phase was zero, and the amplitude was set to 1, 10, and 20, respectively. Figures 28–33 show the response curves of different sinusoidal signals and the local amplification of those response curves. The actual amplitudes of the systems that are controlled by the GA-PID and PSO-PID are larger than the ideal amplitude in all figures, and the actual amplitudes of the systems that were controlled by other PID controllers are lower than the ideal amplitude. The system controlled by the GA-PID has the maximum amplitude overshoot. The system controlled by BA-PID has the maximum amplitude difference between the output amplitude and the ideal amplitude. Amplitude differences increase as angular velocity increase. Local enlarged drawings for the LBAS-PID controller are closest to the ideal amplitude, and the response curves for the LBAS-PID controller have the smallest distance between the ideal amplitude and the actual amplitude. Local amplification figures apparently display that the LBAS-PID controller has the performance of suppressing interference signals and reducing overshoot. For unusual signals combined with interference, the proposed PID controller can keep reliable inspection and eliminate the error. Therefore, we can deduce that the LBAS-PID controller has remarkable shock resistance and interference-resistant property under different signals.Figure 28
Response curves of the sinusoidal signal whose angular velocity is 2π: amplitude is (a) 1, (b) 10, and (c) 20.
(a)
(b)
(c)Figure 29
Response curves of the sinusoidal signal whose angular velocity is 4π: amplitude is (a) 1, (b) 10, and (c) 20.
(a)
(b)
(c)Figure 30
Response curves of the sinusoidal signal whose angular velocity is 5π: amplitude is (a) 1, (b) 10, and (c) 20.
(a)
(b)
(c)Figure 31
Response curves of the sinusoidal signal whose angular velocity is 6π: amplitude is (a) 1, (b) 10, and (c) 20.
(a)
(b)
(c)Figure 32
Response curves of the sinusoidal signal whose angular velocity is 8π: amplitude is (a) 1, (b) 10, and (c) 20.
(a)
(b)
(c)Figure 33
Response curves of the sinusoidal signal whose angular velocity is 10π: amplitude is (a) 1, (b) 10, and (c) 20.
(a)
(b)
(c)The Bode plot is composed of a magnitude characteristic plot and a phase characteristic plot. Bode plot is used to show the stability of a system in the frequency domain by observing its phase margin. The phase margin can be regarded as the phase changing before the system enters the unstable state. If the phase margin is greater than zero, the system is stable. The Bode plot for the proposed approach is shown in Figure34. In Figure 34, the blue dot represents minimum stability margins. In addition, the phase margin and delay margin obtained from Bode analysis are 85.9 and 0.0503, respectively. According to these results, the system is stable and has very good robustness.Figure 34
The Bode plot.


## 7. Conclusions

To enhance the control performance and efficiency of the electrohydraulic force servo control systems, the LBAS-PID controller, whose PID parameters are tuned by the LBAS algorithm, has been incorporated into a force servo control system. The LBAS technique that is inspired by the Lévy flight trajectory, which is an advanced form of the basic beetle antennae search algorithm, has better convergence capabilities and can also compensate for the weak local search problem. Then, the basic mathematical model of an electrohydraulic force servo control system was systematically discussed using theoretical analysis and study. The transfer function model was calculated by system parameter identification. In order to demonstrate the effectiveness of the proposed LBAS-PID controller, the performance of the LBAS-PID controller was assessed using the ITAE and compared with that of BAS, GA, PSO, and BA. Finally, in this paper, the frequency response analysis and the temporal response analysis were carried on for electrohydraulic force servo control system in this paper. The comparative analysis results showed that the system with the LBAS-PID controller had the smallest temporal response indices and frequency response indices, which demonstrated that the system that was controlled by the LBAS-PID controller had the high robustness, the strong antiload disturbance ability, the quick reaction capability, and the strong parameter recognition characteristic when different jamming signals were received. Moreover, the performance of the LBAS computing ability was also tested using 10 benchmark functions and compared with the functional values that were calculated by the BAS, GA, PSO, and BA. The calculation results clearly proved that the LBAS had better identification performance and better global search competency than other metaheuristic algorithms. To obtain a more accurate solution when the LBAS is used to deal with complicated optimization problems, we would like to develop new metahybrid versions that fully use the advantages of other metaheuristic algorithms. In future work, we will study potential applications of the LBAS, including the optimal tuning of the fractional-order PID controllers in dynamic and heterogeneous environments. The FO-PID controller, which has two additional fractional-order operators, is a generalized expansion of the traditional PID controller. The FO-PID controller can be applied to increase the robustness to gain variations in plants and to obtain good output disturbance rejection.



---

*Source: 4287315-2020-03-24.xml*