Computational Analysis of the Active Control of Incompressible Airfoil Flutter Vibration Using a Piezoelectric V-Stack Actuator

: The ﬂutter phenomenon is a potentially destructive aeroelastic vibration studied for the design of aircraft structures as it limits the ﬂight envelope of the aircraft. The aim of this work is to propose a heuristic design of a piezoelectric actuator-based controller for ﬂutter vibration suppression in order to extend the allowable speed range of the structure. Based on the numerical model of a three degrees of freedom (3DOF) airfoil and taking into account the FEM model of a V-stack piezoelectric actuator, a ﬁltered PID controller is tuned using the population decline swarm optimizer P D SO algorithm, and gain scheduling (GS) of the controller parameters is used to make the control adaptive in velocity. Numerical simulations are discussed to study the performance of the controller in the presence of external disturbances.


Introduction
The flutter phenomenon is a dynamic instability that arises in lifting structures. It starts when the velocity becomes greater than a critical speed value, namely the flutter boundary. The fluid does not absorb the structure's energy anymore but increases it as the damping becomes negative in value, causing a condition of fatigue or damage to the structure. The amplitude of the flutter oscillations can grow indefinitely in the linear model, particularly for divergent flutter, or reach a constant amplitude, such as with limit cycle oscillations (LCOs), when there are stiffness or aerodynamic nonlinearities, like in the actual case. By studying the variation of the damping ratio with the airstream velocity, it is possible to recognize different types of flutter [1] such as soft flutter, which is characterized by a slight slope of the damping ratio-speed curve, and hard flutter, which presents an abrupt variation of the damping ratio beyond the critical point and the hump mode, which is characterized by an initial reduction of damping that becomes negative in value and then increases to become positive again, describing a closed flutter region. The danger in reaching flutter instability is not only related to the high amplitude that can be encountered, but also to the fatigue effect and the increase of the control surface's freeplay that it can cause. Over the years, many measures have been taken to prevent the rise of aeroelastic vibration phenomena. These can be mainly divided into two categories: passive and active methods. The study of mass balance belongs to the category of passive methods. The dimension and the position of masses strongly influence the amplitude and the frequency of flutter oscillations and modify the coupling mechanism between the bending and torsional vibration modes. A possible passive solution may be reached by moving the center of mass of each structural element of the wing closer to the axis of elasticity (AE) or by installing some concentrated masses to modify the wing's moment of inertia, with the aim of reducing the amount of energy extracted from the fluid [2]. The modern methods used for flutter suppression are active methods. These use the data recorded from sensors installed on the structure, such as accelerometers, that are elaborated by a control system to properly move some aerodynamic surfaces or to implement other changes in the structure shape in order to obtain the desired dynamic behavior. The first examples of technology for the control of aeroelastic oscillations were based on the movement of aerodynamic surfaces, such as the trailing edge flap, by means of servo-hydraulic systems [3]. However, hydraulic systems present some drawbacks, such as the multiple conversions of the energy from one type to another, (e.g., from mechanical to hydraulic and to mechanical again). Conversions of energy lead to losses which, in turn, require higher nominal power of the system and higher weight as well. Lastly, it has to be said that hydraulic systems usually do not reach high bandwidth values. Thus, considering that aeroelastic phenomena are characterized by high frequencies, it follows that they might not be suitable for flutter suppression. The above-mentioned problems have motivated a growing interest in deformation-induced actuators, such as piezoelectric actuators. These, in fact, allow one to obtain mechanical power directly from electric energy. However, piezoelectric actuation does not allow for obtaining large displacements; thus, amplification mechanisms are to be studied. The first attempts of piezoelectric actuation were related to the use of bimorph beams realized by bonding two piezoelectric layers such that when an electric voltage is applied, one stretches while the other contracts, and bending deformations of the beam are obtained. Applications of this technology have been studied for the active control of helicopter blade twist by Waltz and Chopra [4] and for the flutter suppression of cantilever wings by Heeg [5]. An innovative method for wing oscillation control consists of the integration of actuators into the structure such that the geometry of the structure itself can be modified, usually in terms of the airfoil camber. This technology is called morphing actuation and has been studied since the 1990s, a period in which greater knowledge of piezoelectric materials was acquired, allowing the modeling of these materials in thin patches that can be integrated into the structure. The first studies that aimed at verifying the feasibility of piezoelectric morphing actuation were carried out by Lazarus et al. [6] along with numerous other scholars. However, they focused on the realization of all-movable aerodynamic surfaces not specifically designed for aeroelastic suppression. On the contrary, the technology presented by Fichera et al. [7] was specifically designed for aeroelastic suppression. They proposed a trailing edge morphing actuator capable of a 25 Hz bandwidth and a ±15 • maximum equivalent rotation. It is made of two panels: one is an Macro-Fiber Composite-MFC piezo patch sandwich, which is the active part of the actuator, while the other is made of an aluminum alloy. Both panels are fixed on the wing box on one end and on a slider on the other end. In particular, the slider has airfoil trailing edge geometry and allows the two panels to slide in the deformed configuration. Other works dealing with piezoelectric actuation for aeroelastic suppression are based on the installation of piezo patches on the structure's surface. The actuation of the patches is combined with modern control techniques in order to obtain a robust suppression system [8].
Other simpler configurations of piezoelectric actuation deal with the use of piezoelectric stacks. Piezo stacks are made by overlapping several piezoelectric layers in order to obtain a rod able to elongate or compress after the application of an electric voltage. Since this configuration does not allow one to obtain large displacements, some amplification mechanisms have been studied. Among others, there are the so-called X and V configurations. The mechanism presented by Hall and Prechtl [9] is based on an X configuration that allows amplification of the compression deformation of the stacks, resulting in a vertical displacement. On the other hand, the amplification mechanism proposed by Ardelean and Clark [10] takes advantage of two piezoelectric stacks arranged in a V configuration. The latter configuration is the one selected for the investigation carried out in the present study.
In order to realize the flutter suppression of the airfoil, using the movement of aerodynamic surfaces, many studies have been carried out over the years, and some solutions were reported in the following. In the work of Chen et al. [11], a two-degrees-of-freedom airfoil was considered, while the presence of both the leading edge and trailing edge aero-dynamic surfaces was taken into account only in the expressions of the aerodynamic loads. The control method used by Chen et al. [11] was based on a sliding mode controller (SMC) (i.e., a nonlinear state feedback controller capable of changing its own structure when the system assumes different states). Thus, this controller presents robustness properties to external disturbances and uncertainties. However, the conventional SMC assumes that the controller could switch from one value to another instantaneously, and this is not possible in real-life applications due to the presence of delays. Therefore, the same authors proposed a high-order sliding mode controller (HOSMC) that was introduced after the backstepping procedure. In this way, a continuous input is provided to the controller [12]. In the studies of Na et al. [13] and Yang et al. [14], the model of a three-degrees-of-freedom airfoil was taken into account. Thus, the dynamic of the aerodynamic surface was considered. The control strategy proposed by Na et al. [13] was based on an LQR control implementing full-state feedback. Again, this leads to the necessity of having a system that is completely observable. The authors introduced the presence of measure noises and output disturbances, and thus a linear quadratic Gaussian (LQG) controller associated with a sliding mode observer has been implemented, since this architecture is capable of noise decoupling and shows robustness properties to disturbances. Yang et al. [14] instead specified that the LQR and LQG controllers are based on mathematical knowledge of the disturbances. Thus, when there are some uncertainties, these strategies lose robustness. In order to take into account the uncertainties in the design of the controller, Yang et al. [14] proposed H ∞ loop shaping, which is implemented by defining a generalized plant that encloses all the dynamic features to be controlled, such as sensors and actuator dynamics and so on, and introducing two different input vectors, namely the control input vector and the exogenous input vector. Then, the disturbance-output transfer function is defined, and the minimization problem of its infinity norm is solved. In the works of Wang et al. [15] and Zhang et al. [16], neural networks were associated with a full-state feedback-feedforward controller. This strategy has been studied by Bemelli-Zazzera et al. [17], too. They introduced the use of dynamic neural networks that were capable of changing their own weights through active training in order to adapt to the parameter variation during the normal functioning of the system. In the work of Andrievsky et al. [18], simple adaptive control (SAC) was implemented for the flutter suppression of an airfoil using an implicit reference model (IRM) and the pacification theorem. The authors implemented an adaptive control capable of flutter suppression in both the cases of the presence of actuator delays and the lack thereof.
Recent works dealing with active flutter suppression, carried out in the framework of projects financed by the European Union and the USA, have been presented. In the work carried out by Takarics and Vanek [19], a low-order, control-oriented LPV model of the aircraft was obtained using the bottom-up approach. The order of the nonlinear aeroservoelastic model subsystems is reduced, paying attention to capture accurately the flutter modes. Then, the controller is designed considering a polytopic LPV representation. In the work proposed by Waitman and Marcos [20], a high-fidelity, nonlinear aeroservoelastic model of the FLEXOP aircraft was considered, and three H ∞ optimal controllers were tuned and scheduled for different true airspeed values, increasing the damping of the flutter modes at these velocities. This approach led to a flutter boundary increase of 30%. Another recent work carried out by Theis et al. [21] dealt with the flutter suppression of the mini MUTT aircraft model, built at the University of Minnesota. The H ∞ loop shaping method was used in order to design a controller that presented good robustness to a wide variety of uncertainties. Pusch [22] studied a method based on the blending of control inputs and measurement outputs to isolate individual aeroelastic modes in order to use SISO controllers. In Schmidt et al. [23], the authors presented a study on drone flight dynamics, taking into account aeroelastic effects modeled by using the mean-axis formulation and quasi-steady aerodynamics. The rigid and aeroelastic model of the drone has since been used to study several controllers for autopiloting, considering the flutter instability as well.
In this work, one of the objectives is to present an FEM model of a V-stack piezoelectric actuator for the flutter suppression of an airfoil. Thus, the finite element model of the actuator is first introduced, and the behavior of the FEM model is compared with the experimental data in order to validate it. The second aim is to propose a heuristic method called population decline swarm optimization for tuning a simple control strategy based on a filtered PID controller by looking for the minimization of the integral of time absolute error (ITAE). This approach is used to set gain scheduling of the controller parameters in order to realize a velocity-dependent flutter suppression system and extend the flight envelope of the airfoil. This paper is organized as follows. In Section 2, the numerical model of a threedegrees-of-freedom airfoil is described, and the augmented state space model for the time domain analysis is introduced. In Section 3, the actuator is studied, introducing the finite element formulation of the piezo stacks. Starting from the state space model, the transfer function that relates the actuator tip's vertical displacement to the input voltage is computed. In Section 4, the tuning of the gain scheduling-filtered PID controller is done. The filtered PID controller is tuned using the P D SO algorithm by minimizing the ITAE of the error of the pitch angle. In order to set up an adaptive control in the velocity, the gain scheduling approach is introduced, and the optimization algorithm is used to find the optimal controllers for different velocity values. In Section 5, the open loop analysis of the system for the identification of the flutter boundary is first carried out. Then, the validation of the static and dynamic behavior of the actuator is done. Moreover, the convergence results of the P D SO optimization procedure are also shown. At last, the comparison between open-loop and closed-loop responses is presented, and some critical situations are simulated in order to study the performance of the closed-loop system.

