**Part 5**

**Vehicle Applications** 

406 MATLAB for Engineers – Applications in Control, Electrical Engineering, IT and Robotics

John Chiasson& Marc Bodson (1998) .Differential – Geometric methods for control of lectric motors*. International Journal of Robust and Nonlinear Control.* No.8, pp.923-954. Arnold V.I.( 1983) *Geometric methods in the theory of ordinary differential equations.* Springer-

Kang .W & Krener.A.J. (1992) . Extended quadratic controller normal form and dynamic

Devanathan R. (2001) . Linearisation Condition through State Feedback. *IEEE transactions on* 

Devanathan R. (2004). Necessary and sufficient conditions for quadratic linearisation of a linearly controllable system. *INT.J. Control*, Vol.77, No.7,pp. 613-621

A.K.Parvathy, Aruna Rajan, R.Devanathan (2005) . Complete Quadratic Linearization of PM

A.K.Parvathy, R.Devanathan (2006) Linearisation of Permanent Magnet Synchronous Motor

Pavol Brunovsky (1970) *A classification of linear controllable systems.* Kybernetika,Vol. 6. No.3,

Asriel U Levin & Kumpati S Narendra (1993). Control of nonlinear dynamical systems using

A.K.Parvathy, V.Kamaraj & R.Devanathan (2011) A generalized quadratic linearization

Synchronous Motor Model. , *proceedings of NEPC conference.* IIT Karagpur, pp 49-

neural networks: controllability and stabilization. *IEEE Transactions on neural* 

technique for PMSM. accepted for publication in , *European Journal of Scientific* 

state feedback linearisation of nonlinear systems. *SIAM Journal of Control and* 

Verlag, NewYork, pp 177-188.

52.

pp. 173-188

*Research.*

*optimization* . No.30,pp.1319-1337.

*Automatic Control* .ol 46,no.8,pp.1257-1260

Model.*proc.of IEEE ICIT 2006*,pp. 483-486.

*networks.*Vol-4, No.2,1993, pp.192-206.

Benjamin.C.Kuo(2001) . *Automatic Control Systems.* Prentice-Hall India.

**17** 

*Brazil* 

Anibal Azevedo

*State University of São Paulo* 

**Automatic Guided Vehicle Simulation** 

**in MATLAB by Using Genetic Algorithm** 

The prodigious advances in robotics in recent times made the use of robots more present in modern society. One important advance that requires special attention is the development of an unmanned aerial vehicle (UAV), which allows an aircraft to fly without having a human crew on board, although the UAVs still need to be controlled by a pilot or a

Today's UAVs often combine remote control and computerized automation in a manner that built-in control and/or guidance systems perform deeds like speed control and flightpath stabilization. In this sense, existing UAVs are not truly autonomous, mostly because air-vehicle autonomy field is a recently emerging field, and this could be a bottleneck for

It could be said that the ultimate goal in the autonomy technology development is to replace human pilots by altering machines decisions in order to make decisions like humans do. For this purpose, several tools related with artificial intelligence could be employed such as expert systems, neural networks, machine learning and natural language processing (HAYKIN, 2009). Nowadays, the field of autonomy has mostly been following a bottom-up

One interesting methodology from the hierarchical control systems approach is the subsumption architecture that decomposes complicated intelligent behavior into many "simple" behavior modules, which are organized into layers. Each layer implements a particular goal and higher layers are increasingly more abstract. The decisions are not taken by a superior layer, but by listening to the information that are triggered by sensory inputs (lowest layer). This methodology allows the use of **reinforcement learning** to acquire

Inspired by old behaviorist psychology, **reinforcement learning** (RL) concerned with how an *agent* ought to take *actions* in an *environment,* so as to maximize some notion of cumulative *reward*. Reinforcement learning differs from standard supervised learning in that correct input/output pairs which are never presented (HAYKIN, 2009). Furthermore, there is a focus on an on-line performance, which involves finding a balance between exploration (of uncharted territory) and exploitation (of current knowledge). The reinforcement learning has been applied successfully to various problems, including robot control, elevator scheduling, telecommunications, backgammon and chess (SHIM, 2000).

approach, such as hierarchical control systems (SHIM, 2000).

behavior with the information that comes with experience.

**1. Introduction** 

future UAV development.

navigator.

### **Automatic Guided Vehicle Simulation in MATLAB by Using Genetic Algorithm**

Anibal Azevedo

*State University of São Paulo Brazil* 

#### **1. Introduction**

The prodigious advances in robotics in recent times made the use of robots more present in modern society. One important advance that requires special attention is the development of an unmanned aerial vehicle (UAV), which allows an aircraft to fly without having a human crew on board, although the UAVs still need to be controlled by a pilot or a navigator.

Today's UAVs often combine remote control and computerized automation in a manner that built-in control and/or guidance systems perform deeds like speed control and flightpath stabilization. In this sense, existing UAVs are not truly autonomous, mostly because air-vehicle autonomy field is a recently emerging field, and this could be a bottleneck for future UAV development.

It could be said that the ultimate goal in the autonomy technology development is to replace human pilots by altering machines decisions in order to make decisions like humans do. For this purpose, several tools related with artificial intelligence could be employed such as expert systems, neural networks, machine learning and natural language processing (HAYKIN, 2009). Nowadays, the field of autonomy has mostly been following a bottom-up approach, such as hierarchical control systems (SHIM, 2000).

One interesting methodology from the hierarchical control systems approach is the subsumption architecture that decomposes complicated intelligent behavior into many "simple" behavior modules, which are organized into layers. Each layer implements a particular goal and higher layers are increasingly more abstract. The decisions are not taken by a superior layer, but by listening to the information that are triggered by sensory inputs (lowest layer). This methodology allows the use of **reinforcement learning** to acquire behavior with the information that comes with experience.

Inspired by old behaviorist psychology, **reinforcement learning** (RL) concerned with how an *agent* ought to take *actions* in an *environment,* so as to maximize some notion of cumulative *reward*. Reinforcement learning differs from standard supervised learning in that correct input/output pairs which are never presented (HAYKIN, 2009). Furthermore, there is a focus on an on-line performance, which involves finding a balance between exploration (of uncharted territory) and exploitation (of current knowledge). The reinforcement learning has been applied successfully to various problems, including robot control, elevator scheduling, telecommunications, backgammon and chess (SHIM, 2000).

Automatic Guided Vehicle Simulation in MATLAB by Using Genetic Algorithm 411

Fig. 1. Description of the elements of the developed environment in Matlab.

Vector

Fig. 2. The correspondence between the obstacles detected by a radar vehicle and a vector

Fig. 3. The correspondence between a rule and decision taken by the vehicle.

with binary information.

Situation 1 1 0 0 0 0 0 0 0

Index 0 1 2 3 4 5 6 7 8

Sector 1 2 3 4 5 6 7 8 9

**Genetic algorithms** (GAs) are developed in order to emulate the process of genetic evolution found in nature, and then perform artificial evolution. They were developed by John Holland [11] in the early 70s, and have been successfully applied to numerous large and complex search space problems ever since (MICHALEWICZ, 1996).

In nature, organisms have certain characteristics that affect their ability to survive and reproduce. These characteristics are contained in their genes. Natural selection ensures that genes from a strong individual are presented in greater numbers in the next generation, rather than those from a weak individual. Over a number of generations, the fittest individuals, in the environment in which they live, have the highest probability of survival and tend to increase in number; while the less fit individuals tend to die out. This is the Darwin's principle of the survival-of-the-fittest and constitutes the basic idea behind the GAs. In order to perform computational tests on how the reinforcement learning could cope with genetic algorithms to provide good rules for the navigation of an automatic vehicle (STAFYLOPATIS, 1998), a program that emulates a navigation environment was developed in a Matlab. This Chapter will describe how the ideas developed by (STAFYLOPATIS, 1998) could be employed to study the integration of the GAs and RL to produce rules for automatic guided vehicles searching for a better performance. The main contributions of this Chapter are: the vehicle, its sensors, and also the environment for training are different from the ones presented in (STAFYLOPATIS, 1998); a new equation for the reinforcement learning was proposed; the influence of the parameters that control the production of automatic generation of rules for vehicle control navigation are also tested. Sections 2 and 3 describe how the navigation decision rules could be encoded in a vector. Sections 4 and 5 show how the reinforcement learning and genetic algorithm uses the encoding of Section 2 to cope with the production of vehicle control navigation rules. In Section 6, some results are presented and finally, in Section 7, some conclusions and future works are given. All the proposed approach had been implemented in a Matlab.

#### **2. Navigation problem representation**

The navigation problem could be defined as how a vehicle could be guided through an ambient with many obstacles and barriers using the information available from the information given by the vehicle radar, as it can be seen in Figure 1.

Figure 1 also details how the vehicle radar works. The radar has 9 sectors in order to better acquire the information of how near is a vehicle to an object. The information of existence of an obstacle for each section is stored in vector **v** of nine positions using the following rule: if one object is in one section, then value 1 is attributed. If not, then value 0 is attributed. Figure 2 details a situation where there are obstacles in sector 1 and 2 and the correspondent representation by a vector. Once the objects had been detected by the radar, it is necessary to implement an appropriate action which could be one of the following three: turn 15o degrees to the left, keep the trajectory or turn 15o degrees to the right. For the sake of clarity, but without loss of generality, velocity vc of the vehicle was assumed to be constant. With these three possible decisions, a new concept could be created, which is a rule. A **rule** is a vector that combines the one that describes the situation for the vehicle in terms of obstacles, like the one presented in Figure 2, and a new component that decides the vehicle action in turn to avoid a set of obstacles: 0 – turn left, 1 – keep the trajectory, 2 – turn right. An example of a rule for the situation showed in Figure 2 is described in Figure 3.

**Genetic algorithms** (GAs) are developed in order to emulate the process of genetic evolution found in nature, and then perform artificial evolution. They were developed by John Holland [11] in the early 70s, and have been successfully applied to numerous large

In nature, organisms have certain characteristics that affect their ability to survive and reproduce. These characteristics are contained in their genes. Natural selection ensures that genes from a strong individual are presented in greater numbers in the next generation, rather than those from a weak individual. Over a number of generations, the fittest individuals, in the environment in which they live, have the highest probability of survival and tend to increase in number; while the less fit individuals tend to die out. This is the Darwin's principle of the survival-of-the-fittest and constitutes the basic idea behind the GAs. In order to perform computational tests on how the reinforcement learning could cope with genetic algorithms to provide good rules for the navigation of an automatic vehicle (STAFYLOPATIS, 1998), a program that emulates a navigation environment was developed in a Matlab. This Chapter will describe how the ideas developed by (STAFYLOPATIS, 1998) could be employed to study the integration of the GAs and RL to produce rules for automatic guided vehicles searching for a better performance. The main contributions of this Chapter are: the vehicle, its sensors, and also the environment for training are different from the ones presented in (STAFYLOPATIS, 1998); a new equation for the reinforcement learning was proposed; the influence of the parameters that control the production of automatic generation of rules for vehicle control navigation are also tested. Sections 2 and 3 describe how the navigation decision rules could be encoded in a vector. Sections 4 and 5 show how the reinforcement learning and genetic algorithm uses the encoding of Section 2 to cope with the production of vehicle control navigation rules. In Section 6, some results are presented and finally, in Section 7, some conclusions and future works are given. All the proposed approach had been

The navigation problem could be defined as how a vehicle could be guided through an ambient with many obstacles and barriers using the information available from the

Figure 1 also details how the vehicle radar works. The radar has 9 sectors in order to better acquire the information of how near is a vehicle to an object. The information of existence of an obstacle for each section is stored in vector **v** of nine positions using the following rule: if one object is in one section, then value 1 is attributed. If not, then value 0 is attributed. Figure 2 details a situation where there are obstacles in sector 1 and 2 and the correspondent representation by a vector. Once the objects had been detected by the radar, it is necessary to implement an appropriate action which could be one of the following three: turn 15o degrees to the left, keep the trajectory or turn 15o degrees to the right. For the sake of clarity, but without loss of generality, velocity vc of the vehicle was assumed to be constant. With these three possible decisions, a new concept could be created, which is a rule. A **rule** is a vector that combines the one that describes the situation for the vehicle in terms of obstacles, like the one presented in Figure 2, and a new component that decides the vehicle action in turn to avoid a set of obstacles: 0 – turn left, 1 – keep the trajectory, 2 – turn right. An example of

and complex search space problems ever since (MICHALEWICZ, 1996).

implemented in a Matlab.

**2. Navigation problem representation** 

information given by the vehicle radar, as it can be seen in Figure 1.

a rule for the situation showed in Figure 2 is described in Figure 3.

Fig. 1. Description of the elements of the developed environment in Matlab.

Fig. 2. The correspondence between the obstacles detected by a radar vehicle and a vector with binary information.

Fig. 3. The correspondence between a rule and decision taken by the vehicle.

Automatic Guided Vehicle Simulation in MATLAB by Using Genetic Algorithm 413

(A) The rule that match the acquired information obtained by the vehicle radar.

(B) Employing the rule – Step 1, Situation 1. (C) Employing the rule – Step 2, Situation 1.

(D) Employing the rule – Step 1, Situation 2. (E) Employing the rule – Step 2, Situation 2. Fig. 5. Employing the same rule in different scenarios and the corresponding consequences. Where: *Si* is the measure of how specific is rule *i*, *n* is the number of sectors the vehicle radar has, *ki* is the number of digits marked with value 2, and *R* is the number of rules. Thus *Si*

*[0,1]* and note that variable *Si* becomes 0 when every rule digit related with the sector state is marked with value 2 (lowest specificity) and *Si* becomes 1, if value 2 does not appear in rule

This discussion also shows that the efficiently vehicle navigation control depends on more than one rule, and this leads to the objective of finding the set of rules that could properly

*i* (highest specificity).

∈

#### **3. Navigation control**

The previous binary representation scheme has a main drawback since before every vehicle movement decision demanded to store 3×29 = 1536 vectors (rules) to precisely describe which action should be performed for each scenario detected by the radar information. This implies that before every step performed by the vehicle, a control action should compare, in the worst case, 1536 vectors in order to find the proper decision to be taken. This scheme is illustrated in Fig. 4.

Fig. 4. Decision scheme followed by the automatic guided vehicle.

This computational work could be avoided by introducing a new symbol 2, for the radar information section, which means that "it could have or not an obstacle at this sector". The advantage of the definition of this new symbol is that it will help to reduce the number of necessary information to be kept by the vehicle control. The disadvantage of using this new symbol is that it could group different situations where the decisions should not be the same. Fig. 5 gives an example of this situation.

The example showed in Fig. 5 emphasizes the importance to construct the rules in a manner that the parts of the rule which do not affect the decision should be numbered as 2, and the other parts that have a great influence in the final behavior of the vehicle should be numbered as most specific as possible or, in other words, with the numbers 0 or 1. Fig. 6 illustrates the application of this concept.

Fig. 6 also gives a guideline procedure for another problem that could emerge, which is the appearance of more than one rule that matches the current radar information when symbol 2 is used. One criteria that will be adopted is to select the rules which match the environment situation, but with as much specific information as possible. For this purpose Eq. (1) will be adopted.

$$S\_i = \frac{(n - k\_i)}{n} \qquad i = 1, \cdots, R \tag{1}$$

The previous binary representation scheme has a main drawback since before every vehicle movement decision demanded to store 3×29 = 1536 vectors (rules) to precisely describe which action should be performed for each scenario detected by the radar information. This implies that before every step performed by the vehicle, a control action should compare, in the worst case, 1536 vectors in order to find the proper decision to be taken. This scheme is

Fig. 4. Decision scheme followed by the automatic guided vehicle.

same. Fig. 5 gives an example of this situation.

illustrates the application of this concept.

Eq. (1) will be adopted.

This computational work could be avoided by introducing a new symbol 2, for the radar information section, which means that "it could have or not an obstacle at this sector". The advantage of the definition of this new symbol is that it will help to reduce the number of necessary information to be kept by the vehicle control. The disadvantage of using this new symbol is that it could group different situations where the decisions should not be the

The example showed in Fig. 5 emphasizes the importance to construct the rules in a manner that the parts of the rule which do not affect the decision should be numbered as 2, and the other parts that have a great influence in the final behavior of the vehicle should be numbered as most specific as possible or, in other words, with the numbers 0 or 1. Fig. 6

Fig. 6 also gives a guideline procedure for another problem that could emerge, which is the appearance of more than one rule that matches the current radar information when symbol 2 is used. One criteria that will be adopted is to select the rules which match the environment situation, but with as much specific information as possible. For this purpose

<sup>−</sup> <sup>=</sup> *i R* <sup>=</sup>1, , (1)

( )*<sup>i</sup>*

*i n k <sup>S</sup> n*

**3. Navigation control** 

illustrated in Fig. 4.

(A) The rule that match the acquired information obtained by the vehicle radar.

(B) Employing the rule – Step 1, Situation 1. (C) Employing the rule – Step 2, Situation 1.

(D) Employing the rule – Step 1, Situation 2. (E) Employing the rule – Step 2, Situation 2.

Fig. 5. Employing the same rule in different scenarios and the corresponding consequences.

Where: *Si* is the measure of how specific is rule *i*, *n* is the number of sectors the vehicle radar has, *ki* is the number of digits marked with value 2, and *R* is the number of rules. Thus *Si* ∈ *[0,1]* and note that variable *Si* becomes 0 when every rule digit related with the sector state is marked with value 2 (lowest specificity) and *Si* becomes 1, if value 2 does not appear in rule *i* (highest specificity).

This discussion also shows that the efficiently vehicle navigation control depends on more than one rule, and this leads to the objective of finding the set of rules that could properly

