Minh-Quyen Le

Minh Tu Pham

Mahdi Tavakoli

Richard Moreau

Richard Moreau An

Enhanced Sliding-Mode

An Enhanced Sliding-Mode Control for a Pneumatic-Actuated Teleoperation System

Keywords: Pneumatic actuator, on/off solenoid valve, sliding mode control, haptic teleoperation, stability

published or not. The documents may come

I. INTRODUCTION

Teleoperation, which involves operating one or several robots remotely, has been an active research area over the past decade [START_REF] Hokayem | Bilateral teleoperation: An historical survey[END_REF]. A teleoperation system typically involves a master robot handled by a human operator, and a slave robot operating at location far from the master robot.

The basic control requirements on a bilateral teleoperation system are related to its stability and transparency. Transparency corresponds to the quality of reproduction of the remote environment mechanical properties for the operator. Transparency in conjunction with stability can ensure successful completion of tasks in the teleoperation mode.

Traditionally, electrical direct-current (DC) motors are used as actuators in teleoperation systems as they are simple and capable of providing forces suitable for small haptic devices. In DC motor actuated haptic devices that require reflecting high output forces to the operator, the use of gears becomes necessary. This, however, will result in backlash, high inertia and discontinuity in output forces (cogging), which are undesirable as they distort the reflected forces. In addition, electrical motors are not compatible to applications that have extensive magnetic interference, (e.g., applications involving Magnetic Resonance Imaging; MRI) because they produce their own magnetic fields and contain ferromagnetic materials.

On the other hand, pneumatic actuators whose output force is a function of the compressed air pressure lead to M.Q. Le, M.T. Pham, and R. Moreau are with the Laboratoire Ampère, UMR CNRS 5005, Université de Lyon, INSA-Lyon, F-69621 Villeurbanne Cedex, France (corresponding author's e-mail: minh-quyen.le@insalyon.fr). M. Tavakoli is with the Department of Electrical and Computer Engineering, University of Alberta, Edmonton, AB T6G2V4, Canada.

higher force-to-weigh ratios than DC motors. Moreover, pneumatic actuators are inert to magnetic fields, which is crucial in applications such as robot-assisted surgery under MRI guidance [START_REF] Yu | Comparison of MRI-Compatible Mechatronic Systems with Hydrodynamic and Pneumatic Actuation[END_REF]. Consequently, pneumatic actuators have found use in various applications of haptic teleoperation in recent years [START_REF] Durbha | Passive bilateral tele-operation and human power amplification with pneumatic actuators[END_REF][START_REF] Le | Development of a hybrid control for a pneumatic teleoperation system using on/off solenoid valves[END_REF][START_REF] Tadano | Development of a master slave system with force sensing using pneumatic servo system for laparoscopic surgery[END_REF].

In this study, we investigate the development and control of a certain class of electro-pneumatic actuators in a teleoperation system. Generally, servo-valves, rather than solenoid (on/off) valves, are used to achieve high performances in pneumatic control. However, they are typically expensive due to the requirements of high-precision manufacturing. Therefore, in this paper, fast-switching on/off valves are used due to their advantages in terms of low cost and small size. One of the objectives of this paper is to show that good teleoperation transparency can be obtained with these inexpensive components as actuators of the teleoperation system.

The traditional control approach for pneumatic systems with solenoid valves involves using Pulse Width Modulation (PWM) [START_REF] Shih | Position control of a pneumatic rodless cylinder using fuzzy PWM control method[END_REF], [START_REF] Taghizadeh | Improving dynamic performances of PWM-driven servo-pneumatic systems via a novel pneumatic circuit[END_REF]. A main disadvantage of the PWM control is the chattering phenomenon caused by the high-frequency switching of the valves even in steady state [START_REF] Le | Force Tracking of Pneumatic Servo Systems using on/off Solenoid Valves based on a Greedy Control Scheme[END_REF].

To overcome the drawbacks of PWM-based control of solenoid valve actuated pneumatic devices, Nguyen et al. [START_REF] Nguyen | Accurate sliding-mode control of pneumatic systems using low-cost solenoid valves[END_REF] proposed a nonlinear sliding mode controller, which controlled the actuator by modeling it as a discrete-input, three-mode system. In [START_REF] Le | Sliding Mode Control of a Pneumatic Haptic Teleoperation System with on/off Solenoid Valves[END_REF], we applied Nguyen's controller to the bilateral control problem of a pneumatic teleoperation system. In this paper, we will present an extension of Nguyen's controller and re-apply it to bilateral teleoperation control -this will result in superior performance compared to [START_REF] Nguyen | Accurate sliding-mode control of pneumatic systems using low-cost solenoid valves[END_REF] and [START_REF] Le | Sliding Mode Control of a Pneumatic Haptic Teleoperation System with on/off Solenoid Valves[END_REF]. Specifically, we first extend the three-mode model used in [START_REF] Nguyen | Accurate sliding-mode control of pneumatic systems using low-cost solenoid valves[END_REF] and [START_REF] Le | Sliding Mode Control of a Pneumatic Haptic Teleoperation System with on/off Solenoid Valves[END_REF] into a five-mode model. The two extra modes introduced in the new proposed scheme will offer more possibilities in terms of valve switching choices and will, therefore, improve the response of the system. Next, this new 5-mode sliding control law is used in a two-channel (2CH) bilateral teleoperation architecture with three different schemes: position-position, force-force, and force-position. These architectures are chosen due to their implementation simplicity and efficiency. Then, we provide analyses regarding the tracking performance, switching activity and closed-loop stability using a Lyapunov candidate function.

For simplicity, the master and the slave are identical and one-degree-of-freedom (DOF) pneumatic manipulators (with on/off actuation) in this study. It should be noted that this paper does not deal with compensating for time delays in a teleoperation system's communication channel.

An Enhanced Sliding-Mode Control for a Pneumatic-Actuated Teleoperation System

M.Q. Le, M.T. Pham, M. Tavakoli, R. Moreau

The paper is organized as follows. The modeling of the pneumatic manipulator composed of a cylinder and four solenoid valves is described in Section II. Section III presents a design of the 5-mode open-loop model of the master and slave manipulators. The sliding controller scheme for a closed-loop teleoperation system is then provided in Section IV. Section V shows experimental results that validate the proposed control laws. Finally, concluding remarks are given in Section VI.

II. MODEL OF THE PNEUMATIC ACTUATOR

As mentioned above, the master and the slave manipulators are identical, thus only one pneumatic robot is presented in this section. A schematic of the 1-DOF pneumatic actuation system is shown in Fig. 1. Assuming air is a perfect gas undergoing isothermal process, the behavior of the pressure inside each chamber of the cylinder can be expressed as

P
γ γ   = -       = +     (1)
where the subscripts p and n respectively denote chamber p and n, U 1 , U 2 , U 3 and U 4 are the discrete control voltages (0 or 1) of valve 1, valve 2, valve 3, and valve 4, y and y are the position (m) and velocity (m/s) of the piston, P p and P n are the pressures (Pa), V p and V n are the volumes of chambers (m 3 ), S p and S n are the piston cylinder area (m 2 ), q p and q n are the mass flow rates (kg/s), T a is the temperature of the supply air (K), r is the perfect gas constant (J/(kg.K)) and γ is the polytropic constant. The mass flow rate characteristics of the on/off valves can be expressed as functions of the discrete control voltages and the pressures: q U U P U q P P U q P P q U U P U q P P U q P P = -= -

1 2 1 2 3 4 4 3 ( , , ) ( , ) ( , ) ( , , ) ( , ) ( , )
where P sa and P atm are the pressures of the supply air and atmosphere. The functions in (2) are given by a standard expression in which the mass flow rate of the valve is regulated by the air passage through an orifice:

C P T T C  ≤   =  -    -  -     (3)
In the above, C val is the valve flow rate coefficient, P up and P down are respectively the absolute upstream and downstream stagnation pressures of the valve, T atm is the atmosphere temperature, and T up is the upstream stagnation temperature.

Finally, the dynamics of the piston and the load are

p p n n st ext

My S P S P by F F

= - --+ ( 4
)
where b is the viscous friction coefficient (N.s/m), M is the moving load (kg), F st is the stiction force, and F ext is the external force (N). For simplicity, the stiction force is assumed to be negligible.

III. 5-MODE OPEN-LOOP MODEL AND ENHANCED SLIDING CONTROL OF THE PNEUMATIC ACTUATOR

A. Controller mode selection

In order to facilitate the control law design, a switching scheme for the four solenoid valves of Fig. 1 is defined so that each of the master and slave robots has the five modes of operation shown in Table I. Note that the three first modes are the three classical ones presented in [START_REF] Nguyen | Accurate sliding-mode control of pneumatic systems using low-cost solenoid valves[END_REF], [START_REF] Le | Sliding Mode Control of a Pneumatic Haptic Teleoperation System with on/off Solenoid Valves[END_REF].

u 1 -1 0 0.5 -0.5 U [1 0 1 0] [0 1 0 1] [0 0 0 0] [1 0 0 0] [0 0 0 1]
In Table I, U (the 4 th row) is the input voltage vector of the four valves, i.e.,

U = [U 1 U 2 U 3 U 4 ]
. The '0' state of each input voltage corresponds to a closed valve and the '1' state corresponds to an open valve. All the states where U 1 = U 2 = 1 and U 3 = U 4 = 1 are prohibited to avoid a bypass of the valves. Also, u (the 3 rd row) is a newly introduced discrete control input that has five levels to match the five modes of operation. This new input can be chosen either as u = k.sign(s) or u = -k.sign (s) where k equals 1 or 0.5, and s is a sliding surface, which is a function of tracking error. Whether to choose the switching control u as k.sign(s) or ask.sign(s) depends on the definition of s and the open-loop system. As it will be shown later, this choice is crucial to ensuring the stability of the teleoperation system. Note that mode 3 in Table I is used to "de-actuate" the piston when the tracking error is small enough. Mode 4 allows for moving the piston in the same direction (to the right in Fig. 1) as in mode 1 but with a slower dynamics (because in mode 4 the chamber n is closed as opposed to exhausting as in mode 1; in both cases the chamber p will be filling). Mode 4 may be considered as an intermediate (or average) mode between modes 1 and 3, due to which the control vector u is chosen as 0.5 (Table I). On the other hand, mode 5 whose control vector u equals -0.5 is used to move the piston in the other direction and could be regarded as an intermediate level of actuation (and piston acceleration) between modes 2 and 3.

B. Open-loop models of master and slave

Ignoring the stiction force in (4), the dynamics of the master and slave robots can be written as

= - - + = - - - (5)
where f h and f e are the operator force exerted on the master and the environment force exerted on the slave, and y m and y s are the master and slave positions. Differentiating (5) and using ( 1)-( 2), the dynamics of the master and slave manipulators are obtained after some manipulations as [START_REF] Nguyen | Accurate sliding-mode control of pneumatic systems using low-cost solenoid valves[END_REF]:

, 0 , 0 , 0 , , 0 , 0 , 0 m m h m s s e s m m m h m s s s e s m h m s e s f M u f M u y f M u y f M u f M u f M u α β α β α β α β α α + + - -   + + > + - >     = - + < = - - <     + = - =     (6)
where u m and u s denote the discrete control input as defined in Table I. In the above,

p p i n n i i i i p i n i S P S P b y y M M V y V y γ α   = - - +       2 2 , , ( ) ( )
When u i is positive, it might take two values, i.e., 1 or 0.5, which corresponds to the modes 1 or 4. Since i β + depends on the control choice of the valves, it can be written as

p sa p i n i atm n a i p i n i i p sa p i a i p i S q P P q P P S rT u M V y M V y S q P P rT u M V y γ β γ +    + =         =   =   (8) , , , ( , ) ( , ) , 1 ( ) ( ) ( , ) , 0.5 ( )
Inversely, when u i is negative, it has two possible values, i.e., -1 or -0.5, which relates to the control modes 2 or 5. In this case, i βis calculated as

, , , ( , ) ( , ) , 1 ( ) ( ) ( , ) , 0.5 ( )
p p i atm sa n i n a i p i n i i sa n i n a i n i S q P P q P P S rT u M V y M V y q P P S rT u M V y γ β γ -    + = -         =   = -   (9)
with i = m or s (for master or slave, respectively).

The operation of the 5-mode sliding controller is based on the following three principles: 1. In order to bring the system to the sliding surface s = 0, which corresponds to perfect tracking performance, at steady state we define a neighbourhood of radius ε << 1 around zero. When |s| is within the interval [0, ε], mode 3 (u = 0) is used to conserve energy and somewhat reduce chattering. 2. To be able to switch between modes 1 and 4 or modes 2 and 5, a threshold ε 1 is introduced where ε 1 > ε. When |s| is within the interval [ε, ε 1 ], either mode 4 or mode 5 (u = ± 0.5) is used to provide slower dynamics compared to modes 1 or 2 (u = ± 1), respectively. Still, the piston is actuated to move in the direction that minimizes s. 3. When |s| is within the interval [ε 1 , ∞], either mode 1 or mode 2 (u = ± 1) is used to provide fast dynamics, highly accelerating the piston in the direction that minimizes s. In summary, we get our 5-mode sliding control law as

1 1 sign( ) if 0.5sign( ) if < 0 if s s u s s s ε ε ε ε ± ≥  = ± <   ≤  (10)

IV. CLOSED-LOOP TELEOPERATION SYSTEM

In this section, the five-mode sliding control is applied to a two-channel bilateral pneumatic teleoperation system with various different architectures (i.e., position-position, forceforce, and force-position). In the following, we provide analyses including tracking performance and closed-loop stability using a Lyapunov candidate function.

A. Position error based (PEB) control

A position-error-based, also called position-position, teleoperation system involves the simplest bilateral controller in which no force sensors are required. This architecture involves the transmission of two types of data between the master and the slave: position (or velocity) from the master to the slave and vice versa. The pneumatic-actuated PEB teleoperation system with our proposed sliding mode control is shown in Fig. 2. The sliding surface s = 0 is reached within a finite time if the following condition is satisfied:

V ss s η = < - (13)
for some constant η > 0. Thus, from [START_REF] Tavakoli | Enhanced Transparency in Haptics-Based Master-Slave Systems[END_REF]

m s λ β β η - + + > ( 19
)
where λ is defined in (18).

Note that, from ( 8) and ( 9), and

i i β β +
-are positive, and can be made as large as desired by choosing a sufficiently large valve orifice C val in (3). Thus, to ensure that the conditions ( 17) and ( 19) are satisfied, we only need to show that λ is bounded. For the demonstration of the boundedness of λ, the readers can refer to [START_REF] Le | Sliding Mode Control of a Pneumatic Haptic Teleoperation System with on/off Solenoid Valves[END_REF]. When λ is bounded, the sliding condition in (13) is ensured, which implies that the position tracking error tends to zero (and that the overall system is stable).

Note that a drawback of the PEB method is that it does not guarantee a good transparency in term of force tracking. In order to improve the tracking performance, other schemes are proposed in the next subsections.

B. Force error based (FEB) control

A force-error-based, also called force-force, system is not commonly used in two-channel bilateral teleoperation since two force sensors are required and since position tracking is not good. However, compared to the PEB architecture, this architecture can improve the force tracking performance. The implementation of the proposed sliding mode control on the FEB system was schematically presented in [START_REF] Le | Sliding Mode Control of a Pneumatic Haptic Teleoperation System with on/off Solenoid Valves[END_REF].

Consider the control law u s = u m = k.sign(s) where the sliding surface is defined as

h e s f f = - (20)
Using the Lyapunov function [START_REF] Lau | Implementation of position-force and position-position teleoperator controllers with cable-driven mechanisms[END_REF], we need to show that the sliding condition (13) is satisfied. From the master and models (6), we can calculate h f and e f as

( ), 0 ( ), 0
( ), 0 , ( ), 0
( ) , 0 ( ) , 0 m m m m s m s s h m m m m e s m s s m m m s m s M y u M y u f M y u f M y u M y u M y u α β α β α β α β α α + + - -   - - > -+ + >   = - + < = -+ - <     - = -+ =   (21)
Case 1: s > 0. In this case u m and u s are positive. From (13), we need ( )

h e m m m s s s s f f M y y α β α β η + + = -= - - + -- < - (22)
It was shown in [START_REF] Le | Sliding Mode Control of a Pneumatic Haptic Teleoperation System with on/off Solenoid Valves[END_REF] that m y , s y , m α and s α are bounded. Thus, by choosing a valve with large orifice, m β + and s β + can be made sufficiently large to satisfy (22).

Case 2: s < 0. In this case, u m and u s are negative. We have ( )

h e m m m s s s s f f M y y α β α β η - - = -= - + + -+ >
(23) Similar to Case 1, the stability of the system can be guaranteed by choosing a large enough value of m βand s β -.

Finally, the force tracking error converges to zero and the overall system is stable. However, the FEB method does not guarantee a good position tracking performance. In order to overcome the PEB and FEB architecture drawbacks, we use the DFR scheme described in the following subsection.

C. Direct force reflection (DFR) control

A direct-force-reflection, also called force-position, system has advantages over the position-position and forceforce architectures. Compared to the PEB method, improvements in term of force tracking is achieved due to the measurement of the interaction force between the slave and the environment. Furthermore, its position tracking performance is better than the FEB case thanks to position information. For the implementation of the proposed sliding mode control on the DFR system, please refer to [START_REF] Le | Sliding Mode Control of a Pneumatic Haptic Teleoperation System with on/off Solenoid Valves[END_REF].

In this section, we use a Lyapunov function to prove the stability of the sliding-mode controlled DFR system. First, we will show the stability of the force-controlled master manipulator. Afterwards, we will show the stability of the position-controlled slave manipulator. However, the stability of the overall system is difficult to show due to the complexities introduced by using different sliding surfaces for the master and for the slave.

1) Force convergence of the closed-loop master system

