# Nonlinear Nonconvex Optimization by Evolutionary Algorithms Applied to Robust Control

Authors: Joaquín Cervera; Alfonso Baños
Journal: Mathematical Problems in Engineering
 (2009)

---

## Abstract

This work focuses on the problem of automatic loop shaping in the context of robust control. More specifically in the framework given by Quantitative Feedback Theory (QFT), traditionally the search of an optimum design, a non convex and nonlinear optimization problem, is simplified by linearizing and/or convexifying the problem. In this work, the authors propose a suboptimal solution using a fixed structure in the compensator and evolutionary optimization. The main idea in relation to previous work consists of the study of the use of fractional compensators, which give singular properties to automatically shape the open loop gain function with a minimum set of parameters, which is crucial for the success of evolutionary algorithms. Additional heuristics are proposed in order to guide evolutionary process towards close to optimum solutions, focusing on local optima avoidance.


---

## Body

## 1. Introduction

It is a well-known fact that there is no general procedure to exactly solve nonlinear nonconvex optimization problems when the solutions belong to continuous solution sets ([1]). In the case of solutions belonging to discrete solution sets, it is always possible to find the optimal solution by a branch and bound type exhaustive search ([2]). Branch and bound techniques can also be applied to continuous solution set problems, specially when combined with interval analysis ([3, 4]). The obtained solution is, in general, only close to optimal, according to a certain accuracy factor. A major drawback of branch and bound techniques is that they are a very costly approach in terms of computation time. The only way to quickly obtain an (approximate) solution to this sort of problem is to use some kind of randomized search algorithm, like random algorithms ([5]) or evolutionary algorithms ([6]) according to the particular problem to be solved.Quantitative Feedback Theory is a robust frequency domain control design methodology which has been successfully applied in practical problems from different domains [7]. One of the key design steps in QFT, equivalent to controller design, is loop shaping of the open loop gain function to a set of restrictions (or boundaries) given by the design specifications and the (uncertain) model of the plant. Although this step has been traditionally performed by hand, the use of CACSD tools (e.g., the QFT MATLAB Toolbox [8]) has made the manual loop shaping much more simple. However, the problem of automatic loop shaping is of enormous interest in practice, since the manual loop shaping can be hard for the nonexperienced engineer, and thus it has received a considerable attention, specially in the last three decades.Optimal QFT loop computation is a nonlinear nonconvex optimization problem, for which there is not yet an optimization algorithm which computes a globally optimum solution in a reasonable time, in terms of interactive design purposes. It must be noticed, however, that the work by Nataraj and others on this subject, based on deterministic optimization procedures, combining branch and bound optimization and interval analysis techniques, is very promising (see e.g., [9, 10]).Other typical approaches to solve this problem have tried to find approximate solutions in different ways. For instance, some authors have simplified the problem somehow, in order to obtain a different optimization problem for which there exists a closed solution or an optimization algorithm which does guarantee a global optimum in a shorter computation time. A trade-off between necessarily conservative simplification of the problem and computational solvability has to be chosen. This is the approach in, for instance, [11, 12]. Some authors have investigated the loop shaping problem in terms of particular structures, with a certain degree of freedom, which can be shaped to the particular problem to be solved. This is the case in [13–15]. Another possibility is to use randomized search algorithms, able to directly face nonlinear and nonconvex optimization problems, at the cost of returning an only close to optimal solution, and not guaranteeing that, in general, the returned solution is not far from the optimum. This is the approach adopted in [16, 17], based on evolutionary algorithms.Evolutionary algorithms are computationally demanding, specially as the dimension of the search space increases. In this paper the authors study the use of evolutionary algorithms-based optimization, proposing the addition, with respect to previous work, of some heuristics, very much specific to the particular problem under consideration, which help to improve obtained solutions accuracy and computation time needed to obtain these solutions. In this sense, a good structure for the compensator, in terms of using a reduced set of parameters, but with a rich frequency domain behavior, is of crucial importance. This is the main heuristic proposed in this paper: to use evolutionary algorithms together with a flexible structure, able to get a close to optimum solution, but with a reduced number of parameters. In previous work, the compensator has been fixed to a rational structure, with a finite (but no necessarily small) number of zeros and poles. In this work, the main contribution is to introduce a fractional compensator that, with a minimum number of parameters, gives a flexible structure in the frequency domain regarding automatic loop shaping. In fact, it can be approximated by a rational compensator, but with a considerably large number of parameters. This dramatic reduction in the number of parameters has shown to be of capital importance for the success of evolutionary algorithms in the solution of the automatic loop shaping problem. Other applied heuristics have to do with including some features in the objective function that guide the evolutionary search towards close to optimum solutions, paying special attention to prevent the search from getting stacked in local minima, which is specially likely to happen in the problem under consideration.In this work, following [18–21] it is considered the particular case of minimum phase open loop gain functions, for which the investigated compensators can give a good structure with a reduced set of parameters. Some first steps towards the nonminimum phase case are given in [18].From here onwards, the structure of the paper stands as follows. In Section2 a brief introduction to QFT is given; in Section 3 the proposed solution for the QFT automatic loop shaping problem by evolutionary algorithms is presented; in Section 4 this procedure is applied to a typical benchmark problem. Finally, Section 5 presents the conclusions.


## 2. Introduction to QFT