3DOF Airfoil Model
In order to study the flutter phenomenon, a simplification was introduced by reducing the wing to its representative mean airfoil and, as a consequence, a two-dimensional structural model could be considered. The torsional deformability of the wing, which was considered as a beam with a thin-walled section, was simulated by a spring of stiffness K α placed in the axis of elasticity (AE), which represented the constraint. The flexional stiffness was simulated by the heave spring K y . Due to the presence of the flap a third spring of stiffness, K δ was introduced in order to take into account the hinge behavior. In this way the three-degrees-of-freedom concentrated parameters model of the wing, shown in Figure 1, was obtained. Of course, this model had to be considered as a first approximation of the structure, useful for the preliminary studies carried out in this work. In order to obtain the equations of motion of the three-degrees-of-freedom (3DOF) airfoil, the balance between the inertial, elastic and aerodynamic forces was imposed, writing the system in Equation (1) [24]. The nonlinear stiffness term was also considered in order to take into account that nonlinear effects could arise for large displacements, introducing the limit cycle oscillation phenomenon: In order to obtain the equations of motion of the three-degrees-of-freedom (3DOF) airfoil, the balance between the inertial, elastic and aerodynamic forces was imposed, writing the system in Equation (1) [24]. The nonlinear stiffness term K α3 was also considered Vibration 2021, 4 373 in order to take into account that nonlinear effects could arise for large displacements, introducing the limit cycle oscillation phenomenon: where α and δ are the pitch and flap rotation angles, respectively, while y is the plunge displacement, I α and I δ are the inertia moments, S α and S δ are the static moments, m tot is the total mass of the wing and the support blocks, L is the lift force and M and M are the aerodynamic moments acting on the airfoil and on the flap, respectively.
In order to carry out the stability analysis of the airfoil, it could be useful to take into account the nondimensional equations of motion. Thus, the first and the second equations of the system in (1) were divided by the term mb 2 , while the third equation was divided by the term mb, obtaining the following system: where ω α and ω δ , ω y are the uncoupled natural frequencies, r α and r δ are the gyration radius divided by the semi-chord of the airfoil b and γ is the nonlinear stiffness coefficient, such that K α3 = γK α . In order to obtain the state space model of the system, it was convenient to write the system in matrix form at first (see Equation (3)). The generalized displacement vector X = α δ y T , the generalized force vector F and the matrices M s and K s , namely the generalized inertia matrix and the stiffness matrix that are reported in Appendix A, were defined: M s ..
In order to introduce the structural damping contribution, the procedure described by Liu and Dowell [25] was followed, and the damping factors obtained experimentally were treated in the numerical model as modal damping factors. The structural damping matrix was written as where Λ is the modal matrix, while the modal damping matrix was defined as where m i , ω i and ζ i are the modal masses, the coupled natural frequencies and the measured damping ratios, respectively. Once the matrices were defined, the model could be written as Before defining the state space representation of the analyzed problem, the aerodynamic force model must be introduced.

Aerodynamic Model
The aerodynamic model considered in this work was the one formulated by Theodorsen [26] for an oscillating airfoil. It was based on the potential flow theory and on the hypothesis that the amplitude of the oscillations is so small that they can be treated as perturbations. According to this theory, the aerodynamic force and moments are functions of the Lagrangian parameters and of the Theodorsen constants T i , reported in Appendix B. The aerodynamic loads were written as follows: α + T 10 vδ π + bT 11 . δ 2π (9) where it appears that the Theodorsen function of the reduced frequency K = bω v , where v is the airstream velocity and ω is the frequency of the motion. The Theodorsen function was defined as a function of the Hankel function of the second kind as follows: The structural and aerodynamic formulation discussed is suitable for the analysis of the system in the frequency domain, from which information about the flutter boundary can be obtained. For the aim of this work, it was convenient to introduce time domain analysis. The Duhamel formulation [27] that allows for modeling the arbitrary motion of the airfoil instead of the simple harmonic one was thus introduced by writing the circulatory contribution to the aerodynamic loads as follows: where [27].
By integrating parts of Equation (11), a simplified expression of L c was obtained: In addition, using the Padè approximation, it could be written as a second-order differential equation, where the two aerodynamic augmented states x and .

Augmented State Space Form
The numerical model of the 3DOF airfoil belongs to the class of parameter varying systems, as it depends on the velocity of the airstream. Considering a fixed speed value, the x and writing Equations (6)-(9) in matrix form: Then, the dynamic matrix A = M −1 N was computed. The state input matrix B was defined, taking into account that the input was the flap deflection u = δ c . Thus, B was a column matrix equal to the second column of the dynamic matrix. Lastly, the output state matrix C was imposed to be the identity matrix in the hypothesis of ideal sensors for all the variables. The augmented state space model is written as follows: where the term F (x) = M −1 F (x) takes into account the nonlinear stiffness contribution. All the matrices are reported in Appendix A.

V-Stack Piezoelectric Actuator
In order to actuate the trailing edge flap for the flutter suppression of the airfoil, the V-stack piezoelectric actuator shown in Figure 2 was chosen [10]. The actuator worked by amplifying the piezoelectric stacks' deformation induced by an electric voltage. In fact, the two stacks were mounted in a V configuration, and this made the actuator easily installable inside a typical wing structure. By imposing an input voltage, one stack stretched while the other contracted, producing a vertical displacement of the actuator's tip. The vertical displacement was then transformed in a rotation of the flap through a slider [10]. The piezoelectric stacks were made by overlapping a number of piezoelectric elements with alternate poling, having alternate positive and negative electrodes printed on the external surfaces. When an electric voltage was imposed, the stacks stretched or contracted along their axes, while the other deformations were negligible. The linear electro-mechanical problem can be expressed by the equation system given by ANSI/IEEE Standard 175-1987 (i.e., the strain charge form), written as follows: where is the strain, is the stress, is the compliance matrix, is the piezoelectric matrix, is the electric field, is the electric displacement and is the dielectric constant.

Piezoelectric Stack Finite Element Model
Taking into account the particular functioning of the piezoelectric stack, Equation (16) can be reduced to the equilibrium equation along the poling direction 3 and, considering the ith piezoelectric layer, it can be written as The piezoelectric stacks were made by overlapping a number n of piezoelectric elements with alternate poling, having alternate positive and negative electrodes printed on the external surfaces. When an electric voltage was imposed, the stacks stretched or contracted along their axes, while the other deformations were negligible. The linear electro-mechanical problem can be expressed by the equation system given by ANSI/IEEE Standard 175-1987 (i.e., the strain charge form), written as follows: where ε is the strain, σ is the stress, S is the compliance matrix, d is the piezoelectric matrix, E is the electric field, D is the electric displacement and is the dielectric constant.

Piezoelectric Stack Finite Element Model
Taking into account the particular functioning of the piezoelectric stack, Equation (16) can be reduced to the equilibrium equation along the poling direction 3 and, considering the ith piezoelectric layer, it can be written as Then, in the hypothesis that every layer has the same thickness t, compliance S 33 and piezoelectric behavior d 33 , the total stack strain is where L is the length of the stack and ∆V is the electric voltage imposed. Since nt = L from Equation (18), the stress can be gained: In order to obtain the element stiffness matrix and the equivalent nodal piezoelectric force vector, the virtual work principle was introduced and, under the hypothesis of absence of external forces, it is written as where dΩ is the infinitesimal volume. Substituting Equation (19) into Equation (20) and taking into account that the strain ε and nodal displacement δ are related (i.e., ε = Bδ, where B is the finite element strain displacement matrix obtained by deriving the shape function matrix with respect to the space variable), the virtual work can be expressed as where A is the section of the stack. Due to the particular operating behavior of the stack, the axial behavior was first considered. Substituting into Equation (21) the B matrix of the linear rod finite element, B = 1 L −1 1 , the equilibrium equation along the axial direction was obtained: Thus, the element stiffness matrix k e = EA L 1 −1 −1 1 and the equivalent piezoelec-

Actuator Finite Element Model
In order to obtain the finite element model of the actuator, the simplified structural scheme shown in Figure 3a was considered, and the mesh shown in Figure 3b was introduced. With the aim of taking into account the presence of the internal hinges of the actuator, four elements of very small lengths and moments of inertia were introduced close to points 3 and 4. In this way, equilibrium was maintained while the flexional stiffness was reduced.

