J Hurley
email: jane.hurley@stfc.ac.uk

• N Murdoch

N A Teanby

• N Bowles

T Warren

S B Calcutt

• D Mimoun

• W T Pike

Isolation of Seismic Signal from InSight/SEIS-SP Microseismometer Measurements

Keywords: Mars, InSight, Seismometer, Decorrelation, Noise

The InSight mission is due to launch in May 2018, carrying a payload of novel instruments designed and tested to probe the interior of Mars whilst deployed directly on the Martian regolith and partially isolated from the Martian environment by the Wind and Thermal Shield. Central to this payload is the seismometry package SEIS consisting of two seismometers, which is supported by a suite of environmental/meteorological sensors (Temperature and Wind Sensor for InSight TWINS; and Auxiliary Payload Sensor Suite APSS). In this work, an optimal estimations inversion scheme which aims to decorrelate the shortperiod seismometer (SEIS-SP) signal due to seismic activity alone from the environmental signal and random noise is detailed, and tested on both simulated and Viking data. This scheme also applies a module to identify measurements contaminated by Single Event Phenomena (SEP). This scheme will be deployed as the pre-processing pipeline for all SEIS-SP data prior to release to the scientific community for analysis.

Introduction 1.InSight

The InSight (Interior Exploration using Seismic Investigations, Geodesy and Heat Transport) mission to Mars was selected by NASA in August 2012, under the framework of the DISCOVERY program. With a launch in May 2018, it will deploy the first geophysical observatory on Mars, providing scientific knowledge essential to understand the fundamental processes of telluric planet formation and evolution: an in-situ investigation of the interior of a truly Earth-like planet. To accomplish its objectives, a tightly focussed payload has been assembled: the Seismic Experiment for Interior Structure (SEIS) and the Heat Flow and Physical Properties Package (HP3) [START_REF] Banerdt | The Insight Team), InSight: an integrated exploration of the interior of Mars[END_REF][START_REF] Banerdt | InSight: a discovery mission to explore the interior of Mars[END_REF][START_REF] Lognonne | The InSight/SEIS Team), Science goals of SEIS, the InSight seismometer package[END_REF][START_REF] Mimoun | The InSight SEIS experiment[END_REF].

The InSight mission will use geophysics to explore the Martian interior using seismic and thermal measurements and rotational dynamics, providing information about the initial accretion of the planet, the formation and differentiation of its core and crust, and the subsequent evolution of the interior. Knowledge of Mars' geophysics feeds directly into our understanding of what processes formed the rocky planets; of which the current knowledge is limited to Earth and the moon. Unlike the Moon, Mars is large and complex enough to have undergone most of the processes that affected early Earth but has not undergone extensive plate tectonics or other major reworking that erased evidence of early events.