The basic idea in QFT ([7]) is to define and take into account, along the control design process, the quantitative relationship between the amount of uncertainty to deal with and the amount of control effort to use. Typically, the QFT control system configuration (see Figure 1) considers two degrees of freedom: a controller C(s), in the closed loop, which manages uncertainty and disturbances affecting u(t) or y(t); and a precompensator, F(s), designed after C(s), used to satisfy tracking specifications.Figure 1
Two degrees of freedom control system configuration.Consider an uncertain plant, represented as the set of transfer functionsP={P(s)=P0(s)Q(s),Q∈Q}, with P0(s) being the nominal plant and Q being a set of transfer functions representing uncertainty. For a certain frequency ω, the templateTP(jω) is defined as the Nichols chart (NC) region ([22]) given by
(2.1)TP(ω)=̇{(∠P(jω),|P(jω)|dB)∈NC,P∈P}.The design of the controllerC(s) is accomplished in the Nichols chart, by shaping the nominal open loop transfer function, L0(s)=P0(s)C(s). A discrete set of design frequencies Ω is chosen. Given quantitative specifications on robust stability and robust performance on the closed loop system, the set of boundariesB={Bω,ω∈Ω} is computed. Each Bω defines the allowed region Aω (conversely, forbidden region Fω=NC∖Aω) in NC for L0(jω), so that L0(jω)∈Aω (conversely, L0(jω)∉Fω) forallω∈Ω (boundaries satisfaction) implies specification satisfaction by L(s)=P(s)C(s)forallP(s)∈P. For example, in Figure 9 we have the following.(i)
The open lines horizontally crossing the Nichols plane from-360° to 0° are performance boundaries. In this example, for each ω∈Ω={0.1,0.5,2,15,100}, A(jω) is defined as the region  above Bω.(ii)
The closed line around the point (0 dB,-180°) is a robust stability boundary, defining A(jω) (conversely F(jω)) as the region outside (conversely inside) Bω. In this example, this boundary is a Universal High-Frequency Boundary (UHFB), a stability boundary called universal because the region inside it, FUHFB, is not forbidden only for a particular frequency ω, but for all ω∈ℝ+, so it imposes L0(jω)∉FUHFB for all ω∈ℝ+.The basic step in the design process,loop shaping, consists of the following nonlinear nonconvex constrained optimization problem: design of L0(s) which satisfies boundaries (constraints) and minimizes the high-frequency gain (optimization criterion—[23]), defined as
(2.2)Khf=̇lims→0seL0(s),
where e is L0 excess of poles over zeros. The comparison in terms of this optimization criterion of two nominal open loops L01(s) and L02(s), with respective excess of poles over zeros e1 and e2, only makes sense if e1=e2. The comparison is performed in terms of the loops high-frequency gains, K01 and K02, respectively.Note that this criterion is defined irrespective of the particular structure used forL0(s). It has been shown ([23]) that for minimum phase L0(s) this optimum exists, it is unique and, in general, has an infinite number of poles and zeros, and so it is not implementable. For this reason it is referred along this work as the QFT theoretical optimum, Lqto(s), with optimal high-frequency gain Kqto.SinceLqto(s) is not practical as a loop design, a more practical definition of optimum can be stated. Define X=̇{x=(x1,…,xn)}, xi∈ℝ, i=1⋯n, as the set of controller structure parameters to be instantiated. Consider LX(s) an open loop structure parameterized by X, that is, each x∈X defines Lx(s), a particular instance of LX(s), with high-frequency gain KLx. Define KLO=̇min{KLx,x∈X}≤Kqto, the minimum high-frequency gain that can be achieved by any LX(s) instance. An L-optimal transfer function or LO is defined as an LX(s) instance La such that KLa=KLO.For a givenLX(s) open loop structure, KLO represents the best approach to Kopt that can be achieved with that structure. The ability of a structure LX(s) to achieve Lqto(s) shape is called close-to-optimality, defined as CO(L)=̇KLO-Kqto. Since computing LO is in general still a hard optimization problem, and thus KLO is not known, it is common to compute instead suboptimal transfer functions LsOi with high-frequency gain KLsOi as close to KLO as possible, with i∈ℕ. This is the approach considered in this work, using evolutionary algorithms to compute an LsOi with KLsOi which approaches KLO as much as possible. A measure of how accurately a certain LsO1 approaches KLO can only be established in relative terms, that is, by comparing it to another LsO2 in terms of their high-frequency gains.The meaning of this optimum definition is minimizing the cost of feedback related with sensor noise amplification at high frequencies by minimizing|L0(jω)| at high frequencies, which, due to phase/magnitude relations given by Bode integrals ([24]), is equivalent to minimize ∫ℝ+∠L0(jω)dω, subject to B satisfaction. Figures 2 and 3 show the typical shape of a QFT optimal loop Lqto(jω), with polygonal line simplified versions of Lqto(jω) and boundaries in B. For frequencies ω∈[0,ωa2], ∠Lqto(jω) minimization is limited by performance boundaries. For frequencies ω∈[ωa2,∞], ∠Lqto(jω) minimization is only limited by the UHFB (which makes Lqto(jω) tightly follow the right and bottom sides of the UHFB) and the specified excess of poles over zeros e. Note that ωc is not directly established in the specifications, but indirectly, as a consequence of boundaries in B.Figure 2
Nichols plot ofLqto(jω), QFT theoretical optimal loop, exhibiting typical Bode Step in the frequency range [ωM1,ω2].Bode plot ofLqto(jω), QFT theoretical optimal loop, exhibiting typical Bode Step in the frequency range [ωM1,ω2].
(a)(b)Note the shape ofLqto(jω) in the frequency range [ωM1,ω2]; see Figures 2 and 3. It is called Bode step ([18, 23, 25], first defined in [25]), a frequency band where Lqto(jω) has constant magnitude and a fast phase decrease (according to ∠Lqto(jω) minimization objective, once the UHFB does not constrain this objective), surrounded by frequencies where ∠Lqto(jω) is constant (απ and -eπ/2, resp.) and |Lqto(jω)| decreases. Bode step is very typical of QFT optimal loops, and when designing a suboptimal loop LsO, the designer has to try to make LsO exhibit a Bode step-like shape in order to succeed in accurately approaching Lqto. Note that a Bode step in LsO(jω) can be interpreted as LsO(jω) tightly following the bottom part of the UHFB.


## 3. QFT Automatic Controller Design by Evolutionary Optimization

The used method for automatic QFT controller design consists of using a fixed structure in the compensator with a certain numbern of free parameters xi∈ℝ, i=1⋯n, which constitute the solution space X=̇{x=(x1,…,xn)} of the nonconvex nonlinear global optimization to be performed by evolutionary algorithms. As it was said in the introduction, there are two critical factors which determine the efficiency of such an approach. One is the dimension of the search space, n, which exponentially increases the computational cost in terms of time. The second factor is the use of adequate heuristics which guide the evolutionary search towards close to optimum solutions.The main contribution of this work has to do with the first factor: to use a fractional structure in the compensator is proposed as a key idea to get flexible structures, able to yield close to optimum solutions, but with a reducedn. Several structures from literature have been adapted to solve the QFT loop shaping problem, including some structures proposed by the authors in previous works. These structures are introduced in Sections 3.1.1–3.1.4.The second contribution of this work has to do with the second factor. Some ad hoc heuristics, with features very much specific to the particular problem under consideration, have been developed in order to help evolutionary search to improve obtained solutions accuracy and computation time needed to obtain these solutions, specially in terms of local minima avoidance. These heuristics are presented in Section3.2. Section 3.3 describes the algorithm used for evolutionary optimization, detailing the objective function, which includes these heuristics.In Section4 a design example is solved by using these heuristics and all the fractional structures. The results obtained by each fractional structure are compared.


