# Some Applications of Fractional Calculus in Engineering

Authors: J. A. Tenreiro Machado; Manuel F. Silva; Ramiro S. Barbosa; Isabel S. Jesus; Cecília M. Reis; Maria G. Marcos; Alexandra F. Galhano
Journal: Mathematical Problems in Engineering
 (2010)

---

## Abstract

Fractional Calculus
(FC) goes back to the beginning of the theory of
differential calculus. Nevertheless, the
application of FC just emerged in the last two
decades, due to the progress in the area of
chaos that revealed subtle relationships with
the FC concepts. In the field of dynamical
systems theory some work has been carried out
but the proposed models and algorithms are still
in a preliminary stage of establishment. Having
these ideas in mind, the paper discusses FC in
the study of system dynamics and control. In
this perspective, this paper investigates the
use of FC in the fields of controller tuning,
legged robots, redundant robots, heat diffusion,
and digital circuit synthesis.


---

## Body

## 1. Introduction

The generalization of the concept of derivativeDα[f(x)] to noninteger values of α goes back to the beginning of the theory of differential calculus. In fact, Leibniz, in his correspondence with Bernoulli, L’Hôpital and Wallis (1695), had several notes about the calculation of D1/2[f(x)]. Nevertheless, the development of the theory of Fractional Calculus (FC) is due to the contributions of many mathematicians such as Euler, Liouville, Riemann, and Letnikov [1–3].The FC deals with derivatives and integrals to an arbitrary order (real or, even, complex order). The mathematical definition of a derivative/integral of fractional order has been the subject of several different approaches [1–3]. For example, the Laplace definition of a fractional derivative of a signal x(t) is
(1.1)Dαx(t)=L-1{sαX(s)-∑k=0n-1skDα-k-1x(t)∣t=0},
where n-1<α≤n, α>0. The Grünwald-Letnikov definition is given by (α∈ℜ):
(1.2)Dαx(t)=limh→0[1hα∑k=0∞(-1)k(αk)x(t-kh)],(αk)=Γ(α+1)Γ(k+1)Γ(α-k+1),
where Γ is the Gamma function and h is the time increment. However, (1.2) shows that fractional-order operators are “global’’ operators having a memory of all past events, making them adequate for modeling memory effects in most materials and systems.The Riemann-Liouville definition of the fractional-order derivative is (α>0):
(1.3)aDtαf(t)=1Γ(n-α)dndtn∫atf(τ)(t-τ)α-n+1dτ,n-1<α<n,
where Γ(x) is the Gamma function of x.Based on the proposed definitions it is possible to calculate the fractional-order integrals/derivatives of several functions (Table1). Nevertheless, the problem of devising and implementing fractional-order algorithms is not trivial and will be the matter of the following sections.Table 1
Fractional-order integrals of several functions.
φ(x),x∈ℜ(I+αφ)(x),x∈ℜ,α∈C(x-a)β-1Γ(β)Γ(α+β)(x-a)α+β-1,Re(β)>0eλxλ-αeλx,Re(λ)>0{sin(λx)cos(λx)λ-α{sin(λx-απ/2),cos(λx-απ/2),λ>0,Re(α)>1eλx{sin(γx)cos(γx)eλx(λ2+γ2)α/2{sin(γx-αϕ),ϕ=arctan(γ/λ)cos(γx-αϕ),γ>0,Re(λ)>1In recent years FC has been a fruitful field of research in science and engineering [1–6]. In fact, many scientific areas are currently paying attention to the FC concepts and we can refer its adoption in viscoelasticity and damping, diffusion and wave propagation, electromagnetism, chaos and fractals, heat transfer, biology, electronics, signal processing, robotics, system identification, traffic systems, genetic algorithms, percolation, modeling and identification, telecommunications, chemistry, irreversibility, physics, control systems as well as economy, and finance [7–18].Bearing these ideas in mind, Sections2–6 present several applications of FC in science and engineering. In Section 2, it is presented the application of FC concepts to the tuning of PID controllers and, in Section 3, the application of a fractional-order PD controller in the control of the leg joints of a hexapod robot. Then in Section 4, it is presented the fractional dynamics in the trajectory control of redundant manipulators. Next, in Section 5, it is introduced the fractional characteristics of heat diffusion along a media and, in Section 6 it is shown the application of FC to circuit synthesis using evolutionary algorithms. Finally, the main conclusions are presented in Section 7.


## 2. Tuning of PID Controllers Using Fractional Calculus Concepts

The PID controllers are the most commonly used control algorithms in industry. Among the various existent schemes for tuning PID controllers, the Ziegler-Nichols (Z-N) method is the most popular and is still extensively used for the determination of the PID parameters. It is well known that the compensated systems, with controllers tuned by this method, have generally a step response with a high percent overshoot. Moreover, the Z-N heuristics are only suitable for plants with monotonic step response.In this section, we study a methodology for tuning PID controllers such that the response of the compensated system has an almost constant overshoot defined by a prescribed value. The proposed method is based on the minimization of the integral of square error (ISE) between the step responses of a unit feedback control system, whose open-loop transfer functionL(s) is given by a fractional-order integrator and that of the PID compensated system [7].Figure1 illustrates the fractional-order control system that will be used as reference model for the tuning of PID controllers. The open-loop transfer function L(s) is defined as (α∈ℜ+):
(2.1)L(s)=(ωcs)α,
where ωc is the gain crossover frequency, that is, |L(jωc)|=1. The parameter α is the slope of the magnitude curve, on a log-log scale, and may assume integer as well as noninteger values. In this study we consider 1<α<2, such that the output response may have a fractional oscillation (similar to an underdamped second-order system). This transfer function is also known as the Bode’s ideal loop transfer function since Bode studies on the design of feedback amplifiers in the 1940s [19].Figure 1
Fractional-order control system with open-loop transfer functionL(s).The Bode diagrams of amplitude and phase ofL(s) are illustrated in Figure 2. The amplitude curve is a straight line of constant slope -20α dB/dec, and the phase curve is a horizontal line positioned at -απ/2 rad. The Nyquist curve is simply the straight line through the origin, arg L(jω)=-απ/2 rad.Figure 2
Bode diagrams of amplitude and phase ofL(jω) for 1<α<2.This choice ofL(s) gives a closed-loop system with the desirable property of being insensitive to gain changes. If the gain changes, the crossover frequency ωc will change, but the phase margin of the system remains PM=π(1-α/2) rad, independent of the value of the gain. This can be seen from the curves of amplitude and phase of Figure 2.The closed-loop transfer function of fractional-order control system of Figure1 is given by
(2.2)G(s)=L(s)1+L(s)=1(s/ωc)α+1,1<α<2.The unit step response ofG(s) is given by the expression:
(2.3)yd(t)=L-1{1sG(s)}=L-1{ωcαs(sα+ωcα)}=1-∑n=0∞[-(ωct)α]nΓ(1+αn)=1-Eα[-(ωct)α].For the tuning of PID controllers, we address the fractional-order transfer function (2.2) as the reference system [8]. With the order α and the crossover frequency ωc we can establish the overshoot and the speed of the output response, respectively. For that purpose we consider the closed-loop system shown in Figure 3, where Gc(s) and Gp(s) are the PID controller and the plant transfer functions, respectively.Figure 3
Closed-loop control system with PID controllerGc(s).The transfer function of the PID controller is(2.4)Gc(s)=U(s)E(s)=K(1+1Tis+Tds),
where E(s) is the error signal and U(s) is the controller’s output. The parameters K, Ti, and Td are the proportional gain, the integral time constant, and the derivative time constant of the controller, respectively.The design of the PID controller will consist on the determination of the optimum PID set gains(K,Ti,Td) that minimize J, the integral of the square error (ISE), defined as
(2.5)J=∫0∞[y(t)-yd(t)]2dt,
where y(t) is the step response of the closed-loop system with the PID controller (Figure 3) and yd(t) is the desired step response of the fractional-order transfer function (2.2) given by (2.3).To illustrate the effectiveness of proposed methodology we consider the third-order plant transfer function:(2.6)Gp(s)=Kp(s+1)3
with nominal gain Kp=1.Figure4 shows the step responses and the Bode diagrams of phase of the closed-loop system with the PID for the transfer function Gp(s) for gain variations around the nominal gain (Kp=1) corresponding to Kp={0.6,0.8,1.0,1.2,1.4}, that is, for a variation up to ±40% of its nominal value. The system was tuned for α=3/2 (PM=45°), ωc=0.8 rad/s. We verify that we get the same desired iso-damping property corresponding to the prescribed (α,ωc) values.Bode phase diagrams and step responses for the closed-loop system with a PID controller forGp(s). The PID parameters are K=1.9158, Ti=1.1407, and Td=0.9040.
(a)(b)In fact, we observe that the step responses have an almost constant overshoot independent of the variation of the plant gain around the gain crossover frequencyωc. Therefore, the proposed methodology is capable of producing closed-loop systems robust to gain variations and step responses exhibiting an iso-damping property. The proposed method was tested on several systems revealing good results. It was also compared with other tuning methods showing comparable or superior results [8].


## 3. Fractional PDα  Control of a Hexapod Robot

Walking machines allow locomotion in terrain inaccessible to other type of vehicles, since they do not need a continuous support surface, but at the cost of higher requirements for leg coordination and control. For these robots, joint level control is usually implemented through a PID-like scheme with position feedback. Recently, the application of the theory of FC to robotics revealed promising aspects for future developments [9]. With these facts in mind, this section compares different Fractional PDα robot controller tuning, applied to the joint control of a walking system (Figure 5) with n=6 legs, equally distributed along both sides of the robot body, having each three rotational joints (i.e.,  j={1,2,3}≡{hip,knee,ankle}) [10].Figure 5
Model of the robot body and foot-ground interaction.During this study leg jointj=3 can be either mechanical actuated or motor actuated (Figure 5). For the mechanical actuated case, we suppose that there is a rotational pre-tensioned spring-dashpot system connecting leg links Li2 and Li3. This mechanical impedance maintains the angle between the two links while imposing a joint torque [10].Figure5 presents the dynamic model for the hexapod body and foot-ground interaction. It is considered robot body compliance because walking animals have a spine that allows supporting the locomotion with improved stability. The robot body is divided in n identical segments (each with mass Mbn-1) and a linear spring-damper system (with parameters defined so that the body behaviour is similar to the one expected to occur on an animal) is adopted to implement the intrabody compliance [10]. The contact of the ith robot feet with the ground is modelled through a nonlinear system [11], being the values for the parameters based on the studies of soil mechanics [11].The general control architecture of the hexapod robot is presented in Figure6 [12]. In this study we evaluate the effect of different PDα, α∈ℜ, controller implementations for Gc1(s), while Gc2 is a proportional controller with gain Kpj=0.9(j=1,2,3). For the PDα algorithm, implemented through a discrete-time 4th-order Padé approximation (aij, bij∈ℜ, j=1,2,3), we have
(3.1)Gc1j(z)≈Kpj+Kαj∑i=0i=uaijz-i∑i=0i=ubijz-i,
where Kpj and Kαj are the proportional and derivative gains, respectively, and αj is the fractional order, for joint j. Therefore, the classical PD1  algorithm occurs when the fractional order αj=1.0.Figure 6
Hexapod robot control architecture.It is analysed the system performance of the differentPDα tuning, during a periodic wave gait at a constant forward velocity VF, for two cases: two leg joints are motor actuated and the ankle joint is mechanical actuated and the three leg joints are fully motor actuated [10].The analysis is based on the formulation of two indices measuring the mean absolute density of energy per traveled distance (Eav) and the hip trajectory errors (εxyH) during walking, according to
(3.2)Eav=1d∑i=1n∑j=1m∫0T|τij(t)θ̇ij(t)|dt[Jm-1],εxyH=∑i=1n1Ns∑k=1Ns(ΔixH2+ΔiyH2)[m],ΔixH=xiHd(k)-xiH(k),ΔiyH=yiHd(k)-yiH(k).To tune the different controller implementations we adopt a systematic method, testing and evaluating several possible combinations of parameters, for all controller implementations. Therefore, we adopt theGc1(s) parameters that establish a compromise in what concerns the simultaneous minimisation of Eav and εxyH. Moreover, it is assumed high-performance joint actuators, with a maximum actuator torque of τijMax=400 Nm, and the desired angle between the foot and the ground (assumed horizontal) is established as θi3hd=-15°. We tune the PDα joint controllers for different values of the fractional order αj while making α1=α2=α3.We start by considering that leg joints 1 and 2 are motor actuated and joint 3 is mechanical actuated. For this case we tune the PDα joint controllers for different values of the fractional order αj, with step Δαj=0.1, namely, αj={-0.9,-0.8,…,+0.9}. Afterwards, we consider that joint 3 is also motor actuated, and we repeat the controller tuning procedure versus αj.For the first situation under study, we verify that the value ofαj=0.6 (Figure 7), being the gains of the PDα controller Kp1=2500, Kα1=800, Kp2=300, Kα2=100 and the parameters of the mechanical spring-dashpot system for the ankle actuation K3=1, B3=2, presents the best compromise situation between the simultaneous minimisation of εxyH and Eav.Figure 7
Locus ofEav versus εxyH  for the different values of α in the Gc1(s)  tuning, when establishing a compromise between the minimisation of Eav  and εxyH, with Gc2=0.9, joints 1 and 2 motor actuated and joint 3 mechanical actuated.Regarding the case when all joints are motor actuated, Figure8 presents the best controller tuning for different values of αj. The experiments reveal the superior performance of the PDα controller for αj≈0.5, with Kp1=15000, Kα1=7200, Kp2=1000, Kα2=800, and Kp3=150, Kα3=240.Figure 8
Locus ofEav versus εxyH  for the different values of α in the Gc1(s) tuning, when establishing a compromise between the minimisation of Eav  and εxyH, with Gc2=0.9 and all joints motor actuated.Forαj={0.1,0.2,0.3,0.4} the results are very poor and for αj={-0.9,…,-0.1}∪{+0.9}, the hexapod locomotion is unstable. Furthermore, we conclude that the best case corresponds to all leg joints being motor actuated.In conclusion, the experiments reveal the superior performance of the FO controller forαj≈0.5 and a robot with all motor actuated joints, as can be concluded analysing the curves for the joint actuation torques τ1jm (Figure 9) and for the hip trajectory tracking errors Δ1xH  and Δ1yH  (Figure 10).Plots ofτ1jm  versus t, with joints 1 and 2 motor actuated and joint 3 mechanical actuated and all joints motor actuated, for αj=0.5.
(a)(b)(c)Plots ofΔ1xH and Δ1yH versus t, with joints 1 and 2 motor actuated and joint 3 mechanical actuated and all joints motor actuated, for αj=0.5.
(a)(b)Since the objective of the walking robots is to walk in natural terrains, in the sequel it is examined how the different controller tunings behave under different ground properties, considering that all joints are motor actuated. For this case, and considering the previously tuning controller parameters, the values of{KxF,BxF,KyF,ByF} are varied simultaneously through a multiplying factor Kmult  that is varied in the range [0.1,4.0]. This variation for the ground model parameters allows the simulation of the ground behaviour for growing stiffness, from peat to gravel [11].The performance measureEav versus the multiplying factor of the ground parameters Kmult is presented on Figure 11. Analysing the system performance from the viewpoint of the index Eav, it is possible to conclude that the best PDα  implementation occurs for the fractional order αj=0.5. Moreover, it is clear that the performances of the different controller implementations are almost constant on all range of the ground parameters, with the exception of the fractional order αj=0.4. For this case, Eav presents a significant variation with Kmult. Therefore, we conclude that the controller responses are quite similar, meaning that these algorithms are robust to variations of the ground characteristics [12].Figure 11
Performance indexEav versus Kmult  for the Gc1(s) PDα controller tuning having all joints motor actuated.


## 4. Fractional Dynamics in the Trajectory Control of Redundant Manipulators

A redundant manipulator is a robotic arm possessing more degrees of freedom (dof) than those required to establish an arbitrary position and orientation of the end effector. Redundant manipulators offer several potential advantages over non-redundant arms. In a workspace with obstacles, the extra degrees of freedom can be used to move around or between obstacles and thereby to manipulate in situations that otherwise would be inaccessible [20–23].When a manipulator is redundant, it is anticipated that the inverse kinematics admits an infinite number of solutions. This implies that, for a given location of the manipulator’s gripper, it is possible to induce a self-motion of the structure without changing the location of the end effecter. Therefore, the arm can be reconfigured to find better postures for an assigned set of task requirements.Several kinematic techniques for redundant manipulators control the gripper through the rates at which the joints are driven, using the pseudoinverse of the Jacobian [22, 24]. Nevertheless, these algorithms lead to a kind of chaotic motion with unpredictable arm configurations.Having these ideas in mind, Section4.1 introduces the fundamental issues for the kinematics of redundant manipulators. Based on these concepts, Section 4.2 presents the trajectory control of a three dof  robot. The results reveal a chaotic behavior that is further analyzed in Section 4.3.


### 4.1. Kinematics of Redundant Manipulators

A kinematically redundant manipulator has moredof  than those required to define an arbitrary position and orientation of the gripper. In Figure 12 is depicted a planar manipulator with k∈ℵ rotational (R) joints that is redundant for k>2. When a manipulator is redundant it is anticipated that the inverse kinematics admits an infinite number of solutions. This implies that, for a given location of the manipulator’s gripper, it is possible to induce a self-motion of the structure without changing the location of the gripper. Therefore, redundant manipulators can be reconfigured to find better postures for an assigned set of task requirements but, on the other hand, have a more complex structure requiring adequate control algorithms.Figure 12
A planar redundant planar manipulator withk rotational joints.We consider a manipulator withn degrees of freedom whose joint variables are denoted by q=[q1,q2,…,qn]T. We assume that a class of tasks, we are interested in can be described by m variables, x=[x1,x2,…,xm]T(m<n) and that the relation between q and x is given by
(4.1)x=f(q),
where f is a function representing the direct kinematics.Differentiating (4.1) with respect to time yields
(4.2)ẋ=J(q)q̇,
where ẋ∈ℜm, q̇∈ℜn, and J(q)=∂f(q)/∂q∈ℜm×n. Hence, it is possible to calculate a path q(t) in terms of a prescribed trajectory x(t) in the operational space. We assume that the following condition is satisfied:
(4.3)maxrank{J(q)}=m.Failing to satisfy this condition usually means that the selection of manipulation variables is redundant and the number of these variablesm can be reduced. When condition (4.3) is verified, we say that the degree of redundancy of the manipulator is n-m. If, for some q we have
(4.4)rank{J(q)}<m
then the manipulator is in a singular state. This state is not desirable because, in this region of the trajectory, the manipulating ability is very limited.Many approaches for solving redundancy [25, 26] are based on the inversion of (4.2). A solution in terms of the joint velocities is sought as
(4.5)q̇=J#(q)ẋ,
where J# is one of the generalized inverses of the J  [26–28]. It can be easily shown that a more general solution to (4.2) is given by
(4.6)q̇=J+(q)ẋ+[I-J+(q)J(q)]q̇0,
where I  is the n×n identity matrix and q̇0∈ℜn  is a n×1 arbitrary joint velocity vector and J+ is the pseudoinverse of the J. The solution (4.6) is composed of two terms. The first term is relative to minimum norm joint velocities. The second term, the homogeneous solution, attempts to satisfy the additional constraints specified by q̇0. Moreover, the matrix I-J+(q)J(q) allows the projection of q̇0 in the null space of J. A direct consequence is that it is possible to generate internal motions that reconfigure the manipulator structure without changing the gripper position and orientation [27–30]. Another aspect revealed by the solution of (4.6) is that repetitive trajectories in the operational space do not lead to periodic trajectories in the joint space. This is an obstacle for the solution of many tasks because the resultant robot configurations have similarities with those of a chaotic system.


### 4.2. Robot Trajectory Control

The direct kinematics and the Jacobian of a 3-link planar manipulator with rotational joints (3R robot) has a simple recursive nature according with the expressions:
(4.7)[x1x2]=[l1C1+l2C12+l3C123l1S1+l2S12+l3S123],J=[-l1S1-⋯-l3S123⋯-l3S123l1C1+⋯+l3C123⋯l3C123],
where li is the length of link i, qi…k=qi+⋯+qk, Si…k=Sin(qi…k), and Ci…k=Cos(qi…k).During all the experiments it is consideredΔt=10-3 seconds, LTOT=l1+l2+l3=3, and l1=l2=l3.In the closed-loop pseudoinverse’s method the joint positions can be computed through the time integration of the velocities according with the block diagram of the inverse kinematics algorithm depicted in Figure13, where xref represents the vector of reference coordinates of the robot gripper in the operational space.Figure 13
Block diagram of the closed-loop inverse kinematics algorithm with the pseudoinverse.Based on (4.7) we analyze the kinematic performances of the 3R-robot when repeating a circular motion in the operational space with frequency ω0=7.0 rad s-1, centre at distance r=[x12+x22]1/2  and radius ρ.Figure14 shows the joint positions for the inverse kinematic algorithm (4.5) for r={0.6,2.0} and ρ={0.3,0.5}. We observe that the following hold.The3R-robot joint positions versus time using the pseudoinverse method for   r={0.6,2.0} and ρ={0.3,0.5}.
(a)(b)(c)(d)(i)
Forr=0.6 occur unpredictable motions with severe variations that lead to high joint transients [13]. Moreover, we verify a low-frequency signal modulation that depends on the circle being executed.(ii)
Forr=2.0 the motion is periodic with frequency identical to ω0=7.0 rad s-1.


### 4.3. Analysis of the Robot Trajectories

In the previous subsection we verified that the pseudoinverse-based algorithm leads to unpredictable arm configurations. In order to gain further insight into the pseudoinverse nature several distinct experiments are devised in the sequel during a time window of 300 cycles. Therefore, in a first set of experiments we calculate the Fourier transform of the3R-robot joints velocities for a circular repetitive motion with frequency ω0=7.0 rad s-1, radius ρ={0.1,0.3,0.5,0.7}, and radial distances r∈]0,LTOT-ρ[.Figure15 shows |F{q̇2(t)}| versus the frequency ratio ω0/ω and the distance r where F{} represents the Fourier operator. Is verified an interesting phenomenon induced by the gripper repetitive motion ω0 because a large part of the energy is distributed along several subharmonics. These fractional-order harmonics (foh) depend on r and ρ making a complex pattern with similarities with those revealed by chaotic systems. Furthermore, we observe the existence of several distinct regions depending on r.|F{q̇2(t)}| of the 3R-robot during 300 cycles versus r and ω/ω0, for  ρ={0.1,0.3,0.5,0.7}, ω0=7.0 rad s-1.
(a)(b)(c)(d)For example, selecting in Figure15 several distinct cases, namely for r={0.08,0.30,0.53,1.10,1.30,2.00}, we have the different signal Fourier spectra clearly visible in Figure 16.|F{q̇2(t)}| of the 3R-robot during 300 cycles versus the frequency ratio ω/ω0, for r={0.08,0.30,0.53,1.10,1.30,2.00}, ρ=0.7, ω0=7.0 rad s-1.
(a)(b)(c)In the authors’ best knowledge thefoh  are aspects of fractional dynamics [14, 15, 31], but a final and assertive conclusion about a physical interpretation is a matter still to be explored.For joints velocities 1 and 3 the results are similar to the verified ones for joint velocity 2.


## 4.1. Kinematics of Redundant Manipulators

A kinematically redundant manipulator has moredof  than those required to define an arbitrary position and orientation of the gripper. In Figure 12 is depicted a planar manipulator with k∈ℵ rotational (R) joints that is redundant for k>2. When a manipulator is redundant it is anticipated that the inverse kinematics admits an infinite number of solutions. This implies that, for a given location of the manipulator’s gripper, it is possible to induce a self-motion of the structure without changing the location of the gripper. Therefore, redundant manipulators can be reconfigured to find better postures for an assigned set of task requirements but, on the other hand, have a more complex structure requiring adequate control algorithms.Figure 12
A planar redundant planar manipulator withk rotational joints.We consider a manipulator withn degrees of freedom whose joint variables are denoted by q=[q1,q2,…,qn]T. We assume that a class of tasks, we are interested in can be described by m variables, x=[x1,x2,…,xm]T(m<n) and that the relation between q and x is given by
(4.1)x=f(q),
where f is a function representing the direct kinematics.Differentiating (4.1) with respect to time yields
(4.2)ẋ=J(q)q̇,
where ẋ∈ℜm, q̇∈ℜn, and J(q)=∂f(q)/∂q∈ℜm×n. Hence, it is possible to calculate a path q(t) in terms of a prescribed trajectory x(t) in the operational space. We assume that the following condition is satisfied:
(4.3)maxrank{J(q)}=m.Failing to satisfy this condition usually means that the selection of manipulation variables is redundant and the number of these variablesm can be reduced. When condition (4.3) is verified, we say that the degree of redundancy of the manipulator is n-m. If, for some q we have
(4.4)rank{J(q)}<m
then the manipulator is in a singular state. This state is not desirable because, in this region of the trajectory, the manipulating ability is very limited.Many approaches for solving redundancy [25, 26] are based on the inversion of (4.2). A solution in terms of the joint velocities is sought as
(4.5)q̇=J#(q)ẋ,
where J# is one of the generalized inverses of the J  [26–28]. It can be easily shown that a more general solution to (4.2) is given by
(4.6)q̇=J+(q)ẋ+[I-J+(q)J(q)]q̇0,
where I  is the n×n identity matrix and q̇0∈ℜn  is a n×1 arbitrary joint velocity vector and J+ is the pseudoinverse of the J. The solution (4.6) is composed of two terms. The first term is relative to minimum norm joint velocities. The second term, the homogeneous solution, attempts to satisfy the additional constraints specified by q̇0. Moreover, the matrix I-J+(q)J(q) allows the projection of q̇0 in the null space of J. A direct consequence is that it is possible to generate internal motions that reconfigure the manipulator structure without changing the gripper position and orientation [27–30]. Another aspect revealed by the solution of (4.6) is that repetitive trajectories in the operational space do not lead to periodic trajectories in the joint space. This is an obstacle for the solution of many tasks because the resultant robot configurations have similarities with those of a chaotic system.


## 4.2. Robot Trajectory Control

The direct kinematics and the Jacobian of a 3-link planar manipulator with rotational joints (3R robot) has a simple recursive nature according with the expressions:
(4.7)[x1x2]=[l1C1+l2C12+l3C123l1S1+l2S12+l3S123],J=[-l1S1-⋯-l3S123⋯-l3S123l1C1+⋯+l3C123⋯l3C123],
where li is the length of link i, qi…k=qi+⋯+qk, Si…k=Sin(qi…k), and Ci…k=Cos(qi…k).During all the experiments it is consideredΔt=10-3 seconds, LTOT=l1+l2+l3=3, and l1=l2=l3.In the closed-loop pseudoinverse’s method the joint positions can be computed through the time integration of the velocities according with the block diagram of the inverse kinematics algorithm depicted in Figure13, where xref represents the vector of reference coordinates of the robot gripper in the operational space.Figure 13
Block diagram of the closed-loop inverse kinematics algorithm with the pseudoinverse.Based on (4.7) we analyze the kinematic performances of the 3R-robot when repeating a circular motion in the operational space with frequency ω0=7.0 rad s-1, centre at distance r=[x12+x22]1/2  and radius ρ.Figure14 shows the joint positions for the inverse kinematic algorithm (4.5) for r={0.6,2.0} and ρ={0.3,0.5}. We observe that the following hold.The3R-robot joint positions versus time using the pseudoinverse method for   r={0.6,2.0} and ρ={0.3,0.5}.
(a)(b)(c)(d)(i)
Forr=0.6 occur unpredictable motions with severe variations that lead to high joint transients [13]. Moreover, we verify a low-frequency signal modulation that depends on the circle being executed.(ii)
Forr=2.0 the motion is periodic with frequency identical to ω0=7.0 rad s-1.


## 4.3. Analysis of the Robot Trajectories

In the previous subsection we verified that the pseudoinverse-based algorithm leads to unpredictable arm configurations. In order to gain further insight into the pseudoinverse nature several distinct experiments are devised in the sequel during a time window of 300 cycles. Therefore, in a first set of experiments we calculate the Fourier transform of the3R-robot joints velocities for a circular repetitive motion with frequency ω0=7.0 rad s-1, radius ρ={0.1,0.3,0.5,0.7}, and radial distances r∈]0,LTOT-ρ[.Figure15 shows |F{q̇2(t)}| versus the frequency ratio ω0/ω and the distance r where F{} represents the Fourier operator. Is verified an interesting phenomenon induced by the gripper repetitive motion ω0 because a large part of the energy is distributed along several subharmonics. These fractional-order harmonics (foh) depend on r and ρ making a complex pattern with similarities with those revealed by chaotic systems. Furthermore, we observe the existence of several distinct regions depending on r.|F{q̇2(t)}| of the 3R-robot during 300 cycles versus r and ω/ω0, for  ρ={0.1,0.3,0.5,0.7}, ω0=7.0 rad s-1.
(a)(b)(c)(d)For example, selecting in Figure15 several distinct cases, namely for r={0.08,0.30,0.53,1.10,1.30,2.00}, we have the different signal Fourier spectra clearly visible in Figure 16.|F{q̇2(t)}| of the 3R-robot during 300 cycles versus the frequency ratio ω/ω0, for r={0.08,0.30,0.53,1.10,1.30,2.00}, ρ=0.7, ω0=7.0 rad s-1.
(a)(b)(c)In the authors’ best knowledge thefoh  are aspects of fractional dynamics [14, 15, 31], but a final and assertive conclusion about a physical interpretation is a matter still to be explored.For joints velocities 1 and 3 the results are similar to the verified ones for joint velocity 2.


## 5. Heat Diffusion

The heat diffusion is governed by a linear one-dimensional partial differential equation (PDE) of the form:(5.1)∂u∂t=k∂2u∂x2,
where k is the diffusivity, t is the time, u is the temperature, and x is the space coordinate. However, (5.1) involves the solution of a PDE of parabolic type for which the standard theory guarantees the existence of a unique solution [16].For the case of a planar perfectly isolated surface we usually apply a constant temperatureU0 at x=0 and analyzes the heat diffusion along the horizontal coordinate x. Under these conditions, the heat diffusion phenomenon is described by a noninteger-order model:
(5.2)U(x,s)=U0sG(s)G(s)=e-xs/k,
where x is the space coordinate, U0 is the boundary condition, and G(s) is the system transfer function.In our study, the simulation of the heat diffusion is performed by adopting the Crank-Nicholson implicit numerical integration based on the discrete approximation to differentiation as [16, 17]
(5.3)-ru[j+1,i+1]+(2+r)u[j+1,i]-ru[j+1,i-1]=ru[j,i+1]+(2-r)u[j,i]+u[j,i-1],
where r=Δt(Δx2)-1, {Δx,Δt}, and {i,j} are the increments and the integration indices for space and time, respectively.


### 5.1. Control Strategies

The generalized PID controllerGc(s) has a transfer function of the form
(5.4)Gc(s)=K[1+1Tisα+Tdsβ],
where α and β are the orders of the fractional integrator and differentiator, respectively. The constants K, Ti, and Td are correspondingly the proportional gain, the integral time constant, and the derivative time constant.Clearly, taking(α,β)={(1,1),(1,0),(0,1),(0,0)} we get the classical {PID,PI,PD,P} controllers, respectively.ThePIαDβ controller is more flexible and gives the possibility of adjusting more carefully the closed-loop system characteristics.In the following two subsections, we analyze the system of Figure17 by adopting the classical integer-order PID and a fractional PIDβ, respectively.Figure 17
Closed-loop system with PID controllerGc(s).


### 5.2. PID Tuning Using the Ziegler-Nichols Rule

In this subsection, we analyze the closed-loop system with a conventional PID controller given by the transfer function (5.4) with α=β=1. Usually, the PID parameters (K,Ti,Td) are tuned by using the so-called Ziegler-Nichols open loop (ZNOL) method [17]. The ZNOL heuristics are based on the approximate first-order plus dead-time model:
(5.5)Ĝ(s)=Kpτs+1e-sT.For the heat system, the resulting parameters are{Kp,τ,T}={0.52,162,28} leading to the PID constants {K,Ti,Td}={18.07,34.0,8.5}.A step input is applied atx=0.0 m and the closed-loop response c(t) is analyzed for x=3.0 m, without actuator saturation (Figure 18). We verify that the system with a PID controller, tuned through the ZNOL heuristics, does not produce satisfactory results giving a significant overshoot ov and a large settling time ts, namely, {ts,tp,tr,ov(%)}≡{44.8,27.5,12.0,68.56}, where tp represents the peak time and tr the rise time. We consider two indices that measure the response error, namely, the integral square error (ISE) and the integral time square error (ITSE) criteria defined as
(5.6)ISE=∫0∞[r(t)-c(t)]2dt,ITSE=∫0∞t[r(t)-c(t)]2dt.Figure 18
Step responses of the closed-loop system for the PID controller andx=3.0 m.We can use other performance criteria such as the integral absolute error (IAE) or the integral time absolute error (ITAE); however, in the present case, the ISE and the ITSE criteria have produced the best results and are adopted in the study.In this case, the ZNOL PID tuning leads to the values(ISE,ITSE)=(27.53,613.97). The poor results indicate again that the method of tuning may not be the most adequate for the control of the heat system.In fact, the inherent fractional dynamics of the system lead us to consider other configurations. In this perspective, we propose the use of fractional controllers tuned by the minimization of the indices ISE and ITSE.


### 5.3. PIDβ Tuning Using Optimization Indices

In this subsection, we analyze the closed-loop system under the action of the PIDβ controller given by the transfer function (5.4) with α=1 and 0≤β≤1. The fractional derivative term Tdsβ in (5.4) is implemented through a fourth-order Padé discrete rational transfer function. It used a sampling period of T=0.1 second.The PIDβ controller is tuned by the minimization of an integral performance index. For that purpose, we adopt the ISE and ITSE criteria.A step reference inputR(s)=1/s is applied at x=0.0 m and the output c(t) is analyzed for x=3.0 m, without actuator saturation. The heat system is simulated for 3000 seconds. Figure 19 illustrates the variation of the fractional PID parameters (K,Ti,Td) as function of the order’s derivative β, for the ISE and the ITSE criteria. The dots represent the values corresponding to the classical PID (ZNOL-tuning) addressed in the previous section.The PIDβ parameters (K,Ti,Td) versus β for the ISE and ITSE optimization criteria. The dot represents the PID-ZNOL.
(a)(b)(c)The curves reveal that forβ<0.4 the parameters (K,Ti,Td) are slightly different, for the two ISE and ITSE criteria, while for β≥0.4 they lead to almost similar values. This fact indicates a large influence of a weak-order derivative on system’s dynamics.To further illustrate the performance of the fractional-order controllers a saturation nonlinearity is included in the closed-loop system of Figure17 and inserted in series with the output of the controller Gc(s). The saturation element is defined as
(5.7)n(m)={m,|m|<δ,δsign(m),|m|≥δ.The controller performance is evaluated forδ={20,…,100} and δ=∞ which corresponds to a system without saturation. We use the same fractional-PID parameters obtained without considering the saturation nonlinearity.Figures20 and 21 show the step responses of the closed-loop system and the corresponding controller output, for the PIDβ tuned in the ISE and ITSE perspectives for δ=10 and δ=∞, respectively. The controller parameters {K,Ti,Td,β} correspond to the minimization of those indices leading to the values ISE: {K,Ti,Td,β}≡{3,23,90.6,0.875} and ITSE: {K,Ti,Td,β}≡{1.8,17.6,103.6,0.85}.Step responses of the closed-loop system and the controller output for the ISE and the ITSE indices, with a PIDβ controller, δ=10  and x=3.0 m.
(a)(b)Step responses of the closed-loop system and the controller output for the ISE and the ITSE indices, with a PIDβ controller, δ=∞ and x=3.0 m.
(a)(b)The step responses reveal a large diminishing of the overshoot and the rise time when compared with the integer PID, showing a good transient response and a zero steady-state error. The PIDβ leads to better results than the classical PID controller tuned through the ZNOL rule. These results demonstrate the effectiveness of the fractional algorithms when used for the control of fractional-order systems. The step response and the controller output are also improved when the saturation level δ is diminished.Figure22 depicts the ISE and ITSE indices for 0≤β≤1, when δ={20,…,100} and δ=∞. We verify the existence of a minimum for β=0.875 and β=0.85 for the ISE and ITSE cases, respectively. Furthermore, the higher the δ  the lower the value of the index.ISE and ITSE versus0≤β≤1 for δ={20,…,100} and δ=∞.
(a)(b)Figures23 and 24 show the variation of the settling time ts, the peak time tp, the rise time tr, and the percent overshoot ov(%), for the closed-loop response tuned through the minimization of the ISE and the ITSE indices, respectively.Parametersts, tp, tr, ov(%) for the step responses of the closed-loop system for the ISE indice, with a PIDβ controller, when δ={20,…,100} and δ=∞, x=3.0 m.
(a)(b)(c)(d)Parametersts, tp, tr, ov(%) for the step responses of the closed-loop system for the ITSE indice, with a PIDβ controller, when δ={20,…,100} and δ=∞, x=3.0 m.
(a)(b)(c)(d)In the ISE casets, tp, andtr diminish rapidly for 0≤β≤0.875, while for β>0.875 the parameters increase smoothly. For the ITSE, we verify the same behavior for β=0.85. On the other hand, ov(%) increases smoothly for 0≤β≤0.7, while for β>0.7 it decreases very quickly, both for the ISE and the ITSE indices.In conclusion, for0.85≤β≤0.875 we get the best controller tuning, superior to the performance revealed by the classical integer-order scheme.


## 5.1. Control Strategies

The generalized PID controllerGc(s) has a transfer function of the form
(5.4)Gc(s)=K[1+1Tisα+Tdsβ],
where α and β are the orders of the fractional integrator and differentiator, respectively. The constants K, Ti, and Td are correspondingly the proportional gain, the integral time constant, and the derivative time constant.Clearly, taking(α,β)={(1,1),(1,0),(0,1),(0,0)} we get the classical {PID,PI,PD,P} controllers, respectively.ThePIαDβ controller is more flexible and gives the possibility of adjusting more carefully the closed-loop system characteristics.In the following two subsections, we analyze the system of Figure17 by adopting the classical integer-order PID and a fractional PIDβ, respectively.Figure 17
Closed-loop system with PID controllerGc(s).


## 5.2. PID Tuning Using the Ziegler-Nichols Rule

In this subsection, we analyze the closed-loop system with a conventional PID controller given by the transfer function (5.4) with α=β=1. Usually, the PID parameters (K,Ti,Td) are tuned by using the so-called Ziegler-Nichols open loop (ZNOL) method [17]. The ZNOL heuristics are based on the approximate first-order plus dead-time model:
(5.5)Ĝ(s)=Kpτs+1e-sT.For the heat system, the resulting parameters are{Kp,τ,T}={0.52,162,28} leading to the PID constants {K,Ti,Td}={18.07,34.0,8.5}.A step input is applied atx=0.0 m and the closed-loop response c(t) is analyzed for x=3.0 m, without actuator saturation (Figure 18). We verify that the system with a PID controller, tuned through the ZNOL heuristics, does not produce satisfactory results giving a significant overshoot ov and a large settling time ts, namely, {ts,tp,tr,ov(%)}≡{44.8,27.5,12.0,68.56}, where tp represents the peak time and tr the rise time. We consider two indices that measure the response error, namely, the integral square error (ISE) and the integral time square error (ITSE) criteria defined as
(5.6)ISE=∫0∞[r(t)-c(t)]2dt,ITSE=∫0∞t[r(t)-c(t)]2dt.Figure 18
Step responses of the closed-loop system for the PID controller andx=3.0 m.We can use other performance criteria such as the integral absolute error (IAE) or the integral time absolute error (ITAE); however, in the present case, the ISE and the ITSE criteria have produced the best results and are adopted in the study.In this case, the ZNOL PID tuning leads to the values(ISE,ITSE)=(27.53,613.97). The poor results indicate again that the method of tuning may not be the most adequate for the control of the heat system.In fact, the inherent fractional dynamics of the system lead us to consider other configurations. In this perspective, we propose the use of fractional controllers tuned by the minimization of the indices ISE and ITSE.


## 5.3. PIDβ Tuning Using Optimization Indices

In this subsection, we analyze the closed-loop system under the action of the PIDβ controller given by the transfer function (5.4) with α=1 and 0≤β≤1. The fractional derivative term Tdsβ in (5.4) is implemented through a fourth-order Padé discrete rational transfer function. It used a sampling period of T=0.1 second.The PIDβ controller is tuned by the minimization of an integral performance index. For that purpose, we adopt the ISE and ITSE criteria.A step reference inputR(s)=1/s is applied at x=0.0 m and the output c(t) is analyzed for x=3.0 m, without actuator saturation. The heat system is simulated for 3000 seconds. Figure 19 illustrates the variation of the fractional PID parameters (K,Ti,Td) as function of the order’s derivative β, for the ISE and the ITSE criteria. The dots represent the values corresponding to the classical PID (ZNOL-tuning) addressed in the previous section.The PIDβ parameters (K,Ti,Td) versus β for the ISE and ITSE optimization criteria. The dot represents the PID-ZNOL.
(a)(b)(c)The curves reveal that forβ<0.4 the parameters (K,Ti,Td) are slightly different, for the two ISE and ITSE criteria, while for β≥0.4 they lead to almost similar values. This fact indicates a large influence of a weak-order derivative on system’s dynamics.To further illustrate the performance of the fractional-order controllers a saturation nonlinearity is included in the closed-loop system of Figure17 and inserted in series with the output of the controller Gc(s). The saturation element is defined as
(5.7)n(m)={m,|m|<δ,δsign(m),|m|≥δ.The controller performance is evaluated forδ={20,…,100} and δ=∞ which corresponds to a system without saturation. We use the same fractional-PID parameters obtained without considering the saturation nonlinearity.Figures20 and 21 show the step responses of the closed-loop system and the corresponding controller output, for the PIDβ tuned in the ISE and ITSE perspectives for δ=10 and δ=∞, respectively. The controller parameters {K,Ti,Td,β} correspond to the minimization of those indices leading to the values ISE: {K,Ti,Td,β}≡{3,23,90.6,0.875} and ITSE: {K,Ti,Td,β}≡{1.8,17.6,103.6,0.85}.Step responses of the closed-loop system and the controller output for the ISE and the ITSE indices, with a PIDβ controller, δ=10  and x=3.0 m.
(a)(b)Step responses of the closed-loop system and the controller output for the ISE and the ITSE indices, with a PIDβ controller, δ=∞ and x=3.0 m.
(a)(b)The step responses reveal a large diminishing of the overshoot and the rise time when compared with the integer PID, showing a good transient response and a zero steady-state error. The PIDβ leads to better results than the classical PID controller tuned through the ZNOL rule. These results demonstrate the effectiveness of the fractional algorithms when used for the control of fractional-order systems. The step response and the controller output are also improved when the saturation level δ is diminished.Figure22 depicts the ISE and ITSE indices for 0≤β≤1, when δ={20,…,100} and δ=∞. We verify the existence of a minimum for β=0.875 and β=0.85 for the ISE and ITSE cases, respectively. Furthermore, the higher the δ  the lower the value of the index.ISE and ITSE versus0≤β≤1 for δ={20,…,100} and δ=∞.
(a)(b)Figures23 and 24 show the variation of the settling time ts, the peak time tp, the rise time tr, and the percent overshoot ov(%), for the closed-loop response tuned through the minimization of the ISE and the ITSE indices, respectively.Parametersts, tp, tr, ov(%) for the step responses of the closed-loop system for the ISE indice, with a PIDβ controller, when δ={20,…,100} and δ=∞, x=3.0 m.
(a)(b)(c)(d)Parametersts, tp, tr, ov(%) for the step responses of the closed-loop system for the ITSE indice, with a PIDβ controller, when δ={20,…,100} and δ=∞, x=3.0 m.
(a)(b)(c)(d)In the ISE casets, tp, andtr diminish rapidly for 0≤β≤0.875, while for β>0.875 the parameters increase smoothly. For the ITSE, we verify the same behavior for β=0.85. On the other hand, ov(%) increases smoothly for 0≤β≤0.7, while for β>0.7 it decreases very quickly, both for the ISE and the ITSE indices.In conclusion, for0.85≤β≤0.875 we get the best controller tuning, superior to the performance revealed by the classical integer-order scheme.


## 6. Circuit Synthesis Using Evolutionary Algorithms

In recent decades evolutionary computation (EC) techniques have been applied to the design of electronic circuits and systems, leading to a novel area of research called Evolutionary Electronics (EE) or Evolvable Hardware (EH). EE considers the concept for automatic design of electronic systems. Instead of using human conceived models, abstractions, and techniques, EE employs search algorithms to develop implementations not achievable with the traditional design schemes, such as the Karnaugh or the Quine-McCluskey Boolean methods.Several papers proposed designing combinational logic circuits using evolutionary algorithms and, in particular, genetic algorithms (GAs) [32, 33] and hybrid schemes such as the memetic algorithms (MAs) [34].Particle swarm optimization (PSO) constitutes an alternative evolutionary computation technique, and this paper studies its application to combinational logic circuit synthesis. Bearing these ideas in mind, the organization of this section is as follows. Section6.1 presents a brief overview of the PSO. Section 6.2 describes the PSO-based circuit design, while Section 6.3 exhibits the simulation results.


### 6.1. Particle Swarm Optimization

In literature about PSO the term ‘swarm intelligence’ appears rather often and, therefore, we begin by explaining why this is so.Noncomputer scientists (ornithologists, biologists, and psychologists) did early research, which led into the theory of particle swarms. In these areas, the term “swarm intelligence’’ is well known and characterizes the case when a large number of individuals are able of accomplish complex tasks. Motivated by these facts, some basic simulations of swarms were abstracted into the mathematical field. The usage of swarms for solving simple tasks in nature became an intriguing idea in algorithmic and function optimization.Eberhart and Kennedy were the first to introduce the PSO algorithm [35], which is an optimization method inspired in the collective intelligence of swarms of biological populations, and was discovered through simplified social model simulation of bird flocking, fishing schooling, and swarm theory.In the PSO, instead of using genetic operators, as in the case of GAs, each particle (individual) adjusts its flying according with its own and its companions experiences. Each particle is treated as a point in a D-dimensional space and is manipulated as described in what follows in the original PSO algorithm:(6.1a)vid=vid+c1rand()(pid-xid)+c2Rand()(pgd-xid),(6.1b)xid=xid+vid,
where c1 and c2 are positive constants, rand() and Rand() are two random functions in the range [0,1], Xi=(xi1,xi2,…,xiD) represents the ith particle, Pi=(pi1,pi2,…,piD)  is the best previous position (the position giving the best fitness value) of the particle, the symbol g represents the index of the best particle among all particles in the population, and Vi=(vi1,vi2,…,viD) is the rate of the position change (velocity) for particle i.However, (6.1a) and (6.1b) represent the flying trajectory of a population of particles. Also, (6.1a) describes how the velocity is dynamically updated and (6.1b) the position update of the “flying” particles. Moreover, (6.1b) is divided in three parts, namely the momentum, the cognitive and the social parts. In the first part the velocity cannot be changed abruptly: it is adjusted based on the current velocity. The second part represents the learning from its own flying experience. The third part consists on the learning group flying experience [36].The first new parameter added into the original PSO algorithm is the inertia weigh. The dynamic equation of PSO with inertia weigh is modified to be(6.2a)vid=wvid+c1rand()(pid-xid)+c2Rand()(pgd-xid),(6.2b)xid=xid+vid,
where w constitutes the inertia weigh that introduces a balance between the global and the local search abilities. A large inertia weigh facilitates a global search while a small inertia weigh facilitates a local search.Another parameter, called constriction coefficientk, is introduced with the hope that it can insure a PSO to converge. A simplified method of incorporating it appears in (6.3), where k is function of c1 and c2 as it is presented as follows:
(6.3)vid=k[vid+c1rand()(pid-xid)+c2Rand()(pgd-xid)],xid=xid+vid,(6.4)k=2(2-ϕ-ϕ2-4ϕ)-1,
where ϕ=c1+c2, ϕ>4.There are two different PSO topologies, namely, the global version and the local version. In the global version of PSO, each particle flies through the search space with a velocity that is dynamically adjusted according to the particle’s personal best performance achieved so far and the best performance achieved so far by all particles. On the other hand, in the local version of PSO, each particle’s velocity is adjusted according to its personal best and the best performance achieved so far within its neighborhood. The neighborhood of each particle is generally defined as topologically nearest particles to the particle at each side.PSO is an evolutionary algorithm simple in concept, easy to implement and computationally efficient. Figures25, 26, and 27 present a generic EC algorithm, a hybrid algorithm, more precisely a MA and the original procedure for implementing the PSO algorithm, respectively.Figure 25
Evolutionary computation algorithm.Figure 26
Memetic algorithm.Figure 27
Particle swarm optimization process.The different versions of the PSO algorithms are the real-value PSO, which is the original version of PSO and is well suited for solving real-value problems; the binary version of PSO, which is designed to solve binary problems; and the discrete version of PSO, which is good for solving the event-based problems. To extend the real-value version of PSO to binary/discrete space, the most critical part is to understand the meaning of concepts such as trajectory and velocity in the binary/discrete space.Kennedy and Eberhart [35] use velocity as a probability to determine whether xid (a bit) will be in one state or another (zero or one). The particle swarm formula of (6.1a) remains unchanged, except that now pid and xid are integers in [0.0,1.0] and a logistic transformation S(vid) is used to accomplish this modification. The resulting change in position is defined by the following rule:
(6.5)if[rand()<S(vid)]thenxid=1;elsexid=0,
where the function S(v) is a sigmoid limiting transformation and rand() is a random number selected from a uniform distribution in the range [0.0,1.0].


### 6.2. PSO Based Circuit Design

We adopt a PSO algorithm to design combinational logic circuits. A truth table specifies the circuits and the goal is to implement a functional circuit with the least possible complexity. Four sets of logic gates have been defined, as shown in Table2, being Gset 2  the simplest one (i.e., a RISC-like set) and Gset 6  the most complex gate set (i.e., a CISC-like set). Logic gate named WIRE means a logical no-operation.Table 2
Gate sets.
Gate setLogic gatesGset 2{AND, XOR, WIRE}Gset 3{AND, OR, XOR, WIRE}Gset 4{AND, OR, XOR, NOT, WIRE}Gset 6{AND, OR, XOR, NOT, NAND, NOR, WIRE}In the PSO scheme the circuits are encoded as a rectangular matrixA (row  ×column=r×c) of logic cells as represented in Figure 28.Figure 28
A3×3 matrix representing a circuit with input X and output Y.Three genes represent each cell: <input1><input2><gate type>, where input1 and input2  are one of the circuit inputs, if they are in the first column, or one of the previous outputs, if they are in other columns. The gate type is one of the elements adopted in the gate set. The chromosome is formed with as many triplets as the matrix size demands (e.g., triplets  = 3×r×c). For example, the chromosome that represents a 3×3 matrix is depicted in Figure 29.Figure 29
Chromosome for the3×3 matrix of  Figure 28.The initial population of circuits (particles) has a random generation. The initial velocity of each particle is initialized with zero. The following velocities are calculated applying (6.2a) and the new positions result from using (6.2b). This way, each potential solution, called particle, flies through the problem space. For each gene is calculated the corresponding velocity. Therefore, the new positions are as many as the number of genes in the chromosome. If the new values of the input genes result out of range, then a re-insertion function is used. If the calculated gate gene is not allowed a new valid one is generated at random. These particles then have memory and each keeps information of its previous best position (pbest) and its corresponding fitness. The swarm has the pbest of all the particles and the particle with the greatest fitness is called the global best (gbest).The basic concept of the PSO technique lies in accelerating each particle towards itspbest  and gbest  locations with a random weighted acceleration. However, in our case we also use a kind of mutation operator that introduces a new cell in 10% of the population. This mutation operator changes the characteristics of a given cell in the matrix. Therefore, the mutation modifies the gate type and the two inputs, meaning that a completely new cell can appear in the chromosome.To run the PSO we have also to define the numberP of individuals to create the initial population of particles. This population is always the same size across the generations, until reaching the solution.The calculation of the fitness functionFs in (6.6) has two parts, f1 and f2, where f1 measures the functionality and f2 measures the simplicity. In a first phase, we compare the output Y produced by the PSO-generated circuit with the required values YR, according with the truth table, on a bit-per-bit basis. By other words, f1 is incremented by one for each correct bit of the output until f1 reaches the maximum value f10 that occurs when we have a functional circuit. Once the circuit is functional, in a second phase, the algorithm tries to generate circuits with the least number of gates. This means that the resulting circuit must have as much genes <gate  type>  ≡ <wire>  as possible. Therefore, the index f2, that measures the simplicity (the number of null operations), is increased by one  (zero) for each wire  (gate) of the generated circuit, yielding
(6.6)f10=2ni×no,f1=f1+1if{bitiofY}={bitiofYR},i=1,…,f10,f2=f2+1ifgatetype=wire,Fs={f1,Fs<f10,f1+f2,Fs≥f10,
where ni and no represent the number of inputs and outputs of the circuit.The concept of dynamic fitness functionFd results from an analogy between control systems and the GA case, where we master the population through the fitness function. The simplest control system is the proportional algorithm; nevertheless, there can be other control algorithms, such as the proportional and the differential scheme.In this line of thought, (6.6) is a static fitness function Fs and corresponds to using a simple proportional algorithm. Therefore, to implement a proportional-derivative evolution the fitness function needs a scheme of the type [18]
(6.7)Fd=Fs+KDμ[Fs],
where 0≤μ≤1 is the differential fractional-order and K∈ℜ is the “gain’’ of the dynamical term.


### 6.3. Experiments and Results

A reliable execution and analysis of an EC algorithm usually requires a large number of simulations to provide a reasonable assurance that the stochastic effects are properly considered. Therefore, in this study are developedn=20 simulations for each case under analysis.The experiments consist on running the three algorithms{GA,MA,PSO} to generate a typical combinational logic circuit, namely, a 2-to-1 multiplexer (M2-1), a 1-bit full adder (FA1), a 4-bit parity checker (PC4) and a 2-bit multiplier (MUL2), using the fitness scheme described in (6.6) and (6.7). The circuits are generated with the gate sets presented in Table 2 and P=3000, w=0.5, c1=1.5, and c2=2.Figure30 depicts the standard deviation of the number of generations to achieve the solution S(N) versus the average number of generations to achieve the solution Av(N) for the algorithms {GA,MA,PSO}, the circuits {M2-1,FA1,PC4,MUL2}, and the gate sets {2,3,4,6}. In these figure, we can see that the MUL2  circuit is the most complex one, while the PC4  and the M2-1  are the simplest circuits. It is also possible to conclude that Gset 6  is the less efficient gate set for all algorithms and circuits.S(N) versus Av(N) with P=3000 and Fs for the GA, the MA, and the PSO algorithms.
(a)(b)(c)Figure30 reveals that the plots follow a power law:
(6.8)S(N)=a[Av(N)]ba,b∈ℜ.
Table 3 presents the numerical values of the parameters (a,b) for the three algorithms.Table 3
The parameters(a,b) and (c,d).
AlgorithmabcdGA0.03651.6020.15261.1734MA0.07281.26020.20891.3587PSO0.26771.15280.01411.1233In terms ofS(N) versus Av(N), the MA algorithm presents the best results for all circuits and gate sets. In what concerns the other two algorithms, the PSO is superior (inferior) to the GA for complex (simple) circuits.Figure31 depicts the average processing time to obtain the solution Av(PT) versus the average number of generations to achieve the solution Av(N) for the algorithms {GA,MA,PSO}, the circuits {M2-1,FA1,PC4,MUL2} and the gate sets {2,3,4,6}. When analysing these charts it is clear that the PSO algorithm demonstrates to be around ten times faster than the MA and the GA algorithms.Av(PT) versus Av(N) with P=3000 and Fs for the GA, the MA, and the PSO algorithms.
(a)(b)(c)These plots follow also a power law:(6.9)Av(PT)=c[Av(N)]dc,d∈ℜ.Table3 shows parameters (c,d) and we can see that the PSO algorithm has the best values.Figures32 and 33 depict the standard deviation of the number of generations to achieve the solution S(N) and the average processing time to obtain the solution Av(PT), respectively, versus the average number of generations to achieve the solution Av(N) for the PSO algorithm using Fd, the circuits {M2-1, FA1, PC4, MUL2}, and the gate sets {2,3,4,6}. We conclude that Fd leads to better results in particular for the MUL2 circuit and for the Av(PT).Figure 32
S(N) versus Av(N) for the PSO algorithm, P=3000 and Fd.Figure 33
Av(PT) versus Av(N) for the GA, P=3000 and Fd.Figures34 and 35 present a comparison between Fs and Fd.Figure 34
Av(N) for the PSO algorithm, P=3000 using Fs and Fd.Figure 35
S(N) for the PSO algorithm, P=3000 using Fs and Fd.In terms ofS(N) versus Av(N) it is possible to say that the MA algorithm presents the best results. Nevertheless, when analysing Figure 31, that shows Av(PT) versus Av(N) for reaching the solutions, we verify that the PSO algorithm is very efficient, in particular, for the more complex circuits.The PSO-based algorithm for the design of combinational circuits follows the same profile as the other two evolutionary techniques presented in this paper.Adopting the study of theS(N) versus Av(N) for the three evolutionary algorithms, the MA algorithm presents better results over the GA and the PSO algorithms. However, in what concerns the processing time to achieve the solutions, the PSO outcomes clearly the GA and the MA algorithms. Moreover, applying the Fd the results obtained are improved further in all gate sets and in particular for the more complex circuits.


## 6.1. Particle Swarm Optimization

In literature about PSO the term ‘swarm intelligence’ appears rather often and, therefore, we begin by explaining why this is so.Noncomputer scientists (ornithologists, biologists, and psychologists) did early research, which led into the theory of particle swarms. In these areas, the term “swarm intelligence’’ is well known and characterizes the case when a large number of individuals are able of accomplish complex tasks. Motivated by these facts, some basic simulations of swarms were abstracted into the mathematical field. The usage of swarms for solving simple tasks in nature became an intriguing idea in algorithmic and function optimization.Eberhart and Kennedy were the first to introduce the PSO algorithm [35], which is an optimization method inspired in the collective intelligence of swarms of biological populations, and was discovered through simplified social model simulation of bird flocking, fishing schooling, and swarm theory.In the PSO, instead of using genetic operators, as in the case of GAs, each particle (individual) adjusts its flying according with its own and its companions experiences. Each particle is treated as a point in a D-dimensional space and is manipulated as described in what follows in the original PSO algorithm:(6.1a)vid=vid+c1rand()(pid-xid)+c2Rand()(pgd-xid),(6.1b)xid=xid+vid,
where c1 and c2 are positive constants, rand() and Rand() are two random functions in the range [0,1], Xi=(xi1,xi2,…,xiD) represents the ith particle, Pi=(pi1,pi2,…,piD)  is the best previous position (the position giving the best fitness value) of the particle, the symbol g represents the index of the best particle among all particles in the population, and Vi=(vi1,vi2,…,viD) is the rate of the position change (velocity) for particle i.However, (6.1a) and (6.1b) represent the flying trajectory of a population of particles. Also, (6.1a) describes how the velocity is dynamically updated and (6.1b) the position update of the “flying” particles. Moreover, (6.1b) is divided in three parts, namely the momentum, the cognitive and the social parts. In the first part the velocity cannot be changed abruptly: it is adjusted based on the current velocity. The second part represents the learning from its own flying experience. The third part consists on the learning group flying experience [36].The first new parameter added into the original PSO algorithm is the inertia weigh. The dynamic equation of PSO with inertia weigh is modified to be(6.2a)vid=wvid+c1rand()(pid-xid)+c2Rand()(pgd-xid),(6.2b)xid=xid+vid,
where w constitutes the inertia weigh that introduces a balance between the global and the local search abilities. A large inertia weigh facilitates a global search while a small inertia weigh facilitates a local search.Another parameter, called constriction coefficientk, is introduced with the hope that it can insure a PSO to converge. A simplified method of incorporating it appears in (6.3), where k is function of c1 and c2 as it is presented as follows:
(6.3)vid=k[vid+c1rand()(pid-xid)+c2Rand()(pgd-xid)],xid=xid+vid,(6.4)k=2(2-ϕ-ϕ2-4ϕ)-1,
where ϕ=c1+c2, ϕ>4.There are two different PSO topologies, namely, the global version and the local version. In the global version of PSO, each particle flies through the search space with a velocity that is dynamically adjusted according to the particle’s personal best performance achieved so far and the best performance achieved so far by all particles. On the other hand, in the local version of PSO, each particle’s velocity is adjusted according to its personal best and the best performance achieved so far within its neighborhood. The neighborhood of each particle is generally defined as topologically nearest particles to the particle at each side.PSO is an evolutionary algorithm simple in concept, easy to implement and computationally efficient. Figures25, 26, and 27 present a generic EC algorithm, a hybrid algorithm, more precisely a MA and the original procedure for implementing the PSO algorithm, respectively.Figure 25
Evolutionary computation algorithm.Figure 26
Memetic algorithm.Figure 27
Particle swarm optimization process.The different versions of the PSO algorithms are the real-value PSO, which is the original version of PSO and is well suited for solving real-value problems; the binary version of PSO, which is designed to solve binary problems; and the discrete version of PSO, which is good for solving the event-based problems. To extend the real-value version of PSO to binary/discrete space, the most critical part is to understand the meaning of concepts such as trajectory and velocity in the binary/discrete space.Kennedy and Eberhart [35] use velocity as a probability to determine whether xid (a bit) will be in one state or another (zero or one). The particle swarm formula of (6.1a) remains unchanged, except that now pid and xid are integers in [0.0,1.0] and a logistic transformation S(vid) is used to accomplish this modification. The resulting change in position is defined by the following rule:
(6.5)if[rand()<S(vid)]thenxid=1;elsexid=0,
where the function S(v) is a sigmoid limiting transformation and rand() is a random number selected from a uniform distribution in the range [0.0,1.0].


## 6.2. PSO Based Circuit Design

We adopt a PSO algorithm to design combinational logic circuits. A truth table specifies the circuits and the goal is to implement a functional circuit with the least possible complexity. Four sets of logic gates have been defined, as shown in Table2, being Gset 2  the simplest one (i.e., a RISC-like set) and Gset 6  the most complex gate set (i.e., a CISC-like set). Logic gate named WIRE means a logical no-operation.Table 2
Gate sets.
Gate setLogic gatesGset 2{AND, XOR, WIRE}Gset 3{AND, OR, XOR, WIRE}Gset 4{AND, OR, XOR, NOT, WIRE}Gset 6{AND, OR, XOR, NOT, NAND, NOR, WIRE}In the PSO scheme the circuits are encoded as a rectangular matrixA (row  ×column=r×c) of logic cells as represented in Figure 28.Figure 28
A3×3 matrix representing a circuit with input X and output Y.Three genes represent each cell: <input1><input2><gate type>, where input1 and input2  are one of the circuit inputs, if they are in the first column, or one of the previous outputs, if they are in other columns. The gate type is one of the elements adopted in the gate set. The chromosome is formed with as many triplets as the matrix size demands (e.g., triplets  = 3×r×c). For example, the chromosome that represents a 3×3 matrix is depicted in Figure 29.Figure 29
Chromosome for the3×3 matrix of  Figure 28.The initial population of circuits (particles) has a random generation. The initial velocity of each particle is initialized with zero. The following velocities are calculated applying (6.2a) and the new positions result from using (6.2b). This way, each potential solution, called particle, flies through the problem space. For each gene is calculated the corresponding velocity. Therefore, the new positions are as many as the number of genes in the chromosome. If the new values of the input genes result out of range, then a re-insertion function is used. If the calculated gate gene is not allowed a new valid one is generated at random. These particles then have memory and each keeps information of its previous best position (pbest) and its corresponding fitness. The swarm has the pbest of all the particles and the particle with the greatest fitness is called the global best (gbest).The basic concept of the PSO technique lies in accelerating each particle towards itspbest  and gbest  locations with a random weighted acceleration. However, in our case we also use a kind of mutation operator that introduces a new cell in 10% of the population. This mutation operator changes the characteristics of a given cell in the matrix. Therefore, the mutation modifies the gate type and the two inputs, meaning that a completely new cell can appear in the chromosome.To run the PSO we have also to define the numberP of individuals to create the initial population of particles. This population is always the same size across the generations, until reaching the solution.The calculation of the fitness functionFs in (6.6) has two parts, f1 and f2, where f1 measures the functionality and f2 measures the simplicity. In a first phase, we compare the output Y produced by the PSO-generated circuit with the required values YR, according with the truth table, on a bit-per-bit basis. By other words, f1 is incremented by one for each correct bit of the output until f1 reaches the maximum value f10 that occurs when we have a functional circuit. Once the circuit is functional, in a second phase, the algorithm tries to generate circuits with the least number of gates. This means that the resulting circuit must have as much genes <gate  type>  ≡ <wire>  as possible. Therefore, the index f2, that measures the simplicity (the number of null operations), is increased by one  (zero) for each wire  (gate) of the generated circuit, yielding
(6.6)f10=2ni×no,f1=f1+1if{bitiofY}={bitiofYR},i=1,…,f10,f2=f2+1ifgatetype=wire,Fs={f1,Fs<f10,f1+f2,Fs≥f10,
where ni and no represent the number of inputs and outputs of the circuit.The concept of dynamic fitness functionFd results from an analogy between control systems and the GA case, where we master the population through the fitness function. The simplest control system is the proportional algorithm; nevertheless, there can be other control algorithms, such as the proportional and the differential scheme.In this line of thought, (6.6) is a static fitness function Fs and corresponds to using a simple proportional algorithm. Therefore, to implement a proportional-derivative evolution the fitness function needs a scheme of the type [18]
(6.7)Fd=Fs+KDμ[Fs],
where 0≤μ≤1 is the differential fractional-order and K∈ℜ is the “gain’’ of the dynamical term.


## 6.3. Experiments and Results

A reliable execution and analysis of an EC algorithm usually requires a large number of simulations to provide a reasonable assurance that the stochastic effects are properly considered. Therefore, in this study are developedn=20 simulations for each case under analysis.The experiments consist on running the three algorithms{GA,MA,PSO} to generate a typical combinational logic circuit, namely, a 2-to-1 multiplexer (M2-1), a 1-bit full adder (FA1), a 4-bit parity checker (PC4) and a 2-bit multiplier (MUL2), using the fitness scheme described in (6.6) and (6.7). The circuits are generated with the gate sets presented in Table 2 and P=3000, w=0.5, c1=1.5, and c2=2.Figure30 depicts the standard deviation of the number of generations to achieve the solution S(N) versus the average number of generations to achieve the solution Av(N) for the algorithms {GA,MA,PSO}, the circuits {M2-1,FA1,PC4,MUL2}, and the gate sets {2,3,4,6}. In these figure, we can see that the MUL2  circuit is the most complex one, while the PC4  and the M2-1  are the simplest circuits. It is also possible to conclude that Gset 6  is the less efficient gate set for all algorithms and circuits.S(N) versus Av(N) with P=3000 and Fs for the GA, the MA, and the PSO algorithms.
(a)(b)(c)Figure30 reveals that the plots follow a power law:
(6.8)S(N)=a[Av(N)]ba,b∈ℜ.
Table 3 presents the numerical values of the parameters (a,b) for the three algorithms.Table 3
The parameters(a,b) and (c,d).
AlgorithmabcdGA0.03651.6020.15261.1734MA0.07281.26020.20891.3587PSO0.26771.15280.01411.1233In terms ofS(N) versus Av(N), the MA algorithm presents the best results for all circuits and gate sets. In what concerns the other two algorithms, the PSO is superior (inferior) to the GA for complex (simple) circuits.Figure31 depicts the average processing time to obtain the solution Av(PT) versus the average number of generations to achieve the solution Av(N) for the algorithms {GA,MA,PSO}, the circuits {M2-1,FA1,PC4,MUL2} and the gate sets {2,3,4,6}. When analysing these charts it is clear that the PSO algorithm demonstrates to be around ten times faster than the MA and the GA algorithms.Av(PT) versus Av(N) with P=3000 and Fs for the GA, the MA, and the PSO algorithms.
(a)(b)(c)These plots follow also a power law:(6.9)Av(PT)=c[Av(N)]dc,d∈ℜ.Table3 shows parameters (c,d) and we can see that the PSO algorithm has the best values.Figures32 and 33 depict the standard deviation of the number of generations to achieve the solution S(N) and the average processing time to obtain the solution Av(PT), respectively, versus the average number of generations to achieve the solution Av(N) for the PSO algorithm using Fd, the circuits {M2-1, FA1, PC4, MUL2}, and the gate sets {2,3,4,6}. We conclude that Fd leads to better results in particular for the MUL2 circuit and for the Av(PT).Figure 32
S(N) versus Av(N) for the PSO algorithm, P=3000 and Fd.Figure 33
Av(PT) versus Av(N) for the GA, P=3000 and Fd.Figures34 and 35 present a comparison between Fs and Fd.Figure 34
Av(N) for the PSO algorithm, P=3000 using Fs and Fd.Figure 35
S(N) for the PSO algorithm, P=3000 using Fs and Fd.In terms ofS(N) versus Av(N) it is possible to say that the MA algorithm presents the best results. Nevertheless, when analysing Figure 31, that shows Av(PT) versus Av(N) for reaching the solutions, we verify that the PSO algorithm is very efficient, in particular, for the more complex circuits.The PSO-based algorithm for the design of combinational circuits follows the same profile as the other two evolutionary techniques presented in this paper.Adopting the study of theS(N) versus Av(N) for the three evolutionary algorithms, the MA algorithm presents better results over the GA and the PSO algorithms. However, in what concerns the processing time to achieve the solutions, the PSO outcomes clearly the GA and the MA algorithms. Moreover, applying the Fd the results obtained are improved further in all gate sets and in particular for the more complex circuits.


## 7. Conclusions

Fractional Calculus (FC) goes back to the beginning of the theory of differential calculus. Nevertheless, the application of FC just emerged in the last two decades, due to the progress in the area of chaos that revealed subtle relationships with the FC concepts.Recently FC has been a fruitful field of research in science and engineering and many scientific areas are currently paying wider attention to the FC concepts. In the field of dynamical systems theory, some work has been carried out but the proposed models and algorithms are still in a preliminary stage of establishment. This article presented several case studies on the implementation of FC-based models and control systems, being demonstrated the advantages of using the FC theory in different areas of science and engineering. In fact, this paper studied a variety of different physical systems, namely(i)
tuning of PID controllers using fractional calculus concepts;(ii)
fractional PDα control of a hexapod robot;(iii)
fractional dynamics in the trajectory control of redundant manipulators;(iv)
heat diffusion;(v)
circuit synthesis using evolutionary algorithms.It has been recognized the advantageous use of this mathematical tool in the modeling and control of these dynamical systems, and the results demonstrate the importance of Fractional Calculus and motivate for the development of new applications.



---

*Source: 639801-2009-11-08.xml*