J P Amartin

J P Josien

PREVENTTVE FACILITIES AND EMERGENCY OPERATIONS IN CASE OFFIRES IN CdF COAL MINES

Charbonnages de France are presently wüming in Lorraine coalfield seams presenüng a high level ofspontaneous conibustion risk associated with high methane emissions. The operators have had to establish a strategy' for spontaneous combustion prevention and fight based on several parameters -linütation ofparasiticairleaks.

-rapid Isolation ofdistricts, -üierting with nitrogen. Several conibustions have tumed into gas ignition in the goa£ hopefillly limited to the upper intersection ofthe longwall face with the air return. As in such cases, after evacuation ofpersoimel, rescue leams were used to rapidly isolate the afleeted dtstricts by closing doors previously erected for this purpose, on the intake and retuni roadways. After Isolation and m order losuppress the fire, the districts have beenkqitunder an inert atmosphere for a period of almost two months. They have then been reopened in aecordance 'with a stricl methodology. It has been possjble then to resume coal winning, which has cor.tmued until the nonnal tennination of

l. INTRODUCTION

The coalfield worked by Houilleres du Bassin de Lorraine (HBL), Charbonnages de France Group,. is situated in north-eastem France. This is a dense deposit made up of two groups of over ten seams of coal, separated by a conglomerate bed nearly 200 metres thick (figure l).

The upper group consists of a bituminous soft coal, the lower coke coal.

The field is sharply folded along two parallel axes resulting in varied types of workings within the following two ranges of dip:

• flat seams frorn 8 to 30° or even more locally; • steep seams from 60 to 80°T his configuration, together with the great depth (over 1,000 metres) leads locally to substantial Stresses causing fairly considerable initial fissuring or even öacturing of the rock mass.

The coalfield is also. fairly gassy, with specific einissions of 50 to 60 m^ of CH4 released per tonne of coal extracted being not uncommon, necessitating very high aüflows in the workings (conunonly > 40 rn^/s per face in Operation) for adequate dilution of the gas.

Finally some seams have a high concentration of-pyrite which facilitates oxidation and accelerates the process of spontaneous heating.

The principal mining method is the retreat longwall technique with caving in the flat seam which accouht for more than 80% of the Output (figure 2).

On the whole therefore, HBL coalfield has the following characteristics:

• fissuring related to depth; • easy oxidation related to the presence of pyrite; • high methane content of coal.

• Working the field requires so the use of large quantities ofair.

27th International Conference of Safety in Mines Research Institutes, New Delhi, 20-22 fevrier 1997, p. 443-456 Hence all the factors that can lead to spontaneous heating and methane ignition are present

In fact m the past HBL have had to deal with many spontaneous combustion tuming into fires, often leading to the definitive closure of a district and the abandonment of its equipmeüt

Over the last seven years, besides the preventive measures taken continuously during mining, HBL have defined a pnevention strategy whereby, in the most recent incidents, the district involved has been rapidly closed off, and then inerted. It was possible subsequently to reopen the district without difiiculty, and to rcsume operations which then continued to the end of the original schedule.

For this strategy to be effecäve it must be possible for the district to be isolated very quiddy; this technique is based on two main requirements:

• closure Systems that are easy to operate; • rapid inerting by injecting nitrogen into the old workings.

FAST-ACTING DBTMCT CLOSURE SYSTEMS

Various arrangements for rapidly closing off districts have been tested, but here are described only the two methods that have been adopted for HBL strategy.

Water stopping (figure 3)

To implement a strategy based uppn the rapid closure of a district by creating a water stopping, a particular configuration is required:

• There must be low points in the district access roadways with a high enough change in level. In fact this drop (d) must exceed the height of the gallery (h) if the water stopping is to be able to cope with the pressure surge resulüng either from any explosion that occurs or simpty from the additional gas introduced by nitrogen inerting;

* Large quantities of water must be rapidly available, because the volume of the water stoppings may reach a few thousand cubic metres and must be established äs quickly äs possible.

For this method of closure to be used. arrangements must be made when the district workings are still at the design stage. In fact the choice is either: * To make use of a natural change in level in the seam itself; * Or to create a low point, not stricüy necessary for normal working of the panel, in both bottom and tail gätes. Also, the water supplies to the district must be sufficient to provide the necessary flow, with the control valves positioned so that they remain accessible from safe areas where the atmosphere is clean.

Stoppings with fast-dosing doors

Erection ofthese stoppings must be completed before work begins in the district.

They do not in fact represent additional expenditure because they are erected outside the panel to be worked and will serve äs definitive closure dams when work on the panel is complete and equipment has been removed.

The stoppings should be erected on both access roadways and be identical in principle but with specific arrangements to : * make allowance for the different equipment used in these roadways; * facilitate any subsequent reopening bf the district. a) Stopping in the retum gate (figure 4) This is a solid stopping, built with good quality anchoring, made of concrete or anhydrite about 3 metres thick.

It is provided with apertures whose total area is related to the Ventilation flow rate (of the order of 35 to 45 n?/s) necessary in the district in order to dilute the methane and to maintain a breathable atmosphere for the wölkers. Hence there is: * l door (1.9 x 2.6 m) capable of withstanding a pressure surge of at least 2 atmospheres operated remotely by compressed air, on each face of the stopping. The door dimensions must be compatible with the space taken up by a monorail that is easily removable or against which the doors dose. * 2 explosion panels, 800 nun or even l metre in diameter, together'with their flaps and a mechanical System for remote opening and closing (to facilitate reopening of the district). The electric cables passing through the stopping are bundled inside pipes which are carefülly sealed with a resin material. * 2 or 3 short pipes let open on the side facing "the coalface for subsequent operations of:

checking the atmosphere behind the stopping; * inerting; * injecting sealing materials. •

In certain circumstances it is worth providing a water drain. This must then terminate in a gooseneck on the accessible side ofthe dam. b) Stopping in the intake gate (figure 5) This has the same geometry äs the dam in the retum gate, although since there is usually a conveyor belt in this roadway the stopping has to be provided with a special System for sealing off the conveyor compartment. This System may be: * a pneumatically operated shutter that grips the conveyor belt in such a way äs to limit Ventilation leakages. When this shutter has been closed it is then necessary to finish its sealing offusing expanded foam; * a shutter fitted by hand after partially dismanüing the conveyor. The other Systems on this stopping will be the same äs those on the stopping in the retum gate, i.e.: * door with space for the monorail; * explosion panel (a single panel is suffident); pipes for instruments, water removal, injection or inerting.

As regards the stopping in the air intake roadway, the closure Systems can either be operated remotely (compressed air) or manually.

The pneumatic controls are grouped on control units which can be installed a few tens of metres from the doors in an area where there is no risk ofbad air.

EMERGENCY INTERVENTIONS

Despite the precautions taken during mimng operations such äs: * strict control of Ventilation; conünuous monitoring of CH4, CO and 02 concentrations in the goaf and retum roadway; * limiting or eliminating air leaks; * pressure balancing; * preventive inerting.

HBL during the last six years have had to deal with a number of critical situations (selfheating in old woridngs (°x)), two of which resulted in the ignition of gas in the goaf of active faces. This led us to cany out six emergency district dosures : these were subsequently reopened in accordance with a careftilly planned procedure.

Evacuation of workers

In the event of a rapidly developing incident, for example a localised gas fire with the release oflarge amounts oftoxic gases or even smoke, it is essential to be able to contact all workers äs quickly äs possible and give them precise evacuation instroctions. This is greatly facilitated by three factors : a) The installation of a large number of sensors in the galleries providing regulär data, at very short time intervals, of changes in the concentrations of gases (CH4, CO, O» ) äs well äs any changes in Ventilation flow, all of which are essential for the extent of the incident to be evaluated from the control and monitoring Station at the surface. These instruments also make it possible to select the best route for withdrawing workers to a safe area. Thus in each HBL mine, it is possible to obtain information frcm about 150 sensors in less than four minutes; these data are processed by Computer and easüy retrieved in the form oftrend curves. b) The use of a dedicated telephone network capable of continuously broadcasting prerecorded waming messages or instructions appropriate to the circumstances, via a simple selection procedure on a microcomputer in the control room. Each HBL mine is fitted with such a network (which is in fact a public address System) which carries out a daily seiftest c) Workers are provided with individual seif rescuers (closed circuit oxygen producüon) with a sufficient reserve for them to reach the fresh air muster points or slightly pressurized shelters.

Thus every Charbonnages de France miner always carries bis seif rescuer (Mini-Apeva) capable of supplying him with 30 minutes of oxygen, sufficient for him to reach the pressurised shelters arranged along the roads giving access to the woridngs (figure 6). The locations of these shelters are determined on the basis of a prior study of evacuation pathways using simulations of the Ventilation System, under normal working conditions and in the event offire.

Intervention by rescue teams