### 3.1. Fractional Structures

#### 3.1.1. TID

TID controller [26] is a modified version of PID controller, where the proportional term is replaced by a tilted (fractional) term, with transfer function s-eT, which permits a better approach to theoretical optimum. The resulting controller, including a low pass filter, is given by
(3.1)CTID(s)=k(TseT+1s+qDsq+s)1(1+s/wh)nh.


#### 3.1.2.PIλDμ

PIλDμ controller [27, 28] is a PID generalization in which both integrator and derivative terms have real order, λ and μ, respectively. In this work, the multiplicative version of PIλDμ given in [29] will be used, with transfer function
(3.2)CPIλDμ(s)=kcxμ(λ1s+1s)λ(λ2s+1xλ2s+1)μ.
CRONE 2 (Section 3.1.3) is a particular case of (3.2), where some parameters are linked, and so flexibility is reduced compared to (3.2).


#### 3.1.3. CRONE-Based Controllers

The CRONE approach [30, 31] defines three generations of fractional controllers based on the use of frequency-band noninteger differentiators. First and second generations (CRONE 1 and CRONE 2) use real non-integer differentiation, whereas third generation (CRONE 3) use complex non-integer differentiation. Three CRONE structures-based compensators are studied. The first one uses the transfer function structure common to CRONE 1 and 2. The basic component of this structure is an order n differentiator in the form of the implementable band-defined transfer function. An order nI band-limited integrator and an order nF lowpass filter are added to manage accuracy and robustness and control effort problems, being the open loop structure finally defined as
(3.3)LCR2(s)=k(ωls+1)nI(1+s/ωh1+s/ωl)n1(s/ωh+1)nF.
The application of the CRONE 2 structure to the QFT problem is quite straightforward.The second compensator uses CRONE 3 structure, consisting of the substitution of the (real) ordern integrodifferentiator in CRONE 2 for the real part Dr(s) of a (complex) order n=a+ib integrodifferentiator in the form of the implementable band-defined transfer function
(3.4)D(s)=Dr(s)+iDi(s)=(C01+s/ωh1+s/ωl)a+ib
with C0=ωh/ωu=ωu/ωl and ωu=ωlωh. The open loop structure L(s) in CRONE 3 is finally defined  as
(3.5)LCR3(s)=k(ωls+1)nI(C01+s/ωh1+s/ωl)a×cos[bLog(C01+s/ωh1+s/ωl)]1(s/ωh+1)nF.The third compensator is decoupled CRONE 3 ([32]), a modified version of CRONE 3 structure (3.5), where some parameters are decoupled in order to obtain higher flexibility. Frequencies ωh and ωl are decoupled by defining new frequencies ωh′, ωl′, and ωh4. C0 is decoupled by redefining C0 in cos() as C0′=cC0, where c is a new free parameter, C0=ωh′/ωu=ωu/ωl′, and ωu=ωl′ωh′. The new structure to be shaped  is
(3.6)LdecCR3(s)=k(ωls+1)nI(C01+s/ωh1+s/ωl)a×cos[bLog(cC01+s/ωh′1+s/ωl′)]1(s/ωh4+1)nF.
For both CRONE 3 and decoupled CRONE 3, a detailed study of the parameters involved, its interrelations and their allowed ranks in order to obtain desired behavior, is developed in [32]. In [33] a CRONE-based structure is proposed for the design of open loops based on Bode optimum-like specifications.


#### 3.1.4. FCT Terms

Fractional order Complex Terms (FCTs) [19, 34] are terms
(3.7)Ti=(ωni2s2+2δiωnis+ωni2)ei
with ei∈ℝ, ei>0 corresponding to poles and ei>0 corresponding to zeros. Different combinations of these terms can be used depending on the problem to be solved. For a typical tracking problem with a minimum phase plant, the structure
(3.8)L0=KT1T2T3
with K∈ℝ+, two FCT poles and one FCT zero, achieves very good results.


### 3.2. Heuristics

Once the problem of a large number of parameters to be optimized has been solved by the use of fractional structures, the main problem is the fact that, due to the typically convex nature of the constraints in the solutions space, the evolutionary search can easily get stacked in local minima. The main goal for the heuristics design has been to help the evolutionary search to quickly avoid these local minima when they are detected.The objective function used as the criterion for the natural selection of individuals during the evolutionary optimization process,O(X), is a multiobjective real function which returns a to-be-minimized real value related with nonviolation of boundaries (constraint), stability (constraint), and minimization of Khf (to-be-optimized value). Violation of boundaries and lack of stability are, in fact, constraints, but due to how evolutionary algorithms work, they have to be translated into to-be-optimized values, as Khf originally is. A first approach could be to assign a large value to any solution which violates any boundary or is unstable. But this does not differentiate between solutions which violate a certain boundary completely, or only a little bit, which produces a blind O(X) in the sense that is not able to look for the way to become out of the violation of a certain given constraint. In order to avoid the evolutionary optimization getting stack in such situations, it is necessary to guide it by establishing, by adequately shaping O(X), that the violation of a certain constraint by a certain solution Xa is not as important as the violation produced by another solution Xb, so O(Xa)<O(Xb). This way, the evolutionary algorithm will prefer the survival of Xa instead of Xb, so that better solutions survive, even if they do not satisfy constraints. This scheme, repeated generation after generation, leads to solutions that do respect constraints. These are some of the components of O(X) intended to guide the evolutionary optimization towards solutions satisfying constraints.(i)
For open boundaries, letBω be an open boundary, and consider that it is defined as Bω(p):ℝ→ℝ, a function that assigns to every phase p in NC the magnitude of Bω at that phase. For certain solutions Xa and Xb, it should happen |LXi(jω)|≥Bω(∠LXi(jω)), i={a,b}. Assume that this constrain is not satisfied; that is, Bω is violated, by both Xa and Xb. Xa is considered better than Xb when LXa(jω) is closer to Bω than LXb is, in terms of magnitude. More formally, the penalty for this violation included in O(X) is a monotonically increasing function of Bω(∠LX(jω))-|LX(jω)|. This choice, made generation after generation, contributes to get solutions X such that LX(jω) is over Bω, thus satisfying Bω constraint.(ii)
For closed boundaries, letBω be a closed boundary, and consider that it is defined by Biω(p):ℝ→ℝ, i=u,l, two functions that assign to every phase p in NC the upper and lower magnitudes of Bω at that phase, respectively (bivalued boundary). For certain solutions Xa and Xb, it should happen |LXi(jω)|≥Buω(∠LXi(jω)) and |LXi(jω)|≤Blω(∠LXi(jω)), i={a,b}. Assume that this constrain is not satisfied; that is, Bω is violated, by both Xa and Xb. Xa is considered better than Xb when LXa(jω) is closer to Buω or Blω than LXb is, in terms of NC Euclidean distance. This choice, made generation after generation, contributes to get solutions X such that LX(jω) is out of Bω, thus satisfying Bω constraint.Another important consideration, in order to avoid the evolutionary search getting stacked, is anO(X) which gives more importance to constraints satisfaction than to khf minimization. This way, the evolutionary search first searches for valid solutions, and then, once they have been achieved, it optimizes among them.


