Jean-Thomas Camino
email: jean-thomas.camino@airbus.com

Christian Artigues
email: artigues@laas.fr

Laurent Houssin
email: houssin@laas.fr

Stéphane Mourgues
email: stephane.mourgues@airbus.com

Linearization of Euclidean Norm Dependent Inequalities Applied to Multibeam Satellites Design

Keywords: Mixed Integer Linear Programming, Mixed Integer Nonlinear Programming, Euclidean Norm Linearization, NP-Hardness, Multibeam Satellites

Euclidean norm computations over continuous variables appear naturally in the constraints or in the objective of many problems in the optimization literature, possibly defining non-convex feasible regions or cost functions. When some other variables have discrete domains, it positions the problem in the challenging Mixed Integer Nonlinear Programming (MINLP) class. For any MINLP where the nonlinearity is only present in the form of inequality constraints involving the Euclidean norm, we propose in this article an efficient methodology for linearizing the optimization problem at the cost of entirely controllable approximations. They make it possible to rely fully on Mixed Integer Linear Programming and all its strengths. This methodology is successfully applied to a critical problem in the telecommunication satellite industry: the optimization of the beam layouts in multibeam satellite systems. We provide a proof of the NP-hardness of this very problem along with experiments on a realistic reference scenario.

Introduction

In the wide literature on mathematical optimization, there are several examples of problems involving continuous point variables in R 2 or R 3 with constraints on the Euclidean distance between pairs of such points. In some of these problems, the possibility to rely on convex optimization is preserved: for instance, an upper-bound on a Euclidean distance between two points is a convex constraint. On the other hand, any equality or lower-bound set on a Euclidean distance makes the corresponding optimization problem non-convex. In any case though, these quadratic constraints position the optimization in the Nonlinear Programming (NLP) class. One of the most famous problems that handles such 2-norm computations over continuous variables is the Euclidean Multifacility Location Problem (EMFL): [START_REF] Xue | An efficient algorithm for minimizing a sum of euclidean norms with applications[END_REF], [START_REF] Gl Xue | A polynomial time dual algorithm for the euclidean multifacility location problem[END_REF]. This problem consists in defining positions for n new facilities with respect to m existing facilities. The minimized cost function terms are proportional to distances between pairs of new facilities, and to pairwise distances between old and new facilities. Most algorithms that solve the EMFL rely on secondorder cone programming and interior point techniques for convex optimization since the problem can be equivalently transformed into another one where convex quadratic proximity constraints appear. For works that consider both proximity and separation constraints, we can mention for instance the research on wireless sensor localization (see [START_REF] Mao | Wireless sensor network localization techniques[END_REF] for a thorough survey on the matter). In this problem, we assume that a set of sensors has been deployed on a certain region, and that some of the sensor positions are known while the others are not, the goal being to estimate these unknown positions. The authors of [START_REF] Biswas | Semidefinite programming for ad hoc wireless sensor network localization[END_REF] have indeed to deal with non-convex equality and separation constraints which they choose to relax just enough to reach a semidefinite programming model. In some cases, discrete variables are necessary to model decisions with finite numbers of possibilities. When combined with the Euclidean norm constraints discussed above, these integer variables lead to Mixed Integer Nonlinear Programming (MINLP), which is known to be one of the most difficult optimization problems class ever to be tackled. As an example of such problems, the authors of [START_REF] Sutou | Global optimization approach to unequal global optimization approach to unequal sphere packing problems in 3d[END_REF] worked on the issue of packing unequal spheres in a 3-dimensional polytope with sphere separation constraints and an objective to maximize the volume occupied by the spheres, the application being radiosurgical treatment planning. In that case, the continuous variables are the sphere centers, and the discrete variables correspond for each sphere to the choice of a radius among a finite set of possibilities. This non-convex quadratic problem is solved with an heuristically improved simplicial branch-and-bound method. Another example in the satellite industry is the problem of optimizing the beam layouts of a multibeam telecommunication satellite system (see [START_REF] Camino | A greedy approach combined with graph coloring for non-uniform beam layouts under antenna constraints in multibeam satellite systems[END_REF] for instance). It consists of defining the positions in the Euclidean plane of a certain number of disks, each one representing the spatial extent of a radiofrequency beam that carries telecommunication signals for various applications: television, telephone, radio or internet by satellite for instance. While the number of user ground stations (modelled by points of known coordinates in the Euclidean plane) covered by these disks of discretely varying diameter is maximized, satellite antenna technological constraints force some couples of disks to be sufficiently separated. Note that this last application is the one that motivated this work.