Actuator Finite Element Model
In order to obtain the finite element model of the actuator, the simplified structural scheme shown in Figure 3a was considered, and the mesh shown in Figure 3b was introduced. With the aim of taking into account the presence of the internal hinges of the actuator, four elements of very small lengths and moments of inertia were introduced close to points 3 and 4. In this way, equilibrium was maintained while the flexional stiffness was reduced. Each finite element of the mesh was a beam and had three degrees of freedom for each node that could be collected in the vector = , with and being the nodes of the beam while x and y are the local reference system variables along, and normal to, the beam axis oriented from i to j. The element stiffness matrix was the one computed from the Euler-Bernoulli beam. The stacks' equivalent nodal force vector was obtained, extending the one computed for the piezoelectric rod: In order to obtain the structural discretized model, a coordinate transformation from the local to the global reference frame was introduced, and the structural stiffness matrix was computed.
The structural nodal force vector is defined as follows: Each finite element of the mesh was a beam and had three degrees of freedom for each node that could be collected in the vector δ T = δ ix δ iy ϕ i δ jx δ jy ϕ j , with i and j being the nodes of the beam while x and y are the local reference system variables along, and normal to, the beam axis oriented from i to j. The element stiffness matrix was the one computed from the Euler-Bernoulli beam. The stacks' equivalent nodal force vector was obtained, extending the one computed for the piezoelectric rod: In order to obtain the structural discretized model, a coordinate transformation from the local to the global reference frame was introduced, and the structural stiffness matrix was computed.
The structural nodal force vector is defined as follows: where d , since the piezoelectric stacks work anti-symmetrically, and θ is the angle between the global X and the beam local x axes. The finite element model K ST δ = F was obtained and, by imposing the boundary conditions, the static behavior of the actuator could be studied, resolving the system: where the subscripts u and k stand for unknown and known, respectively. As the study of the dynamic behavior of the actuator was in the aim of this work, the equations of motion of the structure were written. Thus, the consistent mass matrix of the beam was introduced: Vibration 2021, 4

378
With this, the structural mass matrix M ST was computed. The structure equations of motion were written in line with the findings of Paz and Kim [28]: In order to take into account the structural damping contribution, the damping matrix was computed. For the piezo stacks' damping contribution, the dielectric loss factor tan δ was taken into account as it is related to the quality factor Q by the relation tan δ = 1 Q , and with the damping factor being ζ = 1 2Q , it could be obtained that [29] The equations of motion of the stacks in the local frame are where m e , C stack and k stack are the mass, damping and stiffness matrices in the local frame, respectively. It can be seen that there was an analogy between Equation (29) and the classic mass-spring-damper equation of motion. Therefore, the following equality can be imposed: In addition, the damping matrix is computed as where ω n is the natural frequency of the actuator modes, though only the first is retained in this work. The damping of the steel elements of the actuator was instead introduced using the Rayleigh model: where a and b are arbitrary coefficients, computed in this work using a population decline swarm optimization P D SO algorithm, described in Appendix C. The P D SO algorithm was implemented for the identification of the actuator experimental frequency response, presented by [10], choosing as an objective function the percentage error %E M r between the experimental resonance peak M r = 3.95 × 10 −2 mm V and the one computed from the finite element model. The optimization problem is written as . The research space limits were chosen by a trial and error approach, looking for minimization of the percentage error. The algorithm parameters were chosen as follows: c c = c s = 2.025, µ min = 0.4, µ max = 0.9, Λ = 40, ζ = 0.5 and ∆ = 10. The minimum value of the objective function min%E M r = 0.0055% was already obtained for an initial maximum number of particles equal to 20, and the minimum point had the coordinates a = 1.
By introducing the damping contribution and imposing the boundary conditions, the equations of motion of the actuator are written as In order to obtain the theoretical frequency response, the transfer function of the actuator was computed. The state space model was introduced, defining the state vector δ i and the following matrices: where n is the degree of freedom of the actuator. The C matrix was defined in such a way that the model returned as an output the vertical displacement of point 6 (i.e., the tip of the actuator). The transfer function is computed as Lastly, it must be pointed out that to model the internal hinges of the actuator, four finite elements having a length equal to 1% of the membership elements were introduced. The fictitious finite element moments of inertia were obtained by dividing the membership element ones by a factor f, chosen in such a way that the minimum error on the static and dynamic behavior was obtained. For the sake of clarity, looking at Figure 3b, the membership element of the fictitious hinge finite element 10-3 is the finite element 1-10. Figure 4 shows the percentage errors computed between the FE model and the experimental results of the piezoelectric V-stack actuator stroke, namely δ 6Y , of its blocked force F 6Y and of the resonance frequency as a function of the fictitious hinge factor f. Looking at Figure 4, it appears that the optimal dividing factor choice was f = 10 6 .

Gain-Scheduled Filtered PID Active Controller
The control architecture chosen in this work and shown in Figure 5 was based on an output feedback scheme based on a filtered Proportional-Integral-Derivative fPID controller that presented some benefits. This type of controller presented a simplicity in implementation and installation and, moreover, it only requested one output measure, thus asking for the implementation of only one sensor device and allowing us to avoid the complication of the controller scheme introduced by, for instance, state observers. These are some of the reasons that currently still motivate the use of PID controllers in the industry [30]. Moreover, the PID controller had a fixed structure, and the problem reduced to the tuning of its parameters, which were made in this work using a population decline swarm optimization algorithm. The principal issue of choosing a PID controller is that it is not adaptive with the system's change of parameters. In this work, in order to make the controller adaptive in velocity, a gain scheduling (GS) approach was used.

