As described in Ross (2004),S. Sumathi (2010) and Andrade (2010) the commitment of systems based on fuzzy logic both in respect of the facility as implementation of policies that take into account not only the performance of closed loop of industrial process control as well as the

systems are likely to become the largest presence within the control and automation

end

**4.3 Real time control**

**4.4 Implementing intelligent systems**

Fig. 8. Diagram of the General System

**4.4.1 A fuzzy method implementation**

computational effort. The diagram shown in Figure 8

*Ts* variable.

#### Fig. 7. Graphical Interface of the Platform System

```
function Button_01_Callback(hObject, eventdata, handles, varargin)
Serial_Obj=serial('COM1','BaudRate',9600);
opcon2 =get(handles.menu_01,'Value');
opcon = get(handles.menu_01,'String');
switch opcon{opcon2}
    case 'COM1'
        Serial_Obj.Port='COM1';
    case 'COM2'
        Serial_Obj.Port='COM2';
    case 'COM3'
        Serial_Obj.Port='COM3';
    case 'COM4'
       Serial_Obj.Port='COM4';
    case 'COM5'
      Serial_Obj.Port='COM5';
    case 'COM6'
       Serial_Obj.Port='COM6';
    case 'COM7'
      Serial_Obj.Port='COM7';
    otherwise
        errordlg('Select a Valid Serial Port','Error')
end
fopen(Serial_Obj);
if Serial_Obj.Status ==open
```

```
set(handles.con_text,'String','Connected');
end
```
## **4.3 Real time control**

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

function Button\_01\_Callback(hObject, eventdata, handles, varargin)

errordlg('Select a Valid Serial Port','Error')

Fig. 7. Graphical Interface of the Platform System

switch opcon{opcon2} case 'COM1'

case 'COM2'

case 'COM3'

case 'COM4'

case 'COM5'

case 'COM6'

case 'COM7'

otherwise

fopen(Serial\_Obj);

if Serial\_Obj.Status ==open

end

Serial\_Obj=serial('COM1','BaudRate',9600); opcon2 =get(handles.menu\_01,'Value'); opcon = get(handles.menu\_01,'String');

Serial\_Obj.Port='COM1';

Serial\_Obj.Port='COM2';

Serial\_Obj.Port='COM3';

Serial\_Obj.Port='COM4';

Serial\_Obj.Port='COM6';

Serial\_Obj.Port='COM5';

Serial\_Obj.Port='COM7';

Digital Control Systems has its roots in the interface with the analog world and thus the delay time this conversion to occur immediately. In real-time systems is a concern that the execution time of a given instruction does not exceed a predetermined threshold. Converging, the realtime systems are strictly necessary in digital control because it does necessary to guarantee instruction execution control *Uc* within the limit of sampling system *Ts* where a failure of this requirement may lead to instability as the plant design of discrete controllers is directly related *Ts* variable.

## **4.4 Implementing intelligent systems**

The implementation of intelligent systems through *MATLAB* may well become a very profitable it facilitates the testing of new techniques that use features and tools already implemented in this software. Techniques computational intelligence can withdraw from the digital control system the important factor in the *Ts*, however the system still requires real time control. Each execution cycle lets you have the need to perform numerical derivatives and integrals and introducing a context of researching and mapping , which may require less computational effort. The diagram shown in Figure 8

Fig. 8. Diagram of the General System

It´s possible conclude that with the development of both the hardware and software, these systems are likely to become the largest presence within the control and automation

### **4.4.1 A fuzzy method implementation**

As described in Ross (2004),S. Sumathi (2010) and Andrade (2010) the commitment of systems based on fuzzy logic both in respect of the facility as implementation of policies that take into account not only the performance of closed loop of industrial process control as well as the

Fig. 10. The FLT Main Screen

method chosen.

**5. Conclusion**

from a direct single platform.

of rules based on Table 1 decision-making as an example

Table 1. The knowledge from the Rules of Decision System

Where N*egative Large*, *Negative Small, ZeRo*, *Positive Small* and *Positive Large* are linguistics vary based on experience. The system output can be given as shown in Figure 9 with a bank

MATLAB as Interface for Intelligent Digital Control of Dynamic Systems 289

*Error* \ *changingError* **NL NS ZR PS PL**

And the controller output and results in a weighted sum that depends exclusively on the

With the aid of this important tool can facilitate the development of applications to interface with and control devices industries. A overview of these tools has been shown here to that may be developed to greater diversity of applications according to the need of the designer. As an environment of the computational complex, *MATLAB* can provide solutions for engineering, physics and mathematics with the use of its functions basic and therefore leaving the search more efficient and less costly. This platform is an important and necessary tool, it can lead a project, with integration between *Software* and *Hardware*, for a safe convergence of results, thus integrating computational resources for simulation design and implementation

**NL** PL PL PL PL NL **NS** PS PS PS PS NS **ZR** PS ZR ZR NS ZR **PS** NS NS NS NS PS **PL** NL NL NL NS PS

experience acquired human with the process that can be transferred directly to the core of the *Fuzzy Logic Controller*(*FLC*) with its implementation process described in Figure 9

#### Fig. 9. Diagram of the General System

The consolidation of this type of system could provide designers with a feature not found on controllers classical, ie the aforementioned expertise. To use the Fuzzy controller is needed for the project prior knowledge on the implementation of systems based on fuzzy logic and can be found in Andrade (2010) and especially in Zadeh (1965). In *MATLAB* the design begins by checking the feasibility of the controller and can be done by graphical interface *Fuzzy Logic Toolbox* with a significant help from MathWorks (2010).

To open the *FLT (Fuzzy Logic Toolbox)* typed:

>>fuzzy

#### and the following screen will open

However the system to be implemented with *MATLAB* using the serial communication is necessary to implement the intelligent system in *script* and can be done as follows

```
error = data_ctrl;
NG(k) = trimf(error, [a_i b_i c_i]);
NS(k) = trimf(error, [a_i b_i c_i]);
ZR(k) = trimf(error, [a_i b_i c_i]);
PS(k) = trimf(error, [a_i b_i c_i]);
PB(k) = trimf(error, [a_i b_i c_i]);
chang_error = data_ctrl;
NL(k) = trimf(chang_error, [a_i b_i c_i]);
NS(k) = trimf(chang_error, [a_i b_i c_i]);
ZR(k) = trimf(chang_error, [a_i b_i c_i]);
PS(k) = trimf(chang_error, [a_i b_i c_i]);
PL(k) = trimf(chang_error, [a_i b_i c_i]);
```

Fig. 10. The FLT Main Screen

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

experience acquired human with the process that can be transferred directly to the core of the

The consolidation of this type of system could provide designers with a feature not found on controllers classical, ie the aforementioned expertise. To use the Fuzzy controller is needed for the project prior knowledge on the implementation of systems based on fuzzy logic and can be found in Andrade (2010) and especially in Zadeh (1965). In *MATLAB* the design begins by checking the feasibility of the controller and can be done by graphical interface *Fuzzy Logic*

However the system to be implemented with *MATLAB* using the serial communication is

necessary to implement the intelligent system in *script* and can be done as follows

*Fuzzy Logic Controller*(*FLC*) with its implementation process described in Figure 9

Fig. 9. Diagram of the General System

*Toolbox* with a significant help from MathWorks (2010).

To open the *FLT (Fuzzy Logic Toolbox)* typed:

NG(k) = trimf(error, [a\_i b\_i c\_i]); NS(k) = trimf(error, [a\_i b\_i c\_i]); ZR(k) = trimf(error, [a\_i b\_i c\_i]); PS(k) = trimf(error, [a\_i b\_i c\_i]); PB(k) = trimf(error, [a\_i b\_i c\_i]);

NL(k) = trimf(chang\_error, [a\_i b\_i c\_i]); NS(k) = trimf(chang\_error, [a\_i b\_i c\_i]); ZR(k) = trimf(chang\_error, [a\_i b\_i c\_i]); PS(k) = trimf(chang\_error, [a\_i b\_i c\_i]); PL(k) = trimf(chang\_error, [a\_i b\_i c\_i]);

and the following screen will open

chang\_error = data\_ctrl;

error = data\_ctrl;

>>fuzzy

Where N*egative Large*, *Negative Small, ZeRo*, *Positive Small* and *Positive Large* are linguistics vary based on experience. The system output can be given as shown in Figure 9 with a bank of rules based on Table 1 decision-making as an example

Table 1. The knowledge from the Rules of Decision System

And the controller output and results in a weighted sum that depends exclusively on the method chosen.

## **5. Conclusion**

With the aid of this important tool can facilitate the development of applications to interface with and control devices industries. A overview of these tools has been shown here to that may be developed to greater diversity of applications according to the need of the designer. As an environment of the computational complex, *MATLAB* can provide solutions for engineering, physics and mathematics with the use of its functions basic and therefore leaving the search more efficient and less costly. This platform is an important and necessary tool, it can lead a project, with integration between *Software* and *Hardware*, for a safe convergence of results, thus integrating computational resources for simulation design and implementation from a direct single platform.

Angel L. Morales1, Jem A. Rongong2 and Neil D. Sims2

**Control Systems of Structures**

**MATLAB for the Analysis of Uncertain**

Vibration is an important aspect of many engineering systems. In most cases, such vibration is undesirable and requires attenuation or control, which explains the huge quantity of vibration control techniques that can be found in the literature. Nevertheless, sometimes all this knowledge is not enough to guarantee the success in this purpose. Model inaccuracies or parameter uncertainty are unavoidable in all relevant practical application and not only can they degrade the performance of an otherwise well-designed control system, but they may

**A Fuzzy Finite Element Method Programmed in**

For these reasons, robustness is of crucial importance in control-system design. Although always being appreciated, it was not until late 1970s when a theory capable to explicitly handle the robustness issue in feedback design was developed: the *H*∞ optimal control theory (Zames, 1981; Zames & Francis, 1983). This optimization approach is well developed and we can even find a full *Robust Control Toolbox* in MATLAB which, in collaboration with the *Control System Toolbox*, let us identify worst-case scenarios and generate optimal controllers so that

Obviously, the *H*∞ optimal control theory shows some disadvantages, like the high level of mathematical understanding needed to apply them successfully and the fact that it sometimes leads to potentially very conservative results (Gu et al., 2005). Indeed, even when a solution via *H*∞ techniques is eventually found, the process may have been too expensive in terms of designer's time and the sensitivity to the uncertainty sources is not analyzed. These motivations make us think that a straightforward possibilistic methodology for simulating and analyzing the uncertainty propagation in control system would be a very valuable,

This work focuses on how MATLAB can supply suitable tools in order to manage uncertainty propagation in finite element models of structures from the point of view of fuzzy arithmetic. More specifically, we will describe a methodology fully programmed in MATLAB in order to deal with uncertain control systems of vibrating structures and we will illustrate the

In order to develop this work we need to integrate two different problems: the calculation of the dynamic response of a vibrating structure via finite element models and the study of the

the stability and good performance of the plant is ensured (Gu et al., 2005).

**1. Introduction**

also lead to instability and even structural failure.

complementary and time saving tool.

application of this approach to a case study.

1*Universidad de Castilla - La Mancha*

2*The University of Sheffield*

1*Spain*

**15**

2*United Kingdom*

#### **6. References**

## **A Fuzzy Finite Element Method Programmed in MATLAB for the Analysis of Uncertain Control Systems of Structures**

Angel L. Morales1, Jem A. Rongong2 and Neil D. Sims2 1*Universidad de Castilla - La Mancha* 2*The University of Sheffield* 1*Spain* 2*United Kingdom*

#### **1. Introduction**

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

290 MATLAB – A Ubiquitous Tool for the Practical Engineer

Andrade, G.A. Neto, J. L. L. (2010). A framework for modeling, digital control design

Charles L. Philips, H. T. N. (1995). *Digital Control Systems: Analisys and Design*, third edition

Gene F. Franklin, J. David Powell, M. L. W. (1998). *Digital Control of Dynamic Systems*, 3rd edn,

*International Conference on*, IEEE Computer Society Digital Library.

Ibrahim, D. (2006). *Microcontroller Based Applied Digital Control*, John Wiley and Sons. Lewis, C. S. (2001). *Programming Microcontrollers in C*, 2nd edn, LLH Publishing.

Ross, T. J. (2004). *Fuzzy Logic With Engineering Applications*, 2nd edn, John Wiley and Sons. S. Sumathi, S. P. (2010). *Computational Intelligence Paradigms: Theory and Applications Using*

Wilmshurst, T. (2007). *Designing Embedded Systems with PIC Microcontrollers*, Newnes.

and simulations of dynamic, *Computer Modelling and Simulation (UKSim), 2010 12th*

**6. References**

edn, Prentice Hall.

*MATLAB*, CRC Press.

MathWorks, T. (2010). *Fuzzy Logic Tool Box : User's Guide*. Moudgalya, K. (2007). *Digital Control*, John Wiley and Sons.

Zadeh, L. A. (1965). Fuzzy sets, *Information and Control Proceedings* .

Prentice Hall.

Vibration is an important aspect of many engineering systems. In most cases, such vibration is undesirable and requires attenuation or control, which explains the huge quantity of vibration control techniques that can be found in the literature. Nevertheless, sometimes all this knowledge is not enough to guarantee the success in this purpose. Model inaccuracies or parameter uncertainty are unavoidable in all relevant practical application and not only can they degrade the performance of an otherwise well-designed control system, but they may also lead to instability and even structural failure.

For these reasons, robustness is of crucial importance in control-system design. Although always being appreciated, it was not until late 1970s when a theory capable to explicitly handle the robustness issue in feedback design was developed: the *H*∞ optimal control theory (Zames, 1981; Zames & Francis, 1983). This optimization approach is well developed and we can even find a full *Robust Control Toolbox* in MATLAB which, in collaboration with the *Control System Toolbox*, let us identify worst-case scenarios and generate optimal controllers so that the stability and good performance of the plant is ensured (Gu et al., 2005).

Obviously, the *H*∞ optimal control theory shows some disadvantages, like the high level of mathematical understanding needed to apply them successfully and the fact that it sometimes leads to potentially very conservative results (Gu et al., 2005). Indeed, even when a solution via *H*∞ techniques is eventually found, the process may have been too expensive in terms of designer's time and the sensitivity to the uncertainty sources is not analyzed. These motivations make us think that a straightforward possibilistic methodology for simulating and analyzing the uncertainty propagation in control system would be a very valuable, complementary and time saving tool.

This work focuses on how MATLAB can supply suitable tools in order to manage uncertainty propagation in finite element models of structures from the point of view of fuzzy arithmetic. More specifically, we will describe a methodology fully programmed in MATLAB in order to deal with uncertain control systems of vibrating structures and we will illustrate the application of this approach to a case study.

In order to develop this work we need to integrate two different problems: the calculation of the dynamic response of a vibrating structure via finite element models and the study of the

where *M*, *K* and *C* are the mass, the stiffness and the damping matrices, respectively, *x* the vector of displacements of all the degrees of freedom of the structure and *f* a general vector of forces. Matrices *M* and *K* arise from the finite element analysis of the structure whereas *C*

293 A Fuzzy Finite Element Method Programmed

Both in the state space approach and in the transfer function approach, working with physical coordinates is not practical because of the usually high number of degrees of freedom. Thus, a change of variables into modal coordinates is suggested. In addition, in this way the problem can be restricted to the bandwidth of interest, neglecting the high frequency dynamics of the system. After the computation of modes via finite element analysis one obtains the natural frequencies of the modes and the matrix of the mode shapes Φ. After performing the change of variables *x* = Φ*z*, where *z* is the vector of modal coordinates and Φ is normalized in such a way that the modal mass matrix is the identity matrix, the governing equations in modal

where Ω and *ξ* define the diagonal matrices of natural frequencies and damping factors. Finally, the software also provides several commands related to the definition of the actuators and sensors or, strictly speaking, the input force influence vector *Lu* and the sensor influence vector *Ly*, both of them required for the definition of either the transfer function or the state

The finite element models which are obtained with this toolbox can be used to solve both static and dynamic problems, the results of which being ready to be plotted or post-processed under any other MATLAB toolbox. In this work, we blended these commands which those

The *Control System Toolbox* provides a collection of matrix-based functions, usually expressed as M-files, for the topic of control engineering. The control systems may be modelled as transfer functions, zero-pole-gain or even using the state space approach, which allows the designer to apply classic techniques as well as modern techniques. However, it is not our purpose to study thoroughly this topic, so we will concentrate on briefly reviewing, as an example, one of the most common active control techniques of structures (Preumont, 2002): active damping with collocated piezoelectric pairs via Positive Position Feedback, which, in

The transfer function of the multi degree of freedom vibrating structure obtained via modal

where Φ*n*, *ωn* and *ξn* are the mode shape, the natural frequency and the damping ratio for the *n*th mode of vibration, respectively. The same model based upon the first *N* modes of

*ω*2

Φ*n*ΦT *n*

Φ*n*ΦT *n*

*n* + 2j*ωξnωn* − *ω*2 (3)

*n* + 2j*ωξnωn* − *ω*2 (4)

∞ ∑ *n*=1

*z*¨ + 2*ξ*Ω*z*˙ + Ω2*z* = ΦT *f* (2)

is usually obtained under the hypothesis of Rayleigh damping.

in MATLAB for the Analysis of Uncertain Control Systems of Structures

coordinates read

space matrices.

expansion reads

vibration leads to the expression

corresponding to the *Control System Toolbox*.

**2.2 Active damping and the Control System Toolbox**

the end, will be the case study discussed in Section 5.

*G*(*ω*) = *X*(*ω*)

*G*ˆ(*ω*) ≈

*F*(*ω*) =

*N* ∑ *n*=1

*ω*2

where the residual term due to the contributions of the truncated nodes has been neglected. Combining the equations for the surface-bonded actuator and sensor (Dosch et al., 1992; Sims et al., 2005) and Equation 4 which governs the structural dynamics, we obtain the

uncertainty propagation via fuzzy arithmetic. Both problems can be solved using MATLAB code obtaining very good results in terms of accuracy and computational cost.

The motivation behind this work is threefold: firstly, it is meant to show how MATLAB provides suitable functions and specialized toolboxes in order to study the static or dynamic response of a controlled structure via finite elements and control theory; secondly, it demonstrates the suitability of MATLAB to solve uncertainty propagation problems through fuzzy arithmetic minimizing the computational cost; and thirdly, it illustrates the application of a Fuzzy Finite Element procedure, which integrates these two tools, with a specific case study.

The chapter is organized in the following way: Section 1 has summarized the aims and concerns of the chapter; Section 2 deals with dynamic analysis of structures with active damping by means of the *OpenFEM Toolbox* and the *Control System Toolbox* of MATLAB; Section 3 discusses the different types of uncertainty and how they could be analyzed using fuzzy arithmetic programmed in MATLAB; Section 4 describes a Fuzzy Finite Element procedure in MATLAB in order to study the propagation of uncertain parameters in control systems of structures; Section 5 shows a case study in which the Fuzzy Finite Element procedure is tested; Section 6 deals with the main conclusions of this work; and finally the acknowledgements and references included through the text are included.

## **2. Dynamic analysis and control of structures**

#### **2.1 Finite Element analysis and the OpenFEM Toolbox**

The study of the dynamic behaviour of a structure can be carried out by means of finite element analysis. Among all the different possibilities to solve finite elements problems, in this chapter we focus on the use of a specific toolbox of MATLAB: the *OpenFEM Toolbox* (Balmes et al., 2009).

It is an open-source toolbox for finite element analysis within the matrix computing environment MATLAB. The software is developed in collaboration between Macs and the SDTools Company and it is distributed under a GNU Lesser General Public License, that is, a free software license which also allows the user to make modifications in the code. Performing finite element analysis within a matrix computing environment is of considerable interest, in particular as regards the ease of new developments, integration of external software, portability, post-processing, etc. This relatively young software is already quite successful in the finite element community as it is proven by the about 300 downloads per month (OpenFEM, 2011).

Though the toolbox architecture includes graphical user interfaces for visualization, analysis and animation of results, the powerful of this toolbox does not lie in this fact but in the standard and advanced methods which are implemented. Moreover, the provided functions are open and can be easily extended to suit particular needs of the user (Balmes et al., 2009), as it is our case.

Let us consider a vibrating structure. The first step in order to carry out its finite element analysis consists in discretizing the structure. The mesh nodes, elements which join the nodes, material properties, element properties and any additional information such as boundary conditions, loads, etc., are stored in a *struct* class variable which defines the model of the multi degree of freedom vibrating structure which is governed by the next set of equations of motion

$$\mathbf{M}\ddot{\mathbf{x}} + \mathbf{C}\dot{\mathbf{x}} + \mathbf{K}\mathbf{x} = f \tag{1}$$

2 Will-be-set-by-IN-TECH

uncertainty propagation via fuzzy arithmetic. Both problems can be solved using MATLAB

The motivation behind this work is threefold: firstly, it is meant to show how MATLAB provides suitable functions and specialized toolboxes in order to study the static or dynamic response of a controlled structure via finite elements and control theory; secondly, it demonstrates the suitability of MATLAB to solve uncertainty propagation problems through fuzzy arithmetic minimizing the computational cost; and thirdly, it illustrates the application of a Fuzzy Finite Element procedure, which integrates these two tools, with a specific case

The chapter is organized in the following way: Section 1 has summarized the aims and concerns of the chapter; Section 2 deals with dynamic analysis of structures with active damping by means of the *OpenFEM Toolbox* and the *Control System Toolbox* of MATLAB; Section 3 discusses the different types of uncertainty and how they could be analyzed using fuzzy arithmetic programmed in MATLAB; Section 4 describes a Fuzzy Finite Element procedure in MATLAB in order to study the propagation of uncertain parameters in control systems of structures; Section 5 shows a case study in which the Fuzzy Finite Element procedure is tested; Section 6 deals with the main conclusions of this work; and finally the acknowledgements and

The study of the dynamic behaviour of a structure can be carried out by means of finite element analysis. Among all the different possibilities to solve finite elements problems, in this chapter we focus on the use of a specific toolbox of MATLAB: the *OpenFEM Toolbox*

It is an open-source toolbox for finite element analysis within the matrix computing environment MATLAB. The software is developed in collaboration between Macs and the SDTools Company and it is distributed under a GNU Lesser General Public License, that is, a free software license which also allows the user to make modifications in the code. Performing finite element analysis within a matrix computing environment is of considerable interest, in particular as regards the ease of new developments, integration of external software, portability, post-processing, etc. This relatively young software is already quite successful in the finite element community as it is proven by the about 300 downloads per month

Though the toolbox architecture includes graphical user interfaces for visualization, analysis and animation of results, the powerful of this toolbox does not lie in this fact but in the standard and advanced methods which are implemented. Moreover, the provided functions are open and can be easily extended to suit particular needs of the user (Balmes et al., 2009),

Let us consider a vibrating structure. The first step in order to carry out its finite element analysis consists in discretizing the structure. The mesh nodes, elements which join the nodes, material properties, element properties and any additional information such as boundary conditions, loads, etc., are stored in a *struct* class variable which defines the model of the multi degree of freedom vibrating structure which is governed by the next set of equations of

*Mx*¨ + *Cx*˙ + *Kx* = *f* (1)

code obtaining very good results in terms of accuracy and computational cost.

references included through the text are included.

**2. Dynamic analysis and control of structures**

**2.1 Finite Element analysis and the OpenFEM Toolbox**

study.

(Balmes et al., 2009).

(OpenFEM, 2011).

as it is our case.

motion

where *M*, *K* and *C* are the mass, the stiffness and the damping matrices, respectively, *x* the vector of displacements of all the degrees of freedom of the structure and *f* a general vector of forces. Matrices *M* and *K* arise from the finite element analysis of the structure whereas *C* is usually obtained under the hypothesis of Rayleigh damping.

Both in the state space approach and in the transfer function approach, working with physical coordinates is not practical because of the usually high number of degrees of freedom. Thus, a change of variables into modal coordinates is suggested. In addition, in this way the problem can be restricted to the bandwidth of interest, neglecting the high frequency dynamics of the system. After the computation of modes via finite element analysis one obtains the natural frequencies of the modes and the matrix of the mode shapes Φ. After performing the change of variables *x* = Φ*z*, where *z* is the vector of modal coordinates and Φ is normalized in such a way that the modal mass matrix is the identity matrix, the governing equations in modal coordinates read

$$
\ddot{z} + 2\dddot{\xi}\Omega \dot{z} + \Omega^2 z = \Phi^\mathsf{T} f \tag{2}
$$

where Ω and *ξ* define the diagonal matrices of natural frequencies and damping factors.

Finally, the software also provides several commands related to the definition of the actuators and sensors or, strictly speaking, the input force influence vector *Lu* and the sensor influence vector *Ly*, both of them required for the definition of either the transfer function or the state space matrices.

The finite element models which are obtained with this toolbox can be used to solve both static and dynamic problems, the results of which being ready to be plotted or post-processed under any other MATLAB toolbox. In this work, we blended these commands which those corresponding to the *Control System Toolbox*.

#### **2.2 Active damping and the Control System Toolbox**

The *Control System Toolbox* provides a collection of matrix-based functions, usually expressed as M-files, for the topic of control engineering. The control systems may be modelled as transfer functions, zero-pole-gain or even using the state space approach, which allows the designer to apply classic techniques as well as modern techniques. However, it is not our purpose to study thoroughly this topic, so we will concentrate on briefly reviewing, as an example, one of the most common active control techniques of structures (Preumont, 2002): active damping with collocated piezoelectric pairs via Positive Position Feedback, which, in the end, will be the case study discussed in Section 5.

The transfer function of the multi degree of freedom vibrating structure obtained via modal expansion reads

$$\mathcal{G}(\omega) = \frac{\mathcal{X}(\omega)}{F(\omega)} = \sum\_{n=1}^{\infty} \frac{\Phi\_n \Phi\_n^{\mathrm{T}}}{\omega\_n^2 + 2 \mathrm{j}\omega\_{\mathrm{N}}^{\mathrm{T}} \omega\_{\mathrm{N}} - \omega^2} \tag{3}$$

where Φ*n*, *ωn* and *ξn* are the mode shape, the natural frequency and the damping ratio for the *n*th mode of vibration, respectively. The same model based upon the first *N* modes of vibration leads to the expression

$$\hat{\mathbf{G}}(\omega) \approx \sum\_{n=1}^{N} \frac{\Phi\_n \Phi\_n^\mathrm{T}}{\omega\_n^2 + 2 \mathbf{j}\omega \xi\_n \omega\_n - \omega^2} \tag{4}$$

where the residual term due to the contributions of the truncated nodes has been neglected. Combining the equations for the surface-bonded actuator and sensor (Dosch et al., 1992; Sims et al., 2005) and Equation 4 which governs the structural dynamics, we obtain the

On the one hand, we will use the term aleatory uncertainty to describe the inherent variation associated with the physical system or the environment under consideration. These sources of uncertainty can be represented as randomly distributed quantities, usually in terms of probability distribution functions. Typical examples of this kind of uncertainty are manufacturing tolerances, environmental effects, properties of non-uniform materials or any other identifiable disturbances. On the other hand, epistemic uncertainty derives from some level of ignorance of the system or the environment in any activity of the modelling process. In this case, its definition stresses the lack of knowledge as the origin of this uncertainty. Thus, unlike aleatory uncertainty, a mathematical representation of epistemic uncertainty has proven to be much more of a challenge since it is subjective and based on some expert opinion to a certain extent. Typical examples of epistemic uncertainties are some models for boundary conditions, simplified models for joints, models for material damping or unpredictable model

295 A Fuzzy Finite Element Method Programmed

in MATLAB for the Analysis of Uncertain Control Systems of Structures

Uncertainty propagation in dynamic systems has been a very frequent topic in the literature, which has been tackled from two different points of view (De Gersem et al., 2005): *probabilistic* and *possibilistic* approaches. Probabilistic approaches, such as the well-known Monte Carlo Simulation method, investigate the influence of aleatory uncertainties, this method being particularly suitable for finite element models with certain uncertainties. Possibilistic approaches are complementary to probabilistic approaches, since not only are they able to describe problems of aleatory uncertainties to a certain extent, yet without statistical interpretation, but they can also handle problems where incomplete information (epistemic uncertainty) is available. We consider that computationally expensive "brute-force" probabilistic methods are not worthwhile in most cases, where an intuitive worst-case scenario provides enough information. Thus, a possibilistic method in terms of fuzzy arithmetic will be the tool that we develop and propose in order to simulate and analyze

Fuzzy arithmetic has been used by various researchers as a technique for propagating uncertainty or variability through complex engineering models. The origins of this approach can be found in the theory of Fuzzy Sets (Zadeh, 1965). In contrast to the classical set theory, the elements of a fuzzy set are assigned a degree of membership to the set, which is referred to as the membership level *μ*. The core of the set is defined as the subset of elements for which *μ* = 1, whilst the support is the subset for which *μ* > 0. A fuzzy number is a fuzzy set that is

The issue of performing mathematical operations on fuzzy numbers instead of on conventional numbers can be a challenge rather than a trivial problem. In fact, the Fuzzy Arithmetic processes are usually performed by decomposing the fuzzy numbers into a

> *μi* = *i m*

Then, a interval analysis via interval arithmetic can be carried out separately to each

For monotonic problems this process may be trivial, since the maxima or minima of the inputs will lead to the maxima or minima for the output. Conversely, for non-monotonic problems this is no longer the case (Sims et al., 2010). In addition, the Standard Fuzzy Arithmetic may become problematic and lead to different results for the same problem depending on the form in which the solution procedure is applied, as can be seen in the examples gathered in (Hanss, 2002). In short, the application of Standard Fuzzy Arithmetic usually leads to overestimate

(8)

convex and normal, and whose membership function is piecewise continuous.

number of intervals given by the *α*-cuts at the *α*-levels *μi*, (*i* = 0, 1, . . . , *m*) with

results to a lesser or greater extent and alternative methods must be employed.

changes due to ageing, loading, etcetera.

uncertainty propagation in control systems.

membership level *μi*.

frequency response function between the voltage *Va* applied to the actuator and the voltage *Vs* at the output of the charge amplifier of the sensor

$$H(\omega) = \frac{V\_s(\omega)}{V\_a(\omega)} = k\_d k\_s \sum\_{n=1}^{N} \frac{L\_n^{\text{T}} \Phi\_n \Phi\_n^{\text{T}} L\_y}{\omega\_n^2 + 2 \text{j}\omega \xi\_n \omega\_n - \omega^2} \tag{5}$$

where *ks* and *ka* are the sensor gain and the actuator gain, respectively.

As far as the Positive Position Feedback scheme is concerned, it is appropriate for a structure equipped with strain actuators and sensors (Preumont, 2002). The objective of a Positive Position Feedback controller is to use a second order filter to improve the roll-off of the control system, allowing high frequency gain stabilization. The implementation of this control system is given by the following equations

$$\begin{array}{l} M\ddot{\mathbf{x}} + \mathbf{C}\dot{\mathbf{x}} + \mathbf{K}\mathbf{x} = L\_{\mu}\mu\\ y = L\_{\mu}^{T}\mathbf{x} \\ \ddot{v} + \beta\_{f}\dot{v} + \Omega\_{f}^{2}v = y \\ u = -Gv \end{array} \tag{6}$$

where *u* is the control force acting on the structure through the influence vector *Lu*, *y* is the difference of slope between the ends of the sensor, *G* = diag(*gi*) is the positive gain matrix, and *v* is the output of the second order filter defined by the matrices *β f* = diag(2*ξ f ωf*) and Ω2 *f* = diag(*ω*2 *f*). Note that the fact that the same localization vector (*Lu*) rules the actuator and sensors is a consequence of collocation.

The transfer function of the controller (a second order low pass filter) can be written as follows

$$gD(w) = g\frac{\omega\_f^2}{\omega\_f^2 + 2\mathrm{j}\omega\_f^\tau\omega\_f - \omega^2} \tag{7}$$

This compensator has both advantages and drawbacks. On the one hand, the open-loop transfer function has a roll-off of −40dB/decade, larger in comparison to Direct Velocity Feedback (−20dB/decade), which may decrease the risk of destabilizing high frequency dynamics. On the other hand, there is a stability limit which is reached when the open-loop static gain is equal to 1, this fact being independent of the damping (Preumont, 2002).

Once the open loop transfer functions of the plant and the regulator are defined, they can be used in the *Control System Toolbox* as *transfer function* class variables. Then, the desired active damping can be obtained by selecting the proper gain *g* of the regulator. In order to do this, several techniques can be used such as Bode Diagram Design, Root Locus Design, Nichols Plot design, etc.

Finally, the *Control System Toolbox* also provides valuable functions in order to analyze the stability and performance in terms of stability margins, time responses or the position of the closed loop poles and zeros in the pole-zero map. They will be used in our Fuzzy Finite Element methodology as fuzzy outputs in order to analyze the behaviour of the uncertain control system.

#### **3. Uncertainty and fuzzy arithmetic**

In the literature, the use of the term uncertainty is slightly ambiguous. In order to establish the nomenclature for the ensuing description, the reader is referred to the terminology proposed in (Oberkampf et al., 2004).

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

frequency response function between the voltage *Va* applied to the actuator and the voltage

*N* ∑ *n*=1

As far as the Positive Position Feedback scheme is concerned, it is appropriate for a structure equipped with strain actuators and sensors (Preumont, 2002). The objective of a Positive Position Feedback controller is to use a second order filter to improve the roll-off of the control system, allowing high frequency gain stabilization. The implementation of this control system

> *Mx*¨ + *Cx*˙ + *Kx* = *Luu y* = *LT u x*

> > *u* = −*Gv*

where *u* is the control force acting on the structure through the influence vector *Lu*, *y* is the difference of slope between the ends of the sensor, *G* = diag(*gi*) is the positive gain matrix, and *v* is the output of the second order filter defined by the matrices *β f* = diag(2*ξ f ωf*) and

The transfer function of the controller (a second order low pass filter) can be written as follows

This compensator has both advantages and drawbacks. On the one hand, the open-loop transfer function has a roll-off of −40dB/decade, larger in comparison to Direct Velocity Feedback (−20dB/decade), which may decrease the risk of destabilizing high frequency dynamics. On the other hand, there is a stability limit which is reached when the open-loop static gain is equal to 1, this fact being independent of the damping (Preumont, 2002). Once the open loop transfer functions of the plant and the regulator are defined, they can be used in the *Control System Toolbox* as *transfer function* class variables. Then, the desired active damping can be obtained by selecting the proper gain *g* of the regulator. In order to do this, several techniques can be used such as Bode Diagram Design, Root Locus Design, Nichols

Finally, the *Control System Toolbox* also provides valuable functions in order to analyze the stability and performance in terms of stability margins, time responses or the position of the closed loop poles and zeros in the pole-zero map. They will be used in our Fuzzy Finite Element methodology as fuzzy outputs in order to analyze the behaviour of the uncertain

In the literature, the use of the term uncertainty is slightly ambiguous. In order to establish the nomenclature for the ensuing description, the reader is referred to the terminology proposed

*ω*2

*gD*(*w*) = *g*

*v*¨ + *β f v*˙ + Ω2

*ω*2

*f v* = *y*

*f*). Note that the fact that the same localization vector (*Lu*) rules the actuator and

*ω*2 *f*

*L*T *u*Φ*n*ΦT *n Ly*

*n* + 2j*ωξnωn* − *ω*2 (5)

*f* + 2j*ωξ f ωf* − *ω*2 (7)

(6)

*Vs* at the output of the charge amplifier of the sensor

is given by the following equations

sensors is a consequence of collocation.

**3. Uncertainty and fuzzy arithmetic**

in (Oberkampf et al., 2004).

Ω2