In the end, the focus of this article is laid on the particular MINLP problems where the only nonlinearities are quadratic constraints expressed with a Euclidean distance over continuous variables. Both separation and proximity constraints are handled, and they cohabit with discrete variables, in order to optimize a linear objective. For these problems, we detail in section 2 the methodology we devised for reaching a Mixed Integer Linear Programming model (MILP) thanks to controllable approximations, the goal being to take full advantage of all the efficient techniques developped for this specific class of optimization problems. In section 4, the beam layout optimization is defined more in details, its NP-hardness is proven, and the principles of section 2 are applied in order to reach a MILP model which is presented and commented. Experiments on a reference scenario are presented in section 5 before some concluding remarks in section 6.

2 Linearization of Euclidean norm dependent constraints in R 2 Let X ∈ R 2 and α, β ∈ R + , we develop in this section a MILP-compatible linearization process of the constraints of the following form: ||X|| ≤ α and ||X|| ≥ β with || • || denoting the Euclidean norm. The associated inner product will be denoted by < • | • >. The choice of R 2 is directly motivated by the beam layout optimization application, but the principles presented could be generalized to R 3 , or to higher dimensions. Also note that, even if it is not discussed here, these techniques could also be extended to the constraints of the form

||X 1 || + • • • + ||X k || ≤ α and ||X 1 || + • • • + ||X k || ≥ β with k > 1.
On the topic of linearizing the Euclidean norm in the context of mathematical optimization, the examples are rare but can be found in the applications overviewed in the introduction section 5. In the context of radiotherapy equipment configuration, the authors of [START_REF] Liberti | Optimal configuration of gamma ray machine radiosurgery units: the sphere covering subproblem[END_REF] propose to linearize the quadratic terms of the convex proximity constraints with added variables and a notion of approximation points, but without really discussing the error made in the end on the approximated Euclidean distances. Another way of linearizing the Euclidean distances is to discretize the possible positions of the originally continuous variables allowing then to pre-compute all the possible point-topoint distances, as done both in [START_REF] Liberti | Optimal configuration of gamma ray machine radiosurgery units: the sphere covering subproblem[END_REF] and [START_REF] Camino | A greedy approach combined with graph coloring for non-uniform beam layouts under antenna constraints in multibeam satellite systems[END_REF]. Although, our ambition in this work was to preserve this continuity of the position variables so this type of discretization has been discarded. Then, in the context of wireless sensor location, the author of [START_REF] Gentile | Sensor location through linear programming with triangle inequality constraints[END_REF] approximates the Euclidean distances by the L 1 norm and exploits the triangle inequality for reaching a linear programming model. This last technique falls within a more general wave of research in the field of digital distances on how to substitute cleverly the computationally expensive L 2 norm calculations by optimized combinations of the less operation demanding L 1 and L ∞ norms. Note also that these two norms are much more naturally linearizable norms than the Euclidean norm. In the most recent works ([3], [START_REF] Mukherjee | Linear combination of norms in improving approximation of euclidean norm[END_REF]), they showed that in the Euclidean plane, the best maximum relative error (obtained empirically) of an optimized combination of the L 1 and L ∞ norms with respect to the L 2 norm is approximately equal to 5.6%. Here, we are looking for a totally controllable approximation of the Euclidean norm whose maximum relative error could tend towards 0 if one was willing to pay the price in terms of numerical complexity.

Euclidean norm linearization through plane directions discretization