Gain-Scheduled Filtered PID Active Controller
The control architecture chosen in this work and shown in Figure 5 was based on an output feedback scheme based on a filtered Proportional-Integral-Derivative fPID controller that presented some benefits. This type of controller presented a simplicity in implementation and installation and, moreover, it only requested one output measure, thus asking for the implementation of only one sensor device and allowing us to avoid the complication of the controller scheme introduced by, for instance, state observers. These are some of the reasons that currently still motivate the use of PID controllers in the industry [30]. Moreover, the PID controller had a fixed structure, and the problem reduced to the tuning of its parameters, which were made in this work using a population decline swarm optimization P D SO algorithm. The principal issue of choosing a PID controller is that it is not adaptive with the system's change of parameters. In this work, in order to make the controller adaptive in velocity, a gain scheduling (GS) approach was used.
These are some of the reasons that currently still motivate the use of PID controllers in the industry [30]. Moreover, the PID controller had a fixed structure, and the problem reduced to the tuning of its parameters, which were made in this work using a population decline swarm optimization algorithm. The principal issue of choosing a PID controller is that it is not adaptive with the system's change of parameters. In this work, in order to make the controller adaptive in velocity, a gain scheduling (GS) approach was used. The filtered PID controller acted by minimizing the error function ( ) = ( ) − ( ), being ( ) = 0 , because the aim was to stabilize the system. The reference signal sent to the actuator was the voltage input, and it was computed as the sum of three terms: one proportional to the current error, one proportional to the time integral of the error, ensuring the system reached the desired setpoint, and the last term proportional to the filtered error time derivative [30]: where ( ) is the solution of the following first-order differential equation: The choice of a controller with all three components was related to the peculiarities of each term. In fact, a proportional-only controller could not make the system able to reach the setpoint. By increasing the proportional gain, the difference between the actual and desired values could be reduced, but the system would encounter more oscillations The filtered PID controller acted by minimizing the error function e(t) = α re f (t) − α(t), being α re f (t) = 0rad, because the aim was to stabilize the system. The reference signal sent to the actuator was the voltage input, and it was computed as the sum of three terms: one proportional to the current error, one proportional to the time integral of the error, ensuring the system reached the desired setpoint, and the last term proportional to the filtered error time derivative [30]: where e D (t) is the solution of the following first-order differential equation: The choice of a controller with all three components was related to the peculiarities of each term. In fact, a proportional-only controller could not make the system able to reach the setpoint. By increasing the proportional gain, the difference between the actual and desired values could be reduced, but the system would encounter more oscillations for the rapid transient; thus, the integral contribution was added. At least, the derivative terms allowed the system to tackle the rapid transient easily, as the derivative gave information on the evolution of the error. Thus, the controller was able to forecast and fill the eventual delays, leading to a reduction of the oscillation amplitude.
The controller parameters, namely K p , τ i , τ D and τ DF , had to be tuned in such a way that no instabilities due to the control were introduced, and the flap displacement had to be able to stabilize the airfoil in the shortest time. In order to obtain the PID controller at the flutter speed, a P D SO algorithm was used, and the fitness function to be minimized was chosen to be the integral of time absolute error (ITAE).
It is worth noting that that, strictly speaking, heuristic optimizers did not find the actual optimum, but their solution was always the best parameter configuration among all the other configurations that were tested. This means that the controller might be suboptimal because it is just the best among the tested ones. This is particularly true when a flat region of the performance index is found. Thus, considering that the PID tuning was gained via a numeric metaheuristic approach, the obtained controller was to be considered as a suboptimal one. With regard to the objective function, it should be noted that the ITAE was chosen because it had the ability to highlight the divergence error that occurred when incorrect controller parameters were used and the flutter vibration started increasing slowly. The mathematical expression of the fitness function is the following: where t 0 is the instant of perturbation and T w is the time window chosen for the computation of the function, being large enough to include the entire transient response, but not Vibration 2021, 4 381 too much to affect the computational cost of the algorithm. In this study, it was imposed to be t 0 = 0 s and T w = 1 s. Each particle P i at each iteration λ was defined in the algorithm by four coordinates , and the minimization problem was expressed as follows: where the research space limits were chosen by a preliminary trial and error approach.
As the speed can vary, the 3DOF airfoil model is a parameter varying system, and for this family of dynamic systems, the simplest method to obtain adaptive controllers is the gain scheduling (GS) method. The GS method was implemented by defining a set of controllers tuned for different working conditions of the system, which were identified by an appropriate scheduling variable, in this case the airstream speed. However, each controller guaranteed the desired performance only locally, and thus an interpolation of the controller parameters had to be introduced [31]. The gain scheduling control architecture used is shown in Figure 6. In order to compute the scheduling table, the algorithm was used for different velocity values, and the optimal PID controller for each one was found. A set of linear, time-invariant sub-optimal controller parameters was obtained, and a linear interpolation of the PID parameters was introduced to obtain a parameter-varying controller.

Numerical Results
In this section, the numerical results of the aeroelastic model and the FEM actuator model are analyzed. First, the open-loop analysis of the 3DOF airfoil is carried out, and then the validation of the FEM model of the actuator is done. Lastly, some critical case studies for the closed-loop system are analyzed.

Open-Loop Analysis of the Aeroelastic System
The experimental model considered in this work was the one studied by Ardelean et al. [4] at Duke University, and its parameters are listed in Table 1. In order to study the stability of the airfoil, the eigenvalue problem det( − ) = 0 was solved, and the root locus, shown in Figure 7, was obtained, varying the airstream velocity. It can be seen that the system had three pairs of complex eigenvalues and a pair of real eigenvalues. Instability was reached when the real part of one complex eigenvalues pair became In order to compute the scheduling table, the P D SO algorithm was used for different velocity values, and the optimal PID controller for each one was found. A set of linear, time-invariant sub-optimal controller parameters was obtained, and a linear interpolation of the PID parameters was introduced to obtain a parameter-varying controller.

Numerical Results
In this section, the numerical results of the aeroelastic model and the FEM actuator model are analyzed. First, the open-loop analysis of the 3DOF airfoil is carried out, and then the validation of the FEM model of the actuator is done. Lastly, some critical case studies for the closed-loop system are analyzed.

Open-Loop Analysis of the Aeroelastic System
The experimental model considered in this work was the one studied by Ardelean et al. [4] at Duke University, and its parameters are listed in Table 1. In order to study the stability of the airfoil, the eigenvalue problem det(λI − A) = 0 was solved, and the root locus, shown in Figure 7, was obtained, varying the airstream velocity. It can be Vibration 2021, 4 382 seen that the system had three pairs of complex eigenvalues and a pair of real eigenvalues. Instability was reached when the real part of one complex eigenvalues pair became positive in value.  In order to determine the flutter speed, the damping ratio variation with the airstream velocity was studied, obtaining the diagram shown in Figure 8, where it can be seen that the damping ratio related to the unstable eigenvalue became negative in value for a speed equal to = 19 while the flutter frequency was = 4.2 Hz.
These results are in accordance with the ones presented by Ardelean et al. [32]. When the flutter speed was exceeded, the limit cycle oscillation (LCO) phenomenon arose. Thus, the oscillations presented a growing amplitude with the velocity. In Figure 9, the LCO amplitudes on with the airstream velocity are shown. In order to determine the flutter speed, the damping ratio variation with the airstream velocity was studied, obtaining the diagram shown in Figure 8, where it can be seen that the damping ratio related to the unstable eigenvalue became negative in value for a speed equal to v f lutter = 19 m s while the flutter frequency was ω f lutter = 4.2 Hz. These results are in accordance with the ones presented by Ardelean et al. [32]. In order to determine the flutter speed, the damping ratio variation with the airstream velocity was studied, obtaining the diagram shown in Figure 8, where it can be seen that the damping ratio related to the unstable eigenvalue became negative in value for a speed equal to = 19 while the flutter frequency was = 4.2 Hz.
These results are in accordance with the ones presented by Ardelean et al. [32]. When the flutter speed was exceeded, the limit cycle oscillation (LCO) phenomenon arose. Thus, the oscillations presented a growing amplitude with the velocity. In Figure 9, the LCO amplitudes on with the airstream velocity are shown. When the flutter speed was exceeded, the limit cycle oscillation (LCO) phenomenon arose. Thus, the oscillations presented a growing amplitude with the velocity. In Figure 9, the LCO amplitudes on α with the airstream velocity are shown.

