Nonlinear Nonconvex Optimization by Evolutionary Algorithms Applied to Robust Control

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.


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.

Mathematical Problems in Engineering
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 Section 2 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.

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.
Consider an uncertain plant, represented as the set of transfer functions P {P s P 0 s Q s , Q ∈ Q}, with P 0 s being the nominal plant and Q being a set of transfer functions representing uncertainty.For a certain frequency ω, the template T P jω is defined as the Nichols chart NC region 22 given by T P ω ˙ ∠P jω , P jω dB ∈ NC, P ∈ P .

2.1
The design of the controller C s is accomplished in the Nichols chart, by shaping the nominal open loop transfer function, L 0 s P 0 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 boundaries B {B ω , ω ∈ Ω} is computed.Each B ω defines the allowed region A ω conversely, forbidden region F ω NC \ A ω in NC for L 0 jω , so that L 0 jω ∈ A ω conversely, L 0 jω / ∈ F ω for all ω ∈ Ω boundaries satisfaction implies specification satisfaction by L s P s C s for all P 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, F UHFB , is not forbidden only for a particular frequency ω, but for all ω ∈ R , so it imposes L 0 jω / ∈ F UHFB for all ω ∈ R .
The basic step in the design process, loop shaping, consists of the following nonlinear nonconvex constrained optimization problem: design of L 0 s which satisfies boundaries constraints and minimizes the high-frequency gain optimization criterion-23 , defined as where e is L 0 excess of poles over zeros.The comparison in terms of this optimization criterion of two nominal open loops L 0 1 s and L 0 2 s , with respective excess of poles over zeros e 1 and e 2 , only makes sense if e 1 e 2 .The comparison is performed in terms of the loops highfrequency gains, K 0 1 and K 0 2 , respectively.Note that this criterion is defined irrespective of the particular structure used for L 0 s .It has been shown 23 that for minimum phase L 0 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, L qto s , with optimal high-frequency gain K qto .
Since L qto s is not practical as a loop design, a more practical definition of optimum can be stated.Define X ˙ {x x 1 , . . ., x n }, x i ∈ R, i 1 • • • n, as the set of controller structure parameters to be instantiated.Consider L X s an open loop structure parameterized by X, that is, each x ∈ X defines L x s , a particular instance of L X s , with high-frequency gain K L x .Define K L O ˙ min{K L x , x ∈ X} ≤ K qto , the minimum high-frequency gain that can be achieved by any L X s instance.An L-optimal transfer function or L O is defined as an L X s instance L a such that K L a K L O .
For a given L X s open loop structure, K L O represents the best approach to K opt that can be achieved with that structure.The ability of a structure L X s to achieve L qto s shape is called close-to-optimality, defined as CO L ˙ K L O − K qto .Since computing L O is in general still a hard optimization problem, and thus K L O is not known, it is common to compute instead suboptimal transfer functions L sOi with high-frequency gain K L sOi as close to K L O as possible, with i ∈ N.This is the approach considered in this work, using evolutionary algorithms to compute an L sOi with K L sOi which approaches K L O as much as possible.A measure of how accurately a certain L sO1 approaches K L O can only be established in relative terms, that is, by comparing it to another L sO2 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 |L 0 jω | at high frequencies, which, due to phase/magnitude relations given by Bode integrals 24 , is equivalent to minimize R ∠L 0 jω dω, subject to B satisfaction.Figures 2 and 3 show the typical shape of a QFT optimal loop L qto jω , with polygonal line simplified versions of L qto jω and boundaries in B. For frequencies ω ∈ 0, ω a2 , ∠L qto jω minimization is limited by performance boundaries.For frequencies ω ∈ ω a2 , ∞ , ∠L qto jω minimization is only limited by the UHFB which makes L qto 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.
Note the shape of L qto 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 L qto jω has constant magnitude and a fast phase decrease according to ∠L qto jω minimization objective, once the UHFB does not constrain this objective , surrounded by frequencies where ∠L qto jω is constant απ and −eπ/2, resp.and |L qto jω | decreases.Bode step is very typical of QFT optimal loops, and when designing a suboptimal loop L sO , the designer has to try to make L sO exhibit a Bode step-like shape in order to succeed in accurately approaching L qto .Note that a Bode step in L sO jω can be interpreted as L sO jω tightly following the bottom part of the UHFB.

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 number n of free parameters x i ∈ R, i 1 • • • n, which constitute the solution space X ˙ {x x 1 , . . ., x n } 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 reduced n.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 Section 3.2.Section 3.3 describes the algorithm used for evolutionary optimization, detailing the objective function, which includes these heuristics.
In Section 4 a design example is solved by using these heuristics and all the fractional structures.The results obtained by each fractional structure are compared.

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 −e T , which permits a better approach to theoretical optimum.The resulting controller, including a low pass filter, is given by