To find such a convenient linearization process, we relied on two geometrical results. They are both based on a parametrically controlled discretization of the directions of the Euclidean plane, that are otherwise characterized by the continuous domain [0, 2π[. In practice, for a given n directions ∈ N such that n directions ≥ 3, and for all i ∈ {1, • • • , n directions } = U, let us denote by

U i = U i,x U i,y =      cos 2(i -1)π n directions sin 2(i -1)π n directions      ∈ R 2 (1)
the i th discretized direction (running notation throughout the paper). These n directions -th roots of unity provide a regular discretization of the Euclidean plane directions with each resulting direction representing an exclusive sub-interval of [0, 2π[ of size 2π n directions . See for instance Fig. 1(a) for an example with 8 directions. Note that by definition, we

∀i ∈ U, ||U i || = 1 (2)
In practice, the Proposition 1 described below allows to define a process based on linear operations to check whether two points u, v ∈ R 2 are closer than a given distance: one has to check that the projections of the uv vector on the U i directions are all lower than a precise threshold. The fact that u and v will be decision variables (position variables) while the U i directions will be input data is to be kept in mind to understand the linearity of the process proposed.

Proposition 1: Let u, v ∈ R 2 and let d ∈ R + , [ ∀i ∈ U, u -v | U i ≤ d ] =⇒ ||u -v|| ≤ d cos (θ max ) (3)
where

θ max = π n directions Proof: Let us therefore assume that ∀i ∈ U, u -v | U i ≤ d (4)
Since two consecutive directions are separated by an angle of exactly 2π n directions , we necessarily have

|(u -v , U i min )| ≤ π n directions = θ max ( 5
)
where the direction i min ∈ U is defined as the closest direction to uv in terms of angular separation. Therefore, since ||U i min || = 1 and since θ max ∈ [0, π 3 ]:

cos (θ max ) ||u -v|| ≤ u -v , U i min ( 6
)
and since i min ∈ U and therefore verifies the equation ( 4):

||u -v|| ≤ d cos (θ max ) (7)
⊓ ⊔

As a natural complement to the previous proposition, Proposition 2 defines a linear process to check this time whether the two points u, v ∈ R 2 are sufficiently separated, according to the separation distance d.

Proposition 2:

Let once again u, v ∈ R 2 and d ∈ R + . Then, the following implication holds

[ ∃i ∈ U, u -v | U i ≥ d ] =⇒ ||u -v|| ≥ d (8)
On a practical point of view, when we need to make sure that two points are sufficiently seperated, it means that we only need to find, among the n directions discretized directions, one direction for which this inner product is sufficiently high.

Proof:

It is a direct consequence of the Cauchy-Schwarz inequality for the canonical inner product of R 2 . Let i be the direction such that uv | U i ≥ d, then:

||u -v|| = ||u -v|| • ||U i || (since ||U i || = 1) (9) ≥ | u -v | U i | (10) ≥ d (11) ⊓ ⊔

Impact of the linear approximation

In the case of Proposition 1 (the same analysis could be conducted for Proposition 2), say we are trying to check whether the distance

d uv ∈ R + between two points u, v ∈ R 2 is lower than ∆ ∈ R + .
Relying on our previous results, if all the inner products are less or equal to ∆ as dictated by Proposition 1, we will consider here that d uv ≤ ∆. Though, note that the guarantee resulting from all the inner product inequalities is that the distance between the two points is lower than ∆ lim with

∆ lim = ∆ cos (θ max ) (12)
This means that the two points could be actually at a distance comprised between ∆ and ∆ lim , and still be considered as closer than ∆ according to the process of Proposition

1. This is exactly what is represented in Fig. 1(b): D is a disk of radius ∆ centered on a certain v ∈ R 2
, and P is the set of points verifying all the inner product inequalities with respect to the point v, i.e.

P = {u ∈ R 2 / ∀i ∈ U, u -v | U i ≤ ∆}.
Therefore, P \D is what we could call the exterior approximation set, that is the set of points that are considered at a distance from v less than ∆ although they are not. Note that there is another way of exploiting Proposition 1 by comparing all the inner products to cos (θ max ) ∆ instead of directly ∆, then when all the inequalities are true, we have this time the guarantee that the distance between the two analyzed points is less or equal to ∆. However, in that case, it is possible to find situations where the distance between the two points is comprised between cos (θ max ) ∆ and ∆ and where one of the inner products has a value greater than cos (θ max ) ∆, leading to an impossibility to conclude that the two points are closer than ∆ with the process of Proposition 1. This defines the inte-

rior approximation set D\P ′ with P ′ = {u ∈ R 2 / ∀i ∈ U, u -v | U i ≤ cos (θ max ) ∆},
also represented in Fig. 1(b). In the end, we have to chose between two undesirable consequences of our linear approximation: accepting incorrect close points, or not detecting correct close points. Concerning these two effects, note that the two right-hand sides analyzed here (cos (θ max ) ∆ and ∆ that helped define P ′ and P respectively) correspond to extreme situations. Depending on the application considered, one could try to find a convenient trade-off between the two detrimental effects by chosing a right-hand side α ∈ [cos (θ max ) ∆, ∆]. Concerning the amplitude of the error caused by the linear approximation, it is directly linked to the number of directions n directions : the error tends relatively fast towards 0 when n directions increases. One way to quantify this convergence is to compare the area of D,P and P ′ with a varying number of directions, as done in Fig. 2.

Extension of these principles to R 3 and higher dimensions

Another way to interprete this choice we made to rely on the n th roots of unity in R and it indeed admits the uniform distribution defined by the n th roots of unity (A = U) as an optimal solution of optimal value J (U) = θ max (the entire set of optimal solutions can be obtained by rotating the n th roots of unity of an angle α ∈ 0, 2π n directions , each value of α leading to a distinct solution).

While it is trivial in R 2 to discretize uniformly the directions of the plane for a certain number of aimed directions n directions , it is far from obvious to find a pre-determined number of uniformly distributed points at the surface of the unit sphere in R 3 . For some values of n directions , it can even be proven that there is no solution of exactly uniform distributions of the points on the sphere. As a result, this very simple problem motivated a dedicated wave of research and therefore offers a rich literature on the different methods developped to solve it: [START_REF] Mj Peake | The equal spacing of n points on a sphere with application to partition-of-unity wave diffraction problems[END_REF], [START_REF] Edward | Distributing many points on a sphere[END_REF], [START_REF] Harman | On decompositional algorithms for uniform sampling from n-spheres and n-balls[END_REF], [START_REF] Erber | Equilibrium configurations of n equal charges on a sphere[END_REF]. For our linearization process, it means that we have to extend our principles to non-uniform discretized directions in R 3 . One way to do so is to solve at best the problem defined previously in R 3 : min

A⊂{v∈R 3 | ||v||=1} s.t. card(A)=n directions J (A) = max u∈R 3 s.t. ||u||=1 min u ′ ∈A ( u , u ′ ) (14)
In some of the aforementioned literature, this problem is exactly the one tackled, but there are articles that also address close variants of the problem with other criteria inspired by physical phenomena, such as electrostatic equilibrium for instance. From the point of view of our application, solving this optimization problem can be interpreted as finding the most isotropic treatment of all the R 3 directions before approximating linearly the Euclidean distances. Let therefore U be an optimal (or suboptimal) solution of ( 14), then Proposition 1 and Proposition 2 become valid in R 3 simply by using U as the set of discretized directions and by setting θ max = J (U). Since the application that motivated this work (beam layout optimization, detailed further in the article) is set in the Euclidean plane, we did not perform at this point any analysis for R 3 similar to the one presented in the previous paragraph for R 2 where the evolution of the approximation error with the number of discretized direction has been properly quantified. To do so, one would simply have to implement a solution for solving [START_REF] Muller | A note on a method for generating points uniformly on n-dimensional spheres[END_REF] for each number of directions n directions tested.

Finally, note that several articles tackle the question of sampling uniformly n-dimensional spheres (for instance [START_REF] Cook | Rational formulae for the production of a spherically symmetric probability distribution[END_REF] and [START_REF] Muller | A note on a method for generating points uniformly on n-dimensional spheres[END_REF]), allowing us to further extend our principles to dimensions even higher than R 3 by applying the exact same reasoning.

A simple application: the continuous k-center problem

To illustrate the applicability of the proposed linearization process, this section introduces how it could be used to solve a well-known operations research problem that mixes continuous and discrete aspects: the continuous k-center problem. It consists in defining the position of K ≥ 1 centers that are used to cover N ≥ 1 cities of known

positions C c = (X c ,Y c ) ∈ R 2 (c ∈ {1, • • • , N})
in order to provide a certain service: these centers can be fire stations, hospitals, police stations, warehouses... The goal is to place the centers and to allocate the cities to the centers in such a way as to minimize the maximum time needed to provide service to a city. To produce such solutions, most of the literature on the continuous k-center problem proposes to minimize the maximum Euclidean distance between a center and its allocated stations, but other norms can be used, as it is done in [START_REF] Wolgang | The k-center problem[END_REF] with l 1 and l ∞ norms for instance. This is naturally the Euclidean norm we considered here in order to apply the linearization principles presented in the previous section. We use continuous variables for the position of the

K centers ((x k , y k ) ∈ R 2 , k ∈ {1, • • • , K}), that are allowed to vary in a certain bound- ing box B = [X min , X max ] × [Y min , Y max ] ⊂ R 2 ,
and boolean variables to materialize the allocation of the cities to the centers (α c,k ∈ {0, 1}). Finally, λ ∈ R + is the continuous variable that will represent the maximum distance of a city to its center (that is to be minimized).

Minimize λ (15)

under the following constraints

∀k ∈ {1, • • • , K}, x k ≥ X min (16) ∀k ∈ {1, • • • , K}, x k ≤ X max (17) ∀k ∈ {1, • • • , K}, y k ≥ Y min (18) ∀k ∈ {1, • • • , K}, y k ≤ Y max (19) ∀c ∈ {1, • • • , N}, ∑ k∈{1,••• ,K} α c,k = 1 (20) ∀c ∈ {1, • • • , N}, ∀k ∈ {1, • • • , K}, ∀u ∈ {1, • • • , n directions }, cos(θ max )λ ≥ x k -X c y k -Y c • U u,x U u,y -max Z∈B ||Z -C c ||(1 -α c,k ) (21)
Equation ( 16) corresponds to the minimization of the maximum distance between a center and its associated cities. Equations ( 16), ( 17), ( 18), ( 19) define the boundaries of the center positions. The constraints [START_REF] Gl Xue | A polynomial time dual algorithm for the euclidean multifacility location problem[END_REF] force each city to be allocated to one and only center. Finally, the constraints ( 21) allows to lower-bound the continuous distance variable λ by all the center-city distances of all active center-city couples. Note that in equations ( 21), Proposition 1 is used with a "conservative" approach: the unit disk is approximated by P ′ in order to over-estimate the minimized maximum distance, instead of under-estimating it as it would have been the case with an approximation of the disk by P . Fig. 3 provides examples of solutions generated by this MILP respectively for (N = 59, K = 9, n directions = 12) and (N = 144, K = 30, n directions = 12). Note that in this problem, there are no strict proximity constraints, nor separation constraints, but simply a notion of maximum distance to be minimized. This allows to have a wide range of different strategies based on branch-and-bound procedures (see [START_REF] Hatem | A mixed breadth-depth first strategy for the branch and bound tree of euclidean k-center problems[END_REF] for instance) or on metaheuristics ([16]) to produce optimal solutions for very large instances that we could not solve with the MILP model above. However, this example seemed simple enough to illustrate how easily our linearization process allows to reach a direct algorithmic solution that relies on all the powerful solving principles inherent to Mixed Integer Linear Programming. Most importantly, such models are adaptable to variants of the problem that are characterized by added sets of constraints and variables.

4 A more complex application: the beam layout optimization in multibeam satellite systems A well-known reference system in this application (see [START_REF] Maral | Satellite communications systems: systems, techniques and technology[END_REF] for instance) is the satellitecentred (x, y, z) coordinate system, presented on Fig. 4(a). The z axis is in the satellite-Earth centre direction, the x axis is perpendicular to the meridian plane of the satellite (defined by the North, the z axis, and the position of the satellite) and is oriented toward the east, and the y axis is perpendicular to the equatorial plane and oriented in such a way as to complete a right-handed coordinate system (i.e. the south for a geostationary satellite). On the figure, S represents the satellite, O the centre of the Earth, P a point on the Earth's surface and P its projection on the equatorial plane. Then, a useful property is that there exists a subset Z of Θ × Φ and a bijection from A into Z that allows us to use the couple of true view angles (θ, ϕ) ∈ Z, as they are called, to completely identify any point P in A, as shown in Fig. 4(a). In this figure, there is also an illustration of the notion of angular separation of two points on the surface of the Earth from the point of view of the satellite, which is crucial for the problem we addressed in this paper: there is indeed an example of such an angular separation with θ, which is the angle between the sub-satellite point (on the surface of the Earth) and the point P.

There exists a very convenient representation of these true view angles in the xy plane.

To the two true view angles (θ, ϕ) of a given point P in A, we associate the so called projected true view angles θ x ∈ Θ x and θ y ∈ Θ y defined as follows:

θ x = θ cos ϕ (22) θ y = θ sin ϕ (23)
which bijectively defines Θ x and Θ y from Z, as represented in figure (b) of Fig. 4. A well known result on these projected true view angles is that, for any two points P 1 and P 2 in A, the following approximation

( -→ SP 1 , -→ SP 2 ) ≃ (θ x,P 2 -θ x,P 1 ) 2 + (θ y,P 2 -θ y,P 1 ) 2 (24)
is perfectly acceptable in the case of geostationary satellites (see [START_REF] Maral | Satellite communications systems: systems, techniques and technology[END_REF] for instance for more details on this point). This means that the angular distance from the point of view of the satellite between two points on the surface of the Earth can be computed with a simple Euclidean norm in the projection space Θ x × Θ y . For this very reason, these coordinates have been chosen for our study and our models.

Definition of the problem: variables, constraints and objective

A multibeam satellite is particular type of telecommunication satellite that provides service to its users thanks to a plurality of relatively narrow beams, a beam being a zone of significant electromagnetic power on the surface of the Earth for a given radiofrequency source. After receiving the signals from a gateway connected to the terrestrial network, the satellite payload converts in frequency, amplifies, and retransmits the input signals in the different beams through the reflector antennas, as depicted in Fig. 5 where we have 13 beams transmitted by 4 reflector antennas (one per color in the figure). In the coordinate system presented in 4.1, and with the antenna technology For antenna feasibility reasons detailed in [START_REF] Camino | A greedy approach combined with graph coloring for non-uniform beam layouts under antenna constraints in multibeam satellite systems[END_REF], two beams associated to the same reflector must have sufficiently separated beam centers, which makes the link with our previous work on Euclidean separation constraints. The rule adopted in this study is that this separation distance is proportional to the mean of two beam diameters, according to a proportionality coefficient κ ∈ R + (physically realistic values varying in [ 3 2 , √ 3]). Finally, the objective is to maximize the covered traffic ∑ s∈{covered stations} T s .

NP-hardness of the beam layout optimization

The Circle Covering Problem (CCP) defined as the following decision problem: CCP can be reduced polynomially to BLDP. The beam layout optimization problem is therefore NP-hard.

Let x 1 , • • • , x n ∈ R 2 be

Proof:

Let n and p be two positive integers defining an instance I of CCP. We are looking for an instance I ′ of BLDP polynomially defined from I and such that I true ⇔ I ′ true. Let therefore I ′ be defined as follows:

N S = n, T s = 1 for all s ∈ S, ζ = n, N B = p, N R = N B , κ = 0, N W = 1 and W 1 = 2.
Finally, for all s ∈ S, the position of the station s in the plane coincides with the s th point of I. This set of input parameters for I ′ disables several constraints and the problem becomes: can all the N S stations be covered by N B beams of radius 1 (no further constraints) ? Thus, the equivalency with I is clear. ⊓ ⊔

Mixed Integer Linear Programming model

Relying fully on the principles of section 2 allows to reach the following MILP model:

Maximize ∑ (s,b)∈S ×B T s α s,b (25)
under the following constraints for all w ∈ W , r ∈ R , s ∈ S) and corresponding "at most one" and "exactly one" constraints: (26), ( 27) and (28). They help write linearly the objective (25): traffic covered maximization. The proximity constraint that says that a station must be inside the disk of a beam to be covered by it is expressed in (29) thanks to Proposition 1 (approximation by P ′ ).

∀s ∈ S, ∑ b∈B α s,b ≤ 1 (26) ∀b ∈ B, ∑ w∈W ω b,w = 1 (27) ∀b ∈ B, ∑ r∈R ρ b,r = 1 (28) ∀s ∈ S,∀b ∈ B,∀u ∈ U, x b -X stations,s y b -Y stations,s T U u,x U u,y ≤ 1 2 cos (θ max ) ∑ w∈W W w ω b,w + (1 -α s,b )M s (29) ∀b, b ′ ∈ B such that b ′ > b, ∀r ∈ R , β b,b ′ + ρ b,r + ρ b ′ ,r ≤ 2 (30) ∀b, b ′ ∈ B such that b ′ > b, β b,b ′ + ∑ u∈U γ b,b ′ ,u ≥ 1 (31) ∀b, b ′ ∈ B such that b ′ > b, ∀u ∈ U, x b ′ -x b y b ′ -y b T U u,x U u,y ≥ κ 2 ∑ w∈W W w ω b,w + ∑ w∈W W w ω b ′ ,w -N(1 -γ b,b ′ ,u ) (32)

Experiments

The experiments were conducted on instances consisting of N S ∈ {100, 200, 300} user stations of fixed position, with the traffic demand distribution T s being generated ran-domly: each draw defines an instance. The number of beams N B ∈ {10, 20, 30} grows with N S . The number of reflectors has been set to N R = 4, the number of diameters to

N W = 2 (W 1 = 0.3 • and W 2 = 0.5 • ), κ = √ 3.
We tested different values of n directions to assess its impact on the numerical complexity: 10 instances were generated per number of user stations N S considered, the resulting 30 instances being all tested on the numbers of discretized directions n directions in the set 3; 50 . The MILP solver used is Gurobi with a timeout per instance set to 180 seconds. Some other minor industrial constraints were integrated to the model but we chose not to discuss them here. The results are given in Fig. 6(a) in the form of relative gaps between best solution found and best known objective bound for the three types of instances tested. Each point of the three curves is an average gap value obtained over the 10 instances of the corresponding category of instances. As we could have predicted, the more beams and stations, the harder the convergence toward the optimal solution, materialized by higher average gaps in Fig. 6(a). Then, the main observation that can be made from these results is that, for a too low number of directions (3 ≤ n directions ≤ 8), the approximation of the disks by P (separation constraint) and P ′ (proximity constraint) is too rough and does not allow to reach solutions of good quality. On the other hand, for a too high number of directions (20 ≤ n directions ≤ 50), the gain in approximation accuracy becomes so small that the solution quality improvements, if any, do not compensate the increase in numerical complexity due to the growing model size: this explains the degradation of the average gaps in Fig. 6(a) for these values of n directions . This is a general rule to keep in mind for applying this Euclidean norm linearization technique. For illustra-

Conclusion

We introduced a general methodology for linearizing inequality constraints involving the L 2 norm, and validated it on a crucial industrial problem for satellite manufacturers, for which we have reached significantly improved solutions compared to the existing hand-crafted ones. The next steps, that are already a work in progess, will consist in improving the current model to remove symmetries and solve bigger instances, for finally benchmarking the methodology with MINLP solvers.

Fig. 1 .

 1 Fig. 1. (a) Discretization of the directions of the Euclidean plane (n directions = 8) : the 8 th roots of unity (b) Approximation of the Euclidean plane disk D by the regular n directions -sided polygons P and P ′ with the linear approximation of the Euclidean norm

Fig. 2 .

 2 Fig. 2. (a) Evolution of the area of the approximating polygons P and P ′ with respect to the area of the disk D (b) Examples of values appearing in the curves

4. 1

 1 Coordinate system used to define the optimization problemThis paragraph describes the coordinate system used to identify the points on the Earth's surface in the context of satellite communications. It is a necessary information to understand the beam layout optimization problem as it is presented in the next paragraph.

Fig. 3 .

 3 Fig. 3. (a) k-center solution example for N = 59, K = 9, n directions = 12(b) k-center solution example for N = 144, K = 30, n directions = 12

Fig. 4 .

 4 Fig. 4. (a) Satellite-centered coordinate system and true view angles (b) True view representation

Fig. 5 .

 5 Fig. 5. Standard architecture of a multibeam satellite system

  Variables: α s,b , ω b,w , ρ b,r , β b,b ′ , γ b,b ′ ,u ∈ {0, 1}, x b , y b ∈ R For each beam b ∈ B, we introduce beam center variables (x b , y b ) ∈ R 2 , beam diameter, satellite reflector and station allocation variables (respectively ω b,w , ρ b,r , α s,b ∈ {0, 1}

  M s ∈ R + (precisely tuned) relaxes the constraints when b ∈ B does not cover s ∈ S (α s,b = 0). Constraints (30) force the β b,b ′ ∈ {0, 1} variables to be equal to 0 if b, b ′ ∈ B (b = b ′ ) use the same reflector. The antenna separation constraints are activated in that case, first through constraints (31) that force at least one γ b,b ′ ,u ∈ {0, 1} (u ∈ U) to be equal to 1: this is the materialization of the existential quantifier of Proposition 2. Constraints (32) operate the separation according to Proposition 2 and coefficient κ when γ b,b ′ ,u = 1 (approximation by P ). They are relaxed by N ∈ R + when γ b,b ′ ,u = 0.

Fig. 6 .

 6 Fig. 6. (a) Campaign of runs for varying numbers of discretized directions and varying numbers of user stations (b) Example of instance with 200 user stations solved optimally with 12 discretized directions

  2 to discretize the Euclidean plane directions is to observe that they are a solution to the following problem: what subset of size n directions of the R 2 unit circle minimizes the maximum angular distance between a point of the unit circle and its angularly closest point among the selected n directions points ? Mathematically, this problem can be ex-

							Area(P)
		1.6					Area(D) Area(P ′ )
							Area(D)
		1.4					n directions	Area(P) Area(D)	Area(P ′ ) Area(D)
	Relative areas	1 1.2 0.8					4 8 12 16 20	1.273 1.055 1.023 1.013 1.008	0.637 0.900 0.955 0.974 0.984
		0.6					24 28	1.006 1.004	0.989 0.992
		0.4					(b)
		5	10	15	20	25	30
			Number of discretized directions
				(a)
	pressed as follows
		min A⊂{v∈R 2 | ||v||=1} s.t. card(A)=n directions	J (A) = max u∈R 2 s.t. ||u||=1	min u ′ ∈A	( u , u ′ )	(13)

  n points of the Euclidean plane and let p ∈ N * . Can we find p unit disks D 1 , • • • , D p ⊂ R 2 such that each point is contained by at least one disk ? In [12], the authors prove that the Circle Covering Problem (CCP) is NP-complete. The beam layout decision problem (BLDP) associated to the beam layout optimization problem defined in the previous paragraph is the following Let ζ ∈ R + . Can the aggregate covered traffic ∑ T s be greater than or equal to ζ under the constraints of section 4.2 ?

	Proposition 3: