F Němec

M Liška

M Parrot

Seismic influence on the VLF transmitter signal intensity measured by the low-altitude satellite DEMETER

We present results of a statistical study of a possible influence of the seismic activity on the intensity of very low frequency (VLF) transmitter signals observed by a low-altitude satellite. Electric field measurements performed by the Detection of Electro-Magnetic Emissions Transmitted from Earthquake Regions (DEMETER) satellite during its entire mission spanning almost 6.5 years were used. Among various VLF transmitter signals detected, we focused particularly on the NWC and JJI transmitters, because of their favorable locations close to seismically active areas. We evaluate the intensities of the detected transmitter signals at the times when they passed in the vicinity of an imminent earthquake during the propagation in the Earth-ionosphere waveguide, and we statistically compare them with the intensities measured at the times when there was no earthquake present. Only earthquakes with magnitudes larger than or equal to 5 and depths shallower than or equal to 40 km were considered in the analysis. Moreover, due to the low intensity of detected transmitter signals during the day, the analysis is limited exclusively to the nightside. Although the amount of relevant data is rather low, the obtained results show that there is a decrease of the detected intensity shortly (0-3 hours) after the times of the main shocks observed both for the NWC and JJI transmitter signals. The effect is spatially rather limited, observed when the signal passes within about 4 degrees from the earthquake epicenter. The intensity decrease appears to be consistent with acoustic-gravity waves propagating from the earthquake region and influencing the bottom of the ionosphere.

Introduction

A considerable number of electromagnetic effects in the ionosphere possibly related to the seismic activity have been reported over the last few decades [1][2][3]. Related theoretical development suggests that these are due to litosphere-atmosphere-ionosphere coupling [START_REF] Pulinets | Ionospheric Precursors of Earthquakes[END_REF][START_REF] Harrison | [END_REF][6][7][8]. Given that they may occur not only after the earthquake occurrence, but also shortly before the time of the main shock, they might -if confirmed and a e-mail: frantisek.nemec@gmail.com understood -eventually serve as short-time precursors of an imminent earthquake [START_REF] Parrot | Electromagnetic Noise Due to Earthquakes, Handbook of Atmospheric Electrodynamics[END_REF][START_REF] Li | [END_REF].

Possible seismic-related electromagnetic effects may occur in various frequency ranges from ultra low frequency (ULF, from 1 mHz to 3 Hz) through extra low frequency (ELF, from 3 Hz to 3 kHz) up to very low frequency (VLF, from 3 kHz to 30 kHz), and on various time scales, from several minutes up to a few months before an earthquake [11][12][13][14][15][START_REF] Parrot | Atmospheric and Ionospheric Electromagnetic Phenomena Associated with Earthquakes, Statistical Studies with Satellite Observations of Seismogenic Effects[END_REF]. However, there is still not an overall agreement on the nature of the observed effects, as both intensity increases [START_REF] Parrot | Electromagnetic Noise Due to Earthquakes, Handbook of Atmospheric Electrodynamics[END_REF]14,[START_REF] Parrot | [END_REF]18] and wave intensity attenuations [19][20][21][22][23] were reported. In fact, the very existence of this exciting phenomenon is still doubtful [24][25][26][27][28], as -if existing -the seismic related effects are very weak and they can be easily hidden in common natural variations of electromagnetic waves.

This complication can be effectively overcome by large statistical studies including many different earthquakes and using a variation of a superposed epoch analysis.

A systematic analysis of ELF/VLF electromagnetic wave intensity measured by the low-altitude DEMETER satellite revealed that there might be a statistically significant decrease of the wave intensity at frequencies of about 1.7 kHz observed shortly (0-4 hours) before the times of earthquakes during the night [20]. This result was confirmed using a larger data set [21], and eventually re-confirmed using the entirety of the DEMETER data measured during the mission duration [22,23]. It was suggested that the observed effect might be due to changes in the Earth-ionosphere waveguide [29,30] and a related increase of the waveguide critical frequency [START_REF] Budden | The Wave-Guide Mode Theory of Wave Propagation[END_REF]. This seems to be consistent with subionospheric VLF signal perturbations possibly related to earthquakes [1,[START_REF] Hayakawa | [END_REF] and to ion density variations reported in low-altitude satellite data [33][34][35][36][37]. A decrease of the VLF transmitter intensity observed by the DEME-TER spacecraft was also reported in the first 2 years of data, and it was speculated that it might be due to the influence of atmospheric gravity waves induced by the gas-water release from the earthquake preparatory zone [19]. Case studies related to powerful earthquakes were further reported [38][39][40][41][42][43][44][45]. These studies typically reported a decrease of detecting the transmitter signal intensity, which might occur as long as several days before the time of the main shock [46,47]. However, occasional signal enhancements were also reported [48].

We investigate possible seismic-related variations of the intensity of VLF transmitter signals observed by the DEMETER spacecraft using the data measured during the entire duration of the mission (about 6.5 years). A variation of the two-step data processing method, which allows us to account for non-seismic related signal variations, is applied [20]. The used data sets are introduced in Section 2. The way of associating the recorded signals with individual earthquakes and subsequent statistical analyses is described in Sections 3 and 4, respectively. Section 5 contains the results of the statistical analysis of the seismic influence on the transmitter signal intensity. These are discussed in Section 6, and they are briefly summarized in Section 7.

Data sets

DEMETER was a low-altitude satellite operating at a circular polar orbit with an altitude of about 700 km (the original altitude of about 710 km was decreased to about 660 km in December 2005). It was active since June 2004 until December 2010. The orbit of the satellite was almost Sun-synchronous; the north-to-south half-orbits corresponded to the local time of about 10:30 hours, while south-to-north half-orbits corresponded to the local time of approximately 22:30 hours. This allows us to conveniently separate the nighttime and the daytime data. DEMETER measurements took place at geomagnetic latitudes between -65 • and 65 The great circle path between the transmitter location and the spacecraft location at a given time is shown by the blue curve. The distance of the particular measurement from the earthquake considered for the performed analysis is then defined as the shortest distance between the earthquake epicenter and the great circle path. (b) A sketch of the transmitter signal propagation along the magnetic field line to the conjugate hemisphere. The transmitter location is shown by the black triangle at the top. Its signal propagates in the Earth-ionosphere waveguide (dark red arrows) and eventually leaks through the ionosphere to higher altitudes (dark green arrows) where it is measured by the DEMETER spacecraft. It further propagates roughly along the magnetic field line to the opposite hemisphere. The signal observed in the conjugate hemisphere may be significantly modified due to wave-particle interactions taking place in the equatorial region.

field measurements were performed. However, VLF magnetic field data suffer from a significant amount of onboard interferences, and only the electric field data are thus used in the present paper. In a continuously active survey mode, the VLF electric field data cover the frequency range from 20 Hz to 20 kHz with a frequency resolution of about 19.5 Hz and a time resolution of about 2 s. More detailed description of the electric field instrument was given in a dedicated instrument paper [49]. We use data from the entire time period when the DEMETER satellite was active (approximately 6.5 years). The Global Centroid-Moment-Tensor catalog (http://www.ldeo. columbia.edu/ ~gcmt/projects/CMT/catalog/jan76_dec13.ndk) was used as a source of the seismic activity data, providing us with the time, magnitude, location, and depth of hypocenter of earthquakes which occurred all over the world during the time period of interest.

Data processing

In order to evaluate a possible seismic influence on the transmitter signal intensity, we suggest a scheme shown in Figure 1a. The geographic location of a considered transmitter is depicted by the black triangle in the lower left part of the figure. The black curve in the right-hand side of the figure corresponds to the projection of the DEMETER spacecraft orbit. The spacecraft location at a selected time is then depicted by the black square. The blue curve shows the great circle path between the transmitter and the spacecraft location, corresponding to the path of the transmitter signal propagation in the Earth-ionosphere waveguide. Assuming that an earthquake occurs close to this path, it may potentially influence the waveguide properties and thus also the wave intensity detected by the spacecraft. The position of a hypothetical earthquake is marked by the red cross. The short black curve connecting the earthquake location and the great circle path represents the shortest distance between the earthquake epicenter and the signal propagation path, which is considered when evaluating whether a given earthquake is close enough to have a detectable effect or not. We note that another possibility of the transmitter signal propagation would be the propagation nearly along the magnetic field lines to the opposite hemisphere. A sketch of the situation is shown in Figure 1b. It shows a part of the Earth, a transmitter (black triangle at the top), and a corresponding magnetic field line. The dashed circular curve corresponds to the DEMETER orbit. A VLF signal generated by the transmitter propagates in the Earth-ionosphere waveguide, as shown by the dark red arrows at the top. The dark green arrows depict an eventual leakage of the transmitter signal through the ionosphere to higher altitudes. There it is detected by DEMETER, as shown by the yellowish ellipse. The signal, however, further propagates, roughly along the magnetic field line, to the opposite hemisphere (light green arrows). There, close to the magnetically conjugate point, it again reaches the DEME-TER altitudes and can be detected by the spacecraft. Finally, the signal may make it down to the Earth-ionosphere waveguide in the conjugate hemisphere and propagate therein (light red arrows). Importantly, as the waves propagate through the equatorial region, their intensity and spectral features can be significantly altered due to wave-particle interactions [50,51]. This vastly complicates a proper characterization of wave intensities in the conjugate region and limits their usability for the analysis of possible seismic-related effects. Only DEMETER satellite measurements performed at the geomagnetic hemisphere of the transmitter are thus used in the present study.

Possible variations of the VLF transmitter signal intensity related to the seismic activity are expected (if existent) to be very weak. In order to distinguish them from other intensity variations of a comparable amplitude, a statistical analysis including many different events is needed. It is thus desirable to identify VLF transmitters whose signal would often encounter an imminent earthquake on its great circle path to DEMETER. The results of this identification of the transmitters most favorable for the intended statistical analysis are shown in Figure 2. We considered the 12 continuously operating VLF transmitters from Table 1 [52], which are readily detectable by the DEMETER spacecraft. Out of these, we disregard the 3 Russian Alpha transmitters, which alternate between three different frequencies, and are thus not optimally suited for our purpose. This leaves us with 9 transmitters possibly suitable for the analysis. In order to estimate the total number of individual earthquakes ("events") possibly influencing the transmitter signal path over the entire DEME-TER mission duration, we assume that: (i) only earthquakes with magnitudes larger than or equal to 5.0 and depths shallower than or equal to 40 km can influence the transmitter intensity, (ii) a given earthquake affects the transmitter signal propagation at great circle paths within 3 degrees from the earthquake epicenter, (iii) the effect lasts for 1 hour around the time of the main shock, and (iv) the effect occurs only during the local night. Although these assumptions are clearly somewhat arbitrary, they provide us for each transmitter with a rough estimate of the number of relevant events. Most importantly, the obtained results are rather independent of the exact assumptions made, as they remain qualitatively the same.

Figure 2 shows the color coded number of events (i.e., number of DEMETER orbits for which the transmitter-DEMETER great circle path was potentially affected by an earthquake) as a function of the transmitter code name (abscissa) and maximum allowed distance between the transmitter and DEMETER measurements at a given time (ordinate). It can be seen that the number of events is generally rather low, with two transmitters standing out clearly. The first of them is the NWC transmitter (latitude -21.82 • , longitude 114.17 • , operating frequency 19.8 kHz), which is the most powerful VLF transmitter operating in the southern hemisphere. The second transmitter resulting in considerable number of events is the JJI transmitter (latitude 32.04 • , longitude 130.81 • , operating frequency 22.2 kHz), which has a lower output power, but it is favorably located close to seismically active regions. The maximum considered distance between the transmitter location and the location of the DEMETER spacecraft clearly strongly affects the results. This parameter expresses the situation that the spatial extent of the transmitter signal is limited, i.e., only satellite measurements not too far from the satellite location should be considered. We express this distance in the units of angular degrees, i.e., one degree corresponds to about 110 km. The actual value of this parameter is clearly related to the transmitter power output, penetration characteristics of the ionosphere, and transmitter signal attenuation in the Earth-ionosphere waveguide.

Although the statistical analysis presented in Section 4 would in principle allow us to combine signal intensity variations from all relevant transmitters together, we opted not to do so. Considering that these transmitters operate at different frequencies (and locations), the natural background variations for each of them would be slightly different, as might be also underlying physical mechanisms responsible for possible observed effects. Taking into account that the number of suitable events for other than the NWC and JJI transmitters is very low, using only the NWC and JJI transmitter data results in disregarding only small amount of data and provides us with a better picture of the situation. Additionally, analyzing the signals of these transmitters separately allows us for an independent check whether the same intensity variations are observed for the two considered transmitters. Only the NWC and JJI transmitter signals are thus used in the analysis.

The geographic locations of the NWC and JJI transmitters are shown in Figure 3 by the black triangle and black diamond, respectively. The black curve shows the geomagnetic equator. The NWC transmitter signal intensity is thus considered as potentially affected by the seismic activity only when measured at the times when DEMETER is located in the southern hemisphere, and the JJI transmitter signal intensity is considered as potentially affected by the seismic activity only when measured at the times when DEMETER is located in the northern hemisphere. Otherwise, the transmitter signal would propagate between the hemispheres primarily along the magnetic field lines, being influenced by wave-particle interactions in the equatorial region, and a possible seismic influence on the great circle path could thus not be properly evaluated. The color points correspond to the locations of individual earthquakes with magnitudes larger than or equal to 5.0 and depths shallower than 40 km that occurred during the DEMETER mission duration, with the magnitude color coded according to the color scale on the right-hand side. It can be seen that both the NWC and JJI transmitters occur in the vicinity of seismically active regions and, consistently with the results from Figure 2, are thus optimally suited for the intended analysis.

Statistical analyses

Geographic maps of NWC and JJI transmitter signal intensities detected by the DEMETER spacecraft are depicted in Figures 4 and5, respectively. Color coded are median power spectral densities in individual 2 penetration through the ionosphere is thus roughly by an order of magnitude higher [54][55][56]. Additionally, the attenuation of a signal propagating in the Earth-ionosphere waveguide is generally larger during the day than during the night [57]. Consequently, the transmitter signal intensity detected by DEMETER is significantly larger during the night than during the day. Further, the transmitter signal intensity is -as expected -largest close to the transmitter locations and geomagnetically conjugated points. The slight shift of the peak intensity from the transmitter location toward the geomagnetic equator is due to the wave propagation between transmitter and the satellite in the direction along the ambient magnetic field rather than vertically [52]. We note that although the DEMETER altitude slightly changed over the mission duration, the signal attenuation at these high altitudes is extremely low, and this minor altitudinal change thus has principally no effect on the detected transmitter signal intensity [53]. Consequently, the full DEMETER data set is used in the presented analysis. We further note that the time periods when the respective transmitters were not operating were determined by investigating the wave intensity close to the transmitter location [52] and excluded from the analysis. Finally, as discussed above, it is necessary to limit the analysis only to latitude-longitude intervals with a reasonable intensity of the transmitter signals. Specifically, we require the median intensity in a given latitude-longitude bin to be larger than 1 µV 2 m -2 Hz -1 . If it is not, then all the measurements performed in a given spatial bin are excluded from the analysis. Although this intensity threshold is set slightly arbitrary, it roughly corresponds to the intensity safely larger than the background noise (see Figs. 4 and5), but low enough in order not to exclude too much data. In any case, applying this intensity threshold removes the vast majority of the daytime data from the analysis, in particular in the case of the JJI transmitter. However, for both transmitters, the amount of intense-enough daytime data is very limited and effectively prevents us from performing the analysis. The subsequent analysis thus uses exclusively the nighttime data.

The detected transmitter signal intensity varies due to many different factors. In order to identify any possible signal intensity variations related to the earthquake occurrence, we would optimally need to account for these in some way. We thus apply a two-step data processing method based on the data processing formerly developed for natural VLF emissions [20]. The basic idea of the method is that, in the first step, we construct a long-term "map" of the wave intensity parameterized by relevant controlling factors, which is used to express "how intense a signal is expected to be at a given location under given conditions". Then, in the second step of the data processing, we evaluate the wave intensities measured in the vicinity of earthquakes by comparing them with the intensities expected according to the first step of the data processing. The actual application of the method to the specific analyzed situation of the NWC and JJI VLF transmitters is described below.

When constructing the long-term map of the transmitter signal intensity, the following parameters are considered as the most relevant controlling factors: (i) geographic longitude of the spacecraft (2 • resolution, i.e., 180 bins), (ii) geographic latitude of the spacecraft (2 • resolution, i.e., 90 bins), (iii) local time (daytime vs nighttime, i.e., 2 bins), and iv) season of the year (January-March, April-June, July-September, October-December, i.e., 4 bins). We have also investigated a possible influence of the geomagnetic activity represented by the Kp index, but principally no relation between the transmitter signal intensity and the geomagnetic activity was found (not shown, consistent with the statistical results [52]). The level of the geomagnetic activity is thus further not considered in the presented analysis. In order to characterize the transmitter signal intensities measured by DEMETER at given locations under given conditions, a histogram of the transmitter signal intensities measured for each combination of the four parameters was stored with a 0.1 log(µV 2 m -2 Hz -1 ) resolution.

A crucial point of this first step of the performed processing is that it allows us to evaluate whether the intensity measured at a given place under given conditions (i.e., day/night and season) is lower or larger than the long-term average and by how much. This can be conveniently expressed by a corresponding cumulative probability, i.e., a ratio of the measurements which took place at the same place and under the same conditions which are less intense than the actually measured value [20]. The values of the cumulative probability equal to 0 and 1 correspond to the lowest and largest intensities ever measured at a given location under given conditions, respectively. The value of the cumulative probability equal to 0.5 corresponds to a situation of the median intensity, when exactly one half of the intensities measured at a given place and under given conditions was less intense and one half more intense. Importantly, independently of the original intensity distribution, the values of the cumulative probability are uniformly distributed [20]. This allows us not only to compare the measured intensities with the expected long-term distribution, but also to conveniently evaluate the associated statistics.

In the second step of the data processing, we use the map of the transmitter signal intensities constructed in the first step and the associated cumulative probabilities to investigate whether there is a variation of the signal intensity related to the earthquake occurrence or not. The crucial point in this analysis is to determine whether the transmitter signal is located "close to an earthquake" or not. As discussed above, we consider a picture of an imminent earthquake possibly influencing the electron ionization/temperature at the bottom of the ionosphere, modifying thus inherently the Earth-ionosphere waveguide. The transmitter signal propagating therein would be modified, and its intensity change as compared to normal conditions might be detected. In order to possibly influence the transmitter signal intensity, the earthquake epicenter thus has to be close to the transmitter-DEMETER great circle path. This simple picture assumes the transmitter signal to propagate principally all the way to the spacecraft in the Earth-ionosphere waveguide, i.e., we disregard the signals which were ducted along the ambient magnetic field between the hemispheres. Consequently, only DEMETER measurements performed in the same hemisphere as the transmitter are considered for the evaluation of possible seismic-related effects.

In order to evaluate signal intensity variations possibly related to the seismic activity, we considered the time interval from 72 hours before the main shock times until 24 hours after the main shock times. Although seismic related effects might possibly occur more than three days before the main shock, one should consider that the sequence of processes on these longer time scales is surely slightly different for each earthquake. The effect would then occur at a slightly different time for each event, and the superposed epoch analysis would average out the effect. Also, should the considered time window be too long, too much data would be eliminated from the analysis in an attempt not to mix pre-and post-seismic effects (see below). Otherwise, the particular choice of the considered time window does not affect the obtained results. Only earthquakes with magnitudes larger than or equal to 5.0 and with depths of hypocenters shallower than or equal to 40 km were included in the analysis. This choice follows former results [20], and it is consistent with the intuitive expectations that no effect is to be observed for weak/very deep earthquakes. We note that any other (even non-seismic events) possibly affecting the detected transmitter signal intensity, which are not included in the analysis, would be randomly distributed with respect to the considered time intervals, i.e., they only introduce a random "noise" in the performed analysis.

DEMETER data for which an earthquake epicenter is within 10 degrees from the great circle path connecting the transmitter and the DEMETER location are then evaluated. If -for a given spacecraft position -more than one earthquake was close to the transmitter-DEMETER great circle path within the considered 96 hours long time window, we discard the data from the analysis. This is done in order not to mix pre-and post-seismic effects, and, moreover, to prevent a single piece of data to be included in the analysis several times [20]. The relevant data are then organized as a function of the distance between the earthquake epicenter and the transmitter-DEMETER great circle path and the time to/from an earthquake. We use 1 hour time resolution (i.e., 96 time bins) and 1 degree resolution in distance (i.e., 10 distance bins). We then define a "probabilistic intensity" in each 96 × 10 bins as an average value of the respective cumulative probabilities minus 0.5. All relevant DEMETER data are used in the calculation. Finally, we normalize the probabilistic intensities by their standard deviations to obtain "normalized probabilistic intensities" [20]. The resulting normalized probabilistic intensities allow us to directly evaluate whether the intensity in a given bin is less/more intense than average, and how statistically significant this effect is in terms of the number of standard deviations. This normalization benefits from the fact that according to the central limit theorem, the probabilistic intensities have a nearly normal distribution. Although an infinite number of cumulative probabilities need to be averaged for the probabilistic intensities to follow exactly a normal distribution, it turns out that averaging of only three cumulative probabilities results in the distribution of probabilistic intensities very close to normal (with an obvious limitation of extremely low/large probabilistic intensities, which cannot occur in our statistics, albeit the normal distribution predicts nonzero probability density for them).

The validity of the performed analysis was verified using surrogate earthquake catalogs. We started with the real earthquake catalog, and we constructed in total 24 surrogate catalogs by shifting the times of individual earthquakes by ±30, ±60, ±90, ±120, ±150, ±180, ±210, ±240, ±270, ±300, ±330, and ±360 days. Taking into account the considered time scales of possible effects related to the seismic activity, the applied time shifts are large enough to ensure that the data measured close to these surrogate virtual earthquakes should not be in any way affected by the seismic activity. Consequently, the main aim is to verify that the values of normalized probabilistic intensity indeed follow a normal distribution, and to demonstrate the noise level in resulting distance-time plots. The results obtained for these surrogate earthquake catalogs and the NWC transmitter are shown in Figure 6. Figure 6a shows normalized probabilistic intensities obtained for a single surrogate earthquake catalog. The results obtained for other surrogate earthquake catalogs are qualitatively the same. The normalized probabilistic intensities are color coded as a function of the time relative to the times of the (virtual) main shocks (abscissa) and distance of the NWC transmitter-DEMETER great circle paths from the respective earthquake epicenters (ordinate). As expected, one can see a significant scatter of the normalized probabilistic intensity values, both positive and negative, as corresponds to the Gaussian distribution they should follow. Considering that a surrogate earthquake catalog was used to obtain this plot, the observed variations are in no way related to the seismic activity, and they are strictly a statistical fluctuation. Figure 6b demonstrates that the values of normalized probabilistic intensity indeed follow a Gaussian distribution with expected parameters. It was obtained by combining the normalized probabilistic intensities corresponding to all 24 surrogate earthquake catalogs together and plotting the resulting distribution. Overplotted by the red dashed line is a normal distribution with the mean value equal to 0 and the standard deviation equal to 1, demonstrating that the values of normalized probabilistic intensities obtained for surrogate earthquake catalogs indeed follow this expected distribution.

Results

The resulting plot of the probabilistic intensity obtained for the NWC transmitter signal intensity and the real earthquake catalog is shown in Figure 7a. The probabilistic intensity in this plot is color coded according to the color scale at the top as a function of the time to/from an earthquake (abscissa) and distance between the earthquake epicenter and the NWC-DEMETER great circle path (ordinate). Although the values may in principle range between -0.5 and 0.5, the real obtained range is somewhat lower due to the averaging. It can be seen that the behavior is rather noisy, in particular at lower distances where the amount of data in individual bins is quite low. However, focusing on effects lasting longer than one time bin, two apparently most significant effects identified in the figure are: (i) a decrease of the wave intensity at low distances (within about 4 degrees) shortly (0-3 hours) after the times of the main shocks, and (ii) an increase of the wave intensity at low distances observed about 50-57 hours before the times of the main shocks, followed by a subsequent intensity decrease. Figure 7b uses the same format as Figure 7a to represent the dependence of the normalized probabilistic intensity. It can be seen that although the normalized probabilistic intensity values exhibit significant noise-like fluctuations, in agreement with their expected statistical behavior, the aforementioned intensity variations identified in the probabilistic intensity at low distances remain arguably the most pronounced effects. Figure 7c then shows the number of DEMETER half-orbits contributing to individual distance-time bins. This number is typically rather low, in particular at low distances. Unfortunately, this is a principal limitation given by the DEMETER orbit and transmitter/earthquake locations, and it can be in no way overcome by the performed analysis. Altogether, 715 seismic events contributed to the statistics of the NWC transmitter signal intensity variations.

Considering that the observed intensity variations possibly related to seismic activity occur -as it might be expected -at low distances, we further limit the analysis only to the data obtained at distances less than 4 • and use the same processing method but without an explicit dependence on the distance. This allows us to depict the results using line dependencies, as it is done in Figure 8. The black line in Figure 8a shows the time dependence of the probabilistic intensity, while the black line in Figure 8b shows the time dependence of the normalized probabilistic intensity. It can be seen that the normalized probabilistic intensity is generally within about ±2 standard deviations. Nevertheless, both the intensity decrease just after the main shock times and the intensity variations (first increase, then decrease) about 50 hours before the main shock times are quite unique in the sense that they are not observed in single isolated bins, but they rather span over several consecutive time intervals. This correspondingly increases the statistical significance of the effects. In order to express this quantitatively, one might use 3-hours long time bins in place of the 1-hour long time bins, which would result in the dependencies plotted by the red lines in Figures 8a and8b. The two aforementioned effects are thus the most striking features detected. Their statistical significance is, however, only about 2 standard deviations, i.e., there is about 5% chance that any of them occurred randomly. The number of contributing DEMETER half-orbits depicted in Figure 8c shows that they are formed by some 25 events.

Exactly the same data processing as for the NWC transmitter was performed also for the JJI transmitter signal. The obtained results are depicted in Figure 9. Although the position of the JJI transmitter is favorable for the intended analysis, its output power is lower than the output power of the NWC transmitter. Additionally, the frequency of the JJI transmitter is higher than the Nyquist frequency of the DEMETER VLF measurements. The transmitter signal is thus detected only due to the frequency aliasing, and its intensity is further decreased by the low-pass frequency filters applied onboard. All in all, the number of data relevant for the analysis of a possible seismic influence on the JJI transmitter signal intensity is lower than in the case of the NWC transmitter, as can be seen when comparing Figures 9c-8c. In fact, the number of relevant events is so low that the color representation using the format of Figure 7 is nearly impossible, as it contains many distance-time bins with (nearly) no data. Consequently, only the line representation is used for the JJI transmitter results. The total number of earthquakes contributing to the statistics in the case of the JJI transmitter is 316. Figures 9a and9b show that -although the number of events is low -the decrease of the wave intensity just after the times of the main shocks observed previously for the NWC transmitter is well detectable also in the case of the JJI transmitter. This is a significant supporting piece of evidence that this effect is indeed real. The intensity increase well before the main shock times, followed by a subsequent decrease, which was observed for the NWC transmitter is not apparent in the JJI transmitter results. However, a decrease of the JJI transmitter signal intensity is observed about 48-51 hours before the main shock times, roughly 8, but for the results related to the JJI transmitter signal intensity. Note that the occasional gaps in the black curve in the are due to no data in these particular bins.

corresponding to the times when the decrease in the NWC transmitter signal intensity is observed.

Discussion

The observed decreases of the transmitter signal intensity might be consistent with the effect reported formerly [19] using a different data processing method. A comparable decrease of natural signal intensities at frequencies of about 1.7 kHz, i.e., close to the nighttime cut-off frequency of the Earth-ionosphere waveguide, was also reported [20][21][22][23]. These results suggest that around the times of strong shallow earthquakes, there might be significant changes at the bottom of the ionosphere, which, in turn, affect the wave propagation in the Earth-ionosphere waveguide [START_REF] Budden | The Wave-Guide Mode Theory of Wave Propagation[END_REF]. An activation of highly mobile electronic charge carriers in rocks that are subjected to increasing levels of stress causes positive surface potential and field-ionization of air molecules [6]. The observed wave intensity decrease could be then due to the absorption of VLF electromagnetic waves by the O + 2 molecules in the rising plume of ionized air (out of the air molecules, O 2 has the lowest ionization potential, and, moreover, it has a dipole moment and, hence, can couple to the electromagnetic field). However, a considerable difference between our former and recent results is that while the formerly reported intensity decrease at frequencies of about 1.7 kHz occurred shortly (0-4 hours) before the times of the main shocks, the VLF transmitter signal intensity decrease we reveal occurs after the main shock times. The origin of this difference remains unclear at the moment. However, we note that the analyzed effect is quite different. The formerly reported intensity decrease before main shock times was observed for natural electromagnetic emissions in a wide frequency range around the epicenter location, possibly related to an excess of ionization above the epicenter areas. On the other hand, the intensity decrease after main shock times we report on in the present study is observed for monochromatic artificial VLF transmitter signals at higher frequencies propagating in the Earth-ionosphere waveguide above seismic zones.

The ionospheric variations might be possibly related to a complicated chain of processes termed litosphere-atmosphere-ionosphere coupling, which has been theoretically analyzed in several recent papers [START_REF] Pulinets | Ionospheric Precursors of Earthquakes[END_REF][START_REF] Harrison | [END_REF]7]. We note, however, that the number of relevant events entering the performed statistics is rather small, resulting in not entirely convincing statistical significance of the observed effects. We further note that the intensity decrease observed just after the main shock times both in the NWC and JJI transmitter signal data may be possibly explained directly by seismically generated acoustic-gravity waves influencing the bottom of the ionosphere, and thus, in turn, the wave propagation in the Earth-ionosphere waveguide.

The variations of the transmitter signal intensity observed in 1-hour resolution plots possibly related to the earthquake occurrence are rather weak, and, given the limited number of relevant events, statistically significant only at less than 2 standard deviations level. Furthermore, comparably strong intensity variations were found also for other times and distances. Nevertheless, these other intensity variations seem unlikely to be related to the seismic activity, either because they only appear at large distances, or because the intensity in the neighboring time intervals is significantly different. On the other hand, the reported intensity variations possibly related to the seismic activity seem to be limited to low distances of the earthquake epicenters from the transmitter-DEMETER great circle path, and they span over several consecutive time bins. Taking this into account, and considering 3-hours resolution plots, the statistical significance of the observed phenomena is about 2 standard deviations. This means that there is about 5% chance that such intensity variations might have occurred randomly. Additionally, the intensity decrease after the main shocks was revealed independently for two different VLF transmitters, which is an important additional piece of evidence that the effect is indeed real and related to the earthquake occurrence.

Unfortunately, the available data set is rather limited, and it does not allow us to obtain more firm conclusions. While DEMETER acquired a huge amount of data during its mission, our analysis is limited only to the data from a relatively small area close to the transmitters. Given the geometry of the situation, the amount of the usable data is on the order of 10 half-orbits in each of the bins. At the lowest distances between the earthquake epicenters and the transmitter-DEMETER great circle paths, which are arguably most important for our analysis, the number of contributing half-orbits is even lower. Additionally, the amount of data further decreases for earthquakes with larger magnitudes, which makes it very hard to investigate whether the observed effect would be stronger for stronger earthquakes or not. For the same reason, no attempt to divide the earthquakes according to their type and/or hypocenter location below landmass/shoreline/ocean has been done.

The observed decreases of the transmitter signal intensity appear to be consistent with formerly reported transmitter signal variations related to the seismic activity, which also typically revealed some sort of an intensity decrease (e.g., [19,38,41-44, 46,47]). At least for the intensity decrease observed after the main shock times, one may likely readily adopt the explanation based on the influence of atmospheric gravity waves propagating from an earthquake region [19]. On the other hand, the increase of the wave intensity well before the main shock times observed in the NWC transmitter data might be due to a change of the lower ionosphere resulting in a lower attenuation of signals propagating in the Earth-ionosphere waveguide. Such a change could eventually be a result of the litosphere-atmosphere-ionosphere coupling processes (see, e.g., [START_REF] Pulinets | Ionospheric Precursors of Earthquakes[END_REF][START_REF] Harrison | [END_REF]7] and references therein). However, given that the intensity increase is not reproduced in the JJI transmitter results, a possibility that it might be only a statistical fluctuation cannot be excluded.

The identified variations of the VLF transmitter signal intensity are statistically significant only at about 2 standard deviations level, i.e., there is about 5% change that they might have occurred randomly. Moreover, they were obtained by combining together the data measured during the entire DEMETER mission. Considering individual events, the intensity variation may be clearly rather different than the average behavior obtained in the present study. Although our results thus indicate that there might be a relation between the seismic activity and VLF wave intensity measured by satellites in the upper ionosphere, their application to individual events or even to a short-term earthquake prediction does not appear to be feasible at the moment. We believe that the use of the DEMETER data was in this sense pushed nearly to its boundaries. In order to provide a definite conclusion on the existence of seismic related effects on the wave intensity in the upper ionosphere, it would strongly help to have another satellite similar to DEMETER (the more similar to DEMETER, the better -and the more satellites, the better). This would allow us to perform the same analyses using independent data set(s). Only the effects reproduced by these independent analyses should be then considered as real. In fact, one such satellite, China Seismo-Electromagnetic Satellite (CSES), was recently (February 2018) launched [58], and we are looking forward to its data set and results.

Conclusions

We have analyzed a possible influence of the seismic activity on the intensity of VLF transmitter signals measured by the low-orbiting spacecraft. Due to the transmitter signal intensity issues, only the data measured during the nighttime were included in the analysis. Moreover, following former statistical results, we focused exclusively on earthquakes with magnitudes larger than or equal to 5.0 and depths of hypocenters shallower than 40 km. We have used primarily the signal of the NWC transmitter, as it is the most powerful transmitter advantageously located close to seismically active regions. We applied a two-step data processing which allows us to account for intensity variations related to non-seismic factors, and to evaluate the statistical significance of the observed effects. The obtained statistical results suggest that there might be slight (about 2 standard deviations) intensity variations related to the earthquake occurrence observed at the times when the great circle path between the NWC transmitter and the DEMETER spacecraft is within about 4 degrees from the earthquake epicenter. Specifically, an intensity decrease shortly (0-3 hours) after the times of the main shocks significant at about 2 standard deviations level was observed and independently reconfirmed using the JJI transmitter data, although the number of relevant events was very low. Additionally, some NWC transmitter signal intensity variations were observed also before the times of the main shocks. These were, however, only partially reproduced by the JJI transmitter results. The relatively low number of events included in the analysis, along with rather weak magnitudes of the observed effects, are insufficient to unambiguously demonstrate the statistical significance of the obtained results. The observed intensity decrease just after the main shock times is consistent with acoustic-gravity waves propagating from the earthquake region and influencing the bottom of the ionosphere. The transmitter signal intensity variations observed before the main shock times might be a manifestation of pre-seismic litosphere-atmosphere-ionosphere coupling processes. However, they still remain to be confirmed and eventually fully explained.

Fig. 1 .

 1 Fig.1. (a) A scheme of the situation used to evaluate the earthquake distance. The location of a transmitter is shown by the black triangle. The location of the earthquake epicenter is shown by the red cross. The black curve shows the projection of the DEMETER orbit. The great circle path between the transmitter location and the spacecraft location at a given time is shown by the blue curve. The distance of the particular measurement from the earthquake considered for the performed analysis is then defined as the shortest distance between the earthquake epicenter and the great circle path. (b) A sketch of the transmitter signal propagation along the magnetic field line to the conjugate hemisphere. The transmitter location is shown by the black triangle at the top. Its signal propagates in the Earth-ionosphere waveguide (dark red arrows) and eventually leaks through the ionosphere to higher altitudes (dark green arrows) where it is measured by the DEMETER spacecraft. It further propagates roughly along the magnetic field line to the opposite hemisphere. The signal observed in the conjugate hemisphere may be significantly modified due to wave-particle interactions taking place in the equatorial region.

Fig. 2 .

 2 Fig.2. Number of DEMETER half-orbits with data potentially influenced by seismic effects, assuming that these effects occur exclusively during the night, last for 1 hour around the time of the main shock, and they affect the transmitter signal propagation at great circle paths within 3 degrees. The number of relevant half-orbits is color coded according to the scale on the right-hand side as a function of the considered transmitter (abscissa) and maximum allowed distance between DEMETER and the transmitter (ordinate).

Fig. 3 .

 3 Fig. 3. Geographic locations of earthquakes with magnitudes larger than or equal to 5 and depths less than or equal to 40 km that occurred in the region of interest during the duration of the DEMETER mission. The magnitude of individual earthquakes is color coded according to the color scale on the right-hand side. The black curve shows the geomagnetic equator. The geographic locations of the NWC and JJI transmitters are shown by the black triangle and black diamond, respectively.

Figure 4 .Fig. 4 .Fig. 5 .

 445 Fig. 4. Median power spectral density of electric field fluctuations measured by DEMETER at frequencies around 19.8 kHz, corresponding to the NWC transmitter frequency. (a) Data measured during the daytime half-orbits. (b) Data measured during the nighttime halforbits. The spatial resolution is 2 • both in longitude and latitude. The black triangle marks the location of the NWC transmitter. The black curve corresponds to the geomagnetic equator.

Fig. 6 .

 6 Fig. 6. (a) Normalized probabilistic intensities obtained for the NWC transmitter and a surrogate earthquake catalog are color coded as a function of the time relative to the times of the main shocks (abscissa) and distance of the NWC transmitter-DEMETER great circle paths from earthquake epicenters (ordinate). (b) Distribution of values of normalized probabilistic intensity obtained for surrogate earthquake catalogs (obtained by combining 24 plots similar to the one in panel a) together). Overplotted by the red dashed curve is a Gaussian distribution with the mean value equal to 0 and the standard deviation equal to 1, demonstrating that the values of normalized probabilistic intensity follow this distribution.

Fig. 7 .

 7 Fig. 7. Color coded as a function of the time relative to the times of the main shocks (abscissa) and distance of the NWC transmitter-DEMETER great circle paths from earthquake epicenters (ordinate). (a) Probabilistic intensity. (b) Normalized probabilistic intensity. (c) Number of DEMETER half-orbits contributing to the statistics.

Fig. 8 .

 8 Fig. 8. Line representation of the results from Figure 7 at distances between the NWC transmitter-DEMETER great circle paths and earthquake epicenters lower than 4 degrees: (a) Probabilistic intensity. (b) Normalized probabilistic intensity. (c) Number of DEMETER half-orbits contributing to the statistics. The time resolution of the dependencies plotted in black is 1 hour, the time resolution of the dependencies plotted in red is 3 hours.