P I λ 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 is a particular case of 3.2 , where some parameters are linked, and so flexibility is reduced compared to 3.2 .

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 n I band-limited integrator and an order n F lowpass filter are added to manage accuracy and robustness and control effort problems, being the open loop structure finally defined as 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 order n integrodifferentiator in CRONE 2 for the real part D r s of a complex order n a ib integrodifferentiator in the form of the implementable band-defined transfer function with C 0 ω h /ω u ω u /ω l and ω u √ ω l ω h .The open loop structure L s in CRONE 3 is finally defined as

3.5
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 .C 0 is decoupled by redefining C 0 in cos as C 0 cC 0 , where c is a new free parameter, C 0 ω h /ω u ω u /ω l , and ω u ω l ω h .The new structure to be shaped is

3.6
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.

FCT Terms
Fractional order Complex Terms FCTs 19, 34 are terms 7 with e i ∈ R, e i > 0 corresponding to poles and e i > 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 with K ∈ R , two FCT poles and one FCT zero, achieves very good results.

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 K hf 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 K hf 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 X a is not as important as the violation produced by another solution X b , so O X a < O X b .This way, the evolutionary algorithm will prefer the survival of X a instead of X b , 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, let B ω be an open boundary, and consider that it is defined as B ω p : R → R, a function that assigns to every phase p in NC the magnitude of B ω at that phase.For certain solutions X a and X b , it should happen Assume that this constrain is not satisfied; that is, B ω is violated, by both X a and X b .X a is considered better than X b when L X a jω is closer to B ω than L X b is, in terms of magnitude.More formally, the penalty for this violation included in O X is a monotonically increasing function of B ω ∠L X jω − |L X jω |.This choice, made generation after generation, contributes to get solutions X such that L X jω is over B ω , thus satisfying B ω constraint.
ii For closed boundaries, let B ω be a closed boundary, and consider that it is defined by Bi ω p : R → R, 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 X a and X b , it should happen Assume that this constrain is not satisfied; that is, B ω is violated, by both X a and X b .X a is considered better than X b when L X a jω is closer to Bu ω or Bl ω than L X b is, in terms of NC Euclidean distance.This choice, made generation after generation, contributes to get solutions X such that L X jω is out of B ω , thus satisfying B ω constraint.
Another important consideration, in order to avoid the evolutionary search getting stacked, is an O X which gives more importance to constraints satisfaction than to k hf minimization.This way, the evolutionary search first searches for valid solutions, and then, once they have been achieved, it optimizes among them.

Algorithm
This section describes the optimization algorithm which has been implemented for controllers synthesis.It consists on the use of commercial evolutionary algorithm software the Genetic 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 : 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 L ex1 , L ex2 , and L ex3 , respectively, corresponding to the first loop design in Section 4 based on a PID controller .In Figure 4, performance boundaries are violated by L ex1 at design frequencies w 1 0.1 rad/s and w 2 0.5 rad/s, by about 5 dB.The function yields a boundary violation indicator CPBV L ex1 4.1e3.Using BIG VALUE 1e6, this produces an objective function value around O L ex1 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 w 1 0.1 rad/s, by 1 or 2 dB.CPBV function yields a boundary violation index CPBV L ex2 1.5e3, which corresponds to an objective function value O L ex2 1.5e9, still a big penalization, but lower than O L ex1 .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 L ex3 in Figure 6, where no boundary is violated, and so CPBV L ex3 0 and O L ex3 is contributed only by the optimization criterion, with no penalization by constraints violation, yielding O L ex3 152. Figure 7 shows the evolution, for an optimization of the loop structure shown in Figures 4, 5 and 6, of O L i , where i is an index for each generation in the evolutionary  algorithm execution, and L i 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.

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 K hf 136.6 dB, and a third-order one, with K hf 130 dB, with n pe 3 in both cases.A common n pe 3 will be used along this section, so that loops can be compared in terms of K hf .In those structures which cannot get n pe 3 by themselves, a term  For comparison purposes, a classical PID TID with e T 0 based design Figure 9 is also considered, with K hf 152 dB and parameters in 3.1 : k 9.08 * 10 6 , T 9.9 * 10 5 , q 184.9, D 1.97, w h 2091.3,n h 2.
The result obtained with TID controller is shown in Figure 10  CRONE 2 is the same structure as PI λ D μ , but with some parameters linked, which implies less flexibility.
In Figure 13 it is shown a design using CRONE 3 structure, with K hf 126.8 dB, and parameters in 3.5 , are k 0.0084, ω l 153, n I 1.4,C 0 2.3, ω h 825, a 0.85, b 0.34, n F 3. The result is slightly better than using CRONE 2.
In Figure 14 it is shown a design using decoupled CRONE 3 structure, with K hf 105. 3        Figure 15 shows a comparison of the noise amplification at the plant input, T N s −C/ 1 L , achieved by each design.Table 1 summarizes the results obtained in terms of K hf .As it can be easily checked, there is a direct correlation between K hf and T N s .Note how the most flexible the used structure is and so the best K hf it achieves and the better Bode step-like shape its associated loop achieves.

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.

Figure 1 :
Figure 1: Two degrees of freedom control system configuration.

Figure 8 :Figure 9 :
Figure 8: Evolution of O L i with generation index i, reduced scale, showing the minimization of optimization criterion, once constraints are satisfied.