J Lévêque

T Lubin

S Mezani

A Rezzoug

Eddy current in a rotating cylinder in a static field by a stochastic method

This paper deals with the calculation of eddy current in copper cylinder. This cylinder rotates in an imposed static magnetic field. The electromagnetic problem is solved in two-dimension by considering transient motion. Two methods for eddy current computation are compared: stochastic method and classical Finite Element Method. The main goal of this paper is to compare these methods.

Introduction

In this paper we present the calculation of eddy current in a copper cylinder when it moves in an applied static magnetic field. This is a very important problem for thermal treatment of metal [1][2]. For heating any metal by induction field, we have two solutions. The first one is to submit the conductor material to a low value (about 0.7T) and high frequency variable magnetic field. The second method consists in the achievement of a high magnetic field (several Tesla), by using a superconducting magnet for example, in which the material can rotate. So, the control of the body temperature is obtained by adjusting the rotating speed and/or the value of the magnetic field. There are some studies on this kind of heating. In this paper we propose a new numerical method to calculate eddy current distribution in the cylinder. The aim of this paper is to present a stochastic method to calculate induced current in a moving material. This method has some advantage like calculation in only one point or in some points and not in the whole domain. We present these advantages on a 2D transient motion problem with a constant applied magnetic field. The results obtained with the proposed method are compared with those obtained by the classical finite element method.

Equation to be solve

First, we present the equation to solve [3] in order to find the induced current in a cylinder which rotates in a uniform magnetic field. The considered problem is solved in twodimension (2D). From this assumption, the flux density presents only two components in the x-y plane and the potential vector has only one component, this one being orthogonal to this plane (z-component). The induced current density in the cylinder is given by ( )

B v e J × + = e σ (1)
where σ e is the electrical conductivity of the cylinder.

The relation between the magnetic vector potential and the current density is

J A 0 µ - = ∆ (2)
Substituting ( 1) into (2), we obtain ( )