Automatic Guided Vehicle Simulation in MATLAB by Using Genetic Algorithm 415

rules where used, and to punish just the rules mostly related to the collision, only the last

Eq. (3) was considered by (STAFYLOPATIS, 1998), but it represents that the rules that lead the vehicle to perform more than 200 steps without a collision will suffer the same punishment as the rules that conduct the vehicle to perform only 20 steps without a collision, for example. To avoid this problem, this book Chapter proposes, for the first time, a punishment formula

> 1 1 1 ( )/ exp *<sup>p</sup> ns np <sup>k</sup>* − + <sup>=</sup> <sup>+</sup>

<sup>1</sup> /2 *C C it it* <sup>+</sup> = (3)

*C kC it p it* <sup>+</sup><sup>1</sup> <sup>=</sup> (4)

(5)

three employed rules will have their credit updated by Eq. (3).

which is a function on the number of steps as showed in Eq. (4) and (5).

The behavior of Eq. (5) for *np = 20* is showed in Fig. 7.

Fig. 7. Behavior of the parameter penalty factor *kp* for *np = 20*.

The purpose of Eq. (4) and (5) is to punish the rules according with total number of steps (*ns*) performed by the vehicle before the collision, and it was reflected in a different credit actualization. For example, if the three last used rules lead to only 25 steps without a collision, then the new credit value will be 75% of the old credit value (see Fig. 7). The system proposed by Eq. (2) and (3) or (4) and (5) starts by assigning the same value, same credit, for every rule. Then, the rules performance is evaluated in an environment, created in Matlab, where the vehicle control is simulated, and during the simulation, after each step performed by the vehicle, it is verified whether there is a collision (then the Eq. (3) or Eq. (4) and (5) is used) or not (then the Eq. (2) is used) following the scheme described in Fig. 8.

control the vehicle. To achieve this objective the next sections will discuss concepts on how to automatically construct this set of rules by employing computational tools such as Reinforcement Learning and Genetic Algorithm.

Fig. 6. Employing a rule that combines general (number 2) with specific (number 0 or 1) knowledge.

#### **4. Reinforcement learning**

The developed program starts with rules totally generated at random and which quality, in terms of helping the vehicle navigation control, is unknown. Then, it was created a system that evaluates the quality of rules used through the vehicle navigation performed by a computational simulation. As it was seen before, during the navigation, more than one rule could be necessary and the system must punish or reward not only a rule, but also all rules employed that lead the vehicle to a collision or to avoid an obstacle, respectively. One difficulty related with this scheme is to know the correct contribution weight of each rule in a possible vehicle collision. It is also difficult to estimate the reward value for a set of rules that helps the vehicle to avoid an immediate collision, because this set could put the vehicle in a trajectory that leads to a future collision with another obstacle. Besides those problems, the initial credit received for each rule is modified by Eq. (2), if there are no collisions on the current vehicle step *t*.

$$\mathbf{C}\_{it+1} = \mathbf{C}\_{it} + kS\_{it}\mathbf{C}\_{it} \tag{2}$$

Where: *i* is the rule that perfectly matches with the current vehicle situation; in order words, rule *i* is the one which *Sit* value is as greater as possible at the step *t* of the vehicle, *Cit* is the credit value assigned for each rule *i* at the step *t,* and *k* is a constant that indicates the learning rate.

When the vehicle crashes at some obstacle a punishment must be applied to the rules mostly related to the event. To avoid the storage of large amounts of historical record about what

control the vehicle. To achieve this objective the next sections will discuss concepts on how to automatically construct this set of rules by employing computational tools such as

(A) The rule that matches the acquired information obtained by the vehicle radar.

(B) Employing the rule – Step 1. (C) Employing the rule – Step 2.

The developed program starts with rules totally generated at random and which quality, in terms of helping the vehicle navigation control, is unknown. Then, it was created a system that evaluates the quality of rules used through the vehicle navigation performed by a computational simulation. As it was seen before, during the navigation, more than one rule could be necessary and the system must punish or reward not only a rule, but also all rules employed that lead the vehicle to a collision or to avoid an obstacle, respectively. One difficulty related with this scheme is to know the correct contribution weight of each rule in a possible vehicle collision. It is also difficult to estimate the reward value for a set of rules that helps the vehicle to avoid an immediate collision, because this set could put the vehicle in a trajectory that leads to a future collision with another obstacle. Besides those problems, the initial credit received for each rule is modified by Eq. (2), if there are no collisions on the

Where: *i* is the rule that perfectly matches with the current vehicle situation; in order words, rule *i* is the one which *Sit* value is as greater as possible at the step *t* of the vehicle, *Cit* is the credit value assigned for each rule *i* at the step *t,* and *k* is a constant that indicates the

When the vehicle crashes at some obstacle a punishment must be applied to the rules mostly related to the event. To avoid the storage of large amounts of historical record about what

*C C kS C it it it it* <sup>+</sup><sup>1</sup> = + (2)

Fig. 6. Employing a rule that combines general (number 2) with specific (number 0 or 1)

Reinforcement Learning and Genetic Algorithm.

knowledge.

**4. Reinforcement learning** 

current vehicle step *t*.

learning rate.

rules where used, and to punish just the rules mostly related to the collision, only the last three employed rules will have their credit updated by Eq. (3).

$$\mathbf{C}\_{it+1} = \mathbf{C}\_{it} \;/\; \mathbf{2} \tag{3}$$

Eq. (3) was considered by (STAFYLOPATIS, 1998), but it represents that the rules that lead the vehicle to perform more than 200 steps without a collision will suffer the same punishment as the rules that conduct the vehicle to perform only 20 steps without a collision, for example. To avoid this problem, this book Chapter proposes, for the first time, a punishment formula which is a function on the number of steps as showed in Eq. (4) and (5).

$$\mathbf{C}\_{it+1} = k\_p \mathbf{C}\_{it} \tag{4}$$

$$k\_p = \frac{1}{1 + \exp^{-(ns+1)/np}}\tag{5}$$

The behavior of Eq. (5) for *np = 20* is showed in Fig. 7.

Fig. 7. Behavior of the parameter penalty factor *kp* for *np = 20*.

The purpose of Eq. (4) and (5) is to punish the rules according with total number of steps (*ns*) performed by the vehicle before the collision, and it was reflected in a different credit actualization. For example, if the three last used rules lead to only 25 steps without a collision, then the new credit value will be 75% of the old credit value (see Fig. 7). The system proposed by Eq. (2) and (3) or (4) and (5) starts by assigning the same value, same credit, for every rule. Then, the rules performance is evaluated in an environment, created in Matlab, where the vehicle control is simulated, and during the simulation, after each step performed by the vehicle, it is verified whether there is a collision (then the Eq. (3) or Eq. (4) and (5) is used) or not (then the Eq. (2) is used) following the scheme described in Fig. 8.

Automatic Guided Vehicle Simulation in MATLAB by Using Genetic Algorithm 417

*Vm* - The vector with the index of the rules that matches the information extracted

*MRE* - Matrix with only the rules that match the information extracted from vehicle

*Specific* - This function evaluates all the rules contained in matrix MRE according with the equation Eq. (1) in order to measure how specific is each rule.

*max* - This function determines which specific measure is the biggest (value) and

*Records* - Vector that stores the index of the rules used in a certain vehicle step through

*indexp* - Vector with the index of the last three rules applied before the vehicle

*kp* - Penalization factor applied in the credit of the last three rules when a collision

The system showed in Fig. 8 tries to establish a reward and a punishment scheme among the rules and their impact in the vehicle navigation through equations that monitor the performance, although the random generation of rules can also produce sets of rules which lead to a bad control navigation performance. After a pre-specified number of collisions, the set of rules could be changed by a new randomly generated set of rules formed with the insertion of new rules, the exclusion of rules with a bad performance (low credit value), and the maintenance of rules that help to avoid the obstacles (high credit value). To perform the formation of this new set of rules, a genetic algorithm was coupled to the credit evaluation

The complete description of the Genetic Algorithm developed is described in Section 5.

The genetic algorithm keeps a population of individuals, represented by: *A(t) =* {A1

adopted here, the population is stored in a matrix *A(t)* and each column *Ai*

for each generation (iteration) *t*, and each individual represents a rule to guide the vehicle through the environment and to avoid the obstacles. In the computational implementation

and then, the *fitness*, a measure of how this individual is sucessful for the problem, is calculated. The *fitness* is calculated for the entire population and is based on this new population that combines the most capable individuals that will form generation *t+1*.

*<sup>t</sup>* is evaluated according to the number of vehicle movements without a collision,

t , ..., Ant}

*<sup>t</sup>* represents a rule.

*VRM* - Vector with the rule selection to be used in current environment state. *Simulate* - This function simulates the vehicle through environment using VRM rules. *collision* - Variable that indicates whether a vehicle collides (equals to 1) or not (equals

collision and which will be punished by a credit decreasing.

*E* - Matrix with all data about the current simulation environment state. *vE* - Vector with information extracted from the radar, as show in Figure 2.

*MR* - Matrix with a set of rules. Each matrix line represents one rule.

The functions and symbols used in Fig. 8 are defined as follows: *C* - Vector containing the credit associated to a set of rules.

from vehicle radar (contained in vE).

the corresponding line of MRE (indexsmax).

*indexsmax* - Index of the rule with maximum matching value.

*sm* - Vector with all rules specific measure.

radar.

to 0).

occurs.

**5. Genetic algorithm** 

Each rule *Ai*

the environment.

*ns* - Number of steps without a collision.

scheme described in Fig. 8, as could be seen in Fig.9.

```
function EvaluateRules (MR) 
 % Initial learning rate for each step without a collision. 
 k = 0.2; 
 collision = 0; 
 % Initial credit assignment for all rules. 
 for i=1:R 
 C(i) = 0.2; 
 end 
 % Starts the initial conditions of the environment and the objects that appear in the radar 
 % vehicle. 
 [vE, E]=Starts(); 
 ns = 0; % Number of steps. 
 Records = [ ]; % Stores the index of the rules used.
 % Loop that evaluates the rules inside MR while there is no collision. 
 While (collision == 0) 
 % Determining the rules that matches with the actual environment. 
 [vm]=Match(MR, vE); 
 % Verifying the measure of how specific each matching rule is (use Equation (1)). 
 MRE = MR(vm,:); 
 [sm]=Specific(MRE); 
 [value, indexsmax] = max(sm); 
 % Using the matching rule with maximum Sit value and verifying the consequences. 
 VRM = MRE(indexsmax,:); 
 % Update the vehicle in environment by using the selected rule. It also returns if 
 % a collision happened (collision == 1). 
 [collision, vE, E]= Simulate(VRM, E) ; 
 % There is no collision with the rule application, then update the credit according 
 % with the equation (2). 
 If (collision == 0) 
 % Update the credit of the rule applied. 
 C(vm(indexsmax)) = (1 + k*sm(indexsmax))*C(vm(indexsmax)) 
 % Store the index of what rule was applied in a certain vehicle step. 
 Records = [Records vm(indexsmax)]; 
 ns = ns + 1; 
 % There is a collision after the rule application then update the credit of the last three 
 % rules used according with equation (3) or (4) and (5). 
 else 
 indexp = Records(end-2:end); 
 % Or use: kp = 1/(exp(-(ns+1)/20)+1); C(indexp)=C(indexp)*kp; 
 C(indexp) = C(indexp)/2; 
 ns = 0; 
 end 
 end
```
Fig. 8. The scheme to reward or to punish the rules during a simulation in the environment.

% Starts the initial conditions of the environment and the objects that appear in the radar

**function** EvaluateRules (MR)

% Initial credit assignment for all rules.

ns = 0; % Number of steps.

[value, indexsmax] = max(sm);

 % a collision happened (collision == 1). [collision, vE, E]= Simulate(VRM, E) ;

% Update the credit of the rule applied.

Records = [Records vm(indexsmax)];

indexp = Records(end-2:end);

C(indexp) = C(indexp)/2;

VRM = MRE(indexsmax,:);

 % with the equation (2). **If** (collision == 0)

ns = ns + 1;

**else** 

 ns = 0; **end end**

**While** (collision == 0)

[vm]=Match(MR, vE);

 MRE = MR(vm,:); [sm]=Specific(MRE);

 k = 0.2; collision = 0;

 **for** i=1:R C(i) = 0.2; **end** 

 % vehicle. [vE, E]=Starts();

% Initial learning rate for each step without a collision.

Records = [ ]; % Stores the index of the rules used.

% Loop that evaluates the rules inside MR while there is no collision.

% Determining the rules that matches with the actual environment.

 C(vm(indexsmax)) = (1 + k\*sm(indexsmax))\*C(vm(indexsmax)) % Store the index of what rule was applied in a certain vehicle step.

% Or use: kp = 1/(exp(-(ns+1)/20)+1); C(indexp)=C(indexp)\*kp;

% rules used according with equation (3) or (4) and (5).

% Verifying the measure of how specific each matching rule is (use Equation (1)).

% Using the matching rule with maximum Sit value and verifying the consequences.

% Update the vehicle in environment by using the selected rule. It also returns if

% There is no collision with the rule application, then update the credit according

% There is a collision after the rule application then update the credit of the last three

Fig. 8. The scheme to reward or to punish the rules during a simulation in the environment.


The system showed in Fig. 8 tries to establish a reward and a punishment scheme among the rules and their impact in the vehicle navigation through equations that monitor the performance, although the random generation of rules can also produce sets of rules which lead to a bad control navigation performance. After a pre-specified number of collisions, the set of rules could be changed by a new randomly generated set of rules formed with the insertion of new rules, the exclusion of rules with a bad performance (low credit value), and the maintenance of rules that help to avoid the obstacles (high credit value). To perform the formation of this new set of rules, a genetic algorithm was coupled to the credit evaluation scheme described in Fig. 8, as could be seen in Fig.9.

The complete description of the Genetic Algorithm developed is described in Section 5.

#### **5. Genetic algorithm**

The genetic algorithm keeps a population of individuals, represented by: *A(t) =* {A1 t , ..., Ant} for each generation (iteration) *t*, and each individual represents a rule to guide the vehicle through the environment and to avoid the obstacles. In the computational implementation adopted here, the population is stored in a matrix *A(t)* and each column *Ai t* represents a rule. Each rule *Ai <sup>t</sup>* is evaluated according to the number of vehicle movements without a collision, and then, the *fitness*, a measure of how this individual is sucessful for the problem, is calculated. The *fitness* is calculated for the entire population and is based on this new population that combines the most capable individuals that will form generation *t+1*.

Automatic Guided Vehicle Simulation in MATLAB by Using Genetic Algorithm 419

Each genetic algorithm individual is associated to a rule by using a vector *v* with 10 elements and values inside the elements 1 to 9, corresponding to the presence (value 1) or not (value 0) or whatever (value 2) in the corresponding vehicle radar section, respectively (see Fig. 3). The value inside the position 10 indicates the decision that has to be made: 0 turn 15o degrees to the left, 1 – keep the trajectory and 2 - turn 15o degrees to the right. The set of rules is stored in a matrix and each column represents a rule, as showed in Fig. 10. In the Example showed in Fig. 10, the first column contains a vector **v1** with 10 positions which

Fig. 10. Relation between the genetic algorithm individual encoding and a rule. Matrix A has

Once the individual is defined, the population which consists in *numpop* individuals stored

value of *numpop* equals to 100 was adopted. Since each column of matrix *A* represents individuals/solutions, then every element *A(i,j) = k* means what situation (*i* equals 1 to 9) or decision (*i* equals to 10) for the rule *j* should be made if the solution given from individual *j* is chosen. For example, *A(1,1) = 2* means that the individual/solution 1 vehicle radar does

The fitness evaluation or *fitness* is responsible to rank the individual among the population in a generation t. Then, the fitness is constructed in manner that the solutions with less credit measure have the higher fitness value. This explanation justifies why the fitness for

not care (value 2) whether there is or not an object in sector 1 of the vehicle radar.

 *numpop* can also be defined. For the vehicle navigation, the

() () *Fit Ai i* = *f A* (6)

×

each individual *Ai* is calculated with Eq. (6), given as follows:

The implementation details adopted for the Genetic Algorithm are as follows:

**(A) Data structure codification for each individual:** 

each position meaning explained in Fig. 2.

two rules.

in a matrix A of dimension *10* 

**(B) Fitness evaluation and assignment:** 

Fig. 9. Composition of the credit scheme and Genetic Algorithm.

During the formation of the new population, some individuals from generation *t* are submitted to a transformation process by genetic operators in order to form new rules. These transformations include unary operators *mi* (mutation), that allow the creation of new rules by small changes in the individual attributes (mi: A*<sup>i</sup>* → A*i*), and superior order transformation cj (crossover) that produces new individuals by combining one or more individuals (cj: A*<sup>j</sup>* × ... ×A*<sup>k</sup>* → A*j*). This process is carried out until a previous specified maximum number of generations is reached (MICHALEWICZ, 1996).

During the vehicle navigation along the environment different situations may occur, and each one of them could require a proper action implying the necessity of using more than just one rule. For this reason, the **reinforcement leaning** is acopled with a population produced by the genetic algorithm and, when a collision occurs, a proper number of rules used before the collision are penalized (STAFYLOPATIS, 1998). This system consists in firstly give the same value for each rule of the population, and then, evaluate each rule according to the number of times it was used until a collision happened.

The implementation details adopted for the Genetic Algorithm are as follows:

#### **(A) Data structure codification for each individual:**

418 MATLAB for Engineers – Applications in Control, Electrical Engineering, IT and Robotics

Fig. 9. Composition of the credit scheme and Genetic Algorithm.

