Erik Vanhoutte
email: erik.vanhoutte@univ-amu.fr

Stefano Mafrica
email: stefano.mafrica@mpsa.com

Franck Ruffier
email: franck.ruffier@univ-amu.fr

Reinoud J Bootsma
email: reinoud.bootsma@univ-amu.fr

Julien Serres
email: julien.serres@univ-amu.fr

Felipe Gonzalez

Antonios Tsourdos

Time-of-Travel Methods for Measuring Optical Flow on Board a Micro Flying Robot

Keywords: optic flow sensor, sense and avoid, VLSI retina, micro air vehicle (MAV), bionics, bio-inspired robotics, biorobotics

For use in autonomous micro air vehicles, visual sensors must not only be small, lightweight and insensitive to light variations; on-board autopilots also require fast and accurate optical flow measurements over a wide range of speeds. Using an auto-adaptive bio-inspired Michaelis-Menten Auto-adaptive Pixel (M 2 APix) analog silicon retina, in this article, we present comparative tests of two optical flow calculation algorithms operating under lighting conditions from 6 × 10 -7 to 1.6 × 10 -2 W•cm -2 (i.e., from 0.2 to 12,000 lux for human vision). Contrast "time of travel" between two adjacent light-sensitive pixels was determined by thresholding and by cross-correlating the two pixels' signals, with measurement frequency up to 5 kHz for the 10 local motion sensors of the M 2 APix sensor. While both algorithms adequately measured optical flow between 25 • /s and 1000 • /s, thresholding gave rise to a lower precision, especially due to a larger number of outliers at higher speeds. Compared to thresholding, cross-correlation also allowed for a higher rate of optical flow output (99 Hz and 1195 Hz, respectively) but required substantially more computational resources.

Introduction

To detect and avoid obstacles in unpredictable environments, flying insects rely heavily on optical flow (OF) [START_REF] Srinivasan | Honeybees as a model for the study of visually guided flight, navigation, and biologically inspired robotics[END_REF], defined as the vector field of angular velocities of contrasted points, edges or surfaces resulting from the relative motion between the observer and the surrounding objects [START_REF] Gibson | The perception of the visual world[END_REF][START_REF] Nakayama | Optical velocity patterns, velocity-sensitive neurons, and space perception: A hypothesis[END_REF]. Their OF-based strategies therefore provide inspiration for the development of smart autopilots for micro-air-vehicles (MAVs) [START_REF] Srinivasan | Honeybees as a model for the study of visually guided flight, navigation, and biologically inspired robotics[END_REF][START_REF] Franceschini | Optic Flow Based Visual Guidance: From Flying Insects to Miniature Aerial Vehicles[END_REF][START_REF] Serres | Optic Flow-Based Robotics[END_REF] and smart artificial retinas [START_REF] Moeckel | Motion detection chips for robotic platforms[END_REF][START_REF] Expert | A mouse sensor and a 2-pixel motion sensor exposed to continuous illuminance changes[END_REF][START_REF] Roubieu | Two-Directional 1-g Visual Motion Sensor Inspired by the Fly's Eye[END_REF][START_REF] Floreano | Miniature curved artificial compound eyes[END_REF][START_REF] Song | Digital cameras with designs inspired by the arthropod eye[END_REF][START_REF] Duhamel | Biologically inspired optical-flow sensing for altitude control of flapping-wing microrobots[END_REF][START_REF] Mafrica | A bio-inspired analog silicon retina with Michaelis-Menten auto-adaptive pixels sensitive to small and large changes in light[END_REF].

Insect-sized MAVs are increasingly becoming a reality [START_REF] Duhamel | Altitude feedback control of a flapping-wing microrobot using an on-board biologically inspired optical flow sensor[END_REF][START_REF] Kushleyev | Towards a swarm of agile micro quadrotors[END_REF][START_REF] Ma | Controlled flight of a biologically inspired, insect-scale robot[END_REF][START_REF] Dunkley | Visual-inertial navigation for a camera-equipped 25g nano-quadrotor[END_REF][START_REF] Moore | Autonomous MAV guidance with a lightweight omnidirectional vision sensor[END_REF][START_REF] Floreano | Science, technology and the future of small autonomous drones[END_REF] and have to be fitted in the future with sensors and flight control devices enabling them to perform all kinds of aerial maneuvers, including ground and obstacle avoidance, terrain-following and landing.

The fitted sensors should not only be non-emissive to allow an MAV to save energy resources and to be stealth in flight, but must also guarantee a high refresh rate because GPS (GPS stands for the Global Positioning System) signals are limited in both spatial (∼1 m) and temporal resolution (∼7 Hz). Micro cameras have been proposed for MAV applications (e.g., a PAL (PAL stands for Phase Alternating Line) -camera with 720 × 576 pixels at 25 fps [START_REF] Dunkley | Visual-inertial navigation for a camera-equipped 25g nano-quadrotor[END_REF] or a CMOS (CMOS stands for Complementary Metal-Oxide Semiconductor) camera with 752 × 480 pixels at 80 fps [START_REF] Liu | Computationally efficient visual-inertial sensor fusion for Global Positioning System-denied navigation on a small quadrotor[END_REF]), including visual-based simultaneous localization and mapping algorithms (SLAM).

A few cameras deliver fast and reliable OF measurement, but still remain bulky: (i) a PIX4FLOW CMOS camera with 4× binning and 188 × 120 pixels delivers only one OF vector-direction and magnitude-at 250 fps [START_REF] Honegger | An open source and open hardware embedded metric optical flow cmos camera for indoor and outdoor applications[END_REF]; and (ii) the event-based camera called DAVIS (DAVIS stands for Dynamic and Active-pixel VIsion Sensor) with 240 × 180 pixels delivers accurate OF measurement asynchronously [START_REF] Brandli | A 240 × 180 130 dB 3 µs latency global shutter spatiotemporal vision sensor[END_REF][START_REF] Rueckauer | Evaluation of event-based algorithms for optical flow with ground-truth from inertial measurement sensor[END_REF]. The drawbacks of CMOS cameras moreover include their poor robustness to various lighting conditions (e.g., ∼600 lux in [START_REF] Honegger | An open source and open hardware embedded metric optical flow cmos camera for indoor and outdoor applications[END_REF]) and their high (pixel number based) demands on computational resources.

It is essential for MAV applications to use lightweight, low-power consumption and high refresh rate sensors to comply with fast dynamics and to stay reactive in unpredictable environments. Hovering requires lower refresh rates than obstacle avoidance maneuvers while flying at high speed. Inertial measurement units (IMUs) are used to stabilize an MAV in flight or simply to hover at a refresh rate up to 500 Hz. Recently, both an IMU sampled at 500 Hz and OF measurements sampled at 25 Hz were fused using an extended Kalman filter to make an MAV hover robustly [START_REF] Briod | A method for ego-motion estimation in micro-hovering platforms flying in very cluttered environments[END_REF]. The authors' idea was to use the OF direction, but not its scale to correct for inertial sensor drift in particular during changes of direction and to enhance the accuracy of a hovering positioning: this strategy helps in updating both an attitude and a positioning control outputs at 100 Hz [START_REF] Briod | A method for ego-motion estimation in micro-hovering platforms flying in very cluttered environments[END_REF]. Lately, a 40-gram pocket drone avoided obstacles on the basis of OF computed at 20 Hz [START_REF] Mcguire | Efficient Optical flow and Stereo Vision for Velocity Estimation and Obstacle Avoidance on an Autonomous Pocket Drone[END_REF]. The refresh rate of the visual processing appears to be the key parameter to make an MAV achieve "aggressive" maneuvers based on on-board vision (as demonstrated in [START_REF] Falanga | Aggressive Quadrotor Flight through Narrow Gaps with Onboard Sensing and Computing[END_REF]) with sufficient reactive abilities in unpredictable environments. More generally, the criteria for the evaluation of the potential of OF sensors for MAV applications include:

• robustness to light level variations, defined by the number of irradiance decades in which the visual sensor can operate; • range of OF angular speeds (or magnitudes) covered, defined by the minimum and maximum values measured; • accuracy and precision, defined by systematic errors and coefficients of variation; • output refresh rate, defined by the instantaneous output frequency.

A recent OF sensor was based on the M 2 APix (M 2 APix stands for Michaelis-Menten auto-adaptive pixel) retina that can auto-adapt in a seven-decade lighting range and responds appropriately to step changes up to ±3 decades [START_REF] Mafrica | Minimalistic optic flow sensors applied to indoor and outdoor visual guidance and odometry on a car-like robot[END_REF]. The pixels do not saturate thanks to the (normalization based) intrinsic properties of the Michaelis-Menten equation [START_REF] Normann | The effects of background illumination on the photoresponses of red and green cones[END_REF]. A comparison of the characteristics of auto-adaptive Michaelis-Menten and Delbrück pixels under identical lighting conditions demonstrated better performance of the Michaelis-Menten pixels in terms of dynamic sensitivity and minimum contrast detection [START_REF] Mafrica | A bio-inspired analog silicon retina with Michaelis-Menten auto-adaptive pixels sensitive to small and large changes in light[END_REF].

A single local motion sensor (LMS) fitted with two auto-adaptive pixels has been demonstrated [START_REF] Expert | A mouse sensor and a 2-pixel motion sensor exposed to continuous illuminance changes[END_REF] to allow measuring an OF range from 50 to 350 °/s at relatively constant output refresh rate (∼5-7 Hz) despite variations in lighting conditions from ∼50 lux to 10,000 lux. With a larger number of LMSs, output refresh rate increased (64 Hz with 5 LMSs) in an OF range from 25 to 350 °/s (i.e., 1.1-decade) at a constant natural lighting ∼1500 lux [START_REF] Roubieu | Two-Directional 1-g Visual Motion Sensor Inspired by the Fly's Eye[END_REF]. A similar OF range from 50 to 250 °/s (i.e., 0.7-decade) was measured with a semi-panoramic artificial eye, called CurvACE [START_REF] Floreano | Miniature curved artificial compound eyes[END_REF][START_REF] Viollet | Hardware architecture and cutting-edge assembly process of a tiny curved compound eye[END_REF]. The OF range can also be adjusted by tuning the inter-pixel angle, as done for outdoor flights in [START_REF] Sabiron | Low-speed optic-flow sensor onboard an unmanned helicopter flying outside over fields[END_REF]. An ad hoc interpolation-based "time-of-travel" algorithm relying on thresholding was embedded into a small 40 MIPS (MIPS stands for Million Instructions Per Second) dsPIC ® (Microchip Technology Inc., Chandler, AZ, USA) microcontroller to allow a trade-off between OF range [25 °/s; 350 °/s], accuracy (σ error < 15 °/s), the sample rate of the visual signals (200 Hz), computational resources (10.5% on average of the processing time available) and the OF refresh rate (95 Hz with 5 LMSs) [START_REF] Expert | Interpolation based "time of travel" scheme in a Visual Motion Sensor using a small 2D retina[END_REF]. With this limited OF range, the robot's operating OF was adjusted to the middle of the range: 125 °/s in [START_REF] Roubieu | A biomimetic vision-based hovercraft accounts for bees' complex behaviour in various corridors[END_REF], 160 °/s in [START_REF] Roubieu | A fully-autonomous hovercraft inspired by bees: Wall following and speed control in straight and tapered corridors[END_REF], or even 200 °/s in [START_REF] Expert | Flying over uneven moving terrain based on optic-flow cues without any need for reference frames or accelerometers[END_REF] (which was by the way quite high in comparison with related works: ∼60-70 °/s in [START_REF] Moore | Autonomous MAV guidance with a lightweight omnidirectional vision sensor[END_REF][START_REF] Zufferey | Fly-inspired visual steering of an ultralight indoor aircraft[END_REF][START_REF] Zufferey | A 10-gram vision-based flying robot[END_REF][START_REF] Beyeler | Vision-based control of near-obstacle flight[END_REF]).

In this paper, we compare two different methods for the calculation of OF. Both methods are based on time-of-travel algorithms, relying either on thresholding or on cross-correlation of adjacent bandpass-filtered visual signals. The two algorithms are applied on the same visual signals provided by the M 2 APix retina when operating under different combinations of the following conditions:

• in an OF range from 25 °/s to 1000 °/s; • under irradiance conditions varying from 6 × 10 -7 W•cm -2 to 1.6 × 10 -2 W•cm -2 ; • with sampling rates between 100 Hz and 1 kHz; • in real flight when fitted onto a 350-gram MAV.

It will be shown that at low sampling rate (F e < 500 Hz), the cross-correlation method is more precise than the thresholding method over a wide range of OF speeds. Requiring much less computational resources, the thresholding method functions adequately at a high sampling rate (F e close to 1 kHz), with the OF range covered depending mainly on the inter pixel angle ∆ϕ.

In Section 2, optics and front-end pixels of the M 2 APix sensor will be briefly introduced with respect to previous work based on the M 2 APix sensor. In Section 3, the two methods of OF computation based on time-of-travel algorithms will be described in detail. In Section 4, both time-of-travel algorithms applied on the same M 2 APix sensor visual signals will be tested in an OF range from 25 °/s to 175 °/s at different levels of lighting, then in an OF range from 25 °/s to 1000 °/s at a constant lighting by using an experimental setup based on a moving pattern. In Section 5, both time-of-travel algorithms with the M 2 APix sensor will be tested on-board a 350-gram quadrotor during flights over a flat ground textured with a natural scene. Both time-of-travel methods will be compared in terms of computational resources, accuracy, precision and output refresh rate. In Section 6, the benefits of each method will be discussed, leading to the conclusion that time-of-travel algorithms with a M 2 APix sensor can measure a wide OF range at high refresh output frequency despite large lighting variations in real flight.

Optics and Front-End Pixels of the M 2 APix Sensor

The M 2 APix sensor is composed of two groups of light-sensitive pixels: 12 Michaelis-Menten pixels and 12 additional Delbrück pixels [START_REF] Mafrica | A bio-inspired analog silicon retina with Michaelis-Menten auto-adaptive pixels sensitive to small and large changes in light[END_REF] (Figure 1a,b). Each group is composed of two rows of six pixels with rows offset by half the inter-pixel distance (Figure 1c). In this paper, we only used the 12 Michaelis-Menten pixels, sampled at the frequency F e (from 100 Hz to 1 kHz, depending on the experimental setting) by an Overo AirSTORM computer-on-module (COM) (Gumstix, Redwood City, CA, USA) featuring a 1-GHz CPU DM3703 processor (Texas Instruments, Dallas, TX, USA) comprising an ARM Cortex-A8 architecture. From [START_REF] Mafrica | A bio-inspired analog silicon retina with Michaelis-Menten auto-adaptive pixels sensitive to small and large changes in light[END_REF].

The M 2 APix sensor is equipped with a lens taken from a Raspberry Pi camera (focal length = 2 mm). The lens is deliberately defocused (distance between the lens and the retina = 1 mm) to provide a bio-inspired low-pass spatial filter [START_REF] Land | Visual Acuity in Insects[END_REF]. The optical configuration can be described by the inter-pixel angle ∆ϕ = 4.3 • and the full width at half height of the Gaussian sensitivity, the acceptance angle ∆ρ = 3.0 • , thereby respecting the ∆ρ ∆ϕ < 2 rule for proper visual sampling established for diurnal insects. In this configuration, the M 2 APix sensor has a field of view of 32.8 • × 13.4 • in the horizontal and vertical plane, respectively.

In most animals' visual systems, the photoreceptors have been observed to auto-adapt their response to the ambient light level [START_REF] Normann | The effects of background illumination on the photoresponses of red and green cones[END_REF]. The M 2 APix pixels are also auto-adaptive, as their responses are normalized by a signal that depends on the low-pass filtered average of all of the 12 photoreceptors' currents, according to the Michaelis-Menten equation [START_REF] Mafrica | A bio-inspired analog silicon retina with Michaelis-Menten auto-adaptive pixels sensitive to small and large changes in light[END_REF][START_REF] Michaelis | Die kinetik der invertinwirkung[END_REF]. The adaptation time constant of the pixels' response was set here at 0.5 s by using an external capacitor of 47 nF (see [START_REF] Mafrica | A bio-inspired analog silicon retina with Michaelis-Menten auto-adaptive pixels sensitive to small and large changes in light[END_REF] for more details). This allows the amplitude of the output signals to remain within the same constant range notwithstanding variations in light level and, therefore, makes it possible for the pixels to operate in high-dynamic-range lighting conditions.

Optical Flow Computed by Time-Of-Travel-Based Algorithms

Time-of-travel algorithms compute the magnitude of OF (Equation (1)) by measuring the time delay ∆t between the signals of two adjacent photoreceptors (here Michaelis-Menten pixels) constituting a local motion sensor (LMS):

OF = ∆ϕ ∆t (1)
In the work presented here, we used 10 LMSs (two rows of five), each computing the 1D OF in the direction parallel to the rows of pixels, as shown in Figure 2. The optics geometry of the first two-pixel LMS (composed of the red and blue pixels) is described in Figure 3a.

Time-Of-Travel Based on Signal Thresholding

This time-of-travel algorithm is based on opto-electrophysiological studies of the fly's lobula plate tangential cells (H1-neuron) [START_REF] Franceschini | Directionally Selective Motion Detection by Insect Neurons[END_REF] and is described in detail in [START_REF] Franceschini | Optic Flow Based Visual Guidance: From Flying Insects to Miniature Aerial Vehicles[END_REF][START_REF] Expert | A mouse sensor and a 2-pixel motion sensor exposed to continuous illuminance changes[END_REF][START_REF] Roubieu | Two-Directional 1-g Visual Motion Sensor Inspired by the Fly's Eye[END_REF][START_REF] Sabiron | Low-speed optic-flow sensor onboard an unmanned helicopter flying outside over fields[END_REF].

The algorithm is illustrated in Figure 3b and computes the "time-of-travel" ∆t by means of signal thresholding. To this end, the output of each pixel is first digitally band-pass filtered with cut-off frequencies of 3 Hz and 30 Hz. Two adjacent filtered signals are then thresholded (hysteresis thresholding) in order to measure the time delay ∆t between them: a counter is started when the first pixel's signal (the blue signal in Figure 3b) reaches a threshold, and it is stopped when the second pixel's signal (the red signal in Figure 3b) reaches the same threshold. Lastly, the OF magnitude is computed using Equation [START_REF] Srinivasan | Honeybees as a model for the study of visually guided flight, navigation, and biologically inspired robotics[END_REF].

The key parameter in this algorithm is the threshold value. A too large threshold value will allow detection only of high-amplitude pixel signals induced by strong contrasts and is therefore likely to neglect relevant signals. A too small threshold value will allow noise to cross the threshold and induce erroneous OF values. The thresholds' values (hysteresis thresholding in Figure 3b) were therefore adjusted manually with respect to the best trade-off between the maximum refresh rate of the LMSs' output and the minimum number of outliers.

