Marc Louis

Maurice François

M Arc

Fran Coi

-convex f ree energy f or quasi -bri ttl em ateri al s

teaching and research institutions in France or abroad, or from public or private research centers.

A damage model with non-convex free energy for quasi-brittle materials

Introduction

A ttheul ti m atestateofdam age,quasi -bri ttl em ateri al scan beregarded asgranul ar m ateri al s w hose behavi or,due to grai n i nterl ock,i s strongl y dependent on the con nem ent. C l assi cal yi el d cri teri a (M ohr-C oul om b, D rucker-Prager, . . . ) use the con nem ent (the hydrostati c strai n) as a rei nforcem ent factor. B ut there i s experi m entalevi dence that el asti ci ty of pow der m ateri al s, thus thei r free energy,al so depend on i t [ 1] . T hatpoi nti sgeneral l y nottaken i nto account i n dam age m odel s. A s a consequence,the l arge di l antancy observed at these states,for exam pl e on concretes [ 2,3] ,i s poorl y descri bed by cl assi caldam age m odel s as apparent Poi sson' s rati o cannot exceed 0. 5,the l i m i t for l i near el asti ci ty,w hi l e experi m ents exhi bi t greater val ues. A nother speci ci ty ofconcrete behavi ori s the very di erentcrack angl esobserved i n tensi on and com pressi on. In m ostdam agem odel s,theH adam artand R i cecri teri on [START_REF] W R Udni Ckiand J R I Ce | C ondi ti ons for the l ocal i zati on ofdeform ati on i n pressure-sensi ti ve di l atantm ateri al s[END_REF][START_REF] En E B I L L Adon | O n l ocal i zati on i n ducti l e-bri ttl e m ateri al s under com pressi ve l oadi ngs[END_REF]l eadsto thesam e l ocal i zati on angl e i n both cases.Furtherm ore,fora concrete speci m en,even after a rupture i n tensi on,som e carryi ng capaci ty i n com pressi on rem ai ns: thi s i m pl i es a dam age l evelat the onset ofl ocal i zati on m uch l ower i n tensi on than i n com pressi on. T he present m odel descri bes these e ects that are general l y m i ssed by m ost ofdam age m odel s.

T hi swork,i n conti nui ty w i th [START_REF] Le | C ontribution a l a m od el isation du com portem entdesb etons[END_REF] ,consti tutesan attem ptto use a non convex potenti ali n the el d ofdam age m echani cs. For thi s reason,a si m pl e i sotropi c dam agel aw i sconsi dered.A si n the K el vi n' sapproach ofel asti ci ty [START_REF] Ski | O n hooke' s l aw[END_REF] ,pl asti ci ty theory,soi lm echani csand som e dam age m odel s [START_REF] Pi | O n ani sotropi c dam age theory. Fail ure criteria ofStructured m edia[END_REF] ,the i sotropi c and devi atori c decom posi ti on i s used. T he retai ned yi el d cri teri on [START_REF] Coi | A new yi el d cri teri on for the concrete m ateri al s[END_REF]i s sm ooth and convex. T he m odeli s i denti ed w i th respect to the wel lknow n uni axi aland m ul ti axi al testi ngs of [ 3] .

C onstitutive law

H elm holtz free energy. T he dam age l eveli s descri bed by the scal arvari abl e d that ranges from 0 for sound m ateri alto 1 for ful l y dam aged m ateri al [START_REF] Lem Ai Tre | A C ourse on D am age M echanics[END_REF] . T he presentm odeli san i sotropi c one:the hydrostati c and devi atori c parti ti ons ofthe stress = h + d and the strai n " = " h + " d are used (detai l s are gi ven i n A ppendi x A ).T he state vari abl es(" d ;" h ;d)descri be the m ateri al ' sstate and, w i th respectto the general i zed standard m ateri alfram ework [START_REF] Uoc Son N Guyen | Sur l es m at eri aux standards g en eral i s es[END_REF] ,the associ ated therm odynam i c forcesare respecti vel y ( d ; h ;Y )w here Y i sthe energy rel ease rate densi ty. T he proposed free energy i s:

2 (" d ;" h ;d)= 3K " h :" h + 2 1 d(1 + 2'" h ) " d :" d ;
(1) w here " h = tr(")= p 3 (A ppendi x A ). T he new constant i ntroduced i s ', the other ones are the m ass densi ty and the bul k and shear m odul i ,respecti vel y K and . T he rol e of' w i l lbe detai l ed further but i t can be al ready seen that i t i ntroduces a cubi c term i n the free energy and that setti ng ' = 0 l eads to a si m pl e dam age m odeli n w hi ch onl y the devi atori c part i s a ected by dam age. T hi s choi ce has been m ade i n orderto getri d ofthe uni l aterale ectofdam age on bul k m odul us (that com es from crack openi ng and cl osure) because the use of posi ti ve parts of the strai n tensor i nduces di cul ti es associ ated w i th non regul ar free energy [START_REF] Pi Erre B Ernard B Adel ; B Apti Ste Lebl Ond | A ppl icati on ofsom e ani sotropi c dam age m odelto the predi cti on ofthe fai l ure of som e com pl ex i ndustri alconcretestructure[END_REF] . In Eq.1 tensors " h and " d can be repl aced by thei r al gebrai c val ues " h and " d (A ppendi x A ).

Stress to strain relationship. T he hydrostati c and devi atori c stresses,as therm odynam i c forces,are obtai ned by di erenti ati on ofthe free energy w i th respect to " h and " d :

h = 3K " h 2 p 3 'd(" d ) 2 I; (2) d = 2 1 d(1 + 2'" h ) " d : (3)
Setti ng d = 0 l eads to recover the l i near i sotropi c el asti ci ty l aw i n the K el vi n' s decom posi ti on form [START_REF] Ski | O n hooke' s l aw[END_REF] . T he devi atori c stress and strai n rem ai n col l i near together(and col l i nearto the uni tary tensorD ,see A ppendi x A ).O nce projected on the orthogonaltensor base (I;D ),the previ ous expressi on becom es:

h = 3K " h 2 'd(" d ) 2 ; ( 4
) d = 2 1 d 2'd" h " d : (5)
Letussupposean i m posed devi atori cstrai n w hi l ethehydrostati cstrai n rem ai ns equalto zero (i sochori c transform ati on) then a con ni ng pressure ( h < 0) i s necessary to keep the vol um e unchanged:

h (" h = 0;" d )= 2 'd(" d ) 2 : (6)
Letussupposenow thatthedevi atori cstrai n i sjoi ned w i th an hydrostati cstress i m posed equalto zero then an i nduced di l atati on " h > 0 ari ses:

" h ( h = 0;" d )= 2 'd(" d ) 2 1 d : (7)
T hese two aspectsare rel evantto the di l atancy e ecti n concrete-l i ke m ateri al s.

T he m i cro-m echani calpoi ntofvi ew associ ated to these e ectsi sthe surm ounti ng ofconcrete parti cl es,i nduci ng voi ds creati on,that ari ses under i rreversi bl e sheari ng [START_REF] Yazdaniand | A n ani sotropi c dam age m odelw i th di l atati on for concrete[END_REF][START_REF]coi sand G i anniR oyerC arfagni .Structured deform ati on ofdamaged conti nua w i th cohesi ve-fri cti onalsl i di ng rough fractures[END_REF] . E lastic tangent m odulus. T he tangentm odul usi sde ned asH = d =d". In caseofel asti ctransform ati on,d rem ai ns xed,then H reducesto H 0 = @ =@". From the stress to strai n rel ati ons (2,3),we have:

H 0 = 3K 2~ 3 I I + 2~ I I " d + " d I ; (8)
2~ = @ d @" d = 2 (1 d) 4 'd" h ; (9) = 4 'd p 3 : (10)
In thi s expressi on I i s the fourth order i denti ty w hose expressi on i n i ndex form i s [ I] ijkl = ( ik jl + il jk )=2, the sym bol refers to the tensor product and (from Eq.5) ~ represents the apparentshear m odul us. T hen H 0 has the i ndex sym m etri es ofan el asti ci ty tensor: [

H 0 ] ijkl = [ H 0 ] klij = [ H 0 ] ijlk .
Inverse stress to strain relationship. W hen d 6 = 0,(the case d = 0 i s strai ghtforward),repl aci ng " h i n Eq.5 by i ts val ue (Eq.4)gi ves, rst," d asthe sol uti onsofa thi rd orderequati on and,second,the rel ati on between " h and " d :