Fig. 9 .

 9 Fig.9. Same as Figure8, but for the results related to the JJI transmitter signal intensity. Note that the occasional gaps in the black curve in the are due to no data in these particular bins.

  • × 2 • latitude-longitude bins. Seven frequency bins of the DEMETER ICE instrument (19.53 Hz each) centered at the 19.8 kHz frequency of the NWC transmitter were used for the analysis in

DEMETER was operated by CNES and the corresponding staff are deeply acknowledged. We also thank J. J. Berthelier, who was the principal investigator of the electric field instrument. The data used in the study are available at the Centre de Données de la Physique des Plasmas (http://cdpp2.cnes.fr/cdpp). The earthquake data used in the analysis can be downloaded from http://www.ldeo.columbia.edu/ ~gcmt/projects/CMT/catalog/jan76_ dec13.ndk.

Publisher's Note The EPJ Publishers remain neutral with regard to jurisdictional claims in published maps and institutional affiliations.


Mickaël Binois
email: mickael.binois@inria.fr

Abderrahmane Habbal
email: habbal@unice.fr

Victor Picheny
email: victor@secondmind.ai

A game theoretic perspective on Bayesian multi-objective optimization

This chapter addresses the question of how to efficiently solve many-objective optimization problems in a computationally demanding black-box simulation context. We shall motivate the question by applications in machine learning and engineering, and discuss specific harsh challenges in using classical Pareto approaches when the number of objectives is four or more. Then, we review solutions combining approaches from Bayesian optimization, e.g., with Gaussian processes, and concepts from game theory like Nash equilibria, Kalai-Smorodinsky solutions and detail extensions like Nash-Kalai-Smorodinsky solutions. We finally introduce the corresponding algorithms and provide some illustrating results.

Introduction

In the previous Chapter 11, methods for multi-objective (MO) Bayesian optimization (BO) have been introduced. Scaling with many-objectives burst the seams of these existing solutions on several aspects: i) The positive-definiteness requirement for multi-output Gaussian process (GP) covariance matrix severely restricts their modelling ability, letting independent GPs as the default alternative that does not directly exploit possible correlations between objectives (see also Chapter 10); ii) Summing the objectives is a go-to strategy often observed in practice. But picking weights in scalarization methods becomes more arbitrary, with effects that are difficult to apprehend. More generally, specifying preferences for a large number of conflicting objectives can lead to situations where the connection between achievable solutions and the given preferences is not easy to interpret; iii) For Pareto-based infill criteria, the computational cost grows quickly, for instance with hypervolume computations. Plus there are less closed form expressions available than for two or three objectives; iv) The effect of more objectives on the acquisition function landscape is to be studied, but they are expected to be more multimodal due to larger Pareto sets. And these are only a few additional technical problems faced by many-objective optimization (MaOO).

But even the desired outcome is problematic. As the number of objectives increases, the dimension of the Pareto front generally grows accordingly, and representing this latter is more complex. With a limited budget of evaluations, this may not even be achievable: a few hundred design points would not give an accurate representation of a complex ten-dimensional manifold. In addition, even if representing the entire Pareto front was achievable, the decision maker would be left with a very large set of incomparable solutions, with limited practical relevance.

Hence, MaOO is more of an elicitation problem, for which good and principled ways to choose solutions are needed. When the decision maker is involved in the optimization process by giving feedback on pairs of solutions, methods learning this preference information such as [START_REF] John R Lepird | Bayesian preference elicitation for multiobjective engineering design optimization[END_REF][START_REF] Astudillo | Multi-attribute Bayesian optimization with interactive preference learning[END_REF] are available. More diffuse preference information can also be included via the choice of the reference point for hypervolume computation (nullifying the contributions of equivalent or dominated solutions), interactively modifying ranges for objectives [START_REF] Hakanen | On using decision maker preferences with ParEGO[END_REF] or via an aspiration point serving as target (see, e.g., [START_REF] Gaudrie | Targeting solutions in Bayesian multi-objective optimization: sequential and batch versions[END_REF]). Ordering the importance of objectives is another option, exploited e.g., in [START_REF] Majid Abdolshah | Multi-objective Bayesian optimisation with preferences over objectives[END_REF]. With many objectives, possibly of various natures and scales, the burden on the decision maker is heavier, and the effect on the outcome of these preferences is less intuitive.

When the decision is only made a posteriori, knee points, Pareto solutions from which improving slightly one objective would deteriorate much more some others, are appealing to decision makers, or in the design of evolutionary algorithms, see e.g., [START_REF] Zhang | A knee point-driven evolutionary algorithm for many-objective optimization[END_REF]. But their existence and location on the Pareto front may vary, especially in the MaOO case. Other arguments have been directed toward finding a single solution at the "center" of the Pareto front to keep balance between objectives, e.g., by [START_REF] Gaudrie | Budgeted multi-objective optimization with a focus on the central part of the Pareto front-extended version[END_REF][START_REF] Paria | A flexible framework for multi-objective Bayesian optimization using random scalarizations[END_REF]. The rationale is that trade-off solutions at the extremes are not desirable unless prior information has been given. Finding such balanced solutions, without input from the decision maker during the optimization, is the default framework considered in this chapter. The difficulty lies in an appropriate and principled definition of a central solution. A possible approach is to frame the MaOO problem as a game theory one, where each objective is impersonated by a player, and the players bargain until they agree on a satisfying solution for all, referred to as equilibrium.

Game theory is already popular in a variety of applications such as engineering [START_REF] Désidéri | Cooperation and competition in multidisciplinary optimization[END_REF], control [START_REF] Habbal | Neumann-Dirichlet Nash strategies for the solution of elliptic Cauchy problems[END_REF], multi-agent systems [START_REF] Brown | Hierarchical abstraction, distributed equilibrium computation, and post-processing, with application to a champion no-limit Texas hold'em agent[END_REF] or machine learning [START_REF] Lanctot | No-regret learning in extensive-form games with imperfect recall[END_REF]. Game theory has also revealed to be a powerful tool in multi-disciplinary optimization, where decentralized strategies for noncooperative agents takes its full and relevant meaning, as referring to discipline-specific solvers in a parallel, asynchronous and heterogeneous multi-platform setting, see e.g., the pioneering framework developed in [START_REF] Periaux | Evolutionary optimization and game strategies for advanced multi-disciplinary design: applications to aeronautics and UAV design[END_REF]. It comes equipped with principled ways of eliciting solutions (the bargaining axioms of scale invariance, symmetry, efficiency, and independence of irrelevant alternatives), echoing several desirable properties for improvement functions used in MO BO [START_REF] Svenson | Computer Experiments: Multiobjective Optimization and Sensitivity Analysis[END_REF]: reflecting Pareto dominance, no input in the form of external parameters (such as a reference point), or robustness to scale invariance. Indeed, with more objectives of different nature, finding a reasonable common scale for all is much harder and the effects cannot be apprehended.

Compared to the more standard multi-objective BO framework described in the previous chapter, the games framework comes with some specific challenges. The definition of solutions such as Nash equilibrium or Kalai-Smorodinsky solutions can be quite involved, and finding such equilibria generally implies solving several inference tasks at once. For instance, this makes the use of improvement-based approaches difficult. Still, several approaches have emerged using regret or stepwise uncertainty reduction (SUR) frameworks, see e.g., [START_REF] Picheny | A Bayesian optimization approach to find Nash equilibria[END_REF][START_REF] Majid Abdolshah | Multi-objective Bayesian optimisation with preferences over objectives[END_REF] and references therein.

The structure of the chapter is as follows. In Section 2, we briefly review key concepts in game theory such as Nash equilibria. Then in Section 3 we present high-level BO approaches to solve many-objective problems under the Games paradigm. Illustrations are provided in Section 4 before discussing remaining challenges and perspectives in Section 5.

Game equilibria to solution elicitation

The standard multi-objective optimization problem (MOP) corresponds to the simultaneous optimization of all objectives:

(M OP ) min x∈X (f 1 (x), . . . , f m (x)). (1)
Besides considering objectives as individual players' goals, general non-cooperative games, and related solutions, such as Nash ones, do need territory splitting. That is, partitioning of the optimization variables (the input space) among players. Other games and related solution concepts depend on choices of anchor points, as is the case for the Kalai-Smorodinsky (KS) solution, which depends on ideal and disagreement points, as illustrated in Figure 1. We shall introduce in the following the Nash equilibrium concept, then motivated by the generic inefficiency of the latter, move to considering the Kalai-Smorodinsky solution.

Nash games and equilibria

When considering primarily the standard (static, under complete information) Nash equilibrium problem [START_REF] Gibbons | Game Theory for Applied Economists[END_REF], each objective becomes a player's outcome. Compared to a standard MOP [START_REF] Majid Abdolshah | Multi-objective Bayesian optimisation with preferences over objectives[END_REF] where players share the control of the same set of variables (or action space), variables are uniquely allocated to a player in the so-called territory splitting. Denote x i the variables of player i and X i its corresponding action space, where X = i X i . Accordingly, the variable vector x consists of block components x 1 , . . . , x m (x = (x j ) 1≤j≤m ). We shall use the convention f i (x) = f i (x i , x -i ) when we need to emphasize the role of x i , where x -i is the subset of variables controlled by players j = i.

Definition 2.1. A Nash equilibrium problem (NEP) consists of m ≥ 2 decision makers (i.e., players), where each player i ∈ {1, . . . , m} wants to solve its optimization problem:

(P i ) min x i ∈X i f i (x), (2)
and f (x) = [f 1 (x), . . . , f m (x)] : X ⊂ R n → R m (with n ≥ m) denotes a vector of cost functions (a.k.a. pay-off or utility functions when maximized), f i denotes the specific cost function of player i.