The sliding surface s m and the Lyapunov function V m are defined as in (20) and ( 12), respectively. The controller u m is chosen to be similar to the FEB system in subsection IV.B.

To ensure the sliding condition (13), we need

e m m m f My M M α β η + -+ - - < - (25)
Similar to how it was demonstrated in [START_REF] Le | Sliding Mode Control of a Pneumatic Haptic Teleoperation System with on/off Solenoid Valves[END_REF], e f , m y and m α can be shown to be bounded. Thus, the stability condition (25) is satisfied by choosing a large enough value for m β + .

Case 2: s m < 0. In this case, u m < 0. From ( 21) and (20) we need

e m m m f My M M α β η - -+ - + >
(26) Similar to Case 1, it is possible to choose a large enough value of m βto ensure the stability of the master robot.

Consequently, the sliding surface (the force tracking error) tends to zero, i.e., f h tends towards f e .

2) Position convergence of the closed-loop slave system

The sliding surface s s and the Lyapunov function V s are defined as in [START_REF] Tavakoli | Enhanced Transparency in Haptics-Based Master-Slave Systems[END_REF] and ( 12) respectively. The controllers u s is chosen to be similar to the PEB system in subsection IV.A. Case 1: s s > 0. In this case, u s < 0. The sliding condition (13) is equivalent to

Using the expression of s y in (16) leads to (29) The straightforward reasoning described in [START_REF] Le | Sliding Mode Control of a Pneumatic Haptic Teleoperation System with on/off Solenoid Valves[END_REF] allows us to infer that φ is bounded. Therefore, there exists a high value of s βsuch as (28) is satisfied.

s ϕ β η - - < - (28)
Case 2: s s < 0. In this case, u s > 0. Thus, we need

s ϕ β η + + > ( 30
)
where φ is defined in (29) This condition is achieved by choosing a large enough s β + .

Note that for both cases, the convergence of the sliding surface to zero is proved, so x s tends towards x m .

V. EXPERIMENTAL

A. Experimental setup

In this section, experiments with a 1-DOF teleoperation system are reported. As illustrated in Fig. 3, the prototype consists of two identical master and the slave pneumatic manipulators. For more detail about the experimental description, please refer to [START_REF] Le | Development of a hybrid control for a pneumatic teleoperation system using on/off solenoid valves[END_REF]. The controller is implemented using a dSPACE board (DS1104), running at a sampling rate of 500 Hz. This is because, in practice, it is easier if the sliding surface for position control is chosen to have the same dimension as positions as in (31) and not accelerations as in [START_REF] Tavakoli | Enhanced Transparency in Haptics-Based Master-Slave Systems[END_REF]. The first derivative of the position error in (31) is computed through a backward difference method applied on the position signal followed by a second-order Butterworth filter with a cutoff frequency of 70Hz. The second derivative is computed in the same way from the filtered first-derivative signal.

For a good trade-off between the position tracking performance and the chattering problem, the parameters ξ = 0.5 and ω = 70 rad/s are chosen in practice.

In the experiments, the position and force thresholds ε p and ε f are respectively chosen equal to 0.5 mm and 0.1 N in order to achieve good tracking responses without causing too much switching of the valves.

B. Experimental results

Fig. 4 shows the master and the slave force and position tracking profiles in free space and in contact motion for the PEB teleoperation system. Fig. 5 and Fig. 6 illustrate the same profiles for the FEB and the DFR systems, respectively.

As illustrated in Fig. 4, the PEB system provides a good position tracking responses. However, the force response is not as good because no force sensor is used. On the other hand, the force tracking performance of the FEB system is much better, thanks to the knowledge of the force information. Nonetheless, the transparency of the position tracking deteriorates in FEB control -as it can be seen in Fig. 6, the slave's movement does not accurately track the master's movement. Interestingly, the DFR system in Fig. 6 provides an improvement in terms of position tracking response compared to the FEB system. It also displays a superior force tracking performance compared to the PEB, especially under contact motion where good quality force feedback is required. Thus, DFR control provides superior performance compared to both PEB and FEB control. Since the DFR system uses the measurement of slave/environment contact forces, the feeling of contact motion is highly realistic in our experiments. Therefore, the performance of the teleoperation system is improved significantly by feeding the operator with the slave/environment contact force. This result agrees with the previous theoretical work.

Among the three architectures, the DFR scheme seems to be a better choice to obtain a good transparency. Although the various teleoperation controllers have previously been compared from a performance perspective in the literature [START_REF] Tavakoli | Enhanced Transparency in Haptics-Based Master-Slave Systems[END_REF], [START_REF] Lau | Implementation of position-force and position-position teleoperator controllers with cable-driven mechanisms[END_REF], this is the first study to show that it is possible to achieve stability and satisfactory performance using manipulators actuated by low-cost switching on/off valves.

To show the benefits of the 5-mode control system teleoperation, we need to take the transparency analysis beyond only studying the force and position responses. Thus, a frequency analysis of the control signals (controller outputs) is carried out. In the following, due to the space restrictions, we only compare the spectra of the control signal U 1 (instead of U 1 , U 2 , U 3 and U 4 ) of the master side. Similar results can be observed for the other control signals. The frequency responses can be found by using the discrete fast Fourier transform fft command in Matlab. Fig. 7 shows the spectral analysis of the control signal U 1 obtained with the 5mode scheme and the 3-mode scheme of [START_REF] Le | Sliding Mode Control of a Pneumatic Haptic Teleoperation System with on/off Solenoid Valves[END_REF]. As it can be observed, for all three architectures (PEB, FEB, and DFR), the magnitude |U 1 | of the spectra is lower in the 5-mode case than in the 3-mode case over almost all frequencies, especially so at high frequencies. This shows that by using two additional modes, the on/off activity of the valves has been reduced, which allows to provide better tracking responses with smoother dynamics and less oscillations, particularly over high frequencies.

VI. CONCLUSION

In this paper, pneumatic actuators with inexpensive solenoid valves are chosen for the development of a masterslave teleoperation system. To improve the dynamic performance and reduce the switching activities of the valves, a five-mode sliding control scheme has been used. Our study demonstrated that by increasing the number of the possible control actions for the valves, we can reduce the valves' switching activities, hence improving the valve's life times at no cost to teleoperation transparency. In order to evaluate the efficacy of the sliding mode approach, a comparison of the transparency and stability has been performed between three control architectures (PEB, FEB, and DFR) in a 2-channel bilateral teleoperation framework. The drawbacks of the PEB and FEB schemes in terms of less-than-ideal force or position tracking performance are analytically and experimentally demonstrated. Also, the DFR control scheme is shown to be highly transparent thanks to its used of force and position measurements.