In the experiments presented here, the high threshold was set at 5% of maximum signal amplitude (half the value of the high threshold for the low threshold for the hysteresis comparator; see Figure 3b), and all time delays falling outside a specific range (∆t from 4 to 717 ms, i.e., OF from 6 to 1145 • /s) were excluded and considered as outliers. The algorithm was run at the pixels' sampling frequency F e , which varied between 100 Hz and 1 kHz, giving a quantization time T e = 1 F e varying between 1 ms and 10 ms. The OF was then computed by using Equation ( 1), where ∆t is a multiple of T e . The OF resolution ∆OF of this method is therefore non-constant, as it is given by an inverse function of the time delay ∆t. 2 in [START_REF] Mafrica | Minimalistic optic flow sensors applied to indoor and outdoor visual guidance and odometry on a car-like robot[END_REF] for details). Note that for simplicity, the time-of-travel computation is shown here using the raw signals V ph1 and V ph2 , but it is actually computed after band-pass filtering these signals (see [START_REF] Franceschini | Optic Flow Based Visual Guidance: From Flying Insects to Miniature Aerial Vehicles[END_REF][START_REF] Expert | A mouse sensor and a 2-pixel motion sensor exposed to continuous illuminance changes[END_REF][START_REF] Roubieu | Two-Directional 1-g Visual Motion Sensor Inspired by the Fly's Eye[END_REF][START_REF] Mafrica | Minimalistic optic flow sensors applied to indoor and outdoor visual guidance and odometry on a car-like robot[END_REF][START_REF] Sabiron | Low-speed optic-flow sensor onboard an unmanned helicopter flying outside over fields[END_REF] for details).

LMS

