Position Control of the Dielectric Elastomer Actuator Based on Fractional Derivatives in Modelling and Control

Successful control of a dielectric elastomer actuator (DEA) can be a challenging task, especially if no overshoot is desired. The work presents the first use of the PIλDμ control for a dielectric elastomer actuator to eliminate the overshoot. The mathematical model of the dielectric elastomer was established using the fractional Kelvin-Voigt model. Step responses are first tested in the Laplace domain, which gave the most satisfactory results. However, they did not represent the real model. It cannot have negative force acting on the dielectric elastomer actuator. Simulations in Matlab/Simulink were performed to obtain more realistic responses, where output of the PIλDμ controller was limited. Initial parameters for a PID control were obtained by the Wang–Juang–Chan algorithm for the first order plus death time function approximation to the step response of the model, and reused as the basis for the PIλDμ actuator control. A quasi-anti-windup method was introduced to the final control algorithm. Step responses of the PID and the PIλDμ in different domains were verified by simulation and validated by experiments. Experiments proved that the fractional calculus PIλDμ step responses exceeded performance of the basic PID controller for DEA in terms of response time, settling time, and overshoot elimination.


Introduction
Soft actuators are actuators whose performance tries to mimic the behavior of a biological muscle. A biological muscle performs only two simple motions: contraction and relaxation. Even though this movement seems very simple, the reality is quite the opposite. The basic structure of a biological muscle is a bundle of muscle fibers. Muscle fibers are a compound of myofibrils, of which the basics are sarcomeres. A sarcomere has two parts: myosin filaments and actin filaments. They interact with each other with the help of Ca 2+ ions. The presence of Ca 2+ ions contracts the muscle fiber by δL, whereby the duration of the contraction depends on the amount of these ions [1]. Contraction of the muscle generates F force; the product of contraction and force is A = δL · F. If A is measured over time, one can obtain power P = F · δL/δt generated by the contraction of a muscle.
Regarding the activation principle, soft actuators that are most like biological muscles are ionic electroactive polymers [2,3]. They deform their principal shape when there is an internal difference in the ionic state. Even though the principle is similar, its deformation differs from biological muscles. They usually exhibit bending movement with a very slow deformation rate. Activation energy to change ionic state is applied by electrical energy. It is a reverse process, but it has a slow recovery rate.
On the other hand, there is another group of soft actuators called dielectric elastomer actuators, also known as DEAs. Their activation energy is also electricity with 1000 times higher potential then in ionic actuators. The principle of the activation of the structure differs from ionic actuators. DEA has a similar structure as a parallel plated capacitor. It has upper and lower conductive electrodes with elastomer used as the dielectric. Elastomer is considered incompressible. Its Poisson number is 0.5. Only direct current can be applied on a DEA. When DC voltage is applied, attractive force is generated between the positive and negative potentials. This is also called the Maxwell stress. Figure 1 shows the basic structure and principle of activation of a planar DEA. When DC voltage is applied, the elastomer dimensions change from 1 → 1 ', 2 → 2 ′, and 3 → 3 ′ according to Figure 1. The ratio of change in each direction can be described as 1 = 1 / 1 ′, 2 = 2 / 2 ′, and 3 = 3 / 3 ′, respectively. Since elastomer is incompressible, the following applies: 1 = 1 2 3 . The basic principles of DEAs' movement are not the same as for biological muscles, but, in comparison with ionic actuators, they move much faster and exhibit larger movements. Under applied voltage they can also hold the tension or position. As such, they can be used as soft actuators with similar behavior as biological muscles.
The introduction of fractional calculus dates back to the 17th century, where Leibniz addressed a letter to L'Hopital asking him about the existence of a half derivative [4]. This date is marked as the beginning of fractional calculus. Fractional calculus defines integrals and derivatives with non-integer orders, so the name might be a bit confusing, but since then, this term has been generally accepted. Two centuries later the proper definition was developed, with the works of Riemann et al. [5].
Fractional calculus can be used in different areas. Lately it is being used in describing the bursting of the two pancreatic β-cells working in synchronization. With the help of changing only fractional orders, they are able to model different bursting activities of insulin [6]. Studying biological conditions with memory characteristics that can be defined with the fractional calculus is another interesting area of use [7]. Readers interested in areas where fractional calculus is expanding should refer to the work of Sun et al. [8]. Among others, the main areas of research with fractional calculus are in physics, control, signal and image processing, mechanics and dynamic systems, biology, environmental science, materials, viscoelasticity, and economics [8]. Fractional calculus is also being used to describe fractional order chaotic systems, systems that are very sensitive to initial conditions [9]. One such interesting area for using fractional calculus is rheology, part of which is viscoelasticity. The behavior of viscoelastic materials can be simulated using fractional calculus in rheological models, where basic rheological models can be converted into fractional ones [10]. It can also be used in control theory to convert a classical PID controller into a fractional one [5,11,12]. Expressed as λ , the fractional controller has five tuning parameters, which enables a wider fitting range for fine tuning. The basic principles of DEAs' movement are not the same as for biological muscles, but, in comparison with ionic actuators, they move much faster and exhibit larger movements. Under applied voltage they can also hold the tension or position. As such, they can be used as soft actuators with similar behavior as biological muscles.
The introduction of fractional calculus dates back to the 17th century, where Leibniz addressed a letter to L'Hopital asking him about the existence of a half derivative [4]. This date is marked as the beginning of fractional calculus. Fractional calculus defines integrals and derivatives with non-integer orders, so the name might be a bit confusing, but since then, this term has been generally accepted. Two centuries later the proper definition was developed, with the works of Riemann et al. [5].
Fractional calculus can be used in different areas. Lately it is being used in describing the bursting of the two pancreatic β-cells working in synchronization. With the help of changing only fractional orders, they are able to model different bursting activities of insulin [6]. Studying biological conditions with memory characteristics that can be defined with the fractional calculus is another interesting area of use [7]. Readers interested in areas where fractional calculus is expanding should refer to the work of Sun et al. [8]. Among others, the main areas of research with fractional calculus are in physics, control, signal and image processing, mechanics and dynamic systems, biology, environmental science, materials, viscoelasticity, and economics [8]. Fractional calculus is also being used to describe fractional order chaotic systems, systems that are very sensitive to initial conditions [9]. One such interesting area for using fractional calculus is rheology, part of which is viscoelasticity. The behavior of viscoelastic materials can be simulated using fractional calculus in rheological models, where basic rheological models can be converted into fractional ones [10]. It can also be used in control theory to convert a classical PID controller into a fractional one [5,11,12]. Expressed as PI λ D µ , the fractional controller has five tuning parameters, which enables a wider fitting range for fine tuning.
Since most dielectrics used in DEA applications perform viscoelastic behavior, it is challenging to achieve satisfactory position control without overshoot. Babič et al. [13] used a basic Kelvin-Voigt model for material parametrization and self-tuning PID control. They got the first parameters for PID from the Ziegler-Nichols open-loop method. An experiment was conducted with DEA attached to a servo linear motor that maintained a constant force for pre-stretching. They were able to generate displacements up to 2 mm. Position control had 10-20% of overshoot and 0% of stationary error.
In the work of Sarban et al. [14], they used a Kelvin-Voigt model with three springs in series, where two springs had dampers in parallel connections for material parametrization of an adaptive vibrational isolator in cylindrical form. They used feedforward control, a method based on knowing the mechanical behavior of the plant. The feedforward adaptive control approach was used to predict the mechanical behavior of the isolator. They achieved active vibrational isolation but had no measured step responses.
The Gent model is also commonly used for the parametrization of viscoelastic materials, and it has one spring and a pair of spring dampers in parallel connection. It was used in the work of Li et al. [15] for DEA parametrization. They used PID control on the proposed model on a sine wave reference signal. They were able to capture amplitude and frequency, but with a phase shift of 90 • .
In the work of Ngujen et al. [16], they used the Mooney-Rivlin model for the material parametrization of cylindrically configured DEA. They used a PWM-PID control that did not eliminate overshoot and still performed small oscillation around the reference point.
In the work of Mustaza et al. [17], they used a modified Voigt model for the material parametrization of a DEA compound of three equally spaced cylinders disposed at 120 • in a symmetrical radial configuration. The modified Voigt model had been already introduced by Yeoh [18]. Displacement control was implemented with the help of model-based control. The performance of the control depended on the accuracy of the material parametrization. Experiments show that although control is possible, overshoot and static error are not eliminated.
A worm-like soft actuator was created as well as a controller in the work of Cao et al. [19]. For material parametrization they used a multi-layer Kelvin-Voigt model. They used feedforward plus feedback control schemes. They used a PI control in the feedback loop. Its parameters were gained with the help of the Nichols-Ziegler method. Simulations were tested in MATLAB/Simulink, where step responses were given without overshoot or static errors. Once tested on a real model, the performance of the proposed control algorithm showed poor performance regarding overshoot and static error.
In the work of Rizzello et al. [20,21], DEA in the shape of an annular membrane was used for position control. It was pre-stretched on the disc, so there was no need for additional load. However, a spring and an additional load were attached to the actuator. The model was derived with the help of Helmholtz free energy and a nonlinear version of a standard linear solid model. Position control was performed with PID, with the linearization of the model with PID, and with a robust H ∞ control algorithm. Although PID control gave satisfactory results, overshoot was not eliminated. On the other hand, linearized PID and H ∞ control gave excellent performance with little to no overshoot.
As can be seen, it is quite challenging to eliminate overshoot when dealing with position control on viscoelastic materials.
The use of fractional calculus in mathematical modeling and position control of the DEA is presented for the first time to improve the step response of the DEA and to eliminate the overshoot. The planar structure of DEA is used where pre-stretch is maintained with the help of appropriate weight. Since the performance of the PID control did not give satisfactory results, the control was upgraded to fractional control with the help of the PI λ D µ . The step response of the model was compared to the step response of the first order plus death time (FOPDT) function. The parameters for the PID controller were obtained with the help of FOPDT and the Wang-Juang-Chan algorithm. These parameters and λ = µ = 0.5 are starting parameters for a PI λ D µ controller. Nevertheless, the FOMCON toolbox was used in the search for the optimal parameters for the PI λ D µ controller [22]. With its help, different starting parameters were chosen to find the optimal parameters. The step responses of optimal parameters were tested in Laplace, time, and discrete domains, where responses were, surprisingly, not the same. Consequently, all optimized parameters needed to be validated with a real experiment. For the implementation of PI λ D µ to a real controller, its fractional order of integration and derivation needs to be approximated with an Oustaloup filter. The new transfer function of PI λ D µ with fractional orders approximated can be transformed into a discrete domain and implemented into the microcontroller. Implementing the standard PID Simulink block on a microcontroller provides an anti-windup method for limiting the output of the integrator. However, the possibility of integral windup is not eliminated by implementing PI λ D µ in the form of an approximated discrete transfer function. Consequently, a quasi-anti-windup method was developed and implemented on the microcontroller. Real experiments showed the advantages of using PI λ D µ over a PID controller regarding response time, settling time, and overshoot.
In the Materials and Methods section, a summary is given of the fractional calculus applied to the Kelvin-Voigt model and on parameter identification. Some basics are also presented for the use of fractional calculus in control theory. A Laplace transformation was employed to transform the fractional Kelvin-Voigt model into the s-domain. The step response was analyzed and compared with the step response of the first order plus death time (FOPDT) function. After the parameters of the FOPDT function were identified, the Wang-Juang-Chan algorithm was introduced to identify parameters of the PID controller [23]. The identified parameters represented the starting points for parameter identification of the PI λ D µ controller. The optimization of the PI λ D µ controller was explained using the FOMCON toolbox in Matlab [22]. The PI λ D µ controller needed to be transformed into the Laplace and discrete domain from the time domain. Also presented is the implementation of the PI λ D µ controller to a microcomputer. A quasi-anti-windup method was introduced to limit windup output of the approximated controller in the real experiment.
The simulation results were validated by experiments and deviations in different domain models are presented in the Results section.
The Discussion section highlights conclusions and outlines proposals for future work.