Mars seismic activity likely stems from faulting (due to stress-release from the crust as the interior cools; [START_REF] Golombek | A prediction of Mars seismicity from surface faulting[END_REF][START_REF] Knapmeyer | Working models for spatial distribution and level of Mars' seismicity[END_REF][START_REF] Taylor | Seismic exploration of Mars and the NASA InSight mission[END_REF]) and from meteorite impacts [START_REF] Davis | Meteoroid impacts as seismic sources on Mars[END_REF][START_REF] Teanby | Seismic detection of meteorite impacts on Mars[END_REF][START_REF] Teanby | Predicted detection rates of regional-scale meteorite impacts on Mars with the InSight shortperiod seismometer[END_REF]. Although it is expected to be ≈100 times less seismically active than Earth, it is thought to be still seismically active today, as evidenced (for example) by recent tracks in dust leading to boulders in HiRISE images near known faults [START_REF] Roberts | Possible evidence of paleomarsquakes from fallen boulder populations, Cerberus Fossae, Mars[END_REF].

SEIS-SP

SEIS comprises of two independent, three-axis seismometers: an ultra-sensitive very broad band (VBB) oblique seismometer, and a miniature, short-period (SP) seismometer. Both are mounted on the LVL (a precision levelling structure) along with their respective signal preamplifier stages, within sealed, nitrogen-backfilled enclosures. Together, the VBB, SP, and LVL are deployed on the ground as an integrated package, isolated from weather by the WTS (wind and thermal shield). They are connected by a flexible cable tether to the EBOX (a set of electronic cards located inside the lander's thermal enclosure).

The SEIS-SP microseismometer represents the UK contribution to the InSight SEIS payload which is a three-axis microseismometer, using three microfabricated sensors which include monolithic in-plane silicon proof mass and folded-cantilever suspensions with electroplated coils and capacitive sensors which are connected through proximity electronics (via the tether) to the associated analogue feedback circuits within the EBOX. The SEIS-SP electronics are based upon circuits described in [START_REF] Pike | Design, fabrication and testing of a micromachined seismometer with nano-g resolution[END_REF][START_REF] Pike | A silicon microseismometer for Mars[END_REF], [START_REF] Pike | Fabrication process and package design for use in a micromachined seismometer or other device[END_REF], [START_REF] Gu | Solder pump for through-wafer vias[END_REF], [START_REF] Delahunty | Integrating solder bumpers for high shock applications[END_REF], [START_REF] Delahunty | Metal-armouring for shock protection of MEMS[END_REF].

Each sensor head (Fig. 1) consists of the sensor, a micromachined die package, surrounding magnetic assembly, proximity electronics and an encompassing hermetic package with connector to the tether harness: a 25 mm × 25 mm micromachined silicon die bonded to further structures to complete the displacement transducer, and protect the sensor from damage due to shock and vibration. The silicon die contains the moving part of SEIS-SP, a moving proof mass suspended within a stationary frame by a series of flexures. The proof mass is constrained to move in a single, compliant direction in the plane of the silicon die by the suspension design. The motion of the proof mass is measured capacitively using two sets of electrodes, a driving pair on the proof mass and a differential pick-up pair on a separate micromachined die, called the displacement transducer (DT). As the electrodes on the proof mass move with respect to the fixed electrodes on the DT, the overlap between the electrodes changes and the displacement can therefore be transduced. The magnetic assembly provides a static field perpendicular to the coil-current flow which together produces a Lorentz force to reposition the proof mass under electronic feedback control of the coil current. The assembly consists of a series of magnets, pole pieces and yokes to provide a closed magnetic circuit.

SEIS-SP is sensitive to frequencies between 0.05-40 Hz, recording data sampled at a rate of 100 Hz, and detects the acceleration of the ground along the sensitive axis along which the proof mass is free to move. The capacitive sensor detects the motion of the mass, operating under feedback control to increase the bandwidth and linearity, and to provide a velocity output. Thus, SEIS-SP does not model as a driven damped oscillator in the presence of an external force (such as those induced by the wind on the WTS) but rather as shown in Fig. 2.

On Mars, the deployed set-up will include the lander (lander body + solar panels, mounted on 3), the WTS (mounted on 3 feet) and SEIS (mounted on its 3 LVL feet). For the purposes of this work, the geometry of the lander, the WTS and SEIS instrument system are given in Fig. 3; the actual deployment geometry upon arrival on Mars is easily modifiable in the decorrelation pipeline. It is possible that the lander, WTS and SEIS may not be at the same vertical (z) height, as currently assumed.

It is usual that terrestrial seismometers are deployed in seismic vaults, in highly environmentally-controlled installations in which temperature, humidity, pressure and airflow are known and stable, and in this way the recorded seismic measurements have reduced levels of environmental noise. For deployments in which vaults are not possible, the general approach involves installation of the seismometer, recording and analysis of installationspecific noise characteristics, and adjustment of the set-up in order to lower the overall noise levels-for instance, by mounting on sand. It is also common to bury sensors or to enclose sensors with an upturned container, when vaults are not possible, in order to mimic vault-like conditions. However, for SEIS, neither vault nor terrestrial empirical installation routes are possible, meaning that environmental noise (as well as the sought seismic activity) will be a ubiquitous part of whatever measurements are taken by SEIS on the surface of Mars. On Earth, there are also seismometers mounted on the ocean bottom, which similarly to SEIS, are not enclosed by vaults; it is possible, given suitable environmental sensors mounted on future ocean-bottom seismometers, that the decorrelation technique developed here could be applied to such seismometers

TWINS and APSS

As usual for lander and rovers, InSight has a suite of environmental sensors, called the APSS (Auxiliary Payload Sensor Suite) to complement the main scientific payload (SEIS and HP3). The APSS consists of a magnetometer, a pressure sensor, and a wind and temperature sensor (TWINS). The pressure sensor will measure at 2 Hz in nominal mode and 20 Hz in high-rate mode; the wind and temperature sensors will measure at 0.1 Hz and 1 Hz in nominal and high-rate modes, respectively.

Environmental Noise Sources on SEIS

As SEIS will be deployed in the Martian surface, and though it will be protected from the environment by WTS, it is by no means isolated from the environment. Thus, it is expected that SEIS measurements will include "noise" from the environment (namely from pressure, temperature and wind fluctuations) which contaminate the seismic records sought by the mission science drivers. It is usual to capture data in a seismically-and environmentallyquiet period to establish a quiet baseline which can be used to isolate seismic signals alone; however, this method requires accurate knowledge that the baseline is indeed quiet, that the measurements are additive, and does not allow for knowledge of the intrinsic errors in seismometer nor environmental sensor data. The scheme presented here allows for any dataset to be inverted without need for a quiet baseline dataset for comparison, and incorporates measurement errors and correlations into the isolation of the seismic signal from the random noise and environmental noise sources. [START_REF] Murdoch | Evaluating the wind-induced mechanical noise on the InSight seismometers[END_REF] and [START_REF] Mimoun | The SEIS Team), The noise model of the SEIS seismometer of the InSight mission to Mars[END_REF] have formulated and validated a model deriving the forces imparted to the SEIS-VBB (and hence the SEIS-SP) by wind and temperature variations via various possible pathways; their model also includes a module on pressure variations and their pathways through to the seismometer signals from Lognonné and Mosser (1993). This model is implemented in this work as the forward model f, and is described qualitatively in Sects. 2.1-2.3; for details, please see [START_REF] Murdoch | Evaluating the wind-induced mechanical noise on the InSight seismometers[END_REF]. Additionally, the environmental noise forward modelling from [START_REF] Murdoch | Evaluating the wind-induced mechanical noise on the InSight seismometers[END_REF] implemented in this work has also been implemented in [START_REF] Clinton | Preparing for InSight: an invitation to participate in a blind test for Martian seismicity[END_REF].

Additionally, results from [START_REF] Teanby | Seismic coupling of short-period wind noise through Mars' regolith for NASA's InSight Lander[END_REF] on regolith inelastic coupling have been implemented in the modelling here, in conjunction with Murdoch's regolith model; their results show that only a factor of 0.02 of the force transfers through a Mars-analogue regolith, due to inelastic losses.

Finally, a single event detection module has been included and implemented, which will filter spurious readings stemming from single event radiation on the SEIS-SP electronics, which is discussed in more detail in Sect. 2.4.

Wind

The wind field interacts with anything large on the surface of Mars, causing such objects to create disturbances which can be felt by the LVL feet which are then measured by the SEIS instruments; it is not possible to tell from what source the incoming disturbance signal comes when it arrives at the LVL feet and ultimately is measured by SEIS.

Though there may be other large features (such as boulders) on the surface nearby the deployment site, it is certain that the WTS itself and the lander will be sources of windinduced disturbances. Sections 2.1.1 and 2.1.2 detail this pathway for WTS; it is identical for the lander.

Force of Wind on WTS

When wind is incident upon the WTS, it induces drag and lift forces onto the WTS. The WTS feet then exert forces on the regolith due to translation and rotation: translational motion will always happen, and for the purposes of this work, is assumed to be felt at the Centre of Pressure (CoP) of the WTS; rotational motion will result when the line of action of the force does not pass directly through the Centre of Rotation (CoR) of the WTS. These forces are then transferred to the regolith via the WTS feet-for the purposes of this work, we assume that each foot takes one third of the total force on the WTS as a whole.

Transfer Through Regolith

The regolith model assumes a fully elastic half-space (defined only for the positive vertical direction, along with all horizontal directions) with properties of the regolith assumed as in [START_REF] Murdoch | Evaluating the wind-induced mechanical noise on the InSight seismometers[END_REF]. Using a Green's tensor, the force from the disturbance from each WTS foot can be propagated through the regolith to the LVL feet, such that the ground deformation as a function of time can be estimated. The vertical component of this force then results in each LVL foot causing an indentation in the regolith from which it is then possible to calculate the overall change in disturbance in the position of each of the LVL feet which is measured by SEIS-both by VBB and SP-as an acceleration from one timestep to the next registered by the seismometers.

However, results from [START_REF] Teanby | Seismic coupling of short-period wind noise through Mars' regolith for NASA's InSight Lander[END_REF] confirm that the regolith cannot be modelled as fully elastic, and suggests that a fraction f inelas = 0.02 of the disturbance is actually transferred, due to inelastic losses.

Due to Lander

As discussed in Sects. 2.1.1-2.1.2, disturbances due to wind acting on the lander are also included: the mathematical derivation of which can be assessed by substituting lander parameters instead of WTS parameters throughout.

Temperature

Fluctuations in temperature cause disturbances in the registered signals measured by SEIS-SP. Firstly, the WTS and SEIS-SP act as a thermal filters with associated thermal transfer functions: any fluctuation in environmental temperature will pass through the filtering effect of the WTS, and then through the SEIS-SP box (which also acts as a filter to the temperature fluctuation) before encountering the SEIS-SP sensors themselves; this temperature effect is independent of the SEIS-SP transfer function. Secondly, environmental temperature fluctuations, having passed through the filter of the WTS protection, will also cause the individual LVL legs to expand/contract via thermoelastic expansion; this will cause a tilt in the plane of the three SEIS-SP axes, which will cause the components of the signal to be altered primarily as acceleration on the horizontal axes.

Pressure

Pressure fluctuations cause propagated accelerations caused by tilting into the SEIS-SP measurement chain. This is modelled using the low-frequency oscillation solution of Lognonné and Mosser (1993); the range of oscillations expected on Mars have been considered, and even for dust devils, the low-frequency solution should hold here.

Single Events

Single event phenomena (SEP) are caused when very energetic heavy ions (any ion having an atomic number Z > 1) and protons pass through semiconductor materials and deposit a charge. If this direct or indirect (through spallation) charge is sufficient, then a Single Event Effect-either destructive (Single Event Latchup SEL, Single Event Burnout SEB, Single Event Gate Rupture SEGR, etc.) or non-destructive (Single Event Transient SET, Single Event Upset SEU, etc.)-can occur. Galactic cosmic rays and solar particle events are sources of such ions/protons; and hence are concern for the InSight electronics both in cruise to Mars as well as in its deployed phase.

SEIS-SP's electronics are chosen such that it operates within the performance specification during and after exposure to the high-energy radiation environments as detailed by the mission's environmental specifications, with a Radiation Design Factor (RDF) of 1, such that:

-Temporary loss of function or loss of data shall be permitted provided that the loss does not compromise the subsystem/system's health, that full performance can be recovered rapidly, and that there is no time in the mission that the loss is mission-critical. -Normal operation and function shall be restored via internal correction methods without external intervention in the event of an SEU. -Fault traceability shall be provided in the telemetry stream to the greatest extent practical for all anomalies involving SEEs.

A transient single event effect on a particular component will result in a spurious signal in the output from the axes channel in which the component corresponds in the feedback electronics. The problem is whether these transients will cause signals which are differentiable from potential seismic signals. It is expected that SEP will occur a maximum of O(1e -2 ) events/day, and more typically at a rate of O(1e -5 ) events/day-so it is expected to be very rare that spurious signals from SEPs will be registered in SEIS-SP signals.

Unless one seismometer axis is perfectly aligned with the tilt of the seismometer-the probability of which is virtually nil-the transient can only affect one axis at a time. From an operational perspective, then, simple comparison of all three SEIS-SP axis channels together will identify what signals are seismic (i.e. those for which a disturbance in signal is seen in all three channels) and those which are a result of SEU transients/effects (i.e. disturbance in only one channel). Experience in terrestrial seismometry indicates that this is an effective way in which to identify spurious signals in different seismometer channels.

Thus, the SEP filtering module acts simply to run through each timestep of data signal y i measured by each of the three SEIS-SP axes i, comparing the timestep j with that immediately before, and using a simple threshold y thresh , identifying large signals which are only evident in one axis:

if y i (t j ) -y i (t j -t j -1 ) > y thresh and |y k (t j ) -y k (t j -t j -1 )| < y thresh |y m (t j ) -y m (t j -t j -1 )| < y thresh (1)
whereby k, m = i and k = m, then there has been a SEP on axis i at timestep t j . This timestep should then be disregarded from further analysis (Fig. 4). This assertion is dependent upon the effect of all the SEPs being non-damaging to the components considered. It is thus worth assessing and monitoring the axes upon which an SEP has been identified to ensure that there is no change from its historical behaviour from timestep t j onwards.

Fig. 4 Example of simulated SEP overlaid on top of SEIS-SP signal for SEIS-SP axes H1 in blue, H2 in green, V in red, and with SEP-detection module's identification of timing of detected SEP events and attribution of which axis is affected in larger black circles. Percentage of correct diagnoses per case (axis or no SEP event) highlighted by right vertical axis

Decorrelation Scheme: Optimal Estimations Inversion

Using the forward model f described in the previous sections, an iterative optimal estimation method is employed to decorrelate the environmental signals from the total signal measured by the SEIS-SP sensors, in order to isolate the signal attributable only to seismic activity alone. This methodology is an application of that used in atmospheric science, whereby satellite spectra containing information upon temperature, pressure, chemical composition and any atmospheric phenomena such as aerosols and clouds are inverted in order to determine information on any one or more of these individual quantities. Due to the high recorded sampling rate of SEIS-SP measurements f SP (of 100 Hz), inversions are run in discrete timebands t (usually around 2 s) consisting of tf SP individual timestep measurements (usually around 200) in the time-domain in order to keep array sizes manageable. The inversion scheme outlined below is applied for each defined timeband in the input SEIS-SP signal, and the results are concatenated together to provide the whole timeseries. This is based on the schemes in [START_REF] Rodgers | Inverse Methods for Atmospheric Sounding: Theory and Practice[END_REF].

Inverted Parameters

The sought result of the decorrelation is the state vector x in typical optimal estimations terminology: for the case at hand, it is the force attributable to the seismic signal alone, in xyz coordinates, and is defined as

x = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ F(x 1 ) F(x 2 )
. . .

F(x N ) F(y 1 ) F(y 2 ) . . . F(y N ) F(z 1 ) F(z 2 ) . . . F(z N ) ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (2)
for a timeband having N individual time-sampled measurements therein indexed by j , for force F consisting of three perpendicular cartesian x j , y j , and z j components. It should be noted that in the assumed deployment geometry, that these xyz coordinates correspond to north, east and vertical directions. The measurement vector (from which the state vector is derived), then, is the signal S measured directly from the three SEIS-SP axes, in SEIS-SP coordinates (H1, H2, V):

y = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ S(H 1 1 ) S(H 1 2 ) . . . S(H 1 N ) S(H 2 1 ) S(H 2 2 ) . . . S(H 2 N ) S(V 1 ) S(V 2 ) . . . S(V N ) ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (3)
for a timeband having N individual time-sampled measurements therein indexed by j , for signal S consisting of three SEIS-SP (non-perpendicular) axes components.

It should be noted here that whilst forces are additive, signals are not as they have been convolved through the SEIS-SP transfer function.

The Jacobian, K, describes how the forward model f varies with respect to the state vector ∂f ∂x -and essentially advises the step-direction of the iterative inversion processcalculated in this implementation by finite differences. It captures the change in signal at times i through j when each timestep i through j in the retrieval vector x force is perturbed by x, such that ∂x i = (x 1 + x)x i . Here:

K = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ∂f H 1 1 ∂x 1 ∂f H 1 2 ∂x 1 . . . ∂f H 1 N ∂x 1 ∂f H 2 1 ∂x 1 ∂f H 2 2 ∂x 1 . . . ∂f H 2 N ∂x 1 . . . ∂f V 1 ∂x 1 ∂f V 2 ∂x 1 . . . ∂f V N ∂x 1 ∂f H 1 1 ∂x 2 ∂f H 1 2 ∂x 2 . . . ∂f H 1 N ∂x 2 ∂f H 2 1 ∂x 2 ∂f H 2 2 ∂x 2 . . . ∂f H 2 N ∂x 2 . . . ∂f V 1 ∂x 2 ∂f V 2 ∂x 2 . . . ∂f V N ∂x 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ∂f H 1 1 ∂x N ∂f H 1 2 ∂x N . . . ∂f H 1 N ∂x N ∂f H 2 1 ∂x N ∂f H 2 2 ∂x N . . . ∂f H 2 N ∂x N . . . ∂f V 1 ∂x N ∂f V 2 ∂x N . . . ∂f V N ∂x N ∂f H 1 1 ∂y 1 ∂f H 1 2 ∂y 1 . . . ∂f H 1 N ∂y 1 ∂f H 2 1 ∂y 1 ∂f H 2 2 ∂y 1 . . . ∂f H 2 N ∂y 1 . . . ∂f V 1 ∂y 1 ∂f V 2 ∂y 1 . . . ∂f V N ∂y 1 ∂f H 1 1 ∂y 2 ∂f H 1 2 ∂y 2 . . . ∂f H 1 N ∂y 2 ∂f H 2 1 ∂y 2 ∂f H 2 2 ∂y 2 . . . ∂f H 2 N ∂y 2 . . . ∂f V 1 ∂y 2 ∂f V 2 ∂y 2 . . . ∂f V
. . . ∂f V 1 ∂y N ∂f V 2 ∂y N . . . ∂f V N ∂y N ∂f H 1 1 ∂z 1 ∂f H 1 2 ∂z 1 . . . ∂f H 1 N ∂z 1 ∂f H 2 1 ∂z 1 ∂f H 2 2 ∂z 1 . . . ∂f H 2 N ∂z 1 . . . ∂f V 1 ∂z 1 ∂f V 2 ∂z 1 . . . ∂f V N ∂z 1 ∂f H 1 1 ∂z 2 ∂f H 1 2 ∂z 2 . . . ∂f H 1 N ∂z 2 ∂f H 2 1 ∂z 2 ∂f H 2 2 ∂z 2 . . . ∂f H 2 N ∂z 2 . . . ∂f V 1 ∂z 2 ∂f V 2 ∂z 2 . . . ∂f V
. . . ∂f V 1 ∂z N ∂f V 2 ∂z N . . . ∂f V N ∂z N ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (4)

Errors

Errors on the measurements y are accounted for in the measurement covariance matrix S y which is, in this case, a 3N -by-3N diagonal matrix having diagonal elements equal to the square of the standard deviation of the random noise, σ y , believed to occur on each SEIS-SP signal timestep, such that

S y = ⎛ ⎝ σ 2 y H 1 I N . . . . . . . . . σ 2 y H 2 I N . . . . . . . . . σ 2 y V I N ⎞ ⎠ (5)
where I M is the N -by-N identity matrix. The value of σ y H 1 and σ y H 2 , and σ y V for the SEIS-SP design parameters are 1e -9 m/s 2 /Hz 1 2 and 5e -8 m/s 2 /Hz 1 2 , respectively, estimated at 1 Hz.

Iterative Formulation and Error Propagation

An iterative form of optimal estimation is used in order to numerically converge to a solution. Starting at i = 0 with the first guess x 0 = a, the iterative state vector is defined as

x i+1 = x i + (1.0 + γ )S -1 a + K T i S -1 y K i -1 K T i S -1 y (y -f i ) + S -1 a (x i -a) (6)
where S a is the a priori covariance which is described in detail in Sect. 3.4, and γ is the Levenberg-Marquardt parameter which can take any finite value and that acts to adjust the iterative "step-size" in order to maintain monotonic progress towards the solution of lowest least squares error [START_REF] Rodgers | Inverse Methods for Atmospheric Sounding: Theory and Practice[END_REF]. This least squares error is a metric of the difference between the measurements and the estimate of the measurements via the forward model for each iteratively determined inversion vector, and is quantified by the χ 2 error, defined as

χ 2 i+1 = 1 N (y -f i ) T S -1 y (y -f i ) + (x i -a)S -1 a (x i -a) (7)
which should approach unity for a fit that is consistent with the measurements. The error in the inverted parameters themselves thus calculated is expressed as the error covariance matrix

S x S x = K T i S -1 y K i + S -1 a -1 . (8)
Typically, the difference between two iterative values of χ 2 determines when the optimal estimation inversion is said to converge-when χ 2 i+1χ 2 i → 0, the inversion is said to have converged-here when

χ 2 i+1 -χ 2 i χ 2 i × 100% ≤ 0.1%.

A Priori

As it is an optimal estimation scheme that is implemented, an a priori state vector a must be provided to bound the possible inverted values-and in this implementation, this a priori is taken as the first guess of the inverted state vector x in the iterative inversion process.

A pessimistic a priori could be taken that there is no seismic activity (thus a = 0); however, inversion of the measurements can provide a better starting point.

A Priori Preliminary Inversion

When starting a decorrelation, the analyst has the SEIS-SP measured signal and the TWINS/APSS metrological data available. However, he does not know how much of the measured signal is attributable to seismic vs. aseismic components-and in any case, the seismic and aseismic sources are not additive. Unlike the signals, the measured force and attributable environment forces are additive, and can be subtracted in order to yield the a priori force attributable to the seismic activity alone:

a = y force -x env (9)
However, y force and x env are not readily available without some estimation.

Taking the environmental pseudo-a priori data (wind, pressure and temperature) from APSS/TWINS, and running these through the forward model f provides the force attributable to the environment factors alone in xyz coordinates x env .

The measured signal y in SEIS-SP-axes coordinates (H1, H2 and V) can then be used in a preliminary inversion to convert the signal into measured force in xyz coordinates y force . In this preliminary inversion, which is run for each SEIS-SP axes individually, the signal as a function of time for each SEIS-SP axes is converted into the force for each SEIS-SP axes x prel (t) as a function of time. This is accomplished by defining a preliminary inversion forward model f prel such that

y = x(t) * T (t) (10)
whereby y(t) is the signal measured by each SEIS-SP axes, which is defined by the force as a function of time x(t) convolved with the complex SEIS-SP transfer function T (t) which itself has been converted from its frequency domain T (s) into the time domain T (t). This preliminary inversion uses the same iterative formulation defined in Eq. ( 6); it thus requires an preliminary a priori a (t) which is taken here to be 1, and error estimates on both the preliminary a priori and the measurement used in the preliminary retrieval (y prel (t) = y(t) from the main decorrelation inversion) which are taken to be σ x prel = 1, σ y prel = 0.1. It also requires preliminary covariance matrices S a prel = σ 2

x prel 1 and S y prel = σ 2

x prel 1. This iterative preliminary retrieval is then let run until it converges for the inverted state vector x prel , which is said to happen when χ 2 is below 0.001 or the difference between two consecutive iterative χ 2 values is less than 0.001.

Once inverted, the estimated forces x prel (t) are combined into vector form (whereby there are N values for each of the 3 SEIS-SP axes H1, H2 and V), and rotated from SP axes coordinates (H1, H2 and V) to xyz coordinates such that

y force = ⎛ ⎜ ⎝ x prel H 1 x 2 prel H 1 -x 2 prel H 2 cos 2 θ/ sin 2 θ x prel V ⎞ ⎟ ⎠ (11)
whereby θ = 60 • . This y force and x env can then be supplied to Eq. ( 9) to estimate the a prior for the decorrelation.

Implementation of a Priori in Decorrelation

The inversion is allowed to vary from the a priori forces in a manner proportional to the uncertainty σ a in a, which is specified for each of the xyz coordinates. Typically, σ a is taken to be the variance of expected inverted values; in this case, this is defined for each of the three xyz coordinates: σ ax = 3e -8 m/s 2 /Hz 1 2 , σ ay = 6e -8 m/s 2 /Hz 1 2 and σ az = 2e -7 m/s 2 /Hz 1 2 , estimated at 1 Hz. This is accounted for in the 3N -by-3N a priori covariance matrix S a , which includes correlations between different inverted variables-and to aid in obtaining a physically smooth inverted state vector, correlations between each sub-timestep within an inversion timeband are assumed via a linear correlation in time between 0 and 1 (Fig. 5):

S a ij = ⎛ ⎜ ⎜ ⎝ (1 -|j -i| N )σ 2 ax I N . . . . . . . . . (1 -|j -i| N )σ 2 ay I N . . . . . . . . . (1 -|j -i| N )σ 2 az I N ⎞ ⎟ ⎟ ⎠ (12)
Because the correlation is assumed over the timeband-which itself is generally taken to be about 1 s-but the data is taken at 100 Hz, this may limit the frequency content that can be recovered to lower-frequency seismograms. Thus, care should be taken to match the S a to the kind of signal that is being recovered in order to best capture the frequency of the noise and signal.

Application to Simulated Data

The performance and behaviour of optimal estimations decorrelation schemes are best tested using simulated data for which the "true" state is known, and in this case, the various simulated (and hence known) environmental noise sources can be added on in a controlled manner using the forward model. The applicability of the inversion to real data and assumed performance, of course then assumes that the forward model does a sufficient and representative job of representing reality, which cannot be quantified in these simulation tests, but is assessed in [START_REF] Murdoch | Evaluating the wind-induced mechanical noise on the InSight seismometers[END_REF].

To this end, simulated meteorological conditions as shown in Fig. 6 have been generated using representative magnitudes and timescales of variability from [START_REF] Murdoch | Evaluating the wind-induced mechanical noise on the InSight seismometers[END_REF]. The simulated seismic velocity measurement is taken from [START_REF] Teanby | Seismic detection of meteorite impacts on Mars[END_REF]; representative of an explosive source at 50 km depth, rendering a total moment of 10 17 Nm (magnitude 5.3) at a location 60 • (≈ 3000 km) away from the assumed SEIS-SP location. This has then been resampled at the SEIS-SP frequency (100 Hz), and convolved with the SEIS-SP transfer function in order to simulate what SEIS-SP would measure in a noise-free and isolated set-up (Fig. 7 top two rows).

Detectability and Sensitivities to Environmental Parameters

The ability of the decorrelation scheme to identify and remove various components of signal is dependent upon the sensitivity of the measurements to show differential in signal in each axis above the level of expected noise; if the change in signal for a particular parameter is greater than this noise level, then the signal due to the parameter is detectable and removable. In standard optimal estimations theory, this is captured by the Jacobian matrix K.

Figure 8 shows these sensitivities (i.e. change in signal corresponding to a change) for seismic forces in each of the xyz directions independently, for wind speeds in each of the xyz directions independently for both the WTS and lander interactions, for temperature fluctuations as fedthrough both the SEIS-SP transfer function and through thermoelastic tilt routes, and for fluctuations in pressure. Generally:

-SEIS-SP axes are very sensitive to the correlated xyz seismic signal, with some crossover in the horizontal axes as expected; -All SEIS-SP axes are sensitive to the effects of wind on the WTS for wind speeds greater than ≈ 1 m/s; -All SEIS-SP axes are sensitive to the effects of wind on the lander for wind speeds greater than ≈ 2-3 m/s; -The horizontal SEIS-SP axes are sensitive to all expected temperature fluctuations; -Temperature and pressure fluctuations are not detectable through the SEIS-SP transfer function.

For the seismic example shown in Fig. 7 and the environmental series shown in Fig. 6, the components of signal attributable to each environmental noise source as a function of time are shown in Fig. 9. Note that whilst components of force are additive, the components of signal (which are convolved with the complex SEIS-SP transfer function) are not.

Extraction of Seismic Signal from Noisy Simulations

The simulated signal used to assess the overall performance and ability of the optimal estimations decorrelation is obtained by running the forward model with the simulated meteorological conditions from Fig. 6 to obtain the individual force components (xyz) attributable 9 Example simulation: the contribution to the total signal measured by each SEIS-SP axis H2, and due to the wind on the WTS (first row) and lander (second row), the temperature fluctuations (third row), the thermoelastic tilt (fourth row), the pressure variation (last row). The wind, temperature and pressure time-series used in this example simulation are shown in Fig. 6 to each environmental noise source, combining these with the seismic force (xyz) components, and then convolving these forces through the SEIS-SP transfer function to obtain the overall signal registered in each of the SEIS-SP axes (H1, H2, and V). Next random noise is added on top of each of the SEIS-SP axes signal time-series, with magnitude of each random distribution of noise corresponding to σ y for each axis. The blue lines in the first and third rows of Fig. 10 show the simulated "noisy" velocity (i.e. force) in xyz and signal in H1/H2/V axes; the red line shows the "pure" signal-that is, the seismic "truth" signal for which the decorrelation is ultimately aiming from starting the optimal estimations inversion at the "noisy" input signal and the meteorological data.

The optimal estimations scheme is then applied to this "noisy" input signal, and the inversion allowed to iterate until convergence is achieved, over the full 300 s timeseries in discrete timebands of 2 s; the results of which are shown in Fig. 10 in both the velocity and signal time-domains. Looking at the difference plots in the second and fourth row, the impact of the environmental noise (as opposed to the noise-free "pure" signal) is clearly Fig. 10 Result of application of optimal estimations scheme to example seismic signal contaminated with example environmental inputs: top row shows the inverted (black) and noisy (blue) velocity in xyz axes; second row shows the difference between the inverted and pure seismic-signal-along velocity (black), between the noisy and pure seismic-signal-along velocity (blue) in the xyz axes; third and fourth row show the same plots, but in the signal measured in the SEIS-SP H1, H2 and V axes; the final row shows the χ 2 error (left panel), the number of iterations required for convergence (middle row) and the convergence status (right panel; 0 = unconverged, 1 = converged) in each of the 2 s-long timebands in which the total time series has been inverted visible-for instance, it is a nearly uniform offset of 40 nm/s in velocity in the x direction (this is the dominant wind impact), with more structure visible in the H1 SP axes signal plot. Looking at the velocities and signals inverted from the optimal estimations scheme, the difference between the "pure" and the inverted noise-free signal is largely zero, with a couple of notable exceptions when there is information missing on the horizontal axes (x/y vs. H1/H2) which are not always distinguishable. Fig. 11 shows a only the x and SP-H1 axis, over the first 30 s, for ease of viewing. Thus, it is clear that the optimal estimations inversion is capable of differentiating between environmental noise and seismic signals, assuming that the forward model properly represents the environmental forces. Looking at convergence parameters, for this simulated example, the scheme converges within two iterations, with a low cost function value, as quantified by the χ 2 value.

Sensitivities and Operational Considerations

Environmental noise is present continuously, and is ubiquitous on any signal recorded. However, the SEIS-SP signal is recorded at a particular frequency, the various environmental monitoring signals at other frequencies; and the scheme's ability to decorrelate the environmental signal will be influenced by the relationship between these as well as the frequency of fluctuation of the environmental parameters themselves. In theory, the highest-possible Fig. 11 As for Fig. 10, but just for the x and SP-H1 axis for the first 30 s of data frequency data of concurrent SEIS-SP and APSS/TWINS measurements would give the best decorrelation; however, this is not the design case (SEIS-SP has nominal frequency of 100 Hz; APSS has nominal frequency of 2 Hz and high-rate mode of 20 Hz; TWINS has nominal frequency of 0.1 Hz and a high-rate mode of 1 Hz), nor is it operationallyreasonable to expect that such a high data rate is possible continuously. From a processing perspective, it is also advantageous to use the fewest number of datapoints (i.e. the lowest frequency) to increase the speed of the processing pipeline which decorrelates the environmental signal from the seismic signal.

During the commissioning phase when the lander is first deployed on the Martian surface, both high-rate and nominal mode APSS and TWINS measurements will be taken during expected noisy and less-noisy periods (day and night, respectively), which will be used to assess how the decorrelation performs with real data and environmental inputs.

Figure 12 shows the impact of varying the frequency at which the input APSS/TWINS pressure, temperature and wind-speed measurements are provided to the optimal estimation scheme on the inverted performance. In this investigation, environmental noise is applied via the forward model to each 100 Hz SEIS-SP frequency record; this yields the simulated Fig. 12 Difference between the inverted and pure simulated noise-free velocity timeseries for example simulation when environmental information (namely temperature information, pressure information and wind information) is provided at different frequencies. Each panel shows a different combination of frequency of temperature information (F T ) for varying frequency of pressure information (F P , each column) and frequency of wind information (F W , each row), with the nominal operational mode frequencies for each highlighted in blue and high-frequency mode frequencies highlighted in red time-series. Next, the scheme is run with pressure, temperature and wind-speed information provided at varying frequencies between 0.01-100 Hz. The frequency at which pressure is input within this range makes no difference to the performance of the inversion; similarly the rate of temperature information makes little difference; the inversion performance improves with increased frequency of the input wind-speed information, but makes little difference between nominal and high-rate TWINS modes. This implies that for environmental fluctuations having expected frequencies (i.e. not high-frequency perturbations like dust-devils), that both the nominal and high-rate modes should be sufficient to decorrelate the environmental noise from the seismic signal.

Application to Viking Data

The Viking Lander 2 seismometer captured data on Sol 80 at 03h00 of the Viking mission which has been asserted to be seismic in origin (of estimated magnitude 2.7), given interpolation on wind, pressure and temperature data taken around the time period of the event [START_REF] Lorenz | A bump in the night: wind statistics point to Viking 2 Sol 80 seismometer event as a real marsquake[END_REF]. This represents a good test case of a real Marsquake in the presence of environmental signal contamination, upon which to apply the scheme developed here.

The Viking 2 lander seismometer data was captured at a frequency of 1.01 Hz in 3 perpendicular axes; the metrology (wind, pressure and temperature) data was taken much less frequently with the closest data taken before the event at 02h40 and after the event at 03h45.

Fig. 13

Result of application of optimal estimations scheme to Viking Sol 80 signal, which implicitly contains environmental noise sources: top row shows the inverted (black) and noisy (blue) velocity in xyz axes; second row shows the difference between the inverted and Viking velocity (blue); third and fourth row show the same plots, but in the signal measured in the SEIS-SP H1, H2 and V axes; the final row shows the χ 2 error in each of the 2 s-long timebands in which the total time series has been inverted

For this experiment, the seismometer and metrology data is interpolated onto the SEIS-SP frequency, and the seismometer data convolved through the SEIS-SP transfer function. There is no further environmental nor random noise added to the seismometer data since it already includes components of both these errors sources. The optimal estimations scheme is then applied using these interpolated measurements as inputs, and the SEIS-SP forward model in order to isolate the portion of signal attributable to seismic sources alone-it is worth noting that because of the SEIS-SP parameters used here, that this will yield only a rough estimate of whether the signal measured is attributable to environmental sources or rather that it requires a seismic component of input in order to reach the signal measured. Figure 13 shows the results of this experiment; whilst there are clearly artefacts stemming from instrument-specific and regolith assumptions on the inverted values, the decorrelation scheme suggests that the large signal features are indeed attributable to non-environmental sources, taking out small magnitude components of signal which are probably attributable to environmental fluctuations.

Conclusions

An optimal estimations inversion has been developed using the forward model developed in [START_REF] Murdoch | Evaluating the wind-induced mechanical noise on the InSight seismometers[END_REF] in order to decorrelate seismic components of InSIght SEIS-SP signal from environmental and random noise. Using simulated data and expected design noise levels, the optimal estimations scheme is successful in detecting and removing signals from wind and temperature fluctuations, but pressure fluctuations are below the noise level. Application of the scheme to Viking data implies that the Viking 2 seismic event on sol 300 was indeed likely a seismic event, assuming all geophysical parameters and the forward model are representative. Analysis shows that nominal mode measurements are sufficient to decorrelate environmental signals from seismic, and that high-resolution mode measurements are not required except to invert low windspeeds or small temperature variations.

A SEP-identification module has been implement as a preliminary filter to the decorrelation scheme; this has also proven capable of identifying timesteps in continuous data with single event contamination; a requirement of the InSight environmental requirements pertaining to radiation. This will be implemented as part of the operational processing; it is, however, expected that the raw data will also be provided to end users.

Future work includes field testing of test units of SEIS-SP packages modified for terrestrial use (i.e. g earth vs. g mars ) in a representative unisolated installation in order to test the representativeness of the modelling and assumptions made. The scheme will also be applied to the simulated dataset in [START_REF] Clinton | Preparing for InSight: an invitation to participate in a blind test for Martian seismicity[END_REF], which can be upsampled form 2 Hz.

Fig. 2 Fig. 3

 23 Fig. 2 The complex transfer function T (s) for SEIS-SP in frequency space (s) (top = real component, bottom = imaginary component) as a function of frequency

Fig. 5

 5 Fig. 5 Example a priori covariance matrix S a for H1 over a 2 second timeband, showing correlation (0 to 1) between timesteps

Fig. 6

 6 Fig. 6 Simulated meteorological conditions for example: top three panels are components of wind in xyz coordinates, fourth panel is temperature at ground level outside of the WTS, bottom panel is the pressure at ground level. Please note that example is 300 s long, and sampled at 100 Hz

Fig. 7

 7 Fig. 7 Example simulation: the motion of Martian regolith in each xyz axis (top row), and the contribution to the total signal measured by each SEIS-SP axis (H1, H2, and V) due to the seismic motion (bottom row)

Fig. 8

 8 Fig. 8 Sensitivity of the SEIS-SP signal in SP axes (H1, H2, V; vertical axis of each panel) to strength of an imposed seismic signal applied in xyz direction (top row, each panel a different xyz axis), and to varying environmental conditions: second row is due to wind on the WTS (each panel for wind in a different xyz direction), third row is due to wind on the lander (each panel for wind in a different xyz direction), fourth row left panel is due to temperature fluctuations, middle panel is due to thermoelastic tilt and right panel is due to pressure fluctuations

Acknowledgements

We acknowledge the UKSA Aurora framework for supporting the principal and several co-authors.


Yahyaoui Khadidja
email: yahyaouikhadidja@yahoo.fr

Bouri Abdenour
email: bouriabdenour17@yahoo.fr

Scheduling in Real-Time Systems Using Hybrid Bees Strategy

Keywords: Real-time systems, scheduling, HBMO, Optimization, Metaheuristics

In the last decade, stochastic and meta-heuristic algorithms have been extensively used as intelligent strategies to resolve different combinatorial optimization problems. Honey Bee Mating Optimization is one of these most recent algorithms, which simulate the mating process of the queen of the hive. The scheduling algorithm is of paramount importance in a real-time system to ensure desired and predictable behavior of the system. Within computer science real-time systems are an important while often less known branch. Real-time systems are used in so many ways today that most of us use them more than PCs, yet we do not know or think about it when we use the devices in which they reside. Finding feasible schedules for tasks running in hard, real-time computing systems is generally NP-hard. In this work, we are interested in hybridizing this HBMO algorithm with other metaheuristics: Genetic Algorithms (GA), Greedy Random Adaptive Search Procedure (GRASP), Tabu Search (TS) and Simulated Annealing (SA) to resolve a real-time scheduling problem and obtain the optimal tasks schedule with respecting all temporal contraints. This is a complex problem which is currently the object of research and applications. In this scheduling problem, each task is characterized by temporal, preemptive and static periodicity constraints. The quality of the proposed procedure is tested on a set of instances and yields solutions which remain among the best.

Introduction

Finding feasible schedules for tasks running in hard, real-time computing systems is generally NP-hard [START_REF] Azerou | Scheduling Algorithmsfor Real-Time Systems[END_REF]. A vast amount of work has been done in the area of real-time scheduling by both operations research and computer science communities. The proposed approach is based on an optimization algorithm inspired from the Honey Bees Mating Optimization (HBMO). It combines an HBMO algorithm, Genetic Algorithms (GA), Simulated Annealing (SA), Tabu Search (TS) and the Greedy Randomized Adaptive Search Procedure (GRASP). The scheduling algorithms can be rated based on different parameters like makespan, flowtime, communication cost, reliability cost, and makespan [START_REF] Devi | Multiprocessor Scheduling of Dependent Tasks to Minimize Makespan and Reliability Cost Using NSGA-II[END_REF]. The task duplication scheduling scheme has been utilized to minimize the processor idle time and lessen the mean task response time in [START_REF] Singh | An Optimized Task Duplication Based SchedulinginParallelSystem[END_REF]. In this paper, we aim to obtain a schedule called feasible if both of the following conditions are met 1) No precedence constraints are violated and 2) all temporal constraints are satisfied. In order to realize this goal, we have used a new approach based on the hybridization of several algorithms. Our strategy uses the HBMO as a base algorithm which combines a number of different procedures. Each of them corresponds to a different phase of the mating process of the honey bees. In the next section, we see also that the honey bee mating algorithm has never been applied to resolve real-time tasks scheduling problems. This prompted us to use this strategy and hybridizing it with other metaheuristics in order to see its performance to resolve this kind of problem. Our contribution consists in combining the global search represented by HBMO with some local search metaheuristics for the resolution of realtime task scheduling problems in order to intensify the search in promising zones detected by the HBMO exploration process. This hybridization is realized with methods which have been applied in isolation in the resolution of real-time task scheduling problems: GA, TS, SA and GRASP [START_REF] Fleischer | Simulated Annealing: Past, Present, and Future[END_REF][START_REF] Casey | GRASP the examination scheduling problem[END_REF] First, we have to select the population of honey bees that will make up the initial hive. In the proposed approach, this population is created by using GA alone, or GRASP alone, or random selection selected. Then a percentage is selected from these three populations using a hybridization of these processes. The best member of each initial population obtained is selected as the queen of the hive. All the other members of each population are the drones. The two local search strategies SA and TS are used as workers for improving the broods. The paper is organized in the following way. The related works are presented in the second section. The presentation and formulation of real-time tasks scheduling problems is described in Section 3. Section 4 presents the proposed Scheduling Approach. The resolution process of the scheduling real-time tasks problem is detailed in Section 5. The sixth section presents the simulation and experimental results arrived at. The paper finishes with the conclusion and recommendations for future research.

