Yue Wang

Fabio Pistolesi

Sensitivity of the mixing current technique to detect nano-mechanical motion

Keywords: numbers: 85.85.+j, 73.23.-b, 73.23.Hk

Detection of nano-mechanical displacement by transport techniques has reached high level of sensitivity and versatility. In order to detect the amplitude of oscillation of nano-mechanical oscillator a widely used technique consists to couple this motion capacitively to a single-electron transistor and to detect the high-frequency modulation of the current through the non-linear mixing with an electric signal at a slighltly detuned frequency. The method known as current-mixing technique is employed in particular for the detection of suspended carbon nanotubes. In this paper we study theoretically the limiting conditions on the sensitivity of this method. The sensitivity is increased by increasing the response function to the signal, but also by reducing the noise. For these reasons we study systematically the response function, the effect of current-and displacement-fluctuations, and finally the case where the tunnelling rate of the electrons are of the same order or larger of the resonating frequency. We find thus upper bounds to the sensitivity of the detection technique.

I. INTRODUCTION

Nano-electromechanical systems have great potentials as ultra-sensitive detectors for several physical quantities. Recent advances allowed to reach record sensitivity in mass sensing. [1][2][3] This has been possible by the detection of the frequency shift of ultralight oscillators when an additional mass is attached to it. Other exemples concern the detection of the tiny magnetic field generated by nuclear spins. This can be done by the opto-mechanical detection of the force generated by the magnetic dipoles, 4 but also with electro-mechanical means, 5 or by coupling to two-level systes. [6][7][8] The force sensitivity of the device is then the limiting factor for the sensitivity, and again recent advances showed that it is possible to obtain record force sensing with carbon-nanotube oscillators. 9,10 At the same time nano-mechanical oscillators can be so small that interaction between electronic and mechanical degrees of freedom may lead to new and unexpected phenomena [11][12][13][14][15] like the blockade of the current [16][17][18][19] , cooling [20][21][22] or unusual mechanical response. 23,24 In order to exploit nanomechanical resonators, or to study their properties, detection of mechanical motion is crucial. Most detection methods exploiting electronic transport are based on the high sensitivity of singleelectron transistors (SET) to a variation of the gate charge. By coupling capacitively the oscillator to the gate of the SET it is possible to detect the motion of the oscillator with a high accuracy. 25 The method has been used also to cool the oscillator by the back-action of the electronic transport. 26 The main difficulty of the method stems from the high frequency character of the oscillator motion that is typically in the 100 MHz-1 GHz range. Due to the high impedance of the SET, it is more convenient to down-convert the signal to lower frequency before extracting it. This can be achieved by non-linear mixing the mechanically generated modulation with a second high-frequency signal injected between source and drain. The signal at the difference of the two frequencies can be extracted and measured. To our knowledge, for nanomechanical resonators this method was implemented in metallic SET by the group of A. Cleland back in 2003. 27 It was later adapted to the detection of carbon nanotube by the group of P.L. McEuen. 28 It then became the method of choice for carbon nanotubes, leading to several breakthroughs: the observation of the first singleelectron backaction effects in carbon nanotubes, 29,30 ultrasensitive mass detection, 2,3 the detection of the charge response function in quantum dot, 31 the detection of magnetic molecules 32,33 and the observation of decoherence of mechanical motion. 34 The same method can also be implemented by frequency modulation. 35 It is clear that the technique is powerful and that it will continue to be used both for fundamental research and for applications. The question we want to address in this paper is which is the ultimate resolution that can be reached with this kind of detection. In order to do this we investigated three main issues. The first one is how to optimize the response function, that is the quantity ∂I mx /∂x m , where I mx is the measured signal, the mixing current, and x m the amplitude of the mechanical oscillation. The second one is to study the effect of current and mechanical fluctuations. These contribute to the fluctuation of the measured signal and in the end are at the origin of the signal to noise ratio. The third is to consider the case of a mechanical oscillator with a resonating frequency ω m faster than the typical tunneling rate of the electrons Γ. We will develop a theory of transport to obtain the mixing current for any ratio ω m /Γ. The case of a metallic and single-electronic level SET will be considered in details and explicit expressions will be given.

The paper is structured as follows: Section II gives an introduction to the mixing technique and provides the expression of the detector gain λ. Section III analyze current and mechanical fluctuations giving general ex- pressions. Section IV provides a general theory for the detector gain when the resonator frequency is comparable or larger than the typical tunnelling time. Section V and VI gives the explicit expressions for the response function and for the noise in the case of a metallic and a single electron SET. Finally Section VI gives our conclusions.

II. MIXING TECHNIQUE AND RESPONSE FUNCTION

Let us begin by describing the typical system used to measure the oscillation amplitude of a mechanical oscillator by detection of the mixing current. [27][28][29][30] As shown in Fig. 1 a conducting oscillator is capacitively coupled to the central island of a single-electron transistor: its displacement modulates thus the the gate capacitance C g (x), where x(t) is the displacement of the oscillator. We assume the presence of a single mechanical mode whose displacement is parametrized by x, a generalized coordinate with the dimensions of a length. We will consider that the SET is operated in the incoherent transport regime valid for Γ k B T , where Γ is the electron tunnelling rate and T the temperauture ( and k B are the reduced Planck constant and the Boltzmann constant, respectively). This is the standard case for nanomechanical devices. The current I through the device can be obtained by using the Master equation and in general it can be expressed as a function of the source-drain bias voltage V and on the gate charge n g = C g (x)V g /e, where V g is the gate voltage (see Appendix A for a short derivation). The current reads thus:

I = I(V, n g ) .
(1)

In this section we want to obtain the current response of the system when both V and V g are modulated at two slightly different frequencies ω 1 and ω 2 , both much smaller than Γ. We write

V g (t) = V g0 + V g1 (t) , V (t) = V 0 + V 1 (t) , (2)
where

V g1 (t) = V g1 cos(ω 1 t) and V 1 (t) = V 1 cos(ω 2 t).
Choosing ω 1 close to the mechanical resonating frequency ω m allows to drive the resonator, since the modulation of the gate voltage modulates the charge on the suspended part and thus induces an oscillating force [see also Eq. ( 18) in the following]. For small driving amplitude the oscillator responds linearly to the external drive:

x(t) = x m cos(ω 1 t + φ), (3)
where we always measure x from its equilibrium position. (Note that in general x m and φ depend on the driving frequency ω 1 .) The modulation of V g induces thus the following modulation of n g at linear order in the driving:

n g (t) = n g0 + C g V g1 e cos(ω 1 t) + C g x m V g0 e cos(ω 1 t + φ), (4
) where C g ≡ dC g /dx. It is convenient to introduce a length scale by defining L = C g /C g . From geometric considerations L has to be of the order of the distance of the gate from the oscillator, thus typically undreds of nm. The fluctuting part of n g can then be written as

n g1 (t) = n g0 V g1 V g0 cos(ω 1 t) + x m L cos(ω 1 t + φ) , (5)
where n g0 = C g V g0 /e. The mechanical term (x m /L) has a strong frequency dependence close to the mechanical resonance, and can thus be distinguished by the background electrostatic term (V g1 /V g0 ). The two contributions to the modulation of the gate charge can be combined in a single cosine term:

n g1 (t) = n g1 cos(ω 1 t + ϕ) . (6)
Assuming now that the oscillator frequency ω m , and thus also ω 1 and ω 2 , are much smaller than the typical tunneling rate Γ, one can use Eq. (1) to obtain the time dependent current in presence of time-dependent V and n g . For small modulation amplitude we Taylor expand Eq. (1) to second order in V 1 and n g1 obtaining

I(t) = I(V 0 , V g0 ) + ∂I ∂V V 1 (t) + ∂I ∂n g n g1 (t) + 1 2 ∂ 2 I ∂V 2 V 2 1 (t) + ∂ 2 I ∂V ∂n g V 1 (t)n g1 (t) + 1 2 ∂ 2 I ∂n g n 2 g1 (t) + . . . (7)
Only the term proportional to ∂ 2 I/∂V ∂n g has a component that oscillates at the frequency ω ∆ = ω 1 -ω 2 . This signal can be extracted by a standard lock-in technique that essentially allows to measure the quantity I mx :

I c mx = Tm 0 dt T m I(t) cos[ω ∆ t] . (8)
The other quadrature I s mx with sin(ω ∆ t) is defined in a similar way. Averaging over a long measurement time T m 1/ω ∆ one obtains:

I c mx = V 1 4 ∂ 2 I ∂V ∂n g C g V g1 + C g V g0 x m cos φ , (9)
I s mx = - V 1 4 ∂ 2 I ∂V ∂n g C g V g0 x m sin φ . (10)
The detector gain with respect to the two quadrature of x m is thus:

λ = 1 4e ∂ 2 I ∂V ∂n g C g V g0 V 1 . (11)
It measures the sensitivity of the mixing current signal with respect to the two quadratures of x. This quantity depends on the particular bias conditions of the SET, and will be studied in some details in Section V and VI for two explicit models. Note also that in order to obtain λ we need only the static expression for the current. This assumes that the electronic mechanism is much faster than the time dependence of the driving. In order to describe the case of a fast oscillator (to be discussed in Section IV we will need a detailed description of the charge dynamics, and the response function will be no more expressed only in terms of derivatives of the static non-linear current voltage characteristics.

III. EFFECT OF CURRENT AND DISPLACEMENT FLUCTUATIONS

Expression (11) assumes a deterministic evolution of both the current and the displacement of the oscillator x(t). In practice both quantity fluctuate, the first due to shot or thermal noise, and the second due to stochastic fluctuations induced either by the bias voltage or by the thermal fluctuations. In general one can then write the value of I mx in a specific time region as follows:

(I c mx ) n = (n+1)Tm nTm [I(t) + δI(t)] cos(ω ∆ t)dt , (12)
(we write the expression for I c mx , the one for I s mx is similar) where δI(t) and I(t) are the stochastic and deterministic (in phase with the external drive) part, respectively. We can define the time dependent mixing current as

I c mx (t) = (I c mx ) [t/Tm]
, where [α] stands here for the integer part of α. In terms of that the spectral density of the fluctuation of I c mx reads:

S mx (ω) = +∞ -∞ dte iωt I c mx (t)I c mx (0) -I c mx 2 . ( 13
)
We assume that the measuring time is much longer than any correlation time of the quantity δI(t). Different sections of the measurement time are thus uncorrelated and we can write:

S mx (ω) = Tm 0 dte iωt Tm 0 dt 1 T m Tm 0 dt 2 T m cos(ω ∆ t 1 ) cos(ω ∆ t 2 ) δI(t 1 )δI(t 2 ) . (14)
Defining S I (ω = 0) = 2 +∞ ∞ dt δI(t)δI(0) (the numerical factor 2 is conventional for the current-noise spectrum) we have

S mx (ω) = 1 4 S I (ω = 0) (e iωTm -1) iωT m ≈ 1 4 S I (ω = 0) . ( 15
)
Thus the mixing-current low-frequency noise is given simply by the low-frequency current noise spectrum S II . The factor of 4 comes from a different definition of the correlation functions and from the fact that we are collecting a single quadrature. The current noise can have different sources, we consider in the following the two main ones.

A. Shot-noise and thermal current fluctuations

The current fluctuates due to the discrete nature of the the charge. This is characterized by the current-spectral function (for time-independent bias and gate voltages):

S shot I (ω) = 2 dte iωt δI(t)δI(0) , (16)
where δI(t) = I(t) -I . For the case of a SET the current spectral function is well known. 36 As shown there it has a frequency dependent part at low frequency on the scale of the typical tunneling rate Γ. This implies that the correlation function is short ranged with respect to the measuring time T m . Actually it is typically even short ranged with respect to the time dependence of x and of the V or V g potentials. Its value can thus be obtained adiabatically, by assuming these parameters to be static. We only need its low frequency part that can, in general, be expressed in terms of the Fano factor F and the current I:

S shot I (ω = 0) = 2FIe (17)
where F depends on the details of the SET. In the tunnelling limit of uncorrelated tunneling F = 1, in most other cases the Fano factor is typically of the order of 1.

B. Displacement fluctuations

The electrons that cross the structure modify the charge on the gate that in turn modifies the force acting on the oscillator. This stochastic force, that has the same origin of the current-shot noise, induces fluctuations of the displacement, that changes in a much slower way, since the oscillator responds to an external force on the time scale given by its damping coefficient γ. 13,18,37,38 In order to keep the assumption that different averages over the measuring times are uncorrelated one needs T m γ 1. In principle, for very high-Q resonators the approximation should be reconsidered.

Let's begin by considering the force acting on the oscillator as a consequence of a variation of the charge on the gate. A recall of the basic expressions for the electrostatic energy is given in the Appendix A and Fig. 3 there shows the electrical scheme. The force acting on the oscillator is given by the derivative of the electrostatic energy performed at constant charge:

F = -Q 2 g ∂ ∂ x 1 2C g (x) = Q 2 g C g 2C 2 g , (18)
where Q g is the charge on the gate voltage (Fig. 3). The fluctuation of the force δF (t) due to fluctuation of Q g reads thus:

δF (t) = Q g C g C 2 g δQ g (t) . (19)
In general the variation of the charge on the gate is proportional to the variation of the charge on the central island of the SET. By an elementary electrostatic calculation (see Appendix A) δQ g /e = (C g /C Σ )δn, where

C Σ = C g + C L + C R
is the sum of the capacitances of the central island to all the electrodes and -ne is the total charge on the island. In conclusion one finds that

δF (t) = F 0 δn(t) (20)
with

F 0 = Q g eC g C g C Σ = 2 Q g e E C L (21)
the force acting on the oscillator when an electron is added to the dot and with E C = e 2 /(2C Σ ) the Coulomb energy of the SET. Note that F 0 is a crucial parameter, since it constitutes the electro-mechanical coupling constant. 18,23 One can estimate the typical value of F 0 :

Q g /e = 10-100, E C = 1 K, L = 100 nm, thus F 0 ≈ 10 -11 -10 -12 N.
The correlation function of the stochastic force acting on the resonator [S F (t) = δF (t)δF (0) ] is thus simply proportional to the correlation function of the charge on the island [S δn (t) = δn(t)δn(0) ]:

S F (t) = F 2 0 S n (t) , (22)
that can be calculated by the standard method of the master equation. For the case of a metallic dot see for instance Refs. 18 and 39. Its Fourier transform has a Lorentzian form with a width on the scale of Γ. Thus this force act as a white noise on the slow oscillator.

Let's now turn to the displacement correlation function. In order to evaluate it we use a simple Langevin approach. 11,16 We neglect the driving, since we are interested in the low frequency response. The Langevin equation reads

mẍ + mγ ẋ + kx = δF (t), ( 23
)
where m is the (effective mass) of the oscillator mode considered, γ the damping coefficient, and k the effective spring constant. The stochastic force generated by the electrons is also at the origin of the damping coefficient.

In general other effects participate, but close to the degeneracy point of the SET, when the current is maximal, the electronic contribution to the damping can dominate, as observed experimentally in Ref. 32. We will assume thus that γ is due only to the electronic damping. In equilibrium the fluctuation-dissipation theorem gives

S δF (ω = 0) = 2γmk B T. (24)
For finite eV k B T , the system is out of equilibrium and one has to evaluate explicitly γ and S δF from a direct calculation of S F (ω). As shown in Ref. 40 2m γ = dS F (ω)/dω| ω=0 . One can then always define an effective temperature by the relation S F (ω = 0) = 2γmk B T eff , since the oscillator has a very sharp response in frequency and the correlation functions are flat on that scale, one can always interpret the ratio of the fluctuation and the dissipation as an effective temperature. In the case of the SET it has been shown that the typical value of k B T eff is of the order of eV . 13 The Langevin equation ( 23) can then be solved by Fourier transform giving

S x (ω) = x(ω)x(-ω) = F 2 0 S n (ω) m 2 |ω 2 m -ω 2 -iγω| 2 (25)
and in particular in the low-frequency limit:

S x (ω = 0) = F 2 0 S n (ω = 0) m 2 ω 4 m . ( 26
)
We can now use the expansion (7) to find the lowest order contribution of the stochastic fluctuations of x(t) to the current. We denote these fluctuations δx(t) to distinguish them from the time-dependent average induced by the external driving:

δI(t) = ∂I ∂n g V g0 C g e δx(t) + . . . . (27)
The back-action current noise is then

S ba I (ω) = 2 ∂I ∂n g F 0 n g kL 2 S n (ω = 0) . ( 28
)
As discussed in Refs. 19 and 41 the mechanical backaction noise can be very strong and induce effective giant Fano factors.

Finally the measurement added noise can be obtained as is done for the amplifiers, 42 by dividing the fluctuation of the current signal by the amplifier gain squared. This gives:

S add x = S mx λ 2 = S shot I + S ba I 4λ 2 . ( 29
)
This quantity gives the upper bound on the detection sensibility, since the limitations considered are intrinsic to the detection method. We will evaluate explicitly these quantities for two specific models in sections V and VI.

IV. FAST OSCILLATOR

In this section we relax the condition ω m Γ for the calculation of the mixing current. We assume Γ, ω m k B T , the electronic transport is then decribled by sequential transport and we will find the mixing current to lowest non-vanishing order in the amplitude of the oscillating field by making use of a master equation description.

Let's begin by introducing in some details the electron tunnelling description. We assume that the only available charge states on the island are those associated with two charge states N e and (N + 1)e. We will call these two states 0 and 1. The state of the SET is thus fully described by the probabilities of one of these two state to be realized: P n , with n = 0, 1. We define Γ L+(-) as the rate for adding (subtracting) one electron on (from) the central island through the left tunnel junction. Similarly we define Γ R+(-) for the right junction. We define also

Γ α = Γ Lα + Γ Rα , with α = ±, Γ L = Γ L+ + Γ L-, Γ R = Γ R+ + Γ R-
, and Γ T = Γ + + Γ -. The master equation for the the probability reads ( Ṗ ≡ dP/dt):

Ṗ0 = -Γ + P 0 + Γ -P 1 (30) Ṗ1 = Γ + P 0 -Γ -P 1 . ( 31
)
Using the conservation of probability (P 0 + P 1 = 1) we are left with

Ṗ0 = -Γ T P 0 + Γ -. (32)
We consider now that the rate equations are modulated by two oscillating parameters, in our specific case V and n g . We expand in power series of the amplitude of oscillation the rates keeping only the lowest orders:

Γ α (t) = Γ α(0) (t) + Γ α(1) (t) + Γ α(2) (t) + . . . ( 33
)
where α stands for any of the previously introduced labels, and the term into parenthesis indicates the order in the expansion. As far as the driving frequency is smaller than the temperature, ω i k B T , the explicit expression of the time-dependent rates can be obtained by that for the static case by substituting the time-dependent fields: 43

for instance Γ α (t) = Γ α (a(t), b(t)), where a = a 0 + a 1 (t), b = b 0 + b 1 (t), and a 1 (t) = a 1 cos(ω 1 t), b 1 (t) = b 1 cos(ω 2 t).
One can then expand to second order in the time dependent part of the two parameters to obtain:

Γ α (t) = Γ α + ∂Γ α ∂a a 1 (t) + ∂Γ α ∂b b 1 (t) + 1 2 ∂ 2 Γ α ∂a 2 a 2 1 (t) + ∂ 2 Γ α ∂a∂b a 1 (t)b 1 (t) + 1 2 ∂ 2 Γ α ∂b 2 b 2 1 (t) + . . . .
The expansion up to second order can then be rearranged in a Fourier series:

Γ α (t) = Γ α 00 + n=-1,1 Γ α(1) n,0 e inω1t + Γ α(1) 0,n e inω2t + Γ α(2) 1,-1 e i(ω1-ω2)t + cc + . . . ( 34
)
where the static part Γ α 00 has contributions of zero and second order in the driving fields. The notation Γ α(p) n,m indicates a contribution of order p in the driving intensity. Concerning the time dependent second order terms, we keep only the interesting part at the mixing-current frequency ω ∆ .

We look for a solution of the master equation in terms of the stationary Fourier components

P 0 (t) = n,m
A nm e i(nω1+mω2)t .

(

) 35
This gives for each Fourier component the equation:

(inω 1 + imω 2 )A nm + n ,m Γ T n m A n-n ,m-m -Γ - nm = 0 .
(36) We further expand the A coefficients writing:

A nm = ∞ p=0 A (p) nm , (37)
where again p indicates the order in the driving fields. This leads to a set of equations that can be solved recursively. The zeroth-order one reads:

(inω 1 + imω 2 )A (0) nm + Γ T (0) 00 A (0) n,m -Γ -(0) 00 δ n,0 δ m,0 = 0 . (38)
It gives immediately the static solution:

A (0) nm = δ n,0 δ m,0 Γ -(0) 00 Γ T (0) 00 . ( 39
)
For the next two orders we obtain:

A (1) nm = Γ T 00 Γ -(1) nm -Γ T (1) nm Γ -(0) 00 Γ T (0) 00 (inω 1 + imω 2 + Γ T (0) 00 ) , (40)
and

A (2) nm = Γ -(2) nm -Γ T (2) 00 A (0) n,m -n m Γ T (1) n m A (1) n-n ,m-m inω 1 + imω 2 + Γ T (0) 00 . ( 41
)
The non-vanishing terms up to order two are , andA (2) 0,±2 . As usual for the Fourier transform of real functions the following relation holds: A * n,m = A -n,-m . Let us now consider the particle current. It can be expressed in terms of P and Γ, for instance, on the left junction (note that this expression does not include the displacement current):

A (0) 0,0 , A (1) ±1,0 , A (1) 0,±1 , A (2) 0,0 , A (2) ±2,0 , A (2) ±1,±1
I(t)/e = Γ L+ P 0 -Γ L-P 1 = Γ L P 0 -Γ L-. (42)
Substituting the expansion (35) into Eq. ( 42) we obtain for I a similar expansion to Eq. (35). The first three orders read:

I (0) nm /e = Γ L(0) 00 A (0) 00 -Γ L-(0) 00 δ nm δ n0 (43)
I (1) nm /e = Γ L(1) nm A (0) 00 + Γ L(0) 00 A (1) nm -Γ L-(1) nm (44)
I (2) nm /e = Γ L(0) 00 A (2) nm + n m Γ L(1) n-n ,m-m A (1) n m +Γ L(2) nm A (0) 00 -Γ L-(2) nm . (45)
The mixing current is given by

I c mx = ReI 1,-1 /2 , I s mx = -ImI 1,-1 /2 . ( 46
)
In order to simplify the expressions obtained above we use the fact that in general ω 1 ≈ ω 2 ≡ ω D so that even in the fast oscillator limit

|ω 1 -ω 2 | Γ T (0)
00 . This gives the approximate expressions:

A (1) 10 = Γ T 00 Γ -(1) 10
-Γ T (1) 10 Γ -(0) 00 Γ T (0) 00 (iω D + Γ T (0) 00 ) (47)
A (2) 1,-1 = Γ -(2) 1,-1 -Γ T (1) 1,0 A (1) 0,-1 -Γ T (1) 0,-1 A (1) 1,0 Γ T (0) 00 (48)
One can see that the residual ω D -dependence is due to the relaxation time of the charge in the island. As expected it disappears for ω D Γ T (0) 00 . The contribution from

I (1) 1,-1 vanishes since Γ α(1)
1,-1 = 0. The interesting part is the contribution of second order which reads:

I (2) 1-1 = Γ L(0) 00 A (2) 1,-1 + Γ L(1) 1,0 A (1) 0,-1 + Γ L(1) 0,-1 A (1) 1,0 +Γ L(2) 1,-1 A (0) 00 -Γ L-(2) 1,-1 . (49)
One can verify that for ω Γ T (0) 00 expression Eq. (49) reduces to ∂I 2 /∂a∂b recovering the standard results for the mixing-current [cf. expressions (10) and ( 9)].

In the opposite limit of ω Γ T (0) 00

the first order correction to the charge variation vanishes (A (1) 1,0 → 0): the charge has not the time to follow the driving. Only a second order correction survives

A (2) 1,-1 = Γ -(2) 1,-1 /Γ T (0) 00 .
The residual time dependence at the mixing frequency is only due to the direct modulation of the tunneling rates (Γ α(2) 1,-1 ). The final expression for I 1,-1 in the limit ω → ∞ reads:

I (2) 1-1 fast = Γ L(0) 00 Γ -(2) 1,-1 Γ T (0) 00 + Γ L(2) 1,-1 A (0) 00 -Γ L-(2) 1,-1 (50)
In the following two sections we consider explicitly the case of a metallic dot and of a single electronic level dot and we derive explicit expressions for the mixing current, its fluctuation and the response function in the highfrequency regime.

V. THE METALLIC DOT SINGLE-ELECTRON TRANSISTOR

The expression for the tunnelling rate are well known for a metallic dot in the Coulomb blockade regime. [START_REF] Ingold | Single electron tunnelling[END_REF] For convenience of the reader, we report in the appendix a very short derivation of the electrostatic relations. We consider only the two states with N and N + 1 electrons.

A. Low temperature case

We begin by discussing the low temperature case k B T eV E C where E C = e 2 /2C Σ is the Coulomb energy. In this case there are only two non-vanishing rates (for V > 0)

Γ + L (N ) = Γ o (v + ñg )θ(v + ñg ) (51)
Γ - R (N + 1) = Γ o (v -ñg )θ(v -ñg ) (52)
where

Γ o = 1/RC Σ , v = (C + C g /2)
V /e and ñg = C g (x)V g /e -N -1/2, we assume a symmetric device with tunneling resistance R. The stationary solution to the master equation ( 32) and the stationary current ( 42) read

P st 1 = ñg + v 2v , I = eΓ o v 2 -ñ2 g 2v , (53)
both equations valid for |ñ g | < v. The current vanishes continuosly for |ñ g | ≥ v while the probability is 1 for ñg > v and 0 for ñg < -v.

The driving amplitudes in terms of the dimensionless variables introduced read v 1 and n g1 . Note that the dependence of the rates on v and ñg is non-analytic for ñg = ±v, this gives a constraint on the amplitude of the oscillations since the Taylor expansions are not valid if the parameters cross this values. This gives the constraints |ñ g ± n g1 | < v and v -v 1 > ñg , that can be written n g1 , v 1 < v -ñg . Using Eq. (51) and Eq. ( 52) we can readily obtain the non-vanishing coefficients of the expansion (34):

Γ L+ 00 = Γ o (v + ñg ), Γ L+ 10 = Γ o e iϕ n g1 /2, Γ L+ 01 = Γ o v 1 /2, Γ R- 00 = Γ o (v -ñg ), Γ R- 10 = -Γ o e iϕ n g1 /2, Γ R- 01 = Γ o v 1 /2.
For ω 1 ≈ ω 2 = ω D we obtain a very simple expression for the component I 1,-1 :

I 1-1 = eΓ o ñg v 1 n g1 e -iϕ ω2 D + 4v 2 (54)
here we defined ωD = ω D /Γ o . One finds thus a Lorentzian behaviour, the amplification factor decreases quite rapidly for large frequency driving ω D . The main reason for the reduction of sensitivity is the incapacity of the charge in the dot to follow the driving signal. The crossover value for the frequency is ω D ≈ V /Re, above this value one cannot use anymore the adiabatic approximation for the relaxation of the charge on the dot. It simply coincides with the frequency for which one electron per driving period crosses the device. For instance for ω m = 100 MHz, R = 10 5 Ohm, for voltage below a mV the corrections due to the retardation of the charge on the dot becomes relevant This regime has been observed in the experiment presented in Ref. 45, where the crossover from slow to fast oscillator has been investigated by a fine tuning of the tunnelling resistances.

The amplification factor for the mechanical quadratures is thus:

λ = eΓ o L n g ñg v 1 (ω 2 D + 4v 2 ) . ( 55
)
It is maximum for ñg = ±v, but one should also take into account the constraint on the amplitude of v 1 < v -|ñ g |. One way to take that into account is to set v 1 = v -|ñ g |, this is the maximum allowed value for the driving amplitude, and since the signal increases linearly with v 1 , it gives the maximum value for λ. This gives:

λ = eΓ o L n g ñg (v -|ñ g |) (ω 2 D + 4v 2 ) . ( 56
)
The maximum of λ as a function of the gate voltage is obtained for ñg = ±v/2 and its value (for ωD v) is

λ = eΓ o n g 16L (57)
independently of v. For a typical device one has n g ≈ 100 L ≈ 1µm, Γ 0 = 10 11 Hz leading to λ ∼ 0.1 A/m. 9,10

The gain is only a part of the detection, one has also to evaluate the noise. For that we need the two contributions considered in the section III. The Fano factor has been obtained in Ref. 36 (cf. Eq. 41 there):

F = Γ + L 2 + Γ - R 2 (Γ + L + Γ - R ) 2 = v 2 + ñ2 g 2v 2 (58)
it varies between 1/2 and 1. The shot noise becomes thus:

S shot II = e 2 Γ o v 4 -ñ4 g 2v 3 . ( 59
)
To obtain the contribution of the displacement fluctuation we need to calculate the charge noise correlation function: S n (t) = δn(t)δn(0) . This symmetric correlator can be obtained by the conditional probability P (1t|10) that the dot it occupied at time t > 0 with the condition that it was occupied at time 0: S n (t) = [P (1t|10) -P st 1 ]P st 1 . Solving the master equation with the initial condition P 1 = 1 one finds

P (1t|10) = 1 + (P st 1 -1)(1 -e -Γ T t ), (60)
with P st 1 = Γ + /Γ T . By Fourier transforming we obtain:

S n (ω) = P st 1 (1 -P st 1 ) 2Γ T ω 2 + Γ T 2 (61)
As expected the correlation function is flat for ω Γ T , the required low frequency correlator reads then:

S n (ω = 0) = 2Γ + Γ -/(Γ T ) 3 . ( 62
)
In the specific case of low temperature one obtains thus

S n = (v 2 -ñ2 g )/(4Γ o v 3 ).
In the typical working regime of a SET V V g , and n g ≈ N . Using the Eq. (A1) one finds that Q g /e ≈ n g ≈ N . We thus have F 0 = 2N E C /L. Collecting all the terms we can substitute into Eq. ( 28) to obtain:

S ba I = 2e 2 Γ o N 4 E c kL 2 2 (v 2 -ñ2 g )ñ 2 g v 5 . ( 63
)
The ratio of the mechanical to the shot noise is thus:

S ba I S shot I = E c kL 2 2 4N 4 ñ2 g v 2 (v 2 + ñ2 g ) . ( 64
)
For large mechanical coupling (L small and N large) the mechanical noise dominate even if for small ñg it is always suppressed, due to the vanishing of ∂I/∂n g . From Eqs (57), (59), and (63), we obtain the seeked added noise as defined by Eq. ( 65). In order to study its dependence on the different parameters it is convenient to introduce the two dimensionless variables ν ≡ ñg /v and the dimensionless coupling constant δ ≡ (E C /kL 2 )(N 2 /v) = P /eV , where P = F 2 0 /k is the energy scale of the electromechanical coupling. 13,14,18 The added noise then reads:

S add x = E c kΓ 0 f (ν, δ) , (65)
with

f (ν, δ) = 2(1 + ν)[4(δ 2 + 1)ν 2 + 1] ν 2 (1 -ν)δ , (66)
and 0 < ν < 1. The function diverges for ν → 1 due to the fact that we have to limit the amplitude of the voltage modulation and diverges for ν → 0 due to the vanishing of the amplification factor. The minimum added noise is thus always for values of ν between 0 and 1. In the weak coupling limit, for δ 1, one finds that the minimum is at ν ≈ 0.54 and reads

S add x ≈ 14.8 E c /k Γ 0 eV P . ( 67
)
For strong coupling, δ 1, instead the minimum is close to ν = 0 with a value for

S add x ≈ 8 E c kΓ 0 P eV 2 . ( 68
)
In both cases the noise diverges when δ becomes very small or very large. In the weak coupling limit the added noise is dominated by the current noise (imprecision noise), in the strong coupling it is instead given essentially by the back-action noise. As usual 42 the optimal situation is in the middle for δ ≈ 1.

In Fig. 2 we plot S add x [δ, ν m (δ)], where ν m (δ) is the value of ν that minimizes S add xx for given δ. We thus find that the absolute minimum for the added noise is obtained for ν = 0.32 and δ = 1.857 and reads S add xx = 57.61 This is the ultimate sensitivity that can be obtained with this device in ideal conditions, when all other sources of imprecisions have been eliminated. Inserting typical values of E c ≈ 10K, Γ 0 ≈ 10 11 Hz, k = 10 -5 N/m one obtains the value of S add xx ≈ 10 -26 m 2 /Hz. One should regard this value with some caution. Let's consider the value of the coupling that is required to obtain this sensitivity. The optimal value of δ is for eV ≈ P . As discussed in the litterature (see for instance Ref. 18, where this energy is called E E ) this scale determines the value at which the system undergoes a current blockade. It is difficult to reach this limit (since one needs also k B T P ) in metallic SETs. On the other side P of the order of 0.3 K has been observed in suspended carbon nanotubes. 45 The dramatic effects expected at low tempereture on the mechanical resonators have been discussed recently. 23,24 This extreme limit need to be reconsidered, since the resonating frequency of the resonator is renormalized by the coupling, and the added noise induced by the oscillator is expected to be more effective. In particular the oscillator becomes strongly non-linear close to the transition.

E c kΓ 0 ( 69
)
We can estimate in a simple way the effect of the softening of the mechanical resonator following Ref. 18. The correction to the variation of the energy reads 11

∆E ± → ∆E ± ± F 0 x (70)
this changes the form of P st 1 given by Eq. (53) as follows:

P st 1 = ñg + xF 0 /(2E C ) + v 2v . (71)
Substiting into the equation for the average force F 0 P 1 st and taking the derivative with respect to x one obtains the renormalized spring constant:

k = k(1 -δ) . (72)
The instability appears for δ = 1, where two new stable solutions bifurcate. The only change in our previous calculations is the value of k entering Eq. ( 65):

S add xx = E c kΓ 0 1 1 -δ f [ν, δ/(1 -δ)] . (73)
Repeting the minimization procedure we find that the minimum is now for δ = 0.48 holding the value of 132.7(E c /kΓ 0 ). Thus the renormalization of the resonating frequency reduces the precision of a factor of 2, leaving space for high sensitivity detection. The actual limitation in current experiments will be the value of the coupling, since in practice the typical temperature reached in experiments on metallic quantum dot is much larger than P . In the following section we considerthe detection at finite temperature and low voltage.

B. Finite temperature case

Let us now consider the finite temperature case eV k B T E C . In this case we have to take into account the four possible tunnelling processes that change the charge on the dot from the N to the N +1 state (cfr. Appendix). The respective rates read:

Γ + L(R) (N ) = Γ Th h[(∓eV -2ñ g E C )/k B T ], (74)
Γ + L(R) (N + 1) = Γ Th h[(±eV + 2ñ g E C )/k B T ], (75)
with h(y) = -y[1-e y ] and Γ Th = k B T /e 2 R. We consider the low bias voltage limit eV /k B T 1. In this limit the expression for the current Eq. ( 42) becomes:

I = eΓ Th eV 2k B T g(2ñ g E C /k B T ), (76)
where

g(y) = h + h -+ h + h - h + + h - = e y y e 2y -1 (77)
and h ± = h(±y). From the expression of the current we obtain

∂ 2 I ∂V ∂ ñg = - E C Rk B T g (2ñ g E C /k B T ), (78)
(for brevity, we omit in the following the arguments of g and of the other functions of y = 2ñ g E C /k B T ) with the amplification factor:

λ = n g0 Γ 0 eV 1 4Lk B T g ≈ F 0 V 1 8Rk B T g . (79)
The factor g (x) has a maximum for y = 1.16 for which it holds the approximate value 0.154. Thus tuning ñg = 0.58k B T /E C allows to obtain the maximum value of the amplification factor. Comparing this value to Eq. (57), valid for k B T eV , we see that the amplification factor is reduced by the term eV 1 /k B T 1. The shot noise at low frequency reads: 36

S shot I = e 2 Γ + L Γ - R + Γ - L Γ + R Γ T -2 (Γ + L Γ - R -Γ - L Γ + R ) 2 Γ 3 T .
(80) For small V the first term (thermal noise) dominates and gives:

S shot I = e 2 Γ Th h + h - h + + h - . (81)
The charge noise in the same limit reads

S n = 1 Γ Th h + h - (h + + h -) 3 . ( 82
)
From the expression of the back-action noise (28) we see that for V → 0 there is (apparently) no back action of the measurement. It is possible to set V = 0 and exploit its modulation around 0 to detect the motion of the oscillator. But in this case we need to consider the next order in the expansion (7). For V = 0 we have:

δI = ∂I ∂n g ∂V δn g V 1 + . . . . (83)
From this we have for the current-current correlator:

δI(t 1 )δI(t 2 ) = ∂G ∂n g 2 V 1 (t 1 )V 1 (t 2 ) δn g (t 1 )δn g (t 2 ) , (84)
where G = dI/dV for V = 0 is the conductance. The product of the two V 1 terms gives an oscillating term depending on t 1 + t 2 that averages to zero and a second term proportional to cos[ω 2 (t 1 -t 2 )]. Using δn g (t) = (C g V g /e)δx(t) we have

S ba I = 1 2 ∂G ∂n g 2 (V g0 C g V 1 /e) 2 S x (ω 2 ) . (85)
Typically ω 2 ≈ ω m , we thus assume that it is resonant in order to evaluate the case of maximal back-action:

S ba I = g 2 h + h - 16(h + + h -) 3 e 2 Γ 0 2 P Q 2 (eV 1 ) 2 (k B T ) 3 E C , (86)
with h ba = (g ) 3 and Q = ω m /γ the oscillator quality factor.

2 h + h -/(h + + h -)
Adding the two sources of current noise Eq. (86) and Eq. ( 81) we obtain for the added noise:

S add x = E C kΓ 0 α ba P Q 2 k B T + α shot (k B T ) 3 P (eV 1 ) 2 , (87)
with the numerical factors α ba = 4h + h -/(h + + h -) 3 and

α shot = 32h + h -/[(g ) 2 (h + + h -)].
Choosing the value ñg = 1.60 that maximizes λ their values are α ba = 0.23 and α shot = 449.

The minimum of the added noise is obtained for

P = α shot α ba 1/2 (k B T ) 2 QeV 1 , (88)
with a minimum noise of

S add x = 2 E C kΓ 0 α ba α shot 1/2 Q k B T eV 1 . (89)
Since eV 1 /k B T 1, at best we can set this ratio to 0.1. This gives for the optimal value of the coupling

P k B T ≈ 441 Q (90)
and the minimum of the added noise

S add x = 203 QE C kΓ 0 . ( 91
)
Some comments are at order. First we assumed that the frequency driving the voltage bias is resonant with the oscillator. This is un upper limit to the back action, in particular if Q 1 this condition is not fulfilled and the back action will be reduced. For the nonresonant case it is sufficient to use the above results with Q ≈ ω m /ω ∆ , reducing enormously the minimum added noise, to the expenses of finding a much larger coupling constant. The second comment concern the value of the coupling constant P necessary to reach the minimum. One can see that even with the assumption of resonant back action it is relatively large. For a typical Q ≈ 10 4 one finds P /k B T ≈ 0.04. To our knowledge the largest value of the ratio k B T / P is ≈ 0.017 has been reported in Ref. 45. Since as soon as Q 1 it is possible to avoid resonant back-action, in most cases the main limitation is to reach large values of P .

It is interesting to compare the shot-noise contribution of the added noise with the resonant brownian motion fluctuations:

S B x (ω m ) = 2 k B T kγ . (92)
The ratio reads:

S add x S B x = α shot 2 E C γ P Γ 0 k B T eV 1 2 . ( 93
)
Detection of brownian motion can then be done for When the temperature and the voltage bias is much smaller than the electronic level separation the rates for electron transfer reads: 46

P /E C > 2 × 10 4 γ/Γ 0 (
Γ + L(R) = Γ L(R)0 f F [( -µ L(R) )/k B T ], (94)
Γ - L(R) = Γ L(R)0 1 -f F [( -µ L(R) )/k B T ] , (95)
where f F (y) = 1/(1 + e y ) if the Fermi function, is the level position, µ L (R) is the left (right) chemical potential, and Γ L(R)0 are the transfer rates. For simplicity in the following we choose Γ L0 = Γ R0 = Γ 0 . The modulation of the gate voltage leads to the time-dependence (t) = 0 + 1 (t) of the electronic level energy with

0 = d0 + eC g V g0 /C Σ , (96)
1 (t) = e[C g V g0 x(t) + C g V g1 (t)]/C Σ , (97)
and d0 the position of the electronic level for vanishing V g . We assume symmetric bias so that the chemical potential read:

µ L(R) (t) = µ L(R)0 + (-)e(V + V 1 cos ω 2 t)/2 . ( 98
)
Following the steps of the previous section we can calculate the current

I = eΓ 0 2 [f F [( -µ L )/k B T ] -f F [( -µ R )/k B T ]] (99)
from which we obtain for vanishing V the amplification factor:

λ = en g0 Γ 0 4L eV 1 E C (k B T ) 2 f F (y), (100)
where the argument of the Fermi function is y = ( 0µ)/k B T , and will be omitted in the following. The maximum of f F is obtained for y = 1.31 with a value of 0.096. The thermal part of the shot noise and the charge noise read:

S shot I = e 2 Γ 0 f F (1 -f F ), (101)
S n = f F (1 -f F ) Γ 0 . (102)
Using Eq. ( 85) for the back-action noise we obtain

S ba I = f F (1 -f F )f F 2 8 e 2 Γ 0 eV 1 Q P (k B T ) 2 2 . ( 103
)
The added noise has thus the form:

S add x = k B T kΓ 0 α ba Q 2 P k B T + α shot k B T eV 1 2 k B T P (104) with α ba = 2f F (1 -f F ) and α shot = 16f F (1 -f F )/f F 2 .
Their values for y = 1.31 are α ba = 0.34 and α shot = 289.2. We find the same value of P for the minimum of the added noise in the metallic case [cf. Eq. ( 88)], but the minimum of the noise has a different expression:

S add x = 2 α ba α shot 1/2 (k B T ) 2 kΓ 0 eV 1 . (105)
Essentially the energy scale of the Coulomb blockade is substituted by the temperature, in principle reducing the added noise. The conclusion is that the single-level SET should allow a better resolution of the metallic SET by a factor E C /k B T .

VII. CONCLUSIONS

In this work we have studyed theoretically the sensitivity of the mixing-current technique. We first found general expressions valid when the oscillator resonating frequency is comparable or larger of the transfer rate of electrons. We find that a reduction of the amplification factor of the order of (Γ 0 /ω D ) 2 is expected. This effect should be relatively small in most practical experimental realizations. We then analysed the fundamental limitations due to the intrinsic noise present in the (current) signal and the effect of the back-action fluctuations. On general grounds one finds that an optimal value of the electromechanical coupling ( P ) exists that minimizes the added noise. This value is larger than what is realized in the present experiments, showing that increasing the coupling allows to reach higher sensitivity. At finite temperature the relevant parameter is the ratio P /k B T and values of the order of 1 are needed to reach the optimal minimum added noise. At vanishing temperature the relevant parameter is instead P /eV . In all cases the scale of the sensitivity is given by E C /Γ 0 k. Optical means can detect CNTs displacement with good accuracy, even if the small size of the object does not allows to reach the spectacular sensitivity obtained with macroscopic mirrors. A sensitivity of 5 • 10 -22 m 2 /Hz as been reported 47 by cavity-enhanced optical detection of CNTs.

We considered only classical fluctuations. It seems difficult to use the mixing technique to reach the quantum limit of detection, since the effective temperature of the oscillator, even at vanishing temperature, is of the order of eV that typically needs to be larger than ω m . On the other side it may be instructive to compare the sensitivity found at vanishing temperature with the zero point fluctuations spectrum at resonance: S SQL x = 2 ω m /kγ. One sees that the ratio to the typical mixingcurrent technique added noise at zero temperature is 10 -2 ( ω m /E C )(Γ/γ)( P /E C ), since Γ/γ 1, for sufficiently large P the added noise can be of the same order of the zero-point fluctuations.

We conclude that the sensitivity of the mixing technique can still be improved by increasing the electromechanical coupling till reaching P of the order of the temperature or the Coulomb blockade energy where the backaction will be of the same order of the intrinsic current noise of the device. In this appendix we present, mainly for clarifying the notation, a brief derivation of the electrostatic energy variation for the tunnelling of an electron in a single electron transistor. [START_REF] Ingold | Single electron tunnelling[END_REF] The electric scheme is presented in Fig. 3 where the potentials of the left, right, and gate leads are defined as V L , V R , and V g , respectively. In the same way the charge on each capacitance (on the leads side) is indicated with Q i with i = L, R, and g. Defining V I the potential of the island one has

Q i = (V i -V I )C i . (A1)
Summing the three equations one obtains immediately the expression for the potential on the island:

V I = i C i V i + Q /C Σ , (A2)
where C Σ = i C i and Q =i Q i is the total charge on the island. The total electrostatic energy E e (Q) = i Q2 i /2C i = Q 2 /2C Σ + constant, where the constant term does not depend on Q. From Eq. (A1) and Eq. (A2) one then finds that adding a charge q on the island will change the charge on each capacitor plate of

δQ i = - qC i C Σ . ( A3
)
The total electrostatic energy variation (final energy minus initial energy) for the transfer of an electron from the left electrode on the island is then given by the variation of the total electrostatic energy plus the the work done by the voltage sources (i V i δQ i , with δQ i = eC i /C Σ for i = L and δQ L = eC L /C Σ -e):

∆E + L = E e (Q -e) -E e (Q) -e i V i C i C Σ + eV L . (A4)
The general expression reads then:

∆E ± L,R = -e (-e ± 2Q) 2C Σ ∓ e C Σ ( i V i C i -C Σ V L,R
) .

(A5) The variation of the energy depends only on the difference of the three potentials, we can thus choose to express it in terms of V = V R -V L and V g = V g -(V L + V R )/2. For simplicity we write the expressions in the symmetric case of C L = C R = C:

∆E ± L = e 2C Σ (e∓2Q)∓ e C Σ (C V + C g V g ) (A6
)

∆E ± R = e 2C Σ (e∓2Q)∓ e C Σ (-C V + C g V g ) (A7)
with C = C + C g /2. Typically V is very small, while V g can be very large, in particular V g C g /e = n g is normally regarded as finite, while C g → 0 and V g → ∞. For this reasons we can normally neglect the displacement dependence induced by C g (x) in C or C Σ , while it is necessary to keep the x dependence in C g (x) that appears in the expression C g V g .

Let now focus on the four energy variations associated with the change of the number of electrons in the dot between the two states N and N + 1. We need ∆E + L,R (N ) = -∆E - L,R (N +1) that can be explicitly written as:

∆E + (N ) L,R = - e 2 C Σ (n g ± v -N -1/2), (A8)
with n g = C g V g /e and v = C V /e. The expression of the tunneling rate is obtained then by the Fermi golden rule:

Γ α± = k B T e 2 R α h(∆ ± α /k B T ) (A9)
with h(x) = -x/(1 -e x ). In particular for T → 0 the expression for the rate becomes simply Γ α± = -∆ ± α /e 2 R α θ(-∆ ± α ). These expressions allow to obtain the tunneling rates necessary for the calculations presented in the main text of the paper.

Figure 1 .

 1 Figure 1. Schematic of the typical experimental set-up used to measure the displacement of a mechanical oscillator by detection of the mixing current (adapted from Ref. 28)

Figure 2 .

 2 Figure 2. S add xx as a function of δ = P /eV for ν(δ) minimizing the function . In a inset the value of ν that minimizes the function for given δ [νm(δ)].

Figure 3 .

 3 Figure 3. Electric scheme of a single electron transistor

  where we assumed as before eV 1 /k B T = 0.1). For instance in Ref. 10 γ/Γ 0 ≈ 10 -8 allowing the detection of the brownian motion fluctuations even for very weak coupling. For a rough estimate of the coupling in that experiment one can use the expression given in Ref. 23 P /k B T ≈ 2δω m /ω m , where δω m is the modulation of the resonating frequency near the degeneracy point (see Fig.3in Ref. 10). For Ref. 10 one finds P ≈ 16m K to be compared to E C of the order of 10K. Notwithstanding the low value of the coupling constant, the resolution is largely sufficient to detect the Brownian motion of the carbon nanotube.

	VI. THE SINGLE-ELECTRONIC LEVEL SET

K. L. Ekinci, X. M. H. Huang, and M. L. Roukes, Appl. Phys. Lett. 84, 4469 (2002).

B. Lassagne, D. Garcia-Sanchez, A. Aguasca, and


R M Chandima Ratnayake
email: chandima.ratnayake@uis.no

Inspection Interval Estimation: A Fuzzy Logic Based RBI Analysis Approach

Keywords: Risk based inspection analysis, production system, inspection interval, fuzzy inference system

Risk based inspection analysis (RBIA) on offshore oil and gas (O&G) production systems optimize level of in-service inspection. The potential failure risk of a system, a sub-system or a thickness measurement location (TML) of O&G production systems comprise the consequence of failure (CoF) and probability of failure (PoF). A tailor-made risk matrix supports the estimation of maximum inspection intervals. When the inspection intervals are calculated using a risk matrix, suboptimal classification tends to occur as there are no means to incorporate actual circumstances at the boundary of the input ranges or at the levels of linguistic data and risk categories. This manuscript suggests a fuzzy inference system (FIS) to overcome the aforementioned. Membership functions and the rule base development have been carried out in alignment with a tailor-made risk matrix which has been utilized by a production plant owner operator organization. A rule view and a calculation result have been demonstrated to illustrate the methodology.

Introduction

The performance of production and/or process plants are sustained by appropriate inspection planning and scheduling [START_REF] Ratnayake | Maintaining Technical Integrity of Petroleum Flowlines on Offshore Installations: A Decision Support System for Inspection Planning[END_REF][START_REF] Wang | A Model of Multiple Nested Inspections at Different Intervals[END_REF]. However, in this context, it is vital to determine the optimal inspection intervals in terms of a criteria of interest [START_REF] Ratnayake | A Decision Model for Executing Plant Strategy: Maintaining the Technical Integrity of Petroleum Flowlines[END_REF][START_REF] Wang | A Model of Multiple Nested Inspections at Different Intervals[END_REF]. Consequently, risk based inspection analysis (RBIA) has been accepted over the last few years as a method for prioritizing the in-service inspection of a plant as well as for estimating corresponding inspection intervals [START_REF] Ratnayake | Maintenance Integrity: Managing Flange Inspections On Aging Offshore Production Facilities[END_REF]. These methods have been developed nationally (e.g. American Petroleum Institute (API), a number of private organizations (particularly in the petrochemical industry), etc.) [START_REF] Geary | Risk Based Inspection -A Case Study Evaluation of Onshore Process Plant[END_REF] and internationally (e.g. RIMAP -Risk based inspection and maintenance procedures for European industries) [START_REF] Kauer | Plant Asset Management RIMAP (Risk-Based Inspection and Maintenance for European Industries) the European Approach[END_REF]. The importance of risk was recognized principally as an important measure in assuring system safety [START_REF] Brown | Risk-based Hazardous Protection and Prevention by Inspection and Maintenance[END_REF]. However, there is a fundamental challenge in the mathematical modeling of RBIA to perform optimum maintenance as a subject [START_REF] Scarf | On the Application of Mathematical Models in Maintenance[END_REF].

The mathematical modeling enables mitigating subjective judgments based on limited information and also, some of the inherent challenges present in the current RBIA [START_REF] Scarf | On the Application of Mathematical Models in Maintenance[END_REF][START_REF] Ratnayake | Utilization of Piping Inspection Data for Continuous Improvement: a Methodology to Visualize Coverage and Finding Rates[END_REF]. Alternatively, it mitigates the significant variability and discrepancy present in the current inspection interval estimations. For instance, a report published on a case study evaluation of an onshore process plant revealed that "subjec-tive judgments based on limited information did lead to some significant differences in inspection periods" [START_REF] Geary | Risk Based Inspection -A Case Study Evaluation of Onshore Process Plant[END_REF]. It also revealed that although "generally, the inspection periods reflected the assessed risk", "considerable scatter was apparent in the data and some participants exhibited greater conservatism in their assessments than others" [START_REF] Geary | Risk Based Inspection -A Case Study Evaluation of Onshore Process Plant[END_REF]. The same report suggests that "software, expert systems and expert judgment all have merits, greater integration of these elements might be beneficial" [START_REF] Geary | Risk Based Inspection -A Case Study Evaluation of Onshore Process Plant[END_REF]. Hence, it is vital to develop expert systems to support expert judgments and alternatively to develop sophisticated software to minimize the variability present in estimating inspection periods (i.e. inspection intervals).

This manuscript proposes a fuzzy logic based expert system for estimating inservice inspection intervals [START_REF] Ratnayake | Plant Systems and Equipment Maintenance: Use of Fuzzy Logic for Criticality Assessment in NORSOK Standard Z-008[END_REF]. The estimation of in-service inspection intervals is based on PoF, CoF and currently established values of inspection intervals with respect to different risk levels.

Industrial Challenge

Currently, recommended practices, standards (e.g. DNV-RP-G101), operator company procedures, etc., provide decision matrices for estimating inspection intervals (i.e. time to inspect) [START_REF]Risk Based Inspection of Offshore Topsides Static Mechanical Equipment[END_REF]. However, when the classifications are carried out, there is no formal mechanism to incorporate data and information at the boundaries of the risk categories (i.e. alternatively at the boundaries of the ranges and levels of linguistic data). This is mainly due to the fact that there are no means to incorporate real data (qualitative or quantitative) in a consistent manner in estimating the maximum allowable time intervals. For instance, along a boundary, the spontaneous jumps of risk classification together with recommended inspection intervals (e.g. VH to H: the recommended inspection interval changes from 6 months to 48 months) hinder realistic values depending on the estimated PoF and CoF, leading to suboptimal inspection interval estimations (see Fig. 1). Fig. 1 illustrates such a risk matrix along with relevant maximum inspection intervals that have been utilized for piping RBISIA (i.e. in an operator organization which owns a production and process plant).

Due to the lack of a consistent approach, the inspection interval recommendations made are mostly confined to the PoF intervals, CoF intervals, and corresponding inspection interval values in an ad hoc manner dependent on the person who is involved in the analysis. Hence, it is vital to have a consistent approach to incorporate PoF, CoF and inspection intervals.

Methodology

In order to cater for rapid changes (in inspection interval) at the boundaries of each risk level (see Fig. 1), fuzzy membership functions have been introduced for each PoF, CoF and inspection interval. The introduction of fuzzy membership functions enables the inspection interval estimation to be made more realistic. Furthermore, a fuzzy inference system (FIS) enables the mitigation of discrepancies that may occur during the risk assessment process as a result of simultaneous considera-tion of different PoF and CoF ranges of values for estimating maximum inspection intervals.

Risk Matrix

During the detailed risk analysis process, the static mechanical pressure systems are subject to an investigation of PoF and CoF according to the categorization presented in Fig. 1. This is a 5x5 matrix, indicating levels of both PoF and CoF whilst providing five risk levels (i.e. VL, L, M, H and VH). These risk levels represent a combination of PoF and CoF based on the relevant numerical value ranges. The numerical value ranges for PoF and CoF have been retrieved from the documentation pertaining to the case study plant operator's organizations Fig. 1). In addition, the corresponding maximum inspection interval (in months) has been indicated within the parenthesis under each risk level. 6)

H (1E-2 -1E-3) (3 -7) VL (144) L (120) M (72) H (48) H (48) H (48) VH (6) VH (6) VH (6) M (1E-3 -1E-4) (7 -

Probability of Failure Assessment

In essence, the piping equipment is organized into corrosion groups which can contain several degradation mechanisms. Two models have mainly been used for evaluating the probability of failure (PoF) of piping equipment due to degradation: I. a susceptibility model for stainless steels; and II. a rate model for carbon steels.

Susceptibility models are used when the PoF is related to operating conditions. In this context, for a given set of conditions that are constant over time, the PoF also remains constant over time. This indicates that it is not easy to monitor the development of damage mechanisms by using inspection. Hence, actions are related to the monitoring of key process parameters, which are used as a trigger for inspection. DNV-RP-G101, Appendix A [START_REF]Risk Based Inspection of Offshore Topsides Static Mechanical Equipment[END_REF] provides guidance about typical materials and environmental conditions where this model is expected to be applicable and suggest values for PoF for typical conditions. For the susceptibility models there are two governing conditions required for degradation: I. wet environment; and II. temperature. The results from the monitoring of these conditions are the most important in setting the probability.

Rate models have increasing probability over the time, which makes it difficult to express as one probability value. Even though 'time to release' is itself not a probability expression, it explains the speed at which the probability increases, and therefore represents a useful profile of the probability. Hence, the development of degradation is measured by inspection. Then, the PoF is documented as an estimated 'time to release', based on the wall thickness and degradation rate on the area found to have the shortest time to release. In corrosion loops, both stagnant and varying flow conditions indicate very different estimated time to release (leak). Hence, the time to inspection is split among them in order to obtain the optimal time to inspection and this is reflected in the RBI analysis.

This manuscript provides a knowledge based engineering approach with the help of a fuzzy logic based inference system for estimating the inspection interval using a tailor-made risk matrix.

Fuzzy Logic Based Inference System

A 'pure fuzzy logic system' consists of a fuzzy rule base, which comprises a collection of fuzzy IF-THEN rules. These rules are utilized by the fuzzy inference engine to determine a mapping from fuzzy sets in the input universe of discourse U ⊂ R n to fuzzy sets in the output universe of discourse V ⊂ R based on fuzzy logic principles. The fuzzy IF-THEN rules follow the form below: R (1) :

IF x 1 is and … x n is THEN y is G 1 (1)
where and G j are fuzzy sets, x = (x 1 ,x 2 ,…, x n ) T ϵ U and y ϵ V are input and output linguistic variables which belong to the input and output universes, respectively, and j=1, 2,…, m. Practical experience reveals that these fuzzy IF-THEN rules provide a convenient framework to incorporate human expert knowledge. In Eq. ( 1), each fuzzy IF-THEN rule defines fuzzy set , … => G j for i = 1, 2,…,n, in the product space U × V. Expert opinions and data/information retrieved from different sources are taken into the mathematical model using the aforementioned rules. The main focus is to enhance the discerning power in the risk analysis process, whilst minimizing the uncertainties that may occur in dealing with the linguistic variables of the risk levels (i.e. H, VH, VL, etc.) at the boundaries of quantitative ranges. Essentially, membership functions (MFs) are developed with experienced personnel who are familiar with the risk analysis process [START_REF] Lapa | Effect Analysis Fuzzy Inference System in Nuclear Problems using Approximate Reasoning[END_REF]. Then, along with a rule base, MFs provide the possibility of recycling experts' knowledge in a consistent manner.

In 1975, Mamdani built one of the first fuzzy systems which used a set of fuzzy rules supplied by experienced human operators to control a steam engine and boiler combination [START_REF] Mamdani | An Experiment in Linguistic Synthesis with a Fuzzy Logic Controller[END_REF]. To date, Mamdani's approach has been successfully applied to a variety of industrial processes and consumer products [START_REF] Wang | Adaptive Fuzzy Systems and Control-Design and Stability Analysis[END_REF]. Fig. 2 illustrates the work process of the proposed fuzzy inference system. The PoF, time to release (TTR), safety consequence of failure, which is estimated by potential loss of life (PPL) (i.e. CoF PPL ) and economic consequence of failure (i.e. CoF Econ ) have been selected as the input variables. The inspection interval has been selected as the output. These variables consist of quantitative, qualitative and judgmental (i.e. linguistic) data. For each of the aforementioned variables, there is an associated membership function, which is established with the help of data, information and expert opinion [START_REF] Mamdani | An Experiment in Linguistic Synthesis with a Fuzzy Logic Controller[END_REF]. The fuzzification process aids fuzzifiying the inputs by determining the value of the membership functions corresponding to the different inputs. Furthermore, instead of restricting the user to a single, crisp, input value, this process allows an interval of values to be given, with values near the center of the interval being assumed to be 'more certain' than those near the edges, and the width of the interval indicating the amount of uncertainty present in the different input variables. The aforementioned has been achieved by associating appropriate membership functions (MF) for the input variables. Using an appropriate MF, the user has 'more confidence' that the input parameter lies relatively closer to the center of the interval than at the edges. In this study the author has incorporated triangular membership functions [START_REF] Juang | Design of Fuzzy PID Controllers using Modified Triangular Membership Functions[END_REF]. The fuzzy inference system (FIS) parameters were selected as follows: 'And' method with 'minimum', 'Or' method with 'maximum', 'Implication' with 'minimum', 'Aggregation' with 'maximum' and 'Defuzzification' with 'centroid' algorithm. Fuzzy rule bases were developed using the table-look-up approach (see Fig. 1). The toolbox simulator of the MATLAB (R2011a) tool was utilized to implement a fuzzy criticality ranking system [16].

4

Fuzzy Logic Based Modeling

Membership Functions for PoF

In essence, the PoF per year has been related to degradation mechanisms in the susceptibility model. At the same time, the TTR has been utilized in the rate models to express the speed at which the probability increases (note: the probability related to rate models is a variable with time). The PoF is also subjected to an analysis for both external and internal failure mechanisms. The membership functions (MFs) for inputs PoF and TTR have been developed based on data, information and expert opinions (see Fig. 3 and Fig. 4).

Membership Functions for CoF

The CoF PPL values are estimated by quantitative risk analysis. The COF Econ values are determined based on the inputs from relevant operation and process personnel (production loss) and from piping/vessels discipline engineers (repair costs). Fundamentally, COF Econ covers costs related to loss of production and costs for repair in the case of a leak. Consequence of loss of functionality is also considered as an economic loss due to reduced or lost production and/or major repair cost, not necessarily given a leak or shutdown of equipment. For instance, typical loss of functionality is mostly related to failure of internals in vessels and coolers. However, in this manuscript the values of CoF PPL and COF Econ have been retrieved from the tailor-made company specific topside inspection manual of the case study operator company. The MFs for inputs CoF PPL and COF Econ have been developed based on data, information and experts' opinions (see Fig. 5 and Fig. 6).

Analysis and Results

PoF and CoF PLL have been utilized to illustrate the calculation of MII (see Fig. 8). Using the table look-up approach (see Fig. 2) to estimate MII in relation to different PoF and CoF PLL levels, 45 rules have been generated. Fig. 8

Discussion

The suggested approach enables the quantitative values in a range to be distributed in relation to their significance. For instance, the left-hand side and right-hand side of a quantitative range have more relation to the previous and following qualitative category (i.e. VH, H, M, etc.) respectively. In the suggested approach, the aforementioned relationship is consistently established with the help of MFs. Using a rule base, the different MFs are consistently interrelated incorporating the practical significance of the PoF and CoF in calculating the MII.

Conclusion

The suggested method enables gaps present to be mitigated, for instance inbetween H to VH (i.e. 48 to 6 months), L to M (i.e. 120 to 72 months), etc. In addition, the suggested method enables the experts' knowledge to be recycled. Such recycling provides the opportunity to reduce variability present in the analysis due to human error, inconsistency of awareness, lack of experience, etc. Alternatively, the use of the suggested approach facilitates the improvement of the quality of the MII calculation process in terms of time and accuracy. However, special attention should be paid to developing the MFs.

Future research should be carried out to investigate the effect of triangular vs. Gaussian MFs on the accuracy of the MII analysis process.

Fig. 1 .

 1 Fig. 1. RBI matrix (maximum inspection interval in months).

Fig. 2 .

 2 Fig. 2. An expert system for estimating maximum inspection intervals.

Fig. 3 .Fig. 4 .

 34 Fig. 3. MF for PoF.

Fig. 5 .Fig. 6 .Fig. 7 .Fig. 8 .

 5678 Fig. 5. MF for CoF PPL .

  illustrates approximately 29 rule views and the corresponding calculation of MII (i.e. MII = 4.07) for PoF=1E-1.5 and CoF PLL =1E-3. It is possible to use a similar approach for PoF & CoF Econ vs. MII; TTR & CoF PLL vs. MII; and TRR & CoF Econ vs. MII. All the combinations can be merged into a single workspace using the 'Simulink' facility available in 'Matlab' (MATLAB, 2011).


M Sabar
email: msabar@groupeiscae.com

Multi-agent approach for Personnel Scheduling and Rescheduling in Assembly Centers

Keywords: Personnel shift scheduling, Multi-agent systems, Coalition, Kernel stability, Cross-training, Flexible assembly lines

This article presents a multi-agent based algorithm for personnel scheduling and rescheduling in the dynamic environment of a paced multi-product assembly center. Our purpose is, on the one hand, to elaborate daily assignment of employees to workstations to minimize the operational costs as well as the personnel dissatisfactions and, on the other hand, to generate an alternative planning when the first solution has to be rescheduled due to disturbances related to operators' absenteeism. The proposed approach considers the individual competencies, mobility and preferences of each employee, as well as the personnel and competency requirements associated with each assembly activity given both the current master assembly schedule and the line balancing for each product. To benchmark the performance of the multi-agent approach, we use solutions obtained through a simulated annealing algorithm. Experimental results show that our multi-agent approach can produce high-quality and efficient solutions in a short computational time.

Introduction

Personnel scheduling problems are particular cases of resource allocation problems [START_REF] Hao | A Neural Network Application in Personnel Scheduling[END_REF]). They aim to construct a working timetable for each employee by defining start time periods, duration of work, break intervals, as well as the tasks to be fulfilled. The objective is for the timetable to optimize one or several criteria while respecting a set of constraints such as labor requirements, individual preferences, or specific competencies [START_REF] Thompson | Improved implicit optimal modeling of the labor shift scheduling problem[END_REF][START_REF] Ernst | Staff Scheduling and Rostering, A review of applications, methods and models[END_REF]). Typical classifications of personnel scheduling problems tend to separate problems in three categories [START_REF] Ernst | Staff Scheduling and Rostering, A review of applications, methods and models[END_REF]): shift, days-off and tour scheduling problems.

