Nicolas Picot-Clemente

ON BEHALF OF THE Antares Collaboration

Search for magnetic monopoles with the ANTARES underwater neutrino telescope

Keywords: ANTARES, magnetic monopoles, exotic physic, neutrino astronomy Vol. 5, 93

The ANTARES underwater neutrino telescope, located in the Mediterranean Sea at a depth of 2475 m, 40km off the Provencal coast, is composed by an array of 885 photomultipliers distributed on 12 vertical lines. The detector is fully operational since May 2008. Besides the detection of neutrino-induced muons, the telescope is more generally sensitive to particles which emit Cherenkov light, and, thanks to its large volume, offers new opportunities to improve the sensitivity to exotic particles, such as magnetic monopoles. Magnetic monopoles are stable magnetically charged particles first introduced in quantum mechanics by Dirac in 1931, which showed that their existence would give an explanation to the quantization of the electric charge. They would have been produced in the Early Universe, and would bring a first proof to the existence of grand unified models. The paper introduces magnetic monopoles and their signal caracteristics in a neutrino telescope, and then describes the analysis which uses a dedicated algorithm able to reconstruct the velocity of the incoming particles. A new upper limit on the upward going magnetic monopole flux, extracted from ANTARES data taken in 2008, is presented. It is, at present, the best worldwide constraint in the velocity range β = [0.625, 1], and provides, for the first time with a neutrino telescope, a constraint for monopole velocities below the Cherenkov threshold.

Introduction

Magnetic monopoles are hypothetical magnetically charged particles first introduced by Dirac in 1931. He demonstrated that the existence of magnetic monopoles naturally leads to the quantization of the electric charge [1]. In 1974, 't Hooft and Polyakov discovered independently that, in certain spontaneously broken gauge theories, magnetic monopoles are not only a possibility, but a requirement [2,3]. However, despite intensive search efforts, no magnetically charged particles were detected up to now, and stringent experimental flux limits beyond the original Parker bound of flux F P ∼ 10 -15 cm -2 s -1 sr -1 [4] were provided. The development of huge detectors for performing neutrino astronomy gives rise to new expectations for the search for magnetic monopoles, and ANTARES offers new interesting opportunities in this field. In this paper a search for upward going relativistic magnetic monopoles for one year of data taking in 2008 with the ANTARES detector [5] is presented. Section 2 will first introduce magnetic monopoles and their signature in the ANTARES neutrino telescope will be briefly presented. Then the simulation and reconstruction algorithm will be detailled in section 3, and finally the analysis method will be explained in section 4, and the final results will be given in section 5 2 Magnetic monopoles and their signal in the ANTARES neutrino telescope

Magnetic monopoles

Most of the Grand Unified Theories (GUTs) predicts the creation of magnetic monopoles in the early Universe. Indeed, in 1974, 't Hooft and Polyakov showed independently that elements caracterising well a magnetic charge, as introduced by Dirac in 1931 [1], occur as solitons in Gran Unified Theories (GUT) in which a larger gauge group breaks down into a semi-simple subgroup containing the explicit U(1) group of electromagnetism [2,3]. These particles are topologically stable and carry a magnetic charge defined as a multiple integer of the Dirac charge

g D = c 2e
, where e is the elementary electric charge, c the speed of light in vacuum and Planck's constant. Depending on the specific model, the masses inferred for magnetic monopoles range over many orders of magnitude, from 10 4 to 10 20 GeV [6]. Since fast monopoles have a large interaction with matter, they can loose large amounts of energy in the terrestrial environment. The total energy loss of a relativistic monopole with one Dirac charge is of the order of 10 11 GeV after having crossed the full diameter of the Earth [7]. Because magnetic monopoles are expected to be accelerated in galactic coherent magnetic field domains to energies of about 10 15 GeV [8], calculations suggest that monopoles with mass below ∼ 10 14 GeV would be able to cross the Earth and reach the ANTARES detector as upward going signals.

Signal in the ANTARES neutrino telescope

The ANTARES detector is an underwater neutrino telescope immersed in the Western Mediterranean Sea at a depth of 2475 m [9]. In its full configuration the detector consists in 12 mooring lines, each comprising 25 storeys separated by 14.5 m. A storey consists of 3 glass spheres housing a photomultiplier. Finally each line is connected to a Junction Box, which is connected to the shore station at La Seyne-sur-Mer by an electro-optical cable of 40 km length. As for electric charged particles, magnetic monopoles crossing ANTARES with a velocity higher than the speed of light in the sea water of β ∼ 0.74, are expected to emit Cherenkov light. The expected light is however much more intense, with about 8550 times more Cherenkov photons than muons of the same velocity [10]. In addition to the direct photon emission above the Cherenkov threshold, magnetic monopoles can produce indirect light below. Indeed, by ionizing the sea water along its path, a monopole with a velocity β > 0.51 will knock out electrons (δ-rays), which will get enough energy for emitting Cherenkov light [11]. Fig. 1 indicates the amount of light produced during the crossing of a monopole of charge g D as a function of its velocity, compared to the emission of a through going muon. 3 Simulation and reconstruction

Simulation

Upward going magnetic monopoles with one unit of Dirac charge (g = g D ) have been simulated using a Monte Carlo program based on GEANT3 [12] for ten ranges of velocities in the region β = [0.550, 0.995] and with the incom-ing direction distributed isotropically over the lower hemisphere. The simulation of emitted photons is processed inside a cylindrical volume surrounding the instrumented volume. A radius of 480 m, four times larger than that used for the standard ANTARES muon simulation, is chosen in order to take into account the large amount of light emitted by a magnetic monopole. When searching for upward going magnetic monopoles, the main physical backgrounds are upward going atmospheric neutrino-induced muons and downward going atmospheric muons. For simulating atmospheric muons, the CORSIKA program [START_REF] Heck | Technical report FZKA 6019[END_REF] in combination with the QGSJET [START_REF] Kalmykov | [END_REF] code was used, employing the Horandel model for the cosmic ray spectrum, whereas atmospheric neutrinos were simulated according to the Bartol flux [16,17] and the RQPM model [18]. In order to match the real detector conditions, the simulations have been performed using the dead channel information and the optical background rates extracted from the data.

Reconstruction algorithms

Before the reconstruction step, the different trigger algorithms are applied to data before events are stored [19]. These conditions lead to less than 15 % inefficiency for monopoles with a velocity greater than β = 0.58 that produced at least six hits in the detector. The standard track reconstruction assumes that particles travel at the speed of light. In order to improve the sensitivity for magnetic monopoles travelling with lower velocities, the velocity was introduced as a free parameter in the track fit. The modified algorithm performs two independent fits, a track fit and a bright point fit. The former reconstructs the track of a particle crossing the detector at a velocity β, introduced as a free parameter, while the latter reconstructs the event as a point-like light source. Both fits minimize the same quantity χ 2 defined as :

χ 2 = N i=1 (t γ -t i ) 2 σ 2 i + A i (1)
The first term on the right hand side is the sum for N hits of the square of the time residual, where t γ is the expected arrival time of a hit, t i the measured time and σ i the time uncertainty for each hit. The second term A i corresponds to the contribution of hit charge weighted by its minimal approach distance. It brings a penalty to hits with large amplitude combined with a large minimal approach distance. This algorithm yields a resolution δβ on the reconstructed velocity of magnetic monopoles of about δβ 0.025 for velocities lower than the Cherenkov threshold, β C = 0.74, improving to δβ 0.003 for higher velocities.

The analysis has been performed for data taken from December 2007 to December 2008. Quality requirements are applied to the data to ensure low levels of bioluminescent activity and a well calibrated detector. After this selection the data is equivalent to a total of 136 days of live time : 44 days with 12 lines, 46 days with 10 lines and 47 days with 9 lines. For checking the agreement between simulated atmopheric events and real data, 15 % of the whole dataset, corresponding to 20 days of data taking, was taken and used at each step of the analysis. To avoid any bias, the sample of 15 % of data has at the end been removed, and so has not been considered in the final result. Despite the good agreement of shape between data and MonteCarlo atmospheric distributions observed, a normalisation factor of 1.8 is needed to be applied on simulated atmopheric muons, which is consistent with the expected uncertainties on the optical module angular acceptance for downgoing particles and on the input parameters of the model-dependent atmospheric muon flux.

Preliminary cuts

Only magnetic monopoles crossing the Earth are considered in the analysis, and a natural precut is then to only keep particules reconstructed as upward going, i.e. with a reconstructed zenith angle less than 90 • . This cut is motivated by the large amount of background of downward going events coming from atmospheric muons. Moreover a second preliminary cut consists in choosing only particles for which the track was reconstructed by using at least 2 lines, to improve the reconstruction quality. Finally, the last preliminary cut consists in removing events for which the quality factor of the fit is better with the bright point hypothesis than with the track one. Electromagnetic and hadronic showers dominate such events, and a large fraction of poorly reconstructed atmospheric muons are then removed.

Discriminating variables

The event selection for optimising the Model Discovery Factor was performed by using two discriminating variables. As shown fig. 1, the first discriminant between atmospheric events and magnetic monopoles is the number of hits nhit produced in the detector. Fig. 2 represents distributions of the number of hits used by the reconstruction algorithm for events reconstructed in the velocity range β rec = [0.775, 0.825]. A good agreement is observed between simulated atmospheric events and the 15 % data sample, whereas the distribution for magnetic monopoles simulated in the same velocity bin confirms the potential of discrimination for the nhit variable. However, for the lowest monopole velocities below the Cherenkov threshold, using only the number of hits is not sufficient for discrimating monopoles from atmospheric background events, and a new variable has thus been defined by applying two reconstruction algorithms. The first one is the initial muon fit algorithm in which the velocity β is fixed at 1, whereas the second modified algorithm allows β as a free parameter in the procedure. The discriminating parameter λ is then defined as :

λ = log χ 2 t (β = 1) χ 2 t (β free ) (2)
where χ 2 t (β = 1) and χ 2 t (β free ) are the track quality parameters for a fixed and free velocity β respectively. This new variable is then expected to be positive for monopoles and negative for atmospheric muon and neutrino events. In fig. 3 are shown the λ distributions of events reconstructed in the velocity range β = [0.775, 0.825], for atmospheric simulated muons and neutrinos, and for the 15% data sample, as well as for magnetic monopoles simulated in the same range of velocity. Finally, the optimisation of the model discovery factor for a 5σ discovery at 90% probability was performed for 8 bins

Vol. 5, 95 of monopole velocity between β = [0.625, 0.995] by scanning the nhit-λ-plane.

Results

The [nhit, λ] cuts for every velocity bins, defined previously were then applied to the 85 % data sample. The number of atmospheric background events expected after 116 days of data taking is indicated in Table 1 for each bin of reconstructed velocity, as well as the number of observed events when the selection cuts are applied to the 85 % data sample. No significant deviations from the background-only hypothesis are observed by comparing the number of expected events and the number of observed events. All the observed events are therefore considered as background, with no claim for a monopole discovery. The Feldman-Cousins 90 % C.L. upper limit on the upward going magnetic monopole flux is reported in Table 1. The Feldman-Cousins 90 % C.L. flux limit is also reported.

In the flux upper limit calculation, systematic uncertainties were considered. They mainly depend on uncertainties on the detector efficiency. The major contributions to the uncertainties are the OM angular acceptance which contributes to 15 % and the uncertainty in the light absorption length in water which gives an uncertainty of 10 % [21]. To estimate the effect of the loss of detection efficiency, 18 % of hits per event were removed randomly in Monte Carlo monopole simulation. This leads to a deterioration of the upper limit of 3 % above the Cherenkov threshold and between 30 % and 7 % ffor velocities below. The flux upper limit for upward going magnetic monopoles is shown in Fig. 4 as a function of the monopole velocity, and is compared to previous limits.

Conclusion

The analysis performed on 116 days of data taking in 2008 by the ANTARES detector, for the search for magnetic monopoles yields to new constraints on the upward going magnetic monopole flux. This new upper limit improves by a factor of about 3 the AMANDA's one over the range 0.75 ≤ β ≤ 0.995 (γ = 10), whereas for the first time for a neutrino telescope, an analysis is sensitive and improves the flux upper limit below the Cherenkov threshold. going magnetic monopole with g = gD for relativistic velocities 0.625 ≤ β ≤ 0.995 is compared to the published upper limits set by MACRO [22] for an isotropic flux, and Baikal [23] and A-MANDA [24] for up-going monopoles, as well as the theoretical Parker bound [4].

Figure 1 :

 1 Figure 1: Number of Cherenkov photons in the 300 -600 nm wavelength range per centimeter emitted by a monopole with g = gD (solid line), by the δ-rays produced along its path (dashed line) and by a muon (dotted line), as a function of the particle velocity.

Figure 2 :

 2 Figure 2: Distributions of nhit for events reconstructed with βrec = [0.775, 0.825]. Atmospheric down-going muons and upgoing neutrinos (solid histograms) are compared to the 15 % data sample (points). The dashed line stands for magnetic monopoles generated in the same velocity bin.

Figure 3 :

 3 Figure 3: Distributions of λ for events reconstructed with βrec = [0.775, 0.825]. Atmospheric down-going muons and up-going neutrinos (solid histograms) are compared to the 15 % data sample (points). The dashed line stands for magnetic monopoles generated in the same velocity bin.

Figure 4 :

 4 Figure 4: The ANTARES flux upper limit (90 % C.L.) on up-

Table 1 :

 1 Number of background events expected after 116 days, compared to the number of observed events from the 85 % of data.

	βrec range [0.625, 0.675] [0.675, 0.725] [0.725, 0.775] [0.775, 0.825] [0.825, 0.875] [0.875, 0.925] [0.925, 0.975] [0.975, 1.025]	Number of exp. background events 2.2 × 10 -2 1.3 × 10 -1 4.5 × 10 -2 1.1 × 10 -6 8.2 × 10 -7 6.9 × 10 -7 2.3 × 10 -5 1.3 × 10 -2	Number of observed events 0 1 0 0 0 0 0 0	90% C.L. flux u. l. (cm -2 s -1 sr -1 ) 7.5 × 10 -17 8.9 × 10 -17 4.0 × 10 -17 2.4 × 10 -17 1.8 × 10 -17 1.7 × 10 -17 1.6 × 10 -17 1.3 × 10 -17

Vol. 5, 94

Vol. 5, 96


Rivi Ère
email: riviere@in2p3.fr

Carla Distefano

Moon shadow observation with the ANTARES neutrino telescope

Keywords: high energy neutrino, ANTARES, Moon shadow

The ANTARES neutrino telescope is operating in the Mediterranean sea in its full configuration since May 2008. While designed to observe up-going neutrinos, it also records many signals of down-going muons produced by the interaction of high energy cosmic rays in the atmosphere. The shadowing of cosmic rays by the Moon produces a deficit of muons coming from this direction at the ANTARES level. The observation of this deficit of events can be used to check the correct pointing of the detector with respect to a known object.

The strategies to observe this Moon shadow are discussed ant the current results are presented.

Introduction

The ANTARES neutrino telescope is operating in the Mediterranean Sea in full 12 lines configuration since May 2008 [1]. Thanks to the good optical properties of the water of the deep sea, the median angular resolution obtained with the current reconstruction algorithm is estimated to be 0.4 • for cosmic neutrinos with an E -2 flux [START_REF] Adrián-Martínez | Submitted to Astropart. Phys[END_REF]. The pointing performance of the detector relies on the knowledge of parameters, such as the relative delays of the optical sensors within each detector line as well as between the lines, the instantaneous shape of the detector lines which is influenced by the sea current or the absolute orientation of the detector. The values of these parameters have been calibrated before immersion of the detector lines and are regularly measured in situ during operation. The relative time delays are measured using a laser beacon and LED beacons, and also with physical quantities such as the 40 K decay or the atmospheric muon tracks [START_REF] Aguilar | [END_REF]. The alignment measurements are performed continuously with acoustic emitters and transceivers, tiltmetres and compasses [4]. The absolute orientation of the telescope is obtained by the triangulation of acoustic signals between anchors of the lines and the deployment vessel at the sea surface, positioned by GP-S [5]. The uncertainties on these values are small enough to ensure the optimal performances of the telescope, but it remains very important to verify the correct operation of the detector using a calibration source. While no such source exists, some possibilities remain: • The observation of a coincidence between an air shower at sea level with a surface array and the resulting muons at ANTARES level [6].

• The observation of the Moon shadow. This technique, first proposed in 1957 [8], relies on the absorption of cosmic rays by the Moon, as illustrated on figure 1. As the cosmic rays are the source of the down-going muons at the telescope level, this absorption induces a deficit of the number of observed muons in the Moon direction.

2 Data and Monte Carlo

Data set and reconstruction

The apparent Moon radius is on average R M oon = 0.259 • , which is smaller than the resolution of ANTARES, in particular for down-going atmospheric muons; thus the observation of the Moon shadow requires large statistics. For the work reported here, 884 days of data taken between 2007 and 2010 (including periods with 5, 10 and 12 lines configurations) have been analysed. These data have been processed using the standard track reconstruction algorithm of ANTARES. As this algorithm is optimised for up-going tracks of cosmic origin, an algorithm dedicated to the reconstruction of down-going atmospheric muons would improve the results presented in this paper. The event selection criteria used here are the same as presented in [7], i.e. an error estimate β < 1 • and a quality of reconstruction of the track λ > -5.5.

Monde Carlo

The point spread function (PSF) used to generate pseudo experiments and to compute our search strategy is obtained from Monte Carlo simulation (MC). Each data run is simulated as close as possible to the experimental condition (same detector condition, bioluminescence rate, etc.), generating both atmospheric muons and atmospheric neutrinos. An example of data/MC comparison is represented in figure 2.

The PSF is extracted from this MC using all the events reconstructed in a region of 10 • around the Moon position (figure 3). With the chosen selection criteria, the median angular resolution is 0.75 • . Concerning the shape of the event density in the Moon region, we assume the event probability density function (PDF) is uniform in the abscence of Moon1 , we call this hypothesis H 0 . The PSF in the Moon hypothesis (H 1 ) is obtained by computing the 2D convolution product of the PSF with the Moon shape:

P DF (x|H 0 ) ∝ 1 P DF (x|H 1 ) ∝ (1 -P SF ⊗ δ M oon )(x)
3 Search strategy Different search strategies have been studied, a simple on/off search, a χ 2 test of histograms [7] and the one presented here, a likelihood search. The results of these different approaches are compatible.

Likelihood method

For the likelihood search we build a test function based on the likelihood ratio using the expected distributions corresponding to the Moon hypothesis H 1 and no Moon hypothesis H 0 . The test function used is the logarithm of the likelihood ratio:

t = events i log P DF (r i |H 1 ) P DF (r i |H 0 )

Pseudo-experiments

In order to test the search method, to estimate its power and to compute the significance, pseudo-experiments are performed. After generating uniform event distributions in an area around the Moon position, a smearing is applied using the PSF, followed by a shadowing which is or is not applied for the events originating from the Moon direction.

The pseudo experiments presented here are generated with the same event statistics as available in the 2007-2010 data sample.

The distributions of the values of the test function obtained with the pseudo experiments are presented in figure 4. For a given value of the test function, its p-value (probability to obtain a value at least as important from a background fluctuation) is computed using the H 0 distribution (blue). The significance of this p-value can then equivalently be expressed in term of sigma.

The H 1 distribution (red) is used to estimate the power of the test. With the current event statistics, we could expect 2.1 +0.9 -0.8 σ. With 2, 5 and 10 times this statistics, we obtain respectively 2.8 ± 0.9, 4.2 ± 0.9 and 5.9 ± 0.9σ.

Preliminary results

The experimental event density in the Moon region is represented in figure 5 and the density as a function of the distance to the Moon position is shown on figure 6.

Vol. 4, 117 The value of the test function obtained experimentally for the Moon position is Q = 1.17, which has a p-value of 0.007. This corresponds to 2.7σ, which is within the expected range. This is however not significant enough to unambiguously reject the H 0 assumption. Scanning over directions around the Moon position (see figure 7), we find similar hot spots; which are compatible with statistical fluctuations.

Conclusion and outlook

A search for the observation of the Moon shadow with the ANTARES neutrino telescope has been performed on a dataset corresponding to 884 days of livetime. The resulting Moon shadow significance, 2.7σ, is within the expected range and compatible with other analysis [7]. At this point, this result is not strong enough to put interesting constrains on ANTARES pointing capabilities. The increase of statistics or the development of a track reconstruction algorithm optimized for down-going atmospheric muons should improve this results. ANTARES will also perform campaigns of measurement of atmospheric showers by mean of a surface array installed on a boat above the detector, allowing to compare the direction of muons in coincidence with ANTARES. This will provide an additional check of the ANTARES pointing performance.

Figure 1 :

 1 Figure 1: Illustration of the Moon blocking the cosmic rays, thus producing a deficit of muons originating from this direction.

Figure 2 :

 2 Figure 2: Top: Distribution of the reconstructed zenith of the events of the data set used in this analysis: data (black), Monte Carlo (red=yellow+blue). Bottom: Data/MC ratio.

Figure 3 :

 3 Figure 3: Black: Parametrization of the PSF of the events at less than 10 • from the Moon position. Blue: Constant PDF in the H 0 hypothesis. Red: PDF in the H 1 hypothesis.

Figure 4 :

 4 Figure 4: Distributions of the likelihood test function values, for the H 0 hypothesis in blue and for H 1 in red.

Figure 5 :

 5 Figure 5: Event density around the Moon location (0, 0).

Figure 6 :

 6 Figure 6: Event density as a function of the distance to the Moon.

Vol. 4, 115

If we consider a 10 • radius region around the Moon direction, there is actually a counting rate modulation of ±20%, but this is a dipolar modulation which cancels out during the likelihood computation, so it is neglected for simplicity. Only second or higher order modulations could change the likelihood.Vol. 4, 116


Matteo Zecchin
email: matteo.zecchin@eurecom.fr

David Gesbert
email: david.gesbert@eurecom.fr

Marios Kountouris
email: marios.kountouris@eurecom.fr

UAV-Aided Decentralized Learning over Mesh Networks

Keywords: Decentralized learning, unmanned aerial vehicles, machine learning, mesh networks

Decentralized learning empowers wireless network devices to collaboratively train a machine learning (ML) model relying solely on device-to-device (D2D) communication. It is known that the convergence speed of decentralized optimization algorithms severely depends on the degree of the network connectivity, with denser network topologies leading to shorter convergence time. Consequently, the local connectivity of real world mesh networks, due to the limited communication range of its wireless nodes, undermines the efficiency of decentralized learning protocols, rendering them potentially impracticable. In this work we investigate the role of an unmanned aerial vehicle (UAV), used as flying relay, in facilitating decentralized learning procedures in such challenging conditions. We propose an optimized UAV trajectory, that is defined as a sequence of waypoints that the UAV visits sequentially in order to transfer intelligence across sparsely connected group of users. We then provide a series of experiments highlighting the essential role of UAVs in the context of decentralized learning over mesh networks.

I. INTRODUCTION

Decentralized learning is a promising distributed optimization paradigm that enables a network of wireless edge devices to collaboratively train a machine learning (ML) model, and to harness the aggregate computational and data resources of the network. Differently from federated learning (FL) [START_REF] Mcmahan | Communication-efficient learning of deep networks from decentralized data[END_REF], decentralized training is not limited to star network topologies with a central orchestrator, but is more flexible as it allows wireless edge devices to exchange training related information in a peer-to-peer (server-less) manner [START_REF] Tsitsiklis | Distributed asynchronous deterministic and stochastic gradient optimization algorithms[END_REF]. For this reason, they are particularly appealing for future wireless networks based on device-to-device (D2D) communication.

Many different decentralized learning schemes over wireless networks have been recently proposed and analyzed [START_REF] Xing | Decentralized federated learning via SGD over wireless D2D networks[END_REF]- [START_REF] Jeong | Asynchronous decentralized learning over unreliable wireless networks[END_REF]. The major assumption in all these works is that the network topology is strongly connected on average. However, real world mesh networks are characterized by local, rather than global connectivity, and groups of nodes are often isolated or sparsely connected to the rest of the network due to their limited communication range. In these scenarios, decentralized learning is either not possible or its performance is severely hampered.

The work of M. Zecchin is funded by the Marie Sklodowska Curie action WINDMILL (grant No. 813999).

At the same time, unmanned aerial vehicles (UAVs) represent an appealing solution to mitigate limited ground connectivity. UAVs have been used as smart flying relays to improve multi-hop routing capabilities [START_REF] Esrafilian | Autonomous UAV-aided mesh wireless networks[END_REF], to self-organize in flying mesh networks [START_REF] Behnke | Comparison of distributed ad-hoc network planning algorithms for autonomous flying robots[END_REF] and to improve coverage to ground users [START_REF] Sabino | Topology control of unmanned aerial vehicle UAV mesh networks: A multi-objective evolutionary algorithm approach[END_REF]. In this work, we investigate the role of UAVs in aiding decentralized learning protocols over ground mesh wireless networks.

The combination of FL and UAV assisted communication has recently been explored; however, these studies have been limited to scenarios in which the UAV has the role of a parameter server (PS), i.e., aggregating model estimates received from ground nodes and subsequently broadcasting the aggregated model back to the ground [START_REF] Donevski | Federated learning with a drone orchestrator: Path planning for minimized staleness[END_REF], [START_REF] Mrad | Federated learning for UAV swarms under class imbalance and power consumption constraints[END_REF].

Our work differs from these previous works as it considers the UAV serving as a relay and it assumes that ground nodes are able to carry out learning even in the absence of a UAV, by exploiting the already existing ground D2D links. This feature dramatically improves the convergence speed, versatility and fault tolerance of the proposed solution. We propose an optimized trajectory, given as a sequence of waypoints visited by the UAV, which is designed to intelligently provide relaying opportunities to ground nodes and to diffuse locally optimized model across subsets of users with limited connectivity. We provide experiments showing that with the aid of a UAV following the proposed trajectory, it is possible to harness the full potential of the mesh network in spite of sparse and local connectivity, and to accelerate learning compared to UAVaided federated learning algorithms.

II. SYSTEM MODEL

We consider a network of m + 1 devices comprising m ground users plus a UAV serving as a flying relay. We index the ground user by 1, . . . , m and we denote the location of the i-th ground device by p i = [x i , y i , z i ] ∈ R 3 , where x i and y i are the horizontal coordinates while z i denotes the elevation. We assume that ground devices are static, namely their position is not a function of time. On the other hand, the UAV location is denoted p uav = [x, y, z] ∈ R 3 and is assumed to be time-varying in the horizontal coordinates x and y, but not in the vertical one z. Furthermore, the UAV elevation z is set to be larger than a safety altitude z min .

A. Communication Model

Communication among network nodes takes place in rounds. At every communication round t ∈ {τ, 2τ, . . . }, the channel gain coefficient g (t) i,j ∈ R, expressed in dB, between each pair of distinct user (i, j) ∈ [1 : m] 2 is given by

g (t) i,j = g (t) j,i = β g -α g 10 log 10 d i,j + η (t) g (1)
where α g is the path loss exponent, β g is the average channel gain in dB at a reference distance d = 1, d i,j = ∥p i -p j ∥ 2 is the distance between nodes i and j, and η g ∼ N (0, σ 2 g ) models the shadowing effects. For simplicity, we assume that the link parameters α g , β g and σ g are homogeneous across pairs of ground users; however, the proposed solution can easily accommodate heterogeneous channel parameters. At communication round t, the channel gain link between the UAV and a ground node i under Line-of-Sight (LoS) conditions is modeled as

g (t) i,L = β L -α L 10 log 10 d (t) i + η (t) L , (2)
while under Non-Line-of-Sight (NLoS) propagation it follows

g (t) i,N = β N -α N 10 log 10 d (t) i + η (t) N (3)
where

d (t) i = p i -p (t) uav 2
denotes the time-dependent distance between the UAV and user i, α L , β L and η

(t) L ∼ N (0, σ 2
L ) are the channel parameters under LoS, while α N , β N and η

(t) N ∼ N (0, σ 2
N ) describe the channel under NLoS propagation. These parameters are assumed to be homogeneous across users for simplicity.

The LoS probability between the UAV at a position p (t) uav and user i is modeled using the s-model [START_REF] Al-Hourani | Optimal LAP altitude for maximum coverage[END_REF] ρ

(t) i = 1 1 + e -aiθ (t) i +bi (4)
where a i and b i are model coefficients related to the propagation environment, and θ

(t) i
is the elevation angle between the UAV and ground user i. At every communication round t, a link between network nodes is modeled using a simple, yet classical, on-off channel model. Accordingly two nodes communicate if and only if the associated channel gain exceeds a threshold g th . Therefore, the resulting ground connectivity matrix

A (t)
gr ∈ [0, 1] m×m , is symmetric, has diagonal elements being equal to 1 and Bernoulli distributed off-diagonal entries

[A (t) gr ] j,i = [A (t) gr ] i,j ∼ Bern 1 -Φ g th - ḡ(t) i,j √ σ g (5)
where Φ(•) denotes standard Gaussian cumulative distribution function and

ḡ(t) i,j = E[g (t)
i,j ]. Similarly, the connectivity between the UAV and ground users is described by a vector a

(t) uav ∈ [0, 1] 1×m with entries [a (t) uav ] i ∼Bern 1-ρi Φ g th - ḡ(t) i,N √ σ N -ρ i Φ g th - ḡ(t) i,L √ σ L (6)
where ρi = 1

-ρ i , ḡ(t) i,L = E[g (t) i,L ] and ḡ(t) i,N = E[g (t)
i,N ]. Based on the instantaneous connectivity status, determined by the realization of a (t) uav , the UAV serves as a one-hop relay for the communication among grounds users. The connectivity matrix resulting from the relaying opportunities offered by the UAV to ground users is obtained as

A (t) uav = (a (t) uav ) T a (t) uav . (7)
If follows that

A (t)
uav is a symmetric random binary matrix whose entry (i, j) is 1 if an only if there exists a relaying opportunity between ground user i and j, and 0 otherwise.

Overall, the aggregated connectivity matrix, accounting for link existence either by D2D ground communication or thanks to UAV relaying, is given by

A (t) = J m -(J m -A (t) uav ) ⊙ (J m -A (t) gr ) (8)
where J m is the m × m all-one matrix and ⊙ denotes the Hadamard product.

For every realization of the connectivity matrix A (t) , the set of devices connected to node i is

N (t) (i) := {j : [A (t) ] i,j = 1}. (9)
Note that every ground user is connected to itself.

B. Learning procedure

We assume that the goal of ground devices is to collaboratively train a machine learning model in order to benefit from the aggregation of local computational resources and insitu data. In particular, we assume that each ground device is endowed with a local loss function f i : R d → R and local parameter estimate θ i ∈ R d . This corresponds to the distributed empirical risk minimization problem, which entails most decentralized learning problems, whenever {f i } m i=1 are loss terms defined over local datasets {D i } m i=1 and {θ i } m i=1 are local model estimates.

The network objective consists in the minimization of the aggregate network loss subject to a consensus constraint minimize θ1,...,θm

f (θ 1 , . . . , θ m ) := 1 m m i=1 f i (θ i ) (10)
s.t. θ 1 = θ 2 = • • • = θ m .
In the following, we denote the average network estimate as θ = 1/m m i=1 θ i . To solve [START_REF] Sabino | Topology control of unmanned aerial vehicle UAV mesh networks: A multi-objective evolutionary algorithm approach[END_REF], we consider the asynchronous decentralized stochastic gradient descent (DSGD) algorithm proposed in [START_REF] Jeong | Asynchronous decentralized learning over unreliable wireless networks[END_REF]. According to this optimization scheme, ground devices alternate between a local optimization phase based on gradient information (computation phase) and a communication phase to exchange the updated local estimates with one-hop neighbours. To locally optimize the model estimate θ i , we assume that each device i can query a stochastic oracle that is unbiased,

E[g i (θ i )] = ∇ θ f i (θ i ), (11)
and has bounded variance and magnitude

E ∥g i (θ i ) -∇ θ f i (θ i )∥ 2 ≤ σ 2 , (12)
E ∥g i (θ i )∥ ≤ G 2 . ( 13
)
Furthermore, to account for computation impairments and energy constraints, we admit the existence of straggling ground users that can become inactive or postpone the local optimization computation. At each communication round t, the local update rule at device i becomes