In order to isolate the district in the safest possible conditions, it is necessary first -* to expose the smallest possible number of workers; * to take action quickly.

For these two reasons it is important to recroit well-trained rescue teams and to equip them with flameproof suits, high-capacity breathing apparatus and even refrigerated suits if the air is hot and, naturally, portable monitoring equipment.

These rescue teams will go into action simultaneously, being co-ordinated by the unit manager in the surveillance and control room on the suriäce.

The rescue leams working in the intake airway can usually work in a breathable atmosphere with no need to wear breathing apparatus; however the team in the retum airway will face a polluted atmosphere and sometimes have to work in reduced visibility. In these circumstances the rate of progress will be dictated by whoever is working in the more difficult conditions.

3.3Closure of doors and isolation of the district

Once the decision has been taken to intake the district, the task given to the rescue teams is therefore simultaneously to dose the doors in the previousty erected stoppings.

The Job ofthe team in the intake airway will therefore be to clear away the monoraü installation (dismantle a rail and cutting the cable), to dose the shutter on the service conveyor or dismantle it and seal the compartment using expanded foam, to dean up the door surroundings to ensure that it will dose easily, and to check the functioning of the explosion panel.

Besides dealing with the service conveyor, the team in the retum gate will have the same preparatory work to do, occasionally in conditions of reduced visibility, wearing breathing apparatus and continuously monitoring the atmosphere, the results being immediatdy transmitted to the control room.

To facilitate the task of these teams, the equipment necessary. for this preparatory work is placed in Containers which are instaüed near the stoppings together with pumps and Stocks of sealing resin.

The teams make regulär progress reports to the control room so that the manager in Charge of the operations can take whatever This instruction is given by the manager in Charge ofthe operations in such a way that the district is closed off simultaneously at the two stoppings. following an increase in the amount of nitrogen injected into the district if necessary.

The district is isolated by closing the large doors (pneumatic remote Operation), and then finally by closing the explosion panels.

During actual events at HBL, Ihe preparatory work for the closures carried out by the rescue teams required between one hour and five hours (when there was smoke in the return airway); simultaneous closure of the doors took five to ten minutes.

Thus it is possible in less than eight hours to isolate a district in which there has been severe seif heating and buming gases, with minimum exposure to workers and with maximum precautions.

Inerting

Once the district has been isolated by closing the doors and panels in the stoppmgs, the inerting procedure shculd be carried out äs quickly äs possible by : * Continuing to seal off the stoppings (at the anchoring points, door frames, pipe penetrations, etc.) by spraying sealant; * Adjusting the nitrogen flow to maintain an overpressure at the stoppings on both roadways; * Assisting the removal of heat by establishing an extraction flow, if possible, .compensating for the excess nitrogen injected; * Improving the sealing of the district by additional injection sealing, etc.

REOPENING AN ISOLATED DISTMCT

Monitoring the district

When a district has been isolated for the reasons given. above, changes in its atinosphere should be monitored frequently from samples taken from the pipes penetrating the stopping in the retum gate.

Depending on the changes occurring in the difierent parameters monitored, a number of other Jobs, of varying extent, will haye to be carried out while a reopening strategy is prepared.

The reopening strategy

The reopening strategy should be based upon two factors:

a) Ventilation caiculations

These caiculations should be carried out beforehand to give guidance during the reopening Operation. They should also indicate the methodology to be adopted for closing the district again without risk if a fürther uncontrolled incident should occur during the reopening Operation.

They are based upon the need to avoid the effects, first, of any sudden impact by too great a natural pressure surge on the stopping in the retum gate and, secondly, too high an air flow during the reopening process. The caiculations determine in advance : * the extent to which work in the mine should be halted even before operations commence; * the possible flow rate through the district, together with the pressure distribution after every change made at the stoppings (door, panels, etc.); * the size of the optimal open section to achieve the selected flow rate on coropletion ofthe Operation. This will involve : * installing an adequate number of wellplaced sensors to monitor changes in the atmosphere. within the confined zone and in the primary Ventilation System during purging ofthedistnct; * installing a System to facilitate decompression of the stopping in the retum gate. and the purging of the confined zone. This System may consist of a linc of rigid ducts connected to one of the explosion panels, terminatiRg at its other end in a diluting duct operated remotely (figure 9).

c) Reopening the distrkt

District reopening is supemsed from the control room at the surfäce havmg access in real time to all the Information provided by the sensors installed for monitoring the Operation. Two teams are needed, both fülly briefed on the procedure to be followed: one at the stopping in the intake, the other at the stopping in the retum end. Both teams are in telephone communication with the conu-ol centre and act only on its instructions.

The reopening procedure is in five stages l. Passing on the overpressure at the stopping to a purgc System exposed to primary Ventilation This is done by opening the flap in the explosion panel to which the previously installed purging line is connected. In this way the high pressure is fed directiy to the diluting duct exposed to substantial primary Ventilation.

Decomprcssing the confined zone

This involves simply opening the diluting duct at the end of the line. By this stage the overpressure (500 Pa WG) has disappeared, a small amount of gas has been driven out. but äs yet no air is circulating in the distnct.

Controlled purge ofthe distnct

The order is then given IQ the team located at the stopping in the intake airway to open the flap on the explosion panel-Once this is done. air begins to circulate naturally in the district and will drive the gases to the outlet System. The team at the retum side will then adjust the position of the diluting duct so äs to limit the flow in the district and observe the ruies govenung the concentrations linüts in the primary Ventilation downstream of the stopping. When the diluting duct is fülly open this indicates that preliminary degassing of the district is complete.

Increasing the flow

The degassing process must then be continued by gradually increasing the flow of air into the district; the rate ofincrease in flow is öontinuously monitored by a series of sensors located in the Ventilation System in the air return.

The following operations are carried out in the order shown: * the large door in the stopping in the intake airway is opened; * the second explosion panel at the air retum side is opened: * finally the door at the air retum is opened and the flow adjusted to the desired value by installing an air door having the aperture section previously caiculated.

Throughout this Operation, nitrogen injection and methane drainage must be constantly adjusted. It is also important to check regularly that the flow entering the district for a given configuration agrees with the prior caiculations for the same configuration.

Inspection of the district

This can be done once a flow of fresh air is esiablished and the atmosphere in the retum roadway is breathable. This inspection should follow the direction of normal airflow.

CONCLUSIONS

The strategy developed at HBL. based upon the construcüon of stoppings fitted with fastclosing. remolely-conu-olled doors at the ends of the roadways serving areas of longwall faces with'Sublevel recovery has proved ifs worth.

In the last six years, these arrangements have been used on six occasions, for closing off districts in which serious incidents had occurred, induding two gas fires with considerable release oftoxic gases and smoke. In every case, these arrangements have allowed the workings to be reopened easily and operations to restart. This strategy is supplemented by: • a reliable and effecüve warning System; • a proven seif rescuer in which the workers have confidence;

Figure l :

 l Figure l: Section of the Houilleres du Bassin de Lorraine coalHeld

  Figure 2: Diagram of a mcchanised retreating longwall facc

  Water trap made uslhg access roads

  Figure 3: Rapid closure of a distdct by watcr trap

  modifying the flow of inerting nitrogen, and so on) in conjunction with the engineers in the Ventilation Computer centre who carry out the necessary simulations in real time.When the teams have completed their preparatory work. they retreat to points designated by the operations manager, report in and await the closure Order.

Figure 7 :

 7 Figure 7: Schematic diagram: PEGASSING A FACE AREA.

  Figure 7: Schemattc diagra.m: degassmg a face area,


P Bigarre

M Bennani

SEISMIC MONITORING APPLIED TO MINES SAFETY AND OPTIMAL DESIGN OF MINE LAYOUTS IN HARD ROCK MASS SITÜATIONS

The paper intends to show how seismic data can be usefülly mtegrated in daily, mining operations and contribute to make mining safer in tenns of daily decisions ofmine managemenl äs well äs plaoning ofmine layout design. Thus, case examples from French coal mines facing rockbürsting Problems in longwall faces are discussed and detailed.

l. INTRODUCTION

Mining operations in hard rock, tabular situations are usually associated with induced, seismic activity, i.e. the occurrence of dynamic rock fäilures ranging from surrounding audible rock noises to large scale failure phenomena äs: air blasts in tabular shallow nüning conditions, rockbursts along working faces or roadways in deep, longwall extraction, sudden safety pillar fäilures or tectonic fault slips. Cleariy. seismic behaviour of the rock mass and orebody creates hazardous working conditions for people Underground, damage to equipment and delays to production.

Severity of rockbürsting is affected by many parameters. However, high induced stress changes are a basic parameter of large scale disruptions in the rock mass, and this happens to be especially the case in tabular mining. The openings and large voids created by the mining operations cause pre-existing force to be redistributed and largely maglüfied in the Spund rock. When the stress exceed the strength of the rock mass, generally in the vicinity of the openings, failure occurs, generating a seismic ^tremor, noticeable or not. Only few seismic events give place to damages and are then classified äs rockbursts.