The experiment results are encouraging for future work aimed at implementing the proposed sliding bilateral control on a direct-drive, multi-DOF, pneumatic-actuated teleoperation system. Another aspect of this work is to incorporate sliding mode control in more complex teleoperation architectures (e.g., three or four-channel methods) to improve the teleoperation transparency.

Fig. 1 .

 1 Fig. 1. Electro-pneumatic system with four valves

Fig. 2 .

 2 Fig. 2. Position-error-based approach with sliding mode controlIn the PEB scheme, the sliding surface s is defined as 2 2 s e e e ξω ω = + +

Case 1 :

 1 s m > 0. In this case, u m > 0. Using the expression of h f in (21) and the definition of s m as in (20) we have

Fig. 3 .

 3 Fig. 3. Pneumatic master-slave teleoperation experimental setup For the PEB and DFR systems, the sliding surface of the position-controlled slave is normalized as 2 2 e s e e ξ ω ω = + +

Fig. 4 .Fig. 5 .Fig. 6 .

 456 Fig. 4. Position and force profiles for the PEB teleoperation system

Fig. 7 .

 7 Fig.7. Discrete Fourier transforms of the control signal U1 for the master manipulator: 3-mode scheme (Solid), 5-mode scheme (Dotted)

TABLE I :

 I FIVE POSSIBLE CONTROL MODES

		Mode 1	Mode 2	Mode 3	Mode 4	Mode 5
	Ch. p	fills	exhausts	closed	fills	closed
	Ch. n	exhausts	fills	closed	closed	fills


Sean Hodgson

Minh Quyen Le

Mahdi Tavakoli

Minh Tu Pham

Sliding-Mode Control of Nonlinear Discrete-Input Pneumatic Actuators

Keywords: Pneumatic actuator, on/off solenoid valve, robot control, sliding mode control, position tracking, switching activity, stability

This paper proposes a sliding mode law for precise position control with minimal switching activity for a robotic system that uses on/off (solenoid) pneumatic actuators, For a two-chamber pneumatic actuator with four binary solenoid valves, there is a total of sixteen possible input combinations defined directly from the state of the four on/off solenoid valves present in the system. However, only seven of these discrete operating modes are considered both functional and unique. Accordingly, we use a seven-mode sliding controller that minimizes the position error using modes that have both the necessary and sufficient amounts of drive energy and, thus, involve reduced switching activity. An analysis of the closedloop system stability is carried out. The performance of the proposed control design is experimentally verified on a single pneumatic actuator comprising of two chambers driven by four on/off solenoid valves.

I. INTRODUCTION

Recent developments in robotics and telerobotics have allowed working in environments that are normally difficult or hazardous for the human hand to reach or operate in. Robotic systems have been developed for applications ranging from surgery to space exploration. In robot-assisted surgery, for instance, robots may be required to work in areas that have strong magnetic fields. Magnetic Resonance Imaging (MRI) affords images with high resolution and contrast and allows surgeons to access a patient's three-dimensional visualization in real-time [START_REF] Oura | Development of MRI compatible versatile manipulator for minimally invasive surgery[END_REF].

Pneumatic actuators offer many advantages for positioning applications in addition to being inert to magnetic fields: low maintenance cost, high ratio of power to weight, cleanliness and safety are some of these advantages [START_REF] Rosas-Flores | Optimal linearization of the dynamic behavior of an on/off actuated single pneumatic cylinder[END_REF]. However, they suffer from common drawbacks including friction and sensitivity of actuator dynamics to load and piston position along the cylinder stroke [START_REF] Girin | High-order slidingmode controllers of an electropneumatic actuator: Application to an aeronautic benchmark[END_REF]. Also, from a control perspective, controlling a pneumatic actuator is a challenge because the system dynamics are highly non-linear [START_REF] Xing | Tracking control of pneumatic artificial muscle actuators based on sliding mode and non-linear disturbance observer[END_REF].

The non-linear nature of a pneumatics actuator is exacerbated when it uses on/off solenoid valves, which have wide-spread applications due to the high cost of servovalves (due to the precision machining required in their manufacturing) [START_REF] Nguyen | Accurate slidingmode control of pneumatic systems using low-cost solenoid valves[END_REF]. In this case, precise control is difficult due to the discrete-input nature of the system. Most solenoidvalve pneumatic systems utilize a pulse width modulated (PWM) controller. Using time averaging, a PWM input with a sufficiently high frequency can approximate the continuous input properties of a servo-valve [START_REF] Shen | Nonlinear modelbased control of pulse width modulated pneumatic servo systems[END_REF]. Regardless of whether Fig. 1. Electro-pneumatic system with 4 on/off solenoid valves a servo-valve or PWM input is utilized, a control method must be chosen to compensate for the non-linearities of the system. In [START_REF] Rosas-Flores | Optimal linearization of the dynamic behavior of an on/off actuated single pneumatic cylinder[END_REF] the nonlinear system is transformed into an equivalent linear system (controlled by PWM). Because of the highly nonlinear nature of a pneumatic system it is better to address its control with a nonlinear model based control. Other research areas use non-linear averaging techniques based on a sliding mode approach to control pneumatic systems [START_REF] Girin | High-order slidingmode controllers of an electropneumatic actuator: Application to an aeronautic benchmark[END_REF], [START_REF] Xing | Tracking control of pneumatic artificial muscle actuators based on sliding mode and non-linear disturbance observer[END_REF], [START_REF] Shen | Nonlinear modelbased control of pulse width modulated pneumatic servo systems[END_REF], [START_REF] Le | Sliding mode control of a pneumatic haptic teleoperation system with on/off solenoid valves[END_REF].

Sliding mode control is a form of variable structure control that alters the dynamics of a nonlinear system by the application of a high-frequency switching control [START_REF] Utkin | Sliding mode control design principles and applications to electric drives[END_REF]. Sliding mode control can account for the non-linearities of the system as well as its dynamic uncertainties [START_REF] Nguyen | Accurate slidingmode control of pneumatic systems using low-cost solenoid valves[END_REF]. In this paper we consider a pneumatic actuator comprised of two chambers as shown in Figure 1. Each chamber has two solenoid valves, each of which can be either locked or connected to a compressed air supply (source pressure) or to exhaust (atmosphere pressure). Since we cannot have a chamber connected to both pressure and exhaust at the same time, only nine discrete modes exist. We assume the three modes where both chambers are locked, venting, and pressurizing are functionally equivalent, then we can assume the system has a total of seven unique discrete modes.

The sliding control system utilizes discrete operating modes defined directly from the state of the on/off solenoid valves, with the switching between these discrete modes decided based on the current tracking error, thereby bypassing the need for a PWM input. In [START_REF] Nguyen | Accurate slidingmode control of pneumatic systems using low-cost solenoid valves[END_REF], three discrete modes are considered for a two-chamber actuator (see Figure 1): two modes defined by connecting one chamber to exhaust and connecting one chamber to pressure (and the reverse thereof) and one mode defined by closing all valves.

In this paper, we expand the three-mode open-loop model of [START_REF] Nguyen | Accurate slidingmode control of pneumatic systems using low-cost solenoid valves[END_REF] into a seven-mode open-loop model. The new system has four extra modes defined by opening only one of the four solenoid valves at a time, which facilitate appropriate amounts of drive energy for good positioning precision and reduced switching activity.

Limiting the drive power allows us to make smaller and fewer adjustments, improving the positioning accuracy of the controlled system. Compared to the 3 mode paper we are extending to 7 modes and providing guidelines for the parameters τ , β and , which will be discussed later.

We will decide the thresholds for transition from one control mode to another control mode based on the full state of the system: piston position, velocity, and acceleration as well as chamber pressures. By optimizing these transition thresholds, we can decrease the open-close activity of the on/off solenoid valves, which increases their operating lifespan.

The organization of this paper is as follows: The modeling of the pneumatic actuator with its chamber and 4 solenoid valves is reported in Section II. The discrete input model of the actuator is given in Section III. The design and theoretical analysis of the sliding controller are discussed in Section IV. The simulation study to verify the control laws are presented in this paper in Section V The experimental results are presented in this paper in Section VI. Finally, the concluding remarks are presented in Section VII.

II. MODELING OF THE OPEN-LOOP PNEUMATIC ACTUATOR

This section of the paper derives the open-loop model of the pneumatic actuator including its on/off solenoid valves as shown in Figure 1. To describe the air flow dynamics in a cylinder, we assume that: Air is a perfect gas and its kinetic energy is negligible in the chamber, the pressure and the temperature are homogeneous in each chamber, the evolution of the gas in each chamber is polytropic, with the exception of some controller thresholds which were selected using isothermic assumptions. The temperature variation in chambers is negligible with respect to the supply temperature. The mass flow rate leakages are negligible, and the supply and exhaust pressures are constant.

The Simscape simulations in Section V does not assume/model all of these assumptions.

A. Model of the Pneumatic Chambers

If the charging and discharging of the chambers' cylinders are assumed to be polytropic, then the pressure dynamics of these chambers can be approximated [START_REF] Le | Transparency of a pneumatic teleoperation system using on/off solenoid valves[END_REF] as

ṖP = k V P (rT Q P -AP P ẏ) ṖN = k V N (rT Q N + AP N ẏ) (1)
where P P and P N refer to pressures (Pa) inside the chambers P and N, respectively, V P and V N refer to volumes (m 3 ) of the chambers P and N, respectively, Q P and Q N refer to mass flow rates (kg/s) of the chambers P and N, respectively, A refers to the piston cylinder area (m 2 ), T refers to the chamber temperature (K), k refers to the polytropic constant, r refers to the universal gas constant (J/(kgK)), and y refers to the piston position (m) shown in Figure 1. Note the arrows for position y, Force (N) F Ext and F st shown in Figure 1 refer to their positive directions, a negative magnitude to any of these quantities will refer to the opposite direction.

B. Model of the valves

The mass flow rates Q P and Q N can be derived in terms of the discrete voltage inputs U 1 , U 2 , U 3 and U 4 shown in Figure 1 and the continuous pressure inputs P P and P N :

Q P = U 1 Q(P S , P P ) -U 2 Q(P P , P E ) (2a)
Q N = U 3 Q(P S , P N ) -U 4 Q(P N , P E ) (2b)
Here, P S and P E are the pressures of the supply and the exhaust. In general, Q(P U p ,P Down ) used in (2a) and (2b), in which P U p is the upstream pressure and P Down is the downstream pressure, refers to the expression for the mass flow rate through an orifice. This generalized model has two parameters to describe the mass flow rate: the critical pressure ratio b = 0.433 and the sonic conductance C val (mass flow rate constant) [START_REF] Beater | Pneumatic Drives -System Design, Modelling and Control[END_REF]:

Q(P U p , P Down ) = C val P U p T Atm T U p ×      1 - P Down P U p -b 1-b 2 , if P Down P U p > b (subsonic) 1 , if P Down P U p ≤ b (choked) (3)
In the above, T U p is the upstream temperature of air and T Atm is the atmospheric temperature. Also, C V al is a characteristic of the valve.

C. Model of the Piston

Finally, the dynamics of the mechanical actuator involving the applied force on the piston and the resulting piston motion is

A(P P -P N ) -b m ẏ + F Ext -F St = M m ÿ (4)
where b m is the viscose coefficient (N s/m)), M m is the total mass of the load and the piston (Kg), F St is the stiction force (N), and F Ext is the external force (N). For simplicity, the stiction force and external force is assumed to be negligible. The F St was considered to be negligible since the pneumatic actuator used in experiment was an Airpel anti-stiction cylinder (ALL AIR Inc, New York, US).

III. DISCRETE INPUT MODEL OF THE OPEN-LOOP ACTUATOR

It is possible to combine the equations in Section II to write the dynamics of the open-loop pneumatic actuator in a 7-mode discrete-input. Differentiating (4) and substituting

(1) in it, the dynamics of the actuator are obtained as ...

y = F (Z) + krT M m Q P l/2 + y - Q N l/2 -y (5)
where Z = {y, ẏ, ÿ, P P , P N } is the state vector and