Vph2 Vph1 Vph2 V p h 2 ( t -Δ t 1 ) V p h 2 ( t -Δ t 2 ) V p h 2 ( t -Δ t 3 ) V p h 2 ( t -Δ t 4 ) V p h 2 ( t -Δ

Time-of-Travel Based on Signals' Cross-Correlation

This algorithm was inspired by the correlation-based models [START_REF] Hassenstein | Systemtheoretische analyse der zeit-, reihenfolgen-und vorzeichenauswertung bei der bewegungsperzeption des rüsselkäfers chlorophanus[END_REF][START_REF] Albus | depth, and image flow[END_REF] and is based on signals' cross-correlation (Figure 3c), as presented in [START_REF] Mafrica | Minimalistic optic flow sensors applied to indoor and outdoor visual guidance and odometry on a car-like robot[END_REF]. First, the pixels' output signals are digitally filtered with the same band-pass filter used in the thresholding algorithm (see Section 3.1). Then, one of the two signals in each pair of adjacent pixels is delayed by several time delays ∆t i , and the Pearson cross-correlation coefficients i are computed between the delayed (V ph2 in Figure 3c) and non-delayed signals (V ph1 in Figure 3c) within a fixed time window T w . Lastly, the time delay ∆t m giving the maximum coefficient m is obtained and used to compute the OF using Equation (1) as long as m is greater than a given threshold thr . This threshold on the cross-correlation coefficients was set at 0.99 to avoid OF measurement errors due to signals mismatching [START_REF] Mafrica | Minimalistic optic flow sensors applied to indoor and outdoor visual guidance and odometry on a car-like robot[END_REF].

The size of the fixed time window T w is one of the key parameters and requires being tuned as function of the signals' bandwidth in order to have reliable cross-correlation coefficients ( i ). Since the signals' bandwidth is given by the band-pass filter (i.e. 3-30 Hz; see Section 3.1), T w was fixed here at a value of T w = 0.14 s so that 3 ≤ 1 T w ≤ 30. The cross-correlation time window was also defined as T w = N • T e , with N being the number of samples in the time window. Therefore, a sampling rate F e = 500 Hz entails N = 70, F e = 250 Hz entails N = 35, and so on.

To compute the cross-correlation coefficients i , the signal time window T w can be delayed by time steps ∆t i that are calculated from the desired OF values OF i using the inverse function given in Equation ( 1). This can be done by taking time steps ∆t i between two consecutive sampling steps (n T e and (n + 1) T e ) after linearly interpolating the sampled signals. Therefore, by choosing the desired OF values OF i using constant steps, the computed OF will as a result be quantized with a constant resolution ∆OF, contrary to the algorithm based on signal thresholding (see Section 3.1).

Measuring Optical Flow with a Moving Texture

Method

In the following experiments, the M 2 APix sensor was fixed perpendicularly to a moving pattern whose linear speed could be controlled (Figure 4). Lighting was provided by three sources: daylight coming through the top windows (intensity manipulated by means of blinds), fluorescent tubes attached to the ceiling (intensity manipulated by means of a dimmer) and an LED projector that could be positioned in front of the moving pattern for high luminosity conditions. The M 2 APix sensor was thus confronted with a pure translational wide-range OF (from 25°/s to 1000°/s) under various lighting conditions (from 6 × 10 -7 W•cm -2 to 1.6 × 10 -2 W•cm -2 , i.e., from 0.2 lux to 12,000 lux for human vision). The moving-pattern experiments were separated in two parts: in the first part, the sensor was placed at a distance d = 40 cm from the moving pattern, and both the pattern speed and the lighting conditions were varied, whereas in the second part, only the pattern speed was varied at long constant steps, and the sensor was placed closer to the moving pattern (d = 8 cm) to obtain high OF.

During the first part, each of the two algorithms was tested at five different irradiance levels. For each irradiance level, the pattern angular speed OF pattern varied over 10 periods as a sinus function with an angular frequency f = 0.2π rad/s and a magnitude speed V p varying from V pmin = 0.18 m•s -1 to V pmax = 1.22 m•s -1 (Equation ( 2)). With the M 2 APix sensor placed at a distance d = 40 cm (Figure 4), the moving pattern generated an OF range from 25°/s-175°/s (green lines in Figure 5b,c), according to the following equation:

OF pattern = 180 π V pmax -V pmin 2d • cos(2π f t + π) + V pmax + V pmin 2d (2)
During the second part, the level of lighting was fixed at an irradiance of 7 × 10 -4 W•cm -2 , and the two algorithms were used to compute OF in a wide OF range, from 25 °/s to 1000 °/s, by varying both the pattern speed (V p from 0 to 1.4 m•s -1 ) and the sensor's distance d (see Figure 4). In particular, to obtain an OF ranging from 25 °/s to 200 °/s, d was set at 40 cm, and V p varied from 0.18 m•s -1 to 1.34 m•s -1 ; whereas, for an OF ranging from 250 °/s to 1000 °/s, d was set at 8 cm, and

V p varied from 0.35 m•s -1 to 1.4 m•s -1 .
These two experiments allowed us to compare the precision and the output refresh rate of both time-of-travel algorithms when using an auto-adaptive retina, such as the M 2 APix sensor in a wide light-level range (five-decade) and in a wide OF range (1.6-decade). However, the results should be similar with any other auto-adaptive retina implementing a similar auto-adaptation process.

The two algorithms were run at the highest operating rate/frequency that still guaranteed a proper operation for future embedded MAV applications, taking into account the constraint to maintain 40% of CPU (CPU stands for central processing unit) load-free for MAV control and navigation tasks. Therefore, the time-of-travel algorithm based on signals thresholding was run at 1000 Hz, whereas the time-of-travel algorithm based on signals cross-correlation was run at 500 Hz due to the high CPU load of this method.

Results

Figure 5 shows the OF measurements and their instantaneous refresh rates obtained when using the thresholding method (Figure 5b) and the cross-correlation method (Figure 5c) at five different light levels (five colored columns in Figure 5a-c). The OF measurements corresponding to the five light levels shown in Figure 5a were actually obtained during five distinct tests (see the oblique black line segments separating each test at 5, 10, 15 and 20 s in Figure 5b,c). Furthermore, for each light level (each column/sinusoidal period in Figure 5b,c), the OF measurements of the 10 LMSs obtained during 10 sinusoidal periods of the moving pattern (see Section 4.1) were overlaid on only one sinusoidal period in order to show all of the measurements in one figure. The refresh rates given in the lower plots in Figure 5b,c were computed by low-pass filtering the average number of OF measurements obtained at every time step over the last three and the next three steps (six centered point low-pass digital averaging filter). The OF errors were computed as the difference between the OF measurements and their ground-truth values (blue dots and green lines in Figure 5b,c, respectively). The boxplots and the distributions of these errors shown in Figure 5d were computed with 1 °/s beams in the [150 °/s; 150 °/s] range.

The number of outliers is much lower for the cross-correlation method (Figure 5c) than for the thresholding method (Figure 5b). While both error distributions (Figure 5d) have a similar standard deviation, the boxplots reveal the differences in the number of outliers.

The thresholding method running at 1 kHz (Figure 5b) is very light in terms of CPU load (2.8% ± 0.8%) and provides a good OF accuracy and precision until an irradiance of 7 × 10 -6 W•cm -2 (blue column in Figure 5) despite more outliers than the cross-correlation method. However, the refresh rate obtained is relatively low (average of 97 Hz and maximum of 607 Hz for the full sensor). At very low irradiance, i.e., 6 × 10 -7 W•cm -2 , both the OF precision and the refresh rate are strongly deteriorated (see the violet column in Figure 5).

In contrast, the cross-correlation method running at 500 Hz gives a much lower number of outliers (boxplots in Figure 5d) and a much higher refresh rate (average of 2326 Hz and maximum of 5000 Hz for the full sensor), but it is very heavy in terms of CPU load (52.5% ± 2.2%). Similarly to the thresholding method, at very low light levels, i.e., 6 × 10 -7 W•cm -2 (violet column in Figure 5), the refresh rate is deteriorated, especially at very low OF values. With the cross-correlation method, the OF precision is however less deteriorated than with the thresholding method.

The cross-correlation method does not delay signals in constant steps, but with steps following an inverse function to obtain a constant OF-measurement resolution. The five levels of irradiance from 6 × 10 -7 W•cm -2 to 1.6 × 10 -2 W•cm -2 (i.e., from 0.2 lux to 12,000 lux for human vision) used during the first part of the experiments (see the experimental setup shown in Figure 4). (b,c) OF measurements obtained with the M 2 APix sensor using (b) the thresholding method (F e = 1 kHz) and (c) the cross-correlation method (F e = 500 Hz) for the 10 LMSs during 10 laps of the moving pattern. Tests were run separately for each of the five lighting conditions. Each blue dot represents one OF measurement of one LMS; however, the mean output refresh rate is the average of the number of OF measurement from the 10 LMSs in 1 s. The green line represents the OF ground-truth computed using Equation ( 2). (d) Error distribution of each time-of-travel method: the thresholding method is in red and the cross-correlation method in blue.

In the second experiment, various levels of OF were tested from 25 °/s to 1000 °/s with an OF step of 50 °/s (except for the first step) at a constant irradiance level (7 × 10 -4 W•cm -2 ). For each step, the OF mean µ, standard deviation σ and coefficient of variation σ µ were computed for both time-of-travel algorithms, as shown in Figure 6.

Inspection of Figure 6 reveals that the OF mean values µ are very close to the ground-truth values for both methods (blue circles in Figure 6). However, the standard deviation σ (blue vertical bars in Figure 6), and therefore, the coefficient of variation C V = σ µ (red curves in Figure 6), is higher when using the thresholding method (C V > 0.1 in Figure 6a) than the cross-correlation method (C V < 0.1 in Figure 6b), except for very low OF values, i.e., 25°/s, where the precision is nearly the same.

The increase in the coefficient of variation with the thresholding method at low angular speed could be explained by the auto-adaptive response of the M 2 APix sensor that respects a time constant of 0.5 s on average on the whole sensor. Consequently, the lower the angular speed, the higher the number of outliers, as the pixels' adaptation time constant is close to the signals' dynamics, therefore generating matching errors. If one of the two signals does not reach the hysteresis threshold, no OF is measured. This does not happen with the cross-correlation method because the OF is computed by considering a time window comprising N samples (see Section 3.2).

Nevertheless, even with the cross-correlation method, the coefficient of variation increases abruptly around 25°/s. If signals coming from the same contrast on two adjacent pixels are too different due to the auto-adaptation process, the minimum cross-correlation coefficient (set at i = 0.99) will not be reached, and no OF measurement will be generated. Moreover, at low speed, the signal dynamics is close to the fixed time window size set at T w = 0.14 s. Different signals may therefore appear similar in this time window, thereby generating matching errors.

Toward the high angular speeds, the delay ∆t between each photoreceptors' signal approaches the M 2 APix sensor sampling rate F e . The OF calculation therefore becomes less precise, especially for the thresholding method in which precision depends directly on the sampling rate F e .

In the experiments presented in this section, the M 2 APix sensor was fixed while a textured pattern was moving in front of it in a wide OF range and under controlled wide-range light levels. However, real flight conditions can generate disturbances due to vibrations or uncontrolled lighting conditions, which can change abruptly during the flight. Therefore, to test the two algorithms in real flight conditions, one M 2 APix sensor was mounted underneath a 350-gram quadrotor pointing downwards to measure the OF produced by the visual motion of the ground below.

Measuring Optical Flow in Flight

Method

Here, we evaluate the two algorithms running on the M 2 APix sensor under real flight conditions, in the presence of vibrations, movement disturbances (due to the MAV's attitude stabilization) and light variations. To this end, the M 2 APix sensor was attached to the bottom of a 350-gram X4-MaG quadrotor [START_REF] Manecy | X4-MaG: A Low-Cost Open-Source Micro-Quadrotor and Its Linux-Based Controller[END_REF], aligned to the direction of frontal movement of the MAV. In this configuration, OF measurements were only influenced by the linear speed vector's components of the quadrotor and its pitch rotation movements. Light intensity was measured with a photodiode having the same spectral sensitivity, also attached to the bottom of the X4-MaG and oriented in the same downward direction as the M 2 APix sensor. To get the OF ground-truth with high accuracy, the X4-MaG flew inside the Mediterranean Flying Arena (http://www.flying-arena.eu) equipped with 17 motion-capture cameras covering a 6 × 8 × 6 m volume. The X4-MaG followed a predefined 3D trajectory using the VICON™ (Oxford, UK) system outputs as feedback signals.

The trajectory pattern retained for this experiment was a half-moon shape, positioned at a constant height of 0.4 m above the textured ground (cf. Figure 7). The X4-MaG was controlled to move along the half-circle trajectory segment at a constant speed of 0.5 m•s -1 . During the straight trajectory segment, velocity first increased progressively from 0 to 1.5 m•s -1 and then decreased again to 0 m•s -1 . Between these two segments, the X4-MaG hovered in place while making a pure yaw rotation. Since the speed vector of the X4-MaG is tangential to its trajectory, the half-moon pattern confronted the M 2 APix sensor with three distinct types of movement: translation, rotation and a combination of both. Lighting was provided by three sources: daylight coming through the top windows, fluorescent tubes attached to the ceiling and infra-red light emitted by the motion capture system. Manipulation of the window blinds not only modulated the global intensity of the entering daylight, but also gave rise to high-range variations in lighting on the ground texture (see Figure 8). The intensity of the light emitted by the fluorescent tubes could be adjusted manually by means of a dimmer. Because it was necessary for the proper functioning of the tracking, the 9 × 10 -5 W• cm -2 irradiance provided by the motion capture system could not be altered.

In-Flight Results

Two flight laps of the pattern were recorded, during which the light was arbitrarily modulated. During the first lap, the blinds were left open, and daylight projected very luminous bands on the textured ground (Figure 8) that the drone encountered between the third and eighth seconds of the experiment (Figure 9a). As a result, irradiance during this first lap varied between 10 -3 W•cm -2 and 10 -1 W•cm -2 . The blinds were then fully closed (at 10 s in Figure 9). Manual modulation of the fluorescent tubes induced irradiance variations ranging from 9×10 -5 W•cm -2 to 2 × 10 -3 W•cm -2 .

In Figure 9b1,c1, white zones correspond to periods of forward movement of the drone, whereas gray zones correspond to hovering periods with 90°yaw rotation, orienting the drone in preparation of the next trajectory segment. Half-circle segments are identified by larger white zones between 0 and 10 s and between 18 and 28 s, and straight segments are identified by smaller white zones between 13 and 16 s and between 31 and 34 s. Thus, one large white zone and one small white zone represent one lap. For a trajectory height of 0.4 m and a speed range between 0.5 and 1.5 m•s -1 without pitch rotations, the OF theoretically varies from 72 to 215 °/s.

The OF oscillations observed during the half-circle part are mainly due to the drone's attitude stabilization system, with the resulting pitch rotations inducing changes in OF.

Inspection of Figure 9b revealed that the thresholding method provided a virtually uninterrupted stream of OF measurements. Even in the presence of forward drone movement, measurements however showed quite some variability, as captured by the overall (excluding gray areas) standard deviation (σ) of 43 °/s. During forward movement (white zones), the OF-output refresh rate reached an average of 99 Hz for the 10 LMSs. The cross-correlation method (see Figure 9c) occasionally gave rise to measurement gaps, especially during the high acceleration phases (see Figure 9c at 13 s). However, compared to the thresholding method, measurement variability was lower with σ = 16 °/s, and the OF-output refresh rate was much higher, reaching an average of 1195 Hz for the full M 2 APix sensor (10 LMSs). The foregoing results bring out the differential capabilities of the two methods in different situations. Although, overall, the cross-correlation method provides better results in terms of precision, robustness, working range and output refresh rate, its use in MAV applications may be compromised by a crippling CPU load. Indeed, as demonstrated in Section 4 , for a single M 2 APix sensor, the cross-correlation method, operating at F e = 500 Hz, gave rise to an average CPU load of 52.5%, whereas with the thresholding method operating at F e = 1000 Hz, the CPU load was always lower than 3%. We therefore explored OF-output under reduced working frequencies, allowing lessening of the CPU load.

Offline Results at a Low Sampling Rate

The recorded flight data of the online calculations presented in Figure 9 were down-sampled before running the two algorithms again; the results of this offline OF computation are presented in the same format in Figure 10 with F e = 250 Hz and in Figure 11 with F e = 100 Hz.

Comparison of Figures 9 to 11 (see also Table 1) revealed that for the thresholding method, a reduction of the sampling rate resulted in a stronger OF quantization (especially visible with higher OF values), accompanied by a decrease in OF-output refresh rate (from 99 Hz at F e = 1000 Hz to 36 Hz at F e = 100 Hz) and a widening of the error distribution. For the cross-correlation method, a reduction in sampling rate did not affect the error distribution to any noticeable extent, whereas the OF-output refresh rate decreased proportionally with the sampling rate (from 1195 Hz at F e = 500 Hz to 264 Hz at F e = 100 Hz).

Discussion and Conclusions

Experiments were performed to test two time-of-travel algorithms with an auto-adaptive bio-inspired silicon retina, called M 2 APix, intended for MAV applications. The M 2 APix sensor is a prototype sensor [START_REF] Mafrica | A bio-inspired analog silicon retina with Michaelis-Menten auto-adaptive pixels sensitive to small and large changes in light[END_REF] composed of only 12 auto-adaptive pixels forming 10 LMSs. The auto-adaptive pixels allowed us to evaluate each time-of-travel algorithm largely independently of lighting conditions. In the field of robotics, working in a seven-decade range of irradiance is a considerable advantage because robots may then be led to work in environments subject to strong lighting variations, such as obstacle forests, urban canyons or inside buildings.

In the experiments underlying Figure 5, we controlled the lighting in a five-decade range of irradiance. The lower limit was reached at around 10 -6 W•cm -2 irradiance, when the M 2 APix sensor was no longer able to deliver visual signals.

As can be seen in Figure 6, our results revealed that both time-of-travel algorithms are accurate at all OF speeds tested. However, the precision of the cross-correlation method (C V < 0.1 in Figure 6b) is considerably better than the thresholding method (C V > 0.1 in Figure 6a) over a wide OF range. At the lowest OF (25 • /s), both methods strongly lack precision, probably due to the auto-adaption: working at a 0.5-s time constant, the auto-adaptive process may affect OF measurement at low speed when a time delay ∆t approaches 0.5 s. It is in fact difficult to measure both slow and fast OFs with a constant inter-pixel angle ∆ϕ; with ∆ϕ = 1.5 • , the thresholding allowed measuring an OF range of [1.5 • /s; 25 • /s] in real flight conditions [START_REF] Sabiron | Low-speed optic-flow sensor onboard an unmanned helicopter flying outside over fields[END_REF]. Modulating ∆ϕ as a function of the OF required appears to be the only way to measure the OF in the three-decade range (e.g., [1 • /s; 1000 • /s]) with a time-of-travel algorithm.

Experiments with a moving pattern and in real flight showed that even under varying lighting conditions, both time-of-travel algorithms are quite accurate with respect to ground-truth OF. The cross-correlation method however appears to be more precise with a higher output refresh rate than the thresholding method. For MAV applications with fully-embedded computational resources, the cross-correlation method requires 18-times more computational resources than the thresholding method. The cross-correlation was therefore tested offline at a low sampling rate (250 Hz in Figure 10b1, 100 Hz in Figure 11b1) and found to work well despite some brief measurement gaps when OF accelerated beyond 260 • /s 2 (e.g., Figure 10b1). The thresholding method, on the other hand, only worked adequately at a high sampling rate (1 kHz in Figure 9b1); results were poor at a low sampling rate (250 Hz in Figure 10a1 or 100 Hz in Figure 11a1). For future work, it appears worthwhile to explore combinations of both methods by merging their OF measurements (e.g., a cross-correlation method running at 250 Hz with a thresholding method working in parallel and running at 1 kHz), as this would allow profiting from the benefits of each method without too strongly loading the CPU.

Results presented in Table 1 will allow future research to choose between the time-of-travel algorithms for measuring the OF as a function of embedded computational resources. If an algorithm using little computational resources is required (e.g., for computing dozens or hundreds of LMSs, while using one and the same embedded target), the thresholding method may be preferred. By adding a median filter eliminating matching errors [START_REF] Roubieu | Two-Directional 1-g Visual Motion Sensor Inspired by the Fly's Eye[END_REF], the thresholding method moreover becomes more precise. Its precision can also be improved when running at 200 Hz ≤ F e < 1 kHz using a linear interpolation applied to the photoreceptor signals [START_REF] Expert | Interpolation based "time of travel" scheme in a Visual Motion Sensor using a small 2D retina[END_REF]. While the cross-correlation method is inherently both accurate and precise and can provide a high output refresh rate (that is, about proportional to the sampling rate F e ), it requires much more computational resources (also approximately proportional to the sampling rate F e ), which will limit the number of LMSs that can be embedded into the same target. The CPU load of the cross-correlation method may be reduced by using on-board digital signal processing (DSP).

Finally, from a robotics point of view, our 350-gram X4-MaG quadrotor was equipped with an IMU working at 250 Hz to stabilize attitude and with a downward OF sensor, composed of 10 LMSs running at a refresh output frequency close to the IMU sampling rate. A high refresh rate for both inertial and visual sensors will be a perquisite in the near future to endow MAVs with strong reactive abilities in unpredictable environments. Moreover, this study has demonstrated that a time-of-travel algorithm coupled with a M 2 APix sensor can accurately measure the OF in a wide range [25°/s; 1000°/s], allowing MAVs to fly close to surrounding obstacles at high speed in any lighting conditions.

Figure 1 .

 1 Figure 1. (a) The M 2 APix sensor chip with its wire bounding; (b) the silicon retina composed of 12 Michaelis-Menten pixels and 12 additional Delbrück pixels; (c) view of three neighboring pixels.From[START_REF] Mafrica | A bio-inspired analog silicon retina with Michaelis-Menten auto-adaptive pixels sensitive to small and large changes in light[END_REF].

10 Figure 2 .

 102 Figure 2. The 12 Michaelis-Menten pixels (Figure 1b) and the 10 corresponding local motion sensors (LMSs). LMSs are unidirectional and can only measure the OF magnitudes OF i i ∈ {1, ..., 10} in 1D.The optics geometry of the first two-pixel LMS (composed of the red and blue pixels) is described in Figure3a.

Figure 3 .

 3 Figure 3. (a) Optics geometry of a LMS. A defocused lens confers an inter-pixel angle ∆ϕ between the optical axes of two adjacent pixels Ph1 (blue) and Ph2 (red) and an acceptance angle ∆ρ given by the width of the Gaussian angular sensitivity at half height. The output signals V ph1 and V ph2 coming from pixels Ph1 and Ph2, respectively, are naturally delayed by a "time-of-travel" ∆t when a contrast is moving in front of them (Figure 3b,c). The OF can therefore be computed with Equation (1). (b) Principle of the time-of-travel algorithm based on signal thresholding: the time delay ∆t is computed when both signals reach a given hysteresis threshold [4,7,8,29]. (c) Principle of the time-of-travel algorithm based on signals' cross-correlation: the time delay ∆t is computed as the time delay giving the maximum cross-correlation coefficient between the delayed and non-delayed signal (see Figure2in[START_REF] Mafrica | Minimalistic optic flow sensors applied to indoor and outdoor visual guidance and odometry on a car-like robot[END_REF] for details). Note that for simplicity, the time-of-travel computation is shown here using the raw signals V ph1 and V ph2 , but it is actually computed after band-pass filtering these signals (see[START_REF] Franceschini | Optic Flow Based Visual Guidance: From Flying Insects to Miniature Aerial Vehicles[END_REF][START_REF] Expert | A mouse sensor and a 2-pixel motion sensor exposed to continuous illuminance changes[END_REF][START_REF] Roubieu | Two-Directional 1-g Visual Motion Sensor Inspired by the Fly's Eye[END_REF][START_REF] Mafrica | Minimalistic optic flow sensors applied to indoor and outdoor visual guidance and odometry on a car-like robot[END_REF][START_REF] Sabiron | Low-speed optic-flow sensor onboard an unmanned helicopter flying outside over fields[END_REF] for details).

Figure 4 .

 4 Figure 4. Experimental setup using a moving pattern and varying the artificial lighting to test algorithms with the M 2 APix sensor both in various OF ranges and in various lighting conditions. The pattern speed V p varied from 0 to 1.4 m•s -1 ; the sensor was placed perpendicularly in front of the moving pattern at a distance d from 8 cm to 40 cm.

Figure 5 .

 5 Figure 5. (a)The five levels of irradiance from 6 × 10 -7 W•cm -2 to 1.6 × 10 -2 W•cm -2 (i.e., from 0.2 lux to 12,000 lux for human vision) used during the first part of the experiments (see the experimental setup shown in Figure4). (b,c) OF measurements obtained with the M 2 APix sensor using (b) the thresholding method (F e = 1 kHz) and (c) the cross-correlation method (F e = 500 Hz) for the 10 LMSs during 10 laps of the moving pattern. Tests were run separately for each of the five lighting conditions. Each blue dot represents one OF measurement of one LMS; however, the mean output refresh rate is the average of the number of OF measurement from the 10 LMSs in 1 s. The green line represents the OF ground-truth computed using Equation (2). (d) Error distribution of each time-of-travel method: the thresholding method is in red and the cross-correlation method in blue.

Figure 7 .

 7 Figure 7. MAV trajectory: (a) top view; (b) perspective view. The two laps were performed at a flight height of 0.4 m.

Figure 8 .

 8 Figure 8. X4-MaG drone hovering at 1.2 m above the natural texture in the flight arena. The half-moon shaped trajectory is represented in green. The M 2 APix sensor is embedded under the quadrotor looking downwards and measuring the OF coming from the ground below. See Video S1 to watch the experiment.

Figure 9 .

 9 Figure 9. (a) Light intensity in W•cm-2 measured with a photodiode oriented in the same direction as the M 2 APix sensor embedded under the X4-MaG; (b,c) In the left panels (b1,c1), blue points represent single-LMS OF measurements; the green line is the ground-truth provided by high-precision Vicon tracking. Right panels (b2,c2) present the OF error distribution, with 1°/s beams in the [300°/s; 300°/s] range, with respect to the ground-truth during forward drone movement (white zones in left panels). The excluded (gray) zones correspond to periods of pure yaw rotation during the MAV trajectory. See Supplemental Movie 1 to watch the experiment.

Figure 10 .Figure 11 .

 1011 Figure10. OF computed with (a) the thresholding method and (b) the cross-correlation method after down-sampling the data underlying Figure9to F e = 250 Hz.

Table 1 .

 1 Comparison between the results obtained with the two time-of-travel algorithms.

			Thresholding Method		Cross-Correlation Method
	Sampling Rate (F e )	1 kHz 500 Hz 250 Hz 100 Hz	1 kHz	500 Hz 250 Hz 100 Hz
	CPU Load (%)	2.2	1.1 *	< 1 *	< 1 *	overload	52.5	26.3 *	10.5 *
	Precision σ (°/s)	43	44	49	47	-	16	17	20
	Refresh rate (Hz/10 LMSs)	99	51	48	36	-	1195	701	264

* Theoretical value because corresponding tests were made offline.

Acknowledgments:

We thank Julien Diperi for the mechanical design of the flying robot and the test bench and Marc Boyron for his involvement in the overall electronic design of the flying robot. This work was supported by the French Direction Générale de l'Armement (DGA), CNRS, Aix-Marseille University, the Provence-Alpes-Côte d'Azur region and the French National Research Agency for Research (ANR) with the Equipex/Robotex project.

), and the OF produced by the moving pattern was made to vary from 25 °/s to 1000 °/s with steps of 50 °/s (except for the first step). Each graph shows the average µ (blue circles), the standard deviation σ (blue vertical bars) and the coefficient of variation σ µ (red lines) of the OF measurements with respect to the angular speed of the moving pattern (see Figure 4). The ground-truth values of the OF are given by the blue lines.

Supplementary Materials:

The following are available online at http://www.mdpi.com/1424-8220/17/3/571/s1, Video S1: a video of the experiment with the micro flying robot.

Author Contributions: E.V., J.S. and F.R. conceived of the experiments. E.V. and S.M. implemented the algorithms. E.V. conducted the experiments. All authors actively participated in writing the manuscript.

Conflicts of Interest:

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:

OF Optical flow LMS

Local motion sensor M 2 APix Michaelis-Menten auto-adaptive pixel


David Navarro-Alarcon

Omar Zahra

A Depth-Based Algorithm for Manipulating Deformable Objects Using Smooth Parametric Surfaces and Energy Minimisation

In this brief work, we present a new method for controlling deformations of soft objects by using parametric surfaces as a new type of deformation feedback features. This new approach allows us to actively deform objects into complex 3D shapes. A kinematic-based motion controller is derived using an energy minimisation strategy.

I. THE METHOD Consider a set of R robotic manipulators (or fingers) that rigidly grasp a deformable object O. We denote the 3D position of the ith robot (measured at the time instant t) with the vector r i t . The total position vector of the R robots is denoted by

r t = [r 1 t , . . . , r R t ] ⊺ ∈ R 3R .
To measure the object's deformations, we use a static RGB-D sensor (e.g. Kinect) that continuously observes the manipulated object. Let us denote by p a 3D point located over the portion of object's surface that can be observed by the image sensor. It is assumed that this surface is smooth and twice differentiable, therefore, its 3D coordinates can be parameterised as p = [x, y, f (x, y)]. However, since the object's deformation model is unknown, the analytical structure of z = f (x, y) is also not known.

The object's feedback surface can be approximated with the following parametric function:

f (x, y) = o + M -1 j=0 N -1 i=0 a i,j b i,j c i,j d i,j • cos(ixπ/l x ) sin(ixπ/l x ) cos(jyπ/l y ) cos(jyπ/l y ) = g(x, y) • s (1)
where N, M > 0 is the number of harmonics taken into account, s = [a i,j , b i,j , c i,j , d i,j , o]

⊺ represent the vector of constant surface parameters, and l x , l y > 0 are scaling factors. Note that the vectorial function g(x, y) is known and can be constructed with feedback data.

To approximate the unknown parameters s, we first collect L data points of p k = [x k , y k , z k ] at various locations over the surface, and then construct the long vector and matrix

z = z 1 • • • z L ⊺ (2) G = g(x 1 , y 1 ) ⊺ • • • g(x L , y L ) ⊺ ⊺ (3)
With this data, the vector s t of surface parameters can be computed at every time instant t as follows:

s t = G ⊺ G -1 G ⊺ z ( 4
)
This work is supported in part by the HKPolyU under grant 4-ZZHJ, and in part by the HK RGC under GRF grant 14203917.

The authors are with the Department of Mechanical Engineering, Hong Kong Polytechnic University, HKSAR. e-mail: dna@ieee.org The vector s t can be used to characterise the object's feedback surface. To manipulate/deform the object into a target configuration, we must find the relation between the (input) motion of the robots and the (output) surface parameters [START_REF] Navarro-Alarcon | Uncalibrated vision-based deformation control of compliant objects with online estimation of the Jacobian matrix[END_REF], [START_REF] Navarro-Alarcon | Automatic 3D manipulation of soft objects by robotic arms with adaptive deformation model[END_REF]. We approximate this unknown model with the following linearly parametrisable expression:

s t+1 = s t + A t u t (5)
where the mapping matrix A t represents the deformation properties of the object and the vector u t = r tr t-1 stands for the controllable motions of the robots. This unknown deformation matrix A t can be computed online by using the method described in [START_REF] Navarro-Alarcon | Fourier-based shape servoing: A new feedback method to actively deform soft objects into desired 2D image shapes[END_REF]. Given a desired surface f d (x, y) described by the parameters s * , a controller is derived by minimising the following quadratic cost function:

P = s t+1 -s * , s t+1 -s * H + u t , u t Q (6) for H = H ⊺ > 0 and Q = Q ⊺ > 0 as weight matrices.
Computing the extremum ∇P(u t ) = 0 and solving for u t yields the motion control input:

u t = -(A ⊺ t HA t + Q) -1 A ⊺ t H (s t -s * ) (7)
that iteratively minimises P. Note, however, that this method only guarantees local minimisation of the cost function.

Fig. 1 .

 1 Fig. 1. The robotic system manipulating a deformable object.


Sophie Mercier

Michel Roussignol

A block replacement policy for a bivariate wear subordinator

Classical models for wear accumulative indicators are univariate Gamma and compound Poisson processes, which both are univariate subordinators (increasing Lévy processes). Bivariate subordinators are here proposed to model correlated couples of univariate wear indicators, namely bivariate wear indicators. A few properties of these bivariate subordinators are pointed out, and a special simple construction is provided, for application purpose. The use of these bivariate wear subordinors is illustrated through the study of a classical block replacement policy. The influence of the dependence between the marginal indicators on the optimal block replacement policy is pointed out.

INTRODUCTION

In case of a system submitted to an accumulative random damage, classical stochastic models are compound Poisson processes and Gamma processes, which both are increasing Lévy processes, see [START_REF] Van Noortwijk | A survey of the application of Gamma processes in maintenance[END_REF], [START_REF] Nakagawa | Shock and damage models in reliability theory[END_REF] or [START_REF] Abdel-Hameed | A gamma wear process[END_REF] e.g.. Such classical wear models typically are univariate. However, the deterioration level of a system cannot always be synthetized into one single indicator and several indicators may be necessary, see [START_REF] Mercier | Bivariate Gamma wear processes for track geometry modelling, with application to intervention scheduling[END_REF]) for an industrial example. In that case, a multivariate wear model must be used to account for the dependence between the different univariate indicators of the system. Another context where multivariate wear models are required is the case of different systems submitted to common stresses, which make their wear indicators dependent. Multivariate increasing stochastic models hence are of interest in different contexts.

We here propose to use multivariate increasing Lévy processes (or multivariate subordinators) as wear processes. Two main arguments are the following: first, such processes present some positive dependence property (consequence of [START_REF] Bäuerle | Dependence properties and comparison results for lévy processes[END_REF])); secondly, thresholds reaching times may be proved to present aging properties. Such properties are typical of what one may expect from a multivariate wear model: indeed, correlation between wear indicators reveals some common stresses or so, and high values for one of these indicators should entail higher values for the other ones. Roughly speaking, this is the positive dependence property. Also, in a reliability context, thresholds reaching times for increasing wear processes typically represent times to failure and are consequently likely to present aging properties. Multivariate subordinators hence are good candidate to model multivariate deterioration. All the paper is written in the bivariate setting but might easily be translated into the multivariate one.

Considering bivariate subordinators as bivariate wear processes, the objective of this paper is to illustrate their ability to optimize maintenance strategies. With such an aim, we revisit a classical block replacement policy within this new context, with the optimization of a cost function as an objective. In view of the numerical assesment of the cost function, and also of an effective use for the applications, a special construction of a bivariate subordinator is proposed, by so-called trivariate reduction. Though valid in the most general setting, this construction is here proposed in the case of Gamma marginal processes, to make it clearer.

The paper is organized as follows: in Section 2, known and easy facts are summed up for univariate and bivariate subordinators. The construction of a bivariate subordinator by trivariate reduction is given in Section 3, in the case of Gamma marginal processes. The block replacement policy is presented in Section 4 and the results are numerically illustrated in Section 5. Conclusive and prospective remarks end this paper in Section 6.