In this article, the focus is on shift-scheduling problems in a large assembly line environment where the pace setting takt time between individual product units is preset, equal to at least a few minutes.

We consider an assembly line with multiple workstations responsible to sequentially assemble different product-models. Each assembly task requires one or several employees with a specific competency profile. In addition, concerning the manpower pool, we take into account individual competencies, mobility and preferences. Specifically, we consider that:

 Each worker has a specific degree of cross-training, enabling them to carry one or several types of assembly activities during a work shift.

 Workers are allowed to move between workstations in order to fulfill specific assembly activities according to the product assembly schedule.

 Workers can be assigned to secondary activities, which can be either productive, administrative or learning, when the employees are not assigned to an assembly workstation.

 Each worker has a set of individual preferences related to (1) the shift duration, (2) the assignable activities and (3) the number of transfers between activities.

Due to product changeovers and the specific manpower competency requirements associated with each product at each station, there are often large waves of personnel moves among stations. This causes significant disruptions to operations, deterring the overall productivity of the line and causing dissatisfactions among the personnel [START_REF] Sabar | Competency and preference based personnel scheduling in large assembly lines[END_REF][START_REF] Sabar | An agent-based algorithm for personnel shift-scheduling and rescheduling in flexible assembly lines[END_REF]). Also, this paper considers that if a disturbance related to operators' absenteeism occurs at a given time period, the variables representing the personnel' scheduling up to this period are fixed to the matching values from the original plan and, the disturbance' parameters are considered for the remaining time periods. Then, a new rescheduling has to be performed on the global level through complete regeneration of personnel schedule. Concerning the operators' absence, we distinguish full absenteeism (operator absents from work during the entire shift) from partial absenteeism (operator arriving late to work or leaving work during working hours due to sickness or personal affairs). In both cases, the rescheduling aims to generate a new allocation plan which replaces absent operators by transferring activities to the available and on-call workers.

In this context, we present a multi-agent based approach that aims to tackle the complexity of our targeted personnel scheduling / rescheduling problem through distributed problem-solving. The proposed approach is based on cooperation among several rational agents which encapsulate individual competencies and preferences of employees. In this approach, the agents negotiate to form coalitions which allow them to improve their individual schedules, and consequently to iteratively improve the global solution of personnel scheduling problem. To benchmark the performance of the multi-agent approach, we use solutions obtained through a simulated annealing approach.

The remainder of the article is organized as follows. Section 2 defines multi-agent settings and presents in detail the multi-agent architecture and principles retained to sustain our scheduling/rescheduling approach. Section 3 presents the formal description of the proposed multi-agent approach. The experimental setup and results are discussed in section 4. Finally, section 5 presents conclusive remarks.

Multi-agent architecture proposed for personnel scheduling and rescheduling problem

In this research, the real environment to model is an assembly line system. It is made up of several workstations which constantly require a mix of skilled employees. In order to model such a system, we developed a multi-agent system composed of heterogeneous and autonomous agents, which cooperate with one another to produce a personnel schedule. Each agent represents a physical entity of the assembly system, or encapsulates a planning and decision making function. Our multi-agent system includes four categories of autonomous agents: a production-agent; a station-agent for each workstation; a coordinating agent; and an employee-agent representing each employee. These agents are autonomous, rational and able to communicate with each other.