B v e A 0 × + - = ∆ e σ µ (3)
The relation between the electric field and the magnetic vector potential is defined as:

) (V t grad A e - ∂ ∂ - = (4)
Equation ( 4) is then replaced in (3) to give:

B) v ( A A 0 0 × - ∂ ∂ + = ∆ e e t σ µ σ µ (5)
Now, we have to calculate the term (v B) × . To do that, we define the position of a point M and the components of the vector speed in cylindrical coordinates system:

v (7) θ & or ( ) dt dθ
is the angular speed of the copper cylinder Knowing that the magnetic field is given by:

0 A A ) A ( curl B x y ∂ ∂ - ∂ ∂ = = (8)
The product (v B) × is developed and gives:

θ θ θ θ cos A sin A 0 0 B v & & r y r x ∂ ∂ - ∂ ∂ = × (9)
We now have to solve the following equation in the cylinder and in the air regions:

) A A ( A A 0 0 θ θ σ µ σ µ & & x y y x t e e ∂ ∂ + ∂ ∂ - + ∂ ∂ + = ∆ ( 10
)
where σ=0 in the air region.

Boundary condition

Now, we need to set the boundary conditions. The rotating cylinder is placed in a constant magnetic field, as shown in figure n°1.

Fig. 1. Rotating cylinder in a constant magnetic field

The whole considered domain, air and rotating cylinder regions, is a square. On two sides of the domain, we have Neumann's condition and on the other sides, the value of the magnetic vector potential is imposed, either a constant or its opposite. With these boundary conditions, we have a constant and uniform magnetic field created by external sources. In the study, the boundaries are placed far enough from the copper cylinder to not perturb the solution. So, we now have the following problem to solve: We have also a condition on the movement of the cylinder. The speed increase linearly with time from 0 to 3 rd.s -1 during 100 ms. So, the initial value of the speed is null.

             Γ + = Γ = Γ = ∂ ∂ = ∂ ∂ + ∂ ∂ + ∂ ∂ + = ∆ 2 1 0 0 A - A 0 n A 0 ) A A ( A A
In this first part, the system of equations to describe the problem has been developed. In the two following parts, we present first the stochastic method allowing to solve the problem and then the obtained results.

Calculation by stochastic method

Stochastic methods are used to solve various kinds of problem but are very rarely used in electrical engineering. To show the relevance of such method, we applied it to solve our problem [4][5][START_REF] Davey | Proc. COMPUMAG 2003[END_REF][START_REF] Netter | Proc. EUCAS 2003[END_REF].

Stochastic resolution of PDE

The aim of this part is to explain how to solve equation (11). We begin by a short presentation of Monte-Carlo method to describe the resolution of PDE by stochastic process. In a second step, we present the resolution of our problem thanks to the Brownian motion.

Fixed random walk

The resolution of a stationary problem by a Monte Carlo method consists in generating a set of fixed random walks starting from a point of calculation. A virtual particle randomly moves on a mesh. A random walk stops when the particle reaches a border with a Dirichlet's condition. The particle retrogresses when it meets a border with a Newman's condition. The probabilities of displacement are given by the discretization of the partial derivative equation as explain in the following example.

We will generalise the construction of a random walk to solve a time dependent problem. We take for example a very simple one-dimensional diffusion equation:

2 2 x U t U ∂ ∂ = ∂ ∂ (12)
The first step consists of discretizing the PDE. We notice that x = i.∆x and t = n. ∆t, where ∆x and ∆t are the spatial and temporal steps.

² x U U 2 U ² x U ² n , 1 i n , i n , 1 i ∆ + - ≈ ∂ ∂ - + (13) and t U U t u 1 n , i n , i ∆ - ≈ ∂ ∂ - (14) Boundary Air A=constant A=-constant 0 n A = ∂ ∂ 0 n A = ∂ ∂ Rotating cylinder Γ n Γ d1 Γ d2
Substituting (13), ( 14) and ( 15) into (12), we obtain the following equation ( )

1 , , , 1 , , 1 2 -
- + - = + - n i n i n i n i n i U U U U U α (15) α + α + α + + α + = - - + 2 U . 2 U 2 U u 1 n , i n , 1 i n , 1 i n , i (16) with t ² x ∆ ∆ = α (17)
From this equation, we can define the following probabilities:

α + = = - + 2 1 P P x x (18) α α + = -2 t P (19)
Thus, we have:

1 p p p t x x = + + + - . Fig. 2. Spatio temporal random walk. 1 1 ( , ) N n n U Xo To U N = = ∑
So, we define the probabilities which will be used to build the random walks.

1 n , i t n , 1 i x n , 1 i x n , i U p U p U p U - + + - - + + = , (20)
As shown in Fig. 2, the Monte Carlo Markov Chain Method consists in computing a set of random walks starting from the same point P(Xo=i o .∆x, To=n o .∆t) for which we need to approximate the value of U. We consider that dx and dt values are quite chosen for the numerical calculation. To do that let us consider the value of α. If it is very small, whereas the value of the probability P x+ and P x-tends to 0.5, P t-tends to 0 and the convergence process is time consuming. If α is very large, the space probabilities tend toward 0 while P t-tends to 1. A reasonable choice for the value of α is roughly1 which gives an equiprobability for the three walks.

Brownian motion

We saw that the probability of displacement on the mesh comes from the discretization by finite difference of the equation ( 12). Knowing that when the spatial step tends to zero, the random walk tends to a Brownian motion, we use this property to build now a random walk without mesh.

The general expression of a diffusion processes is given by:

∑ ∑ = = ∂ ∂ + ∂ ∂ ∂ = = ∂ ∂ d j i d i i i j i ij x U x x x U x S t x Lu t U 1 , 1
2

) ( ) ( 2
) x ( ). x ( ) x ( S T ij σ σ = (23)
X(t) is a stochastic variable which defines the position of a particle at a time t. At t=0 this particle is at the position where we need to calculate the solution of the equation.

t t X t
)) ( , ( µ is the expectation of the stochastic variable.

)) ( , ( 2 t X t σ is the variance of the stochastic variable. B t is the Brownian motion. It is a stochastic variable whose increment B t+δt -B t follows a normal distribution law N(0, δt).

With this method, a stationary problem is calculated like the limit of a time dependent problem.

We now present a simple example to illustrate the calculation process. We consider the Laplace equation in 2D:

0 y U x U 2 2 2 2 = ∂ ∂ β + ∂ ∂ α (24)
This equation could be rewritten to correspond to the form defined by (21):

0 y U 2 x U 2 2 1 2 2 2 2 =         ∂ ∂ β + ∂ ∂ α (25)
The parameters of the stochastic variable X(t) are now the following :

U(x=0,t) U(Xo,To) U(x=L,t) U(x,t=0) U 2 U n U 1 σ σ =         β α = T ij 2 0 0 2 ) x ( S (26)
where

        β α = σ 2 0 0 2 (27) and 0 = µ .

How many Brownian motion do we have to take?

It is obvious that we cannot take only one Brownian motion to estimate the solution at the point P(X o ,T o ). The question is how many Brownian motion we have to take?

Let us define the stochastic variable U N by:

∑ = = = N n 1 n n N U N 1 U . ( 28
)
The stochastic variable associated to the n th Brownian

displacement is U n . n U ~
is the simulation result for the n th random process starting from the point P(xo,to). The central limit theorem indicates that U N follows a Gaussian stochastic law with U S the mathematical expectation and with We obtain not only one estimate of the solution but also an estimate of the error by the value of the standard deviation S σ ~.

Application to a rotating cylinder

In this part, we use Brownian motion to solve the problem of a copper cylinder rotating in magnetic field as represented in figure 1. To solve equation (11), we first rewrite it in order to separate time derivative term:

) A A ( A A 0 0 θ θ σ µ σ µ & & x y y x t e e ∂ ∂ + ∂ ∂ - ∆ = ∂ ∂ (30) 2 2 2 2 0 1 A A A A A y x t x y x y θ µ σ     ∂ ∂ ∂ ∂ ∂ = + - +     ∂ ∂ ∂ ∂ ∂     & (31)
An elementary displacement of the Brownian motion is given by the previously described equation:

t dB )) t ( X , t ( dt )) t ( X , t ( ) t ( dX σ + µ = (32)
In this case, the parameters are the following characteristics:

0 0 2 0 ( , ( )) 2 0 e e t X t µ σ σ µ σ       =         (33)
and

      θ - = µ x y )) t ( X , t ( & (34)
The main characteristics of the studied problems are summarized in

      + σ = σ (35) d c / ) 2 y 2 x ( e ) x (       + θ = θ & & , (36)
where c and d are parameters adjusting the curves. The variation of the electrical conductivity versus space is shown on figure 3. If the c term is too weak, the variation of the speed is not sharp and the description of the phenomenon is not good. We do not see influence of this term on the computation speed and the precision of the calculus.

-0,2 -0,16 -0,12 -0,08 -0,04 0 0,04 0,08 0,12 0,16

Results

In order to study the effectiveness of the stochastic method, the results have been compared with those obtained from 2D transient finite element simulations. Figure 4 and figure 5 show the flux lines obtained respectively by stochastic method and finite element simulation for the same condition. We can observe that the results are in good agreement but very noisy with stochastic method. Figure 6 shows the time evolution of the potential vector at a given point (x= 0.8, y= 0.8) on the cylinder surface obtained by finite element simulation and Brownian motion. We note again a good agreement between the results obtained by the two methods. The main advantage of the stochastic method is in the choice of the number of Brownian motions and consequently in the precision of the calculus because we have an estimator of the results. If we need very short computer time consuming for the calculus, we can associate the stochastic method with a numeric filter to avoid the noisy part of the numerical results. Figure 7 and 8 show a comparison between two computations results of the magnetic vector potential, one with a low number of Brownian motions and another with finite element method. We also present on Fig. 7 the curve obtained thanks to a noise filtering. It can be seen in Fig. 7 and 8 that the results are much closed for FEM and stochastic method. We show on Fig. 9 a comparison between stochastic method and FEM on a line at y =0.04 and we a good agreement between the two methods.

Conclusion

In this paper, we have solved the problem of induced current in a moving copper cylinder placed in a uniform magnetic field. We have shown that the stochastic method gives results very close to those obtained by finite element simulations. Each method has its own advantages and drawbacks.

Among the advantages of the finite element method, we can quote its popularity and the large number of software packages, some of them being free. One of the problems of this method is that we need to calculate the solution in the whole space and for transient problem, at each space time. Another problem is that we have, when thin meshing is required, a large system to solve with the risk numerical instability.

Among the advantages of stochastic method to solve partial differential equation, we can notice that:

• It is very easy to program it.

• It is not necessary to mesh all space and, so, to use a lot of computer memory. • It is possible to calculate solution in only one point and not in the whole space. • It is also possible to calculate solution at one time and not at each space time.

Because of these advantages, the method can be easily included in a larger design program to study problems of eddy currents and new thermal heating systems.

  correspond to the current Brownian motion and subscript s correspond to the final value of the stochastic variable As N increases, the diameter of the 95 % confidence interval

Fig. 3 .

 3 Fig. 3. Electrical conductivity versus space for c=0.01 and d=5.3

y

  Flux line at t=100 ms with copper

Fig. 4 .Fig. 5 .

 45 Fig. 4. Flux lines obtained with Brownian motion method

Fig. 6 .

 6 Fig. 6. Comparison of the two methods at one point on the border of the cylinder

Fig. 7 .Fig. 8 .Fig. 9 .

 789 Fig. 7. Comparison of two stochastic methods at one point on the border of the cylinder

table 1 Table 1 .

 11 Characteristic of the problem

										Copper	Air
	Angular speed	0 to 3 rd/s	0
	Electrical		37 10 6 S/m	0
	conductivity
	Size					cylinder radius of	Square of 0.16 m
										0.08 m
	Step size		0.004 m
	Simulation time	0.1 s
	Time step	10 -7 s
	The value of the magnetic vector potential is fixed at 0.1
	Wb/m on the boundary Γ d2 and at -0.1 Wb/m on the boundary Γ d1 . These values give a constant magnetic flux density of 1.25
	T in the whole domain when the copper cylinder doesn't
	move.
	Electrical conductivity σ e and angular speed are null in the air
	region and have a non null value in the copper cylinder region.
	To overcome the numerical problem due to any discontinuity
	of the parameters in equation (11), functions	( σ	t	,	X	(	t	))	and
	µ	(	t	,	X	(	t	))	must be continuous. For that, we adopt a sigmoid
	representation of speed and conductivity

& & & &


Yu-Yin Chen

Dejan Vasic

François Costa

Wen-Jong Wu

Chih-Kung Lee

A self-powered switching circuit for piezoelectric energy harvesting with velocity control

The rapid development of low power consumption electronics and the possibility of harvesting energy from environmental sources can make totally autonomous wireless devices. Using piezoelectric materials to convert the mechanical energy into electrical energy for batteries of wireless devices in order to extend the lifetime is the focus in many researches in the recent years. It is important and efficient to improve the energy harvesting by designing an optimal interface between piezoelectric device and the load.

In this paper, a self-powered piezoelectric energy harvesting device is proposed based on the velocity control synchronized switching technique (V-SSHI). Comparing to the standard full bridge rectifier technique, the synchronized switching harvesting on inductor (SSHI) technique can highly improve harvesting efficiency.

However, in real applications when the energy harvesting device is associated with wireless sensor network (WSN), the SSHI technique needs to be implemented and requires being self-powered. The conventional technique to implement self-powered SSHI is to use bipolar transistors as voltage peak detector. In this paper, a new self-powered device is proposed, using velocity control to switch the MOSFET more accurately than in the conventional technique. The concept of design and the theoretical analysis are presented in detail. Experimental results are examined.

Introduction

Nowadays with improvement and rapid growth of lowpower electronics, it is possible to supply portable devices such as mobile phone, MP3 player, wireless sensors and human or animal detecting devices, with harvest energy from ambient. Among these low-power devices, wireless sensor network (WSN) is one of the most important and valuable applications which is highly investigated. Wireless sensor networks can be used to monitor the health of structures, environment, wild animals, tire pressure of running cars, etc. In most of WSN applications, the devices are far from the power line or the devices need to be embedded into the structure to monitor. So, it is hard to use power line to transmit energy to device; battery is the only conventional solution. However, there are lots of disadvantages with using batteries. The major problem is the lifetime: using a 3 V battery a WSN can only be operated for 1 or 2 years. The batteries cannot be a permanent energy supply for WSN. A WSN with selfpowered system can be operated for a longer time without replacing the battery. Harvesting the ambient energy close to the sensor nodes of the WSN is the most likely and suitable solution to extend the its lifetime [1][2][3][4][START_REF] Paradiso | Pervasive Computing[END_REF][START_REF] Roundy | [END_REF][7][8][START_REF] Ottman | Power Electronics[END_REF][10]. There are various methods to convert mechanical ambient energy into electrical energy. According to the conversion mechanism, the power harvesting devices can be classified into three types: piezoelectric, electromagnetic and electrostatic [11]. Among these different types, piezoelectric Send offprint requests to:

a Present address: vasic@satie.ens-cachan.fr transducers have been of great promise and interest.

The electric energy issued from a piezoelectric device is alternating voltage and current, so the simplest interfacing circuit that can be used is a standard DC approach (full bridge diodes rectifier). From studies [12], the efficiency of the standard DC approach is not optimal. There are lots of studies devoted to efforts to design and optimizing the interfacing circuits in order to improve the efficiency and maximizing the output power. A DC-DC converter is studied to match the impedance of the resistor load [START_REF] Ottman | Power Electronics[END_REF]13]. The Synchronized Switch Harvesting on Inductor (SSHI) technique is a very successful and efficient technique to boost the output power from piezoelectrics [14][15][16][17]. This approach is derived from a semi-passive damping technique: Synchronized Switch Damping on Inductor (SSDI) [START_REF] Richard | Proc. SPIE: Smart Structures and Materials 1999: Passive Damping and Isolation[END_REF][START_REF] Richard | Proc. SPIE: Smart Structures and Materials 2000: Damping and Isolation[END_REF]. The SSHI technique consists in adding up a nonlinear switching. This nonlinear process increases the output voltage of the piezoelectric elements that increase the output power. The switching device is triggered at the zero crossing of velocity. In order to realize the synchronized switching technique in real applications without external power source to supply the system, many researches present self-powered supply system for piezoelectric energy harvesting devices [4,[START_REF] Ng | [END_REF][21][START_REF] Junrui | Information and Automation[END_REF]. The design concept of self-powered system proposed by [21] is shown in Figure 1. This conventional self-powered system works thank to a peak voltage detector to control the switching time for SSHI technique. However, the energy supplying the peak detector and the switching control is drawn from the piezoelectric device. The energy losses in the circuit can be accurately controlled by circuit design; the larger excitation levels leading to relatively smaller losses. As the conventional self-powered system uses peak detector, there is always a phase lag between the peak voltage and the actual switching time. Moreover, the phase lag for large excitation level is less than for a small one [START_REF] Junrui | Information and Automation[END_REF].

In this paper, we present a new self-powered piezoelectric energy harvesting system using velocity control SSHI technique, called V-SSHI. The schematic design concept of self-powered velocity control SSHI is shown in Figure 2.

The SSHI used herein is in series type (inductor and switch are in series with piezoelectric patch). Comparing to the conventional design concept, the piezoelectric material is separated into three parts. The main part is dedicated to harvest ambient vibration energy. The second small part is designed to supply energy for switching MOSFET and the last small part is designed for velocity control and for switching on the optimal time. There are two major advantages of this new technique: (1) theoretically, there isn't phase lag by using velocity control signal to determine the switching time; (2) the supply energy for the switching driver can be designed and optimized by the size of the piezoelectric material.

The energy flow chart of the conventional self-powered technique and the self-powered V-SSHI technique is shown in Figure 3. There is a common path for the main stream of energy and for the supply of the self-switching system.

In the V-SSHI technique, they are three energy paths. The energy supplying the self-switching system and the velocity control patches can be designed optimally. The theo-

F = K E P x + αV I = α ẋ -C 0 V (1)
Where:

K E P = c E 11 w.t l (2)
α = e 31 w (3)
C 0 = ε S 33 w.l t (4)
l, w, t denote length, width and thickness, and l, w are much larger than the thickness t. K E P is short-circuit stiffness, c E 11 is elastic constant in 1 direction, α is forcevoltage coupling factor, e 31 is piezoelectric constant, C 0 is clamped capacitance of piezoelectric and ε S 33 is permittivity constant. The governing equations of this model are expressed by [START_REF] Paradiso | Pervasive Computing[END_REF]: D is the damping ratio of the structure, K E is the equivalent stiffness of the structure and piezoelectric materials.

mẍ + D ẋ + K E x = F E -αV I = α ẋ -C 0 V (5)

Standard DC technique

Before talking about the models of the V-SSHI technique, the standard DC technique is proposed to be a reference.

The schematic diagram of piezoelectric energy harvesting transducer with full bridge rectifier connected to a resistor is shown in Figure 6. This is called standard DC technique.

Assuming the structure is excited at the mechanical resonance frequency, the vibration velocity can be modeled by a current source, noted I eq . Figure 6 shows also the key waveforms of the standard DC approach. When the absolute voltage value of the piezoelectric patch V P is less than voltage V C , the diode bridge is in open-circuit. The diodes conduct and piezoelements charge the load only when V P reaches load voltage V C .

According to the relation between force F E , displacement

x and velocity ẋ (5), the output voltage V C can be a function of excitation force FE as [START_REF] Roundy | [END_REF] shown. The output power in the standard DC approach can be calculated using simple equation V 2 C /R as equation (7) shows.

V C = 2αR FE (2RC0ω0+π)D (6) P = V 2 C R = 4α 2 R F 2 E (2RC0ω0+π) 2 D 2 (7)

Self-powered V-SSHI technique

According to the self-powered V-SSHI concept presented in Figure 2, the model can be easily separated into three parts. The details are presented hereunder.

Main patch for SSHI

The main patch concerned with our new concept is designed to act like a classic SSHI technique. The schematic diagram of a SSHI technique is shown in Figure 7. The fundamental concept of SSHI is to use an inductor L and achieving a LC 0 resonance between piezo-element and L.

Through LC 0 resonance and switches to confine the current flow, more power can be harvested from the piezoelement. Assuming the structure is excited at the mechanical resonance frequency, the excitation source can be modeled as a current source I eq . Figure 8 shows the waveform of SSHI-series including voltage across piezoelectric V P , current source I eq and displacement x.

Current I R which flows through R from time period T 1 to T 3 is equal to the resonant current I L . Considering the 8). The output power from the transducer can be calculated using the simple formula V 2 C /R as shown in equation [START_REF] Ottman | Power Electronics[END_REF].

V C = 2αR(1+exp( -π 2Q I )) FE (2RC0ω0(1+exp( -π 2Q I )+π(1-exp( -π 2Q I ))D (8)
P = V 2 C R = 4α 2 R(1+exp( -π 2Q I ) 2 F 2 E (2RC0ω0(1+exp( -π 2Q I )+π(1-exp( -π 2Q I ))) 2 D 2 (9)
Considering the derivate of equation ( 9) relative to load resistance R, the optimal value R opt of R can be calculated, it is expressed by (10).

R opt = π(1-exp( -π 2Q I )) 2C0ω0(1+exp( -π 2Q I )) (10)
Substituting ( 10) into (9), the maximum power value P max can be obtained, it is expressed by (11). (13). Resistor R eq is the equivalent load between +V CC and -V CC .

P max = α 2 (1+exp( -π 2Q I ) F 2 E 2πC0ω0(1-exp( -π 2Q I )D 2 (11)
V CC = αReq FE (C0Reqω0+π)D (12) P = V 2 CC Req = α 2 Req F 2 E (C0Reqω0+π) 2 D 2 (13)

Sensor patch for Velocity Control

The third patch is designed for generating the velocity control signal. The equivalent circuit is shown in Figure 10. A low value resistor R C is connected in parallel with the patch to sense the mechanical current I eq . A passive low-pass filter is used to reduce the high frequency noise.

When SSHI works, the high frequency noise of the velocity signal is very large, so it is impossible to apply directly voltage V P to comparator. The current sensing resistor used herein must be small enough to avoid the effect of the piezoelectric capacitance. The low-pass filter should be carefully designed to guarantee there is no phase lag for the considered frequency. The key waveforms are also shown in Figure 10. The blue line V P is the opencircuit waveform of the piezoelectric patch and the red line V S is the velocity control signal which is in phase with I eq . There is 90 degree phase lag inherently between V S and V P . When the circuit switches by velocity control, the switching time can be accurate; the current is always in phase with voltage when SSHI works. The power output from piezoelectric can be always positive.

3 Experimental results and discussion

Experimental Setup

The experimental structure under test is a cantilever steel beam. Three 31-type PZT-QA piezoelectric ceramic patches provide by the company ELECERAM were bonded on the beam. Table 2 gives the dimensions of the beam and the patches. Figure 11 shows the experimental setup and the self-powered V-SSHI circuit diagram. Figure 12 shows a picture of the experimental setup. In the experimental setup, the SSHI circuit part is a little bit different from the one of Figure 7, but it works identically. The four diodes act like a full bridge rectifier to confine the current flow and the inductor is in series with P 1 .

The cantilever beam is excited by a vibration shaker (LDS-V406). Three piezoelectric patches (P 1 to P 3 ) are bounded close to the fixed end. An accelerometer (PCB-353B03) is situated at the fixed end to measure acceleration; a laser vibrometer (LK-G32) measures the displacement at the free end. P 1 is the main patch for harvesting power. The circuit connected to P 1 for SSHI is composed of several parts: an inductor L for LC 0 resonance to enhance the power; four schottky diodes (D 1 to D 4 ) for confining the current flow; the load composed of a resistor and a capacitor; NMOS (2N7002) and PMOS (NDS0610) pair for positive and negative switching. There are two stages for switching:

-when velocity crosses zero from negative to positive, voltage P 1 is at the maximal positive value, the NMOS is switched at this time and the SSHI process will occur through the path L-D 1 -LOAD-D 3 -NMOS.

-the negative stage works with the same logic through the path L-D 2 -LOAD-D 4 -PMOS.

Patch P 2 is connected to circuit composed of two schottky diodes and three capacitors to generate the positive voltages +V CC and -V CC for supplying comparator (TLV3701).

Voltage V CC should be larger than 2.5V to make sure that the comparator fully works to drive MOSFET. The comparator chosen here is a nano-power comparator from TI and the supplying current is only 560nA/per channel. This nano-power comparator is very easy to drive and suitable for low power circuit design. Patch P 3 is designed for velocity control. It is connected to a current sensing resistor followed by a passive low-pass filter. The resistor used herein is small enough to make sure there is no phase lag. The velocity signal noise (sine wave ideally) is attenuate by the low-pass filter. A comparator is used in order to generate the switching signal (square wave ideally) to drive NMOS and PMOS. The low-pass filter is designed to reduce high-frequency noise without phase lag. Figure 13 shows the experimental waveforms of the self-powered V-SSHI device.

Experimental results

Figure 13 shows three waveforms:

-black line V P is the waveform of the piezoelectric patch

P 1 ,
-blue line V Cout is the output waveform of the comparator, -red line V S is the velocity control signal after the lowpass filter.

Although there is still some high frequency noise in the velocity control signal V S , the comparator work still well;

it is a tradeoff between reducing noise and phase lag. Observing waveform V P , we can note that the switching time occurs almost at the peak value of the voltage.

The model parameters, identified by measurements, are

given in table 3. The experimental and theoretical results of output power are shown in Figure 14. All experimental data are acquired for the same acceleration (a = 2.5m/s 2 ). The theoretical curves for standard DC and standard SSHI are drawn from equations ( 5), (7) and parameters in table 2. The standard DC experiment (measured using Figure 6 circuit) and "SSHI-Experiment"

(measured by power supply and function generator using The conventional self-powered technique proposed by [START_REF] Junrui | Information and Automation[END_REF] is the line with green points. Experimental results show that the maximum output power of self-powered V-SSHI is higher than the conventional technique, essentially due to the efficient phase control.

Conclusion

In this study, a self-powered V-SSHI piezoelectric energy harvesting is proposed, this is a new design concept. Based on the outstanding performance of SSHI technique, the self-powered V-SSHI circuit is fully self-powered, requiring no external power supply and though the velocity control, the switching time can be more accurate than with state-of-the-art techniques. The performance of the conventional self-powered circuit is close to the theoretical values of the SSHI; however, it requires an excitation level high enough to work properly [START_REF] Junrui | Information and Automation[END_REF]. In the self-powered V-SSHI technique, the excitation level doesnt influence the performance and when the supply voltage of the comparator is larger than 2.5V, the whole circuit fully works. The experimental results show better performance and lead to a gain of around 200% compared to the standard DC approach. Of course, the V-SSHI output power is lower than the theoretical SSHI, because the energy is split to supply the auxiliary self-powered circuit.

The architecture proposed in this study is more benefi-cial and represents a new step of the design concept. This circuit is easily used in real applications and may be combined with wireless sensor networks.

Fig. 1 .Fig. 2 .Fig. 3 .

 123 Fig. 1. Schematic design concept of conventional self-switched system

2

 2 Theoretical analysis The mechanical and electrical models of standard DC technique and V-SSHI technique are analyzed based on the following model. The analysis of piezo-element in this study is the same as for the classical model [17]. Considering the piezoelectric patches bounded on a cantilever beam, they can be regarded as a simple energy harvesting device. The dimensions of piezo-elements, schematic and coordinate directions are shown in Figure 4. When the cantilever vibrates, at the first mode, the force which acts on the piezoelectric patches can be simplified to a 1-D model and regarded as a force F acting on the lateral surface as shown in Figure 4. The piezoelectric equations (mechanical and electrical parts) can be expressed as follows (1):

Fig. 4 .Fig. 5 .

 45 Fig. 4. Schematic of piezoelectric harvesting cantilever beam

Figure 4 .

 4 Figure 4. The electrical and mechanical quantities are defined in table 1.

Fig. 6 .

 6 Fig. 6. The schematic diagram of the Standard DC approach and waveforms

Fig. 7 .Fig. 8 .

 78 Fig. 7. Schematic diagram of SSHI piezoelectric energy harvesting device with full bridge rectifier to a resistor load

Fig. 9 .

 9 Fig. 9. The equivalent circuit diagram of the supplying circuit

Fig. 10 .

 10 Fig. 10. The Equivalent circuit diagram of the velocity control circuit and the waveforms

Fig. 11 .Fig. 12 .Fig. 13 .

 111213 Fig. 11. Experimental setup and circuit diagram of V-SSHI device

Figure 7

 7 Figure 7 circuit) are the reference lines compared to theoretical lines; results show good agreement with predictions. Piezoelectric patch P 3 can be replaced by a smaller one. So, in the experiments of this study, the effect of patch P 3 is neglected. The experimental results (blue point) called "SSHI-Experiment" are measured on the conventional SSHI technique powered by external switching signal. The selfpowered V-SSHI technique (red point) is measured by only one patch P 1 . The output power of V-SSHI circuit is lower than the one of "SSHI-Experiment" circuit, because the energy is split to supply the auxiliary self-powered circuit.

Fig. 14 .

 14 Fig. 14. Experimental Results

Table 2 .

 2 Dimension of the electromechanical transducer

	Steel beam
	LengthxWidthxThickness	168.5mm x 94.3mm x 15mm
	First bending mode	41.4Hz
	Piezoelectric pathes (PZT-QA)
	P1	38.1mm x 16.5mm x 0.5mm
	P2	15mm x 5mm x 0.5mm
	P3	15mm x 5mm x 0.5mm

Table 3 .

 3 Measurements and model parameters

	f0	Short circuit resonance frequency	41.41Hz
	f1	Open circuit resonance frequency	41.45Hz
	ξ	Open circuit damping coefficient	0.00105
	C0	Clamped capacitance P1	25nF
	C0	Clamped capacitance P2	3.5nF
	C0	Clamped capacitance P3	3.5nF
	α	Force-Voltage coupling factor	0.00069N/V
	k 2	Electromechanical coupling factor 0.00193
	M	Mass	182g
	K E Equivalent stiffness of the source	12 320N/m
	D	Damping of the structure	0.1N/m/s


between erosion and tectonics are thought to play

a determinant role in mountains evolution. Here, we investigate the interplay in this coupling between the assumed erosion law and the crustal rheology at the margin of a collisional plateau, like the Himalaya of central Nepal.

Lithospheric deformation is calculated over a time scale of 100 kyr by a twodimensional finite elements model that incorporates the rheological layering of the crust and the main features of the convergence across the range including shortening rate and geometry of the Main Himalayan Thrust fault.

For the upper boundary condition, several surface processes were tested: a linear diffusion model and a 1D1/2 integrative model including fluvial incision along the fluvial network and hillslope erosion by landsliding. Model results and their sensitivity to the chosen combinations of erosion law and crustal properties are discussed in light of the constraining geologic and geomorphologic observations (topography, river elevation, denudation and fluvial incision rates). In contrast with the conclusions of [START_REF] Cattin | Modeling moutain building and the seismic cycle in the Himalaya of Nepal[END_REF] where a compliant, quartz-rich crustal rheology and diffusion law were required, we rather propose to use a composite quartz-diabase rheology for the crust associated with fluvial incision to account for erosion and elevation profiles across the Himalaya of central Nepal. More generally, it is proposed that, because of the interplay between the dominant denudation conditions and rheology of the crust, well documented erosion rates and processes D R A F T February 15, 2013, 9:11am D R A F T

Introduction

The land surface is a dynamic interface that results from the combination of tectonic uplift and denudation. Knowledge of the linkages and feedbacks between those two processes is essential for the understanding of the structure and the evolution of mountain belts. Because of several numerical limitations, few 3D studies with full coupling between tectonics and erosion have been conducted so far. Most of the thermomechanical finite element models are bidimensional, where erosion processes are usually reduced to a 1D process acting on the surface profile. Erosion is either resumed to a diffusion law [START_REF] Avouac | Erosion as a driving mechanism of intracontinental mountain growth[END_REF][START_REF] Cattin | Modeling moutain building and the seismic cycle in the Himalaya of Nepal[END_REF], or to a linear relation of the surface slope [START_REF] Beaumont | Himalayan tectonics explained by extrusion of a low-viscosity crustal channel coupled to focused surface denudation[END_REF], or to a simple river incision law where the mean topographic profile is represented by a river profile [Willet, 1999]. These simplified models ignore the respective role of river network and hillslopes in controling the morphology and evolution of the landscape. In this paper we investigate the influence of the assumed erosion law and of the rheological properties of the crust on crustal deformation through the use of a 2-D thermo-mechanical finite element model. We apply our approach to the Himalaya of central Nepal, one of the best documented example of active mountain belt. After a short presentation of the main features of the Himalayan orogenic belt, we describe the modeling approach and the two distinct erosion laws we want to test. We next compare the existing geomorphologic observations (including horizontal shortening, topography, river elevation, denudation and fluvial incision rates) with model-simulated landscape produced by various combinations of erosion law and crustal properties. We finally discuss the sen-

D R A F T February 15, 2013, 9:11am D R A F T X -5
sitivity of our results and try to highlight the most discriminant observations to unravel rheologic and erosional conditions prevailing in a mountain range.

Geodynamical Setting and characteristics of the fluvial network

The Himalayan belt has resulted from the ongoing collision between the Indian and Asian plates. It is one of the most active orogen of the Earth, characterized by a steep topographic front from the 5000m elevated Tibetan Plateau down to the Gangetic plain.

This topographic step goes through four major morphotectonic domains: the rugged South Tibetan plateau, the High Himalaya (HH) with deep gorges and ∼ 8000m summits, the lower relief of the Lesser Himalaya (LH), and the frontal low elevation relief of the Siwaliks Hills. The Himalayan range is affected by an intense ongoing seismicity [e.g. [START_REF] Ni | Seismotectonics of the Himalayan collision zone: Geometry of the underthrusting Indian plate beneath the Himalaya[END_REF][START_REF] Pandey | Interseismic strain accumulation on the Himalayan crustal ramp (Nepal)[END_REF], and displays abundant records of active deformation [START_REF] Nakata | Geomorphic history and crustal movements of the foothills of the Himalayas[END_REF]. The long term shortening rate across the range is ∼ 20 mm.yr -1 [Lyon [START_REF] Lyon-Caen | Gravity anomalies, flexure of the indian plate, and the structure support and evolution of the Himalaya and Ganga basin[END_REF][START_REF] Armijo | Quaternary extension in southern Tibet[END_REF]. During the Holocene, this convergence has been mostly transferred to the southernmost thrust or Main Frontal Thrust (MFT) [START_REF] Lavé | Active folding of fluvial terraces across the Siwaliks Hills, Himalayas of central Nepal[END_REF]. This frontal fault branches on the Main Himalayan Thrust (MHT) rooting at 30-40km deep beneath the South Tibet [Zhao et al., 1993], and displaying a ramp and flat geometry beneath the HH and LH domains [START_REF] Schelling | Thrust tectonics, crustal shortening and the structure of the Far-Eastern Nepal Himalaya[END_REF][START_REF] Lavé | Fluvial incision and tectonic uplift across the Himalayas of central Nepal[END_REF]. Several major north south rivers drain the Himalayas of Nepal intense rainfalls is deeply entrenched in the topography and actively participate to the denudation of the orogen (figure 1). Transhimalayan rivers and fluvial terraces profiles [Lavé andAvouac, 2000, 2001] or thermochronologic data [START_REF] Burbank | Decoupling of erosion and precipitation in the Himalayas[END_REF] suggest in addition that erosion is maximal across the Siwaliks and the HH, lower in the LH and minimal in South Tibet.

Modeling approach

Following [START_REF] Cattin | Modeling moutain building and the seismic cycle in the Himalaya of Nepal[END_REF], our model is based on a 700 km long N18 cross section perpendicular to the range, from the Gangetic Plain to the Tibetan Plateau (see figure 1 for the location). We use a 2-D finite element model, ADELI [Hassani et al., 1997] that accounts for the mechanical layering of the crust, the non-newtonian rheology of rocks, and their dependencies on temperature and pressure. Three lithological layers are distinguished: the upper and lower crusts, and the upper mantle. Each layer has specific mechanical properties. We use the empirical rheological equations and laboratory-derived material properties for quartz, diabase and dry olivine (see supplemented electronic document 1). Those rheologies are dependent on temperature which is prescribed as an initial condition [START_REF] Henry | Kinematic, thermal and petrological model of the Himalayas: Constraints related to metamorphism within the underthrust Indian crust and topographic elevation[END_REF] and do not evolve during the simulation, considering the typical duration of ∼100 kyr (see Appendix 1 (suppl. Electr. Material) for the detail of time scenario). For the different runs, this duration was usually sufficient to reach a stabilized topographic profile, i.e. an average equilibrium between uplift and erosion (in

D R A F T February 15, 2013, 9:11am D R A F T X -7
a Lagrangian sense). The principal geometric characteristics of our model are similar to Cattin and Avouac's model [2000] (figure 2). The geometry of the Moho beneath South Tibet is derived from INDEPTH seismic profile [Zhao et al., 1993] and from gravity data [START_REF] Cattin | Gravity anomalies, crustal structure and thermo-mechanical support of the Himalaya of central Nepal[END_REF]. The boundary conditions applied to the system are constrained by geodynamics data: we apply a 20 mm.yr -1 horizontal velocity on the northern vertical face for depth above 40 km and leave vertical velocities free for the other vertical faces.

The structure is supported by hydrostatic pressure at its base which allows isostatic compensation and thus realize a coupling between uplift and denudation. Due to the duration of our simulations we are not taking into account the seismic cycle and the slip on a low friction MHT is considered as continuous. Our main goal is to test the importance of the upper boundary condition imposed with the introduction of surface processes on the evolution of the system.

Surface processes

We distinguish two domains in term of surface processes: the foreland, south of the MFT

∂ h ∂t x = κ ∂ 2 h ∂x 2 x , (1)
where κ is the mass diffusivity coefficient, h is the elevation of the topography and x is the distance from the MFT. The diffusion model can be applied to describe the shape and evolution of small scale topographic feature. However, it hardly applies at the scale Whipple and Tucker, 2002], in an attempt to develop a simple approach, we have used a detachment limited relation that provides satisfying first order results in the Subhimalaya [START_REF] Lavé | Fluvial incision and tectonic uplift across the Himalayas of central Nepal[END_REF]. This relation states that bedrock incision rate of a river is proportional to the fluvial shear stress τ in excess of some threshold τ c , computed from the river slope S x and flood discharge derived from rainfall profile (figure 2) [START_REF] Lavé | Fluvial incision and tectonic uplift across the Himalayas of central Nepal[END_REF]:

∂h riv ∂t x = K x (τ -τ c ) (2)
with,

τ = k 1 ( Px -P r ) γ (L(X -x)) β S x s 0 α , (3)
where K x is the erodability coefficient (figure 2) depending on rock strength, k 1 a coefficient that depends on the river network geometry, sediment size and flood distribution, L the width of the watershed, s 0 the sinuosity of the river and α, β, γ are exponents consid-

D R A F T February 15, 2013, 9:11am D R A F T X -9
ered as constant in our study area. P is the average precipitation on the watershed and P r some treshold runoff. X is the abscissa of the drainage divide. Despite their leading role, the main rivers do not account for the mean topography, which represent the pertinent variable for the upper boundary condition of mechanical modelling. The elevation profile of the transhimalayan river represents in fact the base level for the network of tributaries which are draining the whole topography, from their sources at the base of the hillslopes to their confluence with the trunk stream. At a given abscissa, the mean elevation of the topography h is therefore the sum of three contributions: (1) the elevation of the main river h riv , (2) the fluvial relief associated to the tributaries ∆h trib that we assume to be controlled by the same incision law as the main river and (3) the relief of the hillslope from the fluvial network to the crest ∆h hill . In active orogens, hillslopes are dominated by landslides [START_REF] Hovius | Sediment flux from a mountain belt derived from landslide mapping[END_REF]: we thus assume that they display a critical slope angle of repose θ c and that they react instantaneously to any local base level drop. A new formalism to integrate the fluvial relief associated to the tributaries network is proposed by Lavé [in prep.] and enable to compute at each time step the changes of the elevation of the trunk stream from equation 2 and of the changes in mean topography according to

∂ h ∂t x = ∂h riv ∂t x + K x (k 2 (P x -P r ) γ ∆h α tribx -τ c ) (4) with, ∆h trib = ∆h total -∆h hill = ∆h total - ∆l 2 tan θ c , (5)
where k 2 is a coefficient which depends like k 1 on the tributary network and the flood distribution, and ∆l the horizontal distance between the crest and the base of the hill-

D R A F T February 15, 2013, 9:11am D R A F T
slope. The different parameters of the two erosion models are detailed in Appendix 2 (supplemented electronic material).

Sensitivity to the erosion law

To investigate the influence of the denudation law, the rheological layering used is quartz (upper crust), diabase (lower crust) and wet olivine (mantle) as on figure 2. The comparison between the two tested erosion models is carried out on 4 different profiles:

mean elevation, horizontal velocity, uplift and erosion rates (figure 3). The river elevation and incision rates are specific to the fluvial incision model, and thus do not appear in the diffusive case. The two models display very similar profiles of horizontal velocity: the southward transfer of the Himalayan shortening is allowed by the low friction condition on the MHT [START_REF] Cattin | Modeling moutain building and the seismic cycle in the Himalaya of Nepal[END_REF]. In contrast, the uplift and erosion rates display distinct profiles suggesting that part of the vertical motion is controlled by deformation in the middle and lower crust in response to erosional unloading. Diffusion processes localize the highest erosion at the southern edge of Tibetan plateau (x∼ 140 km on our profile), associated with the maximum amount of slope variations along the profile. With the fluvial incision model, the erosion of the whole topography acts at a slower rate than incision of the main river: less denudation leads to a faster southward migration or gravitational spreading of the plateau and counteracts against river regressive erosion.

Finally, diffusion based landscape evolution leads to relatively smooth topographic profiles, with continuous slopes variations. By contrast using the incision law allows to preserve a clear slope transition between the LH and HH. This enhances that, in addition of being the most realistic of the two erosion laws, only the fluvial incision model provides a good

D R A F T February 15, 2013, 9:11am D R A F T X -11
agreement with all the available observations: i.e. the highest rates of erosion and incision across the HH (x∼ 100 km in Figure 3) [START_REF] Lavé | Fluvial incision and tectonic uplift across the Himalayas of central Nepal[END_REF][START_REF] Burbank | Decoupling of erosion and precipitation in the Himalayas[END_REF],

positive denudation in the LH, ≤ 1 mm.yr -1 of sedimentation in front of the MFT [START_REF] Lavé | Active folding of fluvial terraces across the Siwaliks Hills, Himalayas of central Nepal[END_REF], the sharp topographic transition between the LH and HH.

Sensitivity to the crustal rheology

In contrast with our results, [START_REF] Cattin | Modeling moutain building and the seismic cycle in the Himalaya of Nepal[END_REF] obtained with a linear diffusion model a good fit to the geophysical observations data including the estimated pattern of river incision. However [START_REF] Cattin | Modeling moutain building and the seismic cycle in the Himalaya of Nepal[END_REF] have been using a homogeneous quartzlike rheology for the crust. We thus suspect a trade off between the assumed denudation condition and the rheological properties of the crust. To test this hypothesis we calculate the deformation field resulting from various combinations of erosion law (diffusion or fluvial incision) and crustal properties (homogeneous or composite). As shown in figure 4, a complex interplay indeed exists between surface processes and crustal properties. A soft rheology for the lower crust (quartz) induces a gravitational collapse of the plateau. This localizes the maximum slope variations at the foothill of the HH, leading to a maximum of erosion by diffusion across the HH. In contrast, with the fluvial incision model, this collapse increases the river embanking, the denudation and uplift rates within the LH, which is found inconsistent with available data in LH [START_REF] Lavé | Fluvial incision and tectonic uplift across the Himalayas of central Nepal[END_REF]. The other end member diabase like crustal rheology implies a stronger coupling between mantle and crustal deformation. This model predicts a buckling of the lithosphere. Using diffusion law with this rheology gives sedimentation in the Lesser Himalaya, and a denudation peak far too the North. Fluvial incision, associated to this unrealistic diabase like upper crust, D R A F T February 15, 2013, 9:11am D R A F T displays also the same inconsistent peak of erosion in South Tibet and provides a poorer fit of the erosion profile than when using a composite rheology.

Conclusion

Among the reduced set of simulations presented in this paper, two combinations provide

  from the southern Tibet down to the Indo-Gangetic plain. In Central and East Nepal, across the HH, those trans-himalayan rivers are distant of 50km on average before joining, in the southern part of the LH, two major rivers systems, the Narayani and Sapt Kosi basins. Both rivers are tributary of the Ganga. Precipitations in Nepal are controlled by D R A F T February 15, 2013, 9:11am D R A F T the barrier of the Himalayas, with a brutal condensation against the HH of the wet air coming from the Indian ocean during monsoon. Whereas a marked rain shadow develops on the north flank of the HH, the prominent fluvial network of the south flank, fed by

(

  x<0), with active sedimentation and the mountain range north of the MFT, dominated by active erosional processes. The first domain, the Indo-gangetic basin is classically described as a low elevation over-filled basin[START_REF] Lyon-Caen | Gravity anomalies, flexure of the indian plate, and the structure support and evolution of the Himalaya and Ganga basin[END_REF]: we thus assume in the following a constant ∼ 0m elevation for it. In the range, two distinct erosion models are explored: a diffusion model and a detachment limited fluvial incision model including an implicit description of the tributaries and hillslope in the mean topography(Lavé, in prep.). The diffusion model is probably one of the most simple self-consistent model, in the sense that erosion, transport and deposition are described by a single relation. The D R A F T February 15, 2013, 9:11am D R A F T basic form of this relation is a 1-D linear diffusion law that simply writes :

  of a whole mountain range, because it can not account for the advective nature of the fluvial processes. Recent studies have, in particular, underlined the key role played by the fluvial incision in leading unglaciated landscape denudation [e.g. Whipple and Tucker, 1999]. Whereas different functional forms have been proposed to model fluvial incision [eg.

aFigure 1 .

 1 Figure 1. Topographic map of the studied area (GTOPO30 digital elevation model),

Acknowledgments. All the illustrations of this were prepared using GMT (Wessel and Smith, 1991). We thank J.-P. Avouac for helpful discussions. We are also grateful to J. Chéry and R. Hassani for providing the finite element code.

X -17

 [START_REF] Lavé | Fluvial incision and tectonic uplift across the Himalayas of central Nepal[END_REF], and to fission track datas from [START_REF] Burbank | Decoupling of erosion and precipitation in the Himalayas[END_REF] for the denudation profile (denudation rates obtained using the thermal model from [START_REF] Henry | Kinematic, thermal and petrological model of the Himalayas: Constraints related to metamorphism within the underthrust Indian crust and topographic elevation[END_REF]).

The control points on horizontal velocity are derived from folding of fluvial terraces in the Siwaliks (1) [START_REF] Lavé | Active folding of fluvial terraces across the Siwaliks Hills, Himalayas of central Nepal[END_REF], progradation of the sediments in the Gangetic Plain and flexure of the indian plate (2) [START_REF] Lyon-Caen | Gravity anomalies, flexure of the indian plate, and the structure support and evolution of the Himalaya and Ganga basin[END_REF] and quaternary grabbens extension in Southern Tibet (3) [START_REF] Armijo | Quaternary extension in southern Tibet[END_REF]. Topography is derived from GTOPO30 DEM, and river elevation from [START_REF] Lavé | Fluvial incision and tectonic uplift across the Himalayas of central Nepal[END_REF].


Jorge Guevara

Roberto Hirata
email: hirata@ime.usp.br

Stéphane Canu
email: scanu@insa-rouen.fr

Kernel Functions in Takagi-Sugeno-Kang Fuzzy System with Nonsingleton Fuzzy Input

Algorithms for supervised classification problems usually does not consider imprecise data, e.g., observed data whose samples can be represented by a collection of intervals, histograms, list of values, fuzzy sets among others. Fuzzy theory is a naturally choice for imprecise data. On the other hand, the state of the art techniques such a kernel methods are still a natural choice for supervised classification problems because of its robustness. Under some assumptions, Takagi-Sugeno-Kang (TSK) fuzzy system rules are equivalent to positive definite kernels (PDK) but such relationship was given considering only singleton fuzzy sets to model crisp input. Imprecise data are better modeled by membership functions of nonsingleton fuzzy sets but the relationship among nonsingleton fuzzy input, TSK fuzzy systems and PDK was up to our knowledge unknown. In this work, we study such relationship. We formulate an extension of TSK fuzzy systems to deal with nonsingleton fuzzy input and then we show that a new class of PDK are derived from it. We give three examples of nonsingleton TSK kernels which are close related to Vapnik's vicinal kernels. Also, based on those TSK induced kernels and on the concept of distance substitution kernels, we formulate two PDK for interval data. Potential applications for the proposed kernels are pattern recognition problems with imprecise data. Experiments conducted with interval datasets show better performance than state of the art approaches.

I. INTRODUCTION

Machine learning algorithms usually work only with crisp data or with crisp representation of imprecise data. Imprecise data is due to measurement errors, missing values, noise values, or it is usually added by instruments or by some preprocessing data transformations that involve sampling, randomness or numerical imprecision.

Some real examples of imprecise data are: [START_REF] Buckley | Fuzzy statistics[END_REF] geneexpression data -because of the complexity of the biological experiments, almost all causes written before are present, in special, the final quantification of the expression for each gene that is due to a statistical model of the many sources of variation; (2) some medical data -in the diagnosis of dyslexia, for instance, patients solve graphical tests which are subjectively evaluated by experts, in the form of linguistic terms (e.g., "near 10", "between 4 and 7"). Clinical data is prune to several problems, in special, missing values; [START_REF] Jeng | Support vector machines for the fuzzy neural networks[END_REF] weather databases -climate indicators along the day are better represented by intervals, for instance, temperature between 3 and 7 degrees, CO 2 between 394.1 and 394.2 parts per million, Wind between 7 and 10 miles per hour, Humidity between 10 and 16 percent and Sea level between 2.3 and 3.5 meters. [START_REF] Sohn | Advantages of using fuzzy class memberships in self-organizing map and support vector machines[END_REF] economics -some data is better represented by intervals or histograms among others, e.g., country's income distribution, histograms of people's systolic pressure, among others, are better described by density-like distributions or fuzzy sets. Some emergent areas such as fuzzy statistics [START_REF] Buckley | Fuzzy statistics[END_REF] and Symbolic Data Analysis (SDA) [START_REF] Billard | Symbolic data analysis: conceptual statistics and data mining[END_REF] search to deal with that kind of data.

From this perspective, machine learning algorithms able to perform some tasks using imprecise data would have a pragmatic importance. However, imprecise data is frequently preprocessed to obtain crisp values. The drawback is that we loose important information about the problem.

Fuzzy theory is the naturally choice to deal with imprecise data. On the other hand, kernel methods such as Support Vector Machines (SVM) are the state of the art methods in Machine Learning. It is worth noticing that fuzzy rule-base classifiers can be trained using SVMs, introducing some of the SVM's advantages in fuzzy classifier design [START_REF] Jeng | Support vector machines for the fuzzy neural networks[END_REF]- [START_REF] Cheng | An incremental support vector machinetrained ts-type fuzzy system for online classification problems[END_REF], such as the nonlinear mapping given by positive definite kernels.

A. kernel methods and fuzzy rules

The research done in machine learning using fuzzy set theory and positive definite kernels for data analysis can be grouped in two main groups. In the first one, fuzzy theories does not involve as a necessary consequence positive definite kernels (or vice versa), some specific fuzzy techniques and positive definite kernels are used to solve particular problems. Some works in clustering [START_REF] Kanzawa | On kernel fuzzy c-means for data with tolerance using explicit mapping for kernel data analysis[END_REF]- [START_REF] Raza | Interval type-2 approach to kernel possibilistic c-means clustering[END_REF], classification problems with outliers or noises [START_REF] Yang | A kernel fuzzy c-means clusteringbased fuzzy support vector machine algorithm for classification problems with outliers or noises[END_REF], feature extraction [START_REF] Wang | Fuzzy maximum scatter discriminant analysis with kernel methods[END_REF], discriminant analysis [START_REF] Heo | Robust kernel discriminant analysis using fuzzy memberships[END_REF] among others are in this group.

In the second group, fuzzy set concepts always imply positive definite kernels (or vice versa), we list some of them

(1) Kernels and fuzzy rules: under a general assumption of T-norm operators and membership functions, positive definite kernels are equivalent to interaction between the fuzzy system input and fuzzy rules [START_REF] Chen | Support vector learning for fuzzy rule-based classification systems[END_REF], [START_REF] Chen | Support vector machines and fuzzy systems[END_REF]. Papers [START_REF] Jeng | Support vector machines for the fuzzy neural networks[END_REF], [START_REF] Lin | Supportvector-based fuzzy neural network for pattern classification[END_REF]- [START_REF] Cheng | An incremental support vector machinetrained ts-type fuzzy system for online classification problems[END_REF] use this characteristic to train fuzzy classifiers with SVM.

(2) Kernels and fuzzy basis functions: using general fuzzy implication operators in fuzzy systems, and by constructing vectors with fuzzy basis function as components, a positive definite linear kernel is obtained by computing the inner product of those vectors [START_REF] Lski | On support vector regression machines with linguistic interpretation of the kernel matrix[END_REF].

(3) Kernels and fuzzy equivalence relations: positive definite kernels that maps to the unit interval with constant one in the diagonal can be seen as fuzzy equivalence relations, and vice versa under some assumptions [START_REF] Moser | On the t-transitivity of kernels[END_REF], [START_REF]On representing and generating kernels by fuzzy equivalence relations[END_REF]. A recent work [START_REF] Liu | Design of natural classification kernels using prior knowledge[END_REF] describe kernels that include some prior knowledge as fuzzy equivalence relations.

B. Contributions

All the precedent works relating positive definite kernels and fuzzy rules only consider crisp inputs in the form of singletons fuzzy sets in their analysis. Because imprecision in data is better modeled by some functions rather that crisp values, i.e., nonsingleton fuzzy sets, we will consider that situation in the context of positive definite kernel methods and Takagi-Sugeno-Kang fuzzy systems.

We claim the following contributions.

• Nonsingleton TSK fuzzy logic system. We present an extension of classical TSK fuzzy systems to deal with nonsingleton fuzzy input. • Nonsingleton TSK positive definite kernels. We show that interaction between nonsingleton TSK fuzzy rules and nonsingleton fuzzy inputs induce a class of positive definite kernel functions. • Kernels for interval data. Using nonsingleton TSK fuzzy kernels, and distance substitution kernels, we give a new general view to construct positive definite kernels for interval data. Experiments are conducted for all the proposed kernels using a Soft C-SVM [START_REF] Chang | LIBSVM: A library for support vector machines[END_REF] in datasets with samples containing missing and interval values.

II. FUZZY SYSTEMS

Fuzzy systems are universal approximators of functions [START_REF] Kosko | Fuzzy systems as universal approximators[END_REF]. They are composed of a set of if-then rules called the rule-base, a fuzzifier and, optionally, an inference algorithm and a defuzzifier.

A. Fuzzy set, Rule base and T-norm

Let U the universe of discourse. A fuzzy set F is the set defined on U and with membership function (MF)

µ F : U → [0, 1].
A rule-base is a set of If-Then rules indexed by l ∈ {1, 2, . . . , L}, of the form:

If u 1 is F l 1 and. . . and If u p is F l p Then v is G l , ( 1
)
where { u j } p j=1 ∪ { v} are called linguistic variables and { F l j } p j=1 ∪ G l are their respective linguistic values [START_REF] Zadeh | The concept of a linguistic variable and its applicationto approximate reasoning -ii[END_REF]. In fuzzy systems, linguistic values are represented by fuzzy sets. Rules can be completely described by fuzzy sets. For example, the rule l given by Equation (1) can be described by the fuzzy sets {F l j } p j=1 ∪ G l defined in their respective universes of discourse {U j } p j=1 ∪ V .

Fuzzy relations, i.e., fuzzy sets defined in the Cartesian product of several universes of discourses, are used to represent the if part of rules or antecedent part. The MF's of such fuzzy relations are obtained by T-norms operators [START_REF] Klement | Triangular norms. position paper i: basic analytical and algebraic properties[END_REF] which are used to implement conjunction in fuzzy logic and consequently fuzzy set intersection.

A triangular norm or T-norm is the function T : [0, 1] 2 → [0, 1], that for all x, y, z ∈ [0, 1] satisfy: T1 commutativity: T (x, y) = T (y, x); T2 associativity: T (x, T (y, z)) = T (T (x, y), z); T3 monotonicity: y ≤ z ⇒ T (x, y) ≤ T (x, z); T4 boundary condition T (x, 1) = x. Using n ∈ N and associativity, a multiple-valued extension

T n : [0, 1] n → [0, 1] of a T-norm T is given by T n (x 1 , x 2 , . . . , x n ) = T (x 1 , T n-1 (x 1 , x 2 , . . . , x n-1 )). (2)
We will use T to denote T or T n .

Definition

is represented by the fuzzy relation

A l ⊂ p j=1 U j with MF µ A l : p j=1 U j → [0, 1] given by µ A l (x) = T (µ F l 1 (x 1 ), µ F l 2 (x 2 ), . . . , µ F l p (x p )), (4)
Following Definition (II.1), fuzzy rules from Equation (1) can be characterized by fuzzy relations

R l ⊂ p j=1 U j × V with MF µ R l : p j=1 U j × V → [0, 1] given by µ R l (x, y) = I(µ A l (x), µ G l (y)), (5)
where I is some fuzzy implication operator (see [START_REF] Shi | On dependencies and independencies of fuzzy implication axioms[END_REF]).

B. Fuzzifier and fuzzy input

Because all the operations performed in fuzzy systems are done using fuzzy sets, it is necessary to turn crisp input into fuzzy sets, this is done by the fuzzifier.

The fuzzifier transform the crisp input

(x 1 , . . . , x p ) ∈ p j=1 U j into fuzzy sets X 1 ⊂ U 1 , . . . , X p ⊂ U p with MF's {µ Xj : U j → [0, 1]} p j=1 1 .
The singleton fuzzifier converts a value x j ∈ R in a singleton fuzzy set whose support is a single value. i.e., the membership function satisfy µ Xj (x j ) = 1 and µ Xj (x) = 0 for all x ∈ R with x = x j .

On the other hand, the nonsingleton fuzzifier converts the value x j ∈ R in a fuzzy set whose support are several values, for example, the Gaussian fuzzifier converts the point x j in a fuzzy set with membership function µ Xj (x) = exp(-γ(x -x j ) 2 ), γ ∈ R + and the triangular fuzzifier in the fuzzy set with membership function µ

Xj (x) = (1 -|x -x j |/b j ), b j ∈ R + . Definition II.2 (Fuzzy input). Let X 1 ⊂ U 1 , . . . , X p ⊂ U p be fuzzy sets with MF's {µ Xj : U j → [0, 1]} p j=1 obtained after the fuzzification of x = (x 1 , x 2 , . . . , x p ) ∈ p j=1 U j . The fuzzy input is the fuzzy relation I ⊂ p j=1 U j with MF µ I : p j=1 U j → [0, 1] given by µ I (x) = T (µ X1 (x 1 ), µ X2 (x 2 ), . . . , µ Xp (x p )). (6)
We will call singleton fuzzy input if it was used the singleton fuzzifier or nonsingleton fuzzy input if it was used the nonsingleton fuzzifier.

C. Inference Algorithm

The third part of a fuzzy system is an inference algorithm whose arguments are the fuzzy input and the rule-base. In Mamdani fuzzy systems, the inference algorithm uses fuzzy composition between the fuzzy input I (Definition II.2) and one element of the rule base given by R l (Equation 5) to obtain the fuzzy set Y l with membership function

µ Y l (y) = sup x∈R p T (µ I (x), µ R l (x, y)).
The system output is a fuzzy set given by computing

⊕ L l=1 Y l .
where ⊕ is a T-conorm operator [START_REF] Klement | Triangular norms. position paper i: basic analytical and algebraic properties[END_REF]. Optionally, it can be used a defuzzifier to convert the output fuzzy set to a crisp value.

D. Takagi-Sugeno-Kang fuzzy system

Takagi-Sugeno-Kang (TSK) fuzzy system does not have a fuzzy set in the consequent part of their rules, it has a function of the p antecedents values instead, then, it is not necessary either fuzzy composition or defuzzification step. In first order TSK fuzzy systems each rule is written as

If u 1 is F l 1 . . . If u p is F l p Then g l (x) = c l 0 + p j=1 c l j x j . (7)
The input space is divided into fuzzy regions by antecedents (if -part) and the system behavior in those regions is described by consequents in the form of functions g l (x) : R p → R given by g l (x) = c l 0 + x, c , x, c ∈ R p , c l 0 ∈ R. Following Definition (II.1) and setting all the universe of discourses equals to R, each rule of a TSK fuzzy system define a function f l : R p → R given by

f l (x) = g l (x)µ A l (x), (8)
The output of the first order TSK fuzzy system is given by the combination of the their L rules in the following manner

h(x) = L l=1 f l (x) L l=1 µ A l (x) . (9)

III. KERNELS AND FUZZY RULES

A Kernel is a function k from X × X to R, where X is an arbitrary finite set. The function k is called positive definite if the n × n matrix Q ij = k(x i , x j ) is positive semidefinite, that is, c ⊤ Qc ≥ 0 for any choice of c ∈ R n and any choice of x i ∈ X (see for instance ref. [START_REF] Schölkopf | Learning with kernels : support vector machines, regularization, optimization, and beyond, ser. Adaptive computation and machine learning[END_REF] or [START_REF] Vapnik | The nature of statistical learning theory[END_REF] for details). In machine learning this property is important because k(x i , x j ) defines implicitly an inner product Φ(x i ), Φ(x j ) in a Reproducing Kernel Hilbert Space H k , using the implicit map Φ : X → H k . Also, positive definite kernels are useful to construct optimization algorithms for machine learning problems, ensuring that kernel algorithms converge to a relevant solution.

A. Takagi-Sugeno kernel

By the form of consequents, there are two positive definite kernels induced from interaction between TSK fuzzy rules and fuzzy inputs.

Lemma 1 (PDFC kernel [START_REF] Chen | Support vector learning for fuzzy rule-based classification systems[END_REF]). If all the MF's of the fuzzy sets of the antecedent part of the l rule (Definition (II.1)) are generated by translating the positive definite functions u : R → [0, 1] defined as:

u(x) = u(-x) if x = 0 1 if x = 0, , x ∈ R, (10)
such that µ F l j (x j ) = u j (x j -z l j ) for some parameter z l j ∈ R, then the kernel k : R p × R p → R given by k(x, z l ) = p j=1 u j (x j -z l j ) = p j=1 µ F l j (x j ), x, z ∈ R p , ( 11
)
is a positive definite invariant translation kernel.

The proof can be found in reference [START_REF] Chen | Support vector learning for fuzzy rule-based classification systems[END_REF] (Theorem 3.11). Lemma (1) assume algebraic product as T-norm operator. By setting g(x) l = c 0 , where c 0 is some constant value. Equation ( 8) can be rewritten as

f l (x) = c 0 k(x, z l ). (12)
It is worth to noting that by using Gaussian functions in Equation ( 10), the gaussian RBF kernel is obtained, but in general PDFC kernels are not RBF kernels [START_REF] Chen | Support vector machines and fuzzy systems[END_REF].

Lemma 2 (First order TSK-kernel [START_REF] Juang | A self-organizing ts-type fuzzy network with support vector learning and its application to classification problems[END_REF], [START_REF] Cheng | An incremental support vector machinetrained ts-type fuzzy system for online classification problems[END_REF]). If all the MF's of the fuzzy sets of the antecedent part are gaussian functions, and the consequent part is the function g l (x) = x, z , then the kernel k : R p × R p → R given by

k(x, z l ) = x, z l exp(-γ||x -z l || 2 ), (13)
is positive definite.

The proof is given in [START_REF] Juang | A self-organizing ts-type fuzzy network with support vector learning and its application to classification problems[END_REF], [START_REF] Cheng | An incremental support vector machinetrained ts-type fuzzy system for online classification problems[END_REF]. Using Lemma (2), Equation (8) can be rewritten as

f l (x) = k(x, z l ). (14)
It is worth to noting that a more general case can be derived using PDFC kernels for the antecedent part and the inhomogeneous polynomial kernel in the consequent part to represent TSK fuzzy rules, supported by the following fact. Lemma 3. Let µ A l (x) in Equation (8) equals to the PDFC kernel k(x, z l ) given in Lema [START_REF] Buckley | Fuzzy statistics[END_REF]. Let the consequent function g l (x) = ( x, z l ) + z 0 ) d for d = 0, 1, . . . , n, then the TSK fuzzy rule (Equation (8)) defines the positive definite kernel

f l (x) = k(x, z l ) = ( x, z l ) + z 0 ) d k(x, z l ). ( 15
)
Proof: The inhomogeneous polynomial kernel is a positive definite kernel, the PDFC is positive definite. The product of two positive definite kernels is also positive definite.

IV. KERNELS AND NONSINGLETON FUZZY INPUT

Positive definite kernels related to TSK rules described in the previous section only consider singleton fuzzy inputs, In this section, we define an extension of TSK fuzzy systems for nonsingleton fuzzy inputs and then we relate the resulting fuzzy rules to positive definite kernel functions.

A. TSK fuzzy system with nonsingleton input

Nonsingleton fuzzy logic system (NFLS) [START_REF] Mouzouris | Nonsingleton fuzzy logic systems: theory and application[END_REF] is a Mamdani fuzzy system that accounts for the uncertainty in the input. NFLS has been applied successfully in engineering applications [START_REF] Mouzouris | Nonsingleton fuzzy logic systems: theory and application[END_REF]- [START_REF] Kim | Modeling corrupted time series data via nonsingleton fuzzy logic system 2 nonsingleton fuzzy logic system[END_REF], with better results that singleton Mamdani fuzzy systems [START_REF] Cara | On comparing non-singleton type-1 and singleton type-2 fuzzy controllers for a nonlinear servo system[END_REF]. Its principal characteristic is the nonsingleton fuzzification step with the effect of noise suppression [START_REF] Mouzouris | Nonsingleton fuzzy logic systems: theory and application[END_REF].

In TSK fuzzy systems only the antecedent part of the rules is represented by fuzzy relations (Definition (II.1)), because the consequent is a function. For TSK fuzzy system can to deal with nonsingleton fuzzy inputs we compute a representative value of the intersection between the fuzzy relations that represent the nonsingleton fuzzy input and rules.

Definition IV.1. Let A l , I ⊂ p j=1 U j two fuzzy relations that represent the antecedent part of the l TSK rule and the nonsingleton fuzzy input with MF's given by Definition (II.1) and Definition (II.2) respectively. Let g l : p j=1 U j → R the function that represent the consequent part of the TSK rule. The function f l : p j=1 U j → R induced by A l , I and g l for the rule l is given by

f l (x) = g l (x)r(µ A∩B (x)), (16)
where µ A∩B :

p j=1 U j → [0, 1]
is the MF of the fuzzy set A ∩ I given by

µ A∩B = T (µ I (x), µ A l (x)), (17)
and r :

[0, 1] → [0, 1].
Using this definition, TSK fuzzy system can naturally extended to deal with nonsingleton fuzzy input.

Definition IV.2 (Output of nonsingleton TSK fuzzy system).

The output of TSK fuzzy system with nonsingleton fuzzy input is given by the function ĥ :

p j=1 U j → R given by ĥ(x) = L l=1 f l (x) L l=1 r(µ A∩B (x)) = L l=1 g l (x)r(µ A∩B (x)) L l=1 r(µ A∩B (x)) . (18)
In the same manner that classical TSK fuzzy systems, the output of the above formulation is a weighted combination of functions {g l } M l=1 by values between [0, 1]. Lemma 4. Let function r (Definition (IV.1)) be the sup function. If the fuzzifier is a singleton fuzzifier, then the output of nonsingleton TSK fuzzy system (Definition (IV.2)) reduces to output of TSK fuzzy system (Equation (9)) Proof:

ĥ(x) = L l=1 g l (x) sup x∈U1×...Up {T (µ I (x), µ A l (x))} L l=1 sup x∈U1×...Up {T (µ I (x), µ A l (x))} = L l=1 g l (x) sup x∈U1×...Up {T (1, µ A l (x))} L l=1 sup x∈U1×...Up {T (1, µ A l (x))} = L l=1 g l (x)µ A l (x) L l=1 µ A l (x)
Then ĥ = h. Let I a nonempty and finite set of indices. Let µ i , µ j ∈ X for i, j ∈ I. Then, the kernel k : X × X → [0, 1] given by

B. Relation between TSK rules for nonsingleton fuzzy inputs and positive definite kernels

k(µ i , µ j ) = sup x∈U T (µ i (x), µ j (x)), (19)
is positive definite.

Proof: By commutativity property of T-norms k is symmetric. Noting that i,j∈I

c i c j k(µ i , µ j ) = i∈I c 2 i k(µ i , µ i )+2 i>j,i,j∈I c i c j k(µ i , µ j )
and

sup x∈U T (µ i (x), µ i (x)) = 1, ∀i ∈ I then i,j∈I c i c j k(µ i , µ j ) = i∈I c 2 i + 2 i>j,i,j∈I c i c j k(µ i , µ j )
Using ( i∈I c i ) 2 = i∈I c 2 i + 2 i>j,i,j∈I c i c j ≥ 0 and by the fact that k(µ i , µ j ) ∈ [0, 1], we have a)

If k(µ i , µ j ) = 0, ∀i, j ∈ I : i > j, then i,j∈I

c i c j k(µ i , µ j ) = i∈I c 2 i ≥ 0 b) If k(µ i , µ j ) = 1, ∀i, j ∈ I : i > j, then i,j∈I c i c j k(µ i , µ j ) = i∈I c 2 i + 2 i,j∈I,i =j.i>j c i c j = ( i∈I c i ) 2 ≥ 0
Thus k is positive definite.

Lemma 6. Function induced by nonsingleton TSK fuzzy rule (Equation ( 16)) can be written as

f l (x) = g l (x)k(µ I , µ l A ) (20)
Proof: Let function r (Definition (IV.1)) be the sup function, substituting µ I = µ i , µ l A = µ j in Equation ( 17) and by Lemma (5) completes the proof By basic properties of kernels, if g l (x) = c 0 where c 0 is some positive constant then

f l (x) = c 0 k(µ I , µ l A ) is positive definite, also if If g l (x) is a positive definite function then f l (x) = g l (x)k(µ I , µ l
A ) is positive definite The above results mean that we could use whatever T-norm operator and MF's of normal fuzzy sets, we could obtain a positive definite kernel k.

Also, kernel given by Lemma ( 5) is a fuzzy equivalence relation [START_REF] Moser | On the t-transitivity of kernels[END_REF] with respect to a given T-norm and can be represented by a fuzzy bi-implication formula [START_REF]On representing and generating kernels by fuzzy equivalence relations[END_REF], because every positive definite kernel that maps to the unit interval with constant one in the diagonal fullfill these requeriments. Details are omitted due to space constraints.

C. Positive definite TSK gaussian kernels

We follows with a classical result of NFLS's [START_REF] Mouzouris | Nonsingleton fuzzy logic systems: theory and application[END_REF].

Lemma 7. Let µ X1 , µ X2 , . . . , µ Xp be the MF's of the fuzzy sets of the fuzzy input with parameters given by (m 1 , . . . , m p ) ⊤ ∈ R p , (σ 1 , . . . , σ p ) ⊤ ∈ R p and let µ l F1 , µ l F2 , . . . , µ l Fp be the MF's of the fuzzy sets of the antecedent part of the rule l with parameters given by

(m l 1 , . . . , m l p ) ⊤ ∈ R p , (σ l 1 , . . . , σ l p ) ⊤ ∈ R p such that for j = 1, 2 . . . p µ Xj (x j ) = exp -0.5(x j -m j ) 2 /σ 2 j , µ F l j (x j ) = exp -0.5(x j -m l j ) 2 /(σ l j ) 2 ,
Let µ l A and µ I given by Definition (II.1) and (II.2) respectively.

If is used algebraic product as T-norm operator, then

sup x∈R p T (µ I (x), µ A l (x)) = p j=1 exp - 1 2 (m j -m l j ) 2 σ 2 j + (σ l j ) 2 , ( 21
)
The proof is in [START_REF] Mendel | Uncertain rule-based fuzzy logic system: introduction and new directions[END_REF].

Lemma 8. Let the set of MF's X = {µ|µ : R p → [0, 1] and ∃x ∈ R p : sup

x∈R p (µ(x)) = 1}. Let µ i , µ j ∈ X for i, j ∈ I = {1, 2, . . . , N }. The kernel k : X × X → [0, 1] given by k(µ i , µ j ) =sup x∈U T (µ i (x), µ j (x)) = p j=1 exp - 1 2 (m j -m l j ) 2 σ 2 j + (σ l j ) 2 , ( 22
)
is positive definite.

Proof: By Lemma (5), k is a positive definite kernel.

Lemma 9. Let k given in Lemma [START_REF] Wee | Efsvmrobus-fcm: Evolutionary fuzzy rule-based support vector machines classifier with fcm clustering[END_REF]. If g l given in Definition (IV.1) is a positive definite function or g l is an affine function, then the function f l given in Definition (IV.1) can be written as the positive definite kernel k :

R p × R p → [0, 1] given by k(x, x l ) = g(x -x l )k(µ I , µ A l ). ( 23
)
Proof: Kernel g is positive by hypothesis and k through the previous lemma. The results product kernel is then positive since the product of two positive kernels is positive, then following Lemma (6) completes the proof.

Note that the kernel k γ : X × X → [0, 1] given by

k γ (µ I , µ A l ) = p j=1 exp - 1 2 (m j -m l j ) 2 σ 2 j + (σ l j ) 2 + γ , (24)
is also positive definite by the change of variables σ ′2 = σ 2 + γ.

Parameter γ plays the same regularization role as it does in the RBF kernel. Kernels given by Equations ( 22) and ( 24) are closed related to Gaussian vicinal kernels defined by Vapnik in [START_REF] Vapnik | The nature of statistical learning theory[END_REF]. These kernels can be viewed as vicinal kernels where the vicinities are defined by the volume given by the spread of each feature.

D. kernels for interval data

Let's assume that imprecise data is given as the interval

[x, x] = [(x 1 , x 2 , . . . , x p ) ⊤ , (x 1 , x 2 , . . . , x p ) ⊤ ] ⊂ R p , where each feature {[x j , x j ]} p j=1 is a closed interval in R with {x j ≤ x j } p j=1 .
In this context, all the induced TSK kernels with nonsingleton fuzzy input can be used, e.g., it is possible to construct a fuzzy data set by setting one fuzzy set X j with MF µ Xj for each interval [x j , x j ].

It is possible to derive new positive definite kernels for fuzzy data using the concept of distance substitution kernels [START_REF] Haasdonk | Learning with distance substitution kernels[END_REF], for example the distance substitution Gaussian kernel:

k rbf d (x, x ′ ) = exp(-γd(x, x ′ ) 2 ), γ ∈ R + , x, x ′ ∈ X (25)
only needs to construct a specific distance measure.

Lemma 10. Let µ = (µ X1 , µ X2 , . . . , µ Xp ) ⊤ and µ ′ = (µ ′ X1 , µ ′ X2 , . . . , µ ′ Xp ) ⊤ two vectors of MF's with param- eters: (m 1 , . . . , m p ) ⊤ ∈ R p , (σ 1 , . . . , σ p ) ⊤ ∈ R p and (m ′ 1 , . . . , m ′ p ) ⊤ ∈ R p , (σ ′ 1 , . . . , σ ′ p ) ⊤ ∈ R p respectively, constructed for the intervals [x, x] ⊂ R p and [x ′ , x ′ ] ⊂ R p respectively. Functions d 1 (µ, µ ′ ) = p j=1 (m j -m ′ j ) 2 p j σ 2 j + p j (σ ′ j ) 2 d 2 (µ, µ ′ ) = p j=1 (m j -m ′ j ) 2 (σ 2 j + (σ ′ j ) 2 )
have zero diagonal and are pseudometrics on X .

Proof: d 1 and d 2 are symmetric and have zero diagonal by construction. They are positive as the product of positive function and subadditive. Thus both functions are pseudometrics distances since they can be set to zero for some distinct variances values.

Lemma 11. The kernels

k γ (µ, µ ′ ) = exp(-γd 1 (µ, µ ′ ) 2 ), γ ∈ R + (26)
k γ (µ, µ ′ ) = exp(-γd 2 (µ, µ ′ ) 2 ), γ ∈ R + (27)
are positive definite Proof: By definition of substitution kernel both kernels are positive definite Of course, many other kernels can be derived using the this concept and this is question of further research.

V. EXPERIMENTS

A. Data and Implementation

Four low quality datasets from the KEEL-dataset repository [START_REF] Alcal-Fdez | Keel data-mining software tool: Data set repository, integration of algorithms and experimental analysis framework[END_REF] were used. Low quality data [START_REF] Palacios | Extending a simple genetic cooperative-competitive learning fuzzy classifier to low quality datasets[END_REF] refers to data with uncertainty about the actual value of a feature. These KEEL datasets contain samples with missing values and interval features. Table I contains the summary of these datasets.

1) Description of the data:

• Dataset Long-4 is an athletic performance dataset, it contains 4 features. The two first features are determined by a coach in form of linguistic values, intervals or numbers; the other two indicators are measured three times producing uncertainty information represent as intervals. • Dataset 100mlI-4 is an athletic performance dataset used to classify whether or not a mark in a 100 meters race is being achieved. The features are given by the weight to height ratio, the reaction time, the starting speed, and 40m. speed. Measurements were obtained by three different observers.

• Dataset 100mlP-4 is the same as Dataset 100mlI-4 but measurements were obtained by a subjective judgment of the coach, in linguistic terms as " reaction time is low". • Dataset Dyslexic-12-4 is a dataset with twelve features and four classes { dyslexia, no dyslexia, control, other} that contain missing values. Table II shows two examples of Long-4 dataset where each feature is an interval. A detailed description of these datasets can be found in [START_REF] Palacios | Extending a simple genetic cooperative-competitive learning fuzzy classifier to low quality datasets[END_REF], [START_REF] Palacios | Future performance modeling in athletism with low quality data-based genetic fuzzy systems[END_REF]. 2) Scaling and fuzzification: Some samples of these datasets belong to two different classes at the same time. Because it is not the purpose of this work, they were removed. Next, intervals were scaled (by a linear transformation) to be in the unit square.

For each (feature) interval [x j , x j ] of the four datasets, a fuzzy set with Gaussian MF µ Xj with m j = (x j + x j )/2 and σ i j = (x jx j )/(2 √ 2 ln(2)) has been built2 . That is, for each sample [x, x] = [(x 1 , x 2 , . . . , x p ) ⊤ , (x 1 , x 2 , . . . , x p ) ⊤ ] ⊂ R p , we have the fuzzy sets X 1 , . . . , X p with MF's µ X1 , . . . , µ Xp with the following parameter vectors (m 1 , m 2 , . . . , m p ) ⊤ and (σ 1 , σ 2 , . . . , σ p ) ⊤ obtained directly from the data.

In the case of missing values, we constructed an interval that spans the whole range of the variable.

B. Kernels Setting

Were used the kernels given by Equations ( 22), ( 23), ( 24), ( 26) and [START_REF] Schölkopf | Learning with kernels : support vector machines, regularization, optimization, and beyond, ser. Adaptive computation and machine learning[END_REF]. For comparison purposes, it was set the RBF kernel as baseline.

Table III show three different setting for experiments. The first one consider model selection (see Section V-B2 for details) over the γ parameter. Kernels with id ∈ {0, 1, 2, 3, 4} are in this group.

The second group is composed by one kernel id = 5 and it does not have model selection.

The last one group is composed of kernels id ∈ {6, 7, 8, 9} and experiments in these kernels consider model selection over γ and multiple kernel setting for fuzzy rules.

1) Multiple kernel setting for fuzzy rules: By considering each fuzzy rule as combination of two (or more) rules and, using positive definite kernels to represent it. The resulting kernel is given by k(x, x ′ ) = r β r k s (x, x ′ ), for β r > 0, and r β r = 1. That permits to use more general approach as the multiple kernel learning [START_REF] Lanckriet | Learning the kernel matrix with semidefinite programming[END_REF] to optimize the β r parameters. For the third group of kernels id = {6, 7, 8, 9} from Table III it was set arbitrarily β = 0.25.

Kernel arguments in

k rbf γ (x, x ′ ) 1 k eq(24) γ (µ X , µ X ′ ) 2 a k eq(23) (x, x ′ ) 3 b βk eq(26) γ (µ X , µ X ′ ) 4 b βk eq(27) γ (µ X , µ X ′ ) No model selection 5 c k eq(23) (x, x ′ )
Multiple kernel setting and model selection over γ 6 βk eq (22) (µ X , µ

X ′ ) + (1 -β)k rbf γ (x, x ′ ) 7 a βk eq(23) (x, x ′ ) + (1 -β)k rbf γ (x, x ′ ) 8 b βk eq(26) γ ′ (µ X , µ X ′ ) + (1 -β)k rbf γ (x, x ′ ) 9 b βk eq(27) γ ′ (µ X , µ X ′ ) + (1 -β)k rbf γ (x, x ′ )
a Function g was setting to the rbf kernel. b The γ ′ parameter in k eq (26) was 1.

c Function g was setting to the one-degree polinomial kernel.

2) Model Selection: It was performed model selection in the kernel parameter γ with values {2 4 , . . . , 2 -10 } and in the parameter C of soft C-SVM with values {2 -1 , . . . , 2 14 } using grid search. For each point in the grid, it was performed 10-fold cross-validation using partitions from KEEL-dataset repository [START_REF] Alcal-Fdez | Keel data-mining software tool: Data set repository, integration of algorithms and experimental analysis framework[END_REF]. Finally, we reported the pair (C, γ) with best cross validation accuracy.

