\$ \$ -. /(0 1


 

  -

 -

 

 -


-


time T\_Clk

 -

 

 


 

 

 

 


'
 

! "

2 " -\$ 
 

 "

! 
" 

 
" 

> !

" # -

% & 
 '

\$ 
 

From Control Design to FPGA Implementation 143

 () \*

one simulation interval of the plant and environment model. This requires the introduction of a *new simulation base sample time* with a resolution high enough to fit the respective controller simulation cycles. In the following, this simulation sample period will be referred to as T\_Clk. The controller design shown in Figure 12 is now simulated with T\_Clk - each unit delay only adds one T\_Clk latency. Still, there are values, that have to be updated only once every T\_Control. Since this can no longer implicitly accomplished by using unit delays based on T\_Control, it has to be expressed by unit delays enabled for one T\_Clk every T\_Control interval by an explicitly modelled trigger signal. The relation of the control interval T\_Control, the execution cycle interval T\_Clk<sup>4</sup> and the trigger signal, representing the beginning of the control cycle in the execution time domain, is displayed in Figure 13.

+ ,

 

Fig. 12. Model of a controller design with registered I/O ports.

 -



Fig. 13. Simulation rate multiplication for computation cycle oriented simulation.

 

<sup>4</sup> Simulation solver restrictions still require T\_Control being an integer multiple of the base sample

 -


 

 

> 

 

will be directly propagated through the circuit, only subject to gate delay - an effect that cannot be expressed in the Simulink model with a time resolution of T\_Control. This makes it very difficult to determine, e.g. when all current input data will be available and stable in relation to the clock edge, when to trigger the Kalman filter component on valid input data and when to submit the result values to the surrounding circuitry.

Fig. 11. Integration of the generated controller and Kalman filter components in an FPGA design.

If the logic embedding the controller component, as depicted in Figure 11, is known to provide the required clock and trigger signals, to guarantee for the stability of the input values for the computation time and to correctly time the take-over of the result values, the utilization of a combinationally integrated component should not pose a problem. Nevertheless, the correct design of the embedding logic and the behavioural validation of both, the component and the overall implementation are efforts, that cannot be undertaken on modelling level. To provide a more robust component implementation and to ease the integration of the component in a chip design a more advanced modelling approach should be considered.

## **4.2 Execution cycle oriented clocked design**

Guidelines for good design practice regarding the creation of reusable IP components have been formulated by Keating & Bricaud (2002). A key feature is the independence from the influence of embedding logic regarding the data propagation timings, especially when the timing behaviour is unknown or cannot be guaranteed at component design time.

The method to ensure the simultaneous submission of input values to the component design and the stability of those inputs for a defined time period is to provide input registers. These registers are triggered to store the input values at defined points in time, when the validity of the values can be assumed. This decouples the timing of the component from any surrounding issues, so the internal behaviour of the component can be designed based on this trigger event.

On modelling level the Simulink controller design has to be equipped with input registers, that are triggered to store the input values at the sample points defined by T\_Control. Additionally, the output has to be provided with a register to submit a stable and valid value to the embedding logic after the computation has been completed. These port registers are modelled by additional *Enabled Unit Delays*, as depicted in Figure 12.

Now the data propagation along the forward path of the data flow has been divided in multiple logic stages and is therefore delayed by multiple simulation cycles. Still, the simulation of the controller has to be accomplished within the period T\_Control, i.e. within 14 Will-be-set-by-IN-TECH

will be directly propagated through the circuit, only subject to gate delay - an effect that cannot be expressed in the Simulink model with a time resolution of T\_Control. This makes it very difficult to determine, e.g. when all current input data will be available and stable in relation to the clock edge, when to trigger the Kalman filter component on valid input data and when

> **- -**

Fig. 11. Integration of the generated controller and Kalman filter components in an FPGA

If the logic embedding the controller component, as depicted in Figure 11, is known to provide the required clock and trigger signals, to guarantee for the stability of the input values for the computation time and to correctly time the take-over of the result values, the utilization of a combinationally integrated component should not pose a problem. Nevertheless, the correct design of the embedding logic and the behavioural validation of both, the component and the overall implementation are efforts, that cannot be undertaken on modelling level. To provide a more robust component implementation and to ease the integration of the component in a

Guidelines for good design practice regarding the creation of reusable IP components have been formulated by Keating & Bricaud (2002). A key feature is the independence from the influence of embedding logic regarding the data propagation timings, especially when the

The method to ensure the simultaneous submission of input values to the component design and the stability of those inputs for a defined time period is to provide input registers. These registers are triggered to store the input values at defined points in time, when the validity of the values can be assumed. This decouples the timing of the component from any surrounding issues, so the internal behaviour of the component can be designed based on this trigger event. On modelling level the Simulink controller design has to be equipped with input registers, that are triggered to store the input values at the sample points defined by T\_Control. Additionally, the output has to be provided with a register to submit a stable and valid value to the embedding logic after the computation has been completed. These port registers are

Now the data propagation along the forward path of the data flow has been divided in multiple logic stages and is therefore delayed by multiple simulation cycles. Still, the simulation of the controller has to be accomplished within the period T\_Control, i.e. within

timing behaviour is unknown or cannot be guaranteed at component design time.

 

 to submit the result values to the surrounding circuitry.

 **! !**

chip design a more advanced modelling approach should be considered.

modelled by additional *Enabled Unit Delays*, as depicted in Figure 12.

**4.2 Execution cycle oriented clocked design**

design.

Fig. 12. Model of a controller design with registered I/O ports.

one simulation interval of the plant and environment model. This requires the introduction of a *new simulation base sample time* with a resolution high enough to fit the respective controller simulation cycles. In the following, this simulation sample period will be referred to as T\_Clk.

The controller design shown in Figure 12 is now simulated with T\_Clk - each unit delay only adds one T\_Clk latency. Still, there are values, that have to be updated only once every T\_Control. Since this can no longer implicitly accomplished by using unit delays based on T\_Control, it has to be expressed by unit delays enabled for one T\_Clk every T\_Control interval by an explicitly modelled trigger signal. The relation of the control interval T\_Control, the execution cycle interval T\_Clk<sup>4</sup> and the trigger signal, representing the beginning of the control cycle in the execution time domain, is displayed in Figure 13.

Fig. 13. Simulation rate multiplication for computation cycle oriented simulation.

<sup>4</sup> Simulation solver restrictions still require T\_Control being an integer multiple of the base sample time T\_Clk

Trajectory generator refence trajectory

**4.3 Remarks**

of this chapter.

RT RT RT RT RT RT RT Convert Convert Convert Convert Convert Convert Convert

under development with regard to the system it is to be embedded in.

as necessary" - and simulation time - "as short as possible".

**5. From simulation model to implementation model**

boolean Convert Convert Trigger Reset Integrator x-axis Enable Kalman x-axis Reset Integrator y-axis Enable Kalman y-axis Enable angle controller lx\_measured vx\_measured ly\_measured vy\_measured

RT RT

RT

output x-axis [A]

From Control Design to FPGA Implementation 145

Fig. 15. Section of the multi-rate model containing the controller design simulated with

This section introduced two different implementable designs of the example controller - first, a largely combinational design with feedback registers to be clocked with T\_Control; and second, an I/O registered design to be clocked with at least T\_Clk and triggered with period T\_Control. Both represent different styles of hardware implementation, depending on the grade of robustness, decoupling, universality or reusability, that is required for the component

The introduction of the sample time T\_Clk allows to simulate the controller execution cycles within the control system sample steps. In principle, T\_Clk does not need to be more fine-grained, than necessary to simulate the *N* designed computation steps - i.e. *T*\_*Clk* ≤ 1/*N* × *T*\_*Control*. Of course, on the chip the component will be clocked with a much higher frequency, but a true cycle-accurate simulation of a multi-MHz clocked controller design is not feasible in the early design and validation stages of a control design with a sample rate of some kHz. Scaling up the resolution of T\_Clk with continuing refinement of the execution path and logic stages is a way to find the balance between behavioural validation - "as accurate

The controller model with registered I/O resembles only the first step into the direction of *synchronous component* design. In principle it is possible to insert registers after each operation step, thus assuming the total control over the execution of the algorithm in terms of timing and determinism. The subdivision of the data path by registers shortens the combinational paths and allows to drive the hardware design with a higher clock - in turn, the overall latency will increase due to the higher cycle count. Aiming on either, the smaller overall latency of combinational designs, or the higher frequency of synchronous designs, is possible by adapting advanced hardware design principles on modelling level, but exceeds the focus

To summarize the activities discussed in this chapter, Figure 16 shows a detailed view of the *PIM-to-PSM transformation* as a part of the Simulink-based development flow from the introduction. Based on a controller partition embedded in a plant and environment model the refinement process towards a platform-specific execution model starts with the *data path design*. This iterative process substitutes the algorithmic blocks with constructs supported

NPMM

input z [A]

positions

velocities

angles

turn rates

Instrumentation

In

input x [A]

input y [A]

Sequence control control signals

T\_Control\_Trigger\_gen

higher time resolution.

This trigger signal is used to enable the following elements for one cycle of T\_Clk in order to store a value at the beginning of each T\_Control interval. First, the controller's input unit delays take over the current input values. Second, the unit delay storing the feedback output value is enabled. Third, the Discrete Integrator block has to store the last accumulated value. For this purpose, the integrator block from Figure 5 has been turned into a triggered subsystem. Until the next enable phase these values will be kept stable by the registers. Within one T\_Clk the results are propagated through the design and on a further trigger submitted to the output register. Since the result depends on the output of the Kalman filter component, the trigger for the output register is derived from the Valid signal of this component.

Fig. 14. Integration of the generated synchronized controller into an FPGA design.

Figure 14 depicts the integration of the I/O-synchronized controller into an FPGA design. The component is independent from surrounding logic, except the clock sources. Thanks to the registered I/Os and the trigger signal, the execution of the controller can be started upon the event of input data availability. This is best signalled by an upstream component like an ADC, which in turn relays its own trigger signal with period T\_Control. The execution flow is continued in a deterministic way further downstream to trigger the output, constituting a real-time behaviour as illustrated in Figure 9.

The introduction of T\_Clk as simulation base sample time allows to simulate all significant controller computation steps within the control system sample period. An decrease of the simulation sample period always increases the simulation complexity and duration, since now also the plant and environment would have to be simulated with increased rate, although they only produce a significant value change with a period of T\_Control. To minimize this effect, the construction of a *multi-rate model* is recommended, allowing to simulate the controller partition with T\_Clk while only executing the plant model once every T\_Control, as in the platform-independent model.

Figure 15 shows a section of the model surrounding the x-axis controller module, illustrating the interface adaptions necessary to accommodate a HDL compatible, cycle-oriented partition. The data lines for values coming from the plant and environment model have been supplemented with *Data type conversions* (Convert) to provide the fixed-point values and *Rate transitions* (RT) to transition from the T\_Control to the T\_Clk time domain. Furthermore, the pulse generator block producing the execution trigger pulse with a width of T\_Clk and period T\_Control is shown.

## **4.3 Remarks**

16 Will-be-set-by-IN-TECH

This trigger signal is used to enable the following elements for one cycle of T\_Clk in order to store a value at the beginning of each T\_Control interval. First, the controller's input unit delays take over the current input values. Second, the unit delay storing the feedback output value is enabled. Third, the Discrete Integrator block has to store the last accumulated value. For this purpose, the integrator block from Figure 5 has been turned into a triggered subsystem. Until the next enable phase these values will be kept stable by the registers. Within one T\_Clk the results are propagated through the design and on a further trigger submitted to the output register. Since the result depends on the output of the Kalman filter component, the trigger for the output register is derived from the Valid signal of this component.


**- -** 

Figure 14 depicts the integration of the I/O-synchronized controller into an FPGA design. The component is independent from surrounding logic, except the clock sources. Thanks to the registered I/Os and the trigger signal, the execution of the controller can be started upon the event of input data availability. This is best signalled by an upstream component like an ADC, which in turn relays its own trigger signal with period T\_Control. The execution flow is continued in a deterministic way further downstream to trigger the output, constituting a

The introduction of T\_Clk as simulation base sample time allows to simulate all significant controller computation steps within the control system sample period. An decrease of the simulation sample period always increases the simulation complexity and duration, since now also the plant and environment would have to be simulated with increased rate, although they only produce a significant value change with a period of T\_Control. To minimize this effect, the construction of a *multi-rate model* is recommended, allowing to simulate the controller partition with T\_Clk while only executing the plant model once every T\_Control, as in the

Figure 15 shows a section of the model surrounding the x-axis controller module, illustrating the interface adaptions necessary to accommodate a HDL compatible, cycle-oriented partition. The data lines for values coming from the plant and environment model have been supplemented with *Data type conversions* (Convert) to provide the fixed-point values and *Rate transitions* (RT) to transition from the T\_Control to the T\_Clk time domain. Furthermore, the pulse generator block producing the execution trigger pulse with a width of T\_Clk and

Fig. 14. Integration of the generated synchronized controller into an FPGA design.

 

 -


real-time behaviour as illustrated in Figure 9.

platform-independent model.

period T\_Control is shown.


 

This section introduced two different implementable designs of the example controller - first, a largely combinational design with feedback registers to be clocked with T\_Control; and second, an I/O registered design to be clocked with at least T\_Clk and triggered with period T\_Control. Both represent different styles of hardware implementation, depending on the grade of robustness, decoupling, universality or reusability, that is required for the component under development with regard to the system it is to be embedded in.

The introduction of the sample time T\_Clk allows to simulate the controller execution cycles within the control system sample steps. In principle, T\_Clk does not need to be more fine-grained, than necessary to simulate the *N* designed computation steps - i.e. *T*\_*Clk* ≤ 1/*N* × *T*\_*Control*. Of course, on the chip the component will be clocked with a much higher frequency, but a true cycle-accurate simulation of a multi-MHz clocked controller design is not feasible in the early design and validation stages of a control design with a sample rate of some kHz. Scaling up the resolution of T\_Clk with continuing refinement of the execution path and logic stages is a way to find the balance between behavioural validation - "as accurate as necessary" - and simulation time - "as short as possible".

The controller model with registered I/O resembles only the first step into the direction of *synchronous component* design. In principle it is possible to insert registers after each operation step, thus assuming the total control over the execution of the algorithm in terms of timing and determinism. The subdivision of the data path by registers shortens the combinational paths and allows to drive the hardware design with a higher clock - in turn, the overall latency will increase due to the higher cycle count. Aiming on either, the smaller overall latency of combinational designs, or the higher frequency of synchronous designs, is possible by adapting advanced hardware design principles on modelling level, but exceeds the focus of this chapter.

## **5. From simulation model to implementation model**

To summarize the activities discussed in this chapter, Figure 16 shows a detailed view of the *PIM-to-PSM transformation* as a part of the Simulink-based development flow from the introduction. Based on a controller partition embedded in a plant and environment model the refinement process towards a platform-specific execution model starts with the *data path design*. This iterative process substitutes the algorithmic blocks with constructs supported

Two possible designs have been introduced along with their properties concerning chip integration. Depending on the number of logic stages incorporated in the design, the platform-specific model has to be simulated with a more fine-grained simulation time resolution. The Mathworks (2011b) HDL Coder provides the possibility to create a cycle-accurate model of the HDL code derived from Simulink blocks configured to more advanced architectural implementations, which might introduce additional execution cycles not reflected in the original model. In those cases the execution control path of platform-specific model has to be further refined, regarding both the logic stages and the,

From Control Design to FPGA Implementation 147

The simulative validation of the platform-specific model including a black-box component can also be conducted using HDL co-simulation, via links to external HDL simulators, like e.g. ModelSim. This method allows to analyse the behaviour of the black-box implementation and its interaction with the surrounding controller design on a cycle-accurate level. HDL co-simulation is also supported by the Xilinx System Generator. Without the HDL co-simulation option the validation of the black-box behaviour can be conducted only after leaving the modelling level by simulating the complete HDL generated design with instantiated black-box component with an HDL simulator. This course of action is aided by the HDL Coder by giving the option to generate an HDL testbench from the recorded Simulink simulation scenario. A method for validation of the implemented chip-design within the environment of the simulated plant model is the hardware-in-the-loop (HIL) simulation, that

In conclusion it can be stated that MATLAB/Simulink supports the seamless top-down development flow for control designs with the intention of implementing the digital real-time controller on a reconfigurable logic platform. Within the extends of the Simulink block set and associated tools, and under consideration of the issues discussed in this chapter, the engineer can conduct the process of data path design by successive refinement and transformation of the controller model under constant validation against the surrounding plant model. The construction of the correct control flow of the algorithm's execution is a decisive engineering step, but it is well supported by the model-based approach. From a certain level of refinement, the usage of the HDL Coder supported block set or the the FPGA vendors' block sets provide the ability to construct and validate a hardware-oriented execution model. While the former allows a more integrated top-down modelling process, and the latter provide more platform-specific design and validation abilities, either tool facilitates in closing the gap

The work presented here is related to the research within the Collaborative Research Centre 622 "Nano-Positioning and Nano-Measuring Machines", funded by the German Research

Altera Corporation (2011). DSP Builder Handbook Volume 1 : Introduction to DSP Builder.

is supported by the vendor-specific tools DSP Builder and System Generator.

between model-based design, automatic code generation and FPGA synthesis.

simulation time resolution, to cycle-accurate simulation.

**6. Conclusion**

**7. Acknowledgements**

**8. References**

Council (DFG) under grant SFB 622 / CRC 622.

by the code generation facilities. Additionally, all signals carrying the data are converted to fixed-point representations with as minimal bit widths as possible.

Fig. 16. The platform-specific controller modelling and validation flow in detail.

The substitutions and conversions necessary to fulfill the HDL-compatibility can produce deviations from the behaviour of the specified controller design. Therefore, a *simulative validation* of each modification towards a more platform-specific representation against the initial platform-independent specification is necessary, which might motivate an iteration of the initial controller design and parametrization. Validation of the fixed-point solution is a critical factor. Both the automated conversion using the Fixed-Point tools as well as the manual conversion rely on the specified data ranges for all operations. For operators with a-priori unknown output data ranges, the tools offer the option to determine the minima and maxima from the simulation results. Obviously, the consistency of this solution is based on the completeness of the simulation scenario - not covering the real minima/maxima of certain values will result in falsely scaled fixed-point representations. As a side note, as pointed out by The Mathworks (2011b), MATLAB has a restriction on simulation data widths of 128 bit, for larger word lengths a bit-accurate simulation is not possible, and the simulation results will deviate from the actual execution results.

The data path design approach also applies to the utilization of the vendor specific tools Xilinx, Inc. (2011) System Generator or Altera Corporation (2011) DSP Builder. But, when using these tools, from some point onwards during refinement process the model will have to be completely (re-)designed using the vendor-specifc block sets to create a platform-specific model, that is compatible to the respective code generation tool chains.

The activity of *behavioural design* applies for all mentioned tools in a similar fashion and concerns the design of the execution behaviour of the algorithm on the hardware. 18 Will-be-set-by-IN-TECH

by the code generation facilities. Additionally, all signals carrying the data are converted to

**--

 **-**

> - -

-

 **-**

**-**

**-**

**--

**-** 

**-** 

**-**

**-**


 **-**

Fig. 16. The platform-specific controller modelling and validation flow in detail.

 **-** The substitutions and conversions necessary to fulfill the HDL-compatibility can produce deviations from the behaviour of the specified controller design. Therefore, a *simulative validation* of each modification towards a more platform-specific representation against the initial platform-independent specification is necessary, which might motivate an iteration of the initial controller design and parametrization. Validation of the fixed-point solution is a critical factor. Both the automated conversion using the Fixed-Point tools as well as the manual conversion rely on the specified data ranges for all operations. For operators with a-priori unknown output data ranges, the tools offer the option to determine the minima and maxima from the simulation results. Obviously, the consistency of this solution is based on the completeness of the simulation scenario - not covering the real minima/maxima of certain values will result in falsely scaled fixed-point representations. As a side note, as pointed out by The Mathworks (2011b), MATLAB has a restriction on simulation data widths of 128 bit, for larger word lengths a bit-accurate simulation is not possible, and the simulation results

The data path design approach also applies to the utilization of the vendor specific tools Xilinx, Inc. (2011) System Generator or Altera Corporation (2011) DSP Builder. But, when using these tools, from some point onwards during refinement process the model will have to be completely (re-)designed using the vendor-specifc block sets to create a platform-specific

The activity of *behavioural design* applies for all mentioned tools in a similar fashion and concerns the design of the execution behaviour of the algorithm on the hardware.

model, that is compatible to the respective code generation tool chains.

**-**

**-**

**-** 

**-** 

**-** 

**-**  

 

 !-

 --

 -

> -

 -

> -

**-**

**-** 

**-**

**-**

fixed-point representations with as minimal bit widths as possible.

#\$%


 " - 

' - --

will deviate from the actual execution results.





> -- -

Two possible designs have been introduced along with their properties concerning chip integration. Depending on the number of logic stages incorporated in the design, the platform-specific model has to be simulated with a more fine-grained simulation time resolution. The Mathworks (2011b) HDL Coder provides the possibility to create a cycle-accurate model of the HDL code derived from Simulink blocks configured to more advanced architectural implementations, which might introduce additional execution cycles not reflected in the original model. In those cases the execution control path of platform-specific model has to be further refined, regarding both the logic stages and the, simulation time resolution, to cycle-accurate simulation.

The simulative validation of the platform-specific model including a black-box component can also be conducted using HDL co-simulation, via links to external HDL simulators, like e.g. ModelSim. This method allows to analyse the behaviour of the black-box implementation and its interaction with the surrounding controller design on a cycle-accurate level. HDL co-simulation is also supported by the Xilinx System Generator. Without the HDL co-simulation option the validation of the black-box behaviour can be conducted only after leaving the modelling level by simulating the complete HDL generated design with instantiated black-box component with an HDL simulator. This course of action is aided by the HDL Coder by giving the option to generate an HDL testbench from the recorded Simulink simulation scenario. A method for validation of the implemented chip-design within the environment of the simulated plant model is the hardware-in-the-loop (HIL) simulation, that is supported by the vendor-specific tools DSP Builder and System Generator.

## **6. Conclusion**

In conclusion it can be stated that MATLAB/Simulink supports the seamless top-down development flow for control designs with the intention of implementing the digital real-time controller on a reconfigurable logic platform. Within the extends of the Simulink block set and associated tools, and under consideration of the issues discussed in this chapter, the engineer can conduct the process of data path design by successive refinement and transformation of the controller model under constant validation against the surrounding plant model. The construction of the correct control flow of the algorithm's execution is a decisive engineering step, but it is well supported by the model-based approach. From a certain level of refinement, the usage of the HDL Coder supported block set or the the FPGA vendors' block sets provide the ability to construct and validate a hardware-oriented execution model. While the former allows a more integrated top-down modelling process, and the latter provide more platform-specific design and validation abilities, either tool facilitates in closing the gap between model-based design, automatic code generation and FPGA synthesis.

## **7. Acknowledgements**

The work presented here is related to the research within the Collaborative Research Centre 622 "Nano-Positioning and Nano-Measuring Machines", funded by the German Research Council (DFG) under grant SFB 622 / CRC 622.

## **8. References**

Altera Corporation (2011). DSP Builder Handbook Volume 1 : Introduction to DSP Builder.

**7** 

**Describing Function Recording with** 

Krunoslav Horvat, Ognjen Kuljaca and Tomislav Sijak

Describing function is an equivalent gain of nonlinear element, defined by the harmonic linearization method of nonlinear static characteristic (Novogranov, 1986, Slotine and Li, 1991, Schwarz and Gran, 2001, Vukic et al. 2003 and many others). It is a known method of analysis and synthesis when nonlinear system can be decoupled into linear and nonlinear parts (Fig. 1). If the linear part of the system has the characteristics of low-pass filter and if we apply periodical signal to the system, output signal will have the same base frequency as

*<sup>F</sup>x <sup>y</sup> t <sup>n</sup> <sup>G</sup> p <sup>L</sup>*

Fig. 1. Nonlinear system represented with decoupled nonlinear F(x, px) and linear parts

, can be expressed by the Fourier expressions:

2

0

2

0

1

1

where *YP*1 and *YQ*1 are first Fourier coefficients.

Im *NP Q*

*N PQ*

*y t Y jY e*

<sup>1</sup> *Y F X t td t P m* sin sin

<sup>1</sup> *Y F X t td t Q m* sin cos

 

If the amplitudes of higher harmonics are relatively small when compared to the amplitude of the first harmonic, i.e. if filter hypothesis holds, output signal can be approximated by it's first harmonic. Mathematically, first harmonic of the output signal, for the sinusoidal input

> 1 1

*y t Y tY t*

1 1

 

 

 *j t*

(2)

(3)

sin cos

**1. Introduction** 

GL(p), p=d/dt.

signal *X t <sup>m</sup>* sin

input signal with damped higher frequencies.

*X t <sup>m</sup>* sin

**Simulink and MATLAB** 

*yt*

(1)

*Brodarski Institute* 

*Croatia* 


## **Describing Function Recording with Simulink and MATLAB**

Krunoslav Horvat, Ognjen Kuljaca and Tomislav Sijak *Brodarski Institute Croatia* 

## **1. Introduction**

20 Will-be-set-by-IN-TECH

148 Technology and Engineering Applications of Simulink

Amthor, A., Zschack, S. & Ament, C. (2008). Position control on nanometer scale based on an

Auger, D. (2008). Programmable hardware systems using model-based design, *2008 IET and*

Carletta, J. & Veillette, R. (2005). A methodology for FPGA-based control implementation,

Caspi, P. & Maler, O. (2005). From control loops to real-time programs, *Handbook of networked*

Krasner, J. (2004). Model-based design and beyond: Solutions for today's embedded systems requirements, *Analyst report, American Technology International* (January): 1–12. Monmasson, E. & Cirstea, M. (2007). FPGA Design Methodology for Industrial Control Systems - A Review, *IEEE Transactions on Industrial Electronics* 54(4): 1824–1842. Müller, M., Fengler, W., Amthor, A. & Ament, C. (2009). Model-driven development and

Pacholik, A., Klöckner, J., Müller, M., Gushchina, I. & Fengler, W. (2011). LiSARD: LabVIEW

Zschäck, S., Amthor, A., Müller, M., Klöckner, J., Ament, C. & Fengler, W. (2010). Integrated

system development process for high-precision motion control systems, *2010 IEEE*

*2011, Cancun (Mexico)*, IEEE Computer Society CPS, pp. 442–447.

*International Conference on Control Applications*, IEEE, pp. 344–350.

Rushton, A. (2011). *VHDL for Logic Synthesis*, 3rd edn, John Wiley & Sons.

The Mathworks (2011a). MATLAB/Simulink Online Documentation. URL: *http://www.mathworks.com/help/toolbox/simulink/* The Mathworks (2011b). Simulink HDL Coder User's Guide R2011b. Xilinx, Inc. (2011). Xilinx System Generator for DSP User Guide.

The Mathworks (2010). Simulink Fixed Point 6 User Guide.

multiprocessor implementation of a dynamic control algorithm for nanopositioning and nanomeasuring machines, *Journal of Systems and Control Engineering* 223: 417–429.

Integrated Softcore Architecture for Reconfigurable Devices, *2011 International Conference on Reconfigurable Computing and FPGAs (ReConFig '11), 30. Nov.- 02. Dec.*

*Electronics Weekly Conference on Programmable Hardware Systems* .

*IEEE Transactions on Control Systems Technology* 13(6): 977–987.

*and embedded control systems*, Birkhäuser, pp. 395–418. Keating, M. & Bricaud, P. (2002). *Reuse Methodology Manual*, 3rd edn, Springer.

*Electronics*, IEEE, Orlando (USA), pp. 2568–2573.

adaptive friction compensation scheme, *2008 34th Annual Conference of IEEE Industrial*

Describing function is an equivalent gain of nonlinear element, defined by the harmonic linearization method of nonlinear static characteristic (Novogranov, 1986, Slotine and Li, 1991, Schwarz and Gran, 2001, Vukic et al. 2003 and many others). It is a known method of analysis and synthesis when nonlinear system can be decoupled into linear and nonlinear parts (Fig. 1). If the linear part of the system has the characteristics of low-pass filter and if we apply periodical signal to the system, output signal will have the same base frequency as input signal with damped higher frequencies.

Fig. 1. Nonlinear system represented with decoupled nonlinear F(x, px) and linear parts GL(p), p=d/dt.

If the amplitudes of higher harmonics are relatively small when compared to the amplitude of the first harmonic, i.e. if filter hypothesis holds, output signal can be approximated by it's first harmonic. Mathematically, first harmonic of the output signal, for the sinusoidal input signal *X t <sup>m</sup>* sin, can be expressed by the Fourier expressions:

$$\begin{cases} y\_N\left(t\right) \approx Y\_{p1} \sin \alpha t + Y\_{Q1} \cos \alpha t\\ y\_N\left(t\right) \approx \text{Im}\left\{ \left(Y\_{P1} + jY\_{Q1}\right)e^{j\alpha t} \right\} \end{cases} \tag{1}$$

$$Y\_{P1} = \frac{1}{\pi} \int\_0^{2\pi} F(X\_m \sin \alpha t) \sin \alpha t \, d\left(\alpha t\right) \tag{2}$$

$$Y\_{Q1} = \frac{1}{\pi} \int\_0^{2\pi} F(X\_m \sin \alpha t) \cos \alpha t \, d\left(\alpha t\right) \tag{3}$$

where *YP*1 and *YQ*1 are first Fourier coefficients.

Describing Function Recording with Simulink and MATLAB 151

An interesting approach to obtaining describing function using Simulink is given by Schwartz and Gran. The approach is very simple, yet effective. In essence, equations (2) and (3) are directly applied in Simulink on nonlinear element outputs. A simulation scheme adapted from Schwartz is given in Fig. 3. In this case dead zone nonlinearity is analyzed.

Triggered

In2

Scope 4

In1

Sine Wave 1

Constant

E

Scope 1

Constant 2

pi ./omega 1

Product 1

> Gain 1

1

Integrator 1

Divide 1

1

s

Pulse

Generator

Subsystem

**3. Describing function by numerical integration** 

To Workspace

Clock

Scope 3

Constant 1

pi ./omega 1

> Scope 2

time

Instead of dead zone, any other nonlinearity can be included in model.

Fig. 3. Simulink model 'csm' for generating describing function (Schwartz)

Sine Wave

Gain

Dead Zone

E'\* u

Product

Integrator

Divide

Compare

>= 0

To Zero

1

Constant 3

2\*pi /omega 1

Add

s

in Fig. 5.

Model is called from Matlab script (adapted form Schwartz and Gran as well) as shown in Fig. 4. simulation parameters are passed to Simulink schematic from Matlab script, as given

Dynamic

Describing function is the ratio between first harmonic of the output signal and input signal in complex form:

$$f\_N\left(X\_m\right) = P\left(X\_m\right) + jQ\left(X\_m\right) = \frac{Y\_{p1}}{X\_m} + j\frac{Y\_{Q1}}{X\_m} \tag{4}$$

where *P X <sup>m</sup>* and *Q X <sup>m</sup>* are coefficients of the harmonic linearization (Novogranov, 1986, Slotine and Li, 1991, Schwarz and Gran, 2001, Vukic et al., 2003 and many others).

Determination of describing function boils down to the determination of integral expressions for the known static characteristic of the nonlinear part of the system.

There are many conventional nonlinearities for which static characteristics and describing functions are theoretically derived and given in analytical form.

The problem arises when the static characteristic of nonlinear system cannot be analytically expressed or integral expressions can not be solved. In that case describing function can be determined by experiment or simulation (Kuljaca et al., Sijak et al.k 2004, 2005, 2007a, 2007b) or some method of numerical integration (Schwarz and Gran, 2001).

## **2. Nonlinear elements in Simulink**

Nonlinear elements are given in Simulink library Discontinuities. Twelve discontinuities given there are shown in Fig. 2.

Fig. 2. Discontinuities (nonlinearities) given in Simulink.

We can see that there are three elements called "Dynamic" (Dead Zone Dynamic, Saturation Dynamic and Rate Limiter Dynamic). The term dynamic is considered at these elements with respect to change of nonlinearity limits (output values limits or rate change limits), not as dynamics in control systems sense (i.e. behavior in time domain). Never the less, we will not deal with such elements here. Describing function method analysis or synthesis is not suitable for systems with varying parameters.

One can see that Simulink gives only basic non-dynamic nonlinearities. More nonlinearities and their describing function can be found in work by Vukic et al., (2003). It is clear that one needs to build them out of basic Simulink models or write them as m or s functions. Since given functions don't have dynamics it is recommended to build them out of given Simulink blocks, or if that is not possible, then using type 2 s-functions. Even in that case some real time extension software for Simulink and Matlab might not be able to handle user designed s-functions.

Describing function is the ratio between first harmonic of the output signal and input signal

*<sup>P</sup>*<sup>1</sup> *<sup>Q</sup>*<sup>1</sup>

*G X P X jQ X j X X*

where *P X <sup>m</sup>* and *Q X <sup>m</sup>* are coefficients of the harmonic linearization (Novogranov, 1986,

Determination of describing function boils down to the determination of integral

There are many conventional nonlinearities for which static characteristics and describing

The problem arises when the static characteristic of nonlinear system cannot be analytically expressed or integral expressions can not be solved. In that case describing function can be determined by experiment or simulation (Kuljaca et al., Sijak et al.k 2004, 2005, 2007a, 2007b)

Nonlinear elements are given in Simulink library Discontinuities. Twelve discontinuities

We can see that there are three elements called "Dynamic" (Dead Zone Dynamic, Saturation Dynamic and Rate Limiter Dynamic). The term dynamic is considered at these elements with respect to change of nonlinearity limits (output values limits or rate change limits), not as dynamics in control systems sense (i.e. behavior in time domain). Never the less, we will not deal with such elements here. Describing function method analysis or synthesis is not

One can see that Simulink gives only basic non-dynamic nonlinearities. More nonlinearities and their describing function can be found in work by Vukic et al., (2003). It is clear that one needs to build them out of basic Simulink models or write them as m or s functions. Since given functions don't have dynamics it is recommended to build them out of given Simulink blocks, or if that is not possible, then using type 2 s-functions. Even in that case some real time extension software for Simulink and Matlab might not be able to handle user designed

Quantizer

Hit Crossing

Wrap To Zero

Saturation Dynamic

y 1

Rate Limiter Dynamic

up u lo

Rate Limiter

up u lo

Relay Saturation

*m m Y Y*

(4)

*Nm m m*

Slotine and Li, 1991, Schwarz and Gran, 2001, Vukic et al., 2003 and many others).

expressions for the known static characteristic of the nonlinear part of the system.

functions are theoretically derived and given in analytical form.

or some method of numerical integration (Schwarz and Gran, 2001).

**2. Nonlinear elements in Simulink** 

Dead Zone Dynamic

suitable for systems with varying parameters.

y

Fig. 2. Discontinuities (nonlinearities) given in Simulink.

up u lo

Coulomb & Viscous Friction

s-functions.

Backlash

Dead Zone

given there are shown in Fig. 2.

in complex form:

## **3. Describing function by numerical integration**

An interesting approach to obtaining describing function using Simulink is given by Schwartz and Gran. The approach is very simple, yet effective. In essence, equations (2) and (3) are directly applied in Simulink on nonlinear element outputs. A simulation scheme adapted from Schwartz is given in Fig. 3. In this case dead zone nonlinearity is analyzed. Instead of dead zone, any other nonlinearity can be included in model.

Fig. 3. Simulink model 'csm' for generating describing function (Schwartz)

Model is called from Matlab script (adapted form Schwartz and Gran as well) as shown in Fig. 4. simulation parameters are passed to Simulink schematic from Matlab script, as given in Fig. 5.

Describing Function Recording with Simulink and MATLAB 153

Triggered subsystem is shown in Fig. 6. The whole file is built to stop simulation after one

Let us now see how Simulink file with Matlab function will perform for a common nonlinearity, saturation with saturation limit 0.5. Saturation type nonlinearity is shown in

F(z)

Saturation has known describing function given by (5). It can be seen that imaginary part of

<sup>2</sup> ( ) arcsin *m m* 1 ;

*tg cc c G GZ A c*

Let us now see the results of running our simulation script in order to obtain describing

the gain is zero. That is a case with all symmetrical nonlinearities.

where Am is amplitude of signal entering saturation.

function for saturation. Simulation results are shown in Fig. 8 and Fig. 9.

c=

Trigger

To Workspace1 DFimag

To Workspace DFreal

z


(5)

*m m m*

*A A A*

full period for a given frequency.

In 2 2

Fig. 6. Triggered subsystem

Fig. 7.

Fig. 7. Saturation

In 1 1

Fig. 4. Matlab script for describing function generation


Fig. 5. Simulation configuration parameters – parameter omega1 is taken from Matlab script

Triggered subsystem is shown in Fig. 6. The whole file is built to stop simulation after one full period for a given frequency.

Fig. 6. Triggered subsystem

152 Technology and Engineering Applications of Simulink

Fig. 5. Simulation configuration parameters – parameter omega1 is taken from Matlab script

Fig. 4. Matlab script for describing function generation

Let us now see how Simulink file with Matlab function will perform for a common nonlinearity, saturation with saturation limit 0.5. Saturation type nonlinearity is shown in Fig. 7.

Fig. 7. Saturation

Saturation has known describing function given by (5). It can be seen that imaginary part of the gain is zero. That is a case with all symmetrical nonlinearities.

$$\mathbf{G} = \mathbf{G}(Z\_m) = \frac{2 \cdot \text{tg}\alpha}{\Pi} \left( \arcsin \frac{c}{A\_m} + \frac{c}{A\_m} \cdot \sqrt{1 \cdot \frac{c^2}{A\_m}} \right); A\_m > c \tag{5}$$

where Am is amplitude of signal entering saturation.

Let us now see the results of running our simulation script in order to obtain describing function for saturation. Simulation results are shown in Fig. 8 and Fig. 9.

Describing Function Recording with Simulink and MATLAB 155

zero. Obtained result from simulation is given complex matrix s. In case like this, where there is no frequency dependency and phase is zero, there is no need for plotting describing function in three dimensions. To make function plot in two dimensions we run the

0 0.5 1 1.5 2 2.5 3

Amplitude

We will see in subsequent sections that this representation of symmetrical nonlinearities as gain dependent on amplitude of nonlinearity input signal is important for stability

This technique can be used to obtain describing functions of real nonlinearities by using Simulink model from Fig. 3 with appropriate real time additions with Real Time Window Target, xPc Target or other real time systems for use with Simulink. Of course, in that case

Fig. 10. Describing function gain – two dimensional representation

amplitudes will have to be given one by one, not as here as one vector.

plot(E,real(gdf)),grid, ylabel('Describing function gain');

following short script:

xlabel('Amplitude');

Two dimensional describing function gain is shown in

gdf=s(1,:);

0.2

analysis.

0.3

0.4

0.5

0.6

Describing function gain

0.7

0.8

0.9

1

Describing Function Amplitude vs. Frequency and Amplitude of Input

Fig. 8. Describing function gain as function of amplitude and frequency

Describing Function Phase vs. Frequency and Amplitude of Input

Fig. 9. Describing function phase as function of amplitude and frequency

It can be seen from Fig. 8 that describing function in this case is not dependent on frequency. Also, phase is zero since saturation is symmetrical nonlinearity. That can be seen from very small values of describing function phase as given in Fig. 9. These values are practically

Describing Function Amplitude vs. Frequency and Amplitude of Input

0

Describing Function Phase vs. Frequency and Amplitude of Input

0

It can be seen from Fig. 8 that describing function in this case is not dependent on frequency. Also, phase is zero since saturation is symmetrical nonlinearity. That can be seen from very small values of describing function phase as given in Fig. 9. These values are practically

0

Fig. 9. Describing function phase as function of amplitude and frequency

1

2

Amplitude

3

0

Fig. 8. Describing function gain as function of amplitude and frequency

5

5

Frequency (rad/sec)

Frequency (rad/sec)

10 0.2

10 0

1

2

Describing Function Phase

3

4

x 10-8

0.4

0.6

Describing Function Gain

0.8

1

1

2

Amplitude

3

zero. Obtained result from simulation is given complex matrix s. In case like this, where there is no frequency dependency and phase is zero, there is no need for plotting describing function in three dimensions. To make function plot in two dimensions we run the following short script:

gdf=s(1,:); plot(E,real(gdf)),grid, ylabel('Describing function gain'); xlabel('Amplitude');

Two dimensional describing function gain is shown in

Fig. 10. Describing function gain – two dimensional representation

We will see in subsequent sections that this representation of symmetrical nonlinearities as gain dependent on amplitude of nonlinearity input signal is important for stability analysis.

This technique can be used to obtain describing functions of real nonlinearities by using Simulink model from Fig. 3 with appropriate real time additions with Real Time Window Target, xPc Target or other real time systems for use with Simulink. Of course, in that case amplitudes will have to be given one by one, not as here as one vector.

Describing Function Recording with Simulink and MATLAB 157

Product

Constant A

Constant 1 pi ./omega 1

> Integrator 1 1 s

Constant 2 pi ./omega 1 Divide 1

Divide

Constant 3 2\*pi /omega 1

Clock

Integrator 1 s

Sine Wave 1

Product 1

Fig. 12. Numerical integration model 'csm\_fuzzy'– describing function recording

Fuzzy Logic Controller

Scope 3

Scope 1

2

Fig. 13. Matlab script for running 'csm\_fuzzy' model

Gain 3 0.001 Gain 2 1.2

Derivative du/dt

Gain 1 1

Sine Wave

Scope 2

Gain A\* u

> Triggered Subsystem

To Workspace time

Compare To Zero >= 0

Scope 4

In1 In2 Add

## **4. Fuzzy control systems and their describing functions in Simulink**

Use of fuzzy control systems is today well known in control systems design. Unfortunately, while there are well developed methods for design or training of fuzzy systems, stability analysis remains problematic. If the fuzzy system is of Mamdani type, than its mathematical description can be very complicated and stability analysis is complicated as well. However, if we can obtain describing function of fuzzy system, then we can use known stability analysis methods developed for linear systems. One type of experimental method of obtaining describing function for fuzzy elements is given in (Kuljaca et al.) and (Sijak et al., 2005, 2007b). Here, we will use numerical integration in Simulink as given by Schwartz. Simulation file will have to be reworked since we cannot use vectors in this simulation due to fuzzy element.

Of course, there are different fuzzy controllers and elements, but describing function for all of them can be achieved using this method. One should note that here we are not talking about adaptive fuzzy controllers, but about fuzzy controllers with fixed membership functions and rule base.

Let us first look on one example of fuzzy controller, as shown in Fig. 11 (Kuljaca et al., Sijak et al. 2007a), where kp and kd are gains.

Fig. 11. The block diagram of fuzzy controller

First of all, it can be seen that fuzzy element here has two inputs. Fuzzy elements in general can have more inputs and outputs than here, but, this structure is most often in control systems (either as given here with derivative of the error or with integral of the error). When more inputs are used it becomes extremely complicated to tune the controller.

Numerical integration Simulink model (Schwartz) with fuzzy controller will look like in Fig. 12, without vector inputs for amplitudes since fuzzy block in Simulink cannot handle such input type. This is also much closer to real measurement, since in real measurement we would not able to use vector inputs. Meaning of this is that we need to use adjusted script to run the Simulink model. Adjusted script is given in Fig. 13.

Fuzzy Logic controller Simulink block is regular block from Simulink Fuzzy Logic toolbox library.

Fuzzy element itself is developed using FIS editor from Matlab Fuzzy Logic Toolbox.

This is an illustrative example; however, it is quite useful in giving an insight in use of describing function for fuzzy controllers. The given method can give a graphic representation of any fuzzy controller based on error and derivatives of error signal. In most cases fuzzy controllers are also symmetrical, thus subsequent stability analysis will be simpler. Amplitudes and frequencies are to be chosen to satisfy expected operational environment of controller.

Gain 3

156 Technology and Engineering Applications of Simulink

Use of fuzzy control systems is today well known in control systems design. Unfortunately, while there are well developed methods for design or training of fuzzy systems, stability analysis remains problematic. If the fuzzy system is of Mamdani type, than its mathematical description can be very complicated and stability analysis is complicated as well. However, if we can obtain describing function of fuzzy system, then we can use known stability analysis methods developed for linear systems. One type of experimental method of obtaining describing function for fuzzy elements is given in (Kuljaca et al.) and (Sijak et al., 2005, 2007b). Here, we will use numerical integration in Simulink as given by Schwartz. Simulation file will have to be reworked since we cannot use vectors in this simulation due to fuzzy element.

Of course, there are different fuzzy controllers and elements, but describing function for all of them can be achieved using this method. One should note that here we are not talking about adaptive fuzzy controllers, but about fuzzy controllers with fixed membership

Let us first look on one example of fuzzy controller, as shown in Fig. 11 (Kuljaca et al., Sijak

error

errordot

First of all, it can be seen that fuzzy element here has two inputs. Fuzzy elements in general can have more inputs and outputs than here, but, this structure is most often in control systems (either as given here with derivative of the error or with integral of the error). When

Numerical integration Simulink model (Schwartz) with fuzzy controller will look like in Fig. 12, without vector inputs for amplitudes since fuzzy block in Simulink cannot handle such input type. This is also much closer to real measurement, since in real measurement we would not able to use vector inputs. Meaning of this is that we need to use adjusted script to

Fuzzy Logic controller Simulink block is regular block from Simulink Fuzzy Logic toolbox

This is an illustrative example; however, it is quite useful in giving an insight in use of describing function for fuzzy controllers. The given method can give a graphic representation of any fuzzy controller based on error and derivatives of error signal. In most cases fuzzy controllers are also symmetrical, thus subsequent stability analysis will be simpler. Amplitudes and frequencies are to be chosen to satisfy expected operational

Fuzzy element itself is developed using FIS editor from Matlab Fuzzy Logic Toolbox.

Fuzzy element y

kp

de dt

more inputs are used it becomes extremely complicated to tune the controller.

run the Simulink model. Adjusted script is given in Fig. 13.

k d

functions and rule base.

library.

environment of controller.

et al. 2007a), where kp and kd are gains.

e

Fig. 11. The block diagram of fuzzy controller

**4. Fuzzy control systems and their describing functions in Simulink** 

Fig. 12. Numerical integration model 'csm\_fuzzy'– describing function recording

Fig. 13. Matlab script for running 'csm\_fuzzy' model

Describing Function Recording with Simulink and MATLAB 159

Describing function phase as function of amplitude and frequency shown in Fig. 17 shows very small values and for all practical purposes can be considered to be zero. Since our fuzzy controller design is symmetrical one that is in compliance with theoretical expectations.

Describing Function Amplitude vs. Frequency and Amplitude of Input

0

Describing Function Phase vs. Frequency and Amplitude of Input

0

Once describing function is obtained, it could be used for stability analysis. We will use a simplified secondary frequency control model of the isolated thermo power system with

0

**5. Stability analysis with describing functions in Simulink – Case study** 

generation rate constraint (Sijak et al.,2001). Control model is shown in Fig. 18.

1

Fig. 17. Describing function phase as function of amplitude and frequency

2

Amplitude

3 -0.005 0 0.005 0.01 0.015 0.02 0.025

Describing Function Phase

2

4

6

Frequency (rad/sec)

8

10

0

1

Fig. 16. Describing function gain as function of amplitude and frequency

2

Amplitude

3 1 1.2 1.4 1.6 1.8 2

Describing Function Gain

2

4

6

Frequency (rad/sec)

8

10

In example given here gains kp and kd are set to 1.2 and 0.001 respectively, min-max and centroid defuzzification Mamdani type fuzzy controller membership functions are given in Fig. 14 and rulebase is given in Fig. 15.

Fig. 14. Membership functions of the fuzzy element: a) proportional input, b) derivative input and c) output


Fig. 15. Rulebase of the fuzzy element

Describing function for given fuzzy controller is given in Fig. 16 and Fig. 17. It can be seen from Fig. 16 that describing function gain is function of amplitude and it really does not depend on frequency in this case.

In example given here gains kp and kd are set to 1.2 and 0.001 respectively, min-max and centroid defuzzification Mamdani type fuzzy controller membership functions are given in

( ) *<sup>p</sup> μ x*

 *<sup>d</sup> x*

1

1

21 0 1

6 5 4 3 2 3 4 5 6

*y*

1

2 1 0 1

*p x d x*

Fig. 15. Rulebase of the fuzzy element

depend on frequency in this case.

2 1 0 1

4 3 2 3 4

3 2 3

Fig. 14. Membership functions of the fuzzy element: a) proportional input, b) derivative

Describing function for given fuzzy controller is given in Fig. 16 and Fig. 17. It can be seen from Fig. 16 that describing function gain is function of amplitude and it really does not

*p x*

*y*

*d x*

Fig. 14 and rulebase is given in Fig. 15.

input and c) output

Describing function phase as function of amplitude and frequency shown in Fig. 17 shows very small values and for all practical purposes can be considered to be zero. Since our fuzzy controller design is symmetrical one that is in compliance with theoretical expectations.

Fig. 16. Describing function gain as function of amplitude and frequency

Fig. 17. Describing function phase as function of amplitude and frequency

## **5. Stability analysis with describing functions in Simulink – Case study**

Once describing function is obtained, it could be used for stability analysis. We will use a simplified secondary frequency control model of the isolated thermo power system with generation rate constraint (Sijak et al.,2001). Control model is shown in Fig. 18.

Describing Function Recording with Simulink and MATLAB 161

*F x i i i i im i im* , *px FX x X t* , , sin 

Fig. 20. The structure of the harmonically linearised nonlinear system and fuzzy controller

1, , 0 *FX G j FX G j* 11 1 22 2 *mL m L*

and nonlinear part of the plant respectively, then the stability analysis of the system can be

Let us now deal with nonlinearities. Generation rate constraint is saturation type nonlinearity shown in Fig. 21. Describing function for such nonlinearities is known and it is

 

output

The self-oscillations (limit cycles) of the system in Fig. 20 are described by solution of:

4 sin 2 cos ( ) , arcsin 2 4 *N m <sup>m</sup> <sup>m</sup> G Z <sup>Z</sup> <sup>Z</sup>*

and *F X* 2 2 *<sup>m</sup>* ,

 ,

<sup>1</sup>*x* <sup>2</sup>*x y j*

Harmonically linearized system is shown in Fig. 20.

 ,*F*<sup>1</sup> *X*1*<sup>m</sup>* ( ) *GL*<sup>1</sup> *j*

conducted using the solution of complex equation (7).

*f j*

Assuming that *F X* 1 1 *<sup>m</sup>* ,

Fig. 21. Generation rate constraint

given in (8).


 (7)

input

 

(8)

are describing functions of the fuzzy controller

*F*<sup>2</sup> *X*2*<sup>m</sup>* ( ) *GL*<sup>2</sup> *j*

(6)

Fig. 18. Power system secondary load-frequency control model

where:

1 1 *<sup>G</sup> G G sT* - the transfer function of the turbine governor

*TCH* - the steam turbine time constant

Pm – the change of the mechanical power of the turbine

PL – the change of the power system load

Pr – the power system active power reference change

f – the power system frequency change

1 *<sup>s</sup> <sup>s</sup> s <sup>K</sup> <sup>G</sup> sT* - the power system transfer function

*R* - the static speed droop of the uncontrolled system

*F z*( ) - the power system generation rate constraint, static characteristic of saturation nonlinearity.

The system given in Fig. 18 consists of nonlinear parts divided by linear parts for which the filter hypothesis is satisfied. Such system can be represented as in Fig. 19.

Fig. 19. The structure of the nonlinear system and fuzzy controller

where:

 , , *ii i <sup>d</sup> F x px p dt* , - nonlinear parts of the system G s - transfer functions of the linear parts of the system Li

Assuming that the filter hypothesis is valid for *GLip* , the nonlinearities can be harmonically linearized and their describing functions used instead (Vukic et al., Netushil et al.).


controller bl

*F z*( ) - the power system generation rate constraint, static characteristic of saturation

The system given in Fig. 18 consists of nonlinear parts divided by linear parts for which the

<sup>1</sup> <sup>1</sup> <sup>1</sup> *F x* , *px* ( ) <sup>1</sup> *G p <sup>L</sup>* <sup>2</sup> <sup>2</sup> <sup>2</sup> *F x* , *px* ( ) <sup>2</sup> *G p <sup>L</sup>*

Assuming that the filter hypothesis is valid for *GLip* , the nonlinearities can be harmonically

linearized and their describing functions used instead (Vukic et al., Netushil et al.).

<sup>1</sup>*x* <sup>2</sup>*x yt*

Fuzzy controller

 

1 s - Pm

PL

1/R

1 TCH


filter hypothesis is satisfied. Such system can be represented as in Fig. 19.

Fig. 19. The structure of the nonlinear system and fuzzy controller

*<sup>d</sup> F x px p dt* , - nonlinear parts of the system

G s - transfer functions of the linear parts of the system Li

z

+ -

Fig. 18. Power system secondary load-frequency control model

Pm – the change of the mechanical power of the turbine

Pr – the power system active power reference change

*sT* - the power system transfer function *R* - the static speed droop of the uncontrolled system

<sup>G</sup> <sup>+</sup> <sup>g</sup>

*TCH* - the steam turbine time constant

PL – the change of the power system load

f – the power system frequency change

where:

1 *<sup>G</sup>*

1 *<sup>s</sup> <sup>s</sup>*

nonlinearity.

where:

, , *ii i*

*<sup>K</sup> <sup>G</sup>*

*s*

*f t*


*G*

1

*sT*

*G*

$$F\_i(\mathbf{x}\_i, \mathbf{p}\mathbf{x}\_i) \approx F\_i(X\_{im}, o\mathbf{o})\_\prime \mathbf{x}\_i = X\_{im} \sin(o\,t) \tag{6}$$

Harmonically linearized system is shown in Fig. 20.

Fig. 20. The structure of the harmonically linearised nonlinear system and fuzzy controller

The self-oscillations (limit cycles) of the system in Fig. 20 are described by solution of:

$$\mathbf{1} + F\_1(X\_{1m'}o\rho)\mathbf{G}\_{11}(j\rho)F\_2(X\_{2m'}o\rho)\mathbf{G}\_{12}(j\rho) = \mathbf{0} \tag{7}$$

Assuming that *F X* 1 1 *<sup>m</sup>* , and *F X* 2 2 *<sup>m</sup>* , are describing functions of the fuzzy controller and nonlinear part of the plant respectively, then the stability analysis of the system can be conducted using the solution of complex equation (7).

Let us now deal with nonlinearities. Generation rate constraint is saturation type nonlinearity shown in Fig. 21. Describing function for such nonlinearities is known and it is given in (8).

Fig. 21. Generation rate constraint

$$\mathcal{G}\_{\mathcal{N}}(Z\_m) = \frac{4}{\pi} \left( \frac{\alpha}{2} - \frac{\sin 2\alpha}{4} + \frac{\cos \alpha}{Z\_m \sqrt{\delta}} \right) , \alpha = \arcsin \frac{\delta}{Z\_m} \tag{8}$$

Describing Function Recording with Simulink and MATLAB 163

0 0.5 1 1.5 2 2.5 3

Amplitude

The static characteristics of the fuzzy regulator, i.e. its describing functions, shows that harmonically linearized fuzzy regulator can be regarded as input dependent variable gain element. Consequently, the stability analysis can be conducted in the plane GNF = f(GNZ),

GNZ = GN(Zm) - describing function of the power system generation rate constraint.

1 1

From (9) the closed loop characteristic equation of the controlled system can be obtained:

( )0

*CH s G G s NZ CH G s*

*T TT s T TG T T T s*

 

*K GG T K bG*

The stability of the equilibrium point for the system given in Fig. 18 can be determined from the characteristic polynomial of the closed loop harmonically linearised equation (Kuljaca et

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

(9)

(10)

*s NZ NZ*

*sT sT s s*

*s G*

3 2

Fig. 23. Describing function of fuzzy controller in two dimensional representation

GN(Xm) = GNF = GN(Fm) - describing function of the fuzzy regulator ,

*T T TG s G K bG KG*

*NZ o l NF s NZ*

By applying Hurwitz stability criterion the following inequality is obtained:

*CH G s NZ*

<sup>0</sup>

*CH l NF*

1

where

where: 0

al., Sijak et al., 2007):

<sup>1</sup> *<sup>K</sup> R*

1.1

1.2

1.3

1.4

1.5

Describing function gain

1.6

1.7

1.8

1.9

2

The parameters of the system are: TG = 0.08 s, TCH = 0.3 s, = 0.0017, Ks = 120 Hz/p.u., Ts = 20 s, R = 2.4 Hz/p.u., bl = 1. With given parameters, describing function of generation rate constraint is shown in Fig. 22.

Fuzzy controller represents a bit more complex problem. We do not have its describing function in analytical form and we need to obtain it by experimental method. The method is described in Section 5 and describing function gain is plotted in Fig. 16 in three dimensional representation. Phase changes are zero, thus there is no complex component of describing function.

Representation in three dimension in this case is not required since there is no dependency of describing function gain on frequency. After running simulation script given in Section 5, one should run the following Matlab code in order to obtain two-dimensional representation of fuzzy controller describing function:

```
gdf=s(:,1); 
plot(E,real(gdf)),grid, ylabel('Describing function gain'); 
xlabel('Amplitude');
```
Describing function in two dimensional representation is given in Fig. 23.

Fig. 22. Generation rate constraint describing function

The parameters of the system are: TG = 0.08 s, TCH = 0.3 s, = 0.0017, Ks = 120 Hz/p.u., Ts = 20 s, R = 2.4 Hz/p.u., bl = 1. With given parameters, describing function of generation rate

Fuzzy controller represents a bit more complex problem. We do not have its describing function in analytical form and we need to obtain it by experimental method. The method is described in Section 5 and describing function gain is plotted in Fig. 16 in three dimensional representation.

Representation in three dimension in this case is not required since there is no dependency of describing function gain on frequency. After running simulation script given in Section 5, one should run the following Matlab code in order to obtain two-dimensional

0 0.5 1 1.5 2 2.5 3

Amplitude

Phase changes are zero, thus there is no complex component of describing function.

plot(E,real(gdf)),grid, ylabel('Describing function gain');

Describing function in two dimensional representation is given in Fig. 23.

representation of fuzzy controller describing function:

constraint is shown in Fig. 22.

gdf=s(:,1);

xlabel('Amplitude');

0

Fig. 22. Generation rate constraint describing function

0.005

Describing function gain

0.01

0.015

Fig. 23. Describing function of fuzzy controller in two dimensional representation

The static characteristics of the fuzzy regulator, i.e. its describing functions, shows that harmonically linearized fuzzy regulator can be regarded as input dependent variable gain element. Consequently, the stability analysis can be conducted in the plane GNF = f(GNZ), where


The stability of the equilibrium point for the system given in Fig. 18 can be determined from the characteristic polynomial of the closed loop harmonically linearised equation (Kuljaca et al., Sijak et al., 2007):

$$\left(T\_{\rm CH} + \left(K\_0 + b\_I G\_{\rm NF}\right) \cdot \frac{K\_s}{1 + s \cdot T\_s} \cdot \frac{1}{1 + s \cdot T\_G} \cdot \frac{G\_{\rm NZ}}{s} + \frac{G\_{\rm NZ}}{s} = 0\right) \tag{9}$$

where: 0 <sup>1</sup> *<sup>K</sup> R*

From (9) the closed loop characteristic equation of the controlled system can be obtained:

$$\begin{aligned} &T\_{CH}T\_sT\_GS^3 + \left[T\_GS\_s\mathbf{G}\_{NZ} + T\_{CH}\left(T\_G + T\right)\_s\right]s^2 + \\ &\left[T\_{CH} + \left(T\_G + T\_s\right)\mathbf{G}\_{NZ}\right]s \\ &+ \mathbf{G}\_{NZ} + \left(\mathbf{K}\_o + b\_I\mathbf{G}\_{NF}\right)\mathbf{K}\_s\mathbf{G}\_{NZ} = \mathbf{0} \end{aligned} \tag{10}$$

By applying Hurwitz stability criterion the following inequality is obtained:

Describing Function Recording with Simulink and MATLAB 165

for obtaining describing function of any single input single output non-dynamic nonlinearity. Method used for conventional static nonlinearities is then adjusted to work with fuzzy systems (and other systems that cannot handle vector based batch simulation in Simulink). Matlab scripts required to run given Simulink model are also given in Chapter. This method can be also used to obtain describing function for single input single output static neural networks. Finally, an example from use describing function for

Described method can be used with practical controllers with Matlab real time tools.

*Systems (in Russian)*, Masinostroenie, Moskow, Russia

Dubrovnik, Croatia, June 2005

0813-X, Warsaw, Poland, September 2007.

Engelwood Cliffs, New Jersey, USA

Kreyszig, E. (1993). *Advanced Engineering Mathematics (7)*, John Wiley and Sons, New York,

Novogranov, B.N. (1986). *Determination of Frequency Domain Characteristics of Nonlinear* 

Kuljaca, O., Tesnjak, S., Vukić, Z. (1999). Describing Function of Mamdani Type Fuzzy

Sijak, T., Kuljaca, O., Antonic, A., Kuljaca, Lj., (2005). Analytical Determination of Describing

Sijak, T., Kuljaca, O., Kuljaca, Lj. (2004). Describing function of generalized static

Sijak, T., Kuljaca, O., Kuljaca, Lj., (2007). Computer Aided Harmonic Linearization of SISO

Sijak, T., Kuljaca, O., Tesnjak, S. (2001). Stability analysis of fuzzy control system using

Slotine, J.J.E., Li, W. (1991). *Applied Nonlinear Control*, Prentice Hall, ISBN 0-13-040890-5,

Schwartz, C., Gran, R. (2001). Describing Function Analysis Using MATLAB and Simulink.

*and Automation*, ISBN 953-6037-35-1, Dubrovnik Croatia, June2001

*on Robotics and Control,* Volume I, HKSRC '99, Hong Kong, 1999.

Controller with Input Signals Derived From Single System Input and Singleton Ouput Membership Functions, *Proceedings of the 1999 IEEE Hong Kong Symposium* 

Function of Nonlinear Element with Fuzzy Logic, *Proceedings of the IEEE International Symposium on Industrial Electronics 2005*, ISBN 0-7803-8739-2,

characteristic of nonlinear element, *Proceedings of REDISCOVER 2004, Southeastern Europe, USA, Japan and European Community Workshop on Research and Education in Control and Signal Processing*, ISBN 953-184-077-6, Cavtat, Croatia, June 2004 Sijak, T., Kuljaca, O., Kuljaca, Lj., (2007). Engineering Procedure for Analysis of Nonlinear

Structure Consisting of Fuzzy Element and Typical Nonlinear Element, *Proceedings of the 15th Mediterranean Conference on Control and Automation*, Athens, Greece, July

Systems Using Linearly Approximated Static Characteristic. *Proceedings of EUROCON 2007 The International Conference on Computer as a Tool*. ISBN 1-4244-

describing function method, *Proceedings of 9th Mediterranean Conference on Control* 

*IEEE Control Systems Magazine*, Vol. 21, No. 4, (August 2001), pp. 19-26, ISSN 1066-

stability analysis is given.

**7. References** 

USA

2007

033X

$$\begin{aligned} \left[T\_G T\_s G\_{NZ} + T\_{CH} \left(T\_G + T\_s\right)\right] \left[T\_{CH} + \left(T\_G + T\_s\right)G\_{NZ}\right] -\\ T\_{CH} T\_s T\_G \left[G\_{NZ} + \left(K\_0 + b\_l G\_{NF}\right)K\_s G\_{NZ}\right] > 0 \end{aligned} \tag{11}$$

The system is stable as long as (11) is positive. Thus, the stability boundary can be derived as:

$$\begin{split} \mathbf{G}\_{NF} &= \frac{\left(T\_{G} + T\_{S}\right)\mathbf{G}\_{NZ}}{b\_{l}T\_{CH}\mathbf{K}\_{s}} + \frac{T\_{CH}\left(T\_{G} + T\_{s}\right)^{2} - T\_{G}T\_{s}T\_{CH}\mathbf{K}\_{0}\mathbf{K}\_{s}}{b\_{l}T\_{CH}T\_{s}T\_{G}\mathbf{K}\_{s}} \\ &+ \frac{T\_{CH}\left(T\_{G} + T\_{s}\right)}{b\_{l}T\_{s}T\_{G}\mathbf{K}\_{s}\mathbf{G}\_{NZ}} \end{split} \tag{12}$$

With the parameters used for simulation, the stability boundary function GNF = f(GNZ) can be evaluated. The function GNF = f(GNZ) is shown in Fig. 24.

Fig. 24. Function GNF = f(GNZ)

With given values system is stable in operating region of interest.

#### **6. Conclusion**

Chapter deals with use of Simulink for obtaining describing function for different elements. Described method of numerical simulation performed in Simulink is suitable for obtaining describing function of any single input single output non-dynamic nonlinearity. Method used for conventional static nonlinearities is then adjusted to work with fuzzy systems (and other systems that cannot handle vector based batch simulation in Simulink). Matlab scripts required to run given Simulink model are also given in Chapter. This method can be also used to obtain describing function for single input single output static neural networks. Finally, an example from use describing function for stability analysis is given.

Described method can be used with practical controllers with Matlab real time tools.

## **7. References**

164 Technology and Engineering Applications of Simulink

 <sup>0</sup> ( )0 *G s NZ CH G s CH G s NZ*

The system is stable as long as (11) is positive. Thus, the stability boundary can be derived as:

*l CH s l CH s G s*

With the parameters used for simulation, the stability boundary function GNF = f(GNZ) can be

0 0.1 0.2 0.3 0.4 0.5

GNZ

Chapter deals with use of Simulink for obtaining describing function for different elements. Described method of numerical simulation performed in Simulink is suitable

With given values system is stable in operating region of interest.

unstable region

*bT K bT TT K*

 

*T TG T T T T T T G*

*CH s G NZ l NF s NZ*

*T TT G K bG KG*

*CH G s l s G s NZ*

*TTT bTT K G*

0

Fig. 24. Function GNF = f(GNZ)

**6. Conclusion** 

5

10

15

G

NF

20

25

30

*G*

*NF*

evaluated. The function GNF = f(GNZ) is shown in Fig. 24.

*G S NZ CH G s G s CH s*

*T T G T T T T TT K K*

2

0

(11)

(12)


**8** 

**Performance Evaluation of a Temperature** 

E.N. Vázquez-Acosta1, S. Mendoza-Acevedo1, M.A. Reyes-Barranca1, L.M. Flores-Nava1, J.A. Moreno-Cadenas1 and J.L. González-Vidal2

Actually, the use of CAD-CAE is a widespread activity for design and evaluation purposes. Hence, several programs can be found for simulation and optimization of systems described with 3D graphic representation, giving the morphology of the model (Muñoz et al., 2008). Also, suites can be found that can simulate and prove some kind of controls being described with transfer functions (TFs), differential equations (DE), space states or ASIC models

In general, the design begins creating the file of the solids that describe the geometry of the device through CAD tools and then exported to simulation and analysis environments (multi-physics). Thereafter, the results obtained are used in a multi-domain platform, where more complete results can be obtained for a system in consideration, taking advantage of modules that can link these environments. Now, one should have in mind that if complex dynamic models are to be considered, a serious limitation can be present if two or more variables of the system are mutually dependant. Therefore, it is suggested to export static models. In this case, the role of the modules is centred only to a basic interchange of data between the environments, having not a real and full interaction between the models considered. Then as a recommendation, it is convenient to include all the possible dynamics in the exported model, even as a part of it or as an external issue, such that the dynamic simulation can give more precise results. On the other hand, when a dynamic simulation is done based on a static model, there could be a dramatic increase in simulation time since a call to the static exported model will be done at each simulation time step. For instance, if a simulation is specified for a 1s lapse with 1ms steps, the exported model will be called 1000

The particular case exposed in this paper is the analysis of a coupled thermo-electrical model designed to operate a micro hotplate (MHP) with a temperature control circuit, used

**1. Introduction** 

(Goering, 2004).

times.

**Control Stage Used on a Semiconductor** 

**Gas Sensor 3D Electro-Thermal Model** 

*1CINVESTAV-IPN, Electrical Engineering Department* 

*2Universidad Autónoma Del Estado de Hidalgo,* 

**Through Simulink®**

*Computing Academic Area* 

*Mexico* 

Vukic, Z., Kuljaca, Lj., Donlagic, D., Tesnjak, S. (2003). *Nonlinear Control Systems*, Marcel Dekker, ISBN 0-8247-4112-9, New York, USA

Netushil at al. *Theory of Automatic Control*, Mir, Moscow 1978

## **Performance Evaluation of a Temperature Control Stage Used on a Semiconductor Gas Sensor 3D Electro-Thermal Model Through Simulink®**

E.N. Vázquez-Acosta1, S. Mendoza-Acevedo1, M.A. Reyes-Barranca1, L.M. Flores-Nava1, J.A. Moreno-Cadenas1 and J.L. González-Vidal2 *1CINVESTAV-IPN, Electrical Engineering Department 2Universidad Autónoma Del Estado de Hidalgo, Computing Academic Area Mexico* 

## **1. Introduction**

166 Technology and Engineering Applications of Simulink

Vukic, Z., Kuljaca, Lj., Donlagic, D., Tesnjak, S. (2003). *Nonlinear Control Systems*, Marcel

Dekker, ISBN 0-8247-4112-9, New York, USA Netushil at al. *Theory of Automatic Control*, Mir, Moscow 1978

> Actually, the use of CAD-CAE is a widespread activity for design and evaluation purposes. Hence, several programs can be found for simulation and optimization of systems described with 3D graphic representation, giving the morphology of the model (Muñoz et al., 2008). Also, suites can be found that can simulate and prove some kind of controls being described with transfer functions (TFs), differential equations (DE), space states or ASIC models (Goering, 2004).

> In general, the design begins creating the file of the solids that describe the geometry of the device through CAD tools and then exported to simulation and analysis environments (multi-physics). Thereafter, the results obtained are used in a multi-domain platform, where more complete results can be obtained for a system in consideration, taking advantage of modules that can link these environments. Now, one should have in mind that if complex dynamic models are to be considered, a serious limitation can be present if two or more variables of the system are mutually dependant. Therefore, it is suggested to export static models. In this case, the role of the modules is centred only to a basic interchange of data between the environments, having not a real and full interaction between the models considered. Then as a recommendation, it is convenient to include all the possible dynamics in the exported model, even as a part of it or as an external issue, such that the dynamic simulation can give more precise results. On the other hand, when a dynamic simulation is done based on a static model, there could be a dramatic increase in simulation time since a call to the static exported model will be done at each simulation time step. For instance, if a simulation is specified for a 1s lapse with 1ms steps, the exported model will be called 1000 times.

> The particular case exposed in this paper is the analysis of a coupled thermo-electrical model designed to operate a micro hotplate (MHP) with a temperature control circuit, used

Performance Evaluation of a Temperature Control Stage Used on a

Fig. 2. Layout of the integrated gas sensor system.

membrane.

shown in Fig. 3.

Semiconductor Gas Sensor 3D Electro-Thermal Model Through Simulink® 169

This kind of gas sensor can be fabricated starting with a silicon substrate and includes a micro hotplate (MHP), whose structure stands as a mechanical support for the layers contained in the thin membrane (SiO2/Polysilicon/SiO2/Polysilicon/Metal Oxide). The thin film located above all these layers is a sensitive metal oxide layer that senses the gaseous species to be detected. Besides, the main purpose of the membrane is to isolate the heat generated in that zone and to uniformly distribute the temperature in the thin

The size of the MHP is 80µm x 80µm and will be suspended by four arms over the microcavity obtained after a bulk micromachining post-process usual in MEMS technology, giving an effective thermal isolation from the substrate. This way, the thermal mass of the MHP is reduced and as a consequence, its thermal inertia is decreased. Therefore, it is easy to attain an operation temperature range between 250°C and 350°C with low voltages applied (< 2.5V). This temperature promotes the chemical reaction between the thin metal oxide layer and the gas species (Afridi et al., 2002). Also, the MHP contains a micro temperature sensor (TS) and a microheater (MH) that is biased to raise its temperature due to the Joule Effect. Both will be fabricated with polysilicon and the proposed configuration is

in an integrated semiconductor gas sensor (SGS). The analysis should consider that the MHP will be physically located in a thin membrane created with MEMS technology. The complexity of this system is rather high since a thermo-electrical coupling must exist among devices that operate as a heat source and the different layers affected by the generated heat. Specifically, the case here studied treats the Joule Effect applied to the polysilicon material used in integrated circuits. Thus, it is necessary to combine somehow the analysis programs to obtain the required electro-thermal simulation results and to verify if the expected performance, i.e. temperature variation in time, is correct.

## **2. Description of the system**

A SGS was designed based on layers available in standard CMOS technology, having two polysilicon layers and two metal layers. The layout of the prototype is shown in Fig. 1.

Fig. 1. Layout of the gas sensor prototype.

One advantage in using this kind of technology is its compatibility with MEMS micromachining processes made with anisotropic etching solutions as tetramethylammonium hydroxide (TMAH), resulting in a preferential etching of the silicon substrate. This way, a very thin membrane is obtained containing the MHP (Fedder, 2005). The layout of the integrated sensor circuit is shown in Fig. 2, where the different blocks implemented for the system, are identified.

in an integrated semiconductor gas sensor (SGS). The analysis should consider that the MHP will be physically located in a thin membrane created with MEMS technology. The complexity of this system is rather high since a thermo-electrical coupling must exist among devices that operate as a heat source and the different layers affected by the generated heat. Specifically, the case here studied treats the Joule Effect applied to the polysilicon material used in integrated circuits. Thus, it is necessary to combine somehow the analysis programs to obtain the required electro-thermal simulation results and to verify if the expected

A SGS was designed based on layers available in standard CMOS technology, having two polysilicon layers and two metal layers. The layout of the prototype is shown in Fig. 1.

One advantage in using this kind of technology is its compatibility with MEMS micromachining processes made with anisotropic etching solutions as tetramethylammonium hydroxide (TMAH), resulting in a preferential etching of the silicon substrate. This way, a very thin membrane is obtained containing the MHP (Fedder, 2005). The layout of the integrated sensor circuit is shown in Fig. 2, where the different blocks implemented

performance, i.e. temperature variation in time, is correct.

**2. Description of the system** 

Fig. 1. Layout of the gas sensor prototype.

for the system, are identified.

Fig. 2. Layout of the integrated gas sensor system.

This kind of gas sensor can be fabricated starting with a silicon substrate and includes a micro hotplate (MHP), whose structure stands as a mechanical support for the layers contained in the thin membrane (SiO2/Polysilicon/SiO2/Polysilicon/Metal Oxide). The thin film located above all these layers is a sensitive metal oxide layer that senses the gaseous species to be detected. Besides, the main purpose of the membrane is to isolate the heat generated in that zone and to uniformly distribute the temperature in the thin membrane.

The size of the MHP is 80µm x 80µm and will be suspended by four arms over the microcavity obtained after a bulk micromachining post-process usual in MEMS technology, giving an effective thermal isolation from the substrate. This way, the thermal mass of the MHP is reduced and as a consequence, its thermal inertia is decreased. Therefore, it is easy to attain an operation temperature range between 250°C and 350°C with low voltages applied (< 2.5V). This temperature promotes the chemical reaction between the thin metal oxide layer and the gas species (Afridi et al., 2002). Also, the MHP contains a micro temperature sensor (TS) and a microheater (MH) that is biased to raise its temperature due to the Joule Effect. Both will be fabricated with polysilicon and the proposed configuration is shown in Fig. 3.

Performance Evaluation of a Temperature Control Stage Used on a

Fig. 4. MH and MS resistance value as a function of temperature.

following sections establish the proposed methodology.

transfer by conduction across the silicon dioxide between them.

**3. Electro-thermal model** 

**3.1 Simulation based on finite elements** 

Semiconductor Gas Sensor 3D Electro-Thermal Model Through Simulink® 171

Based on the later and the geometry drawn (Fig. 3), then the system will have a particular thermal response time for heating and cooling of the MHP. This is precisely the performance that is desired to be evaluated. Given the above background, the objectives of this work are: a) to propose a simulation and analysis strategy to find the dynamic behaviour of the model of the MHP driven by a temperature control circuit, based on the results obtained from the electro-thermal analysis, and b) to optimize the computing resources and simulation time. If the results are the expected, then the system can be fabricated with confidence. The

Using specialized CAD software, a 3D solid model was created in order to propose the geometry of the structure needed for the electro-thermal analysis and simulation. This geometry should include the microcavity that is supposed to be made with the anisotropic etching after the fabrication of the chip. Once the geometry is finished, it can be exported to a multi-physics platform and proper sub-domains were assigned as well as boundary conditions (DC conductive media and heat transfer) to perform the analysis of the electrothermal model. Also, to assure convergence with the finite elements analysis, a convenient mesh density should be selected. The analysis was made using an applied voltage sweep from 1V to 3V with 0.2V steps. The result is shown in Fig. 5, where the electro-thermal behaviour of the 3D geometric model of the gas sensor can be seen after the analysis with finite elements. Here it is confirmed that the heat is focused at the MHP as a result of the micromachining specified before. It should be noted that although the TS is not biased, it will have the same temperature induced by Joule Effect in the MH as a consequence of heat

Fig. 3. Geometry of the microheater resistor and the temperature sensor.

When implemented as a resistor, polysilicon increases its resistance value proportionally with an associated temperature rise. This property is named Thermal Coefficient of Resistance (TCR) (Cerda et al., 2006). The relation of the resistance value with a temperature change is expressed with the following equation:

$$\frac{R - R\_0}{R\_0} = TCR \left( T - T\_0 \right) \tag{1}$$

where T0 is room temperature and R0 is the resistance value at room temperature. TCR and R0 can be experimentally determined. If the relation between temperature and resistance is known, it can be conveniently used to design a temperature control circuit if a constant current is passed through TS. If a DC voltage is applied across the microheater terminals, the resistance value of polysilicon increases as predicted by equation (1). Fig. 4 shows the variation of the resistance of TS and MH as a function of temperature from where TCR can be obtained.

TCR property sets a problem hard to solve from the simulation point of view when using traditional methods or programs like SPICE, since there is no possibility to consider at the same time the Joule Effect and the electro-thermal coupling between MH and TS. Also, it is difficult to properly consider the thermal inertia of the MHP's thermal mass using electrical simulators.

From Fig. 3, it is easy to deduce that the thermal coupling between the temperature sensor and the microheater is possible because of their proximity, so the fact that TS follows the temperature changes in MH can be used favourably. As it is well known, the thermal coupling behaviour depends on the properties of the materials used. For the structure here analysed, the materials involved are crystalline silicon, polysilicon, silicon dioxide and aluminium.

Fig. 4. MH and MS resistance value as a function of temperature.

Based on the later and the geometry drawn (Fig. 3), then the system will have a particular thermal response time for heating and cooling of the MHP. This is precisely the performance that is desired to be evaluated. Given the above background, the objectives of this work are: a) to propose a simulation and analysis strategy to find the dynamic behaviour of the model of the MHP driven by a temperature control circuit, based on the results obtained from the electro-thermal analysis, and b) to optimize the computing resources and simulation time. If the results are the expected, then the system can be fabricated with confidence. The following sections establish the proposed methodology.

## **3. Electro-thermal model**

170 Technology and Engineering Applications of Simulink

Fig. 3. Geometry of the microheater resistor and the temperature sensor.

0

*R*

change is expressed with the following equation:

be obtained.

simulators.

aluminium.

When implemented as a resistor, polysilicon increases its resistance value proportionally with an associated temperature rise. This property is named Thermal Coefficient of Resistance (TCR) (Cerda et al., 2006). The relation of the resistance value with a temperature

<sup>0</sup>

*R R TCR T T*

where T0 is room temperature and R0 is the resistance value at room temperature. TCR and R0 can be experimentally determined. If the relation between temperature and resistance is known, it can be conveniently used to design a temperature control circuit if a constant current is passed through TS. If a DC voltage is applied across the microheater terminals, the resistance value of polysilicon increases as predicted by equation (1). Fig. 4 shows the variation of the resistance of TS and MH as a function of temperature from where TCR can

TCR property sets a problem hard to solve from the simulation point of view when using traditional methods or programs like SPICE, since there is no possibility to consider at the same time the Joule Effect and the electro-thermal coupling between MH and TS. Also, it is difficult to properly consider the thermal inertia of the MHP's thermal mass using electrical

From Fig. 3, it is easy to deduce that the thermal coupling between the temperature sensor and the microheater is possible because of their proximity, so the fact that TS follows the temperature changes in MH can be used favourably. As it is well known, the thermal coupling behaviour depends on the properties of the materials used. For the structure here analysed, the materials involved are crystalline silicon, polysilicon, silicon dioxide and

0

(1)

## **3.1 Simulation based on finite elements**

Using specialized CAD software, a 3D solid model was created in order to propose the geometry of the structure needed for the electro-thermal analysis and simulation. This geometry should include the microcavity that is supposed to be made with the anisotropic etching after the fabrication of the chip. Once the geometry is finished, it can be exported to a multi-physics platform and proper sub-domains were assigned as well as boundary conditions (DC conductive media and heat transfer) to perform the analysis of the electrothermal model. Also, to assure convergence with the finite elements analysis, a convenient mesh density should be selected. The analysis was made using an applied voltage sweep from 1V to 3V with 0.2V steps. The result is shown in Fig. 5, where the electro-thermal behaviour of the 3D geometric model of the gas sensor can be seen after the analysis with finite elements. Here it is confirmed that the heat is focused at the MHP as a result of the micromachining specified before. It should be noted that although the TS is not biased, it will have the same temperature induced by Joule Effect in the MH as a consequence of heat transfer by conduction across the silicon dioxide between them.

Performance Evaluation of a Temperature Control Stage Used on a

Fig. 7. MH cooling response with applied voltage as a parameter.

**3.2 Simulation of the static model with the transfer function at the output** 

typical TCR of polysilicon can be easily calculated from this temperature reading.

no over damping, using a first order TF will simplify the task (Ogata, 1999).

The model created in the multi-physics environment can be exported to a multi-domain environment using the existing modules specially created to link different software suits. Once the model is exported, one arbitrary exciting input and several reading outputs can be declared. In the case here exposed, the input is described as the voltage applied to MH. The value (or even the waveform) of this voltage can be modified along the simulation process in the multi-domain software. This is not possible to do with the multi-physics program. One of the outputs is declared as the reading of the temperature reached by MHP. The

Regardless of the possibilities and advantages that linking modules have, the model presents limitations when it is exported due to the mutual dependence of the variables as described by the electro-thermal coupling, therefore it is exported only as a static model.

Due to the above reason, transfer functions (TFs) must be added at the output of the simulation block of the exported system to consider dynamic features of the model. This makes possible to emulate the dynamic behaviour while heating or cooling, following the results obtained in Figs. 6 and 7. From now on, this model will be named as the exported model. For the case here reported, the TFs of interest can be obtained after making a series of heating and cooling cycles using the voltage range specified in Figs. 6 and 7. With the results obtained, an expression that nearly fits with the heating and cooling curves can be extrapolated. This can be done either with computational methods or by trial and error. Selection of the order of the expression of the TFs highly depends on the response of the system to a step input signal and to the desired accuracy. Since the system here analysed has

Semiconductor Gas Sensor 3D Electro-Thermal Model Through Simulink® 173

Fig. 5. 3D exported model from AutoCAD® after simulated with COMSOL Multiphysics®.

Using this same platform (multi-physics), the thermal response as a function of time having the applied voltage as a parameter was also obtained. Figs. 6 and 7 show the results for heating and cooling, respectively and it should be mentioned that a different response will be obtained for different geometries of MH and TS. This behaviour follows the change in resistance of polysilicon predicted by equation (1).

Fig. 6. MH heating response with applied voltage as a parameter.

Fig. 5. 3D exported model from AutoCAD® after simulated with COMSOL Multiphysics®.

resistance of polysilicon predicted by equation (1).

Fig. 6. MH heating response with applied voltage as a parameter.

Using this same platform (multi-physics), the thermal response as a function of time having the applied voltage as a parameter was also obtained. Figs. 6 and 7 show the results for heating and cooling, respectively and it should be mentioned that a different response will be obtained for different geometries of MH and TS. This behaviour follows the change in

Fig. 7. MH cooling response with applied voltage as a parameter.

## **3.2 Simulation of the static model with the transfer function at the output**

The model created in the multi-physics environment can be exported to a multi-domain environment using the existing modules specially created to link different software suits. Once the model is exported, one arbitrary exciting input and several reading outputs can be declared. In the case here exposed, the input is described as the voltage applied to MH. The value (or even the waveform) of this voltage can be modified along the simulation process in the multi-domain software. This is not possible to do with the multi-physics program. One of the outputs is declared as the reading of the temperature reached by MHP. The typical TCR of polysilicon can be easily calculated from this temperature reading.

Regardless of the possibilities and advantages that linking modules have, the model presents limitations when it is exported due to the mutual dependence of the variables as described by the electro-thermal coupling, therefore it is exported only as a static model.

Due to the above reason, transfer functions (TFs) must be added at the output of the simulation block of the exported system to consider dynamic features of the model. This makes possible to emulate the dynamic behaviour while heating or cooling, following the results obtained in Figs. 6 and 7. From now on, this model will be named as the exported model. For the case here reported, the TFs of interest can be obtained after making a series of heating and cooling cycles using the voltage range specified in Figs. 6 and 7. With the results obtained, an expression that nearly fits with the heating and cooling curves can be extrapolated. This can be done either with computational methods or by trial and error. Selection of the order of the expression of the TFs highly depends on the response of the system to a step input signal and to the desired accuracy. Since the system here analysed has no over damping, using a first order TF will simplify the task (Ogata, 1999).

Performance Evaluation of a Temperature Control Stage Used on a

before, and finally a block that displays the output waveform (Scope).

Fig. 9. Block diagram of Simulink® for the pulsed simulation.

Semiconductor Gas Sensor 3D Electro-Thermal Model Through Simulink® 175

change at the input signal is detected, saving this way time and computational resources. Here, existing modules that detect signal changes (Detect Change) are used together with a Boolean converter and a module that evaluates the exported model (COMSOL Triggered Simulation) followed by a Hold block that keeps the output value until there is a change at the input signal. All these modules are contained in the library of the multi-domain platform. Other supporting blocks included in Simulink® are a temperature unit converter (from K to °C), the Dynamic Adjust block that includes the heating and cooling TFs deduced

The simplest expressions that fit approximately with the dynamics of the original system are the following:


$$G\_c(\mathbf{s}) = \frac{1}{0.00077\mathbf{s} + 1} \tag{2}$$


$$G\_{\rm fl}(\mathbf{s}) = \frac{1}{0.000095\mathbf{s} + 1} \tag{3}$$

Fig. 8 shows the result after the simulation of the static model, together with the deduced TFs output when a voltage pulse was applied to the MH.

Fig. 8. Output of the static model when excited with an input step signal (continuous line) and the transfer function deduced from Figs. 6 and 7 (dashed line).

Once these TFs are included it is possible to apply arbitrary inputs to the system and read their corresponding dynamic outputs storing them as variables of the multi-domain software, with Simulink® of Matlab®, for instance. Also, it is possible to add and simulate some kind of temperature control block together with a signal processing or ASIC circuit. The later is not possible to do within the environment from where the created model was exported. Using the approximation with TFs it is possible to use certain blocks that call the exported model only when there are abrupt changes in the system's input. This helps to decrease the computing processing time (Hsu, 2002). Therefore, the suggested block diagram of the model that can be simulated with this alternative is shown in Fig. 9. This figure shows the block of the simulated electro-thermal system (COMSOL Subsystem), a signal generator block (Signal Builder), a block that calls the exported static model simulation (Pulsed COMSOL Sim), having the task to evaluate the model only when a

The simplest expressions that fit approximately with the dynamics of the original system are

<sup>1</sup> ( ) 0.00077 1 *G s <sup>c</sup> <sup>s</sup>*

<sup>1</sup> ( ) 0.000095 1 *G s <sup>h</sup> <sup>s</sup>*

Fig. 8 shows the result after the simulation of the static model, together with the deduced

Fig. 8. Output of the static model when excited with an input step signal (continuous line)

Once these TFs are included it is possible to apply arbitrary inputs to the system and read their corresponding dynamic outputs storing them as variables of the multi-domain software, with Simulink® of Matlab®, for instance. Also, it is possible to add and simulate some kind of temperature control block together with a signal processing or ASIC circuit. The later is not possible to do within the environment from where the created model was exported. Using the approximation with TFs it is possible to use certain blocks that call the exported model only when there are abrupt changes in the system's input. This helps to decrease the computing processing time (Hsu, 2002). Therefore, the suggested block diagram of the model that can be simulated with this alternative is shown in Fig. 9. This figure shows the block of the simulated electro-thermal system (COMSOL Subsystem), a signal generator block (Signal Builder), a block that calls the exported static model simulation (Pulsed COMSOL Sim), having the task to evaluate the model only when a

and the transfer function deduced from Figs. 6 and 7 (dashed line).

TFs output when a voltage pulse was applied to the MH.

(2)

(3)

the following: - Cooling TF:


change at the input signal is detected, saving this way time and computational resources. Here, existing modules that detect signal changes (Detect Change) are used together with a Boolean converter and a module that evaluates the exported model (COMSOL Triggered Simulation) followed by a Hold block that keeps the output value until there is a change at the input signal. All these modules are contained in the library of the multi-domain platform. Other supporting blocks included in Simulink® are a temperature unit converter (from K to °C), the Dynamic Adjust block that includes the heating and cooling TFs deduced before, and finally a block that displays the output waveform (Scope).

Fig. 9. Block diagram of Simulink® for the pulsed simulation.

Performance Evaluation of a Temperature Control Stage Used on a

Fig. 11. Model with the approximated equation.

arbitrary pulse train as input.

Semiconductor Gas Sensor 3D Electro-Thermal Model Through Simulink® 177

(corresponding to a voltage range of 2V ± 0.42V) was used. The first test was done using a pulse

Fig. 12. Results from the simulation using the three analysing alternatives proposed, with an

As shown in Fig. 12, the three curves are quite similar among them. Deviation of methods 3.2 and 3.3 compared with 3.1 is only 0.22% in the worst case, meaning a variation of 1°C each 461.4°C. Hence it can be concluded that these options give similar results such that each of them can be used with confidence in the voltage range studied. However, it should be remembered that the exported model is limited to tests with simple on-off temperature controls and/or for abrupt inputs or with almost no changes in time. Otherwise, if inputs like a ramp or a sine waveforms are used, for instance, the static model will be invoked at each

train having arbitrary widths and magnitudes. The results are shown in Fig. 12.

#### **3.3 Simulation with the approximated equation and transfer functions**

The model illustrated in Fig. 9 presents great performance advantages if the applied inputs present abrupt changes in time like step signals. However it could be extremely slow if input signals have no abrupt changes during time like ramps or sine waveforms. So, if signals continuously changing in time are to be used, a more useful method to apply is based on the substitution of the block that calls the model (illustrated in Fig. 9), with an approximated equation that closely follows the temperature dependence of MH to the applied DC voltage, as obtained in the multi-physics platform. The kind of equation can be arbitrarily chosen while the parameters of it should be adjusted using the data obtained from the simulation of the exported model. Again, this can be done with the help of computational regression methods or by trial and error. Care should be taken to assure that the equation used is valid at least within the voltage sweep range of interest. The equation found for the case here presented is the following:

$$f(\mathbf{x}) = 64\mathbf{x}^{\prime(1.938 - 0.0605\mathbf{x})} + 20\tag{4}$$

where f(x) is the temperature output in K and x is the voltage applied to the MH in the exported model. Fig. 10 shows good agreement between the curve obtained after the static simulation of the exported model and the plot of equation (4).

Fig. 10. Comparison between the static model and Eq. (4).

Therefore, this different alternative can be used properly in the desired dynamic simulation of the model. Once the approximated equation is obtained, it can be integrated into the simulation model instead of the module that invokes the pulsed exported model, using a function block as is shown in Fig. 11. Then, three options for the model simulation were presented: the exported model (3.1), the exported static model (3.2) and the approximated model (3.3). A series of simulations can be done to compare performance, validate the models and make a convenient choice among these methods. At first glance, a MH temperature range of 172.6 °C – 327.8°C (corresponding to a voltage range of 2V ± 0.42V) was used. The first test was done using a pulse train having arbitrary widths and magnitudes. The results are shown in Fig. 12.

Fig. 11. Model with the approximated equation.

176 Technology and Engineering Applications of Simulink

The model illustrated in Fig. 9 presents great performance advantages if the applied inputs present abrupt changes in time like step signals. However it could be extremely slow if input signals have no abrupt changes during time like ramps or sine waveforms. So, if signals continuously changing in time are to be used, a more useful method to apply is based on the substitution of the block that calls the model (illustrated in Fig. 9), with an approximated equation that closely follows the temperature dependence of MH to the applied DC voltage, as obtained in the multi-physics platform. The kind of equation can be arbitrarily chosen while the parameters of it should be adjusted using the data obtained from the simulation of the exported model. Again, this can be done with the help of computational regression methods or by trial and error. Care should be taken to assure that the equation used is valid at least within the voltage sweep range of interest. The equation

where f(x) is the temperature output in K and x is the voltage applied to the MH in the exported model. Fig. 10 shows good agreement between the curve obtained after the static

Therefore, this different alternative can be used properly in the desired dynamic simulation of the model. Once the approximated equation is obtained, it can be integrated into the simulation model instead of the module that invokes the pulsed exported model, using a function block as is shown in Fig. 11. Then, three options for the model simulation were presented: the exported model (3.1), the exported static model (3.2) and the approximated model (3.3). A series of simulations can be done to compare performance, validate the models and make a convenient choice among these methods. At first glance, a MH temperature range of 172.6 °C – 327.8°C

1.938 0.0605 ( ) 64 <sup>20</sup> *<sup>x</sup> fx x* (4)

**3.3 Simulation with the approximated equation and transfer functions** 

found for the case here presented is the following:

simulation of the exported model and the plot of equation (4).

Fig. 10. Comparison between the static model and Eq. (4).

Fig. 12. Results from the simulation using the three analysing alternatives proposed, with an arbitrary pulse train as input.

As shown in Fig. 12, the three curves are quite similar among them. Deviation of methods 3.2 and 3.3 compared with 3.1 is only 0.22% in the worst case, meaning a variation of 1°C each 461.4°C. Hence it can be concluded that these options give similar results such that each of them can be used with confidence in the voltage range studied. However, it should be remembered that the exported model is limited to tests with simple on-off temperature controls and/or for abrupt inputs or with almost no changes in time. Otherwise, if inputs like a ramp or a sine waveforms are used, for instance, the static model will be invoked at each

Performance Evaluation of a Temperature Control Stage Used on a

(see Fig. 11) begin at 0°C, as shown in Fig. 16.

Fig. 16. Stabilization time to room temperature.

Semiconductor Gas Sensor 3D Electro-Thermal Model Through Simulink® 179

Fig. 15. Comparative with a sine waveform with a 3V amplitude and 100Hz frequency.

Near to real operating conditions can be considered this way. Also, it should be noted that the initial temperature considered in the multi-physics simulation is 20°C, as established in the boundary conditions, but the simulation results obtained with the approximated model

So in order to compare at same temperature points, the multi-domain simulation is first allowed to stabilize at 20°C, then the input signal is applied. Therefore, the considered waveform is started after 0.01s although the temperature stabilization time is shorter as is

simulation step declared, increasing considerably the processing time. Since the objective is to optimize the simulation procedure, method 3.2 will be discarded from now on and only the original model (3.1) and the last one (3.3) will be considered in the following analyses. Comparative simulations between these two last methods were made, using the results of the original model obtained in the multi-physics environment and those obtained using the approximated model in the multi-domain platform. The results are shown in Figs. 13-15.

With these last results using arbitrary waveforms as input to the MH, heating and cooling can be effectively evaluated considering also heat dissipation due to the environment that is surrounding the system.

Fig. 13. Comparative with a saw tooth input signal.

Fig. 14. Comparative with a sine waveform with a 3V amplitude and 10Hz frequency.

simulation step declared, increasing considerably the processing time. Since the objective is to optimize the simulation procedure, method 3.2 will be discarded from now on and only the original model (3.1) and the last one (3.3) will be considered in the following analyses. Comparative simulations between these two last methods were made, using the results of the original model obtained in the multi-physics environment and those obtained using the approximated model in the multi-domain platform. The results are shown in Figs. 13-15.

With these last results using arbitrary waveforms as input to the MH, heating and cooling can be effectively evaluated considering also heat dissipation due to the environment that is

Fig. 14. Comparative with a sine waveform with a 3V amplitude and 10Hz frequency.

surrounding the system.

Fig. 13. Comparative with a saw tooth input signal.

Fig. 15. Comparative with a sine waveform with a 3V amplitude and 100Hz frequency.

Near to real operating conditions can be considered this way. Also, it should be noted that the initial temperature considered in the multi-physics simulation is 20°C, as established in the boundary conditions, but the simulation results obtained with the approximated model (see Fig. 11) begin at 0°C, as shown in Fig. 16.

Fig. 16. Stabilization time to room temperature.

So in order to compare at same temperature points, the multi-domain simulation is first allowed to stabilize at 20°C, then the input signal is applied. Therefore, the considered waveform is started after 0.01s although the temperature stabilization time is shorter as is

Performance Evaluation of a Temperature Control Stage Used on a

system.

Semiconductor Gas Sensor 3D Electro-Thermal Model Through Simulink® 181

From the mentioned electro-thermal characteristics of polysilicon, MH and TS can be used as part of a temperature control circuit. When the temperature sensor is connected as a feedback from the output to the input of the operational amplifier, the voltage drop in it is compared with a reference voltage (VREF). If the voltage difference between both inputs is zero, a power MOS transistor used as a switch connected at the output of the OP-AMP, will interrupt the biasing to the MH. A regulated current source applying pulses with constant amplitude should be used to bias TS (note that this does not cause Joule Effect to TS). Due to heat transmission by conduction, the resistance in the TS will increase (see Fig. 4) and its voltage drop will be modified. This property makes feasible to apply the TS in a temperature control circuit. On the other hand, the MH should be biased with a pulsed voltage source to heat the MHP. The magnitude of the pulse depends on the desired temperature. Besides, a pulsed stimulus will help to reduce the power consumption of the

Hence, when the reference voltage is higher than the feedback voltage drop, the output of the OP-AMP is high making the MOS transistor to conduct allowing a current flow through the MH. So, due to Joule Effect the temperature will gradually rise from room temperature up to a temperature established by the applied voltage, VREF. At the same time the TS will follow the MH's temperature, changing also the voltage drop on it. At the moment when both voltages are equal, the OP-AMP output will go low turning off the MOS transistor and heating will stop. This cycle will be repeated according to the heating and cooling dynamic behaviour of the microheater, operating around the desired temperature. It is expected that the MH will keep this temperature until the set point is modified if a different temperature target within the physical and mechanical limits of the MHP is desired. Next, to ease the dynamic analysis of the controlling circuit, an equivalent circuit was used to consider the thermal coupling between the heater and the sensor as well as the time constant of the electro-thermal system. The MH was replaced with a voltage controlled voltage source and the TS with an RC circuit, where R and C are the thermal resistance and thermal capacitance of polysilicon, respectively. The corresponding time constant of this circuit should be

consistent with the thermal behaviour of the 3D model of the MH considered.

temperature around the operating point as can be seen in Fig. 19.

From the finite element analysis made above, different time constants resulted for heating and for cooling (see Figs. 6 and 7). That is the reason why the equivalent circuit shown in Fig. 18 has two resistance elements corresponding to each value. The magnitude of R for heating is 770 Ω and 950 Ω for cooling. The capacitance value used to fit with each time constants is 1μF; the circuit was biased with VDD=1.5V. The transient behaviour of the system can be evaluated using the circuit shown in Fig. 18. The voltage drop in C is equivalent to the temperature reached by the system and is also the feedback voltage of the OP-AMP. Finally, the MOS transistor switch is replaced by two MOS transistors (PMOS and NMOS) to handle charge and discharge (heating and cooling, respectively) of the capacitor, as the output of the OP-AMP goes high or low. This is the reason why two resistors are used, so heating and cooling can be considered. The results of the transient analysis of the circuit show that there will be control over the current flow of the microheater until the voltage drop on TS reaches the reference voltage (1V for the case illustrated in Fig. 19). As soon as the voltage at the gate of the MOS transistor, VGS, equals the reference voltage, then biasing of the microheater stops. According to the cooling dynamics of the microheater, bias will be switching on and off to hold the

seen in Fig.16. Taking the temperature stabilization to times shorter than 0.001s may affect the simulation dynamics at the beginning.

Based on the above, it may be concluded that the outlined strategies have good performance and therefore they can be used to continue with the system design. Now the temperature control circuit will be included for the next analyses. So, after the validation of the models, the following step is to prove a temperature control, either using blocks from the library of the platform or using integration modules with other simulation programs. Once this is made, the designer can have a global idea about the performance of the system with the possibility to optimize the design.

## **4. Model of the temperature control system and simulation**

For the temperature control of the micro hotplate obtained with post-processes compatible with CMOS technology described before, an analogue proportional control based on an operational amplifier (OP-AMP) is proposed (Barrentino et al., 2004). The schematic used for the basic simulations of the system is shown in Fig. 17.

Fig. 17. Schematic of the temperature control circuit.

seen in Fig.16. Taking the temperature stabilization to times shorter than 0.001s may affect

Based on the above, it may be concluded that the outlined strategies have good performance and therefore they can be used to continue with the system design. Now the temperature control circuit will be included for the next analyses. So, after the validation of the models, the following step is to prove a temperature control, either using blocks from the library of the platform or using integration modules with other simulation programs. Once this is made, the designer can have a global idea about the performance of the system with the

For the temperature control of the micro hotplate obtained with post-processes compatible with CMOS technology described before, an analogue proportional control based on an operational amplifier (OP-AMP) is proposed (Barrentino et al., 2004). The schematic used for

**4. Model of the temperature control system and simulation** 

the basic simulations of the system is shown in Fig. 17.

Fig. 17. Schematic of the temperature control circuit.

the simulation dynamics at the beginning.

possibility to optimize the design.

From the mentioned electro-thermal characteristics of polysilicon, MH and TS can be used as part of a temperature control circuit. When the temperature sensor is connected as a feedback from the output to the input of the operational amplifier, the voltage drop in it is compared with a reference voltage (VREF). If the voltage difference between both inputs is zero, a power MOS transistor used as a switch connected at the output of the OP-AMP, will interrupt the biasing to the MH. A regulated current source applying pulses with constant amplitude should be used to bias TS (note that this does not cause Joule Effect to TS). Due to heat transmission by conduction, the resistance in the TS will increase (see Fig. 4) and its voltage drop will be modified. This property makes feasible to apply the TS in a temperature control circuit. On the other hand, the MH should be biased with a pulsed voltage source to heat the MHP. The magnitude of the pulse depends on the desired temperature. Besides, a pulsed stimulus will help to reduce the power consumption of the system.

Hence, when the reference voltage is higher than the feedback voltage drop, the output of the OP-AMP is high making the MOS transistor to conduct allowing a current flow through the MH. So, due to Joule Effect the temperature will gradually rise from room temperature up to a temperature established by the applied voltage, VREF. At the same time the TS will follow the MH's temperature, changing also the voltage drop on it. At the moment when both voltages are equal, the OP-AMP output will go low turning off the MOS transistor and heating will stop. This cycle will be repeated according to the heating and cooling dynamic behaviour of the microheater, operating around the desired temperature. It is expected that the MH will keep this temperature until the set point is modified if a different temperature target within the physical and mechanical limits of the MHP is desired. Next, to ease the dynamic analysis of the controlling circuit, an equivalent circuit was used to consider the thermal coupling between the heater and the sensor as well as the time constant of the electro-thermal system. The MH was replaced with a voltage controlled voltage source and the TS with an RC circuit, where R and C are the thermal resistance and thermal capacitance of polysilicon, respectively. The corresponding time constant of this circuit should be consistent with the thermal behaviour of the 3D model of the MH considered.

From the finite element analysis made above, different time constants resulted for heating and for cooling (see Figs. 6 and 7). That is the reason why the equivalent circuit shown in Fig. 18 has two resistance elements corresponding to each value. The magnitude of R for heating is 770 Ω and 950 Ω for cooling. The capacitance value used to fit with each time constants is 1μF; the circuit was biased with VDD=1.5V. The transient behaviour of the system can be evaluated using the circuit shown in Fig. 18. The voltage drop in C is equivalent to the temperature reached by the system and is also the feedback voltage of the OP-AMP. Finally, the MOS transistor switch is replaced by two MOS transistors (PMOS and NMOS) to handle charge and discharge (heating and cooling, respectively) of the capacitor, as the output of the OP-AMP goes high or low. This is the reason why two resistors are used, so heating and cooling can be considered. The results of the transient analysis of the circuit show that there will be control over the current flow of the microheater until the voltage drop on TS reaches the reference voltage (1V for the case illustrated in Fig. 19). As soon as the voltage at the gate of the MOS transistor, VGS, equals the reference voltage, then biasing of the microheater stops. According to the cooling dynamics of the microheater, bias will be switching on and off to hold the temperature around the operating point as can be seen in Fig. 19.

Performance Evaluation of a Temperature Control Stage Used on a

**5. Simulation using Simulink®**

Semiconductor Gas Sensor 3D Electro-Thermal Model Through Simulink® 183

Fig. 20. Output of the OP-AMP when a temperature sweep is applied to MH and TS.

used in the electro-thermal system together with the temperature control circuit.

be used in Simulink® to read more than one input or output if this is the case.

It is important to remark that simulating the control circuit with SPICE presents a limitation since it is not possible to directly model the electro-thermal coupling between the microheater and the temperature sensor. Therefore, the simulation result shown in Fig. 19 does not take into account the heat transfer by conduction. As a consequence, the result obtained in that simulation is not fully consistent with the real behaviour of the system. Then to consider the thermal coupling, other tools are required such that an analysis of the system with conditions near to real can be done as a function of the elements and materials

To simulate the dynamic performance of the circuit upon temperature changes over the MHP, it is necessary to follow a strategy that can take the SPICE model to a dynamic environment including random temperature changes in time. This is not a task that can be done exclusively with SPICE. Simulink® can help to achieve this using the SLPS integrating block and handling the electro-thermal feature of the micro hotplate with function blocks. It should be mentioned that the SLPS block allows interaction between models of the circuits described in SPICE and models found in Simulink®, whether they were designed using this platform or imported from some other environment. To make use of the SLPS block, it should be initialized from Matlab® as recommended by the supplier; once this is done functional blocks will be available in the library of Simulink. There, the desired interacting block can be loaded identifying the name of the circuit described in SPICE and the respective inputs and outputs are then selected. On one side, inputs are previously defined in SPICE as input voltage sources of the circuit, and on the other side, the output ports will drive a component or device simulated with SPICE. Nevertheless, the SLPS block has only one input and one output, both multiplexed since the use of MUX or DEMUX blocks must

Fig. 18. Equivalent circuit of the temperature control circuit.

In order to have a complete analysis of the performance of the circuit, it was simulated also in SPICE with a temperature sweep (27°C to 300°C) applied to MH and TS and keeping all other circuit components at room temperature (27°C). The result is shown in Fig. 20, where it can be seen that the output of the OP-AMP goes low when the set point (1V) is reached. This is the expected behaviour for the control circuit but it should be noted that this simulation is independent to time, so it is not a dynamic analysis. Then, from Figs. 19 and 20 it can be considered, as a first approximation, that the system behaves as it should even for heating or cooling of MH.

Fig. 19. Simulation results from the equivalent circuit shown in fig. 18 using SPICE.

Fig. 20. Output of the OP-AMP when a temperature sweep is applied to MH and TS.

## **5. Simulation using Simulink®**

182 Technology and Engineering Applications of Simulink

In order to have a complete analysis of the performance of the circuit, it was simulated also in SPICE with a temperature sweep (27°C to 300°C) applied to MH and TS and keeping all other circuit components at room temperature (27°C). The result is shown in Fig. 20, where it can be seen that the output of the OP-AMP goes low when the set point (1V) is reached. This is the expected behaviour for the control circuit but it should be noted that this simulation is independent to time, so it is not a dynamic analysis. Then, from Figs. 19 and 20 it can be considered, as a first approximation, that the system behaves as it should even for

Fig. 19. Simulation results from the equivalent circuit shown in fig. 18 using SPICE.

Fig. 18. Equivalent circuit of the temperature control circuit.

heating or cooling of MH.

It is important to remark that simulating the control circuit with SPICE presents a limitation since it is not possible to directly model the electro-thermal coupling between the microheater and the temperature sensor. Therefore, the simulation result shown in Fig. 19 does not take into account the heat transfer by conduction. As a consequence, the result obtained in that simulation is not fully consistent with the real behaviour of the system. Then to consider the thermal coupling, other tools are required such that an analysis of the system with conditions near to real can be done as a function of the elements and materials used in the electro-thermal system together with the temperature control circuit.

To simulate the dynamic performance of the circuit upon temperature changes over the MHP, it is necessary to follow a strategy that can take the SPICE model to a dynamic environment including random temperature changes in time. This is not a task that can be done exclusively with SPICE. Simulink® can help to achieve this using the SLPS integrating block and handling the electro-thermal feature of the micro hotplate with function blocks. It should be mentioned that the SLPS block allows interaction between models of the circuits described in SPICE and models found in Simulink®, whether they were designed using this platform or imported from some other environment. To make use of the SLPS block, it should be initialized from Matlab® as recommended by the supplier; once this is done functional blocks will be available in the library of Simulink. There, the desired interacting block can be loaded identifying the name of the circuit described in SPICE and the respective inputs and outputs are then selected. On one side, inputs are previously defined in SPICE as input voltage sources of the circuit, and on the other side, the output ports will drive a component or device simulated with SPICE. Nevertheless, the SLPS block has only one input and one output, both multiplexed since the use of MUX or DEMUX blocks must be used in Simulink® to read more than one input or output if this is the case.

Performance Evaluation of a Temperature Control Stage Used on a

operating temperature of MHP will be considered.

Fig. 22. On-Off Configuration.

temperature variations.

Semiconductor Gas Sensor 3D Electro-Thermal Model Through Simulink® 185

The block diagram is shown in Fig. 24. Random temperature variation affecting the system can be included in the descriptions shown in Figs. 22 and 24 with the objective to observe the response of the control circuit of the gas sensor system when temperature perturbations are added to the surrounding environment. Hence, any external factor that should alter the

When the diagram shown in Fig.22 is simulated, a plot of the temperature variation across time can be obtained showing the behaviour of the MHP connected to a temperature control. As soon as the MH reaches a given temperature established with the set point, the control system maintains this temperature as close as possible, through the control of the switch transistor with the output of the OP-AMP. Fig. 23 shows that a rather good temperature control is achieved around 250°C, as the proposed operating temperature.

Fig. 23. Simulated operating temperature of MHP without SLPS and considering random

Based on this, the circuit described with SPICE must be modified to work properly with the SLPS block. The proposal here reported gives good results if the temperature sensor resistor and its biasing current supply are eliminated, keeping only the OP-AMP, the MOS transistor switch and the microheater resistor, connecting them to voltage sources in SPICE, as is shown in Fig. 21. With this procedure, the SLPS block uses the simulation profile of SPICE and the output files from the corresponding simulation. The power supplies and the MHP model are external to this block in the Simulink® environment. The simulation will contain every element of the micro hotplate and the temperature control circuit must be adjusted to the previous simulation conditions. The values used in the voltage sources are shown in Table 1. They correspond to requirements like low power consumption, VDD=1.5V, a compensation voltage for the OP-AMP of Vnb=0.6V, a feeding voltage for the MH of Vtrans=3V and a reference voltage corresponding to 250°C, VREF=0.971V.

Fig. 21. SPICE circuit description for simulation with Simulink®.

## **6. Results and discussion**

First, a simulation was done for an on-off control circuit using the function blocks that consider the electro-thermal model including TS and MH with their respective features like thermal coupling, temperature sensor morphology and its biasing current. This simulation helps to prove the models developed for the MHP and the dynamic performance. Fig. 22 shows the configuration used for this purpose. On the other hand, it is possible to use the SLPS block with the simulation profile together with the output files from SPICE to integrate the results from both simulations within only one process. Using the conditions specified in Table 1, the performance of the proposed control circuit can be evaluated.


Table 1. Variable values for simulation.

The block diagram is shown in Fig. 24. Random temperature variation affecting the system can be included in the descriptions shown in Figs. 22 and 24 with the objective to observe the response of the control circuit of the gas sensor system when temperature perturbations are added to the surrounding environment. Hence, any external factor that should alter the operating temperature of MHP will be considered.

Fig. 22. On-Off Configuration.

184 Technology and Engineering Applications of Simulink

Based on this, the circuit described with SPICE must be modified to work properly with the SLPS block. The proposal here reported gives good results if the temperature sensor resistor and its biasing current supply are eliminated, keeping only the OP-AMP, the MOS transistor switch and the microheater resistor, connecting them to voltage sources in SPICE, as is shown in Fig. 21. With this procedure, the SLPS block uses the simulation profile of SPICE and the output files from the corresponding simulation. The power supplies and the MHP model are external to this block in the Simulink® environment. The simulation will contain every element of the micro hotplate and the temperature control circuit must be adjusted to the previous simulation conditions. The values used in the voltage sources are shown in Table 1. They correspond to requirements like low power consumption, VDD=1.5V, a compensation voltage for the OP-AMP of Vnb=0.6V, a feeding voltage for the MH of

Vtrans=3V and a reference voltage corresponding to 250°C, VREF=0.971V.

Fig. 21. SPICE circuit description for simulation with Simulink®.

Table 1, the performance of the proposed control circuit can be evaluated.

First, a simulation was done for an on-off control circuit using the function blocks that consider the electro-thermal model including TS and MH with their respective features like thermal coupling, temperature sensor morphology and its biasing current. This simulation helps to prove the models developed for the MHP and the dynamic performance. Fig. 22 shows the configuration used for this purpose. On the other hand, it is possible to use the SLPS block with the simulation profile together with the output files from SPICE to integrate the results from both simulations within only one process. Using the conditions specified in

> Variable Values Vdd 1.5V Vnd 0.6V Vrans 3V VREF 0.971V

**6. Results and discussion** 

Table 1. Variable values for simulation.

When the diagram shown in Fig.22 is simulated, a plot of the temperature variation across time can be obtained showing the behaviour of the MHP connected to a temperature control. As soon as the MH reaches a given temperature established with the set point, the control system maintains this temperature as close as possible, through the control of the switch transistor with the output of the OP-AMP. Fig. 23 shows that a rather good temperature control is achieved around 250°C, as the proposed operating temperature.

Fig. 23. Simulated operating temperature of MHP without SLPS and considering random temperature variations.

Performance Evaluation of a Temperature Control Stage Used on a

Semiconductor Gas Sensor 3D Electro-Thermal Model Through Simulink® 187

either using a simple on-off control block or with a more elaborate control circuit based on an OP-AMP. If a rapid evaluation of the performance is desired, a simple block with Simulink® can be used, which can give an approximated idea if the system is giving the expected results. This can be the basis for a more complete analysis using a more complex temperature control circuit based in the SLPS integrating block of Simulink®. It should be remembered that a single geometry for the MH and TS was studied in this work, but depending on the results obtained in the multi-physics environment, it should be optimized also to avoid heat concentration on borders and corners of MH that will damage the structure and to minimize as well the voltage applied in order to reach a given operating temperature. Besides, as this can be a time consuming process until a highly optimized design is obtained, it is convenient to have a strategy that helps to readily have the system design that can be fabricated with no fatal consequences. Finally, it was possible to couple different software environments to complete in short time the electro-thermal analysis of a 3D structure with its temperature control circuit included and with no huge computational resources. For example, the simulation performed to obtain the results shown in Fig. 12 took approximately from 4 to 6 hours to finish using a laptop loaded with a 3D solids engineering and analysis software based on finite elements, to generate the model. On the other hand, the simulation of the exported static model shown in Fig. 9, took from 1 to 2 hours to give the same results in the same machine. Finally, the approximated model illustrated in Fig. 11, gave also the same results in as much as 15 minutes. This is a dramatic difference and the criteria used in the present

analysis can be used in turn, to extrapolate to other kind of systems or analyses.

The simulation with approximated models have excellent results in repetitive tests, since they reduce the use of hardware requirements, optimize simulation time and avoid interconnection as well as related problems among software programs. There may be cases in which models created in a multi-physics environment can be exported to Simulink® in a dynamic fashion, but the risk of no convergence may be present if there are inconvenient or unexpected changes in the simulation time step. In this case, the exported model will not be useful and the evaluation of the control stage cannot be done. Regardless the approximated model method is faster and allows driving signals as ramps or sine waveforms, results will depend on the equation accuracy and its range of validity. Therefore, the simulation option using approximated models can be the solution to accelerate tests and to apply controls to the model. Also, it can be considered a useful methodology for those users working with rather complex models and having no high capacity resources. However, it is recommended to always verify, if possible, the results with the original program and to try to find models using higher order TFs. There is no doubt that the CAD resources definitely support developing and design activities. Nevertheless, sometimes they may be limited to specific tasks so it is necessary to assess and to explore different alternatives. For some cases, partial results thrown by simulations can be considered good enough, but it will be preferable to have wide information about the behaviour of a system. Coupling compatible simulation environments can be a reliable strategy as was demonstrated for the case presented in this work, where satisfactory results were obtained at different stages. Good results were obtained for the operation of a temperature control circuit applied to a gas sensor system

**7. Conclusions** 

using two different methods.

Fig. 24. Block diagram using Simulink®.

Later, including the SLPS block as indicated in Fig. 24, the same simulation as before was made and the results are similar to those presented in Fig. 23, oscillating too around 250°C, as can be seen in Fig. 25.

Fig. 25. Simulated operating temperature of MHP with the SLPS block included, considering random temperature variations.

Obviously the oscillations shown in Figs. 23 and Fig. 25 are different but there are two reasons for this: first, since random temperature variations were used, every simulation will generate different values that are not coincident each other; second, it should be remembered that the circuit simulated using SLPS, considered a pulsed current driver to decrease power consumption. Although it can be considered that real external temperature variations will not have such abrupt and rapid transitions, the temperature control circuit is able to hold the temperature around the desired set point. After these simulations, it can be said that the objectives of this study were achieved, this is, a proper design of a temperature control for the MHP that considers the electro-thermal coupling of the materials used for the SGS and with a particular geometry. As it was shown, both strategies can be used for evaluation and validation since they give reasonable results either using a simple on-off control block or with a more elaborate control circuit based on an OP-AMP. If a rapid evaluation of the performance is desired, a simple block with Simulink® can be used, which can give an approximated idea if the system is giving the expected results. This can be the basis for a more complete analysis using a more complex temperature control circuit based in the SLPS integrating block of Simulink®. It should be remembered that a single geometry for the MH and TS was studied in this work, but depending on the results obtained in the multi-physics environment, it should be optimized also to avoid heat concentration on borders and corners of MH that will damage the structure and to minimize as well the voltage applied in order to reach a given operating temperature. Besides, as this can be a time consuming process until a highly optimized design is obtained, it is convenient to have a strategy that helps to readily have the system design that can be fabricated with no fatal consequences. Finally, it was possible to couple different software environments to complete in short time the electro-thermal analysis of a 3D structure with its temperature control circuit included and with no huge computational resources. For example, the simulation performed to obtain the results shown in Fig. 12 took approximately from 4 to 6 hours to finish using a laptop loaded with a 3D solids engineering and analysis software based on finite elements, to generate the model. On the other hand, the simulation of the exported static model shown in Fig. 9, took from 1 to 2 hours to give the same results in the same machine. Finally, the approximated model illustrated in Fig. 11, gave also the same results in as much as 15 minutes. This is a dramatic difference and the criteria used in the present analysis can be used in turn, to extrapolate to other kind of systems or analyses.

## **7. Conclusions**

186 Technology and Engineering Applications of Simulink

Later, including the SLPS block as indicated in Fig. 24, the same simulation as before was made and the results are similar to those presented in Fig. 23, oscillating too around 250°C,

Fig. 25. Simulated operating temperature of MHP with the SLPS block included, considering

Obviously the oscillations shown in Figs. 23 and Fig. 25 are different but there are two reasons for this: first, since random temperature variations were used, every simulation will generate different values that are not coincident each other; second, it should be remembered that the circuit simulated using SLPS, considered a pulsed current driver to decrease power consumption. Although it can be considered that real external temperature variations will not have such abrupt and rapid transitions, the temperature control circuit is able to hold the temperature around the desired set point. After these simulations, it can be said that the objectives of this study were achieved, this is, a proper design of a temperature control for the MHP that considers the electro-thermal coupling of the materials used for the SGS and with a particular geometry. As it was shown, both strategies can be used for evaluation and validation since they give reasonable results

Fig. 24. Block diagram using Simulink®.

as can be seen in Fig. 25.

random temperature variations.

The simulation with approximated models have excellent results in repetitive tests, since they reduce the use of hardware requirements, optimize simulation time and avoid interconnection as well as related problems among software programs. There may be cases in which models created in a multi-physics environment can be exported to Simulink® in a dynamic fashion, but the risk of no convergence may be present if there are inconvenient or unexpected changes in the simulation time step. In this case, the exported model will not be useful and the evaluation of the control stage cannot be done. Regardless the approximated model method is faster and allows driving signals as ramps or sine waveforms, results will depend on the equation accuracy and its range of validity. Therefore, the simulation option using approximated models can be the solution to accelerate tests and to apply controls to the model. Also, it can be considered a useful methodology for those users working with rather complex models and having no high capacity resources. However, it is recommended to always verify, if possible, the results with the original program and to try to find models using higher order TFs. There is no doubt that the CAD resources definitely support developing and design activities. Nevertheless, sometimes they may be limited to specific tasks so it is necessary to assess and to explore different alternatives. For some cases, partial results thrown by simulations can be considered good enough, but it will be preferable to have wide information about the behaviour of a system. Coupling compatible simulation environments can be a reliable strategy as was demonstrated for the case presented in this work, where satisfactory results were obtained at different stages. Good results were obtained for the operation of a temperature control circuit applied to a gas sensor system using two different methods.

**1. Introduction**

(WSP) systems.

chapter.

Grazyna Barna ˙

*Poland*

**9**

*Rail Vehicles Institute "TABOR"*

**Matlab Simulink® Model of** 

When a braking force applied to the axle sets of a rail vehicle exceeds a critical value, which depends on the wheel-rail adherence, the wheels start to slide. If no corrective action is taken, in a very short time the wheels can be locked. Locking of the wheels or their excessive slide can result in increased braking distance and damage to the wheel rims (flat spots, also called "flats"). Wheel flats are sources of vibration and noise, they lower riding quality of a vehicle as well as passenger comfort, but first and foremost increasing of the braking distance directly impairs safety of the train staff and passengers and also of people nearby. In order to prevent excessive wheel slide and wheel lock, rail vehicles are equipped with Wheel Slide Protection

**a Braked Rail Vehicle and Its Applications** 

From the point of view of a WSP controller, the rail vehicle is a strongly non-linear and non-stationary plant. From the other hand there exist intuitive expert knowledge concerning the way the slide should be controlled. For these two reasons fuzzy logic controllers (FLC) are

One of the basic problems concerning FLCs is lack of formal methods of design. Designing the rule bases of the fuzzy controllers is usually performed using a trial-and-error method, which in turn requires performing numerous tests of the controlled plant. When a controlled plant is a rail vehicle, possibilities of performing such tests are very limited due to high costs of tests and danger of damaging the wheels. A commonly used solution to this problem is employing a simulator of a braked rail vehicle, which can be used for preliminary design, optimization and tuning of the controllers. Tests on a real object are performed at the last stage of the design

The purpose of this chapter is to present the Matlab Simulink model of a braked rail vehicle and its many applications, which comprise both designing and testing of the WSP systems. In the next two introductory sections the slide phenomenon as well as the structure and principle of operation of WSP systems are described in order to provide the readers who are not familiar with the WSP devices with some useful information which would facilitate reading of this

During braking of a rail vehicle, a braking torque generated by a rail vehicle braking system is applied to the axle sets. When a value of this torque exceeds a boundary value, which

widely used in WSP systems (Caldara et al. (1996), Barna (2009)).

process in order to verify the designed controller.

**2. Slide during braking of a rail vehicle**

## **8. Acknowledgements**

This work was financed by CONACYT through the project number 57429. The authors thanks to Oliverio Arellano Cárdenas for his help in figures edition and formatting.

## **9. References**


## **Matlab Simulink® Model of a Braked Rail Vehicle and Its Applications**

Grazyna Barna ˙ *Rail Vehicles Institute "TABOR" Poland*

#### **1. Introduction**

188 Technology and Engineering Applications of Simulink

This work was financed by CONACYT through the project number 57429. The authors

Afridi, M. Y.; Suehle, J. S.; Zaghloul, M. E.; Berning, D. W.; Hefner, A. R.; Cavicchi, R. E.;

Barretino, D.; Graf, M.; Song, W. H.; Kirstein, K. U.; Hierlemann A. & Baltes, H. (2004).

Cerdà, J.; Manzano, J.; Arbiol, J.; Cirera, A.; Puigcorbé, J.; Vilà, A.; Sabaté, N.; Gràcia, I.;

Fedder, G. K. (2005). CMOS-Based sensors. *Proceedings of the 4th IEEE Conference on Sensors* 

Goering, R. (2004). Matlab® Edges Closer to Electronic Design Automation World. In: *EE* 

Hsu, T. R., (2002), *MEMS & Microsystems* 1st Ed., Mc Graw Hill, ISBN 0-07-239391-2, New

Muñoz, J. S., Valencia, R. and Nieto, C. (2008). COMSOL and MATLAB Integration to

Ogata, K., (1999), *Problemas de Ingeniería de control utilizando Matlab® 1ª Ed.*, Prentice Hall

Optimize Heat Exchangers Using Genetic Algorithms Technique, *Proceedings of* 

http://www.eetimes.com/electronics-news/4050334/Matlab-edges-closer-to-

Semanick, S.; Montgomery, C. B. & Taylor, C. J. (2002). A monolithic CMOS microhotplate-based gas sensor system. *IEEE Sensors Journal,Vol.*2, pp. 644-655,

Hotplate-based monolithic CMOS microsystem for gas detection and material characterization for operating temperatures up to 500ºC. *IEEE Journal of Solid-State* 

Cané, C. & Morante, J. R. (2006). Micromachined twin gas sensor for CO and O2 quantification based on catalytically modified nano-SnO2, *Sens. and Act. B: Chem.*,

*Oct. 30-Nov. 3 2005, IEEE*, pp. 125-128, ISBN 0780390563, Irvine, CA, U.S.A. Oct. 30

thanks to Oliverio Arellano Cárdenas for his help in figures edition and formatting.

*Ciruits*, Vol.39, pp. 1202-1207, ISSN 0018-9200.

*Times.* Last Accessed: Feb. 8, 2011. Available from:

*COMSOL Conference 2008*, Boston, U.S.A, 2008.

Iberia, ISBN 84-8322-046-6, Madrid, España, 1999.

Vol.114, pp. 881-892, ISSN 09254005.

electronic-design-automation-world

York, NY, U.S.A., 2002.

**8. Acknowledgements** 

ISSN 1530437X.

– Nov. 3, 2005.

**9. References** 

When a braking force applied to the axle sets of a rail vehicle exceeds a critical value, which depends on the wheel-rail adherence, the wheels start to slide. If no corrective action is taken, in a very short time the wheels can be locked. Locking of the wheels or their excessive slide can result in increased braking distance and damage to the wheel rims (flat spots, also called "flats"). Wheel flats are sources of vibration and noise, they lower riding quality of a vehicle as well as passenger comfort, but first and foremost increasing of the braking distance directly impairs safety of the train staff and passengers and also of people nearby. In order to prevent excessive wheel slide and wheel lock, rail vehicles are equipped with Wheel Slide Protection (WSP) systems.

From the point of view of a WSP controller, the rail vehicle is a strongly non-linear and non-stationary plant. From the other hand there exist intuitive expert knowledge concerning the way the slide should be controlled. For these two reasons fuzzy logic controllers (FLC) are widely used in WSP systems (Caldara et al. (1996), Barna (2009)).

One of the basic problems concerning FLCs is lack of formal methods of design. Designing the rule bases of the fuzzy controllers is usually performed using a trial-and-error method, which in turn requires performing numerous tests of the controlled plant. When a controlled plant is a rail vehicle, possibilities of performing such tests are very limited due to high costs of tests and danger of damaging the wheels. A commonly used solution to this problem is employing a simulator of a braked rail vehicle, which can be used for preliminary design, optimization and tuning of the controllers. Tests on a real object are performed at the last stage of the design process in order to verify the designed controller.

The purpose of this chapter is to present the Matlab Simulink model of a braked rail vehicle and its many applications, which comprise both designing and testing of the WSP systems. In the next two introductory sections the slide phenomenon as well as the structure and principle of operation of WSP systems are described in order to provide the readers who are not familiar with the WSP devices with some useful information which would facilitate reading of this chapter.

#### **2. Slide during braking of a rail vehicle**

During braking of a rail vehicle, a braking torque generated by a rail vehicle braking system is applied to the axle sets. When a value of this torque exceeds a boundary value, which

subjected to cracks and crumbling. It develops first of all at the surface, and then propagating far into the wheel material. It results in developing the loss of the material of the wheel

Matlab Simulink® Model of a Braked Rail Vehicle and Its Applications 191

Damage to the wheel treads may be caused not only by wheel lock, but also by excessive difference (above 30 km/h) between vehicle velocity and circumferential velocity of a wheel. Vibrations also lower the comfort of the passengers and the trains staff, negatively influencing their mood and ability to work after the journey, and also their health. In Fig. 3 simulation characteristics are shown for both a vehicle with round wheels and a vehicle with one flat spot on one wheel, at speed 140 km/h. The influence of a flat spot is seen the most clearly for frequency of app. 13,5 Hz. Lowering of comfort, defined as increasing of the effective value of vertical acceleration in the middle of the body amounts in this particular case to app. 30 %,

The concern about wheel slide has increased, as rail vehicles began running at higher speeds, because not only braking torque which is necessary to brake a vehicle needs to be bigger, which increases the probability of the wheel slide, but also consequences of such event are

To prevent the described above situations, rail vehicles are equipped with Wheel Slide Protection (WSP) devices, which detect slide of the vehicle wheels and adequately control the braking torque, not only preventing wheels from being locked, but also increasing the

Because a safety critical system is a system the failure of which can result in severe consequences, e.g. death or injury of people, or significant damages to property or environment, thus it is evident that a WSP device is a safety critical system (Barna & Kaluba (2009)). Standard CEN (2009) and leaflet UIC (2005) contain specification of requirements for

adhesion coefficient value, thus making the braking distance as short as possible.

Fig. 3. Swing and frequency characteristics of vertical acceleration in the middle of the passenger car body, respectively without and with a wheel flat (Ofierzy ´nski (2008))

A WSP system consists of the following elements: wheel speed sensors, controller and dump valves. The block diagram of a WSP system for a two-axle bogie with block brake is shown in

A WSP controller on the basis on signals from speed sensors of all wheel-sets, performs calculation of circumferential wheel velocities and determines circumferential wheel

treads (Jergéus (1998); Kwa´snikowski & Firlik (2006)).

and it depends on the vehicle velocity (Ofierzy ´nski (2008)).

more serious when it occurs at high vehicle velocities.

both structure and functions of WSP devices.

**3. Wheel Slide Protection systems**

Fig. 4 ( Barna (2009), UIC (2005), CEN (2009)).

depends on the instantaneous value of wheel-rail adherence, than the circumferential speed of the wheels starts to decrease (which is called sliding). If no corrective action is taken, in a very short time the wheels can be locked.

Locking of the wheels of a rail vehicle, as well as an excessive wheel slide have several negative consequences, out of which two are critical. First and foremost, due to decreasing of the adhesion coefficient, the braking force remains constant at a small value, which makes impossible effective braking of a vehicle and results in significant increasing of the braking distance. As an example, the braking distance of a 150A passenger car braked at 120 km/h amounts, according to both field and simulation tests approximately 480 m, while at braking with all the wheels locked at the beginning of braking it can be increased even up to 800 m. Increasing of the braking distance impairs the safety of passengers, train staff and people in proximity. Secondly, when locked wheels slide along the rails, flat spots (called "flats") can be produced on wheel treads. The depth of a wheel flat can amount up to several millimeters (Pawełczyk (2008)), especially in case of a prolonged slide. In Fig. 1 and Fig. 2 photographs of wheel flats of a passenger car wheels are shown.

Fig. 1. Flat spot on a wheel tread 1

Fig. 2. Flat spot on a wheel tread 2

Increase of the temperature of the sliding wheel and rapid cooling of the outer wheel layer, caused by regaining speed after the wheel slide has ceased, may result in forming of martensite around the flat spot. Martensite is a form of steel, firm, but fragile, which is 2 Will-be-set-by-IN-TECH

depends on the instantaneous value of wheel-rail adherence, than the circumferential speed of the wheels starts to decrease (which is called sliding). If no corrective action is taken, in a

Locking of the wheels of a rail vehicle, as well as an excessive wheel slide have several negative consequences, out of which two are critical. First and foremost, due to decreasing of the adhesion coefficient, the braking force remains constant at a small value, which makes impossible effective braking of a vehicle and results in significant increasing of the braking distance. As an example, the braking distance of a 150A passenger car braked at 120 km/h amounts, according to both field and simulation tests approximately 480 m, while at braking with all the wheels locked at the beginning of braking it can be increased even up to 800 m. Increasing of the braking distance impairs the safety of passengers, train staff and people in proximity. Secondly, when locked wheels slide along the rails, flat spots (called "flats") can be produced on wheel treads. The depth of a wheel flat can amount up to several millimeters (Pawełczyk (2008)), especially in case of a prolonged slide. In Fig. 1 and Fig. 2

Increase of the temperature of the sliding wheel and rapid cooling of the outer wheel layer, caused by regaining speed after the wheel slide has ceased, may result in forming of martensite around the flat spot. Martensite is a form of steel, firm, but fragile, which is

very short time the wheels can be locked.

Fig. 1. Flat spot on a wheel tread 1

Fig. 2. Flat spot on a wheel tread 2

photographs of wheel flats of a passenger car wheels are shown.

subjected to cracks and crumbling. It develops first of all at the surface, and then propagating far into the wheel material. It results in developing the loss of the material of the wheel treads (Jergéus (1998); Kwa´snikowski & Firlik (2006)).

Damage to the wheel treads may be caused not only by wheel lock, but also by excessive difference (above 30 km/h) between vehicle velocity and circumferential velocity of a wheel.

Vibrations also lower the comfort of the passengers and the trains staff, negatively influencing their mood and ability to work after the journey, and also their health. In Fig. 3 simulation characteristics are shown for both a vehicle with round wheels and a vehicle with one flat spot on one wheel, at speed 140 km/h. The influence of a flat spot is seen the most clearly for frequency of app. 13,5 Hz. Lowering of comfort, defined as increasing of the effective value of vertical acceleration in the middle of the body amounts in this particular case to app. 30 %, and it depends on the vehicle velocity (Ofierzy ´nski (2008)).

The concern about wheel slide has increased, as rail vehicles began running at higher speeds, because not only braking torque which is necessary to brake a vehicle needs to be bigger, which increases the probability of the wheel slide, but also consequences of such event are more serious when it occurs at high vehicle velocities.

To prevent the described above situations, rail vehicles are equipped with Wheel Slide Protection (WSP) devices, which detect slide of the vehicle wheels and adequately control the braking torque, not only preventing wheels from being locked, but also increasing the adhesion coefficient value, thus making the braking distance as short as possible.

Because a safety critical system is a system the failure of which can result in severe consequences, e.g. death or injury of people, or significant damages to property or environment, thus it is evident that a WSP device is a safety critical system (Barna & Kaluba (2009)). Standard CEN (2009) and leaflet UIC (2005) contain specification of requirements for both structure and functions of WSP devices.

Fig. 3. Swing and frequency characteristics of vertical acceleration in the middle of the passenger car body, respectively without and with a wheel flat (Ofierzy ´nski (2008))

### **3. Wheel Slide Protection systems**

A WSP system consists of the following elements: wheel speed sensors, controller and dump valves. The block diagram of a WSP system for a two-axle bogie with block brake is shown in Fig. 4 ( Barna (2009), UIC (2005), CEN (2009)).

A WSP controller on the basis on signals from speed sensors of all wheel-sets, performs calculation of circumferential wheel velocities and determines circumferential wheel

**4. Characteristics of the WSP controllers**

• controlled plant is highly non-linear

**–** dump valves behavior is non-linear

(2009)).

most crucial:

measured

the brake cylinders

Wheel Slide Protection (WSP) devices control braking torque, not only preventing wheels from being locked, but also increasing the adhesion coefficient value, making the braking distance as short as possible. Thus the main task of WSP systems is "*to make the best use of available adhesion for all intended-operating conditions by a controlled reduction and restoration of the brake force to prevent axle sets from locking and uncontrolled sliding due to low adhesion*" (CEN

Matlab Simulink® Model of a Braked Rail Vehicle and Its Applications 193

The main difficulties concerning the control algorithm resulting from the characteristics of a braked rail vehicle as a controlled plant are listed below, the first two out of which are the

• instantaneous value of the coefficient of adhesion between wheel and rail *ψ* cannot be

• controlled plant is non-stationary — parameters of characteristics *ψ* = *f*(*s*) change in a

• direct control of a braking torque is not possible; actuators control the pressure values in

The task of a WSP controller is to achieve the highest possible at the current adhesion conditions adhesion forces (forces transferred to the wheels of the axle sets), which is equivalent with achieving the shortest possible braking distance. This goal is achieved, as

In Fig. 5 a trajectory *ψ* = *f*(*s*) for an optimally working controller has been shown, superimposed on a generalized adhesion versus slide curve. A WSP device equipped with such a controller, called in Boiteux (1999) a Second Generation WSP, realizes control in the neighbourhood of the point B of an adhesion curve and maintains a value of a relative slide *s*

A control system such as described above is called an *extremum regulation system*. The task of the extremum regulation system is maintaining a controlled signal as close as possible to the extremum value, which changes depending on disturbance signals (Kaczorek (1977)). In the case of a braked rail vehicle and a WSP controller this extremum value should be the

In Fig. 6 an overall block diagram of a WSP control system is shown, which is a base for

A braked rail vehicle as a control plant is strongly non-linear and non-stationary, therefore the WSP controllers are usually designed using Fuzzy Logic and Expert

• controlling of the brake cylinder pressures is performed indirectly and discretely

mentioned in Section 3 by proper control of the dump valves of each axle set.

instantaneous value of adhesion coefficient between wheel and rail *ψ*.

practically realized control algorithms of WSP controllers.

• translation vehicle speed *vT* is usually not measured due to costs

**–** characteristics *Mb* = *f*(*pc*) is characterized by hysteresis

• instantaneous values of the cylinder pressures are not known.

close to its optimum value *sB* (see subsection 5.2.4).

way which is possible to comprehend in stochastic manner only

**–** characteristics *ψ* = *f*(*s*) is non-linear (see Fig. 13)

Fig. 4. Block diagram of a WSP system: DV1, DV2 - dump valves, V1, V2 - speed measurement signals, EV1, EV2, BV1, BV2 - dump valve control signals

accelerations and wheel slides, performs estimation of vehicle velocity (called the reference velocity), and then, on the basis of the obtained values appropriately controls the dump valves, in order to adjust the braking torque to the instantaneous adhesion conditions (Barna (2010a)). When a slide is detected, the braking torque should be adequately decreased by proper control of the dump valves; after recovery of adhesion the torque should be increased in order to provide effective braking of the vehicle.

Speed sensors make possible measuring the angular speeds of the wheels – they generate square-wave signals from which circumferential wheel speeds can be calculated, and which are main input signals for the WSP controller.

The WSP actuators are dump valves, mounted as close as possible to the brake cylinders. The valves can adopt one of the three states, which makes possible increasing, decreasing or maintaining the pressure value in the brake cylinders:


Both filling and venting of the dump valves can be performed in one of the two following ways (Boiteux (1999)):


In Barna (2009) seven levels of pressure change have been proposed, designated P1, P2, P3, U1, U2, U3, H. The levels are obtained by cyclical applying sequences of the three mentioned above states of the valves **P3**: P, **P2**: P H, **P1**: P **H**: H, **U1**: U H H, **U2**: U H, **U3**: U. Each of the states is applied for a defined period of time, e.g. 100 ms.

## **4. Characteristics of the WSP controllers**

4 Will-be-set-by-IN-TECH

Fig. 4. Block diagram of a WSP system: DV1, DV2 - dump valves, V1, V2 - speed

accelerations and wheel slides, performs estimation of vehicle velocity (called the reference velocity), and then, on the basis of the obtained values appropriately controls the dump valves, in order to adjust the braking torque to the instantaneous adhesion conditions (Barna (2010a)). When a slide is detected, the braking torque should be adequately decreased by proper control of the dump valves; after recovery of adhesion the torque should be increased

Speed sensors make possible measuring the angular speeds of the wheels – they generate square-wave signals from which circumferential wheel speeds can be calculated, and which

The WSP actuators are dump valves, mounted as close as possible to the brake cylinders. The valves can adopt one of the three states, which makes possible increasing, decreasing or

Both filling and venting of the dump valves can be performed in one of the two following

• *graduated*: increasing or decreasing of the pressure are realized by applying repeated

In Barna (2009) seven levels of pressure change have been proposed, designated P1, P2, P3, U1, U2, U3, H. The levels are obtained by cyclical applying sequences of the three mentioned above states of the valves **P3**: P, **P2**: P H, **P1**: P **H**: H, **U1**: U H H, **U2**: U H, **U3**: U. Each of the

measurement signals, EV1, EV2, BV1, BV2 - dump valve control signals

in order to provide effective braking of the vehicle.

maintaining the pressure value in the brake cylinders:

• cutting-off the air supply without venting (H).

sequences of respectively F and H or V and H states.

states is applied for a defined period of time, e.g. 100 ms.

• *continuous*: applying either F or V state

are main input signals for the WSP controller.

• filling of the cylinders (F) • venting of the cylinders (V)

ways (Boiteux (1999)):

Wheel Slide Protection (WSP) devices control braking torque, not only preventing wheels from being locked, but also increasing the adhesion coefficient value, making the braking distance as short as possible. Thus the main task of WSP systems is "*to make the best use of available adhesion for all intended-operating conditions by a controlled reduction and restoration of the brake force to prevent axle sets from locking and uncontrolled sliding due to low adhesion*" (CEN (2009)).

The main difficulties concerning the control algorithm resulting from the characteristics of a braked rail vehicle as a controlled plant are listed below, the first two out of which are the most crucial:

	- **–** characteristics *ψ* = *f*(*s*) is non-linear (see Fig. 13)
	- **–** characteristics *Mb* = *f*(*pc*) is characterized by hysteresis
	- **–** dump valves behavior is non-linear

The task of a WSP controller is to achieve the highest possible at the current adhesion conditions adhesion forces (forces transferred to the wheels of the axle sets), which is equivalent with achieving the shortest possible braking distance. This goal is achieved, as mentioned in Section 3 by proper control of the dump valves of each axle set.

In Fig. 5 a trajectory *ψ* = *f*(*s*) for an optimally working controller has been shown, superimposed on a generalized adhesion versus slide curve. A WSP device equipped with such a controller, called in Boiteux (1999) a Second Generation WSP, realizes control in the neighbourhood of the point B of an adhesion curve and maintains a value of a relative slide *s* close to its optimum value *sB* (see subsection 5.2.4).

A control system such as described above is called an *extremum regulation system*. The task of the extremum regulation system is maintaining a controlled signal as close as possible to the extremum value, which changes depending on disturbance signals (Kaczorek (1977)). In the case of a braked rail vehicle and a WSP controller this extremum value should be the instantaneous value of adhesion coefficient between wheel and rail *ψ*.

In Fig. 6 an overall block diagram of a WSP control system is shown, which is a base for practically realized control algorithms of WSP controllers.

A braked rail vehicle as a control plant is strongly non-linear and non-stationary, therefore the WSP controllers are usually designed using Fuzzy Logic and Expert

obtained by feeding the calculated continuous wheel velocities (see section 5.2.3) through a model of speed sensors. Fuzzy Logic Controllers are of PI type. The Simulink subsystem determining the reference velocity *vref* is shown in Fig. 8. The output of each FLC is an integer in the range from −3 to 3, designating one of seven levels of pressure change (see section 3), which is then transformed in a post-processing block into signals controlling a dump valve.

Matlab Simulink® Model of a Braked Rail Vehicle and Its Applications 195

FLC 4

FLC 3

FLC 2

FLC 1

3.6

Trigger

0.01

v 1

v 2

vref 1

10.0

Trigger

Fig. 7. Fuzzy Logic WSP control system implemented in Matlab Simulink

Fig. 8. Matlab Simulink subsystem determining the reference velocity *vref*

max max

Valve

round 1

vref 1

Fig. 5. The principle of operation of an optimal Wheel Slide Protection (Boiteux (1999))

Fig. 6. Practical block diagram of WSP control system (Barna (2009))

Knowledge techniques (Caldara et al. (1996); Cheok & Shiomi (2000); Mauer (1995); Sanz & Pérez-Rodríguez (1997); Will & Zak (2000)). ˙

The WSP control system using Fuzzy Logic controllers implemented in Matlab Simulink is shown in Fig. 7. The control system is realized as Triggered Subsystem Simulink block, triggered every 100 ms, which corresponds to a real controller program cycle time. The controller blocks for every axle are implemented with Fuzzy Logic Controller Simulink blocks. The inputs for each FLC block are absolute axle slide *σ* [km/h] and axle acceleration *a*, which are calculated on the basis of the measured wheel circumferential velocities and vehicle reference velocity *vref* , as shown in Fig. 6. The measured wheel circumferential velocities are 6 Will-be-set-by-IN-TECH

B

Fig. 5. The principle of operation of an optimal Wheel Slide Protection (Boiteux (1999))

*V O*(*k*)*,*

*ψ*

Δ

Δ

*σ***ˆ**(*k*) **Δ***v*(*k*) **Δ<sup>2</sup>***v*(*k*)

> Determining the reference speed

*vref* (*k*)

Fig. 6. Practical block diagram of WSP control system (Barna (2009))

˙

+ *−*

*σ***ˆ**(*k*)

Pérez-Rodríguez (1997); Will & Zak (2000)).

*α*

A

*s*

1

Rail vehicle

*V L*(*k*) **<sup>Δ</sup>***pc*(*k*) *<sup>v</sup>*(*k*) Controller Dump

valves

Speed sensors

Knowledge techniques (Caldara et al. (1996); Cheok & Shiomi (2000); Mauer (1995); Sanz &

The WSP control system using Fuzzy Logic controllers implemented in Matlab Simulink is shown in Fig. 7. The control system is realized as Triggered Subsystem Simulink block, triggered every 100 ms, which corresponds to a real controller program cycle time. The controller blocks for every axle are implemented with Fuzzy Logic Controller Simulink blocks. The inputs for each FLC block are absolute axle slide *σ* [km/h] and axle acceleration *a*, which are calculated on the basis of the measured wheel circumferential velocities and vehicle reference velocity *vref* , as shown in Fig. 6. The measured wheel circumferential velocities are obtained by feeding the calculated continuous wheel velocities (see section 5.2.3) through a model of speed sensors. Fuzzy Logic Controllers are of PI type. The Simulink subsystem determining the reference velocity *vref* is shown in Fig. 8. The output of each FLC is an integer in the range from −3 to 3, designating one of seven levels of pressure change (see section 3), which is then transformed in a post-processing block into signals controlling a dump valve.

Fig. 7. Fuzzy Logic WSP control system implemented in Matlab Simulink

Fig. 8. Matlab Simulink subsystem determining the reference velocity *vref*

**5.2.2 Model of a braking system**

• model of the lever system • model of the friction elements.

equation (Kaluba (1999)):

equations: for filling

for venting

with:

and for holding

considered as part of this system.

5.2.2.2 Model of the lever system

independently for filling *TF* and venting *TV*.

cylinders of the axle set and its braking torque.

This model can be divided into the following three subsystems:

• model of the pneumatic system and the dump valves

5.2.2.1 Model of the pneumatic system and the dump valves

The model of the braking system describes the relation between pressure in the brake

Matlab Simulink® Model of a Braked Rail Vehicle and Its Applications 197

The pressure at the input of a dump valve can be described with the following

The dump valves are a part of delivery of the WSP system. However, from the modeling point of view, because they are integrated into the pneumatic system of the vehicle, they are

The dump valve has been modeled as the inertia element of adjustable time constant

The pressure in the brake cylinders supplied by the valve are given with one of the following

 1 − *e*

<sup>1</sup> <sup>−</sup> *<sup>e</sup>*−0.75*<sup>t</sup>*

<sup>−</sup> *<sup>t</sup>*−*tz TF* 

. (1)

, (2)

<sup>−</sup> *<sup>t</sup>*−*t*<sup>0</sup> *TV* , (3)

*pc* = *p*0. (4)

*Nk*<sup>1</sup> = (*pcAk* − *Ff*)*ηsηrip*, (5)

*Nk*<sup>2</sup> = (*pcAk* − *Ff*)(2 − *ηsηr*)*ip*. (6)

*pc in* <sup>=</sup> *pcmax*

*pc* = *p*<sup>0</sup> + (*pc in* − *p*0)

*pc* = *p*0*e*

The relation between the pressure of the friction linings of the brake blocks onto the wheel treads of the brake disks, which is exerted by the piston of the brake cylinder and the pressure in the brake cylinder is highly non-linear and is characterized by hysteresis, resulting from friction in the brake cylinder and in the joints of the clamp mechanism (Knorr (2002); Tao & Kokotovic (1996)). This relation is shown in Fig. (9). The two rays limiting the area are given

## **5. Simulator of a braked rail vehicle**

## **5.1 Introduction**

Designing and testing of the FLC WSP controllers demands performing numerous experiments. Performing tests using a real rail vehicle is not feasible, because of high costs and possibility of damaging the object, therefore for designing and testing WSP controllers a test stand is indispensable. One of the most important elements of the test stand is a simulator of a braked rail vehicle.

The mathematical model of a braked rail vehicle, simulation model and various versions of a simulation stand are presented in this subsection.

## **5.2 Mathematical model of a braked rail vehicle**

## **5.2.1 Introduction**

The basis for a rail vehicle simulator is a mathematical model of a braked rail vehicle, taking into consideration the basic phenomena occurring during sliding and omitting the phenomena that are of no or slight significance.

The model has the following parameters:


The inputs to the model are:


The purpose of the model is performing the simulations of braking of a rail vehicle at reduced adhesion. Therefore all the basic phenomena influencing the wheel speeds during braking must be modeled as faithfully as feasible. A special attention must be put to all the main non-linear subsystems of the system. The most important example is the wheel-rail adhesion coefficient.

Some simplifying assumptions have been adopted, which simplify the model, while do not significantly alter its functionality. One of these assumptions is considering the vehicle as a rigid body with the agglomerated mass.

The model consists of the following subsystems:


The subsystems of the model are described in the subsequent subsections.

#### **5.2.2 Model of a braking system**

8 Will-be-set-by-IN-TECH

Designing and testing of the FLC WSP controllers demands performing numerous experiments. Performing tests using a real rail vehicle is not feasible, because of high costs and possibility of damaging the object, therefore for designing and testing WSP controllers a test stand is indispensable. One of the most important elements of the test stand is a simulator

The mathematical model of a braked rail vehicle, simulation model and various versions of a

The basis for a rail vehicle simulator is a mathematical model of a braked rail vehicle, taking into consideration the basic phenomena occurring during sliding and omitting the

• parameters describing properties of the vehicle (e.g. mass, number of axles, inertia of the

The purpose of the model is performing the simulations of braking of a rail vehicle at reduced adhesion. Therefore all the basic phenomena influencing the wheel speeds during braking must be modeled as faithfully as feasible. A special attention must be put to all the main non-linear subsystems of the system. The most important example is the wheel-rail adhesion

Some simplifying assumptions have been adopted, which simplify the model, while do not significantly alter its functionality. One of these assumptions is considering the vehicle as a

The subsystems of the model are described in the subsequent subsections.

• intensiveness of braking defined as the maximum pressure in the brake cylinders.

**5. Simulator of a braked rail vehicle**

simulation stand are presented in this subsection.

**5.2 Mathematical model of a braked rail vehicle**

phenomena that are of no or slight significance.

• vehicle velocity in the moment of commencement of braking

• state of dump valves generated by the WSP controller

The model has the following parameters:

• state of the rail (adhesion coefficient).

rigid body with the agglomerated mass.

• model of a braking system

• model of the adhesion curve

The model consists of the following subsystems:

• model of a rotational motion of a wheel-set

• model of a translation vehicle motion.

**5.1 Introduction**

of a braked rail vehicle.

**5.2.1 Introduction**

axle sets)

coefficient.

The inputs to the model are:

The model of the braking system describes the relation between pressure in the brake cylinders of the axle set and its braking torque.

This model can be divided into the following three subsystems:


5.2.2.1 Model of the pneumatic system and the dump valves

The pressure at the input of a dump valve can be described with the following equation (Kaluba (1999)):

$$p\_{c\_{\rm ini}} = p\_{c\_{\rm max}} \left( 1 - e^{-0.75t} \right). \tag{1}$$

The dump valves are a part of delivery of the WSP system. However, from the modeling point of view, because they are integrated into the pneumatic system of the vehicle, they are considered as part of this system.

The dump valve has been modeled as the inertia element of adjustable time constant independently for filling *TF* and venting *TV*.

The pressure in the brake cylinders supplied by the valve are given with one of the following equations:

for filling

$$p\_{\varepsilon} = p\_0 + (p\_{\varepsilon\_{\rm in}} - p\_0) \left(1 - e^{-\frac{l - l\_{\rm r}}{T\_{\rm F}}}\right),\tag{2}$$

for venting

$$p\_{\mathbb{C}} = p\_0 e^{-\frac{t-t\_0}{T\_V}},\tag{3}$$

and for holding

$$p\_{\mathbb{C}} = p\_0. \tag{4}$$

#### 5.2.2.2 Model of the lever system

The relation between the pressure of the friction linings of the brake blocks onto the wheel treads of the brake disks, which is exerted by the piston of the brake cylinder and the pressure in the brake cylinder is highly non-linear and is characterized by hysteresis, resulting from friction in the brake cylinder and in the joints of the clamp mechanism (Knorr (2002); Tao & Kokotovic (1996)). This relation is shown in Fig. (9). The two rays limiting the area are given with:

$$N\_{k1} = (p\_c A\_k - F\_f) \eta\_s \eta\_r i\_{p\_f} \tag{5}$$

$$N\_{k2} = (p\_c A\_k - F\_f)(2 - \eta\_s \eta\_r) i\_p. \tag{6}$$

5.2.2.3 Model of the friction elements

Because some vehicles are equipped with disk brakes and some with block brakes, models of both types of brakes have been prepared.

For both types of the brakes the braking torque is given with:

$$M\_{\flat} = N\_k \mu r\_{\flat \prime} \tag{7}$$

<sup>20</sup> <sup>0</sup> <sup>60</sup> <sup>40</sup> <sup>100</sup> <sup>80</sup> <sup>140</sup> <sup>120</sup> <sup>160</sup>

Fig. 10. A surface of the friction coefficient *μ* vs. circumferential speed of the vehicle wheel

Matlab Simulink® Model of a Braked Rail Vehicle and Its Applications 199

v [km/h]

<sup>0</sup> <sup>5</sup> <sup>10</sup> <sup>15</sup> <sup>20</sup> <sup>25</sup> <sup>30</sup> <sup>35</sup> <sup>40</sup> <sup>45</sup> <sup>0</sup>

Fig. 11. A surface of the friction coefficient *μ* vs. circumferential speed of the vehicle wheel and the sum of unitary pressure of the friction linings for a disk brake (Kaluba (1999))

Nk [kN]

and the unitary pressure of the friction linings (Sachs (1973))

0.05

0.2 0.25 0.3 0.35 0.4 0.45 0.5

μ [−]

0.1

0.15

0.2

0.25

μ [−]

0.3

0.35

0.4

0.45

0

2 4 6

p x [daN/cm2 ]

50

100

150

v [km/h]

200

where *rh* is the braking radius.

For the block brake, the braking radius *rh* is equal to the wheel radius *r*, for the disk brake it is equal to the radius of the brake disk *rd*.

The friction coefficient between the friction linings of the brake blocks and the wheel tread or the brake disk *μ* depends in a non-linear way on many factors, first of all on wheel velocity and pressure of the disk block.

A simplified value of the friction coefficient *μ* for a block brake is given with the empirical equation. In this model the equation given in Sachs (1973) has been adopted:

$$
\mu = 0,16\sqrt[3]{\frac{40,25}{v}}\sqrt[3]{\frac{2,461}{p\_x}},\tag{8}
$$

while:

$$p\_{\mathbf{x}} = \frac{N\_{\mathbf{k}}}{A\_{\mathbf{x}}}.\tag{9}$$

In Fig. 10 a surface of the friction coefficient *μ* vs. circumferential speed of the vehicle wheel and the unitary pressure of the friction linings for a block brake is shown.

For the disk brake the friction coefficient *μ* has been determined on the basis of the results presented in Kaluba (1999). In Fig. 11 a surface of the friction coefficient *μ* vs. circumferential speed of the vehicle wheel and the sum of unitary pressure of the friction linings for a disk brake is shown.

10 Will-be-set-by-IN-TECH

*filling*

*Nk*<sup>2</sup>

*venting*

*holding*

*p<sup>F</sup><sup>f</sup> pc*<sup>2</sup> *pc*<sup>1</sup>

Fig. 9. Relation between the air pressure in brake cylinders and the brake blocks pressure

Because some vehicles are equipped with disk brakes and some with block brakes, models of

For the block brake, the braking radius *rh* is equal to the wheel radius *r*, for the disk brake it

The friction coefficient between the friction linings of the brake blocks and the wheel tread or the brake disk *μ* depends in a non-linear way on many factors, first of all on wheel velocity

A simplified value of the friction coefficient *μ* for a block brake is given with the empirical

40, 25 *v* 3 2, 461 *px*

*px* <sup>=</sup> *Nk Ax*

In Fig. 10 a surface of the friction coefficient *μ* vs. circumferential speed of the vehicle wheel

For the disk brake the friction coefficient *μ* has been determined on the basis of the results presented in Kaluba (1999). In Fig. 11 a surface of the friction coefficient *μ* vs. circumferential speed of the vehicle wheel and the sum of unitary pressure of the friction linings for a disk

equation. In this model the equation given in Sachs (1973) has been adopted:

*μ* = 0, 16 <sup>3</sup>

and the unitary pressure of the friction linings for a block brake is shown.

*N<sup>k</sup>*

5.2.2.3 Model of the friction elements

where *rh* is the braking radius.

and pressure of the disk block.

while:

brake is shown.

both types of brakes have been prepared.

is equal to the radius of the brake disk *rd*.

For both types of the brakes the braking torque is given with:

*holding*

*Mb* = *Nkμrb*, (7)

*Nk*<sup>1</sup>

*pc*

, (8)

. (9)

Fig. 10. A surface of the friction coefficient *μ* vs. circumferential speed of the vehicle wheel and the unitary pressure of the friction linings (Sachs (1973))

Fig. 11. A surface of the friction coefficient *μ* vs. circumferential speed of the vehicle wheel and the sum of unitary pressure of the friction linings for a disk brake (Kaluba (1999))

#### **5.2.3 Model of a rotational motion of a wheel-set**

A simplified diagram of forces acting upon a braked axle set is shown in Fig. 12.

Fig. 12. A simplified diagram of forces acting upon a braked axle set

The dynamics of the axle set rotary motion is given with the following equation set:

$$J\frac{d\omega\_i}{dt} = F\_{a\dot{l}}r - M\_{b\dot{v}}\tag{10}$$

$$F\_{\rm ai} = \psi\_i(s\_{i\prime}\mathbf{p})Q\_{i\prime} \tag{11}$$

$$s\_i = \frac{v\_T - v}{v\_T},\tag{12}$$

*s*[*−*]

, (13)

. (14)

wheel lock

1

*ψ*

*ψ<sup>α</sup>*

*α*

*ψ*<sup>B</sup>

*ψ*<sup>A</sup>

*ψ*l

0

slide *s* (ORE (1985))

• state of the rail • vehicle velocity *vT* • absolute wheel slide *σ* • relative wheel slide *s*

and

1998); ORE (1985; 1990)):

*s<sup>α</sup> s*<sup>A</sup> *s*<sup>B</sup>

• circumferential deceleration of the wheels *d*

(1985; 1990), is described in detail in Barna (2009).

• slide energy *E* developed in the wheel and a rail contact point.

A

B

Fig. 13. Generalized characteristics of instantaneous coefficient of adhesion *ψ* versus relative

Matlab Simulink® Model of a Braked Rail Vehicle and Its Applications 201

The shape and parameters of the curve depend on many factors, most of which are difficult or impossible to establish. The most important factors are listed below (Boiteux (1987; 1990;

The mathematical model of adhesion coefficient, based on Boiteux (1987; 1990; 1998); ORE

As an example, according to Boiteux (1987), maximum exploitable wheel-rail adhesion coefficient *ψ<sup>B</sup>* and optimum relative slide *sB* can be estimated with the following equations:

One of the most important behavior of the adhesion curve, which has been taken into consideration is, that during the braking at low adhesion with efficiently operating WSP

*Q ψα Eo* 2 *dm*

> 2 *Eo dm Q ψα*

*ψ<sup>B</sup>* = *k*

*sB* <sup>=</sup> <sup>1</sup> *vT*

where (11) defines the adhesion force, and (12) defines the relative slide. The model of value of *ψ* is described in the following subsection.

#### **5.2.4 Model of the adhesion curve**

Generalized characteristics of instantaneous adhesion coefficient *ψ* versus relative slide *s* is shown in Fig. 13. Value of adhesion coefficient at point B (*ψB*) is called *maximum exploitable adhesion* and corresponding slide (*sB*) is called *optimal slide* (Boiteux (1987)).

The parameters of the generalized curve are the following values:


Fig. 13. Generalized characteristics of instantaneous coefficient of adhesion *ψ* versus relative slide *s* (ORE (1985))

The shape and parameters of the curve depend on many factors, most of which are difficult or impossible to establish. The most important factors are listed below (Boiteux (1987; 1990; 1998); ORE (1985; 1990)):

• state of the rail

12 Will-be-set-by-IN-TECH

*Q*

*M<sup>b</sup>*

*vT*

*dt* <sup>=</sup> *Fair* <sup>−</sup> *Mbi*, (10) *Fai* = *ψi*(*si*, **p**)*Qi*, (11)

, (12)

*N*

*si* <sup>=</sup> *vT* <sup>−</sup> *<sup>v</sup> vT*

where (11) defines the adhesion force, and (12) defines the relative slide. The model of value

Generalized characteristics of instantaneous adhesion coefficient *ψ* versus relative slide *s* is shown in Fig. 13. Value of adhesion coefficient at point B (*ψB*) is called *maximum exploitable*

A simplified diagram of forces acting upon a braked axle set is shown in Fig. 12.

*r*

The dynamics of the axle set rotary motion is given with the following equation set:

**5.2.3 Model of a rotational motion of a wheel-set**

*Fa*

of *ψ* is described in the following subsection.

**5.2.4 Model of the adhesion curve**

• relative slide at point *α s<sup>α</sup>*

• relative slide at point *A sA*

• optimum relative slide *sB*

• available adhesion coefficient *ψα*

• adhesion coefficient value at point *A ψ<sup>A</sup>*

• adhesion coefficient value for wheel lock *ψl*.

Fig. 12. A simplified diagram of forces acting upon a braked axle set

*J dω<sup>i</sup>*

*adhesion* and corresponding slide (*sB*) is called *optimal slide* (Boiteux (1987)).

The parameters of the generalized curve are the following values:

• maximum exploitable wheel-rail adhesion coefficient *ψ<sup>B</sup>*

*ω*


The mathematical model of adhesion coefficient, based on Boiteux (1987; 1990; 1998); ORE (1985; 1990), is described in detail in Barna (2009).

As an example, according to Boiteux (1987), maximum exploitable wheel-rail adhesion coefficient *ψ<sup>B</sup>* and optimum relative slide *sB* can be estimated with the following equations:

$$
\psi\_B = k \sqrt{\frac{Q \,\psi\_\alpha \, E\_o}{2 \, d\_m}},
\tag{13}
$$

and

$$s\_B = \frac{1}{v\_T} \sqrt{\frac{2 \, E\_o \, d\_m}{Q \, \psi\_\alpha}}. \tag{14}$$

One of the most important behavior of the adhesion curve, which has been taken into consideration is, that during the braking at low adhesion with efficiently operating WSP

• adhesion variation as occurs in real life

• drag braking test (when a vehicle is hauled at constant speed).

or dead leaves

**5.3.2 Simulation model**

**5.3.3 Simulation stands**

sensors.

laboratory test stand have been shown.

controlled by WSP system • maximum rail slope of 50‰

• vehicle velocity up to 240 km/h

as well as braking force characteristics

(for the braking systems independent of adhesion)

brake independent of adhesion, e.g. track brake

(electrodynamic brake cooperating with friction brake).

• sudden changes of adhesion as the ones occurring when a wheel encounters a spilled oil

Matlab Simulink® Model of a Braked Rail Vehicle and Its Applications 203

• dynamic variation of wheel-rail adhesion coefficient vs. vehicle velocity and the slide

• changes of parameters such as wheel diameter and inertia of axle sets, mass of the vehicle and particular bogies as well as location of the centroid, loading with passenger or cargo

• various brake positions with braking rate values *λ* being within the range 25% do 200%

• braking systems which, apart from the friction brake are additionally equipped with the

It is recommended, that for the traction vehicles it should be possible to simulate blending

The simulation model of the braked rail vehicle integrated with the WSP system has been implemented in Matlab Simulink. A simplified block diagram of the model is shown in Fig. 15. Validation of the model is basically performed by comparing the braking distances obtained from simulations performed with the experimental results (CEN (2009); UIC (2005)).

In the Figures 16, 17 i 18 block diagrams showing various possibilities of realizing the

In Fig. 16 a stand described in UIC (2005) and CEN (2009) is shown, in which the computer simulation is minimized. This stand comprises a hardware model of a brake system with dump valves installed and pressure sensors mounted in brake cylinders. These pressure values are read by the analog inputs of the computer vehicle simulator, which calculates the adhesion forces and resulting wheel velocities. Fast analog output or counter output board of the computer vehicle simulator outputs pulse train simulating impulses from the WSP speed

A real WSP controlling device inputs the simulated speed signals from the computer vehicle simulator and outputs signals controlling dump valves in the hardware model of a brake system. The advantage of this stand is fidelity of simulating the pneumatic part of the vehicle braking system, which makes possible obtaining accurate values of brake cylinder pressures

In Fig. 17 another version of a test stand is shown in which, comparing to the previous stand, also the pneumatic system is computer simulated. A reliable model of a pneumatic system

during WSP operation. The disadvantage is relatively high cost of the stand.

Simulator should make possible modeling the vehicles of the following parameters:

system, the value of the maximum exploitable wheel-rail adhesion coefficient *ψ<sup>B</sup>* increases, so it can double or even triple within 10 s.

Simulated characteristics of adhesion coefficient *ψ* vs. time for a braking from the initial velocity of 120 km/h, showing increase of the maximum exploitable wheel-rail adhesion coefficient in the course of braking is presented in Fig. 14

Fig. 14. Simulated characteristics of adhesion coefficient *ψ* vs. time, where *ψ<sup>i</sup>* – adhesion coefficient for i-th axle set (Barna (2009))

#### **5.2.5 Model of a translation vehicle motion**

The model of a translation vehicle motion is described with:

$$m\frac{dv\_T}{dt} + \sum\_{i=1}^{n} \frac{f\_i}{r} \frac{d\omega\_i}{dt} = -F\_{a\Sigma} \tag{15}$$

where *i* is the number of axle sets of the vehicle.

#### **5.3 Simulation model and a simulation stand**

#### **5.3.1 Requirements**

In UIC (2005) and CEN (2009) there are contained requirements concerning the simulator of a braked rail vehicle. According to these normative documents the simulator should make possible simulating various conditions and scenarios concerning the wheel-rail adhesion as well as vehicle and track parameters:

• adhesion coefficient *ψ* between 0,02 (extremely poor adhesion) and 0,15 (dry rail)


14 Will-be-set-by-IN-TECH

system, the value of the maximum exploitable wheel-rail adhesion coefficient *ψ<sup>B</sup>* increases,

Simulated characteristics of adhesion coefficient *ψ* vs. time for a braking from the initial velocity of 120 km/h, showing increase of the maximum exploitable wheel-rail adhesion

0 5 10 15 20 25 30

t[s]

Fig. 14. Simulated characteristics of adhesion coefficient *ψ* vs. time, where *ψ<sup>i</sup>* – adhesion

*n* ∑ *i*=1

*Ji r dω<sup>i</sup>*

In UIC (2005) and CEN (2009) there are contained requirements concerning the simulator of a braked rail vehicle. According to these normative documents the simulator should make possible simulating various conditions and scenarios concerning the wheel-rail adhesion as

• adhesion coefficient *ψ* between 0,02 (extremely poor adhesion) and 0,15 (dry rail)

*dt* <sup>=</sup> <sup>−</sup>*Fa*Σ, (15)

so it can double or even triple within 10 s.

−0.02

0

coefficient for i-th axle set (Barna (2009))

**5.2.5 Model of a translation vehicle motion**

where *i* is the number of axle sets of the vehicle.

**5.3 Simulation model and a simulation stand**

well as vehicle and track parameters:

**5.3.1 Requirements**

The model of a translation vehicle motion is described with:

*m dvT dt* <sup>+</sup>

0.02

0.04

0.06

ψ [−]

0.08

0.1

0.12

0.14

0.16

coefficient in the course of braking is presented in Fig. 14

ψ1 ψ2 ψ3 ψ4 • drag braking test (when a vehicle is hauled at constant speed).

Simulator should make possible modeling the vehicles of the following parameters:


It is recommended, that for the traction vehicles it should be possible to simulate blending (electrodynamic brake cooperating with friction brake).

### **5.3.2 Simulation model**

The simulation model of the braked rail vehicle integrated with the WSP system has been implemented in Matlab Simulink. A simplified block diagram of the model is shown in Fig. 15. Validation of the model is basically performed by comparing the braking distances obtained from simulations performed with the experimental results (CEN (2009); UIC (2005)).

### **5.3.3 Simulation stands**

In the Figures 16, 17 i 18 block diagrams showing various possibilities of realizing the laboratory test stand have been shown.

In Fig. 16 a stand described in UIC (2005) and CEN (2009) is shown, in which the computer simulation is minimized. This stand comprises a hardware model of a brake system with dump valves installed and pressure sensors mounted in brake cylinders. These pressure values are read by the analog inputs of the computer vehicle simulator, which calculates the adhesion forces and resulting wheel velocities. Fast analog output or counter output board of the computer vehicle simulator outputs pulse train simulating impulses from the WSP speed sensors.

A real WSP controlling device inputs the simulated speed signals from the computer vehicle simulator and outputs signals controlling dump valves in the hardware model of a brake system. The advantage of this stand is fidelity of simulating the pneumatic part of the vehicle braking system, which makes possible obtaining accurate values of brake cylinder pressures during WSP operation. The disadvantage is relatively high cost of the stand.

In Fig. 17 another version of a test stand is shown in which, comparing to the previous stand, also the pneumatic system is computer simulated. A reliable model of a pneumatic system

Fig. 16. Laboratory test stand using a hardware model of a brake system

Fig. 17. Laboratory test stand with computer simulation of the pneumatic system

In Fig. 18 a third test stand is shown, in which a WSP controller is also computer simulated. Such stand can be realized if a control algorithm of a WSP system is known. The biggest difference comparing to the previous stand is, that the fast counters of the WSP system are also simulated. The stand is realized using only a computer with a Matlab Simulink program, thus the advantage is a very low cost of the test stand. The stand is a good choice at the stage of development of WSP control algorithm, because the WSP developer can perform numerous tests of the WSP systems, examining various variants of the control algorithm. However, this stand does not allow for the testing of the ready WSP control device. One of the described above test stands should be used for such purpose. It is this variant of the test stand that is

to UIC (2005) and CEN (2009).

used for simulations described in this chapter.

and of the dump valves is additionally needed. The advantage of the stand is much lower cost. It is also much cheaper in operation, because it does not require supplying with the compressed air. However it cannot be used for homologation tests of WSP systems according

Matlab Simulink® Model of a Braked Rail Vehicle and Its Applications 205

Fig. 15. A simplified block diagram of the braked rail vehicle integrated with the WSP system implemented in Matlab Simulink

16 Will-be-set-by-IN-TECH

100ms

Obliczanie predkosci referencyjnej

4

Model

Model

ukladu hamulcowego

zestawu kolowego

Fig. 15. A simplified block diagram of the braked rail vehicle integrated with the WSP

vzest 1..6 Regulator przeciwposlizgowy 0

4

4

4

4

Przetwornik obr−imp 10ms

4

4

0.06

system implemented in Matlab Simulink

4

4

4 4

4

psi

Losowa zmiana

wspolczynnika

przyczepnosci

Krzywa

4

4

4

Mh1..6

4

4

E0

4

przyczepnosci

lh

1

s

1

vT

s

−1/m

4

4

Tp1..6

4

4

vT

s1..6

STOP 10ms

pmax[kPa] 385

Clock

t

t

t

4

4

4

ZO

pcyl1..6

4

1000

f(u)

p

2

Fcn

Sterowanie zaworow

4

ZL

Stan zaworow

Model zaworow

4

upustowych

Fig. 16. Laboratory test stand using a hardware model of a brake system

and of the dump valves is additionally needed. The advantage of the stand is much lower cost. It is also much cheaper in operation, because it does not require supplying with the compressed air. However it cannot be used for homologation tests of WSP systems according to UIC (2005) and CEN (2009).

Fig. 17. Laboratory test stand with computer simulation of the pneumatic system

In Fig. 18 a third test stand is shown, in which a WSP controller is also computer simulated. Such stand can be realized if a control algorithm of a WSP system is known. The biggest difference comparing to the previous stand is, that the fast counters of the WSP system are also simulated. The stand is realized using only a computer with a Matlab Simulink program, thus the advantage is a very low cost of the test stand. The stand is a good choice at the stage of development of WSP control algorithm, because the WSP developer can perform numerous tests of the WSP systems, examining various variants of the control algorithm. However, this stand does not allow for the testing of the ready WSP control device. One of the described above test stands should be used for such purpose. It is this variant of the test stand that is used for simulations described in this chapter.

120 100 80 60 40 20 0

[km/h]

vT

Matlab Simulink® Model of a Braked Rail Vehicle and Its Applications 207

Fig. 19. Characteristic of vehicle acceleration *aT* vs. vehicle velocity *vT* for simulated braking from 120 km/h at maximum brake cylinder pressure *p*cmax = 385 kPa at good adhesion

wheel accelerations are shown. The numerical data can be used for setting and adjusting the

With a braked rail vehicle simulator it is also possible to perform simulations of braking at poor adhesion conditions aimed at acquiring data concerning the behavior of circumferential wheel speeds and accelerations, which can facilitate the process of developing the WSP control algorithms. In Fig. 21 the characteristics of vehicle velocity and circumferential wheel speed

The rule bases of fuzzy controllers can be designed using expert knowledge and numerical simulation results, provided that such data is available and that the controller which has been used for obtaining the data has controlled the slide in an effective way. In Cheok & Shiomi (2000) a method of designing the fuzzy controller rule base has been proposed, which uses both expert knowledge and chosen numerical results obtained during operation of a PID controller. Basing on the numerical results some rules have been determined, and the lacking rules have been defined on the basis of expert knowledge concerning the slide control. In a simulator of a braked vehicle it is possible to access the signals, which are nor available during normal operation of the systems, i.e. wheel-rail adhesion coefficient and translation velocity of the vehicle. An advantage of having a simulator can be taken to create so called reference

and acceleration during the first cycle of slide detection and suppression.

controllers which can be used for generating data for FLC WSP rule base design.

A block diagram of a WSP control system with a reference controller is shown in Fig. 22.

−1.4

parameters of the WSP control algorithms.

**6.3 Designing the rule bases of FLC WSP**

−1.2

−1

−0.8

a

(Barna (2009))

T [m/s2

]

−0.6

−0.4

−0.2

0

## **6. Application and simulation results**

## **6.1 Introduction**

The simulator of a braked rail vehicle can be used for manifold purposes, which fall into two categories. The first one is developing WSP control algorithms, and it may comprise:


The second category is preliminary testing of the WSP control algorithms or the designed WSP control devices or systems against the requirements given in CEN (2009) and UIC (2005).

In the remaining part of this section the three chosen examples of using the simulation stand for some of the mentioned above purposes are shown.

## **6.2 Acquiring reference data for designing the WSP control algorithms**

By performing simulated brakings at both good and poor adhesion conditions, it is possible to acquire data, which can be used for designing both structure and parameters of WSP controllers. Several examples are given below.

In Fig. 19 the characteristic of vehicle acceleration for simulated braking from 120 km/h at maximum brake cylinder pressure *p*cmax = 385 kPa at good adhesion is shown. Such data, acquired for various brake positions and initial test speed are especially valuable for developing algorithm of calculating the reference speed.

In Fig. 20 the results of simulated braking at poor adhesion conditions (initial value of *ψ<sup>B</sup>* ≈ 0, 06) are presented. The initial test velocity was 160 km/h and the brake position was R (*pcmax* = 385 kPa). In the figure vehicle velocity, circumferential wheel velocity and circumferential 18 Will-be-set-by-IN-TECH

The simulator of a braked rail vehicle can be used for manifold purposes, which fall into two

• acquiring the numerical data concerning the critical values used in the control algorithms

The second category is preliminary testing of the WSP control algorithms or the designed WSP control devices or systems against the requirements given in CEN (2009) and UIC (2005).

In the remaining part of this section the three chosen examples of using the simulation stand

By performing simulated brakings at both good and poor adhesion conditions, it is possible to acquire data, which can be used for designing both structure and parameters of WSP

In Fig. 19 the characteristic of vehicle acceleration for simulated braking from 120 km/h at maximum brake cylinder pressure *p*cmax = 385 kPa at good adhesion is shown. Such data, acquired for various brake positions and initial test speed are especially valuable for

In Fig. 20 the results of simulated braking at poor adhesion conditions (initial value of *ψ<sup>B</sup>* ≈ 0, 06) are presented. The initial test velocity was 160 km/h and the brake position was R (*pcmax* = 385 kPa). In the figure vehicle velocity, circumferential wheel velocity and circumferential

categories. The first one is developing WSP control algorithms, and it may comprise:

Fig. 18. Laboratory test stand with computer simulation of all the systems

• establishing the influence of the changes of the control plant parameters • establishing the influence of the changes of the controller parameters • acquiring the reference data for design of the FLC WSP controllers

**6.2 Acquiring reference data for designing the WSP control algorithms**

**6. Application and simulation results**

• tests of the WSP controllers.

for some of the mentioned above purposes are shown.

developing algorithm of calculating the reference speed.

controllers. Several examples are given below.

**6.1 Introduction**

Fig. 19. Characteristic of vehicle acceleration *aT* vs. vehicle velocity *vT* for simulated braking from 120 km/h at maximum brake cylinder pressure *p*cmax = 385 kPa at good adhesion (Barna (2009))

wheel accelerations are shown. The numerical data can be used for setting and adjusting the parameters of the WSP control algorithms.

With a braked rail vehicle simulator it is also possible to perform simulations of braking at poor adhesion conditions aimed at acquiring data concerning the behavior of circumferential wheel speeds and accelerations, which can facilitate the process of developing the WSP control algorithms. In Fig. 21 the characteristics of vehicle velocity and circumferential wheel speed and acceleration during the first cycle of slide detection and suppression.

#### **6.3 Designing the rule bases of FLC WSP**

The rule bases of fuzzy controllers can be designed using expert knowledge and numerical simulation results, provided that such data is available and that the controller which has been used for obtaining the data has controlled the slide in an effective way. In Cheok & Shiomi (2000) a method of designing the fuzzy controller rule base has been proposed, which uses both expert knowledge and chosen numerical results obtained during operation of a PID controller. Basing on the numerical results some rules have been determined, and the lacking rules have been defined on the basis of expert knowledge concerning the slide control. In a simulator of a braked vehicle it is possible to access the signals, which are nor available during normal operation of the systems, i.e. wheel-rail adhesion coefficient and translation velocity of the vehicle. An advantage of having a simulator can be taken to create so called reference controllers which can be used for generating data for FLC WSP rule base design.

A block diagram of a WSP control system with a reference controller is shown in Fig. 22.

*sB*(*k*) *es*(*k*)

*V O*(*k*)*,*

*<sup>−</sup>* Controller Dump

Measurement of vehicle velocity

*v<sup>T</sup>* (*k*)

a Fuzzy Logic controller, and the other is an expert knowledge based controller.

Fig. 22. Block diagram of a WSP control system with a reference controller (Barna (2009))

In Barna (2010b) a method of designing the rule bases of WSP fuzzy controllers for rail vehicles is presented. In this method results of simulation of two reference controllers are used, both of which use signals which are nor available during normal operation of the systems, i.e. optimal relative slide *sB* and translation velocity of the vehicle *vT*. Two types of reference controllers have been designed in order to obtain data for designing a real-life FLC WSP . One of them is

A block diagram of a Fuzzy Logic reference WSP control system implemented in Matlab Simulink is shown in Fig. 23. It's structure is similar to a regular Fuzzy Logic WSP controller, but the reference value is the optimal relative slide *sB*. A standard Mac-Vicar Whelan rule base

Fig. 23. Fuzzy Logic reference WSP control system implemented in Matlab Simulink

A knowledge based algorithm is based on a MSG1 WSP controller made by Knorr-Bremse

FLC 4

FLC 3

FLC 2

FLC 1

+ *<sup>−</sup> <sup>∗</sup>*

*÷*

*s*(*k*)

has been adopted (Yager & Filev (1995)).

Trigger

and presented in Boiteux (1999).

sB 2 s 1

+ *V L*(*k*) **Δ***pc*(*k*) *v*(*k*)

Rail vehicle

> Valve round 1

valves

Matlab Simulink® Model of a Braked Rail Vehicle and Its Applications 209

Speed sensors

Fig. 20. Results of simulated braking at poor adhesion conditions (initial test velocity 160 km/h, initial value of *ψ<sup>B</sup>* ≈ 0, 06, brake position R (*pcmax* = 385 kPa)) (Barna (2009))

Fig. 21. Characteristics of vehicle velocity and circumferential wheel speed and acceleration for the first cycle of slide detection and suppression (Barna (2009))

20 Will-be-set-by-IN-TECH

0 1 2 3 4 5

v 1 vT

a 1

2.5 <sup>3</sup> −20

−20

−10

0

a[m/s2

]

10

20

−10

0

a[m/s2

]

10

20

t[s]

0 0.5 1 1.5 2 2.5 3 3.5

t[s]

Fig. 21. Characteristics of vehicle velocity and circumferential wheel speed and acceleration

Fig. 20. Results of simulated braking at poor adhesion conditions (initial test velocity 160 km/h, initial value of *ψ<sup>B</sup>* ≈ 0, 06, brake position R (*pcmax* = 385 kPa)) (Barna (2009))

0

120

130

v1 vT

a 1

for the first cycle of slide detection and suppression (Barna (2009))

140

v [km/h]

150

160

50

100

v [km/h]

150

200

Fig. 22. Block diagram of a WSP control system with a reference controller (Barna (2009))

In Barna (2010b) a method of designing the rule bases of WSP fuzzy controllers for rail vehicles is presented. In this method results of simulation of two reference controllers are used, both of which use signals which are nor available during normal operation of the systems, i.e. optimal relative slide *sB* and translation velocity of the vehicle *vT*. Two types of reference controllers have been designed in order to obtain data for designing a real-life FLC WSP . One of them is a Fuzzy Logic controller, and the other is an expert knowledge based controller.

A block diagram of a Fuzzy Logic reference WSP control system implemented in Matlab Simulink is shown in Fig. 23. It's structure is similar to a regular Fuzzy Logic WSP controller, but the reference value is the optimal relative slide *sB*. A standard Mac-Vicar Whelan rule base has been adopted (Yager & Filev (1995)).

Fig. 23. Fuzzy Logic reference WSP control system implemented in Matlab Simulink

A knowledge based algorithm is based on a MSG1 WSP controller made by Knorr-Bremse and presented in Boiteux (1999).

Valve

2−D T[k]

a3

Trigger

tresholds

2−D T[k]

State 4

State 3

0.1

a2

−0.7

a1

−4.4

a

2−D T[k]

10.0

3

Simulink

v

State 2

State 1

2−D T[k]

5.1

a4

State

Matlab Simulink® Model of a Braked Rail Vehicle and Its Applications 211

10

3.6

1

Fig. 26. Block diagram of a knowledge based WSP controller implemented in Matlab

s

>=

>

u(1)\*7/60+3

AND

20

>

u(1)\*11/60+9

<

NOT

15

<

2

>

3.6

2

vref

u(1)\*9/60+6

NOT

AND

3

1

The input controller values for each axle set are:


The output controller values are the same as in the Fuzzy Logic controllers described earlier. Three reference functions are defined, which divide the slide space into four ranges versus the reference velocity *vref* : Δ*σ*1(*vref*), Δ*σ*2(*vref*) and Δ*σ*3(*vref*). These functions are shown in Fig. 24.

Fig. 24. Reference values of absolute slide used in control algorithm of MSG1 WSP system (Boiteux (1999))

Four values of circumferential wheel accelerations are also defined: *a*1, *a*2, *a*<sup>3</sup> oraz *a*4. These values are demonstratively pictured in Fig. 25.

Fig. 25. Values of circumferential wheel accelerations used in control algorithm of MSG1 WSP system pictured demonstratively (Boiteux (1999))

The decision table of MSG1 WSP system control algorithm is presented in Table 1.

A block diagram of a knowledge based WSP control system implemented in Matlab Simulink is shown in Fig. 26.

22 Will-be-set-by-IN-TECH

The output controller values are the same as in the Fuzzy Logic controllers described earlier. Three reference functions are defined, which divide the slide space into four ranges versus the reference velocity *vref* : Δ*σ*1(*vref*), Δ*σ*2(*vref*) and Δ*σ*3(*vref*). These functions are shown in

Fig. 24. Reference values of absolute slide used in control algorithm of MSG1 WSP

Four values of circumferential wheel accelerations are also defined: *a*1, *a*2, *a*<sup>3</sup> oraz *a*4. These

*a*3

Fig. 25. Values of circumferential wheel accelerations used in control algorithm of MSG1

A block diagram of a knowledge based WSP control system implemented in Matlab

The decision table of MSG1 WSP system control algorithm is presented in Table 1.

*a*4

Δ*σ*<sup>1</sup> = *f*(*vref* )

*vref*

*t*

*a*2

Δ*σ*<sup>2</sup> = *f*(*vref* )

Δ*σ*<sup>3</sup> = *f*(*vref* )

The input controller values for each axle set are:

• circumferential wheel acceleration *a*.

• absolute wheel slide *σ*

Δ*σ*

system (Boiteux (1999))

*a*

Simulink is shown in Fig. 26.

values are demonstratively pictured in Fig. 25.

*a*1

WSP system pictured demonstratively (Boiteux (1999))

Fig. 24.

Fig. 26. Block diagram of a knowledge based WSP controller implemented in Matlab Simulink

0 5 10 15 20 25 30 35

Matlab Simulink® Model of a Braked Rail Vehicle and Its Applications 213

v T v1 v2 v3 v4 v ref

vT v1 v 2 v3 v4 v ref

t[s]

0 5 10 15

t[s]

Fig. 29. Characteristics of velocities v, *vT* and *vref* for simulated braking from 50 km/h at maximum brake cylinder pressure (*p*cmax = 385 kPa) at decreased adhesion for a WSP

the following: a passenger coach, a wagon, a locomotive, a train-set or a high speed train. The initial vehicle speed and brake position, as well as additional conditions, are specified for each

Fig. 28. Characteristics of velocities v, *vT* and *vref* for simulated braking from 120 km/h at maximum brake cylinder pressure (*p*cmax = 385 kPa) at decreased adhesion for a WSP

0

controller (Barna (2009))

v [km/h]

controller (Barna (2009))

20

40

60

v [km/h]

80

100

120


Table 1. The decision table of MSG1 WSP system control algorithm (Boiteux (1999))

Simulations have been performed with a simulation model of a braked rail vehicle and controllers implemented in Matlab Simulink. In Fig. 27 exemplary results of simulated braking from 120 km/h with a FLC reference controller at decreased adhesion have been shown.

Fig. 27. Characteristics of velocities v, *vT* and *vref* for simulated braking from 120 km/h at maximum brake cylinder pressure *p*cmax = 385 kPa at decreased adhesion for a FLC reference controller (Barna (2009))

After running a series of tests contained to UIC 541-05 leaflet and EN 15595, simulation data have been used to develop a preliminary rule base using a dedicated C program.

The controller has been tested against the UIC 541-05 requirements with positive results. An exemplary simulated braking is shown in Fig. 28.

#### **6.4 Testing the WSP systems**

When a WSP system has been designed, simulated tests can be performed in order to check whether the WSP meets the requirements of CEN (2009) and UIC (2005). The test program consists of several tests, divided into three groups: slip tests, drag braking test and a test at low adhesion. The exact specification of tests depends on a vehicle type, which can be one of 24 Will-be-set-by-IN-TECH

*σ* a**<sup>1</sup>** → a**<sup>3</sup>** a**<sup>3</sup>** → a**<sup>4</sup>** a**<sup>4</sup>** → a**<sup>4</sup>** a**<sup>4</sup>** → a**<sup>2</sup>** a**<sup>2</sup>** → a**<sup>1</sup>** ↓ σ < **Δ**σ**<sup>1</sup>** H H P3 H P2 **Δ**σ**<sup>1</sup>** <σ< **Δ**σ**<sup>2</sup>** U1 H P2 H P1 **Δ**σ**<sup>2</sup>** <σ< **Δ**σ**<sup>3</sup>** U2 H P1 H H **Δ**σ**<sup>3</sup>** < σ U3 U3 U3 U3 U3

Simulations have been performed with a simulation model of a braked rail vehicle and controllers implemented in Matlab Simulink. In Fig. 27 exemplary results of simulated braking from 120 km/h with a FLC reference controller at decreased adhesion have been

> v T v 1 v2 v3 v 4 v ref

0 5 10 15 20 25 30

t[s]

Fig. 27. Characteristics of velocities v, *vT* and *vref* for simulated braking from 120 km/h at maximum brake cylinder pressure *p*cmax = 385 kPa at decreased adhesion for a FLC reference

After running a series of tests contained to UIC 541-05 leaflet and EN 15595, simulation data

The controller has been tested against the UIC 541-05 requirements with positive results. An

When a WSP system has been designed, simulated tests can be performed in order to check whether the WSP meets the requirements of CEN (2009) and UIC (2005). The test program consists of several tests, divided into three groups: slip tests, drag braking test and a test at low adhesion. The exact specification of tests depends on a vehicle type, which can be one of

have been used to develop a preliminary rule base using a dedicated C program.

Table 1. The decision table of MSG1 WSP system control algorithm (Boiteux (1999))

*a* →

0

exemplary simulated braking is shown in Fig. 28.

controller (Barna (2009))

**6.4 Testing the WSP systems**

20

40

60

v [km/h]

80

100

120

shown.

Fig. 28. Characteristics of velocities v, *vT* and *vref* for simulated braking from 120 km/h at maximum brake cylinder pressure (*p*cmax = 385 kPa) at decreased adhesion for a WSP controller (Barna (2009))

Fig. 29. Characteristics of velocities v, *vT* and *vref* for simulated braking from 50 km/h at maximum brake cylinder pressure (*p*cmax = 385 kPa) at decreased adhesion for a WSP controller (Barna (2009))

the following: a passenger coach, a wagon, a locomotive, a train-set or a high speed train. The initial vehicle speed and brake position, as well as additional conditions, are specified for each

**7. Conclusions**

controllers.

slide and acceleration values.

fitted with I/O boards.

**8. Acknowledgements**

and Higher Education.

*aT* [m/s2] – vehicle acceleration

*e*[−] – control error

*a* [m/s2] – circumferential wheel acceleration

*Ak* [m2] – area of piston of the brake cylinder *Ax* [m2] – area of friction surface of the brake disk

the instant moment *t*

*Fa* [N] – adhesion force for an axle set

the axle set axis

working point (*s*, *ψ*) point B *Ff* [N] – force of the return spring of the brake cylinder

*Fa*<sup>Σ</sup> [N] – sum of adhesion forces for all axle sets of a vehicle *ip* [−] – final ratio of clamp mechanisms of an axle set

**9. List of symbols**

regeneration of adhesion by controlled slide.

In order to design an efficient WSP controller, a simulator of a braked rail vehicle is indispensable. The mathematical model of the wheel-rail adhesion must comprise

Matlab Simulink® Model of a Braked Rail Vehicle and Its Applications 215

The simulator can be used for manifold purposes, including developing and testing of WSP

One of possible applications is using the simulation results of the reference controllers for designing the rule bases for WSP FLC controllers. From the analysis of performance of the WSP fuzzy controller, the rule base of which has been designed in this way, it can be

The test program realized by the simulator must comply with the requirements of CEN (2009) and UIC (2005), thus making possible testing of the controllers in the whole possible range of

• further development of the mathematical model of particular sub-models, especially the

• developing a software/hardware simulator with xPCTarget toolbox using a computer

This paper has been produced as part of the Research Projects "Microprocessor based Anti-slip System for traction rail vehicles meeting the requirements of Technical Specifications of Interoperability" N R10 0046 06/2009 with the financial support of the Ministry of Science

*dm* [m/s2] – average value of circumferential wheel deceleration from the beginning of

*E* [J] – instantaneous value of axle set slide energy, from the start of the braking until

*Eo* [J] – value of axle set slide energy from the start of braking until reaching of the

*J* [kg m2] – moment of inertia of rotating elements associated with an axle set, reduced to

braking and the moment when the slide value reaches 0.4

The futer research concerning the simulator of the braked rail vehicle would comprise:

concluded, that this method makes possible designing efficient controllers.

model of the adhesion coefficient and model of the pneumatic system

test. In order to check whether a WSP meets the requirements of the normative documents, appropriate simulated tests should be performed, and the results assessed.

Testing WSP systems against normative references is described in Barna (2009) and Barna (2010a). In Figures 29, 30 and 31 exemplary test results are shown: in Fig. 29 the results of tests from 50 km/h, and in Fig. 30 the results of tests for rail covered with soap. In Fig. 31 characteristic of pressure in a brake cylinder for simulated braking from 120 km/h is shown.

Fig. 30. Characteristics of velocities v, *vT* and *vref* for simulated braking from 120 km/h at maximum brake cylinder pressure (*p*cmax = 385 kPa) for rail covered with soap for a WSP controller (Barna (2009))

Fig. 31. Characteristics of pressure *pc* for simulated braking from 120 km/h at maximum brake cylinder pressure (*p*cmax = 385 kPa) at decreased adhesion for a WSP controller (Barna (2009))

## **7. Conclusions**

26 Will-be-set-by-IN-TECH

test. In order to check whether a WSP meets the requirements of the normative documents,

Testing WSP systems against normative references is described in Barna (2009) and Barna (2010a). In Figures 29, 30 and 31 exemplary test results are shown: in Fig. 29 the results of tests from 50 km/h, and in Fig. 30 the results of tests for rail covered with soap. In Fig. 31 characteristic of pressure in a brake cylinder for simulated braking from 120 km/h is shown.

0 2 4 6 8 10

t[s]

0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30

t[s]

Fig. 31. Characteristics of pressure *pc* for simulated braking from 120 km/h at maximum brake cylinder pressure (*p*cmax = 385 kPa) at decreased adhesion for a WSP controller (Barna

Fig. 30. Characteristics of velocities v, *vT* and *vref* for simulated braking from 120 km/h at maximum brake cylinder pressure (*p*cmax = 385 kPa) for rail covered with soap for a WSP

40

controller (Barna (2009))

0

1

2

pc1 [kPa]

(2009))

3

4

50

60

70

80

v T v1 v 2 v 3 v 4 vref

v [km/h]

90

100

110

120

appropriate simulated tests should be performed, and the results assessed.

In order to design an efficient WSP controller, a simulator of a braked rail vehicle is indispensable. The mathematical model of the wheel-rail adhesion must comprise regeneration of adhesion by controlled slide.

The simulator can be used for manifold purposes, including developing and testing of WSP controllers.

One of possible applications is using the simulation results of the reference controllers for designing the rule bases for WSP FLC controllers. From the analysis of performance of the WSP fuzzy controller, the rule base of which has been designed in this way, it can be concluded, that this method makes possible designing efficient controllers.

The test program realized by the simulator must comply with the requirements of CEN (2009) and UIC (2005), thus making possible testing of the controllers in the whole possible range of slide and acceleration values.

The futer research concerning the simulator of the braked rail vehicle would comprise:


## **8. Acknowledgements**

This paper has been produced as part of the Research Projects "Microprocessor based Anti-slip System for traction rail vehicles meeting the requirements of Technical Specifications of Interoperability" N R10 0046 06/2009 with the financial support of the Ministry of Science and Higher Education.

## **9. List of symbols**

*a* [m/s2] – circumferential wheel acceleration


Barna, G. (2010b). Simulation based design of fuzzy wheel slide protection controller for

Matlab Simulink® Model of a Braked Rail Vehicle and Its Applications 217

Barna, G. & Kaluba, M. (2009). Safety and reliability of microprocessor brake control systems

Boiteux, M. (1987). Influence de l'énergie de glissment sur l'adhérence exploitable en freinage,

Boiteux, M. (1990). Influence de l'énergie de glissment sur l'adhérence exploitable en freinage,

Caldara, C., Rivera, M. G. & Poma, G. (1996). Software implementation of an anti-skidding

Cheok, A. D. & Shiomi, S. (2000). Combined heuristic knowledge and limited measurement

Kaluba, M. (1999). *Influence of selected factors upon developing of the phenomenon of fluid friction in*

Knorr Bremse AG. (2002). *Bremsen für Schienenfahrzeuge; Handbuch Bremstechnische Begriffe und*

Kwa´snikowski, J. & Firlik, B. (2006). Excessive wear of the wheel treads of the rail bus,

Mauer, G. F. (1995). A fuzzy logic controller for an ABS braking system, *IEEE Transactions on*

Ofierzy ´nski, M. (2008). The method of automatic generation of a model for dynamic-running

ORE (1985). *Adhesion During Braking and Anti-skid Devices, ORE B 164, RP 1, Synthesis of*

ORE (1990). *Adhesion During Braking and Anti-skid Devices, ORE B 164, RP 2, Fundamental Laws*

Pawełczyk, M. (2008). Mathematical model of a system: wheel with a flat spot – rail,

Sachs, K. (1973). *Elektrische Triebfahrzeuge*, Vol. 1, 2 edn, Springer-Verlag, Wien New York. Sanz, M. G. R. R. & Pérez-Rodríguez, J. (1997). An antislipping fuzzy logic controller for a

*Systems, Man, and Cybernetics — Part C: Applications and Reviews* 30(4). Jergéus, J. (1998). Martensite formation and residual stresses around railway wheel flats, *Proc*

Kaczorek, T. (1977). *Theory of Automatic Control Systems (in Polish)*, WNT, Warszawa.

control system for traction electrical drives based on fuzzy-identification techniques, *Symposium on Power Electronics, Industrial Drives, Power Quality, Traction Systems*,

based fuzzy logic antiskid control for railway applications, *IEEE Transactions on*

*the disk brakes of rail vehicles (in Polish)*, PhD thesis, Poznan University of Technology.

*Proceedings of 17th Scientific Conference "Rail Vehicles"*, Kazimierz Dolny, pp. 413–422.

calculations of rail vehicles (in polish), *Proceedings of 18th Scientific Conference "Rail*

*Current Knowledge Concerning Adhesion*, Office for Research and Experiments of the

*of Adhesion in braking*, Office for Research and Experiments of the International Union

*Proceedings of 18th Scientific Conference "Rail Vehicles"*, Vol. 2, Katowice – Ustro ´n,

railway traction system, *Proceedings of the sixth IEEE International Conference on Fuzzy*

Boiteux, M. (1998). Le problème de l'adhérence en freinage, *Chemins de Fer* 5(452): 28–39. Boiteux, M. (1999). Auxiliaires sophistiqués du freinage d'aujourd'hui — les antienrayeurs,

*Automation and Robotics*, Miedzyzdroje.

*Revue Générale des Chemins de Fer* 106(October): 05–15.

*Revue Générale des Chemins de Fer* 109(Juillet - Août): 31–38.

CEN (2009). *EN 15595, Railway applications — Braking — Wheel Slide Protection*.

*and Ecology in Transport)*, Kraków.

*Chemins de Fer* 6(459): 24–35.

Capri, pp. C3–19 – C3–25.

*Werte*, Knorr Bremse AG.

*Fuzzy Systems* 3(4): 381–388.

of Railways, Utrecht.

*Systems*, Barcelona.

pp. 425 – 434.

*Vehicles"*, Vol. 2, Katowice – Ustro ´n, pp. 392 – 412.

International Union of Railways, Utrecht.

*Insts Mech Engrs*, Vol. 212, pp. 69–79.

rail vehicles, *Proceedings of 15th International Conference on Methods and Models in*

for multiple units, *Proceedings of XI International Conference QSEV 2009 (Quality, Safety*


#### **10. References**

Barna, G. (2009). *Control Algorithms of Wheel Slide Protection Systems for Rail Vehicles (in Polish)*, PhD thesis, Poznan University of Technology.

Barna, G. (2010a). Simulation based design and tests of wheel slide protection systems for rail vehicles, *Application of System Science*, Computer Science, Academic Publishing House EXIT, Warsaw, pp. 271–280.

28 Will-be-set-by-IN-TECH

*Nk* [N] – pressure exerted by brake discs, reduced to the braking radius

*pcmax* [Pa] – maximum pressure in a brake cylinder for a given brake position

*pc*<sup>0</sup> [Pa] – outlet dump valve air pressure value at the moment of change of the dump

*η<sup>r</sup>* [−] – coefficient of increasing the efficiency of the lever mechanism of the brake

Barna, G. (2009). *Control Algorithms of Wheel Slide Protection Systems for Rail Vehicles (in Polish)*,

Barna, G. (2010a). Simulation based design and tests of wheel slide protection systems for

rail vehicles, *Application of System Science*, Computer Science, Academic Publishing

*η<sup>s</sup>* [−] – static efficiency of the lever mechanism of the brake clamp mechanism

*Lh* [m] – total braking distance *Mb* [Nm] – braking torque of an axle set

*N* [N] – vertical rail reaction force for an axle set

**p** [−] – vector of parameters of *ψ* vs. *s* characteristics

*pcin* [Pa] – instantaneous air pressure at inlet of a dump valve

*px* [daN/cm2] – unitary pressure of the friction linings of a brake block

*n* [−] – number of axle sets of a vehicle *pc* [Pa] – brake cylinder pressure

valve state

*sB* [−] – optimum relative slide *s<sup>α</sup>* [−] – relative slide for point *α*

*vref* [km/h] – vehicle reference velocity

*γ* [−] – coefficient of rotational mass

*vT* [km/h] – vehicle velocity

*λ* [−] – braking rate

**10. References**

*Q* [N] – instantaneous load of the axle set

*s*[−] – relative slide of wheels against the rails

*Th* [s] – total time of braking until standstill *TV* [s] – time constant of dump valve venting *TF* [s] – time constant of dump valve filling *v* [km/h] – circumferential velocity of axle set wheels

clamp mechanism in motion

*ψ<sup>l</sup>* [−] – adhesion coefficient for locked wheels

*ω* [rad/s] – angular velocity of a vehicle axle set

House EXIT, Warsaw, pp. 271–280.

*σ* [km/h] – absolute slide of wheels against the rails

*ψα* [−] – available adhesion coefficient

*μ* [−] – friction coefficient of the brake block friction lining *ψ* [−] – instantaneous value of wheel-rail adhesion coefficient *ψ<sup>B</sup>* [−] – maximum exploitable wheel-rail adhesion coefficient

*σ*ˆ [km/h] – estimated absolute slide of wheels against the rails

PhD thesis, Poznan University of Technology.

*m* [kg] – vehicle mass

*r* [m] – wheel radius *rb* [m] – brake radius *rd* [m] – brake disk radius

*t* [s] – time


**10** 

*France* 

**Using of Hybrid Supply** 

*Construction Automobile, Cedex 9,* 

Starter-Alternator

Engine Inverter

ECU

Conductor

**for Electric or Hybrid Vehicles** 

*1Ecole Supérieure des Techniques Aéronautiques et de* 

*2Ecole Normale Supérieure (ENS Cachan), Cachan Cedex,* 

Hybrid-Supply

N. Rizoug1, G. Feld2, B. Barbedette1 and R. Sadoun1

For automotive applications, the batteries are sized to ensure many constraints like startup, acceleration, braking and energy recovery. All these constraints give us a very heavy battery with very high energy compared to that required for these applications. To reduce the weight of the storage system, the battery can be associated with high power component like supercapacitors. This last one is used like power booster, and the battery is used just to ensure the energy needed for each application (hybrid or electric vehicles). The use of Matlab-Simulink software allows us to make a modular simulation (Fig. 1). This software resolves the differential equations using several numerical methods (Runge-kutta, Dormand-Prince, Heun, Euler,…). This paper deals with the simulation and conception of hybrid power supply composed with battery and supercapacitors for a micro-hybrid vehicles. In this case, the battery is used as energy tank and supercapacitors as power booster. This design allows to

The hybrid storage system (battery and supercapacitors) supply the starter-alternator through an inverter. Two topologies can be used to associate the battery and the

increase the lifetime of the supply and to downsize this last one.

Fig. 1. Modular simulation of the system using Matlab-Simulink

ICE

**2. Topology of the power system** 

**Starter-Alternato**

**1. Introduction** 

**Load Torque (ICE)**

**Starter-alternator Torque** 


## **Using of Hybrid Supply for Electric or Hybrid Vehicles**

N. Rizoug1, G. Feld2, B. Barbedette1 and R. Sadoun1 *1Ecole Supérieure des Techniques Aéronautiques et de Construction Automobile, Cedex 9, 2Ecole Normale Supérieure (ENS Cachan), Cachan Cedex, France* 

## **1. Introduction**

30 Will-be-set-by-IN-TECH

218 Technology and Engineering Applications of Simulink

Tao, G. & Kokotovic, P. (1996). *Adaptive control of systems with actuator and sensor nonlinearities*,

UIC (2005). *UIC 541-05, Brakes — Specifications for the Construction of Various Brake Parts —*

Yager, R. R. & Filev, D. P. (1995). *Essentials of Fuzzy Modeling and Control*, Wiley-Interscience,

˙ *Int. J.*

Will, A. B. & Zak, S. H. (2000). Antilock brake system modelling and fuzzy control,

John Wiley Sons, Inc., New York.

*Vehicle Design* 24(1): 01–18.

New York.

*Wheel Slide Protection Device (WSP)*, 2 edn.

For automotive applications, the batteries are sized to ensure many constraints like startup, acceleration, braking and energy recovery. All these constraints give us a very heavy battery with very high energy compared to that required for these applications. To reduce the weight of the storage system, the battery can be associated with high power component like supercapacitors. This last one is used like power booster, and the battery is used just to ensure the energy needed for each application (hybrid or electric vehicles). The use of Matlab-Simulink software allows us to make a modular simulation (Fig. 1). This software resolves the differential equations using several numerical methods (Runge-kutta, Dormand-Prince, Heun, Euler,…). This paper deals with the simulation and conception of hybrid power supply composed with battery and supercapacitors for a micro-hybrid vehicles. In this case, the battery is used as energy tank and supercapacitors as power booster. This design allows to increase the lifetime of the supply and to downsize this last one.

Fig. 1. Modular simulation of the system using Matlab-Simulink

## **2. Topology of the power system**

The hybrid storage system (battery and supercapacitors) supply the starter-alternator through an inverter. Two topologies can be used to associate the battery and the

Using of Hybrid Supply for Electric or Hybrid Vehicles 221

Because of their high starting torque, the series DC motors is used like a starter for the automobile applications. The torque developed by this motor can be written as follow:

The relation between the current and the flux of this type of motors complicates the identification of this torque. For that, this motor (starter) is transformed to a separate-wound

Figure 6 presents the transformation of the series DC motor to a separate-wound DC motor.

*T KI starter* (1)

Fig. 4. The transmission between the starter and the ICE

Fig. 5. Electrical diagram of starter (series motor).

This transformation allows us to identify their torque.

**3.1 Formula of the starter Torque** 

DC motor.

supercapacitors: in series or parallel. For this last configuration, the two components are connected to the DC link through two choppers. In the series configuration (fig. 2), the supercapacitors is charged by the battery and discharged during the high power demand (startup, acceleration,…).

Fig. 2. Topology of the series hybrid supply for automotive application

Supercapacitors are chosen like booster because of their very high specific power, which can reach 17kw/kg. The capacitance of this component exceeds the 5000F (Maxwell technology) for a lower voltage (3V). To use this component in power applications, we must make some elements in series.

Fig. 3. Using the hybrid supply with starter or starter-alternator

## **3. Startup torque measurement**

During the starting-up, the ICE develops a load torque applied to the starter or the starteralternator. The computation of this torque allows us to validate the hybrid supply interest. To estimate this torque we need to know the evolution of the starter torque according to its current.

supercapacitors: in series or parallel. For this last configuration, the two components are connected to the DC link through two choppers. In the series configuration (fig. 2), the supercapacitors is charged by the battery and discharged during the high power demand

Supercapacitors are chosen like booster because of their very high specific power, which can reach 17kw/kg. The capacitance of this component exceeds the 5000F (Maxwell technology) for a lower voltage (3V). To use this component in power applications, we must make some

Fig. 2. Topology of the series hybrid supply for automotive application

Fig. 3. Using the hybrid supply with starter or starter-alternator

During the starting-up, the ICE develops a load torque applied to the starter or the starteralternator. The computation of this torque allows us to validate the hybrid supply interest. To estimate this torque we need to know the evolution of the starter torque according to its

**3. Startup torque measurement** 

current.

(startup, acceleration,…).

elements in series.

Fig. 4. The transmission between the starter and the ICE

### **3.1 Formula of the starter Torque**

Because of their high starting torque, the series DC motors is used like a starter for the automobile applications. The torque developed by this motor can be written as follow:

$$T\_{satter} = K \oplus I \tag{1}$$

The relation between the current and the flux of this type of motors complicates the identification of this torque. For that, this motor (starter) is transformed to a separate-wound DC motor.

Fig. 5. Electrical diagram of starter (series motor).

Figure 6 presents the transformation of the series DC motor to a separate-wound DC motor. This transformation allows us to identify their torque.

Using of Hybrid Supply for Electric or Hybrid Vehicles 223

To identify the load torque developed by the ICE during the starting-up phase, tow startingup tests are made. The first it's without the spark plugs and the second it's with

The first test (fig. 8) allows us to estimate the inertia and the friction torque of the ICE. Using

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

Time (S)

Fig. 8. Evolution of the current ant the starter velocity for the test without spark plugs

1.5 3

*J E Kgm f E Nm rd s*

*dt*

0

compression on the ICE motor.

0

*C Nm*

 

*<sup>d</sup> J KI <sup>f</sup> <sup>C</sup>*

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

Time (S)

Using the measured current and starter velocity, we can compute the parameters J, f and C0.

4.2 3 / /

The second test (fig. 9) allows us to estimate the rippled torque (Tcomp) due to the

0

2

(4)

(5)

the mechanic equation we can compute the two parameters (J and f).

**3.2 Formula of the startup torque** 

compression (with spark plug).


Velocity (rd/s)

0

100

200

Starter Current (A)

300

400

Fig. 6. Transformation of the series DC motor to separate-wound DC motor

The acquisition of the armature voltage, the armature current and the starter velocity at deferent wound current allows us to plot the evolution of the ratio K according to the current (fig. 7).

$$E = V - R\,I = K \text{ } \Phi \text{ } \Omega\_{\text{start}} \tag{2}$$

According to the results plotted on the figure 7, the ratio K can be written as follow:

$$K\,\Phi = -2\,\,E-\Psi\,I^3 + 3\,\,E-7\,\,I^2 + 9\,\,E-5\,\,I + 0.0063\tag{3}$$

Fig. 7. Evolution of the ratio K according to the wound current

### **3.2 Formula of the startup torque**

222 Technology and Engineering Applications of Simulink

 Fig. 6. Transformation of the series DC motor to separate-wound DC motor

current (fig. 7).

The acquisition of the armature voltage, the armature current and the starter velocity at deferent wound current allows us to plot the evolution of the ratio K according to the

According to the results plotted on the figure 7, the ratio K can be written as follow:

Fig. 7. Evolution of the ratio K according to the wound current

*E V RI K starter* (2)

3 2 *K EI EI EI* 2 9 3 7 9 5 0.0063 (3)

To identify the load torque developed by the ICE during the starting-up phase, tow startingup tests are made. The first it's without the spark plugs and the second it's with compression (with spark plug).

The first test (fig. 8) allows us to estimate the inertia and the friction torque of the ICE. Using the mechanic equation we can compute the two parameters (J and f).

Fig. 8. Evolution of the current ant the starter velocity for the test without spark plugs

$$
\int \frac{d\Omega}{dt} = K \,\Phi \, I - f \, \, \Omega - \mathbb{C}\_0 \tag{4}
$$

Using the measured current and starter velocity, we can compute the parameters J, f and C0.

$$\begin{cases} f = \text{1.5 E} - \text{3 Kg} \, m^2 \\ f = \text{4.2 E} - \text{3 Nm} / \, rd / s \\ \text{C}\_0 = \text{0 N} m \end{cases} \tag{5}$$

The second test (fig. 9) allows us to estimate the rippled torque (Tcomp) due to the compression on the ICE motor.

Using of Hybrid Supply for Electric or Hybrid Vehicles 225

The supercapacitors sizing is based on the energy and power required for the starting-up of the vehicle. Figure 11 presents the evolution of the current and the voltage of the Buggy starter at the starting-up phase. At the beginning of this phase, the battery voltage deceases to reach 9V. At the same time, the current exceeds the 250A, what gives us a maximum

In order to ensure the functioning of the boost chopper, the minimum voltage of the supercondensator should not decrease lower than 12V (battery voltage). To take a safety margin, the minimum supercapacitors voltage is chosen about 15V. For this component, more than 75% of energy is stored between the maximum voltage of supercapacitors

3

6

9

Voltage (V)

12

15

The power integral on all the starting-up duration allows us to compute the energy required

By the neglecting of the Conv.2 losses, the starting-up energy can be written as follow:

1 3 ( ) <sup>2</sup> <sup>10</sup>

max

min\_ max\_ / 2 max\_ min\_ <sup>2</sup> <sup>30</sup> *VV V V V SuperCaps SuperCaps SuperCaps SuperCaps* (8)

22 2 max\_ min\_ max\_

2 2

The best component for this application can be chosen using the characteristics of all Maxwell cells. The element 2,5V/310F gives us the best compromise between the power and

*<sup>E</sup> C F*

10 10 900 3,34 3 3 <sup>30</sup>

*perCaps Su V perCaps C Vsc SuperCaps* (10)

<sup>250</sup> *E P dt mWh start start* (9)

<sup>0</sup> 0.2 0.4 0.6 0.8 <sup>1</sup> <sup>0</sup>

Time (S)

0

100

200

300

Current (A)

400

500

Current

Voltage

*<sup>V</sup>* (11)

**4. Sizing of the supercapcitors module for the startup** 

Fig. 11. Power and energy required for the starting-up of the vehicle

<sup>0</sup> 0.2 0.4 0.6 0.8 <sup>1</sup> <sup>0</sup>

Time (S)

(Vmax\_SupCaps) and the half of this voltage (Vmax\_SupCaps/2) :

The supercapcitors efficiency reaches 80% for 250A of RMS current.

the energy of the component.

*E CV start sc sc Su*

*start sc*

power of 2400W.

500

1000

1500

power (W)

2000

2500

for this phase.

Fig. 9. Evolution of the current ant the starter velocity for the test with spark plugs

Using the parameters computed at the first test (J, f and C0), the current and the starter velocity presented in the figure 9, we can compute the torque due to the compression (Tcomp):

$$I \int \frac{d\Omega}{dt} = K \,\Phi \, I - f \,\, \Omega - \mathbb{C}\_0 - T\_{comp} \tag{6}$$

Figure 10 presents the wave form of the compression torque Tcomp. This can be written as:

1 5.9sin(0.190\* ) *Tcom starter* (7) <sup>0</sup> 0.2 0.4 0.6 0.8 <sup>1</sup> 1.2 1.4 1.6 1.8 <sup>2</sup> -20 -15 -10 -5 0 5 10 15 20 Time(S) ICE Compresion Torque (Nm)

Fig. 10. Waveform of the compression torque Tcomp

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

Time (S)

Velocity (rd/s)

Starter Current (A)

Fig. 9. Evolution of the current ant the starter velocity for the test with spark plugs

presented in the figure 9, we can compute the torque due to the compression (Tcomp):

Fig. 10. Waveform of the compression torque Tcomp

ICE Compresion Torque (Nm)

*<sup>d</sup> J KIf C T dt*

Figure 10 presents the wave form of the compression torque Tcomp. This can be written as:

1 5.9sin(0.190\* ) *Tcom*

<sup>0</sup> 0.2 0.4 0.6 0.8 <sup>1</sup> 1.2 1.4 1.6 1.8 <sup>2</sup> -20

Time(S)

Using the parameters computed at the first test (J, f and C0), the current and the starter velocity

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

Time (S)

0 *comp*

(6)

*starter* (7)

### **4. Sizing of the supercapcitors module for the startup**

The supercapacitors sizing is based on the energy and power required for the starting-up of the vehicle. Figure 11 presents the evolution of the current and the voltage of the Buggy starter at the starting-up phase. At the beginning of this phase, the battery voltage deceases to reach 9V. At the same time, the current exceeds the 250A, what gives us a maximum power of 2400W.

Fig. 11. Power and energy required for the starting-up of the vehicle

In order to ensure the functioning of the boost chopper, the minimum voltage of the supercondensator should not decrease lower than 12V (battery voltage). To take a safety margin, the minimum supercapacitors voltage is chosen about 15V. For this component, more than 75% of energy is stored between the maximum voltage of supercapacitors (Vmax\_SupCaps) and the half of this voltage (Vmax\_SupCaps/2) :

$$V\_{\text{min\\_SuperCays}} = V\_{\text{max\\_SuperCays}} / \text{\textdegree 2} \implies V\_{\text{max\\_SuperCays}} = \text{2} \; V\_{\text{min\\_SuperCays}} = \text{30}V \tag{8}$$

The power integral on all the starting-up duration allows us to compute the energy required for this phase.

$$E\_{start} = \int P\_{start}dt = 250 \, mV \, \text{lb} \tag{9}$$

The supercapcitors efficiency reaches 80% for 250A of RMS current.

By the neglecting of the Conv.2 losses, the starting-up energy can be written as follow:

$$E\_{\rm start} = \eta\_{sc} \frac{1}{2} \,\mathrm{C\_{sc}} \left( V\_{\rm max\\_SuperCaps}^2 - V\_{\rm min\\_SuperCaps}^2 \right) = \frac{3}{10} \,\mathrm{C\_{sc}} \, V\_{\rm max\\_SuperCaps}^2 \tag{10}$$

$$C\_{sc} = \frac{10}{3} \frac{E\_{start}}{V\_{max}^2} = \frac{10}{3} \frac{900}{30^2} \approx 3,34F \tag{11}$$

The best component for this application can be chosen using the characteristics of all Maxwell cells. The element 2,5V/310F gives us the best compromise between the power and the energy of the component.

Using of Hybrid Supply for Electric or Hybrid Vehicles 227

supply the starter-alternator. In the startup phase, the velocity of the starter-alternator is controlled to reach the set point value (70 rd/S). The load torque is expressed according to

To avoid the angular representation, the starter-alternator is modeled with two phases inputs/outputs (d,q). The model is transformed into three phases inputs/outputs using the Park's transformations. The rotoric current and the load torque generates by the ICE are

considered also like inputs. The starter alternator equations can write as follow:

*di V R i L pL i dt*

*sd sd a sd sd sq sq*

*sq sq a sq sq sd sd af f*

*di V R i L pL i kpM I dt*

<sup>3</sup>

*C p kM I I L L I I*

*em af f sq sd sq sd sq*

(15)

the velocity using the formula 7.

Fig. 13. Test bench developed on the mechatronic laboratory

2

Fig. 14. Model of the starter-alternator


Table 1. Supercapacitors sizing

$$\mathbf{C}\_{component} = \mathbf{310}\mathbf{F}; \quad \mathbf{V}\_{component} = \mathbf{2}, \mathbf{5}\mathbf{V};$$

$$\mathbf{N}^{\odot}\_{elements-en-series} = \mathbf{12}; \quad \mathbf{N}^{\odot}\_{Brandes} = \mathbf{1} \tag{12}$$

$$Volume = \mathbf{11}\mathbf{i}; \quad RMSE = \mathbf{26}\, m\Omega;$$

To reach 30V for the supercapacitors module, we must make 12 elements in series, what give a capacitance module about 26F:

$$C\_{\text{mod}ulle} = \frac{\text{310F}}{\text{12}} = \text{26F} \tag{13}$$

With this component, we can make 8 starting-up without reloading the component:

$$N\_{start} = \frac{\text{26F}}{\text{3,34F}} \approx 8 \text{ starting} - \text{up} \tag{14}$$

#### **5. Simulation of the system**

Before the development of the test bench for this application, a simulation of the operated system is carried out using Matlab-Simulink software. In this case an inverter is used to

Fig. 12. Simulink model of the global system

310 ; 2,5 ; *C FV V component component*

*Volume l RMS m* 1 ; 26 ; To reach 30V for the supercapacitors module, we must make 12 elements in series, what

> 310F 26F 12

mod

3,34 *start*

With this component, we can make 8 starting-up without reloading the component:

<sup>26</sup> <sup>8</sup>

Before the development of the test bench for this application, a simulation of the operated system is carried out using Matlab-Simulink software. In this case an inverter is used to

12 ; 1 *N N elements en series Branches* (12)

*C ule* (13)

*<sup>F</sup> <sup>N</sup> starting up <sup>F</sup>* (14)

Table 1. Supercapacitors sizing

give a capacitance module about 26F:

**5. Simulation of the system** 

Fig. 12. Simulink model of the global system

supply the starter-alternator. In the startup phase, the velocity of the starter-alternator is controlled to reach the set point value (70 rd/S). The load torque is expressed according to the velocity using the formula 7.

Fig. 13. Test bench developed on the mechatronic laboratory

To avoid the angular representation, the starter-alternator is modeled with two phases inputs/outputs (d,q). The model is transformed into three phases inputs/outputs using the Park's transformations. The rotoric current and the load torque generates by the ICE are considered also like inputs. The starter alternator equations can write as follow:

$$\begin{aligned} V\_{sd} &= R\_a \dot{i}\_{sd} + L\_{sd} \frac{d\dot{i}\_{sd}}{dt} - pL\_{sq} \dot{i}\_{sq} \,\Omega\\ V\_{sq} &= R\_a \dot{i}\_{sq} + L\_{sq} \frac{d\dot{i}\_{sq}}{dt} + pL\_{sd} \dot{i}\_{sd} \,\Omega + k \, p \, M\_{af} \, I\_f \, \Omega\\ C\_{cm} &= \frac{3}{2} \, p \Big[ k \, M\_{af} \, I\_f \, I\_{sq} + \left(L\_{sd} - L\_{sq}\right) I\_{sd} \, I\_{sq} \, \right] \end{aligned} \tag{15}$$

Fig. 14. Model of the starter-alternator

Using of Hybrid Supply for Electric or Hybrid Vehicles 229

40,7 28

2

*f pI*

> 1 2 3

(19)

(18)

*f*

(17)

*L H L H*

To calculate the parameter Maf, the starter alternator must turn at constant velocity with the variation of the inductor current. In this case, the value of the voltage between two phases

*sd sq*

32 3 2 *af f ab ab af*

Parameter Value unit Lsd 40 H Lsq 40 H Ra 13.3 m Maf 2.5 mH

The inverter is modeled using the formulas between the starter-alternator voltages (Vas,Vbs

*V K <sup>U</sup> V K V K*

 

The voltage of the Dc link is also one of the inputs for the inverter model. On the other hand, this model imposes the current of the hybrid supply. This last one is composed with two storage systems (battery and supercapacitors) associated with a chopper (fig. 17), which is

Two models are used to represent the battery behaviour. The first one is the simplified R-C model which uses just a resistance and capacitance to modelling the battery. The second representation is based on the frequency behaviour of the battery. This method allows the

211 12 1 <sup>3</sup> 1 12

*U U pM I <sup>M</sup>*

 

The applying of this formula gives us the values of Lsd and Lsq:

Table 2 summarizes the values of the starter-alternator parameters:

and Vcs) and the control of the semiconductors (K1,K2 and K3):

*as*

*bs cs*

representation of the component using the Shepherd model (fig. 18)

connected to the battery through an inductive filter.

Fig. 17. Simulink model of the hybrid supply

*bus*

can expressed as follow:

Table 2. Starter alternator parameters

The starter-alternator parameters ( Lsd, Lsq, Ra, Maf) are identified making some tests. The supplying of the two phases of the starter-alternator by two bridges chopper (fig. 15) at high frequency signal (20 kHz) allows to compute the values of the two inductances Lsd and Lsq. These last ones are computed for two positions of the rotor.

Fig. 15. Identification of the starter-alternator inductances

The first position (fig. 16.a) gives us the value of the direct inductance (Lsd), and The second position of the rotor (fig. 16.b) gives us the value of the quadratic inductance (Lsq). The same formula is used to calculate the two inductances:

$$\begin{cases} L\_{sx} = \mathcal{U}\_{bus} \frac{\Delta t}{2 \Delta i} \\ x = d \text{ or } q \end{cases} \tag{16}$$

Fig. 16. Identification of the machine inductances

The starter-alternator parameters ( Lsd, Lsq, Ra, Maf) are identified making some tests. The supplying of the two phases of the starter-alternator by two bridges chopper (fig. 15) at high frequency signal (20 kHz) allows to compute the values of the two inductances Lsd and Lsq.

The first position (fig. 16.a) gives us the value of the direct inductance (Lsd), and The second position of the rotor (fig. 16.b) gives us the value of the quadratic inductance (Lsq). The

> 2 *sx bus <sup>t</sup> L U*

a) Test to identify the Lsd inductance b) Test to identify the Lsq inductance




0

t i

Voltage (V) / Current (A)

20

40

60

80

 

*x d or q*

*i*

(16)

<sup>0</sup> <sup>20</sup> <sup>40</sup> <sup>60</sup> <sup>80</sup> <sup>100</sup> -80

Time (us)

These last ones are computed for two positions of the rotor.

Fig. 15. Identification of the starter-alternator inductances

same formula is used to calculate the two inductances:

Fig. 16. Identification of the machine inductances

<sup>0</sup> <sup>20</sup> <sup>40</sup> <sup>60</sup> <sup>80</sup> <sup>100</sup> -80

Time (us)

t i

Voltage (V) / Current (A)

The applying of this formula gives us the values of Lsd and Lsq:

$$\begin{cases} L\_{sd} = 40,7 \,\mu H\\ L\_{sq} = 28 \,\mu H \end{cases} \tag{17}$$

To calculate the parameter Maf, the starter alternator must turn at constant velocity with the variation of the inductor current. In this case, the value of the voltage between two phases can expressed as follow:

$$\frac{\mathbf{U}\_{ab}}{\sqrt{3}} = \frac{p \, M\_{af} \, \mathbf{I}\_f \, \mathbf{\varDelta}}{\sqrt{2}} \implies M\_{af} = \frac{\sqrt{2}}{\sqrt{3}} \frac{\mathbf{U}\_{ab}}{2 \, \pi \, f \, p \, \mathbf{I}\_f} \tag{18}$$


Table 2 summarizes the values of the starter-alternator parameters:

Table 2. Starter alternator parameters

The inverter is modeled using the formulas between the starter-alternator voltages (Vas,Vbs and Vcs) and the control of the semiconductors (K1,K2 and K3):

$$
\begin{bmatrix} V\_{as} \\ V\_{bs} \\ V\_{cs} \end{bmatrix} = \frac{\mathcal{U}\_{bus}}{3} \cdot \begin{bmatrix} 2 & -1 & -1 \\ -1 & 2 & -1 \\ -1 & -1 & 2 \end{bmatrix} \cdot \begin{bmatrix} K\_1 \\ K\_2 \\ K\_3 \end{bmatrix} \tag{19}
$$

The voltage of the Dc link is also one of the inputs for the inverter model. On the other hand, this model imposes the current of the hybrid supply. This last one is composed with two storage systems (battery and supercapacitors) associated with a chopper (fig. 17), which is connected to the battery through an inductive filter.

Fig. 17. Simulink model of the hybrid supply

Two models are used to represent the battery behaviour. The first one is the simplified R-C model which uses just a resistance and capacitance to modelling the battery. The second representation is based on the frequency behaviour of the battery. This method allows the representation of the component using the Shepherd model (fig. 18)

Using of Hybrid Supply for Electric or Hybrid Vehicles 231

\_ \_

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

*Ls*

(20)

\_

*I U I II I*

*IU U*

1 *bat bat bus*

*chopper out bus bat SupCaps inverter out chopper out*

The starter-alternator velocity is estimated through Hall sensors. Using the three phases positions (fig. 21), given by these sensors we can estimate the velocity and the electric angle

Two simulations of the system functioning are made using the developed model. The first model it's to validate the startup of the vehicle using the hybrid supply. In this case, the supercapacitor discharge and the voltage of the DC link decreases. The second test is to

To turn-on the vehicle, the velocity of the starter-alternator is controlled to reach the set point value (32rd/s). During this phase, the supercapacitor discharge. Figure 22 presents the waveform of the starter-alternat or velocity (real velocity) and that estimated using the hall sensors. This figure shows also that the starter-alternator velocity reaches the set-point value

Figure 23 presents the waveform of the electric angle and that estimated using the hall signals. This result shows the deference between the real value of the electric angle and that

estimated for the low velocity. After that the estimation gives almost the real value.

The relations between the inputs and outputs of the chopper can write as follow:

 

Fig. 21. Starter-alternator position using Hall sensor

**6.1 Startup of the vehicle using the hybrid supply** 

validate the electric power generation and the charge of the battery.

of the starter-alternator

**6. Simulation results** 

after 50 ms.

Fig. 18. Shepherd model of the battery

For the supercapacitors modeling a transmission line model is used. This approach is based on the representation of the charges propagation along the electrodes surface by RC circuits. So, the supercapacitor is considered like a capacitance distributed along the volume, with accessible areas and less accessible ones. The distribution of the capacitance depends on the electrode material and especially on the geometrical properties of the pores where ions have to access to constitute the double layer capacitance. This approach is well represented by an electrical circuit representing the transmission line and using r-c parallels branches. Then, the proposed model will only include 3 parameters (C, rs and R).

Fig. 19. Transmission line model

The modelling of the chopper is made making the formulas between the inputs and outputs of this converter (fig. 20).

Fig. 20. Chopper inputs/outputs

For the supercapacitors modeling a transmission line model is used. This approach is based on the representation of the charges propagation along the electrodes surface by RC circuits. So, the supercapacitor is considered like a capacitance distributed along the volume, with accessible areas and less accessible ones. The distribution of the capacitance depends on the electrode material and especially on the geometrical properties of the pores where ions have to access to constitute the double layer capacitance. This approach is well represented by an electrical circuit representing the transmission line and using r-c parallels branches. Then,

**rs R/n R/n R/n**

The modelling of the chopper is made making the formulas between the inputs and outputs

**Chopper**

**Ubat Ubus**

**C/n C/n C/n**

**I\_bat I\_chopper\_out I\_inverter\_out**

**SuperCaps**

**I\_Supcaps**

the proposed model will only include 3 parameters (C, rs and R).

**L**

Fig. 18. Shepherd model of the battery

Fig. 19. Transmission line model

Fig. 20. Chopper inputs/outputs

of this converter (fig. 20).

**Battery**

**VSC**

**ISC**

The relations between the inputs and outputs of the chopper can write as follow:

$$\begin{cases} I\_{bat} = \left[ \left( \mathcal{U}\_{bat} - \left( \left( \mathbf{1} - \alpha \right) \mathcal{U}\_{bus} \right) \right) \right] \frac{1}{L\_s} \\\\ I\_{chupper\\_out} = \left( \mathbf{1} - \alpha \right) \mathcal{U}\_{bus} \, I\_{bat} \\ I\_{SupCaps} = I\_{interter\\_out} - I\_{chupper\\_out} \end{cases} \tag{20}$$

The starter-alternator velocity is estimated through Hall sensors. Using the three phases positions (fig. 21), given by these sensors we can estimate the velocity and the electric angle of the starter-alternator

Fig. 21. Starter-alternator position using Hall sensor

## **6. Simulation results**

Two simulations of the system functioning are made using the developed model. The first model it's to validate the startup of the vehicle using the hybrid supply. In this case, the supercapacitor discharge and the voltage of the DC link decreases. The second test is to validate the electric power generation and the charge of the battery.

### **6.1 Startup of the vehicle using the hybrid supply**

To turn-on the vehicle, the velocity of the starter-alternator is controlled to reach the set point value (32rd/s). During this phase, the supercapacitor discharge. Figure 22 presents the waveform of the starter-alternat or velocity (real velocity) and that estimated using the hall sensors. This figure shows also that the starter-alternator velocity reaches the set-point value after 50 ms.

Figure 23 presents the waveform of the electric angle and that estimated using the hall signals. This result shows the deference between the real value of the electric angle and that estimated for the low velocity. After that the estimation gives almost the real value.

Using of Hybrid Supply for Electric or Hybrid Vehicles 233

During the startup phase, the supercapacitors current reaches -90A at the beginning and stabilizes after that at -38A. In this case, the supercapacitors discharge and the voltage of this component deceases. During this phase, the battery charges the supercapacitors at low

<sup>0</sup> 0.1 0.2 0.3 0.4 0.5 0.6 -100

<sup>0</sup> 0.1 0.2 0.3 0.4 0.5 0.6 <sup>27</sup>

Time (s)

<sup>0</sup> 0.1 0.2 0.3 0.4 0.5 0.6 -600

Time (s)

Fig. 24. Current of the starter-alternator during the start-up phase


Phase current


SupCaps voltage (V)

SupCaps Current (A)

0

50

(A)

current (20A) and the SOC of this component decreases from 100% to 99.8%.

Fig. 25. Current and voltage of the supercapacitor module during the start-up phase

Fig. 22. Velocity of the starter-alternator during the start-up phase

Fig. 23. Angle of the starter-alternator during the start-up phase

Figure 24 presents the evolution of starter-alternator current during the startup phase. This result shows that the current for the startup phase reach 560A at the beginning and stabilize after that at 280A.

<sup>0</sup> 0.1 0.2 0.3 0.4 0.5 0.6 <sup>0</sup>

Time (s)

Fig. 22. Velocity of the starter-alternator during the start-up phase

5

1

2

3

4

Angle(rd)

5

6

7

10

15

Velocity (rd/s)

20

25

30

35

Fig. 23. Angle of the starter-alternator during the start-up phase

after that at 280A.

Figure 24 presents the evolution of starter-alternator current during the startup phase. This result shows that the current for the startup phase reach 560A at the beginning and stabilize

<sup>0</sup> 0.1 0.2 0.3 0.4 0.5 0.6 <sup>0</sup>

Time (s)

Fig. 24. Current of the starter-alternator during the start-up phase

During the startup phase, the supercapacitors current reaches -90A at the beginning and stabilizes after that at -38A. In this case, the supercapacitors discharge and the voltage of this component deceases. During this phase, the battery charges the supercapacitors at low current (20A) and the SOC of this component decreases from 100% to 99.8%.

Fig. 25. Current and voltage of the supercapacitor module during the start-up phase

Using of Hybrid Supply for Electric or Hybrid Vehicles 235

<sup>0</sup> <sup>2</sup> <sup>4</sup> <sup>6</sup> <sup>8</sup> <sup>10</sup> -100

0


8 10 12

51

Battery SOC (%)

53

0

Battery current (A)

Battery voltage (V)

20

SupCaps voltage (V)

100

SupCaps current (A)

200

300

Fig. 27. Current and voltage of the supercapacitor module during the power generation phase

<sup>0</sup> <sup>2</sup> <sup>4</sup> <sup>6</sup> <sup>8</sup> <sup>10</sup> <sup>24</sup>

Time (s)

0 2 4 6 8 10

<sup>0</sup> <sup>2</sup> <sup>4</sup> <sup>6</sup> <sup>8</sup> <sup>10</sup> <sup>6</sup>

<sup>0</sup> <sup>2</sup> <sup>4</sup> <sup>6</sup> <sup>8</sup> <sup>10</sup> <sup>49</sup>

Time (s)

Fig. 28. Current, voltage and SOC of the battery during the power generation phase

Fig. 26. Current, voltage and SOC of the battery during the start-up phase

## **6.2 Electric power generation**

When the engine is turned-on, the starter alternator charges the battery and the supercapacitors. Firstly, the supercapacitors is charged with high current to reach 30V (fig. 27). In the same time, the chopper insures the charge of the battery. In this case, the voltage of the battery is controlled to reach the set point value (13V). During this time, the state of charge increases (fig. 28).

The battery current is limited between -20A and +20A. At the beginning, the battery participate to charge the supercapacitor module and after that the battery current reaches 20A, in the same time battery SOC increases during this phase.

The ICE drive the starter-alternator and this last one provide the power needed to charge the battery (fig. 28 and 29). The current in this case is very low compared to that of the startup phase (280A for the sturtup and 4A to charge the battery).

Fig. 26. Current, voltage and SOC of the battery during the start-up phase

20A, in the same time battery SOC increases during this phase.

startup phase (280A for the sturtup and 4A to charge the battery).

When the engine is turned-on, the starter alternator charges the battery and the supercapacitors. Firstly, the supercapacitors is charged with high current to reach 30V (fig. 27). In the same time, the chopper insures the charge of the battery. In this case, the voltage of the battery is controlled to reach the set point value (13V). During this time, the state of

<sup>0</sup> 0.1 0.2 0.3 0.4 0.5 0.6 <sup>0</sup>

0 0.1 0.2 0.3 0.4 0.5 0.6

0 0.1 0.2 0.3 0.4 0.5 0.6

Time (s)

The battery current is limited between -20A and +20A. At the beginning, the battery participate to charge the supercapacitor module and after that the battery current reaches

The ICE drive the starter-alternator and this last one provide the power needed to charge the battery (fig. 28 and 29). The current in this case is very low compared to that of the

**6.2 Electric power generation** 

99.8 99.9 100 100.1

5

10

Battery voltage (V)

Battery current (A)

Battery SOC (%)

15


charge increases (fig. 28).

Fig. 27. Current and voltage of the supercapacitor module during the power generation phase

Fig. 28. Current, voltage and SOC of the battery during the power generation phase

Using of Hybrid Supply for Electric or Hybrid Vehicles 237

the lifetime of the supply. In this paper, hybrid supply for micro-hybrid vehicle has been presented. In this case, two storage components are used: supercapacitors like power supply and battery like energy supply. The power density of supercapacitors makes them very

Two architectures are possible to implement the association of supercapacitors and battery: a parallel architecture by connecting the two sources to DC bus through two choppers. The drawback of this architecture is the very large size of the induction coil in series with the supercapacitor module. In the case of series architecture, the size of the induction coil used to smooth the battery current is less important because of the weakness of the current

After the choose of the architecture of the power supply and the identification of the load torque, the parameters of the starter-alternator are identified. In the next part, the sizing of

Matlab-simulink software is used to simulate the supplying of the starter-alternator by the hybrid supply through an inverter. The association of battery and supercapacitors is made using a boost chopper. For the startup phase, the velocity of the starter-alternator is

N. Rizoug (2006), *Modélisation électrique et énergétique des supercondensateurs et méthodes de* 

P. Bartholomeüs, B. Vulturescu, X. Pierre, N. Rizoug, P. Le Moigne (2001), *A 60V-400A test bench for supercapacitor modules*, EPE'2003, Toulouse, Septembre 2003. S. Buller, E. Karden, D. Kok, R.W. De Doncker (2001), *Simulation Of Supercapacitors In Highly* 

N. Rizoug, P. Bartholomeüs, P. Le Moigne and B. Vulturescu (2004), *Electrical and thermal* 

Y. Hyunjae, S. Seung-Ki, P. Yongho and J. Jongchan (2007), *System Integration and Power-Flow* 

C. Plasse, A. Akémakou, P. Armiroli and D. Sebille (2003), *Lalterno-Demarreur, Du Stop And Go Au Groupe Motopropulsseur Mild Hybride*, Prop'Elec 2003, France 2003. J.M. Dubus, P. Masson and C. Plasse (2003), *The Starter Alternator Reversible Systems Of* 

S. B. Ozzturk, B, Akin, H. A. Toliyat and F. Ashrafzadeh (2006), *Low-Cost Direct Torque* 

D. Feroldi, M. Serra and J. Riera (2009), *Design and Analysis of Fuel-Cell Hybrid Systems* 

W. Gao (2005), *Performance Comparison of a Fuel Cell-Battery Hybrid Powertrain and a Fuel Cell-*

IEEE Trans. Ind. Applic, Vol. 44, N° 1, pp. 108–114, Jan.-feb. 2007.

*en grande puissance ,* Thesis, feb 2006, Ecole centrale de Lille.

*Dynamic Applications*, ESV, 2001, Berlin, Germany.

*caractérisation : Application au cyclage d'un module de supercondensateurs basse tension* 

*behaviour of a supercapacitor module: on-line characterization*, ESSCAP'2004, Belfort,

*Management for a Series Hybrid Electric Vehicle Using Supercapacitors and Batteries*,

*control of Permanent Magnet Synchronous motor using Hall-effect sensors*, IEEE APEC

*Oriented to Automotive Applications*, IEEE Trans. Veh. Technol., vol. 58, no. 9, Nov.

*Ultracapacitor Hybrid Powertrain*, IEEE Trans. Veh. Technol., vol. 54, no. 3, May 2005.

interesting for the applications which need high power during short time.

the two storage system (battery and supercapacitors) is presented.

controlled to follow the set point value (32rd/s).

*Valeo*, Valeo Systèmes Electriques.

charge/discharge.

**8. References** 

France.

2006.

2009.

Fig. 29. Waveform of the phase current during the power generation phase (during 2 seconds)

In the power generation phase, the inverter is used like a synchronous rectifier and the signal period depends to the engine velocity.

Fig. 30. Waveform of the current during the power generation phase (during 0.2 seconds )

## **7. Conclusion**

The batteries used for the micro-hybrid vehicles have very high specific energy. Unfortunately, the power of this technology is very low compared to other technologies (Ni-Mh, Li-Ion,…). The association of battery for this application (micro-hybrid) with powerful component like supercapacitors can decrease the weight of the battery and/or ameliorate the lifetime of the supply. In this paper, hybrid supply for micro-hybrid vehicle has been presented. In this case, two storage components are used: supercapacitors like power supply and battery like energy supply. The power density of supercapacitors makes them very interesting for the applications which need high power during short time.

Two architectures are possible to implement the association of supercapacitors and battery: a parallel architecture by connecting the two sources to DC bus through two choppers. The drawback of this architecture is the very large size of the induction coil in series with the supercapacitor module. In the case of series architecture, the size of the induction coil used to smooth the battery current is less important because of the weakness of the current charge/discharge.

After the choose of the architecture of the power supply and the identification of the load torque, the parameters of the starter-alternator are identified. In the next part, the sizing of the two storage system (battery and supercapacitors) is presented.

Matlab-simulink software is used to simulate the supplying of the starter-alternator by the hybrid supply through an inverter. The association of battery and supercapacitors is made using a boost chopper. For the startup phase, the velocity of the starter-alternator is controlled to follow the set point value (32rd/s).

### **8. References**

236 Technology and Engineering Applications of Simulink

Fig. 29. Waveform of the phase current during the power generation phase (during 2 seconds)

<sup>0</sup> <sup>2</sup> <sup>4</sup> <sup>6</sup> <sup>8</sup> <sup>10</sup> -300

Time (s)

In the power generation phase, the inverter is used like a synchronous rectifier and the

Fig. 30. Waveform of the current during the power generation phase (during 0.2 seconds )

<sup>4</sup> 4.05 4.1 4.15 4.2 -15

Time (s)

The batteries used for the micro-hybrid vehicles have very high specific energy. Unfortunately, the power of this technology is very low compared to other technologies (Ni-Mh, Li-Ion,…). The association of battery for this application (micro-hybrid) with powerful component like supercapacitors can decrease the weight of the battery and/or ameliorate

**7. Conclusion** 

signal period depends to the engine velocity.



0

Starter-Alternator current (A)

5

10

15



Starter-Alternator current (A)

0

100

200

300


**11** 

*Algeria* 

**The Uses of Artificial Intelligence for** 

**Electric Vehicle Control Applications** 

This chapter presents a novel speed control design of electric vehicle (EV) to improve the comportment and stability under different road constraints condition. The control circuit using intelligent fuzzy PI controller is proposed. Parameters which guide the functioning of PI controller are dynamically adjusted with the assistance of fuzzy intelligent control. Actually, electric vehicle (EV) including, full cell and hybrid vehicle have been developed very rapidly as a solution to energy and environmental problem. Driven EVs are powered by electric motors through transmission and differential gears, while directly driven vehicles are propelled by in-wheel or, simply, wheel motors [1, 2]. The basic vehicle congurations of this research has two directly driven wheel motors installed and operated inside the driving wheels on a pure EV. These wheel motors can be controlled independently and have so quick and accurate response to the command that the vehicle chassis control or motion control becomes more stable and robust, compared to indirectly driven EVs. Like most research on the torque distribution control of wheel motor, wheel motors [3, 15] proposed a dynamic optimal tractive force distribution control for an EV driven by four wheel motors,

Research has shown that EV control methods such as, PI control are able to perform optimally over the full range of operation conditions and disturbances and it is very effective with constant vehicle torque, Moreover these non-linear vehicle torque are not fixed and change randomly. However EV with conventional PI control may not have satisfactory performance in such fast varying conditions, the system performance deteriorates. In addition to this, it is difficult to select suitable control parameters Kp and Ki in order to achieve satisfactory compensation results while maintaining the stability of EV traction, due to the highly complex, non-linear nature of controlled systems. These are two of the major drawbacks of the PI control. In order to overcome these difficulties, adaptive PI controller by fuzzy control has been applied both in stationary and under roads constraints, and is shown to improve the overall performance of EV.The Direct Torque Control strategy (DTC) is one kind of high performance driving technologies for AC motors, due to its simple structure and ability to achieve fast response of flux and torque has attracted growing interest in the recent years. DTC-SVM with PI controller direct torque control without hysteresis band can effectively reduce the torque ripple, but its system's robustness will be fur there enhanced. DTC-SVM

thereby improving vehicle handling and stability [4, 5].

**1. Introduction** 

Brahim Gasbaoui and Abdelfatah Nasri

*Department of Electrical Engineering,* 

*Bechar University, Faculty of Sciences and Technology,* 


## **The Uses of Artificial Intelligence for Electric Vehicle Control Applications**

Brahim Gasbaoui and Abdelfatah Nasri *Bechar University, Faculty of Sciences and Technology, Department of Electrical Engineering, Algeria* 

## **1. Introduction**

238 Technology and Engineering Applications of Simulink

J. Bauman and M. Kazerani (2008), *A Comparative Study of Fuel-Cell–Battery, Fuel-Cell–*

S. Lu, Keith A. Corzine and M. Ferdowsi (2007), *A Unique Ultracapacitor Direct Integration* 

Y. Yang, J. Liu, T. Wang, K. Kuo and P. Hsu (2007), *An Electric Gearshift With Ultracapacitors* 

Technol., vol. 57, no. 2, Mar. 2008.

Technol., vol. 56, no. 4, JULY 2007.

Trans. Veh. Technol., vol. 56, no. 5, Sep. 2007.

*Ultracapacitor, and Fuel-Cell–Battery–Ultracapacitor Vehicles*, IEEE Trans. Veh.

*Scheme in Multilevel Motor Drives for Large Vehicle Propulsion*, IEEE Trans. Veh.

*for the Power Train of an Electric Vehicle With a Directly Driven Wheel Motor*, IEEE

This chapter presents a novel speed control design of electric vehicle (EV) to improve the comportment and stability under different road constraints condition. The control circuit using intelligent fuzzy PI controller is proposed. Parameters which guide the functioning of PI controller are dynamically adjusted with the assistance of fuzzy intelligent control. Actually, electric vehicle (EV) including, full cell and hybrid vehicle have been developed very rapidly as a solution to energy and environmental problem. Driven EVs are powered by electric motors through transmission and differential gears, while directly driven vehicles are propelled by in-wheel or, simply, wheel motors [1, 2]. The basic vehicle congurations of this research has two directly driven wheel motors installed and operated inside the driving wheels on a pure EV. These wheel motors can be controlled independently and have so quick and accurate response to the command that the vehicle chassis control or motion control becomes more stable and robust, compared to indirectly driven EVs. Like most research on the torque distribution control of wheel motor, wheel motors [3, 15] proposed a dynamic optimal tractive force distribution control for an EV driven by four wheel motors, thereby improving vehicle handling and stability [4, 5].

Research has shown that EV control methods such as, PI control are able to perform optimally over the full range of operation conditions and disturbances and it is very effective with constant vehicle torque, Moreover these non-linear vehicle torque are not fixed and change randomly. However EV with conventional PI control may not have satisfactory performance in such fast varying conditions, the system performance deteriorates. In addition to this, it is difficult to select suitable control parameters Kp and Ki in order to achieve satisfactory compensation results while maintaining the stability of EV traction, due to the highly complex, non-linear nature of controlled systems. These are two of the major drawbacks of the PI control. In order to overcome these difficulties, adaptive PI controller by fuzzy control has been applied both in stationary and under roads constraints, and is shown to improve the overall performance of EV.The Direct Torque Control strategy (DTC) is one kind of high performance driving technologies for AC motors, due to its simple structure and ability to achieve fast response of flux and torque has attracted growing interest in the recent years. DTC-SVM with PI controller direct torque control without hysteresis band can effectively reduce the torque ripple, but its system's robustness will be fur there enhanced. DTC-SVM

The Uses of Artificial Intelligence for Electric Vehicle Control Applications 241

**3. Direct torque control strategy based space vector modulation (SVM-DTC)**  With the development of microprocessors and DSP techniques, the SVM technique has become one of the most important PWM methods for Voltage Source Inverter (VSI) since it gives a large linear control range, less harmonic distortion, fast transient response, and

simple digital implementation. The induction motor stator flux can be estimated by:

0

*s*

Fig. 2. Bloc diagram for DTC strategy based space vector modulation

And electromagnetic torque Tem can be calculated by:

duty cycle of the switches.

( ) *<sup>t</sup>*

*s ds qs* 

<sup>1</sup> tan *qs*

2 2

*ds* 

> 

<sup>3</sup> ( ) <sup>2</sup> *Tem ds qs qs ds p* 

The SVM principle is based on the switching between two adjacent active vectors and two zero vectors during one switching period. It uses the space vector concept to compute the

*qs qs s qs V R i dt* (5)

(6)

(7)

*i i* (8)

method can improve the system robustness, evidently reduce the torque and flux ripple, and effectively improve the dynamical performance. The DC-DC converter is use with a control strategy to assure the energy require for the EV and the propulsion system. The aim of this chapter is to contribute to understanding the intelligent fuzzy PI controller for utility EV tow rear deriving wheel applied direct torque control based space vector modulation under several scenarios.

## **2. Electric vehicle description**

According to Fig. 1 the opposition forces acting to the vehicle motion are: the rolling resistance force Ftire due to the friction of the vehicle tires on the road; the aerodynamic drag force Faero caused by the friction on the body moving through the air ; and the climbing force Fslope that depends on the road slope [1,2 3].

The total resistive force is equal to Fr and is the sum of the resistance forces, as in (1).

$$\mathbf{F}\_{\mathbf{r}} = \mathbf{F}\_{\text{time}} + \mathbf{F}\_{\text{aero}} + \mathbf{F}\_{\text{slopre}} \tag{1}$$

The rolling resistance force is defined by:

$$\mathbf{F}\_{\text{time}} = \mathbf{r} \mathbf{n} \mathbf{g} \mathbf{f}\_{\text{r}} \tag{2}$$

The aerodynamic resistance torque is defined as follows:

$$\mathbf{F\_{aero} = 1/2 \,\mathrm{\upmu\_{dir}} A\_l C\_d v^2} \tag{3}$$

The rolling resistance force is usually modeled as:

$$F\_{\text{slope}} = \text{rng} \sin(\mathbf{a}) \tag{4}$$

Fig. 1. The forces acting on a vehicle moving along a slope.

Where Mv is the total masse of vehicle r is the tire radius, fr is the rolling resistance force constant, g the gravity acceleration, ρAIR is Air density , Cd is the aerodynamic drag coefficient,Af is the frontal surface area of the vehicle, v is the vehicle speed, is the road slope angle. Values for these parameters are shown in Table1.


Table 1. Parameters of the electric vehicle model

method can improve the system robustness, evidently reduce the torque and flux ripple, and effectively improve the dynamical performance. The DC-DC converter is use with a control strategy to assure the energy require for the EV and the propulsion system. The aim of this chapter is to contribute to understanding the intelligent fuzzy PI controller for utility EV tow rear deriving wheel applied direct torque control based space vector modulation under

According to Fig. 1 the opposition forces acting to the vehicle motion are: the rolling resistance force Ftire due to the friction of the vehicle tires on the road; the aerodynamic drag force Faero caused by the friction on the body moving through the air ; and the climbing force

Fr=Ftire+Faero+Fslope (1)

Ftire=mgfr (2)

Faero=1/2 ρairAfCdv2 (3)

Fslope=mgsin() (4)

Where Mv is the total masse of vehicle r is the tire radius, fr is the rolling resistance force constant, g the gravity acceleration, ρAIR is Air density , Cd is the aerodynamic drag coefficient,Af is the frontal surface area of the vehicle, v is the vehicle speed, is the road

r 0.32 m Af 2.60 m2 m 1300 Kg Cd 0.32 fr 0.01 ρair 1.2 Kg/m3

The total resistive force is equal to Fr and is the sum of the resistance forces, as in (1).

several scenarios.

**2. Electric vehicle description** 

Fslope that depends on the road slope [1,2 3].

The rolling resistance force is defined by:

The aerodynamic resistance torque is defined as follows:

Fig. 1. The forces acting on a vehicle moving along a slope.

slope angle. Values for these parameters are shown in Table1.

Table 1. Parameters of the electric vehicle model

The rolling resistance force is usually modeled as:

## **3. Direct torque control strategy based space vector modulation (SVM-DTC)**

With the development of microprocessors and DSP techniques, the SVM technique has become one of the most important PWM methods for Voltage Source Inverter (VSI) since it gives a large linear control range, less harmonic distortion, fast transient response, and simple digital implementation. The induction motor stator flux can be estimated by:

$$
\sigma\_{qs} = \int\_0^t (V\_{qs} - R\_s i\_{qs}) dt \tag{5}
$$

$$\left| \varphi\_s \right| = \sqrt{\varphi\_{ds}^2 + \varphi\_{qs}^2} \tag{6}$$

$$\varphi\_s = \tan^{-1} \left( \frac{\varphi\_{qs}}{\varphi\_{ds}} \right) \tag{7}$$

And electromagnetic torque Tem can be calculated by:

$$T\_{em} = \frac{3}{2} p (\varphi\_{ds} i\_{qs} - \varphi\_{qs} i\_{ds}) \tag{8}$$

The SVM principle is based on the switching between two adjacent active vectors and two zero vectors during one switching period. It uses the space vector concept to compute the duty cycle of the switches.

Fig. 2. Bloc diagram for DTC strategy based space vector modulation

The Uses of Artificial Intelligence for Electric Vehicle Control Applications 243

Fuzzy controllers have been widely applied to industrial process. Especially, fuzzy controllers are effective techniques when either the mathematical model of the system is nonlinear or not the mathematical model exists. In this paper, the fuzzy control system adjusts the parameter of the PI control by the fuzzy rule. Dependent on the state of the system . The adaptive PI realized is no more a linear regulator according to this principle. In most of these studies, the Fuzzy controller used to drive the PI is defined by the authors

<sup>1</sup> ( ) \*[ ( ) ( ) ]

*e t*( ) : Input of the control. The error of the reference current w (t) \* and the injected speed

<sup>1</sup> ( ) \*[ ( ) ( ) ] *<sup>k</sup>*

*ek ek ek* ( ) ( ) ( 1)

<sup>1</sup> ( ) \*[ ( ) ( ) ] *<sup>k</sup>*

<sup>1</sup> ( ) \* () () *<sup>k</sup> <sup>i</sup> <sup>j</sup> <sup>y</sup> k Kp e t K e j*

*y k Kp et e j T <sup>T</sup>*

*<sup>y</sup> t Kp e t e t dt <sup>T</sup>*

0

1

1

*<sup>j</sup> <sup>i</sup>*

*<sup>y</sup> k Kp e k e j T <sup>T</sup>* (11)

*kp Kp* 20 0.8( 2.5) (12)

*ki* 0.0125 0.003( 2.5) *Ki* (13)

*<sup>j</sup> <sup>i</sup>*

(10)

*i*

*t*

**5. Adaptive fuzzy PI controller** 

from a series of experiments [21, 22, 23, 24].

*y t*( ) : Output of the control.

*kp* : Parameter of the scale *Ti* : Parameter of the integrator.

*y k*( ) : Output on the time of k th sampling. *e k*( ) : Error on the time of k sampling

The on-line tuning equation for kp and ki are show above:

The discrete equation:

*T* : Cycle of the sampling

On-line Tuning:

Where:

*w t*( )

Where:

The expression of the PI is given in the equation (10).

## **4. Conventional PI controller**

The reason behind the extensive use of proportional integral (PI) ..controller is its effectiveness in the control of steady-state error of a control system and also its easy implementation. However, one disadvantage of this conventional compensator is its inability to improve the transient response of the system. The conventional PI controller figure 3 has the form of Eq. (9), where \* *Tem* is the control output. Kp and Ki are the proportional and integral gains respectively, these gains depend on the system parameters. ε is the error signal, which is the difference of the injected voltage to the reference voltage.

Fig. 3. Control of the injected speed using conventional PI controller

$$T\_{em}(t) = Kp\varepsilon(t) + K\dot{t} \int\_{T} \varepsilon(t)dt\tag{9}$$

Equation (9) shows that the PI controller introduces a pole in the entire feedback system, consequently, making a change in its original root locus. Analytically the pole introduces a change in the control system's response. The effect is the reduction of steady-state error. On the other hand, the constants Kp and Ki determine the stability and transient response of the system, in which, these constants rely on their universe of discourses: min max min max *Kp Kp Kp and Ki Ki Ki* [ , ] [,]

Where the values of the minimum and maximum proportional and integral constants (gains) are practically evaluated through experimentation and using some iterative techniques. This makes the design of the conventional PI controller dependent on the knowledge of the expert. When the compensator constants exceed the allowable values, the control system may come into an unstable state. After the determination of the domain of the proportional and integral constants, the tuning of the instantaneous values of the constants takes place. Depending on the value of the error signal, ε, the values of the constants adjusts formulating an adaptive control system. The constants Kp and Ki changes to ensure that the steady-state error of the system is reduced to minimum if not zero.

## **5. Adaptive fuzzy PI controller**

Fuzzy controllers have been widely applied to industrial process. Especially, fuzzy controllers are effective techniques when either the mathematical model of the system is nonlinear or not the mathematical model exists. In this paper, the fuzzy control system adjusts the parameter of the PI control by the fuzzy rule. Dependent on the state of the system . The adaptive PI realized is no more a linear regulator according to this principle. In most of these studies, the Fuzzy controller used to drive the PI is defined by the authors from a series of experiments [21, 22, 23, 24].

The expression of the PI is given in the equation (10).

$$\log y(t) = Kp^\* \left[ e(t) + \frac{1}{T\_i} \right] e(t)dt \tag{10}$$

Where:

242 Technology and Engineering Applications of Simulink

The reason behind the extensive use of proportional integral (PI) ..controller is its effectiveness in the control of steady-state error of a control system and also its easy implementation. However, one disadvantage of this conventional compensator is its inability to improve the transient response of the system. The conventional PI controller figure 3 has the form of Eq. (9), where \* *Tem* is the control output. Kp and Ki are the proportional and integral gains respectively, these gains depend on the system parameters. ε is the error signal, which is the difference of the injected voltage to the

Fig. 3. Control of the injected speed using conventional PI controller

min max min max *Kp Kp Kp and Ki Ki Ki* [ , ] [,]

zero.

() () () *em <sup>T</sup>*

Equation (9) shows that the PI controller introduces a pole in the entire feedback system, consequently, making a change in its original root locus. Analytically the pole introduces a change in the control system's response. The effect is the reduction of steady-state error. On the other hand, the constants Kp and Ki determine the stability and transient response of the system, in which, these constants rely on their universe of discourses:

Where the values of the minimum and maximum proportional and integral constants (gains) are practically evaluated through experimentation and using some iterative techniques. This makes the design of the conventional PI controller dependent on the knowledge of the expert. When the compensator constants exceed the allowable values, the control system may come into an unstable state. After the determination of the domain of the proportional and integral constants, the tuning of the instantaneous values of the constants takes place. Depending on the value of the error signal, ε, the values of the constants adjusts formulating an adaptive control system. The constants Kp and Ki changes to ensure that the steady-state error of the system is reduced to minimum if not

 

*t Ki t dt* (9)

*TtK p*

**4. Conventional PI controller** 

reference voltage.

*y t*( ) : Output of the control.

*e t*( ) : Input of the control. The error of the reference current w (t) \* and the injected speed *w t*( )

*kp* : Parameter of the scale

*Ti* : Parameter of the integrator.

The discrete equation:

$$y(k) = kp^\* \left[ \varepsilon(k) + \frac{1}{T\_i} \sum\_{j=1}^k \varepsilon(j)T \right] \tag{11}$$

Where:

*y k*( ) : Output on the time of k th sampling.

*e k*( ) : Error on the time of k sampling

*T* : Cycle of the sampling

$$\begin{aligned} \Delta e(k) &= e(k) - e(k-1) \\\\ y(k) &= Kp^\* \left[ e(t) + \frac{1}{T\_i} \sum\_{j=1}^k e(j)T \right] \\\\ y(k) &= Kp^\* \left[ e(t) + K\_i \sum\_{j=1}^k e(j) \right] \end{aligned}$$

On-line Tuning:

The on-line tuning equation for kp and ki are show above:

$$kp = 20 + 0.8(Kp - 2.5) \tag{12}$$

$$ki = 0.0125 + 0.003(Ki - 2.5) \tag{13}$$

The Uses of Artificial Intelligence for Electric Vehicle Control Applications 245

**Z S M L**

0 1 2 3 4 5 6

kp

**Z S M L**

0 1 2 3 4 5 6

ki


**de e**


e de


4

0


0 0.2 0.4 0.6 0.8 1

> 0 0.2 0.4 0.6 0.8 1

**kp**

Fig. 9. View plot surface of fuzzy controller for kp

ki

Fig. 10. View plot surface of fuzzy controller for ki

**Degree of membership**

Fig. 8. The Membership function of output ki.

**Degree of membership**

Fig. 7. The Membership functions of output kp.

The frame of the fuzzy adaptive PI controller is illustrated in figure 4.

Fig. 4. PI gains online tuning by fuzzy logic controller

The linguistic variables are defines as {NL,NM,NS,Z,PS,PM,PB}meaning negative large, negative medium ,negative small, zero, positive small, positive medium, positive big.

The Membership function is illustrated in the figures 7, 8, 9 and 10.Figures 9 and 10 shows the view plot of fuzzy controller for kp and ki respectively.

Fig. 5. The Membership function of input e(k).

Fig. 6. The Membership function of input *e k*( ) .

Fig. 7. The Membership functions of output kp.

The linguistic variables are defines as {NL,NM,NS,Z,PS,PM,PB}meaning negative large, negative medium ,negative small, zero, positive small, positive medium, positive big.

The Membership function is illustrated in the figures 7, 8, 9 and 10.Figures 9 and 10 shows

**NL NM ZE PM NS PS PL**


e

**N P Z**


de

The frame of the fuzzy adaptive PI controller is illustrated in figure 4.

Fig. 4. PI gains online tuning by fuzzy logic controller

the view plot of fuzzy controller for kp and ki respectively.

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1

**Degree of membership**

Fig. 6. The Membership function of input *e k*( ) .

**Degree of membership**

Fig. 5. The Membership function of input e(k).

Fig. 8. The Membership function of output ki.

Fig. 9. View plot surface of fuzzy controller for kp

Fig. 10. View plot surface of fuzzy controller for ki

The Uses of Artificial Intelligence for Electric Vehicle Control Applications 247

*v v V Rd V Rd*

*v v V Rd V Rd*

*L d*

The difference between wheel drive angular speeds is then:

And the steering angle indicates the trajectory direction:

*<sup>v</sup>* is the vehicle angular speed according to the center of turn.

0 *Straight ahead*

0 *Turnraight*

*L d*

( 2 ( 2

( 2 ( 2

( 2)tan( )

*mr mr*

 

 

*L*

*ml ml*

\* \* tan( ) *mr ml v d L* 

*L*

( 2)tan( )

 

0 *Turnleft* (19)

(15)

(16)

(17)

(18)

1 2

1 2

And their angular speed by:

Where

Fig. 11. Differential electronic.


Table 2. Fuzzy tuning rules

## **6. Implementation of electronic differential**

The proposed control system principle could be summarized as follows:

A speed control is used to control each motor torque. The speed of each rear wheel is controlled using speed difference feedback. Since the two rear wheels are directly driven by two separate motors, the speed of the outer wheel will need to be higher than the speed of the inner wheel during steering maneuvers (and vice-versa). This condition can be easily met if the speed estimator is used to sense the angular speed of the steering wheel. The common reference speed *ref* is then set by the accelerator pedal command. The actual reference speed for the left drive \* *left* and the right drive \* *right* are then obtained by adjusting the common reference speed \* using the output signal from the DTC speed estimator. If the vehicle is turning right, the left wheel speed is increased and the right wheel speed remains equal to the common reference speed \* . If the vehicle is turning left, the right wheel speed is increased and the left wheel speed remains equal to the common reference speed \* [7, 9, and 11]. Usually, a driving trajectory is adequate for an analysis of the vehicle system model. From the mode show in Fig. 6, the following characteristic can be calculated:

$$R\_{\alpha} = \frac{L\_{\alpha}}{t \text{g}(\mathcal{S})} \tag{14}$$

Where is the steering angle. Therefore, the linear speed of each wheel drive is given by:

$$\begin{cases} V\_1 = \alpha\_v (R - d\_o/2) \\ V\_2 = \alpha\_v (R - d\_o/2) \end{cases} \tag{15}$$

And their angular speed by:

246 Technology and Engineering Applications of Simulink

NL NM NS ZE PS PM PB

N L M S M S M L

Z L M L Z L M L

P L M L Z L M L

N Z S M L M S Z

Z Z S M L M S Z

P Z M L L L M Z

*ref* is then set by the accelerator pedal command. The actual

(14)

*right* are then obtained by

*left* and the right drive \*

[7, 9, and 11]. Usually, a driving trajectory is adequate for an analysis of

*e*( ) 

*e*( ) 

**6. Implementation of electronic differential** 

The proposed control system principle could be summarized as follows:

A speed control is used to control each motor torque. The speed of each rear wheel is controlled using speed difference feedback. Since the two rear wheels are directly driven by two separate motors, the speed of the outer wheel will need to be higher than the speed of the inner wheel during steering maneuvers (and vice-versa). This condition can be easily met if the speed estimator is used to sense the angular speed of the steering wheel. The

adjusting the common reference speed \* using the output signal from the DTC speed estimator. If the vehicle is turning right, the left wheel speed is increased and the right wheel speed remains equal to the common reference speed \* . If the vehicle is turning left, the right wheel speed is increased and the left wheel speed remains equal to the common

the vehicle system model. From the mode show in Fig. 6, the following characteristic can be

*<sup>L</sup> <sup>R</sup> tg* 

Where is the steering angle. Therefore, the linear speed of each wheel drive is given

 ( )

kp and ki

kp

ki

Table 2. Fuzzy tuning rules

common reference speed

reference speed \*

calculated:

by:

reference speed for the left drive \*

$$\begin{cases} V\_1 = \alpha\_v (R - d\_o/2) \\ V\_2 = \alpha\_v (R - d\_o/2) \end{cases} \tag{16}$$

$$\begin{aligned} \alpha\_{mr}^{\*} &= \frac{L\_{oo} - (d\_{oo}/2)\tan(\mathcal{S})}{L\_{oo}} \alpha\_{mr} \\ \alpha\_{ml}^{\*} &= \frac{L\_{oo} - (d\_{oo}/2)\tan(\mathcal{S})}{L\_{oo}} \alpha\_{ml} \end{aligned} \tag{17}$$

Where *<sup>v</sup>* is the vehicle angular speed according to the center of turn.

The difference between wheel drive angular speeds is then:

$$
\Delta oo = o\_{mr}^\* - o\_{ml}^\* = -\frac{d\_{o0}\tan(\delta)}{L\_{o0}} o o\_v \tag{18}
$$

And the steering angle indicates the trajectory direction:

Fig. 11. Differential electronic.

$$
\delta > 0 \implies Term{left}\tag{19}
$$

0 *Straight ahead*

0 *Turnraight*

The Uses of Artificial Intelligence for Electric Vehicle Control Applications 249

In order to characterize the driving wheel system behavior, Simulations were carried using the model of Fig 13.The following results was simulated in MATLAB .The test demonstrate the EV performances using an intelligence fuzzy PI controlled with DTC-SVM strategy

The topology studied in this present work consists of three phases: the first one represent the acceleration phase's beginning with 60 Km/h in straight road, the second phase represent the deceleration one when the speed became 30 Km/h, and finally the EV is moving up the slopped road of 10% under 80 Km/h, the specified road topology is shown

Phases Event information Vehicle Speed [km/h]

Phase 1 Acceleration 60

Phase 2 Bridge, Break 30

Phase 3 Acceleration and climbing a slope 10% 80

**A: Intelligent fuzzy PI controller with space vector modulation** 

in Fig. 14, when the speed road constraints are described in the Table 2.

**7. Simulation results** 

under several speed variation.

Fig. 14. Specify driving route topology

Table 3. Specified driving route topology

Fig. 12. Structure of vehicle in curve

Fig. 13. The riving wheels control system

## **7. Simulation results**

248 Technology and Engineering Applications of Simulink

Fig. 12. Structure of vehicle in curve

Fig. 13. The riving wheels control system

In order to characterize the driving wheel system behavior, Simulations were carried using the model of Fig 13.The following results was simulated in MATLAB .The test demonstrate the EV performances using an intelligence fuzzy PI controlled with DTC-SVM strategy under several speed variation.

## **A: Intelligent fuzzy PI controller with space vector modulation**

The topology studied in this present work consists of three phases: the first one represent the acceleration phase's beginning with 60 Km/h in straight road, the second phase represent the deceleration one when the speed became 30 Km/h, and finally the EV is moving up the slopped road of 10% under 80 Km/h, the specified road topology is shown in Fig. 14, when the speed road constraints are described in the Table 2.

Fig. 14. Specify driving route topology


Table 3. Specified driving route topology

The Uses of Artificial Intelligence for Electric Vehicle Control Applications 251

Figs 17 and 18 explains the variation of phase current and driving force respectively. In the first step and to reach 60 km/h The EV demand a current of 50.70 A for each motor which explained with driving force of 329.30N. In second phase the current and driving forces demand decrees by means that the vehicle is in recharging phase's which explained with the decreasing of current demand and developed driving forces shown in Figs 14 and 15 respectively . The last phases explain the effect of acceleration under the slope on the straight road EV moving. The driving wheels forces increase and the current demand undergo double of the current braking phases the battery use the maximum of his power to satisfy the motorization demand under the slopped road condition which can interpreted physically the augmentation of the globally vehicle resistive torque illustrate in Table 4. In the other hand the linear speeds of the two induction motors stay the same and the road drop does not influence the torque control of each wheels. The results are listed in Table 3.

0 1 2 3 4 5 6

Phase 1 Phase 2 Phase 3

Motor right

**Time [s]**

Phases 1 2 3 the Vehicle resistive torque [N.m] 95.31 68.53 168.00

compared with nominal motor torque of 476 Nm 20.02 % 14.39 % 35.29 %

Fig. 19. Evaluation of the globally vehicle resistive torque compared to nominal motor

Fig. 18. Variation of driving force of the right motor in different phases.

0

the globally vehicle resistive torque Percent

Table 4. Variation of vehicle torque in different Phases.

200

400

600

**Driving forces [N]**

torque in different phases.

800

1000

Fig. 15. Variation of vehicle speeds in different phases.

Refereed to Fig. 15 at time of 2 s the vehicle driver move on straight road with linear speed of 60 km/h, the assumption's that the two motors are not disturbed. In this case the driving wheels follow the same path with no overshoot and without error which can be justified with the good electronic differential act coupled with DTC-SVM performances.

Fig. 16. Evaluation of vehicle and distance travelled in different phases.

Fig. 16, reflect the relationship between vehicle speed's variation and distance traveled in different phases. The distance travelled of 310 m in three electronic differential references acts 60 then break of 30 and acceleration until 80 km/h.

Fig. 17. Variation of phase current of the right motor in different phases.

Refrence speed Linear vehicle speed Linear speed of the right wheel Linear speed of the left wheel

<sup>0</sup> <sup>1</sup> <sup>2</sup> <sup>3</sup> <sup>4</sup> <sup>5</sup> <sup>6</sup> <sup>0</sup>

Refereed to Fig. 15 at time of 2 s the vehicle driver move on straight road with linear speed of 60 km/h, the assumption's that the two motors are not disturbed. In this case the driving wheels follow the same path with no overshoot and without error which can be justified

0 1 2 3 4 5 6

Phase 1 Phase 2 Phase 3

**Time [s]**

Fig. 16, reflect the relationship between vehicle speed's variation and distance traveled in different phases. The distance travelled of 310 m in three electronic differential references

0 1 2 3 4 5 6

Phase 1 Phase 2 Phase 3

**Time [s]**

with the good electronic differential act coupled with DTC-SVM performances.

Distance traveled Vehicle speed

Fig. 16. Evaluation of vehicle and distance travelled in different phases.

Fig. 17. Variation of phase current of the right motor in different phases.

acts 60 then break of 30 and acceleration until 80 km/h.

Phase 1 Phase 2 Phase 3

**Time [s]**


0

Motor right

100

**Distance traveled [m]**

200

300

400

0


0

100

**Phase current[A]**

200

300

20

40

**Linear vehicle speed [Km/h]**

60

80

Fig. 15. Variation of vehicle speeds in different phases.

**Vehicle speed [Km/h]**

Figs 17 and 18 explains the variation of phase current and driving force respectively. In the first step and to reach 60 km/h The EV demand a current of 50.70 A for each motor which explained with driving force of 329.30N. In second phase the current and driving forces demand decrees by means that the vehicle is in recharging phase's which explained with the decreasing of current demand and developed driving forces shown in Figs 14 and 15 respectively . The last phases explain the effect of acceleration under the slope on the straight road EV moving. The driving wheels forces increase and the current demand undergo double of the current braking phases the battery use the maximum of his power to satisfy the motorization demand under the slopped road condition which can interpreted physically the augmentation of the globally vehicle resistive torque illustrate in Table 4. In the other hand the linear speeds of the two induction motors stay the same and the road drop does not influence the torque control of each wheels. The results are listed in Table 3.

Fig. 18. Variation of driving force of the right motor in different phases.


Table 4. Variation of vehicle torque in different Phases.

Fig. 19. Evaluation of the globally vehicle resistive torque compared to nominal motor torque in different phases.

The Uses of Artificial Intelligence for Electric Vehicle Control Applications 253

Mag (% of Fundamental)

(a) (b)

Referring to figures 20 (a) and (b) we show harmonic analyses of stator current. DTC-SVM with intelligent fuzzy PI controller present 29.96%, DTC-SVM with classical PI controller give 33.82 %. The first controller offer an reduction of 13.84% .This remarkable change obtained enables us to say that the current inject by voltage source inverter in DTC-SVM classical PI controller is harmonics current polluting what to justify the great oscillations of the torque and the attraction force .as a consequence this ripple present negative effects on the autonomy of the battery and heating of the both motors and

Designation PI controller intelligent fuzzy PI

THD [%] 33.82 29.96

Table 5. Comparative between PI and intelligent fuzzy PI controller

Comportment of electric vehicle in sloped road Less adaptive More adaptive

Driving forces and electromagnetic torque More oscillation Less oscillation

The research outlined in this paper has demonstrated the feasibility of an improved vehicle stability which utilizes two independent back drive wheels for motion by using DTC-SVM controls. DTC-SVM with intelligent Fuzzy control is able to adapt itself the suitable control parameters which are the proportional and integral gains kp and ki to the variations of vehicle torque. This method was Improved EV steering and stability during different trajectory. The advantage DTC-SVM controller is robustness and performance, there capacity to maintain ideal trajectories for two wheels control independently and ensure good disturbances rejections with no overshoot and stability of vehicle perfected ensured with the speed variation and less error speed. The DTC-SVM with intelligent fuzzy PI controller is more adaptive for propelled systems. The electric vehicle was proved best comportment and stability during different road path by

Fig. 22. (a) Total harmonic distortion using DTC with intelligent fuzzy PI controller, (b)

<sup>0</sup> <sup>5</sup> <sup>10</sup> <sup>15</sup> <sup>20</sup> <sup>25</sup> <sup>30</sup> <sup>35</sup> <sup>40</sup> <sup>45</sup> <sup>50</sup> <sup>0</sup>

Total Harmonic Distortion THD= 33.82%

Harmonic order

controller

<sup>0</sup> <sup>5</sup> <sup>10</sup> <sup>15</sup> <sup>20</sup> <sup>25</sup> <sup>30</sup> <sup>35</sup> <sup>40</sup> <sup>45</sup> <sup>50</sup> <sup>0</sup>

Total Harmonic Distortion THD= 29.96%

Harmonic order

Total harmonic distortion using DTC-SVM with classical PI.

increase power losses.

**8. Conclusions** 

Mag (% of Fundamental)

According to the formulas (2),(3),(4) and Table. 3, the variation of resistive vehicle torques in different cases as depicted in Table 4. , the vehicle resistive torque was 95.31 N.m in the first case (acceleration phase) when the power propulsion system resistive one is only 68.53 Nm in the breaking phases (phases 2) , the back driving wheels develop more and more efforts to satisfy the traction chain demand which impose an resistive torque equal to 168.00 N.m .The result prove that the traction chain under acceleration demand develop the double effort comparing with the breaking phase case's by means that the vehicle needs the half of its energy in the deceleration phase's compared with the acceleration one's as it specified in table 2.

#### **B. Comparative study of two method of controlling**

In simulations the two different methods to control the EV were used .Because of the sweeping of the kp on the interval [15 43 ] and the ki on the interval [228 243 ] as shown in the figure 20 and 21.The DTC with Fuzzy Adaptive PI Control method improves EV performance. The intelligent fuzzy PI controller was proved in efficiency adaptation for stability of the vehicle. The results obtained by simulation show that this structure permits the realization of the robust control based on intelligent fuzzy inference system, with good dynamic and static performances for the multi-converters/multi-machines propelled system.

Fig. 20. Variation gain kp of intelligent fuzzy PI controller

Fig. 21. Varaition gain ki of adaptive fuzzy PI controller

Fig. 22. (a) Total harmonic distortion using DTC with intelligent fuzzy PI controller, (b) Total harmonic distortion using DTC-SVM with classical PI.

Referring to figures 20 (a) and (b) we show harmonic analyses of stator current. DTC-SVM with intelligent fuzzy PI controller present 29.96%, DTC-SVM with classical PI controller give 33.82 %. The first controller offer an reduction of 13.84% .This remarkable change obtained enables us to say that the current inject by voltage source inverter in DTC-SVM classical PI controller is harmonics current polluting what to justify the great oscillations of the torque and the attraction force .as a consequence this ripple present negative effects on the autonomy of the battery and heating of the both motors and increase power losses.


Table 5. Comparative between PI and intelligent fuzzy PI controller

## **8. Conclusions**

252 Technology and Engineering Applications of Simulink

According to the formulas (2),(3),(4) and Table. 3, the variation of resistive vehicle torques in different cases as depicted in Table 4. , the vehicle resistive torque was 95.31 N.m in the first case (acceleration phase) when the power propulsion system resistive one is only 68.53 Nm in the breaking phases (phases 2) , the back driving wheels develop more and more efforts to satisfy the traction chain demand which impose an resistive torque equal to 168.00 N.m .The result prove that the traction chain under acceleration demand develop the double effort comparing with the breaking phase case's by means that the vehicle needs the half of its energy in the deceleration phase's compared with the acceleration one's as it

In simulations the two different methods to control the EV were used .Because of the sweeping of the kp on the interval [15 43 ] and the ki on the interval [228 243 ] as shown in the figure 20 and 21.The DTC with Fuzzy Adaptive PI Control method improves EV performance. The intelligent fuzzy PI controller was proved in efficiency adaptation for stability of the vehicle. The results obtained by simulation show that this structure permits the realization of the robust control based on intelligent fuzzy inference system, with good dynamic and static performances for the multi-converters/multi-machines propelled

0 1 2 3 4 5 6

**Phase 1 Phase 2 Phase 3**

**Time [s]**

**Phase 1 Phase 2 Phase 3**

0 1 2 3 4 5 6

**Time [s]**

specified in table 2.

system.

**B. Comparative study of two method of controlling** 

210

220

230

**Gain[ki]**

240

250

Fig. 20. Variation gain kp of intelligent fuzzy PI controller

Fig. 21. Varaition gain ki of adaptive fuzzy PI controller

**Gain[kp]**

The research outlined in this paper has demonstrated the feasibility of an improved vehicle stability which utilizes two independent back drive wheels for motion by using DTC-SVM controls. DTC-SVM with intelligent Fuzzy control is able to adapt itself the suitable control parameters which are the proportional and integral gains kp and ki to the variations of vehicle torque. This method was Improved EV steering and stability during different trajectory. The advantage DTC-SVM controller is robustness and performance, there capacity to maintain ideal trajectories for two wheels control independently and ensure good disturbances rejections with no overshoot and stability of vehicle perfected ensured with the speed variation and less error speed. The DTC-SVM with intelligent fuzzy PI controller is more adaptive for propelled systems. The electric vehicle was proved best comportment and stability during different road path by

The Uses of Artificial Intelligence for Electric Vehicle Control Applications 255

[8] Lin Chen j, Kang-Ling Fang chen . A Novel Direct Torque Control for Dual-Three-Phase

[9] P. Vas. Sensorless Vector and Direct Torque Control, Oxford University Press, 1998. [10] [10] A. Schell., H. Peng., D. Tran., E.Stamos. Modeling and Control Strategy

[11] A. Haddoun., M. Benbouzid., D. Diallo., R. Abdesseme.,, J. Ghouili., K. Srairi. Modeling

[12] A. Nasri., A .Hazzab .,I.K. Bousserhane., S. Hadjeri., P. Sicard., Two Wheel Speed

[13] K. Hartani., M.bourahla,. Y.miloud, M.sekour., Electronic Differential with Direct

[14] L.T. Lam., R. louey. Developpement of ultra-battery for hybrid-electric vehicle

[15] Larminie. Electric Vehicle Technology Explained, *Edited by John Wiley and John Lowry,* 

[16] A. Haddoun et al. Analysis, modeling and neural network traction control of an electric

[17] M. Vasudevan., R. Arumugam ., New direct torque control scheme of induction

[18] M. E. H. Benbouzid et al. Advanced fault-tolerant control of inductionmotor drives

[19] A. Gupta., A. M. Khambadkone. A space vector pwm scheme for multilevel inverters

[20] T. G. Habetler., F. Profumo., M.Pastorelli., L. Tolbert. Direct torque control of induction

[21] Mingqian Gao., Shanghong He. Self-adapting Fuzzy-PID Control of Variable Universe

[22] J.-Y. Chen, P.-S. Tsai and C.-C. Wong .Adaptive design of a fuzzy cerebellar model

*vol.* 28, no. 5, pp. 1045-1053, septembre/October 1992.

applications, *Elservier, power sources, vol*. 158, pp. 1140-1148, 2006.

*Cybernetics,* pp. 876-88, 2003.

vol.17, no.1, 2009, TUBITAK.

*industriel electronic vol*. 55,N 6 June 2008.

Electrical Engineering , vol. 5, no.2, pp. 199-216, 2008.

pp. 159-168, 2005.

*England*, 2003.

pp. 854–859.

1383, 2004.

53, October 2006.

*Technology and Automation*.

152, no.2, pp. 133-137, 2005.

2007.

Induction Motor, *Conf. Rec. IEEE International Conference on Machine Learning and* 

development for Fuel Cell Electric Vehicle, *Annual Review in control Elseiver, vol*. 29,

,Analysis and neural network control of an EV Electrical Differentiel, *Transaction on* 

Robust Sliding Mode Control For Electric Vehicle Drive, Serbian Journal of

Torque Fuzzy Control for Vehicle Propulsion System, Turk J Elec Eng & Comp Sci,

vehicle without differential gears, *in Proc. IEEE IEMDC, Antalya, Turkey, May* 2007,

motor for electric vehicles, 5th Asian Control Conference, vol. 2, pp. 1377 –

for EV/HEV traction applications, *From conventional to modern and intelligent control techniques, IEEE Trans. Veh. Technol., vol.* 56, no. 2, pp. 519–528, Mar.

based on two-level space vector pwm, *IEEE Transaction on Industrial Electronics, vol.*

machines using space vector modulation*, IEEE Transaction on Industry Applications,* 

in the Non-linear System, *2008 International Conference on Intelligent Computation* 

arithmetic controller neural network, *IEE Proc.-Control Theory and pplications, vol.*

maintaining the motorization error speed equal zeros and gives a good distribution for deriving forces. The electric vehicle was proved efficiency comportment in the different road constraints.


Table 6. Vehicle Parameters.

## **9. References**


maintaining the motorization error speed equal zeros and gives a good distribution for deriving forces. The electric vehicle was proved efficiency comportment in the different

Te Motor traction torque 238 Nm

Je Moment on inertia of the drive train 7.07Kgm2

Rw Wheel radius 0.32m

M Vehicle mass 1300Kg

fe Bearing friction coefficient 0.32

Kd Aerodynamic coefficient 0.32

A Vehicle frontal area 2.60 m2

fv Vehicle friction coefficient 0.01

Lw Distance between two wheels and axes 2.5m

dw Distance between the back and the front wheel 1.5m

[1] Yee-Pien Yang., Chun-Pin Lo. Current Distribution Control of Dual Directly Driven

[2] P. He., Y. Hori., M. Kamachi ., K. Walters., H. Yoshida. Future Motion Control to be

[3] Jun-Koo Kang ., Seung-Ki Sul. New Direct Torque Control of Induction Motor for

[4] C. C. Chan . Electric vehicles charge forward*, IEEE Power Energy Mag, vol*. 2, no. 6, pp.

[5] Z. Q. Zhu ., David HoweZ. Zhu . Electrical machines and drives for electric, hybrid, and

fuel cell vehicles, *Proc. IEEE, vol. 95, no. 4,* pp. 764–765, 2007.

[6] P. Vas. Sensorless Vector and Direct Torque Control, *Oxford University Press*, 1998. [7] K. Itoh ., H. Kubota. Thrust ripple reduction of linear induction motor with direct torque

Wheel Motors for Electric Vehicles, *Control Engineering Practice, vol.* 16, no. 11, pp.

Realized by In-wheel Motored Electric Vehicle, *In Proceedings of the 31st Annual Confer- ence of the IEE Industrial Electronics Society, IEEE Press, Raliegh South Carolina,*

Minimum Torque Ripple and Constant Switching Frequency, *IEEE Trans. Ind.* 

control, *Proceedings of the Eighth International Conference on Electrical Machines and* 

road constraints.

Table 6. Vehicle Parameters.

1285-1292, 2008.

24–33, 2004.

USA, pp. 2632-2637, 2005.

*Applicat., vol*. 35, no. 5, p. 1076-1082, 1999.

*Systems, ICEMS 2005, vol*. 1, pp. 655-658, 2005.

**9. References** 


[23] Chih-Min Lin., Ya-Fu Peng., Adaptive CMAC-Based Supervisory Control for Uncertain

[24] Hugang Han., Chun-Yi Su., Yury Stepanenko. Adaptive control of a class of nonlinear

cybernetics, vol. 34, no. 2, april 2004.

*fuzzy systems, vol.* 9, no. 2, april 2001.

Nonlinear Systems, IEEE Transactions on systems, man, and cybernetics-part b:

systems with nonlinearly parameterized fuzzy approximators, *IEEE transactions on* 

## *Edited by Subhas Chakravarty*

Building on MATLAB (the language of technical computing), Simulink provides a platform for engineers to plan, model, design, simulate, test and implement complex electromechanical, dynamic control, signal processing and communication systems. Simulink-Matlab combination is very useful for developing algorithms, GUI assisted creation of block diagrams and realisation of interactive simulation based designs. The eleven chapters of the book demonstrate the power and capabilities of Simulink to solve engineering problems with varied degree of complexity in the virtual environment.

Photo by NiPlot / iStock

Technology and Engineering Applications of Simulink

Technology and Engineering

Applications of Simulink

*Edited by Subhas Chakravarty*