*f* = diag(*ω*2

Plot design, etc.

control system.

*H*(*ω*) = *Vs*(*ω*)

*Va*(*ω*) = *kaks*

where *ks* and *ka* are the sensor gain and the actuator gain, respectively.

On the one hand, we will use the term aleatory uncertainty to describe the inherent variation associated with the physical system or the environment under consideration. These sources of uncertainty can be represented as randomly distributed quantities, usually in terms of probability distribution functions. Typical examples of this kind of uncertainty are manufacturing tolerances, environmental effects, properties of non-uniform materials or any other identifiable disturbances. On the other hand, epistemic uncertainty derives from some level of ignorance of the system or the environment in any activity of the modelling process. In this case, its definition stresses the lack of knowledge as the origin of this uncertainty. Thus, unlike aleatory uncertainty, a mathematical representation of epistemic uncertainty has proven to be much more of a challenge since it is subjective and based on some expert opinion to a certain extent. Typical examples of epistemic uncertainties are some models for boundary conditions, simplified models for joints, models for material damping or unpredictable model changes due to ageing, loading, etcetera.

Uncertainty propagation in dynamic systems has been a very frequent topic in the literature, which has been tackled from two different points of view (De Gersem et al., 2005): *probabilistic* and *possibilistic* approaches. Probabilistic approaches, such as the well-known Monte Carlo Simulation method, investigate the influence of aleatory uncertainties, this method being particularly suitable for finite element models with certain uncertainties. Possibilistic approaches are complementary to probabilistic approaches, since not only are they able to describe problems of aleatory uncertainties to a certain extent, yet without statistical interpretation, but they can also handle problems where incomplete information (epistemic uncertainty) is available. We consider that computationally expensive "brute-force" probabilistic methods are not worthwhile in most cases, where an intuitive worst-case scenario provides enough information. Thus, a possibilistic method in terms of fuzzy arithmetic will be the tool that we develop and propose in order to simulate and analyze uncertainty propagation in control systems.

Fuzzy arithmetic has been used by various researchers as a technique for propagating uncertainty or variability through complex engineering models. The origins of this approach can be found in the theory of Fuzzy Sets (Zadeh, 1965). In contrast to the classical set theory, the elements of a fuzzy set are assigned a degree of membership to the set, which is referred to as the membership level *μ*. The core of the set is defined as the subset of elements for which *μ* = 1, whilst the support is the subset for which *μ* > 0. A fuzzy number is a fuzzy set that is convex and normal, and whose membership function is piecewise continuous.

The issue of performing mathematical operations on fuzzy numbers instead of on conventional numbers can be a challenge rather than a trivial problem. In fact, the Fuzzy Arithmetic processes are usually performed by decomposing the fuzzy numbers into a number of intervals given by the *α*-cuts at the *α*-levels *μi*, (*i* = 0, 1, . . . , *m*) with

$$
\mu\_i = \frac{\dot{i}}{m} \tag{8}
$$

Then, a interval analysis via interval arithmetic can be carried out separately to each membership level *μi*.

For monotonic problems this process may be trivial, since the maxima or minima of the inputs will lead to the maxima or minima for the output. Conversely, for non-monotonic problems this is no longer the case (Sims et al., 2010). In addition, the Standard Fuzzy Arithmetic may become problematic and lead to different results for the same problem depending on the form in which the solution procedure is applied, as can be seen in the examples gathered in (Hanss, 2002). In short, the application of Standard Fuzzy Arithmetic usually leads to overestimate results to a lesser or greater extent and alternative methods must be employed.

Fuzzy inputs Fuzzy outputs

Iterative FE calculations

Fig. 1. Scheme of the numerical procedure to perform a Fuzzy Finite Element analysis using

variables. Note that if *μ* = 1, the interval analysis degenerates into the deterministic analysis. Finally, the membership functions of the required output variables are assembled from the

The Fuzzy Finite Element analysis is based on the results of a previous deterministic design of the control system. The control engineer must set the control parameters and gains of the regulator in order to accomplish all the features required for the active damping in terms of both stability and performance. After this process, which is usually easier and faster than other complex control techniques such as those derived from *H*∞ control, the Fuzzy Finite Element analysis will simulate the behaviour of the control system when some parameters are considered uncertain. This method is not a robust designing technique but a possibilistic method which tests whether the deterministic design is robust enough to keep the closed loop behaviour of the structure inside the allowable bounds for its stability and performance

The set of calculations which are carried out at each permutation of input values is compound

*Step 3:* Analysis of stability and performance via fuzzy variables such as stability margins or

The method chosen for carrying out the interval fuzzy arithmetics was the efficient implementation of the Transformation Method proposed in (Klimke, 2003). Nevertheless, we

The first modification lies in the fact that Klimke's MATLAB script is made under the assumption that the fuzzy output is the result of evaluating an analytical expression. In our case, we need to replace the evaluation of such analytical function (*inline* class) with a list of MATLAB commands which develops the aforementioned steps. This also allows us to include the possibility of achieving several fuzzy outputs instead of only one. The second modification was not necessary but it endows the method with enhance capabilities. It consists on using *cell* class arrays instead of *double* class arrays to define the fuzzy outputs, so an extra field may

*Step 2:* Closed loop simulation of the control system via the *Control System Toolbox*

have included some modifications which will be introduced in the ensuing.

*Step 1:* Finite Element calculations via the *OpenFEM Toolbox*

Interval analysis at α1

Deterministic analysis at α N

297 A Fuzzy Finite Element Method Programmed

in MATLAB for the Analysis of Uncertain Control Systems of Structures

Interval analysis at α N-1

...

Interval analysis at αn

Interval analysis at α2

...

α1 α2

y1

αn

μ( ) y1

. . .

. . .

αN-1 α N

α1 α2 αn

αN-1 α N

> ...

...

μ( ) x 2

α1 α2 αn

four *α*-cuts.

features.

of these steps:

active damping.

x 1

x 2

interval results at all sublevels.

αN-1 α N

> ...

μ( ) x 1

...

The Transformation Method, introduced in (Hanss, 2002), represents a special implementation of Fuzzy Arithmetic that avoids the undesired overestimation which usually arises when Fuzzy Arithmetic is reduced to interval computation. A later efficient implementation of this method was carried out in (Klimke, 2003), taking advantage of the matrix computing environment of MATLAB and providing enhanced features such as fast processing of discretized fuzzy numbers through multi-dimensional arrays, elimination of recurring permutations, automatic decomposition of models, treatment of single occurrences of variables through interval arithmetic and a monotonicity test based on automatic differentiation. Nevertheless, we have to note that the fuzzy parameters are restricted to convex fuzzy numbers when using the Transformation Method.

Specially relevant for an efficient implementation of the Transformation Method is the elimination of recurring permutations. The decomposition scheme of this method sometimes produces recurring points which in the end generate recurring permutations. By removing these recurring permutations from the evaluation procedure one could save a very valuable computation time. Indeed, we could even try to reuse as many points for different *α*-cuts as possible by selecting only inner points which have already occurred in a higher-level *α*-cut. For symmetric triangular membership functions, the obtained discretization is identical to the original formulation, but for other arbitrarily shaped membership functions, the distribution of the points is less regular, but of similar density (see Figure 1). However, for the same number of *α*-cuts, the less regular distribution of the inner points results in less accurate results compared to the original formulation.

The number of permutations *r* by using this efficient implementation of the Transformation Method is given by the following equation

$$r = m^{\text{\textquotedblleft}} + (m+1)^{\text{\textquotedblright}} \tag{9}$$

where *m* + 1 is the number of *α*-cuts and *n* is the number of fuzzy parameters. As proven in (Klimke, 2003), the new scheme provides a valuable improvement in terms of computational cost unless the number of *α*-cuts is too large compared to the number of fuzzy inputs. If that was the case, one could disregard the removing of recurring permutations in favour of more accurate results.

#### **4. Fuzzy Finite Element procedure in MATLAB**

The fuzzy theory has led to the development of Fuzzy Finite Element methods as that described in (De Gersem et al., 2005; De Munck et al., 2008; Moens & Vandepitte, 2005) for analyzing the propagation of uncertainties in the Frequency Response Functions of a dynamic system. When using this methodology, the uncertainties in the model parameters can be taken into account by fuzzy numbers with their shape obtained from experimental data or assumed from expert knowledge. Then, the numerical procedure for the implementation of these Fuzzy Finite Element methods consists of a sequence of Interval Finite Element procedures based on a general strategy for computing the fuzzy result from operations on fuzzy operands, i.e. the *α*-sublevel or *α*-cuts technique.

Figure 1 clarifies this procedure. This technique subdivides the membership range into a number of *α*-cuts. At each level, the intersection with the membership function of the input uncertainties results in an interval. Based on these input intervals for all uncertain parameters, an interval analysis can be performed obtaining the interval of the output variables at the considered *α*-cut. Strictly speaking, the interval analysis for a specific *α*-cut consists in running all the finite elements calculations corresponding to all the possible permutations among input variables at said level and choosing the lower and upper bound of the required output 6 Will-be-set-by-IN-TECH

The Transformation Method, introduced in (Hanss, 2002), represents a special implementation of Fuzzy Arithmetic that avoids the undesired overestimation which usually arises when Fuzzy Arithmetic is reduced to interval computation. A later efficient implementation of this method was carried out in (Klimke, 2003), taking advantage of the matrix computing environment of MATLAB and providing enhanced features such as fast processing of discretized fuzzy numbers through multi-dimensional arrays, elimination of recurring permutations, automatic decomposition of models, treatment of single occurrences of variables through interval arithmetic and a monotonicity test based on automatic differentiation. Nevertheless, we have to note that the fuzzy parameters are restricted to

Specially relevant for an efficient implementation of the Transformation Method is the elimination of recurring permutations. The decomposition scheme of this method sometimes produces recurring points which in the end generate recurring permutations. By removing these recurring permutations from the evaluation procedure one could save a very valuable computation time. Indeed, we could even try to reuse as many points for different *α*-cuts as possible by selecting only inner points which have already occurred in a higher-level *α*-cut. For symmetric triangular membership functions, the obtained discretization is identical to the original formulation, but for other arbitrarily shaped membership functions, the distribution of the points is less regular, but of similar density (see Figure 1). However, for the same number of *α*-cuts, the less regular distribution of the inner points results in less accurate results

The number of permutations *r* by using this efficient implementation of the Transformation

where *m* + 1 is the number of *α*-cuts and *n* is the number of fuzzy parameters. As proven in (Klimke, 2003), the new scheme provides a valuable improvement in terms of computational cost unless the number of *α*-cuts is too large compared to the number of fuzzy inputs. If that was the case, one could disregard the removing of recurring permutations in favour of more

The fuzzy theory has led to the development of Fuzzy Finite Element methods as that described in (De Gersem et al., 2005; De Munck et al., 2008; Moens & Vandepitte, 2005) for analyzing the propagation of uncertainties in the Frequency Response Functions of a dynamic system. When using this methodology, the uncertainties in the model parameters can be taken into account by fuzzy numbers with their shape obtained from experimental data or assumed from expert knowledge. Then, the numerical procedure for the implementation of these Fuzzy Finite Element methods consists of a sequence of Interval Finite Element procedures based on a general strategy for computing the fuzzy result from operations on fuzzy operands, i.e. the

Figure 1 clarifies this procedure. This technique subdivides the membership range into a number of *α*-cuts. At each level, the intersection with the membership function of the input uncertainties results in an interval. Based on these input intervals for all uncertain parameters, an interval analysis can be performed obtaining the interval of the output variables at the considered *α*-cut. Strictly speaking, the interval analysis for a specific *α*-cut consists in running all the finite elements calculations corresponding to all the possible permutations among input variables at said level and choosing the lower and upper bound of the required output

*r* = *mn* + (*m* + 1)*n* (9)

convex fuzzy numbers when using the Transformation Method.

compared to the original formulation.

accurate results.

*α*-sublevel or *α*-cuts technique.

Method is given by the following equation

**4. Fuzzy Finite Element procedure in MATLAB**

Fig. 1. Scheme of the numerical procedure to perform a Fuzzy Finite Element analysis using four *α*-cuts.

variables. Note that if *μ* = 1, the interval analysis degenerates into the deterministic analysis. Finally, the membership functions of the required output variables are assembled from the interval results at all sublevels.

The Fuzzy Finite Element analysis is based on the results of a previous deterministic design of the control system. The control engineer must set the control parameters and gains of the regulator in order to accomplish all the features required for the active damping in terms of both stability and performance. After this process, which is usually easier and faster than other complex control techniques such as those derived from *H*∞ control, the Fuzzy Finite Element analysis will simulate the behaviour of the control system when some parameters are considered uncertain. This method is not a robust designing technique but a possibilistic method which tests whether the deterministic design is robust enough to keep the closed loop behaviour of the structure inside the allowable bounds for its stability and performance features.

The set of calculations which are carried out at each permutation of input values is compound of these steps:

The method chosen for carrying out the interval fuzzy arithmetics was the efficient implementation of the Transformation Method proposed in (Klimke, 2003). Nevertheless, we have included some modifications which will be introduced in the ensuing.

The first modification lies in the fact that Klimke's MATLAB script is made under the assumption that the fuzzy output is the result of evaluating an analytical expression. In our case, we need to replace the evaluation of such analytical function (*inline* class) with a list of MATLAB commands which develops the aforementioned steps. This also allows us to include the possibility of achieving several fuzzy outputs instead of only one. The second modification was not necessary but it endows the method with enhance capabilities. It consists on using *cell* class arrays instead of *double* class arrays to define the fuzzy outputs, so an extra field may

*p*˜1 = �2.5, 2.5, 2.5�TFN and *p*˜2 = �3, 2, 2�TFN, in accordance to the nomenclature proposed in (Dubois & Prade, 1980). The fuzzy output obtained with the modified version clearly matched the results obtained with the original version. Nevertheless, an unavoidable loss

299 A Fuzzy Finite Element Method Programmed

Our modified version results more time-consuming because of the different way in which the simulated function is handled, as it can be seen in Table 1. In the original version, the *inline* class analytical expression can be executed only once with the array containing all the required permutations of fuzzy inputs at a time. In the modified version, the expression is evaluated separately at each permutation as it would be the case when a control system is simulated with several MATLAB commands. The elapsed time during these operations is called *tf cn*. As far as the obtention of the fuzzy output as a *cell* array instead of as a *double* array, this fact does not add any noticeable increase of computational cost but it does lead to heavier variables and higher memory requirements. The elapsed time during these operations is called *tmax*. In all the cases, the loss of efficiency is acceptable with regard to the benefits they provide, specially considering that the elapsed time in any finite element simulation may be on the order of

Consider an aluminium plate in free conditions, the dimensions and material properties of which are detailed in Table 2. The finite element model of this plate was carried out by means of the *OpenFEM Toolbox* of MATLAB (Balmes et al., 2009). After defining the mesh, the matrices of mass (*M*) and stiffness (*K*) which define the undamped model were obtained. The eigenvalue problem can then be solved, the natural frequencies and mode shapes being shown in Figure 4. The damping matrix (*C*) is constructed considering a structural damping

The piezoelectric actuator and sensor, the properties of which are detailed in Table 3, are collocated in the centre of the plate. As described in (Dosch et al., 1992; Sims et al., 2005), the piezoelectric actuator and sensor may be modelled as a beam which provides bending moments at its ends and measures a voltage proportional to the difference of slope of its ends. Due to collocation, the localization vectors of actuator (*Lu*) and sensor (*Ly*) are identical and

> *ttot* (s) 0.047 0.344 *ttot* (%) 100.00 100.00 *tf cn* (s) 0.016 0.313 *tf cn* (%) 34.04 90.99 *tmax* (s) 0.000 0.000 *tmax* (%) 0.000 0.000 *toth* (s) 0.031 0.031 *toth* (%) 65.96 9.01 Weight (bytes) 336 8568

implementation. The total time (*ttot*), the time for evaluating the analytical expression (*tf cn*), the time for calculating the maxima and minima of the output fuzzy number (*tmax*) and the time for other calculations (*toth*) are shown. The weight means the size in bytes of the output

Table 1. Efficiency comparison between the original and the modified Klimke's

Original Modified

*y* = (0, . . . , −1, 0, . . . , 1, 0, . . . , 0) (10)

seconds and it has much more restrictive memory requirements.

in MATLAB for the Analysis of Uncertain Control Systems of Structures

**5. A case study of uncertainty propagation**

*L*T *u* = *L*T

of *ξ*i = 0.6% in all the modes.

equal to

fuzzy variable.

of computational efficiency arose.

be included where information about each iteration is stored (for instance, the exact value of the fuzzy inputs at each permutation). Figures 2 and 3 show the MATLAB implementation of the efficient Transformation Method after modifying the Klimke's implementation. Note that the MATLAB implementation of the external function in Figure 3 only shows an example corresponding to the evaluation of a test function. In our case, this function would contain all the required commands which carry out the three foregoing calculating steps.

In order to prove the validity of our modified version of the Klimke's Transformation Method, a comparative fuzzy analysis was run. We took as a reference one example found in (Klimke, 2003) where the test function *f*(*x*1, *x*2) = cos(*πx*1)*x*2 is subjected to the triangular fuzzy inputs

```
function [fzr]=gtrmrecur_modif(varargin)
% Angel L. Morales, 2011.
% FZR=GTRMRECUR_MODIF(FZ1,FZ2,...,FZN) calls the general transformation
% method implemented by Klimke in 2003 with extended capabilities.
% The function allows N fuzzy numbers in alpha-cut representation.
% An accompanying function called GTRMRECUR_EXTERNAL_FCN is required.
fzin=varargin;
n = nargin;
m = size(fzin{1},1) - 1;
indices.subs = {};
indices.type = '()';
fzr = cell(m+1,2);
for k = 1:2
t = (floor((m+2-k)/2)+rem(m+2-k,2))*2-2+k;
for i = 1:n
c = reshape(fzin{i}(end+1-k:-2:1,:)',t+2-k,1);
c = c(3-k:end);
repvec = t*ones(1,n);
repvec(i) = 1;
x{i} = repmat(shiftdim(c,1-i),repvec);
end
s = k;
z_cell_matrix=gtrmrecur_external_fcn(x);
for j = m+2-k:-2:1
indices.subs(1:n) = {1:s};
w = reshape(subsref(z_cell_matrix,indices),s^n,1);
wc = zeros(size(w));
for cont=1:length(wc) wc(cont)=w{cont}{2}; end;
[none,pos_min]=min(wc); fzr{j,1} = w{pos_min};
[none,pos_max]=max(wc); fzr{j,2} = w{pos_max};
s = s + 2;
end
end
for j = m:-1:1
if fzr{j+1,1}{2}>=fzr{j,1}{2}; fzr{j,1}=fzr{j,1}; else fzr{j,1}=fzr{j+1,1}; end
if fzr{j+1,2}{2}>=fzr{j,2}{2}; fzr{j,2}=fzr{j+1,2}; else fzr{j,2}=fzr{j,2}; end
end
```
#### Fig. 2. MATLAB implementation of the efficient Transformation Method after modifying Klimke's script.

```
function [z_cell_matrix]=gtrmrecur_external_fcn(x)
% Angel L. Morales, 2011.
% Z_CELL_MATRIX=GTRMRECUR_EXTERNAL_FCN(FZ1,FZ2,...,FZN) calls the function
% which executes the commands in order to obtain the fuzzy output.
% It is called during the execution of FZR=GTRMRECUR_MODIF(FZ1,FZ2,...,FZN)
func = inline('cos(pi.*x1).*x2');
x_vector=cell(size(x));
for cont=1:length(x)
    x_vector{cont} = reshape(x{cont},numel(x{cont}),1);
end
iterations=length(x_vector{1});
z_cell=cell(1,2); z_cell_array=cell(1,iterations);
for iteration=1:iterations
    x1=double(x_vector{1}(iteration));
    x2=double(x_vector{2}(iteration));
    z_cell{1}=[x1 x2];
    z_cell{2}=feval(func,x_vector{1}(iteration),x_vector{2}(iteration));
    z_cell_array{iteration}=z_cell;
end
z_cell_matrix=reshape(z_cell_array,size(x{1}));
```
Fig. 3. Example of the MATLAB function called by the modified Transformation Method which contains the commands for the obtention of the fuzzy output

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

be included where information about each iteration is stored (for instance, the exact value of the fuzzy inputs at each permutation). Figures 2 and 3 show the MATLAB implementation of the efficient Transformation Method after modifying the Klimke's implementation. Note that the MATLAB implementation of the external function in Figure 3 only shows an example corresponding to the evaluation of a test function. In our case, this function would contain all

In order to prove the validity of our modified version of the Klimke's Transformation Method, a comparative fuzzy analysis was run. We took as a reference one example found in (Klimke, 2003) where the test function *f*(*x*1, *x*2) = cos(*πx*1)*x*2 is subjected to the triangular fuzzy inputs

the required commands which carry out the three foregoing calculating steps.

function [fzr]=gtrmrecur\_modif(varargin)

t = (floor((m+2-k)/2)+rem(m+2-k,2))\*2-2+k;

x{i} = repmat(shiftdim(c,1-i),repvec);

z\_cell\_matrix=gtrmrecur\_external\_fcn(x);

c = reshape(fzin{i}(end+1-k:-2:1,:)',t+2-k,1);

w = reshape(subsref(z\_cell\_matrix,indices),s^n,1);

function [z\_cell\_matrix]=gtrmrecur\_external\_fcn(x)

z\_cell=cell(1,2); z\_cell\_array=cell(1,iterations);

z\_cell\_matrix=reshape(z\_cell\_array,size(x{1}));

x1=double(x\_vector{1}(iteration)); x2=double(x\_vector{2}(iteration));

z\_cell\_array{iteration}=z\_cell;

x\_vector{cont} = reshape(x{cont},numel(x{cont}),1);

for cont=1:length(wc) wc(cont)=w{cont}{2}; end; [none,pos\_min]=min(wc); fzr{j,1} = w{pos\_min}; [none,pos\_max]=max(wc); fzr{j,2} = w{pos\_max};

% FZR=GTRMRECUR\_MODIF(FZ1,FZ2,...,FZN) calls the general transformation % method implemented by Klimke in 2003 with extended capabilities. % The function allows N fuzzy numbers in alpha-cut representation. % An accompanying function called GTRMRECUR\_EXTERNAL\_FCN is required.

if fzr{j+1,1}{2}>=fzr{j,1}{2}; fzr{j,1}=fzr{j,1}; else fzr{j,1}=fzr{j+1,1}; end if fzr{j+1,2}{2}>=fzr{j,2}{2}; fzr{j,2}=fzr{j+1,2}; else fzr{j,2}=fzr{j,2}; end

% Z\_CELL\_MATRIX=GTRMRECUR\_EXTERNAL\_FCN(FZ1,FZ2,...,FZN) calls the function % which executes the commands in order to obtain the fuzzy output. % It is called during the execution of FZR=GTRMRECUR\_MODIF(FZ1,FZ2,...,FZN)

z\_cell{2}=feval(func,x\_vector{1}(iteration),x\_vector{2}(iteration));

which contains the commands for the obtention of the fuzzy output

Fig. 2. MATLAB implementation of the efficient Transformation Method after modifying

Fig. 3. Example of the MATLAB function called by the modified Transformation Method

% Angel L. Morales, 2011.

m = size(fzin{1},1) - 1; indices.subs = {}; indices.type = '()'; fzr = cell(m+1,2); for k = 1:2

fzin=varargin; n = nargin;

for i = 1:n

end s = k;

c = c(3-k:end); repvec = t\*ones(1,n); repvec(i) = 1;

for j = m+2-k:-2:1 indices.subs(1:n) = {1:s};

wc = zeros(size(w));

Klimke's script.

% Angel L. Morales, 2011.

func = inline('cos(pi.\*x1).\*x2'); x\_vector=cell(size(x)); for cont=1:length(x)

iterations=length(x\_vector{1});

for iteration=1:iterations

z\_cell{1}=[x1 x2];

s = s + 2; end end for j = m:-1:1

end

end

end

*p*˜1 = �2.5, 2.5, 2.5�TFN and *p*˜2 = �3, 2, 2�TFN, in accordance to the nomenclature proposed in (Dubois & Prade, 1980). The fuzzy output obtained with the modified version clearly matched the results obtained with the original version. Nevertheless, an unavoidable loss of computational efficiency arose.

Our modified version results more time-consuming because of the different way in which the simulated function is handled, as it can be seen in Table 1. In the original version, the *inline* class analytical expression can be executed only once with the array containing all the required permutations of fuzzy inputs at a time. In the modified version, the expression is evaluated separately at each permutation as it would be the case when a control system is simulated with several MATLAB commands. The elapsed time during these operations is called *tf cn*. As far as the obtention of the fuzzy output as a *cell* array instead of as a *double* array, this fact does not add any noticeable increase of computational cost but it does lead to heavier variables and higher memory requirements. The elapsed time during these operations is called *tmax*. In all the cases, the loss of efficiency is acceptable with regard to the benefits they provide, specially considering that the elapsed time in any finite element simulation may be on the order of seconds and it has much more restrictive memory requirements.

### **5. A case study of uncertainty propagation**

Consider an aluminium plate in free conditions, the dimensions and material properties of which are detailed in Table 2. The finite element model of this plate was carried out by means of the *OpenFEM Toolbox* of MATLAB (Balmes et al., 2009). After defining the mesh, the matrices of mass (*M*) and stiffness (*K*) which define the undamped model were obtained. The eigenvalue problem can then be solved, the natural frequencies and mode shapes being shown in Figure 4. The damping matrix (*C*) is constructed considering a structural damping of *ξ*i = 0.6% in all the modes.

The piezoelectric actuator and sensor, the properties of which are detailed in Table 3, are collocated in the centre of the plate. As described in (Dosch et al., 1992; Sims et al., 2005), the piezoelectric actuator and sensor may be modelled as a beam which provides bending moments at its ends and measures a voltage proportional to the difference of slope of its ends. Due to collocation, the localization vectors of actuator (*Lu*) and sensor (*Ly*) are identical and equal to

$$L\_u^\mathrm{T} = L\_y^\mathrm{T} = (0, \dots, -1, 0, \dots, 1, 0, \dots, 0) \tag{10}$$

Table 1. Efficiency comparison between the original and the modified Klimke's implementation. The total time (*ttot*), the time for evaluating the analytical expression (*tf cn*), the time for calculating the maxima and minima of the output fuzzy number (*tmax*) and the time for other calculations (*toth*) are shown. The weight means the size in bytes of the output fuzzy variable.

the maximum damping of the first vibration mode (*ξ*1,max). In this particular case, a gain of *g* = 0.200 entails a selected damping coefficient of *ξ*1,sel = 88.30%. Figure 5 also shows the characteristic stability limit of the Positive Position Feedback control systems and the fact that only the first and the fifth modal shapes seem to be controllable in the view of the wide loops

301 A Fuzzy Finite Element Method Programmed

g = 0.200 ξ 1,sel = 88.30% 1 2

*ξ*1,max is the highest damping achievable in the first mode

*ξ* (structural

*ξ*1,sel is the

g = 0.257 Stability limit 3 4 5

6

Real part

Fig. 5. Root locus of the control system after tuning the second order filter parameters to the

In the proposed case study the collocated pair of actuator and sensor are subjected to uncertainty. More precisely, the position of the centre of both the actuator and sensor are considered to be uncertain parameters. Additionally, another source of uncertainty corresponding to the structural damping of the plate is included as this is known to be uncertain in practice. Thus, we have a total of five different fuzzy numbers: *x*˜*a* (position "x" of the centre of the actuator), *y*˜*a* (position "y" of the centre of the actuator), *x*˜*s* (position

damping). The membership functions of these fuzzy parameters, assuming that they are both

After applying the Fuzzy Finite Element algorithm described in Section 4 one can obtain several fuzzy outputs which describe the stability and performance of the control system. In fact, it is possible to choose those variables which better fit the effects which one is interested

The gain margin *G*˜ is the increase in the system gain when the phase is −180 that will result in a marginally stable system with intersection of the −1 + *j*0 point on the Nyquist diagram. The phase margin *P*˜ is the amount phase shift of the system at unity magnitude that will result in a marginally stable system with intersection of the −1 + *j*0 point on the Nyquist diagram.

loop independently of the selected gain, whilst the selected damping coefficient ˜

"x" of the centre of the sensor), *y*˜*s* (position "y" of the centre of the sensor) and ˜

going from the poles to the zeros (Preumont, 2002).

in MATLAB for the Analysis of Uncertain Control Systems of Structures

Imaginary part

symmetric and triangular, can be seen in Figure 6.

in. In this particular case, the following are calculated:

*ξ*1,max Maximum damping coefficient of the first mode.

*ξ*1,sel Selected damping coefficient of the first mode.

The maximum damping coefficient ˜

*G*˜ Gain margin. *P*˜ Phase margin.

˜

˜

first flexible mode of vibration (*ξ*f = 0.9 and *ω*f = 900 rad/s).

Table 2. Geometries and material properties of the plate.

Table 3. Geometries and material properties of the piezoelectric (Preumont, 2002).

where the non-zero values are placed in the position corresponding to the active rotational degree of freedom of the nodes in which the piezoelectric is located.

As described in the Section 2.2, a Positive Position Feedback control strategy designed with the *Control System Toolbox* of MATLAB is to be used. Figure 5 presents the root locus of the control system after tuning the second order filter parameters (*ξ*f and *ω*f) so that the first flexible mode of vibration become actively damped. After the tuning process (*ξ*f = 0.9 and *ω*f = 900 rad/s), the next step consists of selecting the control gain which most approaches

Fig. 4. Mode shapes and natural frequencies for the plate in free conditions.

the maximum damping of the first vibration mode (*ξ*1,max). In this particular case, a gain of *g* = 0.200 entails a selected damping coefficient of *ξ*1,sel = 88.30%. Figure 5 also shows the characteristic stability limit of the Positive Position Feedback control systems and the fact that only the first and the fifth modal shapes seem to be controllable in the view of the wide loops going from the poles to the zeros (Preumont, 2002).

Fig. 5. Root locus of the control system after tuning the second order filter parameters to the first flexible mode of vibration (*ξ*f = 0.9 and *ω*f = 900 rad/s).

In the proposed case study the collocated pair of actuator and sensor are subjected to uncertainty. More precisely, the position of the centre of both the actuator and sensor are considered to be uncertain parameters. Additionally, another source of uncertainty corresponding to the structural damping of the plate is included as this is known to be uncertain in practice. Thus, we have a total of five different fuzzy numbers: *x*˜*a* (position "x" of the centre of the actuator), *y*˜*a* (position "y" of the centre of the actuator), *x*˜*s* (position "x" of the centre of the sensor), *y*˜*s* (position "y" of the centre of the sensor) and ˜ *ξ* (structural damping). The membership functions of these fuzzy parameters, assuming that they are both symmetric and triangular, can be seen in Figure 6.

After applying the Fuzzy Finite Element algorithm described in Section 4 one can obtain several fuzzy outputs which describe the stability and performance of the control system. In fact, it is possible to choose those variables which better fit the effects which one is interested in. In this particular case, the following are calculated:

*G*˜ Gain margin.

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

Width (*mm*) 300.00 Length (*mm*) 500.00 Thickness (*mm*) 3.00 Young's modulus (GPa) 72.00 Poisson's coefficient 0.30 Density (kg/m3) 2700.00 Structural damping (%) 0.60

Width (*mm*) 40 Length (*mm*) 100 Thickness (*mm*) 0.250 *d*33 (m/V) 0.30·10−9 *d*31 (m/V) -0.15·10−9 *e*31 (N/Vm) -7.5 *k*33 0.7 Young's modulus (GPa) 50 Maximum traction (MPa) 80 Maximum compression (MPa) 600 Maximum electric Field (V/mm) 2000 Density (kg/m3) 7600 Maximum strain Brittle Maximum temperature (C) 80°C - 150°C

Table 3. Geometries and material properties of the piezoelectric (Preumont, 2002).

degree of freedom of the nodes in which the piezoelectric is located.

Mode 2 f = 66.46 Hz 2

Mode 5 f = 183.42 Hz 5

Fig. 4. Mode shapes and natural frequencies for the plate in free conditions.

Mode 1 f = 64.01 Hz 1

Mode 4 f = 176.53 Hz 4

where the non-zero values are placed in the position corresponding to the active rotational

As described in the Section 2.2, a Positive Position Feedback control strategy designed with the *Control System Toolbox* of MATLAB is to be used. Figure 5 presents the root locus of the control system after tuning the second order filter parameters (*ξ*f and *ω*f) so that the first flexible mode of vibration become actively damped. After the tuning process (*ξ*f = 0.9 and *ω*f = 900 rad/s), the next step consists of selecting the control gain which most approaches

> Mode 3 f = 150.29 Hz 3

> Mode 6 f = 231.02 Hz 6

Table 2. Geometries and material properties of the plate.

*P*˜ Phase margin.

˜ *ξ*1,max Maximum damping coefficient of the first mode.

˜ *ξ*1,sel Selected damping coefficient of the first mode.

The gain margin *G*˜ is the increase in the system gain when the phase is −180 that will result in a marginally stable system with intersection of the −1 + *j*0 point on the Nyquist diagram. The phase margin *P*˜ is the amount phase shift of the system at unity magnitude that will result in a marginally stable system with intersection of the −1 + *j*0 point on the Nyquist diagram. The maximum damping coefficient ˜ *ξ*1,max is the highest damping achievable in the first mode loop independently of the selected gain, whilst the selected damping coefficient ˜ *ξ*1,sel is the

specifically analyzed via fuzzy variables, we can note that its zero crosses the imaginary axis and becomes unstable for some values of the gain. Note that zeros are essentially related to mode shapes of the structure, which tend to be more sensitive to the parameter uncertainty than the natural frequencies. This fact warns us that a higher gain level may lead to instability of the fourth mode when actively damping the first, which can be detected easily by means of including new fuzzy variables in the analysis such as the maximum real part of the closed

303 A Fuzzy Finite Element Method Programmed

in MATLAB for the Analysis of Uncertain Control Systems of Structures

Fig. 8. Root locus of the control system for two different *α*-cuts: *μ* = 1 (black) and *μ* = 0

methodology for the analysis and simulation of uncertain control of structures.

accuracy results together with low computational costs.

one method stand out from the other.

In this chapter we have shown the suitability of the matrix computing environment MATLAB to solve control problems of structures via finite element analysis and root locus design and how these two issues have been blended and integrated into a Fuzzy Finite Element

The Fuzzy Finite Element procedure has been described. The approach relies on the *α*-cut technique which subdivides the membership range into a number of *α*-levels. The fuzzy arithmetic is based upon an efficient implementation in MATLAB code for the Transformation Method. This implementation has been modified and improved with some enhanced features. This methodology has been shown to be applicable to a realistic case study related to structural control using a Positive Position Feedback technique and piezoelectric actuators and sensors. In the light of the obtained results, the proposed procedure may prove useful to other analysis involving any other control technique or uncertainty source. In addition, the efficient implementation obtained due to MATLAB programming allows us to achieve high

Finally, further research may be carried out in this topic. Firstly, a more efficient implementation of the Transformation Method can be found when using *cell* arrays instead of *double* arrays or when considering multiple fuzzy outputs. Secondly, the Transformation Method also allows a sort of sensitivity analysis which should be included in the MATLAB implementation in order to endow our method with extended capabilities. Thirdly, and in-depth comparison between this possibilistic methodology and the robust *H*∞ optimal control theory must be carried out in order to identify more precisely the situations in which

loop poles.

Imaginary part

**6. Conclusion**

(grey).

Fig. 6. Membership functions of the input fuzzy numbers: (a) coordinates of the centre of the actuator (*x*˜a and *x*˜s), (b) coordinates of the centre of the sensor (*y*˜a and *y*˜s) and (c) structural damping of the plate ( ˜ *ξ*).

damping that the control systems provides at the value of gain selected in the deterministic case. Clearly, the gain and phase margins can be used to quantify controller stability, whilst the damping ratios serve to quantify performance and reliability.

The membership functions of these variables are shown in Figure 7. We observe that the uncertainty considered may lead to a both better and worse behaviour of the control system in comparison to the considered deterministic situation.

Another possibility for obtaining a global view of the whole control system lies in plotting its root locus for different *α*-sublevels or *α*-cuts. In particular, we show in Figure 8 the root locus for two different values of the membership. Although the fourth mode has not been

Fig. 7. Membership function of the output fuzzy numbers: (a) stability indicators (*G*˜ in black color and *P*˜ in gray color) and (b) performance indicators ( ˜ *ξ*1,max in black color and ˜ *ξ*1,sel in gray color).

specifically analyzed via fuzzy variables, we can note that its zero crosses the imaginary axis and becomes unstable for some values of the gain. Note that zeros are essentially related to mode shapes of the structure, which tend to be more sensitive to the parameter uncertainty than the natural frequencies. This fact warns us that a higher gain level may lead to instability of the fourth mode when actively damping the first, which can be detected easily by means of including new fuzzy variables in the analysis such as the maximum real part of the closed loop poles.

Fig. 8. Root locus of the control system for two different *α*-cuts: *μ* = 1 (black) and *μ* = 0 (grey).

## **6. Conclusion**

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

(a) (b) (c)

245 247 249 251 253 255

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

86.0 86.8 87.6 88.4 89.2 90.0

ξ1,max (%), ξ1,sel (%)

*ξ*1,max in black color and ˜

*ξ*1,sel in

ξ (%) 0.1 0.3 0.5 0.7 0.9 1.1

μ(ξ)

y , y (mm) a s

Fig. 6. Membership functions of the input fuzzy numbers: (a) coordinates of the centre of the actuator (*x*˜a and *x*˜s), (b) coordinates of the centre of the sensor (*y*˜a and *y*˜s) and (c) structural

damping that the control systems provides at the value of gain selected in the deterministic case. Clearly, the gain and phase margins can be used to quantify controller stability, whilst

The membership functions of these variables are shown in Figure 7. We observe that the uncertainty considered may lead to a both better and worse behaviour of the control system

Another possibility for obtaining a global view of the whole control system lies in plotting its root locus for different *α*-sublevels or *α*-cuts. In particular, we show in Figure 8 the root locus for two different values of the membership. Although the fourth mode has not been

(a) (b)

μ(ξ ) 1,sel

μ(ξ ), 1,max 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

Fig. 7. Membership function of the output fuzzy numbers: (a) stability indicators (*G*˜ in black

145 147 149 151 153 155

x , x (mm) a s

*ξ*).

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

the damping ratios serve to quantify performance and reliability.

in comparison to the considered deterministic situation.

0.2565 0.2570 0.2575 0.2580 0.2585 0.2590

G 101.70 102.00 102.30 102.60 102.90 103.20

P (deg)

color and *P*˜ in gray color) and (b) performance indicators ( ˜

μ(y ), μ(y ) a s

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

gray color).

μ(G), μ(P)

damping of the plate ( ˜

μ(x ), μ(x ) a s

> In this chapter we have shown the suitability of the matrix computing environment MATLAB to solve control problems of structures via finite element analysis and root locus design and how these two issues have been blended and integrated into a Fuzzy Finite Element methodology for the analysis and simulation of uncertain control of structures.

> The Fuzzy Finite Element procedure has been described. The approach relies on the *α*-cut technique which subdivides the membership range into a number of *α*-levels. The fuzzy arithmetic is based upon an efficient implementation in MATLAB code for the Transformation Method. This implementation has been modified and improved with some enhanced features. This methodology has been shown to be applicable to a realistic case study related to structural control using a Positive Position Feedback technique and piezoelectric actuators and sensors. In the light of the obtained results, the proposed procedure may prove useful to other analysis involving any other control technique or uncertainty source. In addition, the efficient implementation obtained due to MATLAB programming allows us to achieve high accuracy results together with low computational costs.

> Finally, further research may be carried out in this topic. Firstly, a more efficient implementation of the Transformation Method can be found when using *cell* arrays instead of *double* arrays or when considering multiple fuzzy outputs. Secondly, the Transformation Method also allows a sort of sensitivity analysis which should be included in the MATLAB implementation in order to endow our method with extended capabilities. Thirdly, and in-depth comparison between this possibilistic methodology and the robust *H*∞ optimal control theory must be carried out in order to identify more precisely the situations in which one method stand out from the other.

**Part 4**

**Modelling, Identification and Simulation**

#### **7. Acknowledgment**

The authors are grateful for the support received from the "Ministerio de Educación" of the Spanish government through the "Programa Nacional de Movilidad de Recursos Humanos del Plan Nacional de I-D+i 2008-2011" and from the EPSRC through grant reference EP/D078601/1.

#### **8. References**

URL: *http://support.sdtools.com/gf/project/openfem*

OpenFEM (2011).

## **Part 4**

## **Modelling, Identification and Simulation**

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

304 MATLAB – A Ubiquitous Tool for the Practical Engineer

The authors are grateful for the support received from the "Ministerio de Educación" of the Spanish government through the "Programa Nacional de Movilidad de Recursos Humanos del Plan Nacional de I-D+i 2008-2011" and from the EPSRC through grant reference

Balmes, E., Bianchi, J.-P. & Leclére, J.-M. (2009). *Structural Dynamics Toolbox for use with*

De Gersem, H., Moens, D., Desmet, W. & Vandepitte, D. (2005). A fuzzy finite element

De Munck, M., Moens, D., Desmet, W. & Vandepitte, D. (2008). A response surface based

Dosch, J. J., Daniel, D. J. & Garcia, E. (1992). A self-sensing piezoelectric actuator for collocated control, *Journal of Intelligent Material Systems and Structures* 3: 166–185. Dubois, D. & Prade, H. M. (1980). *Fuzzy sets and systems : theory and applications*, Mathematics

Gu, D.-W., Petkov, P. H. & Konstantinov, M. M. (2005). *Robust control design with MATLAB*,

Hanss, M. (2002). The transformation method for the simulation and analysis of systems with

Klimke, A. (2003). An efficient implementation of the transformation method of fuzzy

Moens, D. & Vandepitte, D. (2005). A fuzzy finite element procedure for the calculation of

Oberkampf, W. L., Helton, J. C., Joslyn, C. A., Wojtkiewicz, S. F. & Ferson, S. (2004). Challenge

Preumont, A. (2002). *Vibration control of active structures: An introduction*, Kluwer Academic,

Sims, N. D., Bayly, P. V. & Young, K. A. (2005). Piezoelectric sensors and actuators for milling

Sims, N. D., Manson, G. & Mann, B. (2010). Fuzzy stability analysis of regenerative chatter in

Zames, G. (1981). Feedback and optimal sensitivity: Model reference transformations,

Zames, G. & Francis, B. A. (1983). Feedback, minimax sensitivity and optimal robustness,

multiplicative seminorms and approximate inverses, *IEEE Transactions on Automatic*

tool stability lobes, *Journal of Sound and Vibration* 281: 743–762.

milling, *Journal of Sound and Vibration* 329: 1025–1041. Zadeh, L. A. (1965). Fuzzy sets, *Information and Control* 8: 338–353.

*IEEE Transactions on Automatic Control* 28: 585–601.

arithmetic, (IANS), extended preprint 2003/009, *Technical report*, University of

uncertain frequency-response functions of damped structures: Part 1 - procedure,

problems: uncertainty in system response given uncertain parameters, *Reliability*

uncertain properties, *Computers & Structures* 86: 1080–1092.

in science and engineering, Academic Press, New York.

uncertain parameters, *Fuzzy Sets and Systems* 130: 277–289.

*Journal of Sound and Vibration* 288: 431–462.

URL: *http://support.sdtools.com/gf/project/openfem*

*Engineering & System Safety* 85: 11–19.

procedure for the calculation of uncertain frequency response functions of damped structures: Part 2 - numerical case studies, *Journal of Sound and Vibration* 288: 463–486.

optimisation algorithm for the calculation of fuzzy envelope frfs of models with

**7. Acknowledgment**

EP/D078601/1.

**8. References**

*MATLAB: User's guide*.

Springer, London.

Stuttgart.

Dordrecht.

*Control* 26: 301–320.

OpenFEM (2011).

**16**

*India*

**Visual and Thermal Image**

*1Division of Avionics, Department of Aerospace Engineering,*

*2Department of Electronics and Communication Engineering,*

*3Division of Avionics, Department of Aerospace Engineering,*

Unmanned aerial vehicles (UAVs) are aircrafts which have the capability of flight without an onboard pilot. These vehicles are remotely controlled, semi-autonomous, autonomous, or have a combination of these capabilities. UAV's has its applications in a whole lot of domains. Image processing applications with specific importance to surveillance and

UAVs are equipped with imaging sensor platform, which operates remotely controlled, semi-autonomously or autonomously, without a pilot sitting in the vehicle. The platform may have a small or medium size still-video or video camera, thermal or infrared camera systems, airborne light detection and ranging (LIDAR) system, or a combination thereof. All these different kinds of cameras are an effective sensor tool which is portable, light weight

Thermal images have a valuable advantage over the visual images. Thermal images do not depend on the illumination, the output is the projection of thermal sensors of the emissions of heat of the objects. This unique merit gives rise for effective segmentation of objects.

With the development of new imaging sensors arises the need of a meaningful combination of all employed imaging sources. Image fusion of visual and thermal sensing outputs adds a new dimension in making the target tracking more reliable. Target tracking at instances of smoke, fog, cloudy conditions gets improved. With conditions of same background colour perception of target unnoticed getting eliminated with thermal image inclusion, image fusion gives complementary information. A holistic system which represents the combined

**1. Introduction**

reconnaissance is of immense interest.

and airborne in a platform on the UAV.

Ultimately, surveillance measure using an UAV gets improved.

fused data is perceived at the control level of the UAV's.

*Madras Institute of Technology, Anna University,*

*Madras Institute of Technology, Anna University*

*Madras Institute of Technology, Anna University 4National Aerospace Laboratories (NAL) Bangalore,*

**Fusion for UAV Based**

K. Senthil Kumar1, G. Kavitha2, R. Subramanian3 and G. Ramesh4

**Target Tracking**

## **Visual and Thermal Image Fusion for UAV Based Target Tracking**

## K. Senthil Kumar1, G. Kavitha2,

R. Subramanian3 and G. Ramesh4

*1Division of Avionics, Department of Aerospace Engineering, Madras Institute of Technology, Anna University, 2Department of Electronics and Communication Engineering, Madras Institute of Technology, Anna University 3Division of Avionics, Department of Aerospace Engineering, Madras Institute of Technology, Anna University 4National Aerospace Laboratories (NAL) Bangalore, India*

## **1. Introduction**

Unmanned aerial vehicles (UAVs) are aircrafts which have the capability of flight without an onboard pilot. These vehicles are remotely controlled, semi-autonomous, autonomous, or have a combination of these capabilities. UAV's has its applications in a whole lot of domains. Image processing applications with specific importance to surveillance and reconnaissance is of immense interest.

UAVs are equipped with imaging sensor platform, which operates remotely controlled, semi-autonomously or autonomously, without a pilot sitting in the vehicle. The platform may have a small or medium size still-video or video camera, thermal or infrared camera systems, airborne light detection and ranging (LIDAR) system, or a combination thereof. All these different kinds of cameras are an effective sensor tool which is portable, light weight and airborne in a platform on the UAV.

 Thermal images do not depend on the illumination, the output is the projection of thermal sensors of the emissions of heat of the objects. This unique merit gives rise for effective segmentation of objects. Ultimately, surveillance measure using an UAV gets improved.

With the development of new imaging sensors arises the need of a meaningful combination of all employed imaging sources. Image fusion of visual and thermal sensing outputs adds a new dimension in making the target tracking more reliable. Target tracking at instances of smoke, fog, cloudy conditions gets improved. With conditions of same background colour perception of target unnoticed getting eliminated with thermal image inclusion, image fusion gives complementary information. A holistic system which represents the combined fused data is perceived at the control level of the UAV's.

Visual and Thermal Image Fusion for UAV Based Target Tracking 309

Segmentation is the key and the first step to automatic target recognition, which will directly affect the accuracy of the following work. As a result, the division methods and its precision degree are essential. Infrared heat wave image is different from the visible light images. It reflects the distribution of the object surface temperature and latent characteristics of

The infrared heat radiation, due to the imperfections of the system, will bring a variety of noise in the imaging process. The noise of complex distribution of infrared images makes

The two dimensional Otsu algorithm is given as follows. Suppose an image pixel size is M × N, gray-scale of the image ranges from 0 to L-1. The neighborhood average gray g (m, n) of

Calculating the average neighbourhood gray of each pixel point, a gray binary group (i, j) may form. We use C*ij to* represent the occurrence frequency of (i, j). Then the probability P*ij*

Assuming the existence of two classes *C0* and *C1* in Two dimensional form, the histogram represents their respective goals and background, and with two different probability density distribution function. If making use of two-dimensional histogram threshold vector (s, t) to segment the image (of which 0 ≤ s, t <L), then the probability of two classes are respectively:

��� �

��� �����

σB = ω0 [(µ0-µτ)( µ0-µτ)T] + ω1 [(µ1-µτ)( µ1-µτ)T] (5)

ω1 = P(C1) = ∑ ∑ P��

corresponding threshold of segmentation is the optimal threshold (S, T), namely:

���

When the track of the above-mentioned dispersion matrix gets the maximum, the

 tr(σB (S,T)) = ������,��� {tr(σB (S,T))} (6) We know that 2-D thermal images with noise segmented by Otsu way may get better results compared to one dimensional threshold segmentation methods. However, the computation cost gets huge, which is because the determination of the optimal threshold need to travel all the s and t, of which 0 ≤ s, t < L. That is to say, the more gray scale value of images is, the

(���)��

��� ∑ ∑ f(m + i, n + j) (���)�� ���(���)��

���(���)�� (1)

����� (2)

��� = ω0 (s,t) (3)

����� = ω1 (s,t) (4)

**3. Image segmentation**

**3.1 2D OTSU algorithm**

the signal to noise ratio lower than visible light images.

the coordinate definition (m, n) pixel point is as follows:

g (m,n) = �

��� ���

ω0 = P(C0) = ∑ ∑ P�� �

��� =1.

of vector (i, j) may be determined by the formula:

Pij = ���

���

The probability of background occurrence is:

The probability of object occurrence is:

The definition of dispersion matrix:

longer choice time of the threshold is.

Here, 0 ≤ I, j < L, and ∑ ∑ ���

material form.

## **2. Thermal imaging**

Thermography which uses black body radiation law makes it to have information gathering without visible illumination. Thermal imaging cameras detect radiation in the infrared (IR) range of the electromagnetic spectrum (3-6 µm and 8–14 µm).

The charge coupled device (CCD) and complementary metal oxide semiconductor (CMOS) sensors are used for visible light cameras. These can detect only the non thermal part of the infrared spectrum called near-infrared (NIR). On the other hand, thermal imaging cameras make use specialized focal plane arrays (FPAs) that respond to longer wavelengths (midand long-wavelength infrared).

There is also a difference between how far one can see with a cooled and with an uncooled thermal imaging camera. Cooled camera systems are more expensive, but generally have a longer range than uncooled systems under many conditions. Extremely long range thermal imaging applications are best served by cooled camera systems. This is particularly true in the midwave band in humid atmospheric conditions.

The heat radiation is focused onto special receptors in the camera which convert it into a format which is displayed on a monitor in monochrome which is recognisable by the human eye. The objects emitting the greatest intensity of heat are usually presented as the darkest (black) in the greyscale, i.e. known as 'black-hot'. Many cameras have a function whereby the functionality can be switched from 'blackhot' to 'white-hot' and back again at the operator's wish. Probably the greatest 'enemy' of thermal imaging is extended rainfall since that has the effect of cooling all inanimate objects and severely reducing contrast.

Thermal imaging makes it possible of real time target tracking. Detection of targets in dark and low light conditions can be done. All weather operation and dull, dirty and dangerous (DDD) roles are possible.

Thermal imaging cameras produce a clear image in the darkest of nights, in light fog and smoke and in the most diverse weather conditions. There has also been an increased interest in thermal imaging for all kinds of security applications, from long-range surveillance at border crossings, truck and shipping container inspection, to the monitoring of high-security installations such as nuclear power stations, airports, and dams. But thermal imaging has a lot more to offer than just a night vision solution for security applications.

Car manufacturers are integrating night vision modules for driver vision enhancement into cars. By helping drivers to see at night, accidents can be avoided. Boats and yachts are being equipped with thermal imaging cameras for night time navigation and other maritime applications like man overboard searches.

Often the thermal imager is just a small part of the complete system as in an UAS, so it needs to be as small, light and inexpensive as possible. Low-cost thermal imager is used as a pilot's night vision enhancement. It helps pilots by enhancing the ability to see terrain and other aircraft at long ranges, even in total darkness, light fog, dust and smoke. Thermal imaging is a technology that enables detection of people and objects in total darkness and in very diverse weather conditions.

A typical application for thermal imaging is border security, where most threats occur at night. Thermal imaging allows the aircraft to fly in total darkness and detect targets through smoke. The same aircraft can also be used to detect such things as forest fires. Areas which are hotter than the surroundings can indicate the start of a fire and can clearly be seen on a thermal image.

#### **3. Image segmentation**

308 MATLAB – A Ubiquitous Tool for the Practical Engineer

Thermography which uses black body radiation law makes it to have information gathering without visible illumination. Thermal imaging cameras detect radiation in the infrared (IR)

The charge coupled device (CCD) and complementary metal oxide semiconductor (CMOS) sensors are used for visible light cameras. These can detect only the non thermal part of the infrared spectrum called near-infrared (NIR). On the other hand, thermal imaging cameras make use specialized focal plane arrays (FPAs) that respond to longer wavelengths (mid-

There is also a difference between how far one can see with a cooled and with an uncooled thermal imaging camera. Cooled camera systems are more expensive, but generally have a longer range than uncooled systems under many conditions. Extremely long range thermal imaging applications are best served by cooled camera systems. This is particularly true in

The heat radiation is focused onto special receptors in the camera which convert it into a format which is displayed on a monitor in monochrome which is recognisable by the human eye. The objects emitting the greatest intensity of heat are usually presented as the darkest (black) in the greyscale, i.e. known as 'black-hot'. Many cameras have a function whereby the functionality can be switched from 'blackhot' to 'white-hot' and back again at the operator's wish. Probably the greatest 'enemy' of thermal imaging is extended rainfall since

Thermal imaging makes it possible of real time target tracking. Detection of targets in dark and low light conditions can be done. All weather operation and dull, dirty and dangerous

Thermal imaging cameras produce a clear image in the darkest of nights, in light fog and smoke and in the most diverse weather conditions. There has also been an increased interest in thermal imaging for all kinds of security applications, from long-range surveillance at border crossings, truck and shipping container inspection, to the monitoring of high-security installations such as nuclear power stations, airports, and dams. But thermal imaging has a lot more to offer than just a night vision solution for

Car manufacturers are integrating night vision modules for driver vision enhancement into cars. By helping drivers to see at night, accidents can be avoided. Boats and yachts are being equipped with thermal imaging cameras for night time navigation and other maritime

Often the thermal imager is just a small part of the complete system as in an UAS, so it needs to be as small, light and inexpensive as possible. Low-cost thermal imager is used as a pilot's night vision enhancement. It helps pilots by enhancing the ability to see terrain and other aircraft at long ranges, even in total darkness, light fog, dust and smoke. Thermal imaging is a technology that enables detection of people and objects in total darkness and in

A typical application for thermal imaging is border security, where most threats occur at night. Thermal imaging allows the aircraft to fly in total darkness and detect targets through smoke. The same aircraft can also be used to detect such things as forest fires. Areas which are hotter than the surroundings can indicate the start of a fire and can clearly be seen on a

that has the effect of cooling all inanimate objects and severely reducing contrast.

range of the electromagnetic spectrum (3-6 µm and 8–14 µm).

the midwave band in humid atmospheric conditions.

**2. Thermal imaging**

and long-wavelength infrared).

(DDD) roles are possible.

security applications.

applications like man overboard searches.

very diverse weather conditions.

thermal image.

Segmentation is the key and the first step to automatic target recognition, which will directly affect the accuracy of the following work. As a result, the division methods and its precision degree are essential. Infrared heat wave image is different from the visible light images. It reflects the distribution of the object surface temperature and latent characteristics of material form.

The infrared heat radiation, due to the imperfections of the system, will bring a variety of noise in the imaging process. The noise of complex distribution of infrared images makes the signal to noise ratio lower than visible light images.

### **3.1 2D OTSU algorithm**

The two dimensional Otsu algorithm is given as follows. Suppose an image pixel size is M × N, gray-scale of the image ranges from 0 to L-1. The neighborhood average gray g (m, n) of the coordinate definition (m, n) pixel point is as follows:

$$\log \text{(m,n)} = \frac{1}{\text{kxk}} \sum\_{\mathbf{l}=\mathbf{\bar{\cdot}}=\mathbf{\bar{\cdot}}=\mathbf{\bar{\cdot}}=\mathbf{\bar{\cdot}}=\mathbf{\bar{\cdot}}=\mathbf{\bar{\cdot}}=\mathbf{\bar{\cdot}}=\mathbf{\bar{\cdot}}}{} \text{f} \mathbf{\bar{\cdot}} (\mathbf{m} + \mathbf{\bar{\cdot}}, \mathbf{n} + \mathbf{j}) \tag{1}$$

Calculating the average neighbourhood gray of each pixel point, a gray binary group (i, j) may form. We use C*ij to* represent the occurrence frequency of (i, j). Then the probability P*ij* of vector (i, j) may be determined by the formula:

$$P\_{\text{ij}} = \frac{\text{Cl}\,\text{j}}{\text{M} \,\text{X} \,\text{N}} \tag{2}$$

Here, 0 ≤ I, j < L, and ∑ ∑ ��� ��� ��� ��� ��� =1.

Assuming the existence of two classes *C0* and *C1* in Two dimensional form, the histogram represents their respective goals and background, and with two different probability density distribution function. If making use of two-dimensional histogram threshold vector (s, t) to segment the image (of which 0 ≤ s, t <L), then the probability of two classes are respectively: The probability of background occurrence is:

$$\alpha\_0 = \mathbf{P}(\mathbf{C}\_0) = \sum\_{l=0}^{S} \sum\_{j=0}^{U} \mathbf{P}\_{l} = \alpha\_0 \text{ (s,t)}\tag{3}$$

The probability of object occurrence is:

$$\alpha\_1 = \mathbf{P(C\_1)} = \sum\_{l=s+1}^{L-1} \sum\_{l=t+1}^{L-1} \mathbf{P\_{l[}} = \alpha\_1 \text{ (s,t)} \tag{4}$$

The definition of dispersion matrix:

$$\mathbf{a} \cdot \mathbf{a} \mathbf{B} = \mathbf{a} \mathbf{0} \left[ (\mu \mathbf{0} \text{-} \mu \text{t}) (\mu \mathbf{0} \text{-} \mu \text{t}) \mathbf{I} \right] + \mathbf{a} \mathbf{1} \left[ (\mu \mathbf{1} \text{-} \mu \text{t}) (\mu \mathbf{1} \text{-} \mu \text{t}) \mathbf{I} \right] \tag{5}$$

When the track of the above-mentioned dispersion matrix gets the maximum, the corresponding threshold of segmentation is the optimal threshold (S, T), namely:

$$\text{tr}(\text{o}\_{\text{B}}(\text{S}, \text{T})) = \max\_{0 \le s, t \le L} \{ \text{tr}(\text{o}\_{\text{B}}(\text{S}, \text{T})) \}\tag{6}$$

We know that 2-D thermal images with noise segmented by Otsu way may get better results compared to one dimensional threshold segmentation methods. However, the computation cost gets huge, which is because the determination of the optimal threshold need to travel all the s and t, of which 0 ≤ s, t < L. That is to say, the more gray scale value of images is, the longer choice time of the threshold is.

Visual and Thermal Image Fusion for UAV Based Target Tracking 311

When u = 4, the system is the chaos state of the biggest ergodicity. Genetic algorithm is a kind of random search algorithm drawn from natural selection and natural genetic mechanisms of biology, which is particularly well suited to deal with the complex and

The segmentation results are shown in Fig. 1 to Fig. 4 with the segmented image of the

Image fusion is the process by which two or more images are combined into a single image retaining the important features from each of the original images. The fusion of images is often required for images acquired from different instrument modalities or capture techniques of the same scene or objects. Important applications of the fusion of images include medical imaging, microscopic imaging, remote sensing, computer vision, and

nonlinear problems that traditional search methods are difficult to resolve.

Fig. 3. Thermal Image

Fig. 4. Segmented Image

**4. Image fusion process**

thermal image.

The segmentation of various thermal images is illustrated here. The segmentation results can be made more efficient in identification of targets and optimized by a number of methods. These segmentation results are used in fused image target tracking. One such method will be determining optimum threshold using histogram analysis. The method of using chaos based genetic algorithm makes the process time lower. The chaos based genetic algorithm uses Otsu algorithm as fitness function and proceeds for segmentation.

The so-called chaos refers to the uncertainty in the system appearing in seemingly without rules. From the mathematical sense, for determined series initial values, it is possible to predict the long-term behaviour of the system, even know its past behaviour and state by the power system. It is a unique phenomenon of non-linear systems. Chaos has randomicity, ergodicity and regularity. Based on these characteristics of chaos, using chaotic variables to optimize search is no doubt better than random search. Logistic map is the most basic chaotic map. The chaos equation of the Logistic definition can be described as follows:

$$\chi\_{\{\mathfrak{n}+1\}} = \mathfrak{u} \propto\_{\mathfrak{n}} (1 \cdot \chi\_{\mathfrak{n}}) \; ; \; 0 \le \mathfrak{u} \le 4 \tag{7}$$

xn lies between 0 and 1. n varies from 0,1,2, etc.

Fig. 1. Thermal Image

Fig. 2. Segmented Image

When u = 4, the system is the chaos state of the biggest ergodicity. Genetic algorithm is a kind of random search algorithm drawn from natural selection and natural genetic mechanisms of biology, which is particularly well suited to deal with the complex and nonlinear problems that traditional search methods are difficult to resolve.

Fig. 3. Thermal Image

310 MATLAB – A Ubiquitous Tool for the Practical Engineer

The segmentation of various thermal images is illustrated here. The segmentation results can be made more efficient in identification of targets and optimized by a number of methods. These segmentation results are used in fused image target tracking. One such method will be determining optimum threshold using histogram analysis. The method of using chaos based genetic algorithm makes the process time lower. The chaos based genetic

The so-called chaos refers to the uncertainty in the system appearing in seemingly without rules. From the mathematical sense, for determined series initial values, it is possible to predict the long-term behaviour of the system, even know its past behaviour and state by the power system. It is a unique phenomenon of non-linear systems. Chaos has randomicity, ergodicity and regularity. Based on these characteristics of chaos, using chaotic variables to optimize search is no doubt better than random search. Logistic map is the most basic chaotic map. The chaos equation of the Logistic definition can be described as follows:

x(n+1) = u xn (1- xn) ; 0 < u ≤ 4 (7)

algorithm uses Otsu algorithm as fitness function and proceeds for segmentation.

xn lies between 0 and 1. n varies from 0,1,2, etc.

Fig. 1. Thermal Image

Fig. 2. Segmented Image

Fig. 4. Segmented Image

The segmentation results are shown in Fig. 1 to Fig. 4 with the segmented image of the thermal image.

## **4. Image fusion process**

Image fusion is the process by which two or more images are combined into a single image retaining the important features from each of the original images. The fusion of images is often required for images acquired from different instrument modalities or capture techniques of the same scene or objects. Important applications of the fusion of images include medical imaging, microscopic imaging, remote sensing, computer vision, and

Visual and Thermal Image Fusion for UAV Based Target Tracking 313

The Haar wavelet is a certain sequence of rescaled "square-shaped" functions which together form a wavelet family or basis. Wavelet analysis is similar to Fourier analysis in that it allows a target function over an interval to be represented in terms of an orthonormal

Neural Network and Fuzzy Logic approach can be used for sensor fusion. Such a sensor fusion could belong to a class of sensor fusion in which case the features could be input and decision could be output. The system can be trained from the input data obtained from the

The basic concept is to associate the given sensory inputs with some decision outputs. The following algorithm for pixel level image fusion using Fuzzy Logic illustrate the process of defining membership functions and rules for the image fusion process using FIS (Fuzzy Inference System) editor of Fuzzy Logic toolbox in MATLAB. The process flow of the

1. The visual and thermal image forms the inputs for the fusion system. The inputs must

2. With a gray level conversion applied, the two images are transformed to a column

3. The number and type of membership functions for both the input images are given to the FIS. The rule base for FIS decides how the output fused image should be.

The fuzzy system considered here is Mamdani type and its' FIS model in MATLAB is shown in Fig. 6. The Mamdani rule base is a crisp model of a system, i.e. it takes crisp inputs and produces crisp outputs. It does this with the use of user-defined fuzzy rules on user-defined

4. The fused image is converted back to matrix format from column form.

function basis. The Haar wavelet is also the simplest possible wavelet.

Fig. 5. Wavelet Multi Resolution Analysis

**4.2 FIS based image fusion**

sensors.

process is as follows:

form.

fuzzy variables.

be with the same field of view.

robotics. Fusion techniques include the simplest method of pixel averaging to more complicated methods such as principal component analysis (PCA) and wavelet transform (WT) based fusion.

Image fusion is the process of relevant information extraction from two or more images. The resulting image will encompass more information than any of the given input images. A multi sensor, multi temporal and multi view information techniques are required which overcome the limitations of using single senor information. The benefits of image fusion are extended range of operation, extended spatial and temporal resolution and coverage, reduced uncertainty, higher accuracy, reliability and compact representation.

Various kinds of image fusion for visual, Infrared and Synthetic Aperture Radar (SAR) images exist. Some of the primitive fusion process algorithms also have disadvantages. Direct fusion method makes the image blurred. The pixel based image fusion is computationally complex for large and high resolution images. The image averaging method produces a reduced contrast of information. The requirement is thus a novel and efficient information fusion process. Image fusion has a significant role of recognition of targets and objects. Target identification, localisation, filtering and data association forms an important application of the fusion process. Thus an effective surveillance and reconnaissance system can be formed.

There is information which is redundant and at the same time complementary too. The following summarize several approaches to the pixel level fusion of spatially registered input images. A generic categorization of image fusion methods in the following:

• Linear Superposition, Nonlinear Methods, Expectation Maximisation, Image Pyramids, Wavelet Transform, Generic Multiresolution Fusion Scheme, Optimization Approaches, Artificial Neural Networks, Fuzzy Techniques

Some of the prominent applications of image fusion are:

• Concealed Weapon Detection, Night Time Surveillance, Automatic Landing System, Digital Camera Applications, Medical Diagnosis, Defect Inspection and Remote Sensing

#### **4.1 Wavelet transform**

A signal analysis method similar to image pyramids is the discrete wavelet transform. The main difference is that while image pyramids lead to an over complete set of transform coefficients, the wavelet transform results in a nonredundant image representation. The discrete two dimensional wavelet transform is computed by the recursive application of low pass and high pass filters in each direction of the input image (i.e. rows and columns) followed by sub sampling.

These basis functions or baby wavelets are obtained from a single prototype wavelet called the mother wavelet, by dilations or contractions (scaling) and translations (shifts). They have advantages over traditional Fourier methods in analyzing physical situations where the signal contains discontinuities and sharp spikes.

Image fusion process is achieved by multiresolution decomposition at fourth level. The multiwavelet decomposition coefficients of the input images are appropriately merged and a new fixed image is obtained by reconstructing the fused multiwavelet coefficients. The theory of multiwavelets is also based on the idea of multiresolution analysis (MRA) as shown in Fig. 5. During a single level of decomposition using a scalar wavelet transform, the two-dimensional (2-D) image data is replaced with four blocks corresponding to the subbands representing either low-pass or high-pass filtering in each direction.

robotics. Fusion techniques include the simplest method of pixel averaging to more complicated methods such as principal component analysis (PCA) and wavelet transform

Image fusion is the process of relevant information extraction from two or more images. The resulting image will encompass more information than any of the given input images. A multi sensor, multi temporal and multi view information techniques are required which overcome the limitations of using single senor information. The benefits of image fusion are extended range of operation, extended spatial and temporal resolution and coverage,

Various kinds of image fusion for visual, Infrared and Synthetic Aperture Radar (SAR) images exist. Some of the primitive fusion process algorithms also have disadvantages. Direct fusion method makes the image blurred. The pixel based image fusion is computationally complex for large and high resolution images. The image averaging method produces a reduced contrast of information. The requirement is thus a novel and efficient information fusion process. Image fusion has a significant role of recognition of targets and objects. Target identification, localisation, filtering and data association forms an important application of the fusion process. Thus an effective surveillance and

There is information which is redundant and at the same time complementary too. The following summarize several approaches to the pixel level fusion of spatially registered

• Linear Superposition, Nonlinear Methods, Expectation Maximisation, Image Pyramids, Wavelet Transform, Generic Multiresolution Fusion Scheme, Optimization Approaches,

• Concealed Weapon Detection, Night Time Surveillance, Automatic Landing System, Digital Camera Applications, Medical Diagnosis, Defect Inspection and Remote Sensing

A signal analysis method similar to image pyramids is the discrete wavelet transform. The main difference is that while image pyramids lead to an over complete set of transform coefficients, the wavelet transform results in a nonredundant image representation. The discrete two dimensional wavelet transform is computed by the recursive application of low pass and high pass filters in each direction of the input image (i.e. rows and columns)

These basis functions or baby wavelets are obtained from a single prototype wavelet called the mother wavelet, by dilations or contractions (scaling) and translations (shifts). They have advantages over traditional Fourier methods in analyzing physical situations where the

Image fusion process is achieved by multiresolution decomposition at fourth level. The multiwavelet decomposition coefficients of the input images are appropriately merged and a new fixed image is obtained by reconstructing the fused multiwavelet coefficients. The theory of multiwavelets is also based on the idea of multiresolution analysis (MRA) as shown in Fig. 5. During a single level of decomposition using a scalar wavelet transform, the two-dimensional (2-D) image data is replaced with four blocks corresponding to the

subbands representing either low-pass or high-pass filtering in each direction.

input images. A generic categorization of image fusion methods in the following:

reduced uncertainty, higher accuracy, reliability and compact representation.

(WT) based fusion.

reconnaissance system can be formed.

**4.1 Wavelet transform**

followed by sub sampling.

Artificial Neural Networks, Fuzzy Techniques Some of the prominent applications of image fusion are:

signal contains discontinuities and sharp spikes.

#### Fig. 5. Wavelet Multi Resolution Analysis

The Haar wavelet is a certain sequence of rescaled "square-shaped" functions which together form a wavelet family or basis. Wavelet analysis is similar to Fourier analysis in that it allows a target function over an interval to be represented in terms of an orthonormal function basis. The Haar wavelet is also the simplest possible wavelet.

## **4.2 FIS based image fusion**

Neural Network and Fuzzy Logic approach can be used for sensor fusion. Such a sensor fusion could belong to a class of sensor fusion in which case the features could be input and decision could be output. The system can be trained from the input data obtained from the sensors.

The basic concept is to associate the given sensory inputs with some decision outputs. The following algorithm for pixel level image fusion using Fuzzy Logic illustrate the process of defining membership functions and rules for the image fusion process using FIS (Fuzzy Inference System) editor of Fuzzy Logic toolbox in MATLAB. The process flow of the process is as follows:

The fuzzy system considered here is Mamdani type and its' FIS model in MATLAB is shown in Fig. 6. The Mamdani rule base is a crisp model of a system, i.e. it takes crisp inputs and produces crisp outputs. It does this with the use of user-defined fuzzy rules on user-defined fuzzy variables.

Visual and Thermal Image Fusion for UAV Based Target Tracking 315

Fig. 8. Rule Base of FIS

Fig. 9. Surface View of FIS

Fig. 6. FIS Model

Fig. 7. FIS Membership Function

Fig. 8. Rule Base of FIS

Fig. 6. FIS Model

Fig. 7. FIS Membership Function

Fig. 9. Surface View of FIS

Visual and Thermal Image Fusion for UAV Based Target Tracking 317

Fig. 10. Process Flow

The idea behind using a Mamdani rule base to model crisp system behavior is that the rules for many systems can be easily described in terms of fuzzy variables. Thus we can effectively model a complex non-linear system, with common-sense rules on fuzzy variables.

The operation of the Mamdani rule base can be broken down into four parts: 1) mapping each of the crisp inputs into a fuzzy variable (fuzzification), 2) determining the output of each rule given its fuzzy antecedents, 3) determining the aggregate output of all of the fuzzy rules; 4) mapping the fuzzy output to crisp output (defuzzification).