Related works

The task scheduling problem is defined as one of the many popular academic NPhard problems. In [START_REF] George | Dynamic Task Scheduling Methods in Heterogeneous Systems-A Survey[END_REF], different aspects in scheduling and issues in various levels of real time systems are discribted. The use of metaheuristic methods shows their efficiency and effectiveness to solve these categories of complex problems [7]. In the literature, many metaheuristics have been proposed based on methods and approaches to task scheduling.. In [8], E. G. Talbi et al. have hybridized Simulated Annealing (SA) with Genetic Algorithms (GA). The authors in [9] resolve the scheduling task problem through hybridization of Particle Swarm Optimization (PSO) with GA.

The comparison between various scheduling algorithms and simulations results is given in [START_REF] Yoo | Multiobjective GA for Real Time Task Scheduling[END_REF]. An scheduling algorithm has multiobjective to minimize the total tardness and total number of processors is used in [START_REF] Konar | An Improved Hybrid Quantum-Inspired Genetic Algorithm (HQIGA) for Scheduling of Real-Time Task in Multiprocessor System[END_REF]. The authors combine Adaptive Weight Approach with genetic algorithm and simulated annealing method. In order to explore the computing power of quantum computation, an approach based on the hybrid Quantum genetic algorithm is used to resolve the real-time scheduling problem in multiprocessor environment In [12]. In [START_REF] Koudil | Using Artificial Bees to Solve Partitioning and Scheduling Problems in Codesign[END_REF], the authors applied Bee Colony Optimization (BCO) as a heuristic algorithm to solve the problem of static scheduling of independent tasks on identical machines. This strategy is based on the intelligent behavior of honey bees in the foraging process. Koudil et al. adapted the MBO algorithm to solve integrated partitioning/scheduling problems in co-design in [START_REF] Nirmala | Fuzzy Scheduling Algorithm for Real -Time multiprocessor system[END_REF]. This algorithm gives good results in terms of solution quality and execution time. The Fuzzy logic is also used to solve the real time tasks scheduling problem. In [START_REF] Monnier | A genetic algorithm for scheduling tasks in a real-time distributed system[END_REF] , the autors propose a model witch in the input stage consists of three lin-guistic variables i.e. CPU time, deadline and communication overhead.