(" d ) 3 + p" d + q = 0; p = 2'd h 3K (1 d) 4 ' 2 d 2 ; q = 3K d 8 2 ' 2 d 2 ; " h = h + 2 'd " d 2 3K : (11)
T he di scri m i nant = (p=3) 3 + (q=2) 2 l eads to anal yse from one to three real sol uti ons. T hi s non uni vocali nverse rel ati onshi p woul d not be used i n ni te el em ent cal cul ati on,but we shal lhave to consi der these three possi bi l i ti es for the sem i -anal yti calresol uti ons i n Sec.5 and 6.

3 Free energy analysis

T he proposed H el m hol tz non convex free energy exhi bi tscom pl ex behavi orthat has to be careful l y studi ed. Fi g.1 show s the contour pl ot of the free energy de ned by Eq. 1. T he transform ati on A PC Z i s the com pressi on studi ed i n Sec.5 and the dam age val ue retai ned for the draw i ng i s d(C ), the m axi m um reached duri ng thi s transform ati on.

Fi gure 1: Iso-val ues ofthe free energy i n the D rucker-Pragerpl ane for d= 0. 64 P ositivity ofthe free energy. Itsi gni esthatthe m ateri alcannotrestore m oreenergy than hasbeen stored i nsi de.In l i nearel asti ci ty,i ti sassoci ated w i th the posi ti veness of the K el vi n' s m odul i(3K ;2 ) and i nfers cl assi c bounds for the Poi sson' srati o ( 1 6 6 0: 5). In the generalcase,from Eq.1 we have:

< 0 , " d > s 3K 2~ " h : (12)
T hi s condi ti on corresponds to the forbi dden dom ai n i nsi de the curve = 0 i n Fi g.1. It onl y exi sts w hen ~ < 0. It w i l lbe show n i n Sec.3 that i t cannot be reached duri ng a transform ati on because l ocal i zati on occurs at l east w hen ~ = 0. P articular lines. Li nesofi nterestare the l ociof d = 0 and h = 0. From Eq.5 we have:

" h = 1 d 2'd ) d = 0: (13)
From Eq.9,thi s corresponds to ~ = 0. O n the ri ght si de ofthi s l i ne on Fi g.1, ~ < 0. T he second sol uti on for d = 0 i s m ore cl assi cal l y " d = 0. T he l ocus of h = 0 (Fi g.1) i s gi ven by Eq.4:

" d = s 3K " h 2 'd ) h = 0: (14)
B etween the devi atori c axi s " h = 0 and thi s l i ne, " h > 0 and h < 0. T hi s represents the di l atancy e ect ofthe m odel : the devi atori c strai n i nduces di l atancy even i fthe m ateri ali s under (m oderated) pressure. T he saddl e poi nt H i s de ned as the i ntersecti on ofthe two l i nes h = 0 and d = 0,we have:

" h (H )= 1 d 2'd ; " d (H )= 1 2'd s 3K (1 d); (H )= 3K (1 d) 2 8' 2 d 2 : (15)
A s (H ) = 0, i t corresponds to an unstabl e free stress state w i th "(H ) 6 = 0. Fi nal l y,i tcan be easi l y show n thatthe poi ntN ,correspondi ng to the m i ni m um devi atori c strai n for = 0 i s such as:

" h (N )= 2" h (H )," d (N )= p 2" d (H
). P ositivity of the energy release rate. T he energy rel ease rate Y i s associ ated w i th the dam age d:

Y = d dd = (1 + 2'" h )(" d ) 2 : (16)
T he therm aldi ssi pati on i s _

q d i = Y _ d.
A s dam age d cannot physi cal l y decrease, i ts posi ti veness i m pl i es that i fY i s negati ve,the evol uti on ofthe dam age m ust stop:

" h < 1 2' ) _ d = 0: (17)
In the exam pl e ofSec.5,' = 160 then thi s condi ti on i s " h < 3: 1210 3 . T hat regi on, denoted as _ d = 0 on Fi g. 1 corresponds to hi gh con nem ents out of m ost practi calappl i cati ons (far to be reached i n the presented exam pl es); i t i s strongl y due to the i ni ti alchoi ce ofdam age acti ng onl y on the hydrostati c part. O ne can rem ark that som e experi m ents tend to con rm the exi stence of a dam age l ocki ng at hi gh con nem ents [START_REF] Url I On | C om paction des b etons : el em ents de m od el isation et caract erisation exp erim ental e. Phd thesi s[END_REF] .

Localization. T he H adam art and R i ce cri teri on [START_REF] W R Udni Ckiand J R I Ce | C ondi ti ons for the l ocal i zati on ofdeform ati on i n pressure-sensi ti ve di l atantm ateri al s[END_REF][START_REF] En E B I L L Adon | O n l ocal i zati on i n ducti l e-bri ttl e m ateri al s under com pressi ve l oadi ngs[END_REF]ofl ocal i zati on authori zes the exi stence ofa l ocal i zati on pl ane (a concentrati on ofl arge strai ns) orthogonalto the uni t vectorñ as soon as: det(ñ: H : ñ)6 0:

A l though we sti l lconsi der here el asti c transform ati ons,we assum e thi s l ocali zati on to correspond to the appari ti on ofa m acroscopi c crack.O nce appeared, the dam age m odelceases to appl y as the body i s spl i t i nto pi eces: i t i nfers a restri cti on to the dom ai n ofde ni ti on ofthe m odel . A m ong possi bl e vectorsñ, we consi derñ I ,an ei genvector ofthe strai n devi ator tensor,then " d : ñI = " d I ñI , w here " d I i s the correspondi ng ei genval ue of" d (a pri nci paldevi atori c strai n). A s H = H 0 i n thi s case,Eq.8 gi ve:

ñI : H 0 : ñI = 3K + ~ 3 2 " d I ñI ñI + ~ I: (19)
T he determ i nant ofthi s expressi on i s:

det(ñ I : H 0 : ñI )= ~ 2 3K + 4~ 3 2 " d I : (20)
T hi s show s that ~ = 0 i s a su ci ent condi ti on of l ocal i zati on. T he dom ai n (grayed) on the ri ght si de of thi s l i ne i n Fi g. 1 w i l lnot be concerned by the possi bl e evol uti ons. From Eq.12,i t al ways contai ns the dom ai n 6 0. From Eq.63," d I > " d = p 6 > 0,consi deri ng that K , ~ and are posi ti ve and usi ng Eq.9 and 10,we obtai n a second su ci ent condi ti on for l ocal i zati on:

" d > 3K + 4 (1 d) 4 p 2 'd p 2" h : (21)
T he correspondi ng l i ne i s denoted as on the Fi g.1. T he l ocal i zati on occurs at l east w hen a transform ati on reaches the l i nes or ~ = 0;the dom ai n above them (grayed)i s not concerned by possi bl e evol uti ons.

O ther rem arks about stability. T he currentl i ne ! passes by the saddl e poi nt H .From i ts de ni ti on _

" h = h = _ " d = d com es: 2 (1 d)" d _ " h 4 'd" h " d _ " h = 3K " h _ " d 2 'd(" d ) 2 _ " d : (22)
T hi sexpressi on hasno si m pl e anal yti c sol uti on and ! hasbeen draw n on Fi g.1 w i th a steepest sl ope al gori thm . B el ow !, the stress tends to bri ng back the system to the stabl e ori gi nalstate 0;above ! i ttendsto m ake the system reach the dom ai n above the l i ne w here the l ocal i zati on occurs. T hen the dom ai n between ! and can be consi dered as unstabl e, l eadi ng to an i nstantaneous evol uti on towards the l i nes and ~ = 0 w here the l ocal i zati on occurs. by [START_REF] Coi | A new yi el d cri teri on for the concrete m ateri al s[END_REF] . Its shape,sm ooth and convex,can be regarded as a softened approxim ati on ofthe Von M i ses and R anki ne cri teri ons. It i s i n good agreem ent w i th the el asti c l i m i t i denti ed by [ 3]i n bi axi altesti ngs (Fi g.2). T he rst m em ber correspondsto the Von M i sesexpressi on;the second one,due to a tensorexponenti al ,grow thsqui ckl y forposi ti ve val uesofany pri nci palstress.T he m ateri al constants are y ,that pri nci pal l y rul es the l i m i t stress i n com pressi on and 0 , that strongl y i n uences the tensi on to com pressi on stress rati o. T he functi ons g(d) and h(d),w i th g(0) = h(0) = 1,are hardeni ng functi ons (w hose detai l ed expressi ons are gi ven by Eq.47 and 48).

D issipative behaviour

f( ;d)= j j d j j g(d) + 0 h(d) exp 0 p 3 y : (23)
N am i ng ( I ; II ; III ) the ei genval uesofthe stress tensorand consi deri ng sound m ateri al(d = 0),thi s cri teri on w ri tes:

f( I ; II ; III ;d = 0) = r ( II III ) 2 + ( III I ) 2 + ( I II ) 2 3 + 0 s exp 2 I 0 + exp 2 II 0 + exp 2 III 0 p 3 ! y : (24)
T he i denti cati on of( 0 ; y ) i s done w i th respect to the m easured el asti c l i m i t stresses c = 9: 98 M Pa i n pure com pressi on and t = 1: 29 M Pa i n pure tensi on. T he di erence f( t ;0;0;0) f( c ;0;0;0)gi ves:

0 s exp 2 t 0 + 2 + t r 2 3 = 0 s exp 2 c 0 + 2 c r 2 3 : (25)
T hi s equati on accepts a num eri c resol uti on gi ven by the i ntersecti on ofthe two functi ons of 0 that correspond to l eft and ri ght m em bers. T hi s i ntersecti on i s uni que because each one i s m onotoni c. T he val ue of y i s obtai ned from f( t ;0;0;0)= 0 or f( c ;0;0;0)= 0. N um eri calval ues of( 0 ; y ) are reported i n tabl e (1).

Tangent m odulus and evolution equation. In case ofdi ssi pati ve transform ati on,the consi stency equati on df( ;d) = 0 and the di erenti alof (";d) gi ve the dam age evol uti on dd:

dd = @f @d + @f @ : @ @d 1 @f @ :H 0 :d": (26)
and the expressi on ofthe tangent m odul us H = d =d":

H = " I+ @f @d 1 @ @d @f @ # 1 :H 0 : (27)
T he i nverse of the fourth rank tensor can easi l y be obtai ned i n a conveni ent tensori albase [START_REF] Ski | O n hooke' s l aw[END_REF]asthe i nverse ofa 6x6 square m atri x.T he i nvol ved deri vati ves are obtai ned w i thout di cul ty from Eq.2,3 and 23:

@f @d = j j d j j g 0 (d) g 2 (d) 0 h 0 (d) h 2 (d) exp 0 p 3 ; (28) @f @ = d g(d)j j d j j + 1 h(d) exp(2 = 0 ) j jexp( = 0 )j j ; (29)
@ @d = 2 '(" d ) 2 I p 3 2 (1 + 2'" h )" d : (30)
T he el asti ci ty equati ons(2,3),the yi el d condi ti on [START_REF] Truski Novsky | T herm odynam i cs of ratei ndependant pl asti ci ty[END_REF]and the evol uti on equati ons(26 to 30)consti tute thesetthati srequi red forthe useofthem odelw i thi n a FEM code.

Tension and com pression sim ulation

Stress to strain relationship. T he stress tensor expresses as = ẽ1 ẽ1 . D ue to the i sotropy ofboth the m ateri aland the m odel ,the strai n tensorw ri tes:

" = 2 4 " 0 0 0 ~ " 0 0 0 ~ " 3 5 : (31)
w here " i sthe uni -di m ensi onalstrai n and ~ the apparentPoi sson' srati o.U si ng the projecti on equati ons(57,58)and denoti ng D = " d =j j " d j j(see A ppendi x A ), we have: (32)

D = si gn("(1 + ~ )) 2
" h = 1 2~ p 3 "; " d = r 2 3 j "(1 + ~ )j: (33)
U si ng these equati ons i n Eq.4 and 5 l eads to two expressi ons for the stress to strai n rel ati onshi p:

= 3K (1 2~ )" 4 'd p 3 (1 + ~ ) 2 " 2 ; (34) = 2 1 d 2'd 1 2~ p 3 " (1 + ~ )": (35)
El i m i nati ng between them gi ves a second order equati on i n ~ :

~ 2 12 'd p 3 "+ ~ 2 (1 d)+ 6K + 12 'd" p 3 + 2 (1 d) 3K = 0: (36)
Eq. 34 or 35 gi ves from zero to two sol uti ons for a gi ven strai n " (these sol uti ons correspond to di erent val ues of~ ,then to di erent strai n tensors ").

D issipative transform ations equations. T he yi el d functi on (23) becom es:

f ( ;d)= r 2 3 j j g(d) + 0 h(d) s exp 2 0 + 2 p 3 ! y : (37)
A s i n the general3D case,from the di erenti alof (";d) and the consi stency equati on df ( ;d)= 0 i s obtai ned the stress to strai n di erenti alrel ati onshi p:

d @f @d + @f @ @ @d = @f @d @ @" d": (38)
From Eq.37 we have:

@f @ = r 2 3 si gn( ) g(d) + exp 2 0 h(d) r 2 + exp 2 0 ; (39)
@f @d = r 2 3 j j g 0 (d) g 2 (d) 0 h 0 (d) h 2 (d) s 2 + exp 2 0 p 3 ! : (40)
From the stress to strai n rel ati onshi p (Eq.34, 35 and 36),the fol l ow i ng deri vati ves,i n w hi ch @ =@" i sthe tangentYoung m odul usand the di scri m i nant(i n Eq.36),are obtai ned:

@ @" = 3K (1 2~ ) 8 p 3 'd(1 + ~ ) 2 " 3K + 4 p 3 'd(1 + ~ )" 2" @~ @" ; (41) 3 2 "(1 + ~ ) @ @d = 2' p 3 2" @~ @" 1 + 2~ (1 + ~ )" 1 + 2' 1 2~ p 3 " 2(1 + ~ )+ 2" @~ @" ; (42) 2" @~ @" = (1 + 2~ ) 2 (1 d)+ 12 'd"= p 3 p : (43)
Localization in tension and com pression. In thetensi on orcom pressi on state,the axi s 1 i s the sym m etry axi s. A s axi s 2 and 3 have an equi val ent rol e, we can de ne w i th no restri cti on ñ i n the pl ane [ ẽ1 ;ẽ 2 ] as: ñ = cos( )ẽ 1 + si n( )ẽ 2 . For di ssi pati ve transform ati ons,the generalequati ons ofSec.4 al l ow to com putedet(ñ: H : ñ)num eri cal l y by searchi ng i tsm i ni m um val uew i th respect to . T he non convex potenti alauthori zes possi bl e l ocal i zati on duri ng el asti c transform ati ons.From Eq.8 we have:

ñ: H 0 : ñ = 2 4 b+ (a + 4c)cos 2 ( ) (a + c)si n( )cos( ) 0 (a + c)si n( )cos( ) b+ (a 2c)si n 2 ( ) 0 0 0 b 3 5 ; (44)
a = K + 3 (1 d) 2 3 p 3 'd(1 2~ )"; b= (1 d) 2 p 3 'd(1 2~ )"; c = 4 3 p 3 'd(1+ ~ )": (45)
Identi cation. T hei denti cati on oftheYoung M odul usE = 31: 2 G Pa and the Poi sson' s rati o = 0: 162 bei ng very cl assi c,they w i l lnot be detai l ed;the correspondi ng bul k and shear m odul i(K ; ) are gi ven i n tabl e (1). T he procedure ofi denti cati on of( 0 ; y )i sgi ven i n Sec.4.El i m i nati ng d between Eq.34 and 35 gi ves the fol l ow i ng expressi on that can be used w i th the experi m ental data i n order to determ i ne ':

' = p 3 2 =" 3K (1 2~ ) 3K "(1 2~ ) 2 2 "(1 + ~ ) 2 + ~ : (46)
T he hardeni ng functi ons g(d) and h(d) are i denti ed w i th respect to the transform ati on A C (see Fi g.3). T hei r structure has been chosen i n order to have a sm ooth transi ti on at the yi el di ng poi nt A . U niaxialcom pression sim ulation. T hesetofequati ons(31 to 43)de nes the strai n dri ven com pressi on and tensi on curves (Fi g.3 and 4),al l ow i ng the com pari son w i th K upfer' s data (ci rcl es). T he com pressi on starts by the l i near el asti c transform ati on O A ,w i th d = 0. T he m ateri albegi ns to yi el d atpoi nt A w here (A )= c . T he transform ati on A B C correspondsto the dam age grow th, up to poi nt C w here d(C )= 0: 64. T he l ocal i zati on cri teri on i s never reached duri ng A B C .A t poi nt P, ~ = 0: 5: the m ateri albegi ns to di l ate. U nl oadi ngs such asO B O ford = 0: 5 and O C O ford = 0: 64 revealweak nonl i nearel asti ci ty: i tcan benoti ced thatthecurvatureofO C (thetangentm odul usdecreasesasthe l oad i ncreases)i squal i tati vel y i n good agreem entw i th theexperi m entsof [START_REF] Sal Ah R Am | C ontribution a l a m od el isation du com portm entm ul tiaxial du b eton endom m ag e avec description du caract ere unil at eral[END_REF]for exam pl e. O fcourse hysteresi s l oops or perm anent strai ns are not descri bed by the presentsi m pl e form ofthe m odelthatdoesnottakesi nto accountpl asti ci ty e ects. Poi ntC correspondsto a hori zontalsl ope for the transform ati on.T hen d = 0 and,asthe val ue of@f=@d 6 = 0 atthi spoi nt,Eq.38 show sthat@ =@" = 0: the el asti c transform ati on O C Z has a al so an hori zontaltangent at poi nt C . A fterC ,@ =@" < 0,the transform ati on C Z i sa non l i nearel asti csofteni ng.T he dam age rem ai ns xed at d(C ). From som e poi nt just before Z,an unl oadi ng fol l ow s the unreal i sti c el asti c transform ati on ZC O (thi s aspect of the m odel rem ai ns to be enhanced). A t poi nt Z the l ocal i zati on condi ti on rel ati ve to the el asti c transform ati ons (Eq.45), i s reached and thi s poi nt i s cl ose to the Fi gure 3: U ni axi alcom pressi on si m ul ati on experi m ental breakdow n poi nt i denti ed by K upfer (al though thi s val ue has notbeen used duri ng the i denti cati on process).T he conti nuati on ofthe el asti c transform ati on after Z i s, for i nform ati on,draw n on Fi g.1; as " h and " d are dependent upon " by Eq.33,they becom e si m ul taneousl y equalto zero at the saddl e poi ntH .T he l ocal i zati on ari sesatpoi ntZ forthe angl e c = 56 degrees: the crack pl ane i s rather cl ose to the l oad axi s ẽ1 , as cl assi cal l y observed i n experi m ents. T he apparent Poi sson' s rati o ~ reaches 0. 70 at the poi nt Z and thi sl argeval ue,thatcannotbereached by a cl assi cdam agem odel ,i sresponsi bl e ofthe good descri pti on ofthe transverse strai n.

g(d) = 1 + g 0 p d; (47) h(d) = 1 d + h 0 ( p d d): (48
U niaxial tension sim ulation. T he tensi on curve (Fi g.4) i s al so referred to poi nts A , B , C ' , Z w hose rol es are si m i l ar to the ones used i n com pressi on. T he yi el di ng occurs at poi nt A ,at the i denti ed val ue t that i s,accordi ng to K upfer and con rm ed by [START_REF]A cousti c em i ssi on and post-cri ti calm echani calbehavi our ofa concrete under tensi l e stress[END_REF] ,l ess than the hal fofthe peak stress: thi s resul t i s di erent from the cl assi c vi si on ofa fragi l e behavi or i n tensi on. A t poi nt C ' , @f=@d = 0 then (from Eq.38)the sl ope i s hori zontal .Iti s notreached because the l ocal i zati on occursbefore,atpoi ntZ,duri ng the di ssi pati ve transform ati on. If [START_REF]A cousti c em i ssi on and post-cri ti calm echani calbehavi our ofa concrete under tensi l e stress[END_REF]descri besa post-peak evol uti on,K upferdoesnot,nei therthi sm odelw i th the retai ned constants. A gai n,al though not used i n the i denti cati on process, the stress (Z) i s cl ose to the experi m entalval ue. T he angl e ofl ocal i zati on at poi ntZ i s t = 27 degrees,l eadi ng to a crack pl ane rathercl ose to be orthogonal to the tensi l e axi sẽ 1 ,ascom m onl y observed (m oreoverm ostofdam age m odel s suppose a crack at 0 degree). T he l ateralstrai n, due to l ow dam age val ues, rem ai ns very sm al l , consi stentl y w i th experi m ents. T he dam age val ue at the l ocal i zati on poi nt d(Z)= 0: 22 i s weak: al though i sotropi c,thi s m odelal l ow s to have a rupture i n tensi on and to keep a l oadi ng capaci ty i n com pressi on.

Fi gure 4: U ni axi altensi on si m ul ati on

E qui-biaxial sim ulations

In thi s secti on we consi der an equi -bi axi all oadi ng (strai n dri ven) i n order to val i date the m odel w i th the constants i denti ed i n tensi on and com pressi on. T he stresstensori s = (ẽ 2 ẽ2 + ẽ3 ẽ3 ). D ue to the i ni ti ali sotropy ofboth the m ateri aland the m odel ,the strai n tensor w ri tes:

" h = 2 ^ p 3 "; " d = r 2 3 j "(1 + ^ )j: (50)
U si ng these equati ons i n Eq. 4 and 5 gi ves two expressi ons for the stress to strai n rel ati onshi p:

= 3K 2 (2 ^ )" 2 'd p 3 (1 + ^ ) 2 " 2 ; (52) = 2 1 d 2'd 2 ^ p 3 " (1 + ^ )": (53)
El i m i nati ng gi vesthe fol l ow i ng second orderequati on w hi ch acceptsfrom zero to two sol uti ons i n ^ :

^ 2 2 p 3 'd"+ ^ 2 (1 d)+ 3K 2 + 2 (1 d) 3K 2 p 3 'd" = 0: (54)
O ther cal cul iare si m i l ar to the previ ous case. E qui-biaxial com pression sim ulation. D i rectl y dependent upon the Fi gure 5: Equi -bi axi alcom pressi on si m ul ati on yi el d surface,the yi el d poi nt A (Fi g.5) i s wel ldepi cted. T he dam age i ncreases duri ng transform ati on A Z up to d(Z)= 0: 70,thi s val ue bei ng greater than the m axi m alone obtai ned i n uni axi alcom pressi on.A tpoi ntZ,(w hi ch i sbefore the poi ntC 'w here@f=@d = 0)the l ocal i zati on condi ti on fordi ssi pati vetransform ati ons(Eq.18)i sreached.T hel ocal i zati on angl ei s cc = 40 degrees,consi stentl y w i th K upfer' s postm ortem pi cture. T he coni calcrack shape i n Fi g.5 has been chosen i n order to respect the sym m etry axi s ofthe probl em . W i thout contradi cti on,the observed structure i spyram i dal ,thi s shape al l ow i ng the ki nem ati cs ofthe concrete pi eces. T he l ocal i zati on appears too qui ckl y and the l arge di l atancy i sa l i ttl e overesti m ated.M oreover,thi spostpeak regi on di ersaccordi ng to di erentexperi m entalm ethodsand concretes:forexam pl e the hi gh strength concrete tested by [START_REF] Arzouk | B ehavi or of hi gh-strength concrete under bi axi alstresses[END_REF]exhi bi ts no post peak. E qui-biaxial tension sim ulation. T he part before the peak stress i s i n good agreem entw i th the experi m entsforboth axi aland l ateralstrai ns(Fi g.6). C ontrary to experi m ents,a di ssi pati ve post-peak evol uti on C ' Z i s descri bed (i t represents a m i nor probl em com pared to the el asti c transform ati on C Z i n the uni axi alcom pressi on case). T he dam age d(Z)= 0: 67 at ul ti m ate poi nt reaches about the sam e val ue than i n bi axi altensi on. T he l ocal i zati on angl e tt = 12: 6 degreesobtai ned i n thi s case i nvol vescracksapproxi m ati vel y orthogonalto the l oad di recti ons,very consi stentw i th K upfer' spostm ortem pi ctures.A gai n,the Fi gure 6: Equi -bi axi altensi on si m ul ati on coni calshape i s onl y the representati on ofthe sym m etry ofthe probl em as the m odelonl y de nes the angl e ofthe crack,not i ts shape.

C onclusions

T hi s study show s that the si m pl e coupl i ng i ntroduced between the hydrostati c and devi atori c partsofthe strai n i n the free energy l eadsto a num ber ofconsequences,even i fthe pri m ary e ect,the non l i near el asti ci ty,sl i ghtl y i n uences the shape of the el asti c di scharges. T he i m portant di l atancy of the concrete at ul ti m ate stages i s wel ldepi cted. T he l ocal i zati on cri teri on i s used as an i ndi cator of the creati on of a m acroscopi c crack,i t appl i es at stress and strai n l evel s that are qui te consi stent w i th experi m ents. T he dam age val ue predi cted at l ocal i zati on i n pure tensi on rem ai ns sm al l ,l etti ng to the m ateri alsom e l oad beari ng capaci ty,forexam pl e forfurthercom pressi on.T he crack angl e gi ven by the l ocal i zati on condi ti on i sdi erenti n tensi on and com pressi on and consi stent w i th observati ons.

T he descri pti on ofthese e ects consti tutes an i m provem ent w i th respect to m ost ofavai l abl e dam age m odel s ofcom parabl e com pl exi ty ( ve constants are used). T hey are parti cul arl y rel evant w hen,as i n earthquake engi neeri ng,the response of a concrete structure under severe condi ti ons i s consi dered. T he gl obalstructure of the m odel and i ts em beddi ng therm odynam i c fram ework, opens the way to num eri cali m pl em entati on w i thi n a ni te el em ent code.

T he i ntroducti on of a pl asti ci ty form al i sm [START_REF] Fr Ed Eri C R Agueneau | D am age m odelfor concrete l i ke m ateri al s coupl i ng cracki ng and fri cti on,contri buti on towards structuraldam pi ng: Fi rst uni axi alappl i cati on[END_REF][START_REF] Saeb | A coupl ed el astopl asti c dam age m odelfor geom ateri al s. C om p[END_REF]woul d hel p to descri be m i ssi ng e ects such as hysteresi s l oops and perm anent strai ns and to correct som e i rrel evant post-peak responses. T he hydrostati c part of the free energy coul d be a ected by dam age i n order to extent the el d ofappl i cati on ofthe m odelto con ned statesand avoi d the dam age l ocki ng condi ti on presenti n the actualform . T he three-di m ensi onnalstress to strai n rel ati onshi p i s uni vocal , however the i nverse one i s not: i nteresti ng theoreti calanal ysi s can be done on stressdri ven l oadi ngsasthey m ay i nvol ve m ul tiphased states(thi srefersto the work of [START_REF] Ksen | Equi l i bri um ofbars[END_REF]i n el asti ci ty and [START_REF] Arfagni | A m echani calm odel for the el asti c-pasti c behavi or of m etal l i c bars[END_REF][START_REF] Truski Novsky | T herm odynam i cs of ratei ndependant pl asti ci ty[END_REF]i n pl asti ci ty).

A ppendices A Tensorialform alism

For any sym m etri c second order tensor A the hydrostati c and devi atori c parts are denoted respecti vel y by A h and A d . T he hydrostati c part i s obtai ned by the projecti on onto the norm ed hydrostati ctensorI= p 3 (w here I i sthe i denti ty tensor).

A = A h + A d (55)
A h = A : I p 3

I p 3 (56)
T he sym bol ": " denotes the doubl e contracti on i.e. A : I = A ij ji = tr(A ).

In thi s m odel ,the devi atori c parts d and " d (ofthe C auchy stress and the i n ni tesi m al strai n ") w i l l rem ai n col l i near but m ay be a priori of opposi te di recti ons.W e choose to de ne thi sdi recti on w i th respectto the norm ed tensor D gi ven by: D = " d j j " d j j (57)

T he norm used i s the eucl i dean naturalone i.e. j j A j j= p A ij A ij . T he al gebrai c val ues (" h ;" d ; h ; d ) are de ned as: O ne can rem ark that " d i s posi ti ve but d can a priori ei ther be posi ti ve or negati ve. T hese projecti ons are know n under the nam es p and q w here p = h = p 3 and q = p 3=2j d ji n soi l s m echani cs.

B U sefulproperty of the deviator

" d I > " d II > " d III (61)
T hi s l eads to the fol l ow i ng expressi ons: [ 2] P Jam et, A l ai n M i l l ard, and G N ahas. Tri axi al behavi our of a m i croconcrete com pl ete stress-strai n for con ni ng pressures rangi ng from 0 to 100 M Pa.In Proc.ofInt.C onf.on C oncrete under M ul tiaxialC onditions, vol um e 1,pages133{140,Toul ouse,France,1984.U ni versi t ePaulSabati er.

[ 3] H el m utK upfer,H ubertK H i l sdorf,and H ubertR usch.B ehavi orofconcrete under bi axi alstresses.A C I Journal ,66: 656{666,1969.

Y

  ield surface. T hi s yi el d surface i s based on the el asti c cri teri on proposed Fi gure 2: Ini ti alyi el d surface (M Pa)

  ) (G Pa) (M Pa) (M Pa)

  " 6 0. T he projecti on equati ons (57,58) gi ve:D = si gn("(1 + ^ ))

" h = " : I p 3 ; h = : I p 3 ;

 33 " d = " :D ; d = :D (58)

6 ;" d = p 6 ]

 66 d ) 2 3(" d II ) 2 (62) A s soon as " d II evol ves between i ts bounds [ " d = p (gi ven respecti vel y by " d II = " d III and " d II = " d I ),these functi ons are m onotoni c and: A nna O ugl ova,Y vesB erthaud,M arcFran coi s,and Fran coi sFoct.M echani cal properti es of ferri c oxi de form ed by corrosi on i n rei nforced concrete structures. C orrosion Sci.,48: 3988{4000,2006.

  T he sorted ei genval ues of" d are (" d I ;" d II ;" d III ). From Eq.58,we have:

	(" d I ) 2 + (" d II ) 2 + (" d III ) 2 = (" d ) 2	(59)
	" d I + " d II + " d III = 0	(60)


Nicolas Debarsy
email: ndebarsy@fundp.ac.be.

Vincenzo Verardi
email: vverardi@fundp.ac.be.

Estimating Nonlinearities in Spatial Autoregressive Models *

Keywords: Spatial econometrics, semiparametric estimations JEL Classification: C14, C21

In spatial autoregressive models, the functional form of autocorrelation is assumed to be linear. In this paper, we propose a simple semiparametric procedure, based on Yatchew's (1998) partial linear least squares, that relaxes this restriction. Simple simulations show that this model outperforms traditional SAR estimation when nonlinearities are present. We then apply the methodology on real data to test for the spatial pattern of voting for independent candidates in US presidential elections. We find that in some counties, votes for "third candidates" are non-linearly related to votes for "third candidates" in neighboring counties, which pleads for strategic behavior.

Introduction

In spatial econometrics, autoregressive models (or SAR) have been developed to estimate how changes in a given variable spread to the neighborhood (see Anselin, 1988 or LeSage and Pace, 2009, for further details). This effect is generally assumed to be linear, which is restrictive. The objective of this note is to propose a simple estimation method, based on Yatchew's (1998) difference estimator, that allows for non-linear spatial interdependence.

The structure of the paper is the following: after this short introduction we propose in section 2 a procedure to estimate a semiparametric spatial autoregressive model. Furthermore we recommend a test that allows to understand if assuming linearity in the autoregressive component is legitimate. In section 3 we present some simple simulations to assess the performance of the procedure and we present an empirical example based on US presidential elections in section 4. Section 5 concludes.

Estimation method 2.1 Nonlinear spatial autoregressive model

The general form of a linear first order spatial autoregressive model is

y i = x i θ + ρW y i + ε i , i = 1, ..., N (1)
where y i is the value taken by the dependent variable and x i is the row vector of characteristics of individual i. W y i measures the value of y in the neighborhood N (i) of individual i and is defined as

W y i = j∈N (i)
w ij y j where w ij models spatial interactions between i and j. Column vector θ and the autoregressive spatial parameter ρ are the coefficients to be estimated. Given the linearity of (1), a unit change in W y is associated to a ρ units change in the conditional expectation of y, whatever the value of W y. This assumption could be relaxed by considering a more general model of the type

y i = x i θ + f (W y i ) + ε i , i = 1, ..., N (2)
where f is any bounded and differentiable function. This model can be easily estimated extending Yatchew's (1998) semiparametric procedure to the case of a spatial autoregressive model by assuming (i) that W y is drawn from a continuous distribution with a convex finite support, (ii) that the first derivative of f is bounded by a constant L and (iii) that the ε i are independent and identically distributed (i.i.d) with mean 0 and variance σ 2 ε . It is well-known in the spatial econometric literature that W y is endogenous. To deal with this, we suggest using Newey et al. (1999) procedure which consists in introducing the estimated residuals (η i ) of a first step equation where the troublesome variable W y is regressed on the exogenous variables in (2) and a set of instruments. The instruments we use are standard in the spatial econometric literature and are the aggregate value of the exogenous variables measured in the neighborhood (i.e. W x and W 2 x). To simplify notation, we consider ηi as an additional exogenous variable in the matrix x in (2). If the coefficient associated to ηi is not significant, the endogeneity bias is negligible.

The rationale underlying the semi-parametric method is simple: suppose that we rearrange the observations by sorting them in increasing order according to variable W y. By taking first differences, we have:

y i -y i-1 = (z i -z i-1 )δ dif f + [f (W y i ) -f (W y i-1 )] + (ε i -ε i-1 ) , i = 2, ..., N (3) with z i = [x i η i ], δ ′ dif f = [θ ′ dif f γ dif f ] and γ dif f the coefficient associated to (η i -η i-1
). Increasing the number of observations (which, broadly speaking, means filling the finite support interval of W y with values) will cause the difference W y i -W y i-1 to shrink at a rate of 1/N. Since the first derivative of f is assumed bounded by L,

|f (W y i ) -f (W y i-1 )| ≤ L |W y i -W y i-1 |. The shrinkage of (W y i -W y i-1 ) will induce f (W y i-1
) to cancel out with f (W y i ). Reordering and taking differences thus allows to estimate the θ vector of parameters consistently with LS, whatever the functional form of f , as soon as its first derivative is bounded.

Note that this simple estimator is inefficient (it has a Gaussian efficiency of 66.7%).

To increase efficiency, Yatchew (1998) suggests to use higher order differences. In our setup, this generalization is written as:

Dm j=0 d j y t-j = Dm j=0 d j z i-j δ dif f + Dm j=0 d j f (W y i-j ) + Dm j=0 d j ε i-j , i = D m + 1, ..., N (4
d 2 j = 1,
guarantees that the variance of the transformed residual in ( 4) is σ 2 ε . Yatchew (1998) shows that, with D m sufficiently large, the estimator approaches asymptotic efficiency. [START_REF] Anselin | Spatial Econometrics, Methods and Models[END_REF] As far as inference is concerned, Yatchew (1998) shows that θdiff has the approxi-

mate sampling distribution θdiff ∼ N θ, 1 + 1 2D m σ 2 ε Nσ 2 u (5)
where σ 2 u is the conditional variance of z given W y. The standard errors of the estimated parameters can thus be easily computed. Concerning the relevance of the nonlinear effect of the variable W y , Yatchew (1998) developed a simple test based on the comparison of the scale of the residuals of the difference equation (s 2 dif f ) with that of the LS regression where the function f is assumed to be linear (s 2 lin ). The underlying idea of the test is that if nonlinearities exist, a linear approximation of the tested relation will lead to an overestimation of the variance of the residuals.

The proposed test statistic is:

V = √ D m N (s 2 lin -s 2 dif f ) s 2 dif f (6) 1 When D m = 1, d 0 = 1 √ 2 and d 1 = -1 √ 2 , equation (4) boils down to equation (3).
which is asymptotically distributed as a N(0, 1). A rejection of the null would suggest a non-linear relation between W y and y.

In our setup, the estimation of s 2 lin cannot be based on LS since it is biased and inconsistent when an autoregressive term is present. We therefore estimate it using the residuals of ( 1), the SAR model, estimated by maximum likelihood.

Finally, since θdiff is a consistent estimator of θ, the relation between y and W y can be assessed by running a non-parametric estimation of the partialled out ỹi (which is ỹi = y ix i θdiff ) on W y. This variable still contains the information on the spatial dependence of interest, but is filtered of the influence of the control variables. The nonparametric estimator we consider is cubic spline. Alternative such as kernel regression methods can be used. They lead to very similar results.

To assess the performance of the proposed methodology, we present some simple simulations in the following section.

Simulations

The four following data generating processes (DGP) are considered:

a) y i = x i θ + ε i b) y i = 0.75W y i + x i θ + ε i c) y i = 0.75W y i -0.4 (W y i ) 2 + x i θ + ε i d) y i = 1 1+exp(-2W y i ) -0.5 + x i θ + ε i
where x i is a 1 × 3 vector whose elements are drawn from three independent N(0, 4), θ is a 3 × 1 vector of ones and ε i ∼ N(0, 0.1). The simulated sample size is 300. The xcoordinates are generated from a U[0, 20] and the y-coordinates from a U[0, 50]. Spatial weights are

w ij =          1/b ij j 1/b ij if b ij < b 0 otherwise
where b ij are all pairwise distances. Parameter b (the threshold value above which the interaction between i and j is assumed to be negligible) is set to 5. By convention, w ii = 0. All the models are fitted assuming D m = 1.

To illustrate the fitting performance of the proposed estimation procedure, we generate four samples according to the DGPs discussed above and present the scatter plots, the non-parametric fit (thick plain line) and the true DGP (thin dashed line) in Figure 1. As expected, the results are unambiguous.

[INSERT FIGURE 1 HERE]

In the case of no spatial autocorrelation (panel a), no clear pattern emerges and the non-parametric curve lies close to the horizontal line (the true DGP). In the three other cases (panels b, c and d), the nonparametric estimation of the autocorrelation matches the true functional form quite well. The last two panels (c and d) shed doubt on the appropriateness of a linear approximation for the spatial component.

To check the performance of the V-statistic in detecting nonlinearities, we replicate the four DGPs described above 1000 times. Each time a new error term is randomly drawn and a new dependent variable is generated. The design space is kept unchanged.

We then compute the percentage of rejection of the null (at 5%). Results are presented in Table 1. The null is always rejected in the quadratic and the sigmoid configurations. The size of the test is about 5% for the absent and linear cases. Though we are aware that these simulations are too simple to provide a clear assessment of the quality and power of the test, it seems that it could be an interesting complementary tool to assess the linearity of potential spatial autocorrelation.

Application

In this section, we present an illustrative example of the nonlinear SAR model. The objective of the analysis is to study the voting behavior for independent candidates (focusing on US presidential elections of 2000) in a given county as a function of the votes cast by this candidate in neighboring counties (which are assumed to be well anticipated by electors). The hypothesis tested is that voters will not vote for the third candidate if they believe that it might help the candidate they dislike the most win the elections (as occurred, for instance, in the first round of the French presidential elections of 2002 when Le Pen, the extreme right wing candidate, obtained 16.86% of the votes and qualified for the second round at the expense of the socialist candidate).

Hence, if interested voters anticipate that the third candidate will collect a limited number of votes, they will vote for him to declare their dissatisfaction with the political establishment. Furthermore, the resulting share of votes is expected to increase jointly with the share of votes in the neighborhood as the message sent will be stronger.

However, if voters perceive that the third candidate will obtain so many votes that it results in jeopardizing the political scenario they will stop voting for him. We therefore expect to observe a concave-shaped spatial autoregressive component in the vote for the third candidate. To test for this, we estimate the relation

y i = x i θ + f (W y i ) + ε i for i = 1, ..., N (7)
where variable y is the log of the vote share cast by outsider candidates and i indexes counties. The control variables (x) are those generally considered in this type of regressions i.e. (i) the log of the votes share of independents in the previous elections, (ii) the total population, (iii) the log of the vote share for Republicans, (iv) the average per capita income, (v) the log of the ratio between democrats and republican votes for the previous elections and (vi) the ethnic composition (i.e. the proportion of blacks, whites, asians, native indians others and pacific island which is the reference). [START_REF] Lesage | Introduction to Spatial Econometrics[END_REF] The order of difference considered is 10 and the weighting matrix is defined as follows:

counties located in different states do not interact. [START_REF] Newey | Nonparametric estimation of triangular simultaneous equation models[END_REF] Within states, we assign a spatial weight proportional to the inverse of the distance between counties' centroids which implicitly assumes that individuals are better informed on closer counties. Data come from Polidata, a national demographic and political data consulting firm in the US.

We concentrate our analysis only on swing states (since the described theory only holds in these) [START_REF] Yatchew | Nonparametric Regression Techniques in Economics[END_REF] and focus on Oregon, Pennsylvania, Tennessee and Wisconsin where the concaved-shaped relation is clear. 5

[INSERT FIGURE 2 HERE]

For the state of Pennsylvania, we observe that the sincere voting behavior occurs as long as W y is smaller than 2.7%. Indeed, in this situation, we observe that an increase in the vote share for independents in the neighborhood induces an increase in the vote share for independents in the considered county. This is probably because voters think that their vote will strengthen the message conveyed by the neighbors on the dissatisfaction with the political establishment. However, when the votes for outsiders become more numerous in the neighborhood (around 2.7%), the "votes for change" start decreasing. This could be explained by the fact that voters realize that they should vote strategically to prevent the candidate they dislike the most to win the elections. This theory seems to hold well in Oregon and Tennessee, where the concavedshaped relation is evident. For Wisconsin, we observe a similar behavior except that 2 Including the log of the vote share for the Democrats instead does not affect results. [START_REF] Newey | Nonparametric estimation of triangular simultaneous equation models[END_REF] This assumption relies on the majoritarian system in place. [START_REF] Yatchew | Nonparametric Regression Techniques in Economics[END_REF] In 2000 swing states were Tennessee, Nevada, Ohio, Missouri, New Hampshire, Florida, New Mexico, Wisconsin, Iowa, Oregon, Minnesota and Pennsylvania. 5 The graphs for all the other States are available from the authors upon request. 8 a first threshold is observed when W y is around 2.3% and a decrease takes place at 2.45%.

As far as inference is concerned, the V-statistic is larger than 2 (in absolute value) in all States. Its value is 4.19 for Oregon, 10.10 for Pennsylvania, 25.53 for Tennessee and 40.39 for Wisconsin. This clearly rejects the linearity assumption in all the cases considered, which confirms the impression given by the graphs.

Conclusion

In spatial econometrics, the SAR is one of the most commonly used models. In this paper, we propose a simple generalization of it for the case of a non-linear spatial autoregressive component. We present some simulations and a simple empirical application to show the usefulness of the procedure.

  ) where D m (∈ N + 0 ) is the order of differencing and d 0 , ..., d m are differencing weights. Two conditions are imposed on d 0 , ..., d m . The first is m j=0 d j = 0, which ensures that the nonparametric spatial component part is partialled out. The second, m j=0

Table 1

 1

	: Test for linearity
	Spatial autocorrelation absent linear quadratic sigmoid
	% Rejection	5.2%	4.8%	100%	100%


E Côme

M Cottrell

M Verleysen

J Lacaille

Self Organizing Star (SOS) for health monitoring

Self Organizing Maps (SOM) have been successfully applied in a lot of real world hard problems since their apparition. In this paper we present new topologies for SOM based on a planar graph. The design of a specific graph to encode prior information on the dataset topology is the central question addressed in this paper. In this context, star-shaped graphs are advocated for health monitoring applications, leading to a new kind of SOM that we denote by Self Organizing Star (SOS). Experiments using aircraft engine measurements show that SOS lead to meaningful and natural dataset representation.

Introduction

SOM are well-known data mining tools which are frequently used to project high-dimensional data into a low-dimensional discrete space for visualization and exploratory purpose. In the context of dimensionality reduction, SOM can be thought as a kind of nonlinear but discrete PCA. SOM are classically defined by two key elements: a set of prototype vectors {m 1 , . . . , m K }, which are ddimensional vectors, d being the dimensionality of the dataset {x 1 , . . . , x n } and a lattice or grid which defines a neighborhood relation between the prototypes. This lattice plays the role of the embedding space. Each unit or prototype is associated with a fixed position in this lattice or grid. Self-organization means that after convergence prototypes which are close in the lattice (with a small distance between them) will also be close in the data space. New data can be projected after learning of the SOM by affecting them to the location on the map of the best matching unit (BMU): the closest prototype vector to the datapoint.

In this paper we investigate the opportunity of using other lattices different from regular rectangular or hexagonal grids defined thanks to undirected graphs, which will have the additional interesting property of being planar. By using such graphs all the interesting properties of SOM are kept, in particular easy data visualization and data exploration with the additional freedom of designing lattices more suited to the expected dataset topology. We also show how such graph can be designed in the context of a health monitoring application. The classical lattices for SOM are rectangular or hexagonal grids. Extensions of these classical topologies have already been investigated by several authors. Toroidal or spherical maps [START_REF] Wu | Spherical self-organizing map using efficient indexed geodesic data structure[END_REF][START_REF] Andreu | Selecting the toroidal self-organizing feature maps (tsofm) best organized to object recognition[END_REF], growing grid with free topology [START_REF] Fritzke | Growing cell structures -a self-organizing network for unsupervised and supervised learning[END_REF], and tree structured SOM (to speed up computational time required for BMU search) [START_REF] Pakkanen | The elvoving tree -analysis and applications[END_REF][START_REF] Campos | S-tree: self organizing tree for data clustering and online vector quantization[END_REF], among others. We follow here a path more similar to [START_REF] Barsi | Neural self-organization using graphs[END_REF] which focuses on the use of graphs to define SOM topologies. Our solution differs from [START_REF] Barsi | Neural self-organization using graphs[END_REF] by the nature of the graph being used for defining the SOM topology. We promote in particular the star topology, which leads to the method proposed here and called Self Organizing Star (SOS).

The paper is organized as follow. In Section 2 the basics of SOM will be reminded and the question of how the neighborhood relationship induced by a lattice influences the results of the algorithm will be addressed. The use of undirected graphs to define new lattices will also be investigated. In Section 3, the SOS solution proposed to deal with health monitoring problems will be advocated and described. Finally, some experimental results on real and simulated datasets will be given using such a topology; the paper ends with some concluding remarks.

Self-Organizing Maps Background

In its classical presentation [START_REF] Kohonen | Self-Organizing Maps. Information siences[END_REF][START_REF] Cottrell | Theoretical aspects of the som algorithm[END_REF], the SOM algorithm is an iterative algorithm, which iterates the two following steps over training patterns x j for computing the set of prototypes m i , i ∈ {1, . . . , K} which define the map:

• Competitive step, this step aims at finding the best matching unit (BMU)

for sample

x j : c = arg min i∈{1,...,K} ||x j -m i ||. (1)
• Cooperative step, this step aims at moving the BMU and close units, closer to the training pattern:

m i (t + 1) = m i (t) + α(t)h ci (t) [x j -m i (t)] , (2)
with t the time step, α(t) the learning rate of the algorithm and h ci (t) the neighborhood function between units c and i at time t.

Several neighborhood functions are commonly used such as:

h ci (t) = 1 (dci<σt) (3)
h ci (t) = exp(-d 2 ci /2σ 2 t ). (4)
All of them depend on a radius σ t which is classically decreasing during the learning process. These neighborhood function also depends on d ci : the distance between neuron c and i, which is determined by the map topology. During the cooperative step the coordinates of the best matching unit c are updated in order to move it closer to the training pattern. Other prototypes m i are also moved towards the training pattern according to their distance d ci from the BMU defined by the lattice. Closer units from the BMU will be more affected than the others. The two steps are iterated over the dataset during several epochs until convergence. Convergence is being monitored by moves of prototype vectors. The learning rate is classically decreased during the learning process using an inverse linear adaptation. Thanks to the cooperative steps, auto-organization is reached at the end of the algorithm. Looking at this algorithm, one can clearly see that the distance between units in the map plays a key role in the auto-organization property of the algorithm; therefore modifications of this distance will clearly have an impact on the results of the algorithm. The lattice structure can be a way to incorporate prior information concerning the dataset topology into the dimensionality reduction process. In the context of classical SOM, one assumes that the dataset topology can effectively be represented by a grid or a straight line but other hypotheses can be interesting and can advantageously be investigated. The key element that implicitly or explicitly defines the lattice is the distance between map units. It has been already noted that graph theory can be used to define such distance [START_REF] Barsi | Neural self-organization using graphs[END_REF][START_REF] Pakkanen | The elvoving tree -analysis and applications[END_REF]. In this case, map units are the nodes of a graph, and distances between them are defined as the minimum number of edges needed to reach one node starting from the other, that is the so called shortest path distance.We therefore propose to modify the SOM algorithm only by taking as input an adjacency matrix which specifies the graph topology that the user desires. All undirected graphs can theoretically be used, but a specific class is of interest: the class of planar graph, because such graphs can easily be represented in a 2-dimensional setting allowing us to supply the SOM with visualization tools, as seen in the next section.

Star topology

In health monitoring applications, there exists a particular state which corresponds to the normal operation of the component under monitoring. The aim of health monitoring is to look at departures from this normal state, because faults lead to deviation of state variables from this normal state. Furthermore, we may hope that different faults are characterized by different shapes of drifts. The classical lattices of SOM do not take into account this prior information on Fig. 1: Planar graph representation of the star topology with 7 rays and 3 units by rays the dataset topology: there is no specific unit which takes a central place on the map. Moreover the users do not easily understand the meaning of the map since the normal states can for example be represented by codebooks vectors located on an edge of the map, instead of a central position that would be more natural. To prevent this difficulty, we propose to use a star-shaped graph as shown in Figure 1. This graph has a clear natural center from which different arms or rays grow. Such graph can easily be understood by users if different rays correspond to different faults and distances to the center describe the seriousness of the problem. These graphs can be characterized by their number of rays and by the length of the rays. SOM using such topologies will be called Self-Organizing Stars (SOS).

Results

Toy dataset

The first experiment aims at demonstrating the ability of such topology to represent a star-shaped dataset. For this purpose a toy dataset was generated in a 3-dimensional space. One hundred data points have been generated using a normal distribution with mean (0, 0, 0) and standard deviation 2.I, where I denotes the identity matrix of size 3 × 3. Three other sets of one hundred data points have been drawn, each corresponding to a noisy version of the segment [0, 15] along one of the three vectors of the R 3 canonical basis. These sets were obtained by first sampling the location along the vectors using a uniform distribution between 0 and 15 and by adding isotropic noise with standard deviation equal to 1. A realization of such a dataset is presented in Figure 2 the SOS codebooks positions after learning. These results were obtained using a SOS with 3 rays of 15 units each. The training was performed using the classical inverse linear learning rate scheme of SOM algorithms and a decreasing radius between 5 and 1; the number of epochs of the algorithm was fixed to 5 and the initialization was performed using random assignments in the data range for all units. Figure 2 (b) gives a representation of map results in the form of a component plane. We can also see on this figure that one ray corresponds to the part of the dataset simulated along the first canonical vector. This two representations clearly show that such dataset can be effectively described by a SOS; moreover using the planar representation of the star graph, such a SOM can be effectively represented in a 2 dimensional way.

Health monitoring dataset

The second dataset comes from an health monitoring application which concerns aircraft engines. The data are multi-dimensional measurements on the engines (9 variables describing engines state and 7 variables describing the environment). The dataset was simulated using a thermo-dynamical model of the engine; six different sets of points were simulated, each one starting from a normal engine and gradually moving to a defective state. Such sets will be denoted as failure trajectory. Five sets concern progressive failures of different engine macro-components : fan, booster, high pressure core, high pressure turbine, low pressure turbine; the last one was generated with progressive degradation of all macro-components together. Random environmental conditions (external temperature, altitude, Mach, ...) were also simulated for all samples and finally realistic noise was added to all variables. Each one of the six sets contains 2000 points, leading to a dataset of 12000 multi-dimensional measurements. Rough engine measurements are inappropriate for direct use; they are indeed influenced by external conditions. We therefore first processed the data by a Linear Model (LM), to eliminate the effect of measured environmental conditions (external temperature, altitude, Mach, rotation speed, air levy, engine variable geometries positions) on the other 9 variables (temperatures and pressures on different locations and fuel consumption). Finally, the LM residuals are filtered using a simple moving average (15 samples) in order to reduce noise. The filtered residuals are then used as inputs to a SOS for the easy visualization of fault degradations. The map was trained using an inverse linear scheme for the learning rate, 5 epochs and a neighborhood radius decreased from 5 to 1. The topology of the map was fixed manually: it contains 6 rays (one for each trajectory) of 20 prototypes each. Figure 3 presents the projection of two trajectories on the map. As expected each fault corresponds to a ray; similar results were obtained for the other faults. Once learned, the map is therefore a very interesting tool for health monitoring as it enables the diagnosis of specific sub-component failures and the estimation of the failure gravity.

Conclusion

The proposed SOM extension (SOS) based on a star shaped planar graph gives interesting results in the health monitoring application investigated in this paper.

booster deterioration fan deterioration Fig. 3: Star-shaped map for turbo fan engine health monitoring; projection of the 2 first failure trajectories on the map; color (from white to black) and size (from small to big) of the squares indicates the position in the trajectory (from normal state to defective state).

This approach is interesting when prior information on the dataset topology is available as in health-monitoring applications. However, several questions are still open, first including the choice of the number of rays and their lengths which is left up to now to the user. Classical model selection tools can probably help for such a task. In addition, in the context of health monitoring applications we may face multiple simultaneous faults; this is not addressed by the proposed algorithm and is therefore also a topic for future work.

Fig. 2 :

 2 Fig. 2: (a) toy dataset in R 3 and SOS codebook's position after learning depicted by big black dots, (b) Variable 1 component plane; map units colors are proportional to the first value of their codebooks vectors.


Mabrouka El Guedri
email: mabrouka.elguedri@supelec.fr

Guy D ' Urso

Christian Lajaunie
email: christian.lajaunie@ensmp.fr

Gilles Fleury
email: gilles.fleury@supelec.fr

TIME-FREQUENCY CHARACTERISATION FOR ELECTRIC LOAD MONITORING

Electric utilities and consumers are increasingly interested in energy monitoring for economic and environmental reasons. A non-intrusive solution may rely on information extracted from the electric consumption measured at a centralized part of a distribution network. The problem at hands consists in the separation of the electric load into its major components. This problem of source separation from one sensor is quite tractable under certain conditions. In this work, the focus is made on the most consuming household appliance in France: the space-heating. It is a sum of an unknown number of pseudo-periodic signals embedded in the global active power. An unsupervised algorithm to determine the space-heating schedule from the global consumption based on the interpretation of the spaceheating signature in the time-frequency domain is proposed. The proposed method conjoins a time-frequency detector and a frequent itemsets extraction. First results on real data are quite satisfying.

INTRODUCTION AND BACKGROUND

The electric power industry and consumers recently face many challenges such as energy saving and greenhouse gas emissions reducing. Accurate and reliable information about the nature and the state of the electric systems will undoubtedly be helpful to meet these challenges. Actually, a good knowledge of the electric load and the targeted appliances help consumers understanding their bills and better control their consumption A non-intrusive and economical solution may rely on information extracted from the electric consumption measured at a centralized easily accessible part of a distribution network. Non-intrusive electric load monitoring has been subject to several approaches over the last twenty years. General overviews can be found in [1, 2 , 3]. The available solutions require measurements of the active and the reactive power, which carry out the finger-prints of the electric appliances. They are mostly made up of three steps. Event detection determines the appliances operating schedule. Load identification uses steady state powers and transient patterns, if available, to recognize the elementary components. In this paper, the focus is made on the most consuming enduse in France: the space-heating. The observed signal (the active power) is a sum of pseudo-periodic square waves characterized by slowly time-varying duty cycles. The problem statement is given in section 2. In section 3, the detection of locally-stationary signals embedded in a nonstationary observation in the time-frequency (T-F) domain is explained. A novel procedure of interpretation of the timefrequency content is investigated. Experimental tests have been carried out on real data provided by Electricité de France.

PROBLEM STATEMENT

In this paper, we are interested in the breakdown of the whole electric consumption given the active power y(t) available at the electricity meter (sampled at a sampling rate T e =2s) in France. In fact, this measurement is available and does not need any installation at home customers. Moreover, the recorded measurement could be sent to a centralized platform through radio transmission. Notice that only one measurement is required compared to existent solutions based on three measurements of current and voltages at higher sampling frequency. The electric load decomposition consists in the detection and the estimation of its major components. The observed signal is the sum of an unknown number of loads ݕ of the electric appliances as detailed in equation 1,

,ݐ∀ )ݐ(ݕ = ∑ ݕ )ݐ( + )ݐ(ܾ ୀଵ (1)
where b is an additive Gaussian noise. Classical methods of sources separation [START_REF] Hyvärinen | Independent component analysis: Algorithms and Applications[END_REF] are not suited for this problem. Recent techniques [START_REF] Benaroya | Wiener based source separation with HMM/GMM using a single sensor[END_REF] developed for audio signals are based on the temporal and spectral parsimony of the elementary signals. In this work, we propose a dedicated method to electric load characterization. An analysis and a characterization of the individual signals are crucial to deal with this problem. In this paper, the focus is made on the space-heating by convectors for two reasons. Firstly, it is the most consuming household appliance in France (~70% of the national consumption in the residential sector). Secondly, a convector load obviously differs from the other signals. When it is on, it generates periodic square signal with constant period (40s or 80s) and a slowly time-varying duty cycle. Figure 1 illustrates three signals of three (electronic) convectors operating simultaneously at a customer house. During this study, we analyzed the elementary signals recorded for eleven different convectors during three weeks. This data basis is certainly not sufficient, but necessary for a feasibility study. In fact, gathering such data is a costly and difficult task. Some common properties of the convectors loads have been established thanks to this analysis. Actually, we noted that a convector signal is locally stationary, where the local stationarity means that the covariance function of the studied signal can be considered as the time modulation of a stationary covariance function [START_REF] Flandrin | Temps-Fréquence[END_REF]. Let us underline that more data were gathered during the winter 2008-2009 from two different sites. The active power of each convector (almost thirteen convectors) has been measured. The signals analysis have the same properties of those used during this study. The detection of space-heating from the active power is formulated as the detection of a mixture of an unknown number of periodic locally stationary signals, where the set of possible values for periodicities is known. Since the observed signal y is non-stationary, time-frequency representations [START_REF] Flandrin | Temps-Fréquence[END_REF] are suitable tools to characterize its components. In this space, the signal components are described by structures called spectral patterns. Given that targeted signals are piecewise stationary, the appropriate TFR is the Short Time Fourier Transform (STFT) [START_REF] Flandrin | Temps-Fréquence[END_REF]. A sample of an electric load, where three convectors are simultaneously in use, and the magnitude of the corresponding STFT, namely the spectrogram, are illustrated by Figure 2. Equidistant horizontal lines correspond to the space-heating signature in the time-frequency domain, because of the periodic property of the convectors signature. The interpretation of the space-heating signature obviously becomes straightforward in this domain. Note that other time-frequency features may also be distinguished. The interpretation of the time-frequency representation content mostly consists of two processes: the TFR segmentation and the classification of the extracted patterns. Several parametric or non parametric, supervised or unsupervised methods have already been proposed in literature [START_REF] Pierson | Watershed segmentation of timefrequency images[END_REF][START_REF] Vincent | Nonstationary signals classification using time-frequency distributions[END_REF]. This work investigates a new unsupervised and semiparametric procedure, dedicated to the characterization of an unknown mixture of periodic locally-stationary signals embedded in a non-stationary observation. The segmentation of the TFR is tuned with a time-frequency detector based on a statistical model of the TFR. The classification process aims at separating the extracted spectral patterns into two families: the space-heating features and the other devices' ones. This operation is performed by a data-mining method, namely the frequent itemsets extraction. The proposed method is based on priors extracted from convectors measurements with electronic regulation. As for other convection systems, such as mechanics convectors, we will show that our method might be still appropriate since the space-heating load is locally stationary. The proposed approach to characterize the space-heating load is described in Figure 3.

DETECTION AND EXTRACTION OF PIECEWISE STATIONARY SIGNALS FROM THE SPECTROGRAM

Spectrogram segmentation based on a series of hypothesis tests

A TFR might be considered in two different ways. In [START_REF] Pierson | Watershed segmentation of timefrequency images[END_REF][START_REF] Leprettre | Extraction of pertinent subsets from time-frequency representations for detection and recognition purposes[END_REF] it is considered as an image. The segmentation is performed by active contours and mathematical morphology. These tools do not provide an automatic detection and characterization of the time-frequency features. In [START_REF] Coates | Time-frequency modeling[END_REF][START_REF] Coates | Regionally optimised timefrequency distributions using finite mixture models[END_REF], the TFR is modelled with a mixture of probability densities.

Its segmentation is based on the statistical properties of a specific time-frequency representation: the spectrogram, which is the magnitude of the STFT [START_REF] Flandrin | Temps-Fréquence[END_REF].

As the observation y is a sequence of N samples, each spectrogram coefficient is given as following (h is the analysis window of length N f < N):

ߩ ௬ ሾ݇, ߥሿ = ଵ ே ฬ∑ ℎሾ݉ -݇ሿ݉|ݕሿ݁ ିଶగ ഌ ಿ ே ୀଵ ฬ (2)
The segmentation procedure consists in labelling each T-F coefficient using a statistical model of the chosen TFR. In this work, our purpose is to provide two regions: a region which contains the spectral patterns and a region defined as the background of the TFR. The segmentation task might be performed through a hypothesis test. The method described below has been developed in the GIPSA-lab as part of studies conducted in this lab. The segmentation procedure consists in labelling each T-F coefficient using a statistical model of the chosen TFR. In this work, our purpose is to provide two regions: a region which contains the spectral patterns and a region defined as the background of the TFR. The segmentation task might be performed through a hypothesis test. The method described below has been developed in the GIPSA-lab as part of studies conducted in this lab.

Statistical hypothesis test

The detection of abrupt temporal changes is formulated as following: H 0 : The observed signal r is stationary, possibly added with a deterministic signal. It follows a known probability density 0 p , whose parameters are unknown. H 1 : Non-stationary signal of unknown probability density p 1 Lack of information about ruptures makes the definition of an optimal detector not feasible. In [START_REF] Hory | Mixtures of Chi2 distributions for the interpretation of a time frequency representation[END_REF], the authors propose a recursive detector based on the knowledge of the 0 p 's family. They proposed a new detector in a particular TFR: the spectrogram, and reformulated the previous hypothesis test in the new space. As the prior density of observations and the cost function are unknown, the decision rule's threshold might be determined for a given false alarm Pfa

Consequently, if ] [ ] , [ ν η ν ρ Pfa y k > ,
k K I I card I I k p p L k L y µ µ ν ρ ν ν ν ρ < Ε ⊂ = ∈ = ∈ = ) ( 1 ν H L is the complementary subspace of ) ( 0 ν H L
. Notice that the decision rule requires the knowledge of the statistical distributions of the T-F coefficients.

Statistical model of the spectrogram coefficients

For better explanation, the statistical test is presented in the particular case of a stationary white Gaussian signal. χ law [START_REF] Martin | Décision temps-fréquence[END_REF] with 2 degrees of freedom. Therefore, the spectrogram might be modelled as a mixture of central and non-central chi-square laws [START_REF] Martin | Décision temps-fréquence[END_REF]. Non-central parameters are proportional to the signal to noise ratio evaluated for the spectrogram coefficients. The estimation of the distributions' parameters helps to make a decision as for the category of each spectrogram coefficient.

Decision Rule

The threshold value is determined for a fixed false alarm Pfa . The decision rule is:

If ] [ ] , [ ν η ν ρ Pfa y k > , the alternate hypothesis 1 H is true, and ) ( ] , [ 1 ν ν ρ H y L k ∈
. Otherwise, the null hypothesis is satisfied and we have:

) ( ] , [ 0 ν ν ρ H y L k ∈ 1 H
L is the set of all the spectral patterns.

Algorithm

The key idea of this detector is to construct the subspace of interest iteratively. More precisely, let ) (

0 ν i H L be the current estimation of ) ( 0 ν H L at iteration i, { } ) ( 0 ν i H i L card P = and ] [ν µ i be the empirical mean of ) ( 0 ν i H L
. For more details, we refer to these works [START_REF] Hory | Mixtures of Chi2 distributions for the interpretation of a time frequency representation[END_REF][START_REF] Martin | Décision temps-fréquence[END_REF]. The stopping criterion is the stability of 0 H L . As for initialisation, one affects to each ) (

0 ν H L
the p% lowest values of the observed data. The parameter p is chosen a priori. This procedure is used to detect all the spectral and the temporal changes in the T-F domain. It provides a mapping of the spectral patterns including those associated to the piecewise-stationary periodic signals. The next step aims at automatically extracting these patterns so that one can select the convectors patterns.

A novel method to extract "some" time-frequency patterns: Frequent Itemsets extraction

The subspace 1 H L is a time-frequency map of two families of features: the first one is associated to periodic signals, the second one to other signals. Our objective is to extract the first category of patterns. The extraction algorithms of compact subsets [START_REF] Martin | Décision temps-fréquence[END_REF] are not suited to our problem because we aim at extracting only some itemsets having some specific properties. Actually, the features associated to the convectors have three properties. Each pattern is nearly rectangular. Besides, for a given temporal support, the convectors structures have almost the same size. Finally, whenever the space-heating is "on", the patterns of each convector are equidistant. These priors are used to model the targeted structures and to develop an automatic and unsupervised procedure of extraction. A specific method namely the Frequent Itemsets Extraction selects these itemsets. Early works on this method [START_REF]Fast Discovery of association rules[END_REF] led to the APRIORI algorithm, where the processed data are binary matrix. The rows are called objects and the columns are the attributes. Definition: A frequent itemset is an itemset for which the number of objects including this itemset, called support, is greater than a threshold named cut frequency (f c ) [START_REF] Crémilleux | Optimisation d'extraction de motifs : une nouvelle méthode fondée sur la transposition de données[END_REF].

Frequent Itemsets Extraction looks for the frequent itemsets in a large data basis. This procedure is carried on a condensed representation: a Galois lattice. The lattice construction is illustrated in Figure 5 for a non realistic binary matrix (Figure 4). The frequent itemsets for a cut frequency 2 = c f are presented in Figure 6. The APRIORI algorithm proceeds repeatedly as following: 1 srt iteration: evaluate the itemsets' frequencies k th iteration:

-Frequent itemsets of iteration (k-1) are used to determine candidates for k-frequent itemsets, with a greater frequency than the cut frequency. -Update data: pruning the non-frequent itemsets. In our case, time instants are the objects (transactions). Frequencies are the attributes. A spectral pattern could be considered as an association between an itemset of attributes and the objects (instants) containing these attributes (frequencies). Note that the matrix data must be binarised and transposed for computing reasons. Moreover, a morphological filtering is used to eliminate all the patterns with duration inferior to the minimum duration of a convector. The patterns extraction is then performed by specific software: MvMiner [START_REF] Rioult | Mining strong emerging patterns in wide sage data[END_REF]. As the space-heating spectrum covers the entire frequential axis, the extraction of the spectral signature of the convectors might be performed by selecting the longer itemset of attributes. As for the cut frequency, it could be defined regarding the minimal duration of a convector's operating 0 ch t given consumer habits an operating constraints as following: The output of this module is a map including mostly the binarised signature of the convectors in the spectrogram. We can then deduce the schedule of this end-use from a daily electrical load.

EXPERIMENTAL RESULTS

Real data have been gathered only at a customer house, because collecting real data is a difficult and costly task. The results given herein are not complete yet. We only demonstrate the efficiency of the proposed method on real data in the particular case of electronics convectors. The evaluation of the detector on large data basis for different customers will be performed in future works, using data recorded at tow houses.

Figure 1 -Figure 3 -

 13 Figure 1 -Active power of three convectors

* 0 H 1 H

 01 : Random signal b , stationary, with probability density ) : The alternate hypothesis. The statistical properties of the spectrogram coefficients under the null hypothesis have to be determined in order to evaluate the detection threshold. For each frequency, of zero mean independent Gaussian random variables. It follows asymptotically a central2

  then the alternate hypothesis H 1 is accepted. Otherwise, one accepts the null hypothesis.Let us now give some definitions and notations. For a given signal y, the observation space is ܮ = {(݇, ߥ), ݇ ∈ ,ܭ ߥ ∈ }ܨ where ܭ × ܨ is the time-frequency domain. For each normalised frequency, we define:

	where ( 0 H	: )	s	{ (	{ ( max ) arg ,	)	/	( ) ,	/	H	0	,	( 0	[ and y	,	) ]	s	} s } )
	s	is	a	threshold	chosen	a	priori.