Rockbursts and air-blasts are one of the major problem for the mining industry all around the worid. These are special seismic events and their mechanisms are not clearly identified wheri compared to those öf other seismic events. In other words, exact rock-bürsting prediction in the time-space domain still remains a challenge open to the research area. However, better understanding of those phenomena can be approached within details using available methods like field stress measurements, numerical modelling, and moreover seismic monitoring.

In situ stress measurements

Natural stress measurements are important in the way that induced Stresses (that triggers eventually failure mechanisms in the rock mass initially at equilibrium) around the mine layout are directiy dependent of the initial, in situ stresses in the virgin rock mass; depending on the tectonics of the site and the vicinity of old workings. natural stresses may change considerably in the same zone (presence of folding or faulting zones, geological anomaties). Overcöring technique is a stress relief method that call be applied on a quasi routine basis in Order to map with accuracy Stresses from place to place following the mine layout Underground and its tectonic changes.

Overcöring technique is based on the measurements recorded from a cell equipped with 'strain gauges. After being placed in a smäil borehole away from a mine opening (Fig. l) overcoring ofthe cell induces ftill destressing of the core and strains recorded are used to back caiculate natural stresses (CTJ, c^, 03), in terms ofmodule and orientation.

Numerical modelling

Numerical modelling is a basic approach' when stress prediction is needed inside a large volume ofrock. As soon äs voids and openings are large, or dose together with ancient mine workings overiying or underlying the volume of interest, potential instabilities become very difScult to assess.

When field stress measurements are available and elastic parameters of the bedding rocks are known, numerical modelling offers a low-cost, efficient means to study and identify stability conditions following a lot of situations encountered in practice in mining engineeiing concemed with mine layouts and critical stability problems: stability of pillars in pillarand-room extraction, stability of upper roof in shallow tabular mining method, influence of older panels lying in the vicinity of new mined out openings.

However, if numerical modelling techniques are nowadays widely used, most of the codes marketed present specific linütations when dealing with tabular mine layouts in 3 dimensions, especially when superimposed seam levels are considered by the mining engineer. These geometries are bidimensional, that is ratio of height over length/width of the opening is at least two Orders of magnitude, and stress redistribution over a large volume of rock all aroimd the voids modelled usually provides for non consistent solutions in stress or displacement prediction at least in one direction.

In order to address the particular problem of multilevel tabular orebodies, INERIS has developed the code SUIT3D, based on the boundary element method, and incorporating algorithms and graphic facilities fülly suited to this kind of application. Running on PC based Computer, SUIT3D pennits to model panels extracted by longwall or room-and-pillar methods ofdifferent geometries, distributed on different seam levels with different dip and dip directions. Numbers of panels and seams modelled are restricted to the Computer capacities. Eventually, SürT3D has been optimized to set up sünulations ofdifferent mining sequences.

Fig. 3 shows an application of SUIT3D applied to the modelling of the stress field around and inside a stiff pillar at U.E. Reumaux, Lorraine Collieries ( [START_REF] Heib | Modelisation numerique de l'exploitation du panneau IRMA 1250 -U.E. Reurnaux. Numerical modelling of panel IRMA[END_REF].

Here, the study is based on the project of mining out a future panel crossing from south to north this pillar surrounded by 15 old workings lying in 4 different seams. Analysis of the caiculated stress field pennits to the mine management to decide proper safety countermeasures (local monitoring, destressing action, support design, etc.) to be taken in regards to the past experience and previous modelling results conducted in similar conditions. Seismic tremors, whatever is the ränge of magnitude considered, are not predictable. Thus, seismic monitoring Systems need to be desigried in order to :

(^ provide information äs fast and detailed äs possible on each new, recorded event: reliable real tune processing of wavefonns must include at least caiculated ground motion parameters: peak paräcle velocity, acceleration, duration, frequency, etc., and source parameters: location, energy and magnitude; cf^ pemut to back analyze äs quickly äs possible recent sequences of seismicity m Order to study the space-time patterns related to the different zones of interests äs far äs they are covered by the seismic network.

SEISMIC MONITORING IN FRENCH COLLIERIES

Presentation

Novradays. Charbonnages de France operates two collieries facing rockbursting problems. At Houilleres du Bassin de Lorraine, four mining nnits mine coal down to 1.250 meters deep, following different methods depending on the dip of the seams: caving longwall, mechanised working face or .hydraulically stowed, rising faices: Thickness ranges from 3 to 5 meters. At U.E Provence, the coal deposit lies sub-horizontally, with a thickness varying between 2.2 and 3 metres, $nd is presently mined at a depth of 1.300 meters by the caving longwall face method, equipping panels of 220 meters of face width. Adiustic and seismic monitoring has been a majoi aim since 1985 with research programs to improve rockburst control.

Description of the monitoring neiw'orks

In both Provence and Lorraine coalfields, real-time seismic monitoring is operated in order to analyse and track seismic activity generated by the Underground workings. Mine scale networks (4 and 10 km extension respectively) are equipped äs follows: Local magnitude coverage ranges in both cases from magnitude 1.3 up to 3.2 with sufficient, reliable location estimates (Semäute, 1994). Lower magnitude events are generally poorly located. Most surface seismic stations are equipped with dual gain, increasing largely the dynamic of the System. Digital data acquisnion and processing architectures are showed in Fig. 5 and6.

31
The large extension of the mine scale seismic network deployed in Lorraine (>10km) impose to run local sismo-acousüc networks concentrated around working faces. In this case, stations are equipped with 14 Hz geophones and magnitude covered comes down to 0.

Seismic data interpretation on a routine basis

On a real time basis, seismic data are processed and edited in order to get ground motions and source parameters äs soon äs needed.

4-

• 'iuisrc 4 : l H/, intrinsically sale conditioned cophonc lor coal ininc : Underground horizontal omponcnt lying on ils base. l-'igure 5 : description of thc inonitoring systein of minc scale seismic network operatcd at Provence and Lorraine C'ollierics. Masicr/Slavc prolocol intcgralcd in two parallel acquisition units öfters füll safety in case of inicrvcntion or failurc. Ncwly rccordcd scisniic cvcnt on inastcr acquisition unit is immednrtcly senl t(i Ihc ccniral i.'iiiiipiiti.'i-llin'us;li llic l-nc.il Arca Nclwork l»r real-tinic. auloinutic data procc-ssing, dataliasu u|'idaling ;>nd ri;|ioil ediliiii;. liUL'i'ani\'L'. processing of thc data or inore in dcplli sludics may bc-donc from ccntral Computer or any L-onncL-ied on ihc LAN (Bigarrc. 1994).

-5-Figurc 6: duscriplion Öl monitoring System opei-ating a sisrno-acouslic network in Lorraine Collicry, at cach shal'l cntrance (Fell, 1996).

In practice, at Provence Colliery for example, dozens of seismic'events are recorded daily when mining face ädvance is normal. Special attention is then brought only to the major events, i.e. for which local magnitude is approxünately greater thän 2.5. In this particular case, wavefonns may be interactively checked up and hypocenter location is attentively detailed while information is coUectßd rrom Underground. Figure 7 shows an application of such a routine while rockburst MS-5 was not feit by miners working on T21 face. The mapping detection of large seismic tremors in unexpected, active seisnuc areas is an effident mean to lead to appropriate measures äs quick inspection and fast decisions like floor slotting, with optimal time planning.

UEP HBCM SYTWIS

Secousses du 95 /10 11 ä 5h57 au 95 /10 / 2 ä 23M4 Every month, special seismic report is edited for each currently mined panels, including updated parameters such äs [START_REF] Bigarre | Analyse de ta sismicite miniere pour la planification et le controle des exploitations affectees par des coups de terrains. Data analasis ofmine induced seismic events applied to mine layout design subject to rockburst occurrences[END_REF]: b-value: this widely used parameter is the slope of the linear tendency of the distribution ofthe number of seismic events per class of magnitude (Fig. 8). K merely quantifies how much large seisroic tremors are present in the current seismic activity ofthe zone ofinterest; a-value: this parameter is caiculated äs the linear tendency ofthe daily seismic energy released versus the face ädvance. It quantifies the rate of seismic energy re-lease per meter of ädvance of the panel (Fig. 9). These parameters and others are routinely used to quantify and characterise the seismie activity of the working faces relatively to the coal extraction rate. They are examined in terms of classification of mining areas and panels and for a same panel during its whole nrining phase. .100 600 MO

CumiUulrd Mitviuicp (m)

-F,ST17

-E'i r; i

• EHf:17 -EI'T'';.', ,wrri

Airr".'

"-Atn')'

Fiyurc 10 : plol i)l' i.-umulalivc seismic encrgy rclease vcrsus cumulaied face advance conccrning 9 panels l'runi 4 ilillcrcni.s' (lisuicis .'.pai'iciJ in ilic niinc l.iynut. Lcvcis, ircnds and changes ot seismic rcgiine arc roulinely aiid roliiparalively cxaluincd in order lo flassily currenl woiking laces.

of seismic energy release and then rockbursting potential. Related to Fig. 10, Fig. 11 shows the ratio of total seismic energy release E^ over total face advance LI per panel versus number of recorded rockburst. When correlated with other Parameters, such a basic approach provides for a rough but fast and easy mean to quantüy the susceptibilhy of currently mined panels to rockburst occurrences. Fig. 10 shows the cumulative seismic energy release versus cumulated face advance conceming panel T71, mined out southwards in direction to the shaft pillar (Fig. 12). Rßcorded Small seismic impacts with small associated damages observed Underground, in relation with increase of seismic regime during the previous 140 meters of face advance. permitted to conclude that Stresses were shifting at a higher rate than expected inside the safety zone. These seismic observations did conduct the mine management to the stopping ofthe panel, in order to preserye a sub-critical stress regime in the northem part of the shaft pillar.

/ ^aft

Fig. 12 : localization of sesmic evente inside shaft safety pillar related to panel T71 geometry, following a plane view

Considerations described here-above intend to illustrate the importance of seismic monitoring applied. to mine induced seismic activity and instability control. Seismic routine analysis offers large means to investigate response of Underground structures without considerable efforts and help to save time in terms of safety countermeasures to be taken by mine management.

Seismic prediction

Rockburst are seismic events associated with damages but withoüt other clear distinction and thus very difficult to predict, whether studying intrinsic source mechanisms or patterns of seismic seqüences. [START_REF] Senfaute | La. surveillance microsisnüque de l'exploitation souterraine du charbon aux Houilleres de Provence Microseismic monitoring of Underground coal mining, at the Provence coal mines[END_REF] did study prediction of large seismic tremors (Es ^O^ J) äs a data bank that contain most of the times rockbursts associated with large damages. study based on data analysis technique.s.

The discriminant analysis technique undertaken over the seismic activity induced by the mining of a particular panel has been applied successfully, determining a fünction Y=f(Xl,X2. X3,..) (linear discriminant fünc-tion of N variables); able to classify days with large tremors occurrences and days withoüt, where XI, X2, X3,...X12 are variables äs: face advance the day before and two days before, Logarithm of seismic energy released the day before and two days before, number of seismic tremors per class of magnitude the day before and two days before, etc. More Details can be found in [START_REF] Senfaute | La. surveillance microsisnüque de l'exploitation souterraine du charbon aux Houilleres de Provence Microseismic monitoring of Underground coal mining, at the Provence coal mines[END_REF].

Fig. 13 iUustrates the discriminant fünction «Y» determiried with its value along x axis and number of days along y axis, following two classes (days with and withoüt large tremors). The graph shows cleariy that all the days with large tremors are predicted with positive Y value, except for l day (arrow pointed), i.e. with 97% of rate of success. Days withoüt large tremors are successfully predicted äs such with 76% of rate of success.

-n

Days with large tremor(s) \ Fiyure 13: discriiniiiaiil l'unclion plollcd l'or llic days willi large Iremors (34 days) und withoin larae u'emors (S7 Uays). panel 'l'()7 sKoiiieli'y 'l. U.E. Provence. Success öl prediction öl' days wilh large Irciiiors lo appcar is W/i (Scnlaulc. 19 1 )?).