Partitioning and scheduling

A real-time application is normally composed of multiple tasks with different levels of criticality : Missing deadlines is not desirable. In this context, three categories of realtime applications are defined [START_REF] Devi | Multiprocessor Scheduling of Dependent Tasks to Minimize Makespan and Reliability Cost Using NSGA-II[END_REF]16]:

 Soft real-time applications : The system could still work correctly Although the tasks could miss some deadlines. However, missing some deadlines for soft realtime tasks will lead to paying penalties.  Hard real-time applications : The tasks can not miss any deadline, otherwise, undesirable or fatal results will be produced in the system.  Firm real-time applications : The tasks, which are such that the sooner they finish their computations before their deadlines, the more rewards they gain.

In this paper, a real-time task 𝑡 𝑖𝑗 graph is defined by a set of temporal constraints such as { 𝑠𝑡(𝑡 𝑖𝑗 𝑘 ), 𝑐(𝑡 𝑖𝑗 𝑘 ), 𝑓𝑡1(𝑡 𝑖𝑗 𝑘 ), 𝑑𝑙(𝑡 𝑖𝑗 𝑘 ) } as well as preemption and static periodicity constraints. The table 1 reports the temporal contraints or parameters description of real-time system to study.

Table 1 Description of the real-time tasks adopted Temporal parameters

Temporal

Parameters description

Temporal parameters definition

𝑠𝑡(𝑡𝑖𝑗 ) )

Start time Time at which a j ieme task of i eme graph starts its execution 𝑐(𝑡𝑖𝑗)

Computing time Time necessary for completion of a j eme task of i eme graph i instance 𝑓𝑡(𝑡𝑖𝑗) ) Completion time Time at which a j ieme task of i eme graph instance finishes its execution 𝑑𝑙(𝑡𝑖𝑗)

Deadline time Time by which execution of the j ieme task of i eme graph instance should be completed, after the task is released. τij the period of 𝜏 𝑖𝑗 The periodic behavior of task, the study be restricted to the. hyperperiod H which is the least common multiple of the periods. 𝜏 𝑖

The task model real-time system is usually described by a set of task graphs. The related concepts are defined as follows : System: A real-time system is modeled as a set G = {Gi, i = 1, . . g} of periodic task graph.

Task graph: A task graph Gi is a directed acyclic graph (DAG) of period τi and defined by a 2-tuple Gi = {Ti, Ei, τi} where :

-𝑇𝑖 = {𝑡𝑖𝑗, 𝑗 = 1, . . 𝑛𝑖} represents the set of tasks, -𝐸𝑖 = {𝑒𝑖𝑗𝑙, 𝑗 = 1, . . 𝑛𝑖} represent s the set of communications.

Task : A task 𝑡𝑖𝑗 is a node of task graph and defined by 𝑡 𝑖𝑗 𝑘 denote the processor implementation ( instance of task). The problem consists of three parts : Part 1 : Partitioning P ,which determines the assignement of each task to p u , u ∈ {1, … . p} , processors :

∀t ij , P(t ij ) = p u , u ∈ {1, … . p} (1)
Part 2 : Scheduling on processors SC1 u , which associates each task instance assigned on processors {1, … . p} a start time st.

∀P(t ij ) = p u , t ij → SCh 1 u (t ij k ) = st(t ij k ), (2)
u ∈ {1, … . p} Part 3 : Scheduling on interconnexion system 𝑆𝐶ℎ_2 𝐼 , Which associates each interprocessors communication a start time st. The communication defined by 𝑒 𝑖𝑗𝑙 is a directed edge, which means task 𝑡 𝑖𝑗 precedes 𝑡 𝑖𝑙 . Each communication is associated with a data transmission time on a specific communication channel. It is commonly assumed, in multiprocessor systems research, that communication between tasks assign to the same processor is effectively instantaneous, compared with inter processors communication.

e ijl k → SCh_2I(e ijl k ) = st(e ijl k ), {e ijl k ∈ E i k ∕ P(t ij ) ≠ P(t il )} (3)
Scheduling in real-time systems requires that the order of execution be stipulated. A schedule is called feasible if both of the following conditions are met : 1) All precedence constraints are respected and 2) All timing constraints are satisfied such as the deadline defined by 𝑑𝑙(𝑡 𝑖𝑗 𝑘 ).

which means the task must be finished before the deadline, otherwise the system will fail. In the multiprocessor system The task scheduling is a mapping process of all tasks to p processors. The goal of the scheduler is to schedule the system's tasks on these p processors, so that every task is completed before the expiration of the deadline 𝑑𝑙(𝑡 𝑖𝑗 𝑘 ).

Here we use the similar cost function as [START_REF] Abbass | A Single Queen Single Worker Honey-Bees Approach to 3-SAT. The genetic and Evolutionary computation conference[END_REF]. The cost function formulated by equation ( 4) serves to calculate the tardiness scheduling of the real time system.

𝑇𝑎𝑟𝑑𝑖𝑛𝑒𝑠𝑠

The Real-time tasks complete their computing time before their deadlines as early as possible. Differents comportments are showed in Figure 1.

Hybrid HBMO for real time scheduling resolution problem

The HBMO algorithm was proposed by Abbass in [START_REF] Abbass | A Single Queen Single Worker Honey-Bees Approach to 3-SAT. The genetic and Evolutionary computation conference[END_REF]. Since then it has been used on a number of different applications [START_REF] Sahoo R R | Navigational Path Planning of Multi-Robot using Honey Bee Mating Optimization Algorithm (HBMO)[END_REF][19] [START_REF] Sabar | Solving Examination Timetabling Problems using Honey-bee Mating Optimization (ETP-HBMO)[END_REF][21] [START_REF] Abedinia | Robust LFC in deregulated environment: Fuzzy PID using HBMO[END_REF]. HBMO was found to outperform some better known algorithms. However, it has not been applied to realtime scheduling problems. A honey-bee colony consists of the queen(s), drones, worker(s) and broods. The Honey-bee Mating Optimization algorithm mimics the natural mating behavior of the queen bee when she leaves the hive to mate with drones in the air. After each mating, the genetic pool of the queen is enhanced by adding sperm to her spermatheca. Before the mating flight begins, the queen is initialized with a certain amount of energy and only ends her mating flight when the energy level drops below a threshold level (which is close to zero) or when her spermatheca is full. The probability of a drone mating with a queen obeys the following annealing function:

𝑃𝑟𝑜𝑏(𝐷) = 𝑒 [-∆(𝑓∆)/𝑠𝑝𝑒𝑒𝑑(𝑡)] ( 6 )
where (f) is the absolute difference between the fitness of D and the fitness of the queen and Speed(t) is the speed of the queen at a given time t. After each flight, the queen's speed and energy evolve according to the following equations:

𝑆𝑝𝑒𝑒𝑑(𝑡 + 1) =∝× 𝑆𝑝𝑒𝑒𝑑(𝑡) (7) 𝐸𝑛𝑒𝑟𝑔𝑦(𝑡 + 1) = 𝐸𝑛𝑒𝑟𝑔𝑦(𝑡) (8)
where factor   (0,1) is the amount of speed and energy reduction after each flight.

The workers are presented as heuristics whose functionality is to improve the broods produced during the mating process.

To diversify the initial HBMO population and select the best solution as the queen, Genetic Algorithms (GA) and GRASP method are used. The former, GA, works very well on mixed (continuous and discrete) combinatorial problems. Genetic Algorithms create a population of solutions and apply genetic operators such as mutation and crossover to improve the solutions in order to find the best one(s) [24]. The latter, GRASP, is a multi-start or iterative process, in which each GRASP iteration consists of two phases, a construction phase, in which a feasible solution is produced, and a local search phase, in which a local optimum in the neighborhood of the constructed solution is sought. GRASP has proved its efficiency in research and computer science applications as well as in industrial applications.

To improve the broods produced during the queen's mating, we use two workers based on local search metaheuristics: Tabu Search (TS) and Simulated Annealing (SA).

PRESENTATION OF THE RESOLUTION PROCESS

In this section, we present the approach adapted to resolve the issues addressed.

To start with, a set of parameters must be defined:

 The bee population represents the set of scheduling plans of real-time tasks with temporal constraints,  The queen represents the tasks best scheduling plan in the population generated by GA, by GRASP, randomly or hybridization of a percentage from each of the three populations obtained using the three aforementioned procedures and gradually improved by the implementation of the neighborhood generations iterative procedures (crossover in the HBMO) and two workers SA and TS.



The drones make up the remaining task scheduling plans Algorithm 2 describes The HBMO approach for scheduling

SIMULATIONAND EXPERIMENTAL RESULTS

This part is devoted to the implementation and testing of the algorithmic methods developed above. In order to assess the effectiveness of the approaches under study, we constructed a testbench consisting three task systems, which are randomly generated. The systems were composed of 40, 60 and 80 task graphs. Each task graph had its own temporal characteristics. The parameters of the algorithm have been selected after thorough testing. A number of different values were tested and the ones selected are those that yielded the best results in both solution quality and computational time. The best values for these parameters appear in the table 2 and table 3. In the first time, we generate the initial population of the proposed algorithm randomly and with GA. After 10 runs , The figure 2 illustrate the graphic results obtained through the application of HBMO, GA and GRASP approaches. The best real-time task scheduling approach is to be selected taking into consideration the Cost function criteria. The graphics results show that the HBMO approach converge to the optimal results comparing with the GA and GRASP methods.

Fig. 2. Graphic results obtained through different hybridization processes example of 60 tasks

In order to ameliorate scheduling real-time problem results obtained through the application of basic HBMO (HBMO_ Random) , we have hybridized the HBMO algorithm with other metaheuristics: Genetic algorithms GA, GRASP and hybridization (three sub-populations made up of one percentage from random, one percentage from AG, one percentage from GRASP) in the stage of initial population generation. The figure 3 presents all graphic results obtained through different hybridization processes. The graphics results show that the results given by HBMOhybridization approach were best than the results given by the HBMO with GA and HBMO with GRASP algorithms. We notice also that the results given by the HBMO-GRASP and HBMO-hybridization were nearly the same. (in the iterations number 34, 24, 16 and 9) the values obtained by HBMO_ hybridization and the HBMO_GRASP are the same ).

CONCLUSION

In this paper, an approach is proposed to solve real-time task scheduling problems. This approach is based on the hybridization of several algorithms inspired from nature such as Genetic, HBMO, Tabu Search and Simulated Annealing and GRASP algorithms. In the stage of generating initial solutions, the adopted approach uses Genetic Algorithms alone, random selection alone, GRASP alone, and hybridization of the three processes selecting a percentage of the population from each. In the stage of solution improvement by the workers, HBMO uses Tabu Search and Simulated Annealing algorithms. The proposed approach was tested on three real-time task scheduling instances with 40 tasks, 60 tasks and 80 tasks where each task has its own temporal constraints. The gain coefficient from the optimal solution is used as the standard to evaluate the quality of the solutions obtained. Results confirm the positive impact of using a hybrid strategy with regard to function objective quality and computing time in comparison with the hybridization algorithms object of this paper. The results arrived at show that for real-time scheduling problems (with a number of tasks 80), the best approach is HBMO-Hybridization. In the future works, we will focus our searches to resolve the problem of minimizing the energy consumption in the realtime embedded systems by using the proposed algorithms.

Fig. 1 .

 1 Fig. 1. Comportments of the real time instances of the task

Algorithm 2 .. Initialization 2 .

 22 Hybrid HBMO scheduling approach 1Generate the initial population of the bees using separately GA, GRASP, Random or through Hybridization of percentages selected from these procedure; killed All broods will be used in the next mating flighFitness function Fitness functionObjective function 𝑙𝑎𝑡𝑒𝑛𝑒𝑠𝑠(𝐺) of the real time system

Table 2 .Table 3 .

 23 Parameters and their values for the HBMO approach adopted Parameters and their values for the used approaches

Fig. 3 .

 3 Fig. 3. Graphic results obtained through different hybridization processes example of 80 tasks).

A summary of differences and similarities between the HBMO approach applied for the resolution of the real time scheduling problem and the original algorithm are showed in table 1


Farrukh Jamal

Christophe Chesneau

Mohammed Elgarhy

Type II General Inverse Exponential family of distributions

Keywords: T-X transformation, Inverse exponential distribution, Maximum likelihood estimation. 2000 MSC: 60E05, 62E15, 62F10

come

Introduction

In recent years, numerous methods have been introduced to increase the flexible properties of classic probability distributions. One of the most famous of them is the so called T-X transformation introduced by [START_REF] Alzaatreh | A new method for generating families of continuous distributions[END_REF]. The related T-X family of distributions is characterized by a cumulative distribution function (cdf) described below. Let (a, b) ∈ R 2 with a < b, p(t) is a probability density function (pdf) with support [a, b], G(x) a cdf and W (x) a function such that W [G(x)] satisfies the three following conditions :

i) W [G(x)] ∈ [a, b], ii) W [G(x)
] is differentiable and monotonically non-decreasing, iii) W [G(x)] → a when x → -∞ and W [G(x)] → b when x → +∞. Then the cdf of the T-X family is given by

F (x) = W [G(x)] a p(t)dt,
x ∈ R.

(1)