### 3.3. Algorithm

This section describes the optimization algorithm which has been implemented for controllers synthesis. It consists on the use of commercial evolutionary algorithm software (theGenetic and evolutionary algorithm toolbox for use with MATLAB (GEATbx, [35])) together with an ad hoc objective function, programmed by the authors, implementing the heuristics described in Section 3.2.The evolutionary algorithm is, in particular, a multiple subpopulations evolutionary search algorithm. In this kind of evolutionary search, each subpopulation evolves in an isolated way for a few generations (as it happens in a single population evolutionary algorithm). After that, one or more individuals are exchanged between subpopulations. The way this process models species evolution is more similar to nature, compared to single population evolutionary algorithms, which helps to avoid local minima. The basic structure of this kind of algorithm is the following (adapted from [35]):PROCEDURE Multipopulation_Evolutionary_Search_Algorithm (search_space)BEGIN PROCEDUREINITIALIZATION, consisting on:creation of initial population // set of randomly chosen// points from search space,// grouped in speciesevaluation_of_individuals =OBJECTIVE_FUNCTION(initial_population);WHILE NOT(any termination criteria is met)BEGIN WHILENEW POPULATION GENERATION: // consisting on...fitness assignment selection;recombination;mutation;evaluation_of_offspring = OBJECTIVE_FUNCTION(offspring);reinsertion;migration;competition;END WHILE;RESULT = best individual in last population;END PROCEDURE;This is the code which implements the objective functionO(X) described in Section 3.2:PROCEDURE Objective_Function (set_of_individuals)BEGIN PROCEDUREFOR EACH individual IN set_of_individuals, indexed as iSBV = Compute_Stability_Boundary_Violation(individual);SRV = Compute_Stability_Ray_Violation(individual);PBV = Compute_Performance_Boundary_Violation(individual);NBV = Compute_Noise_Boundary_Violation(individual);PSP = Compute_Phases_Separation_Penalization(individual);MPP = Compute_Maximum_Phase_Penalization(individual);HFG = Compute_High_Frequency_Gain(individual);Constraints_Violation =BIG_VALUE * (weights[1] * SBV+ weights[2] * SRV+ weights[3] * PBV+ weights[4] * NBV);Technical_Penalizations =BIG_VALUE * (weights[5] * PSP+ weights[6] * MPP);Optimization_Criterion =weights[7] * HFG;RESULT[i] = Constraints_Violation+ Technical_Penalizations+ Optimization_Criterion;END FOR;END PROCEDURE;ValuesSBV, SRV, PBV and NBV are related with a direct violation of any of the problem constraints. This is the reason they are amplified by  BIG_VALUE  so that any solution satisfying constraints is preferred (has a lower objective value) compared to any other which does not satisfy any of the constraints.ValuesPSP and MPP are related with technical issues, which are necessary to consider in order to avoid nonreal solutions. Nonreal solutions are solutions that, due to intrinsic limitations of discrete computation, such as the need to consider a finite number of points to represent boundaries and loops, would satisfy the algorithm discrete translation of continuous constraints, but not the original continuous constraints. These values are also amplified by BIG_VALUE.ValueHFG is directly related with the optimization criterion, high-frequency gain minimization. It is not amplified by BIG_VALUE, and so it is only considered, in practice, when no constraint is violated.Weights are conceived to give more importance to a certain penalization compared to others, but have not been used for the moment; that is, all of them are equal.As explained in Section3.2, SBV, SRV, PBV, and NBV values correspond to a gradient, related with how much a certain boundary is violated, so that the objective function distinguishes different degrees of violation of a certain constraint, so that it helps the evolutionary algorithm to approach solutions which do not violate that constraint. For instance, consider the following MATLAB code, corresponding to the implementation of the function Compute Performance Boundary Violation (CPBV).PROCEDURE Compute_Perform_Boundary_Violation(loop_points,boundaries)loop_phases = loop_points.phases;loop_magnitudes = loop_points.magnitudes;bnd_phases = boundaries.phases;bnd_magnitudes = boundaries.magnitudes;nphases = length(bnd_phases);bnd_phases_indexes =mod(round(loop_phases*(nphases/360)+nphases),nphases);bnd_phases_indexes = nphases - bnd_phases_indexes;minimum_allowed_magnitudes = bnd_magnitudes(bnd_phases_indexes);magnitude_differences =loop_magnitudes - minimum_allowed_magnitudes;negative_differences =magnitude_differences .* (magnitude_differences < 0);squared_differences = negative_differences.̂2;RESULT = sum(squared_differences);END PROCEDURE;For each design frequencyω, this procedure checks whether the loop is below the performance boundary at that frequency, Bω, and in that case quantifies how much below it is, which is called magnitude difference in the code. These magnitude differences are squared, so that bigger magnitude differences are considered much worse than small differences. After that, these values are summed up to obtain the final boundary violation indicator, so that this indicator doubles when the loop violates boundaries at two frequencies, compared to a single violation, is triple when there are three violating frequencies, and so on. For an example of CPBV function behaviour, consider Figures 4, 5, and 6, representing loops Lex1, Lex2, and Lex3, respectively, corresponding to the first loop design in Section 4 (based on a PID controller). In Figure 4, performance boundaries are violated by Lex1 at design frequencies w1=0.1rad/s and w2=0.5 rad/s, by about 5 dB. The function yields a boundary violation indicator CPBV(Lex1)=4.1e3. Using BIG_VALUE  = 1e6, this produces an objective function value around O(Lex1)=4.1e9, which is a big penalization. The evolutionary algorithm will prefer individuals (loops) with lower objective function values. For instance, in Figure 5 performance boundaries are violated only at design frequency w1=0.1 rad/s, by 1 or 2 dB. CPBV function yields a boundary violation index CPBV(Lex2)=1.5e3, which corresponds to an objective function value O(Lex2)=1.5e9, still a big penalization, but lower than O(Lex1). This way, the evolutionary algorithm tends to move the loop, at each frequency ω∈Ω, out of F(jω) (forbidden area) and towards A(jω) (allowed area), producing a design like Lex3 in Figure 6, where no boundary is violated, and so CPBV(Lex3)=0 and O(Lex3) is contributed only by the optimization criterion, with no penalization by constraints violation, yielding O(Lex3)=152.Figure 4
Lex1 loop, with CPBV(Lex1)=4.1e3 and O(Lex1)=4.1e9.Figure 5
Lex2 loop, with CPBV(Lex2)=1.5e3 and O(Lex2)=1.5e9.Figure 6
Lex3 loop, with CPBV(Lex3)=0 and O(Lex3)=152.Figure7 shows the evolution, for an optimization of the loop structure shown in Figures 4, 5 and 6, of O(Li), where i is an index for each  generation in the evolutionary algorithm execution, and Li is the best individual (loop) in generation i, best in the sense of minimizing the objective function.   During the first three generations there are some boundary violations, which produce big objective function values. These values decrease as generation index i increases, due to the effect of the implemented heuristics. From fourth generation there is no boundary violation, so BIG_VALUE  is not applied any more, and so from that point the task of the evolutionary algorithm is reduced to get better and better values for the optimization criterion. This process can be better visualized in Figure 8, where scale has been adapted for this purpose. This figure permits checking how the objective function converges to a certain value, in this case around 152, from a certain generation, in this case around generation 130.Figure 7
Evolution ofO(Li) with generation index i, big scale, showing the effect of boundaries being violated at first generations.Figure 8
Evolution ofO(Li) with generation index i, reduced scale, showing the minimization of optimization criterion, once constraints are satisfied.Figure 9
PID-based loop design.