F (Z) = -b m M m ÿ - Ak M m P P l/2 + y + P N l/2 -y ẏ
where l is the total length of the chamber. We find that there are a total of nine discrete modes for the solenoid valves [START_REF] Le | Transparency of a pneumatic teleoperation system using on/off solenoid valves[END_REF]. These modes are shown in Table I. This paper will outline a controller that only utilizes modes M 1 through M 7 of the open-loop system. For the seven discrete modes, dynamic

M 1 M 2 M 3 M 4 M 5 M 6 M 7 M 8 M 9 U 1 0 1 0 0 0 1 0 0 1 U 2 0 0 1 0 0 0 1 1 0 U 3 0 0 0 0 1 0 1 0 1 U 4 0 0 0 1 0 1 0 1 0
equation can be obtained by substituting (2a) and (2b) into [START_REF] Nguyen | Accurate slidingmode control of pneumatic systems using low-cost solenoid valves[END_REF]. We obtain ...

y = F (Z) , mode M 1 F (Z) + (-1) i B i (Z) , mode M i = M 1 (6)
B 2 (Z) = krT M m Q(P S , P P ) (l/2 + y) B 3 (Z) = krT M m Q(P P , P E ) (l/2 + y) B 4 (Z) = krT M m Q(P N , P E ) (l/2 -y) B 5 (Z) = krT M m Q(P S , P N ) (l/2 -y) B 6 (Z) = B 2 (Z) + B 4 (Z) B 7 (Z) = B 5 (Z) + B 3 (Z)
Note that because P E ≤ P P ≤ P S , P E ≤ P N ≤ P S and -l/2 ≤ y ≤ l/2 and mass flow rates are non-negative functions B 2 (Z) through B 7 (Z) are all positive or equal to zero. The value i takes a range of integer values from 2 to 7.

IV. THE SLIDING MODE CONTROLLER DESIGN

For a position-controlled system, we can define the following sliding surface s = 0 where s is defined as:

s = ë ω 2 + 2ξ ė ω + e ( 7
)
where e is the position error y -y d , y is the actual position, y d is the desired position, and ξ and ω are constant and positive numbers. We use this function s and invoke the seven different modes of the open-loop system based on five different regions of the function s. These regions of s and the selected operating mode of the system are illustrated in Table II. Please note that the input voltages (i.e., control actions) for each mode is listed in Table I.

V is a positive-valued function, therefore if V < 0, V and |s| will be decreasing. Assuming s is initially bounded, s will asymptotically approach zero if we control the system so that

V = ṡs < -η|s| (9)
for some constant η > 0 [8] [START_REF] Nguyen | Accurate slidingmode control of pneumatic systems using low-cost solenoid valves[END_REF]. Take the derivative of ( 7) and substitute ( 6) to obtain

ṡ = λ , mode M 1 λ + (-1) i B i (Z)/ω 2 , mode M i = M 1 (10)
where λ = (F (Z) -... y d )/ω 2 + 2ξë/ω + ė. λ will be bounded if ẏ, ÿ, ẏd , ÿd , ... y d , P P and P N are bounded. P P and P N are bounded between P S and P E . y d is a controlled input to the system thus y d , ẏd , ÿd , ... y d are assumed to be bounded. If we rewrite the (4) as

M m ÿ + b m ẏ = A(P P -P N ) (11)
we see that the right side of ( 11) is always bounded. Therefore velocity is a 1 st order differential equation, which means velocity is a decaying exponential. Therefore if ẏ is initially bounded, then ẏ will always be bounded. If we rewrite (4)

ÿ = 1 M m (A(P P -P N ) -b m ẏ) (12)
We find that it is defined purely in terms of bounded functions, thus it too must always be bounded. Therefore λ is also bounded. Therefore if the positive-valued functions B i (Z) are sufficiently large then modes M 2 , M 4 and M 6 can ensure ṡ > η and the modes M 3 , M 5 and M 7 can ensure ṡ < -η.

All six functions B i (Z) are linearly proportional to C val and P S , the valve's mass flow rate constant in (3), thus choosing a large enough valve or supply pressure will ensure that these scalar functions will be sufficient in magnitude. Thus, using the modes M 2 , M 4 and M 6 when s < 0, and using M 3 , M 5 and M 7 when s > 0. Will ensure [START_REF] Le | Transparency of a pneumatic teleoperation system using on/off solenoid valves[END_REF], and thus the convergence of s to zero over time.

if s converges to zero and y d , ẏd , ÿd are bounded the output y, ẏ, ÿ will also be bounded and thus the system will be BIBO stable.

B. Controller Mode Selection

The 7-mode controller requires knowledge of the current chamber pressures to pick the appropriate operating modes. This requires additional sensors as compared to the 3-mode controller. These sensors are required for deciding which control mode to use when β ≥ |s| > (see Table II).

We will first study the negative region -β ≥ s > -, if we evaluate (4) under mode M 2 and assume that the filling chamber has a sufficient amount of time to be fully pressurized we can see that the ÿ ∝ (P S -P N ). If we evaluate (4) under mode M 4 and assume that the venting chamber has a sufficient amount of time to be fully vented we can see that the ÿ ∝ (P P -P E ) If we define

E 1 = (P S -P N ) -(P P -P E ) = (P S + P E ) -(P P + P N ) ( 13
)
the magnitude of E 1 is positive when the pressure difference P S -P N is greater than the pressure difference P P -P E . Therefore, when E 1 is positive, the appropriate operating mode for the region -≥ s > -β is M 2 as it will result in a higher piston acceleration compared to the mode M 4 .

Conversely when the magnitude of E 1 is negative M 4 will result in a higher piston acceleration compared to M 2 . In the positive region β ≥ s > , comparing modes M 5 and M 3 , from ( 4) we see that the magnitudes of ÿ are based on the pressure differences P S -P P and P N -P E respectively. Let us define

E 2 = (P S -P P ) -(P N -P E ) = E 1 (14)
The magnitude of E 2 is positive when the pressure difference P S -P P is greater then the pressure difference P N -P E .

Therefore when E 2 is positive the appropriate operating mode for the region ≥ s > β is M 5 as it will accelerate the piston more compared to mode M 3 . Conversely when the magnitude of E 2 is negative mode M 3 will result in a higher piston acceleration compared to M 5 In summary, the magnitude of s in ( 7) and the magnitude of E 1 in (13) can be used by the controller to select the best of the seven operating modes see Figure 2.

C. Selecting Parameters τ , β and

This section proposes appropriate ways for selecting the controller parameters of τ , β and for smoothest motions and least switching activity.

1) Selecting τ : The purpose of the timeout parameter τ is to reduce switching between the modes used in the region β ≥ |s| > by enforcing a minimum amount of time between modes transitions in this region. A larger τ value will reduce switching. A smaller τ value will give a more accurate output. There is a trade-off between these goals, therefore an appropriate value of τ should be determined by evaluating the open-loop responses M 2 and M 4 . For initial conditions (at t = 0) of y = 0 and V P = V N = A l 2 (i.e., the piston positioned in halfway along the cylinder length), using (1) -( 4) we can simulate the open-loop system response in each of these modes (Figure 3). These initial conditions where selected because they are symmetric in nature.

For the pressurizing case M 2 the intial conditions P P = P N = P E was selected. For the venting case M 4 the intial conditions P P = P N = P S was selected.

We define t f ill as the time it takes to reach the maximum velocity max( ẏ) for the M 2 Chamber P filling profile and t vent as the time it takes to reach the maximum velocity max( ẏ) for the M 4 Chamber N venting profile. If we compare t f ill and t vent obtained from simulation we see that t f ill is a shorter span of time. For our system tests we have selected τ = t f ill .

2) Selecting β: The magnitude β is the transition threshold between a mode connecting the two chambers to supply and exhaust pressures (i.e., M 6 or M 7 ), and the alternating modes for opening only a single valve (i.e., M 2 /M 4 or M 3 /M 5 ).

For our purposes β was chosen empirically. We used the following equation as a starting point. β = 3A(P S -P E ) 8ωb m (15)

A higher β was found to be more optimal for y ref signals with large variation and a smaller β was found to be more optimal for y ref signals with smaller variation.

For our experiments the we selected a β = 3mm. This selection of β is not unique however it has been demonstrated in our results as being appropriate.

3) Selecting : When the magnitude of s is less than a positive valued but small , the controller removes the actuation from the system by closing all valves. The magnitude of should ideally be selected so that when this occurs and the effect of past actuations settles out, the position difference, y -y d , will be less than some desired small amount e min . Therefore, is selected as = e min

V. SIMULATION STUDY

To analyze the performance of the 7-mode controller described in this paper we will be comparing it against the original 3-mode controller it was based on in [START_REF] Nguyen | Accurate slidingmode control of pneumatic systems using low-cost solenoid valves[END_REF].

A. Simulation Parameters and Test Inputs

For our simulation, we selected the model parameters listed in Table III. These model parameters correspond to the experimental setup that will be used in Section VI. The 0.9 kg Total Mass of load following controller parameters were selected: ω = 100 rad/s, ξ = 0.5, τ = 40ms, β = 3mm and = 1mm. To model our system we utilized the Simulink Simscape toolbox.

B. Simulation Results

The simulation was run utilizing a 40mm peak-peak sine wave test input with frequencies varying from 0.1Hz -3.0Hz, the results of the simulations are charted in Figure 4. From these results, we find that for both the 3-mode and the 7mode systems increasing the input frequency leads to an increasing RMS tracking error In the 7 mode system, for moderate differences between y d and y, this difference may not be enough to switch from the |s| > β region. As we can see this can cause an increase in tracking error in the 2-3hz range for the 7-mode system. This is one of the costs of using the 7-mode sliding control algorithm. The advantage of the 7-mode controller algorithm is the reduction in the solenoid valves switching as observed in Figure 4. At 1.5Hz, this is equal to 63% reduction in switching activity and 0.45mm improvement in tracking error.

To test the system with a more complicated signal, we used an input that is a summation of 8 sine waves with 8 different frequencies. As we can see from Figure 5, both controllers have no trouble tracking the multi-sine wave; however the switching activity for the 7-mode controller is greatly reduced compared to the 3-mode controller. From the 3-mode controller to the 7-mode controller, there is a 64% reduction in solenoid switching activity.

VI. EXPERIMENTAL TESTING

In this section, experiments with a 1-DOF system are reported. This experimental setup is the same as the one described in [START_REF] Le | Development of a hybrid control for a pneumatic teleoperation system using on/off solenoid valves[END_REF] for a single actuator. For this experiment, the following controller parameters were selected were the same as in simulation except ω = 60 rad/s.

A. Experimental Results

The experiment was run utilizing a 40mm peak-peak sine wave test input with frequencies varying from 0.1Hz -3.0Hz. The recorded results from these experiments are charted in Figure 6. When we compare these results for the 3-mode controller and the 7-mode controller, we can see that there was a notable improvement to tracking performance and switching activity for the 7-mode controller (compared to the 3-mode controller).

To test the system with a more complicated signal, we used an input that is a summation of 8 sine waves with 8 different frequencies (exact same as in simulation). As we can see from Figure 5, both controllers have no trouble tracking the multi-sine wave; however the switching activity for the 7-mode controller is reduced compared to the 3mode controller. From the 3-mode controller to the 7-mode controller, there is a 49% reduction in solenoid switching activity.

To test the load disturbance rejection of the algorithm a 1kg weight was attached to the actuator using a rope and a pulley. The experiment was run utilizing the sine wave test input with a 1.5Hz frequency. The recorded results are shown in Figure 8. Tracking error for both 3 mode and 7 mode increased by 0.2mm RMS. Switching was found to increase for both 3 mode and 7 mode by ∼ 10%.

VII. CONCLUDING REMARKS

This paper proposed a sliding mode law for precise position control with minimal switching activity designed for use on a pneumatic actuator. The two-chamber actuator with four on/off solenoid valves had a total of sixteen possible modes only seven of which were considered functional and unique and were used for this particular control law. This controller uses four more additional modes compared to [START_REF] Nguyen | Accurate slidingmode control of pneumatic systems using low-cost solenoid valves[END_REF] that decrease the coarseness in the drive force for lower

Fig. 2 . 7 -

 27 Fig. 2. 7-mode controller diagram

Fig. 4 .Fig. 5 .

 45 Fig. 4. Sine-wave simulation results

Fig. 6 .Fig. 7 .

 67 Fig. 6. Sine-wave experimental results

Fig. 8 .

 8 Fig. 8. Sine Wave 1.5Hz with 1kg load attached to actuator

TABLE II SELECTION

 II OF THE OPERATING MODE BASED ON POSITIONING ERROR s

	Region	Discrete	Magnitude of	ṡ
		operating modes	from (10)
	s > β	M 7			Large negative
	β ≥ s >	M 3 and M 5	Modest negative
	≥ s > -	M 1			Minimal
	-≥ s > -β	M 2 and M 4	Modest positive
	-β ≥ s	M 6			Large positive
	A. Stability
	To be able to analyze stability, consider the Lyapunov
	function candidate	V =	1 2	s 2

  val 3.4 × 10 -9 kg/(s Pa) Mass Flow Rate Const.

		TABLE III
		SYSTEM PARAMETERS TABLE
	Var.	Value	Label
	l	0.1m	Chamber Length
	T	23 C	Chamber Temperature
	C P S	300, 000 Pa	Supply Air Pressure
	P E	100, 000 Pa	Exhaust Air Pressure
	k	1.2	polytropic constant
	A	1.814 cm 2	Piston Cylinder Area
	b m	50 (N s)/m	Viscosity Coefficient
	M m


Élise Janvresse
email: elise.janvresse@univ-rouen.fr

Thierry De

L A Rue

AVERAGING ALONG UNIFORM RANDOM INTEGERS

Keywords: 2000 Mathematics Subject Classification. 11A63, 11B05, 60G50, 60G55 Benford's law, log-density, H∞-density, uniform random integers

Motivated by giving a meaning to "The probability that a random integer has initial digit d", we define a URI-set as a random set E of natural integers such that each n ≥ 1 belongs to E with probability 1/n, independently of other integers. This enables us to introduce two notions of densities on natural numbers: The URI-density, obtained by averaging along the elements of E, and the local URI-density, which we get by considering the k-th element of E and letting k go to ∞. We prove that the elements of E satisfy Benford's law, both in the sense of URI-density and in the sense of local URI-density. Moreover, if b 1 and b 2 are two multiplicatively independent integers, then the mantissae of a natural number in base b 1 and in base b 2 are independent. Connections of URI-density and local URI-density with other well-known notions of densities are established: Both are stronger than the natural density, and URI-density is equivalent to log-density. We also give a stochastic interpretation, in terms of URI-set, of the H∞-density.

1. Introduction 1.1. Benford's law and Flehinger's theorem. Benford's law describes the empirical distribution of the leading digit of everyday-life numbers. It was first discovered by the astronomer Simon Newcomb in 1881 [START_REF] Newcomb | Note on the frequency of use of the different digits in natural numbers[END_REF] and named after the physicist Franck Benford who independently rediscovered the phenomenon in 1938 [START_REF] Benford | The law of anomalous numbers[END_REF]. According to this law, the proportion of numbers in large series of empirical data with leading digit d ∈ {1, 2, . . . , 9} is log 10 (1 + 1/d). More generally, defining the mantissa M (x) of a positive real number x as the only real number in [1, 10[ such that x = M (x)10 k for some integer k, Benford's law states that for any 1 ≤ α < β < 10, the proportion of numbers whose mantissa lies in [α, β] is log 10 βlog 10 α.

Giving Benford's law a mathematical meaning requires to formalize the notion of "everyday-life numbers", which is far from obvious. However there have been many attempts to explain mathematically the ubiquity of this distribution in empirical datasets. One of them is Betty J. Flehinger's theorem, published in 1965 in an article entitled On the probability that a random integer has initial digit A [START_REF] Flehinger | On the probability that a random integer has initial digit A[END_REF]. It occurred to Flehinger that the most natural set of numbers on which we should verify Benford's distribution is the whole set of positive integers. Unfortunately, defining

P 1 n (d) := 1 n n j=1 ½ [d,d+1[ M (j)
(that is the proportion of integers between 1 and n with leading digit d), we see that the sequence P 1 n (d) has no limit as n → ∞: It oscillates over longer and longer periods. Flehinger's idea was then to seek the limit by iteration of the process of Cesaro averaging: She inductively set which is the proportion predicted by Benford's law. (Donald Knuth generalized Flehinger's theorem to the distribution of the whole mantissa in 1981 [START_REF] Knuth | The art of computer programming[END_REF].)

P
In spite of its title, Flehinger's article has no probabilistic content. A good reason is that there is no way of picking an integer uniformly at random in the set of all natural numbers. The first motivation for the present work was nevertheless to translate Flehinger's theorem in the context of probability theory: How can we interpret the probability that a (random) integer has a given initial digit? Our purpose is thus to give a meaning to the sentence "An integer picked uniformly at random has such a property". 1.2. Roadmap of the paper. We construct in Section 2 a random infinite set E of integers, which we call a URI-set, such that averaging along the elements of E reflects the expected behaviour of a random integer. This random set enables us to introduce two notions of densities on natural numbers: The URI-density (see Section 3.1), obtained by averaging along the elements of E, and the local URIdensity (see Section 5.2), which we get by considering the k-th element of E and letting k go to ∞. We prove that the elements of E satisfy Benford's law, both in the sense of URI-density (Theorem 3.5), and in the sense of local URI-density (Theorem 5.1). Our point of view also enables to consider simultaneously the mantissae of a number in different bases, and in particular we prove a result which can be interpreted as follows: If b 1 and b 2 are two multiplicatively independent integers, then the mantissae of a natural number in base b 1 and in base b 2 are independent (Theorem 4.1). Connections of URI-density and local URI-density with other well-known notions of densities are established: We prove that both are stronger than the natural density (Theorems 3.3 and 5.2), and that in fact URI-density is equivalent to log-density (Theorem 3.6). We finish in Section 6 by giving a stochastic interpretation, in terms of URI-set, of the H ∞ -density used in Flehinger's theorem, and by raising some open problems.

The construction of the random set of integers is inspired by a previous article by the same authors [START_REF] Janvresse | From uniform distributions to Benford's law[END_REF], where a probabilistic proof of Flehinger's theorem was provided. We summarize this proof in Section 2.1

Acknowledgements. The authors are grateful to Bodo Volkmann for stimulating questions.

Uniform random set of integers

P(M k ∈ S) -µ B (S) ≤ 9 10 k .
A connection is made between the quantities P k n (d), n ≥ 1, and the k-th step of our Markov chain: It is established in [START_REF] Janvresse | From uniform distributions to Benford's law[END_REF] that for all a ∈ [1, 10[ and all k ≥ 1,

lim j→∞ P k ⌊a10 j ⌋ (d) = P M k ∈ [d, d + 1[ M 0 = a .
A proof of (1), with an estimation of the speed of convergence, follows: We get for all k ≥ 1 lim inf

n→∞ P k n (d) -µ B [d, d + 1[ ≤ 9 10
k and lim sup n→∞

P k n (d) -µ B [d, d + 1[ ≤ 9 10 k . (4)
2.2. Construction of a random set of integers. We can interpret the Markov chain (M k ) as the sequence of mantissae of positive random variables X k , where the sequence (X k ) is itself a Markov chain such that, given X 0 , . . . , X k , X k+1 is uniformly distributed in ]0, X k [. Then X := {X k , k ≥ 0} is a discrete random set of real numbers which satisfies the following property: For any t > 0, conditionally to the fact that X∩]t, ∞[ = ∅, max(X∩]0, t]) is uniformly distributed in ]0, t], and independent of X∩]t, ∞[. Our idea is thus to imitate the structure of this random set of reals, but inside the set of natural numbers. We are looking for a random infinite set of integers E satisfying (U) for all n ≥ 1, max(E ∩ {1, . . . , n}) is uniformly distributed in {1, . . . , n};

(I) for all n ≥ 1, max(E ∩ {1, . . . , n}) is independent of E ∩ {n + 1, n + 2, . . .}.

For such a random set E, we must have by (U), for each n ≥ 1,

P(n ∈ E) = P max(E ∩ {1, . . . , n}) = n = 1/n,
and (I) implies that all events (n ∈ E) are independent.

Conversely, picking elements of E using independent Bernoulli random variables, with P(n ∈ E) = 1/n for each n ≥ 1 gives a random set satisfying the required conditions. Indeed, for each j ∈ {1, . . . , n}, we get

P max(E ∩ {1, . . . , n}) = j = P j ∈ E, j + 1 / ∈ E, . . . , n / ∈ E = 1 j j j + 1 • • • n -1 n = 1 n .
Observe also that, with probability 1, the cardinality of E is infinite. Because of the uniformity property (U), such a random set E appears as a good way to modelize the uniform distribution in the set of natural numbers and will therefore be referred to as a set of uniform random integers, or URI-set.

URI-density and Benford's law

From now on, E denotes a URI-set, and we denote its ordered elements by

E = {N 1 = 1 < N 2 < . . . < N k < . . .}.
For each n ≥ 1, we set E n := E ∩ {1, . . . , n}.

It will be useful to give the following estimation of

|E n |. Lemma 3.1. |E n | ln n a.s.
----→ n→∞ 1.

We recall Theorem 12 page 272 in Petrov's book [START_REF] Valentin | Sums of independent random variables[END_REF]:

Theorem 3.2. Let (Z n ) be a sequence of independent centered real-valued random variables Z n . If t n ↑ ∞ and n≥1 E [Z p n ] t p n < ∞, for some p, 1 ≤ p ≤ 2, then n j=1 Z j t j a.s.
----→ n→∞ 0.

Proof of Lemma 3.1. We consider the independent centered random variables

Z n := ½ E (n) -1/n. Since n≥1 E Z 2 n (ln n) 2 ≤ n≥1 1 - 1 n 1 n(ln n) 2 < ∞, we get by Theorem 3.2 (5) n j=1 Z j ln n = 1 ln n n j=1 ½ E (j) - 1 ln n n j=1 1/j a.s. ----→ n→∞ 0 .
This concludes the proof of the lemma.

3.1. URI-density. We say that a subset A of the set of natural numbers has

URI-density α if 1 n n k=1 ½ A (N k ) a.s.
----→ n→∞ α.

Note that an equivalent formulation is 1

|E n | j∈En ½ A (j) a.s. ----→ n→∞ α.
As we expect, the URI-density generalizes the notion of natural density.

Theorem 3.3. Let A ⊂ Z + . If 1 n n j=1 ½ A (j) ----→ n→∞ α, then A has URI-density α.
Proof. When considering the elements of E n , it will be convenient to order them backwards:

E n = E ∩ {1, . . . , n} = Y (n) 1 > Y (n) 2 > . . . > Y (n) |En| = 1 .
For each n ≥ 1, let a n := ½ A (n). We are going to prove the result in the form

(6) 1 |E n | j∈En a j a.s.
----→ n→∞ α.

We split the above average as

(7) 1 |E n | j∈En a j = 1 |E n | |En| i=1 a Y (n) i -K n i + 1 |E n | |En| i=1 K n i ,
where

K n i :=    E a Y (n) i Y (n) i-1 , . . . , Y (n) 1 (2 ≤ i ≤ |E n |), E a Y (n) 1 (i = 1).
We first deal with the second term of [START_REF] Hernández | Markov chains and invariant probabilities[END_REF]. By Properties (U) and (I), Y

is uniformly distributed in {1, . . . , n}, and conditionally to Y

(n) i-1 , . . . , Y (n) 1 , Y (n) i is uniformly distributed in 1, . . . , Y (n)
i-1 -1 , as long as Y (n) i-1 > 1. Hence, K n 1 = 1 n n j=1 a j and K n i = 1 Y (n) i-1 -1 Y (n) i-1 -1 j=1 a j , (2 ≤ i ≤ |E n |).
By hypothesis, as soon as

Y (n) i-1 is large, K n i is close to α. For any fixed ε > 0, the number of i's such that |K n i -α| > ε is bounded independently of n. Since |E n | → ∞ a.s., it follows that 1 |E n | |En| i=1 K n i a.s. ----→ n→∞ α .
Let us now turn to the first term of [START_REF] Hernández | Markov chains and invariant probabilities[END_REF]. Define for any n ≥ 1 and i ≥ 1

A n i := a Y (n) i -K n i (1 ≤ i ≤ |E n |), 0 (i > |E n |).
In order to prove [START_REF] Flehinger | On the probability that a random integer has initial digit A[END_REF], it remains to show that

(8) 1 |E n | |En| i=1 A n i a.s. ----→ n→∞ 0.
Using a standard method, we first prove that (8) holds along a subsequence, then we control the oscillations to conclude. Consider n m := ⌊exp(m 2 )⌋. By lemma 3.1, convergence along the subsequence (n m ) amounts to

A(n m ) := 1 ⌊ln n m ⌋ ⌊ln nm⌋ i=1 A nm i a.s. ----→ m→∞ 0 .
Observe now that the variance of A nm i is bounded, and that for i = j, E A nm i A nm j = 0. Therefore, the variance of A(n m ) is of order (ln n m ) -1 = 1/m 2 . By Tchebychev's inequality,

m≥1 P A(n m ) > m -1/4 ≤ m≥1 Var(A(n m )) m -1/2 = m≥1 m 1/2 O 1 m 2 < ∞.
Hence, by Borel-Cantelli, with probability 1 we have A(n m ) ≤ m -1/4 for m large enough. This proves that

1 |E nm | |En m | i=1 A nm i a.s.
----→ m→∞ 0.

Consider now an integer

n ∈]n m , n m+1 [. We can write (9) 1 |E n | |En| i=1 A n i = 1 |E n |   |En| i=1 A n i - |En m | i=1 A nm i   + 1 |E n | - 1 |E nm | |En m | i=1 A nm i + 1 |E nm | |En m | i=1 A nm i Since A nm i is bounded and |E nm | -|E n | = o(|E n |), the second term on the RHS vanishes as m → ∞. Moreover, E nm ⊂ E n . Therefore each A nm i in the first term (except A nm 1
which has a slightly different definition) is annihilated by some A n j , and the first term of (9) reduces to

1 |E n | |En|-|En m |+1 i=1 A n i - 1 |E n | A nm 1 ,
which goes to zero as m → ∞. Since we already know that the third term goes to zero, this proves (8).

3.2.

Benford's law. We say that a sequence of positive real numbers (x n ) follows Benford's law if for all 1 ≤ t < 10

1 n n j=1 ½ M (xj )<t ----→ n→∞ log 10 t.
Remark 3.4. We recall that this is equivalent to the uniform distribution mod 1 of the sequence (log 10 x k ) (see e.g. [START_REF]The distribution of leading digits and uniform distribution mod 1[END_REF]). Therefore, this is also equivalent to the fact that the sequence (1/x k ) follows Benford's law.

The following theorem shows that the elements of the URI-set E almost surely follow Benford's law. Theorem 3.5. For any 1 ≤ t ≤ 10, the URI-density of {n ≥ 1 : M (n) < t} is log 10 t.

Proof. By Duncan's work [START_REF] Luce | Note on the initial digit problem[END_REF], this result can be viewed as a corollary of the equivalence between URI-density and log-density (see Theorem 3.6 below). However we provide a direct proof of it, in which useful ideas will be presented.

It is convenient to consider a coupling of the URI-set E and its continuous analog defined as follows: Let ξ be a Poisson process on R * + with intensity 1/x. It can be viewed as a random set of points. For any interval I, let ξ I denote the number of points in I ∩ ξ: ξ I is Poisson distributed with parameter I dx x . From ξ, define the random set E as the set of integers n ≥ 1 such that ξ ]n-1,n] ≥ 1. Since the random variables (ξ ]n-1,n] ) n≥1 are independent and

P(n ∈ E) = 1 -P ξ ]n-1,n] = 0 = 1 n , E is a URI-set.
The process ξ satisfies a property analogous to Property (U): For any a > 0, the largest point of [0, a] ∩ ξ is uniformly distributed in [0, a]. Indeed, for any t < a

P max([0, a] ∩ ξ) ≤ t = P ξ [t,a] = 0 = t a .
Let us order backwards the points of ξ ----→ n→∞ log 10 t.

∩ [0, 1]: 1 > Y 1 > Y 2 > . . .. Conditionally to Y k , Y k+1 is uniformly distributed in [0, Y k ]. By Proposition 3.1 of [8], the mantissae (M (Y k )) constitute a Markov chain whose unique invariant distribution is µ B . Moreover, Y 1 being uniformly distributed in [0, 1], the distribution of M (Y 1 )

Now observe that

n≥1 P(ξ ]n-1,n] ≥ 2) = n≥1 1 n + n -1 n log n -1 n < ∞,
hence by Borel-Cantelli, with probability one there is only a finite number of n's

such that ξ ]n-1,n] ≥ 2. Coming back to the URI-set E = {N 1 = 1 < N 2 < • • • }, this implies the almost-sure existence of R such that, for all large enough k, 0 ≤ N k -X k+R < 1. Since N k → ∞ a.s., we have |M (N k ) -M (X k+R )| → 0 unless N k be of
P k n = α,
where the P k N 's are inductively defined by P 0 n := ½ A (n) and P k+1 n

:= (1/n) n j=1 P k j .
Obviously, H ∞ -density is stronger than natural density, in the sense used by Diaconis in [START_REF] Diaconis | Weak and strong averages in probability and the theory of numbers[END_REF]: If A has natural density α, then A has H ∞ -density α. The example of the set A of integer whose initial digit (when written in base 10) is 1 shows that H ∞ -density is strictly stronger than natural density.

Still stronger than H ∞ -density is the notion of log-density.

Recall that A ⊂ Z + has log-density α if 1 ln n n j=1 1 j ½ A (j) ----→ n→∞ α.
A proof that log-density is stronger than H ∞ -density can be found in [START_REF]The distribution of leading digits and uniform distribution mod 1[END_REF], together with an example of a set A with a log-density, but for which the H ∞ -density fails to exist.

Theorem 3.6. Let A ⊂ Z + . 1 |E n | n j=1 ½ A (j)½ E (j) - 1 ln n n j=1 ½ A (j)/j a.s. ----→ n→∞ 0 .
In particular, A has URI-density α if and only if A has log-density α.

Proof. We apply Theorem 3.2: We consider the independent random variables

Z n := ½ A (n)(½ E (n) -1/n) which are centered. Since n≥1 E Z 2 n (ln n) 2 ≤ n≥1 1 - 1 n 1 n(ln n) 2 < ∞, we get (11) n j=1 Z j ln n = 1 ln n n j=1 ½ A (j)½ E (j) - 1 ln n n j=1 ½ A (j)/j a.s. ----→ n→∞ 0 . Since, by Lemma 3.1 |E n | ln n a.s. ----→ n→∞ 1,
we can conclude the proof of Theorem 3.6.

Independence of mantissae in different bases

All the previous results concerned numbers written in the base-10 numeral system, but extend straightforwardlly to any integer base b. We denote by M b (x) ∈ The purpose of this section is to prove the following theorem, which states that the mantissae in different bases of the elements of the URI-set E are independent, under some algebraic condition on the bases. Theorem 4.1. Let (b i ) 1≤i≤ℓ be positive integers, satisfying [START_REF] Valentin | Sums of independent random variables[END_REF] ∀a 1 , . . . , a ℓ ∈ Z,

a 1 ln b 1 + • • • + a ℓ ln b ℓ = 0 =⇒ a 1 = • • • = a ℓ = 0.
Then for any

1 ≤ t i < b i (1 ≤ i ≤ ℓ), {n ∈ Z + : M bi (n) ≤ t i ∀1 ≤ i ≤ ℓ} has URI-density equal to ℓ i=1 log bi t i .
Recall that the positive integers (b i ) 1≤i≤ℓ are said to be multiplicatively independent if b s1 1 . . . b s ℓ ℓ = 1 where (s i ) 1≤i≤ℓ ⊂ Z implies that s i = 0 for all i. Note that, in the case ℓ = 2, property [START_REF] Valentin | Sums of independent random variables[END_REF] exactly means that b 1 and b 2 are multiplicatively independent. To our knowledge, it is unknown whether, in the general case, multiplicative independence of b 1 , . . . , b ℓ implies property [START_REF] Valentin | Sums of independent random variables[END_REF]. This question is related to the so-called Schanuel's conjecture in transcendental number theory (see [START_REF] Lang | Introduction to transcendental numbers[END_REF], p. 30-31 or [START_REF] Waldschmidt | Séminaire sur les nombres transcendants 1972-73[END_REF]). Lemma 4.2. Let (Z k ) k≥1 be i.i.d. random variables taking values in (R/Z) ℓ with common distribution ν. Assume that the only probability distribution µ such that µ * ν = µ is the Lebesgue measure on (R/Z) ℓ . Then the random walk

(P k ) k≥1 := (Z 1 + • • • + Z k ) k≥1 is uniformly distributed on (R/Z) ℓ . In other words, it satisfies: For all cylinder C = [u 1 , v 1 ] × • • • × [u ℓ , v ℓ ] where 0 ≤ u i < v i < 1, (13) 1 n n k=1 ½ C (P k ) a.s. ----→ n→∞ ℓ i=1 (v i -u i ).
Proof. Denote by µ the Lebesgue measure on (R/Z) ℓ . Let M 0 be a random variable with law µ, independent of (Z k ) k≥1 . Setting

M k := M 0 + Z 1 + • • • + Z k = M 0 + P k ,
we get a stationary random walk (M k ) k≥0 . Since µ is the unique invariant measure under convolution by ν, the stationary process (M k ) k≥0 is ergodic, and by Birkhoff ergodic theorem, we get that for all cylinder

C = [u 1 , v 1 ] × • • • × [u ℓ , v ℓ ] where 0 ≤ u i < v i < 1, 1 n n k=1 ½ C (M k ) a.s. ----→ n→∞ ℓ i=1 (v i -u i ).
(See e.g. [START_REF] Hernández | Markov chains and invariant probabilities[END_REF], Corollary 2.5.2 page 38.) Therefore, we can find some m 0 ∈ (R/Z) ℓ such that, with probability 1, for all cylinder

C = [u 1 , v 1 ] × • • • × [u ℓ , v ℓ ] where 0 ≤ u i < v i < 1 are rational numbers, 1 n n k=1 ½ C (m 0 + P k ) ----→ n→∞ ℓ i=1 (v i -u i ).
We thus obtain that, for all cylinder C with rational endpoints,

1 n n k=1 ½ C-m0 (P k ) a.s. ----→ n→∞ ℓ i=1 (v i -u i ) = µ(C -m 0 ).
By density of the rationals, ( 13) is satisfied for any cylinder C.

Proof of Theorem 4.1. As in the proof of Theorem 3.5, we consider the coupling of the URI-set with the Poisson process ξ, and we denote by

• • • > X 2 > X 1 > 1 > Y 1 > Y 2 > • • • the points of ξ. Define, for all k ≥ 1, U k := Y k /Y k-1 (
where

Y 0 := 1): Then (U k ) k≥1 is a sequence of i.i.d. uniform random variables in [0, 1],
and

Y k = U 1 U 2 . . . U k . Set Z k := log b1 (U k ) mod 1, . . . , log b ℓ (U k ) mod 1 ∈ (R/Z) ℓ
, and let ν be the common law of the Z k 's. We claim that the only probability measure µ which is invariant under convolution by ν is the Lebesgue measure on (R/Z) ℓ . Indeed, for such an invariant measure, the Fourier coefficients must satisfy

∀(m 1 , . . . , m ℓ ) ∈ Z ℓ , µ(m 1 , . . . , m ℓ ) = µ(m 1 , . . . , m ℓ ) ν(m 1 , . . . , m ℓ ).
We just have to check that ν(m 1 , . . . , m ℓ ) = 1 when (m 1 , . . . , m ℓ ) = (0, . . . , 0).