Among the numerous members of the T-X family, all involving new configurations for p(t), W (x) and g(x), there are the Exp-G (Kw-G type 2) family by Cordeiro et al. (2013), the Weibull-X family by [START_REF] Alzaatreh | A new method for generating families of continuous distributions[END_REF], the Gamma-X family by [START_REF] Alzaatreh | The gamma-normal distribution: Properties and applications[END_REF], the Exponentiated T-X family by [START_REF] Alzaghal | Exponentiated T-X family of distributions with some applications[END_REF], the Weibull-G family by [START_REF] Bourguignon | The Weibull-G family of probability distributions[END_REF], the Logistic-G family by Torabi and Montazari (2014), The T-Burr family by [START_REF] Nasir | A New Generalized Burr Family of Distributions Based on Quantile Function[END_REF], the Ofr family by [START_REF] Haq | The Odd Frêchet-G Family of Probability Distributions[END_REF] and the TIIGE family by [START_REF] Hamedani | Type II General Exponential Class of Distributions[END_REF]. All these families demonstrate nice theoretical and practical properties, offering new solutions in terms of statistical models for the practitioners.

Let us now present and motivate the considered distribution in this paper. Adopting the notations related the cdf of the T-X family given by (1), we consider for p(t) the pdf of the inverse exponential distribution defined by

p(t) = λ t 2 e -λ t , t ∈ (0, +∞),
(observe that it is also a special case of the Frêchet distribution), implying that a = 0 and b = +∞, and the function W (x) given by

W (x) = 1 (1 -x) θ -1.
By taking any cdf G(x), one can observe that W [G(x)] satisfies the three required conditions, i.e. i), ii) and iii). Note that W [G(x)] can be written as the odds function of a cdf G * (x) belonging to the Lehmann type II distribution with cdf G(x) and parameter θ, i.e. W [G(x)] = G * (x)

1-G * (x) , with G * (x) = 1-[1-G(x)] θ . Further details on the Lehmann type II distribution can be found in [START_REF] Gupta | Modeling failure time data by Lehman alternatives[END_REF]. Another point of view is to view W (x) as the quantile function of the Lomax distribution with parameters α = 1 θ and β = 1, i.e. with cdf H(x) = 1 -(1 + x) -1 θ . With the presented functions, the cdf F (x) given by (1) becomes

F (x) = W [G(x)] 0 λ t 2 e -λ t dt = e -λ t W [G(x)] 0 = e -λ W [G(x)] = e -λ [1-G(x)] θ 1-[1-G(x)] θ . (2)
For the purpose of this paper, we call the related family of distributions the TIIGIE family (for Type II Generalized Inverse Exponential). The pdf of the general Frêchet distribution for p(t) has ever been considered in [START_REF] Haq | The Odd Frêchet-G Family of Probability Distributions[END_REF] but with a completely different function W (x). On the other side, the considered function W (x) has ever been considered (as a quantile function) in [START_REF] Nasir | A New Generalized Burr Family of Distributions Based on Quantile Function[END_REF] but with a completely different pdf p(t). To the best of our knowledge, the TIIGIE family is new in the literature. This paper deals with a complete mathematical and practical studies of this family, with discussions. We first present the expressions of the crucial functions related to the family, then obtain useful expansions of the pdf, cdf and hrf, a general expression for the quantile function, probability weighted, ordinary and incomplete moments, generating function, entropies and order statistics, with moments. Practical investigations are done for special cases, showing that the TIIGIE family is a quite flexible family of distributions to fit real data from several fields. To be more specific, from Figures 1, 2, and 3 presented in a next section, we can see that the possible pdf shapes of TIIGIE family are J shape, reverse J shape, right skewed, left skewed and symmetrical. This means that the TIIGIE family can show suitable fit to those data sets, whose histograms are similar to the TIIGIE family pdfs shapes. Further, the TIIGIE family exhibits monotone [increasing (IFR) and decreasing (DFR)], non-monotone [bathtub (BT) and upside-down bathtub (UBT)] and decreasing-increasing-decreasing (DID) hrf shapes to cope with all types of lifetime data sets. This potentiality is illustrated with applications to two real life data sets, by considering the Bur III distribution as baseline and the maximum likelihood method for the estimations of the parameters. The rest of the paper is organized as follows. Section 2 studies some immediate properties of the functions characterizing the TIIGIE family. Some special cases with plots are presented in Section 3. A comprehensive account of some of its mathematical properties is proposed in Section 4. Section 5.1 provides the necessary to the estimation of the unknown parameters with the maximum likelihood method and the least square method, with a short simulation study. Applications to two real life data sets are given in Section 6.

On the crucial functions of the TIIGIE family

Some expressions

We now give the general expressions of the reliability function (rf), the pdf and the hrf of the TIIGIE family. The rf is given by

R(x) = 1 -F (x) = 1 -e -λ [1-G(x)] θ 1-[1-G(x)] θ , x ∈ R.
Let us now denote by g(x) the pdf related to G(x). By the derivation of F (x), the pdf is given by

f (x) = λθg(x)[1 -G(x)] θ-1 (1 -[1 -G(x)] θ ) 2 e -λ [1-G(x)] θ 1-[1-G(x)] θ , x ∈ R. ( 3
)
The hrf is given by

h(x) = f (x) R(x) = λθg(x)[1 -G(x)] θ-1 (1 -[1 -G(x)] θ ) 2 1 -e -λ [1-G(x)] θ 1-[1-G(x)] θ e -λ [1-G(x)] θ 1-[1-G(x)] θ , x ∈ R.

Shapes of the pdf and hrf

The shapes of the pdf and the hrf of the TIIGIE family can be described analytically. As usual, the critical points of the pdf f (x) are the roots of the equation given by

∂ ∂x ln(f (x)) = 0 with ∂ ∂x ln(f (x)) = g (x) g(x) -(θ -1) g(x) 1 -G(x) -2θ [1 -G(x)] θ-1 g(x) 1 -[1 -G(x)] θ + λθ [1 -G(x)] θ-1 g(x) (1 -[1 -G(x)] θ ) 2 .
This equation can have more than one root. If x * is a root of this equation, then it corresponds to a local maximum if

∂ 2 ∂x 2 ln(f (x * )) < 0, a local minimum if ∂ 2 ∂x 2 ln(f (x * )) > 0 and a point of inflection if ∂ 2 ∂x 2 ln(f (x * )) = 0.
In a similar way, the critical points of the hrf h(x) are the roots of the equation given by ∂ ∂x ln(h(x)) = 0 with

∂ ∂x ln(h(x)) = g (x) g(x) -(θ -1) g(x) 1 -G(x) -2θ [1 -G(x)] θ-1 g(x) 1 -[1 -G(x)] θ + λθ [1 -G(x)] θ-1 g(x) (1 -[1 -G(x)] θ ) 2 + λθg(x)[1 -G(x)] θ-1 (1 -[1 -G(x)] θ ) 2 1 -e -λ [1-G(x)] θ 1-[1-G(x)] θ e -λ [1-G(x)] θ 1-[1-G(x)] θ .
Again, this equation can have more than one root. If x * is a root of this equation, then it corresponds to a local maximum if

∂ 2 ∂x 2 ln(h(x * )) < 0, a local minimum if ∂ 2 ∂x 2 ln(h(x * )) > 0 and a point of inflection if ∂ 2 ∂x 2 ln(h(x * )) = 0.

Asymptotic results

In order to show the effect of the parameters on tails distributions, we now study the asymptotic properties of F (x), f (x) and h(x). First of all, let us remark that W (y) ∼ θy when y → 0 and W (y) ∼ 1

(1-y) θ when y → 1. Therefore

F (x) ∼ e -λ θG(x) , G(x) → 0, F (x) ∼ e -λ[1-G(x)] θ , x → +∞.
Similarly, we have

f (x) ∼ λθg(x) (1 -[1 -G(x)] θ ) 2 e -λ θG(x) , G(x) → 0 and f (x) ∼ λθg(x)[1 -G(x)] θ-1 , x → +∞.
Concerning the hrf, the following asymptotic results hold :

h(x) ∼ λθg(x) (1 -[1 -G(x)] θ ) 2 e -λ θG(x) , G(x) → 0 and h(x) ∼ θ g(x) 1 -G(x) , x → +∞.
We observe that, when x tends to +∞, the hrf h(x) is proportional to the hrf related to G(x) ; the constant of proportionality is given by θ.

Linear representation

Let h > 0. We now provide useful linear representations for several functions, including [F (x)] h f (x). It follows from the formula of exponential power series that

F (x) = +∞ k=0 1 k! (-1) k λ k [1 -G(x)] θk (1 -[1 -G(x)] θ ) -k .
The generalized binomial series gives

(1 -[1 -G(x)] θ ) -k = +∞ =0 -k (-1) [1 -G(x)] θ and [1 -G(x)] θ(k+ ) = +∞ m=0 θ(k + ) m (-1) m [G(x)] m .
Combining these equalities, we obtain the following series expansion:

F (x) = +∞ m=0 a m Π m (x),
where

a m = (-1) m +∞ k=0 +∞ =0 (-1) k+ 1 k! λ k -k θ(k+ ) m and Π m (x) = [G(x)
] m is the cdf of the well-known Exp-G distribution with power m and Π 0 (x) = 1 (see [START_REF] Nadarajah | The Exponentiated-type Distributions[END_REF]). By derivation of F (x), f (x) can be expressed as

f (x) = +∞ m=0 a m+1 π m+1 (x), (4)
where

a m+1 = (-1) m+1 +∞ k=0 +∞ =0 (-1) k+ 1 k! λ k -k θ(k+ ) m+1 and π m+1 (x) = (m+1) [G(x)] m g(x) is the pdf of the Exp-G distribution with power m + 1. Now let us observe that [F (x)] h f (x) = 1 h+1 f * (x)
, where f * (x) denotes a pdf of the TIIGIE with parameter (h + 1)λ instead of λ. Then, by analogy with (4), we can write

[F (x)] h f (x) = +∞ m=0 b m,h π m+1 (x), ( 5
) with b m,h = (-1) m+1 +∞ k=0 +∞ =0 (-1) k+ 1 k! λ k (h + 1) k-1 -k θ(k+ )
m+1 . The series expansion (5) can be used to derive most of the mathematical properties of the TIIGIE family. This is developed in Section 4.

Some special cases

Among the numerous possible distributions arising from the TIIGIE family, we now present three special cases using classics distributions as baselines.

TIIGIE-Uniform distribution

We define the TIIGIE-U distribution with parameters (α, λ, θ) the TIIGIE family with the uniform distribution on the interval [0, α] as baseline distribution, i.e. with cdf

0 1 2 3 4 0.0 0.5 1.0 1.5 α = 4 x pdf λ = 3 θ = 0.1 λ = 3 θ = 0.5 λ = 3 θ = 2 λ = 3 θ = 6 λ = 1.5 θ = 6 λ = 0.5 θ = 6 0 1 2 3 4 0.0 0.5 1.0 1.5 2.0 2.5 3.0 α = 4 x hrf λ = 3 θ = 0.1 λ = 3 θ = 0.4 λ = 1 θ = 0.1 λ = 2 θ = 2 λ = 0.4 θ = 0.1
G(x) = x/α, x ∈ [0, α], G(x) = 1, x > α and g(x) = 1/α, x ∈ [0, α]. The pdf of the TIIGIE-U distribution is given by f (x) = λθ[1 -x/α] θ-1 α(1 -[1 -x/α] θ ) 2 e -λ [1-x/α] θ 1-[1-x/α] θ , x ∈ [0, α].
The associated hrf is given by

h(x) = λθ[1 -x/α] θ-1 α(1 -[1 -x/α] θ ) 2 1 -e -λ [1-x/α] θ 1-[1-x/α] θ e -λ [1-x/α] θ 1-[1-x/α] θ , x ∈ [0, α].
Plots of the pdf and hrf of the TIIGIE-U distribution for some parameter values are displayed in Figure 1.

TIIGIE-Weibull distribution

We define the TIIGIE-W distribution with parameters (λ, θ, α, β) the TIIGIE family with the Weibull distribution with parameters (α, β) as baseline distribution, i.e. with cdf G(x) = 1 -e -αx β and g(x) = αβx β-1 e -αxβ , x > 0. The pdf of the TIIGIE-W distribution is given by

f (x) = λθαβx β-1 e -αxβ e -α(θ-1)x β (1 -e -αθx β ) 2 e -λ e -αθx β 1-e -αθx β , x > 0.
The associated hrf is given by