(1) Production-agent: elaborates and manages the production planning. It determines the dynamic sequence and quantity of the product models to be assembled. The production plan is then commu-nicated to stations-agents. The Production agent uses a set of priority rules to decide which job orders are to be planned. Its objective is to optimize criteria such as the maximization of the workstations' utilization or the minimization of delays. In this article, we consider that the production planning decisions are independent of the influence of human resources management.

(2) Station-agent: manages and controls the assembly activities of a workstation. Based on the production planning, this agent defines the needs of the workstation concerning number of employees and their required competencies, which it sends to the coordinator-agent.

(3) Coordinator-agent: is responsible for coordinating the employee-agents. First, it elaborates an initial solution of personnel scheduling. Then, it takes the active role of mediator in the negotiation process among the employee-agents who will try to improve their initial work plans through activities swapping.

(4) Employee-agent: represents the individual interests of an employee. It encapsulates the state as well as the main characteristics of the matching employee, in particular his competencies, his preferences and his allocation history. These agents can negotiate and cooperate among them in order to maximize their profit and their satisfaction. In the proposed architecture, they are coordinated by the coordinator-agent which plays a mediator's role.

In our approach, the personnel scheduling process is supported by the Coordinator-agent and the Employee-agents. In the first step, the coordinator-agent uses priority rules to produce an initial solution taking into account employees 'needs and competencies. Next, the initial work plan of each employee is transmitted to the corresponding Employee-agent. Considering that the initial solution is often mediocre quality in terms of total cost and employees' satisfaction, we use coalitions to improve the performance of the schedule. The issue of coalition formation has been studied in the game theory literature in the context of cooperative N-person games [START_REF] Rapoport | N-person game theory: concepts and applications[END_REF].

In our approach, coalitions are formed among Employee-agents who negotiate a potential mutual agreement on what activities to swap in order to increase their individual profits and ultimately improve the global personnel scheduling solution. Each Employee-agent is rational and self-interested. It has interest in forming coalitions which releases him from less satisfying assembly activities and/or allows him to get a more satisfying set of activities. The proposed coalition approach is round based.

In each round at most one coalition is formed. Each round involves two phases:

-Phase 1 aims to generate a stable coalitional configuration, which consists of a partition of the set of employee-agents into disjoint stable coalitions. To maintain polynomial complexity of the formation process, we restrict our research on coalitions of size two. At the beginning of each round, each employee-agent contacts other employee-agents with whom it has common competencies. It seeks to identify whether there are assembly activities that can be swapped and sends coalition offers. In the scheduling case, all activities' exchange possibilities are tested. However, in the rescheduling case, only activities which are not carried out yet are taken into consideration. A coalition offer contains the list of tasks that may be switched, and the corresponding payoff. For each received coalition-offer by an employee-agent, it uses the Kernel stability concept [START_REF] Davis | The Kernel of a cooperative game[END_REF] to test the coalition equilibrium. The Kernel is based on the idea that the members of a coalition should be in an equilibrium concerning their power to object to each other's payoff. If an employee-agent finds that it can outweigh the other according to the initial payoff, it uses the Transfer Scheme pro-posed in Streans (1968) to demand a side payments transfer. Using a stable payoff distribution, the employee-agents can compare different coalition structures. In fact, each employee-agent compares the payoff of all received proposals. It chooses the coalition which is most beneficial for the employee it represents (i.e. the one that maximizes his utility function), and informs the Coordinator-agent about the accepted coalitions.

