Home

(non) linear Dynamical Systems using Genetic Programming

1. 22 7 GP amp controller design 7 1 Introduction Application like pick and place units robot arms printers etc are created to perform spe cific tasks The movable part is driven by a mo tor and the position is measured The field of motion control practises the realtime creation of the motor power supply according to measure ments considering user requirements and con sidering knowledge about the system 5 The suggested applications above all have slightly different requirements A pick and place unit should stand still at the desired place as fast as possible A robot arm must fol low a trajectory as good as possible at every time And a printer should follow a trajectory with a steady velocity The power of this im plementation of the GP algorithm is that ar bitrary requirements can be specified in time and frequency domain Because the GP pro gram only optimizes dynamical models with a fixed number of inputs only the fitness function has to be changed to make the step from system identification to controller design 7 2 Fitness function Creating a fitness function for controller design is a difficult task Each model has to be judged properly else GP will create wrong controllers with an apparently good fitness GP optimizes what the user specifies no more no less So the user has to consider carefully what he really de sires of a controller without using assumptions valid for PID like controllers For
2. disp simulink models disp system system model md1 1 disp gp model gp_model md1 disp 7 dispC controller Ep sys zpk gp_z gp_p gp_k gp_sys minreal gp_sys gp_z gp_p gp_k zpkdata gp_sys gp_z gp_z gp_p gp_p disp poles real_gp_p index_gp_p disp gp p index gp p disp zeros real_gp_z index_gp_z disp gp_z index_gp_z disp gain num2str gp_k sort real gp_p sort real gp_z 42 B FITNESS FUNCTION M FILE disp if f BadSimPenalty disp performance if settling time disp settling time gt num2str settling time s error lt num2str accuracy 100 gt else disp settling time num2str settling time s error lt num2str accuracy 100 end disp bandwidth num2str bandwidth Hz displ displ constraints disp max sensitivity num2str maxS dB lt num2str SensLimit dB 1 disp max controller output num2str maxU N lt num2str ActuatorLimit N else disp The gp model is unstable or simulation problems have occurred return end determining transfer function of plant and controller Tplant Terror EstimTF system model fs sens n sens n sens parts sens overlap S in Amount s in Controller ut s out PlantOut MST Tcontroller Terror EstimTF gp model fs sens n sens n sens pa
3. nodeproperties cpp nodeproperties h settings cpp settings h fitness h eogaoperator h eogencontinue2 h eomutationoperators h eoparsetreedepthinit2 h eosgatransform2 h 55 D GP PROGRAM D 2 C code description D 2 2 Description gp4so cpp This is the main application During exploration of the EO library and the creation of the first application a sample application was used which is dedicated to symbolic regression This was quite a mess which has been partly reduced in the current application It is still no optimal architecture or fully user friendly programmed but it has improved a lot It just takes quite a time the EO library begins to live Besides the EO library is still in development and not finished engine cpp engine h This is the class which makes the Matlab engine available It reads the Matlab path from the windows registry and uses libeng d11 and libmx d11 from Matlab to connect to some necessary functions to operate the Matlab engine At the moment it is only compatible with Matlab 6 1 and Matlab 6 5 inireader cpp inireader h This class enables reading of the data types double int bool and string from an ini file The ini file should consist of sections with and keys belonging to a section The value of a key is separated with With these parameters the program can be tuned without recompiling the program logdata cpp logdata h This class collects data about settin
4. static void mdlTerminate SimStruct S ifdef MATLAB MEX FILE Is this file being compiled as a MEX file include simulink c MEX file interface mechanism Helse include cg_sfun h Code generation registration function fiendif 35 A S FUNCTIONS 4 2 protecteddiv A 2 protecteddiv protecteddiv dll is an S function for Simulink which is called a protected division It works the same as a normal division except when the denominator approaches zero The denominator saturates if its absolute value is smaller than 10 9 The sign is unaffected If the denominator is zero it is set to 1076 The source code of this S function is shown below define S_FUNCTION_NAME protecteddiv define S_FUNCTION_LEVEL 2 include simstruc h static void mdlInitializeSizes SimStruct S ssSetNumSFcnParams S 0 if ssGetNumSFcnParams S ssGetSFcnParamsCount S return Parameter mismatch will be reported by Simulink if ssSetNumInputPorts S 2 return ssSetInputPortWidth S 0 1 ssSetInputPortWidth S 1 1 ssSetInputPortDirectFeedThrough S 0 1 ssSetInputPortDirectFeedThrough S 1 1 if ssSetNumDutputPorts S 1 return ssSet utputPortWidth S 0 1 ssSetNumSampleTimes S 1 Take care when specifying exception free code see sfuntmpl_doc c ssSetOptions S SS_OPTION_EXCEPTION_FREE_CODE SS_OPTION_USE_TLC_WITH_ACCELERATOR static void mdlInitializeSampleTimes SimStruct
5. B FITNESS FUNCTION M FILE end end grid title plant gain xlabelC f Hz ylabel dB plant phase subplot 234 semilogx freq 2 end 180 pi unwrap angle Tplant 2 end k axis tight grid title C plant phase xlabel f Hz ylabel degree controller gain subplot 232 semilogx freq 2 end 20 logi0 abs Tcontroller 2 end k axis tight grid title controller gain xlabel f Hz ylabel C dB controller phase subplot 235 semilogx freq 2 end 180 pi unwrap angle Tcontroller 2 end k axis tight grid title controller phase xlabel f Hz ylabel degree closed loop gain subplot 233 semilogx freq 2 end 20 log10 abs 1 S 2 end k axis tight grid title closed loop gain xlabel f Hz ylabel dB Y closed loop phase subplot 236 semilogx freq 2 end 180 pi unwrap angle 1 S 2 end k axis tight grid title closed loop phase xlabel f Hz ylabel degree show the closed loop system open_system system_model plant open_system gp_model GPmodel 45 C GA IMPLEMENTATION C GA implementation Like GP GA also uses a population The major difference is that the structure of an individual is a fixed amount of numbers on a row The working principle is written next The words between brackets are parameters of GA which can be specified
6. D GPPROGRAM D 1 User manual D 1 2 Description of files gp4so exe This is the executable file which starts a GP run Before execution the settings in gp ini should be checked and the fitness function must be set up correctly After that the only thing left is running gp4so exe and waiting for a good solution which minimizes the fitness function This can take some hours to a few days During a GP run you can view intermediate results The easiest way to do that is to press the Pause button on the keyboard and then go to the Matlab command window which has been opened by the GP program Some ways to examine intermediate results are listed next To examine the current model where gp4so is working on type fitness gp_model To only view the current model type open system gp mode1 GPmode1 To examine the best model so far type fitness temp_model md1s To examine a model which was the best model in an earlier state type fitness temp_modelxxx mdls xxx is the number of that model When evaluating a temp_modelxxx mdl also a plot of the fitness progression is drawn This is only drawn until the best model so far so not until the present time You can determine how long ago that best model was found by looking at the amount of fitness evaluations at the end of the current generation This is useful when considering the convergence of the GP run Every generation some statistics are shown like the best fit
7. 1 e m 2r 3F 4 1 1 1 0 0 5 1 1 5 2 t s b Error between output and trajectory Figure 7 5 The output compared with the tra jectory Figure 7 6 shows that stability and robustness are preserved and in Figure 7 7 the open loop system is plotted with a bandwidth of 158 Hz 26 7 8 Results Magnitude dB Phase degree 10 Frequency Hz Figure 7 7 Open loop transfer function 7 3 2 Controller design for a fourth or der system System The second system is the Pato setup the smaller type located in the Lab of Control Sys tems Technology at the University of Eindhoven It consists of two rotating masses with a flexi ble beam in between see Figure 7 8 At one mass an electro motor is connected and the an gular position of both masses are measured by encoders with a resolution of 2000 increments per revolution This system behaves mainly as a fourth or der system A controller can be designed for two possible cases measuring at the motor side and measuring at the load side The transfer function of V Zmotor is different than V 2icaa because of the flexible beam in between Both transfer functions have been determined by mea suring the sensitivity in closed loop A weak 7 GP amp CONTROLLER DESIGN 7 8 Results Sy electro motor X motor Xload Figure 7 8 Schematic view of the Pato setup PD controller is used to stabilize the closed loop Furt
8. 10 10 10 10 10 10 Frequency Hz 2 1 10 Figure 7 15 Controller transfer function GP has created a PID controller with a low pass filter and a notch which compensates the resonance peak The characteristics of the con troller are listed below In Figure 7 16 the response of the output is shown together with the trajectory During movement of the trajectory the output lags lit tle bit behind the trajectory The output does not have any overshoot Figure 7 17 shows that stability and robustness are preserved The open loop system is drawn in Figure 7 18 7 3 Results 0 05 15 2 1 t s Figure 7 16 Time response Figure 7 17 Nyquist plot o Magnitude dB Phase degree I 1 1 1 Q GON Q Oo 8 eoo eo o C Ari 10 10 10 Frequency Hz Figure 7 18 Open loop transfer function 7_GP amp CONTROLLER DESIGN 7 4 Summary Simulations GP is able to generate controllers which are optimized for the specified performance stabil ity and robustness Three simple SISO systems derived from experimental setups where used to create controllers for Every GP run is differ ent so a variety of controllers have been created Slight adjustments can be made to the fitness function settings so the user is able to balance the importance of various parts The mass damper system used a relative slow trajectory After a while GP gene
9. 100000 Ymaximum jerk of trajectory in m s3 amax 1000 Ymaximum acceleration of trajectory in m s2 vmax 50 ymaximum velocity of trajectory in m s Y frequency analysis parameters sens n 2 11 Amount of data points for determining the sensitivity sens parts 2 Division of total signal in smaller parts sens overlap 0 5 Overlap of smaller parts fs 1000 sampling frequency fitness penalty parameters BadSimPenalty 1e10 Penalty for incomplete simulations ConstraintPenalty 1e5 Penalty if a constraint has been exceeded requirement parameters accuracy 001 settling time is determined using a max error of accuracy 100 SensLimit 6 constraint maximum sensitivity in dB ActuatorLimit 1 constraint maximum actuator amplitude in N MST 10 constraint maximum simulation time in s disturb 0 amplitude of disturbance signal noise noise 0 Yamplitude of sensor noise signal noise bandwidth factor 0 bandwidth factor 0 disabled else f f exp bandwidth factor bandwidth 3 pre checks Seek E I Id I AARAA AAA Ahhh hhh hhhhhkhhhhhhhhhhhhhhhhhhhhhhhhh THE REST OF THE M FILE DOES NOT HAVE TO BE CHANGED FOR CONTROLLER DESIGN IIITIII ZII NAR AAA N NAR RIA Y TAA ATA A AA RT Y AL Mr NN A Y RAN AL ND NANA AANG Ahi PRE CHECKS Y this part is only used by gp4so It is used to get the system_model without simulating it if nargout f BadSimPenalty return end check if anoth
10. 6 3 6 4 6 5 6 6 6 7 6 8 6 9 7 1 7 2 7 3 7 4 7 5 7 6 7 7 7 8 7 9 7 10 7 11 7 12 7 13 7 14 7 15 7 16 7 17 7 18 An example of a tree structure y ue u2 2 0 0200 a An overview of the GP algorithm o e eee eee ee An illustration of a crossover operation 2222s An illustration of a mutation operation 2 6 oe e a a ee ee An illustration of a reproduction operation 2 2 eee a Overview of GP research eee A global view of the architecture of the implemented GP algorithm Properties of a second order transfer function eee e Plant to be identified 222222 Best GP model si taart Arent aen eee E ep ae enen Se ed e Y me Comparison of the plant and the best model using the same input data as used in the fitness functlON sven ern eee Be eu ah GU bie pene mU X doe Roe Validation of the best model using different input data Plant to be identified llle Best GP model m nu nsnm are A edd le eh aod doe eus Comparison of the plant and the best model using the same input data as in the fitness FUNCION Poder eee Vg re di ated th da as Ds Validation of the best model using different input data Transfer function gain phase and coherence function o ooo Closed loop systemi 22e els Interpretation of a time domain simulation eee Transfer function o
11. GP only mutates whole nodes so the chance it returns the same block with slightly different parameters is very small However having more parameters per node is unsurmountable for controller design in the s domain using Simulink and a tree representa tion of an individual The advantage of using the s domain is that such a model is directly in terpretable for control engineers Disadvantages are that discretisation is not accounted for and that a transfer function can not be made with el ementary blocks when using a tree structure A feedback coupling is needed to create a transfer function using elementary blocks In the z domain elementary blocks can be used which is preferable for GP A network structure could be used to generate models No conversion from continuous to discrete is needed 5 MULTIDIMENSIONAL OPTIMIZATION OF A NONLINEAR FUNCTION when performing a simulation Another issue is that now a days most controllers are imple mented in a discrete form for example on a DSP So using the z domain already incorporates dis cretisation In short the z domain provides some nice ad vantages but the s domain has been chosen be cause of the interpretability Due to this choice and the choice of using a tree structure repre sentation for an individual transfer functions have to be included to the function set and a parameter optimization algorithm is indispens able Which additional parameter optimization routine
12. The meaning of these inputs is determined by the fitness function Section Matlab Simulink Settings Name of the fitness function Work directory location of fitness function dll s and other needed files Space between blocks in simulink 50 D GPPROGRAM D 1 User manual Section GA optimization set to true to return the best individual of all generation set to false to return the best individuals of the last generation fitness m file The fitness function is an m file which determines how good a model is The input is a Simulink model with a fixed amount of inputs The output is a number which indicates how good a model is the fitness value A small number is good and a large number is bad 0 is the minimum and BadSimPenalty 1010 is the maximum This fitness value is used in the GP and GA algorithms for minimization So it is crucial to define a good fitness function It should always return a valid number and your demands have to cover all possibilities in the search space fitness m plant 2nd order mdl fitness pato m pato motor mdl pato load mdl These are working example files to get a better insight how to set up such files Before using them check the gp ini file if the correct fitness function has been selected gp_model The GP program creates a Simulink model every time it wants to determine a fitness of a particular individual This model is created du
13. a transfer function is called a Bode plot Besides the stability and robustness gain and phase margin also the bandwidth can be determined The bandwidth is defined as the frequency where the open loop crosses the 0 dB line The sensitivity function is shown in Equa tion 7 3 It represents the suppression of the disturbance signal ug s effected on the plant input signal u s Further the reciprocal of the absolute value is the distance to the point 1 in a Nyquist plot 7 2 u s _ 1 Wc worn 1 23 7 2 Fitness function Performance Direct performance requirements are time domain based like settling time overshoot static error etc Since the GP program can handle time and frequency domain criteria the performance is specified in the time domain The fitness function is based on a pick and place like task CO iin DI For the trajectory u a 3th or der set point is chosen see Figure 7 2 T is the time that the trajectory becomes constant again T is the total simulation time Ts is the settling time which is specified as the time after which the difference between the plant output and the trajectory stays within an error band until the end of the simulation time T emaz is the maximum absolute error between Te and T position m time s Figure 7 2 Interpretation of a time domain sim ulation if T gt T then f gt f if T 20 9T then f fm 7 4 The share of the performance crite
14. be returned and it should be a good indication how good a model is Furthermore the fitness function should always return a value else undefined behavior can occur 3 Run the gp4so exe and wait Better think in days instead of minutes Every fitness evaluation a simulink model has to be created and simulated so that just takes a lot of time Intermediate results can be checked as described before 4 If a termination criterion has been reached or you have pressed a key you will be asked if you Nr really want to abort the GP run If you type y you can save the best model including additional data 5 Finally the saved result can be evaluated and validated 54 D GP PROGRAM D 2 C code description D 1 5 Troubleshooting If the Matlab engine does not start check if you have Matlab version 6 1 or 6 5 The application does not work with other versions at the moment When you get a License Manager Error from Matlab you do not have a correct license dat file in the matlabdir bin win32 directory This problem occurs for example when this file is accessed via a network for multiple use Contact the system administrator about this to copy the file to your local computer Oui OVAAL COMMPuLel If at the startup some warnings occur about missing items in gp ini you should check if you spelled those names correctly If strange fitness numbers occur or the GP program throws an error or just quits immediately probably the fi
15. bution of probability For each offspring this standard deviation of the step size is increased or decreased by a factor 1 5 Parents pass on this standard deviation of the step size to its offspring By selecting the best individuals the standard deviation of the step size will adapt to the shape of the fitness landscape This method has shown to be very robust for a variety of problems 5 4 Choice of method Some consideration for the numerical optimiza tion of the parameters of an unknown model are as follows The amount of evaluations has to be re duced to a minimum because each fitness evaluation contains one to several simula tions in Simulink which costs time Model parameters are randomly initialized so compared to the optimum values a huge deviation is likely to be present ini tially Because of this a global optimiza tion method is preferred instead of a local optimization method The landscape of the fitness function is un known and probably highly nonlinear Also the sensitivity of all parameters can vary a lot The algorithm should be somehow self adjustable to the right order of magnitude to change the parameters The solution does not have to be the ex act global minimum after one optimization GP keeps on optimizing models by random selection of individuals from the population so it is possible that a model is optimized more than once 6 GP amp SYSTEM IDENTIFICATION The fitness landscapes
16. ceil n 10 n 2 semilogy fitness plot data 2 n 1 1 fitness plot data 2 n 1 1 1 6 4x33 SAAI 7 fitness plot data 2 n i 2 miny ones n 2 1 Color 7 N M hold on semilogy x y k LineWidth 3 axis min x max x miny maxy grid xlabel fitness evaluations ylabel fitness title fitness progression end figure SECOND FIGURE Kopen loop gain subplot 221 ol_dB 20 log10 abs 1 S 2 end 1 semilogx bandwidth bandwidth min ol_dB max ol_dB color 7 7 7 hold on semilogx freq 2 end ol_dB k axis tight grid title open loop gain xlabelC f Hz ylabel dB open loop phase subplot 223 semilogx freq 2 end 180 pi unwrap angle 1 S 2 end 1 k axis tight grid title open loop phase xlabel f Hz ylabel degree nyquist subplot 222 plot real 1 S 2 end 1 imag 1 S 2 end 1 k 1 0 k hold on linspace 0 2 pi 100 plot 1 sin w 107 SensLimit 20 cos w 107 SensLimit 20 color 7 7 7 axis 5 5 5 5 grid title nyquist xlabel Re ylabel Im sensitivity subplot 224 semilogx freq 2 end 20 logi0 abs S 2 end k axis tight grid title sensitivity xlabel f Hz ylabel dB figure THIRD FIGURE plant gain subplot 231 semilogx freq 2 end 20 logi0 abs Tplant 2 end k axis tight 44
17. disturbance rejection at the lower frequencies below the bandwidth al ways results in a worse disturbance rejec tion in the higher frequencies above the bandwidth Better would be to insert low frequent noise instead of white noise to stimulate the disturbance rejection at the lower frequencies Adding noise to the plant input will also pre vent that GP creates a pure feedforward con troller The added disturbance signal can be seen as a model uncertainty However the in tensity of the noise may not be too high causing a bigger error than the error band The settling time can not be optimized in that case Error handling Finally erroneous controllers controllers which cause the closed loop to be unstable and will often cancel the simulation get the maximum penalty of 1 10 Controllers with stiff differen tial equations often consume a lot of simulation time After a user specified time the simulation will be aborted and the fitness function will re turn the maximum penalty Some simulations throw an error before this time limit because of unsolvable differential equations This is also in tercepted and the maximum penalty is returned The total fitness function consists of the next elements if Ts gt To then f f if T gt 0 9T then f f Emar A 7 GP amp CONTROLLER DESIGN 7 3 Results E A A eeoa if the maximum sensitivity gt 6 dB then f 10 if the actuator limit has be
18. example the gain increases with decreasing frequency in the low frequency region because of an I action For a randomly generated controller containing lots of transfer function blocks this may definitely not be assumed A dip somewhere in the low frequency gain could lead to bad tracking per formance or bad disturbance rejection Dynamical systems Before presenting all elements of the fitness function a brief introduction is given about the dynamics used for controller design The scheme of a closed loop system is shown in Figure 7 1 7 GP amp CONTROLLER DESIGN Figure 7 1 Closed loop system The closed loop transfer function is given in Equation 7 1 for a linear SISO system It shows the relation between the user defined tra jectory and the actual output of the plant The controller has been divided into a feedback and a feedforward part Often both parts are de signed separately The GP program will be of fered the error the velocity and the acceleration signal to create the controller COP u s T G s P s The open loop transfer function is given in Equation 7 2 It represents the relation between the error e and the plant output y T s 7 1 y s d C s P s The open loop transfer function plotted in the complex plane is called a Nyquist plot It is use ful to determine stability passing side of point 1 and robustness distance to point 1 A plot of the absolute value and the angle of
19. genetic operations Results crossover 7096 branch mutation 10 branch addition 596 point mutation 15 GA optimization 20 10 per generation using a tournament selection of 2 The GP run was aborted manually because the fitness did not improve anymore for awhile Properties of the best individual are The best model is shown in Figure 6 2 The structure is exactly correct Only the numerical values vary a little bit 768 476 A 0 265861 s2 3 99296s 100 067 n2 ni InputO Figure 6 2 Best GP model output In Figure 6 3 the plant and the best model are compared for all fitness cases The error left 6 GP amp SYSTEM IDENTIFICATION mainly originates from the inserted noise in the emory y plant input u output y grey y black 1 0 5 o 3 0 1 2 3 1 05 0 0 5 1 0 1 2 3 1 0 5 0 0 5 1 1 2 3 tls t s 0 1 2 tis o Figure 6 3 Comparison of the plant and the best model using the same input the fitness function Validation data as used in A good result is obtained when applying differ ent data sets see Figure 6 4 It can be con cluded that the GP model is a good approxima tion of the plant input u output y grey y black eror y Y 0 0 1 0 1 t e N o da o Nn o o m o o 1 o o o D 1 4 05 2 0 o 0 5 2 4 4 0 1 2 3 o 1 2 3 y 5 0 o 2 5 o 1 2 3 0 1 2 3
20. have been normal ized to a static gain of 1 with an extra static gain parameter multiplied with the transfer function This parameter is initialized with a loguniform 14 4 5 Discussion probability distribution so the chance of ampli fying or reducing the gain is equal This could also be solved by GP as GP can add a gain block in front of a transfer function block The reason that a transfer function is given a static gain parameter is to give the trans fer function more freedom so it can have any shape for that type of transfer function A dis advantage is that it adds a parameter extra in one node so the pressure is put a little bit more at the numerical optimization of the parameters see Chapter 5 instead of GP In this case it does not make a big difference because there are already more than one parameters in a transfer function node s domain vs z domain The main purpose of this research was the ex ploration of the suitability of GP for controller design considering easy plants Results from conventional methods are known so the per formance of GP can be compared with these benchmark problems For interpretability the s domain has been chosen to create models In the future it could be of interest to switch to the z domain This is explained later on Having more parameters per node is not pro motable when using only the GP algorithm for parametric and structural optimization because GP lacks efficiency
21. not have to edit the fitness function by hand anymore The fitness function can be divided into the following parts 1 function description This part describes the usage of the fitness function by the GP program and by the user 2 adjustable parameters These parameters can be adjusted to specify your demands The rest of the fitness function does not need to be changed for controller design 3 pre checks These checks are inserted for the user s convenience It makes the fitness function more flexible For example if you want to evaluated the best model so far during a GP rum you can simply type fitness temp model mdls which evaluates the best model without having to look what the most recent model is 4 setup simulation model This part inserts the GP model in another Simulink model which can then be simulated 5 fitness determination Here simulation are run and the fitness value is determined 6 validity check of the fitness value After the fitness function has been determined a check is made if the fitness value is valid The fitness function should always return a vaild fitness value 7 visualisation of the results If a directory name was specified the results of the fitness evaluation are plotted Because the GP program never specifies the directory name it uses the GP model in the memory and files in the same directory as the fitness function the results are not shown if the GP program evaluates the fitness fun
22. other So it is pos sible that none or more than one genetic operations are applied to the same model of a single generation 4 4 Conversion of a node to a Simulink block GP builds a tree from nodes Each type of node has different properties like the amount of in puts and the amount of parameters it contains If GP creates a tree using these nodes this has to be converted to a Simulink model to be able to evaluate it Each node is converted separately and connected as indicated by the tree In Ap pendix D 1 3 the Simulink blocks are presented belonging to all available node types which are Add Subtract Multiply Divide Integrate Dif ferentiate TF 1st order TF 2nd order Gain Constant Time Delay Input0 Inputl In put7 This can easily be expanded All available Simulink blocks in the Simulink Library could be used if they would be needed Also Simulink makes it possible to define your own block by creating a so called S function An example of a second order transfer func tion is shown in Figure 4 2 Explanation of the chosen structure can be found in Ap pendix D 1 3 For each node type the Simulink block and the formula is given Parameters of a block are initialized with different probability distributions to create more suitable blocks n is the user specified maximum initialization value 4 IMPLEMENTATION OF GP TF 2nd order C2 ci cl s242 c0 cls ci ci TF2 2 y s ora Ys co uniform 0 1 c
23. tfl tf2 gain GP builds trees using these two sets A termi nal does not have any inputs so it finalizes a branch Functions have one or more inputs so they cause a tree to grow With these nodes GP will be able to identify linear systems tfl and tf2 are respectively a first order and a second order transfer function Both blocks are needed to enable creation of both real complex poles and zeros so every possible proper system can be build Result of a GP run A GP run can be terminated several ways a maximum number of generations a certain fitness level or manually Then the GP program will create a Simulink file of the best GP model This GP model can be validated by applying other inputs than the one used for determining the fitness Also comparing the transfer function of the best model and the plant is very useful to assess the result 6 GP SYSTEM IDENTIFICATION 6 2 Fitness function As explained in Chapter 2 a fitness has to be assigned to each GP individual to indicate its performance Only time domain evaluations are used during system identification while exciting the GP model with a sine a step or an impulse Identification will be based on a set of m inputs The output of the GP model should resemble The output of the GP model should resemble the plant as much as possible for all m input cases Therefor a fitness function proportional to the squared output difference between the GP model and the plant
24. used genetic op eration with a common usage of 30 It has a less destructive effect than crossover because on average it makes smaller steps through the search space It also keeps the diversity of node types in a population high because it randomly inserts new nodes from the terminal and func tion sets 3 GP RESEARCH Reproduction Reproduction is simply creating a new individ ual by making an exact copy of the selected par ent Figure 2 5 y uet u 2 y ue u 2 parent child Figure 2 5 An illustration of a reproduction op eration Reproduction is useful to preserve the best individual also called Elitism Research has pointed out that the reproduction operation does not add a lot of performance gain to GP because it does not change the behavior of an individual a low usage of 2 is common Other genetic operations Other genetic operations which have been used during this research are expansion mutation numeric mutation and branch addition Ex pansion mutation replaces a terminal node by a randomly created branch Numeric muta tion changes all parameters of an individual ran domly Branch addition creates a new tree and combines that with the old tree by joining them with a 2 input node The first two have been re moved at a later stage because they appeared to be obsolete Branch addition is still present be cause it stimulates linear combinations of parts which is often beneficial in the
25. used in a simulink model It terminates the simulation if a maximum running time is reached The parameter of the simulink block is specified in seconds This block should be present in every system model the fitness function uses for simulations Without this complex models will take a lot of time to finish The reason is often that the model contains stiff differential equations which are hard to solve Such models have to be thrown away and are given a bad fitness protecteddiv dll An S function which is used instead of a pure division block When dividing by zero an error occurs This is prevented by using this block The denominator saturates if the absolute value is smaller than 1 1078 so it will not become 0 which prevents division by zero gp ini An ini file to set parameters which are used by the GP program gp4so exe The file is divided into the sections GP Settings Genetic Operations Nodes Node Settings Matlab Simulink Settings and GA optimization To understand and get a feeling how to use these settings you have to explore the working principle of GP first The sections and keys in the ini file will be described next Section GP Settings i fault value 0 TournamentSelect Amount of candidates used for Tournament Selection when selecting parents erent cos when potting ln back in tho population er 1 TerminationFitness Fitness value when the GP run stops Generation when the GP run stops 50000
26. 0 TE 49 D GP PROGRAM D 1 User manual Section Genetic Operations fault value Description Exchanges a subtree from two parents BranchMutation Creates a new tree at a randomly chosen point BranchAddition Adds a tree by inserting a block with two inputs at the top and let a new branch grow at the free input PointMutation Exchanges a node by a new node GA Optimization Optimization of the parameters of an individual The key values indicate the chance that it is used The chances are independent so it is possible that more genetic operations are applied to the same individual or even none The GA optimization algorithm is also seen as a genetic operation it only optimizes the values of an individual Section Nodes 10 Constant Time delay node Input Inputi Input The key values indicate the relative chance it is used relative with respect to other node types Setting it to zero it is not used For a more detailed description of these node see Section 4 4 Section Node Settings Maximum initialization value of parameters in a new created node 0 Set to true to enforce negative poles in the transfer functions and differentiator Set to true to enforce negative zeros in the transfer functions Be sure you added AmountOfInputs input nodes to the section Nodes starting with InputO Inputi Input7 At the moment a maximum of 8 inputs is possible
27. 1 uniform 0 n ca loguniform i n Figure 4 2 Properties of a second order transfer function 4 5 Discussion Conditioning and normalization The reason of designing blocks in a special form is to create useful models more frequently If this was not done initially too much initial ized models would result in stiff or even unstable models It is obvious that this has a negative ef fect on the performance of GP During the first tests of system identification only a single uni form probability distribution between n and n was used to initialize the parameters Further more the parameters where not positioned in the transfer function to take account of normal ization and the difference in order of magnitude between all coefficients of transfer function For example the parameters of a second order transfer function like p nr uem are often set to obtain an under critical damped transfer func tion In that case the order of magnitude of cz is about the order of magnitude of cj squared In Appendix D 1 3 the special form of each block is shown The tests after these adjust ments have shown the expected improvements Actually designing the blocks in such a way is equal to inserting knowledge about the problem and lighten the task of GP a little bit The parameters have been put in such a way that they are conditioned well Most of them have the order n and a few the order 1 Further more all transfer functions
28. 10 40 5 20 0 o 5 20 10 40 o 1 2 3 0 1 2 tis tis o N a ts Figure 6 4 Validation of the best model using different input data 20 6 3 Results 6 3 2 System identification of a fourth order system System Two damped masses with a spring and a damper in between will be identified next see Fig ure 6 52 A Simulink model has been made which describes the relation between F and zi see Figure 6 5b In Figure 6 7 the 3 input cases are shown m 0 1 kg k 100 N m d 0 001 Ns m d 0 3 Ns m 152 38 100 0184 060153 20 0952 59 9s 4th order system Random Number var 0 0001 b Simulink model F zi Figure 6 5 Plant to be identified Settings data points 400 ae 0 5 Hz sine impulse at 0 3 s erate ein integrate Ist order TF 2nd order TF genetic operations crossover 70 branch mutation 10 branch addition 5 Results point mutation 15 GA optimization 20 10 per generation using a tournament selection of 2 6 1 hours 0 045 319 impulse high frequency component has not been caught 6 GP amp SYSTEM IDENTIFICATION 6 3 Results The GP run was aborted manually The best GP model is shown in Figure 6 6 In Figure 6 7 the plant and the best model are compared for all fitness cases The GP model performs good for these 3 fitness cases but the structure is not the same It main
29. 2 2 cO c1s cl ci 2 Us omm le Co uniform 0 1 c uniform 0 n c loguniform n TF2 c3 c1 c1 c2 5 03 c1 c1 8242 c0 cts ci ci TF2 2 s co41l y s esc Tons q Us Co uniform 0 1 9 uniform 0 n 53 D GP PROGRAM D 1 User manual c3 loguniform i nl 04 c1 c1 03 03 2 04 01 c1 03 c28 04 01 01 s2 2 c0 c1s c1 c1 2 2 2 TF2 cy 8 2e2038C3 EM u y s Ca E 52 2c9c15 2 s Co C2 uniform 0 1 C1 3 uniform 0 n ca loguniform n The second order transfer functions are scaled differently with a similarity to the normalized form s 2Ewns w2 The parameter at the place of is initialized between 0 and 1 to produce under critical damped systems which are mostly used in the field of system identification and controller design Constant Constant y t c c loguniform n Time Delay Time Delay y t u t eR c c uniform 0 10 n Assuming n lies near the sampling frequency then the maximum initial time delay of this block is 10 samples Inputx Inputx y t Us t D 1 4 Performing a GP run Steps for a successful GP run are 1 Check all settings of gp ini and adjust them where necessary 2 Adjust the fitness function according to your desires This part is the most important and most difficult to do so spend enough time for this You need to specify what you really want which sounds much easier than it actually is One value has to
30. Node types corr eux keers OR C9 Ben Ton e dd RAE WES ue 52 D 1 4 Performing a GP rUdM e 54 Dildo Troubleshooting weent a AD A RE E ub A A 55 D 2 C code description 2 1n AE oe a AA Wm mee old dede 55 DIT C pac ae Roy ete RETE tee nd 55 D 2 2 Deseripulon at arte A NE oR AS RE ERR IEN EE A 56 D 2 3 Known bugs possible improvements experience eee 57 Abstract Ever since digital controllers became both available and cheap a complex controller can be imple mented with only a little bit of more effort However the design of such a controller is more difficult Often a mere PID controller is used because it is easy well understood and the performance is known In this thesis a controller synthesis method has been developed which can be applied to a wide range of plants It is based on Genetic Programming GP which creates controllers with variable structures Nonlinear dynamics can be included in the controller and or the plant Requirements for the controller are specified in time and or frequency domain The requirements are translated in a fitness function which is problem dependant This fitness function returns a number indicating how good a certain controller meets the requirements This fitness function will be minimized by the evolution process of GP and thereby finding the optimal solution for the problem GP is investigated in detail with respect to how it works what kind of problems can be ta
31. S ssSetSampleTime S 0 INHERITED SAMPLE TIME ssSetOffsetTime S 0 0 0 static void mdl0utputs SimStruct S int T tid InputRealPtrsType uPtrs ssGetInputPortRealSignalPtrs S 0 real T y ssGetOutputPortRealSignal S 0 if uPtrs 1 gt 0 amp amp uPtrs 1 lt 1le 6 y uPtrs 0 1e6 else if uPtrs 1 gt 1e 6 amp amp GruPtrs 1 lt 0 y uPtrs 0 1e6 else y uPtrs 0 uPtrs i static void mdlTerminate SimStruct S ifdef MATLAB MEX FILE Is this file being compiled as a MEX file include simulink c MEX file interface mechanism else include cg_sfun h Code generation registration function endif 36 B FITNESS FUNCTION M FILE B Fitness function m file The fitness function returns a value which indicates how good a GP model satisfies user specified requirements It is built into one file so all elements of the fitness function are hold together This ensures that the fitness function is the same as the fitness function used during the GP run Enough info lines are added to explain what is done in the fitness function A suggestion for further research is to subdivide this fitness function because it has grown quite large which is not very user friendly The only part to be edited by the user considering controller design is part 2 Another suggestion is to make a GUI to set all parameters of the fitness function and set the settings of GP Then you do
32. Structural and Parametric Optimisation of non linear Dynamical Systems using Genetic Programming A tool for automated controller design Master s thesis by D J H Bruijnen DCT 2003 75 CTB 555 03 1379 Engineering Thesis Committee prof dr ir M Steinbuch supervisor ir LA C Soute coach at Philips CFT dr ir M G J van de Molengraft coach at TU e dr ir G Z Angelis Philips CFT dr ir D A van Beek TU e Eindhoven University of Technology TU e Department of Mechanical Engineering Control Systems Technology Group Philips CFT Mechatronic Equipment Motion Control Group Eindhoven August 2003 Preface This research was a project proposed by and carried out at the Motion Control Group of Philips CFT It is a first exploration whether Genetic Programming is suitable for system identification and controller design Preceding research has been carried out by Iris Soute who has used Genetic Programming for finding Lyapunov functions to prove stability My coaches have been Iris Soute Philips CFT and Ren van de Molengraft TU e Many thanks to them for their good advise and support My project can be divided into three main elements 1 Genetic Programming Algorithms 2 OO Programming 3 system identification controller design I would like to thank the next people from Philips CFT for their support to the project Dani l Vangheluwe 1 Jos Onokiewicz 2 Georgo Angelis 1 3 Felix Peeters 3
33. als are selected to produce offspring with accompanying step sizes and parameter values This will go on until the amount of generations MaxGen has been reached KeepBest is set to true if you want GA to return the best individual of a GA optimization when the maximum amount of generations is reached When set to false it only returns the best of the last generation An advantage of this is that GP will not get stuck in a local minimum In the GP algorithm the current best individual is quite often selected If GA does not change its values because it is in a local minimum and can not get out the amount of identical individuals will grow in the GP population The diversity decreases and it will be impossible to optimize So keeping the best of all generations could lead to stagnation of the GP progress resulting in a GP population of all the same individuals with identical parameters A more detailed description of the algorithm can be found in 4 The GA is implemented in eoGaOperator h as a genetic operation of GP See the GP CDROM 2 to view the source code 46 D GPPROGRAM D 1 User manual D GP program D 1 User manual To become familiar with the software for controller design with GP gp4cd first all supplied files are listed After that some more detail is given about these files Finally an instruction list is formulated to perform a GP run successfully considering the user s desires D 1 1 Gp4cd files Here
34. and Emiel Botermans 3 Furthermore 1 would like to thank some other students doing also their traineeships graduation project here for having interesting discussions about their projects my project and other stuff My neighbor Rick Scholte who knocks on the wall sometimes and my roommate Angelique Kessels who taught me the art of juggling during the rest breaks of WorkPace CONTENTS Contents 1 2 Introduction Genetic Programming 27 Tntrod ction era e teren Sa OR RI LEE S gas 2 2 Representation e 22 Structure s a doae a A Mos ae A a 22 9 Nodes Los enol due Bae WEN ede io varden ee ae e 2d 2 3 GP algorithm ac oen sneer ae Oe A a A ee x d 23 1 Initialization ee 232 Eval astiOn ve este dedo Rode mE BS S uu E UA 2 3 8 Stop criterion vs ono mmo a ee 2 84 Selection 6050 4 xoxo sek de won e boe eu ex ud ig 2 8 5 Genetic operations leen GP research 3 1 System identification with GP o ee 3 2 Controller design with GP 3 3 Other applications eA Implementation of GP 4T EOUDray a u n as A odo Ia a 4 9 Matlab Simnni o 4 48 eten Sae ARA Se 4 3 Features of the GP program 22 4 4 Conversion of node to a Simulink block 4 5 DISCUSS so ds AWARD RUE UESTRE ee Ar pl ee eee Multidimensional optimization of a nonlinear function 5 1 A need for numerical optimization lees 5 2 Det
35. been done by others There are some libraries available on the inter net which are dedicated to Evolutionary Algo rithms A few examples are GPC created by A Fraser 1994 1 at the University of Sal ford and EO Evolving Objects by the Geneura Team at University of Granada 1998 10 In this thesis the EO library is chosen be cause it appears to be the most complete and flexible library of the considered libraries EO is a template based ANSEC compliant evolu tionary computation library description of EO 18 reference manual of EO 19 It was devel oped to implement evolutionary algorithms like Neural Networks NN Evolutionary Strategies ES Simulated Annealing SA Genetic Algo rithms GA and Genetic Programming GP 4 2 Matlab Simulink System identification and controller design will be performed with GP in this thesis A good choice has to be made for the representation of a solution lt is preferred to create models in a user familiar form so the result can be easily in terpreted and used for further application Also the evaluation of models should be easy and ro bust That is why it is decided to use Simulink models The so called Matlab Engine is shipped with Matlab This enables a C application to communicate with Matlab Simulink 3 7 This approach has a number of advantages The GP algorithm can be implemented in C instead of implementing everything in Matlab This will improve spe
36. can be unknown non linear noisy or contain narrow valleys so a ro bust optimization method is needed Suitable candidates are the stochastic methods Simu lated Annealing and a Genetic Algorithm Both are able to find global optima for a difficult fitness function A Genetic Algorithm has less problem dependent settings than Simulated An nealing Furthermore a useful example of a ge netic algorithm is already available and is easy to implement so this algorithm is selected for optimizing the parameters of a model This Genetic Algorithm has been imple mented in the GP program see Appendix C In the next chapter results of system identifica tion are shown 18 6 GP amp system identification 6 1 Introduction To examine the performance of our GP pro gram see Chapter 4 system identification is explored Defining a fitness function for system identification in the time domain is relatively easy so all effort can be put into the optimiza tion of the computer program A part of the field of system identification is trying to find the system dynamics by examining input output data black box identification From now on the system to be identified is called the plant Terminal and function sets The available nodes which represent Simulink blocks in a Simulink model can be divided into two groups the terminal set constant input and the function set add subtract multiply divide integrate differentiate
37. ce and stability problems Disturbance rejection With the above criteria GP is able to create con trollers which are optimized considering perfor mance robustness and the actuator limit The only thing left which is not considered is distur bance rejection The sensitivity transfer func tion is a measure for disturbance rejection This can be added in several ways Decreasing the sensitivity at low frequen cies will improve disturbance rejection in that region This can be added in the fitness function e g by stimulating con troller gain at low frequencies Stimulation of a higher bandwidth causes a larger frequency range to have an improved disturbance rejection This can be done by adding the term e bandwidth Inserting a disturbance signal will stimu late disturbance rejection because GP opti mizes the controller to be as fast as possible 24 7 2 Fitness function within the specified error band If the error is too big because of the disturbance signal the only possibility is to suppress it by in creasing the controller gain This is imme diately a problem for the optimization pro cess If other criteria in the fitness function do not allow just that then GP is unable to optimize the settling time which would be devastating for the whole optimization process Adding white noise is not smart because a controller can not provide a good distur bance rejection at all frequencies Improve ment of the
38. ckled with it and how to implement it For this a computer program has been written which is actually a structural and parametric optimization algorithm for dynamical systems using a GP GA hybrid GA Genetic Algorithm First some system identification cases are studied to get familiar with the GP program After that the GP program is transformed to perform controller design which is the main objective of this thesis Some simple problems have been studied to test and optimize the GP program Some promising results are achieved but there is still a lot of research left to be done e g with respect to nonlinearities discretisation and MIMO systems This thesis is a first exploration whether GP is suitable to create controllers for a user provided plant A futuristic idealistic goal is to have an application which can generate an optimal controller for an arbitrary plant However in achieving this goal THE main challenge is to define our control objective into a GP fitness function After that let evolution do the job Pressing the start button would be the only skill needed Samenvatting Na de opkomst van digitale regelaars zijn deze drastisch in prijs gedaald Een complexe regelaar structuur kan hierop eenvoudig geimplementeerd worden Echter het ontwerpen van een dergelijke regelaar kost meer inspanning Vaak wordt er voor een eenvoudige PID regelaar gekozen omdat deze eenvoudig en begrijpelijk is met een goede prestatie als resultaat T
39. ction This is desired because hundred of thousands of evaluations are done We are only interested in the best models 37 B FITNESS FUNCTION M FILE Fitness function for controller design The fitness m file is shown next The previously described parts are also indicated 1 function description function f system_model fitness gp model dirname XFITNESS returns a fitness value which indicates the performance of a simulink model according to user specified requirements Jusage by GP f FITNESSC gp model f system_model FITNESS usage by user f FITNESS gp_model dirname X gp model A simulink model with the name lt gp_model gt mdl It contains m inputs and 1 output block and can be used by the command SIM T X Y simCmodel t f t ul inputs 4 t timevector n by 1 matrix u inputvector n by m matrix outputs T timevector n by 1 matrix X respons of all states in model n by matrix Y outputvector n by 1 matrix If gp_model is set to temp_model the most recent numbered file gt temp_modelxxx ndl will be selected if it exists Useful for displaying intermediate results during a GP run dirname optional the gp model file will be search in this directory also the input output data is plotted f positive real value which indicates the performance of a simulink model A larger value of f
40. current fitness function has not been changed you can simply use fitness model7 test If model7 used system mdl exists in that directory it will be used else the system model defined in the fitness function will be used If you want to abort the GP run you do not want to wait for the whole generation and you do not want to save the best result you can easily press Ctrl C This immediately quits the GP run The Matlab command window will not be closed If you do another run it will be used again This does not cause any troubles 48 D GP PROGRAM D 1 User manual estimtf m An m file which determines the transfer function between an input and output of a simulink model by applying white noise to the input channel The transfer function is determined using the Matlab command TFE See estimtf m for required inputs and created outputs or look at a sample fitness file fitness m how it can be used simmodel m An m file which can run multiple simulations of a Simulink model All inputs of a model can be specified See simmodel m for all required inputs or look at a sample of a fitness file how it can be used setp_3c m An m file which generates a third order setpoint The adjustable options are amplitude maximum velocity maximum acceleration maximun jerk start time end time and sample time A time position velocity and acceleration signal is returned See setp_3c m for more info calcstop dll An S function which is
41. d uals have been created by using parts of both parents Figure 2 3 The usage rate of genetic operations can be set for each type Crossover is the most often used genetic operation in practice A usage of about 70 is common The impact on the individual s fitness is huge Big steps are made through the search space Although crossover has a negative effect on the fitness in over 75 of the cases it is useful in making leaps in the search space thereby avoiding ending in a local optimum Mutation Two mutation operations are point mutation and branch mutation Point mutation se lects one randomly chosen node from the struc ture of a parent and replaces it by a new node Figure 2 4 This node should have the same 2 3 GP algorithm y 3 u 2u 7 y w 2u 5 Figure 2 3 An illustration of a crossover opera tion amount of inputs to avoid connectivity errors Branch mutation selects one randomly chosen node and creates a new branch at that place These are only two examples of commonly used mutation operations Much more types of mutation operations exist The main condition is that they should replace some part of a struc ture with new randomly selected nodes A mu tation operation can change the size of the tree but not necessarily and furthermore it pre vents extinction of certain node types y u e u 2 Figure 2 4 An illustration of a mutation opera tion Mutation is the second most
42. determine settling time absE abs y 1 ur spec find absE A gt accuracy if length spec within spec all the way settling_time 0 elseif spec end n settling_time T else interpolate settling time spec spec end spec_interp absE spec A accuracy absE spec absE spec 1 settling time t spec 1 spec_interp t spec 1 spec_interp end a scaled settling_time penalty is added to the fitness if the position is not within the spec at the constant part constant find dur 0 constant constant end if settling_time gt t constant f f settling time T end add an extra scaled amplitude penalty if the position is not within the spec at time T if settling time gt T 9 f f max absE constant end A end controller output constraint maxU max abs y s out ControllerQut if maxU gt ActuatorLimit f f ConstraintPenalty end calculate sensitivity S sim_error EstimTF system model fs sens n sens n sens parts sens overlap s_in Amount s_in Controller0ut s_out Controllerdut MST if sim_error f BadSimPenalty else sensitivity constraint maxS max 20 log10 abs S if maxS gt SensLimit f f ConstraintPenalty maxS SensLimit end bandwidth optimization if bandwidth factor O CPabs abs 1 8 1 bw find CPabs lt 1 if any bw gt 1 amp CPabs 1 gt 1 bw bw 1 bw bw log CPabs bw log CPabs bw 1 CPabs bw Z
43. earch space it does not always find the desired solution after one run of the GP pro gram GP produces different solution every run because of its stochastic nature so it is possible that a better solution is found in another run 2 3 1 Initialization The GP algorithm starts with creating a new population The individuals are created ran domly using nodes of the function set and the terminal set The amount of individuals in a population is a tunable parameter Using more individuals gives a more diverse popula tion but is also more computational demanding The maximum allowed complexity can be set by limiting the tree depth Increasing the maxi mum tree depth increases the search space The search space should cover the complexity of the given problem but increasing it too much can cause convergence problems 2 3 2 Evaluation The performance of an individual has to be de termined This is done by defining a fitness function The fitness of all new individuals in a generation has to be evaluated This part is by far the most computational demanding com pared to the rest of GP algorithm The fitness function can be chosen freely A high level speci fication of the requirements can be implemented here A well defined fitness function is crucial in obtaining a good result with GP It shapes the landscape of the search space 2 3 3 Stop criterion A stop criterion is necessary or else the algo rithm will continue searching
44. ection to generate programs Dept of Computer Science University College London 22 W Langdon R Poli Better Trained Ants for Genetic Programming 1998 22 g gr 23 W Press S Teukolsky W Vetterling B Flannery Numerical recipes in C The Art of Scientific Computing Second Edition Cambridge University Press 2002 34 A S FUNCTIONS Al calcstop A S functions A 1 calcstop calcstop dll is an S function to be used in a Simulink model It tracks the simulation time and terminates the simulation if the simulation time is larger than the parameter value of the S function block The source code is shown below It was compiled in Matlab by typing gt gt mex calcstop c define S FUNCTION NAME calcstop define S FUNCTION LEVEL 2 include simstruc h include time h time_t begin const int request 1 static void mdlInitializeSizes SimStruct 8 ssSetNumSFcnParams S 1 if ssGetNumSFcnParams S ssGetSFcnParamsCount S return if ssSetNumInputPorts S 0 return if ssSetNumDutputPorts S 0 return ssSetNumSampleTimes S 1 ssSetOptions S SS OPTION EXCEPTION FREE CODE SS OPTION USE TLC WITH ACCELERATOR begin time NULL static void mdlInitializeSampleTimes SimStruct S static void mdl utputs SimStruct S int T tid time t now mxArray y ssGetSFcnParam S 0 double k mxGetPr y now time NULL if now begin gt k ssSetStopRequested S request
45. ed Also the available EO library can be used to imple ment GP which saves programming time All available tools of Matlab and Simulink can be used to evaluate models and evalu ate fitness cases Implementing algorithms in C for example to solve differential equations is not needed any more Mat lab Simulink has a lot of sophisticated tools 12 4 2 Matlab Simulink for that which are optimized during the years The result is a ready to run Simulink model which can be easily interpreted and can be used for further application Most controller engineers are familiar with the Matlab Simulink environment When combining C the EO library and the Matlab Engine for access to Matlab and Simulink all ingredients are available to cre ate a GP program which optimizes dynami cal systems according to a user defined fitness function This fitness function has to be defined in Matlab so the GP program will perform a specific task like system identification and con troller design The fitness function alone decides what task is carried out the GP program only optimizes dynamical systems according to this fitness function A global view of the architecture of the im plemented GP algorithm is shown in Figure 4 1 The source code can be found on a cdrom see 2 an example of a fitness function m file can be found in Appendix B Gp4so C apply genetic operators GP model l Far select simulink model par
46. en exceeded then f f 10 if the simulation time is too long then f 1019 if a simulation error occurs then f 1019 7 3 Results Some simple systems have been taken to let GP create controllers for This is done to test which combination of elements belongs in the fitness function A Pentium IV 2 4 GHz with 256 MB RAM is used for performing the GP runs The GP settings are equal to the settings used for system identification in Chapter 6 except for the population size which is set to 500 to have a more diverse population The other difference is that there are three inputs available instead of one 7 3 1 Controller design for a mass damper system System A mass damper system is examined which orig inates from a real setup a pick and place unit placed at an air bearing The transfer function of the system is written in Equation 7 5 and the bode plot is shown in Figure 7 3 _ 2 6 2 652 19s The system has got an actuator limit of 63 N Considering the purpose of the real setup a displacement of 0 3 m has to be made using a third order trajectory with a final accuracy of 5 um The maxima of the jerk acceleration and velocity of the trajectory are respectively 1000 m s 6 3 m s and 0 8 m s For these conditions a controller has to be designed 0 0006s P s 7 5 Settings A time of 2 seconds will be simulated with a sampling frequency of 5 kHz This time should be long enough so that GP can optimize the
47. enerations can be set A minimum fitness level can be set and the GP run can be aborted manually The GP program has got several tuning pa rameters These can be changed in a so called ini file see Appendix D The GP program is divided into three sec tions initialization performing the GP run and saving useful data such as the fitness progression the GP settings and the best model A GP run can take days to evolve If a fit ness improvement occurs then that model together with its fitness is saved to a file One reason is that the current best result can be shown during a GP run Another reason is that if the application crashes the most useful information is kept 13 4 4 Conversion of a node to a Simulink block Sometimes it is not possible to run a ran domly created Simulink model The main reason is numeric problems like singulari ties and stiff or even unstable models This leads to slow progression of the Simulink simulation or even abortion of it To catch these faults an S function has been written which stops the simulation after a speci fied number of seconds code is presented in Appendix A Further also a check is made after each simulation for errors us ing the Matlab command lasterr which returns the last occurred error If the sim ulation was not successful a high fitness is assigned meaning a bad individual The chances of selecting genetic operations are independent of each
48. ents evaluate fitness create population fitness function m file save best y stopcriteria Matlab Simulink Figure 4 1 A global view of the architecture of the implemented GP algorithm The implementation of the GP program and how it should be used is described in detail in Appendix D From now on one run of the GP program is called a GP run 4 IMPLEMENTATION OF GP 4 3 Features of the GP program A brief list of important features of the GP pro gram is listed next More information about the GP program can be found in Appendix D The representation used is a tree structure for both system identification Chapter 6 and controller design Chanter 7 LA DAL LA ULUL UL AL Yii UNL 1 The fitness function returns a value related to the performance of an individual To achieve this an individual in the GP pro gram has to be converted to a Simulink file which is then evaluated in Matlab Simulink resulting in a fitness value Available genetic operations are crossover branch mutation point mutation and branch addition A Genetic Algorithm has been implemented to boost the numerical optimization It is implemented as a genetic operation of GP The effect and motivation of this GA is dis cussed in Chapter 5 Parents are selected by Tournament Selec tion Both Steady State GP and Gener ational GP are possible The stop criterion is threefold The maxi mum number of g
49. er dir is requested dir_exists exist dirname if dir_exists set path current_dir cd addpath current_dir necessary enables the use of the plant and dl s cd dirname check fitness data file if exist cd gp_model _fitness log fitness_plot_data load cd gp model _fitness log elseif exist cd fitness log fitness plot data load cd fitness log else fitness_plot_data 0 end check if there is a plant saved during finalization of a GP run if exist cd gp_model used system mdl system model gp model used system use plant of saved result of a GP run end Xcheck if gp model temp model the most recent numbered gp_model will be selected if isequal gp_model temp_model D dir temp_model md1 B FITNESS FUNCTION M FILE PA A eee D D name if exist gp_model int2str length D mdl1 gp_model gp model int2str length D end end check if the gp model exists if not exist gp model return to the initial directory cd current dir error gp model does not exist end end 4 setup simulation model AAA I ANAANAANAYKAA AA ANA GET AND SETUP CONTROLLER try load_system gp_model catch end try add_block built in Inport gp_model ui 1 add_block built in Inport gp_model u2 add_block built in Inport gp_model u3 add block built in utport gp model y add line gp
50. er setpoint Ti s Settling time u m Plant input signal Ud m Disturbance signal Un m Sensor noise signal Ur m Reference signal w Weighting factor z m Position y m Plant output signal Ym m Output data of a simulation of a Simulink model Up m Output data of a simulation of a plant Abbreviations AD Analog Digital ADF Automatically defined function CFT Centre For Industrial Technology DA Digital Analog EA Evolutionary Algorithms EO Evolving Objects ES Evolutionary Strategies GA Genetic Algorithm GUI Graphical User Interface ITAE Integral of Time weighted Absolute Error m file Matlab file MIMO Multiple Output Multiple Input ML Machine Learning NN Neural Networks PID Proportional Integral Differential S function Simulink function SA Simulated Annealing SISO Single Output Single Input 33 REFERENCES REFERENCES References 1 A Fraser Genetic Programming in C University Of Salford 1994 2 D Bruijnen Cdrom of Automated Controller Design using GP 2003 8 D Hanselman B Littlefield Mastering Matlab 6 A Comprehensive Tutorial and Reference Prentice Hall inc 2001 4 D Vangheluwe Genetic Algorithm applied to Parameter Design A Case study on a Reliability Problem 5 G Franklin J Powell A Emami Naeini Feedback Control Of Dynamical Systems third edition Addison Wesley Publishing Company inc 1994 6 G Gray Y Li D Murray Smith K Sharman S
51. erministic methods a ee eee ee 5 3 Stochastic methods 222r 5 4 Choice of method eee GP amp system identification 6 1 Introduction i numo bem rue EROS wore qoe mes 6 2 Fitness funetion e 6 3 RESTS see hn da a ae A Ber hades ates eagle te de 6 3 1 System identification of a mass spring damper system 6 3 2 System identification of a fourth order system 64 Summary oues alen te eh eee Ga a a eo ea a ded es EARS Pee GP amp controller design TE Introduction 3 ad aen bei ee deep een ha dre 7 2 Fitness funetionis e e pem d doe ar OR eit Bom ah ar lele de ee aes MM cn PC 7 3 1 Controller design for a mass damper system 7 3 2 Controller design for a fourth order system TA SUMA a a A a Se TP ee A Conclusion and recommendations CONTENTS e SO o 0 00 OO OO OO A AA A 10 11 11 11 12 12 12 13 13 14 15 15 16 17 17 18 18 19 19 19 20 22 22 22 22 25 25 26 30 31 CONTENTS CONTENTS List of figures 32 List of symbols 33 Abbreviations 33 Bibliography 34 A S functions 35 A 2 Soo pow aee AAA 35 ADE Drobeeteddiy a o dcs vos exo d eode Be ORGS A TERCER DE ac RR 36 B Fitness function m file 37 C GA implementation 46 D GP program 47 DI User manuals aoe dae veelen or Br er le oe ee Ee ak un e f 47 Dit Gp4cd files ann the Ae eg ur mue ie Helen Arta Se Bey ad es 47 D 1 2 Description of files e 48 DIS
52. esis Although it would not be such a big step to transform the GP program to a nonlinear system identifier lt is just adding some nonlin ear blocks in the function set and choosing the input cases well to ensure observability 3 2 Controller design with GP Automatic Creation of both the Topology and Parameters for a Robust Controller by Means of Genetic Programming 16 In this paper a robust controller has been de veloped with GP for a second order system with no zeros Controllers are generated and exam ined using the SPICE simulator tool to solve differential equations The fitness function con sists of 10 components 8 components are based on an ITAE Integral of Time weighted Abso lute Error For these 8 components 4 differ ent parameter sets and 2 different trajectories 11 3 3 Other applications are used to obtain a robust controller The 9th component is the examination of the closed loop determined using an AC sweep A gain bigger than 3 dB in the low frequent area is penalized The 10th component is a stability test A small pulse is applied at the start A penalty is added to the fitness function when the output exceeds the pulse s amplitude at the end of the simula tion This is a stability test The major differences with this thesis are the fitness function and the simulation environment The same authors of the paper have produced some more papers about controller design with GP The applicati
53. f a mass damper system with time delay Controller transfer function The output compared with the trajectory lt lt ele Nyquist plot 2d in eee A 9 ue ld ata Open loop transfer function ee ee Schematic view of the Pato setup eee ee ee ee Transfer function fit of a Pato setup measured at the motor Transfer function fit of a Pato setup measured at the load Controller transfer function Time response seraa haw See te ae dae ee Bee oo ees were ss Nyquist plots uie a Poe deere leder e eae ee Bind v Hoe eed tg bee a Open loop transfer function 2 22222 ellen Controller transfer function ere ume response usa a a Rr CREE que Bou v eode rhe eS NSyquist plot nan seen ate bee ur E Weten EUR eae RE eee E RIS RUN Ss Open loop transfer function oeoa c e 32 LIST OF FIGURES 25 LIST OF SYMBOLS LIST OF SYMBOLS List of Symbols Symbol Unit Description A ma Amplitude of third order setpoint C N m Feedback controller model Cr N m Feedforward controller model e m Error signal Emas m Maximum error f Fitness value m H Amount of inputs n Amount of data points n C Maximum initialization value P m N Plant model s rad s Complex variable used by Laplace transforms S Sensitivity t s Time T H Complementary sensitivity or closed loop transfer function T s Total simulation time Te s End of third ord
54. fer function fit of a Pato setup measured at the load 3 002 108 0 004s st 3 94353 1 08 10552 P s 7 7 Both systems have a time delay of 4 ms this is because of the use of a DA converter and the time needed to compute the controller out put This is clearly visible in the transfer func tions and limits the performance to control the system Furthermore measuring at the load side has an extra bandwidth limitation it jumps from 180 to 360 Both models have been given to the GP program to create a controller for it Settings The objective is to move as fast as possible to another angular position with a maximum error of 0 1 A third order set point is used A time interval of 2 seconds will be simulated with a sampling frequency of 2 kHz The actuator limit is set to 2 V 7 GP amp CONTROLLER DESIGN 7 3 Results Results of measuring at the motor side GP returned a Simulink model with 8 transfer function blocks 4 gain blocks and 2 subtract blocks It uses the error signal and the velocity signal so GP has created a feedback velocity feedforward controller In Figure 7 11 the trans fer function of the feedback part is shown Magnitude dB Phase degree EN 2 is Frequency Hz Figure 7 11 Controller transfer function The resulting characteristics of the controller are listed below 20 days 0 12 gt 012 Settlin
55. field of linear con troller design Simulations have shown that it has a positive effect on the performance of GP More information about genetic operations can be found in 20 In the next chapter recent research is presented in the field of GP 3 GP research GP has existed only for about a decade Since the foundation was laid in the early 1990 s the popularity of GP has been increasing enor mously because of its promising results Also the fast speed increase of computing power is an important reason because GP is very com putational intensive A lot of books articles and papers have been dedicated to the subject Be cause of the flexibility of GP to tackle an ar bitrary problem it has been applied to a wide range of problems Already a lot of results have been achieved which are equal or better than results produced by experts in the field References that have been used here to gather information about GP are 15 20 and 21 Also a lot of information can be found on the internet e g 11 12 and 13 GP research can be broadly categorized into two groups 1 GP Techniques and Theory 2 GP Applications An overview is given from the latter only see Figure 3 1 because in this re search GP is used as a tool for system identifi cation and controller design For both topics an example is given in the next two sections They are compared with the approach in this thesis GP Research ARA ee GP Applica
56. for better results forever Two commonly used criteria are stop if a user defined number of generations has been reached or stop if the best of run fitness meets a user defined value Manual termination of the program is also possible 2 GENETIC PROGRAMMING 2 3 4 Selection If the stop criterion is not met a new generation has to be created Individuals parents are se lected to produce offspring An often used selec tion method is selecting parents probabilistically in proportion to their fitness Another popu lar selection method is Tournament Selection This is simply selecting the best individual out of a few randomly selected individuals from the population This is repeated until enough par ents are selected Two breeding policies are often used Gener ating a whole new generation using genetic oper ations is called Generational GP Replacing only a few individuals by newly generated in dividuals is called Steady State GP Other breeding policies can be found in 21 2 3 5 Genetic operations Genetic operations are used to create new in dividuals offspring out of existing individuals parents These operations are based on mech anisms seen in nature during evolution Com monly used genetic operations are crossover mutation and reproduction Crossover The GP structures of two selected parents are taken A randomly selected branch of both par ents will be exchanged Now two new indivi
57. g time Ts Tc 0 056 s Bandwidth 55 dB 0 52 V Maximum controller output Figure 7 12 shows that an overshoot of 8 occurs No penalty was included in the fitness function so if a better settling time can be ob tained with such an overshoot GP will do just that cm Oo M xim o NW or 1 t s Figure 7 12 Time response Figure 7 13 shows that stability and robust ness are preserved Figure 7 13 Nyquist plot The bandwidth can be determined from Figure 7 14 It also shows that the anti resonance resonance peak of the plant is sup pressed by a slant notch visible in Figure 7 11 Keep in mind that only the settling time is opti mized The notch suppresses the disturbance of the load at the motor It is obvious that suppressing this disturbance reduces the settling time of the first mass It is nice to see that GP came up with this without explicitly telling it to do that Magnitude dB Phase degree 10 10 10 10 Frequency Hz Figure 7 14 Open loop transfer function 7 GP amp CONTROLLER DESIGN Results of measuring at the load side GP returned a Simulink model with 3 transfer function blocks 1 integrator block 2 subtract blocks and 1 addition block It uses only the error signal so GP has created a feedback con troller In Figure 7 15 the transfer function of the feedback part is shown Magnitude dB Phase degree l piliui
58. gs results run info etc and will be saved at the end of a GP run matlabcom cpp matlabcom h This class provides all the functions needed for communication between the GP program and Mat lab Simulink It also converts a GP individual to a Simulink model node cpp node h This class represents a node in a tree or individual Such a node corresponds with a Simulink block But it is not exactly one on one For example a tree can have more of the same input but a Simulink model can only have one of each input which will then split into more lines from one inport block nodeproperties cpp nodeproperties h This class contains all characteristics of all node types All nodes written in gp ini should be defined in this class to be effective Defined characteristics are type name in gp ini amount of inputs if it is an input usage chance amount of parameters initialization of the parameters and needed Matlab code to generate the block in a Simulink model If you want a new Simulink block this class only needs to be changed and the gp ini should contain the block settings cpp settings h This class is more like a collection of most settings retrieved from gp ini and the CNodeProperties class At a lot of places the same settings are needed By concentrating all these settings it is easily portable and you only have to define a setting once 56 D GP PROGRAM D 2 C code description fitness h This template class determine
59. he result is harder to interpret GP s main task is optimizing a structure of blocks A dynamical model can contain a lot of parameters These parameters are randomly created Because GP is not very good at the optimization of those parameters it must be as sisted by a dedicated algorithm for optimizing the parameters of a model In this thesis a GA is used to perform this task The performance of this algorithm is crucial to obtain a good result within reasonable computing time The param eters should first be optimized to be able to opti mize the structure This is very time consuming and every gain at this is desirable Another opti mization algorithm could be tested to compare it with GA A suitable candidate is Simulated Annealing combined with the Nelder Mead Sim plex method In literature it has achieved some promising results which indicates that it also could be useful for our problem In this research only linear system identifica tion and linear controller design have been ex plored The GP program is the same for both It is a parametric and structural optimization al gorithm for non linear dynamical models The only difference is the fitness function Other fields of application are Nonlinear control e g gain scheduling Identification of nonlinear systems Signal processing e g synthesis of dedi cated filters LIST OF FIGURES List of Figures 2 1 2 2 2 3 2 4 2 5 3 1 4 1 4 2 6 1 6 2
60. hermore a constant speed trajectory is used to avoid nonlinearities like coulomb friction and cogging which would influence the measure ments Noise is added to the controller output and the motor input is measured The trans fer function of this input output is called the sensitivity The system may not change its ro tating direction because of this noise else the nonlinearities will disturb the measurements Assuming that the implemented controller ap proximates the theoretical continuous time con troller the transfer function of the plant is de termined by calculating P C71 S71 1 see Equation 7 3 A sampling frequency of 500 Hz is used In Figure 7 9 the transfer function is plotted measured at the motor side grey line The black line is a fit of the transfer function which is shown in Equation 7 6 Dod Measurement lel Magnitude dB 10 10 10 Frequency Hz b Ge C o o 300 I R e e Phase degree Frequency Hz Figure 7 9 Transfer function fit of a Pato setup measured at the motor Ps 5000s 1 47 10 s 3 002 10 0 004s si 3 94353 1 08 105s 7 6 In Figure 7 10 the transfer function is plot ted measured at the load side grey line The black line is a fit of the transfer function which is shown in Equation 7 7 Measurement Model Magnitude dB Frequency Hz 10 10 10 Frequency Hz Figure 7 10 Trans
61. ijdens dit afstudeeronderzoek is een programma ontworpen waarmee regelaars gegenereerd kunnen worden voor een grote diversiteit aan systemen Het is gebaseerd op Genetisch Programmeren GP dat regelaars cre ert met een variabele structuur Niet lineaire dynamische elementen kunnen toegevoegd worden aan de regelaar en of het te regelen systeem Verder is het mogelijk om eisen te specificeren in het tijd en frequentie domein De eisen worden vertaald in een fitness functie die probleem specifiek is Deze fitness functie retourneert een waarde die aangeeft hoe goed aan de gestelde eisen voldaan wordt Deze fitness functie wordt vervolgens geoptimaliseerd door GP wat dus het probleem oplost GP is uitvoerig onderzocht betreffende de werking de problemen die ermee aangepakt zijn door an dere onderzoekers en hoe GP geimplementeerd kan worden Vervolgens is er een computer programma geschreven wat in feite een structuur en numerieke optimalisatie algoritme is voor dynamische syste men met als achterliggend algoritme een combinatie van GP en GA GA Genetisch Algoritme Ten eerste is er systeem identificatie mee uitgevoerd om vertrouwd te raken met het programma Daarna zijn enkele wijzigingen aangebracht zodat het geschikt is voor het ontwerpen van regelaars Dit is het hoofddoel van dit afstudeeronderzoek Enkele eenvoudige problemen zijn voorgelegd aan het GP programma om het te testen en aan te passen waar nodig is 4 CONTENTS CONTENTS Goede re
62. in a so called ini file gp ini PopSize is the size of the population Every generation the best parents are selected with an amount of Parents For each new individual an amount of Recombi parents Recombined parents are randomly selected out of the best parents A new individual is created by applying a crossover operator at these recombined parents An amount of PopSize new individuals are created The standard deviation of the search step is multiplied or divided equal chance by Alpha for each value of each individual This enables adjustments to the order of magnitude because the standard deviation of the search step can get a factor Alpha bigger or a factor Alpha smaller A choice of 1 5 for Alpha is the default A step is taken for each value of each individual by drawing a number from a normal probability distribution using its new standard deviation of the search step There is a standard deviation parameter for each value of an individual because the order of magnitude of those values can vary a lot The standard deviation of the initial step is determined by multiplying the absolute value of the parameters in a model by InitStepSizeFactor Do not set this value too high else it makes too big jumps in the first generations so it is not able to optimize locally This is the case when it is set to 1 A value somewhere between 0 01 and 0 1 works fine Every generation the best individu
63. interpolation on dB scale bandwidth bw 1 Clength S 1 fs 2 elseif CPabs 1 gt 1 bandwidth fs 2 else bandwidth 0 end f f exp bandwidth_factor bandwidth end end 41 B FITNESS FUNCTION M FILE 6 validity check of the fitness value hhhhhbbbhhhhhh FINALIZATION check validity of f if not exist f f BadSimPenalty elseif not isfinite f f gt BadSimPenalty f lt O not imag f 0 f BadSimPenalty end 7 visualisation of the results plot input output data if a a directory is specified if dir exists Kreturn to the initial directory cd current dir remove GPmodel from closed loop delete line system model ddu r 1 GPmodel 1 delete line system model du r 1 GPmodel 2 delete line system model e 1 GPmodel 3 delete line system model GPmodel i D 2 delete block system model GPmodel determine bandwidth if it isn t calculated before if f BadSimPenalty CPabs abs 1 S 1 bw find CPabs lt 1 if any bw gt 1 CPabs 1 gt 1 bw bw 1 bw bw log CPabs bw log CPabs bw 1 CPabs bw interpolation on dB scale bandwidth bw 1 length S 1 fs 2 elseif CPabs 1 gt 1 bandwidth fs 2 else bandwidth 0 end end plot info to screen disp Genetic Programming for controller design disp FITNESS FUNCTION 1 disp Created by DJH Bruijnen I disp TUE Philips CFT 1 disp
64. is a list with all supplied files They are listed relative to the root dir of the application The default of the simworkdir is and the default of the modelsavedir is mdls Files which do not need to be altered are gp4so exe lt simworkdir gt estimtf m lt simworkdir gt simmodel m lt simworkdir gt setp_3c m lt simworkdir gt calcstop dli lt simworkdir gt protecteddiv d1ll Files which can be altered by the user according to his desires are gp ini lt simworkdir gt fitness m file name is specified in gp ini Some working sample files are lt simworkdir gt fitness m lt simworkdir gt plant_2nd_order md1 lt simworkdir gt fitness_pato m lt simworkdir gt pato_motor mdl lt simworkdir gt pato_load mdl Created files during a GP run are gp_model a temporary Simulink model only in memory which may not be saved lt simworkdir gt lt modelsavedir gt temp_modelxxx mdl lt simworkdir gt lt modelsavedir gt fitness log Created files after a GP run are lt simworkdir gt lt savedir gt lt savename gt mdl lt simworkdir gt lt savedir gt lt savename gt _fitness log lt simworkdir gt lt savedir gt lt savename gt _fitness m lt simworkdir gt lt savedir gt lt savename gt _used_system md1 lt simworkdir gt lt savedir gt lt savename gt _logdata txt lt simworkdir gt loghistory txt 47
65. is equal to a worse performance X system model optional This will return the system model name without simulating anything Koninklijke Philips Electronics N V Centre for Industrial Technology CFT Created 01 2003 By DJH Bruijnen 2 adjustable parameters LALAAAAAAA Ahh hhhhhhihhhhhhkhhhhhhhhhhh PARAMETERS TO BE EDITED BY THE USER system parameters system model pato_load user specified plant must be in the same directory s_out PlantOut 1 koutport at plant output s out ControllerOut 2 koutport at controller output s_in PlantOut 1 Z inport with sensor noise signal s_in ControllerOut 2 Zinport with disturbance signal s_in RefAcc 3 inport with reference acceleration signal s_in RefVel 4 finport with reference velocity signal s_in RefPos 5 Jinport with reference position signal s_in Amount length struct2cell s in amount of inports of the system controller parameters c_out ControllerOut 1 outport at controller output c_in RefAcc 1 jinport with reference acceleration signal c_in RefVel 2 Linport with reference velocity signal c_in PosError 3 inport with position error signal c_in Amount length struct2cell c_in amount of inports of the controller 38 B FITNESS FUNCTION M FILE A NN E aT eee ftrajectory parameters dt 002 sample time for 3rd order trajectory T 2 simulation time A 2 pi amplitude of trajectory jmax
66. is used A weighting fac tor w is applied to scale the order of magni tude of the fitness contribution for input case j with respect to the other input cases The total fitness function is the sum of all input cases In equation 6 1 the fitness function is shown yp is vector with simulation data of the plant for input case j This is a discrete time measure ment with n samples ym is the accompanying output of the GP model ml ypj i 6 1 i l 6 3 Results Some simple systems have been taken to be iden tified by the GP program An Intel Pentium IV 2 4 GHz with 256 MB RAM is used to perform the GP runs 6 3 1 System identification of a mass spring damper system System The Simulink model shown in Figure 6 1 is used to create input output data It is a mass spring damper system with a damping coefficient of 0 2 and a natural frequency of 10 rad s In Fig ure 6 3 the three input cases are shown The measured output is contaminated with random noise to make the problem more realis tic 19 6 3 Results s2 2 0 2 10s 10 2 input mass spring damper output Random Number var 0 005 Figure 6 1 Plant to be identified Settings The main settings of the GP program are shown in the table below data points time span inputs step at 0 3 s 0 8 Hz sine 3 2 Hz sine gain add subtract integrate 1st order TF 2nd order TF maximum tree depth population size function set
67. ly behaves like a second order system The three second order transfer func tion do not exert a lot of influence on the result because of their low gains s2 6 149925 36 3865 ni4 ni3 n12 Figure 6 6 Best GP model input u output y grey y black emory y 1 2 3 2 1 2 t s t s t s Figure 6 7 Comparison of the plant and the best model using the same input data as in the fitness function Validation The other input data is also applied to the GP model see Figure 6 8 For this input data the output of the GP model and the plant are also similar Unfortunately the anti resonance resonance peaks were not caught This is clearly visible when comparing the trans fer function of the plant and the model see Fig ure 6 9 The GP run has been performed sev 21 eral times but a good model structure includ ing the anti resonance rosonance peaks has not been found input u output y grey y black I 1 05 0 02 o 05 0 05 ud E o 0 1 2 3 o 1 o 1 2 3 2 15 0 02 o 1 o 05 0 02 0 1 1 2 3 o 1 2 3 2 3 o 1 2 3 0 2 3 0 1 2 3 eror Y Y 10 3 0 05 2 0 1 10 o 0 05 1 2 3 o t s tls s Figure 6 8 Validation of the best model using different input data e 50 Magnitude dB i Phase degree i 8 8 10 Frequency Hz Figure 6 9 Transfer function gain phase and coherence function The reason that it does not find the co
68. model u1 1 GPmodel 1 add line gp model u2 1 GPmodel 2 add line gp model u3 1 GPmodel 3 add line gp model GPmodel 1 y 1 addterms gp model catch end EDIDI Ahh hhhhkhhhhhhhhhhhh CREATE CLOSED LOOP SYSTEM try bdclose system model catch end try load system system model catch end add block gp model GPmodel system model GPmodel add line system model ddu r 1 GPnodel 1 add line system model du r 1 GPmode1 2 add line system model e 1 GPmode1 3 add line system model GPmodel 1 D 2 5 fitness determination hhhhbbbhhhhhhietbht tb DETERMINING FITNESS f 0 Ydetermine zeroes poles and gain of the controller only linear models gp_sys linmod gp model gp z gp_p gp k ss2zp gp_sys a gp_sys b gp_sys c gp_sys d c_in PosError check stability if any gp z gt 0 any gp_p gt 0 f BadSimPenalty else create 3rd order trajectory ur dur ddur setp_3c A vmax amax jmax dt 0 T n m size ur Y determining size of trajectory matrix t O m 1 dt time vector n by 1 matrix sensor noise and disturbances un 2 rand n m 1 noise ud 2 rand n m 1 disturb 40 B FITNESS FUNCTION M FILE B e en end Xcalculate response of the closed loop ty sim error simmodel system model t un ud ddur dur ur m MST if sim_error f BadSimPenalty else end
69. nals are user provided inputs and constants The function set and the terminal set can he chosen freely provided that USL ittiias DOU CALL VU CLODO de LY pEBEIVVANVA V no syntax errors occur when nodes are connected arbitrarily 2 3 GP algorithm The GP algorithm starts with creating the first generation of a new population containing a user defined amount of randomly generated pro grams also called individuals The fitness of these individuals is evaluated The GP algo rithm stops if a termination criterion is satisfied If not parents are selected in a probabilistic way according to their fitness A new generation is created by applying genetic operations The fit ness of the individuals in the new generation is evaluated The loop is now closed and is con tinued until the termination criterion has been reached An overview of this GP algorithm is shown in Figure 2 2 Figure 2 2 An overview of the GP algorithm To set up GP for a particular problem the main five items which have to be defined are Function set Terminal set Fitness function Evolutionary parameters Termination criterion 2 3 GP algorithm When these items have been defined well and a solution for the particular problem exists then there is a realistic chance that the GP algorithm finds a good solution Mainly it is a matter of time and resources Because GP is a stochastic search technique to find an optimum solution in a wide s
70. nd only for the new vertice the function is evaluated Then ordering takes place as above and the algorithm searches until no further improvement is gained The algorithm performs well for finding a lo cal minimum Also narrow valleys are no prob lem the simplex first moves to the valley and then stretches itself in the direction of the nar row valley 5 MULTIDIMENSIONAL OPTIMIZATION OF A NONLINEAR FUNCTION 5 3 Stochastic methods Simulated Annealing Simulated Annealing is inspired by the cooling process of metals When adjusting the cooling scheme of a metal in the right way properties of the metal can be determined As a metal cools the atoms organize themselves into an ordered minimum energy structure At a high temper ature the atoms move a lot At a lower tem perature they move less which results in a fixed structure below some temperature level If the temperature decreases slow enough the atoms have time to arrange themselves in a fine struc ture low energy Decreasing the temperature very fast results in a coarse structure high en ergy Simulated Annealing tries to imitate this pro cess It tries to find the global optimum of a mul tidimensional function by decreasing the tem perature slow enough When evaluating the fitness with a different state this temperature determines the rate of accepting the state if it has a worse fitness A better fitness is always accepted At a high temperature a worse fi
71. ndomly creates a new individual It builds a tree by selecting available nodes using the usage chances Furthermore it initializes all existing parameters of an individual according to their probability distribution defined in CNodeProperties eosgatransform2 h This template class applies all genetic operations to produce new offspring This class has been changed from the original one because it only could handle two genetic operations D 2 3 Known bugs possible improvements experience The class FitnessType is not used as it is meant to be Its purpose is to define whether GP should minimize or maximize the fitness value In the whole program minimizing the fitness value has been assumed The fitness is represented as a double data type and not as a more sophisticated class like FitnessType The authors experience is that the EO library is very big and looks like spaghetti for a non OO programmer the first month So it is very difficult to use it as supposed to A result is that some parts of the code are not very nicely coded because you have to stick to the way the EO library is coded The author was previously not used to program object oriented At some places this is visible like concentrating all settings in CSettings instead of localizing those settings in the places where they are actually needed 57
72. ness function The fitness function represents the user s requirements For controller design the main items are stability performance and ro bustness All parts of the application have to be constructed properly else the chance to get a good result decreases drastically Recommendations Until now only controllers have been designed for a second and a fourth order system includ ing time delay This research can be contin ued by refining the fitness function e g in corporating discretisation such that the con troller will perform well on an experimental setup Also some more complex systems should be tested to explore the generality of the fitness function We think of systems with coulomb friction and higher order systems The fitness function should be defined in such a way that it can reflect the control objectives for a wide range of systems Another possible implementation is to design controllers in the discrete z domain using a net work structure This has some advantages with respect to the continuous s domain Blocks like TF1 and TF2 are not needed anymore These can be constructed using elementary blocks As 31 seen before GP prefers the use of elementary blocks so the genetic operations can operate bet ter Furthermore the design of a discrete con troller is closer to the implementation on an ex perimental setup It does not have to be discreti sized like continuous controllers A disadvantage is that t
73. ness in the population the average fitness of the population the amount of fitness evaluations and the generation number Another possibility for examining intermediate results is using another Matlab command window Then you do not have to press the Pause button however you also can not examine the current model Only the best models which are saved in the directory mdls can be examined If you see that the GP run has converged to a certain level or if you just want to abort the GP run then you can press any key on the keyboard After a generation is completed you will be asked if you really want to quit By typing y you get the option to save the model by typing the filename and sub directory If the file already exists a message appears if you want to overwrite it Separetely from the model other data is saved so the settings of GP are kept Also the fitness function and the system model are saved so the performance of the generated GP model can be evaluated at any time independently from the current state of the fitness function and the plant system model Finally the fitness progression is saved With all this data it is possible to initialize the GP program equally to an earlier state to reproduce the result For example if you have saved the model as model in the directory test you can examine the model with the original fitness function and system model as follows cd test model7_fitness model7 If the
74. oblem GP and GA are about telling the computer what needs to be done and not how to do it GA and GP are inspired by natural evolution The body of research of algorithms inspired by natural evolution is called Evolutionary Al gorithms EA 9 which is a part of Ma chine Learning The idea is that a population evolves using genetic operations like crossover mutation and reproduction The selection of in dividuals to produce offspring depends on their fitness The fitness is a number which indicates how good an individual satisfies a problem Roughly the fittest individuals survive and new possibly better individuals will be created Fitness weighted selection and genetic opera tions make it possible to create offspring with a better fitness This mechanism can be continued until the performance of the fittest individual is satisfying Just like in nature it is survival of the fittest 2 2 Representation GP produces computer programs To be able to solve problems the computer program needs one or more inputs and an output The con 2 2 Representation nection between the inputs and output must be syntactically correct Furthermore small parts of that structure have to be exchangeable so ge netic operations can be applied 2 2 1 Structure Koza 17 proposed to use a tree structure con sisting of functions and terminals This way ar bitrary programs can be built Modifying a pro gram is easily done by subs
75. ons are different but the fit ness functions is about the same 3 3 Other applications GP has also been applied to a lot of other appli cations A few examples about totally different applications are given next Using Genetic Programming to find Lya punov functions 15 This paper is about finding Lyapunov func tions which are used to prove stability for non linear systems Better Trained Ants for Genetic Program ming 21 In this paper an artificial ant which has to follow the Santa Fe trail is programmed The Santa Fe trail is a virtual path on a check board with food at specific locations The ant should follow the optimal path and eat all the food Genetic Programming Applied to Traffic Control 14 In this paper a set of crossing roads with traf fic lights at the intersections have to be con trolled The problem of controlling the traffic lights optimal in terms of traffic rate and wait ing time of a car at an intersection is presented to GP The working of GP has been studied and the fields where GP has been applied is investigated In the next chapter the implementation of the generalized GP program is discussed This ap plication can be used for both system identifica tion and controller design 4 IMPLEMENTATION OF GP 4 Implementation of GP 4 1 EO library It is possible to build the whole GP program from scratch This would result in a lot of work which has already
76. ontroller is more difficult In this thesis a method is devel oped to let an evolutionary computer program generate such an atypical controller The user only provides the plant in this case a Simulink model and the desired requirements and the computer program will generate a suitable con troller The mechanism used in this computer pro gram is called Genetic Programming GP It is an algorithm that optimizes a variable structure according to a user specified fitness function From an engineering point of view this fitness function is the most challenging part of GP The fitness function and the rest of GP are explained in Chapter 2 The first objective of the graduation project is to investigate how GP works The second objec tive is to implement it and the last objective is to assess the GP potential by performing system identification and controller design The system Philips CFT Mechatronic Equipment Motion Control Group identification examples serve as a test case for GP whereas controller design is our main ob jective G Naus 8 has continued the research of our GP based system identification As this is our first exploration of GP our focus is on linear systems Therefore it will be easy to judge GP results in view of the well known results of linear system identification and con trol Information about GP is gathered to exam ine its working principle see Chapter 2 and 3 Next a generic GP computer p
77. ontrollers which can be implemented immediately on a real setup because during the 30 7 4 Summary tests of GP sensor noise was not included non linear effects like coulomb friction were not in cluded and the discretisation on a real setup is not taken into account This leads to controllers that emphasize the sensor noise too much and sometimes place poles and zeros above half the sampling frequency resulting in aliasing Dur ing tests on a real setup this was observed For future research this should also be penalized in the fitness function In the fitness function the sensitivity function is determined numerically However for linear systems it is also possible to determine the sensi tivity function analytically using the linmod command in Matlab This would improve speed and accuracy This was not done because some evaluated plants already contained a non linear element namely a time delay block When using linmod this time delay would not be included in the analytical model Idea s for future research One type of controller specifications was used namely follow a 3th order trajectory and be as fast as possible within a certain error band at the constant part of the trajectory Further more stability and robustness have to be pre served Such a demand is typical for a pick and place unit Other types like following a constant speed trajectory printer or following an arbi trary trajectory rob
78. ot arm as good as possible could be examined Also other plants which are more complex higher order nonlinearities etc can be con sidered This should be done to test the per formance of the GP program and where neces sary adjustments should be made to the fitness function to improve its generality In this thesis all models are continuous time models An interesting variation is to define all models as discrete models An advantage is that a discrete controller is more close to the imple mentation of it on a real setup Furthermore GP performs better with elementary blocks 8 CONCLUSION AND RECOMMENDATIONS 8 Conclusion and recommen dations A GP program which optimizes a dynami cal model according to user specified require ments has been implemented successfully It has been tested by performing system identifi cation and controller design with it for a few simple systems Some promising results have been achieved by the GP program It can be concluded that Genetic Programming in com bination with a numerical optimization method for the parameters of the GP model is suitable for controller design In this thesis a Genetic Al gorithm has been used to optimize the parame ters of a GP model The application is still under development Improvements have been achieved by fixing some minor bugs by tuning the GP algorithm for better convergence to a good solution and most important by defining a suitable fit
79. r rect structure in this case is observability For all three fitness cases the time responses of the plant and the GP model coincide so GP is not able to find the correct structure in the first place To overcome this problem the fitness function has to be defined better A possibil ity for linear models is comparing the transfer functions of the plant and the GP model 7 GP amp CONTROLLER DESIGN 6 4 Summary GP is able to perform system identification for simple systems The structure of second order systems presented in Section 6 3 1 was exactly reproduced though small numerical variations were present lhe fourth order system could merely be approximated The exact structure was never found The best GP model mainly be haved like a second order system which did not catch the anti resonance resonance peaks The cause of this is that the dynamics of the anti resonance resonance peaks is barely observable using an impulse a step and a sine Looking at the frequency response that peak is only 196 of the low frequent gain so it is obvious that it can not be caught definitely not with contami nating the output with noise A solution would be to compare the models in the frequency do main This research has been continued by G Naus 8 In the next chapter results of controller design with GP are described The only difference with system identification is the fitness function Its composition will be described in detail
80. ram And secondly all important settings are saved so in a later stage a similar GP run can be performed D 1 3 Node types Below the Simulink blocks are presented with their characteristics uniform 0 n means that the parameter is initialized by drawing a random number out of a uniform probability distribution between 0 and n n is the MaxInitValue specified in gp ini loguniform i n means that a number is drawn between 4 and n such that the chance that c1 gt co is equal to c1 lt a with c and cz drawings from a loguniform probability distribution between i and n Add 5 u t u t ual Subtract SERE y t wi t valt Multiply ius y t ux tJua t Divide protecteddiv Function t S Functio y t E a A protected division S function is introduced because a normal division block would produce division by zero problems in case of u t 0 Integrate Integrator y s ier tu s Differentiate ze 46 Differentiator y s u s c uniform 0 n 52 D GP PROGRAM D 1 User manual This block and also all TF blocks use such an initialization so the largest initial pole is n It is wise to chose n somewhere below the sampling frequency in rad s Note that this does not restrict creation of poles or zeros above a frequency n due to numerical optimization The standard differention block in Simulink du dt was not used because of numerical problems The cause of this is
81. rates a con troller which is immediately within the error band at the constant part of the trajectory This results in a fitness value of zero so the GP run is finished by fulfilling all requirements After this a criterion to improve disturbance rejec tion was added to the fitness function namely a bandwidth criterion so GP kept on improving the controller Furthermore when only applying the feedfor ward part the output of the mass damper plant came only halfway the prescribed setpoint It is not fully optimized because GP has created a feedback and feedforward controller simulta neously Separating these parts in the fitness function could be considered in future research Several attempts to incorporate disturbance rejection have been tried but no good result was obtained When inserting white noise or low frequent noise GP was not able to optimize the settling time because the error band could not be reached Reducing the amplitude of the noise did not work because then it did not have any effect The produced controllers were worse than without adding a disturbance signal The bandwidth criterion only works if GP is able to fulfil all other requirements so this is also no good general criterion A recommendation for further research is that the user should constrain the sensitivity or the controller gain for a chosen frequency to incorporate disturbance rejection The fitness function is not yet well defined to produce c
82. ring runtime and will not be saved on the hard disk Anyhow it may not be saved else a problem of equally named models occurs So look out when checking intermediate results Do not save the Simulink models when closing them temp_modelxxx mdl fitness log When a better model has been found it will be saved in the mdls subdirectory The model is saved as temp_modelxxx mdl with xxx a successive number Also its fitness and fitness evaluation count is appended to fitness log With this the fitness progression can be examined later together with the structure of accompanying models when fitness improvements occur lt savedir gt lt savename gt mdl lt savedir gt lt savename gt _fitness log lt savedir gt lt savename gt _fitness m lt savedir gt lt savename gt _used_system mdl lt savedir gt lt savename gt _logdata txt loghistory txt 51 D GPPROGRAM D 1 User manual If you save the best model at the end of the GP program the files above are saved The model itself and the fitness progression file are copied from the mdls subdirectory The GP model is the temp modelxxx mdl with the highest number Furthermore the fitness file and the sys tem model are saved from the work directory Finally some settings and results are saved into lt savename gt _logdata txt and this is also appended to loghistory txt All these files are saved for two reasons The result can be evaluated independently of the current state of the GP prog
83. rion to the fitness function is shown in Equation 7 4 The first line adds the scaled settling time to the fit ness f The second line is added to add more nu ance to the fitness function for less good models Minimizing this fitness results in an optimum settling time The way the fitness landscape is shaped determines the evolution capabilities of GP 7 GP g CONTROLLER DESIGN Stability and robustness Aside from the performance criterion there are some more criteria desirable like a stabil ity robustness criterion and the actuator limit Both criteria are implemented as a constraint If a constraint is exceeded a penalty of 1 105 will be added to the fitness The robustness criterion is included as to 6 dB This prevents that the open loop does not come too close to the point 1 in a Nyquist plot This criterion solely is not enough to guar antee stability because the line in a Nyquist plot could pass the point 1 on the wrong side This would result in an unstable closed loop behavior Together with the previous time domain perfor mance criterion stability is guaranteed because such an unstable model would blow up and be caught by the performance criterion Ritin 1 Max AMIS the maximum sensitivity Actuator limit Every actuator can provide a finite amount of power Limiting the plant input is needed to prevent creation of unrealistic controllers with too high gains Applying such controllers would result in performan
84. rogram has been written see Chapter 4 which optimizes dy namical systems The first tests with system identification pointed out that GP had diffi culties with optimizing the parameters of a model To improve this a numerical optimiza tion method was added to the GP program see Chapter 5 In Chapter 6 results of perform ing system identification with GP are shown Finally controller design has been carried out which is described in Chapter 7 2 GENETIC PROGRAMMING 2 Genetic Programming 2 1 Introduction Genetic Programming GP 20 is a part of a very large body of research called Machine Learning ML This term was first introduced by Samuel Samuel 1963 who meant comput ers programming themselves with it A good contemporary definition of ML is The study of computer algorithms that improve automat ically through experience Mitchell 1996 GP aspires to do precisely that inducing a popula tion of computer programs that improve auto matically as they experience the data on which they are trained GP was introduced by Koza in 1992 17 The underlying basis is Genetic Algorithms GA 4 The major difference between GA and GP is the representation of the solution GA returns a fixed length sequence of binary numbers which has to be interpreted GP returns a ready to run program with a variable length This is the power of GP it is able to adapt itself to the complexity of the given pr
85. rts sens overlap c in Amount c in PosError c out Controller ut MST freq linspace 0 fs 2 length S plot results if the simulation was successful close all if not sim_error figure FIRST FIGURE Zur and y subplot 221 plot t constant t constant O A color 7 7 7 hold on plot t ur color 7 7 7 plot t y 1 k axis tight title r and y xlabel t s ylabel x m wy ur subplot 222 plot t constant T A accuracy A accuracy Araccuracy A accuracy color 7 7 7 hold on plot t y 1 ur k axis tight title error y r xlabelCt s ylabelCx m plant input subplot 223 plot t y 2 k hold on plot 0 Tl ActuatorLimit ActuatorLimit ActuatorLimit ActuatorLimit color 7 7 71 axis tight title plant input xlabelCt s ylabel F N plot fitness data subplot 224 if not isequal fitness_plot_data 0 fitness plot data find fitness plot data 2 0 2 1e 15 n size fitness plot data 1 x zeros 2 n 1 1 y zeros 2 n 1 1 Zsetup plot data 43 B FITNESS FUNCTION M FILE for i 1 n 1 x 2 i 1 2 i fitness plot data i i 1 1 0 1 y 2 i 1 2 i fitness plot data i 2 end x 2 n 1 fitness plot data n 1 y 2 n 1 fitness plot data n 2 plot results miny min fitness plot data 2 2 maxy max fitness plot data
86. s the fitness of an individual So it firstly creates a model then it simulates that model and finally returns the fitness value To do all this it makes use of the CMatlabCom class eogaoperator h This template class implements the GA optimization algorithm It is implemented as a genetic operation The only difference is that it creates models and does fitness evaluations The other genetic operations do not do that They only change a part of the individual and after all genetic operations are applied the new fitness of all offspring is calculated eogencontinue2 h This template class is a changed file from the EO library because some more features are required Its mainly about adding more termination criteria and that a confirmation will be asked if a termination criterion has been reached eomutationoperators h This header file contains some modified genetic operations eoPointMutation2 eoSubtreeXOver2 eoBranchMutation2 and a new genetic operation eoBranchAddition A reason that the first three operators are changed is that the same node type does not have to be identically because it can contain numerical values Another reason is that there was an inconsistency in the EO library When initializing a tree its depth is limited this has always been assumed and can be set in gp ini When the original genetic operation was applied the amount of nodes was limited instead of the depth eoparsetreedepthinit2 h This template class ra
87. s thrown away Now another parabola is fitted through these three points and so on This method works good for sufficiently smooth functions Golden Section Search is designed to handle in effect the worst possible case of function op timization It starts with finding a bracketing triplet This means finding three points with the middle point a lower fitness as the other two Then another point is chosen at the golden sec tion fraction 3 V5 2 0 38197 between the outer points This is the most optimal choice considering the worst case of function optimiza tion A derivation of this can be found in 23 Both Golden Section Search and Parabolic In terpolation are only designed for local optimiza tion near the starting point The mini ANAL Conjugate Gradient method A common error is to assume that any reason able way of incorporating gradient information should be about as good as any other This line of thought leads to the following not very good algorithm the Steepest Descent method Opti mize in the direction of the gradient and do the same at the found optimum For long narrow valleys it will perform many small steps which is not very efficient 16 5 2 Deterministic methods A smarter way of choosing new directions is the concept of non interfering directions more conventionally called conjugate directions A mathematical description can be found in 23 The conjugate directions are chosen to
88. set tling time An extra element has been added to 25 Magnitude dB Frequency Hz Figure 7 3 Transfer function of a mass damper system with time delay the fitness function f f e endwidth which stimulates a higher bandwidth The provided spec was relatively easy to reach Without the bandwidth criterion the fitness would be zero and the controller would not be optimized any more although the controller could still be op timized a lot with respect to disturbance rejec tion Results The main characteristics of the produced con troller are listed below Settling time Ts Te os GP returned a Simulink model with 2 trans fer function blocks 5 gain blocks 5 subtract blocks and 6 addition blocks It uses the error signal the velocity signal and the acceleration signal so GP has created a feedback veloc ity acceleration feedforward controller In Fig ure 7 4 the transfer function of the feedback part is shown When the feedforward part is applied solely the output moves about halfway the trajectory 7 GP CONTROLLER DESIGN Magnitude dB Phase degree o o 10 107 10 Frequency Hz Figure 7 4 Controller transfer function Figure 7 5 a 7 5 b show that the trajectory is tracked good The output is within the spec immediately at the constant part of the trajec tory 1 15 2 t s a Time respons of output and trajectory
89. sultaten zijn behaald maar er is nog veel onderzoek wat nog uitgevoerd moet worden Denk hierbij b v aan niet lineariteiten discretisatie en MIMO systemen Dit afstudeeronderzoek is een eerste aanzet om te onderzoeken of GP geschikt is om regelaars te genereren voor een willekeurig systeem Een futuristisch idealistisch doel is om een programma te hebben dat een optimale regelaar cre ert voor een willekeurig systeem Alvorens dit doel te bereiken is de grootste uitdaging het defini ren van de regelaar eisen in een GP fitness functie Daarna zal de evolutie van GP de rest doen De enige benodigde vaardigheid zou dan op de startknop drukken zijn Structural and Parametric Optimisation of non linear Dynamical Systems using Genetic Programming Dennis Bruijnen Eindhoven University of Technology TU e Department of Mechanical Engineering Control Systems Technology Group 1 Introduction The design of a controller depends on the dy namics of the plant and the desired require ments They are designed by control engineers who mostly implement a PID like SISO con troller It is a simple and understandable control scheme easily tuned and implemented How ever incase of highly coupled and non linear system behavior this SISO PID approach is not always possible Since digital controllers are available and cheap a more complex controller can be im plemented with only a little bit of more effort However the development of such a c
90. t ness is often accepted which causes the state to travel from one optimum to another Decreas ing the temperature this will become less until it stagnates in an optimum This is only a brief explanation of the method in literature there is more information available 23 How to alter the state is not defined by Sim ulated Annealing Often a combination of the Nelder Mead Simplex method with Simulated Annealing is used It has shown to be a very effective method to optimize a multidimensional nonlinear function Genetic Algorithm A Genetic Algorithm GA is highly related to GP The main difference is the representation of the solution GA only optimizes a fixed length representation of parameters Because of this it is much more efficient in optimizing parameters than GP A lot of variants have been developed since the 60s when GA was developed A description of such a variant of GA can be found in 4 and will be described briefly here Like GP GA also generates a population of in 17 5 4 Choice of method dividuals which is modified every generation us ing genetic operations An individual is simply a collection of a fixed amount of numbers Its fitness is determined by a fitness function Each generation the individuals with the best fitness are chosen to reproduce These parents produce offspring using the crossover operation Then the values of each offspring are mutated with a random number using a normal distri
91. t generation The result is that GP probably stagnates in a local optimum The diversity of the population decreases so the chance that the right structure with the correct parameters will be found be comes smaller A solution to this problem is to use a dedi cated parameter optimization routine Optimiz ing a dynamic model according to a calculated fitness function is actually optimizing a multidi mensional nonlinear function In literature a lot of methods are presented Mainly they can be divided into two parts deterministic methods and stochastic methods There is no best of all method Which method is the best to use is problem dependent 5 MULTIDIMENSIONAL OPTIMIZATION OF A NONLINEAR FUNCTION 5 2 Deterministic methods Local gradient Deterministic optimization methods can be di vided into two types methods which use only function evaluations for example the Nelder Mead Simplex method and methods which use also derivative information for example the Conjugate Gradient method These two exam ples are explained briefly in this section More information can be found in 23 1D optimization All gradient methods are based on optimizing a function in one direction and then choose an other direction and so on Two commonly used methods are Parabolic Interpolation and Golden Section Search The first one chooses three points and fits a parabola through it mum of this parabola is added to the points and the outer point i
92. that this block is not proper and will cause problems when solving the differential equation numerically The gain goes to infinity with increasing frequency Instead of this a proper transfer function with a pole is used Below that pole frequency it is equal to a pure differentiator Gain S y t cu t c loguniform 4 n The parameter of a Gain block is initialized using the loguniform probability distribution in order to keep the chance of a gain larger than 1 equal to the chance of a gain smaller than 1 This is done to avoid the assumption that a gain is mostly used for amplification TF 1st order e1 co s c0 y s circio co uniform 0 n c loguniform 1 n 8 00 1 TFI y s cos u s Co C1 uniform 0 n ca loguniform 2 n The TF blocks are scaled in such a way that the static gain is 1 An extra gain with a loguniform probability distribution is added to give the transfer function more freedom so it can have any shape possible for such a transfer function Further TF 1st order is one block with two possibilities At the initialization of the block one of the two possibilities will be chosen randomly This is also done for TF 2nd order which has three possibilities The reason of this construction is that only proper transfer functions can be made with Simulink Without such a construction it is not possible to place real and complex zeros and poles arbitrarily TF 2nd order c2 c1 c1 L s
93. tions GP Techniques and Theory Artificial Life Autonomous Agents Control Financial Trading Neural Networks Art Image amp Signal Processing Prediction amp Classification Optimisation Identification Figure 3 1 Overview of GP research 3 GP RESEARCH 3 1 System identification with GP Structural system identification using genetic programming and a block diagram oriented sim ulation tool 6 In this paper nonlinear system identification using GP is investigated Models are created in Matlab Simulink using continuous time blocks s domain A plant is created containing a sec ond order transfer function a saturation and a time delay Time domain input output data is created by simulating the model The chosen input involved steps and ramps The same in put data is given to generated models by GP The fitness function is defined as the sum of the squared error between the output of the plant and the GP model A combination of simulated annealing and nelder simplex optimization was used to support GP with the parameter opti mization This research is very similar to the sys tem identification approach in this thesis Mat lab Simulink is used for simulations and com parison of models in the time domain The ma jor difference is a different parameter optimiza tion routine Further only linear system identi fication has been carried in this thesis because our main objective was to perform controller synth
94. tituting branches or replacing nodes Maintaining the tree structure will prevent syntax errors An example of a tree structure is shown in Figure 2 1 y P m root node C Q function terminal Figure 2 1 An example of a tree structure y ue u 2 A tree consists of nodes connected with each other The root node of the tree returns the out put of the computer program A node can have an arbitrary number of inputs If a node does not have any inputs a so called terminal it will stop that branch from growing Starting with the root node randomly selected nodes can be attached until all branches have stopped grow ing because of terminal nodes The tree depth is defined as the number of layers in a tree e g the tree in Figure 2 1 has a tree depth of 4 Now a genetic operation can be applied These are described in Section 2 3 5 More detailed infor mation about genetic operations can be found in 20 A tree structure is not the only way to rep resent a function In literature a lot of differ ent architectures have been proposed Roughly the main categories are tree structures linear structures and graph network structures see Section 5 2 of 20 2 GENETIC PROGRAMMING 2 2 2 Nodes There are two types of nodes functions and ter minals Examples of functions are boolean func tions and or not arithmetic functions X transfer functions Ede and many more Examples of termi
95. tness function was not setup correctly The fitness function should always return a number which must be between 0 and 101 BadSimPenalty in every possible situation For example simulation errors which are realistic to occur because of the randomly created models can be caught by a try catch construction If the performance of the model is bad after running the GP program for a long time check the elements of the fitness function and the settings in gp ini The fitness function shapes the search space and the settings of GP determine how to explore the search space Both items are crucial for convergence to a good solution D 2 C code description The C application gp4so Genetic Programming For System Optimization is based on the EO library Evolving Objects Library for C version 0 9 2 This is an object oriented and template based library which is developed for implementing Evolutionary Algorithms It provides a lot of classes that can be of use It is not a ready made application so making it dedicated for a specific task it needs to be extended with self made classes To get familiar with using the EO library has proven to be quite difficult and will take some time D 2 1 C files Besides of the EO library files some additional files are needed to create the gp4so executable gp4so cpp engine cpp engine h inireader cpp inireader h logdata cpp logdata h matlabcom cpp matlabcom h node cpp node h
96. to use is discussed in the next chapter 15 5 Multidimensional optimiza tion of a nonlinear function 5 1 A need for numerical optimiza tion During the first tests with the GP program some problems occurred GP has been implemented using a tree representation for an individual see Chapter 4 The problem is that transfer func tions can not be made in a tree structure us ing elementary blocks The solution is to use first order and second order transfer function blocks so poles and zeros can be placed arbi trarily by the GP program However GP is not able to optimize the parameters of those blocks with the genetic operations GP only mutates whole nodes blocks A block like a second order transfer function can contain up to 5 parameters which is far from elementary The chance that it mutates to the same block with slightly dif ferent parameters is very small so optimization of these blocks is almost impossible with the GP program For simple models it quickly finds the right structure but it takes a long time before it finds the optimal parameters GP then keeps on try ing different models because it does not know that the right structure has already been found For more complex models GP occasionally produces the right structure during the run but the parameters are not yet optimized As the chance is very small that GP also finds the op timal parameters the model gets a bad fitness and it will not survive in the nex
97. tructural system identification using genetic programming and a block diagram oriented simulation tool 1996 7 J Dabney T Harman Mastering Simulink 4 Prentice Hall inc 2001 pum 8 G Naus System Identification using GP 2003 9 http www revolutionaryengineering com EA 1 html 10 11 http eodev sourceforge net http www geneticprogramming com um 3 2 http www genetic programming org http www ieee org mm 4 http asd bbn com papers traffic traffic html KZ an e AS ro ee pa 15 I Soute Genetic Programming TU Eindhoven mechanical engineering Dynamic Systems Design report number 2000 25 september 2000 16 J Koza M Keane F Bennett J Yu W Mydlowec O Stiffelman Automatic Creation of both the Topology and Parameters for a Robust Controller by Means of Genetic Programming 1999 17 J Koza Genetic Programming On the Programming of Computers by Natural Selection MIT press Cambridge MA 1992 18 M Keijzer J J Merelo G Romero M Schoenauer Evolving Objects a general purpose evo lutionary computation library 2001 19 M Keijzer EO Reference Manual 0 9 2 2002 20 W Banzhaf P Nordin R Keller and F Francone Genetic Programming an Introduction on the Automatic Evolution of Computer Programs and its Applications Morgan Kaufman Publishers Inc 1998 21 W Langdon A Qureshi Genetic Programming Computers using Natural Sel
98. try to en sure that the directions that have already been minimized stay minimized Contrary to intu ition this does not mean following the line of steepest descent each time The conjugate directions are actually calcu lated on the assumption that the error surface is quadratic which is not generally the case How ever it is a fair working assumption and if the algorithm discovers that the current line search direction is not actually downhill it simply cal culates the line of steepest descent and restarts the search in that direction Once a point close to a minimum is found the quadratic assump tion holds true and the minimum can be located very quickly Nelder Mead Simplex method The Nelder Mead Simplex method is an opti mization method that does not need derivatives It does not make use of line optimization The main concept used by the Nelder Mead Simplex method is the geometrical concept of a simplex A simplex is the geometrical figure consisting in n dimensions of n 1 points or vertices and all their interconnecting line segments polygo nal faces etc Examples are a triangle n 2 or a tetrahedron n 3 What the Nelder Mead Simplex algorithm ba sically does is to calculate the function at each of these vertices and order them according to the function value The worst is discarded and a new vertice is proposed using rules of reflec tion expansion and contraction The other ver tices are retained a

(non) linear Dynamical Systems using Genetic Programming

Contents

Download Pdf Manuals

Related Search

Related Contents