Days withoüt large tremors n discriniinant fünction

These results show cleariy that data analy-.sis techniques can provide a valuable mean to predict large tremors to occur for the current starting 24 hour window time, äs soon äs others factors (nüning method, geological conditions) remain constant. Research is going on äs to real-time implementation of such indicators.

CONCLUSIONS

Many advances have been made these last years in stress prediction techniques related to tabular mining äs well äs in seismic iristmmentation and processing techniques.

Presently. in French, deep coal basins, stress prediction rely essentially on in situ stress measurements • based on stress relief techniques and then large and/or local scale numerical modelling of the mine layouts, in order to optimise safety and econonücal mining. These two stages are systematically undertaken for each füture panel in ,order to assess and design better countermeasures and mining phases if needed. Results and data are always 9-confronted with the past experience and regulär observational procedures.

As stress prediction from numerical modelling is inherently based on a simplified Problem, real-time seismic monitoring is used to overall control the effective, rock mass response to mining operations. Seismic data are integrated at different time window scales into the set of data collected and interpreted following pre-defined procedures. Objectives are lo improve knowledge äs far and fast äs possible on the working conditions in order to rationalise nüne management decisions confronted with instability occurrences Underground.

Some examples of bn site routine analysis commented here above answer clearly to the need and contribution of rock mechanics engineering techniques, included applied geophysics, in making mining safer in terms of daily decisions of nüne management äs well äs planning of mine layout design. Research carried on focuses on source mechanisms of seismic tremors and fractal analysis of seismic pattems [START_REF] Lanchy | Etüde spatio-temporelle des seismes induits par l'exploitation de la mine de Gardanne -mise en oeuvre d'outils d'analyse et approche fractale. Spacetime distribution ofmine induced seismic events -fractal approach[END_REF] and is hoped to improve in the very near füture reliability ofindications available in terms of major instability prevention.

  Fig. 2 shows plotted strain curves versus drilling advance from data monitored during a measurement trial in Lorraine Collieries, France, 27th International Conference of Safety in Mines Research Institutes, New Delhi, 20-22 fevrier 1997, p. 793-802Accurate stress knowledge represents a key Parameter in order to achieve correct numerical modelling of a mining zone.

FigFig. 2 :

 2 Fig. l: (right) Principle ofovercoring technique

Fig. 3

 3 Fig. 3 : a) orthogonal projection ofthe mine layout around the zone of interest b) contours ofthe ratio of total vertical stress to initial vertical stress a^z /o zzo

  suriäce seismic stations, tridirectional or verücal, with radio transmission of the preamplified modulated Signals towards central site, reception and demodulation, analogue filtering and then digital conversion on an acquisition unit based on PC Computer. All seismic sensors are l Hz geophones (Fig.4); ^ Underground seismic stations, tridirectional or vertical, with cable transmission up to the shaft entrance. Here again, l Hz geophones and preamplifiers are used based on intrinsically safe equipment.

Fig. 7 :

 7 Fig. 7: seisniic activity of 02/10/95, panel T21 district arbois Sud. U.E. Provence.

  Fig 10 shows cumulative seismic energy versus face ädvance plotted for the last 9 panels mined out at Provence between April 92 and August 96. This representation offers large means of classification of panels in terms -6-

Figure 8 :

 8 Figure 8 : Ihc paramctcr h-valuc : 1.94, Islopcl öl' tlic linear Ireiul. indicatcs proporlion öl lui'iic scismic Ircmdrs ovcr ovcrall .scisniii; .iclivily f'or a givcn aclivc niining l'acc.

Figure 9

 9 Figure 9 : a-valuc : l .21 R+05 l/m has bcen ticfined a.s Ihc püsitivc slopc of ihe ilisirihnluin scisinif cncrgy rclc;is>; 011 a il.iily hiisis vcrsus face ativancc. It quaiilil'ics seismic energy releasc rale.

  Fig. 11 : ratio of total sesmic energy release E^ over total face advance LI per panel versus number of recorded rock bürst at UE Provence.

ACKNOWLEDGEMENTS

Author wish to thank management of Houilleres du Bassin .de Lorraine and Houilleres du Bassin de Centre et Midi for the technical support into research programs and permission to publish this paper.


Emmanuel Gobet

Khushboo Surana

A new sequential algorithm for L 2 -approximation and application to Monte-Carlo integration ✩

Keywords: Sequential approximation, Monte-Carlo simulations, Stochastic algorithms, Spectral decomposition 2010 MSC: 65C05, 62Lxx, 74S25

We design a new stochastic algorithm (called SALT) that sequentially approximates a given function in L 2 w.r.t. a probability measure, using a finite sample of the distribution. By increasing the sets of approximating functions and the simulation effort, we compute a L 2 -approximation with higher and higher accuracy. The simulation effort is tuned in a robust way that ensures the convergence under rather general conditions. Then, we apply SALT to build efficient control variates for accurate numerical integration. Examples and numerical experiments support the mathematical analysis.

Introduction

The problem

Given a d-dimensional random variable Y and a measurable function f : R d → R such that E(f 2 (Y )) < +∞, we aim at computing the L 2 decom-✩ The first author research is part of the Chair Financial Risks of the Risk Foundation and the Finance for Energy Market Research Centre. This work has been partly done during the visit of the second author at Ecole Polytechnique during the spring 2013, with the support of International Exchange Programs of Ecole Polytechnique.