θ (t+ 1 2 ) i = θ (t) i , if device i is straggler at round t θ (t) i -η (t) i g i (θ (t-τi) ) otherwise (14) where η (t) i
is a local learning rate and the delay τ i ≥ 0 accounts for the staleness of the gradient information at device i. Subsequently, each device i shares its updated local estimate θ (t+ 1 2 ) i with its neighbours N (t) (i) using either a digital or analog communication protocol [START_REF] Ozfatura | Decentralized SGD with overthe-air computation[END_REF], [START_REF] Xing | Federated learning over wireless device-to-device networks: Algorithms and convergence analysis[END_REF]. The received estimates are then averaged to obtain the new local estimate

θ (t+1) i = j∈N (i) w i,j θ (t+ 1
2 ) j [START_REF] Koloskova | A unified theory of decentralized SGD with changing topology and local updates[END_REF] where w i,j are the entries of the mixing matrix W (t) obtained using a Metropolis-Hastings weighting scheme [START_REF] Xiao | Distributed average consensus with timevarying Metropolis weights[END_REF].

In [START_REF] Jeong | Asynchronous decentralized learning over unreliable wireless networks[END_REF], it has been shown that the performance of the asynchronous DSGD optimization procedure depends both on the activity of users and on the degree of wireless network connectivity. In particular, with more connected network topologies converging faster than sparser ones. This motivates the use of a UAV to facilitate the diffusion of locally optimized models, and to render the decentralized learning protocol more efficient in spite of sparse and local ground connectivity.

III. TRAJECTORY OPTIMIZATION

At every optimization round t, the connectivity matrix A (t) associated to the network of ground users depends on the UAV location p (t) uav and it can be enhanced thanks to the the relaying opportunities it provides to ground devices. In the following, we propose to optimize the UAV trajectory during the optimization process so that the distributed learning procedure is facilitated.

A key quantity that is used to measure the information diffusion capabilities of a network is the expected consensus rate [START_REF] Koloskova | A unified theory of decentralized SGD with changing topology and local updates[END_REF]. While this quantity can be used to characterize the rate of convergence of DSGD procedure, it does not provides a tractable optimization objective to derive the UAV trajectory. For this reason we define a more tractable surrogate objective that yields the optimized trajectory as a sequence of waypoints {w i } n i=1 . In particular, we assume that the initial way point w 0 is equal to the initial UAV location p (0) uav , and that the sequence of waypoints is determined on-the-fly, with the waypoint w i+1 being computed when the UAV reaches the location specified by the previous waypoint w i . The waypoints are designed so as to hover the UAV on a position that maximizes the probability of creating a relaying opportunities between users that, up to communication round t, have not been able to communicate. To this end, we recursively define an link activity rate matrix R (t) as

R (0) = 0 (16) R (t+1) = γR (t) + (1 -γ)E[A (t) ], (17)
for γ ∈ (0, 1). Denoting by t i the communication round in which the UAV reaches the waypoint w i , the subsequent waypoint is obtained solving the following optimization problem maximize

pt (1 -R (ti) ) ⊙ E[A uav ] 1 (18)
where ∥•∥ 1 denotes the entry-wise 1-norm; namely, the sum of absolute values of the matrix entries.

The optimization problem (18) determines the next waypoint so as to maximize the relaying opportunities between pair of users associated to links with a low activity rate. The main challenges in solving (18) are the lack of a close form expression for Φ(•) and the non-convexity of the objective. In order to make the objective differentiable, we approximate Φ(•) using the sigmoid function

S(x) := 1 1 + e -αx ( 19
)
where α is a fitting parameter set to α = -1.702 as proposed in [START_REF] Lieberman | Introduction to operations research[END_REF]. This approximation of the objective allows us to employ efficient gradient based solvers to generate the sequence of waypoints. Nonetheless, the optimization objective (18) remains nonconvex. In order to reduce the probability of obtaining a waypoint associated to poor local maxima, we employ gradient descent with restarts. The number of restart points is chosen to meet the UAV computation constraints and the restart points are sampled uniformly at random inside the convex hull determined by ground user locations.

IV. SIMULATIONS

To test the proposed solution, we consider a network deployment of 30 × 60m 2 with 23 ground devices deployed at the ground level (z = 0m) as depicted in Figure 1a. The propagation parameters describing the ground links channel gain are set to α g = 3, β g = -30 dB and σ 2 g = 1, and the channel gain threshold determining active/inactive links is fixed to g th = -60 dB. In the considered deployment, ground users are naturally clustered together in 3 distinct groups and the path exponent is such that communication within each cluster is possible, but links between users belonging to different clusters are active with negligible probability. Furthermore, we consider an obstruction (gray vertical line) that amounts to a 35 dB attenuation for the links between users residing on opposite side of the line. A UAV flying at a fixed altitude of 10m serves as a relay to enhance ground connectivity. The channel gain parameters describing the link between ground users and the UAV under LoS propagation are α L = 2.5, β L = -30 dB and σ 2 L = 1, while under NLoS are α N = 3, β N = -30 dB and σ 2 N = 1. We assume that the devices store only 10 data samples from the FashionMNIST dataset, which alone would not guarantee good inference capabilities. Therefore, they wish to harness the distributed dataset to jointly train a machine learning model. In the following experiments we consider a fully connected neural network with one hidden layers comprising 25 neurons. Ground devices update the model employing a gradient descent optimizer and a geometrically decaying learning rate η

(t) i = 0.1 • (0.995) t .
In this setting, the role of the UAV is to intelligently create relaying opportunities so to promote collaborative learning, and to facilitate the ground users to harness the entire distributed dataset by global diffusion of the locally optimized models, in spite of sparse and local connectivity.

We benchmark the proposed solution, corresponding to a UAV visiting the sequence of waypoints returned by (18) and serving users for 20 optimization rounds each time a waypoint is reached, and compare it with alternative trajectory optimization schemes.

In particular, we consider the cluster mid-points traversal trajectory according to which the UAV first runs a k-means clustering algorithm to detect natural clusters of the ground devices positions, and then visits sequentially the mid-points between each pair of cluster centers. Once the mid-point is reached, the UAV serves ground users for 20 optimization rounds and then it flies to the next location.

We consider the barycenter placement, in which the UAV hovers at a fixed location p bar uav determined by the mean ground user location.

Finally, we consider the maximum connectivity placement in which the UAV location is fixed and set equal to the coordinate that maximizes the probability of creating a relay links. This is obtained solving the following maximization problem

maximize p ∥(1 -A gr ) ⊙ E[A uav ]∥ 1 . (20)
For both the barycenter and the maximum connectivity placements, the UAV serves ground users at every communication round as it hovers at a fixed location for the entire learning procedure. We consider the centralized learning solution, corresponding to a fully connected topology, as a performance upper bound. For all listed approaches we run the distributed optimization protocol for 500 rounds and we track the testing accuracy of the mean network estimate θ. We also consider the network consensus error metric

ε(θ 1 , . . . , θ m ) = 1 dm m i=1 θ i -θ 2 2 , (21)
which measures the degree of disagreement among ground nodes estimates and it is a proxy to assess how effective is the UAV trajectory in satisfying the consensus constraint in (18).

In Figure 1a we plot the ground user locations (black dots) and we overlay the UAV trajectories during the training process. Specifically, in the top left corner, we report in green the UAV trajectory returned by the proposed scheme. The UAV frequently hovers between the disconnected components in order to relay information between the clusters and to diffuse model estimates between groups of interconnected ground nodes that rarely communicate using D2D ground links. In the top right corner, we provide the trajectory of the cluster mid-points traversal solution (gray). The UAV successfully identifies the clusters and it sequentially visits the mid-points. This strategy enhances the ground connectivity, but it fails at providing relaying opportunities across the two bottom components that are disconnected due to the propagation obstacle. In the bottom row, we plot the barycenter (left plot) and maximum connectivity (right plot) placementes. Both solutions yields a static UAV placement that is fixed throughout the entire training phase.