ν(m 1 , . . . , m ℓ ) = (R/Z) ℓ e -i2π(m1t1+•••+m ℓ t ℓ ) dν(t 1 , . . . , t ℓ ) = [0,1] e -i2π(m1 log b 1 u+•••+m ℓ log b ℓ u) du = [0,1]
e -i2πθ ln u du, where θ := m1 ln b1 + • • • + m ℓ ln b ℓ = 0 for (m 1 , . . . , m ℓ ) = (0, . . . , 0) by [START_REF] Valentin | Sums of independent random variables[END_REF]. Hence,

ν(m 1 , . . . , m ℓ ) = 1 1 -i2πθ = 1,

Local URI-density

For the sake of simplicity, we now return to the classical base 10 numeration system (but obviously the following results also hold in any integer base).

5.1.

A single element of a URI-set satisfies Benford's law. According to Theorem 3.6, Theorem 3.5 turns out to be weaker than Flehinger's theorem. However similar ideas to those developed in the proof can lead to somewhat stronger results than Theorem 3.5.

Theorem 5.1. For all 1 ≤ α < β < 10

lim k→∞ P M (N k ) ∈ [α, β] = µ B ([α, β]).
Proof. As in the proof of Theorem 3.5, we consider the URI-set E constructed from the Poisson process ξ. For a fixed integer m, we number the points of ξ∩]m, +∞[: We now consider the event

ξ∩]m, +∞[= {X m 1 < X m 2 < . . .}. We also set X m 0 := m. Observe that the process (1/X m k ) n≥0 is again a Markov chain such that, given (1/X m 0 , . . . , 1/X m k ), 1/X m k+1 is uniformly distributed in ]0, 1/X m k [. It
A m := n>m ξ ]n-1,n] ≤ 1 ∩ ℓ:10 ℓ >m ξ ]10 ℓ -1,10 ℓ ] = 0 .
Defining the random variable J m as the largest index such that N Jm ≤ m, the realization of A m ensures that the shift of index between the process (X m k ) and the integers N k that are larger than m + 1 remains constant, equal to J m . Therefore, for any k ≥ 1 0 ≤ N Jm+k -X m k < 1. Moreover, since A m also forbids that any N k > m be of the form 10 ℓ , we get (conditionally to A m )

(15) ∀k ≥ 1, 0 ≤ M (N Jm+k ) -M (X m k ) < 10 N Jm+k < 10 m .
Note also that P(A m ) → 1 as m → ∞, so that choosing m large enough will enable us to condition with respect to A m without affecting too much the probability of any event. Indeed, we will make use of the following inequality, valid for any events A and B with P(A) > 0:

(16) |P(B | A) -P(B)| ≤ P(A c ) P(A) .
Let us fix an arbitrary ε > 0. We choose m large enough so that ( 17)

P(A c m ) P(A m ) < ε and 10 m < ε.
Conditioning with respect to J m which takes values in {1, . . . , m}, we get

P M (N k ) ∈ [α, β] -µ B ([α, β]) ≤ m j=1 P(J m = j) P M (N k ) ∈ [α, β] | J m = j -µ B ([α, β]) .
Then we write, for any 1 ≤ j ≤ m,

P M (N k ) ∈ [α, β] | J m = j -µ B ([α, β]) ≤ D 1 + D 2 + D 3 + D 4 ,
where

D 1 := P M (N k ) ∈ [α, β] | J m = j -P M (N k ) ∈ [α, β] | A m , J m = j D 2 := P M (N Jm+k-j ) ∈ [α, β] | A m , J m = j -P M (X m k-j ) ∈ [α, β] | A m , J m = j , D 3 := P M (X m k-j ) ∈ [α, β] | A m , J m = j -P M (X m k-j ) ∈ [α, β] | J m = j D 4 := P M (X m k-j ) ∈ [α, β] | J m = j -µ B ([α, β]
) . Observe that A m , which is measurable with respect to the Poisson process on ]m, +∞[, is independent of (J m = j), which is measurable with respect to the Poisson process on ]1, m]. Hence, using ( 16) and ( 17), we can bound D 1 + D 3 by 2ε.