Definition 2.2. A Nash equilibrium x * ∈ X is a strategy such that:

(N E) ∀i, 1 ≤ i ≤ m, x * i = arg min x i ∈X i f i (x i , x * -i ). (3)
In other words, when all players have chosen to play a NE, then no single player has incentive to move from his x * i .

The above definition shows some attractive features of NE: it is scale-invariant, does not depend on arbitrary parameters (such as the reference point), and convey a notion of firstorder stationarity. On the other hand, it requires partitioning the inputs between players, which may come naturally for some problems, but may be problematic for many MaOO problems. Plus an arbitrary partitioning can have an unknown outcome on the solution. One option proposed in [START_REF] Désidéri | Multiobjective Design Optimization Using Nash Games[END_REF] is to define the partitioning based on sensitivity analysis of one main objective. Hence a perspective would be to define an optimal partitioning in the sense that it allocates variables (or a linear combination of variables) to the player it has most influence on.

Besides the need to define a territory splitting, it is important to notice that, generically, Nash equilibria are not efficient, i.e., do not belong to the underlying set of best compromise solutions, the so-called Pareto front, of the objective vector (f i (x)) x∈X . Indeed, NE efficiency may happen when one of the players has control on all the optimization variables while the others have control on nothing (a somehow degenerate territory splitting). Still, NE have a couple of advantages: a nice notion of stationarity, they are generally well balanced and, of interest for many objectives, scale invariant.

Remark. When block components x 1 , . . . , x m of the variable vector x overlap due to, e.g., constraints, then one should modify (2) and (3) accordingly, so that x i ∈ X i (x -i ). In the literature, the Nash equilibrium problem is then referred to as a generalized Nash equilibrium problem (GNEP).

Remark. The NE solution is scale invariant, and more generally invariant under any strictly increasing transformation

Ψ i (1 ≤ i ≤ m) since, in this case, (3) is equivalent to ∀i, 1 ≤ i ≤ m, Ψ i (f i (x * i , x * -i )) ≤ Ψ i (f i (x i , x * -i ))
Interestingly, it could be shown at least in the convex case that if players share probabilistic control on the optimization variable components (i.e., any player j may control the same component x i with probability p ij ) then there exists a probability matrix (p ij ) such that the associated Nash equilibrium (with an ad hoc definition) lies in the Pareto front, see [START_REF] Aboulaich | A new algorithm for approaching Nash equilibrium and Kalai Smoridinsky solution[END_REF] for a sketch of the approach.

From a practical viewpoint, computing NE for games stated in continuous variable settings (as opposed to discrete or finite games, e.g., in vector spaces with Banach or Hilbert structure) can be based on variational analysis, e.g., with the classical fixed-point algorithms to solve NEPs [START_REF] Uryas | On relaxation algorithms in computation of noncooperative equilibria[END_REF][START_REF] Başar | Relaxation techniques and asynchronous algorithms for on-line computation of noncooperative equilibria[END_REF][START_REF] Li | Distributed algorithms for the computation of noncooperative equilibria[END_REF]. A modified notion of Karush-Kuhn-Tucker (KKT) points, adapted to (generalized) NE, is developed by [START_REF] Kanzow | Augmented Lagrangian methods for the solution of generalized Nash equilibrium problems[END_REF], to propose a dedicated augmented Lagrangian method. Yet these methods require too many evaluations to tackle directly expensive black-boxes, requiring specific BO algorithms that are presented in Section 3.

The Kalai-Smorodinsky solution

The Kalai-Smorodinsky solution was first proposed by Kalai and Smorodinsky in 1975 as an alternative to the Nash bargaining solution in cooperative bargaining. Differently from Nash equilibria, no partitioning of the decision variable among players is required by the KS solution concept.

The problem is as follows: starting from a disagreement or status quo point d in the objective space, the players aim at maximizing their own benefit while moving from d toward the Pareto front (i.e., the efficiency set). The KS solution is of egalitarian inspiration [START_REF] John | The bargaining problem without convexity: Extending the egalitarian and Kalai-Smorodinsky solutions[END_REF] and states that the selected efficient solution should yield equal benefit ratio to all the players. Indeed, given the utopia (or ideal, or shadow) point u ∈ R m defined by

u i = min x∈X f i (x),
selecting any compromise solution y = [f 1 (x), . . . , f m (x)] would yield, for objective i, a benefit ratio

r i (x) = d i -y i d i -u i .
Notice that the benefit from staying at d is zero, while it is maximal for the generically unfeasible point u. The KS solution is the Pareto optimal choice y * = [f 1 (x * ), . . . , f m (x * )] for which all the benefit ratios r i (x) are equal:

x * ∈ {x ∈ X s.t. r 1 (x) = • • • = r m (x)} ∩ {x ∈ X s.t. f (x) ∈ P} .
Geometrically, y * is the intersection point of the Pareto front and the line (d, u) (see Figure 1), which may be empty. For m = 2, y * exists if the Pareto front is continuous. For the general case, we use here the extension of the KS solution proposed by [START_REF] Leth | Nonconvex n-person bargaining: efficient maxmin solutions[END_REF] under the name efficient maxmin solution. Since the intersection with the (d, u) line might not be feasible, there is a necessary trade-off between Pareto optimality and centrality. The efficient maxmin solution is defined as the Pareto-optimal solution that maximizes the smallest benefit ratio among players, that is:

x * * ∈ arg max y∈P min 1≤i≤m r i (x).
It is straightforward that when the intersection is feasible, then y * and y * * := f (x * * ) coincide.

Figure 1 shows y * * in the situation when the feasible space is nonconvex. y * * is always on the Pareto front (hence not necessarily on the (d, u) line). In the following, we refer indifferently to y * (if it exists) and y * * as the KS solution. Note that this definition also extends to discrete Pareto sets. For m = 2, the (non-extended) KS solution can be axiomatically characterized as the unique solution fulfilling all the bargaining solution axioms, which are: Pareto optimality, symmetry, affine invariance, and restricted monotonicity [START_REF] Kalai | Other solutions to Nash's bargaining problem[END_REF]. For m ≥ 3, there is no axiomatic setting, the KS solution being required to fulfill Pareto optimality, affine invariance, and equity in benefit ratio. Still, KS is particularly attractive in a many-objective context since it scales naturally to a large number of objectives and returns a single solution, avoiding the difficulty of exploring and approximating large m-dimensional Pareto fronts-especially with a limited number of observations.

Closely related to the Kalai and Smorodinsky solution is the reference point methodology, developed by [START_REF] Wierzbicki | The use of reference objectives in multiobjective optimization[END_REF]. The reference point methodology uses achievement functions, which refer to aspiration and reservation references. Ideal and nadir points are particular instances of such reference points, respectively, and achievement functions play the same role as the benefit ratio of our present setting. Moreover, the notion of neutral compromise and max-min (of achievement functions) were introduced by Wierzbicki in the cited reference, yielding a framework objectively close to the one of Kalai and Smorodinsky. One important difference between the reference point and KS approaches is that the latter relies on a game theoretic axiomatic construction, in view of palliating the Nash equilibrium inefficiency.

Disagreement point choice

Clearly, the KS solution strongly depends on the choice of the disagreement point d. A standard choice is the nadir point N given by N i = max x∈Pareto set f (i) (x). Some authors introduced alternative definitions, called extended KS, to alleviate the critical dependence on the choice of d [START_REF] Bozbay | Bargaining with endogenous disagreement: The extended Kalai-Smorodinsky solution[END_REF], for instance by taking as disagreement point the Nash equilibrium arising from a previously played non-cooperative game. But such a choice would need a pre-bargaining split of the decision variable vector x among the m players. When some relevant territory splitting is agreed on, then Nash yields KS an interesting and non arbitrary disagreement point, and in return, KS yields the Nash game an interesting and non arbitrary efficient solution, starting from the non efficient equilibrium. The resulting Nash-Kalai-Smorodinsky solution is denoted NKS.

The nadir point is very useful to rescale the Pareto front in case objectives are of different nature or scale, see e.g., [START_REF] Bechikh | Estimating nadir point in multiobjective optimization using mobile reference points[END_REF]. As such, it makes a natural disagreement solution. Still, finding the nadir point is a complex task, especially under a limited budget, as it involves exploring regions good for a few objectives with worse values on others. The many objective context makes this issue more prominent. It is related to the problem of Pareto resistance, i.e., regions of the space with good values on a few objectives but not Pareto optimal (see, e.g., [START_REF] Fieldsend | Enabling dominance resistance in visualisable distance-based many-objective problems[END_REF] and references therein). Unless some objectives are more important than others, the corresponding extremal regions of the Pareto front are of little interest for selecting a single good solution.

A simpler disagreement point is via the pseudo-nadir point, defined as the worst objective solution over designs achieving the minimum for one objective: Ñi = max x∈{x (1) * ,...,x (m) * } f i (x), where x (j) * ∈ arg min x∈X f j (x). It can be extracted directly from the pay-off table, but it is not robust as it can under-or over-estimate the true nadir point, see e.g., [START_REF] Miettinen | Nonlinear multiobjective optimization[END_REF] for details. While for two objectives, nadir and pseudo-nadir coincide, this is not the case anymore in higher dimensions. Using this pseudo nadir point removes the need to search for extremal non-dominated solutions on the Pareto front, and aligns with finding the ideal point.

The egalitarian inspiration of the KS solution is based on the assumption that all the objectives have equal importance. However, in many practical cases, the end user may want to favor a subset of objectives which are of primary importance, while still retaining the others for optimization. One way of incorporating those preferences [START_REF] Thiele | A preferencebased evolutionary algorithm for multi-objective optimization[END_REF][START_REF] Wang | A mini-review on preference modeling and articulation in multi-objective optimization: current status and challenges[END_REF][START_REF] Robin C Purshouse | A review of hybrid evolutionary multiple criteria decision making methods[END_REF] is to discard solutions with extreme values and actually solve the following problem (constrained version of the original MOP (1)):