Actuator Validation
The curves of the tip's vertical displacement and blocked force versus the electric voltage obtained from Equation (25) were compared with the theoretical ones presented by Ardelean et al. [10]. This is shown in Figure 10, where it can be seen that the FEM and theoretical curves almost overlapped. In fact, errors of 0.15% for the stroke and 0.055% for the blocked force were obtained. The actuator's data are listed in Table 2.

Actuator Validation
The curves of the tip's vertical displacement and blocked force versus the electric voltage obtained from Equation (25) were compared with the theoretical ones presented by Ardelean et al. [10]. This is shown in Figure 10, where it can be seen that the FEM and theoretical curves almost overlapped. In fact, errors of 0.15% for the stroke and 0.055% for the blocked force were obtained. The actuator's data are listed in Table 2.

Actuator Validation
The curves of the tip's vertical displacement and blocked force versus the electric voltage obtained from Equation (25) were compared with the theoretical ones presented by Ardelean et al. [10]. This is shown in Figure 10, where it can be seen that the FEM and theoretical curves almost overlapped. In fact, errors of 0.15% for the stroke and 0.055% for the blocked force were obtained. The actuator's data are listed in Table 2.
In order to validate the dynamic behavior of the actuator, the bode diagram of the transfer function computed in Equation (35) was compared with the experimental frequency response of the actuator presented by Ardelean et al. [10]. From Figure 11, it can be seen that in the bandwidth region extending up to 200 Hz, the two curves overlapped, while some differences appeared beyond the first natural frequency of the actuator. These differences were due to the actuator's parameters' uncertainty and the difference between the theorical and real values of the materials' constants. However, it is worth noting that such a discrepancy would not affect the closed-loop system, while the flutter frequency of the airfoil fell in the bandwidth region. Therefore, the FEM model was acceptable for the frequency range of interest.  In order to validate the dynamic behavior of the actuator, the bode diagram o transfer function computed in Equation (35) was compared with the experime frequency response of the actuator presented by Ardelean et al. [10]. From Figure 1 can be seen that in the bandwidth region extending up to 200 Hz, the two cu overlapped, while some differences appeared beyond the first natural frequency o actuator. These differences were due to the actuator's parameters' uncertainty and difference between the theorical and real values of the materials' constants. However worth noting that such a discrepancy would not affect the closed-loop system, whil flutter frequency of the airfoil fell in the bandwidth region. Therefore, the FEM m was acceptable for the frequency range of interest.

GS-fPID Controller Tuning
The minimization problem defined in Equation (39) was solved using the algorithm. First of all, the research space limits were defined as 0.10 ≤ ≤ 2 For a larger research spa higher number of particles was needed to achieve the convergence and computational cost of the procedure's increases. The algorithm parameters were ch as was done before for the damping contribution in the FEM model of the actuator, four optimizations were done with different maximum numbers of particles 5 10 20 50 in order to study the convergence of the procedure. For optimization, the algorithm was run 10 times, and the final values of the obje

GS-fPID Controller Tuning
The minimization problem defined in Equation (39) was solved using the P D SO algorithm. First of all, the research space limits were defined as 0.
For a larger research space, a higher number of particles was needed to achieve the convergence and the computational cost of the procedure's increases. The algorithm parameters were chosen as was done before for the damping contribution in the FEM model of the actuator, and four optimizations were done with different maximum numbers of particles P max = 5 10 20 50 in order to study the convergence of the procedure. For each optimization, the algorithm was run 10 times, and the final values of the objective function were analyzed to compute the minimum, maximum, median and standard deviation. These convergence indexes are shown in Table 1, while in Figure 12, the ITAE trends for any initial swarm size are shown.   Table 4. From the results reported in Table 4, it can be seen that in the speed range higher than 19 m/s, i assumed the research space upper limit value. It has been noted that, by increasing the upper limit, i tended to always assume its higher value, while the other parameters remained unchanged. This behavior indicates that, in such a range of velocities, the integral contribution to the controller was negligible when compared with the proportional and derivative ones. Thus, the heuristic tuning procedure suggests that a PD scheme would be more suitable From the results shown in Table 1, it can be seen that the minimum objective function value was ITAE min = 4.51E −5 rad * s, and it was already obtained for a maximum number of particles P max = 20. The coordinates of the point of minimum were K p = 7.8E 3 V deg ; τ i = 1E 4 s; τ D = 0.255 s; τ DF = 1 s. The values of the time constants were the same for each initial swarm size, so they are not reported in Table 3. It can be seen that, by increasing the initial swarm size, the possibility to lock into local minima reduced (see Figure 12), and the values of the convergence indexes reduced because of the higher research capability. However, the CPU time increased, as expected. In order to tune the gain scheduling adaptive controller, the P D SO algorithm was used for each speed value of the vector v = 15 18 18.5 19.5 20 25 28 , obtaining the scheduling table reported in Table 4. From the results reported in Table 4, it can be seen that in the speed range higher than 19 m/s, τ i assumed the research space upper limit value. It has been noted that, by increasing the upper limit, τ i tended to always assume its higher value, while the other parameters remained unchanged. This behavior indicates that, in such a range of velocities, the integral contribution to the controller was negligible when compared with the proportional and derivative ones. Thus, the heuristic tuning procedure suggests that a PD scheme would be more suitable in such speed range. Moreover, looking at Table 4, it can also be noted that τ D and τ DF reached, in some cases, the limiting values of their research spaces. By carrying out numerical tests, it was noted that if a wider research space was selected for τ D and τ DF , several pairs of τ D and τ DF values would give the same minimum value of the objective function. This indicates a flat region of the objective function with respect to these variables. Moreover, regarding the tendencies of the obtained scheduling values, it must be noted that the heuristic approaches were useful because they allowed for finding good solutions despite the complexity of the objective function. Moreover, they also allowed for finding a good solution even when the shape of the objective function was not known. This was the case when a performance index, such as the ITAE, was selected as an objective function. The physical meaning of the performance index must, of course, be clear, but its behavior with respect to both the tuning parameters (such as the PID gains in the present work) and the system's parameters (such the velocity) is generally unknown or hard to obtain analytically (particularly when nonlinearities are also present). The metaheuristic algorithms allow for treating such optimization problems in an easy way, but to know the reason why some solutions are obtained using the metaheuristic optimization is hard, and this is what it takes. The designer just has to ensure that the obtained solution is meaningful from physical and technological points of view. In particular, one of the points to be addressed is the stability of the controlled plant. In this work, by means of an engineering approach, it was verified by carrying out numerical simulations of the nonlinear closed-loop system and looking for PID parameter sets that ensure the local linearized plant stability (aware of the fact that the local stability may not be straightforward in ensuring the stability of the time-varying systems that can ask for ad hoc adaptive approaches). Figure 13 shows the damping ratio and the natural frequencies of the closed-loop system as a function of the advancing velocity. The modes associated with the structural V-stack mechanism were omitted for the sake of clarity of the diagram. In particular, when looking at the closed-loop modal damping ratio trend, it can be seen that a pair of conjugate poles became unstable for velocities higher than 32 m s . It is worth noting that to obtain the modal characteristics shown in Figure 13, the closed-loop system was linearized. However, when the actuator saturation was taken into account, the controller was no longer able to stabilize the system for speed values higher than 28 m s , as will be shown in the next numerical studies. Thus, the flutter boundary of the closed loop system was 28 m s . In Figure 14, the pole map of the PID transfer function is shown. The PID transfer function had two poles; one pole was always equal to zero, while the other changed with the airstream velocity, according to the change of the PID parameters. In fact, starting from a speed of 15 m s , the normalized frequency of the real non-zero pole assumed a value of ω n = 0.074 Hz, and when the open-loop flutter boundary was reached, it assumed the value ω n = 0.159 Hz, which was maintained until a speed of 25 m s . This stemmed from the absence of poles that were too fast, allowing its implementation into a flight computer which could have a sampling of 5 ms (see, for instance, Waitman and Marcos [20]).