Again, X m k-j is measurable with respect to the Poisson process on ]m, +∞[, hence is independent of (J m = j). By ( 14), the contribution of D 4 can be bounded by (9/10) k-j , hence by (9/10) k-m .

It remains to deal with D 2 . Since everything is conditioned on A m , we can use (15) to get

D 2 ≤ P M (X m k-j ) ∈ [α -10/m, α] | A m , J m = j + P M (X m k-j ) ∈ [β -10/m, β] | A m , J m = j .
Using again (16), ( 17) and (14) yields

D 2 ≤ 2(9/10) k-m +2ε+µ B ([α -10/m, α])+µ B ([β -10/m, β]) ≤ 2(9/10) k-m +4ε.
Remark that the statement of Theorem 5.1 would not hold if we replace the interval [α, β] by any Borel set S. Indeed, since N k is an integer, the probability that its mantissa belong to the set of irrational numbers is zero. In other words, the convergence of the distribution of M (N k ) to Benford's law is only a weak convergence. However, if we denote by Xk the largest point of the Poisson process which is smaller than N k , then Xk ∈]N k -1, N k ], and the distribution of M ( Xk ) converges to Benford's law in total variation norm.

5.2.

Local URI-density is stronger than natural density. In view of Theorem 5.1, it is natural to introduce the local URI-density of A ⊂ Z + as the limit, when k → ∞, of P(N k ∈ A) (whenever the limit exists). The purpose of this section is to prove that local URI-density is stronger than natural density: Theorem 5.2. If A ⊂ Z + possesses a natural density, then the local URI-density of A exists and coincides with its natural density.

In fact, local URI-density turns out to be strictly stronger than natural density, since Theorem 5.1 proves the existence of sets A without natural densities but for which the local URI-density exists.

The proof of Theorem 5.2 is based on the following lemma. Lemma 5.3. Let (P k ) k≥1 be a sequence of probability distributions on Z + satisfying: For any k ≥ 1, there exists n k such that

• The map n → P k (n) is non-decreasing on {1, . . . , n k } and non-increasing on {n k , n k + 1, . . . }; • For any integer m ≥ 1,

lim k→∞ P k (mn k ) P k (n k ) = 1.
Then, if A ⊂ Z + possesses a natural density α,

lim k→∞ P k (A) = α. mn k n k min (t) n k max (t) t θP k (n k ) n k P k (n) n Figure 1. Profile of P k
Proof. Let us fix ε > 0. Let θ ∈]0, 1[, close enough to 1 so that (1θ)/θ < ε. Let m ∈ Z + be such that m > 1/ε and such that, for any n ≥ m,

1 n n i=1 ½ A (i) ∈ ]α -ε, α + ε[.
We choose k large enough such that (18)

P k (mn k ) P k (n k ) > θ.
By a Fubini argument, we can write P k (A) as

(19) P k (A) = P k (n k ) 0 |{n ∈ A : P k (n) > t}| dt.
We split the integral into two terms

I 1 := θP k (n k ) 0 |{n ∈ A : P k (n) > t}| dt and I 2 := P k (n k ) θP k (n k ) |{n ∈ A : P k (n) > t}| dt.
Observe that

θP k (n k ) |{n ∈ Z + : P k (n) > θP k (n k )}| ≤ P k (|{n ∈ Z + : P k (n) > θP k (n k )}|) ≤ 1.
Therefore

I 2 ≤ (1 -θ)P k (n k ) |{n ∈ Z + : P k (n) > θP k (n k )}| ≤ 1 -θ θ < ε.
Let us turn to the estimation of I 1 . By the hypothesis on the variations of P k (n), for any 0 < t < P k (n k ), there exist n k min (t) ≤ n k ≤ n k max (t) such that {n : P k (n) ≥ t} = n k min (t), . . . , n k max (t) . (See Figure 5.2.) We can rewrite I 1 as

θP k (n k ) 0 n k max (t) -n k min (t) + 1 ϕ(t) dt,
where

ϕ(t) := 1 n k max (t) -n k min (t) + 1 n k max (t) i=n k min (t) ½ A (i).
We prove that, for 0 < t < θP k (n k ), ϕ(t) is close to the natural density of A: By (18), for any 0 < t < θP k (n k ), we have n k max (t) > mn k . Thus

1 ≤ n k max (t) n k max (t) -n k min (t) + 1 ≤ m m -1 ≤ 1 1 -ε , and
1 n k max (t) n k min (t)-1 n=1 ½ A (n) ≤ n k min (t) -1 n k max (t) ≤ 1 m < ε. Since n k max (t) > m, 1 n k max (t) n k max (t) n=1 ½ A (n) ∈ ]α -ε, α + ε[, It follows that for 0 < t < θP k (n k ), α -2ε < ϕ(t) < (α + ε)/(1 -ε).
Hence we get the following estimation:

(α -2ε)I 3 < I 1 < α + ε 1 -ε I 3 ,
where

I 3 := θP k (n k ) 0 n k max (t) -n k min (t) + 1 dt.
Using (19) with A = Z + , we get

P k (n k ) 0 n k max (t) -n k min (t) + 1 dt = P k (Z + ) = 1,
and by the same argument as for the estimation of I 2 , we have

P k (n k ) θP k (n k ) n k max (t) -n k min (t) + 1 dt < 1 -θ θ < ε. hence 1 -ε < I 3 ≤ 1.
Observe that the second condition in the lemma is crucial. Indeed, the sequence of binomial distributions of parameter p ∈]0, 1[ defined by

B k (n) := k n p n (1 -p) k-n , 0 ≤ n ≤ k
satisfies the first assumption of the lemma. However there exists a set A possessing a natural density, but for which B k (A) fails to converge to this natural density as k → ∞ (see [START_REF] Diaconis | Weak and strong averages in probability and the theory of numbers[END_REF], Theorem 3 page 25).

Proof of Theorem 5.2. Defining P k (n) := P(N k = n), we have to check the hypotheses of Lemma 5.3. We start by establishing an induction formula for P k (n).

Recall that P(N k = n) = 0 if n < k, and that P(N

1 = n) = ½ n=1 . For 2 ≤ k ≤ n,
we decompose

P(N k = n) = P(N k = n, N k-1 = n -1) + P(N k = n, N k-1 ≤ n -2).
The first term in the RHS is (1/n)P(N k-1 = n-1). The second term can be written as

P(N k = n, N k-1 < n -1) = 1 n 1 - 1 n -1 P(|E n-2 | = k -1) = n -2 n 1 n -1 P(|E n-2 | = k -1) = n -2 n P(N k = n -1).
This yields, for any 2

≤ k ≤ n, (20)
P k (n) = n -2 n P k (n -1) + 1 n P k-1 (n -1).
For 2 ≤ k ≤ n -1, dividing by P k (n -1) we obtain (21)

P k (n) P k (n -1) = 1 + 1 n f n (k) -2 ,
where

f n (k) := P k-1 (n -1) P k (n -1)
.

We prove by induction on n ≥ 4 that k ∈ {2, . . . , n-1} → f n (k) is a non-decreasing function. Observe that f 4 (2) = 0 and f 4 (3) = 1, hence f 4 is non-decreasing. Assume that f n-1 is non-decreasing for some n ≥ 5. Using (20), we get for 3

≤ k ≤ n -1 f n (k) = n -3 n -1 P k-1 (n -2) + 1 n -1 P k-2 (n -2) n -3 n -1 P k (n -2) + 1 n -1 P k-1 (n -2)
.

If 3 ≤ k ≤ n -2, we get (22) f n (k) = n -3 + f n-1 (k -1) n -3 f n-1 (k) + 1
.

Hence, by induction, f n is non-decreasing on {2, . . . , n-2}. Moreover, for k = n-1, since P n-1 (n -2) = 0,

f n (n -1) = n -3 + f n-1 (n -2) ≥ n -3 + f n-1 (n -3) ≥ f n (n -2).
Observe also that f n (2) = 0 and f n (n -1) = (n -3)(n -2)/2 for all n ≥ 4. Hence, for all n ≥ 5, there exists an integer k n such that f

n (k) ≤ 2 for 2 ≤ k ≤ k n and f n (k) > 2 for k > k n . Since f n-1 (k -1) ≤ f n-1 (k) for any 3 ≤ k ≤ n -2, we get by (22) that f n (k) ≤ f n-1 (k), which proves that n → k n is non-decreasing.
For any fixed k ≥ 3, let n k be the smallest integer n such that k n ≥ k. By (21), n → P k (n) is non-decreasing upto n k and non-increasing after n k . Note that n k exists, otherwise n → P k (n) would be non-decreasing, which is obviously impossible. Therefore, the first hypothesis of Lemma 5.3 is satisfied.

For all k ≥ 3, observe that n k is characterized by the following:

(23) f n k (k) ≤ 2, and f n k -1 (k) > 2.
To check that (P k ) satisfies the second hypothesis, we need precise estimations of f n (k). We start by establishing a formula for P k (n). Observe that for all 1 ≤ j < n, P(N k = n|N k-1 = j) is equal to

P(j + 1 / ∈ E, . . . n -1 / ∈ E, n ∈ E) = j j + 1 . . . n -2 n -1 1 n = j n(n -1)
.

Hence, by conditioning, P k (n) = P(N k = n) can be rewritten as

2≤j2<•••<j k-1 ≤n-1 P(N k = n|N k-1 = j k-1 )P(N k-1 = j k-1 |N k-2 = j k-2 ) . . . P(N 3 = j 3 |N 2 = j 2 )P(N 2 = j 2 )
which yields

P k (n) = 2≤j2<•••<j k-1 ≤n-1 j k-1 n(n -1) j k-2 j k-1 (j k-1 -1) . . . j 2 j 3 (j 3 -1) 1 j 2 (j 2 -1) = 1 n(n -1) 2≤j2<•••<j k-1 ≤n-1 1 j k-1 -1 . . . 1 j 3 -1 1 j 2 -1 = 1 n(n -1) 1≤j2<•••<j k-1 ≤n-2 1 j 2 j 3 . . . j k-1
.

We use this formula to estimate the denominator in the definition of f n (k):

P k (n -1) = 1 (n -1)(n -2) 1≤j2<•••<j k-1 ≤n-3 1 j 2 j 3 . . . j k-1 = 1 (n -1)(n -2) 1≤j2<•••<j k-2 ≤n-3 1 j 2 j 3 . . . j k-2 C(j 2 , . . . , j k-2 ), where C(j 2 , . . . , j k-2 ) := 1 k -2 1≤j≤n-3 j / ∈{j 2 ,...,j k-2 } 1 j . Observe that 1 k -2 k-2≤j≤n-3 1 j ≤ C(j 2 , . . . , j k-2 ) ≤ 1 k -2 1≤j≤n-3 1 j ,
which gives the following estimation

(24) k -2 1≤j≤n-3 1 j ≤ f n (k) = P k-1 (n -1) P k (n -1) ≤ k -2 k-2≤j≤n-3 1 j • Let m ≥ 1. Applying this estimation to f n k -1 and f mn k , we get 0 ≤ f n k -1 (k) -f mn k (k) ≤ k -2 k-2≤j≤n k -4 1 j - k -2 1≤j≤mn k -3 1 j

•

The RHS of the above inequality can be written as a product AB, where

A := k -2 k-2≤j≤n k -4 1 j
and B :=

1≤j≤mn k -3 1 j -k-2≤j≤n k -4 1 j 1≤j≤mn k -3 1 j
• By ( 24) and ( 23), we get

(25) k -2 1≤j≤n k -3 1 j ≤ 2,
which ensures that A is bounded (say, by 4). Moreover, an easy computation shows that B ∼ ln k ln n k , which goes to 0 as k → ∞ by (25). Recalling that f n k -1 (k) > 2 by (23), the above estimations prove the following property: For any ε > 0, for k large enough, f n (k) > 2ε for each n k ≤ n ≤ mn k . For such k, we get by (21)

P k (mn k ) P k (n k ) = mn k n=n k +1 P k (n) P k (n -1) = mn k n=n k +1 1 + 1 n (f n (k) -2) ≥ 1 - ε n k (m-1)n k .
On the other hand, we know that P k (n k ) ≥ P k (mn k ), which proves that

lim k→∞ P k (mn k ) P k (n k ) = 1.
6. Open problems and discussion 6.1. Connection with other densities. We conjecture that the existence of local URI-density implies the existence of URI-density (and, in this case, that both coincide). It is not clear either whether local URI-density is equivalent to the H ∞ -density used by Flehinger. However we can provide the following interpretation of the H ∞ -density in terms of our URI-set. Recall that in the proof of Theorem 3.3, we ordered the elements of E n = E ∩ {1, . . . , n} backwards: .

E n = Y (n) 1 > Y (n) 2 > . . . > Y (n) |En| = 1 .
It remains to prove that, for any fixed k, P(D n k ) → 1 as n → ∞. Fix ǫ > 0 and choose δ > 0 such that (1 -3δ) k > 1ǫ. Observe that whenever n > 1/δ, the proportion of integers i ∈ {1, . . . , n} such that δ < i/n < 1δ is larger than 1 -3δ. Now, if n > 1/δ k , we have

P Y (n) 1 n ∈]δ, 1 -δ[, Y (n) 2 Y
(n) 1 ∈]δ, 1 -δ[, . . . , Y (n) k Y (n) k-1 ∈]δ, 1 -δ[ > (1 -3δ) k > 1 -ǫ.
Hence, P(D n k ) > 1ǫ. 6.2. Conditional URI-density. Let P be a subset of Z + with p∈P 1/p = ∞, so that the cardinality of P ∩ E be almost surely infinite. We have two ways to define the URI-density of A conditioned on P . First, by averaging over P ∩ E, and consider (whenever it exists) the almost-sure limit of Question: are these two definitions equivalent? 6.3. Asymptotic independence of successive elements in E. Another question concerning the URI-set E is the following: consider A, B ⊂ Z + , and assume that for both A and B we can define the density d(A) and d(B) (these could be the natural densities, the URI-densities, the H ∞ -densities or maybe some other notions of densities). Under which condition do we have

[ 1 ,

 1 b[ the mantissa in base b of a positive real number x and by µ B b the probability distribution over [1, b[ defined by µ B b ([1, t]) := log b t for all 1 ≤ t < b.

  which proves the claim. It follows by Lemma 4.2 that the sequence log b1 (Y k ) mod 1, . . . , log b ℓ (Y k ) mod 1 is uniformly distributed in (R/Z) ℓ , and the same is true if we replace Y k by X k . The end of the proof goes with similar arguments as for Theorem 3.5. Cassels-Schmidt-Benford sequences. Given two multiplicatively independent positive integers b 1 and b 2 , Bodo Volkmann defined a Cassels-Schmidt number of type (b 1 , b 2 ) as a number which is normal in base b 1 but not in base b 2 . By analogy, he also proposed to define a Cassels-Schmidt-Benford (CSB) sequence of type (b 1 , b 2 ) as a sequence of positive numbers (x k ) which follows Benford's law with respect to base b 1 but not with respect to base b 2 . It turns out that it is far easier to find explicit examples of CSB sequences. Indeed, take x k := (b 2 ) k , then (x k ) does certainly not follow Benford's law with respect to base b 2 . But since ln b 2 / ln b 1 is not a rational number, the sequence log b1 x k mod 1 = (k ln b 2 / ln b 1 mod 1) is uniformly distributed in [0, 1]. Hence (x k ) follows Benford's law with respect to base b 1 . As an application of Theorem 4.1 in the case ℓ = 2, we get that in almost every URI-set, we can find a CSB sequence (x k ) of type (b 1 , b 2 ), such that the sequence (M b2 (x k )) follow any probability distribution prescribed in advance on [1, b 2 [.

  follows by (3) and Remark 3.4 that, for any Borel set S ⊂ [1, 10[, (14) ∀k ≥ 1, P M (X m k ) ∈ Sµ B (S) ≤ 9 10 k .

Proposition 6 . 1 (For½

 61 Stochastic interpretation of H ∞ -density). ForA ⊂ Z + , A has H ∞ -density α if and only if lim For each n ≥ 1, we introduce a non-increasing sequence of random integers ( Y (n) i ) i≥1 , with the following distribution: in {1, . . . , n}, • conditionally to Y A ⊂ Z + , we write P Y (n) k ∈ A as 1≤y k ≤y k-1 ≤•••≤y1≤n ½ A (y k )P Y (n) k = y k | Y A (y k ).We recognize P k n used in the definition of the H ∞ -density (see Section 3.3). Now,

n k=1 ½

 k=1 P ∩A (N k ) n k=1 ½ P (N k ).Second, by numbering the elements ofP = {p 1 < p 2 < • • • < p n < • • • } and averaging over the random subset of P {p N1 < p N2 < • • • < p n k < • • • },that is by considering (whenever it exists) the almost-sure limit of1 n n k=1 ½ A (p N k ).

½

  A (N k )½ B (N k+1 ) a.s. ----→ n→∞ d(A)d(B) ?

  2.1. Flehinger's theorem through Markov chain. We introduce a homogeneous Markov chain (M k ) k≥0 taking values in [1, 10[, defined by its initial value M 0 (which can be deterministic or random) and the following transition probability: for any Borel set S ⊂[1, 10[,

	(2) P M k+1 ∈ S M where U is a uniform random variable in [0, 1].
	Let us denote by µ B the probability distribution on [1, 10[ given by Benford's
	law: For any 1 ≤ t < 10,
	µ B ([1, t]) := log 10 t.
	It is proved in [8] by a standard coupling argument that
	• µ B is the only probability distribution on [1, 10[ which is invariant under
	the probability transition (2);
	• Whatever choice we make for the initial condition M 0 , we have for any
	Borel set S ⊂ [1, 10[ and for all k ≥ 1
	(3)

k = a := P M (aU ) ∈ S ,

  the form 10 p (which, again by Borel-Cantelli, happens almost surely for only finitely many k's). Equivalence with log-density. To deal with the problem of non-existence of natural densities, several alternative densities have been introduced. Flehinger's theorem amounts to considering the so-called H ∞ -density, obtained by iteration of Cesaro averages: A subset A of Z + has H ∞ -density α if

	3.3. lim k→∞	lim sup n→∞	P k n = lim k→∞	lim inf n→∞
	It follows from (10) that
	1 n	n k=1	½ M (N k )<t	a.s. ----→ n→∞	log 10 t.

We conjecture that it is true when both d(A) and d(B) are natural densities. However this is certainly not true for all A and B with URI-densities: As a counterexample, consider the set A of integers with leading digit 1 and the set B of integers with leading digit 9. But what happens if for example d(B) is the natural density whereas d(A) is only the URI-density? 6.4. Iterated URI-density. Diaconis defined in [START_REF] Diaconis | Weak and strong averages in probability and the theory of numbers[END_REF] the iterated log-density: For a subset A of Z + , set

and inductively for all ℓ ≥ 2:

The set A has ℓ-th log-density α if the limit of the above exists as n → ∞ and is equal to α. In fact, this notion does not yield a new density, since Diaconis proved that A has an ℓ-th log-density if and only if A has a log-density (and then, of course, both coincide). Then he proposed to define the L ∞ -density, which extends the logdensity in much the same way as H ∞ -density extends natural density: Consider lim ℓ→∞ lim inf n→∞ L(A, n, ℓ) and lim ℓ→∞ lim sup n→∞ L(A, n, ℓ). If the two limits are equal, call their common value the L ∞ -density of A. As shown in [START_REF]Examples for the theory of infinite iteration of summability methods[END_REF], L ∞density is strictly stronger than log-density.

It is natural in this context to study iterations of the URI-density, which can be defined as follows: Let E (ℓ) ℓ≥1 be a sequence of independent URI-sets, and denote the random elements of E (ℓ) by N (ℓ)

We define the URI-density of order ℓ in the same way, averaging along the subsequence N

(1) N [START_REF] Diaconis | Weak and strong averages in probability and the theory of numbers[END_REF] . . .

We can also introduce the infinite iteration of the URI-method, considering almost-sure limsup and liminf in the above expressions, and see if they converge to the same limit as ℓ → ∞.

Although we have shown that URI-density and log-density coincide, it is not obvious if there are connections between iterated URI-density and iterated logdensity. Can URI-densities of finite order ℓ be strictly stronger than URI-density? Can we compare the infinite iteration of both methods?