Using Fractional Calculus in Rheology and Control Theory
To apply a dielectric elastomer actuator (DEA) as a controlled soft actuator one needs to develop the governing equations of the system. The DEA uses the elastomer as the main platform on which conductive electrodes are applied. The elastomer VHB 4910 from the company 3M was used in this specific case [24]. The planar structure of the DEA was used as shown in Figure 1. The elastomer was cut to the initial non stretched dimensions of (50 × 49 × 1) mm according to the coordinate system used in Figure 1. After inserting the elastomer into the gripping jaws the dimensions between the jaws was (10 × 49 × 1) mm. Conductive paste from Bare Conductive [25] was applied to the stretched elastomer, where the dimensions changed to (60 × 49 × 0.16) mm. Since high voltage can spark over the edge, creating a short circuit, the active surface of the DEA was reduced to (60 × 40) mm. Figure 2 shows the planar configuration of the stretched elastomer and active surface of the DEA, as well as the applied coordinate system. When DC voltage is applied, the elastomer expands in the longitudinal or x direction. However, small movement is also seen in the y direction, which represents the loss of energy in the system. Actuators 2021, 10, x FOR PEER REVIEW 5 of 18

Fractional Calculus in Rheology
Dynamic analysis of the used elastomer revealed viscoelastic behavior, which means that the elastomer performed with both elastic and viscous behavior. According to [10,26,27], one can model material behavior with fractional calculus and rheology. Using the Kelvin-Voigt model and fractional calculus, one can modify the basic Kelvin-Voigt model to a fractional Kelvin-Voigt model, which consists of two parts: a spring and a springpot. The springpot is an element that is a combination of a spring and a dashpot. Its characteristic is determined with the springpot constant and the derivative order which is its main advantage over the basic Kelvin-Voigt model [5,10,27]. In Equation (1) the Grünwald-Letnikov definition of fractional derivative is used ( The Grünwald-Letnikov definition is chosen since it can easily handle fractional derivation (p) and integration (−p), and it is also suitable for the numerical calculation of fractional derivatives and does not lead to difficulties when applying it to real-world problems [5,28].
If the derivative order of a springpot equals 1, the element represents a dashpot; if the derivative order equals 0, the element represents a spring. With fractional calculus the derivative order is no longer limited to integers and can be any real number between 0 and 1.

Fractional Calculus in Control Theory
To show that the fractional Kelvin-Voigt model is stable it must be transformed into the Laplace domain and written in the form of commensurate order.
A Laplace transform in fractional calculus transforms fractional derivation into fractional multiplication and fractional integration into fractional division. The following equation, represents the Laplace transform of a fractional derivative with zero initial conditions.
The following equation,

Fractional Calculus in Rheology
Dynamic analysis of the used elastomer revealed viscoelastic behavior, which means that the elastomer performed with both elastic and viscous behavior. According to [10,26,27], one can model material behavior with fractional calculus and rheology. Using the Kelvin-Voigt model and fractional calculus, one can modify the basic Kelvin-Voigt model to a fractional Kelvin-Voigt model, which consists of two parts: a spring and a springpot. The springpot is an element that is a combination of a spring and a dashpot. Its characteristic is determined with the springpot constant and the derivative order which is its main advantage over the basic Kelvin-Voigt model [5,10,27]. In Equation (1) the Grünwald-Letnikov definition of fractional derivative is used The Grünwald-Letnikov definition is chosen since it can easily handle fractional derivation (p) and integration (−p), and it is also suitable for the numerical calculation of fractional derivatives and does not lead to difficulties when applying it to real-world problems [5,28].
If the derivative order of a springpot equals 1, the element represents a dashpot; if the derivative order equals 0, the element represents a spring. With fractional calculus the derivative order is no longer limited to integers and can be any real number between 0 and 1.

Fractional Calculus in Control Theory
To show that the fractional Kelvin-Voigt model is stable it must be transformed into the Laplace domain and written in the form of commensurate order.
A Laplace transform in fractional calculus transforms fractional derivation into fractional multiplication and fractional integration into fractional division. The following equation, represents the Laplace transform of a fractional derivative with zero initial conditions. The following equation, represents the transfer function of the fractional Kelvin-Voigt model. The dynamic characteristics of the elastomer VHB 4910 were identified in our previous work [29,30], where the parameters were described in the relationship of δ − ε and can be used for any specific dimensions. The constants for the fractional Kelvin-Voigt model were calculated for the specified DEA according to and In Equations (5) and (6), E = 1.06106 N/mm 2 , η = 0.086904 Ns/mm 2 , and L = 60 mm and is the extended length in the x direction, and A = 49 × 0.16 mm 2 and is the cross-section of the elastomer.
After recalculation, the parameters were k = 0.138645 N/mm, c = 0.086904 Ns/mm, and α = 0.47, where α is the same as in the δ − ε relationship. To transform it into commensurate order parameter α needs to be written in the form of α = 47/100 and α = β × 47, where β = 1/100 [31]. If complex variable s β = σ, one can form the commensurate order transfer function of the fractional Kelvin-Voigt model: Since the order of the transfer function changed, the stability criterion also changes from As can be seen in Figure 3, all poles of the commensurate order transfer function of the fractional Kelvin-Voigt model satisfied the stability criterion in Equation (9), which shows that the fractional order Kelvin-Voigt model is mathematically stable.
represents the transfer function of the fractional Kelvin-Voigt model. The dynamic characteristics of the elastomer VHB 4910 were identified in our previous work [29,30], where the parameters were described in the relationship of − and can be used for any specific dimensions. The constants for the fractional Kelvin-Voigt model were calculated for the specified DEA according to In Equations (5) and (6), and is the extended length in the x direction, and A = 49 × 0.16 2 and is the cross-section of the elastomer.
After recalculation, the parameters were = 0.138645 / , = 0.086904 / , and = 0.47, where is the same as in the − relationship. To transform it into commensurate order parameter needs to be written in the form of = 47/100 and = × 47, where = 1/100 [31]. If complex variable = , one can form the commensurate order transfer function of the fractional Kelvin-Voigt model: Since the order of the transfer function changed, the stability criterion also changes from As can be seen in Figure 3, all poles of the commensurate order transfer function of the fractional Kelvin-Voigt model satisfied the stability criterion in Equation (9), which shows that the fractional order Kelvin-Voigt model is mathematically stable.
where m = 0.0001 kg and represents the mass of the elastomer.
Using the Laplace transform on Equation (10) one obtains the following equation: while using the fractional derivative plant model, a logical conclusion was to use fractional PID control, which is written in the form of PI λ D µ , since more adjustable parameters are available for the fine tuning of the controller. The fractional PI λ D µ has the same parameters for gains as the basic PID, except that it has two additional parameters: derivation and integration order. Since fractional calculus is used, orders are limited to real numbers between 0 and 1. It should be noted that derivation and integration are reciprocal operations.
According to the Grünwald-Letnikov definition of fractional calculus, positive order describes derivation, whereas negative order describes integration. In the case of the PI λ D µ , in the control theory it is commonly known that the integral part should be negative, even though, in most literature, it is positive. When using fractional calculus in the plant model as well as for PI λ D µ control, one can expect different slopes in the Bode diagram. Fractional derivatives impose a phase shift of −απ/2, where α is the derivation order. It also changes the slope to −20α dB/dec. Figure 4 shows Bode diagrams for the fractional (k = 0.138645 N/mm, c = 0.086904 Ns/mm, and α = 0.47) and the basic Kelvin-Voigt model (k = 0.138645 N/mm and c = 0.086904 Ns/mm). It can be seen how the slope of the fractional Kelvin-Voigt changes. To create Bode plots it is recommended to use the FOTF Matlab Toolbox [32] or the FOMCON Toolbox [22].
The equation of motion expressed with the fractional Kelvin-Voigt model is given as where = 0.0001 and represents the mass of the elastomer. Using the Laplace transform on Equation (10) (11) while using the fractional derivative plant model, a logical conclusion was to use fractional PID control, which is written in the form of λ , since more adjustable parameters are available for the fine tuning of the controller. The fractional λ has the same parameters for gains as the basic PID, except that it has two additional parameters: derivation and integration order. Since fractional calculus is used, orders are limited to real numbers between 0 and 1. It should be noted that derivation and integration are reciprocal operations. According to the Grünwald-Letnikov definition of fractional calculus, positive order describes derivation, whereas negative order describes integration. In the case of the λ , in the control theory it is commonly known that the integral part should be negative, even though, in most literature, it is positive.
When using fractional calculus in the plant model as well as for λ control, one can expect different slopes in the Bode diagram. Fractional derivatives impose a phase shift of − /2, where α is the derivation order. It also changes the slope to −20 / . Figure 4 shows Bode diagrams for the fractional ( = 0.138645 / , = 0.086904 / , and = 0.47) and the basic Kelvin-Voigt model ( = 0.138645 / and = 0.086904 / ). It can be seen how the slope of the fractional Kelvin-Voigt changes. To create Bode plots it is recommended to use the FOTF Matlab Toolbox [32] or the FOMCON Toolbox [22].    Comparing the step responses of both functions, one can obtain the parameters of FOPDT as k = 5.6652, L = 0.0055 s, and T = 0.7891 s. Using the Wang-Juang-Chan algorithm [23] for PID parameters identification, the following parameters were obtained: k p = 7.179, k i = 5.665, and k d = 0.03107. After the parameters of the PID controller were obtained, the parameters of the PI λ D µ gains were taken from the PID, and derivative and integral orders were chosen to be both at 0.5 [33,34]. Figure 5 shows the step responses of the plant model with the PID and the PI λ D µ controller in the Laplace domain. Both controllers gave fast responses. The PID controller produced no steady-state error, whereas the PI λ D µ produced a small steady-state error of about 0.5% of the desired value, which makes sense since only half of the integrator was used.

 
process time constant s T  Comparing the step responses of both functions, one can obtain the parameters of FOPDT as = 5.6652, = 0.0055 , and = 0.7891 . Using the Wang-Juang-Chan algorithm [23] for PID parameters identification, the following parameters were obtained: = 7.179, = 5.665, and = 0.03107. After the parameters of the PID controller were obtained, the parameters of the λ gains were taken from the PID, and derivative and integral orders were chosen to be both at 0.5 [33,34]. Figure 5 shows the step responses of the plant model with the PID and the λ controller in the Laplace domain. Both controllers gave fast responses. The PID controller produced no steady-state error, whereas the λ produced a small steady-state error of about 0.5% of the desired value, which makes sense since only half of the integrator was used.

Optimization of Parameters for the Controller by the Use of the FOMCON Toolbox in Matlab and Implementation of the Controller on a Microprocessor
Optimization of Parameters for the λ Controller The initial parameters for gains of the λ controller were calculated with the Wang-Juang-Chan algorithm, whereas derivation and integration orders were chosen freely. Optimization of λ parameters can be performed according to the work of Tepljakov [22,33], who developed the FOMCON toolbox, which stands for Fractional Order Modeling and Control. Figure 5 shows the FOMCON toolbox for the λ parameter optimization. According to Tepljakov, the best way to find optimized parameters is to fix the gains of λ and to vary the integration and derivation orders [33]. The plant model, which is defined as a fractional order transfer function, can be optimized according to the gain and phase margin of the open-loop control model. Different optimization algorithms are available, as well as different performance metrics. In this specific case, the Nelder-Mead optimization method was used, since it is a direct search method and is therefore well suited to optimizing a function with unknown derivatives [33,35]. A gain margin of 10 dB and a phase margin of 45° were selected as the optimal responses of the closed-loop plant model. Different initial parameters were taken to be optimized with the FOMCON. It can optimize only gains, only orders, or all parameters at once.
The methodology of the optimization is such that fractional transfer function of the plant is first approximated with the Oustaloup filter with the selected order of approximation and frequency range. Open-loop frequency response is then performed with an approximated transfer function where the gain and phase margins are calculated. The Step response of plant model with the PID and PI λ D µ controller.

Optimization of Parameters for the PI λ D µ Controller by the Use of the FOMCON Toolbox in Matlab and Implementation of the PI λ D µ Controller on a Microprocessor
Optimization of Parameters for the PI λ D µ Controller The initial parameters for gains of the PI λ D µ controller were calculated with the Wang-Juang-Chan algorithm, whereas derivation and integration orders were chosen freely. Optimization of PI λ D µ parameters can be performed according to the work of Tepljakov [22,33], who developed the FOMCON toolbox, which stands for Fractional Order Modeling and Control. Figure 5 shows the FOMCON toolbox for the PI λ D µ parameter optimization. According to Tepljakov, the best way to find optimized parameters is to fix the gains of PI λ D µ and to vary the integration and derivation orders [33]. The plant model, which is defined as a fractional order transfer function, can be optimized according to the gain and phase margin of the open-loop control model. Different optimization algorithms are available, as well as different performance metrics. In this specific case, the Nelder-Mead optimization method was used, since it is a direct search method and is therefore well suited to optimizing a function with unknown derivatives [33,35]. A gain margin of 10 dB and a phase margin of 45 • were selected as the optimal responses of the closed-loop plant model. Different initial parameters were taken to be optimized with the FOMCON. It can optimize only gains, only orders, or all parameters at once.
The methodology of the optimization is such that fractional transfer function of the plant is first approximated with the Oustaloup filter with the selected order of approximation and frequency range. Open-loop frequency response is then performed with an approximated transfer function where the gain and phase margins are calculated. The Nelder-Mead optimization method was selected to optimize gains and order parameters. It varies gains and order parameters inside constraints for each parameter, where gains were constrained between 0 and 100, and orders were constrained between 0 and 1. The integral square error (ISE) method was selected to evaluate the error between the desired and calculated phase and gain margins. The Nelder-Mead method calculates the centroid of all points regarding the optimized parameters. If the calculated error ISE is greater than in the previous step, then the method contracts points for finding the optimum. If the calculated error ISE is smaller than in the previous step, then the method expands points for finding the optimum [33]. The first optimization was performed and tested with a step function in the Laplace domain with initial parameters k p = 1, k i = 1, k d = 1, λ = 0.5, and µ = 0.5. Response time T r [s], settling time T s [s], and overshoot P [%] were measured and observed for both controllers. The fastest response time, settling time, and no overshoot were demanded. The second optimization was performed with initial parameters that were calculated with the help of the Wang-Juang-Chan algorithm, namely, k p = 5.665, k i = 7.179, k d = 0.031, λ = 0.5, and µ = 0.5. Once again, the best response was obtained when optimizing all parameters at once. The observed differences in the step responses of the controlled plant optimized with different initial parameters were small, and may be neglected in the Laplace domain.

Continuous and Discrete Simulink Models of DEA
To obtain a functional controlled DEA actuator, the optimized parameters need to be tested in Simulink, since the real DEA is controlled with a Simulink model uploaded to the DSP controller.
PI λ D µ in the Laplace domain must be approximated with a filter for implementation into Simulink. The Oustaloup filter was used for the PI λ D µ approximation [23,36]. Its mathematical expression is represented as It must be pointed out that the Oustaloup filter only approximates the integration and derivation orders, rather than the complete transfer function. The final transfer function of the approximated PI λ D µ controller is then computed with individual approximation. The combined order of the approximated transfer function unfortunately becomes very high. The order of transfer function depends on the N order of approximation. The higher the order of approximation, the more accurate it is.
A continuous Simulink plant model of the DEA with the PI λ D µ approximation was built with the use of the FOTF toolbox. To implement the PI λ D µ control to the DSP microcontroller, a continuous time approximation with an Oustaloup filter must be converted into a discrete model. This is achieved with the Matlab function c2d and by defining the sampling time.
The discrete Simulink model was built with the FOMCON toolbox, as it includes blocks for discrete PI λ D µ controller implementation. The output of the controller is the force with which it is applied on the DEA. In the continuous and in the discrete Simulink model the controller output is saturated between 0 and 1.5. Negative values of the controller output must not have any impact on the DEA model in the Simulink simulation, since only internal connections of polymers can retract elastomer. The output represents a generated force of 1.5 N, which is the maximum force that can be generated on the real DEA model with the applied voltage source. The maximum force was calculated with the help of the Maxwell stress equation, where stress was converted into force with the specified dimensions of the elastomer in use [3]. The Maxwell stress generated on the DEA when voltage is applied is expressed as Force can be calculated from multiplying Equation (11) with the active surface as In Equation (13), ε r = 4.7 is relative permittivity according to [37], and l 1 , l 2 and z are the dimensions of the active surface and the thickness of the DEA, respectively.

Implementation of PI λ D µ on the DSP Microcontroller
In the discrete Simulink model, the discrete PI λ D µ block from the FOMCON toolbox was used, which performs an Oustaloup approximation of individual parameters automatically. However, as this block does not enable resetting the PI λ D µ controller outputs, an integral windup problem can occur. In order to implement the PI λ D µ controller into the DSP microcontroller TMS320F28069M from Texas Instruments [38], the PI λ D µ needed to be approximated with an Oustaloup filter, written as a transfer function and then translated into a discrete transfer function with the ability to reset it. This is called the quasi-antiwindup method. Figures 6 and 7 show the Simulink model for DSP implementation with serial communication to Matlab/Simulink, and real time monitoring. Figure 8 shows the structure of the ADC_ISR block, where ADC and control are performed. Figure 9 shows the approximated transfer function of fractional control updated with the quasi-anti-windup method implemented on DSP.
DEA model with the applied voltage source. The maximum force was calculated with the help of the Maxwell stress equation, where stress was converted into force with the specified dimensions of the elastomer in use [3]. The Maxwell stress generated on the DEA when voltage is applied is expressed as (14) Force can be calculated from multiplying Equation (11) with the active surface as   (15) In Equation (13), = 4.7 is relative permittivity according to [37], and 1 , 2 and are the dimensions of the active surface and the thickness of the DEA, respectively.

Implementation of on the DSP Microcontroller
In the discrete Simulink model, the discrete λ block from the FOMCON toolbox was used, which performs an Oustaloup approximation of individual parameters automatically. However, as this block does not enable resetting the λ controller outputs, an integral windup problem can occur. In order to implement the λ controller into the DSP microcontroller TMS320F28069M from Texas Instruments [38], the λ needed to be approximated with an Oustaloup filter, written as a transfer function and then translated into a discrete transfer function with the ability to reset it. This is called the quasi-anti-windup method. Figure 6 and Figure 7 show the Simulink model for DSP implementation with serial communication to Matlab/Simulink, and real time monitoring. Figure 8 shows the structure of the ADC_ISR block, where ADC and control are performed. Figure 9 shows the approximated transfer function of fractional control updated with the quasi-anti-windup method implemented on DSP.      The quasi-anti-windup method resets the discrete transfer function when the output of the approximated discrete λ controller exceeds two times the maximum output from saturation, which represents a force of 4 N. The method proved to be very successful, which becomes more evident in the Results section, where the experimental results are shown.

Results
The step responses of optimized λ parameters on average had the same rise time, the settling time was approximately 5.5 times faster, and overshoot was nearly eliminated in comparison to its PID counterpart, as can be seen in Table 1. However, the Laplace domain can easily handle fractional orders since derivative is translated into multiplication and integration into division. Even more, the output of the regulator in the Laplace domain is not limited, which means its output can be infinite. Configuring the model    The quasi-anti-windup method resets the discrete transfer function when the output of the approximated discrete λ controller exceeds two times the maximum output from saturation, which represents a force of 4 N. The method proved to be very successful, which becomes more evident in the Results section, where the experimental results are shown.

Results
The step responses of optimized λ parameters on average had the same rise time, the settling time was approximately 5.5 times faster, and overshoot was nearly eliminated in comparison to its PID counterpart, as can be seen in Table 1. However, the Laplace domain can easily handle fractional orders since derivative is translated into multiplication and integration into division. Even more, the output of the regulator in the Laplace domain is not limited, which means its output can be infinite. Configuring the model The quasi-anti-windup method resets the discrete transfer function when the output of the approximated discrete PI λ D µ controller exceeds two times the maximum output from saturation, which represents a force of 4 N. The method proved to be very successful, which becomes more evident in the Results section, where the experimental results are shown.

Results
The step responses of optimized PI λ D µ parameters on average had the same rise time, the settling time was approximately 5.5 times faster, and overshoot was nearly eliminated in comparison to its PID counterpart, as can be seen in Table 1. However, the Laplace domain can easily handle fractional orders since derivative is translated into multiplication and integration into division. Even more, the output of the regulator in the Laplace domain is not limited, which means its output can be infinite. Configuring the model in Simulink where the output of the controller must be limited is vital to testing optimized parameters before performing the real experiment. The most promising parameters for the real experiment are highlighted with yellow and the best parameters are highlighted with green in Table 1.
The step responses of the continuous and the discrete Simulink model were compared with the step responses of the Laplace domain model for all optimized parameters. Table 2 shows the results of the step responses of the same controller plant model built and tested in different domains with optimized parameters. Table 1. Initial, optimal, and comparing parameters of the PI λ D µ and PID and their step responses in the Laplace domain.

K p K i K d λ µ T r [s] T s [s] P [%]
PI λ D µ Initial param. The best results were obviously obtained in the Laplace domain. There are more explanations for this. Fractional derivatives and integrals can be transformed easily into the Laplace domain. This could be understood as having an exact plant model with a PI λ D µ controller. The second reason is that the output of the PI λ D µ controller in the Laplace domain is not saturated between 0 and 1.5 N to limit its output to the maximum force that could be applied on the DEA. Not only is its output not limited, but it also enables the output to be negative. Negative output from the controller means negative force acting on the DEA model, which reduced the settling time and the overshoot. On the other hand, the continuous Simulink model uses a high order Oustaloup filter for PI λ D µ approximation since direct use is not possible in Simulink. The output of the approximated PI λ D µ controller was saturated between 0 and 1.5, which was responsible for the longer settling time and the overshoot of the DEA model. For the simulation of the continuous Simulink model, a variable simulation step was chosen. The simulation step in the discrete model was set to 0.1 µs, which was the same as the step used in the microprocessor. When a grater step was chosen, discrete simulation in Simulink did not converge successfully or it had large settling time. If a small enough step size is chosen for the discrete simulation then the differences between the continuous and the discrete Simulink models are significantly lowered. Otherwise, the differences can be obvious.
Optimal parameters for the model built in the Laplace domain were also optimal for the continuous and discrete Simulink model. The best parameters for all domains are highlighted in green, and the worst parameters used for PID control in all domains are highlighted in red in Table 2. After the DEA was set up for the experimental part it was then mounted on the platform in a vertical position and loaded with 315 g weight. According to previously conducted experiments, this is the appropriate weight to hold the DEA in the desired position, otherwise its internal polymers' connections would bring it back to the initial unloaded state [29]. Laser position sensor Wenglor YP05MGV80 [39] was mounted on the DEA. A laser was used since it offers contactless measurements. A high voltage supply was generated from DC to the HVDC converter E1-P600/Y [40]. Figure 10 shows the experimental setup for the real DEA.
Actuators 2021, 10, x FOR PEER REVIEW 13 of 18 highlighted in green, and the worst parameters used for PID control in all domains are highlighted in red in Table 2. After the DEA was set up for the experimental part it was then mounted on the platform in a vertical position and loaded with 315 g weight. According to previously conducted experiments, this is the appropriate weight to hold the DEA in the desired position, otherwise its internal polymers' connections would bring it back to the initial unloaded state [29]. Laser position sensor Wenglor YP05MGV80 [39] was mounted on the DEA. A laser was used since it offers contactless measurements. A high voltage supply was generated from DC to the HVDC converter E1-P600/Y [40]. Figure 10 shows the experimental setup for the real DEA.

Step Responses of the Real DEA with PID and Control
The step responses of the DEA were separated into step responses with PID control and step responses with λ control. Figure 11 shows the step responses of individual sets of PID parameters. Eighteen experiments were performed, where different parameters were used. Promising parameters were experimented multiple times. Response time and settling time were measured as well as overshoot (P[%]) and steady-state error. Table  3 shows the results of the step responses of the DEA with PID control. It can be seen from Figure 11 that not all successful step responses performed in the Laplace domain and in

Step Responses of the Real DEA with PID and PI λ D µ Control
The step responses of the DEA were separated into step responses with PID control and step responses with PI λ D µ control. Figure 11 shows the step responses of individual sets of PID parameters. Eighteen experiments were performed, where different parameters were used. Promising parameters were experimented multiple times. Response time and settling time were measured as well as overshoot (P[%]) and steady-state error. Table 3 shows the results of the step responses of the DEA with PID control. It can be seen from Figure 11 that not all successful step responses performed in the Laplace domain and in the Simulink model gave good results with the real DEA. The best results were obtained with the experiments of numbers 4, 6, 8, 10, and 17. The average rise time was 1.9 s. The average settling time was 10.4 s. The average overshoot was 8.1% of the desired value. Some experiments also produced a steady-state error. The best performers in Table 3 are marked with the same color as in Figure 11a. Some experiments also produced a steady-state error. The best performers in Table 3 are marked with the same color as in Figure 11a.
Step responses of a real DEA using (a) a PID control and (b) a λ control. Table 3.
Step responses of a real DEA using a PID control with different parameters.

Exp. Number
[  Figure 11b shows the step responses of the λ control of the real DEA. An integral windup problem can be seen at the first four experiment numbers. In some experiments (4 to 9) the integral windup problem did not occur, but this did not guarantee the stability of the controller. Experiments with numbers 10 to 15 show step responses with the quasi-anti-windup method applied. The results were greatly improved, as rise times and settling times were reduced, and overshoots were reduced to an average 2.6% of the desired value. Table 4 shows the parameters used for individual experiments, as well as the results. The best results are highlighted with the same color as in Figure 11b. Step responses of a real DEA using (a) a PID control and (b) a PI λ D µ control. Table 3.
Step responses of a real DEA using a PID control with different parameters.  Figure 11b shows the step responses of the PI λ D µ control of the real DEA. An integral windup problem can be seen at the first four experiment numbers. In some experiments (4 to 9) the integral windup problem did not occur, but this did not guarantee the stability of the controller. Experiments with numbers 10 to 15 show step responses with the quasianti-windup method applied. The results were greatly improved, as rise times and settling times were reduced, and overshoots were reduced to an average 2.6% of the desired value. Table 4 shows the parameters used for individual experiments, as well as the results. The best results are highlighted with the same color as in Figure 11b.  Figure 12 shows the comparison of the step responses of both controllers for a real DEA. PID controller results are marked with experiment numbers 0 to 4, and PI λ D µ controller results are marked with experiment numbers 5 to 9. The best parameters for PID control were obtained by the Wuan-Juang-Chan algorithm. Average rise time among step responses with the PID controller was 1.7 s, the average settling time was 10.4 s, and the average overshoot was 8% of the desired value. On the other hand, these gain parameters were also the best for the PI λ D µ control where integration and derivation orders were optimized with help of the FOMCON toolbox. Table 5 shows the comparison of the results of the step responses obtained with both controllers. It can be seen that the PI λ D µ control improved step responses compared to the PID. For the PI λ D µ control the average rise time was reduced to 1 s, settling time to 5.4 s, and overshoot to 4.3%. With the implementation of the PI λ D µ control a small steady-state error was introduced, which was not greater than 1% and did not affect the controller's performance significantly. When the order of integration in the PI λ D µ controller was less than 1, steady-state error was introduced to the step response of the plant. Figure 12b shows the comparison of step responses with the same parameters of PID and PI λ D µ control in different domains. These parameters were K p = 5.66, K i = 7.17, K d =0.03, λ = 0.6, and µ = 0.8. The best results were obtained for experiments 0 and 1, representing the Laplace domain, where the rise and settling times were the smallest and there was no overshoot. The reason is that the output of the controller was not saturated between 0 and 1.5 as it was in the Simulink models. The discrete Simulink models had slightly greater response and settling time as well as overshoot, as seen for experiments 2 and 3. Experiments 4 and 5 show the experimental results of the PI λ D µ and the PID control on the real DEA, respectively. Table 6 shows the response time, settling time, overshoot, and steady-state error of the step responses of controllers in different domains.

Comparing Results of Step Responses for Both Controllers and the Best Parameters in All Domains
improved step responses compared to the PID. For the λ control the average rise time was reduced to 1 s, settling time to 5.4 s, and overshoot to 4.3%. With the implementation of the λ control a small steady-state error was introduced, which was not greater than 1% and did not affect the controller's performance significantly. When the order of integration in the λ controller was less than 1, steady-state error was introduced to the step response of the plant.   Figure 12. Comparison of the step responses of (a) the PID (0-4) and PI λ D µ (5-8) control of a real DEA and comparison of step responses of (b) the PID and PI λ D µ control in different domains. Table 5. Comparison of the step responses of the PID (0-4) and PI λ D µ (5-8) control of a real DEA.

Exp.
between 0 and 1.5, which represent the force in N. The introduction of the saturated output of the controller had some negative effects on the step response. It reduced the response time and the settling time and increased the overshoot in the Simulink simulations.
Approximation of the PI λ D µ controller with high order Oustaloup filter is needed to perform a Simulink simulation. Continuous simulation was performed with a variable step solver in Simulink. For discrete simulation, a fixed step was chosen, corresponding to the microprocessor step. The size of the fixed step for discrete Simulink simulation must be properly chosen. Greater step size caused larger response time, settling time, and overshoots in discrete Simulink simulation.
To implement the PI λ D µ to the microprocessor it must be approximated with an Oustaloup filter in transfer function form and transformed into the z-domain with the same step as defined in the microprocessor. The first position control experiments of the DEA actuator resulted in large overshoot with a large steady-state error. This implied an integral windup problem. A quasi-anti-windup method was introduced to reset the transfer function of the PI λ D µ when the output exceeded a threshold value. As shown by the results, successful implementation and realization of the PI λ D µ position control was achieved for the selected DEA. Rise time was reduced by nearly half to around 1 s, as was settling time, which was reduced to around 5 s. Overshoot was also much lower, averaging around 4% of the final value. The main disadvantage of the PI λ D µ controller was the appearance of a steady-state error, even though the error was smaller than 1% of the desired value and did not influence the performance of the controller.
Observing the transition from the Laplace to the discrete domain, one can conclude that overshoot in step response in the Laplace domain predicts even greater overshoot in step response in continuous and discrete domains for the DEA plant. Controller output in the Laplace domain is not limited as it is in continuous and discrete domains between 0 and 1.5 and can also have negative values, which reduces the overshoot. If overshoot cannot be eliminated with unlimited output of the controller in the Laplace domain, then one can expect overshoot in continuous and discrete domains.
In future work, an anti-windup method with a transfer function for PI λ D µ approximation is going to be realized, separated into individual blocks of the z-transform. Each block will be limited with the output. Hopefully this will bring in even better results for position control.
Funding: Authors acknowledges the financial support from Ministry of Education, Science and Sport and European Regional Development Fund Investing in Your Future. Part of this research was done under ROBKONCEL OP20.03530 project.