Closed-Loop Analysis
The PID parameters obtained from the optimization done in Section 5.3 were then used to run a simulation of the closed-loop system at the flutter speed. The initial state vector for the simulation was = 0.05 0.025 0 0 0 0 0 0 , and the responses shown in Figure 15 were obtained. It can be seen that the system was able to suppress the flutter oscillations on in less than 0.1 seconds without any overshoot, while for the

Closed-Loop Analysis
The PID parameters obtained from the optimization done in Section 5.3 were then used to run a simulation of the closed-loop system at the flutter speed. The initial state vector for the simulation was = 0.05 0.025 0 0 0 0 0 0 , and the responses shown in Figure 15 were obtained. It can be seen that the system was able to suppress the flutter oscillations on in less than 0.1 seconds without any overshoot, while for the

Closed-Loop Analysis
The PID parameters obtained from the optimization done in Section 5.3 were then used to run a simulation of the closed-loop system at the flutter speed. The initial state vector for the simulation was x T 0 = 0.05 0.025 0 0 0 0 0 0 , and the responses shown in Figure 15 were obtained. It can be seen that the system was able to suppress the flutter oscillations on α in less than 0.1 seconds without any overshoot, while for the plunge displacement and the flap deflection, the maximum values were y max = 0.0011 m and δ max = 14.9 • , respectively.
Vibration 2021, 4, 24 3 plunge displacement and the flap deflection, the maximum values were = 0.0011 and = 14.9°, respectively. As in the first case study for the closed-loop system with a gain-scheduling adaptiv controller, a linear variation of the speed, starting from 15 and with a 2 slope, w As in the first case study for the closed-loop system with a gain-scheduling adaptive controller, a linear variation of the speed, starting from 15 m s and with a 2 m s 2 slope, was considered, and the responses shown in Figures 16 and 17 were obtained. It can be seen that the initial perturbation was rapidly suppressed, as clearly shown in Figure 17, where the time window was limited to one second of analysis. Then, at the time instant t = 2 s, the flutter boundary of the open-loop system was reached. Thus, it started oscillating, with increasing LCOs as the speed increased. On the other hand, the closed-loop system did not show oscillations until the time instant t = 6.5 s, when flap-unstable oscillations arose, suggesting that the closed-loop flutter boundary was exceeded. Looking at the velocity diagram in Figure 16, the closed loop-flutter boundary is represented by the speed of 28 m s , as previously stated. It can be said that the adaptive controller was able to stabilize the system in the new flight envelope assigned. As in the first case study for the closed-loop system with a gain-scheduling adaptive controller, a linear variation of the speed, starting from 15 and with a 2 slope, was considered, and the responses shown in Figure 16 and Figure 17 were obtained. It can be seen that the initial perturbation was rapidly suppressed, as clearly shown in Figure 17, where the time window was limited to one second of analysis. Then, at the time instant = 2 s, the flutter boundary of the open-loop system was reached. Thus, it started oscillating, with increasing LCOs as the speed increased. On the other hand, the closed-loop system did not show oscillations until the time instant = 6.5 s, when flap-unstable oscillations arose, suggesting that the closed-loop flutter boundary was exceeded. Looking at the velocity diagram in Figure 16, the closed loop-flutter boundary is represented by the speed of 28 , as previously stated. It can be said that the adaptive controller was able to stabilize the system in the new flight envelope assigned. Looking at Figure 8, it can be noted that the instability arose after an abrupt variation of the damping ratio, starting from the velocity value = 15 . Thus, it was useful to study the controller performance in such a critical speed range. The controller stress case was chosen to be a speed ramp from 15 to 28 with different acceleration values, while a constant disturbance on the flap deflection equal to 5° occurred. From the responses shown in Figure 18 , it can be seen that the transient response was almost equal for all the accelerations studied, and the unstable oscillations that occurred after reaching the closed-loop flutter boundary had the same amplitudes. Therefore, the controller reaction was considered good for both the considered velocity variation and the disturbance. Looking at Figure 8, it can be noted that the instability arose after an abrupt variation of the damping ratio, starting from the velocity value v = 15 m s . Thus, it was useful to study the controller performance in such a critical speed range. The controller stress case was chosen to be a speed ramp from 15 m s to 28 m s with different acceleration values, while a constant disturbance on the flap deflection equal to 5 • occurred. From the responses shown in Figure 18, it can be seen that the transient response was almost equal for all the accelerations studied, and the unstable oscillations that occurred after reaching the closed-loop flutter boundary had the same amplitudes. Therefore, the controller reaction was considered good for both the considered velocity variation and the disturbance. values, while a constant disturbance on the flap deflection equal to 5° occurred. From the responses shown in Figure 18 , it can be seen that the transient response was almos equal for all the accelerations studied, and the unstable oscillations that occurred after reaching the closed-loop flutter boundary had the same amplitudes. Therefore, the controller reaction was considered good for both the considered velocity variation and the disturbance. As in the last closed loop case study, a sinusoidal velocity variation within the range ∈ 15 28 was imposed, while a random disturbance on the flap deflection ∈ −5°5° occurred. From the responses shown in Figure 19 , it can be seen that the controller was able to face up to the simultaneous variation of speed and the disturbances, avoiding instability conditions and maintaining the state variable peaks a acceptable values from a mechanical point of view. In fact, it can be seen that during the transient response, the rotation of the airfoil assumed the maximum value = −1.1° Figure 18. Controller stress case.
As in the last closed loop case study, a sinusoidal velocity variation within the range v ∈ 15 m s 28 m s was imposed, while a random disturbance on the flap deflection δ d ∈ −5 • 5 • occurred. From the responses shown in Figure 19, it can be seen that the controller was able to face up to the simultaneous variation of speed and the disturbances, avoiding instability conditions and maintaining the state variable peaks at acceptable values from a mechanical point of view. In fact, it can be seen that during the transient response, the rotation of the airfoil assumed the maximum value α max = −1.1 • , the flap deflection was limited to δ max = 15 • , and the plunge displacement presented a pick equal to y max = −10 −3 m. In addition, all three responses were stabilized to the equilibrium point in almost 0.7 seconds, presenting slight fluctuations due to the presence of the disturbance.