h(x) = λθαβx β-1 e -αxβ e -α(θ-1)x β (1 -e -αθx β ) 2 1 -e -λ e -αθx β 1-e -αθx β e -λ e -αθx β 1-e -αθx β , x > 0.
Plots of the pdf and hrf of the TIIGIE-W distribution for some parameter values are displayed in Figure 2.

= 2 α = 3 β = 0.5 λ = 2 θ = 5 α = 4 β = 0.6 λ = 3 θ = 2 α = 2 β = 1 λ = 2 θ = 2 α = 3 β = 2

TIIGIE-Burr III distribution

We define the TIIGIE-BIII distribution with parameters (λ, θ, c, k) the TIIGIE family with the Burr III distribution, also known as exponentiated Lomax distribution, with parameters (c, k) as baseline distribution, i.e. with cdf

G(x) = (1 + x -c ) -k and g(x) = ckx -c-1 (1 + x -c ) -k-1 , x > 0. The pdf of the TIIGIE-BIII distribution is given by f (x) = λθckx -c-1 (1 + x -c ) -k-1 [1 -(1 + x -c ) -k ] θ-1 (1 -[1 -(1 + x -c ) -k ] θ ) 2 e -λ [1-(1+x -c ) -k ] θ 1-[1-(1+x -c ) -k ] θ , x > 0.
The associated hrf is given by

h(x) = λθckx -c-1 (1 + x -c ) -k-1 [1 -(1 + x -c ) -k ] θ-1 (1 -[1 -(1 + x -c ) -k ] θ ) 2 1 -e -λ [1-(1+x -c ) -k ] θ 1-[1-(1+x -c ) -k ] θ e -λ [1-(1+x -c ) -k ] θ 1-[1-(1+x -c ) -k ] θ , x > 0.
Plots of the pdf and hrf of the TIIGIE-BIII distribution for some parameter values are displayed in Figure 3. We observe shapes of different natures for the functions, which motivates us to used it as statistical model in a next section.

Mathematical properties

We now investigate some mathematical properties of the TIIGIE family. We denote by X a random variable having the cdf F (x) (and the pdf f (x)) and by Y a random variable having the cdf G(x) (and the pdf g(x)).

Quantile function

The quantile function of X is given by In particular, the median is given by

x p = G -1 1 - ln(p) ln(p) -λ 1 θ , p ∈ (0, 1]. 0.0 0.5 1.0 1.5 2.0 2.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 x pdf λ = 3 θ = 0.5 c = 4 k = 0.5 λ = 3 θ = 2 c = 2 k = 1 λ = 3 θ = 3 c = 2 k = 2 λ = 5 θ = 3 c = 0.5 k = 0.7 λ = 1 θ = 3 c = 5 k = 8 0.0 0.5 1.0 1.5 2.0 2.5 0 2 4 6 8 x hrf λ = 3 θ = 0.1 c = 4 k = 0.5 λ = 3 θ = 0.5 c = 2 k = 0.1 λ = 0.5 θ = 0.1 c = 6 k = 0.5 λ = 3 θ = 3 c = 2 k = 2 λ = 3 θ = 0.1 c = 4 k = 0.2 λ = 3.2 θ = 0.5 c = 7.3 k = 0.1
x 1/2 = G -1 1 -ln(2)
ln( 2)+λ 1 θ . We can also use

x p to define well-known quantile measures, as the Bowley skewness based on quartiles and the Moors kurtosis based on octiles, respectively defined by

B = x 3/4 + x 1/4 -2x 1/2 x 3/4 -x 1/4 , M = x 7/8 -x 5/8 + x 3/8 -x 1/8 x 6/8 -x 2/8 .
Further details can be found in [START_REF] Kenney | Mathematics of Statistics[END_REF] and [START_REF] Moors | A quantile alternative for kurtosis[END_REF].

Finally, if U is a random variable having the uniform distribution on the interval (0, 1), then the random variable X = x U has the cdf F (x).

Probability weighted moments

Let (r, h) ∈ (0, +∞) 2 . The (r, h)-th probability weighted moments of X and Y are respectively defined by

τ r,h = E X r [F (X)] h = +∞ -∞ x r [F (x)] h f (x)dx, τ * r,h = E Y r [G(Y )] h .
Using the series expansion given by ( 5), we obtain

τ r,h = +∞ m=0 b m,h +∞ -∞ x r π m+1 (x)dx = +∞ m=0 c m,h τ * r,m , (6)
where c m,h = (m + 1)b m,h . The term τ * r,m can be computed numerically by using the quantile function

Q G (u) = G -1 (u) as τ * r,m = 1 0 [Q G (u)] r u m du.

Moments

The r-th ordinary moment of X is defined by µ r = E (X r ) = +∞ -∞ x r f (x)dx. Using (4), we can express µ r as

µ r = +∞ m=0 a m+1 +∞ -∞ x r π m+1 (x)dx.
The integral terms can be computed numerically by noticing that

+∞ -∞ x r π m+1 (x)dx = (m + 1) 1 0 [Q G (u)] r u m du. Alternatively, using (6), one can notice that µ r = τ r,0 = +∞ m=0 c m,0 τ * r,m .
The mean of X is obtained by taking r = 1. Let r be an integer. The r-th central moments of X can be expressed as

µ r = E[(X -µ 1 ) r ] = r k=0 r k (-1) k (µ 1 ) k µ r-k
and the r-th cumulants of X can be obtained via the equation :

κ r = µ r - r-1 k=1 r-1 k-1 κ k µ r-k with κ 1 = µ 1 . The skewness of X is given by γ 1 = κ 3 /κ 3/2 2
and the kurtosis of X is given by γ 2 = κ 4 /κ 2 2 .

Moment generating function

The moment generating function of X is given by M X (t) = E e tX = +∞ -∞ e tx f (x)dx. Using (4), we can express M X (t) as

M X (t) = +∞ m=0 a m+1 +∞ -∞ e tx π m+1 (x)dx.
The integral terms can be computed numerically by noticing that +∞ -∞ e tx π m+1 (x)dx = (m + 1) 1 0 e tQ G (u) u m du. An alternative expression using (6) is given by

M X (t) = E e tX = +∞ r=0 t r r! µ r = +∞ r=0 +∞ m=0 t r r! c m,0 τ * r,m .

Incomplete moments

The r-th incomplete moment of X is defined by m r (t

) = E X r 1 {X≤t} = t -∞ x r f (x)dx. Using (4), we can express m r (t) as m r (t) = +∞ m=0 a m+1 t -∞ x r π m+1 (x)dx.
Again, the integral terms can be expressed as

t -∞ x r π m+1 (x)dx = (m+1) G(t) 0 [Q G (u)] r u m du.
The incomplete moments of X are useful tools in the definitions of several important quantities. Among them, let us mention the Lorenz curve, the Bonferroni curve, the mean deviation about the mean, the mean deviation about the median, the residual life function and the reversed residual life function.

Entropies

The Rényi entropy of X is defined by

I δ (X) = 1 1 -δ log +∞ -∞ [f (x)] δ dx , with δ > 0 and δ = 1. It follows from (3) that [f (x)] δ = λ δ θ δ [g(x)] δ [1 -G(x)] (θ-1)δ (1 -[1 -G(x)] θ ) 2δ e -λδ [1-G(x)] θ 1-[1-G(x)] θ .
Using the exponential series expansion, we can write k+2δ) .