-Phase 2 proceeds to select the coalition to enact. Once all accepted coalitions have been received by the Coordinator-agent, it randomly selects a coalition among the group of coalitions that have a bilateral acceptance of the two members. In a rescheduling case, the priority is given to a coalition among those which disengage completely or partially an agent-employee of his activities during his absence periods. Next, the Coordinator-agent informs the two coalition's members about the agreement. These two agents complete the process by exchanging tasks. Based on the new task distribution, the employee-agents start a new round of coalition formation.

These two phases are structured within the framework of an anytime algorithm. Such a type of algorithm improves gradually the quality of its solution as computation time increases and can be interrupted at any time during computation to provide a solution [START_REF] Russell | Composing Real-Time Systems[END_REF].

Multi-agent approach for personnel scheduling and rescheduling

In this section, we formally present our multi-agent approach for personnel scheduling. We describe on one hand the algorithm for initial solution generation, and on the other hand, sequentially the algorithm to generate coalitions with K-stable payoff distributions and the algorithm for coalition selection. It should be noted that whenever the process of solution improvement by coalition formation is stopped, it produces the best currently available solution. This process is generally interrupted when the desired state is reached or when having to answer an immediate need of assembly lines about some employee's allocation, for example when an employee becomes out-of-kilter.

The initial solution of employee allocation is performed at the coordinator-agent level based on a priority dispatching rule. Using the production planning, each station-agent has a view of local requirements concerning the number of employees and their competency profile. The used dispatching rule involves the selection, period by period, of the workstation with the least extra number of employees that have the required competency profile. At a given period, the extra number of employees is equal to the difference between the number of available employees and the required number of employees. For each selected workstation, the coordinator-agent assigns the least cross-trained employee available among those who have the required competency profile.