maximum number of generations is reached (MICHALEWICZ, 1996).

according to the number of times it was used until a collision happened.

During the formation of the new population, some individuals from generation *t* are submitted to a transformation process by genetic operators in order to form new rules. These transformations include unary operators *mi* (mutation), that allow the creation of new rules by small changes in the individual attributes (mi: A*<sup>i</sup>* → A*i*), and superior order transformation cj (crossover) that produces new individuals by combining one or more individuals (cj: A*<sup>j</sup>* × ... ×A*<sup>k</sup>* → A*j*). This process is carried out until a previous specified

During the vehicle navigation along the environment different situations may occur, and each one of them could require a proper action implying the necessity of using more than just one rule. For this reason, the **reinforcement leaning** is acopled with a population produced by the genetic algorithm and, when a collision occurs, a proper number of rules used before the collision are penalized (STAFYLOPATIS, 1998). This system consists in firstly give the same value for each rule of the population, and then, evaluate each rule Each genetic algorithm individual is associated to a rule by using a vector *v* with 10 elements and values inside the elements 1 to 9, corresponding to the presence (value 1) or not (value 0) or whatever (value 2) in the corresponding vehicle radar section, respectively (see Fig. 3). The value inside the position 10 indicates the decision that has to be made: 0 turn 15o degrees to the left, 1 – keep the trajectory and 2 - turn 15o degrees to the right. The set of rules is stored in a matrix and each column represents a rule, as showed in Fig. 10. In the Example showed in Fig. 10, the first column contains a vector **v1** with 10 positions which each position meaning explained in Fig. 2.


Fig. 10. Relation between the genetic algorithm individual encoding and a rule. Matrix A has two rules.

Once the individual is defined, the population which consists in *numpop* individuals stored in a matrix A of dimension *10* × *numpop* can also be defined. For the vehicle navigation, the value of *numpop* equals to 100 was adopted. Since each column of matrix *A* represents individuals/solutions, then every element *A(i,j) = k* means what situation (*i* equals 1 to 9) or decision (*i* equals to 10) for the rule *j* should be made if the solution given from individual *j* is chosen. For example, *A(1,1) = 2* means that the individual/solution 1 vehicle radar does not care (value 2) whether there is or not an object in sector 1 of the vehicle radar.

#### **(B) Fitness evaluation and assignment:**

The fitness evaluation or *fitness* is responsible to rank the individual among the population in a generation t. Then, the fitness is constructed in manner that the solutions with less credit measure have the higher fitness value. This explanation justifies why the fitness for each individual *Ai* is calculated with Eq. (6), given as follows:

$$Fit(A\_i) = f(A\_i) \tag{6}$$

Where: ( )*<sup>i</sup> f A* is the fitness evaluation as defined by the credit value Cit, and it corresponds to the number of steps performed with rule *i* without a collision, according to the rule stored in vector *Ai*, through the environment.

#### **(C) Individual Selection for the next generation:**

The population formation process employed is the "*Roulette Wheel"*, in order words, a random raffle where the probability *Q(Ai)* to choose the individual *Ai* for the next generation, in a population with *b* individuals, is used. The value *Q(Ai)* can be obtained by using Eq. (7).

$$Q(A\_i) = \text{Fit}(A\_i) \ne \sum\_{i=1}^{numpop} \text{(Fit}(A\_i))\tag{7}$$

Automatic Guided Vehicle Simulation in MATLAB by Using Genetic Algorithm 421

Table 1 sumarizes the description of the parameters necessary to perform tests with the

*pb0* percentance of digits equals to 0 for the initial population. *pb1* percentance of digits equals to 1 for the initial population. *pb2* percentance of digits equals to 2 for the initial population. *numpop* number of individuals in the genetic algorithm population.

*pg* size of the subpopulation in terms of the A matrix (5%).

*N* number of genetic algorithm aplication(14). Table 1. Parameters used for the developed approach and their meanings.

*M* number of collisions before applying the genetic algorithm (10).

It was also tested two environments whose obstacles initial position are shown in Fig. 12. As shown in Fig. 12, the obstacles could be clustered in two groups: the obstacles that delimits the environment bounds have their position fixed and the obstacles that are in the interior of the environment, whose positions are actualized in a different manner for each

• Environment 1: The moving obstacles actualize their x-positions according to Eq. (8).

*t t i i x x* = +

t is the x-position for obstacle *i* in step *t* of the vehicle,

Fig. 12. Detailed description of the obstacles initial position in the environment.

ε

(8)

is a random uniform

ε

**6. Tests and results** 

environment:

Where: xi

variable in the interval [0, 1].

developed approach described in sections 4 and 5.

*Ci0* Initial credit assigment for a set of rules. *k* reinformecent learning rate parameter.

**Parameter Description**

*pc* Crossover rate. *pm* Mutation rate.

The better individual of the actual generation is always kept to be on the next generation.

#### **(D) Crossover OX:**

The Crossover operators try to generate new individuals to form the next population by combining information from past generations that is present in the individuals. This work used two crossover operators which is described as follows.

Two individuals, A1 and A2 , with *N* elements are randomly chosen from a population in generation t. Then, an integer number δ in the interval [1*, N-1*] is drawn and elements of A1 that are in positions δ until *N* are exchanged with the elements of A2 that are in positions δ until *N*. This exchange will produce two new individuals, nA1 and nA2, that can appear in the next generation.

#### **(E) Mutation operator:**

The mutation operator modifies *10\*numpop\*pm* elements of matrix A, where *pm* is the percentage of total bits to be muted (mutated). The selection of which element *Aij* to be mutated consist in randomly select the line index and the column index, and then change it. Fig. 11 shows how the genetic algorithm components described before are combined.

One important observation is that the best individual of generation At is always kept for the next generation A*t+1*. Furthermore, the size of subpopulations *As1 <sup>t</sup>* and *As2 <sup>t</sup>* are the same and is equal to 5% of the total population (*numpop*).

Fig. 11. How the genetic algorithm elements are combined into a unique algorithm.

#### **6. Tests and results**

420 MATLAB for Engineers – Applications in Control, Electrical Engineering, IT and Robotics

Where: ( )*<sup>i</sup> f A* is the fitness evaluation as defined by the credit value Cit, and it corresponds to the number of steps performed with rule *i* without a collision, according to the rule stored

The population formation process employed is the "*Roulette Wheel"*, in order words, a random raffle where the probability *Q(Ai)* to choose the individual *Ai* for the next generation, in a population with *b* individuals, is used. The value *Q(Ai)* can be obtained by

> ( ) ( ) / ( ( )) *numpop ii i i*

*Q A Fit A Fit A*

The better individual of the actual generation is always kept to be on the next generation.

The Crossover operators try to generate new individuals to form the next population by combining information from past generations that is present in the individuals. This work

Two individuals, A1 and A2 , with *N* elements are randomly chosen from a population in

until *N*. This exchange will produce two new individuals, nA1 and nA2, that can appear in

The mutation operator modifies *10\*numpop\*pm* elements of matrix A, where *pm* is the percentage of total bits to be muted (mutated). The selection of which element *Aij* to be mutated consist in randomly select the line index and the column index, and then change it. Fig. 11 shows how the genetic algorithm components described before are combined. One important observation is that the best individual of generation At is always kept for the next generation A*t+1*. Furthermore, the size of subpopulations *As1t* and *As2t* are the same and

Fig. 11. How the genetic algorithm elements are combined into a unique algorithm.

δ

1

<sup>=</sup> (7)

in the interval [1*, N-1*] is drawn and elements of A1

δ

until *N* are exchanged with the elements of A2 that are in positions

=

in vector *Ai*, through the environment.

generation t. Then, an integer number

δ

is equal to 5% of the total population (*numpop*).

using Eq. (7).

**(D) Crossover OX:** 

that are in positions

the next generation. **(E) Mutation operator:** 

**(C) Individual Selection for the next generation:** 

used two crossover operators which is described as follows.

Table 1 sumarizes the description of the parameters necessary to perform tests with the developed approach described in sections 4 and 5.


Table 1. Parameters used for the developed approach and their meanings.

It was also tested two environments whose obstacles initial position are shown in Fig. 12. As shown in Fig. 12, the obstacles could be clustered in two groups: the obstacles that delimits the environment bounds have their position fixed and the obstacles that are in the interior of the environment, whose positions are actualized in a different manner for each environment:

• Environment 1: The moving obstacles actualize their x-positions according to Eq. (8).

$$\mathbf{x}\_{i}^{t} = \mathbf{x}\_{i}^{t} + \mathbf{c} \tag{8}$$

Where: xi t is the x-position for obstacle *i* in step *t* of the vehicle, ε is a random uniform variable in the interval [0, 1].

Fig. 12. Detailed description of the obstacles initial position in the environment.

• Environment 2: The user could also specify a positive value to a *radius* parameter and the moving obstacles actualize their x and y-positions according to Eq. (9).

$$\begin{cases} \mathbf{x}\_i^t = \mathbf{x}\_i^t + radius \ast \cos(\mathbf{x}\_i^t) \\ \mathbf{y}\_i^t = \mathbf{y}\_i^t + radius \ast \sin(\mathbf{y}\_i^t) \end{cases} \tag{9}$$

Automatic Guided Vehicle Simulation in MATLAB by Using Genetic Algorithm 423

The new proposed reinforcement learning formula proposed by Eq. (4) and (5) was also

**Set Test 1 Test 2 Test 3 Test 4 Test 5 Best** 

20 7 15 7 11 20 7 121 7 7 105 121 36 20 27 7 41 41 37 279 339 296 15 339 15 205 100 126 15 205 84 34 31 617 29 617 32 34 16 37 12 37 31 33 165 31 77 165 Table 4. The number of maximum steps without a collision performed by the best set of rules produced by each set of parameters (described in Table 2) and Eq. (4) and (5).

Fig. 13 shows the simulation of the vehicle control navigation with the best set of rules

The data in Table 4 also attest the better performance of the new reinforcement learning scheme proposed by Eq. (4) and (5) (see Table 4) since almost all set of parameters (expections are 3 and 7) created a set of rules with a number of steps without a collision

Fig. 13. The simulation of the vehicle navigation with best rules obtained with set 4.

This Chapter discussed mathematical techniques that could help the ultimate unmanned aerial vehicle (UAVs) objective, which is the development of rules to replace human pilots by altering machines decisions in order to make decisions like humans do. For this purpose,

**Number**

tested for *np = 20* for the parameter values showed in Table 2.

which is bigger than the set of rules created by Eq. (3) (see Table 3).

obtained from the fourth set of parameters.

**7. Conclusions and future works** 

For all tests performed in environment 2, the radius value was fixed in 5. Five runs with the values listed in Table 2 had been performed in order to carry the best set of values for all the parametes listed in Table 1. Each row in Table 2 corresponds to a new configuration.


Table 2. List of values assigned for each parameter.

Each row of Table 3 corresponds to the number of the steps that could be taken with the final set of rules obtained after the application of the scheme described in Fig. 9 using the parameters described in Table 2 and the **reinforcement learning** that follows Eq. (3).


Table 3. The number of maximum steps without a collision performed by the best set of rules produced by each set of parameters (described in Table 2) and Eq. (3).

In Table 3, the first column indicates the set of parameters and the corresponding obtained set of rules which lead to the number of maximum steps without a collision described in the second column until the sixth column. The results for the set of parameters 1 to 6 show that for environment 1 the best choices of parameters are sets 4 and 5. Those two sets were also tested in environment 2 in order to verify set parameters robustness in order to produce adequate rules for different environments. Configuration 7 had the best performance and three tests were sucessfull in producing rules which could make the vehicle able to take more than 80 steps without a collision.

• Environment 2: The user could also specify a positive value to a *radius* parameter and

*t t t i i i t t t i i i*

For all tests performed in environment 2, the radius value was fixed in 5. Five runs with the values listed in Table 2 had been performed in order to carry the best set of values for all the

**Set Ci0 k Pb0 Pb1 Pb2 numpop pc pm pg M N Environment**  20 2 2.5% 2.5% 95% 100 30 1 5 10 14 1 20 2 2.5% 2.5% 95% 100 1 30 5 10 14 1 20 2 2.5% 2.5% 95% 100 30 1 5 20 5 1 20 2 2.5% 2.5% 95% 100 1 30 5 20 5 1 20 2 10% 10% 80% 100 30 1 5 20 5 1 20 2 10% 10% 80% 100 1 30 5 20 5 1 20 2 2.5% 2.5% 95% 100 1 30 5 20 5 2 20 2 10% 10% 80% 100 30 1 5 20 5 2

Each row of Table 3 corresponds to the number of the steps that could be taken with the final set of rules obtained after the application of the scheme described in Fig. 9 using the

**Set Test 1 Test 2 Test 3 Test 4 Test 5 Best** 

16 7 7 7 7 16 7 110 7 7 80 110 28 57 78 16 7 78 177 32 48 21 18 177 151 18 15 32 116 151 34 32 19 32 106 106 32 163 20 106 84 163 38 15 31 108 31 108 Table 3. The number of maximum steps without a collision performed by the best set of

In Table 3, the first column indicates the set of parameters and the corresponding obtained set of rules which lead to the number of maximum steps without a collision described in the second column until the sixth column. The results for the set of parameters 1 to 6 show that for environment 1 the best choices of parameters are sets 4 and 5. Those two sets were also tested in environment 2 in order to verify set parameters robustness in order to produce adequate rules for different environments. Configuration 7 had the best performance and three tests were sucessfull in producing rules which could make the vehicle able to take

parameters described in Table 2 and the **reinforcement learning** that follows Eq. (3).

rules produced by each set of parameters (described in Table 2) and Eq. (3).

*x x radius x y y radius y*

\* cos( ) \* sin( )

(9)

**Number**

the moving obstacles actualize their x and y-positions according to Eq. (9).

parametes listed in Table 1. Each row in Table 2 corresponds to a new configuration.

Table 2. List of values assigned for each parameter.

more than 80 steps without a collision.

= + = +


The new proposed reinforcement learning formula proposed by Eq. (4) and (5) was also tested for *np = 20* for the parameter values showed in Table 2.

Table 4. The number of maximum steps without a collision performed by the best set of rules produced by each set of parameters (described in Table 2) and Eq. (4) and (5).

Fig. 13 shows the simulation of the vehicle control navigation with the best set of rules obtained from the fourth set of parameters.

The data in Table 4 also attest the better performance of the new reinforcement learning scheme proposed by Eq. (4) and (5) (see Table 4) since almost all set of parameters (expections are 3 and 7) created a set of rules with a number of steps without a collision which is bigger than the set of rules created by Eq. (3) (see Table 3).

Fig. 13. The simulation of the vehicle navigation with best rules obtained with set 4.

#### **7. Conclusions and future works**

This Chapter discussed mathematical techniques that could help the ultimate unmanned aerial vehicle (UAVs) objective, which is the development of rules to replace human pilots by altering machines decisions in order to make decisions like humans do. For this purpose,

**18** 

*México* 

**Robust Control of Active Vehicle** 

**Suspension Systems Using Sliding** 

Esteban Chávez Conde1, Francisco Beltrán Carbajal2,

*3Universidad Politécnica de la Zona Metropolitana de Guadalajara* 

*1Universidad del Papaloapan, Campus Loma Bonita* 

*4Instituto Tecnológico de Cd. Guzmán* 

**Modes and Differential Flatness with MATLAB** 

Antonio Valderrábano González3 and Ramón Chávez Bracamontes3

*2Universidad Autónoma Metropolitana, Unidad Azcapotzalco, Departamento de Energía* 

The main control objectives of active vehicle suspension systems are to improve the ride comfort and handling performance of the vehicle by adding degrees of freedom to the system and/or controlling actuator forces depending on feedback and feedforward

Passenger comfort is provided by isolating the passengers from undesirable vibrations induced from irregular road disturbances, and its performance is evaluated by the level of acceleration which vehicle passengers are exposed. Handling performance is achieved by maintaining a good contact between the tire and the road to provide guidance along the track. The topic of active vehicle suspension control system, for linear and nonlinear models, in general, has been quite challenging over the years and we refer the reader to some of the fundamental work in the area which has been helpful in the preparation of this chapter. Control strategies like Linear Quadratic Regulator (LQR) in combination with nonlinear backstepping control techniques are proposed in (Liu et al., 2006). This strategy requires information about the state vector (vertical positions and speeds of the tire and car body). A reduced order controller is proposed in (Yousefi et al., 2006) to decrease the implementation costs without sacrificing the security and the comfort by using accelerometers for measurements of the vertical movement of the tire and car body. In (Tahboub, 2005) a controller of variable gain that considers the nonlinear dynamics of the suspension system is proposed. It requires measurements of the vertical positions of the car body and the tire, and

On the other hand, many dynamical systems exhibit a structural property called differential flatness. This property is equivalent to the existence of a set of independent outputs, called flat outputs and equal in number to the control inputs, which completely parameterizes every state variable and control input (Fliess et al., 1995). By means of differential flatness

**1. Introduction** 

information of the system obtained from sensors.

the estimation of other states and of the road profile.

several tools related with artificial intelligence could be employed such as expert systems, neural networks, machine learning and natural language processing (HAYKIN, 2009).

Two methodologies, reinforcement learning and genetic algorithms, are described and combined. The methodologies had been implemented in a Matlab and were successfully applied in order to create rules that keep the vehicle away from the obstacles from random initial rules. This Chapter also described how the ideas developed by (STAFYLOPATIS, 1998) could be modified, particularly the setting of parameters and a new reinforcement learning methodology, to provide a better integration of the Genetic Algorithms (GAs) and Reinforcement Learning (RL) which produce rules for automatic guided vehicles with a better performance.