## 3.1. Fractional Structures

### 3.1.1. TID

TID controller [26] is a modified version of PID controller, where the proportional term is replaced by a tilted (fractional) term, with transfer function s-eT, which permits a better approach to theoretical optimum. The resulting controller, including a low pass filter, is given by
(3.1)CTID(s)=k(TseT+1s+qDsq+s)1(1+s/wh)nh.


### 3.1.2.PIλDμ

PIλDμ controller [27, 28] is a PID generalization in which both integrator and derivative terms have real order, λ and μ, respectively. In this work, the multiplicative version of PIλDμ given in [29] will be used, with transfer function
(3.2)CPIλDμ(s)=kcxμ(λ1s+1s)λ(λ2s+1xλ2s+1)μ.
CRONE 2 (Section 3.1.3) is a particular case of (3.2), where some parameters are linked, and so flexibility is reduced compared to (3.2).


### 3.1.3. CRONE-Based Controllers

The CRONE approach [30, 31] defines three generations of fractional controllers based on the use of frequency-band noninteger differentiators. First and second generations (CRONE 1 and CRONE 2) use real non-integer differentiation, whereas third generation (CRONE 3) use complex non-integer differentiation. Three CRONE structures-based compensators are studied. The first one uses the transfer function structure common to CRONE 1 and 2. The basic component of this structure is an order n differentiator in the form of the implementable band-defined transfer function. An order nI band-limited integrator and an order nF lowpass filter are added to manage accuracy and robustness and control effort problems, being the open loop structure finally defined as
(3.3)LCR2(s)=k(ωls+1)nI(1+s/ωh1+s/ωl)n1(s/ωh+1)nF.
The application of the CRONE 2 structure to the QFT problem is quite straightforward.The second compensator uses CRONE 3 structure, consisting of the substitution of the (real) ordern integrodifferentiator in CRONE 2 for the real part Dr(s) of a (complex) order n=a+ib integrodifferentiator in the form of the implementable band-defined transfer function
(3.4)D(s)=Dr(s)+iDi(s)=(C01+s/ωh1+s/ωl)a+ib
with C0=ωh/ωu=ωu/ωl and ωu=ωlωh. The open loop structure L(s) in CRONE 3 is finally defined  as
(3.5)LCR3(s)=k(ωls+1)nI(C01+s/ωh1+s/ωl)a×cos[bLog(C01+s/ωh1+s/ωl)]1(s/ωh+1)nF.The third compensator is decoupled CRONE 3 ([32]), a modified version of CRONE 3 structure (3.5), where some parameters are decoupled in order to obtain higher flexibility. Frequencies ωh and ωl are decoupled by defining new frequencies ωh′, ωl′, and ωh4. C0 is decoupled by redefining C0 in cos() as C0′=cC0, where c is a new free parameter, C0=ωh′/ωu=ωu/ωl′, and ωu=ωl′ωh′. The new structure to be shaped  is
(3.6)LdecCR3(s)=k(ωls+1)nI(C01+s/ωh1+s/ωl)a×cos[bLog(cC01+s/ωh′1+s/ωl′)]1(s/ωh4+1)nF.
For both CRONE 3 and decoupled CRONE 3, a detailed study of the parameters involved, its interrelations and their allowed ranks in order to obtain desired behavior, is developed in [32]. In [33] a CRONE-based structure is proposed for the design of open loops based on Bode optimum-like specifications.


### 3.1.4. FCT Terms

Fractional order Complex Terms (FCTs) [19, 34] are terms
(3.7)Ti=(ωni2s2+2δiωnis+ωni2)ei
with ei∈ℝ, ei>0 corresponding to poles and ei>0 corresponding to zeros. Different combinations of these terms can be used depending on the problem to be solved. For a typical tracking problem with a minimum phase plant, the structure
(3.8)L0=KT1T2T3
with K∈ℝ+, two FCT poles and one FCT zero, achieves very good results.


## 3.1.1. TID

TID controller [26] is a modified version of PID controller, where the proportional term is replaced by a tilted (fractional) term, with transfer function s-eT, which permits a better approach to theoretical optimum. The resulting controller, including a low pass filter, is given by
(3.1)CTID(s)=k(TseT+1s+qDsq+s)1(1+s/wh)nh.


## 3.1.2.PIλDμ

PIλDμ controller [27, 28] is a PID generalization in which both integrator and derivative terms have real order, λ and μ, respectively. In this work, the multiplicative version of PIλDμ given in [29] will be used, with transfer function
(3.2)CPIλDμ(s)=kcxμ(λ1s+1s)λ(λ2s+1xλ2s+1)μ.
CRONE 2 (Section 3.1.3) is a particular case of (3.2), where some parameters are linked, and so flexibility is reduced compared to (3.2).


## 3.1.3. CRONE-Based Controllers

The CRONE approach [30, 31] defines three generations of fractional controllers based on the use of frequency-band noninteger differentiators. First and second generations (CRONE 1 and CRONE 2) use real non-integer differentiation, whereas third generation (CRONE 3) use complex non-integer differentiation. Three CRONE structures-based compensators are studied. The first one uses the transfer function structure common to CRONE 1 and 2. The basic component of this structure is an order n differentiator in the form of the implementable band-defined transfer function. An order nI band-limited integrator and an order nF lowpass filter are added to manage accuracy and robustness and control effort problems, being the open loop structure finally defined as
(3.3)LCR2(s)=k(ωls+1)nI(1+s/ωh1+s/ωl)n1(s/ωh+1)nF.
The application of the CRONE 2 structure to the QFT problem is quite straightforward.The second compensator uses CRONE 3 structure, consisting of the substitution of the (real) ordern integrodifferentiator in CRONE 2 for the real part Dr(s) of a (complex) order n=a+ib integrodifferentiator in the form of the implementable band-defined transfer function
(3.4)D(s)=Dr(s)+iDi(s)=(C01+s/ωh1+s/ωl)a+ib
with C0=ωh/ωu=ωu/ωl and ωu=ωlωh. The open loop structure L(s) in CRONE 3 is finally defined  as
(3.5)LCR3(s)=k(ωls+1)nI(C01+s/ωh1+s/ωl)a×cos[bLog(C01+s/ωh1+s/ωl)]1(s/ωh+1)nF.The third compensator is decoupled CRONE 3 ([32]), a modified version of CRONE 3 structure (3.5), where some parameters are decoupled in order to obtain higher flexibility. Frequencies ωh and ωl are decoupled by defining new frequencies ωh′, ωl′, and ωh4. C0 is decoupled by redefining C0 in cos() as C0′=cC0, where c is a new free parameter, C0=ωh′/ωu=ωu/ωl′, and ωu=ωl′ωh′. The new structure to be shaped  is
(3.6)LdecCR3(s)=k(ωls+1)nI(C01+s/ωh1+s/ωl)a×cos[bLog(cC01+s/ωh′1+s/ωl′)]1(s/ωh4+1)nF.
For both CRONE 3 and decoupled CRONE 3, a detailed study of the parameters involved, its interrelations and their allowed ranks in order to obtain desired behavior, is developed in [32]. In [33] a CRONE-based structure is proposed for the design of open loops based on Bode optimum-like specifications.


## 3.1.4. FCT Terms

Fractional order Complex Terms (FCTs) [19, 34] are terms
(3.7)Ti=(ωni2s2+2δiωnis+ωni2)ei
with ei∈ℝ, ei>0 corresponding to poles and ei>0 corresponding to zeros. Different combinations of these terms can be used depending on the problem to be solved. For a typical tracking problem with a minimum phase plant, the structure
(3.8)L0=KT1T2T3
with K∈ℝ+, two FCT poles and one FCT zero, achieves very good results.


## 3.2. Heuristics

Once the problem of a large number of parameters to be optimized has been solved by the use of fractional structures, the main problem is the fact that, due to the typically convex nature of the constraints in the solutions space, the evolutionary search can easily get stacked in local minima. The main goal for the heuristics design has been to help the evolutionary search to quickly avoid these local minima when they are detected.The objective function used as the criterion for the natural selection of individuals during the evolutionary optimization process,O(X), is a multiobjective real function which returns a to-be-minimized real value related with nonviolation of boundaries (constraint), stability (constraint), and minimization of Khf (to-be-optimized value). Violation of boundaries and lack of stability are, in fact, constraints, but due to how evolutionary algorithms work, they have to be translated into to-be-optimized values, as Khf originally is. A first approach could be to assign a large value to any solution which violates any boundary or is unstable. But this does not differentiate between solutions which violate a certain boundary completely, or only a little bit, which produces a blind O(X) in the sense that is not able to look for the way to become out of the violation of a certain given constraint. In order to avoid the evolutionary optimization getting stack in such situations, it is necessary to guide it by establishing, by adequately shaping O(X), that the violation of a certain constraint by a certain solution Xa is not as important as the violation produced by another solution Xb, so O(Xa)<O(Xb). This way, the evolutionary algorithm will prefer the survival of Xa instead of Xb, so that better solutions survive, even if they do not satisfy constraints. This scheme, repeated generation after generation, leads to solutions that do respect constraints. These are some of the components of O(X) intended to guide the evolutionary optimization towards solutions satisfying constraints.(i)
For open boundaries, letBω be an open boundary, and consider that it is defined as Bω(p):ℝ→ℝ, a function that assigns to every phase p in NC the magnitude of Bω at that phase. For certain solutions Xa and Xb, it should happen |LXi(jω)|≥Bω(∠LXi(jω)), i={a,b}. Assume that this constrain is not satisfied; that is, Bω is violated, by both Xa and Xb. Xa is considered better than Xb when LXa(jω) is closer to Bω than LXb is, in terms of magnitude. More formally, the penalty for this violation included in O(X) is a monotonically increasing function of Bω(∠LX(jω))-|LX(jω)|. This choice, made generation after generation, contributes to get solutions X such that LX(jω) is over Bω, thus satisfying Bω constraint.(ii)
For closed boundaries, letBω be a closed boundary, and consider that it is defined by Biω(p):ℝ→ℝ, i=u,l, two functions that assign to every phase p in NC the upper and lower magnitudes of Bω at that phase, respectively (bivalued boundary). For certain solutions Xa and Xb, it should happen |LXi(jω)|≥Buω(∠LXi(jω)) and |LXi(jω)|≤Blω(∠LXi(jω)), i={a,b}. Assume that this constrain is not satisfied; that is, Bω is violated, by both Xa and Xb. Xa is considered better than Xb when LXa(jω) is closer to Buω or Blω than LXb is, in terms of NC Euclidean distance. This choice, made generation after generation, contributes to get solutions X such that LX(jω) is out of Bω, thus satisfying Bω constraint.Another important consideration, in order to avoid the evolutionary search getting stacked, is anO(X) which gives more importance to constraints satisfaction than to khf minimization. This way, the evolutionary search first searches for valid solutions, and then, once they have been achieved, it optimizes among them.


## 3.3. Algorithm

This section describes the optimization algorithm which has been implemented for controllers synthesis. It consists on the use of commercial evolutionary algorithm software (theGenetic and evolutionary algorithm toolbox for use with MATLAB (GEATbx, [35])) together with an ad hoc objective function, programmed by the authors, implementing the heuristics described in Section 3.2.The evolutionary algorithm is, in particular, a multiple subpopulations evolutionary search algorithm. In this kind of evolutionary search, each subpopulation evolves in an isolated way for a few generations (as it happens in a single population evolutionary algorithm). After that, one or more individuals are exchanged between subpopulations. The way this process models species evolution is more similar to nature, compared to single population evolutionary algorithms, which helps to avoid local minima. The basic structure of this kind of algorithm is the following (adapted from [35]):PROCEDURE Multipopulation_Evolutionary_Search_Algorithm (search_space)BEGIN PROCEDUREINITIALIZATION, consisting on:creation of initial population // set of randomly chosen// points from search space,// grouped in speciesevaluation_of_individuals =OBJECTIVE_FUNCTION(initial_population);WHILE NOT(any termination criteria is met)BEGIN WHILENEW POPULATION GENERATION: // consisting on...fitness assignment selection;recombination;mutation;evaluation_of_offspring = OBJECTIVE_FUNCTION(offspring);reinsertion;migration;competition;END WHILE;RESULT = best individual in last population;END PROCEDURE;This is the code which implements the objective functionO(X) described in Section 3.2:PROCEDURE Objective_Function (set_of_individuals)BEGIN PROCEDUREFOR EACH individual IN set_of_individuals, indexed as iSBV = Compute_Stability_Boundary_Violation(individual);SRV = Compute_Stability_Ray_Violation(individual);PBV = Compute_Performance_Boundary_Violation(individual);NBV = Compute_Noise_Boundary_Violation(individual);PSP = Compute_Phases_Separation_Penalization(individual);MPP = Compute_Maximum_Phase_Penalization(individual);HFG = Compute_High_Frequency_Gain(individual);Constraints_Violation =BIG_VALUE * (weights[1] * SBV+ weights[2] * SRV+ weights[3] * PBV+ weights[4] * NBV);Technical_Penalizations =BIG_VALUE * (weights[5] * PSP+ weights[6] * MPP);Optimization_Criterion =weights[7] * HFG;RESULT[i] = Constraints_Violation+ Technical_Penalizations+ Optimization_Criterion;END FOR;END PROCEDURE;ValuesSBV, SRV, PBV and NBV are related with a direct violation of any of the problem constraints. This is the reason they are amplified by  BIG_VALUE  so that any solution satisfying constraints is preferred (has a lower objective value) compared to any other which does not satisfy any of the constraints.ValuesPSP and MPP are related with technical issues, which are necessary to consider in order to avoid nonreal solutions. Nonreal solutions are solutions that, due to intrinsic limitations of discrete computation, such as the need to consider a finite number of points to represent boundaries and loops, would satisfy the algorithm discrete translation of continuous constraints, but not the original continuous constraints. These values are also amplified by BIG_VALUE.ValueHFG is directly related with the optimization criterion, high-frequency gain minimization. It is not amplified by BIG_VALUE, and so it is only considered, in practice, when no constraint is violated.Weights are conceived to give more importance to a certain penalization compared to others, but have not been used for the moment; that is, all of them are equal.As explained in Section3.2, SBV, SRV, PBV, and NBV values correspond to a gradient, related with how much a certain boundary is violated, so that the objective function distinguishes different degrees of violation of a certain constraint, so that it helps the evolutionary algorithm to approach solutions which do not violate that constraint. For instance, consider the following MATLAB code, corresponding to the implementation of the function Compute Performance Boundary Violation (CPBV).PROCEDURE Compute_Perform_Boundary_Violation(loop_points,boundaries)loop_phases = loop_points.phases;loop_magnitudes = loop_points.magnitudes;bnd_phases = boundaries.phases;bnd_magnitudes = boundaries.magnitudes;nphases = length(bnd_phases);bnd_phases_indexes =mod(round(loop_phases*(nphases/360)+nphases),nphases);bnd_phases_indexes = nphases - bnd_phases_indexes;minimum_allowed_magnitudes = bnd_magnitudes(bnd_phases_indexes);magnitude_differences =loop_magnitudes - minimum_allowed_magnitudes;negative_differences =magnitude_differences .* (magnitude_differences < 0);squared_differences = negative_differences.̂2;RESULT = sum(squared_differences);END PROCEDURE;For each design frequencyω, this procedure checks whether the loop is below the performance boundary at that frequency, Bω, and in that case quantifies how much below it is, which is called magnitude difference in the code. These magnitude differences are squared, so that bigger magnitude differences are considered much worse than small differences. After that, these values are summed up to obtain the final boundary violation indicator, so that this indicator doubles when the loop violates boundaries at two frequencies, compared to a single violation, is triple when there are three violating frequencies, and so on. For an example of CPBV function behaviour, consider Figures 4, 5, and 6, representing loops Lex1, Lex2, and Lex3, respectively, corresponding to the first loop design in Section 4 (based on a PID controller). In Figure 4, performance boundaries are violated by Lex1 at design frequencies w1=0.1rad/s and w2=0.5 rad/s, by about 5 dB. The function yields a boundary violation indicator CPBV(Lex1)=4.1e3. Using BIG_VALUE  = 1e6, this produces an objective function value around O(Lex1)=4.1e9, which is a big penalization. The evolutionary algorithm will prefer individuals (loops) with lower objective function values. For instance, in Figure 5 performance boundaries are violated only at design frequency w1=0.1 rad/s, by 1 or 2 dB. CPBV function yields a boundary violation index CPBV(Lex2)=1.5e3, which corresponds to an objective function value O(Lex2)=1.5e9, still a big penalization, but lower than O(Lex1). This way, the evolutionary algorithm tends to move the loop, at each frequency ω∈Ω, out of F(jω) (forbidden area) and towards A(jω) (allowed area), producing a design like Lex3 in Figure 6, where no boundary is violated, and so CPBV(Lex3)=0 and O(Lex3) is contributed only by the optimization criterion, with no penalization by constraints violation, yielding O(Lex3)=152.Figure 4
Lex1 loop, with CPBV(Lex1)=4.1e3 and O(Lex1)=4.1e9.Figure 5
Lex2 loop, with CPBV(Lex2)=1.5e3 and O(Lex2)=1.5e9.Figure 6
Lex3 loop, with CPBV(Lex3)=0 and O(Lex3)=152.Figure7 shows the evolution, for an optimization of the loop structure shown in Figures 4, 5 and 6, of O(Li), where i is an index for each  generation in the evolutionary algorithm execution, and Li is the best individual (loop) in generation i, best in the sense of minimizing the objective function.   During the first three generations there are some boundary violations, which produce big objective function values. These values decrease as generation index i increases, due to the effect of the implemented heuristics. From fourth generation there is no boundary violation, so BIG_VALUE  is not applied any more, and so from that point the task of the evolutionary algorithm is reduced to get better and better values for the optimization criterion. This process can be better visualized in Figure 8, where scale has been adapted for this purpose. This figure permits checking how the objective function converges to a certain value, in this case around 152, from a certain generation, in this case around generation 130.Figure 7
Evolution ofO(Li) with generation index i, big scale, showing the effect of boundaries being violated at first generations.Figure 8
Evolution ofO(Li) with generation index i, reduced scale, showing the minimization of optimization criterion, once constraints are satisfied.Figure 9
PID-based loop design.