Conclusions
In this paper, the numerical model of a three-degrees-of-freedom airfoil was studied in order to predict the flutter boundary of the system. With the aim of taking into account the piezoelectric V-stack actuator, a finite element model was introduced.  Figure 19. Responses to the sinusoidal velocity variation and random disturbance.

Conclusions
In this paper, the numerical model of a three-degrees-of-freedom airfoil was studied in order to predict the flutter boundary of the system. With the aim of taking into account the piezoelectric V-stack actuator, a finite element model was introduced. The validation of the static and dynamic behavior of the actuator with the experimental data given in the literature was carried out. A heuristic swarm method named P D SO was employed for the tuning of a gain-scheduling adaptive controller based on a filtered PID controller. An increase of the flutter boundary of the system of about 47% was obtained. Lastly, in order to study the performance of the closed-loop system, several simulations were carried out, considering critical values of acceleration and speed and taking into account the presence of disturbances. The actuator FEM model computed in this work can represent a starting point for the study of the influence of materials and geometric configurations in the V-stack actuator flutter suppression capability. Moreover, the airfoil flutter suppression P D SO tuning approach presented in this work could be used, for instance, as a reference approach for a less simple model or hardware in the loop control tuning of an aerodynamic surface, using the local sensors' output time histories to define the objective function to be minimized.

Conflicts of Interest:
The authors declare no conflict of interest.

Population Decline Swarm Optimization: P D SO
The P D SO is a variant of particle swarm optimization (PSO). It is a stochastic optimization procedure based on a social model for which a swarm of particles is defined, and for each one, a set of variables is assigned. The aim of the algorithm is to find the values of the parameters associated with the minimum point of the objective function [30]. At the first iteration step, a set of coordinates inside the research space are defined, and a set of random directional velocity values are assigned to each particle. Then, the objective function is computed for each particle and it is assigned as an individual best value, while the minimum objective function value among all the particles is searched for to define the global minimum. At the end of the iteration, the new coordinates and velocities of the particles are computed in order to make them follow the best solution. During successive iterations, it is verified that the new coordinates fall in the research space defined, and the individual and global best solutions are updated. The algorithm stops when a convergence condition is achieved or when the maximum iterations Λ are reached. The flowchart of the algorithm is shown in Figure A1.
The P D SO is a PSO algorithm, where the reduction of the number of particles is considered during its evolution. In fact, a decline factor ζ < 1 is defined, and when a certain number of iterations ∆ is reached, the population is reduced (P = ζP max ). While a higher maximum number of particles is related to greater research capability and a lower probability of the local minimum lock problem, it leads to a higher computational cost of the procedure. A proper setting of parameters allows for a best compromise between the initial global search capability and final faster exploitation ability of the swarm. A reinitialization of some particle coordinates and velocities can be introduced together with the population reduction. In this way, new trajectories of the particles are explored, reducing the local minimum lock events.
The mathematical formulation of the algorithm is explained in the following. First of all, a number n of parameters ξ is assigned to each particle.

Population Decline Swarm Optimization:
The is a variant of particle swarm optimization (PSO). It is a stochastic optimization procedure based on a social model for which a swarm of particles is defined, and for each one, a set of variables is assigned. The aim of the algorithm is to find the values of the parameters associated with the minimum point of the objective function [30]. At the first iteration step, a set of coordinates inside the research space are defined, and a set of random directional velocity values are assigned to each particle. Then, the objective function is computed for each particle and it is assigned as an individual best value, while the minimum objective function value among all the particles is searched for to define the global minimum. At the end of the iteration, the new coordinates and velocities of the particles are computed in order to make them follow the best solution. During successive iterations, it is verified that the new coordinates fall in the research space defined, and the individual and global best solutions are updated. The algorithm stops when a convergence condition is achieved or when the maximum iterations Λ are reached. The flowchart of the algorithm is shown in Figure A1.
The is a PSO algorithm, where the reduction of the number of particles is considered during its evolution. In fact, a decline factor 1 is defined, and when a certain number of iterations Δ is reached, the population is reduced ( = ). While a higher maximum number of particles is related to greater research capability and a lower probability of the local minimum lock problem, it leads to a higher computational cost of the procedure. A proper setting of parameters allows for a best compromise between the initial global search capability and final faster exploitation ability of the swarm. A reinitialization of some particle coordinates and velocities can be introduced together with the population reduction. In this way, new trajectories of the particles are explored, reducing the local minimum lock events.
The mathematical formulation of the algorithm is explained in the following. First of all, a number of parameters is assigned to each particle.  where λ = 1, 2, . . . , Λ is the actual iteration step. Thus, after the objective function is computed for each particle, their positions are updated: where v i λ+1 is the updated velocity, defined as where c c and c s are the cognitive and social acceleration coefficients, respectively; r 1 and r 2 are random coefficients inside the range 0 1 ; and µ λ is the inertia factor introduced to balance the local and global searches, varying linearly with the iterations according to the following relation: At last, χ is the constriction factor introduced to avoid the explosion of the swarm, computed as where φ = c c + c s > 4 is imposed in order to ensure the stability of the algorithm [33]. The population reduction model is defined as suggested by Orlando and Alaimo [30]: where the notation x stands for the integer part of x, such that x = − − x.