The main contributions of this Chapter were: the vehicle, its sensors, and also the environment for training which are different from the ones presented in (STAFYLOPATIS, 1998); a new equation for the reinforcement learning was proposed; the influence of the parameters that control the production of automatic generation of rules for vehicle control navigation are also tested. The results presented attested the better performance of the new reinforcement learning scheme proposed, implemented in Matlab. For future works, it could be applied a genetic algorithm for finding the better configuration of the parameters for the combined scheme of reinforcement learning and the genetic algorithm used to find the rules.

For future works more tests could be performed by using other equations for the reinforcement learning process and also other methods like Beam Search (SABUNCUOGLU and BAVIZ, 1999; DELLA CROCE and T'KINDT, 2002) could replace the Genetic Algorithm. A more detailed and complex decisions could also be incorporated such as the possibility of increasing or reducing the velocity of the vehicle.

#### **8. References**


## **Robust Control of Active Vehicle Suspension Systems Using Sliding Modes and Differential Flatness with MATLAB**

Esteban Chávez Conde1, Francisco Beltrán Carbajal2,

Antonio Valderrábano González3 and Ramón Chávez Bracamontes3

*1Universidad del Papaloapan, Campus Loma Bonita* 

*2Universidad Autónoma Metropolitana, Unidad Azcapotzalco, Departamento de Energía 3Universidad Politécnica de la Zona Metropolitana de Guadalajara 4Instituto Tecnológico de Cd. Guzmán México* 

#### **1. Introduction**

424 MATLAB for Engineers – Applications in Control, Electrical Engineering, IT and Robotics

several tools related with artificial intelligence could be employed such as expert systems, neural networks, machine learning and natural language processing (HAYKIN, 2009). Two methodologies, reinforcement learning and genetic algorithms, are described and combined. The methodologies had been implemented in a Matlab and were successfully applied in order to create rules that keep the vehicle away from the obstacles from random initial rules. This Chapter also described how the ideas developed by (STAFYLOPATIS, 1998) could be modified, particularly the setting of parameters and a new reinforcement learning methodology, to provide a better integration of the Genetic Algorithms (GAs) and Reinforcement Learning (RL) which produce rules for automatic guided vehicles with a

The main contributions of this Chapter were: the vehicle, its sensors, and also the environment for training which are different from the ones presented in (STAFYLOPATIS, 1998); a new equation for the reinforcement learning was proposed; the influence of the parameters that control the production of automatic generation of rules for vehicle control navigation are also tested. The results presented attested the better performance of the new reinforcement learning scheme proposed, implemented in Matlab. For future works, it could be applied a genetic algorithm for finding the better configuration of the parameters for the combined scheme of reinforcement learning and the genetic algorithm

For future works more tests could be performed by using other equations for the reinforcement learning process and also other methods like Beam Search (SABUNCUOGLU and BAVIZ, 1999; DELLA CROCE and T'KINDT, 2002) could replace the Genetic Algorithm. A more detailed and complex decisions could also be incorporated such as the

Della Croce, F.; T'Kind, V., A Recovering Beam Search Algorithm for the One-Machine

Stafylopatis, A., Blekas, K., Autonomous vehicle navigation using evolutionary

Shim, D. H, Kim, H. J., Sastry, S., Hierarchical Control System Synthesis for Rotorcraft-based

Holland, J. H., Adaptation in natural and artificial systems. The University of Michigan

Michalewicz, Z. Genetic Algorithms + Data Structures = Evolution Programs, 3rd edition,

Sabuncuoglu, I., Baviz, M., Job Shop Scheduling with Beam Search, European Journal of

Haykin, , S. Neural Networks and Learning Machines, 3rd edition, Prentice Hall, 2009. Sutton, R.S., Barto, A.G., *Reinforcement Learning: An Introduction*, MIT Press, 1998.

Dynamic Total Completation Time Scheduling Problem, Journal of the Operational

reinforcement learning, European Journal of Operational Research, Vol. 108(2), p.

Unmanned Aerial Vehicles, AIAA Guidance, Navigation, and Control Conference

possibility of increasing or reducing the velocity of the vehicle.

Research Society, vol. 54, p. 1275-1280, 2002.

Operational Research, vol. 118, pp. 390-412.

better performance.

used to find the rules.

**8. References** 

306-318, 1998.

Press, 1975.

and Exhibit, v. 1, p. 1-9, 2000.

Springer-Verlag, 1996.

The main control objectives of active vehicle suspension systems are to improve the ride comfort and handling performance of the vehicle by adding degrees of freedom to the system and/or controlling actuator forces depending on feedback and feedforward information of the system obtained from sensors.

Passenger comfort is provided by isolating the passengers from undesirable vibrations induced from irregular road disturbances, and its performance is evaluated by the level of acceleration which vehicle passengers are exposed. Handling performance is achieved by maintaining a good contact between the tire and the road to provide guidance along the track.

The topic of active vehicle suspension control system, for linear and nonlinear models, in general, has been quite challenging over the years and we refer the reader to some of the fundamental work in the area which has been helpful in the preparation of this chapter. Control strategies like Linear Quadratic Regulator (LQR) in combination with nonlinear backstepping control techniques are proposed in (Liu et al., 2006). This strategy requires information about the state vector (vertical positions and speeds of the tire and car body). A reduced order controller is proposed in (Yousefi et al., 2006) to decrease the implementation costs without sacrificing the security and the comfort by using accelerometers for measurements of the vertical movement of the tire and car body. In (Tahboub, 2005) a controller of variable gain that considers the nonlinear dynamics of the suspension system is proposed. It requires measurements of the vertical positions of the car body and the tire, and the estimation of other states and of the road profile.

On the other hand, many dynamical systems exhibit a structural property called differential flatness. This property is equivalent to the existence of a set of independent outputs, called flat outputs and equal in number to the control inputs, which completely parameterizes every state variable and control input (Fliess et al., 1995). By means of differential flatness

Robust Control of Active Vehicle Suspension Systems

system.

Using Sliding Modes and Differential Flatness with MATLAB 427

Fig. 1. Schematic diagram of a quarter-vehicle suspension system: (a) passive suspension system, (b) electromagnetic active suspension system and (c) hydraulic active suspension

**2.2 Linear mathematical model of the electromagnetic active suspension system**  A schematic diagram of a quarter-car active suspension system is illustrated in Fig.1 (b). The electromagnetic actuator replaces the damper, forming a suspension with the spring. The friction force of an electromagnetic actuator is neglected. The mathematical model of the

electromagnetic suspension system, presented in (Martins et al., 2006), is given by:

**2.3 Linear mathematical model of hydraulic active suspension system** 

mathematical model of the hydraulic suspension system is given by

between the hydraulic force *FA* and the friction force *Ff* .

where *ms* , *mu* , *<sup>s</sup> k* , *<sup>t</sup> k* , *<sup>s</sup> z* , *uz* and *<sup>r</sup> z* represent the same parameters and variables shown in the passive suspension system. The electromagnetic actuator force is represented by *FA* .

A schematic diagram of an active quarter-car suspension system is shown in Fig. 1(c). The

where *ms* , *mu* , . *<sup>s</sup> k* ., *<sup>t</sup> k* , *<sup>s</sup> z* , *<sup>u</sup> z* and *<sup>r</sup> z* represent the same parameters and variables shown in the passive suspension system. The hydraulic actuator force is represented by *FA* , and *Ff* represents the friction force generated by the seals of the piston with the cylinder wall inside the actuator. This friction force has a significant magnitude (> 200 ) *N* and cannot be ignored (Martins et al., 2006; Yousefi et al., 2006). The net force given by the actuator is the difference

*mz k z z F ss s s u A* + − ( )= (3)

*mz k z z k z z F uu s s u t u r A* − −+ − − ( ) ( )= (4)

*mz c z z k z z F F ss s s u s s u f A* + − + − −+ ( ) ( )= (5)

*mz c z z k z z k z z F F uu s s u s s u t u r f A* − −− −+ − − ( ) ( ) ( )= (6)

the analysis and design of controller is greatly simplified. In particular, the combination of differential flatness with sliding modes, which is extensively used when a robust control scheme is required, e.g., parameter uncertainty, exogenous disturbances and un-modeled dynamics (see Utkin, 1978), qualifies as an adequate robust control design approach to get high vibration attenuation level in active vehicle suspension systems. Sliding mode control of a differentially flat system of two degrees of freedom, with vibration attenuation, is presented in (Enríquez-Zárate et al., 2000).

This chapter presents a robust active vibration control scheme based on sliding modes and differential flatness for electromagnetic and hydraulic active vehicle suspension systems. Measurements of the vertical displacements of the car body and the tire are required for implementation of the proposed control scheme. On-line algebraic estimation of the states variables is used to avoid the use of sensors of acceleration and velocity. The road profile is considered as an unknown input disturbance that cannot be measured. Simulation results obtained from Matlab are included to show the dynamic performance and robustness of the active suspension systems with the proposed control scheme. This chapter applies the algebraic state estimation scheme proposed by Fliess and Sira-Ramírez (Fliess & Sira-Ramírez, 2004a, 2004b; Sira-Ramírez & Silva-Navarro, 2003) for control of nonlinear systems, which is based on the algebraic identification methodology of system parameters reported in (Fliess & Sira-Ramírez, 2003). The method is purely algebraic and involves the use of differential algebra. This method is applied to obtain an estimate of the time derivative from any signal, avoiding model reliance of the system at least in the estimation of states. Simulation and experimental results of the on-line algebraic estimation of states on a differentially flat system of two degrees of freedom are presented in (García-Rodríguez, 2005).

This chapter is organized as follows: Section 2 presents the linear mathematical models of vehicle suspension systems of a quarter car. The design of the controllers for the active suspension systems are introduced in Sections 3 and 4. Section 5 divulges the design of the algebraic estimator of states, while Section 6 shows the use of sensors for measuring the variables required by the controller. The simulation results are illustrated in Section 7. Finally, conclusions are brought out in Section 8.

#### **2. Dynamic model of quarter-car suspension systems**

#### **2.1 Linear mathematical model of passive suspension system**

A schematic diagram of a quarter-car suspension system is shown in Fig. 1(a). The mathematical model of the passive suspension system is given by

$$m\_s \ddot{z}\_s + c\_s(\dot{z}\_s - \dot{z}\_u) + k\_s(z\_s - z\_u) = 0 \tag{1}$$

$$k m\_u \ddot{z}\_u - c\_s (\dot{z}\_s - \dot{z}\_u) - k\_s (z\_s - z\_u) + k\_t (z\_u - z\_r) = 0 \tag{2}$$

where *ms* represents the mass of a quarter car, *mu* represents the mass of one wheel with the suspension and brake equipment, *<sup>s</sup> c* is the damper coefficient of suspension, *<sup>s</sup> k* and *<sup>t</sup> k* are the spring coefficients of the suspension and the tire, *<sup>s</sup> z* and *uz* are the displacements of car body and the wheel and *<sup>r</sup> z* is the terrain input disturbance. This simplified linear mathematical model of a passive suspension system has been widely used in many previous works, such as (Liu et al., 2006; Yousefi et al., 2006).

the analysis and design of controller is greatly simplified. In particular, the combination of differential flatness with sliding modes, which is extensively used when a robust control scheme is required, e.g., parameter uncertainty, exogenous disturbances and un-modeled dynamics (see Utkin, 1978), qualifies as an adequate robust control design approach to get high vibration attenuation level in active vehicle suspension systems. Sliding mode control of a differentially flat system of two degrees of freedom, with vibration attenuation, is

This chapter presents a robust active vibration control scheme based on sliding modes and differential flatness for electromagnetic and hydraulic active vehicle suspension systems. Measurements of the vertical displacements of the car body and the tire are required for implementation of the proposed control scheme. On-line algebraic estimation of the states variables is used to avoid the use of sensors of acceleration and velocity. The road profile is considered as an unknown input disturbance that cannot be measured. Simulation results obtained from Matlab are included to show the dynamic performance and robustness of the active suspension systems with the proposed control scheme. This chapter applies the algebraic state estimation scheme proposed by Fliess and Sira-Ramírez (Fliess & Sira-Ramírez, 2004a, 2004b; Sira-Ramírez & Silva-Navarro, 2003) for control of nonlinear systems, which is based on the algebraic identification methodology of system parameters reported in (Fliess & Sira-Ramírez, 2003). The method is purely algebraic and involves the use of differential algebra. This method is applied to obtain an estimate of the time derivative from any signal, avoiding model reliance of the system at least in the estimation of states. Simulation and experimental results of the on-line algebraic estimation of states on a differentially flat system of two degrees of freedom are presented

This chapter is organized as follows: Section 2 presents the linear mathematical models of vehicle suspension systems of a quarter car. The design of the controllers for the active suspension systems are introduced in Sections 3 and 4. Section 5 divulges the design of the algebraic estimator of states, while Section 6 shows the use of sensors for measuring the variables required by the controller. The simulation results are illustrated in Section 7.

A schematic diagram of a quarter-car suspension system is shown in Fig. 1(a). The

where *ms* represents the mass of a quarter car, *mu* represents the mass of one wheel with the suspension and brake equipment, *<sup>s</sup> c* is the damper coefficient of suspension, *<sup>s</sup> k* and *<sup>t</sup> k* are the spring coefficients of the suspension and the tire, *<sup>s</sup> z* and *uz* are the displacements of car body and the wheel and *<sup>r</sup> z* is the terrain input disturbance. This simplified linear mathematical model of a passive suspension system has been widely used in many previous

+ −+ − (1)

( ) ( ) ( )=0 *mz c z z k z z k z z uu s s u s s u t u r* − −− −+ − (2)

presented in (Enríquez-Zárate et al., 2000).

in (García-Rodríguez, 2005).

Finally, conclusions are brought out in Section 8.

works, such as (Liu et al., 2006; Yousefi et al., 2006).

**2. Dynamic model of quarter-car suspension systems 2.1 Linear mathematical model of passive suspension system** 

mathematical model of the passive suspension system is given by

( ) ( )=0 *mz c z z k z z ss s s u s s u*

Fig. 1. Schematic diagram of a quarter-vehicle suspension system: (a) passive suspension system, (b) electromagnetic active suspension system and (c) hydraulic active suspension system.

#### **2.2 Linear mathematical model of the electromagnetic active suspension system**

A schematic diagram of a quarter-car active suspension system is illustrated in Fig.1 (b). The electromagnetic actuator replaces the damper, forming a suspension with the spring. The friction force of an electromagnetic actuator is neglected. The mathematical model of the electromagnetic suspension system, presented in (Martins et al., 2006), is given by:

$$m\_s \ddot{z}\_s + k\_s(z\_s - z\_u) = F\_A \tag{3}$$

$$m\_u \ddot{z}\_u - k\_s(z\_s - z\_u) + k\_t(z\_u - z\_r) = -F\_A \tag{4}$$

where *ms* , *mu* , *<sup>s</sup> k* , *<sup>t</sup> k* , *<sup>s</sup> z* , *uz* and *<sup>r</sup> z* represent the same parameters and variables shown in the passive suspension system. The electromagnetic actuator force is represented by *FA* .

#### **2.3 Linear mathematical model of hydraulic active suspension system**

A schematic diagram of an active quarter-car suspension system is shown in Fig. 1(c). The mathematical model of the hydraulic suspension system is given by

$$m\_s \ddot{z}\_s + c\_s(\dot{z}\_s - \dot{z}\_u) + k\_s(z\_s - z\_u) = -F\_f + F\_A \tag{5}$$

$$m\_u \ddot{z}\_u - c\_s(\dot{z}\_s - \dot{z}\_u) - k\_s(z\_s - z\_u) + k\_t(z\_u - z\_r) = F\_f - F\_A \tag{6}$$

where *ms* , *mu* , . *<sup>s</sup> k* ., *<sup>t</sup> k* , *<sup>s</sup> z* , *<sup>u</sup> z* and *<sup>r</sup> z* represent the same parameters and variables shown in the passive suspension system. The hydraulic actuator force is represented by *FA* , and *Ff* represents the friction force generated by the seals of the piston with the cylinder wall inside the actuator. This friction force has a significant magnitude (> 200 ) *N* and cannot be ignored (Martins et al., 2006; Yousefi et al., 2006). The net force given by the actuator is the difference between the hydraulic force *FA* and the friction force *Ff* .

Robust Control of Active Vehicle Suspension Systems

**3.2 Sliding mode and differential flatness control** 

following form:

where 1 = *<sup>u</sup>*

*t <sup>m</sup> <sup>d</sup>*

The design gains 2 0

polynomial 3 2 + ++ β

switching surface

shown in Section 5.

where μ , γ

*<sup>k</sup>* , 2 = 1 + + *su s ts t km k <sup>d</sup> km k*

Then, the error dynamics restricted to

β

σ

sliding surface is globally attractive,

following sliding-mode controller:

 ββ

 β

Now, consider a linear switching surface defined by

The input *u* in terms of the flat output and its time derivatives is given by

and 3 = *<sup>s</sup>*

(3)

σ

σ

*<sup>k</sup> <sup>d</sup> <sup>m</sup>* .

Using Sliding Modes and Differential Flatness with MATLAB 429

(4) = 1 *u su s <sup>s</sup> t ts t s m km k k u F F F k km k m* + ++ +

(4) = 1 *u su s <sup>s</sup> t ts t s m km k k u F F F k km k m* + ++ +

where (4) *F v* = defines an auxiliary control input. This expression can be written in the

*s*

 = *F FFF* +++ βββ210

(3) *F FFF* +++ βββ

(8)

(4) *u dF dF dF* = 1 23 + + (9)

= 0 is governed by the linear differential equation

, , are selected to verify that the associated characteristic

2 10 *sss* be Hurwitz. As a consequence, the error dynamics on the

denote positive real constants and "sign" is the standard signum function. The

σ

σ γ ( ) σ

= 0 is globally asymptotically stable. The sliding surface

globally attractive with the continuous approximation to the discontinuous sliding mode

< 0 for

condition for the existence of sliding mode presented in (Utkin, 1978). One then obtains the

[ ] (3) = 2 10 *v F F F sign* − −−− +

This controller requires the measurement of all the state variables of the suspension system, *<sup>s</sup> z* , *<sup>s</sup> z* , *<sup>u</sup> z* and *<sup>u</sup> z* , corresponding to the vertical positions and velocity of the body of the car and the wheel. The variables *<sup>s</sup> z* and *<sup>u</sup> z* are calculated through an online algebraic estimator,

ββμ

controller as given in (Sira-Ramírez, 1993), i.e., by forcing to satisfy the dynamics,

 = − + μ [σ γ *sign*( ) σ

σσ

σ

β

<sup>210</sup> = 0 (11)

*u dv dF dF* = 12 3 + + (13)

(10)

σ

] (12)

≠ 0 , which is a very well known

= 0 is made

#### **3. Control of electromagnetic suspension system**

The mathematical model of the electromagnetic active suspension system illustrated in Fig. 1(b) is given by the equations (3) and (4). Defining the state variables <sup>1</sup> = *<sup>s</sup> x z* , <sup>2</sup> = *<sup>s</sup> x z* , 3 = *<sup>u</sup> x z* and <sup>4</sup> = *<sup>u</sup> x z* for the model of the equations mentioned, the representation in the state space form is,

$$\dot{\mathbf{x}}(t) = Ax(t) + Bu(t) + Ez\_r(t); \qquad x(t) \in \mathbb{R}^4, A \in \mathbb{R}^{4 \times 4}, B \in \mathbb{R}^{4 \times 1}, E \in \mathbb{R}^{4 \times 1}$$

$$\begin{bmatrix} \dot{\mathbf{x}}\_1 \\ \dot{\mathbf{x}}\_2 \\ \dot{\mathbf{x}}\_3 \\ \dot{\mathbf{x}}\_4 \end{bmatrix} = \begin{bmatrix} 0 & 1 & 0 & 0 \\ -\frac{k\_s}{m\_s} & 0 & \frac{k\_s}{m\_s} & 0 \\ 0 & 0 & 0 & 1 \\ \frac{k\_s}{m\_u} & 0 & -\frac{k\_s + k\_t}{m\_u} & 0 \end{bmatrix} \begin{bmatrix} \mathbf{x}\_1 \\ \mathbf{x}\_2 \\ \mathbf{x}\_3 \\ \mathbf{x}\_4 \end{bmatrix} + \begin{bmatrix} 0 \\ \frac{1}{m\_s} \\ 0 \\ 0 \\ -\frac{1}{m\_u} \end{bmatrix} u + \begin{bmatrix} 0 \\ 0 \\ 0 \\ \frac{k\_t}{m\_u} \end{bmatrix} z\_r \tag{7}$$

with *u F* = *<sup>A</sup>* , the force provided by the electromagnetic actuator as control input.

#### **3.1 Differential flatness**

The system is controllable and hence, flat (Fliess et al., 1995; Sira-Ramírez & Agrawal, 2004), with the flat output being the positions of body car and wheel, *F mx mx* = *s u* 1 3 + , in (Chávez-Conde et al., 2009). For simplicity in the analysis of the differential flatness for the suspension system assume that = 0 *t r k z* . In order to show the parameterization of all the state variables and control input, we firstly compute the time derivatives up to fourth order for *F* , resulting in

$$\begin{aligned} F &= m\_s \mathbf{x}\_1 + m\_u \mathbf{x}\_3 \\ \dot{F} &= m\_s \mathbf{x}\_2 + m\_u \mathbf{x}\_4 \\ \ddot{F} &= -k\_t \mathbf{x}\_3 \\ F^{(3)} &= -k\_r \mathbf{x}\_4 \\ F^{(4)} &= \frac{k\_t}{m\_u} u - \frac{k\_s k\_t}{m\_u} (\mathbf{x}\_1 - \mathbf{x}\_3) + \frac{k\_t^2}{m\_u} \mathbf{x}\_3 \end{aligned}$$

Then, the state variables and control input are differentially parameterized in terms of the flat output as follows

$$\begin{aligned} \boldsymbol{\chi}\_1 &= \frac{1}{m\_s} \left( \boldsymbol{F} + \frac{m\_u}{k\_\varepsilon} \ddot{\boldsymbol{F}} \right) \\ \boldsymbol{\chi}\_2 &= \frac{1}{m\_s} \left( \ddot{\boldsymbol{F}} + \frac{m\_u}{k\_\varepsilon} \boldsymbol{F}^{(3)} \right) \\ \boldsymbol{\chi}\_3 &= -\frac{1}{k\_\varepsilon} \ddot{\boldsymbol{F}} \\ \boldsymbol{\chi}\_4 &= -\frac{1}{k\_\varepsilon} \boldsymbol{F}^{(3)} \end{aligned}$$

The mathematical model of the electromagnetic active suspension system illustrated in Fig. 1(b) is given by the equations (3) and (4). Defining the state variables <sup>1</sup> = *<sup>s</sup> x z* , <sup>2</sup> = *<sup>s</sup> x z* , 3 = *<sup>u</sup> x z* and <sup>4</sup> = *<sup>u</sup> x z* for the model of the equations mentioned, the representation in the state space

4 44 41 41 ( ) = ( ) ( ) ( ); ( ) , , , , *<sup>r</sup> x t Ax t Bu t Ez t x t A B E* ××× ++ ∈∈ ∈ ∈

<sup>=</sup> <sup>0</sup> 00 0 1 0

*s s s*

<sup>−</sup> + + <sup>+</sup> − −

*uu u*

*mm m*

*<sup>k</sup> x k kk x*

<sup>1</sup> 0 0

The system is controllable and hence, flat (Fliess et al., 1995; Sira-Ramírez & Agrawal, 2004), with the flat output being the positions of body car and wheel, *F mx mx* = *s u* 1 3 + , in (Chávez-Conde et al., 2009). For simplicity in the analysis of the differential flatness for the suspension system assume that = 0 *t r k z* . In order to show the parameterization of all the state variables and control input, we firstly compute the time derivatives up to fourth order for *F* , resulting in

*s u s u t t*

+ +

*F mx mx F mx mx*

4

( )

*u*

 + + 

*u*

*t st t u u u*

− −+

*k kk k F u xx x mm m*

Then, the state variables and control input are differentially parameterized in terms of the

*s t*

*<sup>m</sup> x FF m k*

*<sup>m</sup> x FF m k*

*s t*

(3)

2

13 3

(3)

01 0 0 0 <sup>0</sup> <sup>1</sup> 0 0 <sup>0</sup>

*r*

(7)

*t*

*m*

*u z*

*u*

1 1 2 2 3 3 4 4

*x x mm m*

with *u F* = *<sup>A</sup>* , the force provided by the electromagnetic actuator as control input.

*s st*

(3)

  = = = =

*F kx F kx*

− −

1

<sup>1</sup> <sup>=</sup>

<sup>1</sup> <sup>=</sup>

<sup>1</sup> <sup>=</sup>

−

*x F k*

<sup>1</sup> <sup>=</sup>

−

*x F k*

*t*

*t*

2

3

4

(4)

=

*s s*

*x x*

*x x k k*

**3. Control of electromagnetic suspension system** 

 

**3.1 Differential flatness** 

flat output as follows

form is,

$$
\mu\_{\mathcal{U}} = \frac{m\_{\mu}}{k\_{\iota}} F^{(4)} + \left(\frac{k\_{s} m\_{\mu}}{k\_{\iota} m\_{s}} + \frac{k\_{s}}{k\_{\iota}} + 1\right) \ddot{F} + \frac{k\_{s}}{m\_{s}} F^{(4)}
$$

#### **3.2 Sliding mode and differential flatness control**

The input *u* in terms of the flat output and its time derivatives is given by

$$
\Delta \mu = \frac{m\_u}{k\_t} F^{(4)} + \left(\frac{k\_s m\_u}{k\_t m\_s} + \frac{k\_s}{k\_t} + 1\right) \vec{F} + \frac{k\_s}{m\_s} F \tag{8}
$$

where (4) *F v* = defines an auxiliary control input. This expression can be written in the following form:

$$
\mu = d\_1 F^{(4)} + d\_2 \ddot{F} + d\_3 F \tag{9}
$$

where 1 = *<sup>u</sup> t <sup>m</sup> <sup>d</sup> <sup>k</sup>* , 2 = 1 + + *su s ts t km k <sup>d</sup> km k* and 3 = *<sup>s</sup> s <sup>k</sup> <sup>d</sup> <sup>m</sup>* .

Now, consider a linear switching surface defined by

$$
\sigma = F^{(3)} + \mathcal{B}\_2 \dot{F} + \mathcal{B}\_1 \dot{F} + \mathcal{B}\_0 F \tag{10}
$$

Then, the error dynamics restricted toσ= 0 is governed by the linear differential equation

$$F^{(3)} + \mathcal{B}\_2 \ddot{F} + \mathcal{B}\_1 \dot{F} + \mathcal{B}\_0 F = 0 \tag{11}$$

The design gains 2 0 β β , , are selected to verify that the associated characteristic polynomial 3 2 + ++ β ββ 2 10 *sss* be Hurwitz. As a consequence, the error dynamics on the switching surface σ = 0 is globally asymptotically stable. The sliding surface σ = 0 is made globally attractive with the continuous approximation to the discontinuous sliding mode controller as given in (Sira-Ramírez, 1993), i.e., by forcing to satisfy the dynamics,

$$\dot{\sigma} = -\mu \left[ \sigma + \gamma \text{sign}(\sigma) \right] \tag{12}$$

where μ , γ denote positive real constants and "sign" is the standard signum function. The sliding surface is globally attractive, σσ < 0 for σ ≠ 0 , which is a very well known condition for the existence of sliding mode presented in (Utkin, 1978). One then obtains the following sliding-mode controller:

$$
\mu \equiv d\_1 \upsilon + d\_2 \dot{F} + d\_3 F \tag{13}
$$

$$
\upsilon = -\beta\_2 F^{(3)} - \beta\_1 \ddot{F} - \beta\_0 \dot{F} - \mu \left[ \sigma + \chi \text{sign}(\sigma) \right]
$$

This controller requires the measurement of all the state variables of the suspension system, *<sup>s</sup> z* , *<sup>s</sup> z* , *<sup>u</sup> z* and *<sup>u</sup> z* , corresponding to the vertical positions and velocity of the body of the car and the wheel. The variables *<sup>s</sup> z* and *<sup>u</sup> z* are calculated through an online algebraic estimator, shown in Section 5.

Robust Control of Active Vehicle Suspension Systems

**4.2 Sliding mode and differential flatness control** 

<sup>2</sup> = + *su s ts t cm c km k* ,

Then, the error dynamics restricted to

σ

sliding surface is globally attractive,

following sliding-mode controller:

estimator, shown in Section 5.

β

 ββ

 β

The design gains 2 0

polynomial 3 2 + ++ β

switching surface

where μ , γ

Now, consider a linear switching surface defined by

following form:

where η<sup>1</sup> = *<sup>u</sup> t m k* , η

The input *u* in terms of the flat output and its time derivatives is given by

η

σ

η

η

<sup>3</sup> = 1 + + *su s ts t km k km k* ,

(3)

σ

controller as given in (Sira-Ramírez, 1993), i.e., by forcing to satisfy the dynamics

 = − + μ [σ γ *sign*( ) σ

σσ

η

σ

η

β

Using Sliding Modes and Differential Flatness with MATLAB 431

(4) (3) = 1 *u su s su s s s t ts t ts t s s m cm c km k c k uF F FFF k km k km k m m* 

(3) = 1 *u su s su s s s t ts t ts t s s m cm c km k c k uv F FFF k km k km k m m*

where (4) *F v* = defines the auxiliary control input. The expression can be written in the

(3) *uvF FFF* =

 = *F FFF* +++ βββ210

(3) *F FFF* +++ βββ

12 3 45 + +++

ηηη

 and η<sup>5</sup> = *<sup>s</sup> s k <sup>m</sup>* .

, , are selected to verify that the associated characteristic

2 10 *sss* be Hurwitz. As a consequence, the error dynamics on the

denote positive real constants and "sign" is the standard signum function. The

σ

ηηη

σ γ ( ) σ

= 0 is globally asymptotically stable. The sliding surface

globally attractive with the continuous approximation to the discontinuous sliding mode

< 0 for

condition for the existence of sliding mode presented in (Utkin, 1978). One then obtains the

(3) *uvF FFF* =

[ ] (3) = 2 10 *v F F F sign* − −−− +

This controller requires the measurement of all the variables of state of suspension system, *<sup>s</sup> z* , *<sup>s</sup> z* , *<sup>u</sup> z* and *<sup>u</sup> z* corresponding to the vertical positions and velocity of the body of the car and the tire, respectively. The variables *<sup>s</sup> z* and *<sup>u</sup> z* are calculated through an online algebraic

ββμ

12 3 4 5 + +++

η<sup>4</sup> = *<sup>s</sup> s c m*

+ + + ++ + +

+ + + ++ + +

(15)

σ

] (19)

≠ 0 , which is a very well known

(20)

= 0 is made

(16)

(17)

= 0 is governed by the linear differential equation

<sup>210</sup> = 0 (18)

#### **4. Control control of hydraulic suspension system**

The mathematical model of the hydraulic active suspension system shown in Fig. 1(c) is given by the equations (5) and (6). Defining the state variables 1 = *<sup>s</sup> x z* , 2 = *<sup>s</sup> x z* , 3 = *<sup>u</sup> x z* and <sup>4</sup> = *<sup>u</sup> x z* for the model of the equations mentioned, the representation in the state space form is, 4 44 41 41 ( ) = ( ) ( ) ( ); ( ) , , , , *<sup>r</sup> x t Ax t Bu t Ez t x t A B E* ××× ++ ∈∈ ∈ ∈

$$
\begin{bmatrix}
\dot{\boldsymbol{x}}\_{1} \\
\dot{\boldsymbol{x}}\_{2} \\
\dot{\boldsymbol{x}}\_{3} \\
\dot{\boldsymbol{x}}\_{4}
\end{bmatrix} = \begin{bmatrix}
0 & 1 & 0 & 0 \\
0 & 0 & 0 & 1 \\
\frac{k\_{s}}{m\_{u}} & \frac{c\_{s}}{m\_{u}} & -\frac{k\_{s} + k\_{t}}{m\_{u}} & -\frac{c\_{s}}{m\_{u}}
\end{bmatrix} \begin{bmatrix}
\boldsymbol{x}\_{1} \\
\boldsymbol{x}\_{2} \\
\boldsymbol{x}\_{3} \\
\boldsymbol{x}\_{4}
\end{bmatrix} + \begin{bmatrix}
0 \\
\frac{\boldsymbol{m}}{m\_{s}} \\
0 \\
0 \\
\end{bmatrix} \boldsymbol{u} + \begin{bmatrix}
0 \\
0 \\
0 \\
\frac{k\_{t}}{m\_{u}} \\
m\_{u}
\end{bmatrix} \tag{14}
$$

with *uF F* = *A f* − , the net force provided by the hydraulic actuator as control input (the net force provided by the actuator is the difference between the hydraulic force *FA* and the frictional force *Ff* ).

#### **4.1 Differential flatness**

The system is controllable and hence, flat (Fliess et al., 1995; Sira-Ramírez & Agrawal, 2004), with the flat output being the positions of body car and wheel, *F mx mx* = *s u* 1 3 + , in (Chávez-Conde et al., 2009). For simplicity in the analysis of the differential flatness for the suspension system assume that = 0 *t r k z* . In order to show the parameterization of all the state variables and control input, we firstly compute the time derivatives up to fourth order for *F* , resulting in

$$\begin{aligned} F &= m\_s \mathbf{x}\_1 + m\_u \mathbf{x}\_3 \\ \dot{F} &= m\_s \mathbf{x}\_2 + m\_u \mathbf{x}\_4 \\ \dot{F} &= -k\_r \mathbf{x}\_3 \\ F^{(3)} &= -k\_t \mathbf{x}\_4 \\ F^{(4)} &= \frac{k\_t}{m\_u} u - \frac{c\_s k\_t}{m\_u} (\mathbf{x}\_2 - \mathbf{x}\_4) - \frac{k\_s k\_t}{m\_u} (\mathbf{x}\_1 - \mathbf{x}\_3) + \frac{k\_t^2}{m\_u} \mathbf{x}\_3 \end{aligned}$$

Then, the state variables and control input are parameterized in terms of the flat output as follows

$$\begin{aligned} \boldsymbol{x}\_{1} &= \frac{1}{m\_{\boldsymbol{s}}} \Big( \boldsymbol{F} + \frac{m\_{\boldsymbol{u}}}{k\_{\boldsymbol{t}}} \ddot{\boldsymbol{F}} \Big) \\ \boldsymbol{x}\_{2} &= \frac{1}{m\_{\boldsymbol{s}}} \Big( \ddot{\boldsymbol{F}} + \frac{m\_{\boldsymbol{u}}}{k\_{\boldsymbol{t}}} \boldsymbol{F}^{(3)} \Big) \\ \boldsymbol{x}\_{3} &= -\frac{1}{k\_{\boldsymbol{t}}} \ddot{\boldsymbol{F}} \\ \boldsymbol{x}\_{4} &= -\frac{1}{k\_{\boldsymbol{t}}} \boldsymbol{F}^{(3)} \end{aligned}$$

The mathematical model of the hydraulic active suspension system shown in Fig. 1(c) is given by the equations (5) and (6). Defining the state variables 1 = *<sup>s</sup> x z* , 2 = *<sup>s</sup> x z* , 3 = *<sup>u</sup> x z* and <sup>4</sup> = *<sup>u</sup> x z* for the model of the equations mentioned, the representation in the state space form

01 0 0 0 <sup>0</sup>

<sup>=</sup> <sup>0</sup> 00 0 1 0

*ss s s s*

− − + + <sup>+</sup> −− −

*uu u u u*

with *uF F* = *A f* − , the net force provided by the hydraulic actuator as control input (the net force provided by the actuator is the difference between the hydraulic force *FA* and the

The system is controllable and hence, flat (Fliess et al., 1995; Sira-Ramírez & Agrawal, 2004), with the flat output being the positions of body car and wheel, *F mx mx* = *s u* 1 3 + , in (Chávez-Conde et al., 2009). For simplicity in the analysis of the differential flatness for the suspension system assume that = 0 *t r k z* . In order to show the parameterization of all the state variables and control input, we firstly compute the time derivatives up to fourth order for *F* , resulting in

= ( )( )

*k ck kk k F u xx xx x mm m m*

Then, the state variables and control input are parameterized in terms of the flat output as

*s t*

*<sup>m</sup> x FF m k*

*<sup>m</sup> x FF m k*

*s t*

(3)

*t st s t t u u u u*

− −− −+

*u*

 + + 

*u*

*mm m m m*

*<sup>k</sup> x k c kk c x*

1

1

2

24 13 3

(3)

0

*u z*

*t*

*m*

*u*

*r*

(14)

**4. Control control of hydraulic suspension system**

 

frictional force *Ff* ).

follows

**4.1 Differential flatness** 

is, 4 44 41 41 ( ) = ( ) ( ) ( ); ( ) , , , , *<sup>r</sup> x t Ax t Bu t Ez t x t A B E* ××× ++ ∈∈ ∈ ∈

1 1 2 2 3 3 4 4

*x x kc k c*

*x x*

*ss s s*

*x x mm m m m*

*s s st s*

1 3 2 4

1

<sup>1</sup> <sup>=</sup>

<sup>1</sup> <sup>=</sup>

<sup>1</sup> <sup>=</sup>

−

*x F k*

<sup>1</sup> <sup>=</sup>

−

*x F k*

*t*

*t*

2

3

4

*s u s u*

+ +

*F mx mx F mx mx*

3

*t t*

*F kx F kx*

− −

4

(3)

  = = = =

(4)

$$\mathcal{U}u = \frac{m\_u}{k\_t} F^{(4)} + \left(\frac{c\_s m\_u}{k\_t m\_s} + \frac{c\_s}{k\_t}\right) F^{(3)} + \left(\frac{k\_s m\_u}{k\_t m\_s} + \frac{k\_s}{k\_t} + 1\right) \ddot{F} + \frac{c\_s}{m\_s} \dot{F} + \frac{k\_s}{m\_s} F^2$$

#### **4.2 Sliding mode and differential flatness control**

The input *u* in terms of the flat output and its time derivatives is given by

$$\dot{\mathbf{u}} = \frac{m\_u}{k\_l} \mathbf{v} + \left(\frac{c\_s m\_u}{k\_l m\_s} + \frac{c\_s}{k\_t}\right) \mathbf{F}^{(3)} + \left(\frac{k\_s m\_u}{k\_l m\_s} + \frac{k\_s}{k\_t} + 1\right) \ddot{\mathbf{F}} + \frac{c\_s}{m\_s} \dot{\mathbf{F}} + \frac{k\_s}{m\_s} \mathbf{F} \tag{15}$$

where (4) *F v* = defines the auxiliary control input. The expression can be written in the following form:

$$
\mu = \eta\_1 \upsilon + \eta\_2 F^{(3)} + \eta\_3 \ddot{F} + \eta\_4 \dot{F} + \eta\_5 F \tag{16}
$$

where η<sup>1</sup> = *<sup>u</sup> t m k* , η<sup>2</sup> = + *su s ts t cm c km k* , η<sup>3</sup> = 1 + + *su s ts t km k km k* , η<sup>4</sup> = *<sup>s</sup> s c m* and η<sup>5</sup> = *<sup>s</sup> s k <sup>m</sup>* .

Now, consider a linear switching surface defined by

$$
\sigma = F^{(3)} + \mathcal{B}\_2 \dot{\bar{F}} + \mathcal{B}\_1 \dot{F} + \mathcal{B}\_0 F \tag{17}
$$

Then, the error dynamics restricted to σ= 0 is governed by the linear differential equation

$$
\dot{F}^{(3)} + \mathcal{B}\_2 \dot{F} + \mathcal{B}\_1 \dot{F} + \mathcal{B}\_0 F = 0 \tag{18}
$$

The design gains 2 0 β β , , are selected to verify that the associated characteristic polynomial 3 2 + ++ β ββ 2 10 *sss* be Hurwitz. As a consequence, the error dynamics on the switching surface σ = 0 is globally asymptotically stable. The sliding surface σ = 0 is made globally attractive with the continuous approximation to the discontinuous sliding mode controller as given in (Sira-Ramírez, 1993), i.e., by forcing to satisfy the dynamics

$$\sigma = -\mu \left[ \sigma + \gamma \text{sign}(\sigma) \right] \tag{19}$$

where μ , γ denote positive real constants and "sign" is the standard signum function. The sliding surface is globally attractive, σσ < 0 for σ ≠ 0 , which is a very well known condition for the existence of sliding mode presented in (Utkin, 1978). One then obtains the following sliding-mode controller:

$$
\mu = \eta\_1 \upsilon + \eta\_2 F^{(3)} + \eta\_3 \ddot{F} + \eta\_4 \dot{F} + \eta\_5 F \tag{20}
$$

$$
\upsilon = -\beta\_2 F^{(3)} - \beta\_1 \dddot{F} - \beta\_0 \dot{F} - \mu \left[ \sigma + \gamma \text{sign}(\sigma) \right]
$$

This controller requires the measurement of all the variables of state of suspension system, *<sup>s</sup> z* , *<sup>s</sup> z* , *<sup>u</sup> z* and *<sup>u</sup> z* corresponding to the vertical positions and velocity of the body of the car and the tire, respectively. The variables *<sup>s</sup> z* and *<sup>u</sup> z* are calculated through an online algebraic estimator, shown in Section 5.

#### **5. On-line algebraic state estimation of active suspension system**

#### **5.1 First time derivative algebraic estimation**

The algebraic methodology proposed in (Fliess & Sira-Ramírez, 2003) allows us to estimate the derivatives of a smooth signal considering a signal model of *n th* − order, thus it is not necessary to design the time derivative estimator from a specific dynamic model of the plant. Through valid algebraic manipulations of this approximated model in the frequency domain, and using the algebraic derivation with respect to the complex variable *s* , we neglect the initial conditions of the signal. The resulting equation is multiplied by a negative power *<sup>n</sup>*−<sup>1</sup> *s* and returned to the time domain. This last expression is manipulated algebraically in order to find a formula to estimate the first time derivative of *y*( )*t* .

Consider a fourth order approximation of a smooth signal *y*( )*t* ,

$$\frac{d^4y(t)}{dt^4} = 0\tag{21}$$

Robust Control of Active Vehicle Suspension Systems

−

*t*

where ( ,) *<sup>i</sup> t t* is the estimation period.

calculated periodically as follows,

Using Sliding Modes and Differential Flatness with MATLAB 433

This formula is valid for *t* > 0 . Since (25) provides an approximated value of the first derivative, this is only valid during a period of time. So the state estimation should be

> 3 321 4

In order to obtain a better and smoother estimated value of the vertical velocity, we have implemented simultaneously two algebraic estimators for each velocity to estimate. Proceeding with an out-of-phase policy for one of these algebraic estimators, the outputs of

The only variables required for the implementation of the proposed controllers are the vertical displacements of the body of the car *<sup>s</sup> z* and the vertical displacement of the wheel *<sup>u</sup> z* . These variables are needed to measure by some sensor. In (Chamseddine et al., 2006) the use of sensors in experimental vehicle platforms, as well as in commercial vehicles is presented. The most common sensors used for measuring the vertical displacement of the body of the car and the wheel are laser sensors. This type of sensors could be used to measure the variables *<sup>s</sup> z* and *<sup>s</sup> z* needed for controller implementation. An accelerometer or another type of sensor is not needed to measure the variables *<sup>s</sup> z* and *<sup>u</sup> z* , these variables are estimated with the use of algebraic estimators from knowledge of the variables *<sup>s</sup> z* and *<sup>u</sup> z* . Fig. 2 shows a schematic diagram of the instrumentation for the active suspension system.

Fig. 2. Schematic diagram of the instrumentation of the active suspension system.

−

λ λ λλ

λ

*tt z z dd*

*t <sup>i</sup> t t i i*

2 2 21

*dt t t* , ∀ − ( )>0 *<sup>i</sup> t t* (26)

3 1

λ λλλ

12( ) 96 ( )

1 2

24 ( ) <sup>=</sup> ( )

*z ddd dz*

− +

λ λ

*tt t ii i <sup>t</sup> <sup>i</sup> <sup>i</sup>*

both are combined properly to obtain a final estimated value.

**6. Instrumentation of the active suspension system** 

*t*

This model indicates that *y*( )*t* is a signal whose behavior can be approximated by a family of polynomials of third order, thus the fourth time derivative is assumed close to zero. The order of this approximation can be increased to enhance the estimation quality of the algebraic estimator. From (21) it is possible to design a time derivative algebraic estimator. First, we apply Laplace transform to (21),

$$s^4Y(s) - s^3Y(0) - s^2\dot{Y}(0) - s\ddot{Y} - Y^{(3)} = 0\tag{22}$$

Now, taking successive derivatives until a number of three with respect to the complex variable *s* , we obtain a expression which is free of initial conditions,

$$\frac{d^4\left(\mathbf{s}^4\mathbf{Y}\right)}{ds^4} = 0\tag{23}$$

Expanding this expression and multiplying by <sup>−</sup><sup>3</sup> *s* ,

$$24s^{-3}Y + 96s^{-2}\frac{dY}{ds} + 72s^{-1}\frac{d^2Y}{ds^2} + 16\frac{d^3Y}{ds^3} + s\frac{d^4Y}{ds^4} \tag{24}$$

Returning to the time domain,

$$\begin{split} &\frac{d}{dt}(t^4 z(t)) - 16t^3 z(t) + 72 \int\_0^t \lambda\_1^2 z(\lambda\_1) d\lambda\_1 \\ & - 96 \int\_0^t \int\_0^{\lambda\_1} \lambda\_2 z(\lambda\_2) d\lambda\_2 d\lambda\_1 + 24 \int\_0^t \int\_0^{\lambda\_1} \int\_0^{\lambda\_2} z(\lambda\_3) d\lambda\_3 d\lambda\_2 d\lambda\_1 = 0 \end{split}$$

From the last equation is possible to obtain the following algebraic estimator,

$$\frac{dz}{dt} = \frac{12t^3 z + 96 \int\_0^t \int\_0^{\lambda\_1} \lambda\_2 z(\lambda\_2) d\lambda\_2 d\lambda\_1 - 24 \int\_0^t \int\_0^{\lambda\_1} \int\_0^{\lambda\_2} z(\lambda\_3) d\lambda\_3 d\lambda\_2 d\lambda\_1}{t^4} \tag{25}$$

This formula is valid for *t* > 0 . Since (25) provides an approximated value of the first derivative, this is only valid during a period of time. So the state estimation should be calculated periodically as follows,

$$12(t - t\_i)^3 z + 96 \int\_{t\_i}^t \int\_{t\_i}^{t\_1} \lambda\_2 z(\lambda\_2) d\lambda\_2 d\lambda\_1$$

$$\frac{dzl}{dt}\Big|\_{t\_i}^t = \frac{-24 \int\_{t\_i}^t \int\_{t\_i}^{t\_1} \int\_{t\_i}^{t\_2} z(\lambda\_3) d\lambda\_3 d\lambda\_2 d\lambda\_1}{\left(t - t\_i\right)^4}, \ \forall (t - t\_i) \ge 0 \tag{26}$$

where ( ,) *<sup>i</sup> t t* is the estimation period.

432 MATLAB for Engineers – Applications in Control, Electrical Engineering, IT and Robotics

The algebraic methodology proposed in (Fliess & Sira-Ramírez, 2003) allows us to estimate the derivatives of a smooth signal considering a signal model of *n th* − order, thus it is not necessary to design the time derivative estimator from a specific dynamic model of the plant. Through valid algebraic manipulations of this approximated model in the frequency domain, and using the algebraic derivation with respect to the complex variable *s* , we neglect the initial conditions of the signal. The resulting equation is multiplied by a negative power *<sup>n</sup>*−<sup>1</sup> *s* and returned to the time domain. This last expression is manipulated

**5. On-line algebraic state estimation of active suspension system** 

algebraically in order to find a formula to estimate the first time derivative of *y*( )*t* .

4 4 ( ) = 0 *dyt*

This model indicates that *y*( )*t* is a signal whose behavior can be approximated by a family of polynomials of third order, thus the fourth time derivative is assumed close to zero. The order of this approximation can be increased to enhance the estimation quality of the algebraic estimator. From (21) it is possible to design a time derivative algebraic estimator.

Now, taking successive derivatives until a number of three with respect to the complex

( ) 4 4 <sup>4</sup> = 0

2 34 24 96 72 16 −− − + + ++ *dY d Y d Y d Y sY s s <sup>s</sup>*

1 1 2

3 1 1 2

12 96 ( ) 24 ( )

*t t tz z d d z ddd dz*

*t*

<sup>111</sup> <sup>0</sup>

λλλ

2 2 21 3 321 0 0 00 0

4

96 ( ) 24 ( ) = 0

*z dd z ddd*

λ λ

2 2 21 3 321 0 0 00 0

*d sY*

*dt* (21)

*ds* (23)

*ds ds ds ds* (24)

<sup>432</sup> (3) *s Y s s Y s Y sY Y* ( ) (0) (0) = 0 − − −− (22)

2 34

λ λλλ

λ λ

λ λλλ

(25)

Consider a fourth order approximation of a smooth signal *y*( )*t* ,

variable *s* , we obtain a expression which is free of initial conditions,

32 1

43 2

− +

− +

λ

λ

λ λ λλ

*dt t*

( ( )) 16 ( ) 72 ( )

*<sup>d</sup> tzt tzt z d*

*t t*

From the last equation is possible to obtain the following algebraic estimator,

+ − λ λ λλ

**5.1 First time derivative algebraic estimation** 

First, we apply Laplace transform to (21),

Expanding this expression and multiplying by <sup>−</sup><sup>3</sup> *s* ,

Returning to the time domain,

*dt*

=

In order to obtain a better and smoother estimated value of the vertical velocity, we have implemented simultaneously two algebraic estimators for each velocity to estimate. Proceeding with an out-of-phase policy for one of these algebraic estimators, the outputs of both are combined properly to obtain a final estimated value.

#### **6. Instrumentation of the active suspension system**

The only variables required for the implementation of the proposed controllers are the vertical displacements of the body of the car *<sup>s</sup> z* and the vertical displacement of the wheel *<sup>u</sup> z* . These variables are needed to measure by some sensor. In (Chamseddine et al., 2006) the use of sensors in experimental vehicle platforms, as well as in commercial vehicles is presented. The most common sensors used for measuring the vertical displacement of the body of the car and the wheel are laser sensors. This type of sensors could be used to measure the variables *<sup>s</sup> z* and *<sup>s</sup> z* needed for controller implementation. An accelerometer or another type of sensor is not needed to measure the variables *<sup>s</sup> z* and *<sup>u</sup> z* , these variables are estimated with the use of algebraic estimators from knowledge of the variables *<sup>s</sup> z* and *<sup>u</sup> z* . Fig. 2 shows a schematic diagram of the instrumentation for the active suspension system.

Fig. 2. Schematic diagram of the instrumentation of the active suspension system.

Robust Control of Active Vehicle Suspension Systems

with *p* = 100 ,

differential flatness controller with algebraic state estimation.

ζ

Fig. 4. Simulink model of the sliding mode and differential flatness controller.

hydraulic suspension controller are depicted in Fig. 7.

those gotten by the controllers using the real velocities.

In Fig. 6 is depicted the robust performance of the electromagnetic suspension controller. It can be seen the high vibration attenuation level of the active vehicle suspension system compared with the passive counterpart. Similar results on the implementation of the

In Fig. 8 is presented the algebraic estimation process of the velocities of the car body and the wheel. There we can observe a good and fast estimation. In Figs. 9 and 10 are shown the simulation results on the performance of the electromagnetic and hydraulic suspension controllers using the algebraic estimators of velocities. These results are quite similar to

2 2 ( )( 2 ) *n n s ps s* ++ + ζω

 ω

Using Sliding Modes and Differential Flatness with MATLAB 435

It is desired to stabilize the suspension system at the positions = 0 *<sup>s</sup> z* and = 0 *<sup>u</sup> z* . The gains of both electromagnetic and hydraulic suspension controllers were obtained by forcing their closed loop characteristic polynomials to be given by the following Hurwitz polynomial:

> = 0.5 , = 90 ω*<sup>n</sup>* ,

The Simulink model of the sliding mode and differential flatness controller of the active suspension system is shown in Fig. 4. For the electromagnetic active suspension system, it is assumed that *cz = 0*. In Fig. 5 is shown the Simulink model of the sliding mode and

μ= 95 y γ= 95 .

#### **7. Simulation results**

The simulation results were obtained by means of MATLAB/Simulink ® , with the Runge-Kutta numerical method and a fixed integration step of 1*ms* . The numerical values of the quarter-car model parameters (Sam & Hudha, 2006) are presented in Table 1.


Table 1. Quarter-car model parameters

In this simulation study, the road disturbance is shown in Fig. 3 and set in the form of (Sam & Hudha, 2006):

$$z\_r = a \frac{1 - \cos(8\pi t)}{2}$$

with *a* = 0.11 [m] for 0.5 0.75 ≤ ≤*t* , *a* = 0.55 [m] for 3.0 3.25 ≤ ≤*t* and 0 otherwise.

Fig. 3. Type of road disturbance.

The simulation results were obtained by means of MATLAB/Simulink ® , with the Runge-Kutta numerical method and a fixed integration step of 1*ms* . The numerical values of the

> Parameter Value Sprung mass ( *ms* ) 282 [ ] *kg* Unsprung mass ( *mu* ) 45 [ ] *kg*

Spring stiffness ( *<sup>s</sup> <sup>k</sup>* ) 17900 *<sup>N</sup>*

Damping constant ( *<sup>s</sup> <sup>c</sup>* ) <sup>1000</sup> *N s*

Tire stiffness ( *<sup>t</sup> <sup>k</sup>* ) 165790 *<sup>N</sup>*

In this simulation study, the road disturbance is shown in Fig. 3 and set in the form of (Sam

2 *<sup>r</sup>*

=

*z a* <sup>−</sup>

with *a* = 0.11 [m] for 0.5 0.75 ≤ ≤*t* , *a* = 0.55 [m] for 3.0 3.25 ≤ ≤*t* and 0 otherwise.

1 (8 )

*cos t*

π

*m* 

*m* <sup>⋅</sup>

> *m*

quarter-car model parameters (Sam & Hudha, 2006) are presented in Table 1.

**7. Simulation results** 

Table 1. Quarter-car model parameters

Fig. 3. Type of road disturbance.

& Hudha, 2006):

It is desired to stabilize the suspension system at the positions = 0 *<sup>s</sup> z* and = 0 *<sup>u</sup> z* . The gains of both electromagnetic and hydraulic suspension controllers were obtained by forcing their closed loop characteristic polynomials to be given by the following Hurwitz polynomial: 2 2 ( )( 2 ) *n n s ps s* ++ + ζω ω with *p* = 100 , ζ = 0.5 , = 90 ω*<sup>n</sup>* , μ = 95 y γ= 95 .

The Simulink model of the sliding mode and differential flatness controller of the active suspension system is shown in Fig. 4. For the electromagnetic active suspension system, it is assumed that *cz = 0*. In Fig. 5 is shown the Simulink model of the sliding mode and differential flatness controller with algebraic state estimation.

Fig. 4. Simulink model of the sliding mode and differential flatness controller.

In Fig. 6 is depicted the robust performance of the electromagnetic suspension controller. It can be seen the high vibration attenuation level of the active vehicle suspension system compared with the passive counterpart. Similar results on the implementation of the hydraulic suspension controller are depicted in Fig. 7.

In Fig. 8 is presented the algebraic estimation process of the velocities of the car body and the wheel. There we can observe a good and fast estimation. In Figs. 9 and 10 are shown the simulation results on the performance of the electromagnetic and hydraulic suspension controllers using the algebraic estimators of velocities. These results are quite similar to those gotten by the controllers using the real velocities.

Robust Control of Active Vehicle Suspension Systems

differential flatness based controller.

differential flatness based controller.

Using Sliding Modes and Differential Flatness with MATLAB 437

Fig. 6. Electromagnetic active vehicle suspension system responses with sliding mode and

Fig. 7. Hydraulic active vehicle suspension system responses with sliding mode and

Fig. 5. Simulink model of the sliding mode and differential flatness controller with state estimation.

Fig. 5. Simulink model of the sliding mode and differential flatness controller with state

estimation.

Fig. 6. Electromagnetic active vehicle suspension system responses with sliding mode and differential flatness based controller.

Fig. 7. Hydraulic active vehicle suspension system responses with sliding mode and differential flatness based controller.

Robust Control of Active Vehicle Suspension Systems

Using Sliding Modes and Differential Flatness with MATLAB 439

Fig. 9. b. Electromagnetic active vehicle suspension system responses with sliding mode and

Fig. 10. a. Hydraulic active vehicle suspension system responses with sliding mode and

differential flatness based controller using algebraic state estimation.

differential flatness based controller using algebraic state estimation.

Fig. 8. On-line algebraic state estimates of the hydraulic active suspension system.

Fig. 9. a. Electromagnetic active vehicle suspension system responses with sliding mode and differential flatness based controller using algebraic state estimation.

Fig. 8. On-line algebraic state estimates of the hydraulic active suspension system.

Fig. 9. a. Electromagnetic active vehicle suspension system responses with sliding mode and

differential flatness based controller using algebraic state estimation.

Fig. 9. b. Electromagnetic active vehicle suspension system responses with sliding mode and differential flatness based controller using algebraic state estimation.

Fig. 10. a. Hydraulic active vehicle suspension system responses with sliding mode and differential flatness based controller using algebraic state estimation.

Robust Control of Active Vehicle Suspension Systems

Germany, October 4-6, 2006.

IFAC Nolcos Conference. Germany, 2004.

CINVESTAV-IPN, México, D.F., México, 2005.

pp. 123-125, Harbin, Heilongjiang, August 7-11, 2006.

an Active Suspension System", IEEE ICIEA, 2006.

Limassol, Cyprus, June 27-29, 2005.

Vehicular Technology. Vol. 55, No. 1, pp. 86-94, January 2006.

**9. References** 

8-10, 2009.

December 2000.

1361, 1995.

96, 2004.

2003.

620, April 1993.

2004.

Using Sliding Modes and Differential Flatness with MATLAB 441

Chamseddine, Abbas; Noura, Hassan; Raharijaona, Thibaut "Control of Linear Full

Chávez-Conde, E.; Beltrán-Carbajal, F.; Blanco-Ortega, A.; Méndez-Azúa, H. "Sliding

Enríquez-Zárate, J.; Silva-Navarro, G.; Sira-Ramírez, H. "Sliding Mode Control of a

Flies, M.; Lévine, J.; Martin, Ph.; Rouchon, P. "Flatness and defect of nonlinear systems:

Fliess, M.; Sira-Ramírez, H. "Reconstructeurs d'états". C.R. Acad. Sci. Paris, I-338, pp. 91-

Fliess, M.; Sira-Ramírez, H. "Control via state estimations of some nonlinear systems", 4th

Fliess, M.; Sira-Ramírez, H. "An Algebraic Framework For Linear Identification", ESAIM

García-Rodríguez, C. "Estimación de estados por métodos algebraicos", Master of Thesis.

Liu, Zhen; Luo, Cheng; Hu, Dewen, "Active Suspension Control Design Using a

Martins, I.; Esteves, J.; Marques, D. G.; Da Silva, F. P. "Permanent-Magnets Linear

Sam, Y. M.; Hudha, K. "Modeling and Force Tracking Control of Hydraulic Actuator for

Sira-Ramírez, H. "A dynamical variable structure control strategy in asymptotic output

Sira-Ramírez, H.; Silva-Navarro, G. "Algebraic Methods in Flatness, Signal Processing and State Estimation", Innovación Editorial Lagares de México, November 2003. Sira-Ramírez, H.; Agrawal, Sunil K. "Differentially Flat Systems", Marcel Dekker, N.Y.,

Tahboub, Karim A. "Active Nonlinear Vehicle-Suspension Variable-Gain Control", 13th

Vehicle Active Suspension System Using Sliding Mode Techniques", 2006 IEEE International Conference on Control Applications. pp. 1306-1311, Munich,

Mode and Generalized PI Control of Vehicle Active Suspensions", 2009 IEEE International Conference on Control Applications. Saint Petersburg, Russia, July

Differentially Flat Vibrational Mechanical System: Experimental Results", 39th IEEE Conference on Decision and Control. pp. 1679-1684, Sydney, Australia,

introductory theory and examples", Int. Journal of Control. Vol. 61, pp. 1327-

Control, Optimization and Calculus of Variations. Vol 9, pp. 151-168, January.

Combination of LQR and Backstepping", 25th IEEE Chinese Control Conference,

Actuators Applicability in Automobile Active Suspensions", IEEE Trans. on

tracking problems", IEEE Trans. on Automatic Control. Vol. 38, No. 4, pp. 615-

IEEE Mediterranean Conference on Control and Automation, pp. 569-574,

Fig. 10. b. Hydraulic active vehicle suspension system responses with sliding mode and differential flatness based controller using algebraic state estimation.

#### **8. Conclusions**

The stabilization of the vertical position of the quarter of car is obtained in a time much smaller to that of the passive suspension system. The sliding mode based differential flatness controller requires the knowledge of all the state variables. Nevertheless the fast stabilization with amplitude in acceleration and speed of the body of the car very remarkable is observed. On-line state estimation is obtained successfully, however when it is used into the controller one can observe a deterioration of the control signal. This can significantly improve with a suitable interpolation between the estimated values at each restart of the integrations. In addition, the simulations results show that the stabilization of the system is obtained before the response of the passive suspension system, with amplitude of acceleration and speed of the body of the car very remarkable. Finally, the robustness of the controllers is observed to take to stabilize to the system before the unknown disturbance.

#### **9. References**

440 MATLAB for Engineers – Applications in Control, Electrical Engineering, IT and Robotics

Fig. 10. b. Hydraulic active vehicle suspension system responses with sliding mode and

The stabilization of the vertical position of the quarter of car is obtained in a time much smaller to that of the passive suspension system. The sliding mode based differential flatness controller requires the knowledge of all the state variables. Nevertheless the fast stabilization with amplitude in acceleration and speed of the body of the car very remarkable is observed. On-line state estimation is obtained successfully, however when it is used into the controller one can observe a deterioration of the control signal. This can significantly improve with a suitable interpolation between the estimated values at each restart of the integrations. In addition, the simulations results show that the stabilization of the system is obtained before the response of the passive suspension system, with amplitude of acceleration and speed of the body of the car very remarkable. Finally, the robustness of the controllers is observed to take to stabilize to the system before the

differential flatness based controller using algebraic state estimation.

**8. Conclusions** 

unknown disturbance.


**19** 

*Tunisia* 

**Thermal Behavior of IGBT** 

Mohamed Amine Fakhfakh, Moez Ayadi,

*F fM g* <sup>1</sup> = *<sup>v</sup>* (1)

Ibrahim Ben Salah and Rafik Neji

 *University of Sfax/Sfax* 

**Module for EV (Electric Vehicle)** 

EVs are divided into three categories: the pure EV, the hybrid EV, and the fuel cell EV. Although these three types of electric vehicle have different system configuration, one (or more) motor drive system is always needed to convert electrical power into mechanical ones. Among the drive systems used for EV, induction motor system and permanent magnet motor systems are mostly used for their high power density, high

The motor drive system for electric vehicle (EV) is composed of a battery, three phase inverter, a permanent magnet motor, and a sensor system. The inverter is a key unit important among these electrical components which converts the direct current of the battery into the alternating current to rotate the motor. Therefore, for predicting the dynamic power loss and junction temperature, the electro-thermal coupling simulation techniques to estimate the power loss and to calculate the junction temperature become

This paper describes a compact thermal model suitable for the electro-thermal coupling simulation of EV inverter module for two current control methods. We can predict the dynamic temperature rise of Si devices by simulating the inverter operation in accordance

As shown in Figure 1 and table 1, there are six forces acting on the electric vehicle: the rolling resistance force, the aerodynamic force, the aerodynamic lift force, the gravity force, the

Rolling resistance is due the tires deforming when contacting the surface of a road and varies depending on the surface being driven on. It can be model using the following

**1. Introduction** 

efficiency.

important.

equation:

with the real EV running.

**2. Dynamic model of the EV** 

normal force, and the motor force.

**2.1 Rolling resistance force** 


## **Thermal Behavior of IGBT Module for EV (Electric Vehicle)**

Mohamed Amine Fakhfakh, Moez Ayadi, Ibrahim Ben Salah and Rafik Neji  *University of Sfax/Sfax Tunisia* 

#### **1. Introduction**

442 MATLAB for Engineers – Applications in Control, Electrical Engineering, IT and Robotics

Utkin,V. I. "Sliding Modes and Their Applications in Variable Structure Systems".

Yousefi, A.; Akbari, A and Lohmann, B., "Low Order Robust Controllers for Active

pp. 693-698, Munich, Germany, Octuber 4-6, 2006.

Vehicle Suspensions", IEEE International Conference on Control Applications,

Moscow: MIR, 1978.

EVs are divided into three categories: the pure EV, the hybrid EV, and the fuel cell EV. Although these three types of electric vehicle have different system configuration, one (or more) motor drive system is always needed to convert electrical power into mechanical ones. Among the drive systems used for EV, induction motor system and permanent magnet motor systems are mostly used for their high power density, high efficiency.

The motor drive system for electric vehicle (EV) is composed of a battery, three phase inverter, a permanent magnet motor, and a sensor system. The inverter is a key unit important among these electrical components which converts the direct current of the battery into the alternating current to rotate the motor. Therefore, for predicting the dynamic power loss and junction temperature, the electro-thermal coupling simulation techniques to estimate the power loss and to calculate the junction temperature become important.

This paper describes a compact thermal model suitable for the electro-thermal coupling simulation of EV inverter module for two current control methods. We can predict the dynamic temperature rise of Si devices by simulating the inverter operation in accordance with the real EV running.

#### **2. Dynamic model of the EV**

As shown in Figure 1 and table 1, there are six forces acting on the electric vehicle: the rolling resistance force, the aerodynamic force, the aerodynamic lift force, the gravity force, the normal force, and the motor force.

#### **2.1 Rolling resistance force**

Rolling resistance is due the tires deforming when contacting the surface of a road and varies depending on the surface being driven on. It can be model using the following equation:

$$F\_1 = f \, M\_v \, \mathcal{g} \tag{1}$$

Thermal Behavior of IGBT Module for EV (Electric Vehicle) 445

*mv a p r dV F M FFF dt*

The power that the EV must develop at stabilized speed is expressed by the following

( ) *vehicle m r a p v dV P FV F F F M V*

We deduce the expression of the total torque by multiplying equation (5) with the wheel

Neglecting the mechanical losses in the gearbox, the t electromagnetic torque Cem developed by the motor is obtained by dividing the wheels torque Cvehicle by the ratio reduction rd.

*vehicle r a p v dV C C C C MR*

*em r a p v*

*dV C C C C MR r dt*

Figure 2 presents the dynamic model of the EV load, implemented under Matlab/simulink.

rm <sup>1</sup>

1/(Mv\*Rwheels)

Gain2

Control of permanent magnet synchronous motor is performed using field oriented control. The stator windings of the motor are fed by an inverter that generates a variable frequency variable voltage. The frequency and phase of the output wave are controlled using a

In our studie, we have used two types of current control, Hysteresis and PWM.

1

*d*

r

Rwheels Gain1

Fig. 2. SIMULINK dynamic model of electric vehicle

**3. Electric motor control** 

position sensor as shown in figure 3.

Mv\*g\*sin(u) F3

0.5\*1.202\*Sf\*Cx\*u^2 F2

f\*Mv\*g F1

= +++ (5)

*dt* = = +++ (6)

*dt*

=+++ (7)

= +++ (8)

dV/dt V(m/s)

r/Rwheels

Gain4 3.6 Gain3

s Integrator

3 Wm

2 V(km/h)

By projection on the (O, x) axis, we obtain:

equation:

radius R:

1 Cr

1 Cem

2 pte\_alpha

Fig. 1. Diagram of forces applied to the EV


Table 1. Applied forces to EV

#### **2.2 Aerodynamic force**

Aerodynamic drag is caused by the momentum loss of air particles as they flow over the hood of the vehicle. The aerodynamic drag of a vehicle can be modeled using the following equation:

$$\mathcal{F}\_2 = \frac{1}{2} \rho \mathcal{S}\_f \mathcal{C}\_\mathbf{x} V^2 \tag{2}$$

#### **2.3 Gravity force**

The gravity force can be calculated as follows:

$$F\_3 = M\_v \, g \sin \theta \,\tag{3}$$

#### **2.4 Motor force**

Using Newton's Second Law, we can deduce the motor force; it can be obtained by the following equation:

$$M\_v \frac{dV}{dt} = \sum \overrightarrow{F\_{ext}} = \overrightarrow{F\_m} + \overrightarrow{F\_p} + \overrightarrow{F\_a} + \overrightarrow{F\_r} \tag{4}$$

By projection on the (O, x) axis, we obtain:

444 MATLAB for Engineers – Applications in Control, Electrical Engineering, IT and Robotics

*Fm*

*Fpy*

*Fpx*

*Far*

*Fprop*

*Fp*

*Fr2x*

Fr1x Rolling resistance force Fr2x Rolling resistance force

θ Slope angle with the horizontal

Aerodynamic drag is caused by the momentum loss of air particles as they flow over the hood of the vehicle. The aerodynamic drag of a vehicle can be modeled using the following

<sup>3</sup>*F M*= *<sup>v</sup> g*sin

Using Newton's Second Law, we can deduce the motor force; it can be obtained by the

*v ext m p a r dV <sup>M</sup> F F FFF*

2 1 2 *F SCV* = ρ

*dt*

2

θ

→ → →→→

*f x* (2)

= = +++ (4)

(3)

Fav Normal force Far Normal force Fa Aerodynamic force Fprop Thrust force Fp Gravity force Fm Motor force

*Fa*

Fig. 1. Diagram of forces applied to the EV

*Fr1x*

*Fav*

*Fprop*

Table 1. Applied forces to EV

The gravity force can be calculated as follows:

**2.2 Aerodynamic force** 

equation:

**2.3 Gravity force** 

**2.4 Motor force** 

following equation:

$$F\_m = M\_v \frac{dV}{dt} + F\_a + F\_p + F\_r \tag{5}$$

The power that the EV must develop at stabilized speed is expressed by the following equation:

$$P\_{\text{vehicle}} = F\_m V = \left(F\_r + F\_a + F\_p + M\_v \frac{dV}{dt}\right) V \tag{6}$$

We deduce the expression of the total torque by multiplying equation (5) with the wheel radius R:

$$\mathbf{C}\_{\text{vehicle}} = \mathbf{C}\_r + \mathbf{C}\_a + \mathbf{C}\_p + M\_v R \frac{dV}{dt} \tag{7}$$

Neglecting the mechanical losses in the gearbox, the t electromagnetic torque Cem developed by the motor is obtained by dividing the wheels torque Cvehicle by the ratio reduction rd.

$$\mathbf{C}\_{em} = \frac{1}{r\_d} \left( \mathbf{C}\_r + \mathbf{C}\_a + \mathbf{C}\_p + M\_v \, R \frac{dV}{dt} \right) \tag{8}$$

Figure 2 presents the dynamic model of the EV load, implemented under Matlab/simulink.

Fig. 2. SIMULINK dynamic model of electric vehicle

#### **3. Electric motor control**

Control of permanent magnet synchronous motor is performed using field oriented control. The stator windings of the motor are fed by an inverter that generates a variable frequency variable voltage. The frequency and phase of the output wave are controlled using a position sensor as shown in figure 3.

In our studie, we have used two types of current control, Hysteresis and PWM.

Thermal Behavior of IGBT Module for EV (Electric Vehicle) 447

The studied module is the Semikron module SKM 75GB 123D (75A/1200V) which contains two IGBTs and with two antiparallel diodes. The structure of the module contains primarily eight layers of different materials, each one of it is characterized by its thickness Li, its thermal conductivity Ki, density ρi and its heat capacity Cpi. Table 2 show the materials properties of the various layers of module as shown in figure 6. These values are given by the manufacturer and/or of the literatures [Dorkel et al., 1996; Uta et

> Material *L (mm) K (W/mK) ρCp (J/Kcm3)*  Silicium 0.4 140 1.7 Solder 1 0.053 35 1.3 Copper 0.35 360 3.5 Isolation 0.636 100 2.3 Copper 0.35 360 3.5 Solder 2 0.103 35 1.3 Base plate 3 280 3.6 Grease 0.1 1 2.1

Solder 2

Copper Isolation

Copper

Grease

Base plate

Solder 1

Silicium

In the power module, the heating flow diffuses vertically and also laterally from the heating source. So, a thermal interaction happens inside the module between the adjacent devices

This thermal interaction depends from [Kojima et al., 2006; Ayadi et al., 2010; Fakhfakh et

Figure 7 shows the thermal influence between the different components of the module. We notice that each component has a thermal interaction with the others and we supposed that

**4. Thermal model of IGBT module** 

al., 2000; Thoams et al., 2000].

Fig. 6. Example of the module structure

Table 2. Thermal parameters of a power module



each module have zero interaction with other modules.

when they operate together.

al., 2010]:

Fig. 3. Drive system schematic

#### **3.1 PWM current controller**

PWM current controllers are widely used. The switching frequency is usually kept constant. They are based in the principle of comparing a triangular carrier wave of desire switching frequency and is compared with error of the controlled signal [Bose, 1996].

Fig. 4. PWM current controller

#### **3.2 Hysteresis current controller**

Hysteresis current controller can also be implemented to control the inverter currents. The controller will generate the reference currents with the inverter within a range which is fixed by the width of the band gap [Bose, 1996; Pillay et al., 1989].

Fig. 5. Hysteresis current controller

#### **4. Thermal model of IGBT module**

446 MATLAB for Engineers – Applications in Control, Electrical Engineering, IT and Robotics

PWM current controllers are widely used. The switching frequency is usually kept constant. They are based in the principle of comparing a triangular carrier wave of desire switching

Hysteresis current controller can also be implemented to control the inverter currents. The controller will generate the reference currents with the inverter within a range which is fixed

frequency and is compared with error of the controlled signal [Bose, 1996].

Fig. 3. Drive system schematic

**3.1 PWM current controller** 

Fig. 4. PWM current controller

**3.2 Hysteresis current controller** 

Fig. 5. Hysteresis current controller

by the width of the band gap [Bose, 1996; Pillay et al., 1989].

The studied module is the Semikron module SKM 75GB 123D (75A/1200V) which contains two IGBTs and with two antiparallel diodes. The structure of the module contains primarily eight layers of different materials, each one of it is characterized by its thickness Li, its thermal conductivity Ki, density ρi and its heat capacity Cpi. Table 2 show the materials properties of the various layers of module as shown in figure 6. These values are given by the manufacturer and/or of the literatures [Dorkel et al., 1996; Uta et al., 2000; Thoams et al., 2000].

Fig. 6. Example of the module structure


Table 2. Thermal parameters of a power module

In the power module, the heating flow diffuses vertically and also laterally from the heating source. So, a thermal interaction happens inside the module between the adjacent devices when they operate together.

This thermal interaction depends from [Kojima et al., 2006; Ayadi et al., 2010; Fakhfakh et al., 2010]:


Figure 7 shows the thermal influence between the different components of the module. We notice that each component has a thermal interaction with the others and we supposed that each module have zero interaction with other modules.

Thermal Behavior of IGBT Module for EV (Electric Vehicle) 449

are introduced between solder 2 and base plate because all module components have the base plate as a common material. So the thermal circuit network of IGBT1 becomes as the

P1 P2 P3

The PM motor drive simulation was built in several steps like abc phase transformation to dqo variables, calculation torque and speed, and control circuit [Ong, 1998; Roisse et

Silicon Copper Isolating Base plate Heat-Sink

Parks transformation used for converting Iabc to Idq is shown in figure 10 and the reverse

The inverter is implemented in Simulink as shown in figure 12. The inverter consists of the "universal bridge" with the parameters of the IGBT module studied. All the voltages and the currents in the motor and the inverter can be deducted. The following figure shows the

(2/3)\*(sin(u(4))\*u(1)+sin(u(4)-2\*pi/3)\*u(2)+sin(u(4)-4\*pi/3)\*u(3)) Fcn1

(2/3)\*(cos(u(4))\*u(1)+cos(u(4)-2\*pi/3)\*u(2)+cos(u(4)-4\*pi/3)\*u(3)) Fcn

> 2 Iq

> 1 Id

For proper control of the inverter using the reference currents, current controllers are implemented generate the gate pulses for the IGBT's. Current controllers used are shown in

transformation for converting Idq to Iabc is shown in figure 11.

figure 9.

Fig. 9. Thermal model of IGBT module

**5. Simulation and results** 

Fig. 10. Iabc to Idq bloc

model of the inverter used.

figure 13 and 14.

al., 1998].

4 theta

3 Ic

2 Ib

1 Ia

Fig. 7. Different thermal influences between the module components

Literature proposes some thermal circuit networks for electrothermal simulation for the semiconductor device. For example the finite difference method (FDM) and the finite element method (FEM). In our study we have used the FEM technique to model our inverter module. Figure 8 shows the thermal circuit example obtained by the FEM of IGBT1 without thermal interaction.

Fig. 8. Thermal circuit obtained by the FEM

Where:


In order to introduce the thermal interaction between the different components of the module, we inserted three other current sources P1, P2 and P3. These sources are deduced from the structure of IGBT module [Drofenik et al., 2005; Hamada et al., 2006; Usui et al., 2006].

The source P1 is the power loss of DIODE1; it is introduced at the interface between the silicon and the copper materials because the IGBT1 and the DIODE1 ships are bounded on the same copper area. The source P2 and P3 are power loss of IGBT2 and DIODE2, they

DIODE 1

DIODE 2

Literature proposes some thermal circuit networks for electrothermal simulation for the semiconductor device. For example the finite difference method (FDM) and the finite element method (FEM). In our study we have used the FEM technique to model our inverter module. Figure 8 shows the thermal circuit example obtained by the FEM of IGBT1 without

In order to introduce the thermal interaction between the different components of the module, we inserted three other current sources P1, P2 and P3. These sources are deduced from the structure of IGBT module [Drofenik et al., 2005; Hamada et al., 2006; Usui et

The source P1 is the power loss of DIODE1; it is introduced at the interface between the silicon and the copper materials because the IGBT1 and the DIODE1 ships are bounded on the same copper area. The source P2 and P3 are power loss of IGBT2 and DIODE2, they

Fig. 7. Different thermal influences between the module components

IGBT 1

IGBT 2

thermal interaction.

Where:

al., 2006].

Fig. 8. Thermal circuit obtained by the FEM



are introduced between solder 2 and base plate because all module components have the base plate as a common material. So the thermal circuit network of IGBT1 becomes as the figure 9.

Fig. 9. Thermal model of IGBT module

#### **5. Simulation and results**

The PM motor drive simulation was built in several steps like abc phase transformation to dqo variables, calculation torque and speed, and control circuit [Ong, 1998; Roisse et al., 1998].

Parks transformation used for converting Iabc to Idq is shown in figure 10 and the reverse transformation for converting Idq to Iabc is shown in figure 11.

Fig. 10. Iabc to Idq bloc

The inverter is implemented in Simulink as shown in figure 12. The inverter consists of the "universal bridge" with the parameters of the IGBT module studied. All the voltages and the currents in the motor and the inverter can be deducted. The following figure shows the model of the inverter used.

For proper control of the inverter using the reference currents, current controllers are implemented generate the gate pulses for the IGBT's. Current controllers used are shown in figure 13 and 14.

Thermal Behavior of IGBT Module for EV (Electric Vehicle) 451

NOT

opérateur logic

1 MLI

double

Cem

V(km/h)

w

pte\_alpha

0 Constant

modèle dynamique du VE

The complete system used for simulation and implemented in MATLAB / Simulink, is shown in Figure 15. This system was tested with two current controls, hysteresis and PWM control. The motor used is an axial flux Permanent Magnet Synchronous Motor (PMSM).

Tem

teta

Idq

onduleur+moteur

NOT

NOT

opérateur logic1

opérateur logic2

For the simulation, we controlled the speed of EV at 30km / h.

g

oolea

oolea

oolea

Fig. 14. Hysteresis controller

Relay2

3 Xc

2 Xb

1 Xa

Relay1

Relay

Fig. 15. PMSM in a traction chain

vitesse <sup>w</sup>

Pulses

Vref

V

idq

the

commande

Fig. 11. Idq to Iabc bloc

Fig. 12. Inverter model

Fig. 13. PWM current controller

Fig. 14. Hysteresis controller

g A B C

Universal Bridge

+


oolea

oolea

oolea

(cos(u(3)-2\*pi/3)\*u(1) +sin(u(3)-2\*pi/3)\*u(2)) Fcn2

(cos(u(3))\*u(1) +sin(u(3))\*u(2)) Fcn

(cos(u(3)-4\*pi/3)\*u(1) +sin(u(3)-4\*pi/3)\*u(2)) Fcn1

> 1 signal de commande

opérateur logic2

NOT opérateur logic1

NOT opérateur logic

3 Ic

> 1 MLI

double

2 Ib

1 Ia

Fig. 11. Idq to Iabc bloc

3 teta

2 Iq

1 Id

Fig. 12. Inverter model

3 Xc

2 Xb

1 Xa Batterie

Fig. 13. PWM current controller

signal triangulaire NOT

>=

>=

>=

The complete system used for simulation and implemented in MATLAB / Simulink, is shown in Figure 15. This system was tested with two current controls, hysteresis and PWM control. The motor used is an axial flux Permanent Magnet Synchronous Motor (PMSM). For the simulation, we controlled the speed of EV at 30km / h.

Fig. 15. PMSM in a traction chain

Thermal Behavior of IGBT Module for EV (Electric Vehicle) 453

0 1 2 3 4 5 6 7

Temps (s)

Time (s)

IGBT1 DIODE1

3.015 3.02 3.025 3.03 3.035 3.04

Temps (s)

Time (s)

Fig. 18. Iabc currents with PWM control


Température de jonction (°C)

Power losses (W)




0

Courant (A)

Current (A)

20

40

60

80

Fig. 19. IGBT1 and DIODE1 power losses with PWM control

Figure 16 shows the EV speed regulated at 30km / h for the two types of control. We note that with the hysteresis control, we reach faster the steady state.

Fig. 16. EV speed; (1): with PWM controller; (2): with hysteresis controller

The stator phase currents corresponding to this regulation are represented by figure 17 and 18 Figure 19 and 20 show the IGBT1 and DIODE1 power losses for hysteresis and PWM current control respectively.

Fig. 17. Iabc currents with hysteresis control

Figure 16 shows the EV speed regulated at 30km / h for the two types of control. We note

(2)

(1)

that with the hysteresis control, we reach faster the steady state.

Consigne

Fig. 16. EV speed; (1): with PWM controller; (2): with hysteresis controller

current control respectively.

0

5

10

15

Speed (km/h)

20

25

30

35

Fig. 17. Iabc currents with hysteresis control





0

Courant (A)

Current (A)

20

40

60

80

The stator phase currents corresponding to this regulation are represented by figure 17 and 18 Figure 19 and 20 show the IGBT1 and DIODE1 power losses for hysteresis and PWM

0 1 2 3 4 5 6 7

Temps (s)

Time (s)

0 1 2 3 4 5 6 7

Time (s)

Fig. 18. Iabc currents with PWM control

Fig. 19. IGBT1 and DIODE1 power losses with PWM control

Thermal Behavior of IGBT Module for EV (Electric Vehicle) 455

A detailed dynamic model for EV was studied using two current control systems. MATLAB / Simulink were chosen from several simulation tools because of its flexibility in working with analog and digital devices, it is able to represent real-time results with the simulation time reduced. A comparative study was carried out in terms of switching frequency for power dissipated by the components of the inverter and junction temperature. The hysteresis current control has a variable switching frequency that depends on the hysteresis band, this type of control allows for fast simulations with a shorter time. The PWM current control has a fixed frequency switching and allows having junction temperatures lower than

0 1 2 3 4 5 6 7

Temps (s)

Time (s)

B. K. Bose, Power Electronics and Variable Frequency Drives. (1996). 1 ed: Wiley, John &

P. Pillay & R. Krishnan. (1989). Modeling, simulation, and analysis of permanent-magnet

Jean-Marie Dorkel, Patrick Tounsi, & Philippe Leturcq. (1996). Three-Dimensional thermal

Uta Hecht & Uwe Scheuermann. (2000). Static and Transient Thermal Resistance of

*Applications, IEEE Transactions on*, vol. 25, pp. 265-273

motor drives. I. The permanent-magnet synchronous motor drive. *Industry* 

Modeling Based on the Two-Port Network Theory for Hybrid or Monolithic Integrated Power Circuits. IEEE Transaction on Electronics Devices, vol. 19, NO. 4,

Advanced Power Modules. *Semikron Elektronik GmbH*, Sigmundstr. 200, 90431

Fig. 22. DIODE1 junction temperature

32

34

36

38

Température de jonction (°C)

Junction tem

perature (°C)

40

42

44

46

**6. Conclusion** 

the hysteresis control.

Sons

pp. 501-507

Nürnberg (Germany).

**7. References** 

Fig. 20. IGBT1 and DIODE1 power losses with hysteresis control

Figure 21 and 22 show the IGBT1and DIODE1 junction temperature obtained by the two types of current control. It is very clear that the junction temperature of IGBT1 and DIODE1 is higher for the hysteresis control; this is due by the increase of power dissipation of the module components this type of control.

Fig. 21. IGBT1 junction temperature

Fig. 22. DIODE1 junction temperature

#### **6. Conclusion**

454 MATLAB for Engineers – Applications in Control, Electrical Engineering, IT and Robotics

IGBT1 DIODE1

Figure 21 and 22 show the IGBT1and DIODE1 junction temperature obtained by the two types of current control. It is very clear that the junction temperature of IGBT1 and DIODE1 is higher for the hysteresis control; this is due by the increase of power dissipation of the

0 1 2 3 4 5 6 7

Temps (s)

Time (s)

3.01 3.015 3.02 3.025 3.03 3.035

Time (s)

Temps (s)

Fig. 20. IGBT1 and DIODE1 power losses with hysteresis control

module components this type of control.

Puissance dissipées (W)

Power losses (W)

Fig. 21. IGBT1 junction temperature

T em

pérature de jonc

Junction temperature (°C)

 tion (°C ) A detailed dynamic model for EV was studied using two current control systems. MATLAB / Simulink were chosen from several simulation tools because of its flexibility in working with analog and digital devices, it is able to represent real-time results with the simulation time reduced. A comparative study was carried out in terms of switching frequency for power dissipated by the components of the inverter and junction temperature. The hysteresis current control has a variable switching frequency that depends on the hysteresis band, this type of control allows for fast simulations with a shorter time. The PWM current control has a fixed frequency switching and allows having junction temperatures lower than the hysteresis control.

#### **7. References**


**Part 6** 

**Robot Applications** 

