Implementation of a New Quasi-Optimal Controller Tuning Algorithm for Time-Delay Systems

Systems and models with dead time or aftereffect, also called hereditary, anisochronic or time-delay systems (TDS), belonging to the class of infinite dimensional systems have been largely studied during last decades due to their interesting and important theoretical and practical features. A wide spectrum of systems in natural sciences, economics, pure informatics etc., both real-life and theoretical, is affected by delays which can have various forms; to name just a few the reader is referred e.g. to (Gorecki et al., 1989; Marshall et al., 1992; Kolmanovskii & Myshkis, 1999; Richard, 2003; Michiels & Niculescu, 2008; Pekař et al., 2009) and references herein. Linear time-invariant dynamic systems with distributed or lumped delays (LTI-TDS) in a single-input single-output (SISO) case can be represented by a set of functional differential equations (Hale & Verduyn Lunel, 1993) or by the Laplace transfer function as a ratio of so-called quasipolynomials (El’sgol’ts & Norkin, 1973) in one complex variable s, rather than polynomials which are usual in system and control theory. Quasipolynomials are formed as linear combinations of products of s-powers and exponential terms. Hence, the Laplace transform of LTI-TDS is no longer rational and socalled meromorphic functions have to be introduced. A significant feature of LTI-TDS is (in contrast to undelayed systems ) its infinite spectrum and transfer function poles decide except some cases of distributed delays, see e.g. (Loiseau, 2000) about the asymptotic stability as in the case of polynomials. It is a well-known fact that delay can significantly deteriorate the quality of feedback control performance, namely stability and periodicity. Therefore, design a suitable control law for such systems is a challenging task solved by various techniques and approaches; a plentiful enumeration of them can be found e.g. in (Richard, 2003). Every controller design naturally requires and presumes a controlled plant model in an appropriate form. A huge set of approaches uses the Laplace transfer function; however, it is inconvenient to utilize a ratio of quasipolynomials especially while natural requirements of internal (impulse-free modes) and asymptotic stability of the feedback loop and the feasibility and causality of the controller are to be fulfilled. The meromorphic description can be extended to the fractional description, to satisfy requirements above, so that quasipolynomials are factorized into proper and stable meromorphic functions. The ring of stable and proper quasipolynomial (RQ)


