Dr Jean-Michel André

Philippe Jonnard

X-ray spontaneous emission control by 1-dimensional photonic bandgap structure

Keywords: 42.25.Bs Wave propagation, transmission and absorption 42.70.Qs Photonic bandgap materials 78.20.-e Optical properties of bulk materials and thin films 78.67.Pt Multilayers, superlattices, photonic structures, metamaterials 78.70.En X-ray emission spectra and fluorescence

The possibility of controlling the x-ray spontaneous emission of atoms embedded in a 1dimensional photonic bandgap structure by the so-called Purcell effect, is studied.

Calculations of the spontaneously emitted power are presented from Fermi's golden rule in the framework of the Wigner-time approach extended to absorbing media. Numerical simulations are compared to experimental results for the case of the K emission from silicon atoms excited by electrons within a Mo/Si multilayer Bragg reflector. The inhibition or enhancement of x-ray emission from such structures appear to be feasible.

I. INTRODUCTION

Alteration of the spontaneous atomic emission rate by non-trivial boundary conditions on the radiation was pointed out in a classic paper by Purcell [1]. From Fermi's golden rule it appears that this rate is proportional to the density of mode (DOM) of the radiation in the medium. Model of this effect valid in the Wigner-Weisskopf regime, has been done in a classical formalism [2,3]. The advent of Photonic Band-Gap (PGB) structures has permitted to manipulate the DOM and then to control the spontaneous emission (SE) of atoms. Both enhancement and inhibition have been performed in the long wavelength domains by means of a one-dimensional PGB (1D-PGB) structure equivalent to a distributed Bragg reflector (DBR): it was possible to enhance the emission provided that the emission frequency of the atomic transition coincides with the frequencies of the edges of the band-gap and to stop it if the emission frequency of the atomic transition is located inside the band-gap [4,5].

Indeed the DOM of a DBR may present important peaks at the band-edges that allow the high extraction efficiency of the emission. The ability to control the spontaneous emission has profound consequence on many optoelectronics devices [6].

On another side, the progresses in the techniques of thin film deposition have permitted to fabricate periodic multilayered stack on nanometric scale that work as efficient Bragg reflector in the short wavelength domain of the electromagnetic spectrum (uv and x-rays) [7,8]. Then it appeared that these structures can be implemented to monitor the atomic emission in this spectral region. Thus Kossel-like structures have been observed under electron excitation in several Bragg mirrors for the Mα emission of tungsten in W/C multilayers [9,10] or the Kα emission of silicon in Mo/Si multilayers [11]. These works were, to our knowledge, the first experimental evidence of the possibility to enhance atomic SE in the x-ray domain by means of a 1D-PBG. The implementation of 1D-PBG as Distributed Bragg Reflector in the x-ray domain should play an important role in the development of xray laser [12] as emphasized by Yariv and Yeh. In this paper, the feasibility of controlling the SE in the x-ray domain by using Bragg multilayer structure is studied through the concept of DOM. First we give the formula deduced from Fermi's golden rule to calculate the SE rate from the DOM and the normal modes (NM). Then we explain how to compute the DOM in the Wigner time approach (WTA) for absorbing media and the NMs for a Bragg reflector. We discuss the feasibility of controlling the SE through numerical examples. Finally the results of the theoretical model are compared to experimental data obtained for the Si-K emission line emitted from a Mo/Si Bragg reflector under electron excitation.

II. SPONTANEOUS EMISSION RATE IN A 1D-PBG STRUCTURE

The spontaneous emission rate (SER) of an emitting atom is given by Fermi's golden rule :

! (R) = 2" ! 2 f H i 2 # ($ i % $ f ) (1)
in which i and f are the initial and final states of energies !! i / f . H is the interaction Hamiltonian of the whole system including emitting atoms, radiation and the optical « box » ;

it is given by :

H = ! ! µ(R). ! E(R) (2)
where ! µ(R) and ! E(R) are the dipole moment operator of the atom positioned at R and the electric field operator respectively. In the Wigner-Weisskopf regime, the mode density is the same classically or quantum electrodynamically ; then the SER can be derived from the classical model of a point dipole inside a medium with an inhomogeneous dielectric constant ε(r). In this approach the emitting atom is modeled by an harmonically oscillating dipole with frequency ω 0 and dipole moment ! µ located at R ; the corresponding current density J(r,t) to is taken :

J(r,t) = ! ! µ(R)cos(" 0 t) #(t) $ (r % R) (3)
The steady-state rate of power emission after all transients have vanished is given by :

P(t) = ! 2 2 " 2 µ 2 a k (R). ! µ # 2 $ (% 0 &% k ) d 3 k (4)
For sake of consistency the calculation are summarized in the Appendix. By integrating this expression over all wave-numbers k, the DOM !(" ) = dk d" appears as a result of changing variables k to ω k . It yields the power spectrum P ω of the emitted radiation as follows :

P ! = C "(! ) a k(! ) 2 ( 5
)
where C is some constant. Since The DOM is independent of position inside either structure, the averaging of the power spectrum is effected by replacing the values of a k(!)

2

by the corresponding averaged value :

a k(! ) (z) 2 = 1 Z max " Z min a k(! ) (z) 2 dz Z min Z max # ( 6
)
The integration is performed along the z-axis in the region occupied by the emitting medium between z = Z min and z = Z max .

III. DOM FOR ABSORBING MEDIA IN THE WIGNER-TIME APPROACH

III.1 Matrix transfer formalism for absorbing media