α ⋆ Centre de Mathématiques Appliquées, Ecole Polytechnique and CNRS, Route de Saclay, 91128 Palaiseau Cedex, France α ⋆ α ⋆ Indian Institute of Technology, Kanpur, Kalyanpur, U.P., India Email addresses: emmanuel.gobet@polytechnique.edu (Emmanuel Gobet), khushboosurana21@gmail.com (Khushboo Surana) 1 Corresponding author position of f on the vector space spanned by the functions (φ k ) k≥0 according to the law induced by Y , i.e. finding coefficients (α ⋆ k ) k≥0 (from now on, we assume they exist) such that

f (Y ) = k≥0 α ⋆ k φ k (Y ) in L 2 (1.1)
(assuming that E(φ 2 k (Y )) < +∞ for any k). To achieve this decomposition, we use a finite number of independent simulations of Y , say (Y m ) 1≤m≤M , which serve to build a stochastic algorithm to compute (α ⋆ k ) k≥0 . We aim at designing a scheme with robust convergence properties (robust w.r.t. f and the distribution of Y ). Replacing f by an approximation on the basis functions is useful for instance in applications where the evaluation of f is particularly costly (e.g. the output of a complex computer program) and for some reasons, in a further step, we need to evaluate it many times which justifies the use of simpler and cheaper representations. The main application we develop below is the Monte-Carlo computation of E(f (Y )) (numerical integration) using the approximative decomposition (1.1) as control variating. Anther potential field of application is Uncertainty Quantification [START_REF] Le Maître | Spectral Methods for Uncertainty Quantification[END_REF] where the emulation of f serves to build statistics about the uncertainty f (Y ).

Throughout the paper, we assume that (φ k ) k≥0 are orthonormalized basis functions in the following sense:

(H1) E (φ k (Y )φ l (Y )) = δ k,l , for any k, l ≥ 0.

Therefore from (1.1), α ⋆ k = E(f (Y )φ k (Y )) and a naive Monte-Carlo estimator is

α M k = 1 M M m=1 f (Y m )φ k (Y m ). (1.2)
Since f (Y )φ k (Y ) is integrable (because each factor is square integrable), the strong law of large numbers yields the a.s.-convergence of α M k towards α ⋆ k for every k, as M → +∞. In other words

K k≥0 α M k φ k (.)
a.s.

-→

M →+∞ K k≥0 α ⋆ k φ k (.) L 2 (P•Y -1 ) -→ K→+∞ f (.).
Nevertheless, the global L 2 -convergence of the left hand side to the right one is a non-trivial issue since: i) there is no guarantee that each coefficient α M k belongs to L 2 (except under stronger assumptions ensuring that f (Y )φ k (Y ) is square-integrable); ii) as usually in statistics with the analysis of the bias-variance trade-off, we have to appropriately tune the joint convergence of K and M to infinity, in order to let both estimation and statistical errors converge to 0.

In this work we design a sequential algorithm which n-th step provides an estimation of the coefficients (α ⋆ k ) 0≤k≤Kn with non-decreasing order K n . Knowing the values of the former step coefficients helps in improving the evaluation at the next step: the algorithm takes the form of a sequential learning algorithm that we call SALT (Sequential Approximation in L-Two). The convergence is studied in L 2 : therefore, in view of the item i) above, we strengthen the hypothesis on the functions basis by assuming

(H2) for any k ≥ 0, there is a constant c k such 2 that |φ k (Y )| ≤ c k < +∞ a.s. .

Literature background and applications

In the last fifteen years, several important algorithms based on Monte-Carlo L 2 -approximations have been developed, they are popular owing to the flexibility of the Monte-Carlo simulations. Among them, we mention empirical regression methods that are now a standard approach to solve optimal stopping problems [START_REF] Longstaff | Valuing American options by simulation: A simple least squares approach[END_REF][START_REF] Egloff | Monte Carlo algorithms for optimal stopping and statistical learning[END_REF][START_REF] Bouchard | Monte-Carlo valuation of American options: facts and new algorithms to improve existing methods[END_REF], Backward Stochastic Differential Equations [START_REF] Gobet | A regression-based Monte Carlo method to solve backward stochastic differential equations[END_REF][START_REF] Lemor | Rate of convergence of an empirical regression method for solving generalized backward stochastic differential equations[END_REF][START_REF] Gobet | Linear regression MDP scheme for discrete backward stochastic differential equations under general conditions[END_REF]... They are several attempts to improve numerically L 2projections, in order to reduce statistical inaccuracies: see for instance [START_REF] Bender | Least-squares Monte Carlo for BSDEs[END_REF][START_REF] Ben Zineb | Preliminary control variates to improve empirical regression methods[END_REF]. Our current work is in this vein. However, a major difference is that our algorithm SALT works sequentially, in several steps, giving at each step a valuable approximation, which accuracy increases step after step: the computational effort increases coherently according to the number of steps, so that the algorithm can be stopped at any time to produce a solution, expected to be reasonably accurate given the computational effort. For the computation of E(f (Y )) with smooth functions f , we show in Section 3 that our algorithm achieves optimal convergence rates (about optimal numerical integration, see [START_REF] Atanassov | A new Monte Carlo method for calculating integrals of smooth functions[END_REF] or [START_REF] Atanassov | What Monte Carlo models can do and cannot do efficiently?[END_REF] and references therein). Some ideas related to the iteration of our algorithm are exposed in [START_REF] Maire | An iterative computation of approximations on Korobov-like spaces[END_REF][START_REF] Maire | Reducing variance using iterated control variates[END_REF] where iterative control variates are built to speed-up Monte-Carlo integration methods (extensions to stochastic processes are designed in [START_REF] Gobet | Sequential control variates for functionals of Markov processes[END_REF][START_REF] Gobet | Solving BSDE with adaptive control variate[END_REF]). But in [START_REF] Maire | An iterative computation of approximations on Korobov-like spaces[END_REF][START_REF] Maire | Reducing variance using iterated control variates[END_REF], essentially only the numerical integration error (i.e. about computing the expectation) is investigated and the number of control variates is fixed. Here, we consider the full L 2 -approximation and we aim at achieving an infinite number of basis functions, with a robust strategy ensuring the convergence in a rather general setting.

Organization of the paper

In Section 2, we present the algorithm and state the main convergence results. An application to numerical integration (computation of E(f (Y ))) is then developed. In Section 3, examples of basis functions are discussed, with explicit convergence rates in the case of smooth functions f and random variables Y taking values in a cube. Despite these specific examples, we argue that our algorithm applies to more general settings. These examples show that we achieve rate-optimality in the computation of E(f (Y )). Section 4 is devoted to numerical tests showing the performance of the algorithm SALT, in particular for the Monte-Carlo evaluation of E(f (Y )). Intermediate results are proved in Appendix.

Extra notations.

• The residual at order K ≥ 0 of the approximation of f on the basis functions is defined by

r K (y) := f (y) - K k=0 α ⋆ k φ k (y). (1.3)
The identity (1.1) means

E(r 2 K (Y )) = +∞ k=K+1 [α ⋆ k ] 2 → 0 as K → +∞.
The more appropriate the functions basis, the faster the convergence rate of the truncation error E(r 2 K (Y )). Examples are given later.

• The notation A ≤ c B means A ≤ CB for a generic constant C > 0 possibly changing from line to line, independent of the algorithm step number n. Similarly, A ≷ c B means A ≤ c B and B ≤ c A.

Algorithm SALT

Heuristics

The principle of our algorithm is to take advantage of the coefficients computed at steps i = 0, . . . , n -1 to speed-up those at step n, like a sequential learning algorithm.

We first expose the intuition of the algorithm. Consider that (non exact) coefficients (α 0 , . . . , α K ) have been already computed, then write (owing to (H1))

α ⋆ k := α k + β k with β k := E f (Y ) - K l=0 α l φ l (Y ) φ k (Y ) , 0 ≤ k ≤ K,
and replace the expectation by an empirical mean over M independent simulations of Y : this provides an approximated correction to the pre-computed value α k . The variance of this estimator is proportional to 1/M and to the variance of

Z k,K := (f (Y )-K l=0 α l φ l (Y ))φ k (Y ).
In particular, if K l=0 α l φ l (.) approximates well f (.), the random variable Z k,K and its variance are expected to be relatively small. Consequently the new estimation of α ⋆ k will be particularly accurate, and likely much more accurate compared to the previous step. Iterating the procedure and incorporating more and more coefficients to compute, we obtain the algorithm SALT described below.

Detailed algorithm

Let (K n ) n≥1 be a non-decreasing sequence of integers, related to the number of basis functions used at each algorithm step. At step n, we compute the coefficients for the indices k = 0, . . . , K n using M n simulations: the k-th coefficient is to be denoted by α n,Mn k . The algorithm is initialized at n = 0 with α n,Mn k = 0 for any k. The step n ≥ 1 works as follows.