Introduction
Systems and models with dead time or aftereffect, also called hereditary, anisochronic or time-delay systems (TDS), belonging to the class of infinite dimensional systems have been largely studied during last decades due to their interesting and important theoretical and practical features. A wide spectrum of systems in natural sciences, economics, pure informatics etc., both real-life and theoretical, is affected by delays which can have various forms; to name just a few the reader is referred e.g. to (Górecki et al., 1989;Marshall et al., 1992;Kolmanovskii & Myshkis, 1999;Richard, 2003;Michiels & Niculescu, 2008;) and references herein. Linear time-invariant dynamic systems with distributed or lumped delays (LTI-TDS) in a single-input single-output (SISO) case can be represented by a set of functional differential equations (Hale & Verduyn Lunel, 1993) or by the Laplace transfer function as a ratio of so-called quasipolynomials (El'sgol'ts & Norkin, 1973) in one complex variable s, rather than polynomials which are usual in system and control theory. Quasipolynomials are formed as linear combinations of products of s-powers and exponential terms. Hence, the Laplace transform of LTI-TDS is no longer rational and socalled meromorphic functions have to be introduced. A significant feature of LTI-TDS is (in contrast to undelayed systems ) its infinite spectrum and transfer function poles decideexcept some cases of distributed delays, see e.g. (Loiseau, 2000) -about the asymptotic stability as in the case of polynomials. It is a well-known fact that delay can significantly deteriorate the quality of feedback control performance, namely stability and periodicity. Therefore, design a suitable control law for such systems is a challenging task solved by various techniques and approaches; a plentiful enumeration of them can be found e.g. in (Richard, 2003). Every controller design naturally requires and presumes a controlled plant model in an appropriate form. A huge set of approaches uses the Laplace transfer function; however, it is inconvenient to utilize a ratio of quasipolynomials especially while natural requirements of internal (impulse-free modes) and asymptotic stability of the feedback loop and the feasibility and causality of the controller are to be fulfilled. The meromorphic description can be extended to the fractional description, to satisfy requirements above, so that quasipolynomials are factorized into proper and stable meromorphic functions. The ring of stable and proper quasipolynomial (RQ) www.intechopen.com meromorphic functions (R MS ) is hence introduced (Zítek & Kučera, 2003;Pekař & Prokop, 2010). Although the ring can be used for a description of even neutral systems (Hale & Verduyn Lunel, 1993), only systems with so-called retarded structure are considered as the admissible class of systems in this contribution. In contrast to many other algebraic approaches, the ring enables to handle systems with non-commensurate delays, i.e. it is not necessary that all system delays can be expressed as integer multiples of the smallest one. Algebraic control philosophy in this ring then exploits the Bézout identity, to obtain stable and proper controllers, along with the Youla-Kučera parameterization for reference tracking and disturbance rejection. The closed-loop stability is given, as for delayless systems, by the solutions of the characteristic equation which contains a quasipolynomial instead of a polynomial. These infinite many solutions represent closed-loop system poles deciding about the control system stability. Since a controller can have a finite number of coefficients representing selectable parameters, these have to be set to distribute the infinite spectrum so that the closed-loop system is stable and that other control requirements are satisfied. The aim of this chapter is to describe, demonstrate and implement a new quasi-optimal pole placement algorithm for SISO LTI-TDS based on the quasi-continuous pole shiftingthe main idea of which was presented in (Michiels et al., 2002) -to the prescribed positions. The desired positions are obtained by overshoot analysis of the step response for a dominant pair of complex conjugate poles. A controller structure is initially calculated by algebraic controller design in R MS . Note that the maximum number of prescribed poles (including their multiplicities) equals the number of unknown parameters. If the prescribed roots locations can not be reached, the optimizing of an objective function involving the distance of shifting poles to the prescribed ones and the roots dominancy is utilized. The optimization is made via Self-Organizing Migration Algorithm (SOMA), see e.g. (Zelinka, 2004). Matlab m-file environment is utilized for the algorithm implementation and, consequently, results are tested in Simulink on an attractive example of unstable SISO LTI-TDS. The chapter is organized as follows. In Section 2 a brief general description of LTI-TDS is introduced together with the coprime factorization for the R MS ring representation. Basic ideas of algebraic controller design in R MS with a simple control feedback are presented in Section 3. The main and original part of the chapter -pole-placement shifting based tuning algorithm -is described in Section 4. Section 5 focuses SOMA and its utilization when solving the tuning problem. An illustrative benchmark example is presented in Section 6.

Description of LTI-TDS
The aim of this section is to present possible models of LTI-TDS; first, that in time domain using functional differential equations, second, the transfer function (matrix) via the Laplace transform. Then, the latter concept is extended so that an algebraic description in a special ring is introduced. Note that for the further purpose of this chapter the state-space functional description is useless.

