A novel geo-inspired earthquake optimization algorithm

A novel geo-inspired earthquake optimization algorithm Efrain Mendez 1,†,‡, Alexandro Ortiz 1,‡∗, Pedro Ponce 2,‡ and Arturo Molina 1,†,‡ 1 Tecnologico de Monterrey, School of Engineering and Science, Mexico; a01336094@itesm.mx 2 Tecnologico de Monterrey, School of Engineering and Science, Mexico; a00787745@itesm.mx 3 Tecnologico de Monterrey, School of Engineering and Science, Mexico; pedro.ponce@itesm.mx 4 Tecnologico de Monterrey, School of Engineering and Science, Mexico; armolina@itesm.mx * Correspondence: a00787745@itesm.mx ‡ These authors contributed equally to this work. Academic Editor: name Version September 11, 2018 submitted to Abstract: A novel metaheuristic optimization method is proposed based on an earthquake that is 1 a geology phenomenon. The novel Earthquake Algorithm (EA) proposed, adapts the principle of 2 propagation of geology waves P and S through the earth material composed by random density 3 to ensure the dynamic balance between exploration and exploitation, in order to reach the best 4 solution to optimization complex problems by searching for the optimum into the search space. The 5 performance and validation of the EA are compared against the Bat Algorithm (BA) and the Particle 6 Swarm Optimization (PSO) by using 10 diverse benchmark functions. In addition, an experimental 7 engineering application is implemented to evaluate the proposed algorithm. Early results show a 8 feasibility of the proposed method with a clearly constancy and stability. It is important highlight the 9 fact that the main purpose of this paper is to present a new line of research, which is opened from the 10 novel EA. 11


Introduction
Optimization techniques are used to solve problems choosing the best option between all possible solutions, the metaheuristic algorithms are becoming powerful methods tough optimization problems [1].
According to [2], optimization algorithms can be divided into deterministic and stochastic classes.Deterministic methods are classified into evaluation and gradient-based methods.Evaluation methods, where calculating gradient functions is not necessary but they are extremely slow and ineffective.Gradient-based methods utilize gradients or derivations of the objective function to direct the exploration.Nevertheless, there is not security that the convergence to the global optimum by these methods is realized in objective functions with rough and complex forms.Otherwise, stochastic methods are not required gradient information to try to find the global optimum but they need to execute the objective function in order to find the optimum point [3].
In real and complex design problems, optimization algorithms may fail or result a local optima.Such problems cannot be resolved by conventional methods that locate local optima when the objective is to find the global optimum [4].
Metaheuristic algorithms offer an effective solution strategies due their distinct mechanisms to avoid from locally optimal solutions [5].The randomization procedure creates unsystematic solutions, which objective is to find a global solution through search space.By the local search, convergence and focus on realize good solutions in a region is determined [6].
In many instances, metaheuristics optimization algorithms have the same characteristics: (1) nature-inspired (with principles from physics, biology, geography, ethology), (2)  is used (including random variables), (3) gradient or Hessian matrix of the objective function is not used and (4) some parameters need to be adjusted on line with the problem to resolve [7].Another important characteristic is the strategy for exploring search spaces.A dynamic balance from diversification and intensification must be given.The first term allude to the examination of the search space, whereas the second term refers to exploitation according to search knowledge [8].
The last category is based from human behavior in different situations or activities such as: education behavior (Teaching-learning-based Optimization,TLBO) [33], neighborhood searching (Tabu-search) [34], socio-politically motivated strategy (Imperialist Competitive Algorithm, ICA) [35], musicians looking for the best melody (Harmony Search, HS) [36], mimic sport teams competition (League Championship Algorithm, LCA) [37].Fig. 1 depicts the broadly classification from metaheuristic optimization algorithms.Nevertheless, it is important highlight the fact that the main purpose of this paper is to present a new line of research, which is opened from the novel algorithm presented throughout this document.

Metaheuristic
In order to achieve the objective, Section 2 explains the basic concepts of an Earthquake behavior, then Section 3 shows the algorithm inspired by it.Later, in order to prove the viability of the proposed algorithm Section 4 shows its performance against some commonly used benchmark functions, and Section 5 its behavior against a real implementation (Speed controller optimization of a DC motor).
Finally, Section 6 opens the discussion of the obtained results and Section 7 explains the conclusions.