min x∈X {f 1 (x), . . . , f m (x)} s.t. f i (x) ≤ c i , i ∈ J ⊂ [1, . . . , m] ,
with c i 's predefined (or interactively defined) constants. Choosing a tight value (i.e., difficult to attain) for c i may discard a large portion of the Pareto set and favor the i-th objective over the others. Incorporating such preferences in the KS solution can simply be done by using c as the disagreement point if all objectives are constrained, or by replacing the coordinates of the nadir with the c i values. In a game theory context, this would mean that each player would state a limit of acceptance for his objective before starting the cooperative bargaining. From geometrical considerations, one may note that d (and hence c) does not need to be a feasible point. As long as d is dominated by the utopia point u, the KS solution remains at the intersection of the Pareto front and the (d, u) line. An example where preferences from a decision maker are integrated this way is provided in [START_REF] Binois | The Kalai-Smorodinsky solution for many-objective Bayesian optimization[END_REF].

As shown for instance in [START_REF] Binois | The Kalai-Smorodinsky solution for many-objective Bayesian optimization[END_REF], by applying a monotonic rescaling to the objectives, the corresponding KS solution can move along the Pareto front. An example of a monotonic value function applied to a raw objective corresponds to taking the logarithm of one objective for a physical or modeling reason. Since Nash equilibrium are invariant to such rescaling, they are more stable than nadir or pseudo-nadir based KS solutions, but still affected. One proposed option to alleviate this effect is to use ranks instead of objective values directly. Interestingly, it induces a dependency on the parameterization of the problem. Studying these points can include the use of a measure on the input space and the use of copulas in the objective one, [START_REF] Binois | The Kalai-Smorodinsky solution for many-objective Bayesian optimization[END_REF][START_REF] Binois | On the estimation of Pareto fronts from the point of view of copula theory[END_REF], with further work needed to quantify the effect on the Pareto front. The combination with more involved preference learning methods such as [START_REF] John R Lepird | Bayesian preference elicitation for multiobjective engineering design optimization[END_REF][START_REF] Astudillo | Multi-attribute Bayesian optimization with interactive preference learning[END_REF] could also be investigated.

Bayesian optimization algorithms for games

Solving expensive black-box games is an emerging topic in Bayesian optimization [START_REF] Aprem | A Bayesian optimization approach to compute Nash equilibrium of potential games using bandit feedback[END_REF][START_REF] Al-Dujaili | Approximating Nash equilibria for black-box games: A Bayesian optimization approach[END_REF][START_REF] Picheny | A Bayesian optimization approach to find Nash equilibria[END_REF][START_REF] Binois | The Kalai-Smorodinsky solution for many-objective Bayesian optimization[END_REF]. Similarly to the single and multi-objective cases, BO leverages the probabilistic information available from the Gaussian process models to balance exploration and exploitation in searching for the game equilibrium. Other metamodelling options could be entertained, as long as they also provide mean and variance predictions at any location, plus the ability to get realizations (posterior draws) at arbitrary locations. We present here several algorithms to solve such equilibrium or solution finding problems.

For the Nash equilibrium (Section 3.1), [START_REF] Aprem | A Bayesian optimization approach to compute Nash equilibrium of potential games using bandit feedback[END_REF] proposed an upper confidence bound-like acquisition function based on a regret function approximating the game-theoretic regret. The case of potential games [START_REF] González | A survey of static and dynamic potential games[END_REF], when a single function can summarize the game, commonly arising in problems with shared resources, is considered by [START_REF] Al-Dujaili | Approximating Nash equilibria for black-box games: A Bayesian optimization approach[END_REF], with a dedicated expected improvementlike acquisition function. In [START_REF] Picheny | A Bayesian optimization approach to find Nash equilibria[END_REF][START_REF] Binois | The Kalai-Smorodinsky solution for many-objective Bayesian optimization[END_REF], the flexible and general step-wise uncertainty reduction is applied to estimate different classes of equilibria, including Nash and KS, which we present in Section 3.2. Although not proposed in the literature (up to our knowledge), a related family of acquisition functions based on Thompson sampling is also applicable, as we describe briefly in Section 3.3.

Fixed point approaches for the Nash equilibrium

Fixed points methods, see an example in Algorithm 8, are one way to get Nash equilibrium. In the expensive case, a naïve version is simply to replace the true objectives by their GP mean prediction: fi , 1 ≤ i ≤ m. While computationally efficient and possibly sufficient for a coarse estimation, such an approach does not address the exploration-exploitation trade-off, and may lead to poor estimates if the models are not accurate representations of the objectives.

[3] provides an acquisition function based on the upper confidence bound to estimate Nash equilibrium for continuous games. They focus on the game theoretic regret, that is the most any player can gain from deviating from the current solution. Following small deviations for each player in turns is also the underlying principle in fixed point methods. The corresponding acquisition function aims to select as next design to evaluate one achieving the minimum of the regret of the GP surrogate. But rather than optimizing directly to get the minimum objective values, these are replaced by an approximation. That is, for a player i at x i , the minimum value is replaced by the sum of the mean objective value over x -i with the corresponding standard deviation, scaled by a factor γ. These can be computed analytically for GPs for some separable covariance functions.

In the special case of potential games [START_REF] González | A survey of static and dynamic potential games[END_REF], also related to fixed point methods, the strategy of improving as much as possible one player leads to the Nash equilibrium (and not any fixed point). This is exploited by [START_REF] Aprem | A Bayesian optimization approach to compute Nash equilibrium of potential games using bandit feedback[END_REF] who directly considers the potential. By definition, the potential Φ :

X → R is a function such that f i (x , x -i ) -f i (x , x -i ) = Φ(x , x -i ) - Φ(x , x -i ), ∀i, x , x .
As the available measurements are on the f i 's and not the potential directly, constructing the GP of the potential Φ involves taking into account integral and gradient operators. Once the GP is fitted, they proposed an acquisition function similar to expected improvement [START_REF] Mockus | On bayesian methods for seeking the extremum[END_REF].

Stepwise uncertainty reduction

Step-wise uncertainty reduction (SUR) methods, see, e.g., [START_REF] Geman | An active testing model for tracking roads in satellite images[END_REF][START_REF] Villemonteix | An informational approach to the global optimization of expensive-to-evaluate functions[END_REF][START_REF] Chevalier | Fast parallel kriging-based stepwise uncertainty reduction with application to the identification of an excursion set[END_REF], require an uncertainty measure Γ about a quantity of interest, here the game equilibrium. Then new evaluations are selected sequentially in order to reduce this uncertainty, leading to an accurate identification of the quantity of interest. This framework is well-adapted to games as it allows to integrate various learning tasks into a single goal (for instance, the definition of KS type solutions involves several other unknown quantities, such as the ideal and nadir points).

For any multivariate function f defined over X , with image Y, denote by Ψ : Y → R m the mapping that associates its equilibrium (Nash, KS or other). When the multivariate function is replaced by a GP emulator conditioned on n observations f (x (1) ), . . . , f (x (n) ) , Ψ is defined on the corresponding distribution Y n (). Then Ψ(Y n ) is a random vector (of unknown distribution). Loosely speaking, the spread of the corresponding distribution in the objective space reflects the uncertainty on the solution location, given by the uncertainty measure Γ(Y n ). One measure of variability of a vector is the determinant of its covariance matrix [START_REF] Vadimovich | Theory of Optimal Experiments[END_REF]:

Γ(Y n ) = det [cov (Ψ(Y n ))],
while other information theoretic measures could be used. For instance, this could be the conditional entropy as in [START_REF] Villemonteix | An informational approach to the global optimization of expensive-to-evaluate functions[END_REF].

The SUR strategy greedily chooses the next observation that reduces the most this uncertainty:

max x∈X {Γ(Y n ) -Γ(Y n,x )} ,
where Y n,x is the GP conditioned on {f (x (1) ), . . . , f (x (n) ), f (x)}. However, such an ideal strategy is not tractable as it would require evaluating all f (x) while maximizing over X .

A more manageable strategy is to consider the expected uncertainty reduction:

Γ(Y n ) - E Yn(x) [Γ(Y n,x )]
, where E Yn(x) denotes the expectation taken over Y n (x). Removing the constant term Γ(Y n ), the policy is defined with

x ∈ arg min x∈X J(x) = E Yn(x) [Γ(Y n,x )] .
Informally speaking, the variation of Ψ is partly caused by not knowing the precise localization of the solution of interest, plus eventually from not knowing the values of the anchor points.

The realizations corresponding to these specific points would either be for designs close to the solution of the GP means, or for designs with large variance. Hence J(x) defines a trade-off between exploration and exploitation, as well as a trade-off between the different learning tasks (d and u points, Nash equilibrium, Pareto front). Hence the BO optimization loop is fully defined (see Chapter 11, Algorithm 1) and the difficulty boils down to evaluating J efficiently.

A crucial aspect in enabling the use of the SUR strategy is the ability to generate realizations or estimate Y n efficiently. When X is discrete (or can be discretized), the approaches proposed in [START_REF] Picheny | A Bayesian optimization approach to find Nash equilibria[END_REF][START_REF] Binois | The Kalai-Smorodinsky solution for many-objective Bayesian optimization[END_REF] relies on the use of conditional simulations, i.e., joint posterior samples on X or a well chosen subset of it, coupled with fast update formulas for the resulting ensembles, following [START_REF] Chevalier | Fast update of conditional simulation ensembles[END_REF]. This approach hinges on the discretization size used, which needs to remain in the thousands due to the cubic computational complexity in the number of samples. The targeted solution is then computed on each realization, resulting in samples from Y n . Note that for discrete Nash equilibrium computations, all combinations of the selected strategies for each player must be simulated at. An alternative solution, applicable to continuous X would rely on approximated sample paths with a closed form expression as used, e.g., by [START_REF] Mutny | Efficient high dimensional Bayesian optimization with additivity and quadrature Fourier features[END_REF]. Again, the appropriate solution can be obtained from the samples, which is inexpensive compared to using the expensive black-box.