State-space model
A LTI-TDS system with both input-output and internal (state) delays, which can have point (lumped) or distributed form, can be expressed by a set of functional differential equations where ∈ x  n is a vector of state variables, ∈ u  m stands for a vector of inputs, ∈ y  l represents a vector of outputs, are lumped (point) delays and convolution integrals express distributed delays (Hale & Verduyn Lunel, 1993;Richard, 2003;Vyhlídal, 2003 Convolution integrals in (1) can be numerically approximate by summations for digital implementation; however, this can destabilize even a stable system. Alternatively, one can integrate (1) and add a new state variable to obtain derivations on the right-hand side only.
In the contrary, the model can also be expressed in more consistent functional form using Riemann-Stieltjes integrals so that both lumped and distributed delays are under one convolution. For further details and other state-space TDS models the reader is referred to (Richard, 2003).

Input-output model
This contribution is concerned with retarded delayed systems in the input-output formulation governed by the Laplace transfer function matrix (considering zero initial conditions) as in (3). Hence, in the SISO case (we are concerning about here), the transfer function is no longer rational, as for conventional delayless systems, and a meromorphic function as a ratio of retarded quasipolynomials (RQ) is obtained instead.
A (retarded or neutral) quasipolynomial of degree n has the generic form However, the transfer function representation in the form of a ratio of two quasipolynomials is not suitable in order to satisfy controller feasibility, causality and closed-loop (Hurwitz) stability (Loiseau 2000;Zítek & Kučera, 2003). Rather more general approaches utilize a field of fractions where a transfer function is expressed as a ratio of two coprime elements of a suitable ring. A ring is a set closed for addition and multiplication, with a unit element for addition and multiplication and an inverse element for addition. This implies that division is not generally allowed.

Plant description in R MS ring
A powerful algebraic tool ensuring requirements above is a ring of stable and proper RQmeromorphic functions (R MS ). Since the original definition of R MS in (Zítek & Kučera, 2003) does not constitute a ring, some minor changes in the definition was made in . Namely, although the retarded structure of TDS is considered only, the minimal ring conditions require the use of neutral quasipolynomials at least in the numerator as well.
An element () MS Ts∈ R is represented by a proper ratio of two quasipolynomials

MS
As Bs∈ R are coprime, i.e. there does not exist a non-trivial (non-unit) common factor of both elements. Note that a system of neutral type can induce problem since there can exist a coprime pair () () , As Bs which is not, however, Bézout coprimewhich implies that the system can not be stabilized by any feedback controller admitting the Laplace transform, see details in (Loiseau et al., 2002).

Controller design in RMS
This section outlines controller design based on the algebraic approach in the R MS ring satisfying the inner Hurwitz (Bounded Input Bounded Output -BIBO) stability of the closed loop, controller feasibility, reference tracking and disturbance rejection. For algebraic controller design in R MS it is initially supposed that not only the plant is expressed by the transfer function over R MS but a controller and all system signals are over the ring. As a control system, the common negative feedback loop as in Fig. 1    The following basic transfer functions can be derived in the control system in general and () Qs , () Ps are from R MS and the fraction (13) is (Bézout) coprime (or relatively prime). The numerator of () M s ∈ R MS agrees to the characteristic quasipolynomial of the closed loop. Following subsections describes briefly how to provide the basic control requirements.

Reference tracking and disturbance rejection
The question is how to select () MS Ts∈ R in (16) where D ⋅ means that the output is influenced only by the disturbance, and symbol W ⋅ expresses that the signal is a response to the reference. Limit (17)  can be found e.g. in

Pole-placement shifting based controller tuning algorithm
In this crucial section, the idea of a new pole-placement shifting based controller tuning algorithm (PPSA) is presented. Although some steps of PPSA are taken over some existing pole-shifting algorithms, the idea of connection with pole placement and the SOMA optimization is original.

Overview of PSSA
We first give an overview of all steps of PPSA and, consequently, describe each in more details. The procedure starts with controller design in MS R introduced in the previous section. The next steps are as follows: is higher then den n , try to move the rest of dominant (rightmost) poles to the left. The same rule holds for shifted zeros, analogously. 5. If all prescribed poles and zeros are dominant, the procedure is finished. Otherwise, select a suitable cost function reflecting the distance of dominant poles (zeros) from prescribed positions and distances of spectral abscissas of both, prescribed and dominant poles (zeros). 6. Minimize the cost function, e.g. via SOMA. Now look at these steps of the algorithm at great length.

Characteristic quasipolynomial and characteristic entire function
Algebraic controller design in the MS R ring introduced in Section 3 results in a controller owning the transfer function () R Gs containing a finite number of unknown (free, selectable) parameters. The task of PPSA is to set these parameters so that the possibly infinite spectrum of the closed loop has dominant (rightmost) poles located in (or near by) the prescribed positions. If possibly, one can prescribe and place dominant zeros as well. Note www.intechopen.com that controller design in MS R using the feedback system as in Fig. 1 results in infinite spectrum of the feedback if the controlled plant is unstable.
If the (quasi)polynomial numerator and denominator of () Gs have no common roots in the open right-half plane, the closed-loop spectrum is given entirely by roots of the numerator () ms of () M s , the so called characteristic quasipolynomial. In the case of distributed delays,

()
Gs has some common roots with

{ }
Re 0 s ≥ in both, the numerator and denominator, and these roots do not affect the system dynamics since they cancel each other. In this case, the spectrum is given by zeros of the entire function () () / U ms m s , i.e. the characteristic entire function, where () U ms is a (quasi)polynomial the only roots of which are the common unstable roots. The (quasi)polynomial denominator of () WY Gs agrees with () ms . Its role is much more important than the role of the numerator of () WY Gs since the closed-loop zeros does not influence the stability. In the light of this fact, the setting of closed-loop poles has the priority. Therefore, one has to set den l free denominator parameters first. Free (selectable) parameters in the numerator of () WY Gs are to be set only if there exist those which are not contained in the denominator. The number of such "additional" parameters is num l .

Closed-loop model and step response overshoots
The task now is how to prescribe the closed-loop poles appropriately. We choose a simple finite-dimensional model of the reference-to-output transfer function and find its maximum overshoots and overshoot times for a suitable range of the model poles.

Direct pole placement
This subsection extends step 3 of PPSA from Subsection 4.1. The goal is to prescribe poles and zeros of the closed-loop "at once". The drawback here is that the prescribed poles (zeros) might not be dominant (i.e. the rightmost). The procedure was utilized to LTI-TDS e.g. in (Zítek & Hlava, 2001 , see (Ben Israel & Greville, 1966). Contrariwise, whenever nl > , it is not possible to place roots exactly and the pseudoinverse provides the minimization of squares of the left-hand sides of (30)-(34). The methodology described in this subsection is utilized on both, the numerator and denominator.

Continuous poles (zeros) shifting
Once the poles (zeros) are prescribed, it ought to be checked whether these roots are the rightmost. If yes, the PPSA algorithm stops; if not, one may try to shift poles so that the prescribed ones become dominant. There are two possibilities. First, the dominant roots move to the prescribed ones; second, roots nearest to the prescribed ones are shifted -while the rest of the spectrum (or zeros) is simultaneously pushed to the left. The following describes it in more details. We describe the procedure for the closed-loop denominator and its roots (poles); the numerator is served analogously for all its free parameters which are not included in the denominator. Recall that den l is the number of unknown (selectable) parameters, den n stands for the number of model (prescribed) poles (including their multiplicities), den n represents the number of real poles and conjugate pairs of prescribed poles and sp n is the number of currently shifted real poles and conjugate pairs. Generally, it holds that den sp den nn l ≤≤ The idea of continuous poles shifting described below was introduced in (Michiels et al., 2002). Similar procedure which, however, enables to shift less number of poles since sp den nl ≤ includes every single complex pole instead of a conjugate pair, was investigated in (Vyhlídal, 2003). Roughly speaking, the latter is based on solution of (30) - (34)  Δ are increments of poles and controller parameters, respectively. In case of a pmultiple pole, the following term is inserted in (36) and (37)  The continuous shifting starts with sp den nn = . Then, one can take the number of den n rightmost poles and move them to the prescribed ones. The rightmost closed-loop pole moves to the rightmost prescribed pole etc. Alternatively, the same number of dominant poles (or conjugated pairs) can be considered; however, the nearest poles can be shifted to the prescribed ones. If two or more prescribed poles own the same dominant pole, it is assigned to the rightmost prescribed pole and removed from the list of moved poles. The number

{ }
, sp den den nn l ∈ is incremented whenever the approaching starts to fail for any pole. If sp den nn > , the rest of dominant poles is pushed to the left. More precisely, shifting to the prescribed poles is described by the following formula If sp den nl = and all prescribed poles become the rightmost (dominant) ones, PPSA is finished. Otherwise, do the last step of PPSA introduced in the following subsection.

Minimization of a cost function via SOMA
This step is implemented whenever the exact pole assignment even via shifting fails. In the first part of this subsection we arrange the cost function to be minimized. Then, SOMA algorithm (Zelinka, 2004) belonging to the wide family of evolution algorithms is introduced and briefly described. Again, the procedure is given for the pole-optimization; the zerooptimization dealing with the closed-loop numerator is done analogously.

Cost function
The goal now is to rearrange feedback poles (zeros) so that they are "sufficiently close" to the prescribed ones and, concurrently, they are "as the most dominant as possible". This requirement can be satisfied by the minimizing of the following cost function Alternatively, one can include both, the zeros and poles, in (46), not separately.

www.intechopen.com
Poles can be found e.g. by the quasipolynomial mapping root finder (QPMR) implemented in Matlab, see (Vyhlídal & Zítek, 2003). Hence, the aim is to solve the problem We use SOMA algorithm based on genetic operations with a population of found solutions and moving of population specimens to each other. A brief description of the algorithm follows.

SOMA
SOMA is ranked among evolution algorithms, more precisely genetic algorithms, dealing with populations similarly as differential evolution does. The algorithm is based on vector operations over the space of feasible solutions (parameters) in which the population is defined. Population specimens cooperate when searching the best solution (the minimum of the cost function) and, simultaneously, each of them tries to be a leader. They move to each other and the searching is finished when all specimens are localized on a small area. In SOMA, every single generation, in which a new population is generated, is called a migration round. The notion of specific control and termination parameters, which have to be set before the algorithm starts, will be explained in every step of a migration round below. First, population { } 12 , ,..., PopSize P = vv v must be generated based on a prototypal specimen. For PPSA, this specimen is a vector of controller free parameters, v , of dimension den Dl = . The prototypal specimen equals the best solution from Subsection 4.5. One can choose an initial radius (Rad) of the population in which other specimens are generated. The size of population (PopSize), i.e. the number of specimens in the population, is chosen by the user. Each specimen is then evaluated by the cost function (46). The simplest strategy called "All to One" implemented here then selects the best specimenleader, i.e. that with the minimal value of the cost function ( ) arg min where L denotes the leader, i is i-th of specimen in the population and mr means the current migration round. Then all other specimen are moved towards the leader during the migration round. The moving is given by three control parameters: PathLength, Step, PRT. PathLength should be within the interval [1.1,5] and it expresses the length of the path when approaching the leader. PathLength = 1 means that the specimen stops its moving exactly at the position of the leader.
Step represents the sampling of the path and ought to be valued [0.11, ] PathLength . E.g. a pair PathLength = 1 and Step = 0.2 agrees with that the specimen makes 5 steps until it reaches the leader.
[ ] 0,1 PRT ∈ enables to calculate the perturbation vector PRTVector which indicates whether the active specimen moves to the leader directly or not. PRTVector is defined as www.intechopen.com where [ ] 0,1 i rnd ∈ is a randomly generated number for each dimension of a specimen. Although authors of SOMA suggest to calculate PRTVector only once in migration round for every specimen, we try to do this in every step of the moving to the leader. Hence, the path is given by Step , the active specimen goes to the leader directly without "zig-zag" moves. For every specimen of the population in a migration round, the cost function (i.e. value of the specimen) is calculated in every single step during the moving towards the leader. If the current position is better then the actual best, it becomes the best now. Hence, the new position of an active specimen for the next migration round is given by the best position of the specimen from all steps of moving towards the leader within the current migration round, i.e.

( )
where MinDiv is the selected minimal diversity. The final value opt v is equal to L v from the last migration round. We implemented the whole PPSA with SOMA in two Matlab m-files.

Illustrative example
In this closing session, we demonstrate the utilization of the PPSA and the methodology described above in Matlab on an attractive example. Consider an unstable system describing roller skater on a swaying bow (Zítek et al., 2008) given by the transfer function using the generalized Euclidean algorithm, see , where p 2 , p 1 , p 0 , q 3 , q 2 , q 1 , q 0 ∈  are free parameters. In order to provide reference tracking and load disturbance rejection, use parameterization (16) and the reference-to-output function reads = − . Hence l den = 7. Now, there are two possibilitieseither set zero exactly to obtain constrained controller parameter (then l den = 6) or to deal with the numerator and denominator of (61) together in (46) -we decided to utilize the former one. Generally, one can obtain e.g.
The concrete quasipolynomial which roots are being set, thus, reads and poles locations in the vicinity of the origin are displayed in Fig. 9. Fig. 9. Initial poles locations The process of continuous roots shifting is described by the evolution of controller parameters, the spectral abscissa (i.e. the real part of the rightmost pole ,1 d σ ) and the distance of the dominant pole from the prescribed one ,1 ,1 dp σσ − , as can be seen in Fig. 10  -Fig. 12, respectively. Note that p 0 is related to shifted parameters according to (64).  When shifting, it is suggested to continue in doing this even if the desired poles locations are reached since one can obtain a better poles distribution -i.e. non-dominant poles are placed more left in the complex space. Moreover, one can decrease the number of shifted poles during the algorithm whenever the real part of a shifted pole becomes "too different" from a group of currently moved poles.
and the poles location is pictured in Fig. 13. www.intechopen.com As can be seen, the desired prescribed pole is reached and it is also the dominant one. Thus, optimization can be omitted. However, try to perform SOMA to find the minimal cost function (16) with this setting: Rad = 2, PopSize = 10, D = 6, PathLength = 3, Step = 0.21, PRT = 0.6, Migration = 10, MinDiv = 10 -6 , . Yet, the minimum of the cost function remains in the best solution from continuous shifting, i.e. according to (67)

Conclusion
This chapter has introduced a novel controller design approach for SISO LTI-TDS based on algebraic approach followed by pole-placement-like controller tuning and an optimization procedure. The methodology has been implemented in Matlab-Simulink environment to verify the results. The initial controller structure design has been made over the ring of stable and proper meromorphic functions, R MS , which offers to satisfy properness of the controller, reference tracking and rejection of the load disturbance (of a nominal model). The obtained controller has owned some free (unset) parameters which must have been set properly.
In the crucial part of the work, we have chosen a simple finite-dimension model, calculated its step-response maximum overshoots and times to the overshoots. Then, using a static pole placement followed by continuous pole shifting dominant poles have been attempted to be placed to the desired prescribed positions. Finally, optimization of distances of dominant (the rightmost) poles from the prescribed ones has been utilized via SOMA algorithm. The whole methodology has been tested on an attractive example of a skater on a swaying bow described by an unstable LTI TDS model. The procedure is similar to the algorithm introduced in (Michiels et al., 2010); however, there are some substantial differences between them. Firstly, the presented approach is made in input-output space of meromorphic Laplace transfer functions, whereas the one in (Michiels et al., 2010) deals purely with state space. Second, in the cited literature, a number poles less then a number of free controller parameters is set exactly and the rest of the spectrum is pushed to the left as much as possible. If it is possible it is necessary to choose other prescribed poles. We initially place the poles exactly; however, they can leave their positions during the shifting. Anyway, our algorithm does not require reset of selection assigned poles. Moreover, we suggest unambiguously how the prescribed poles (and zeros) positions are to be chosen -based on model overshoots. Last but not least, in (Michiels et al., 2010), the gradient sampling algorithm (Burke et al., 2005) on the spectral abscissa was used while SOMA together with more complex cost function is considered in this chapter. The presented approach is limited to retarded SISO LTI-TDS without distributed delays only. Its extension to neutral systems requires some additional conditions on stability and existence of a stabilizing controller. Systems with distributed delays can be served in similar way as it is done here, yet with the characteristic meromorphic function instead of quasipolynomial. Multivariable systems would require deeper theoretic analysis of the controller structure design. The methodology is also time-comsupting and thus useless for online controller design (e.g. for selftuners).
In the future research, one can solve the problems specified above, choose other referenceto-output models and control system structures. There is a space to improve and modify the optimization algorithm.