In Figure 1b we report the testing accuracy attained by the mean network estimate θ when decentralized learning is assisted by a UAV flying according to the different trajectories. The testing accuracy allows us to quantify the extent to which the UAV is beneficial to the collaborative learning process. The proposed solution (in green) is able to take full advantage of the distributed dataset, and it successfully enables fast distributed training with a final accuracy level that matches the accuracy of the centralized solution (in black). The barycenter solution (in red) converges slowly to a lower accuracy level, highlighting the necessity of a dynamic UAV placement to take full advantage of the network resources. The cluster mid-points traversal solution, despite enhancing ground connectivity, it is not able to connect all the disjoint components and therefore it converges to suboptimal solution. Similarly, the maximum connectivity placement is not able to connect all the disjoint network components and to transfer intelligence across different clusters.

In Figure 1c we report the network consensus error evolution attained by the different UAV trajectories. While the proposed trajectory is able to reduce the consensus error during training and it eventually ensures that the edge devices reach a common learning goal, the other baselines are not able to drive network nodes to a globally shared model estimate.

Finally, we propose a comparison between the proposed decentralized learning scheme and a UAV-aided federated learning protocol, as in [START_REF] Donevski | Federated learning with a drone orchestrator: Path planning for minimized staleness[END_REF]. In particular, for the federated learning algorithm, we assume that the UAV serves as a PS, and it orchestrates the learning procedure by collecting locally optimized models by the network devices and broadcasting aggregated estimates back to the ground users. On the other hand, in case of decentralized learning, the UAV serves as a relay and the ground devices can also exploit the available D2D ground links to exchange model estimates, in principle being able to perform learning without the presence of a UAV. As a result, the proposed protocol is more flexible with respect to the communication topology, it can easily accommodate multiple assisting UAVs, and converges faster. To compare these two approaches we study the same deployment as in Figure 1a. We assume that the relaying UAV follows the proposed trajectory, while the UAV serving as orchestrator follows the trajectory obtained solving (18) setting A (t) gr = 0, trying to serve large groups of users prioritizing stale ones, akin to [START_REF] Donevski | Federated learning with a drone orchestrator: Path planning for minimized staleness[END_REF]. In Figure 2 we report the testing accuracies attained by the protocols. The proposed approach drastically reduces the training time, halving the number of iterations required to reach the final performance obtained by the federated learning protocol.

V. CONCLUSION

In this work, we studied the benefits that a flying relay can bring to a network of wireless devices that are jointly training a machine learning model. We proposed a trajectory optimization scheme that enhances the ground connectivity so as to facilitate the diffusion of locally optimized model estimates, and that enables ground users to take full advantage of network computational and data resources. We provided a series of experiments highlighting how a properly designed UAV trajectory can greatly promote decentralized training and outperform UAV-aided federated learning protocols.

  Different UAV trajectory and placements. Black dots represent ground users and the gray vertical line is a propagation obstacle that corresponds to a 35 dB attenuation. Testing accuracy averaged over 5 runs, obtained by the mean network estimate, using different UAV-Aided decentralized learning protocols. Evolution of the average consensus error (21) attained by the benchmarked trajectories. Smaller consensus error corresponds to disagreement between network nodes.

Fig. 2 :

 2 Fig.2: Testing accuracy, averaged over 5 experiments, obtained by the mean network estimate when training is aided by a UAV serving as a relay to assist the decentralized learning protocol (green), or as an orchestrator to perform federated learning (gray dashed).


S L Avila
email: avila@grucad.ufsc.brcarpes@grucad.ufsc.brbergmann@vrac.puc-rio.br

W P Carpes

J R Bergmann

L Krähenbühl

Conception and Comprehension of the Contour Beam Antenna in a Satellite System by Multi-objective Genetic Algorithm with Sensitivity Analysis

Keywords: Multi-objective Optimization, Sensitivity Analysis, Contour Beam Antenna, Satellite System

To obtain the maximum performance of any device the designer needs to have a large experience in the entire universe where it is inserted. Consequently, more important than to find an adequate optimized solution is to comprehend the whole problem, giving to the designer opportunity to make choices and also to learn about the problem. A complete optimization procedure with sensitivity analysis is shown in this paper, which the chief aim is to increase the comprehension of complex engineering conception problems. The multi-objective genetic algorithm and two different approaches of sensitivity studies are shown running together. The contourbeam-shaped reflector antenna for specific illumination in a satellite system is solved as application problem.

2.

Introduction

Contour beam antennas are used on geostationary satellite to provide an efficient coverage of the service area. In many cases, their design has also to satisfy copolar and crosspolar isolation between different areas, gain slope at edges of the coverage area, and multiple frequency bands. As satellite expected life service increases, it creates a demand for reconfigurable contour beam antennas to cope with the changes in the service area. It defines a scenario where the antenna design has to simultaneously satisfy multiple objectives. Among the possible configurations used for onboard antennas, shaped reflector antennas fed by a single feed have been considered. The reflector shaping techniques employ an optimization method to adjust the reflector surface in order to satisfy the desired specifications. The usual optimization methods are deterministic [1,2], leading to an optimum solution at the end of iteration and leaving no freedom for the designer to accommodate further adjustments. Multi-objective Genetic Algorithm (MGA) is a stochastic-based method that yields a set of efficient solutions (Pareto-front), that simultaneously satisfies several goals [3,4]. The analysis of the Pareto-front behavior permits to understand the tradeoff between the different objectives. Compared with the deterministic optimization methods, which lead to unique solution, MGA offers the possibility to the designer to make the final choice among the set of solutions by considering additional constraints not included in the initial steps. In this work, we apply the MGA in the synthesis of a shaped reflector antenna that is designed to provide Brazilian territory illumination and simultaneously minimize the energy level on the adjacent area. Beyond the search for the coverage specifications, the sensitivity analysis is included. This analysis consists in studying the real world inevitable deviations from project specifications, like constructing inaccuracies or different operational circumstances. We work here with two different sensitivity approaches. The first one consists in a classical approach for optimization parameters [START_REF] Vieira | Multi-objective sensitivity analysis in finite domains of a Yagi-Uda antenna[END_REF], but using only and directly the data from the stochastic MGA process, to study variations in the construction devices step mainly. In the second approach, the sensitivity study concerns the model specifications [START_REF]the 8. ference 9. problems solved by stochastic methods[END_REF]. It means that a complementary analysis is applied to fixed parameters, which are not used as optimization parameters. For instance, this approach could to verify the effects of satellite position changes on the geostationary orbit. The results obtained via the multi-objective optimization including sensitivity analysis help to improve the comprehension of the problem.

Multi-objective Genetic Algorithm

The Genetic Algorithm (GA) is a stochastic procedure based on the concepts of natural selection and genetics. There are many papers showing the effectiveness of the GA to solve engineering optimization problems. In most real-world problems, several goals must be satisfied simultaneously in order to obtain an optimal solution. As these objectives are usually conflicting, no single solution may exist that is best regarding all considered criteria. Multi-objective optimization (also called multicriteria, multiperformance or vector optimization) seeks to optimize the components of a vector-valued cost function. Unlike single objective optimization, the solution to this problem is not a single point, but a family of efficient points. Each point on this surface is optimal in the sense that no improvement can be achieved in a cost vector component that does not lead to degradation in at least one of the remaining components. Each element in the efficient set constitutes a non-dominated (non-inferior or non-superior) solution to the multiobjective problem. The main action of the multi-objective optimization is to determine the efficient front. With this set of solutions, it is possible to understand the dependence between each objective, which allows making efficient choices for the final solution decision. Fundamentally, in a multi-objective optimization procedure, given a set of n objective functions (f i ) to be maximized (>) or minimized (<), one should find the set P of solutions defined as:

P = { it does not have to exist x ∈ S | 1) exist i such that f i (x) > f i (x*) and 2) ∀ i≠j f j (x) ≥ f j (x*) } , ( 1
)
where S is the feasible search space. The set P contains the efficient solutions (x*) of the problem and it is called the Pareto-optimal set [3,4]. With this done, we can choose one solution according to some preference (decision step).

The multi-objective genetic algorithm (MGA) used in this work is based on three current populations. Figure 1 shows the procedure.