The membership function used is Guassian as shown in Fig. 7. The rule base is formed for the fusion process as shown in Fig. 8. The surface view of the two inputs and one output is represented in Fig. 9.

## **5. Target tracking of objects**

Automatic detection and tracking of interested targets from a sequence of images obtained from a reconnaissance platform is an interesting area of research for defence related application. The video images are obtained from an unmanned aerial vehicle (UAV) with on-board guidance and navigation system. The aircraft carries a multispectral camera which acquires images of the territory and sends the information to a ground control station (GCS) in real time.

During flight, the pilot in the ground control station may identify a region of interest as a target. This identification can be click and target type or an intelligent perception type.

The target which appears on a small window could be tracked by engaging track mode. Optical flow is an approximation of the local image motion based upon local derivatives in a given sequence of images. That is, in 2D it specifies how much each image pixel moves between adjacent images while in 3D in specifies how much each volume voxel moves between adjacent volumes. The 2D image sequences used are formed under perspective projection via the relative motion of a camera and scene objects.

The differential methods for determining the optical flow are Lucas-Kanade, Horn-Schunck, Buxton-Buxton, Black-Jepson and variational methods.

**Lucas–Kanade method** is a widely used differential method for optical flow estimation developed by Bruce D. Lucas and Takeo Kanade. The assumption in this method is that the flow is essentially constant in a local neighbourhood of the pixel under consideration, and solves the basic optical flow equations for all the pixels in that neighbourhood, by the least squares criterion.