At the end of the first stage, each employee-agent , is the activity k to execute on station l in the period t by the employee-agent i . To evaluate the utility of each employee-agent

  n N i ,..., 1  
  n N i ,..., 1  
, we use a linear function:

) ( ) ( i i i i f S V    
, where S is a constant which corresponds to an initial amount allocated to each employee-agent, it represents an artificial gain that each employee earns if it suc-ceeds to totally release himself from duty. i f is an increasing linear function of work duration and dissatisfaction of employee i . In fact, 8 7

6 5 4 3 2 1 F F F F F F F F f i        
, where 1, :   8 F : Penalty cost (positive or negative) associated to the dissatisfaction or satisfaction of employee i for its assignment to a set of activities.

The algorithm to generate coalitions with K-stable payoff distributions (Phase 1) is as follows:

Each agent-employee i :

1. Maintains a register concerning the references of the employee-agents with whom it can exchange certain activities. 2. For each employee-agent j indexed on its register, i tests all possible permutations of activities and calculates the value of each potential coalition

) ( ) ( ) ( j new j i new i ij V V C v    
. In case of several permutation possibilities with an agent j , i retains the one which generates the highest coalition value. In rescheduling case, the coalitions are based on the exchange of the activities which are not carried out yet (i.e. activities planned between the periods of absence' notification and the shift-end). Activities carried out are considered fixed and cannot be changed. However, they are taken into account in assessing the total cost of personnel staff scheduling.

If

) ( ) ( ) ( j j i i ij V V C v     , then
i sends a coalition proposal ij PR to j . The proposal encapsulates the set of activities to be permuted; the coalition value

) ( ij C v and the initial proposed payoff 2 )) ( ) ( ( ) ( ) ( j j i i ij j j j V V C v V u        (i.

The utility function i

V is designed so as to generate more profit for an employee who succeeds to release himself from duty or acquire a set of activities which creates a higher satisfaction.

Given a pair of employee-agents ( i , j ) with the activity vectors i  and j  , we define the potential value of the coalition j i C , as:

) ( ) ( ) ( , j new j i new i j i V V C v    
, where i new  and j new  are the new activity vectors of i and j if they agree to form the coalition by permuting a part of their initial activities. To accept a coalition, the payoff of each agent after the redistribution of the coalition value must be at least equal to its initial self-value, i.e.

; ) (

j i ij u u C v   ) ( i i i V u   and ) ( j j j V u  
. Each employeeagent uses the Kernel concept to evaluate the offered payoff and to assess its power to object to its partner's payoff. A general strategy used by employee-agents for coalition formation and payoff distribution is defined as follows. 3.1. Initialization of the regression coefficient :

1   3.2. Each employee agent   n i ,..., 1  : -Elaborates the list  , i
 of K-stable coalitions that give him a payoff at least equal to

max , i u   ;
-Sends the list to the coordinator-agent.

3.3. Based on all the received lists, the coordinator-agent selects the set of coalitions which have a bilaterally acceptance from the two members i.e. the coalitions

); ( ) ( |   j ij i ij ij C C C        n j i ,..., 1 ,  
. At this level, there are two possible scenarios :

   BC
: several coalitions have a bilaterally acceptance of their two members:

-In scheduling cases, the Coordinator-agent randomly selects a coalition from BC . In rescheduling cases, it selects in priority a coalition among those which release completely or partially an agent-employee of duties during his absence periods. Then, informs the two coalition's members about the agreement. -These two agents finalize the process by exchanging tasks. Based on the new tasks distribution, the employee-agents start a new round of coalition formation (return to Stage 2).

   BC

: no consensus is reached, then the regression coefficient will be decreased

     : -If 0   return to 3.2.
-If 0   the global solution has reached a local optimum (i.e. given the current activities distribution, employee-agents have no benefit by forming coalitions), then we introduce an artifice for fictitious payoffs distribution. This artifice randomly generates and attributes factitious profits to a certain number of employee-agents in such manners as to incite them to form coalitions. Return to stage 2 in order to generate new K-stable coalitions. At the end of this stage, we obtain a set of potential coalitions with stable payoff distributions. Each employee-agent may have several offers of coalitions with various profits. Since each employeeagent i is rational, it tries to form the coalition, among all possibilities, in which it earns the greatest payoff max , i u . However, if an agent i chooses to form a coalition with the agent j , nothing guaran- tees that agent j will accept because j may earn more by forming another coalition with a third agent k . In case of conflicts of interest between employees-agents, we introduce a regression func- tion reg f which allows agents to reduce the value of their aimed payoff in order to reach a consensus. . After that, it communicates the results to the coordinator-agent which randomly selects a coalition among the group of coalitions that have a bilaterally acceptance of the two members. The detailed procedure for coalition's selection is defined as detailed on Figure 1.

BC

Computational Experiments

In this section, we present experimental results concerning several shift scheduling problems in the context of a paced multi-product assembly center. A set of five problems is conducted to test the performance of the proposed multi-agent approach. For each problem, the production planning is spread out over 60* 9-hour shift. Experiments have been performed in an assembly line consisting of 40 workstations. The takt time between two product units is preset equal to 15 minutes. For each workstation, the employee requirements in a given period are determined according to the assembly activities to be fulfilled on the scheduled product according to the preset line balancing. Concerning the staff, we consider that the offer and the demand per shift for employees vary between 150 and 200 employees. The daily absenteeism rate varies according to the shift number, from a minimum of 1% to a maximum of 5% of total employees.

To evaluate the quality and the efficiency of the proposed multi-agent approach, we report solutions obtained on the same problems through the simulated annealing approach (SA). A detailed description of the proposed simulated annealing algorithm and its use for resolution of personnel scheduling problems can be found in [START_REF] Sabar | A multi-agent based approach for personnel scheduling in assembly centers[END_REF][START_REF] Sabar | An agent-based algorithm for personnel shift-scheduling and rescheduling in flexible assembly lines[END_REF]. For each shift s, we report the cumulated deviation CD s between the best solutions founded by these two approaches for a computation time equal to 10 minutes for scheduling and 3 minutes if rescheduling is required. Figure 2 exhibits the evolution of this deviation between SMA and SA results. It shows clearly that for the five test problems the proposed SMA approach leads to high quality solutions in comparison with the SA approach. It is interesting to observe that SMA systematically outperforms SA for all shift results. Indeed, we notice that the deviations of the SMA approach solutions from SA range between -4.2 % and -0.7 %. These results demonstrate that the proposed multi-agent approach for personnel scheduling is effective and can generate high-quality solutions fast and reliably.

Conclusion

In this article, we developed a multi-agent approach for the personnel scheduling/rescheduling problem in the context of a paced multi-product assembly center. The proposed approach is based on cooperation among several rational agents which encapsulate individual competencies and preferences of workers. The experiments we have performed demonstrate that the multi-agent approach can produce high-quality and efficient solutions in comparison with simulated annealing approach.

Our future research will focus on the impact of dynamic random events such as product quality issues on the line, and probabilistic operation times potentially depending on the operator's skill level.

In addition, we will investigate the impact of modeling employee preferences on the quality of the scheduling solutions obtained.

1 F: 2 F 3 F 4 F 5 F : Transfer cost of employee i ;  6 F;  7 F

 1234567 Salary cost of employee i ; : Activity assignment cost of employee i;  : Idleness penalty cost of employee;  : Cost savings generated by the assignment of employee i to secondary activities;  : Penalty cost associated to the deviation from the number of transfers preferred by employee i : Penalty cost associated to the deviation from the total work duration preferred by employee i ;

1 For

 1 e. Dividing the profit generated by the coalition into two equal parts). 4. Receive coalition proposals from the other employees-agents. 5. Evaluate the received coalition proposals: (a) Use the Kernel concept to test the proposals coalitions equilibrium; (b) If employee-agenti dominates any other agent, it uses the Streans' transfer scheme to evalu- ate the side-payment demand and informs the concerned agent. 6. For each instable coalition, send or receive a part of payoff equal to the side-payment demand. more details concerning the function i f , see[START_REF] Sabar | Competency and preference based personnel scheduling in large assembly lines[END_REF]

Fig. 1 .

 1 Fig.1. Algorithm for coalition selection (Phase 2)

Fig. 2 .

 2 Fig. 2. Deviation between MAS and SA results

  For an employee-agent i , this function is defined as

								f	reg	( i	, 	)	:	u	i	,	max			u	i	,	max	, where			  1 , 0	rep-
	resents the rate of payoff's decreasing. Considering the reduced payoff, each employee-agent
	i		  n ,..., 1	chooses among its K-stable coalitions those which give him a payoff at least equal to
			u
					i	,	max