Basically, the algorithm starts with a set of solutions randomly created. These solutions are evaluated and the Pareto-optimal condition (1) is tested, giving two groups of solutions: one formed by efficient solutions, called non-dominated population (NDOM); and another by non-efficient solutions, called dominated population (DOM). An index (IDOM) indicating how many times each solution is dominated by others is created. This is important to execute the selection process only with solutions around the Pareto set, which makes the convergence faster. After the Pareto's check, it is time to apply the Clearing technique, whose purpose is to obtain a sparse and well-established Pareto front. If similarities among individuals are detected (in parameters or/and objectives spaces [7]), one or some of them are punished. The penalty consists in moving the penalized individual to DOM (by changing IDOM from 0 to 1).

To better control the whole evolutionary process we work with a fixed number of individuals (nbind) in the genetic operators actions (crossover and mutation). We call this fixed minimal size group "real population" (REAL), which is recreated at each generation. REAL is composed by all solutions of NDOM (after Clearing); plus 1/4*nbind solutions of DOM to maintain a certain diversity (chosen among those having the minor IDOM indexes). If the number of REAL's individuals is still lesser than nbind, we complete with individuals of DOM (chosen another ones having the minor indexes once again). The opposite will be handled in the selection process. The number of selected individuals picked up from REAL is always nbind. The selection is done by pulling together a Deterministic Sampling and Tournament methods. The parents are in part selected by Deterministic Sampling (based on the population fitness average, which gives emphasis to the Pareto center area) and the remaining by Tournament (based on each objective individually n objectives times, which gives emphasis to the Pareto edges areas). This approach makes easier the attainment of a well-established Pareto set. After selection we apply the crossover and mutation operators for real coding schemes [8].

Generation after generation, these operators create new individuals (children) from existing ones (parents), driving to a good "search space" investigation. Attention, children do not automatically respect all the constraints (fixed maximal dimensions, for example). In this case, classical techniques are applied: for instance, values above a given limit are reduced using a simple "adjustment by saturation" [4].

The new individuals (children) are evaluated and directly inserted in REAL, whose all solutions are submitted to Pareto-optimal condition (1). This results in a modified NDOM, whose size fluctuates (increases and decreases) continually, while the DOM can only grow. Similar to a mono-objective GA, new good solutions can appear at any step in the MGA methodology. However, these new individuals can be lost or destroyed during the evolutionary process. The concept of Global Elitism [8] guarantees the permanence of efficient solutions. In the multi-objective algorithm proposed here, the global elitism is implicitly incorporated in the use of all NDOM to compound REAL. Finally, the evolutionary process is restarted with the new DOM and NDOM. To reduce the computational cost which is necessary to evaluate IDOM, the size of DOM population is controlled (to keep this maximum size some dominated solutions -with a bad IDOM index -are moved to an outside database (DOMold)). The end of the evolutionary process can be done by some criterion (e.g., maximum number of generations).

Generation process: convergence by number of non-dominated solutions or/and generation maximum number.

Initial Solutions randomly created

Evaluation (problem solutions analysis) Global Elitism "Real" + children (no solution is lost) Pareto's condition

Sensitivity Analysis

In engineering design, the aim of the sensitivity analysis is to verify how stable a device performance is, since deviations from project specifications are inevitable in real world (imprecision of the construction, change of operational conditions due to deformation forces, …). Two sorts of studies are done here. The first one is the sensitivity analysis in the parameters space and the second one is the sensitivity analysis to the model (predetermined parameters not present in the optimization process).

Sensitivity Analysis in the Parameters Space

Among several possible sensitivity estimators, the Lipschitz constant has been used. It is defined as the maximum variation rate in a finite domain, has been used [START_REF] Vieira | Multi-objective sensitivity analysis in finite domains of a Yagi-Uda antenna[END_REF]. Another method, based on containing ellipsoids, has been proposed [START_REF] Krähenbühl | Radiation pattern control by subreflector shaping in a dual-reflector antenna," 2. technique for high performance reflector antennas[END_REF] in order to take profit from the large amount of data generated in the stochastic search. In all cases, the precise evaluation of sensitivity quantities is intrinsically very much time consuming. The procedure used here is based on a metric directly evaluated from the data provided by MGA. This methodology also permits to identify which solution parameter degrades more the performance when deviated. In order to validate the method, we apply it to an analytical function minimization and following to the antenna problem.

The main purpose of the MGA is to search for, among all the generated data, well-established samples of the Pareto-optimal front. If the optimization method convergence is not premature, the feasible space is well represented [7]. Many samples fall around Pareto's region while the MGA is searching for the solution, as shown in Figure 2a, so that it is reasonable to use this data in the sensitivity analysis.

In a first step, we search for the maximum variation on the values of each objective function related to deviations from the optimal parameter values (V * ). Following, we select some samples of the Pareto-set are selected to perform the sensitivity analysis. Both metrics are evaluated using only the already acquired samples inside a perturbation domain around V * . Figure 2a shows the data generated by MGA in the simultaneous minimization of two quadratic objective functions (f) with two design variables (x):

( )

2 2 2 1 1 1 x x f + - = and ( ) 2 2 2 1 2 1 x x f + + = .
(2)

The optimal region is a line connecting the parabolas vertices. In this analytical example, we choose the Pareto sample V * =[-0.003387, -0.004788], which is the closest solution to fitness [1, 1] found by MGA. We defined the perturbation domain D * as a deviation of [-0.2, 0.2] 2 around V * , and then took the sample with the worst fitness inside it. In this example, the worst cases for both objective functions are on the domain corners. The designer considers the solution stable if this worst fitness is acceptable. To compare how stable is a Pareto sample in relation to others, it is essential to associate the fitness values with the distance in the parameter space. We define the behavior associated with an optimal point V * for each objective function f k (⋅) inside the perturbation domain D * as:

j i j k i k j i D V V k V V V f V f V B j i - - = ≠ ∈ | ) ( ) ( | max ) ( , , * * , (3)
where V i and V j are samples inside D * . Note that D * can be any domain containing V * . The evaluation of the behavior is straightforward, since we use only points already sampled. Furthermore, it is a value asymptotically close to the Lipschitz constant as the number of samples tends to infinite all over D * . That means it is an approximation to the maximum mean rate between two points inside D * . The greater B k is, the more unstable the function f k at V * will be. Since distance and objective function variation are positive values in B k evaluation using (3), we may normalize them according to the optimal solution in analysis V *  normalize each parameter by its respective V * component absolute value, and function variation by |f k (V * )|  in order to connect directly sensitivity to significant digits.

It is enough to perform the behavior analysis to each variable separately in order to point out which solution parameter is more relevant to maintain a good device performance. In this case, to create the domain D * we achieve a large perturbation to the parameter in analysis (like for all variables sensitivity analysis case) and small perturbations (a tiny percentage) to the others ones. The procedure shown is not exact, but it gives us a satisfactory notion of sensitivity. Moreover, no additional computation effort is necessary if the optimization stochastic method works well (allowing a good exploration of the feasible search space). The analysis method described here uses direct variables. If using indirect parameters like expansion series, the proposed analysis could not have a real meaning. In that case, it may be necessary to create a parameterization for the perturbation, which must be linked to that used in the optimization process, in order to keep on using the samples. The antenna problem solved in this paper is like that: the shaped surface is described by an expansion global series. We show that changes in results section ( §7.1).

Sensitivity Analysis to the Model

Besides the search for the main objectives, it is also possible to include the analysis of the effects of changes in the system specifications. It is a sensitivity analysis of the model specifications [START_REF]the 8. ference 9. problems solved by stochastic methods[END_REF], which means, differently from the classical parameters sensitivity study (previous section), a complementary study applied to parameters that are not used as optimization parameters. As examples, in our case, displacements of the feed and reflector positions, the effects of changes of satellite position on the geostationary orbit, etc. Hence, these are details that the designer does not control but he can estimate possible alterations.

The sensitivity analysis concerning the model can be done after the efficient solutions set is found. The procedure is: all Pareto solutions and some dominated solutions (those who are the closest to the Pareto front -IDOM index) are re-evaluated with the changes in the model, getting thus a new disposition of the front. To estimate the solution stability, we use:

( )

1 1 1 , y y x x d - - = , (
)

2 2 2 , y y x x d - - = , 2 1 d d d + = (4)
where d is the deviation imposed on the solution (for instance, two perturbations: d 1 and d 2 ) and x and y are two objectives. Of course, it is not recommended to add d for non-correlated perturbations. In this case, the analyses must be done separately or by max(|d 1 |,|d 2 |), for example. There is one point to pay attention: it is important to know the difference between fixed specifications and possible optimization variables or objectives. System demands must be objectives. If a fixed specification is very important to maintain the system performance, may be it is better to put this parameter as optimization variable. Model perturbations have as goal to improve the final solution quality. The sensitivity analyses can be optimization objectives, but this increases the complexity and the computational cost of the problem resolution.

In the results section ( §7), both sorts of sensitivity analysis will be applied, helping in the comprehension of them.

Satellite System Arrangement

The satellite antenna configuration is a single offset reflector illuminated by a corrugated feed horn, operating at 12 GHz (Ku Band). The reflector projects a circular aperture A with diameter D, offset point H and equivalent focal distance F, like shown in Figure 3. The cos q θ model describes the feed radiation pattern where q is chosen to provide a 12 dB attenuation at the reflector edge. The satellite beam points to the center of the Brazilian territory, as illustrated in Figure 3. The reflector shaping has to maximize the minimum directive gain (GT) in the coverage area and simultaneously minimize the maximum directive gain (GI) in the adjacent area (gray), in order to avoid interference. In the synthesis, these gain specifications are attained by employing MGA to adjust the expansion coefficients of the Fourier-Jacobi series [START_REF] Vieira | Multi-objective sensitivity analysis in finite domains of a Yagi-Uda antenna[END_REF]. In order to get information on how good is a given radiation pattern, N sampling points are placed on the whole Brazilian territory and isolation area. The distance between the sampling stations must agree with Nyquist rate and Gibbs' phenomenon [10]. The fitness function is just the search for the worst evaluation for all sampling points (to coverage area: min(GT), to isolation area: max(GI)). The directive gain is obtained by radiated far field computation [START_REF] Silver | Microwave Antenna Theory and Design[END_REF]. It can be evaluated by using Physical Optics (PO):

inc PO H n J r r × ≈ 2 , ( 6
)
where inc H v is the incident magnetic field. The radiated far field of the PO can be written as

[ ] ∫ - - ≈ - ' '. ' ) . ( 4 s r r jk PO PO jkr PO ds e r r J J r e k j E r r r r π η , (7)
where k and η are, respectively, the wave number and the intrinsic impedance of the medium, s' is the reflector surface, rˆ is the unit position vector of the observation point and ' r r is the vector of the source point. The PO analysis is not accurate when used for predicting far-angular regions or cross-polarized fields. In order to improve the precision of ( 7), edge diffraction can be taken into account [START_REF] Rego | Equivalent edge currents for the tim antennas[END_REF]. The total field is then given by the direct (feed field

) and indirect (PO field

f E r PO E r
and fringe field fr E r ) components:

fr PO f E E E E r r r r + + = . ( 8
)
The antenna chosen and the application problem are simple to facilitate the comprehension of the optimization procedure, not being necessarily the most correct for a real communication system.

Optimization Results

Several MGA executions had been carried out to guarantee the Pareto-front found. Figure 4 displays the set of solutions where it is possible to observe that the increase in the value of GT causes the increase of GI. This is simple to explain: to augment GT, it is required to spread the energy, which results in an increase of GI. This comprehension comes from MGA -Pareto front. Understood the tradeoff between the objectives, the designer can proceed with the decision step to choose one solution according to some preference. Otherwise, beyond the attainment of the efficient solutions, we can start the sensitivity analysis complementary studies.

Application of Sensitivity Analysis in the Parameters space

We have chosen two solutions to perform this sensitivity analysis: S1, with GT = 25.50dBi and GI = 7.08dBi; and S2, with GT = 26.50dBi and GI = 13.53dBi (see Figure 5). It is required to change the proposed sensitivity procedure to do this study here. In this case, the parameters adjusted are not direct variables but rather a global expansion series. As a result, little changes in only one parameter could not have real meaning, for example. To control real-world perturbations on the reflector shape, like thermal deformations or mechanical efforts, we can work in a reverse approach (without using complex thermal models) like that: accepting a variation of 0.25dBi in GT objective and 0.5dBi in GI objective as tolerable, a finite domain of solutions is defined around the selected S1 and S2. We search into this domain the solution that has the biggest disparity of its surface in comparison with the selected solution shape. That solution having the biggest difference, with the same domain size, is considered the more stable solution. Of course it is necessary to certify if the perturbation occurred is real and possible to happen.

Figure 5. Selected solutions and respectively finite domains (in dBi).

Figure 6 shows the difference of the extremes solutions to S1 and S2 (only one extreme for each solution as illustration). In this case, considering the four extremes for each solution, we can say that S2 is more stable than S1. It means that S1 changes more yours directive gains when little perturbations are done in yours parameters than in S2. This can be a kind of complementary study to understand the problem behavior. Figure 6. Differences in "z" coordinate between the selected solutions and theirs respective extremes.

Application of Sensitivity Analysis to the Model

For testing the procedure, we chose two perturbations:

Satellite displacement

Due to losses, satellites position in the geo-stationary orbit shows small variations that require further correction procedure. To take into account these orbital oscillations, this sensitivity study considers ±1.0% perturbation (∆r) in the satellite position. Table 1 shows the minimum Table I. Sensitivity Analysis concerning the model. CHANGES 7 (a) show that the ∆r = -1% in the satellite orbit defocuses the beam bringing, as consequence, the reduction of GT and the increase of GI. While, for the ∆r = -1% the opposite is verified. Figure 7 (b) shows that the GT are largely affected by lowering the initial specifications (F, D, H and θ o ). With the Pareto-fronts found, we can calculate the deviation (3) to find the most stable deviation among the solutions for each kind of perturbation and also the two alterations together (for that it is required to re-evaluate with both alterations simultaneously). Figure 8 shows the radiation pattern of the most stable solution found with a minimum of 26dBi direct gain in the coverage area. Figure 9 shows differences in "z" coordinate between the final solution and a perfect parabolic shape. Although the optimized reflector shows not big deviations from a parabola, the difference is more intense at the edge (0.11m maximum -approximately 3λ) of the reflector and it defocus the reflected rays in order to provide the adequate coverage of the desire far-field region. Conclusions lti-objective optimization with sensitivity analysis has been presented in this work. In order to explain the

8.

A study of mu methodology, we considered the design of a shaped reflector antenna used in a satellite communication system. The main goal of this work is to improve the comprehension of engineering problems optimization. MGA features allow a large examination of its results and this is the biggest difference from deterministic methods. Based on all analyses and results in hand, the designer can either choose the final optimized solution in the Pareto set, or restart the optimization process considering different specifications and/or parameters. In that case, he will benefit from his better comprehension of the problem.

Figure 1 .

 1 Figure 1. Multi-objective Genetic Algorithm procedure.

  (a) Parameters space (a') Pareto-front -Fitness space (b) Parameters space zoom (b') Fitness space zoom Figure 2. Two parabolas -analytical function test (2).

Figure 3 .

 3 Figure 3. Geometry of the single-feed offset reflector antenna and coverage problem proposed.The shaped reflector surface is described by an expansion of Fourier-Jacobi series [2]:

Figure 4 .

 4 Figure 4. Pareto-front (scales in dBi).

  Figure 7 (a) shows the solutions behavior (GI and GT are given in dBi).Feeder × Reflector SystemAlterations due to thermal gradient, mechanical forces and other external agents can cause changes in the distances and angles of the feeder × reflector system. In this exercise we work with the focal distance (F=1.524m), diameter reflector (D=1.803m), feeder displacement offset (H=1.278m) and feed offset angle (θ o =40°). The perturbation was set as ±0.25% in the mentioned parameters (e.g., θ o =40°±0.1°). Figure7(b) shows the solutions behavior.

Figure 7 .

 7 Figure 7. Pareto-fronts (in dBi) on Sensitivity Analysis to the Model.

7. 3

 3 Final Decision Step (Pareto-front behavior and both sensitivity studies), the designer can choose the final optimized solution or al deviation in Based in all these analysis restart the optimization process with different and/or more specifications to use as advantage a better knowledge of the problem.

Figure 8 .

 8 Figure 8. Footprint radiation pattern of the final solution (COPOL in dBi).

Figure 9 .

 9 Figure 9. Differences in "z" coordinate (m) between the final solution and a perfect parabolic shape.