C. Results

Table IV shows the cross validation accuracy, number of support vectors and the best parameters C and γ from the model selection step. A soft C-SVM [START_REF] Chang | LIBSVM: A library for support vector machines[END_REF] were used to test all the proposed kernels. The Results show that using the proposed kernels,in most of the cases, outperform the RBF kernel for those datasets. Results in the datasets differ, showing that the kernel choice is a crucial part of the supervised classification problem.

VI. CONCLUSION

Real world applications of supervised classifiers involve imprecise data. This is due to the fact that no measurement instrument can produce an exact result, usually data is gathered by subjective opinions, and is better expressed by linguistic, interval or fuzzy values. In kernels methods, the choice of a kernel function is crucial because it constitutes prior knowledge about a task. No free lunch theorem states that learning algorithms with better generalisation capabilities are obtained by using some prior information of the domain.

In this work, we show that the relationship among TSK fuzzy systems, nonsingleton fuzzy input and positive definite kernels a new class of positive definite kernels can be derived. Also we give a formulation to construct kernels for interval data using the concept of distance substitution kernels. Experiments performed in interval datasets show promising results.

Lemma 5 (

 5 Positivity of the nonsingleton TSK fuzzy kernel). Let U = p j=1 U j the Cartesian product of universes of discourse. Let the set of MF's of normal fuzzy sets. X = {µ|µ : U → [0, 1] and ∃x ∈ U : sup x∈U (µ(x)) = 1}.

  If u 1 is F l 1 and. . . and If u p is F l p ,

	II.1 (Rule antecedent part representation by fuzzy relation and T-norm). Let p j=1 U j the Cartesian product of universes of discourse. Let x = (x 1 , x 2 , . . . , x p ) ∈ p j=1 U j . Let {F l j ⊂ U j } p j=1 fuzzy sets with their respective MF's {µ F l j } p j=1 . The antecedent part of the l rule given by

  Table III denoted by x and x ′ correspond to crisp values, we set those values to be the means of the intervals, i.e., x = (m 1 , m 2 , . . . , m p ) ⊤ .

The fuzzifier can also be applied to transform sets into fuzzy sets or to increase the fuzziness of fuzzy sets, but we are only interested in the above definition.

This measurement for σ parameter is called the full width at half maximum.

ACKNOWLEDGMENT

The authors would like to thank CNPq and CAPES (Brazilian Agencies) for their financial support.