(Step n)-i) Generate a sample (Y n,m ) 1≤m≤Mn of M n independent r.v.
with same distribution as Y , and independent of other simulations.

(Step n)-ii) Define the corrections (β n,Mn k
) k and the coefficients (α n,Mn k ) k as follows:

                     β n,Mn k := 1 M n Mn m=1 f (Y n,m ) - K n-1 l=0 α n-1,M n-1 l φ l (Y n,m ) φ k (Y n,m ) for 0 ≤ k ≤ K n , α n,Mn k := α n-1,M n-1 k + β n,Mn k for 0 ≤ k ≤ K n , β n,Mn k := α n,Mn k := 0 for k > K n .

Convergence results

For n ≥ 1, denote by E n-1 (•) and Var n-1 (•) the expectation and variance conditionally to the sigma-field generated by the random variables

(Y i,m : 1 ≤ m ≤ M i , i ≤ n -1).
To justify the algorithm convergence, we preliminarily analyse the conditional mean and variance of the coefficients. The proof of the above lemma is postponed to Appendix.

Lemma 2.1. Let n ≥ 1 and 0 ≤ k ≤ K n . We have E n-1 (α n,Mn k ) = α ⋆ k , Var n-1 (α n,Mn k ) ≤ c 2 k M n E[r 2 K n-1 (Y )] + K n-1 l=0 (α n-1,M n-1 l -α ⋆ l ) 2 .
At step n of the algorithm, the function f (.) is approximated by

f n (.) := Kn k=0 α n,Mn k φ k (.) (2.4)
and the resulting quadratic error is given by

E n := E f n (Y ) -f (Y ) 2 = E(r 2 Kn (Y )) + Kn k=0 E (α n,Mn k -α ⋆ k ) 2 . (2.5) Note that E 0 = E(f 2 (Y )). From Lemma 2.1, α n,Mn k estimates α ⋆ k without bias and the conditional expectation of α n,Mn k is constant: this implies E (α n,Mn k -α ⋆ k ) 2 = Var(α n,Mn k ) = E Var n-1 (α n,Mn k ) ≤ c 2 k M n E n-1 . (2.6)
Combining (2.5) and (2.6), we establish Theorem 2.1 (Error propagation along algorithm steps). For any n ≥ 1, we have

E n ≤ E(r 2 Kn (Y )) + 1 M n Kn k=0 c 2 k E n-1 .
In case there exists an order K ⋆ < +∞ such that r K ⋆ (.) ≡ 0 (i.e. f can be linearly represented by a finite number of (φ k ) k ), it is enough to take K n = K ⋆ to directly obtain Corollary 2.1 (Finite-dimensional approximation). In the case of finite K ⋆ , we have

E n ≤ 1 M n K ⋆ k=0 c 2 k E n-1 .
In particular, choosing a simulation effort rate which is constant and large enough (in the sense

M n = M with ρ := 1 M K ⋆ k=0 c 2 k < 1)
, the convergence to 0 is geometric along steps, at rate ρ:

E n ≤ ρ n E(f 2 (Y )), n ≥ 0.
The global simulation effort after n steps being C = M n, the error is decreasing at most like exp(-αC) for some α > 0: therefore, when the L 2approximation problem is only finite-dimensional, the convergence is exponentially fast w.r.t. the computational effort. This is fully different from the naive Monte-Carlo procedure. Similar features have been reported in some Sequential Monte-Carlo algorithms, see [START_REF] Booth | Exponential convergence for Monte Carlo particle transport?[END_REF][START_REF] Baggerly | Exponential convergence of adaptive importance sampling for Markov chains[END_REF]. We do not elaborate further in this direction since the focus of this work is rather on the general case K * = +∞.

The next result easily follows from Theorem 2.1, we leave the proof to the reader.

Corollary 2.2 (Infinite-dimensional approximation, general convergence).

The choice

M n := ⌈2 Kn k=0 c 2 k ⌉ yields E n ≤ E(r 2 Kn (Y )) + 1 2 E n-1 , (2.7)
E n ≤ n k=1 2 k-n E(r 2 K k (Y )) + 2 -n E(f 2 (Y )). (2.8)
Thus, for any choice of non-decreasing unbounded sequence (K n ) n (in particular K n → +∞ as n → +∞), we have E(r 2 Kn (Y )) → 0 and thus E n → 0. This specification of M n is very simple and depends only on the a priori knowledge of (c k ) k and the choice of the sequence (K n ) n : this yields robustness concerning the tuning of the local simulation effort.

Application to integration

We now turn to the important application to the computation of E(f (Y )). For this, assume that φ 0 ≡ 1: in view of (H1), it implies that φ k (Y ) is a centered r.v. for k ≥ 1, and

α ⋆ 0 = E(f (Y )). Algorithm 1. A first possibility is to estimate E(f (Y )) simply by I (1)
n := α n,Mn 0 : it is unbiased (see Lemma 2.1) and leads to a quadratic error controlled by E |I (1) n -

E(f (Y ))| 2 = E |α n,Mn 0 -α ⋆ 0 | 2 ≤ E n . (2.9)
Algorithm 2. A second possibility is to generate M other independent samples of Y serving to build a Monte-Carlo estimator of E(f (Y )) with control variates given by f n (Y ) (removing the constant term): it writes

I (2) n := 1 M M m=1 f (Y m ) - Kn k=1 α n,Mn k φ k (Y m ) . (2.10)
As for the first algorithm, it is an unbiased estimator of E(f (Y )) and its quadratic error is given by

E |I (2) n -E(f (Y ))| 2 = 1 M E Var n f (Y ) - Kn k=1 α n,Mn k φ k (Y ) = 1 M Kn k=1 E (α n,Mn k -α ⋆ k ) 2 + +∞ k=Kn+1 [α ⋆ k ] 2 ≤ E n M . (2.11)
Hence, once computed a "good" L 2 -approximation of f , it gives the possibility of more accurately evaluating its integral via Monte-Carlo simulations.

Comparison with a naive Monte-Carlo algorithm

A naive Monte-Carlo algorithm would consist in truncating the decomposition (1.1) at order K and in using M independent copies (Y m ) 1≤m≤M with the same distribution as Y to approximate each coefficient α ⋆ k by α M k defined in (1.2). Then, the approximation of f writes

f Naive MC K,M (.) := K k=0 α M k φ k (.).
(2.12)

The related quadratic error is equal to

