GENETIC ALGORITHMS IN SYSTEM IDENTIFICATION AND CONTROL By Kristinn Kristinsson B. Sc. Electrical Engineering University of Iceland, 1986 A T H E S I S S U B M I T T E D I N P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F M A S T E R O F A P P L I E D S C I E N C E in T H E F A C U L T Y O F G R A D U A T E S T U D I E S E L E C T R I C A L E N G I N E E R I N G We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A August 1989 © Kristinn Kristinsson, 1989 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Electrical Engineering The University of British Columbia 2075 Wesbrook Place Vancouver, Canada V6T. 1W5 Date: Abstract Current online identification techniques are recursive and involve local search tech-niques. In this thesis, we show how genetic algorithms, a parallel, global search tech-nique emulating natural genetic operators can be used to estimate the poles and zeros of a dynamical system. We also design an adaptive controller based on the estimates. The algorithms are shown to be useful for continuous time parameter identifications and to be able to identify directly physical parameters of a system. Simulations and an experiment show the technique to be satisfactory and to provide unbiased estimates in presence of colored noise. ii Table of Contents Abstract ii List of Tables vi List of Figures vii Acknowledgement xi 1 Introduction 1 1.1 General introduction 1 1.2 "Standard" identification methods 1 1.3 Motivation for this work 2 1.4 Outline of this thesis 3 2 Genetic Algorithms 4 2.1 History of Genetic Algorithms 4 2.2 The algorithm 4 2.3 Coding '. 7 2.4 Reproduction 7 2.4.1 Ranking 9 2.5 Crossover 10 2.6 Mutation 11 2.7 Mathematical Foundations 11 2.8 Example 13 iii 2.9 Summary 14 3 System Identification 19 3.1 Background 19 3.2 Discrete time identification 21 3.2.1 Parameter identification 25 3.2.2 Pole-zero identification 25 3.2.3 Results 27 3.3 Recursive Least Squares estimation 31 3.3.1 Results 32 3.4 Continuous time identification 38 3.4.1 Results 39 3.5 Friction compensation 41 3.5.1 Results 43 3.6 Summary 44 4 Controller design 46 4.1 Controller 46 4.2 Parameter based design 48 4.3 Pole-zero based design 50 4.3.1 Multiple poles and zeros 53 4.3.2 Complex poles and zeros 53 4.3.3 Implementation 54 4.4 Results 54 4.4.1 Minimum phase plant 56 4.4.2 Nonminimum phase plant 64 4.4.3 Unrnodeled dynamics 65 iv 4.4.4 Persistently exciting signal 73 4.4.5 Recursive Least Squares 76 4.5 Summary 79 5 Experiment 82 5.1 Water level in a tank 82 5.2 Simulation results 83 6 Conclusions 87 Bibliography 89 A Genetic algorithms procedures 92 v List of Tables 2.1 Fitness dependent reproduction 8 2.2 Reproduction 15 2.3 Crossover and mutation 15 3.4 Search space for a stable minimum phase system 28 3.5 Search space for continuous time parameters 39 4.6 Search space for nonminimum phase 65 4.7 Search space for unmodeled dynamics 68 vi List of Figures 2.1 Genetic Algorithm 6 2.2 Ranking 10 2.3 Function with 11 local maxima 14 2.4 G A generations 16 2.4 G A generations (continued) 17 3.5 Window size = 5 22 3.6 Window size = 10 22 3.7 Window size = 20 23 3.8 Window size = 30 23 3.9 Number of trials = 1 24 3.10 Number of trials = 2 24 3.11 Number of trials = 3 24 3.12 Number of trials = 30 24 3.13 Parameters 26 3.14 Poles and zeros in the complex plane 26 3.15 Reparameterized complex plane 26 3.16 P R B S input and output of a system without noise 29 3.17 Pole-Zero estimate of a system without noise 29 3.18 Estimation of gain and delay of a system without noise 30 3.19 Pole zero locations 30 3.20 P R B S input and output with noise 33 vii 3.21 Parameter identification using RLS 33 3.22 Parameters locations 34 3.23 G A , pole zero identification 34 3.24 G A , parameters identification calculated from the pole zero identification 35 3.25 G A , parameters locations calculated from the pole zero identification . . 35 3.26 G A , parameters identification 36 3.27 G A , parameters locations 36 3.28 Continuous-time-system input and output 39 3.29 Continuous-time parameter estimates 40 3.30 Actual output and the output using the final estimates 40 3.31 Friction model 42 3.32 Motor input (/) and output (w) 43 3.33 Friction parameters identification 44 4.34 Two-degree of freedom controller 48 4.35 G A adaptive controller 49 4.36 Parameter Controller 50 4.37 Factorized Controller 55 4.38 Ladder 55 4.39 Reference input and output of a minimum phase system without noise . 56 4.40 Pole-Zero estimates for a minimum phase system without noise 57 4.41 Estimates of gain and delay for a minimum phase system without noise 57 4.42 Pole zero locations 58 4.43 Reference input and output of a minimum phase system with noise using pole-zeros estimates 59 4.44 Pole-zero estimates for a minimum phase system with noise 60 v i n 4.45 Estimates of gain and delay for a minimum phase system with noise . . 60 4.46 Pole zero locations 61 4.47 Reference input and output of a minimum phase system with noise using parameter estimates 62 4.48 Parameter estimates for a minimum phase system with noise 62 4.49 Estimates of gain and delay for a minimum phase system with noise . . 63 4.50 Parameters locations 63 4.51 Reference input and output for a nonminimum phase system 66 4.52 Parameter estimate for a nonminimum phase system 66 4.53 Gain and delay estimate for a nonminimum phase system 67 4.54 Pole zero locations for a nonminimum phase system 67 4.55 Input-Output for 3 parameters estimate with unmodeled dynamics . . . 68 4.56 Parameter estimate for unmodeled dynamics 69 4.57 Gain and delay estimate for unmodeled dynamics 69 4.58 Input-Output with dead beat control 70 4.59 Parameter estimate for unmodeled dynamics with dead beat control . . 70 4.60 Gain and delay estimate for unmodeled dynamics with dead beat control 71 4.61 Input-Output with desired pole = 0.7 71 4.62 Parameter estimates for unmodeled dynamics with desired pole = 0.7 . 72 4.63 Gain and delay estimate for unmodeled dynamics with desired pole = 0.7 72 4.64 Reference input and output of a system with window size = 30 73 4.65 Parameters for a window size = 30 74 4.66 Gain and delay for a window size = 30 74 4.67 Reference input and output of a system with window size = 60 75 4.68 Parameters for a window size = 60 75 4.69 Gain and delay for a window size = 60 76 ix 4.70 Reference input and output of a system using RLS 77 4.71 Parameter estimates for RLS 77 4.72 Parameter locations for RLS 78 4.73 Reference input and output of a system using G A to compare to RLS 79 4.74 Parameter estimates for G A to compare with RLS 80 4.75 Parameter locations for G A 80 5.76 Tank Input-Output 83 5.77 Parameter estimate for a tank 84 5.78 Estimated gain for a tank 84 5.79 Pole zero locations for a tank 85 x Acknowledgement I would like to use this opportunity for thanking all those that have made the completion of this thesis possible. I would especially like to thank Prof. Guy A . Dumont, my thesis advisor, for introducing the Genetic Algorithms to me and for his advice and guidance in my research. I would also like to thank Ye Fu for providing the RLS routines and Rob Ross for his assistance in operation the Pulp & Paper Centre's fxVAX computer. At last special thanks go to Dr. K . Natarajan who read the final draft of this thesis. xi Chapter 1 Introduction 1.1 General introduction The area of system identification has been given a lot of attention over the years. Many methods have been used and many extended versions exist, but all of them are based upon eighteenth century mathematics which assumes smooth search space with ever present derivatives. In the last few years Artificial Intelligence and learning have been gaining lot of popularity and have been entering many fields, but little has been done to apply them in the field of system identification and control. 1.2 "Standard" identification methods On-line system identification methods used to date are based on recursive implemen-tation of off-line methods such as least-squares, maximum-likelihood or instrumental variable. All those methods are based on the same principle and a unified description exists [23]. Those recursive schemes are in essence local search techniques that search for zero gradient by going in a direction suggested by the local gradient. They go to the nearest point that gives zero gradient and stay there. Nothing will get the methods to search further as long as the gradient stays zero. It is therefore very difficult for those methods to find a global maximum and they often fail in the search for global maxi-mum if the search space is not differentiable or linear in the parameters. Because of 1 Chapter 1. Introduction 2 the linearity condition they have difficulty locating directly poles and zeros or physical parameters of a system. Another aspect is that these methods are all serial. They go from one point in the search space to another at every sampling instant, as a new input-output pair becomes available. They are not capable of iterating more than once on each data they receive, they need new data to direct the search. 1.3 Motivation for this work Genetic algorithms are a parallel, global search technique that emulates natural genetic operators. They search many points simultaneously and thus have the potential to converge more rapidly. In every generation new artificial chromosomes are created by taking parts of the fittest chromosomes of the previous generation and combining them to make a highly fit chromosome. They do not need to assume that the search space is differentiable or continuous because they go from one generation to another with transition rules that are probabilistic. This means the algorithms do not have to wait for new data, but can iterate a few times on each data they receive. They work with a population of binary coded strings so they can explore the search space in each generation and then direct the search to regions where there is a high probability of finding improved performance. Genetic algorithms have also been shown to excel in multimodal optimization [10], and thus have the potential to give unbiased estimates in presence of coloured noise. In this thesis, a Genetic Algorithm (GA) is implemented as an estimator for dis-crete time systems. The results obtained employing this new identification method are particularly favourable and they are considered to be well suited to the adaptive Chapter 1. Introduction 3 control problem. Although the use of G A has been gaining popularity, its use in adap-tive control has not been investigated. The algorithm is used on few discrete time systems, both minimum and nonminimum phase and with or without colored noise. It is used to identify either parameters or poles-and-zeros. The encouraging results are then compared with Recursive Least Squares. 1.4 Outline of this thesis This thesis is organised as follows: Chapter 2. Genetic Algorithms are described and some of the simple genetic opera-tors are explained. The algorithm is then used to find the maximum for a function with eleven local maxima. Chapter 3. A G A for system identification is implemented both in discrete and con-tinuous time and simulations results are shown. The algorithm is also compared to Recursive Least Squares algorithm. Chapter 4. The pole placement controller design is outlined and simulations results are shown using the G A to identify plants with either minimum phase or non-minimum phase characteristic and unmodeled dynamics. One simulation is then done using Recursive Least Squares for the identification for comparison. Chapter 5. A n experiment with a water tank is explained and identification results are shown. Chapter 6. Conclusions and suggestions for further work are given. Chapter 2 Genetic Algorithms 2.1 History of Genetic Algorithms The algorithms come out of work done by John H . Holland and his students at the University of Michigan. The underlying principles of genetic algorithms were first published by Holland in 1962 [18]. The mathematical framework was developed in the late 1960's and in 1975, Holland's pioneering book, Adaptation in Natural and Artificial Systems was published [19]. The same year it was shown by one of his student that Genetic Algorithms (GAs) are very useful in function optimization even on "difficult" domains that are multimodal, noisy and high-dimensional [10]. In 1983, Goldberg used G A to minimize power consumption in gas-pipelines and then combined a learning classifier system with a G A to detect leakage in the system [13]. In the last four years a lot of research has been devoted to G A , two conferences have been held [16,17] and two books have been written on the subject [15,9]. Genetic Algorithms have proven to be useful in many different applications [15], like function optimization, computer network design, travelling salesman problem, pattern recognition and many more. 2.2 The algorithm Genetic algorithms differ from other search techniques by the use of concepts taken from natural genetics and evolution theory. They are different in four ways: 1. GAs work with coding of the parameters, not the parameters themselves. 4 Chapter 2. Genetic Algorithms 5 2. GAs search from a population of points, not a single point. 3. GAs only need fitness values. There is no requirement for derivatives or other auxiliary knowledge. 4. GAs use probabilistic transition rules, not deterministic rules. The parameters to be found by the G A , need to be coded as a finite length string over a finite search space. As an example, consider a stable real pole, (magnitude less than one) the search space would be on the interval [0,1] or if we want a resolution of 1/1000 the search space would be in the integer interval [0,1000] and with a binary coding this would be coded as a 10 bit string. The algorithm works with a population of strings, searching many peaks in parallel, hence reducing the possibility of ending at a local minimum and missing the global minimum. The only available feedback from the system is the value of the performance measure (fitness). Although transition rules are probabilistic, the algorithm is not simply a random search. It is a randomized search that is guided by the fitness values of each string. The algorithm uses information, already in the population, about things that have worked well in the past. By use of operators taken from population genetics the algorithm efficiently explores part of the search space where the probability of finding improved performance is high. A genetic algorithm in its simplest form consists of 3 steps: (see Figure 2.1) 1. Reproduction 2. Crossover 3. Mutation In the next three sections these will be described in details. Chapter 2. Genetic Algorithms 6 Randomly generated initial population r " O n e - ~| generation Reproduction Crossover Mutation L I J Figure 2.1: Genetic Algorithm Chapter 2. Genetic Algorithms 7 2.3 C o d i n g It has been shown that binary coding is in a certain sense the optimal coding [19]. Suppose we have the binary strings 1010010111 and 1110100110, by comparing them we can see some similarity, 1 * 10* *011*, where * is a don't care. We call 1 * 10* *011*, a schema (plural, schemata). A string can be an instance of 2 1 0 = 1024 schemata which can be found by replacing any of the bits in the string by a don't care. The number of possible schemata on the alphabet {*,0,1} for the binary coding is 3 1 0 = 59094 so by carefully selecting the string all schemata could be represented by 58 strings. If the coding is decimal we would need 3 decimal number to represent the same search space as a 10 bit binary string would. A three bit decimal string would be an instance of 2 3 = 8 schemata but the number of possible schemata on the alphabet {*,0,1, . . . , 9} would be l l 3 = 1331 so all schemata would have to be represented by 166 strings. So binary coding would need 3 times fewer strings to explore the search space. Therefore binary coding is chosen with each parameter corresponding to a fixed length binary substring of j bits [0 , . . . , 2d — 1]. The value, (x), of the binary substring is mapped to an interval of the real numbers [/6,u6] to give y = -^—(Ub - lb) + lb (2.1) With n parameters, the final string consists of n concatenated substrings. I ai | • • • | a n | (2.2) | 10 - -01 | ••• | 01---11 I 2.4 R e p r o d u c t i o n In the reproduction part of the algorithm it is decided which strings are going to survive and which ones are going to disappear, based on what in biological terms, is known as Chapter 2. Genetic Algorithms 8 number of F(i) offsprings 100 0.50 1 10 0.05 0 200 1.00 1 300 1.50 2 210 1.05 1 290 1.45 1 310 1.55 1 280 1.40 2 120 0.60 1 180 0.90 0 Table 2.1: Fitness dependent reproduction. the survival of the fittest principle. It is done by assigning a positive number, fitness F(i), to each individual in the population. It must be positive because high fitness individuals should receive more offsprings than low fitness individuals. Based on the normalized fitness Fn(i), the number of offsprings for each individual is calculated. The fitness function tells us how well the system, to optimize or control, is behaving under a certain string. The fitness function can be any nonlinear, nondifferentiable, discontinuous function, because the algorithm only needs a fitness value assigned to each string. The number of offsprings is chosen according to the string normalized fitness. The fitness is normalized with the average value of the fitness, £ . ( 0 = • F}'] (2-3) 1=1 so the strings with above average fitness will have more than 1 offspring and those with below average fitness will have less than 1 offspring on the average (see Table 2.1). The strings are selected according to the expected value of the normalized fitness or what has become known as Stochastic Remainder Selection without Replacement, Chapter 2. Genetic Algorithms 9 [15]. That means the strings will receive number of offsprings equal to the integer value of their normalized fitnesses and then the population is filled up by choosing another offspring for each of the strings with probability equal to their fractional part until the total number of offsprings are equal to the population size TV. The algorithm keeps track of the best string in the population and if it is not in the new population (because some other G A operators destroyed it) it randomly replaces another string in the new population. 2.4.1 Ranking It is important to regulate the number of offsprings an individual can get, to maintain a diversity in the population. Especially for the first few generations, when a few "super" individuals can potentially take over a large part of the population, thereby reducing the diversity of the population. The presence of super individuals can be sensed by monitoring the number of individuals that are going to receive 0 offsprings. It is somewhat a better way than limiting the number of offsprings an individual can get, because it could be desirable to give a good individual many offsprings as long as the diversity is maintained. To control the reproduction, ranking can be introduced [3]. Whenever a certain ratio of the normalized fitness is going to receive 0 offsprings, the strings are sorted according to their fitness values. Then, instead of calculating the normalized fitness as in Equation 2.3, the normalized fitness is given to each string according to „ ... 2 ( m a x - l ) w . , , x iV + l , F ^ = / y _ i + 1 - (ma* ~ VjfZi <2-4) where max as shown in Figure 2.2 is a user defined value, 1 < max < 2, and N is the population size. The range of the normalized fitness will then be [2-max,max]. This means that no matter how big the fitness is for the best string its normalized fitness Chapter 2. Genetic Algorithms 10 max 1 2-max - -1 N rank(i) Figure 2.2: Ranking will never be more than max when ranking is in effect. The lowest ranking string will similarly always be guaranteed 2 — max as its normalized fitness. 2.5 Crossover Reproduction directs the search towards the best but does not create any new indi-viduals. The offsprings are identical to their parent. In nature, the offsprings are not exact copy of the parents, they usually have two parents and then inherit their char-acteristic from both parents to make up a new individual. The main operator to work on the parents is crossover, the main searching operator. This operator takes valuable information from both parents and combines it to find a highly fit individual. To apply this operator, two strings from the reproduced population are mated at random and they are cut once randomly between two bits. The new strings are then created by interchanging the tails. It means that parent A will get the tail cut from parent B as it tail and vice versa. This can best be explained by an example. Suppose there are Chapter 2. Genetic Algorithms 11 two strings 00000000 and 11111111 and assume a random number generator comes up with a 3 as the cutting place or crossover site. Then the new strings will be 11100000 and 00011111 Reproduction and crossover give genetic algorithms much of their power. The search is emphasized towards the best and new regions are explored by using information about things that have worked well in the past. 2.6 Mutation Even though reproduction and crossover come up with many new strings they do not introduce any new information into the population at the bit level. They work with the bits that are already in the population and do not get any new bits into the population. The bits can only reproduce or die, so if at certain position all the bits have the same value, there is no way that crossover and reproduction can get the lost bit back. To insure against such a loss and as a source of new bits, mutation is introduced. In the case of binary coding, the mutation operator simply flips the state of a bit from 0 to 1 or vice versa. But it should be used sparingly because it is a random search operator that searches the space randomly and the algorithm is intended to be a randomized searching algorithm, not a random search. 2.7 Mathematical Foundations The theoretical properties of genetic algorithms can be studied using the theory of schemata1 proposed by Holland [19]. The defining length of a schema, 6(h), is the 1see definition of schemata in Section 2.3 Chapter 2. Genetic Algorithms 12 length between its outermost denning positions, for example 6(0 *0) = 3 — 1 = 2 and S(* — 3 — 3 = 0. The defining length is a measure of how often crossover may be destructive for a particular schema. For the schema 0*0 there are two ways to destruct it by cutting it, but the schema * * 0 can not be destructed by crossover, providing both offsprings created by crossover are kept. The order of a schema, o(h), on the other hand, is a measure of how often mutation will be destructive for a schema. The order of a schema is the number of defining positions for a string, for example o(0 * 0) = 2 and o(* * 0) = 1, or in other words, mutation can possibly destruct schema 0 * 0 in two places but schema * * 0 in one place. In other words, schemata with short defining length and low order, stands the biggest chance of surviving into the next generation. This can be written as the Schemata Theorem [19] T h e o r e m 1 Consider a GA using both crossover and mutation. The expected pro-portion of each schema represented in the population changes in one generation from m(h,t) to m(h,t + 1) > m ( M ) : 7 ^ ( l - - ro(M)) (X - *")° < f c ) Where pc is the probability of a particular mating to undergo a crossover and pm is the probability of a single bit to mutate during a generation. The average fitness of the strings at time t representing the schema h is denoted by F(h,t) and F(t) is the average fitness of the population. What it means is that the number of schemata at time t + 1 is greater or equal to the number at time t multiplied by the expected number of offsprings less those schemata that are destructed by crossover or mutation. In other words the schemata theorem states that the algorithm is going to converge towards the best, but there is no guarantee that it is converging to the optimum. As Goldberg [15] puts it (pp.74): Chapter 2. Genetic Algorithms 13 "Convergent behavior without guarantee of optimality bothers many people who approach genetic algorithms from other, more traditional, optimization backgrounds. . . .the fact of the matter is that genetic algorithms have no convergence guarantees in arbitrary problems. They do sort out interesting areas of a space quickly, but they are a weak method, without the guarantees of more convergent procedures. This does not reduce their utility. Quite the contrary, more convergent methods sacrifice globality and flexibility for their convergence." G A have been shown to behave well on multimodal functions, although there is no known necessary and sufficient condition under which a function is genetically opti-mizable. However, numerous studies have shown that functions on which G A fail are pathological, and generally fail to be optimized by any other known technique except exhaustive search [4]. In a recent study by Goldberg [14] it has been shown that even though the algorithm is misled, it will converge for a wide range of starting conditions (initial population) and under unfavorable conditions. 2.8 Example Suppose we have the function (1 - cos2nt) sin2lint (2.5) and wish to find the maximum on the time interval [0,1]. The function has 11 local maxima, the global one being in the middle as shown in Figure 2.3. By not knowing the underlying function itself, but only the values of it, it is very difficult to locate the maximum. If the direction of steepest gradient is followed the maximum will be the one closest to the starting point. The G A on the other hand should be able to find the maximum by climbing more than one peak at a time. Chapter 2. Genetic Algorithms 14 rt rt 1.6 -o. o.e 0 .4 o.e O .B 1.0 Figure 2.3: Function with 11 local maxima Assume the interval [0,1] is coded as 10 bit binary string and the population size is 10. The initial population is chosen randomly and the binary string then mapped onto the time interval. The fitness is read from Figure 2.3. The average fitness is calculated and the number of offsprings for each individual found (see Table 2.2). Crossover and mutation are done by choosing mates and crossover site, both at random. Mutation is applied by mutating every bit of the new population with probability equal to pm (approx. 1/1000), see Table 2.3. After three generations the algorithm is able to find a solution within 3.5% of the maximum, as can be seen in Figure 2.4. It is not until after 12 generations it finds The Maximum, but that is one of the underlying principles of the algorithm that it does its best while learning to do better. 2.9 Summary Because the algorithm works with a population of strings, it is given more chance to locate the global maximum in a multimodal search space. It is in fact searching many points (peaks) in parallel and exchanging information between the peaks. The initial Chapter 2. Genetic Algorithms 15 para- fitness normalized off-meters fitness springs 1000001101 0.513 1.61 2.18 2 0011110101 0.239 0.78 1.05 2 0000001110 0.014 0.00 0.00 0 1100010111 0.773 0.85 1.16 1 1000110111 0.554 0.17 0.24 0 0101011011 0.339 0.80 1.16 1 1111010111 0.961 0.03 0.04 0 1001100101 0.599 1.67 2.26 2 0111010011 0.457 0.01 0.01 0 0101001011 0.324 1.40 1.89 2 av. 0.74 Table 2.2: Reproduction reproduction mate x-site new generation 10-00001101 3 2 1011110101 100-0001101 8 3 1001100101 00-11110101 1 2 0000001101 00111101-01 10 8 0011110111 110001-0111 7 6 1100110101 01010110-11 9 8 0101011011 100110-0101 5 6 1001100111 100-1100101 2 3 1000001101 01010010-11 6 8 0101001011 01010010-11 4 8 0101001001 Table 2.3: Crossover and mutation Generation 3 Generation 1 CTl Generation 14 Generation 12 Generation 11 Generation 9 Chapter 2. Genetic Algorithms 18 population is generated randomly and the population size is kept constant throughout the process. The algorithm only requires payoff information (fitness) for each of the string, without the need for assumptions such as differentiability, thus making it very useful for discontinuous surfaces. Genetic algorithms are inherently parallel. Indeed, all strings or individuals in a population evolve simultaneously without central coordination. To realize their full potential, they must be implemented on parallel computer architectures. Chapter 3 System Identification 3.1 Background Although a variety of techniques have been developed for system identification, none has proven to be effective in all domains. It would be nice to have a method that is sufficiently robust, that is, could be used on a broad class of problems. GAs have been used on a variety of problems as have been reported in [16,17,15,14]. In this chapter the algorithms are going to be applied to both discrete and continuous time systems. But first we look at some previous work in this area. Etter et. al. [11] studied the system below and modelled it as having only two poles. y(t) = 1 + 1 0 g 1 u(t) (3.6) They identified o,\ and a 2 with the input as a white noise and used population size of 11. They showed that the G A did better in locating the true values than a random search did. Das and Goldberg [8] worked with system of the form M O = ^ £ ^ £ ( « ( < ) + <(«)) (3.7) 1 + a^q 1 + a2q with 60 = —0.2, bi — 0.1, b2 — 0.4, at = —1.6 and a 2 — 0.95. The system they used was non-minimum phase and very oscillating (£ = 0.04). They successfully identified the five parameters with the input as a P R B S signal and e(t) as Gaussian noise with variance equal to 10% of the input. 19 Chapter 3. System Identification 20 Smith and DeJong [22] used G A to calibrate a nonlinear model of US. migration patterns. There has been one application of G A for continuous time systems. Goldberg [12] identified mass-spring system with small damping (£ = 0.05) mx(t) + cx{t) + kx{t) = f(t) (3.8) where m — 1.0, c = 0.1 and k = 1.0. The force function, /(<), was a two step staircase function and he identified directly the parameters of the continuous time system, m,c and k. A l l the applications so far have been on open loop systems and for the discrete time systems have identified the parameters of the models which RLS can easily do. Nobody has seen the ability of the G A to identify directly the poles and zeros. When estimating poles and zeros with conventional estimation methods the problem is that the system is no longer linear in the parameters. Standard algorithms do not identify directly the poles and zeros. They change the system into a concatenation of second order systems and then calculate the poles and zeros for each 2nd order block [24]. GAs on the other hand can directly identify the poles and zeros. There is really no difference from GA's point of view whether it is identifying the poles and zeros or the parameters. A l l it needs is a fitness value to assign to each string. The advantage of knowing the poles and zeros is simpler controller design as can be seen in Chapter 4. GAs could also be used to identify physical parameters, like Goldberg did in [12] for a mass-spring system. For instance, using this method the friction coefficient in a motor drive could be identified directly. Traditionally discrete time estimation is used which results in coefficients that are nonlinearly dependent upon the sampling time. Because of GAs ability to deal with nonlinearity they can be used to identify continuous time systems. Chapter 3. System Identification 21 The simulations were performed on Pulp and Paper Centre /nVAX. Programs were written in P A S C A L for the Genetic Algorithm and in F O R T R A N for the RLS part. 3.2 Discrete time identification Consider the system A{q-')y{t) = 2?( ( « - d) + C{q-X)e{i) (3-9) Where A,B and C are polynomials in the backward shift operator, q'1, i . e. y(t — 1) = q~1y(t) and y,u and e are the output, input and noise respectively. The noise e(t) is a normally distributed random sequence with zero mean and a unit variance (cr^). The polynomials A and C are assumed to be monic. The objective is to estimate A(q~1), B(q~1) and the delay d, when given the input u(t) and the output y(t). The estimates are denoted by Two sequences e(t) and n(t), can be defined, for calculating how well the estimates fit the system, as: Mq-'M*) = Biq-'Ht - d) + e(t) (3.10) or v(t) = y(t)-m (3 . i i ) with Aiq-^yit) = Biq-'Ht - d) Then we try to minimize i?[e2(f.)] or E[n2(t)}. The first case corresponds to the least-squares case and has a search space which is quadratic, the second is akin to the Instrumental Variable (IV) case and has a highly nonlinear search space. Depending on the method used, the fitness function is chosen as = X > "(*(«-*))' (3-12) i=0 Chapter 3. System Identification 22 G e n e r a t i o n s Figure 3.5: Window size = 5 Figure 3.6: Window size = 10 or as F(t) = Y,M-(v(t-i)Y (3.13) t = 0 where M is a bias term needed to ensure a positive fitness as explained in Chapter 2 and w is the window size or the number of time steps the fitness is accumulated over, with a effect akin to that of the forgetting factor in R L S . The effect of different window sizes can be seen in Figures 3.5 to 3.8 where the algorithm is run on a system with P R B S input and colored noise. For the moment just assume that the figures show some parameters of a second order system. From these figures it can be seen that the variance of the parameter estimates reduces as the window size increases. But there is a price to pay for increasing the window size. Implementing these fitness functions is expensive in terms of C P U time. Because of the nature of the algorithm, (i. e. coding, probabilistic transfer rules, etc.) no recursive version of the fitness function exists. So at every generation the algorithm has to calculate the estimated output for the whole window which makes the execution time proportional to the window size. There is also an advantage of not having a recursive fitness function. That means that Chapter 3. System Identification 23 B < ~ •, '<) f i r . ;"i ! : ; ;; i n . j ; , . : •-. : : <\X i i i : rt\:r' i r : . ! ! ! ; ! :IL 0.7170 : ' ! ! ij ' i ' 1 , i ' i " 'I '' i *|; i: i i ;i ' I'i'ii*-| i !?? 1 '* :| J. ^ ?• ':: i i j ' ^ ' • •; ': * -' ; : : : \i t(- •'. i; : ; • "•: • i ( *• • • i * • 3 "' i ' value 'o i' 1 ! J f M h ; Estimated i j i i i i l i i t j i | W^il 1 rljUlUir-Kl •0.7010 Irtj Li L J | f_|j ; :" i •.0OOOEO3 G e n e r a t i o n s Figure 3.7: Window size = 20 1B0. SOO. G e n e r a t i o n s Figure 3.8: Window size = 30 the algorithm does not have to wait for new input-output data before coining up with new estimates. It can actually iterate as often as one likes for each sample but there is some upper limit because of time constraints and in the window the input has to be persistently excited (see Chapter 4, Section 4.4.4) for the algorithm to converge to a certain value. Figures 3.9 to 3.11 show parameter estimates for different number of trials (generations per sample) for a system with P R B S input and no noise. It is seen that the algorithm actually needs fewer generations (200,175,150) to converge as the number of trials increases from 1 to 3 and hence number of data points (200,87,50). But there should be some limit on number of trials as can be seen from Figure 3.12 where the algorithm uses 30 trials for each data it needs about 2250 generations to converge or 75 samples. Chapter 3. System Identification 24 G e n e r a t i o n s Figure 3.9: Number of trials = 1 Figure 3.10: Number of trials = 2 - -7.1IMOC-03 G e n e r a t i o n s G e n e r a t i o n s Figure 3.11: Number of trials = 3 Figure 3.12: Number of trials = 30 Chapter 3. System Identification 25 3.2.1 Parameter identification The system of Equation 3.9 can be described by the following polynomials A{q~l) = 1 + a l 9 - 1 +--- + anq-n B(q-') = 6 0 (l + 6 i g - 1 + --- + M - n ) ( 3 - 1 4 ) C{q-1) = 1 + c1q-1 + • • • + cnq~n The G A can be used to identify the parameters in A and B and the delay, using either Equation 3.12 or 3.13 as a fitness function. For a second order stable system it gives a search space for ai and a2 of the form seen in Figure 3.13. If the system is also inversely stable the search space for f>i and b2 is of the same form too. 3.2.2 Pole-zero identification Because they do not require linearity in the parameters, genetic algorithms can directly identify the poles and zeros of the system. In pole-zero form, the plant can be written as: (3.15) Biq-1) = b0{l-ziq-1){l-z1q-1)-..{l-zmq-l){l-zmq-1) Where m = n/2 if n is even and m = (n + l) /2 if n is odd. The parameters, p m and zm will be zero if n is odd. It can also be reparameterized so that a complex conjugate poles or two real poles will be represented by two parameters. A{q-i) = ( l - ( a 1 ± / 3 1 ) g - 1 ) - . . ( l - ( a m ± / ? m ) g - 1 ) (3.16) B{q~') = M l - ( 7 l ± t f l ) g - 1 ) - " ( l - ( 7 m ± * n x ) ? - 1 ) The parameters Pi and Si can be either imaginary (complex conjugate poles) or real (two real poles). Because the signs on (3 and 8 are of no importance we can use the signs to decide if the numbers are imaginary or real, negative will mean complex number and Chapter 3. System Identification 26 a 2 •2 -2 \ -2 2 i?e _ , _ _ , Figure 3.14: Poles and ze- Figure 3.15: Reparameter-rigure 3.13: Parameters . ., , , . , , , ros in the complex plane lzed complex plane (3.17) positive will mean real numbers. As an example l + 2 9 - 1 - f 0 g - 2 = ( l - ( _ l + l ) g - 1 ) ( l - ( - l - l ) 9 - 1 ) = [-1,+1] l + 2 g - 1 + 2 g - 2 = {l-{-l+jl)q-i){l-(-l-ji)q-i) = [-1,-1] That gives search space for a stable system of the form seen in Figure 3.15. Where the lower half plane excluding the real axes, represents the complex conjugate poles and the upper half plane represents the real axes. If the parameters for a second order system are given by [23] (see Figure 3.13) : Aiq-1) = l.O-l.bq'1 + 0.7q~2 B(q-X) = 1.0(1.0+ 0.5g- 1 + 0.0g- 2) C(q-X) = 1 .0- 1.0?"1 + 0.2g- 2 The poles and zeros are (see Figure 3.14) : A\q~x) = 1.0 - (0.75 i j O . 3 7 ) ? - 1 Biq-1) = 1.0(1.0- (-0.25 i-0.25)?- 1 ) or (see Figure 3.15 ) : ' 0.75 + j'0.37 (3.18) (3.19) pit2 = [0.75,-0.37] = p x , 2 = 0.75 - j'0.37 (3.20) Chapter 3. System Identification 27 and the zeros f -0.5 z l i 2 = [-0.25,+0.25] = zi,2 = < (3.21) ( 0.0 For a stable minimum phase plant, the poles and zeros are inside the unit circle (Figure 3.14), therefore the search space can be limited to be the unit circle or a box enclosing the unit circle. For a nonminimum phase plant some of its zeros will be outside the unit circle, so one has to decide how big the search space is going to be depending on a priori knowledge of the system. 3.2.3 Results Parameter settings The crossover rate is chosen so as to give some of the population the opportunity to survive into the next generation without any changes. The mutation rate is chosen such that on the average one string in the population is mutated. Unless stated otherwise the genetic parameters have therefore been chosen as follows [8] : Pc = 0.8 pm = 0.01 (3.22) population — 100 Second order systems were used so six parameters needed to be identified, that is d and 60 and then either parameters, bi,&2> ai and a 2 or poles-zeros, a i , / ? i , 7 i and h\. The delay is coded as two bit string to give 4 choices for the delay and the other parameters are coded as 7 bit strings, making totally a 37 bit string, which leaves the search space with 2 3 7 = 1.37 1 0 n alternatives. The parameters have been concatenated as follows M|o a |o 2 |6 1 |6 2 |6 0 | (3-23) Chapter 3. System Identification 28 or M |ai|/?i|7il*ilM (3.24) depending on whether poles-zeros or parameters were identified. Upper and lower bound on the parameters are defined (see Figures 3.13 and 3.14) and the resolution of the coding is calculated using Equation 2.1. lower bound upper bound # of bits resolution d 1 4 2 1 b0 0.0 2.0 7 0.016 CLi,bi -2.0 2.0 7 0.032 0-2, b2 -1.0 1.0 7 0.016 ai,/?i»7i»*i -1.0 1.0 7 0.016 Table 3.4: Search space for a stable minimum phase system Ident i f i ca t ion w i t h PRBS A P R B S signal is used as an input for the system of Equation 3.18. The P R B S input has a period of 127 with the bit interval equal to four times the sampling interval. The G A is run for 600 generations and 3 trials are used for each input-output data so 200 samples will be used. The window (forgetting factor) has been set as 30, i . e. the fitness function is calculated for the current input-output and the 30 previous samples. Figure 3.16 shows the output for the P R B S signal when there is no noise in the system (cr2 = 0). Figures 3.171 and 3.18 show the estimated parameters for each generation using instrumental variable criterion as given in Equation 3.13. The values of the estimates of the last generation are written at the right hand side of the graphs. After about 150 generations or 50 samples the algorithm comes up with unbiased estimates for all parameters except the zeros. Parameters —ji and 8i should both be 0.250 so a bias of 0.321 for 6\ seems rather big. But if bi and 62 are calculated we Estimation of — Q i , / ? i , — 71,6X Chapter 3. System Identification 29 3 3 O C3 O.OO 6 0 . l O O . 160 . N u m b e r o f s a m p l e s Figure 3.16: P R B S input and output of a system without n o i s e in T3 - • 0.2130 •7 . I000E-02 ISO. 3 0 0 . 4 6 0 . G e n e r a t i o n s Figure 3.17: Pole-Zero estimate of a system without noise Chapter 3. System Identification 30 Chapter 3. System Identification 31 get 0.426 and 0.050 respectively (true values 0.50, 0.0). The steady state gain of the system is 7.5 so Equation 3.13 is less sensitive to changes in the zeros than the poles and it should also be emphasized that the algorithm does not necessarily converge to T H E optimum. "It does its best while learning to do better." Figure 3.19 show how the poles and zeros move around in the complex plane for each generation where the initial generation is at the back and the final generation is in front. The unit circle is also plotted every 100 generations together with the estimates at that point. 3.3 Recurs i ve Least Squares es t imat ion The actual process, or the system, is assumed to be described by an equation of the form where e(t) is a white noise with zero mean and a variance o~\ and the polynomials A(q~1), B(q~1) and C(q~1) are given as A(q->)y(t) = B(q-*)u(t) + Ciq-'Mt) (3.25) Aiq-1) = 1.0 + a jg- 1 + • • • + anaq-n« Biq-1) = q~l( 1.0 + b i q - l + -.- + bnbq-n>) Ciqr1) = i.o (3.26) Now introduce the parameter vectors 9 and 6 and a vector, 50. 100. 100. N u m b e r of s a m p l e s Figure 3.21: Parameter identification using R L S Chapter 3. System Identification - 2 J Figure 3.22: Parameters locations > -o CO s 7 . 1 0 0 0 E - 0 2 Figure 3.23: G A , pole zero identification Chapter 3. System Identification 35 Figure 3.25: GA, parameters locations calculated from the pole zero identification Chapter 3. System Identification 36 Figure 3.27: GA, parameters locations Chapter 3. System Identification 37 steps (0.93 0 = 0.04) and a 1 ; a 2 , &i and b2 are then identified.2 Figure 3.21 shows the result using the input shown in Figure 3.20. It can be seen that the estimates have rather large variance especially b\ and b2 (the long dashed and dotted line respectively). Their value after 200 samples is 1.3957 and 0.6088 respectively but 10 samples before their values were 0.4214 and 1.2141 respectively so it is difficult to say what value they are converging to. The estimates for ai and a2 have also a variance but much smaller and they converge to biased estimates with the final estimates as -1.1543 and 0.4974 respectively. Figure 3.22 shows plot of a2 as a function of and b2 as a function of bx with the time axis running out of the page and the triangle for a stable estimates plotted every 50 samples. To compare those results with the G A , the G A is run identifying the same parame-ters as those identified by the R L S . That means that the gain, b0 and the delay, d, are assumed to be known so there are only four parameters to be identified. Using same parameter settings as in Table 3.4 it gives a total string length of 28 bits, so the popu-lation size has been set to 50. The G A is run twice, first identifying the poles and zeros and secondly identifying the parameters. The results of the pole-zeros identification is shown in Figure 3.23, it is then converted into the parameters, Figure 3.24 and a 3-D figure is plotted, Figure 3.25. The parameters estimation is shown in Figure 3.26 and the corresponding 3-D figure is shown in Figure 3.27. It can be seen that in both cases the poles have almost zero bias, they are only limited by the resolution of the search space. They converge in about 50 generations for the pole-zero identification but in about 100 generations for the parameters identification or about twice as fast for the pole-zero identification than for the parameters. The zeros converge slowly for both cases but the final estimates are close to the true values (0.5 and 0.0) in both cases. If G A is then compared to the RLS it can be seen that the RLS needs more than 50 2True values from Equation 3.18 are -1.5, 0.7, 0.5 and 0.0 respectively. Chapter 3. System Identification 38 samples for the poles to converge but the zeros do not converge, whereas the G A needs between 50 and 100 generations for the poles to converge which means that with 3 trials per sample it needs between 17 and 33 samples and the zeros are slowly converging. So in terms of number of samples the G A converges faster. But as mentioned earlier the fitness function for the G A can not be calculated recursively so the algorithm has to calculate the outputs for all the window and to calculate every output it involves (?ia + rif, + 1) multiplications and (n a + rif,) additions. The difference in bias of the estimates are mostly caused by different objective (cost) function, the RLS uses a simple least square whereas the G A uses IV alike objective function. 3.4 Continuous time identification Consider n-th order system with a differential operator s = ^ and unknown coefficients a; and bi y(t) = V " + 6 n ' ii»i»»!! I'! J ' 11' ' • " 1 1 . 1 ' I « " ' I I 11, n;;II r i . . . . i ' I. i i i i i i i i' !• ii 'i I iiI, n I " i" 'i i' i' ' i " j , i , i' i i i " i i' *i H I' II «i"iii i "I'II n mi »'•.'! I'll1!! II Mil II H ' ! i l .' , 11' '11' 1 !i i'",'' ii ' i .1 Ii ii i ' ' i ' ' i , i ' ' 1 1 ' ' 1 ' 1 1 ' ' ' ' i i ' i ' i ' i H .1 i pi I. i i i' s ii i ' i n ii ', " i' ' i i i h ii »'!!!!-':'!!!!!..!!!!.!!!!! !! ' ! ' II1 » i ' i ! ! : •' 'I 1 / i n i! i i 1" ii n1 ll 'I " I • I'l". " ' ' ! ' ' 111 " i 1 " i' I' ii 'j i' i"'i:i'!i'i ji "in1:1;;«: •'« I f ' I1 'l I V 'i ,I'II'I ,^,VIII|I'I1' 'i ^'!::i;'!!'.'i!':>•'!'.' I I " :i ;i nr.i >: ii ' i I * *!! l! \i)'t '•!!! "> *)V< 51" 1 i " If!'"!: i1 1 . 1 1 1 /1 O.O 10 . SO. 3 0 . 4 0 . 6 0 . 6 0 . TO. 8 0 9 0 . 1 0 0 . T i m e Figure 3.28: Continuous-time-system input and output discrete time case, the algorithm can be run several times for each sampling time. 3.4.1 Results This has been implemented in ACSL using the PASCAL subroutines from previous implementations of the algorithm. A second order system has been used with the transfer function: y(t) 0.0/1+ 1.0 u{t) s2 +0.5s+ 1.0 [ ' The GA is used to find all four parameters of the plant. The search space has been defined as in Table 3.5. The total string length is 36 bits so the GA parameters have lower bound upper bound # of bits precision oi, a2,h,b2 0.0 12.775 9 0.025 Table 3.5: Search space for continuous time parameters Chapter 3. System Identification 40 12. O. ' ' 1 1 • 1 ' ' 1 1 ' ' 0.0 60. lOO. ISO. 200. BBO. 300. 360. 400. 460. BOO. B60. 600. G e n e r a t i o n s Figure 3.29: Continuous-time parameter estimates 0.4 0.3 0.2 0.1 O.O — . 2 -- .3 I 1 1 1 1 1 1 1 1 . 1 o.o to. eo. 30. 40. BO. eo. 70. eo. 90. 100. T i m e Figure 3.30: Actual output and the output using the final estimates Chapter 3. System Identification 41 been set the same as in Equation 3.23. The input u(t) has been chosen as (see the dashed line in Figure 3.28) : u(t) = 0.08 sin 0.5* + sin 3.0* (3.37)' The input is chosen so as to excite the system both above and below the natural frequency and they should have about the same amplitude in the output for the pa-rameters to converge to the true value. The sampling time T, is 0.5 seconds and the window is chosen as 8 seconds. The parameter estimates are shown in Figure 3.29. The convergence is slow but they approach the actual values and the estimate after 450 generations, i. e. 150 sampling intervals is 0.000. + 1.050 s2 + 0.6005 + 1.025 ' j which has a natural frequency, a;n = 1.01 and a damping ratio, £ = 0.30 instead of 1.0 and 0.25 respectively. If both the actual output and the output using the final estimates are plotted (Figure 3.30) one can see that the response is almost identical excluding the transient. 3.5 Friction compensation Now look at the estimation of physical parameters, namely the friction of a motor. A motor with friction torque T/ and a load disturbance torque TJ can be described by the following model, J ^ = KI{t)-Tf(i) + Tl(t) (3.39) where J is the total moment of inertia, K is the current constant and / is the motor current. Neglecting the load disturbance and introducing m = «(0 + ^ (3-40) Chapter 3. System Identification 42 Figure 3.31: Friction model the motor model can be written as J ^ = Ku(t) + {ff(u>)-Tf(u,)} (3.41) If the estimates are good the term inside the bracket is going to vanish and the model looks like a frictionless motor. Many models for the friction have been suggested but a model used in [5] has been adopted. The model is: aiu> + / ? i iv > 0 z>(o,) = { (3.42) a2iv + j32 iv < 0 Therefore the estimation of Tf(iv) requires estimation of four parameters, on, a2 and f32. Only two of them can be estimated at a given time depending on whether the angular velocity iv is greater than zero or less than zero. G A can be applied to this problem ones the objective function has been defined. Assuming it is required to minimize the error between the actual angular velocity and the estimated one, the fitness function becomes minj (u>(t) - co -»-> ea s Figure 3.33: Friction parameters identification sample which is sampled at the rate of 5 per second. The input is a square wave with period equal to 12 seconds, which can be seen together with the output in Figure 3.32. Figure 3.33 shows the estimated parameters using G A and the final value that the G A comes up with is: dj = 0.100 0X = 0.400 (3.45) d 2 = 0.205 -32 = 0.180 which have almost zero bias. The algorithm takes less than 100 generations to find approximate estimates for each one of the populations and further refinements are then found along the way. 3.6 S u m m a r y In this chapter it has been demonstrated how GAs can be used to estimate both con-tinuous and discrete time systems and for identifying parameters, poles and zeros or Chapter 3. System Identification 45 physical parameters of a system. The G A has proven to be a. robust algorithm, whereas in all the applications the basic algorithm (procedures) stays the same. The only dif-ference between these different applications is a different routine to calculate the fitness function. In all the applications shown in this chapter, the G A has been able to converge towards the actual values of the parameters. In most cases it gives unbiased estimates but in some cases it takes time, like the identification of the zeros, primarily because the objective function is not as sensitive to changes in the zeros as changes in the poles and the G A is only looking for a good solution not necessarily the best. In comparison to some widely known identification technique, RLS, the G A perform as well or even better in terms of number of samples required to converge. But G A can also easily be used on problem where RLS is difficult or even not possible to use because of the requirement of linearity in the parameters of the system. C h a p t e r 4 Con t ro l l e r design 4.1 C o n t r o l l e r There are a large variety of adaptive design techniques. In this thesis I have chosen to use an indirect scheme and I am not identifying the stochastic part, so an adaptive pole placement design has been chosen. It is a simple design method that makes use of the knowledge about the poles and zeros to obtain a desired transfer function or desired response. A SISO plant is described by an A R M A X 1 model: *>=S-(<)+3£(<) (4-46) The control law for a two-degree of freedom pole-placement controller (Figure 4.34) can be written as R(q)u(t) = -S(q)y(t) + T(q)yT(t) (4.47) where yT is the reference signal and R is assumed to be monic. To simplify the writing in the analysis that follows, the arguments of polynomials are suppressed. If Equation (4.47) is written in terms of the system input u(t) and then put into Equation (4.46) the closed loop system becomes TB . . qdRC , . / xAuto Regressive A(q), Moving Average C(q), eXternal signal B(q) 46 Chapter 4. Controller design 47 The desired closed loop transfer function is given by H(q) = f=g (4,9) or TB B m q*AR + BS Am ( 4 - 5 0 ) By choosing the desired closed loop transfer function and estimating the plant, the controller design has been reduced to finding the polynomials R, S and T that satisfy Equation 4.50. Some of the process zeros can be cancelled in the design. If B is factorized as B = B+ B~, where B+ is cancelled stable process zeros and B~ is uncancelled process zeros, Bm must be written as, Bm = B~B'm The parts of B that are not factors of Bm must be factors of qdAR + BS so B+ must be a factor of R. Therefore R must be written as R = B+R. Generally the degree of qdAR + BS is higher than Am so there must be some cancellation on the transfer function. So the Diophantine equations becomes qdAR + B~S = A0Am (4.51) and T can be found from T = B'mAot0 (4.52) where tf0 is to ensure proper gain and A0 is the observer polynomial that is cancelled in the transfer function [1]. In order to design the controller we need to solve the Diophantine Equation 4.51 for R and S. In order to find a unique solution the degree of S has to be less than the sum of d and degree of A. We choose degS = deg A + d — 1. Using the causality conditions (degR > degS) the degree of A0 can be found from Equation 4.51 to be degA0 > 2(degA + d) - degAm - degB+ - 1 (4.53) Chapter 4. Controller design 48 Vr T + o 1_ R U B A Figure 4.34: Two-degree of freedom controller The degree of R can be found from Equation 4.51 as degR = degA0 + degAm — deg A — d (4.54) This procedure assumes that the true plant A, B is known. When the true plant is not known, one can use the G A to do the estimation of the plant A, B and then design an indirect adaptive control scheme as shown in Figure 4.35. 4.2 Parameter based design By assuming that we have knowledge of A and B~ and that A0 is chosen we can solve Equation 4.51 for R and S. Define n = deg A, m = degB~, fc = degA0 and j = degAm, then A(q) = qn + a1qn~1 + • • • + o n B~{q) = b0(qm+b1qm-1 + --- + bm) Ao{q) = g f c + a i 0 g f c _ 1 + ••• + «fc0 Am(q) = q3 + aiwq3"1 H h ajw S(q) = s0qn*d-1 •\-s1qn+d-2 + • • • + J n + r f - i R(q) = q><+J-»- 0} Figure 4.40: Pole-Zero estimates for a minimum phase system without noise •a 900. G e n e r a t i o n s Figure 4.41: Estimates of gain and delay for a minimum phase system without noi Chapter 4. Controller design 58 2-1 -2-J Figure 4.42: Pole zero locations Chapter 4. Controller design 59 fx o ex -2. ' ' 1 1 " ' — BO. 0.00 BO. 100. 160. 200. N u m b e r o f s a m p l e s Figure 4.43: Reference input and output of a minimum phase system with noise using pole-zeros estimates Figure 4.39 shows the reference input and output of the plant when there is no added noise in the system ( & o p. N u m b e r o f s a m p l e s Figure 4.51: Reference input and output for a nonminimum phase syst em CO -B .0000E-03 O.O ISO. 300 . 4 B 0 . eoo . Generat ions Figure 4.52: Parameter estimate for a nonminimum phase system Chapter 4. Controller design 67 Figure 4.54: Pole zero locations for a nonminimum phase system Chapter 4. Controller design 68 3 -60. 0.00 60. 100. 160. 800. 260. 3O0. 360. 400. N u m b e r o f s a m p l e s Figure 4.55: Input-Output for 3 parameters estimate with unmodeled dynamics samples and the search space is defined as in Table 4.7. The desired pole was set lower bound upper bound # of bits precision d 1 4 2 1 bo 0.0 2.0 7 0.016 ai -1.0 1.0 7 0.016 -2.0 20.0 9 0.043 Table 4.7: Search space for unmodeled dynamics at 0 (dead-beat). For the first simulation the delay d is assumed to be known so only the gain, the pole and the zero (b0, a x and bx) are identified. That gives a response that has about 70% overshoot and has a low damping as shown in Figure 4.55 and the parameter estimates are as shown in Figures 4.56 and 4.57 with the estimates for the zero not converging to a certain value. But when all the four parameters (b0, 6 1 } d and aj) are identified (population size = 50, total string length = 23) the overshoot is reduced to about 30% and the damping is higher, Figure 4.58, and the system is Chapter 4. Controller design 69 BO. > to B 'V' * >' S U V " ! ! !l Generat ions Figure 4.56: Parameter estimate for unmodeled dynamics 8 . 4 0 0 0 E 0 2 Figure 4.57: Gain and delay estimate for unmodeled dynamics Chapter 4. Controller design =3 1=5 o -»-> CL, 0.00 60. lOO. 160. 200. 260. 300. 360. N u m b e r o f s a m p l e s Figure 4.58: Input-Output with dead beat control 13 ro 6 Generat ions Figure 4.59: Parameter estimate for unmodeled dynamics with dead beat control Chapter 4. Controller design 71 82. M —» to 1 1 1 1 G e n e r a t i o n s Figure 4.60: Gain and delay estimate for unmodeled dynamics with dead beat control 3 ex o Q L , a - e o . o o o e o . 100 . i s o e o o . e s o . a o o . 3GO. 4 0 0 . N u m b e r o f s a m p l e s Figure 4.61: Input-Output with desired pole = 0.7 Chapter 4. Controller design C D "(0 t> -O C O -4-> s Generat ions Figure 4.62: Parameter estimates for unmodeled dynamics with desired pole •si •B5 eoo. G e n e r a t i o n s Figure 4.63: Gain and delay estimate for unmodeled dynamics with desired pole Chapter 4. Controller design 73 3 QH H-» 3 o "3 QH -60. 0.00 60. 100. 160. ZOO. Z60. 30O 360. 30 N u m b e r o f s a m p l e s Figure 4.64: Reference input and output of a system with window size estimated to be (see Figures 4.59 and 4.60) 0.599g- 2(l + 0.798g-1) 1 - 0.953?-1 Further improvement is obtained if a slower response is chosen (desired pole at 0.7), Figure 4.61. The estimated model becomes (see Figures 4.62 and 4.63) (4.84) 0.662g- 3(l + 0.712g-1) 1 - 0.937?-1 (4.85) 4.4.4 Persistently exciting signal To show the effect of a persistent excitation, the minimum phase system of Equation 4.79 without noise ( ca s co E d 800. Generat ions Figure 4.65: Parameters for a window size = 30 -a -65 CM G e n e r a t i o n s Figure 4.66: Gain and delay for a window size = 30 Chapter 4. Controller design 3 OH -l-> o 3 QH — lOO. —60. O.OO 60. lOO. 160. 200. 260. 300. 360. 400. N u m b e r o f s a m p l e s Figure 4.67: Reference input and output of a system with window size > G e n e r a t i o n s Figure 4.68: Parameters for a window size = 60 Chapter 4. Controller design 76 fz. eoo. G e n e r a t i o n s Figure 4.69: Gain and delay for a window size = 60 starts to deteriorate and the algorithm is not able to bring them back consequently the output does not follow the input very well (Figure 4.64). When the window is increased to 60, to be able to include a step change in the window at all time, the estimates are shown to converge to the true value (Figures 4.68 and 4.69). Because of a small bias in the estimates the output has a small overshoot at every setpoint change (Figure 4.67). 4.4.5 Recursive Least Squares To compare the G A to some method that is widely known, a standard RLS algorithm is used on the same system as before with noise variance o\ = 0.1. The same pole placement controller design is used. A deadbeat observer is chosen and the closed loop poles and zeros are set at zero (deadbeat). The forgetting factor is set to 0.9 to resemble a window for the G A of 30 steps (0.9 3 0 = 0.04) and a 1,a 2,fe 1 and 62 are then identified. Figure 4.70 shows the reference input and the output of the system. It can be seen that Chapter 4. Controller design 77 3 ex -4-» a o a OH T i m e s t e p s Figure 4.70: Reference input and output of a system using RLS CD J 3 !> T 3 CD CO s 2.8200E-02 e o . i o o . 160. N u m b e r o f s a m p l e s Figure 4.71: Parameter estimates for RLS Chapter 4. Controller design Figure 4.72: Parameter locations for RLS Chapter 4. Controller design 79 3 -BO. 0.00 60. lOO. 160. 200. N u m b e r o f s a m p l e s Figure 4.73: Reference input and output of a system using G A to compare to RLS the overshoots becomes larger as the parameter estimates deteriorates (Figure 4.71). Figure 4.72 shows plot of a 2 as a function of a x and b2 as a function of b\. To compare those results with G A , the G A has been run for same parameter estimates (only pole and zeros, not gain and delay) using IV criterion and the results are shown in Figures 4.73, 4.74 and 4.75. Comparing Figure 4.70 and Figure 4.73 one can see that there is not much of a difference in their response. Both have transients while searching for good parameters and after few steps output follows the input. The RLS is quicker to converge to a value whereas the G A is satisfied with suboptimal estimates. 4.5 Summary It has been shown how knowledge of the plant (.4 and B), either its parameters or poles-and-zeros, can be used to design a pole-placement controller. Simulations results Chapter 4. Controller design 80 3 cu .*-> eo 0 m 1BO. 300. 400. Generat ions Figure 4.74: Parameter estimates for G A to compare with R L S Chapter 4. Controller design 81 show that the G A is as well fit for doing the identification as the RLS. The GA has proven to be able to handle both minimum and nonminimum phase systems and has also shown its ability to control when there is unmodeled dynamics. Chapter 5 Experiment 5.1 Water level in a tank A n experiment was carried out on a tank system at the Pulp and Paper Centre. The tank has a sensor to measure the height h of the water and a pump to pump the water, given a drive voltage u, into the tank. The outflow of the tank is a function of the tank level so the dynamics will be nonlinear. Therefore the tank can be described by a nonlinear differential equation of the form [2]: ^ = -AVh + Bf(u(t)) (5.86) at Where A depends on gravity and the ratio between the effective outlet area and the cross section of the tank and B depends on the cross section of the tank and also relates the pump flow to the drive voltage u of the pump motor electronics. A linear model of the tank is given in A s t r o m and Ostberg [2] as: KT His) = (5.87) K ' Ts + l K 1 where T depends on A and the initial height and K depends on the sensor and the constant B. A Z O H is used together with the A/D converter to read the water height so the Z-transform would be: „ , , KT(1 - e-%)z~1 H ( z ) = — r IT-T— (5-88) (1 — e T z~1) where T, is the sampling time. 82 Chapter 5. Experiment 83 10. 80. 30. 40. 60. 60. 70. N u m b e r o f s a m p l e s Figure 5.76: Tank Input-Output 5.2 Simulation results To collect the data the sampling time is chosen as 0.55 sec. and 80 samples are obtained using P R B S as the input (see figure 5.76). Because of the prohibitive time it takes to run the G A on an I B M A T we were not able to do any online control on the tank. The G A is therefore run offline using the IV fitness function (see equation 3.13) to identify directly the poles and zeros. The algorithm assumes that there are two poles, one zero and it also identifies the gain b0. il+Piq-W+Piq-1) [ j The poles are decoded as in chapter 3. Both the poles and the zero are assumed to be stable so they are assumed to lie between -1 and +1. The gain is assumed to be in the range [0,10] The length of each string is chosen as 11 which give resolution of about 1/1000 for the poles and the zero and about 5/1000 for the gain. With four Chapter 5. Experiment 84 cu > cu -•-> 0 co G e n e r a t i o n s Figure 5.77: Parameter estimate for a tank .s 0.8000E-02 G e n e r a t i o n s Figure 5.78: Estimated gain for a tank Chapter 5. Experiment -2-" Figure 5.79: Pole zero locations for a tank Chapter 5. Experiment 86 parameters to identify, the total string length is 44 bits so the population size is set to 100. The probability of crossover and mutation is chosen as before to be 0.80 and 0.01 respectively and 6 generations were generated each sampling interval. Figure 5.77 and 5.78 show the estimated parameters for each generation using the input output data of figure 5.76. The estimated system after 300 generations is: o.oesg-Mi-o.eeig-1) o.oesg-1 (1 - 0.637g-1)(l - 0.965-r1) (1 - 0.965-?-1) [ ' ' So the zero cancels one of the poles and the system is a first order system with a delay and a pole close to the unit circle. The time it takes to fill up the tank is about 10 seconds so with a sampling time close to half a second the pole should be according to equation 5.88 about -0.95 so the estimates seems good. In figure 5.76 the output of the tank is shown if the final estimated parameters were used (dashed line). It can be seen that the estimated output is not far from the actual one and it should be pointed out that the estimated output is put equal to the actual output at the beginning of every window, which means that the two output are much closer than suggested in figure 5.76. In figure 5.79 the locations of the poles and the zero are shown in the complex plane for every generation where the cancellation of one of the poles with the zero can be clearly seen. Chapter 6 Conclusions In this thesis a new approach was taken to the identification problem. The usual hill climbing algorithms that follow the steepest gradient were abandoned for a method that uses concepts from evolutionary theory called Genetic Algorithms. They proved to be able to identify both discrete time and continuous time systems and could give unbiased estimates in the presence of colored noise. They showed some advantage over R L S and could be used in cases were the system is not linear in the parameters were RLS can not be used, for example identify physical parameters, delays and pole-zeros. They were used to design an adaptive pole-placement controller and gave good control for a variety of problems as demonstrated by simulations. A n experiment was presented. The algorithm was tested on a real data from a tank system. The algorithm behaved well but because of the prohibitive time it takes for the algorithm to run on an I B M - A T , no real time online control was attempted. Genetic Algorithms have proven to be useful on wide range of applications without any changes in the basic algorithm. The only interface with the system the algorithm is working on is through the objective or fitness function. That function is the only thing that needs to be changed from one application to another. Because the GAs search within a population not from a single value they are insensitive to noise. As with the R L S , there must be some priori knowledge about the system to identify, the search space that the parameters are likely to lie within must be specified and also the resolution. For a proper choice of the resolution the algorithm will prevent the 87 Chapter 6. Conclusions 88 estimates from jumping around and hence could be used to filter some noise. It should also be remembered that the algorithm is a randomized search technique, so there is no guarantee of optimality, the algorithm does its best while learning to do better. A n area for further research is the exploitation of more than the best string in the population for the design of a robust controller. For example the average of the ten best strings in the population could be used for preventing abrupt changes in the estimates. Also dominance could be used for a changing plant and for a multimodal search space, like the example from chapter 2, some sort of distribution among the peaks could be maintained by introducing sharing, that is the individuals are prevented from crowding around one particular peak by punishing them for being too close together. That could be particularly useful in changing environment where by maintaining diversity in the population the algorithm does not put all of its effort into searching around a particular peak. Finally, GAs are parallel algorithm, so every attempt to run the algorithm on non-parallel computer is bound to be slow. For our case the algorithm uses little bit less than 1 second of C P U time for each generation, on a /A-VAX (1 MIPS) for population size of 100, string length of 37 and window size of 30 steps. Once parallel computer architectures become readily available, G A will become more attractive. Bibliography strom, K . J . and B . Wittenmark, (1984). Computer controlled systems, Prentice Hall Inc., Englewood Cliffs N . J . [2] Ast r6m, K . J . and A . B . Ostberg, (1985). "A teaching laboratory for process control," Proceedings of the American Control Conference, pp. 1380-1385, vol. 3. [3] Baker, J . E . , (1985). "Adaptive Selection Methods for Genetic Algorithms," Pro-ceedings of an International conference on Genetic Algorithms and Their Applica-tions, pp. 101-111. [4] Bethke, A . D. , (1980). Genetic algorithms as function optimizers. Doctoral disser-tation (CCS), University of Michigan, Ann Arbor, MI . [5] Canudas, C , K . J . Astrom and K . Braun, (1987). "Adaptive compensation in DC-motor drives," IEEE Journal of robotics and automation, pp. 680-685, vol. RA-3 , No. 6, December. [6] Clarke, D. W. , (1981). "Implementation of self tuning controllers," Self-tuning and adaptive control: Theory and applications, (eds. C . J . Harris and S.A. Billings), pp. 144-165. [7] Clarke, D. W. , (1984). "Self-tuning control of nonminimum-phase systems," Au-tomatica, vol. 20, pp. 501-517. [8] Das, R. and D. E . Goldberg, (1988). "Discrete-Time parameter estimation with 89 Bibliography 90 Genetic algorithms," Proceedings of the 19l annual Pittsburgh conference on mod-elling and simulation. [9] Davis, L . (Ed.) (1987). Genetic Algorithms and Simulated Annealing, London: Pitman. [10] DeJong, K . A . , (1975). An analysis of the behavior of a class of genetic adaptive systems. Doctoral dissertation (CCS), University of Michigan, Ann Arbor. [11] Etter, D. M . , M . J . Hicks, and K . II. Cho, (1982). "Recursive adaptive filter design using an adaptive genetic algorithm." Proceedings of the IEEE International conference on Acoustics, Speech and Signal Processing, pp. 635-638, vol. 2. [12] Goldberg, D. E . , (1981). System identification via genetic algorithm, Unpublished manuscript, University of Michigan, Ann Arbor, MI . [13] Goldberg, D. E . , (1983). Computer-aided gas pipeline operation using genetic al-gorithms and rule learning. Doctoral dissertation (civil engineering), University of Michigan, Ann Arbor. [14] Goldberg, D. E . , (1987). "Simple genetic algorithms and the minimal deceptive problem," Genetic algorithms and simulated annealing, Lawrence, D. (Ed.), Pit-man Publishing, pp. 74-88. [15] Goldberg, D. E . , (1989). Genetic algorithms in Search, optimization and Machine Learning, Addison-Wesley. [16] Grefenstette, J . J . (Ed.), (1985). Proceedings of an International conference on genetic algorithms and their applications, Hillsdale, N J : Lawrence Erlbaum Asso-ciates. Bibliography 9 1 [17] Grefenstette, J . J . (Ed.), (1987). Genetic algorithms and their applications: Pro-ceedings of the second international conference on genetic algorithms, Hillsdale, N J : Lawrence Erlbaurn Associates. [18] Holland, J . H . , (1970). "Outline for a logical theory of adaptive systems," A. W. Burks (Ed.), Essays on cellular automata, pp. 297-319, University of Illinois Press. [19] Holland, J . H . , (1975). Adaptation in Natural and Artificial Systems, University of Michigan Press, Ann Arbor. [20] Johansson, R., (1986). Identification of continuous time dynamic systems, Tech-nical Report, Department of Automatic Control, Lund Institute of Technology, Lund, Sweden. [21] Kristinsson K . and G. A . Dumont (1988). "Genetic algorithms in system identi-fication", Third I E E E international symposium on intelligent control, Arlington, V A , U S A . [22] Smith, T. and K . A . DeJong, (1981). "Genetic Algorithms applied to the calibra-tion of information driven models of US migration patterns," Proceedings of the 12th Annual Pittsburgh conference on Modelling and Simulation, pp. 955-959. [23] Soderstrom, T., L . Ljung, and I. Gustavsson, (1974). A comparative study of re-cursive identification methods, Report 7427, Lund Institute of Technology, Lund, Sweden. [24] Wittenmark, B . and B . J . Evans, (1988). " A n adaptive pole placement controller based on pole-zero parameterization," Preprint from the 8th IFAC/IFORS sympo-sium on identification and system parameter estimation, pp.98-103, Beijing, China. Appendix A Genetic algorithms procedures Program P P C G A P R O C E D U R E Select a population ; B E G I N F O R all the population size DO F O R all the bits DO IF random > 0.5 T H E N bit := 1 E L S E bit := 0 ; E N D ; P R O C E D U R E Schemata ; B E G I N { Count how many 1 there are in each bit position using one counter for each bit position. Then count lost bits by counting number of counters that have value equal to 0 or population size. Then count converged bits by counting number of counters that have value less than converged % or those with value greater than (1 - converged) %. } E N D ; P R O C E D U R E System ; 92 Appendix A. Genetic algorithms procedures B E G I N { Choose either P R B S input (no control) or setpoint changes } IF setpoint changes T H E N Every bitinterval multiply the input yref(t) with -1.0 ; IF P R B S input T H E N Every bitinterval make the P R B S sequence using shift register as given Eykhoff [1974] ; Find e(t), the normally distributed random sequence ; Call the controller design with the best estimate to find the next controller output, u(t) ; E N D ; P R O C E D U R E Convert the strings into the parameters ; B E G I N F O R all the population DO F O R all the substrings DO B E G I N x := decimal value of the binary substring ; resolution := maximum value - minimum value 2length oi substring _ j ' y := x * resolution + minimum value ; E N D ; E N D ; P R O C E D U R E Fitness evaluation ; Appendix A. Genetic algorithms procedures B E G I N F O R all the population DO B E G I N F O R t-window size T O t-0 DO B E G I N IF RLS T H E N B E G I N e:=y-ju; fitness := bias - e2 + fitness ; E N D E L S E B E G I N y := f u • v-=y-y ; fitness := bias - r/2 + fitness ; E N D E N D ; E N D ; E N D ; P R O C E D U R E Quicksort; Use algorithm given in the book D A T A S T R U C T U R E T E C H N I Q U E S by Thomas A . Stan dish page 25-27 or any other sorting algorithm Appendix A. Genetic algorithms procedures PROCEDURE Rank; BEGIN max — 1 a := 2 popsize — 1 ' max — 1 b := 1 - : -(popsize + 1) ; popsize — 1 FOR all the popsize DO fitness := a * rank + b ; END ; PROCEDURE Offspring ; BEGIN FOR all the population DO BEGIN Normalized fitness := fitness / meanfit ; Int fit := integer value of normalized fitness ; Offspring count := Intfit ; Throw a dice to decide if the string gets one offspring for the fractionalpart or not ; END ; IF there are too few or too many in the population THEN choose the one that were most likely to get additional offspring and did not or the one that were most unlikely to get additional offspring and did, depending on whether the population size is too small or too large respectively. ; END ; Appendix A. Genetic algorithms procedures P R O C E D U R E Copies ; B E G I N W H I L E population size is not full DO B E G I N Next string ; W H I L E offspring count < 0 DO B E G I N Copy string ; Offspring count := offspring count - 1 E N D ; E N D ; E N D ; P R O C E D U R E Crossover ; B E G I N F O R half the population size DO B E G I N Find the first string that has not been used Choose another string randomly ; Apply crossover with pc probability ; IF crossover T H E N B E G I N Choose crossover point randomly ; Appendix A. Genetic algorithms procedures 97 Copy first half of first string up to crossover point and second half of second string from crossover point to the end into the first offspring ; Copy first half of second string up to crossover point and second half of first string from crossover point to the end into the second offspring ; E N D E L S E B E G I N Copy first string into first offspring ; Copy second string into second offspring ; E N D ; E N D ; E N D ; P R O C E D U R E Mutation ; B E G I N F O R all the population size DO F O R all the bits in each string DO Mutate each bit with pm probability ; E N D ; B E G I N Get all parameters ; Initialize random generators ; Appendix A. Genetic algorithms procedures Initialize system input and output ; Print initial values ; Select a population ; Start with the initial estimate as 1, 0, ..., 0 ; F O R kids := 1 T O number of kids DO B E G I N Schemata ; IF (kids-1) / trial = integer T H E N system ; Convert the strings into the parameters ; Fitness evaluation ; Calculate the average fitness ; Offspring ; Count how many receive 0 offsprings ; IF (receive 0 offsprings) > (fitOpct * popsize) T H E N B E G I N Quicksort ; Rank ; Offspring ; E N D ; Copies ; Crossover ; Mutation ; Find the best string ; IF it is not in the new population T H E N Replace a string randomly chosen with the best one Appendix A. Genetic algorithms procedures 99 Print report ; Make the new generation the current one ; E N D ; E N D ;