e -λδ [1-G(x)] θ 1-[1-G(x)] θ = +∞ k=0 1 k! (-1) k λ k δ k [1 -G(x)] θk (1 -[1 -G(x)] θ ) -k . Therefore [f (x)] δ = λ δ θ δ [g(x)] δ +∞ k=0 1 k! (-1) k λ k δ k [1 -G(x)] θ(k+δ)-δ (1 -[1 -G(x)] θ ) -(
The generalized binomial series gives

(1 -[1 -G(x)] θ ) -(k+2δ) = +∞ =0 -(k + 2δ) (-1) [1 -G(x)] θ and [1 -G(x)] θ(k+ +δ)-δ = +∞ m=0 θ(k + + δ) -δ m (-1) m [G(x)] m .
Combining these equalities, we obtain the following series expansion:

[f (x)] δ = λ δ θ δ +∞ m=0 d m,δ [G(x)] m [g(x)] δ , where d m,δ = (-1) m +∞ k=0 +∞ =0 (-1) k+ 1 k! λ k δ k -(k+2δ) θ(k+ +δ)-δ m
. Therefore the Rényi entropy of TIIGIE family is given by

I δ (X) = 1 1 -δ δ log(λ) + δ log(θ) + log +∞ m=0 d m,δ +∞ -∞ [G(x)] m [g(x)] δ dx .
The integrals terms can be evaluated numerically for a given G(x).

The δ-entropy is defined by

H δ (X) = 1 δ -1 log 1 - +∞ -∞ [f (x)] δ dx ,
with δ > 0 and δ = 1. So we have

H δ (X) = 1 δ -1 log 1 -λ δ θ δ +∞ m=0 d m,δ +∞ -∞ [G(x)] m [g(x)] δ dx .
Finally, the Shannon entropy of X is defined by S(X) = E(-log[f (X)]). It is in fact a particular case of the Rényi entropy when δ tends to 1 + .

Order statistics

Order statistics are fundamental in many areas of statistical theory and practice. Let X 1 , . . . , X n be a random sample from the TIIGIE family. The pdf of the ith order statistic, say X i:n , is given by

f i:n (x) = n! (i -1)!(n -i)! n-i j=0 (-1) j n -i j [F (x)] j+i-1 f (x).
Using the series expansion (5), we have

f i:n (x) = n! (i -1)!(n -i)! n-i j=0 (-1) j n -i j +∞ m=0 b m,j+i-1 π m+1 (x).
The r-th ordinary moment of X i:n is defined by µ i:n,r = E (X r i:n ). It follows from the above series expansion, the definition of the probability weighted moments of Y and (6) that

µ i:n,r = n! (i -1)!(n -i)! n-i j=0 (-1) j n -i j τ r,j+i-1 = n! (i -1)!(n -i)! n-i j=0 (-1) j n -i j +∞ m=0 c m,j+i-1 τ * r,m ,
(with c m,j+i-1 = (m + 1)b m,j+i-1 ).

Estimation inference

Maximum likelihood estimation (MLE)

Let x 1 , . . . , x n be n observed values from the TIIGIE family and ξ = (λ, θ, ξ * ) be the vector of unknown parameters, ξ * denoting the vector of parameters related to the distribution characterized by the cdf G(x). The log likelihood function is given by

(ξ) = n log(λ) + n log(θ) + n i=1 log(g(x i )) + (θ -1) n i=1 log(1 -G(x i )) -2 n i=1 log(1 -[1 -G(x i )] θ ) -λ n i=1 [1 -G(x i )] θ 1 -[1 -G(x i )] θ .
The maximum likelihood estimators of the parameters are obtained by maximizing the log likelihood function. They can be obtained by solving the non-linear equations :

∂ ∂λ (ξ) = 0, ∂ ∂θ (ξ) = 0 and ∂ ∂ξ * (ξ) = 0 with ∂ ∂λ (ξ) = n λ - n i=1 [1 -G(x i )] θ 1 -[1 -G(x i )] θ . ∂ ∂θ (ξ) = n θ + n i=1 log(1 -G(x i )) + 2 n i=1 [1 -G(x i )] θ log(1 -G(x i )) 1 -[1 -G(x i )] θ -λ n i=1 [1 -G(x i )] θ log(1 -G(x i )) (1 -[1 -G(x i )] θ ) 2 . Let us set g ξ * (x i ) = ∂ ∂ξ * g(x i ) and G ξ * (x i ) = ∂ ∂ξ * G(x i ). Then we have ∂ ∂ξ * (ξ) = n i=1 g ξ * (x i ) g(x i ) -(θ -1) n i=1 G ξ * (x i ) 1 -G(x i ) -2θ n i=1 [1 -G(x i )] θ-1 G ξ * (x i ) 1 -[1 -G(x i )] θ + λθ n i=1 [1 -G(x i )] θ-1 G ξ * (x i ) (1 -[1 -G(x i )] θ ) 2 .
They can be solved using the Newton method or fixed point iteration methods. As usual, for determining related mathematical quantities as the variance co-variance matrix and the confidence interval for parameters, we need the information matrix which can be generated by the taking the expectation of the second order derivative.

Least square method (LSE)

Adopting the notations above, let x (1) , . . . , x (n) be x 1 , . . . , x n in increasing order. Least square estimates are obtained by minimizing the following function

S(ξ) = n i=1 F (x (i) ) - i n + 1 2 = n i=1 e -λ [1-G(x (i) )] θ 1-[1-G(x (i) )] θ - i n + 1 2 .
Minimizing S(Θ) with respect to λ, θ and ξ, we have following system of non linear equations:

∂S(ξ) ∂λ = -2 n i=1 e -λ [1-G(x (i) )] θ 1-[1-G(x (i) )] θ - i n + 1 [1 -G(x (i) )] θ 1 -[1 -G(x (i) )] θ e -λ [1-G(x (i) )] θ 1-[1-G(x (i) )] θ = 0, ∂S(ξ) ∂θ = -2λ n i=1 e -λ [1-G(x (i) )] θ 1-[1-G(x (i) )] θ - i n + 1 [1 -G(x (i) )] θ log(1 -G(x (i) )) 1 -[1 -G(x (i) )] θ 2 e -λ [1-G(x (i) )] θ 1-[1-G(x (i) )] θ = 0, ∂S(ξ) ∂ξ * = 2λθ n i=1 e -λ [1-G(x (i) )] θ 1-[1-G(x (i) )] θ - i n + 1 [1 -G(x (i) )] θ-1 G ξ * (x (i) ) (1 -[1 -G(x (i) )] θ ) 2 e -λ [1-G(x (i) )] θ 1-[1-G(x (i) )] θ = 0.
This system of non-linear equations can be solved numerically by any software to obtained the estimates.

A simulation study

In order to assess the performance of the maximum likelihood and least square methods for estimating the parameters, a small simulation is carried out. For such purposes, a TIIGIE-BIII distribution is considered using Monte Carlo simulations. The process is carried out as follow:

• The number of Monte Carlo replications was made 1000 times each with sample sizes n = 30, 50, 100 and 300.

• Initial values for the parameters are selected as given in Tables 1, 2, 3 and4, • The formula for the mean squared error (MSE) of the estimate λ of λ is given by 1 1000

1000 i=1 ( λ -λ) 2 .
• The third step is also repeated for the other parameters.

These numerical results show that the considered estimation methods perform quite well in estimating the model parameters of the TIIGIE-BIII distribution.

Table 1: Estimates and MSEs of TIIGIE-BIII distribution for ML and LS estimate, Set 1 : (λ, θ, c, k) = (0.5, 0.5, 0.5, 0.5) and Set 2 : (λ, θ, c, k) = (0.5, 0.5, 0.5, 1.5). Tables 6 and8 provide the values of goodness-of-fit measures for the TIIGIE-BIII model and other fitted models, whereas the MLEs and their corresponding standard errors (SEs) (in parentheses) are listed in Tables 7 and9, respectively. The plots of the fitted model is shown in Figures 4 and5. These plots and Tables 6 and8 indicate that the TIIGIE-BIII model yields the best fit to the others competitive models.

P-P plot

Theoretical probabilities Empirical probabilities

Appendix

The pdfs of the statistical models used in Section 6 are presented below.

• The pdf of the ELx distribution introduced by Abdul-Moniem and Abdel-Hameed ( 2012) is given by

f (x) = aα β 1 + x β -(α+1) 1 -1 + x β -α a-1 , x > 0.
• Tahir, M.H., Cordeiro, G.M., Mansoor, M and Zubair, M. (2015).The Weibull-Lomax distribution: properties and applications, Hacettepe Journal of Mathematics and Statistics, 44, 2, 461-480. Torabi, H. and Montazari, N.H. (2014). The logistic-uniform distribution and its application. Communications in Statistics-Simulation andComputation, 43, 2551-2569.

Figure 1 :

 1 Figure 1: Plots for pdfs and hrfs of the TIIGIE-U distribution.

  θ = 0.5 α = 3 β = 2 λ = 3 θ = 6 α = 3 β = 2.5 λ = 2 θ = 5 α = 4 β = 5 λ = 3 θ = 0.5 α = 0.3 β = 0

Figure 2 :

 2 Figure 2: Plots for pdfs and hrfs of the TIIGIE-W distribution.

Figure 3 :

 3 Figure 3: Plots for pdfs and hrfs of the TIIGIE-BIII distribution.

Figure 4 :

 4 Figure 4: PP, QQ, epdf and ecdf plots of the TIIGIE-BIII distribution for data set 1.

Figure 5 :

 5 Figure 5: PP, QQ, epdf and ecdf plots of the TIIGIE-BIII distribution for data set 2.

••-•

  The pdf of the BLx distribution introduced by Lemonte and Cordeiro (a-1 (1 -x) b-1 dx. • The pdf of the KwLx distribution introduced by Lemonte and Cordeiro (The pdf of the TIIGE-Lx distribution introduced by Hamedani et al. The pdf of the OFr-Lx distribution introduced by Haq and Elgarhy (2018) is given byf (x) = αθ[1 + (x/β)] -(αθ+1) β [1 -[1 + (x/β)] -α ] θ+1 eThe pdf of the WLx distribution introduced by Tahir et aL (2015) is given by f (

Table 2 :

 2 Estimates and MSEs of TIIGIE-BIII distribution for ML and LS estimate, Set 3 : (λ, θ, c, k) = (0.5, 0.5, 1.5, 0.5) and Set 4 : (λ, θ, c, k) = (0.5, 0.5, 1.5, 1.5).

		Set 1				Set 2
		MLE		LSE		MLE		LSE
	n	Estimates MSEs Estimates	MSEs	Estimates MSEs Estimates	MSEs
	30	0.5215	0.0100	0.499743	0.000236	0.5124	0.0092	0.499918 0.000233
		0.5242	0.0224	0.500959	0.000412	0.5350	0.0355	0.500669 0.000411
		0.5415	0.0277	0.497127	0.000047	0.5275	0.0161	0.498370 0.000189
		0.5411	0.0224	0.498763	0.000255	1.5844	0.1707	1.496670 0.002238
	50	0.5148	0.0061	0.500492	0.000138	0.5064	0.0052	0.499829 0.000139
		0.5127	0.0104	0.499697	0.000241	0.5228	0.0119	0.500617 0.000246
		0.5262	0.0140	0.498323	0.000025	0.5181	0.0084	0.499235 0.000114
		0.5274	0.0130	0.499869	0.000147	1.5530	0.1020	1.497740 0.001352
	100	0.5055	0.0026	0.500618	0.000067	0.5080	0.0027	0.499670 0.000070
		0.5074	0.0045	0.499349	0.000114	0.5044	0.0043	0.500584 0.000120
		0.5122	0.0063	0.499296	0.000011	0.5050	0.0033	0.499720 0.000056
		0.5111	0.0055	0.500341	0.000071	1.5414	0.0500	1.497980 0.000680
	300	0.5024	0.0008	0.499964	0.000023	0.5001	0.0009	0.499920 0.000023
		0.5009	0.0012	0.500108	0.000039	0.5044	0.0015	0.500166 0.000039
		0.5015	0.0018	0.499727 3.473×10 -6	0.5016	0.0012	0.499932 0.000019
		0.5028	0.0016	0.499873	0.000024	1.5030	0.0163	1.499490 0.000218

Table 7 :

 7 MLEs and SEs (in parentheses) for data set 1.

	Model		Estimates
	TIIGIE-BIII 1.3035	16.0723	0.5361	11.3254
	(λ, θ, c, k)	(0.3901) (2.8563) (0.2819) (0.0610)
	TIIGE-Lx	0.1892	0.0083	1.4577	10.3473
	(λ, α, a, b)	(0.0155) (0.0032) (0.0128) (0.0225)
	OFr-Lx	1.6698	66.9329	5.2782
	(θ, α, β)	(0.2798) (9.3477) (0.1088)
	WLx	0.1782	1.9231	17.3346	3.7849
	(α, β, a, b)	(0.0916) (2.7434) (3.2247) (1.4321)
	KwLx	83.5491 129.4375 51.6285	0.2297
	(α, β, a, b)	(2.8877) (3.3168) (16.1895) (0.0347)
	BLx	9.6404	5.6777	29.3087	0.6645
	(α, β, a, b)	(0.9320) (2.6266) (0.3286) (0.4529)
	ELx	3.2130	4.2546	41.9873
	(α, β, a)	(1.4229) (6.5601) (0.6488)

Table 8 :

 8 Goodness-of-fit measures for data set 2.

	Model	-	AIC	BIC	CVM	AD	KS	PV
	TIIGIE-BIII -12.5338 -16.8677 -11.9389 0.0333 0.2513 0.0911 0.9642
	TIIGE-Lx	-10.5754 -13.1508 -7.8219 0.0516 0.3876 0.1310 0.7226
	OFr-Lx	-10.6272 -15.2544 -11.2577 0.0650 0.4817 0.1189 0.8231
	WLx	-11.8302 -15.6604 -10.3316 0.0597 0.4468 0.1281, 0.7478
	KwLx	-12.0054 -16.0809 -11.3521 0.0408 0.3142 0.1241 0.7816
	BLx	-12.4135 -16.8270 -11.4982 0.0398 0.2633 0.0999 0.9426
	ELx	-10.9943 -15.9886 -11.9920 0.0575 0.4291 0.1144 0.8567

Table 9 :

 9 MLEs and SEs (in parentheses) for data set 2.

	Model		Estimates
	TIIGIE-BIII	2.2471	15.0705	1.6547	41.7735
	(λ, θ, c, k)	(0.0852) (0.5356) (0.8421) (6.7835)
	TIIGE-Lx	40.0767	40.0138	25.0314	15.1712
	(λ, α, a, b)	(0.4884) (0.4900) (5.1094) (0.0391)
	OFr-Lx	14.2219	27.2773	7.2345
	(θ, α, β)	(3.2679) (0.8324) (2.02942)
	WLx	0.8702	2.7291	15.8424	20.9831
	(α, β, a, b)	(0.1116) (1.4434) (0.0287) (1.9301)
	KwLx	92.56841 67.6157 145.5810 77.3110
	(α, β, a, b)	(1.3827) (2.9101) (1.9805) (5.0347)
	BLx	87.7305	39.5196	8.0644	155.9390
	(α, β, a, b)	(0.4855) (1.9286) (1.8289) (2.1959)
	ELx	288.9114 129.9128 145.7873
	(α, β, a)	(12.9229) (8.5911) (6.9498)

Empirical and theoretical dens.

Email address: a farrukhjamalmphil@gmail.com, b christophe.chesneau@unicaen.fr (Farrukh Jamal a , Christophe Chesneau b , Mohammed Elgarhy c )

Applications

In this section, we prove empirically the flexibility of the TIIGIE-BIII distributions by means of three real data sets. The TIIGIE-BIII distribution will be compared with some competitive models listed in Table 5. The pdfs of these models are given in Appendix. We consider the -(where the maximized log-likelihood), AIC (Akaike information criterion), BIC (Bayesian information criterion), CVM (Cramér-Von Mises), AD (Anderson-Darling) and KS (Kolmogorov Smirnov with its p-value (PV)) statistics to compare the fitted distributions. The results in this section are obtained using the R PROGRAM. Data set 1 : Mead (2016) used actual taxes data set. The data consists of the monthly actual taxes revenue in Egypt from January 2006 to November 2010. The distribution is highly skewed to the right. The actual taxes revenue data (in 1000 million Egyptian pounds) are: 5.9, 20.4, 14.9, 16.2, 17.2, 7.8, 6.1, 9.2, 10.2, 9.6, 13.3, 8.5, 21.6, 18.5,5.1,6.7, 17, 8.6, 9.7, 39.2, 35.7, 15.7, 9.7, 10, 4.1, 36, 8.5, 8, 9.2, 26.2,21.9,16.7, 21.3, 35.4, 14.3, 8.5, 10.6, 19.1, 20.5, 7.1, 7.7, 18.1, 16.5, 11.9, 7,8.6,12.5, 10.3, 11.2, 6.1, 8.4, 11, 11.6, 11.9, 5.2, 6.8, 8.9, 7.1, 10.8. Data set 2 : The data contain Failure stresses of bundles of 1000 impregnated carbon fibers Length 20 mm from Crowder et al. (1991). 2.526, 2.546, 2.628, 2.628, 2.669, 2.669, 2.71, 2.731, 2.731, 2.731, 2.752, 2.752, 2.793, 2.834, 2.834, 2.854, 2.875, 2.875, 2.895, 2.916, 2.916, 2.957, 2.977, 2.998, 3.06, 3.06, 3.06, 3.08.