## 4. Design Example

To illustrate the behavior of the proposed optimization method, the QFT Toolbox for MATLAB [8] Benchmark Example number 2 is used. It has also been used, for instance, in [16, 17]. In [16], two rational loops are designed, a second-order one, with Khf=136.6 dB, and a third-order one, with Khf=130 dB, with npe=3 in both cases. A common npe=3 will be used along this section, so that loops can be compared in terms of Khf. In those structures which cannot get npe=3 by themselves,  a term
(4.1)H(s)=1(1+s/ωh)nh
is added to fix npe=3.For comparison purposes, a classical PID (TID witheT=0) based design (Figure 9) is also considered, with Khf=152 dB and parameters in (3.1): k=9.08*106, T=9.9*105, q=184.9, D=1.97, wh=2091.3, nh=2.The result obtained with TID controller is shown in Figure10, with Khf=140 dB, improving PID. Parameters in (3.1) are k=3.86*10-5, T=4.4*105, eT=0.36, q=5685.4, D=0.77, wh=7625.2, nh=2.Figure 10
TID-based loop design.PIλDμ improves TID result by another 12 dB, with Khf=128 dB, Figure 11. Parameters in (3.2), with term (4.1), are kc=2.79, x=1010, μ=0, λ=0.38, λ1=0.0083, λ2=1010, ωh=951.87, nh=2.Figure 11
PIλDμ-based loop design.CRONE 2 structure yields the loop shown in Figure12, with corresponding Khf=129.5 dB, and parameters in (3.3), are k=2.98*106, ωl=2588, nI=1.28, ωh=1.65*107, n=3.16, nF=3. This result is slightly worse than PIλDμ's, which could be expected, since CRONE 2 is the same structure as PIλDμ, but with some parameters linked, which implies less flexibility.Figure 12
CRONE-2-based loop design.In Figure13 it is shown a design using CRONE 3 structure, with Khf=126.8 dB, and parameters in (3.5), are k=0.0084, ωl=153, nI=1.4, C0=2.3, ωh=825, a=0.85, b=0.34, nF=3. The result is slightly better than using CRONE 2.Figure 13
CRONE-3-based loop design.In Figure14 it is shown a design using decoupled CRONE 3 structure, with Khf=105.3 dB, and parameters in (3.6), are k=1.46, ωl=7, nI=1.07, C0=11.2, ωh=256.9, ωl=7.1, a=1.5, b=0.45, c=0.7, ωh′=2*104, ωl′=166, ωh4=446, nF=3. This result significantly improves the original CRONE 3 design (in more than 20 dB).Figure 14
Decoupled CRONE-3-based loop design.Finally, in Figure15 it is shown the result obtained with the FCT structure (3.8), with Khf=94.2 dB, which improves decoupled CRONE 3 result by more than 10 dB. Parameters for (3.8), are K=8.3, ωn1=2.02, e1=0.16, δ1=2.32, ωn2=743.6, e2=1.4, δ2=0.28, ωn3=7145, e3=-0.57, δ3=2.55.Figure 15
FCT-based loop design.Figure 16
TN(s) comparison. From left to right, in terms of ωcg, FCT, decoupled CRONE 3, PIλDμ, CRONE 3, CRONE 2, TID and PID.Figure15 shows a comparison of the noise amplification at the plant input, TN(s)=-C/(1+L), achieved by each design. Table 1 summarizes the results obtained in terms of Khf. As it can be easily checked, there is a direct correlation between Khf and TN(s). Note how the most flexible the used structure is (and so the best Khf  it achieves) and the better Bode step-like shape its associated loop achieves.Table 1
Khf comparison.
CONTROLLERKhf (dB)PID152TID140PIλDμ128CRONE 2129.5CRONE 3126.8decoupled CRONE 3105.3FCT94.2


## 5. Conclusions

An automatic QFT controller design procedure, based on evolutionary algorithms optimization on the parameters of a fixed structure, has been proposed. The key idea behind this proposal is the introduction of a structure with few parameters (a must in order to get good results from evolutionary optimization) but, at the same time, flexible enough, thanks to its fractional nature, to get results which are close to the optimum. Fractional structures have been proposed as ideal candidates. Additional heuristics, focused on guiding the evolutionary search to prevent it from getting stacked in local minima, have been proposed. These structures and heuristics have achieved very good results in terms of QFT classical optimization criterion.



---

*Source: 671869-2009-09-27.xml*