The resulting criterion for KSE estimation automatically interleaves improving the estimation of the specific quantities involved, by balancing between estimation of the Nash equilibrium, of the objective-wise minima for the ideal point and of the intersection with the Pareto front. A different use of the posterior distribution is also possible, with Thompson sampling.

Thompson sampling

Over the past decade, Thompson sampling (TS) [START_REF] William R Thompson | On the likelihood that one unknown probability exceeds another in view of the evidence of two samples[END_REF] has become a very popular algorithm to solve sequential optimization problems. In a nutshell, TS proceeds by sequentially sampling from the posterior distribution of the problem solution, allowing for efficiently addressing exploration-exploitation trade-offs. In a GP-based setting, this is simply achieved by sampling from the posterior distribution of the objectives and choosing as the next sampling point the input that realizes the equilibrium on the sample.

As it would not directly try to pinpoint the ideal point (say), it cannot be directly transposed for estimating KS solutions. Still, it may serve to estimate Nash equilibria. However, up to our knowledge this solution has not been studied yet (the closest proposition is an algorithm in [START_REF] Picheny | A Bayesian optimization approach to find Nash equilibria[END_REF] based on the probability of realizing the optimum). Note that the same approximation of the samples as for SUR can be used. Next we illustrate the results of this latter on a practical application.

Application example: engineering test case

We consider the switching ripple suppressor design problem for voltage source inversion in powered systems, which was proposed by [START_REF] Zhang | Switching ripple suppressor design of the grid-connected inverters: A perspective of many-objective optimization with constraints handling[END_REF], and also tested in [START_REF] He | A repository of real-world datasets for data-driven evolutionary multiobjective optimization[END_REF]. The device is composed of three components, one with three inductors, one with parallel LC resonant branches and a capacitor branch. Arguing that this type of problem is naturally many objective, [START_REF] Zhang | Switching ripple suppressor design of the grid-connected inverters: A perspective of many-objective optimization with constraints handling[END_REF] describes several types of conflicting objectives: harmonic attenuations of the switching ripple suppressor under different frequencies, power factor, resonance suppression and inductor cost.

In the proposed problem, where n r is the number of resonant branches, there are n r + 4 variables and n r + 1 objectives: variables are inductors L 1 , L 2 , L 3 that can be related to the inductor cost objective (f nr+1 ), and variables C 1 , . . . , C nr related to the goal of attenuating the harmonics at resonant frequencies (f 1 , . . . , f nr ). The last variable is C f . Note that five additional design constraints can also be added. The cost of evaluation is here negligible, but it showcases a realistic application example where the true Pareto front is unknown.

We start by estimating the Nash-Kalai-Smorodinsky equilibrium in the n r = 4 case; with the following partitioning: (f 1 : C 1 , C f ), (f 2 : C 2 ), (f 3 : C 3 ), (f 4 : C 4 ), (f 5 : L 1 , L 2 , L 3 ). For the SUR strategy, we discretize each input subspace using Latin hypercube samples of sizes [START_REF] Hakanen | On using decision maker preferences with ParEGO[END_REF][START_REF] Brown | Hierarchical abstraction, distributed equilibrium computation, and post-processing, with application to a champion no-limit Texas hold'em agent[END_REF][START_REF] Brown | Hierarchical abstraction, distributed equilibrium computation, and post-processing, with application to a champion no-limit Texas hold'em agent[END_REF][START_REF] Brown | Hierarchical abstraction, distributed equilibrium computation, and post-processing, with application to a champion no-limit Texas hold'em agent[END_REF][START_REF] Wierzbicki | The use of reference objectives in multiobjective optimization[END_REF] and use a budget of 80 evaluations for the initial design followed by 70 sequential iterations. Note that this discretization of the input space with the partitioning results in ≈ 2 × 10 6 possible combinations for the Nash equilibrium. To keep it tractable for the SUR procedure, filtering is performed to select 1000 solutions to actually draw posterior realizations, before evaluating the criterion on the 200 most promising candidates. More details on the procedure can be found in [START_REF] Picheny | A Bayesian optimization approach to find Nash equilibria[END_REF][START_REF] Binois | The Kalai-Smorodinsky solution for many-objective Bayesian optimization[END_REF]. The result obtained using the GPGame package [START_REF] Picheny | GPGame: Solving Complex Game Problems using Gaussian Processes[END_REF] is described in Figures 2 and3. It appears that the first four objectives are highly correlated, which does not affect the estimation of the solution. Most of the initial designs (67/80) are dominated, while this is the case for two thirds (47/70) of the added ones, possibly due to the search of the inefficient Nash equilibrium. These new designs are all concentrated in the same region of the output space. The NE is relatively close to the Pareto front in this case, with the obtained solution dominating it.

Then we consider the regular Kalai-Smorodinsky solution, whose results are in Figures 4 and5. This time the input space is discretized with 10 6 possible solutions, uniformly sampled. Again, most initial design are dominated (72/80), but this is the case for only a third (28/70) of the sequentially added ones. Figure 5 shows that the KS solution is almost a straight line, an indication of the centrality of this solution. It can be seen that, in the sequential procedure, new points are added at the extremities of the Pareto front to locate the nadir, as well as in the central part.

Choosing between the obtained NKS and KS solutions now depends on the viewpoint. The spread of the latter is much broader in the objective space than for NKS. Note that the same range is used on all the Figures. To further ease comparison, we added the KS solution on Figure 2 and the Nash and NKS ones on Figure 4, but there were not part of the solutions evaluated in their corresponding BO loops, nor their targets. In case the range of objectives is relevant to the decision maker, the centrality and equity of benefits on all objectives may be appealing to the decision maker. If instead the territory splitting for the objectives is important, the relative scaling of the objectives irrelevant, then the NKS solution would be preferred (also over staying on the Nash equilibrium). As a perspective, it would be pertinent to additionally take into account the constraints defined in [START_REF] Zhang | Switching ripple suppressor design of the grid-connected inverters: A perspective of many-objective optimization with constraints handling[END_REF], for which the Pareto front becomes more complex.

What is done and what remains

We reviewed the enrichment of many-objective optimization with concepts from game theory. In the context of limited budgets of evaluations when the extensive search of the entire Pareto front is unreachable, game theory can guide the estimation of a single solution chosen for its centrality and possible additional game theoretic properties. Several BO-based method- ologies have been developed for estimating the corresponding solution, which can scale to many-objective setups. As this is an emerging topic in BO, it seems that many algorithmic developments could be proposed, for instance based on Thompson sampling or on entropy search.

For the Nash and NKS solutions, a partitioning of the input space is necessary, while it is seldom included in the multi-objective context. A perspective is to randomize this partitioning when it does not come directly from the problem definition. The possible advantage would be to focus on a set of different Nash equilibria, all with stationarity properties, providing more than a single solution.

As with the engineering example we presented, there are usually additional constraints in the problem definition. In the SUR procedure described above, dealing with constraints could amount to filtering unfeasible solutions on the realizations. Other options could include the use of the probability of feasibility as in [START_REF] Schonlau | Global versus local search in constrained optimization of computer models[END_REF]. For more details on the handling of constraints in BO, we refer to [START_REF] Letham | Constrained bayesian optimization with noisy experiments[END_REF] and references therein. Then the difference between an objective and a constraint is sometimes fuzzy. Especially when the constraints can be modified, taking them as objectives instead, akin to multiobjectivization [START_REF] Knowles | Reducing local optima in single-objective problems by multi-objectivization[END_REF], could help finding feasible solutions (or work in the case the constrained problem is unfeasible). Conversely, prioritizing some objectives [START_REF] Desideri | Platform for prioritized multi-objective optimization by metamodel-assisted Nash games[END_REF] while defining constraints on others can help control the trade-offs when moving on the Pareto front.

Figure 1 :

 1 Figure 1: Illustration of solutions of two different elicitation problems with two different disagreement points d: either the Nadir point or a Nash equilibrium, resulting in the KS and NKS points (stars) on the Pareto front P, respectively. The shaded area shows the feasible objective space.

5 NKSFigure 2 :

 52 Figure 2: Scatter plot matrix of the objectives obtained when estimating the NKS solution by the SUR strategy. Black circles: dominated initial designs, black +: nondominated initial designs, green crosses: sequentially added solutions (circles if dominated), blue triangles: estimated Nash equilibrium, and red diamond: estimated NKS solution. The KS solution (yellow square) is marked for reference only.

Figure 3 :

 3 Figure 3: Parallel coordinates plot corresponding to the objectives of Figure 2. Colors are the same, dominated solutions are in dotted lines while the Nash equilibrium is marked by the large blue dashed line and the estimated NKS solution by the large red line.

5 KSFigure 4 :

 54 Figure 4: Scatter plot matrix of the objectives obtained when estimating the KS solution by the SUR strategy. The legend is as for Figure 2 except the yellow square indicates the KS solution.Here the Nash solution (blue diamond) and NKS solution (red triangle) are marked for reference only.

Figure 5 :

 5 Figure 5: Parallel coordinates plot corresponding to the objectives of Figure 4. Colors are the same, dominated solution are depicted with dotted lines while the estimated KS solution is marked by the large yellow line.

Ensure: For all i = 1, ..., m, x * i = arg min x i ∈X i J i (x * -i , x i )

Algorithm 1: Pseudo-code for the fixed-point approach [START_REF] Uryas | On relaxation algorithms in computation of noncooperative equilibria[END_REF] 1 Require: m: number of players, 0 < α < 1 : relaxation factor, n max : max iterations; 2 Construct initial strategy x (0) ; k) small enough then 7 exit;