Earthquake background
Earthquakes are located in the "top list" of the most catastrophic and destructive natural disasters provoking fatalities in human life and economic lost [38].In seismic zones, elastic energy collects around those regions in advance of friction is outdo and unexpected shear movements occurs in form as an earthquake, decreasing surrounding strain energy and generating seismic waves that shake the external layer of the Earth [39].
From [40], seismological community has been discussing in the last decade that the first few seconds of the P-wave provides information about the final magnitude of an earthquake.This allows to use P-wave velocity to make a good exploration by delivering information from search space to find the optimum solution.
In [41], mentioned that seismic waves are elastic waves and each earth material must behave this characteristic to transmit them.Elasticity degree dictates the transmission capacity.By an earthquake, earth material is submitted to stress (compression, tension and/or shearing).Earth material has a ductile characteristic at very slow strain rates, as movements in the order of mm or cm/year.
Conversely, the earth acts elastically to fast and small deformations caused by earthquakes.When high amplitude of seismic deformations and long-period free-oscillations modes are present, the inelastic response of seismic waves must be considered.
According to [42], earth material behavior can be defined by Hooke's Law [43] in elastic range.
Above its elastic limit, earth material can be respond with brittle fracturing (e.g.fault by earthquake) or ductility.
Depending on the type of deformation, elastic material endures to stress and can be quantified by some elastic moduli [41]: • bulk modulus κ is described as the ratio of the hydrostatic pressure change to the resulting relative volume change (i.e.κ = ∆P/(∆V/V).
• Young's modulus E is defined by the response of a cylinder with Length which is pulled on both ends.Its value is given by E = (F/A)/(∆L/L).
• Poisson's ratio σ is the ratio between cylinder lateral contraction being pulled on its ends to its relative longitudinal extension, i.e. σ = (∆W/W)/(∆L/L).
• Lamé parameter λ does not have explanation physical but it can be defined in terms of the elastic moduli already mentioned λ = σE/((1 + σ)(1 − 2σ)) • density ρ of earth material.Denser rocks have faster wave propagation because rigidity increases with density.

P and S waves velocities
When a earthquake occurs, rocks break and generate waves through interior to surface of the Earth, Fig. 2.There are two types of waves P and S that are generated.The P-wave is the fastest and depends on the earth material compressibility.They are transmitted by compression and tension of the medium with volume changes.The S-wave is slower than P-wave, depends on rocks elasticity and move the epicenters up and down, perpendicular to the wave direction.
Eqs. ( 1) and (2) describe the estimation of P and S waves velocities.
where v p and v s are the P and S waves speed, λ and µ the Lamé parameters, and ρ the density of earth material.Records from ground-acceleration (Fig. 3) shows the waves P and S when an earthquake occurred in Chiapas, Mexico (09.15.2010) with an epicenter in 15.59 N, 93.52 W at 95 km of profundity [44].
P-waves are the first to appear (due to faster velocity) and are detected by seismograph accelerometers.
Seconds after, S-waves are observed.In the Fig. 4 can be visualized the movement of the P-waves through a medium by compression and dilation with volume changes.These kind of waves can be transmitted by any medium (solid, liquid and gas).Particle motion of earth material are represented by spheres.
Representation of the S-waves movement is shown in the Fig. 5, the propagation only occurs in solid medium with shearing deformation (perpendicular movements to the wave direction).

Earthquake optimization algorithm
The Earthquake algorithm initially works with a random population of solutions, which are called epicenters.Besides, it also uses the implementation of a function for the random motion in order to produce some members of the epicenters in distant regions of the searching space.The first version of this Earthquake algorithm was proposed in [45].Then, as will be explained later in this section, some of the parameters used for this novel algorithm can be described using a Poisson ratio, reason that inspired the implementation of an Exponential Distribution, from the relation between a Poisson and an Exponential distribution taken from [46], for the random generation already mentioned.
As seen in Section 2, the motion of an earthquake can be parametrized with the velocity equations of the P and S waves, where the principle of operation of the proposed algorithm lies in those parameters.Knowing then, that the transmission of the P-wave is faster than that of the S-wave, the P-wave is used for a quick exploration and the other one for a more detailed one.In order to determine when to use a wave or the other, it is essential for the algorithm to define an operation range for the S-wave, which will be referred to in this document as the S-range or Sr.Also, as seen in Fig. 7, the Sr is defined around the best solution.The Sr, should be assigned with the previous knowledge of the problem needs, nevertheless it is recommended to implement the range in function of the percentage of error, between the obtained solutions and the expected ones.Thus, in Fig. 6 can be seen a population of three epicenters, where the black epicenters are out of Sr, reason why they use the velocity of P-waves (black dotted lines), to map where the solution is (red triangle).The blue epicenter is located in smaller distance from the objective, it is also in Sr and that is why it uses the S-wave (blue dotted lines) for the mapping.
Knowing that the velocity of the P and S waves are given by Eqs. ( 1) and ( 2) respectively, it is also known that the Lamé parameters (λ and µ) and the density (ρ), are needed in order to determine the current transmission speed.
To contextualize property values analyzed in this work, Table 1 shows the principle values of density, Poisson's ratio and ratio of seismic wave velocities from some earth materials (average and/or ranges).Now, in order to validate the proposed constant for the Lamé parameters, [41] explains that most of the rocks have a Poisson ratio between 0.2 and 0.3, leaving a mean optimal value on 0.25 for the ratio.
As well, also taken from [41], the following equation shows the relation between the Poisson ratio and the Lamé parameters: Randomly select positive or negative v i . 19: if (rand > v p ) then 21: end if

23:
Constrain x i if needed.

24:
Obtain the new fitness values.
where σ is the Poisson's ratio.
According to [41], the Lamé parameters can be the same under some circumstances, so for the current algorithm it is taken that λ = µ.In that case, in order to find the optimal Lamé parameters to be used, several tests were performed with different Lamé values, finding that the only real constant that worked was 1.5; that is the reason why, for the algorithm purpose, the wave transmission parameters can be defined as the constant 1.5, taken from Table 1, giving to: Where substituting Eq. ( 4) on Eq. ( 3): Which validates the premise developed to obtain Eq. ( 4).On the other hand, for the implementation of the algorithm, the density of the solids (ρ) is used as a random value, selected from a range between 2200 and 3300 Kg/m 3 , also according to Table 1.
Shown the initial position of the epicenters population are randomly selected, however their speed are initially 0. On the other hand, the update of the current position of the epicenters, is given by Eq. ( 6): where X t i and X t−1 i are the current and the previous positions, meanwhile the V i is the current speed.
Finally, the equation that incorporates the Exponential distribution, to reduce the probability of visiting points already visited for the epicenters, or epicenters is described in [46].Where the random value is generated with the distribution, in a range of ± the maximum value of V P /V S , that is ± where X best is the global best solution, and the Exp µ (s) is the random value generated with the exponential distribution from the value of µ.
Whence, the random update for the position with the Exponential distribution for this algorithm, is given by Eq. ( 7), and the diagram that describes the proposed architecture for the earthquake algorithm is shown by Fig. 8.
Hence, the Sr is recommended to be ±10% from the best solution, but after selecting if the current epicenter is going to use the v s or v p as v i (see Fig. 8), it results very important for the algorithms performance to understand that since both speeds are calculated by a square root, the final result is a positive number, but it is also known that the result contemplates a ±v i .
The above is the reason why the flowchart, that describes the proposed architecture of the algorithm, contemplates the use of a random selection of a positive or a negative v i , just to give to the EA another degree of freedom.To clarify how the proposed algorithm works, Fig. 9 shows the pseudocode of the Earthquake Algorithm.Meanwhile Fig. 10 shows the generalities of the Earthquake Algorithm behavior, which in Then, in Fig. 10(b), the searching behavior is shown, where it can be seen that after some iterations the epicenters star to converge, with a fine search around the current global best.In spite of this, it can be observed that three epicenters took distant routes from the epicenters set, because of the random generation of positions using the exponential distribution (previously explained).
As already said, that final generation of some epicenters allowed the algorithm to "escape" the

Benchmark functions
To validate the proposed optimization method mentioned in Section 3, the performance of the Earthquake Algorithm is tested in this section by 10 benchmark functions, and its performance compared against the Particle Swarm Optimization and the Bat algorithm.The surfaces of the benchmark functions used, are shown in the Fig. 11.Additionally, their mathematical representation for each one is given by Table 2, where it is also shown the global solutions and the search domains for each one.Those were the parameters use for the benchmark of the optimization methods.
To have the comparison of the optimization methods, every algorithm was implemented with a population of 40 particles for PSO, 40 bats for BA and 40 epicenters for the EA.Also, it is important to mention that the results of the three of them were evaluated after 100 iterations, and repeated 100 times.
To quantify the results of the algorithms, Table 3 shows the mean values obtained after taking the average of the 100 results of each algorithm after the 100 iterations.
Additional to that, the second group of results are for the calculated standard deviation, and finally the best solution found of each algorithm against every function is also shown, to validate their convergence property.
Also as it can be seen in Table 3, the standard deviation of the algorithms solutions are very constant, because except against the egg holder function, the STD of the algorithm emphasizes the constant performance of the algorithm.
The reason that increases the standard deviation of the algorithm against the egg golder function, is that searching area and the multiple local solutions, make the algorithm to sometimes look for a solution in an area that is out of the possible Sr, that is why above is recommended to search an adaptive way to improve the searching area to different problems.
On all cases presented in this work (benchmark functions, model optimization and controllers optimizations), the freedom grade given by the exponential distribution random generator, helped the algorithm to send a couple of epicenters to a fast search when it seemed that the group of epicenters started to look for a convergence in a local minimum.
Those heuristic "moves" are also fundamental for the performance of the algorithm, because the capabilities to "escape" from a local solution are improved.Therefore, the paper proves that the proposed algorithm works and that it can be used as an optimization method, knowing its sturdiness against different functions or applications.

Case study
As demonstrated in Section 4, the proposed Earthquake algorithm is able to reach solutions to different benchmark functions.However, in order to discuss the capability of the algorithm against a real application, a PID speed controller for a DC motor was implemented.Additionally, in order to ensure the quality of the RPMs measurements of the DC motor a test environment was designed, where the motor could be fixed to a base, with a disc with a notch that together with an optical sensor of horseshoe type H21A1, allows to obtain a pulse counter per revolution.Fig. 13 shows the diagram of the testbed designed.

Front view Lateral view
Rear view Thus, Fig. 13 shows the testbed designed to set it (motor in red color and components of the base in white), being that the shaft was secured to a disk of two millimeters thick (disk in blue), which passes through an optical horseshoe sensor (gray device) fix to the same base, for the measurement of its RPMs.The dimensions of the motor and optical sensor are found in [50] and [51] respectively.
As far as measurements are concerned, a low pass filter was implemented into the circuit according to input and output maximum values, 5 V and 5290 RPM, respectively: where f max is the maximal frequency expected.Finally, the circuit schematic and board designs are given by Figs. 14 and 15 using the components in [50][51][52][53][54][55].The complete tests environment designed, is shown in Fig. 16, made for experimental speed analysis and control, implemented with the PPN13KB DC motor [50].

Model optimization using EA, BA and PSO
The data acquisition system and the controllers implemented, where deployed in a FPGA NI cRio-9068 [56] integrated with an I/O module NI 9381 [57], with a sample time of 6.85µSec.The principle benefits of using an FPGA, relies on its performance because of the hardware parallelism that it uses, moreover the reliability granted for avoiding the continually risk of a time-critical issue of tasks preempting one another (which is a constant on software tools).
To obtain and later optimize the model of the DC motor, a step input from 0 to 5 volts was taken to acquire the RPM response of the system, obtaining then the behavior shown Fig. 17, where can be seen on the left axis the step input in percentage of duty cycle, and on the right the speed in RPMs, both against time in seconds.Taking the data from Fig. 17, the a first experimental the continuous transfer function to be optimized (obtained by the Analytical Method (AM) [58]), is given by: And with a sample time T = 6.85µS, its discrete representation by: To obtain a quantitative analysis of the control system performance (plant model and closed loop system) to evaluate and compared the improvement of the optimization algorithms implementation, performance indexes were selected to measure the dynamic and stationary error [59].The improvement of the implementation were analyzed with four performance indexes: 1. Integral of the square of the error, 2. Integral of the absolute magnitude of the error, IAE 3. Integral of time multiplied by absolute error, 4. Integral of time multiplied by squared error, To improve the plant model obtained from analytic method, optimization algorithms are used and compared.Such algorithms are: (1) earthquake algorithm (EA), (2) bat algorithm (BA) and (3) particle swarm optimization (PSO).
The optimization process of the DC motor model is represented in Fig. 18, it shows how the EA is implemented to find the best parameters (gain and pole) with the objective to reduce the error from experimental data.
Also its important to mention that the cost function, that evaluates the performance index is based in the model of the Integral of Time multiplied by the Absolute Error (ITAE), which properties lead the optimization to an absolute adjustment.The index is calculated with the discrete model of Eq. ( 10), which is calculated on every iteration with the c2d function of Matlab TM .
For the model optimization, each algorithm was implemented with a population of 40 particles for PSO, 40 bats for BA and 40 epicenters for EA, the three for 50 iterations, and repeated 10 times.
The input model for every algorithm is given by the Eq. ( 15), and the results quantified in Table 4, where the best solution obtained from each algorithm is evaluated.Table 4 shows the values obtained for the tests.
Fig. 19 depicts the graphic behavior of the optimized transfer functions, where can be seen a comparison between the approximated models against the real measured data, all of them reacting to a step input.Besides, Table 5 compares the performance indexes obtained for every algorithm, and also shows the indexes for the analytic method.

PID optimization using EA, BA and PSO
Knowing that the proposed study case, is the speed control for a DC motor, the controller implemented and optimized by EA, BA and PSO, is a parallel PID controller.
Then, from Fig. 20 the fitness function to optimize is given by its transfer function [58]: where k p , k i and k d are the proportional, integral and derivative gains, and e(t) the error.
Thus the non optimized values, obtained by the Ziegler-Nichols Method [58], are:  • k p = 8e −6 • Representation of the optimization process for tuning PID is visualized in Fig. 21 where the Eqs.( 15) and ( 16) in closed loop are evaluated from each iterations obtaining the best parameters (k p , k i and k d ) reducing the dynamic error (ITAE) between the response and reference.
In order to strengthen the controller, and have another aspect to compare with the optimization algorithms, an adaptive PID controller using fuzzy logic was implemented.The topology implemented for the controller is given by Fig. 22.
As explained in [60], plotting the possible system inputs against the obtained outputs of the fuzzy sets, results in a control surface that can represent the entire set; reason why from the input and output sets (Fig. 23), evaluating with the proposed relation matrix (Table 6), the control surface of the system (Fig. 25) was obtained and implemented as show in Eq. ( 17).
The Fig. 24, shows how the mapping of the inputs against the outputs leads to the surface, which is finally obtained in Fig. 25.After obtaining the control surface, the curve is modeled to facilitate its implementation, adjusting it (as already mentioned) to the Eq.(17).where U(e) is the control surface, evaluated on the error e.
The general structure implemented for the PID and the Fuzzy-PID control system, can be resumed as shown in Fig. 26, where in the HMI Host was implemented the Fuzzy configuration to estimate the surface between inputs/outputs, to optimize the PID controller embedded in the cRio FPGA.Finally, the optimization algorithms already mentioned where implemented using the same population, iterations and repetitions as in Section 5.2.Nevertheless, in this section the algorithms were used first to obtain randomly a transfer function, and then its corresponding PID constants.
Taking the same criteria as in the model optimization, the selected constants were taken from the best set of solutions.Table 7, shows the values obtained for the tests.
The Fig. 27 shows a comparative between the behaviors of the implemented methods, where the PID graph represents the non-optimized PID, and the reference plot shows the input step, where the reference is placed on 4000 RPM.On the other hand, the quantification is given by Table 8, that shows the performance indexes of every method.
Evaluating the obtained parameters on a adaptive controller, Fig. 28 shows a comparative between the behaviors of the implemented methods in a fuzzy PID Mamdani type controller, where the PID graph also represents the non-optimized fuzzy PID, and the reference plot shows the input step, where the reference also for this case is placed on 4000 RPM.On the other hand, the quantification is given by Table 9, that shows the performance indexes of every method.Nevertheless, the cost function for this section is also evaluated with the performance index based in the model of the Integral of Time multiplied by the Absolute Error (ITAE), which properties lead the optimization to an absolute adjustment.Resuming in every iteration the global best as the epicenter with the minimal ITAE.

Discussion
The Earthquake Algorithm presented in this paper, was tested against two of the most implemented optimization algorithms (PSO and BA), in order to prove and compare how the EA solves optimization problems.
As most of the modern metaheuristic optimization methods, the EA is inspired on a behavior that exists in the nature, its strengths are found in the capability of using two kinds of velocity to find a solution.In section Section 4, the method was tested with different benchmark functions,   demonstrating to have the best performance against the Keane function, Rastrigin function and the Holder table, compared to the BA and PSO.
Those functions share that they are continuous, non-separable, multimodal, and somehow that their global minimums are surrounded by local solutions.The difference between them and functions like Schaffer n.2 or Levy, is in the separation of the possible solutions, because the Earthquake Algorithm proved to have a better performance that the other two, in cases where it can exploit more its versatility to adapt speed for distant searches and subsequently nearby.
Talking about convergence, the PSO showed that for functions like the Sphere function it has the best performance, but after adding relief (like in the Keane function) the PSO reduces its performance and the EA shows a better or equal performance.That property, is improved in the EA by the v s parameter, that allows the algorithm to do a finer search when needed.
However, against the other benchmark functions, the algorithm also had acceptable performances, seeing that the EA never had the worse mean result against any of the functions.Actually, the egg holder function was the most difficult test for the algorithm, though even that function was perfectly solved by the algorithm a couple of times.
Analyzing the Egg holder function, the BA clearly has the beast mean solution, but the EA also found great solutions but with a bigger standard deviation.That could be improved, with an adaptive range for the random generation with the exponential distribution, to help the algorithm prevent falling into a local minimum on a bigger range.Here (in the egg holder function) or against the table holder, the Keane or the Rastrigin functions, the v p parameter shows its strength, because that is the stochastic components Preprints (www.preprints.org)| NOT PEER-REVIEWED | Posted: 11 September 2018 doi:10.20944/preprints201809.0182.v1© 2018 by the author(s).Distributed under a Creative Commons CC BY license.

Figure 4 .
Figure 4. Travel of P-wave through an earth material by compression and dilation.

Figure 5 .
Figure 5. Travel of S-wave through an earth material by shearing deformation.

Figure 7 .
Figure 7. Procedure to find particles in and out the defined Sr.

Figure 8 .
Figure 8. Proposed architecture for the Earthquake Algorithm.

Figure 9 .
Figure 9. Pseudocode of the Earthquake Algorithm σ

Figure 10 .
Figure 10.Behavior of the Earthquake algorithm using ten epicenters.

Fig. 10 (
Fig. 10(a) appears an example of how could the epicenters randomly be placed on a function, same as when evaluated give the first fitness values.After the first rank of the best solutions, it can be known which of the epicenters are in or out the Sr and hence, which of them are going to use the v s or the v p speeds to update their positions.

Figure 11 .Table 2 .
Figure 11.Surfaces of the Benchmark functions used to compare the performance of the EA against the PSO and BA [47,48].

Preprints
(www.preprints.org)| NOT PEER-REVIEWED | Posted: 11 September 2018 doi:10.20944/preprints201809.0182.v15.1.Implementation systemThe experimental system designed for the tests consists of a DC motor with a 5V of nominal voltage, controlled by an H-Bridge[49] like the one shown in Fig.12, which also shows the resultant circuit after activating the transistor Q 5 , where the red segmented line represents the activation current flow, and the green line the resultant energy leading to the motor movement.

Figure 12 .
Figure 12.Schematic design of the implemented H-Bridge.

Fig. 15
Fig.15shows the implemented circuit board design, taken from the schematic seen in figure Fig.14, which also clearly shows a space without components, same that is used to fix the motor base of figure Fig.13to the circuit board.

Figure 17 .
Figure 17.Experimental behavior of the DC motor.

Figure 23 .
Figure 23.The fuzzy sets implemented for the Mamdani fuzzy PID controller.

Figure 26 .
Figure 26.General structure of the control System implemented

Figure 27 .
Figure 27.Comparisons of transient responses to step change in the reference with conventional PID structure.

Figure 28 .
Figure 28.Comparisons of transient responses to step change with Fuzzy PID structure.

Table 3 .
Optimization algorithms performances from different benchmark functions.

Table 4 .
Parameters obtained for plant models using each optimization algorithm.Transient response of DC motor models.

Table 5 .
Comparative table from indexes of the obtained plant models.

Table 6 .
Relation matrix for the input and output fuzzy sets.e NG out NP out Z out PP out PG out

Table 7 .
PID constants obtained for plant models using each optimization algorithm.

Table 8 .
Comparative table from indexes on conventional PID structure.

Table 9 .
Comparative table from indexes on Fuzzy PID structure.