E f Naive MC K,M (Y ) -f (Y ) 2 = E(r 2 K (Y )) + K k=0 E (α M k -α ⋆ k ) 2 = E(r 2 K (Y )) + K k=0 Var (φ k (Y )f (Y )) M ≤ E(r 2 K (Y )) + E (f 2 (Y )) M K k=0 c 2 k using (
E f Naive MC K,M (Y ) -f (Y ) 2 ≤ C ξ(K).
This result serves as a benchmark for the theoretical performance of sequential algorithms, see the discussion in Section 3.1.

Examples of Y and vector spaces

Chebyshev polynomials

Let us consider the case of d-dimensional random variable Y which distribution has the density

p(x) = d i=1 1 π 1 -x 2 i 1 (-1,1) d (x) (3.13)
with respect to the Lebesgue measure on R d . Its coordinates (Y 1 , . . . , Y d ) are independent and each one can be simulated owing to the inversion method:

Y i d = sin(πU i -π/2
) where U i is uniformly distributed on (0, 1).

Chebyshev polynomials provide natural basis functions orthonormalized w.r.t. the law of Y . The one-dimensional polynomial of degree k ∈ N is defined by

T 0 (x) = 1 and T k (x) = √ 2 cos(k arccos(x)) for k ≥ 1, (3.14)
so that E(T k (Y )T l (Y )) = δ k,l for any k, l ≥ 0. For usual properties related to Chebyshev polynomials, we refer the reader to [11, Sections 1.5 and 2.2]. These polynomials have the advantage of being bounded (useful to check (H2)). The approximation of f on polynomial basis is often referred to as spectral method [START_REF] Funaro | Polynomial approximation of differential equations[END_REF][START_REF] Le Maître | Spectral Methods for Uncertainty Quantification[END_REF].

The multidimensional case is achieved by tensorization: for a multi-

integers k = (k 1 , . . . , k d ) ∈ N d , set T k (x) = d i=1 T k i (x i ). (3.15)
We readily check that the basis functions (T k ) k∈N d satisfy (H1)and (H2) (with c k = 2 d/2 ).

First strategy for K n . Define |k| = max i=1,...,d k i : at the n-th step of the algorithm, we aim at approximating f by

k:|k|≤n α ⋆ k T k ,
i.e by taking all the polynomials with degree at most n in each variable: therefore,

K n = (n + 1) d . ( 3
(P • Y -1 ), we have E r 2 Kn (Y ) ≤ Cn -2L , n ≥ 1, (3.17)
for a constant C independent on n. Then, the choice M n = ⌈2 k:|k|≤n c 2 k ⌉ = 2 d+1 (n + 1) d leads to the error estimate

E n ≤ c n k=1 2 k-n k -2L + 2 -n ≤ c n -2L
(the second inequality follows by easy computations exploiting the fast increase of 2 k compared to k 2L as k → +∞). The number of simulations used up to step n (representing the simulation effort) is equal

M n := n k=1 M k ≷ c n d+1 , therefore E n ≤ CM -2L d+1 n . (3.18)
We recover the curse of dimensionality effect, with a competition between the smoothness L and the dimension d. We mention that the effective algorithm complexity may be different from M n : it depends on the way the user encodes the evaluation of Chebyshev polynomials and the related finite approximation.

Second strategy for K n . Actually the factor d+1 in (3.18) can be improved into a factor d by choosing another numerical strategy with Exponentially Growing approximation space (denoted by EG or SALT-EG). For γ > 0, set

K n+1 = max(K n + 1, ⌊K n e γ ⌋)
for a given K 1 : the ordering of indices k of the polynomials that we incorporate at each step is not really important in our analysis, for instance in dimension 2 we can take 1,

x 1 , x 2 , x 1 x 2 , x 2 1 , x 2 1 x 2 , x 2 1 x 2 2 , x 1 x 2 2 , x 2 2 , x 3 1 ,
. . . and similarly in higher dimension. On the one hand, it is clear that K n ≤ c e γn . On the other hand, since K n+1 ≥ K n e γ -1 and K n+1 ≥ n + K 1 → +∞, we easily deduce e γn ≤ c K n . Hence, to summarize,

K n ≷ c e γn . ( 3

.19)

Denote by ñ the integer such that (ñ + 1) d ≤ K n < (ñ + 2) d : then, from (3.17) the residual at step n is bounded as

E r 2 Kn (Y ) ≤ C ñ-2L ≷ c e -2L γ d n . (3.20)
Moreover, the number of simulations used up to step n can be estimated as

M n := n k=1 M k ≷ c n k=1 2 d+1 K k ≷ c e γn .
At last, after simple computations, we can show that the quadratic error (2.8) is estimated as follows:

E n ≤ c n k=1 2 k-n e -2L γ d k + 2 -n ≤ c      e -2L γ d n if γ < γ ⋆ := d log(2) 2L , n2 -n if γ = γ ⋆ , 2 -n if γ > γ ⋆ , (3.21) ≤ C        M -2L d n if γ < γ ⋆ , M -2L d n log(M n ) if γ = γ ⋆ , M -log(2)/γ n if γ > γ ⋆ . (3.22)
A better convergence rate is thus obtained for γ small enough, smaller than the threshold γ ⋆ := d log(2) 2L ≈ 0.347 d L related to the smoothness of f and the dimension of Y . Comparing (3.18) and (3.22), observe that SALT-EG yields an improved convergence rate w.r.t. the total simulation effort M n , compared to SALT-PG.

Let us briefly compare the above accuracy result with the one obtained using a naive Monte-Carlo approach (see Section 2.5). From Proposition 2.1 and (3.17), we have ξ(K) = K -2L/d and M = ⌈2 d (K + 1)K 2L/d ⌉, which gives

E f Naive MC K,M (Y ) -f (Y ) 2 ≤ CM -2L d+2L .
The above rate is obviously worse than M -2L d in (3.22).

Application to the computation of E(f (Y )). Applying Algorithm 2 of Section 2.4 using the algorithm SALT-EG to the current setting (see Equality (2.10)) with M = M n extra simulations (which only doubles the total number of simulations), we get an accuracy in the evaluation of E(f (Y )) controlled by (see Inequality (2.11))

E |I

(2)

n -E(f (Y ))| 2 ≤ CM -L d -1 2 n . (3.23)
Of course, this is a very significant improvement compared to the naive Monte-Carlo algorithm for which the above estimate is simply CM

-1 2 n . Remind that for the computation of α ⋆ 0 = E(f (Y ))
when f is a smooth function, the optimal error using M n random evaluations of f is M

-L d -1
(see [START_REF] Atanassov | A new Monte Carlo method for calculating integrals of smooth functions[END_REF] or [START_REF] Atanassov | What Monte Carlo models can do and cannot do efficiently?[END_REF] for a review). Our SALT-EG achieves this optimal rate. We mention that in [START_REF] Atanassov | A new Monte Carlo method for calculating integrals of smooth functions[END_REF], the authors provide a different optimal algorithm in the case where Y is uniformly distributed on (-1, 1) d . Nevertheless, not only our SALT-EG algorithm is rate-optimal for the computation of E(f (Y )), but additionally we do compute all the coefficients sequence (α ⋆ k ) k∈N d and not only α ⋆ 0 , which is a significant improvement.

Legendre polynomials

If Y is now uniformly distributed on the cube (-1, 1) d , a similar analysis can be derived using Legendre polynomials instead of Chebyshev ones. For related properties and approximation results, we refer to [11, Sections 1.4 and 2.2] and [9, Inequality (5.8.11), Section 5.8.2]. The main difference is that the orthonormalized polynomials are no more bounded independently of the degree. Indeed, we can easily check that in dimension 1, the orthonormalized Legendre polynomial L k with degree k is uniformly bounded by

√ k + 1 := c k . In dimension d, L k (x) = L k 1 (x 1 ) . . . L k d (x d ) is bounded by c k = (k 1 + 1) . . . (k d + 1
). As for (3.16) with SALT-PG, we take K n = (n + 1) d (projecting on Legendre polynomials with degree at most n in each variable): then, for a function f with L-derivatives that are square integrable on (-1, 1) d , the rate in (3.17) is unchanged. Moreover, following Corollary 2.2, we take

M n = ⌈2 k:|k|≤n c 2 k ⌉ = 2 (n+1)(n+2) 2 d
, which gives

M n ≷ c n 2d+1
. The quadratic error of the algorithm is then

E n ≤ CM -2L 2d+1 n . (3.24)
Similar computations with the SALT-EG version leads to

E n ≤ c M -L d n .
The application to the computation of E(f (Y )) can be analyzed similarly, we leave the details to the reader.

Observe that the unboundedness of Legendre polynomials has a negative effect on the convergence rate which is similar to doubling the dimension: at least from theoretical viewpoint, the values of (c k ) k in (H2) are an important concern of our algorithm.

Other distributions and approximations

We mention further extensions that we do not fully detail here. Our objective is rather to give, for the reader convenience, alternative ways to explore, that depend much on the examples to handle. Obtaining sparse polynomial approximations (i.e. a reduced size of the vector space while maintaining the same accuracy) is possible for functions f belonging to Korobov spaces, see for instance [START_REF] Maire | An iterative computation of approximations on Korobov-like spaces[END_REF]. Some transformations of the distribution of Y are also proposed to reduce to that case.

In the case of other distributions (in particular supported on non compact sets of R d ), we can use damped polynomials (of Laguerre or Hermite type, the damping ensures that (H2) is met) when we have to handle exponential or Gaussian distributions for Y . Alternatively, one can reduce to one of the previous cases by a change of variable Y = ϕ(Z) (and then compute the L 2 projection of f • ϕ w.r.t. the probability measure P • Z -1 ).

It is also standard to use a one-to-one mapping ϕ : (-1, 1) d → R d , to reduce to the cube; for related discussions on ϕ, see [9, Section 2.7]. Besides, we argue that the cube case (like (-1, 1) d or (0, 1) d ) is quite frequent in probabilistic applications, since usually a simulation algorithm maps a r.v. uniformly distributed on a cube into the target random variable Y .

If the distribution of Z = ϕ -1 (Y ) has a bounded density p Z (.) on (-1, 1) d , then, for any g,

E(g 2 (Y )) = (-1,1) d g 2 (ϕ(z))p Z (z)dz ≤ |p Z | ∞ π d (-1,1) d g 2 (ϕ(z)) dz π d (1 -z 2 1 ) . . . (1 -z 2 d )
, meaning that the L 2 (P • Y -1 )-norm is controlled (up to a constant) by the L 2 -norm w.r.t. Chebyshev weights. Hence, computing the approximation of f • ϕ in the latter norm is sufficient to derive a L 2 approximation of f w.r.t. the law of Y . So far, we have mainly mentioned global polynomial approximations which take well advantage of global smoothness properties, but the use of local approximations (local polynomials for instance) may be relevant in practice. The combination with sparse grids [START_REF] Bungartz | Sparse grids[END_REF] may be a source of significant improvement as well.

Numerical experiments

We present numerical results performed when the distribution of Y has the density (3.13), in dimension d = 1 and d = 2. We use the Chebyshev 2-3-4-5-6) Table 1 Table 2 Table 3 Table 4 Table 5 0.134 0.0727 0.144 0.0338 0.337 polynomials as exposed in Section 3.1, with various SALT-steps n and different strategies for increasing the approximation space (either PG or EG). In the results reported below (Tables 2-3-4-5-6), we indicate:

• The size of the vector space K n at the final step n: for PG it is equal to K n = (n + 1) d , and for EG it is given by K 0 = 1, K n+1 = max(K n + 1, ⌊K n e γ ⌋) for a parameter γ specified later.

• The total amount of simulations M n = n k=1 M k used up to step n. • The quadratic error E n , which is computed using the formula (2.5) as an empirical average of

f n (Y ) -f (Y ) 2 .
The integration over Y is made with M n extra independent simulations, and the one over f n with 10 independent runs of SALT.

• The Alg.2 Error is defined by the half-width of the 95%-confidence interval using the Algorithm 2 for the computation of E(f (Y )), i.e. using the SALT-approximation as a control variate, see the estimator Mn . • The Improvement Factor (IF) is defined by the ratio between the halfwidth of the 95%-confidence interval of the standard Monte-Carlo method using the same number of simulations M n , and Alg.2 Error, whence

IF = 1.96 Var(f (Y )) M n /Alg.2 Error = Var(f (Y ))
Var(I

. Restated in terms of numerical efficiency, obtaining for instance IF = 2 is equivalent to a requirement of 2 2 = 4 times fewer simulations for the same accuracy. The variances Var(f (Y )) have been computed numerically, their values are reported in Table 1.

We first illustrate the phenomenon of geometric convergence stated in Corollary 2.1, by taking the function f as a polynomial (see Tables 2 and5). As soon as the size K n is large enough so that the approximation space includes f , we have r Kn (Y ) = 0 and the estimates of Corollary 2.2 lead to a geometric convergence of the error to 0. This is clearly observed in Tables 2 and5.

Second, when the function f is in none of the approximation spaces, the convergence is presumably not geometric. In Table 3, we take f (x) = (x + ) 2 , for which the regularity may be considered equal to L = 2, yielding γ * ≈ 0.173 (see definition (3.21)). We experiment both strategies PG and EG (γ = 0.05 < γ * ) for increasing the size of approximation space. PG and EG coincide for n small since γ is rather small, and they differ only for n large. Observe that the accuracy is similar after n = 50 steps with PG and n = 40 steps with EG. The advantage of EG here is that we have used a reduced number of simulations, since by increasing faster K n , we increase faster M n and the overall computational effort is smaller. In Table 4 where f (x) = (sin(10x)) + , the observation is similar. We have experienced larger γ, with analogous features, but there is no clear picture to us of which values of γ are better.

Third, we observe that the convergence rate seemingly worsens as the regularity of the function f deteriorates, see Table 4 with f (x) = (sin(10x)) + which is only once (almost-everywhere) differentiable with large derivative (due to the factor 10) and Table 6 with f (x, y) = (x + y) + . SALT still enables to speed-up the computation of E(f (Y )) but the improvement is less significant compared to smooth cases. This is coherent with the previous theoretical analysis, in particular that leading to the estimate (3.23) which illustrates the effect of regularity on the convergence order.

Although the significance of the improvement varies from one example to another, the SALT appears as an interesting scheme for L 2 -approximation and variance reduction, since the tuning of its parameters can be made in a robust way, leading to convergence in quite general situations (see Corollary 2.2). Because of independence between steps,

E n-1 (β n,Mn k ) = E(f (Y )φ k (Y ))- K n-1 l=0 α n-1,M n-1 l E(φ l (Y )φ k (Y )) = α ⋆ k -α n-1,M n-1 k , thus E n-1 (α n,Mn k ) = α n-1,M n-1 k + α ⋆ k -α n-1,M n-1 k = α ⋆ k .
We now turn to the variance analysis. By bounding the variance by the second moment and using (H2), we obtain

  H2) and bounding the variance by the second moment. By appropriately tuning M in order to balance both contributions, we obtain Proposition 2.1. Let K ≥ 0 and assume that E(r 2 K (Y )) ≤ ξ(K) for an explicit positive function ξ(.) that measures the accuracy of the L 2 -approximation w.r.t. the order K. For M := ⌈ K k=0 c 2 k ξ(K) ⌉, the quadratic error of the naive Monte-Carlo estimator (2.12) is at most equal (up to a constant C independent of K) to the truncation error:

( 2 .

 2 10); thus we set Alg.2 Error = 1.96 Var(I (2) n )

Table 6 :

 6 Results in dimension d = 2 for f (x, y) = (x + y) + , using SALT-PG 294 × 10 -2 3.487 × 10 -3 5,13 20 441 26488 1.274 × 10 -3 3.782 × 10 -4 18,5 30 961 83328 3.437 × 10 -4 1.258 × 10 -4 31,3 40 1681 190568 1.569 × 10 -4 5.625 × 10 -5 46,3 Appendix A. Proof of Lemma 2.1

Var n- 1 (

 1 φ l (Y ))φ k (Y )

Table 1 :

 1 Variance of f (Y ) in each case (Tables

Table 2 :

 2 Results in dimensiond = 1 for f (x) = (1x 2 )4 using SALT-PG 10 -12 9.936 × 10 -8 2.38 × 10 5 30 31 1984 2.914 × 10 -17 2.375 × 10 -10 6.78 × 10 7 40 41 3444 1.399 × 10 -20 3.950 × 10 -12 3.10 × 10 9 50 51 5304 2.210 × 10 -22 4.001 × 10 -13 2.46 × 10 10 70 71 10224 3.294 × 10 -25 1.112 × 10 -14 6.38 × 10 11 100 100 20604 3.458 × 10 -28 2.539 × 10 -16 1.97 × 10 13

	n 10 11 K n 20 21	M n 264 3.667 × 10 -5 7.305 × 10 -4 Alg. 2 Error E n 924 2.374 ×	IF 6.05 × 10 1

Table 3 :

 3 Results in dimension d = 1 for f (x) = (x + ) 2 = 0.05) 11 264 9.812 × 10 -43.779 × 10 -3 8,60 20 PG 21 924 1.539 × 10 -4 7.999 × 10 -4 21,7 EG (γ = 0.05) 21 924 1.253 × 10 -4 7.217 × 10 -4 24,1 30 PG 31 1984 4.426 × 10 -5 2.927 × 10 -4 40,5 EG (γ = 0.05) 32 1988 4.741 × 10 -5 3.027 × 10 -4 39,1 40 PG 41 3444 1.508 × 10 -5 1.297 × 10 -4 69,4 EG (γ = 0.05) 53 3712 8.086 × 10 -6 9.147 × 10 -5 94,8 50 PG 51 5304 7.485 × 10 -6 7.118 × 10 -5 102 EG (γ = 0.05) 87 6532 1.868 × 10 -6 3.315 × 10 -5 197

	n 10	SALT PG EG (γ	K n M n 11 264 9.430 × 10 -4 3.704 × 10 -3 8,78 Alg. 2 Error IF E n

Table 4 :

 4 Results in dimension d = 1 for f (x) = (sin(10x)) +

	n 10 20 30 40 50	SALT PG EG (γ = 0.05) 11 264 2.761 × 10 -1 6.338 × 10 -2 0,72 K n M n Alg.2 Error IF E n 11 264 3.107 × 10 -1 6.724 × 10 -2 0,68 PG 21 924 2.193 × 10 -2 9.543 × 10 -3 2,57 EG (γ = 0.05) 21 924 2.500 × 10 -2 1.019 × 10 -2 2,40 PG 31 1984 5.812 × 10 -3 3.355 × 10 -3 4,98 EG (γ = 0.05) 32 1988 6.261 × 10 -3 3.478 × 10 -3 4,80 PG 41 3444 1.898 × 10 -3 1.455 × 10 -3 8,72 EG (γ = 0.05) 53 3712 1.283 × 10 -3 1.115 × 10 -3 11,0 PG 51 5304 1.367 × 10 -3 9.951 × 10 -4 10,3 EG (γ = 0.05) 87 6532 2.411 × 10 -4 3.765 × 10 -4 24,5

Table 5 :

 5 Results in dimensiond = 2 for f (x, y) = (1x 2 ) 4 (1y 2 ) 4 using SALT-PG 10 -51.186 × 10 -4 4.77 × 10 1 20 441 26488 6.604 × 10 -14 3.090 × 10 -9 7.16 × 10 5 30 961 83328 9.212 × 10 -18 2.061 × 10 -11 6.06 × 10 6 40 1681 190568 3.656 × 10 -20 8.585 × 10 -13 9.61 × 10 8

	n 10 121 K n	M n 4048 1.529 × E n	Alg.2 Error	IF

observe that c k ≥ 1 since E(φ 2 k (Y )) = 1; in some cases, c k does not depend on k, see Section

n

using (1.3) and (H2) at the last equality.