The frequent itemsets map and the longer frequent itemsets obtained for a weekly electric load (Figure 2) are presented in Figure 7. Notice that some itemsets are split up resulting in nondetected intervals . The space-heating detector in the temporal domain (Figure 8) is thereafter obtained given the patterns matrix and a model of a convector. Besides the non-detection due to split up of some itemsets, Figure 8 (black ellipses) illustrates some non-detection due to the cut frequency chosen and the length of extracted itemsets. A post-processing might be added to the global chain, so that adjacent regions are interpolated. Thereafter, the spaceheating detector from the selected itemsets is deduced. The temporal support of non-detection intervals is around 10% of the space-heating support, without any post-processing. We also stress that the chosen method of extraction does not lead to any false alarms. The detector described above is based on the spectral properties of a convector load and should be applicable to other convection system since the corresponding patterns are frequent compared to the other devices ones.

CONCLUSIONS AND FUTURE WORKS

In this work, a new method of interpretation of a particular time-frequency representation, namely the spectrogram in the case study of electric load monitoring, is investigated. The method conjoins a time-frequency detector and a data mining method to extract some pertinent patterns. The T-F detection is based on a statistical model of the spectrogram. An existing detector of non-stationarities is implemented and tested on real data. A model of the convectors given the priors on the temporal and the spectral properties of the electric signature is used. The itemsets extraction method operates on the binarised time-frequency map and provides the frequent itemsets. Finally, the detection of the periodic signals from the unique mixture is performed by selecting the largest patterns of attributes. The method presented herein is dedicated to the detection and the extraction of periodic locally stationary signals embedded in a non-stationary observation. An application of this algorithm is the estimation of the spaceheating operating schedule given the unique measurement of the active power, without any intrusion in the house and any knowledge on the electric installation. The designed detector has no false alarm, which is an important criterion in the point of view of the customer. The detector would be evaluated on a larger data basis, where each device active power is recorded. If the performances of the proposed detector are satisfying, it would be evaluated in real life through tests on an industrial platform. Other spaceheating system would be thereafter studied.