In the Wigner-time approach the DOM, ! ("), is considered to be the reciprocal of the group velocity v :

!(" ) = 1 v = dk(" ) d" ( 7
)
Since the complex transmission coefficient t(ω) contains phase information from which the dispersion relation k(ω) can be extracted, one understands easily that the group velocity v, hence the DOM ! (") can be deduced from t(ω) as explained in a paper [13] by Bendickson et al. These authors show that :

!(" ) = 1 D y' x # x ' y x 2 + y 2 (8)
where D is the total length of the structure, and x and y the real and imaginary part of t, respectively; the prime denotes differentiation with respect to ω. It is common to present a dimensionless DOM ρ(ω) obtained by multiplying ! (") by the group velocity in the bulk medium v bulk .

Our purpose now is to calculate the transmission coefficient of a N-period Bragg reflector made up with absorbing materials. As usual when dealing with the theory of electromagnetic propagation in linear stratified structure, we apply the transfer matrix formalism. Let us recall that the transfer matrix describes the transfer properties of the stratified system : it relates the amplitude of the incoming and outgoing electric fields and its expression depends on the geometry under consideration. In the Left-To-Right propagation (LTR) case considered in this paper and described in Figure 1, one has :

1 ! r ! " # $ % & = ! M ! t 0 ! " # $ % & (9)
with [14,15] :

! M = 1 ! t # ! r # ! t ! r ! t 1 # ! t ! " # # # # $
where ! t and ! r are the complex transmission and reflection coefficients of the whole system, respectively. The arrow → used as superscript expresses the LTR propagation.

This expression of the matrix transfer is obtained from the linearity of the problem and from the time-reversal symmetry that is the fact that Maxwell's equations are symmetric under time reversal. To describe the time reversal in our context of absorbing media, a special complex conjugation is required. We denote this operation by the superscript #. It means a standard complex conjugation plus the reversal of all the signs of the imaginary part of the complex optical indices. The requirement to this special conjugation can be justified as follows: in a time reversed process for propagation in absorbing media, the fields are actually amplified ; converting absorption to amplification requires the reversal of the sign of the imaginary part of the refractive index. This point is crucial in our calculation. Generally the time reversal for lossless systems is got by the mere complex conjugation [14]. The calculation of the DOM of absorbing composite quasi-periodic stratified structure requires the implementation of generalized Stokes reciprocity, where #-conjugation associated with the time reversal and parity reversal are involved. This generalization is presented in

1 0 Stratified medium Z ! t ! r
Ref. [15]. For convenience we recall that parity reversal is a process where a field of unit amplitude incident from the right gives rise to new Right-To-Left (RTL) complex transmission and reflection coefficients ! t and ! r , respectively. The RTL process is displayed in Figure 2.

III.2 DOM of a N-period Bragg reflector

The N-period Bragg reflector is made of a stack of N bilayers with two materials denoted by 1 and 2. The LTR transfer matrix of the N-period system can be expressed by means of the Chebychev polynomials of second kind U N (x) and of the coefficients C ij of the bilayer transfer matrix as follows [16] :

! M = C 11 U N !1 (Kd) ! U N ! 2 (Kd) C 12 U N !1 (Kd) C 21 U N !1 (Kd) C 22 U N !1 (Kd) ! U N ! 2 (Kd) " # $ % & ' (11a)
where

C 11 = 1 1 ! r 1,2 2 (exp(i") ! r 1,2 2 exp(!i#)) (11b)
C 22 = 1 1 ! r 1,2 2 (exp(!i") ! r 1,2 2 exp(i#)) (11c)
C 12 = r 1,2 1 ! r 1,2 2 (exp(!i") ! exp(i#)) (11d)
C 21 = r 1,2 1 ! r 1,2 2 (exp(i") ! exp(!i#)) (11e)
with

! = k "1 d 1 + k " 2 d 2 and ! = k "1 d 1 # k " 2 d 2 (11f)
The term k !j is the perpendicular component of the wave-vector in the layer j = 1,2. It obeys the Snell's law. The coefficient r 1,2 is the Fresnel coefficient associated with the reflection at the interface between the layer 1 and the layer 2 ; its expression is given by the Fresnel formula [16,[START_REF] Yeh | Optical Waves in Layered Media[END_REF]. K is the so-called Bloch-Floquet wave-number ; it corresponds to the wave-number associated with an infinite periodic system, i.e. where N → !. In fact the main problem in the determination of the transfer matrix comes from the calculation of K. Different ways can be used to calculate it [16,[START_REF] André | [END_REF]. It can be obtained from the trace S of the bilayer transfer matrix as :

K = 1 d cos !1 (S) = 1 d cos !1 (C 11 + C 22 ) (12)
The transmission coefficient ! t (!) and the reflection coefficient ! r (!) can be deduced by identification of Eqs. ( 10) and (11). Thus,

! t (! ) = (C 11 U N "1 (Kd) " U N " 2 (Kd)) "1 (13) and ! r (! ) = ! t (! ) C 21 U N "1 (Kd) (14)

III.3 Numerical for a Mo/Si Bragg reflector with N periods

The thickness d of the bi-layer is obtained from the Bragg condition :

d = c p! " r sin# r ( 15
)
where ! r is the Bragg frequency, ! r the Bragg angle and p the order of diffraction. We shall see that the value of the Bragg frequency and of the Bragg angle depends on the problem that one has to solve. If one wishes to obtain an inhibition of the emission corresponding to the resonant frequency ω 0 , ! r will be chosen to be equal to ω 0 so that the resonant frequency lies in the forbidden gap. On the contrary if one wishes to obtain an exaltation of the emission at resonant frequency ω 0 , ! r will be chosen so that the frequencies ω + and ω -of the band-edges coincide with ω 0 ; the choice is not obvious and can be done from a a priori estimation of the bandwidth Δω; a strategy to reach this goal will be presented hereafter.

As an example, we consider a molybdenum (Mo) and silicon (Si) Bragg reflector. This Mo/Si system is standard in the soft-x-ray domain as Bragg mirrors. We deal with the Silicon K (Si-K) emission, which means the Bragg reflector is designed to be able to diffract at the energy E 0 = 1740 eV, corresponding to the resonant frequency ω 0 = 2643 10 15 s -1 ; as explained previously one obtains the period d by applying the condition (15) with ! r = ω 0 .

One has d = 4.8 nm for a Bragg angle ! r = 4.26°. As expected, one observes peaks at the band-edges (Ω -, Ω + ) but the peak corresponding to the reduced frequency Ω -is more intense than the other one at the reduced frequency Ω + .

This effect corresponds to an asymmetry in the transmittance and arises from the occurrence of absorption in the materials. Figure 4 shows that quasi-symmetry is obtained in the noabsorption case in the case of a Bragg reflector made of 100 Mo/Si bilayers. It is observed that the curves with and without absorption are slightly shifted. This is probably due to the fact that the calculations are made without taking into account the variation of the indices as a function of the wavelength. When the number N increases, the peak at the band-edge Ω - becomes narrower and more intense. This result is in agreement with the well-known fact that for an infinite periodic system, the limit of the Brillouin zones in the (k,ω) domain have a step-like shape and consequently, the DOM !("

)

IV. CHOICE OF THE PERIOD FOR INHIBITION AND ENHANCEMENT OF THE SPONTANEOUS EMISSION IV.1 Inhibition

The goal is to inhibit the SE of frequency ω 0 (energy E 0 ) by locating this frequency between the band-edge ω -(or ω + if absorbing is very small) of the DOM. In this case, one may fix the Bragg frequency ! r that is located practically in the middle of the bandgap (i.e. at ω + Δ ω/2 = ω + -Δ ω/2) just at the resonant frequency ω 0 ; the choice of Bragg angle ! r is generally governed by practical considerations (geometry of the experiment, minimum practical value of the d-spacing, …). Once ! r and ! r are given, one obtains d from Eq. ( 15).

Let us note that in the soft-x-ray domain where the effect of refraction can be important, it is valuable to implement of Bragg law that takes into account this effect; for further details, see for instance [16].

IV.2 Enhancement

This case is more difficult to deal with than the previous one. Indeed to get enhancement it is necessary to locate the resonant frequency ω 0 of the SE at one of the bandedge frequencies ω + or ω -; the problem is that the value of these frequencies depends on the value of the d-spacing. Thus, it is useful to have an approximate value of the bandwidth Δω to deduce an approximate value of the band-edge frequencies from the relationship :

! ± = ! r ± "! 2 (16)
The two possible values d + and d -of the d-spacing are then obtained by combining Eqs. (15) and ( 16) :

d ± = cp! (" 0 ! #" 2 ) sin$ r (17)
In absence of absorption, analytical expressions are available to estimate Δω. From Eq. (A.8) of Ref. [16], one deduces that the following equation can be used to estimate the bandwidth :

!" # 4 $ ! r (" 0 ) sin c( p$ % ) " 0 (18)
where sinc denotes the cardinal sine function, sin(x)/x.

V. CALCULATION OF THE NORMAL MODES V.1 Generalities

Let us recall that the normal modes a k (r) can be defined from the harmonic solution :

A k (r,t) = a k (r)exp(!i" k t) (19)
of the homogeneous form (i.e. in absence of current) of the wave equation for the vector potential A(r,t) :

! " ! " A(r,t) + 1 c 2 #(r) !! A(r,t) = 0 ( 20
)
in the transverse jauge condition. In fact since our problem is one-dimensional (1D), the 1D NM a k (z) obeys the Helmholtz equation :

a k "(z) + ! k 2 c 2 "(z) a k (z) = 0 (21)
and must obey a completeness relationship of the form :

a k * (z ! ) a k (z 0 ) dk = " t (z # z 0 ) (22)
where δ t stands for the ε-transverse Dirac function. They must also fulfill the boundary conditions, that is, their continuity and the continuity of their derivative with respect to the space variable z at the boundary as a result of the continuity of the parallel components of the magnetic and electric fields.

V.2 Case of the N-period Bragg reflector

Let us now consider the two homogeneous regions denoted 1 and 2 of optical indices n 1 and n 2 respectively, forming the unit cell of the stratified medium, that is the bilayer in a Bragg reflector. In each of the region, the NM consists of a superposition of Left-and Rightgoing waves of wave-vector normal component k !j , j = 1 or 2 as shown in figures 1 and 2.

Hence, one writes in the first unit cell with index n = 0 :

a j (n = 0) (z) = A j exp(i k ! j z) + B j exp("i k ! j z) ( 23
)
By virtue of the Bloch-Floquet theorem, it follows that the form of the mode in the p-th unit cell (n = p) can be generated from the one in the first unit cell (n = 0) by :

a j (n= p) (z) = exp(i p K d) A j exp(i k ! j (z " p d)) + Bj exp("i k ! j (z " p d)) # $ % & , pd < z < ( p + 1)d ( 24
)
where K is the Bloch wave-number corresponding to the overall solution a K (z) for a infinite structure ( N ! "). This wave-number will be determined as shown hereafter.

The continuity condition of the mode and its first derivative denoted by prime at z = pd and z = d (p + γ) leads to:

a 2 (n = p) (z = pd) = a 1 (n = p +1) (z = pd) (25) a ' 2 (n = p) (z = pd) = a ' 1 (n = p +1) (z = pd) (26) a 1 (n = p +1) (z = d( p + ! )) = a 2 (n = p +1) (z = d( p + ! )) (27)
and

a ' 1 (n = p +1) (z = d( p + ! )) = a ' 2 (n = p +1) (z = d( p + ! )) (28)
By inserting the solutions Eq. ( 24) in Eqs. [START_REF] Jackson | Classical Electrodynamics[END_REF][START_REF] Tai | Dyadic Green Functions in Electromagnetic Theory[END_REF](27)(28), it follows that the four coefficients A j and Bj (j = 1,2) obey the following 4 X 4 matrix equation :

M(K, k 1 , k 2 , d,! ) A 1 B 1 A 2 B 2 " # $ $ $ $ % & ' ' ' ' = 0 0 0 0 " # $ $ $ $ % & ' ' ' ' (29)
where the subscript ! has been removed for convenience, with :

M(K, k 1 , k 2 , d,! ) = exp(i" 1 + d) exp(i" 1 # d) #1 #1 i k 1 exp(i" 1 + d) #i k 1 exp(i" 1 # d) #i k 2 i k 2 # exp(ik 1 ! d) # exp(#ik 1 ! d) exp(ik 2 ! d) exp(#ik 2 ! d) #ik 1 exp(ik 1 ! d) ik 1 exp(#ik 1 ! d) ik 2 exp(ik 2 ! d) #ik 2 exp(#ik 2 ! d) $ ' ( ) ) ) ) (30)
where ! j + /" = K ± kj .

The condition for obtaining a nontrivial solution of Eq. ( 30) is that :

det[ M(K, k 1 , k 2 , d,! )] = 0 (31)
that is,

cos(Kd) = ! (k 1 ! k 2 ) 2 4 k 1 k 2 cos d(k 1 (" ! 1) + k 2 " ) [ ] + (k 1 + k 2 ) 2 4 k 1 k 2 cos d(k 1 (1 ! " ) + k 2 " ) [ ] (32)
Eq. ( 32) is the dispersion relationship whose solution gives the value of the Bloch wavenumber K. To obtain fully the expression of the mode in the n-th cell, one must determine the four coefficients: A 1 , B 1 , A 2 and B 2 . This can be done by calculating the kernel of the submatrix M (K, k 1 , k 2 , d,! ) . One gets :

A 1 = 2 k 2 ! exp id " 1 ! ( ) k ! exp id # k + ( ) + k + exp id # k ! ( ) $ % & ' (33)
A 2 = exp(!i d " k ! ) [k ! exp id# 1 + ( ) + k + exp id (# 1 + + 2" k 2 ) $ % & ' ! 2 k 2 exp(i d " k + )] (34)
B 1 = exp(!i d " k 2 ) { k + exp id(# 1 + + 2" k 2 ) $ % & ' ! 2 k 1 exp i d (2K + " k 1 ) [ ] ! k ! exp id(# 1 ! + " (2 k 1 + k 2 )) $ % & ' } (35) B 2 = 2 k 1 exp i d (2 K +! k + ) [ ] + k " exp i d # 1 + ( ) " k + exp i d (# 1 " + 2! k 1 ) $ % & ' (36)
with

k ! = k 2 ! k 1 and k + = k 2 + k 1 (37)
These results are a generalization of the Kronig-Penney model [14].

V.3 Kronig-Penney model

By vanishing the ratio γ (γ → 0), one obtains the results corresponding the Kronig-Penney (KP) model; let us recall that in this model the thickness of the layer 1 tends to zero and the optical potential is given by a Dirac comb.

From Eqs. (33-37) one obtains the coefficients A KP and B KP for the KP system :

A 1 KP = 2 k 2 (1 ! exp(id" 1 ! )) (38)
B 1 KP = 2 k 2 (exp(id! 1 + ) " 1) (39)
A 2 KP = k + exp(id! 1 + ) " k " exp(id! 1 " ) " 2 k 1 exp(2idK ) (40)
B 2 KP = 2 k 1 exp 2 i d K ( ) + k ! exp i d " 1 + ( ) ! k + exp i d " 1 ! ( ) (41)
in agreement with the results of Ref. [14]. The KP case for an infinite PBG , named "Diraccomb superlattice" has been studied in details for both TE and TM polarization in [21] and [22].

VI. POWER SPECTRUM IN THE INHIBITION AND ENHANCEMENT CASES

We have calculated from Eq. ( 5) the power spectrum versus the reduced frequency Ω of the Si K emission coming from the Mo/Si multilayer system under study in §. III.3. The number of bilayers is 40. The power is normalized so that, far away from the band-gap, it corresponds to the free-space power, equal to unity. This normalized power is plotted in

VII. EXPERIMENT

The details of the experiment can be found in Ref. [11], so we briefly recall its main characteristics. The multilayer is made of 40 Mo/Si bi-layers, the thickness of the Mo and Si layers being 1.6 and 3.2 nm respectively. An electron gun is used to produce the Si 1s ionizations necessary for the generation of the Si Kα emission. The electron energy is 6 keV and in this case all the bi-layers contribute to the emission. The rotation axis of the sample is perpendicular to the plane defined by the directions of the electrons and photons. The sample is rotated around 4.26°, i.e. the Bragg angle of the multilayer at the Si Kα wavelength (0.713 nm), with a precision of 0.09°. This angle is between the sample surface and the detection direction and we call it detection angle. During an experiment, the spectrometer, whose acceptance angle is 0.23°, is fixed at the Si Kα wavelength whose intensity is monitored as a function of the detection angle. The scan is made about ±2° around the 4.26° detection angle. It has been verified [11] that a similar experiment performed on a silicon single crystal gives a monotonous behavior as a function of the detection angle and a lower intensity than in the case with the multilayer.

We show in Figure 6 the intensity of the Si Kα emission as a function of the detection angle. The angular scale is shifted so that the angle of the maximum is zero. The curve increases with the angle until its reaches a maximum. Then a dip occurs, about 0.6° wide, and afterwards the intensity increases again. This behavior is well reproduced by the calculation of the DOM and the emitted power but on a narrower angular scale : in this case the width of the dip a slightly less than 0.2° for the DOM and about 0.3° for the power. The experimental and calculated curves are scaled vertically in order that the intensities of their maximum and minimum coincide. In the power calculation, the in-depth distribution of the Si K ionizations within the multilayer, calculated with a Monte-Carlo code, is used to weight the distribution of the NM in Eq. ( 6). To take into account the instrumental broadening introduced by the finite angular acceptance of the spectrometer, we apply a 0.23° moving average on the calculated power. The discrepancy between the width of the experimental and calculated curves is not well understood, but could be ascribed to a bad account of the effect of the acceptance angle of the spectrometer.

VIII. DISCUSSION AND CONCLUSION

The x-ray SER of a 1D-PBG-embedded emitter can be enhanced provided that the emission energy coincide with the band edge energies of the DOM. On the contrary the spontaneous emission will be partially extinguished if the emission energy is within the forbidden band of the DOM. To obtain enhancement of the SE, one has to center the forbidden band on an appropriate energy chosen so that the emission energy coincide with band edge energies. To do it, it is necessary to choose the correct d-spacing or detection angle θ, as demonstrated by the presented θ-scan experiment on a Mo/Si multilayer.

Nevertheless one has to keep in mind that the 1D-PBG is not a total omnidirectional band gap, which means that a complete suppression or enhancement of the total SE is not possible in all directions with this kind of device. It has been shown that in a high contrast 1-D PGB, even in mid-gap emission, the SE can be enhanced parallel to the stack layers where the high-index layers act as resonant waveguides [23]. To achieve overall SE suppression, a full 3D-PBG is required. One has to emphasize that this task is not obvious and evidence of crystals exhibiting a full photonic band-gap has not been reported to date even in the visible light. A step towards an extension of the emission inhibition to further directions is the increase of the "dimensionality'" of the PBG. In the x-ray regime, recent development of the so-called lamellar multilayer gratings [24] could offer an opportunity to achieve this goal.

We plan to extend the present study to these gratings both theoretically and experimentally.

The control of the spontaneous emission can also be obtained in more sophisticated systems such as a Fabry-Pérot (FP) resonators with the emitters located within the FP spacer.

In this case one can possibly draw advantages from the resonant defect states located within the forbidden band of the DOM [14,[START_REF] André | [END_REF].

It is important to be aware that the PBG does not act as a filter that would yield only a relative enhancement; the PBG gives rise to an absolute enhancement that could be profitably implemented for low-threshold lasing. The inhibition of the spontaneous emission could be used to favor the stimulated emission in an x-ray lasing system.

Figure 1 :

 1 Figure 1: Left-to-Right (LTR) propagation scheme. Z is the direction of stratification of the medium.

Figure 2 :

 2 Figure 2: Right-To-Left (RTL) propagation scheme.

Figure 3

 3 Figure 3 shows the DOM ! (") and the transmittance T versus the reduced frequency Ω = ω/ω 0 , calculated in the WTA for the Bragg angle ! r = 4.26°, a duty cycle γ (ratio of the Mo layer thickness d 1 to the bilayer thickness d) of 0.33 and for different values of the number of bilayers N. The calculations are performed for the s-polarization case and we use

Figure 3 :

 3 Figure 3: Calculation taking into account the absorption of the DOM (a) and the transmittance (b) of a Mo/Si Bragg reflector versus the reduced frequency as a function of the number of bilayers: 40 (solid line), 100 (dotted line) and 300 (dashed line).

  photonic band-gap ΔE PBG is close to 102 eV. Considering that the natural width of the Si-K line is close to 0.5 eV[20], that is much narrower than ΔE PBG , the control of the SE by matching the frequency ω 0 of the SE to the DOM shape seems feasible. Practically the tuning can be obtained by varying the d-spacing or the Bragg angle.

Figure 4 :

 4 Figure 4: Calculation with (dotted line) and without (solid line) absorption of the DOM (a) and the transmittance (b) of a Mo/Si Bragg reflector with 100 bilayers versus the reduced frequency.

Figure 5 .Figure 5 :

 55 Figure 5. For inhibition, Ω = 1 have to be situated within the band-gap, which means that the propagation of the radiation at the frequency ω 0 of the SE is forbidden. For enhancement, one obtains quite similar figures but with an offset of the abscissa scale : due to the absorption the enhancement cannot be achieved at the upper band-edge (Ω + ) but only at the lower band-edge (Ω -). Assuming a Bragg angle equal to 4.26°, inhibition can be obtained for d = 4.8 nm, while for enhancement at Ω + , d = d + = 4.65 nm, and for enhancement at Ω -, d = d -= 4.95 nm. One notes that these different cases are achieved for values of d-spacing which are different

Figure 6 :

 6 Figure 6: Intensity of the Si Kα emission coming from a Mo/Si mulilayer excited by 6 keV electrons (solid line) as a function of the detection angle and comparison with the calculation of the DOM (dotted line) and the emitted power broadened by the angular acceptance of the spectrometer (dashed line).

& & & & (10)

& & & &

Appendix

The power radiated by the current J(r,t) contained in the volume V and associated with the radiating dipole can be written [START_REF] Jackson | Classical Electrodynamics[END_REF] :

Since the potential A(r,t) can be expressed via the dyadic propagator ! D(r,t; R,! ) which obeys the homogeneous wave Eq. ( 22), as follows [START_REF] Tai | Dyadic Green Functions in Electromagnetic Theory[END_REF] :

and since ! D (r,t;R, !) can be expanded in the NMs a(r) by :

where sinc stands for the cardinal sine, the total energy

is given by :

Assuming that the current density is given by Eq. ( 3), the total energy is in this case :

Eq. ( 4) is deduced from Eq. (A.4) by considering the steady-state rate of emission power :


Ana Paula Estrada-Vargas

Ernesto López-Mellado

Members IEEE Jean-Jacques Lesage
email: jacques.lesage@lurpa.ens-cachan.fr

J-J Lesage

An Identification Method for PLC-based Automated Discrete Event Systems

Keywords: Automated DES, Identification, Interpreted Petri Nets, Programmable Logic Controllers I

This paper deals with identification of automated Discrete Event Systems (DES). A method for processing sequences of input/output signals from PLC-based controlled DES is proposed; it yields an interpreted Petri net model describing the closed-loop behavior of the automated DES. This new method extends a previously presented approach by taking into account the technological characteristics of industrial controllers and the data collection requirements. Based on polynomial-time algorithms, the method is implemented as an efficient software tool; its use is illustrated through a case study dealing with an automated manufacturing system.

In recent years, the scientific community has proposed identification approaches (based on Petri net or automata) for obtaining approximated models of DES whose behavior is unknown or ill-known. In the context of automated DES, identification methods can be complementary to established modeling techniques; identification builds a closed-loop controller-plant model, which is generally obtained by a synchronous composition operation. Three main approaches for identifying DES have been proposed in literature [START_REF] Estrada-Vargas | López-Comparative Analysis of Recent I Discrete-Event Systems[END_REF].

The incremental synthesis approach, proposed in [START_REF] Meda | A passive meth discrete event systems[END_REF][3] deals with unknown partially measurable concurrent DES exhibiting cyclic behavior. Several algorithms have been proposed allowing the on-line building of interpreted Petri net (PN) models from the DES outputs. Although the techniques are efficient, the obtained models may represent more sequences than those observed.

In [START_REF] Giua | Identification integer programming[END_REF] a method to build a free labeled PN from a finite set of transitions strings is presented. This method is based on the resolution of an Integer Linear Programming (ILP) problem; the obtained PN generates exactly the observed language. This approach has been extended to other IPN classes [START_REF] Cabasino | Petri nets from their coverability gra Conference on Decision & Control[END_REF], [START_REF] Dotoli | Mangi discrete event systems using Petri nets[END_REF], [START_REF] Fanti | Fault d discrete event systems using Petri nets Workshop on Discrete Event Systems 435[END_REF]; however, issues regarding applications to actual DES have not yet been addressed in these works.

Another recent off-line method [START_REF] Klein | Fau Systems using an identification ap Congress[END_REF] allows building a non deterministic finite automaton (FA) from a set of input/output sequences, experimentally measured from the DES to be identified. Under several hypotheses, the constructed FA generates exactly the same input/output (I/O) sequences of given length than observed ones. The method was conceived for fault detection in a model-based approach [START_REF] Roth | An F Systems Based on an Identified Mode on Information Control Problems in M Moscow[END_REF]. Extensions to this work propose an identification method performing optimal partitioning of concurrent subsystems for distributed fault detection purposes [START_REF] Roth | Black event systems with optimal partitioni[END_REF].

In a previous paper we proposed a method for synthesizing interpreted PN (IPN) for coping with concurrent partially observable DES [START_REF] Estrada-Vargas | López-Me Identification of Concurrent Discret Cyclic Behavior[END_REF]; it processes a set of cyclic sequences of binary output signals yielding models including silent transitions and non labeled places.

In this paper that method has been extended for dealing with sequences of I/O signals captured during the closedloop operation of PLC-based controlled DES. Several technological characteristics of industrial controllers are taken into account in data collection and processing. The obtained model is a safe (1-bounded) IPN describing the controller-plant concurrent behavior, including those that are not measurable directly from the PLC.

The paper is organized as follows. In section II the background on Petri nets and languages is outlined. Based on these definitions the problem of DES identification is stated in section III in terms of language associated to an IPN. Several constraints inherent of real controlled DES are analyzed in section IV. The identification algorithm is given in section V. Finally, a case study is included in section VI.

II. BASICS ON PETRI NETS AND LANGUAGES

This section presents the basic concepts and notation of PN and IPN used in this paper. The symbol • t j (t j

•

) denotes the set of all places p i such that I(p i ,t j ) 0 (O(p i ,t j ) 0). Such places are called input (output) places of t j . Analogously, • p i (p i

•

) denotes the set of input (output) transitions of p i .

The incidence matrix of G is ܥ ൌ ܥ ା െ ܥ ି , where

ܥ ି ൌ ൣܿ ି ൧; ܿ ି ൌ ሺܫ ǡ ݐ ሻ; and ܥ ା ൌ ൣܿ ା ൧Ǣ ܿ ା ൌ ܱሺ ǡ ݐ ሻ are
the pre-incidence and post-incidence matrices respectively.

A marking function M : P→= represents the number of tokens residing inside each place; it is usually expressed as an |P|-entry vector. = + is the set of nonnegative integers. Definition 2: A Petri Net system or Petri Net (PN) is the pair N = (G,M 0 ), where G is a PN structure and M 0 is an initial marking.

In a PN system, a transition t j is enabled at marking M k if ∀p i ∈ P, M k (p i ) I(p i ,t j ); an enabled transition t j can be fired reaching a new marking M k+1 which can be computed as

M k+1 = M k + Cv k , where v k (i) = 0, ij, v k (j) =
(G,M 0 ). A PN is called safe if ܯ א ܴሺܩǡ ܯ ሻǡ  א ܲǡ ܯ ሺ ሻ ͳ.
Now it is defined IPN, an extension to PN that allows associating input and output signals to PN models.

Definition 3 : An IPN (Q, M 0 ) is a net structure Q = (G,Σ,Φ,λ,ϕ) with an initial marking M 0 . • G is a PN structure, Σ = {α 1 , α 2 , .
.., α r } is the input alphabet, and Φ = {φ 1 , φ 2 ,..., φ q } is the output alphabet.

• λ : T→ Σ∪{ε} is a labeling function of transitions, where ε represents a system internal event externally uncontrollable; it is not allowed that the symbol ε is associated to more than one t j ∈ p i • .

• ϕ : R(Q,M 0 )→(= + ) q is an output function, that associates to each marking in R(Q,M 0 ) a q-entry output vector;

q=|Φ| is the number of outputs. ϕ is represented by a q×|P| matrix, such that if the output symbol φ i is present (turned on) every time that M(p j ) 1, then ϕ (i,j) = 1, otherwise ϕ(i,j) = 0.

When an enabled transition t j is fired in a marking M k , then a new marking M k+1 is reached. This behavior is represented as

1 + ⎯→ ⎯ k j t k M M
; the state equation is completed with the marking projection y k = ϕ M k , where y k ∈ (= + ) q is the k-th output vector of the IPN.

According to functions λ and ϕ, transitions and places of an IPN (Q,M 0 ) can be classified as follows.

Definition 4: If λ(t i ) ε the transition t i is said to be controllable (t i can be fired when the associated input symbol is presented). Otherwise it is uncontrollable (t i is autonomously fired). A place p i ∈P is said to be measurable if the i-th column vector of ϕ is not null, i.e. ϕ(•,i) 0.

Otherwise it is non-measurable. P = P m ∪ P u where P m is the set of measurable places and P u is the set of non-measurable places.

Definition 5: The l-length I/O language of an IPN (Q,M 0 ) is defined as:

͉ ሺ ǡ ሻ ൌ ሼቂ ሺ௧ శభ ሻ ఝሺெ శభ ሻ ቃ ǡ ቂ ሺ௧ శమ ሻ ఝሺெ శమ ሻ ቃ ǡ ǥ ǡ ቂ ሺ௧ శ ሻ ఝሺெ శ ሻ ቃ ȁܯ ௧ శభ ሱۛሮ ܯ ାଵ ௧ శమ ሱۛሮ ܯ ାଶ ǥ ܯ ାିଵ ௧ శ ሱሮ ܯ ା Ǣ ܯ א ܴሺܳǡ ܯ ሻǢ ൌ Ͳǡͳǡ Ǥ Ǥ ሽ.

III. IDENTIFICATON PROBLEM STATEMENT

Consider the following language definitions used in the statement and solution of the identification problem.

Definition 6: The set of observed input/output (I/O) words of a DES S with m inputs and n outputs is

Ȟሺ ሻ ൌ ሼݓ ଵ ǡ ݓ ଶ ǡ ǥ ሽ , such that ݓ ൌ ሺቂ ூ ሺଵሻ ை ሺଵሻ ቃ ǡ ቂ ூ ሺଶሻ ை ሺଶሻ ቃ ǡ ǥ ቂ ூ ሺȁ௪ ȁሻ ை ሺȁ௪ ȁሻ ቃሻ , where ቂ ூ ሺሻ ை ሺሻ
ቃ is the j-th observed input/output vector of size m+n in sequence ݓ and ȁݓ ȁ is the length of the I/O word ݓ .

Definition 7:

The observed k-length I/O language of a DES S is defined as ࣦ ሺ ሻ ൌ ሼݓ ሺ݆ ͳሻݓ ሺ݆ ʹሻ ǥ ݓ ሺ݆ ݈ሻȁݓ א Ȟሺ ሻǡ ݆ ݈ ȁݓ ȁܽ݊݀ ݈ ݇ሽ.

Now the identification problem can be defined as follows: given a set of observed I/O words Ȟሺܵሻ generated by a real DES during its operation, the aim of the proposed method is to construct a safe IPN model (Q, M 0 ) such that ͉ ሺܳǡ ܯ ሻ ൌ ࣦ ሺܵሻ; such a language approximates the actual behavior of the DES.

In our approach, ࣦ ሺ ሻ is computed for inferring the states of the system according to the parameter ߢ, which is used to adjust the accuracy of the identified IPN, similarly as proposed in [START_REF] Klein | Fau Systems using an identification ap Congress[END_REF].

IV. IDENTIFYING REAL AUTOMATED DES

The systems considered in this work are closed loop controlled DES (Figure 1); they consist of a plant and its industrial controller (in many cases a Programmable Logic Controller: PLC). The behavior of such systems (i.e. the PLC-plant compound system behavior) can be observed by collecting the signals exchanged between controller and plant.

Figure 1. Closed loop controller-plant DES

Several phenomena, due to the interaction between plant and controller, increase the complexity of the identification process; however they must be taken into account when real controlled DES have to be identified: • An input evolution (signal emitted by the plant through a sensor) does not always provoke an output evolution (signal emitted by the PLC to an actuator). In practice, few of input changes provoke output evolutions; • Non simultaneous I/O events are often simultaneously observed; • When output changes are provoked by input changes, this causal relationship is not necessarily captured simultaneously; Now, we are going to explain these phenomena.

A. PLC treatment cycle

A PLC cyclically performs three main steps: "input reading" (I) where it reads the signals from the sensors, "program execution" (PEX) to determine the new outputs values for the actuators, and "output writing" (O) where the newly determined commands are sent to the plant actuators.

B. Experimental constraints

In this section, we use SFC only for precisely describe diverse behaviors. Recall that in the identification problem, the PLC program is assumed to be unknown.

Due to the PLC cycle, some situations between inputs and outputs could arise. Consider a situation described in These two signals will be simultaneously read during the next I phase and observed as simultaneous events in the identification data base. In DES theory events cannot occur simultaneously; so an observed event vector will therefore be defined as a change of value in an entry of an I/O vector. changes its value from 0 to 1, the correspondent change in "B" is not provoked immediately, since it is necessary first a change in output "A". In this case cause and effect cannot be captured simultaneously but will be detected only if we observe a sequence of 4 consecutive events. These three scenarios show that the implementation of a controller and its interaction with the plant introduces phenomena that must be taken into account by the identification algorithm.

Some approximations of the causality between inputs and outputs, allowing a more compact representation of the system's behavior can now be introduced. Consider the following I/O vector sequence σ involving one input x and two outputs A, B:

ݔ ܣ ܤ V ൌ Ͳ ͳ Ͳ ǡ ͳ ͳ Ͳ ǡ ͳ Ͳ ͳ
This sequence can be represented as: ܣ

V. IDENTIFICATION METHOD

A. General strategy

The identification procedure consists of several stages (Algorithm 1) that build systematically a safe IPN representing exactly the sampled output language of length ț+1 of the DES. The stages of this procedure are described below.

B. Sample processing B.1 Event sequences

As stated before, the data obtained from the system to be identified is a set of sequences of I/O vectors ݓ ଵ ǡ ݓ ଶ ǡ ǥ, such that ݓ ൌ ሺݓ ሺͳሻǡ ݓ ሺʹሻǡ ǥ ሻ, where ݓ ሺ݆ሻ refers to the j-th Definition 8: An observed event vector ɒ ୧ ሺሻ is the variation between two consecutive I/O vectors ୧ ሺሻ , ୧ ሺ ͳሻ; it is computed as ɒ ୧ ሺሻ ൌ ୧ ሺ ͳሻ െ ୧ ሺሻ. An input event vector Ⱦሺɒ ୧ ሺሻሻ is the variation between two consecutive input vectors ୧ ሺሻ, ୧ ሺ ͳሻ; it is computed as Ⱦሺɒ ୧ ሺሻሻ ൌ ୧ ሺ ͳሻ െ ୧ ሺሻ. A symbolic input event ɉԢሺɒ ୧ ሺሻሻ is a string representation of the input event vector Ⱦሺɒ ୧ ሺሻሻ; it is computed as:

ߣԢሺ߬ ሺ݆ሻሻ ൌ ෑ ቐ ܫ ̴ͳ ݂݅ ܫ ሺ݆ ͳሻ െ ܫ ሺ݆ሻ ൌ ͳ ܫ ̴Ͳ ݂݅ ܫ ሺ݆ ͳሻ െ ܫ ሺ݆ሻ ൌ െͳ ߝ ݂݅ ܫ ሺ݆ ͳሻ െ ܫ ሺ݆ሻ ൌ Ͳ ୀଵ
Then for every sequence ݓ , a sequence of observed event vectors ߬ ൌ ߬ ሺͳሻ߬ ሺʹሻ ǥ ߬ ሺȁ߬ ȁሻ is obtained. The maximum number of possible event vectors is ͵ ሺାሻ െ ͳ. However, in practice, only a small subset of them is observed.

ݓ ଵ ൌ ۉ ۈ ۈ ۈ ۈ ۇ ۏ ێ ێ ێ ێ ێ ێ ۍ Ͳ Ͳ Ͳ ͳ Ͳ Ͳ Ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ǡ ۏ ێ ێ ێ ێ ێ ێ ۍ ͳ Ͳ Ͳ Ͳ Ͳ Ͳ Ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ǡ ۏ ێ ێ ێ ێ ێ ێ ۍ ͳ Ͳ Ͳ Ͳ ͳ Ͳ Ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ی ۋ ۋ ۋ ۋ ۊ ǡ ݓ ଶ ൌ ۉ ۈ ۈ ۈ ۈ ۇ ۏ ێ ێ ێ ێ ێ ێ ۍ Ͳ Ͳ Ͳ ͳ Ͳ Ͳ Ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ǡ ۏ ێ ێ ێ ێ ێ ێ ۍ ͳ Ͳ Ͳ Ͳ Ͳ Ͳ Ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ǡ ۏ ێ ێ ێ ێ ێ ێ ۍ ͳ ͳ ͳ Ͳ Ͳ ͳ Ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ǡ ۏ ێ ێ ێ ێ ێ ێ ۍ ͳ Ͳ ͳ Ͳ Ͳ ͳ ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ǡ ۏ ێ ێ ێ ێ ێ ێ ۍ ͳ Ͳ Ͳ Ͳ Ͳ Ͳ ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ǡ ۏ ێ ێ ێ ێ ێ ێ ۍ Ͳ Ͳ Ͳ ͳ Ͳ Ͳ Ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ی ۋ ۋ ۋ ۋ ۊ ǡ ݓ ଷ ൌ ۉ ۈ ۈ ۈ ۈ ۇ ۏ ێ ێ ێ ێ ێ ێ ۍ Ͳ Ͳ Ͳ ͳ Ͳ Ͳ Ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ǡ ۏ ێ ێ ێ ێ ێ ێ ۍ ͳ Ͳ Ͳ Ͳ Ͳ Ͳ Ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ǡ ۏ ێ ێ ێ ێ ێ ێ ۍ ͳ ͳ ͳ Ͳ Ͳ ͳ Ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ǡ ۏ ێ ێ ێ ێ ێ ێ ۍ ͳ Ͳ Ͳ Ͳ Ͳ Ͳ Ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ǡ ۏ ێ ێ ێ ێ ێ ێ ۍ ͳ Ͳ Ͳ Ͳ Ͳ Ͳ ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ǡ ۏ ێ ێ ێ ێ ێ ێ ۍ Ͳ Ͳ Ͳ ͳ Ͳ Ͳ Ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ی ۋ ۋ ۋ ۋ

ۊ

According to definition 8, for every sequence ݓ , a sequence of observed event vectors ߬ ൌ ߬ ሺͳሻ߬ ሺʹሻ ǥ ߬ ሺȁ߬ ȁሻ is computed. During the process, if the difference has not been observed before, a new event vector ݁ is created and stored (߬ ሺ݆ሻ ൌ ݁ ). Now, for the example 1, the sequences τ i of the detected event vectors e j associated to I/O changes are obtained:

ݓ ଵ ൌ ۉ ۈ ۈ ۈ ۈ ۈ ۇ ۏ ێ ێ ێ ێ ێ ێ ۍ Ͳ Ͳ Ͳ ͳ Ͳ Ͳ Ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې భ ՜ ۏ ێ ێ ێ ێ ێ ێ ۍ ͳ Ͳ Ͳ െͳ Ͳ Ͳ Ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ۏ ێ ێ ێ ێ ێ ێ ۍ ͳ Ͳ Ͳ Ͳ Ͳ Ͳ Ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې మ ՜ ۏ ێ ێ ێ ێ ێ ێ ۍ Ͳ Ͳ Ͳ Ͳ ͳ Ͳ Ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ۏ ێ ێ ێ ێ ێ ێ ۍ ͳ Ͳ Ͳ Ͳ ͳ Ͳ Ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ی ۋ ۋ ۋ ۋ ۋ ۊ ǡ ߬ ଵ ൌ ݁ ଵ ݁ ଶ ߣ ᇱ ሺ݁ ଵ ሻ ൌ ̴ܽͳǡ ߣ ᇱ ሺ݁ ଶ ሻ ൌ ߝ ݓ ଶ ൌ ۉ ۈ ۈ ۈ ۈ ۈ ۇ ۏ ێ ێ ێ ێ ێ ێ ۍ Ͳ Ͳ Ͳ ͳ Ͳ Ͳ Ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې భ ՜ ۏ ێ ێ ێ ێ ێ ێ ۍ ͳ Ͳ Ͳ െͳ Ͳ Ͳ Ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ۏ ێ ێ ێ ێ ێ ێ ۍ ͳ Ͳ Ͳ Ͳ Ͳ Ͳ Ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې య ՜ ۏ ێ ێ ێ ێ ێ ێ ۍ Ͳ ͳ ͳ Ͳ Ͳ ͳ Ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ۏ ێ ێ ێ ێ ێ ێ ۍ ͳ ͳ ͳ Ͳ Ͳ ͳ Ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ర ՜ ۏ ێ ێ ێ ێ ێ ێ ۍ Ͳ െͳ Ͳ Ͳ Ͳ Ͳ ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ۏ ێ ێ ێ ێ ێ ێ ۍ ͳ Ͳ ͳ Ͳ Ͳ ͳ ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ఱ ՜ ۏ ێ ێ ێ ێ ێ ێ ۍ Ͳ Ͳ െͳ Ͳ Ͳ െͳ Ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ۏ ێ ێ ێ ێ ێ ێ ۍ ͳ Ͳ Ͳ Ͳ Ͳ Ͳ ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ల ՜ ۏ ێ ێ ێ ێ ێ ێ ۍ െͳ Ͳ Ͳ ͳ Ͳ Ͳ െͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ۏ ێ ێ ێ ێ ێ ێ ۍ Ͳ Ͳ Ͳ ͳ Ͳ Ͳ Ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ی ۋ ۋ ۋ ۋ ۋ ۊ ǡ ߬ ଶ ൌ ݁ ଵ ݁ ଷ ݁ ସ ݁ ହ ݁ ߣ ᇱ ሺ݁ ଷ ሻ ൌ ̴ܾͳ ̴ܿͳǡ ߣ ᇱ ሺ݁ ସ ሻ ൌ ̴ܾͲǡ ߣ ᇱ ሺ݁ ହ ሻ ൌ ̴ܿͲǡ ߣ ᇱ ሺ݁ ሻ ൌ ̴ܽͲ ݓ ଷ ൌ ۉ ۈ ۈ ۈ ۈ ۈ ۇ ۏ ێ ێ ێ ێ ێ ێ ۍ Ͳ Ͳ Ͳ ͳ Ͳ Ͳ Ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې భ ՜ ۏ ێ ێ ێ ێ ێ ێ ۍ ͳ Ͳ Ͳ െͳ Ͳ Ͳ Ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ۏ ێ ێ ێ ێ ێ ێ ۍ ͳ Ͳ Ͳ Ͳ Ͳ Ͳ Ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې య ՜ ۏ ێ ێ ێ ێ ێ ێ ۍ Ͳ ͳ ͳ Ͳ Ͳ ͳ Ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ۏ ێ ێ ێ ێ ێ ێ ۍ ͳ ͳ ͳ Ͳ Ͳ ͳ Ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ఱ ՜ ۏ ێ ێ ێ ێ ێ ێ ۍ Ͳ Ͳ െͳ Ͳ Ͳ െͳ Ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ۏ ێ ێ ێ ێ ێ ێ ۍ ͳ ͳ Ͳ Ͳ Ͳ Ͳ Ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ర ՜ ۏ ێ ێ ێ ێ ێ ێ ۍ Ͳ െͳ Ͳ Ͳ Ͳ Ͳ ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ۏ ێ ێ ێ ێ ێ ێ ۍ ͳ Ͳ Ͳ Ͳ Ͳ Ͳ ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ల ՜ ۏ ێ ێ ێ ێ ێ ێ ۍ െͳ Ͳ Ͳ ͳ Ͳ Ͳ െͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ۏ ێ ێ ێ ێ ێ ێ ۍ Ͳ Ͳ Ͳ ͳ Ͳ Ͳ Ͳ ے ۑ ۑ ۑ ۑ ۑ ۑ ې ی ۋ ۋ ۋ ۋ ۋ ۊ ǡ ߬ ଷ ൌ ݁ ଵ ݁ ଷ ݁ ହ ݁ ସ ݁

B.2 Sequences of ț-length event vector traces

From every sequence ߬ ൌ ߬ ሺͳሻ߬ ሺʹሻ ǥ ߬ ሺȁ߬ ȁሻ we compute sequences of event vector traces

߬ ൌ ߬ ሺͳሻǡ ߬ ሺʹሻǡ ǥ ߬
ሺȁ߬ ȁሻ such that every ߬ ሺ݆ሻ is the ߢlength substring of ߬ that finishes with ߬ ሺ݆ሻ. For the first ߢ െ ͳ elements of the trace sequence the event vector ߝ (zero vector) is used. Such traces are used to determine equivalent states, according to the following equivalence notion. Definition 9: Two states of the identified system are κequivalent if their I/O vectors are the same and if the ߢ last observed event vectors that lead to these states are the same. Following with the previous example, the sequences of traces using ߢ ൌ ʹ are:

߬ ଵ ଶ ൌ ߝ݁ ଵ ǡ ݁ ଵ ݁ ଶ for ߬ ଵ ߬ ଶ ଶ ൌ ߝ݁ ଵ ǡ ݁ ଵ ݁ ଷ ǡ ݁ ଷ ݁ ସ ǡ ݁ ସ ݁ ହ ǡ ݁ ହ ݁ for ߬ ଶ ߬ ଷ ଶ ൌ ߝ݁ ଵ ǡ ݁ ଵ ݁ ଷ ǡ ݁ ଷ ݁ ହ ǡ ݁ ହ ݁ ସ ǡ ݁ ସ ݁ for ߬ ଷ

C. Building the basic structure C.1 Representing event traces

Once the sequences of event vector traces have been obtained, every trace ߬ ሺ݆ሻ is related to a transition in the IPN through a function γ:T→ሼ߬ ሺ݆ሻሽ; the firing of a transition implies that ߢ consecutive event vectors related to such a transition have been observed. In order to preserve firing order between transitions, dependencies are created between them and associated with an observed marking through the function ߤǣ ܲ ௨ ՜ ሼMܯ ȁܯ א ܴሺܰǡ ܯ ሻሽ, which relates every non-measurable place with an observed marking, such that every transition has only one input place and one output place ( ݐ א ܶǡ ȁ ݐ x ȁ ൌ ȁݐ

x ȁ ൌ ͳ ). Notice that the number of nonmeasurable places is not predefined. When an event vector trace ߬ ሺ݆ሻ is found again in a ߬ the associated dependency must be used if it leads to the same observed marking.

Let e j be the last event vector in the trace ߬ ሺ݆ሻ ; the associated transition will be denoted as ݐ ೕ (more than one transition may have associated the same e j ). This strategy can be systematically performed following the next procedure which is an adaptation from that included in [START_REF] Estrada-Vargas | López-Me Identification of Concurrent Discret Cyclic Behavior[END_REF].

Algorithm 2. Building the basic IPN structure

Input: The set Ȟ ൌ ሼ߬ ሽ

Output: A PN structure G composed by p∈P u 1. ܶ ՚ Ǣ ܶܧ ՚ Ǣ ܲ ՚ ሼ ሽǢ ܯ ሺ ሻ ՚ ͳǢ ߤሺ ሻ ՚ ߬ ሺͳሻ; 2 ߬ א Ȟ 2.1 ݐ݊݁ݎݎݑܿ ՚  Ǣ 2.2 ߬ ሺ݆ሻ א ߬ ǡ ͳ ݆ ȁ߬ ȁ 2.2.1 If ߬ ሺ݆ሻ ב ܶܧ then ܶܧ ՚ ܶܧ  ሼ߬ ሺ݆ሻሽǢ ܶ ՚ ܶ  ൛ݐ ೕ ൟǢ J൫ݐ ೕ ൯ ՚ ߬ ሺ݆ሻǢ  א ܲ ൫ܫ ǡ ݐ ೕ ൯ ՚ ͲǢ ܱ൫ ǡ ݐ ೕ ൯ ՚ ͲǢ ǡݐ݊݁ݎݎݑܿ൫ܫ ݐ ೕ ൯ ՚ ͳ Ǣ If ݆ ൌ ȁ߬ ȁ and ሻݐ݊݁ݎݎݑܿߤሺ ൌ ߤሺ ሻ then ܱ൫ ǡ ݐ ೕ ൯ ՚ ͳǢ else ܲ ՚ ܲ  ሼ ௨௧ ሽǢ ݐ א ܶ ሺܫ ௨௧ ǡ ݐ ሻ ՚ ͲǢ ሺ ௨௧ ǡ ݐ ሻ ՚ ͲǢ ߤሺ ௨௧ ሻ ՚ ሻݐ݊݁ݎݎݑܿߤሺ ݁ Ǣ ܱ൫ ௨௧ ǡ ݐ ೕ ൯ ՚ ͳǢ ݐ݊݁ݎݎݑܿ ՚  ௨௧ Ǣ 2.2.2 If ߬ ሺ݆ሻ א ܶܧ then If  ȁ  ൌ xݐ ೕ ǡ ݐ ೕ א ܶǡ J൫ݐ ೕ ൯ ൌ ߬ ሺ݆ሻ ܽ݊݀ ߤሺ ሻ ൌ ሻݐ݊݁ݎݎݑܿߤሺ then ǡݐ݊݁ݎݎݑܿ݃݁ሺݎ݁݉  ሻǢ ݐ݊݁ݎݎݑܿ ՚ ሺݐ ೕ ሻ x Ǣ
else go to step 2.2.1, and take ߬ ሺ݆ሻ ב .ܶܧ If ݆ ൌ ȁ߬ ȁ and ሻݐ݊݁ݎݎݑܿߤሺ ൌ ߤሺ ሻ then ǡݐ݊݁ݎݎݑܿ݃݁ሺݎ݁݉  ሻǢ Since search operation has linear complexity and the algorithm implies the addition of a transition for every computed trace that has not been yet observed, then Algorithm 2 is executed in polynomial time on the number of observed output sequences and their maximum length .

Proposition 1.

The IPN G built through algorithm 2 represents all and only all the trace sequences (subsequences) in Ȟ .

The proof is similar to that included in [START_REF] Estrada-Vargas | López-Me Identification of Concurrent Discret Cyclic Behavior[END_REF].

Using the previous algorithm, the obtained PN corresponding to the three sequences of event vector traces of the example 1 is shown in Figure 6. Notice that one of the sequences is not cyclic.

ݐ ೕ ǡ ݐ ೕ א  x ȁܽ ് ܾ ݐ݃݁൫ݎ݁݉ ೕ ݐ ೕ ൯݉݁݁݃ݎሺݐ ೕ x ǡ ݐ ೕ x ሻ Rule 2: ݐ ೕ ǡ ݐ ೕ א  x ȁܽ ് ܾ ݐ݃݁൫ݎ݁݉ ೕ ݐ ೕ ൯݉݁݁݃ݎሺ ݐ ೕ x ǡ ݐ ೕ x ሻ Proposition 2. G' preserves the trace sequences in G.
The proof for is similar to that included in [START_REF] Estrada-Vargas | López-Me Identification of Concurrent Discret Cyclic Behavior[END_REF].

In the example, the application of the rules leads to the model of Figure 7(a). Since the initial place has two input transitions associated to e 6 , they can merge and also their input places.

C.3 Concurrent transitions

Other transformations may be performed when there are transitions that appear in the sequences in different order describing their interleaved firing; this behavior is exhibited by concurrent transitions. The analysis can be performed on a model component comprised between two transitions relied by several paths containing the concurrent transitions. If there are m! paths, we can explore if there exists m different transitions in the paths and every sequence is a permutation from each other. When it is verified, the subnet can be transformed into a concurrent component of G' preserving the same behavior.

In Figure 7(a), notice that between the transition associated to e 3 and the new transition associated to e 6 there are paths with all possible permutations of e 4 and e 5 ; then, we can transform this into a concurrent component and we obtain the net showed in Figure 7 Notice that this model preserves the same event vector sequences of the previous one.

The simplification by analysis of concurrency is not strictly necessary for representing the event vector sequences; however the equivalent model with concurrent transitions may be simpler. Although the analysis could be inefficient when the number of paths in the subnet is large, usually this number is reduced.

The aim of the simplification strategies given above is obtaining fairly reduced models useful for analyzing the DES behavior, rather than minimizing the number of nodes in the obtained models.

D. Adding interpretation to the PN model D.1Representing outputs changes

Once the event vector sequences are represented in the basic model, it must be completed by adding the output and input changes. Recall that event vectors are computed from the difference of consecutive I/O vectors; thus an e j relates measurable places representing the outputs changes. This procedure is detailed in [START_REF] Estrada-Vargas | López-Me Identification of Concurrent Discret Cyclic Behavior[END_REF]. The number of observable places is then equal to the number of outputs. The net with measurable places for the example is showed in Figure 7(c).

D.2 Model simplifying

Implicit non-measurable places can be removed; if there is a non-measurable place p k whose input and output transitions are exactly the same than any measurable place, then delete p k and its input and output arcs.

D.3 Representing input changes

Once the output adding and implicit places deleting has been performed, it only remains to add input information to complete the IPN model. Input information is associated with labels for transitions in a natural way given by the symbolic event input function. Algorithm 4 describes a systematic way to do it. The final model for the illustrative example is showed in Figure 8. The associated inputs for transitions are given by: λ'(e 1 )=a_1, λ'(e 2 )= ε, λ'(e 3 )=b_1 c_1, λ'(e 4 )=b_0, λ'(e 5 )=c_0, λ'(e 6 )=a_0.

Figure 8. Simplified IPN model

Since every one of the transitions in the net actually represents a sequence of event vectors of length κ, the I/O language of length κ + 1 of the net is equal to the observed I/O language. Even, for the above example, the I/O language of the IPN is equal to the observed I/O language, i.e. only the observed cyclic I/O sequences are represented by the evolution of the net.

Proposition 3. For a DES S and an identification parameter κ, Algorithm 1 yields an IPN model (Q,M 0 ) which represents exactly ࣦ ାଵ ሺܵሻ.

Proof. Since the deletion of implicit places does not change ͉ ሺܳǡ ܯ ሻ , the proof is made with the model obtained before this procedure. The firing of a transition ݐ in the system is not affected by the addition of arcs to and from ,ݐ since these arcs have been computed from differences of vectors in Ȟሺܵሻ. Then, also in this model, every event vector sequence ߪ of length less or equal than ߢ belongs to the language of the net iff it was observed.

The sequences of transitions of length less or equal than ߢ that can be fired, lead to markings in the measurable places that also belong to Ȟሺܵሻ (since the marking change provoked in the measurable places was obtained from the difference of observed vectors). Then, we have that sequences of observed output vectors of length less or equal than ߢ ͳ correspond to sequences of marking vectors in the net and ͉ ାଵ ሺܳሻ ൌ ࣦ ାଵ ሺܵሻ.ז

VI. CASE STUDY

Based on the algorithms presented in this paper, a software tool has been developed to automate the IPN model synthesis. Many examples on diverse size and complexity have been tested. Below we present a small size case study from [START_REF] Roth | An F Systems Based on an Identified Mode on Information Control Problems in M Moscow[END_REF] whose layout is showed in Figure 9. The function of this system is to sort parcels according to their size. It has 9 inputs (signal sensors from the system): a 0 ,a 1 , a 2 ,b 0 ,b 1 ,c 0 ,c 1 ,k 1 ,k 2 , and 4 outputs (signals to the actuators): A+,A-,B,C. After the identification process and before concurrence transformations, the model shown in Figure 10(b) is obtained. After some transformations considering phenomena described in section IV, the net in Figure 11 is obtained. These models have been automatically obtained using the developed software tool.

From the obtained figure, we can infer the following behavior:

-When a large parcel arrives (k1 and k2 rise), cylinder A is extended (A+ rise). After cylinder A pushes the parcel to the conveyor 3 (a0 falls, k1 falls, k2 falls, a1 rises, a1 falls and a2 rises), cylinder C extends (C rises) and cylinder A retracts to its initial position (A-falls). When parcel is completely pushed, cylinders retract to their initial position.

-When a short parcel arrives (only k1 rise), cylinder A is extended (A+ rise) until conveyor 2 is reached (a0 falls, k1 falls and a1 rises). Then, cylinder B extends (B rises) and cylinder A retracts (A-falls). When parcel is completely pushed, cylinder B retracts until its initial position.

This interpretation is valid only for the part of the behavior we have observed. If we consider that data collection has been made for a long time, we can state that the model is almost exact.

Definition 1 :

 1 An ordinary Petri Net structure G is a bipartite digraph represented by the 4-tuple G = (P,T,I,O) where: P = {p 1 ,p 2 ,...,p |P| } and T = {t 1 ,t 2 ,...,t |T| } are finite sets of vertices named places and transitions respectively; I(O) : P × T → {0, 1} is a function representing the arcs going from places to transitions (from transitions to places).

Figure 2 .

 2 Figure 2. PLC cycle and data collection At the end of the PEX phase the current values of inputs and outputs (I/O) are sent from the PLC to a computer and stored for a later treatment by the identification algorithm.

Figure 3 (

 3 current active step is #10; a and b are two input signals to the PLC; A and B are two output signals).

Figure 3 .

 3 Figure 3. A single input is the condition for state evolution Changes in the state and outputs will occur when signal b is active; however other input signals may evolve without consequence in the outputs. This must be considered in the identification algorithm. Consider now the time diagram in Figure 4. Two signals are asynchronously emitted by sensors of the plant between two successive "input reading" phases (I) of the PLC cycle.These two signals will be simultaneously read during the next I phase and observed as simultaneous events in the identification data base. In DES theory events cannot occur simultaneously; so an observed event vector will therefore be defined as a change of value in an entry of an I/O vector.

Figure 4 .

 4 Figure 4. Apparent simultaneous evolution of several inputs

Figure 5 .

 5 Figure 5. I/O causality and sequences of events

  ܤ Since such a behavior can be a consequence of the PLC cycle, it can be compacted as: ܣ ௫̴ଵ ሱሮ ܤ . Analogously, we can generalize:

Algorithm 1 .

 1 Global identification procedure Inputs: A DES and an accuracy parameter κ Output: (Q,M 0 ): an IPN model 1. Obtain the cyclic sequences w i of observed I/O vectors. 2. Compute event vector sequences τ i and symbolic input events λ' from the observed vectors. 3. For every sequence of event vectors τ i , create event vector traces τ i κ of length κ. 4. Create the non-observable behavior of the IPN and simplify it. 5. Complete the IPN adding the observed behavior and deleting implicit places.

  sequence ݓ . Such sequences may have different length. From these sequences, strings of observed event vectors are first computed.

Example 1 .

 1 Consider a DES with n = 4 output signals, Φ = {A, B, C, D}, and m = 3 input signals Σ = {a, b, c}. Three I/O sequences have been observed; vector entries correspond to distribution [a, b, c | A, B, C, D ] T

Figure 6 .Algorithm 3 .

 63 Figure 6. Basic model with sequences of event vector traces C.2 Simplifying the basic structure Additional node merging operations can be performed on the basic structure in order to obtain an equivalent trace model. Now we can take into account the event vector e j associated to transitions. Consider the following transformation rules. Algorithm 3. Simplifying the basic structure. Input: G Output: G': an equivalent IPN Apply the following rule on the initial place and iteratively on the merged places Rule 1: ݐ

  (b).

Figure 7 .

 7 Figure 7. (a) Model after merging (b) Simplified basic model (c) IPN model including measurable places

Algorithm 4 .

 4 Representing input changes Input: G', ߣԢሺ݁ ሻ Output; Q : the final model of the identification process Step 1. ݐ ೕ א ܶǡ ߣሺݐ ೕ ሻ ՚ ߣԢሺ݁ ሻ //Associate to each ݐ ೕ the symbolic input event registered at the detection of e j .

Figure 9 .

 9 Figure 9. Manufacturing system layout Due to the lack of room, we only consider nine observed cycles of production (a cycle consist in the arrival of a parcel and its sorting) given below by the nine I/O observed sequences of different length. The vector entries correspond to the following distribution: [A+ A-B C k1 k2 a0 a1 a2 b0 b1 c0 c1].After the identification process and before concurrence transformations, the model shown in Figure10(b) is obtained. After some transformations considering phenomena described in section IV, the net in Figure11is obtained. These models have been automatically obtained using the developed software tool.From the obtained figure, we can infer the following behavior:-When a large parcel arrives (k1 and k2 rise), cylinder A is extended (A+ rise). After cylinder A pushes the parcel to the conveyor 3 (a0 falls, k1 falls, k2 falls, a1 rises, a1 falls and a2 rises), cylinder C extends (C rises) and cylinder A retracts to its initial position (A-falls). When parcel is completely pushed, cylinders retract to their initial position.-When a short parcel arrives (only k1 rise), cylinder A is extended (A+ rise) until conveyor 2 is reached (a0 falls, k1 falls and a1 rises). Then, cylinder B extends (B rises) and cylinder A retracts (A-falls). When parcel is completely pushed, cylinder B retracts until its initial position.This interpretation is valid only for the part of the behavior we have observed. If we consider that data collection has been made for a long time, we can state that the model is almost exact.

Figure 10 .Figure 11 .

 1011 Figure 10. (a) Cyclic sequences (b) identified study VII. CONCLUSIONS Identification of automated concurrent Systems (DES) has been addressed. A pr method has been extended and adapte identification of closed-loop PLC-based co from vector sequences including both input order to cope with technological issue automation based on PLC, simplifying trans obtained model have been proposed, a compact representation of the system behavi Current research deals with the infe sequences from a single sequence and co priori known relations among DES identification.

  1, this equation is called the PN state equation. The reachability set of a PN is the set of all possible reachable marking from M 0 firing only enabled transitions; this set is denoted by R