A univariate subordinator is a univariate nondecreasing Lévy process, namely a process (Y t ) t 0 with range R + , starting from Y 0 = 0, which is non decreasing, right continuous with left-side limits and with stationary independent increments. Only subordinators with null drifts are here envisionned, both in the univariate and bivariate cases. The Laplace transform of a univariate subordinator (Y t ) is given by:

L Yt (x) = E e xYt = exp 0 @ t Z R + e xy 1 Y (dy) 1 A (1)
for all x; t 0, where Y (dy) is a non-negative measure on R + such that

Z R + (jyj ^1) Y (dy) < 1;
where symbol ^stands for minimum. This measure Y (dy) is the so-called Lévy measure of the process (Y t ) t 0 . The Laplace transform and the distribution of the process are characterized by its Lévy measure Y (dy). An alternate characterization is provided by the tail integral function U Y , with

U Y (y) = Z [y;+1[ Y (dx) for y 2 R + , U Y (0) = +1 and U Y (+1) = 0.
A well-known example of a univariate subordinator is the Gamma process. Its Lévy measure is equal to Y (dy) = y exp( y) dy, where ; > 0. A great interest of the Gamma process is the explicit knowledge of the probability distribution of Y t (Gamma distributions). For other subordinators, a close form is generally not available. Nevertheless, recent researches in statistics provide tools for estimating the Lévy measure of general (univariate) subordinators, see [START_REF] Comte | Nonparametric adaptive estimation for pure jump Lévy processes[END_REF] e.g.. Just as in the univariate case, a bivariate subordinator is a bivariate non-decreasing Lévy process (here with null drift), namely a process

X t = X (1) t ; X (2) t t 0
with range R 2 + , starting from X 0 = (0; 0), which is component-wise nondecreasing, right continuous with left-side limits and with stationary independent increments. To avoid trivialities, P X

(1) t > 0; X (2) t > 0 is assumed to be non zero.

The Laplace transform of X t is

L Xt (x 1 ; x 2 ) = E e x 1 X (1) t +x 2 X (2) t = exp 0 B @ t ZZ R 2 + e (x 1 y 1 +x 2 y 2 ) 1 X (dy 1 ; dy 2 ) 1 C A (2) for all x = (x 1 ; x 2 ) 2 R 2
+ and all t 0. Here again the Laplace transform and the distribution of the process are characterized by the Lévy measure X which fulfills: ZZ

R 2 + (kyk ^1) X (dy) < 1;
where k:::k stands for any standard norm on R 2 + . The Lévy measure can be characterized by the bivariate tail integral function U X : [0; +1] 2 ! [0; +1] provided by:

U X (x 1 ; x 2 ) = Z [x 1 ;+1[ Z [x 2 ;+1[ X (dy) (3)
for all (x 1 ; x 2 ) 2 R 2 + n f(0; 0)g, U X (x 1 ; 1) = U X (1; x 2 ) = 0 and U X (0; 0) = +1.

The bivariate tail integral function measure the dependence between the two marginal processes in the sense that if X and X are bivariate subordinators with the same marginal processes and respective tail integral functions U X (x) and U X (x) such that U X (x) U X (x) for all x in R 2 + , then the components X

(1) t and X

(2) t are "more dependent" than the components X(1) t and X(2) t (with clear notations). As an illustration, it is easy to check that

U X (x 1 ; x 2 ) U X k (x 1 ; x 2 ) (4) for all (x 1 ; x 2 ) 2 [0; +1] 2 , where U X k (x 1 ; x 2 ) = U 1 (x 1 ) ^U2 (x 2 )
stands for the tail integral function in the case where the components are completely dependent (case of maximal dependence). The maximal dependence is hence reached for the maximal tail integral function.

Given a bivariate subordinator X

(1) t ; X

(2) t t 0

, for each i = 1; 2, the marginal process

X (i) t t 0
is a univariate subordinator with associated Lévy measure i where 1 (dx 1 ) = X (dx 1 R + ) (the same for 2 ) and (univariate) tail integral function denoted by U i .

The linear correlation coefficient Xt of X t is a measure of the dependence between X

(1) t and X

(2) t . Under technical assumptions ensuring its existence, a direct calculation of covariance and standard deviations with the help of Laplace transforms (1) and (2) gives:

Xt = RR R 2 + U X (u 1 ; u 2 ) du 1 du 2 p v 1 v 2 (5)
with

v i = var X (i) 1 = 2 Z R + u U i (u) du for i = 1; 2. Expression (5)
shows that Xt is independent of t (just as for any Lévy process). It is denoted by X in the following. Using 0 U X (x 1 ; x 2 ) and (4), we derive the following bounds for X :

0 X max with max := RR R 2 + (U 1 (u 1 ) ^U2 (u 2 )) du 1 du 2 p v 1 v 2 (6)
The lower bound X = 0 corresponds to the case of independent indicators. The upper bound max is reachable and corresponds to the case of complete (maximal) dependence. We have max = 1 in case U 1 = U 2 , but in general max < 1, see Figure 1 later on. This result shows that, in case of non equal marginal subordinators, high dependence between the marginal processes cannot always be modeled by a bivariate subordinator.

We finally recall a comparison result from [START_REF] Bäuerle | Dependence properties and comparison results for lévy processes[END_REF], to be used later on. This result roughly says that a bivariate subordinator is "increasing" with dependence. More specifically, if X and X are bivariate subordinators with the same marginal processes and respective tail integral functions U X (x) and U X (x) such that U X (x) U X (x) for all x in R 2 + , then X is smaller than X in the sense of the concordance order. This means that:

P X (1) t > s 1 ; X (2) t > s 2 P X(1) t > s 1 ; X(2) t > s 2 ; P X (1) t s 1 ; X (2) t s 2 P X(1) t s 1 ; X(2) t s 2 ;
(7) for all t; s 1 ; s 2 0. As already told, there is a special interest for bivariate subordinators with Gamma marginal processes, because the marginal probability distributions are explicitly known (Gamma distributions). However, no closed form is generally available for the process bivariate distribution. Nevertheless the construction given in the next section provides a special case where such a distribution can be obtained.

CONSTRUCTION BY TRIVARIATE REDUCTION

We here propose a model of bivariate subordinators constructed by trivariate reduction (TR subordinator), in the special case of Gamma marginal processes. This method has already been introduced in [START_REF] Buijs | Time dependent reliability analysis of flood defences using gamma processes[END_REF] or Ebrahimi, and used in [START_REF] Mercier | Modelling track geometry by a bivariate gamma wear process, with application to maintenance[END_REF]. Results on bivariate Gamma distributions constructed by trivariate reduction may also be found in Devroye [START_REF] Devroye | Non-uniform random variate generation[END_REF], section XI.3).

Let us first recall that the probability distribution function (p.d.f.) of a Gamma distribution (a; b) is

f a;b (u) = b a (a) u a 1 e bu 1 R + (u)
and that the distribution at time t of a univariate Gamma process with parameters (a; b) is (at; b).

The tail integral function of a Gamma process with parameters (a; b) is:

U (x) = Z +1 x a e bu u du = aE 1 (bx)
where

E 1 (x) = R +1 x e y
y dy is the exponential integral function.

For a general subordinator with univariate Gamma processes as margins and marginal parameters (a 1 ; b 1 ) and (a 2 ; b 2 ), the upper bound (6) is then provided by:

max ( ) = ZZ R 2 + min p E 1 (u 1 ) ; 1 p E 1 (u 2 ) du 1 du 2 (8)
where = a 1 a 2 . Let us now describe the trivariate construction: starting from three independent univariate Gamma processes Y (i) t t 0 with parameters ( i ; 1) for i 2 f1; 2; 3g and from b 1 ; b 2 > 0, we set: 8 < :

X (1) t = Y (1) t + Y (3) t =b 1 X (2) t = Y (2) t + Y (3) t =b 2 (9)
The interest of such a construction is clear for the applications to reliability, where the term Y

(3) t

will correspond to some common damage accumulation coming from some common wear source, whereas Y

(1) t and Y

(2) t will correspond to some damage accumulations coming from independent wear sources.

The process (X t ) t 0 = X (1) t ; X

(2) t t 0 constructed in (9) is a bivariate subordinator with univariate Gamma processes as margins. Its marginal parameters (a i ; b i ) are provided by: A direct computation provides:

a 1 = 1 + 3 a 2 = 2 + 3
X = 3 p a 1 a 2
and setting = X :

( 1 = a 1 p a 1 a 2 2 = a 2 p a 1 a 2 3 = p a 1 a 2 (10)
Two equivalent alternate parametrizations hence are available for TR subordinators with Gamma marginal processes: either ( 1 ; 2 ; 3 ; b 1 ; b 2 ) or (a 1 ; a 2 ; ; b 1 ; b 2 ). This entails that all the dependence of a TR subordinator with Gamma marginal processes is contained in the linear correlation coefficient . Equations (10) also entail that:

min (a 1 ; a 2 ) p a 1 a 2 = min p ; 1 p := max;T R ( ) (11)
where := a 1 a 2 . In case a 1 = a 2 , we get max; T R = max = 1. In the general case, we however get max; T R < max , as shows Figure 1 where the functions 7 ! max ( ) and 7 ! max; T R ( ) are plotted. This proves that a trivariate construction cannot model as high a dependence as general subordinators, in case of different marginal processes.

For the practical use of TR subordinators in Section 5, we now provide the joint p.d.f. f Xt (u 1 ; u 2 ) of X t : when 3 6 = 0, it is equal to:

Z min(u 1 ;u 2 ) 0 f 1 t;1 (u 1 u 3 )f 2 t;1 (u 2 u 3 )f 3 t;1 (u 3 )du 3 (12) When 3 = 0 (or = 0), both X (1) t and X (2) t are inde- pendent and f Xt (u 1 ; u 2 ) = f a 1 t (u 1 )f a 2 t (u 2 ).

BLOCK REPLACEMENT POLICY

In this section we illustrate the ability of using bivariate subordinator to optimize maintenance strategy with a block replacement policy: a system is considered, with degradation measured by a bivariate subordinator (X t ) t 0 (with null drift). For each univariate subordinator

X (i) t t 0
, there is a threshold s i above which the deterioration level is considered as too high and we set:

(i) s i = inf t 0 such that X (i) t s i
to be the corresponding univariate reaching time (i = 1; 2). The random variable (i) hence stands for a univariate time to failure, only taking indicator X (i) t into consideration. As for the bivariate case, we envision two possibilities for the system failure: either the system is failed as soon as one (univariate) indicator is too high or when both indicators are. If we interpret X

(1) t and X

(2) t as the wear indicators of two dependent components, the first case corresponds to a seriessystem and the second one to a parallel system. The respective times to failure are:

= inf t 0jX (1) t s 1 or X (2) t s 2 = min (1) s 1 ; (2)
s 2 and = inf t 0jX

(1) t s 1 and X

(2) t

s 2 = max (1) s 1 ; (2) s 2 :
The block replacement policy is then applied: the system is not continuously observed. It is perfectly and instantaneously repaired each T time unit (periodically) with a cost c r . If the system is down at repair time, a unitary cost c u is induced per unit time for the down period just before the repair. After a repair, the future evolution of the system is similar as from the beginning and is independent of its past. We are interested in the asymptotic unitary cost (per unit time) denoted by C (T ). We first recall classical results for sake of completeness, see [START_REF] Nakagawa | Block replacement[END_REF]) e.g.. Using standard renewal theory, the asymptotic unitary cost C (T ) is provided by:

C (T ) = E (C ([0; T ])) T
where E (C ([0; T ])) is the mean cost on a generic cycle [0; T ]. The asymptotic unitary cost then is:

C (T ) = c r + c u R T 0 F (t) dt T ( 13
)
where F (t) is the cumulative distribution function (c.d.f.) of the failure time .

A finite optimal T opt exists if and only if E ( ) > cr cu . In that case, this optimal point is characterized by

Z T opt 0 t f (t) dt = c r c u
and the optimal cost is:

C T opt = c u F T opt : (14)
In case, E ( ) cr cu , the cost function C (T ) is non increasing with T and the best is not to ever repair the system (T opt = +1).

We now complete such results in the present context and we first consider the case where = min (1)

s 1 ; (2) s 2 .
In that case, we have:

F (t) = P ( t) = 1 P X (1) t s 1 ; X (2) t s 2
for all t 0 and

E ( ) = Z +1 0 P ( > t) dt = Z +1 0 P X (1) t s 1 ; X (2) t s 2 dt
As a consequence, it can be seen from (7) that F (t), C (T ) and C (T opt ) are non increasing with dependence (see ( 13) and ( 14)), whereas E ( ) is non decreasing with dependence. This entails that

E ? E ( ) E k
where E ? and E k stand for the mean times to failure in case of independent and completely dependent wear indicators, respectively (with fixed margins).

In case E ? > cr cu , we hence have E ( ) > cr cu and a finite T opt exists for any kind of dependence.

If E k cr cu , a finite T opt never exists. If E ? cr cu < E k , a finite T opt exists for high dependence but does not exist for low dependence.

This shows the clear influence of the dependence on the optimal maintenance policy.

Similar results are valid in case = max (1)

s 1 ;
(2) s 2 . In that case, we get:

F (t) = P X (1) t > s 1 ; X (2) t > s 2 ; E ( ) = Z +1 0 1 P X (1) t > s 1 ; X (2) t > s 2 dt
and all conclusions are reversed. In this way, F (t), C (T ) and C (T opt ) are non decreasing with dependence, whereas E ( ) is non increasing with dependence. In case E k cr cu < E ? , a finite T opt exists for low dependence but does not exist for high dependence.

NUMERICAL EXPERIMENTS

This section illustrates the previous results for a bivariate subordinator with Gamma marginal processes, constructed by trivariate reduction (see Section 3).

Numerical computations

Based on the explicit form (12) of the p.d.f. of

X (1) t ; X (2) t
, it is possible to compute all the quantities analytically provided in Section 4.

As an example, in case = min (1)

s 1 ;
(2) s 2 and > 0, we have:

F (t) = 1 Z s 1 0 Z s 2 0 f Xt (u 1 ; u 2 ) du 1 du 2 = 1 Z min(s 1 ;s 2 ) 0 F 1 t;1 (s 1 u 3 ) F 2 t;1 (s 2 u 3 ) f 3 t;1 (u 3 ) du 3
where F a;b is the cumulative distributive function of the Gamma distribution with parameters a and b.

This provides:

Z T 0 F (t) dt = T Z T 0 Z min(s 1 ;s 2 ) 0 F 1 t;1 (s 1 u 3 ) F 2 t;1 (s 2 u 3 ) f 3 t;1 (u 3 ) du 3 dt
Due to numerical difficulties, we actually prefer write:

Z T 0 F (t) dt = T Z T 0 Z min(s 1 ;s 2 ) 0 1 F 1 t;1 (s 1 u 3 ) 1 F 2 t;1 (s 2 u 3 ) f 3 t;1 (u 3 ) du 3 dt = Z T 0 F 3 t;1 (min (s 1 ; s 2 )) dt + Z T 0 Z min(s 1 ;s 2 ) 0 F 1 t;1 (s 1 u 3 ) f 3 t;1 (u 3 ) du 3 dt + Z T 0 Z min(s 1 ;s 2 ) 0 F 2 t;1 (s 2 u 3 ) f 3 t;1 (u 3 ) du 3 dt Z T 0 Z min(s 1 ;s 2 ) 0 F 1 t;1 (s 1 u 3 ) F 2 t;1 (s 2 u 3 ) f 3 t;1 (u 3 ) du 3 dt
and compute the four involved quantities, which provides more accurate results. When = 0, the components are independent and, in case = min (1)

s 1 ;
(2) s 2 , we get:

F (t) = 1 F 1 t;1 (s 1 )F 2 t;1 (s 2 ):
Different numerical experiments are next presented in order to illustrate the behaviour of C(T ), T opt and C(T opt ) with respect of T and/or . As for the influence of the other parameters, C (T ) is clearly increasing with respect of c r , c u , a 1 , a 2 , b 1 and b 2 , and decreasing with respect of s 1 and s 2 , so that no further study is required.

Examples 1 to 3 correspond to the case where = min (1)

s 1 ; (2) s 2 .

Example 1

The parameters of the process are a 1 = 50, a 2 = 70, b 1 = b 2 = 1, = 0:7. The failure levels are s 1 = 50 and s 2 = 40. The costs are c r = 1 and c u = 50. Figure 2 gives C (T ) with respect to T and the asymptotic unitary costs C (i) (T ) for each univariate indicator (i = 1; 2). Such unitary costs are computed taking into account one single component X (i) t , which means that the time to failure is considered to be (i) . The induced costs are assumed to be the same as for the bivariate case. Both components are renewed at the same time when the system is renewed and this induces the same renewal cost for the bivariate and univariate cases. In Figure 2, one can observe that C (T ) and C (2) (T ) are nearly superimposed and are clearly higher than C (1) (T ). This comes from the fact that the different parameters are such that E (1) s 1 ' 1:01 and E

(2) s 2 ' 0:58. As a consequence, the second component X

(2) t will fail before the first one X

(1) t with a high probability and in that case, it is useless to use both components X

(1) t and X

(2) t . In case of very different components, the bivariate model consequently does not bring much and one can simply consider the "worse" one.

Example 2

The parameters of the marginal processes are a 1 = 50, a 2 = 74, b 1 = b 2 = 1. The failure levels are s 1 = 20 and s 2 = 30. The costs are c r = 10 and c u = 500. In this case a 1 < a 2 , s 1 < s 2 but the mean univariate times to failures are close:

E (1) s 1 ' E (2) s 2
' 0:41. The function C (T ) is computed for = 0, and = T R max . The results are presented in Figure 3, as well as C (i) (T ) for i = 1; 2. In that case, the univariate costs C (i) (T ) are very close one from each other. The bivariate cost C (T ) is observed to be higher. In that case, not taking into account both components may lead to an under-estimated cost fonction.

Example 3

The parameters of the marginal processes are a 1 = a 2 = 50, b 1 = b 2 = 1. The failure levels are s 1 = s 2 = 30. The costs are c r = 10 and c u = 500. The two components have the same distribution and the same failure level, so that the marginal times to failure are identical. The cost function C (T ) is computed for different values of and the results are presented in Figure 4. As expected, C (T ) is observed to be decreasing with dependence. In case of an unknown correlation, the most conservative attitude hence is to take into account both components and to consider them as independent. Also, the univariate costs C (i) (T ) for i = 1; 2 coïncides with the bivariate cost function C (T ) in case = 1. Taking into consideration one single indicator hence lead to an under-estimated cost fonction. The optimal cost C(T opt ) and period T opt are next plotted against the linear correlation coefficient in Figures 5 and6. As expected, C(T opt ) is nonincreasing with dependence. The optimal period T opt is observed to be non-decreasing with dependence. That point has not been proved from a theoretical point of view.

Example 4

The numerical parameters are the same as for example 3, but the time to failure now is = max (1)

s 1 ; (2) s 2 .
The results are provided in Figures 7 and8. As expected, the monotony of C(T opt ) and T opt with dependence is reversed. The optimal cost is here a little lower than in case = min (1) s 1 ;

(2) s 2 , and the optimal period a little higher. This is conform to what is expected, because a parallel-system is more reliable than a series-system.

Conclusion, Prospective remarks

A block replacement policy has been studied in the new context of a bivariate deterioration modelled by a bivariate subordinator. The influence of the dependence between the marginal wear indicators on the optimal policy has been pointed out. In the case where the system time to failure is the minimum or the maximum of the two univariate times to failure, the optimal cost function has been proved to be monotone with respect to the dependence, with opposite monotonicities for the two cases. When the correlation is unknown, the safest attitude hence is to consider the correlation which leads to the highest cost, namely consider the components as independent in the first case (minimum) and completely dependent in the second case (maximum).

As for the optimal period, it has been numerically observed to be monotone with dependence for all the examples provided here, as well as for all the other examples we may have looked at. It may then be conjectured that this property is always true in the studied context. However, we have not been able to prove it theoretically. This consequently remains on open question.

In the present paper, the choice has been made to consider a simplified model where the system time to failure is the minimum or the maximum of the two univariate times to failure. Thinking of the two components X

(1) t and X

(2) t as the width and heigth of a girder for instance, one may think that a more adequate criterion for high deterioration might be linked to the grid section X

(1) t :X

(2) t e.g.. The time to failure might then has the following shape:

= inf t 0 : X (1) t :X (2) t s = inf t 0 : X (1) t ; X (2) t 2 S
where s is some fixed threshold and S = (x 1 ; x 2 ) 2 R 2 + : x 1 x 2 s . More generally, failure times may correspond to reaching times of failure zones, which may be quite general (S should be an upper set namely such that if (x 1 ; x 2 ) 2 S and (x 1 ; x 2 ) (y 1 ; y 2 ), then (y 1 ; y 2 ) 2 S). Further studies should then be performed for such general failure zones.

Also, the choice has here been made to renew both indicators at each maintenance time. In case of indi-cators linked to two different components, a wiser attitude might be to renew a component only if it is too deteriorated (X (i) t M i , where M i is some maintenance threshold). Such a maintenance policy also requires further investigation.

Figure 1 :

 1 Figure 1: max ( ) and max; T R ( ) with respect of .

Figure 2 :

 2 Figure 2: C (T ) with respect of T .

Figure 3 :

 3 Figure 3: C (T ) with respect of T for = T R max , = 0 and C (i) (T ) for i = 1; 2.

Figure 4 :

 4 Figure 4: C (T ) with resp. of T for = 0, 0:25, 0:5, 0:75 and 1.

Figure 5 :Figure 6 :

 56 Figure 5: C opt ( ) with respect of .

Figure 7 :

 7 Figure 7: C opt ( ) with respect of .

Figure 8 :

 8 Figure 8: T opt ( ) with respect of .