Augustin Parret-Fréaud

Pierre Gosselet
email: gosselet@lmt.ens-cachan.fr

Christian Rey
email: rey@lmt.ens-cachan.fr

Error estimation for sub-structured problems

published or not. The documents may come

Error estimation for substructured problems

INTRODUCTION

Virtual testing has become a short term aim of industrials which will to replace expensive experimental studies and validations by numerical simulations, even in order to certify large structures as planes and bridges. So, one keypoint of the new numerical methods to develop is the verication of computations which enables to warranty that the computed solution is suciently close to the original continuum mechanics model. This topic of numerical analysis has been the subject of many studies for the last decades.

Another key point of new numerical methods is their ability to quickly provide solutions to large (nonlinear) systems. The most classical answer to this issue is to use domain decomposition methods in order to take advantage of the parallel hardware architecture of recent clusters and grids. In engineering, non-overlapping domain decomposition methods are mostly employed.

Our aim is to provide a fully integrated adaptative computational strategy to compute large structural mechanical problems with certied quality. Our method is based on the error in constitutive relation [START_REF] Ladevèze | Mastering calculations in linear and nonlinear mechanics[END_REF] to measure the quality of our results, and on a generic vision of non-overlapping domain decomposition methods [START_REF] Gosselet | Non-overlapping domain decomposition methods in structural mechanics[END_REF].

1

In this paper, we focus especially on the computation of elementary contribution to the global error in a fully parallel context. First, we give basic notions about non-overlapping domain decomposition and constitutive error for elastic mechanical problem. Then, the present a strategy to compute elementary contribution tted to sub-structuration. At last, we give some results and conclusions.

FORMULATION AND BASIC NOTIONS

We study the static equilibrium of a linear elastic structure occupying the domain Ω under the small perturbation hypothesis. The problem is solved using nite element discretization.

Non-overlapping domain decomposition

Domain Ω is decomposed into a set of subdomains (Ω (s) ), in the following superscript (s) stands for data associated to domain Ω (s) . Assuming conforming decomposition of the structure (one-to-one correspondance of nodes on the interface between subdomains), the equilibrium equation Ku = f (K is the stiness matrix, u and f the displacement and force nodal elds), can be written in a substructured-fashion:

K (s) u (s) = f (s) + t (s) T λ (s) b ∑ s A (s) u (s) b = 0 ∑ s A (s) λ (s) b = 0 (1)
First equation of (1) stands for equilibrium of substructure (s), λ

(s) b
is the nodal reaction of neighbouring subdomains and t (s) is the trace operator (which extracts the boundary nodes from the subdomain nodes). Second equation of (1) sets the continuity of displacements between substructures (u (s) b = t (s) u (s) is the boundary displacement, and A (s) is a signed boolean assembly operator like in FETI method [START_REF] Farhat | The dual schur complement method with well-posed local neumann problems[END_REF]). Third equation of (1) sets the equilibrium of reactions (action-reaction principle) (A (s) is a boolean assembly operator like in BDD method [START_REF] Mandel | Balancing domain decomposition[END_REF]).

Whatever the (most often iterative) resolution strategy to system (1), in the end one obtains boundary elds u

u KA : u |∂1Ω = u d σ SA : ∫ Ω σ : ε(u * )dΩ = ∫ Ω f d • u * dΩ + ∫ ∂2Ω F d • u * dS, ∀u * KA 0 CR : σ = K : ε(u) ( 2
)
where σ is Cauchy stress tensor, ε(u)

e CR ( u h , σ h ) = σ h -K : ε( u h ) σ,Ω with x 2 σ,Ω = ∫ Ω ( x : K -1 : x ) dΩ (3)
An important feature is that whenever (u h , σ h ) satises (2, KA&SA), the associated estimator overestimates the true error.

RECOVERING SA STRESS WITHIN SUB-STRUCTURED CONTEXT

Most classical nite element formulations provide displacement eld u h satisfying kinematic constraints so that u h is chosen equal to u h . However, σ h is seldom statically admissible so that σ h satisfying (2, SA) has to be computed from σ h . A systematic two steps process (illustrated g. 1(a)) consists in: rst, building equilibrated force densities F set on the edges of each element of the mesh, using star-patches of elements and previous eld σ h ; second, recovering equilibrated stress eld σ h from an element-by-element procedure with F applied as Neumann boundary condition, using when possible analytical solution or richer approximation basis ( [START_REF] Babuska | Validation of a posteriori error estimators by numerical approach[END_REF] recommands p + 3 interpolation).

Since the second step is fully local, it can be easily driven in a sub-structured context. Unfortunatly, the rst step involves, for each node, systematic computations on a patch of elements centered on it (cf g. 1(b)). In a substructured # sub-domains 1 2 3 5 8 9

e CR (%) 9.40 9.31 13.59 9.23 9.32 8.77

Table 1: Global error for dierent substructurations These rst results enable to conclude that in most cases our strategy provides satisfying SA elds since the estimated error is smaller when the problem is substructured. The 3-subdomain substructuration is an exception, in this case, an interface between subdomain passes through the stress concentration zone making the linear interpolation of reaction density irrelevant.

Future work will try to overcome this problem, so that no a priori knowledge of stress concentrations will be necessary to provide decompositions leading to correct fully-parallel estimation of constitutive relation error.

2. 2

 2 Constitutive relation error in elasticityA posteriori error estimation, which has been widely studied for the last decades, contains a large number of techniques. The associated estimators are usually classied into three main categories: ux projection techniques, residual on equilibrium equation or error in constitutive relation[START_REF] Ladevèze | Mastering calculations in linear and nonlinear mechanics[END_REF].In the present study, we have chosen to use an estimator based on the constitutive relation error whose principle is recalled below. The usual continuum formulation of any mechanical problem can be rewritten in term of kinematic admissibility (2, KA) , static admissibility (2, SA) through the principle of virtual work and constitutive relation (2, CR):

Figure 1 :

 1 Figure 1: Recovering of SA stress

  is the symmetric part of the gradient of displacement eld u, u d is given displacement on the part ∂ 1 Ω of the boundary, F d is given surface eort on the part ∂ 2 Ω of the boundary, f d is the given volume force.

Usually, a nite-element solution pair (u h , σ h ) satises (2, KA&CR) but not (2, SA). Then, the starting point of the measure of the constitutive relation error is to build a new displacement-stress pair ( u h , σ h ) from (u h , σ h ), which satises (2, KA&SA). Thereafter, the constitutive relation error is computed through the residual on the constitutive relation (2, CR) with the use of following energetic norm (3):


Pablo Adasme
email: pablo.adasme@usach.cl

Abdel Lisser
email: abdel.lisser@centralesupelec.fr

A Stochastic Geometric Programming Approach for Power Allocation in Wireless Networks

Keywords: Geometric and signomial programming, sequential and piece-wise linear approximations, wireless networks, power allocation

In this paper, we consider the power allocation problem for 5G wireless networks using massive multiple input multiple output technologies. Two non-linear optimization models are considered to maximize the worst user signal-to-interference noise ratio and the total capacity of the network subject to power constraints. In particular, we transform the first one into a geometric programming (GP) problem. Whereas the second one leads to a signomial programming formulation. The main contributions of the paper are first to propose novel formulations for power allocation in wireless networks while using stochastic, geometric, and signomial programming frameworks altogether. We derive stochastic formulations for each GP model to deal with the uncertainty of wireless channels. Secondly, since solving optimally the stochastic models represents a challenging task, we obtain tight bounds using approximation solution methods. In particular, the piece-wise linear programming and the sequential approximation methods allow us to obtain tight intervals for the objective function values of the stochastic models. Notice that these intervals contain the optimal solutions. In particular, we propose an approximated GP model that allows obtaining lower bounds for the signomial problem. This is achieved by using the arithmetic-geometric mean inequality. Finally, we compare the deterministic and stochastic models and prove the robustness of the stochastic models. Notice that we solve all the instances and obtain near-optimal solutions for most of them.

Introduction

The fifth-generation (5G) of wireless communications brought increased attention from both academia and industry communities as it offers the potential development of unprecedented future network applications. It is expected that these networks will impact the operation of many industries and society. Example applications include mobile health, autonomous vehicles, smart cities and homes, manufacturing and entertainment, education, smart grid, data analytics, and networks to be developed under the Internet of Things (IoT) paradigm, just to name a few [START_REF] Chiasserini | IEEE 5G for the Automotive Domain[END_REF][START_REF] Somisetty | Big Data Analytics in 5G[END_REF][START_REF] Markakis | [END_REF][START_REF] Orsino | Ericsson Research, Factories of the Future Enabled by 5G Technology[END_REF][START_REF] Rao | The Path to 5G for Health Care[END_REF][START_REF] Mangra | Smart Cities: Connected Ecosystem of Ecosystems[END_REF][START_REF]IEEE 5G and Beyond Technology Roadmap[END_REF]. Notice that 5G technology, as a revolutionary approach, is envisioned to eliminate access bounds to wireless networks, limitations of bandwidth capacity and latency on connectivity [START_REF]IEEE 5G and Beyond Technology Roadmap[END_REF]. Unfortunately, current network infrastructures cannot be adapted straightforwardly to 5G technology. Thus, it is mandatorily required to update previous ones. These new updates will play a critical role since novel components of 5G networks will handle the huge increased data usage, coverage, security, and with low latency. In this paper, we consider the problem of optimal power allocation on these networks by using Massive Multiple Input Multiple Output (MaMIMO) technologies. Notice that MaMIMO appears as an extension of a traditional MIMO system. The main advantage of MaMIMO is that we can use a significantly larger number of antennas than in a classical MIMO system. It has been shown that MaMIMO allows one to achieve better performance in terms of signal paths, link reliability, coverage, and security. As such, it emerges as a serious candidate for future 5G based-networks.

Our main contributions in this paper can be enumerated as follows. First, we propose novel formulations for power allocation in wireless networks while using stochastic, geometric, and signomial programming frameworks altogether. Secondly, since solving optimally the stochastic models represents a challenging task, we obtain tight bounds (intervals) using approximation solution methods. In particular, the piece-wise linear programming and the sequential approximation methods allow us to obtain tight intervals for the objective function values of the stochastic models. Notice that these intervals contain the optimal solutions. Next, we propose an approximated GP model to obtain lower bounds for the signomial model. For this purpose, we use the arithmetic-geometric mean inequality. Finally, we compare the deterministic and stochastic models and prove the robustness of the stochastic models. Notice that we solve all the instances and obtain near-optimal solutions for most of the instances. More precisely, we consider two non-linear programming optimization problems for the optimal power allocation where the aim is to maximize the worst user signal-to-interference noise ratio (SINR) and the total capacity of the network subject to power constraints. In particular, we transform the first one into equivalent geometric programming (GP) problem. Whereas the second one leads to an equivalent signomial programming (SP) formulation. Notice that a GP model is a non-convex problem. However, it can be transformed into a convex one by using standard logarithmic transformations leading to sum-log-exp convex functions. Consequently, all GPs can be solved to global optimality by using interior point algorithms in polynomial time complexity [START_REF] Boyd | A Tutorial on Geometric Programming[END_REF]. Unlike a GP problem, an SP problem cannot be transformed into an equivalent convex one. Thus it is significantly harder to solve optimally. However, in this case, we also derive an approximated GP formulation that allows obtaining lower bounds for the original SP problem. The latter is achieved by using the arithmetic-geometric mean inequality [START_REF] Boyd | A Tutorial on Geometric Programming[END_REF][START_REF] Duffin | Geometric programming-theory and application[END_REF]. Subsequently, we derive stochastic formulations for each GP model to deal with the uncertainty of wireless channels. In particular, we propose individual and joint chance constraints for each GP model. Finally, we obtain lower and upper bounds with piece-wise tangent linear and sequential convex approximation methods. We conduct substantial numerical experiments to compare all the proposed models and algorithms.

A GP model is a type of mathematical optimization problem characterized by objective and constraint functions that have a special form [START_REF] Boyd | A Tutorial on Geometric Programming[END_REF]. In general, a GP model can be written as min

{x∈R m ++ } i∈I0 c i m j=1 x dij j (1)
s.t.

i∈I k c i m j=1 x dij j ≤ 1, k ∈ {1, . . . , K} (2)
where the set {I k , k ∈ K = {0, 1, . . . , K}} is the disjoint index sets of {1, . . . , Q}. Usually, the term c i m j=1 x dij j

for each i ∈ I k is called a monomial function where each term c i for all i ∈ I k must be a nonnegative real number. If at least one of them is negative, then the problem (1)-( 2) is no longer a GP model, but an SP problem. The exponent parameters d ij for all i ∈ I k , k ∈ K and j ∈ M = {1, . . . , m} can be real numbers. Finally, each variable x j for j ∈ M must be strictly positive. Hereafter, we denote by R m ++ the set of positive real numbers of dimension m. A sum of monomials is called a posynomial function. Constraints where a monomial is strictly equal to one, are also allowed in a GP model. However, constraints involving a posynomial function strictly equal to one also lead to a non-GP problem. We will see that the arising power allocation problems considered in this paper can be written as equivalent GP or SP problems. As far as we know, individual and joint chance constraints for GP and SP problems have not been investigated in the literature for power allocation in wireless networks. Individual and joint chance constraints for the constraints in (2) can be formulated respectively as [START_REF] Liu | Stochastic geometric optimization with joint probabilistic constraints[END_REF]

P ξ i∈I k c i (ξ) m j=1 x dij j ≤ 1 ≥ 1 -α, k ∈ {1, . . . , K} (3)
and

P ξ i∈I k c i (ξ) m j=1 x dij j ≤ 1, k ∈ {1, . . . , K} ≥ 1 -α (4)
where we assume that each input parameter (c i (ξ)) for all i ∈ I k , k ∈ K behaves as a random variable which is distributed according to a certain probability distribution function. Notice that each chance constraint (3), for each k ∈ K, ensures that at least (1 -α) percentage of the uncertain constraints must be satisfied separately where α is an input parameter chosen arbitrarily from the interval [0; 0.5). Whereas the joint chance constraint (4) ensures that (1 -α) percentage of the constraints must be satisfied according to a joint probability distribution function. This paper is organized as follows. In Section 2, we present and discuss some related works which are closer to our power allocation problem. We also discuss a few relevant works related to probabilistic individual and joint chance constraints when applied to GP problems. Next, in Section 3, we present two non-linear optimization problems for power allocation in 5G wireless networks and explain how they can be transformed into equivalent GP models. Subsequently, we propose individual and joint chance constraints for each GP model and derive its equivalent deterministic formulation. Next, in Section 4, we present and explain the sequential approximation methods used to solve the GP models involving joint probabilistic constraints and also for solving the approximated GP models that allow computing lower bounds for the SP problem. In this section, we also present GP models using tangent piece-wise linear functions that we use to calculate lower and upper bounds for the two non-linear optimization problems. Subsequently, in Section 5, we conduct substantial numerical experiments and compare all the proposed models and algorithms. Finally, Section 6 concludes the paper.

Related work

Interference is a major problem in wireless networks affecting the broadcast nature of radio transmission channels. It affects data rate metrics including user signal-to-interference noise ratio as well as the maximum achievable capacity of the networks. In particular, in communication systems where users transmit in frequency bands that are non-orthogonal, it is significantly harder to deal with the interference problem. A common strategy that allows minimizing it consists of handling efficiently the amount of power to be assigned to different users. Unfortunately, several service metrics related to the interference are non-linear functions of SINR which is a non-linear and non-convex function of power. As a consequence, the arising power control optimization problems are difficult to solve optimally since they belong to the complexity class of NPhard problems. However, some of these problems can be transformed into GP problems, and hence they can be convexified.

Some recently published papers dealing with power allocation problems in wireless networks where the geometric programming approach is utilized can be consulted for instance in [START_REF] Parsaeefard | User Association in Cloud RANs with Massive MIMO[END_REF][START_REF] Azizi | Profit Maximization in 5G+ Networks with Heterogeneous Aerial and Ground Base Stations[END_REF][START_REF] Bashar | Energy Efficiency of the Cell-Free Massive MIMO Uplink With Optimal Uniform Quantization[END_REF][START_REF] Amani | Energy-Efficient Resource Allocation in Heterogeneous Cloud Radio Access Networks via BBU Offloading[END_REF][START_REF] Tweed | Outage-Constrained Resource Allocation in Uplink NOMA for Critical Applications[END_REF][START_REF] Parsaeefard | Dynamic Resource Allocation for Virtualized Wireless Networks in Massive-MIMO-Aided and Fronthaul-Limited C-RAN[END_REF][START_REF] Lakani | Optimal Design and Energy Efficient Binary Resource Allocation of Interference-Limited Cellular Relay-Aided Systems With Consideration of Queue Stability[END_REF]. More general works can be found in references [START_REF] Gao | Antenna selection and power allocation design for 5G massive MIMO uplink networks[END_REF][START_REF] Benmimoune | Joint Transmit Antenna Selection and User Scheduling for Massive MIMO Systems[END_REF][START_REF] Pedram | Energy Efficiency in 5G Cellular Network Systems[END_REF][START_REF] Zhang | Robust Energy-Efficient Transmission for Wireless-Powered D2D Communication Networks[END_REF][START_REF] Le | Enabling User Grouping and Fixed Power Allocation Scheme for Reconfigurable Intelligent Surfaces-Aided Wireless Systems[END_REF][START_REF] Lu | Distortion-Aware Cross-Layer Power Allocation for Video Transmission Over Multi-User NOMA Systems[END_REF]. For the sake of brevity, we describe a few of them where the GP approach was taken into account for power allocation problems in wireless networks. In [START_REF] Parsaeefard | User Association in Cloud RANs with Massive MIMO[END_REF], the authors consider a resource allocation problem which is formulated as a non-convex problem that allows connecting a set of users to a cloud radio access network structure. To solve the problem, the authors propose an iterative algorithm that requires applying different transformations and convexification techniques and solving a sequence of GP problems with a sequential convex approximation method. By simulation results, they show that their proposed algorithm allows for increasing the total throughput of the network. Similarly in [START_REF] Azizi | Profit Maximization in 5G+ Networks with Heterogeneous Aerial and Ground Base Stations[END_REF], the authors propose a novel framework for 5G and beyond (5G+) heterogeneous wireless networks which require the access of a technology referred to in the literature as power domain nonorthogonal multiple access (PD-NOMA). The main goal of their work is to maximize the total network profit under some practical technological network constraints in addition to power maximum limits, and isolation of the virtualized wireless network. They formulate the problem as a mixed integer non-linear optimization problem. To solve the model, they propose a practical approach with reduced computational complexity that consists of solving an alternating method where the optimization problem is broken down into three sub-problems. The latter is achieved by using the sequential convex approximation method, the GP approach, and a mesh adaptive direct search method. Their numerical experiments reveal that the proposed approach can improve the overall network profit. Also, in [START_REF] Tweed | Outage-Constrained Resource Allocation in Uplink NOMA for Critical Applications[END_REF] the authors consider a resource allocation problem for uplink non-orthogonal multiple access networks for health and public safety applications. The authors apply a chance-constrained robust optimization approach leading to a joint resource allocation problem that allows minimizing user transmit power subject to rate and outage constraints. Since their resulting non-linear problem is non-convex, they apply variable relaxation and complementary geometric programming approaches and develop a two-step iterative algorithm based on successive convex approximations. By simulation results, the authors demonstrate that their proposed algorithm outperforms the traditional orthogonal multiple access transmission schemes in terms of user transmit power and overall system density using fewer subcarriers. Finally, we mention the work proposed by [START_REF] Chiang | Power Control By Geometric Programming[END_REF] where the authors present several power allocation models that can be equivalently written as GP problems. Their examples include maximizing the total system throughput or the worst user throughput, subject to the quality of service constraints, delay constraints on data rates, and outage probabilities.

As can be observed from the literature, several relevant problems related to future wireless networks utilized the GP approach to solve specialized resource allocation problems. However, we notice that there have been no attempts yet to include stochastic individual or joint chance constraints in these types of problems to handle efficiently the uncertainty of the input parameters. The chance-constrained approach allows imposing a probability of occurrence for some or all of the constraints of a mathematical programming problem. This means that some of the constraints will be satisfied, at least for a given percentage of the occurrences of each random variable. Probabilistic constraints can be considered either separately or jointly. Notice that the stochastic programming framework is a relevant approach in the literature and is frequently considered in optimization problems. It allows for obtaining optimal solutions while taking into account the uncertainty of the input parameters of a mathematical programming problem [START_REF] Cheng | Distributionally robust stochastic knapsack problem[END_REF][START_REF] Cheng | A second-order cone programming approach for linear programs with joint probabilistic constraints[END_REF][START_REF] Lobo | Applications of second-Order cone programming[END_REF][START_REF] Luedtke | An integer programming approach for linear programs with probabilistic constraints[END_REF][START_REF] Prékopa | Stochastic Programming[END_REF][START_REF] Prékopa | On probabilistic constrained programming[END_REF][START_REF] Shapiro | Lectures on stochastic programming: Modeling and theory[END_REF][START_REF] Adasme | Stochastic and semidefinite optimization for scheduling in orthogonal frequency division multiple access networks[END_REF][START_REF] Adasme | Uplink scheduling for joint wireless orthogonal frequency and time division multiple access networks[END_REF]. In particular, the authors in [START_REF] Dupacová | Stochastic geometric programming: approaches and applications[END_REF][START_REF] Rao | Engineering Optimization: Theory and Practice[END_REF] show that the probabilistic constraint (3) can be equivalently written employing two deterministic constraints involving posynomials and slack variables. To the best of our knowledge, there are only a few works in the literature including stochastic GP problems subject to joint probabilistic constraints. In [START_REF] Liu | Stochastic geometric optimization with joint probabilistic constraints[END_REF], the authors impose joint probabilistic constraints on a generic GP problem and assume that the input stochastic parameters are normally distributed and independent of each other. They approximate the problem by using piece-wise linear functions which allow transforming the resulting problem into an equivalent convex geometric program. They also prove that their approximation model allows obtaining lower bounds for the original problem. Finally, they design a sequential convex optimization algorithm that allows obtaining upper bounds for the original problem too. Similarly, the authors in [START_REF] Liu | Rectangular chance constrained geometric optimization[END_REF] discuss joint rectangular probabilistic constrained GPs. For this purpose, they assume that the input random parameters are elliptically distributed and pairwise independent. The authors obtain a non-convex reformulation of the joint rectangular chance-constrained problem and propose convex approximations using variable transformation and piece-wise linear functions. Finally, they provide a theoretical bound for the number of segments in the worst-case scenario and show by numerical experiments that their approximations are asymptotically tight.

A few other recently published works dealing with resource allocation of power in wireless networks where the approaches proposed in this paper could be applied as part of future research can be briefly described as follows. In [START_REF] Mohajer | Heterogeneous Computational Resource Allocation for NOMA: Toward Green Mobile Edge-Computing Systems[END_REF], the authors propose a dynamic optimization model to maximize the total uplink and/or downlink energy efficiency while satisfying quality of service constraints using as a viable solution the Mobile Edge Computing (MEC) paradigm to satisfy the growing demand for broadband for the next generation heterogeneous systems. The authors divide the optimization problem into two separate sub-problems, a computational carrier scheduling one, and a resource allocation problem. In particular, they propose a sub-gradient method for the computational resource allocation and a successive convex approximation together with a dual decomposition method. Their simulation results show considerable improvements for various traffic models guaranteeing fairness requirements. It also improved the total throughput for mobile computing services. Similarly, in [START_REF] Mohajer | Mobility-aware load Balancing for Reliable Self-Organization Networks: Multi-agent Deep Reinforcement Learning[END_REF] the authors present a deep learning-based mobility robustness optimization solution that learns the required parameter values for the mobility pattern for self-organizing networks. Their simulation experiments show that the function of mobility robustness optimization not only learns to optimize its performance but also it learns how to distribute the excess load of the network. Finally, they prove that the solution minimizes the number of unsatisfied users guaranteeing a more balanced network. Finally, in reference [START_REF] Mohajer | Energy-Aware Hierarchical Resource Management and Backhaul Traffic Optimization in Heterogeneous Cellular Networks[END_REF], the authors consider a dynamic optimization model to minimize the total energy consumption of fifth-generation (5G) heterogeneous networks while providing the required coverage and capacity. This is achieved by optimizing carrier allocation and power utilization. They do also propose a multi-hop back-hauling strategy to effectively use the existing infrastructure of small-cell networks for simultaneous transmissions. Their numerical results demonstrated considerable power savings for different traffic models. They further show that energy efficiency and system data rates can be significantly improved.

Mathematical formulations

In this section, first, we present the two non-linear optimization problems. Then, for each one of them, we derive equivalent GP formulations. In particular, for the second one, we obtain an equivalent SP formulation. Thus, we derive an approximated GP model that allows obtaining lower bounds for it. Subsequently, we impose individual and joint chance constraints on each GP model obtained and derive for each one of them an equivalent deterministic GP formulation.

Maximizing the worst user signal to interference noise ratio

The first model we consider allows maximizing the worst user signal-tointerference noise ratio for a MaMIMO system where each receiver has perfect channel state information [START_REF] Haleem | On the Capacity and Transmission Techniques of Massive MIMO Systems[END_REF]. Hereafter, we assume that the MaMIMO network is located inside a single cell area and it is composed of a set of K = {1, . . . , K} users and a unique base station (BS). We also assume that each user is using only one antenna to receive the data from the BS which is equipped with a predefined number of antennas T . We formulate the optimization problem as

(M 1 ) : max {p∈R K ++ } min {i∈K} p i | g H i g i | 2 j∈K,(j̸ =i) p j | g H i g j | 2 + | σ i | 2 (5) s.t. P min ≤ p i ≤ P max , ∀i ∈ K (6)
where p = (p i ) denotes the amount of power to be assigned for each user i ∈ K. The input parameters g i ∈ C T ×1 , g H i ∈ C 1×T and σ i denotes the beam domain channel vector associated to user i ∈ K, its Hermitian transpose and Additive White Gaussian Noise (AWGN), respectively. The AWGN is assumed to behave according to an independent complex Gaussian distribution function with zero mean and unit variance (σ i ∼ CN (0, 1)). Finally, we assume that each component of each vector g i , i ∈ K is a complex number generated according to a quasi-static independent and identically distributed Rayleigh fading channel. Notice that from the solution of model (M 1 ), a power allocation mechanism can be designed by simply using the power allocation vector obtained from the model. Similarly, the worst user SINR can be easily obtained from the value of the objective function since we use a unique variable for all users in K. This objective function value will coincide with the user or users having the maximum worst SINR.

In (M 1 ), the objective function ( 5) denotes the worst user SINR [START_REF] Adasme | Maximizing Signal to Interference Noise Ratio for Massive MIMO: A Mathematical Programming Approach[END_REF][START_REF] Valduga | Low complexity heuristics to beam selection and rate adaptation in sparse massive MIMO systems[END_REF]. Whilst the constraints (6) ensure that each power variable is greater than or equal to P min and less than or equal to P max . Notice that by introducing an additional variable t and defining the parameters

a ij =| g H i g j | 2 | g H i g i | -2 > 0 and b i =| σ 2 i || g H i g i | -2
> 0 for all i, j ∈ K, we can write (M 1 ) equivalently as min

{(p,t)∈R K+1 ++ } t -1 (7) s.t.
j∈K,(j̸ =i)

a ij p j p -1 i t + b i p -1 i t ≤ 1, ∀i ∈ K (8)
P min ≤ p i ≤ P max , ∀i ∈ K
We note that the coefficients a ij and b i in ( 8) for all i, j ∈ K, may not be known precisely. Thus, we suppose that each one of them is normally distributed and independent of the other, i.e., a ij (ξ) ∼ N (ā ij , σ 2 aij ) and b i (ξ) ∼ N ( bi , σ 2 bi ). This allows replacing the deterministic constraints in [START_REF] Boyd | A Tutorial on Geometric Programming[END_REF] with the following individual chance constraints

P ξ    j∈K,(j̸ =i) a ij (ξ)p j p -1 i t + b i (ξ)p -1 i t ≤ 1    ≥ (1 -α), ∀i ∈ K (9)
According to [START_REF] Dupacová | Stochastic geometric programming: approaches and applications[END_REF], the constraints (9) can be written in equivalent deterministic form as j∈K,(j̸ =i)

āij p j p -1 i t + bi p -1 i t + ϕ -1 (1 -α)    j∈K,(j̸ =i) σ 2 aij p 2 j p -2 i t 2 + σ 2 bi p -2 i t 2    ≤ 1, ∀i ∈ K (10)
where ϕ -1 (1 -α) is the quantile of the standard normal distribution function N (0, 1). The parameter α can take values in the interval α ∈ [0; 0.5). Notice that the constraints [START_REF] Parsaeefard | User Association in Cloud RANs with Massive MIMO[END_REF] are generalized posynomial constraints. Thus, they can be transformed into standard GP format. The latter can be achieved by introducing additional variables θ i > 0 for all i ∈ K. This allows imposing an upper bound on each term inside the root square j∈K,(j̸ =i)

σ 2 aij p 2 j p -2 i t 2 + σ 2 bi p -2 i t 2 ≤ θ 2 i , ∀i ∈ K (11)
which in turn allows obtaining the equivalent set of constraints j∈K,(j̸ =i)

āij p j p -1 i t + bi p -1 i t + ϕ -1 (1 -α)θ i ≤ 1, ∀i ∈ K (12)
j∈K,(j̸ =i)

σ 2 aij p 2 j p -2 i t 2 θ -2 i + σ 2 bi p -2 i t 2 θ -2 i ≤ 1, ∀i ∈ K (13)
Notice that the constraints ( 12)-( 13) prove that the resulting problem is a GP problem. Consequently, we can obtain the global optimal solution using individual chance constraints by solving the following GP model

(IM 1 ) : min {(p,t)∈R K+1 ++ } t -1 s.t. j∈K,(j̸ =i) āij p j p -1 i t + bi p -1 i t + ϕ -1 (1 -α)    j∈K,(j̸ =i) σ 2 aij p 2 j p -2 i t 2 + σ 2 bi p -2 i t 2    ≤ 1, ∀i ∈ K P min ≤ p i ≤ P max , ∀i ∈ K
Similarly, we can replace the deterministic constraints (8) by the following joint probabilistic constraints

P ξ    j∈K,(j̸ =i) a ij (ξ)p j p -1 i t + b i (ξ)p -1 i t ≤ 1, ∀i ∈ K    ≥ (1 -α) (14)
As these constraints are equivalent, we have

i∈K P ξ   j∈K,(j̸ =i) a ij (ξ)p j p -1 i t + b i (ξ)p -1 i t ≤ 1   ≥ (1 -α) (15)
Next, by introducing auxiliary variables y i ∈ R ++ for all i ∈ K, we get the following equivalent deterministic GP problem using joint chance constraints [START_REF] Liu | Stochastic geometric optimization with joint probabilistic constraints[END_REF][START_REF] Dupacová | Stochastic geometric programming: approaches and applications[END_REF][START_REF] Rao | Engineering Optimization: Theory and Practice[END_REF] (JM 1 ) : min

{(p,t,y)∈R 2K+1 ++ } t -1 s.t. j∈K,(j̸ =i) āij p j p -1 i t + bi p -1 i t + ϕ -1 (y i )    j∈K,(j̸ =i) σ 2 aij p 2 j p -2 i t 2 + σ 2 bi p -2 i t 2    ≤ 1, ∀i ∈ K (16) i∈K y i ≥ (1 -α) (17) 0 ≤ y i ≤ 1, ∀i ∈ K (18)
P min ≤ p i ≤ P max , ∀i ∈ K
Notice that the quantile function now depends on the variable y i for all i ∈ K.

In the next section, we obtain lower and upper bounds for (JM 1 ) using tangent piece-wise linear functions as well as a sequential approximation algorithm [START_REF] Liu | Stochastic geometric optimization with joint probabilistic constraints[END_REF]. This allows reporting an interval where the global optimal solution lies.

Maximizing the total capacity of the wireless network

The second model we consider maximizes the total capacity of the network in presence of signal-to-interference noise ratio for a MaMIMO system where each receiver has a perfect channel state information [START_REF] Haleem | On the Capacity and Transmission Techniques of Massive MIMO Systems[END_REF]. In this case, the non-linear optimization problem we consider can be stated as

(M 2 ) : max {p∈R K ++ } i∈K log 2   1 + p i | g H i g i | 2 j∈K (j̸ =i) p j | g H i g j | 2 + | σ i | 2    (19)
s.t. P min ≤ p i ≤ P max , ∀i ∈ K Similarly as for (M 1 ), we redefine the nonnegative parameters

a ij =| g H i g j | 2 > 0 and b i =| σ i | 2 > 0 for all i, j ∈ K.
Claim 1 (M 2 ) can be equivalently written as a signomial programming problem.

Proof First, notice that

i∈K log 2   1 + p i | g H i g i | 2 j∈K (j̸ =i) p j | g H i g j | 2 + | σ i | 2    = log 2    i∈K   1 + p i | g H i g i | 2 j∈K (j̸ =i) p j | g H i g j | 2 + | σ i | 2      
Next, we introduce a lower bound variable t i for each i ∈ K in order to obtain an optimal solution for (M 2 ) by solving the following equivalent problem max

{(p,t)∈R 2K ++ } i∈K t i s.t. t i ≤ 1 + a ii p i j∈K,(j̸ =i) a ij p j + b i , ∀i ∈ K (20)
P min ≤ p i ≤ Pmax, ∀i ∈ K
In particular, we note that the constraints (20) can be written as

j∈K,(j̸ =i) a ij p j t i + b i t i j∈K,(j̸ =i) a ij p j + b i + a ii p i ≤ 1, ∀i ∈ K (21)
or as j∈K,(j̸ =i)

a ij b -1 i p j t i + t i - j∈K a ij b -1 i p j ≤ 1 ∀i ∈ K (22)
Notice that the left-hand side of constraint [START_REF] Le | Enabling User Grouping and Fixed Power Allocation Scheme for Reconfigurable Intelligent Surfaces-Aided Wireless Systems[END_REF] has negative coefficients which show that (M 2 ) leads to a signomial programming problem. □

Because of Claim 1, we see that finding a global optimal solution for (M 2 ) is significantly harder than solving a GP problem. However, we can still approximate each posynomial function in the denominator of the constraints [START_REF] Zhang | Robust Energy-Efficient Transmission for Wireless-Powered D2D Communication Networks[END_REF] with a monomial function to obtain a lower bound for (M 2 ) [START_REF] Boyd | A Tutorial on Geometric Programming[END_REF]. This can be achieved by using the arithmetic-geometric mean inequality. Thus, we have the following proposition.

Proposition 1 For fixed and nonnegative parameter values δ ij and β i for all i, j ∈ K such that j∈K δ ij + β i = 1, for all i ∈ K, the following GP model gives a lower bound for (M 2 )

Q(δ, β) = max {(p,t)∈R 2K ++ } i∈K t i (23)
s.t. l∈K p -δ il l   j∈K,(j̸ =i) c ij p j t i + d i t i   ≤ 1, ∀i ∈ K (24)
P min ≤ p i ≤ Pmax, ∀i ∈ K
where

c il = a il j∈K a ij δ ij -δij (b i β -1 i ) -βi , ∀i, l ∈ K, (i ̸ = l)
and

d i = b i j∈K a ij δ ij -δij (b i β -1 i ) -βi , ∀i ∈ K
Proof Notice that by applying the arithmetic-geometric mean inequality we have that

j∈K a ij p j + b i ≥ j∈K a ij p j δ ij δij (b i β -1 i ) βi , ∀i ∈ K (25)
for values of δ ij ≥ 0, β i ≥ 0 for all i, j ∈ K with respect to j∈K δ ij + β i = 1 for all i ∈ K. We do also have

j∈K,(j̸ =i) a ij p j t i + b i t i j∈K a ij p j + b i ≤ j∈K,(j̸ =i) a ij p j t i + b i t i j∈K aij pj δij δij (b i β -1 i ) βi ≤ 1, ∀i ∈ K (26)
Next, by performing simple algebraic manipulations in the latter inequalities we obtain

l∈K p -δ il l   j∈K,(j̸ =i) c ij p j t i + d i t i   ≤ 1, ∀i ∈ K, (27)
Finally, a lower bound for (M 2 ) can be obtained by calculating the logarithm of the objective function in [START_REF] Lu | Distortion-Aware Cross-Layer Power Allocation for Video Transmission Over Multi-User NOMA Systems[END_REF]. Thus, concluding the proof. □

Notice that the inequalities [START_REF] Cheng | Distributionally robust stochastic knapsack problem[END_REF] have zero gaps when all the involved terms in the sum and the product are equal. Furthermore, we observe from Proposition 1 that we can adjust the parameter δ ij and β i for all i, j ∈ K to tighten the lower bound obtained with Q(δ, β). More precisely, we see that the optimal values for the parameters δ ij and β i for all i, j ∈ K can be obtained by solving the following optimization problem

Q(δ * , β * ) = max {(δ,β)∈R K 2 +K + } Q(δ, β) (28)
s.t. j∈K δ ij + β i = 1, ∀i ∈ K (29)
Unfortunately, finding the global optimal solution of problem ( 28)-( 29) is not trivial as it is a non-convex problem. However, in an attempt to find its global optimal solution, we can still solve the problem with an iterative algorithm referred to as the single condensation method in the literature which is essentially a sequential approximation method [START_REF] Chiang | Power Control By Geometric Programming[END_REF]. In general, with this method, we can obtain a locally optimal solution to the problem. However, a common strategy may consist of generating different initial solutions for model ( 28)-( 29) to increase the chances to obtain the global optimal solution for model (M 2 ). In the next section, we explain how we apply this method to obtain a lower bound for (M 2 ). Notice that analogously as for the model (M 1 ), we can assume that the input parameters c ij and d i for all i, j ∈ K in model Q(δ, β) are random variables normally and independently distributed, i.e., c ij (ξ) ∼ N (c ij , σ 2 cij ) and

d i (ξ) ∼ N ( di , σ 2
di ), respectively. This allows writing the following stochastic version of Q(δ, β) with individual chance constraints (IM 2 ) : max

{(p,t)∈R 2K ++ } i∈K t i s.t. l∈K p -δ il l   j∈K,(j̸ =i) cij p j t i + di t i   + ϕ -1 (1 -α)      l∈K p -2δ il l   j∈K,(j̸ =i) σ 2 cij p 2 j t 2 i + σ 2 di t 2 i        ≤ 1, ∀i ∈ K (30)
P min ≤ p i ≤ P max , ∀i ∈ K
Observe that model (IM 2 ) is a GP problem and then, it can be solved optimally. Finally, we can also arrive at the following stochastic optimization problem with joint chance constraints (JM 2 ) : max

{(p,t,y)∈R 3K ++ } i∈K t i s.t. l∈K p -δ il l   j∈K,(j̸ =i) cij p j t i + di t i   + ϕ -1 (y i )      l∈K p -2δ il l   j∈K,(j̸ =i) σ 2 cij p 2 j t 2 i + σ 2 di t 2 i        ≤ 1, ∀i ∈ K (31) i∈K y i ≥ (1 -α) 0 ≤ y i ≤ 1, ∀i ∈ K P min ≤ p i ≤ P max , ∀i ∈ K
Similarly as for (JM 1 ), in the next section, we compute bounds for (JM 2 ) using piece-wise linearization and sequential approximation methods.

Piece-wise linear and sequential approximation methods

In this section, we present two GP models using a tangent piece-wise linear approximation method for computing lower and upper bounds for (JM 1 ) and (JM 2 ), respectively. Then, we present a sequential approximation algorithm to obtain upper and lower bounds for (JM 1 ) and (JM 2 ), respectively. Finally, we present a single condensed sequential approximation method used to compute lower bounds for model Q(δ, β) (See proposition 1). We mention that in order to solve the models (IM 2 ) and (JM 2 ), we use the best parameter values (δ * , β * ) obtained with the single condensation method.

GP models obtained using tangent piece-wise linear functions

According to Theorem 1 in [START_REF] Liu | Stochastic geometric optimization with joint probabilistic constraints[END_REF], we can compute lower and upper bounds for (JM 1 ) and (JM 2 ) by solving the following GP models (JM 1 ) Lb : min

{(p,t,y)∈R 2K+1 ++ } t -1 s.t. j∈K,(j̸ =i) āij p j p -1 i t + bi p -1 i t + e fs 2 y gs 2 i    j∈K,(j̸ =i) σ 2 aij p 2 j p -2 i t 2 + σ 2 bi p -2 i t 2    ≤ 1, ∀i ∈ K, s ∈ S ( 32
) i∈K y i ≥ (1 -α) 0 ≤ y i ≤ 1, ∀i ∈ K P min ≤ p i ≤ P max , ∀i ∈ K and (JM 2 ) U b : max {(p,t,y)∈R 3K ++ } i∈K t i s.t. l∈K p -δ il l   j∈K,(j̸ =i) cij p j t i + di t i   + e fs 2 y gs 2 i      l∈K p -2δ il l   j∈K,(j̸ =i) σ 2 cij p 2 j t 2 i + σ 2 di t 2 i        ≤ 1, ∀i ∈ K, s ∈ S ( 33
) i∈K y i ≥ (1 -α) 0 ≤ y i ≤ 1, ∀i ∈ K P min ≤ p i ≤ P max , ∀i ∈ K
respectively where the set S = {1, . . . , S} corresponds to an index set associated with a set of linear functions. Notice that to obtain a convex reformulation of the constraints ( 16), one can use the standard variable transformations x i = log(y i ), r i = log(t i ) and h i = log(p i ) for each i ∈ K. Thus, the constraints ( 16) can be equivalently written as j∈K,(j̸ =i)

e hj -hi+r+log(āij ) + e -hi+r+log( bi) + e 2 log(ϕ -1 (e x i ))

  j∈K,(j̸ =i) e 2hj -2hi+2r+log(σ 2 a ij ) + e -2hi+2r+log(σ 2 b i )   ≤ 1, ∀i ∈ K (34)
Notice that we can replace the terms e fs 2 y gs 2 i for all s ∈ S, i ∈ K in (JM 1 ) Lb and (JM 2 ) U b by tangent linear functions. However, this can be performed once the log transformations have already been applied to convexify the problem [START_REF] Boyd | A Tutorial on Geometric Programming[END_REF]. This leads to the constraints [START_REF] Dupacová | Stochastic geometric programming: approaches and applications[END_REF]. Thus, we approximate the function 2 log(ϕ -1 (e xi )) ≥ F s (x i ) = f s x i + g s for all s ∈ S, i ∈ K. For this purpose, one can choose tangent lines of 2 log(ϕ -1 (e xi )) at different points in [log(1 -α), 0), say χ 1 , . . . , χ S . Then, for each s ∈ S, we have that

g s = 2e χs (ϕ -1 ) (1) (e χs ) ϕ -1 (e χs )
and

f s = -g s χ s + 2 log(ϕ -1 (e χs ))
This is possible to achieve since the function 2 log(ϕ -1 (e xi )) is convex at the required interval. Consequently, we have that the feasible sets of (JM 1 ) and (JM 2 ) are contained in the feasible sets of (JM 1 ) Lb and (JM 2 ) U b , respectively. Finally, by reversing the log transformations, one can easily verify that ϕ -1 (y i ) ≥ e for all s ∈ S, and i ∈ K.

Sequential approximation methods

Now we explain the sequential approximation method used to compute upper and lower bounds for (JM 1 ) and (JM 2 ), respectively [START_REF] Liu | Stochastic geometric optimization with joint probabilistic constraints[END_REF]. For the sake of brevity, we only explain the method when applied to (JM 1 ). The idea is simple and consists of fixing the variable y = y n at iteration n while satisfying the feasibility of the constraints ( 17)-( 18). This allows solving the GP problem (JM 1 )(y n ) : min

{(p,t)∈R K+1 ++ } t -1 s.t. j∈K,(j̸ =i) āij p j p -1 i t + bi p -1 i t + ϕ -1 (y n i )    j∈K,(j̸ =i) σ 2 aij p 2 j p -2 i t 2 + σ 2 bi p -2 i t 2    ≤ 1, ∀i ∈ K P min ≤ p i ≤ P max , ∀i ∈ K Let (p n , t n , θ n , v n
) denote an optimal solution of (JM 1 )(y n ), an optimal solution of the Lagrangian dual variable θ and the optimal objective function value, respectively [START_REF] Liu | Stochastic geometric optimization with joint probabilistic constraints[END_REF]. We then compute the search direction parameters

φ i = θ n i (ϕ -1 ) ′ (y n i ) j∈K,(j̸ =i) σ 2 aij (p n j ) 2 (p n i ) -2 (t n ) 2 + σ 2 bi (p n i ) -2 (t n ) 2
for each i ∈ K and update the y variable by solving the following GP problem

(JM 1 )(p n , t n ) : min {y∈R K + } i∈K φ i y i s.t. y i ≤ ϕ   1 -j∈K,(j̸ =i) āij p n j (p n i ) -1 t n -bi (p n i ) -1 t n j∈K,(j̸ =i) σ 2 aij (p n j ) 2 (p n i ) -2 (t n ) 2 + σ 2 bi (p n i ) -2 (t n ) 2   , ∀i ∈ K i∈K y i ≥ (1 -α) 0 ≤ y i ≤ 1, ∀i ∈ K
where ϕ(•) denotes the normal cumulative distribution function. The iterative sequential method is depicted in Algorithm 1 Finally, we present a single condensed sequential approximation method used to compute lower bounds for Q(δ, β) according to [START_REF] Chiang | Power Control By Geometric Programming[END_REF]. This procedure is depicted in Algorithm 2 as follows. The procedure is simple and starts by initializing the input parameter values δ and β while ensuring that the conditions (29) are satisfied. Then, we iterative until the difference between the objective function values obtained with model Q(δ, β) in two consecutive iterations is less than or equal to a small positive value ϵ. Notice that at each iteration, the input parameters δ and β are updated according to the monomials approximations [START_REF] Cheng | Distributionally robust stochastic knapsack problem[END_REF]. Finally, the best solution obtained for Q * = Q(δ * , β * ) is returned by the algorithm together with the optimal values of parameters δ * and β * . Recall that to solve models (IM 2 ) and (JM 2 ), we use these parameter values. Notice that Algorithm 2 can be sensitive to the initial values of δ and β.

Algorithm 1 Sequential convex approximation method

Require: An instance of problem (JM 1 ) Ensure: An upper bound for (JM 1 )

1: Initialization step: 2: Choose an initial point y 0 of y satisfying the constraints ( 17)-( 18).

3: Solve (JM 1 )(y 0 ) and use the optimal solution obtained for solving (JM 1 )(p 0 , t 0 ). 4: Set n = 1. Denote by y n the optimal solution of (JM 1 )(p 0 , t 0 ) 5: Iterative step: 6: while (∥y n-1 -y n ∥ > ε and number of iterations <= M axIter) do 7:

Solve (JM 1 )(y n ) and use the optimal solution obtained for solving (JM 1 )(p n , t n ).

8:

Let ȳ denote the optimal solution of (JM 1 )(p n , t n ) 9:

y n+1 = y n + τ (ȳ -y n ).
Here, τ ∈ (0, 1) is the step length. Solve the GP problem Q(δ, β).

14:

Let (p, t) ∈ R 2K and W Q denote the optimal solution and objective function value obtained.

15: if (| W Q -V Q |≤ ϵ
W Q = V Q 19: end if 20: end while 21: return best solution obtained (p * , t * , W * Q , δ * , β * )
Consequently, to increase the chances of obtaining a global optimal solution for Q = Q(δ, β), an intuitive approach is to run Algorithm 2 several times by using different initial values for the parameter δ and β. Although this approach does not guarantee finding the global optimal solution to the problem, we can still measure the difference between the stochastic and deterministic models derived from Q(δ * , β * ).

Numerical experiments

In this section, we conduct substantial numerical experiments to compare all the proposed models and the sequential approximation algorithms. For this purpose, we implement python codes using CVXPY version 1.1.12 interfaced with Mosek solver version 9.2.42 [START_REF] Diamond | CVXPY: A Python-embedded modeling language for convex optimization[END_REF][START_REF] Agrawal | A rewriting system for convex optimization problems[END_REF][START_REF]MOSEK solver version 9.2.42[END_REF]. The numerical experiments were carried out on an Intel(R) 64 bits core (TM) with 2.50 GHz and 12GB of RAM under Windows 10. We present numerical results for all the proposed models derived for (M 1 ) and (M 2 ) where the objectives are to maximize the worst user signal-to-interference noise ratio and the total capacity of the network, respectively. In particular, for (M 1 ), we consider instances dimensions ranging from K = 10 to K = 50 users. Whereas for (M 2 ), we generate instances ranging from K = 10 to K = 22 users. We notice that the objective functions of both models (M 1 ) and (M 2 ) deteriorate significantly when the number of users increases. This is clear since the more users the network has, the higher the interference generated. In particular, for (M 2 ) we did not generate instances with more than K = 22 users because doing this leads to infeasible solutions. We arbitrarily set the input parameters values to P min = 0.1, P max = 0.5 and generate the complex vectors g i ∈ C T ×1 , and g H i ∈ C 1×T for each i ∈ K according to an independent complex Gaussian distribution function with zero mean and variance equal to one. Then, we multiply each of these vectors by a factor of 2.5. Finally, the parameter σ i for each i ∈ K is also generated according to an independent complex Gaussian distribution function with zero mean and variance equal to one. The parameters a ij and b i are calculated as explained in Section 3. Without loss of generality, we assume in our stochastic GP models that āij = a ij and bi = b i for all i, j ∈ K coincide with the means of the random variables. In all our experiments, we vary the parameter values of α ∈ {0.1; 0.25}, and σ bi , σ aij ∈ {10 -3 ; 10 -2 ; 10 -1 ; 1; 2} for all i, j ∈ K. In order to solve the GP models (JM 1 ) Lb and (JM 2 ) U b , we consider the set S = {5, 10, 20} to reference the number of tangent linear segments. Finally, the parameter values in Algorithms 1 and 2 are set to ε = 10 -4 , M axIter = 50, τ = 0.5, y i = 0.9999999 for all i ∈ K as initial values, and ϵ = 10 -6 , M axIt = 1000, respectively. Notice that the initial values of y i for all i ∈ K give a value of 5.199337582290661 for the inverse of the standard normal cumulative distribution function. Thus, avoiding the value of 1 which gives infinite for this function. Moreover, these values satisfy all the required constraints of the GP models of Algorithm 1.

Worst user SINR

In Table 1 of Annex A of the supplementary material, we present average numerical results for the GP models (M 1 ), (IM 1 ), and (JM 1 ) over 20 samples. More precisely, in column 1 we present the number of users of each particular instance. Whilst in columns 2-5, 6-9 and 10-13 we present the optimal objective function value, the CPU time in seconds required by Mosek and CVXPY solvers and the standard deviation obtained with the models (M 1 ), (IM 1 ) and (JM 1 ), respectively. In particular, model (JM 1 ) is solved with Algorithm 1. Consequently, in column 14, we present the number of iterations required by Algorithm 1 to obtain the solution.

From Table 1, first, we observe that the optimal objective function values obtained with model (JM 1 ) are larger than those reported for (IM 1 ). In turn, the objective values obtained with (IM 1 ) are larger than the ones obtained with (M 1 ). Next, we see that these differences increase when using higher values of σ and α. Notice that by using higher values for the parameters σ and α, we can obtain more conservative solutions to the problem. Next, we observe that the CPU time required by Mosek solver is significantly smaller than the one reported for CVXPY. This can be explained by the parser property of CVXPY solver which is required to convexify the GP problem [START_REF] Boyd | A Tutorial on Geometric Programming[END_REF]. Next, we observe that the standard deviation values are large and nearly the same for all the instances. This fact gives an idea of how sensitive the solutions obtained with the GP models are depending on the input data. Finally, we observe that the number of iterations is on average less than 22 which evidences the effectiveness achieved by the Algorithm 1.

In Tables 2,3, 4, and 5 (See Annex A), we report numerical results obtained with model (JM 1 ) using randomly generated instances for different values of α and σ. We solve the same instances in these tables. The main purpose of these tables is to show how tight the GP models using tangent piece-wise linear functions are when compared to the sequential approximation Algorithm 1. The legend of these tables is the same and is as follows. In column 1, we present the number of users considered in the network. Next, in columns 2-6 we present the number of tangent linear segments used in model (JM 1 ) Lb and its number of constraints, its optimal objective function value, and the CPU times in seconds obtained with Mosek and CVXPY solvers, respectively. From columns 7 to 10, we present the optimal objective function value of (JM 1 ) obtained with Algorithm 1, the CPU times in seconds required by Mosek and CVXPY solvers, and the number of iterations required by Algorithm 1 to solve (JM 1 ). Finally, in column 11 we respectively report the gaps that we compute by v((JM1))-v((JM1) Lb ) v((JM1)) * 100 where v(•) denotes the optimal objective function value of the respective model. Notice that we report these gap values for a different number of tangent linear segments used in model (JM 1 ) Lb . From Tables 2, 3, 4, and 5, we mainly observe that the objective function values obtained with (JM 1 ) are slightly higher than those obtained with (JM 1 ) Lb . In particular, we see that the gap values reported in column 11 in each of these tables become smaller when using more tangent linear functions in (JM 1 ) Lb . Similarly, we observe that the gap values obtained with σ = 0.01 are smaller than those obtained when using σ = 0.1. Next, we also see that the optimal objective function values reported in Tables 4 and5 for both models (JM 1 ) Lb and (JM 1 ) are larger than in Tables 2 and3, respectively. We notice that this is a consequence of using a larger value of parameter α which is equivalent to imposing a stronger probabilistic constraint in each model. As a consequence, we obtain more conservative solutions with each model.

Another relevant observation is that the CPU times reported in columns 5 and 6, and in columns 8 and 9 for all the instances in each of these tables are significantly higher for CVXPY than for the Mosek solver. Recall that the CPU time of CVXPY includes the CPU time of the Mosek solver. This shows that CVXPY spends a considerable amount of time transforming the problem from GP format to a convex problem that is solvable by the Mosek solver. Finally, we observe that the number of iterations required by Algorithm 1 to solve all the instances remains nearly the same in Tables 2, 3, 4, and 5.

To give more insights regarding the conservatism level achieved with the solutions obtained with models (M 1 ), (IM 1 ), (JM 1 ) and (JM 1 ) Lb , we report in Figures 1234the number of violated constraints and amount of constraint violations for an instance with K = 30 users while using parameter values of α = 0.25 and σ = {0.1, 1}. More precisely, we save the optimal solutions obtained with the models (M 1 ), (IM 1 ), (JM 1 ), and (JM 1 ) Lb for a particular instance randomly generated. And then, we randomly generate another 100 instances (scenarios) and count the number of violated constraints in [START_REF] Boyd | A Tutorial on Geometric Programming[END_REF] and their violated amounts obtained by replacing the saved solutions.

From Figures 1 and2, we mainly observe that the optimal saved solution obtained with (IM 1 ) is more robust than the one obtained with (M 1 ) in the sense that it allows obtaining less violated constraints and a lower amount of violation. Similarly, we observe that the optimal saved solution obtained with (JM 1 ) using Algorithm 1 is more conservative than the one obtained with (IM 1 ). Finally, we see that the solution obtained with (JM 1 ) Lb is less conservative than the one obtained with (JM 1 ) using Algorithm 1. In conclusion, we see that the solution obtained with the joint chance-constrained model (JM 1 ) is significantly more robust than the ones obtained with the other models. In Figures 3 and4, we observe similar trends as in Figures 1 and2. Ultimately, we notice that the higher the value of parameter σ, the higher the level of robustness achieved with the obtained solutions.

Maximum capacity

Now, we report numerical results for the second power allocation problem which is aimed to maximize the total capacity of the wireless network. In particular, in Table 6 (See annex B in the supplementary material), we report numerical results for Q(δ, β). These results are obtained with the condensed sequential approximation Algorithm 2. In this table, we do also report numerical results for the model (IM 2 ). Recall that the latter model is solved using the best values obtained for the parameters δ and β when solving Q(δ, β) with Algorithm 2. In Table 6, column 1 shows the number of users of each instance. Columns 2 to 6 present the objective function value of Q(δ, β), the capacity achieved which is computed by taking the logarithm of the objective function value, the CPU time in seconds required by the Mosek and CVXPY solvers to obtain the solutions and the number of iterations required by Algorithm 2. Columns 7-10 report the optimal objective function value of model (IM 2 ), its capacity which is again computed by taking the logarithm of the objective value, and the CPU times in seconds required by the Mosek and CVXPY solvers, respectively. The instances are randomly generated and solved for different input values for the parameters α and σ. From Table 6, first, we observe that the capacity values decrease significantly when the number of users increases. This is evident since the higher the number of users, the higher the interference in the network. Next, we also see that the CPU times required by CVXPY increase considerably with the number of users. Similarly, the number of iterations required by Algorithm 2 to converge also increases with the number of users. Next, we observe that the optimal objective function values obtained with (IM 2 ) decrease when the parameter σ increases. On the opposite, these objective values increase with an increase of parameter α. Finally, we observe that the CPU time values required by CVXPY are significantly higher than those required by the Mosek solver.

In Tables 7 and8 in annex B, we present upper bounds obtained with model (JM 2 ) U b and lower bounds obtained with model (JM 2 ) while using the sequential approximation Algorithm 1 for given and fixed values of δ and β. The values of δ and β are obtained by solving Q(δ, β) with Algorithm 2 for the instances in Table 6 of annex B. In Tables 6,7, and 8, we solve the same instances. In Tables 7 and8, the legends are the same. In column 1, we present the number of users of each instance. Next, in columns 2 and 3, we present the number of tangent linear segments and constraints of model (JM 2 ) U b . Subsequently, in columns 4, 5, and 6, we report the optimal objective function values of (JM 2 ) U b and the CPU times in seconds required by the Mosek and CVXPY solvers, respectively. Similarly, in columns 7-10 we present the optimal objective function values of model (JM 2 ), the CPU times in seconds required by the Mosek and CVXPY solvers, and the number of iterations obtained when solving (JM 2 ) with Algorithm 1, respectively. Finally, in column 11 we present gap values which are computed by v((JM2) U b )-v((JM2)) v((JM2)) * 100 where v(•) denotes the optimal objective function value of the respective model. From Tables 7 and8, we mainly observe that the optimal objective function values obtained with (JM 2 ) U b and (JM 2 ) get closer when using a higher number of piece-wise linear segments (S = 20). We further notice that these objective values deteriorate when using lower and higher values of α and σ, respectively. Next, we see that the objective function values decrease when more users are present in the network. Notice that we do not consider more than 22 users. Otherwise, we may obtain infeasible solutions with negative capacity values. We also observe that the CPU time values obtained with Algorithm 1 are larger than those obtained with (JM 2 ) U b for all the tested instances using less than S = 20 linear segments. For the remaining ones, we observe the opposite situation. Finally, we see that the number of iterations required by Algorithm 1 to solve all the instances is very stable. To give more insights concerning the solutions obtained with the condensed sequential approximation Algorithm 2, we plot in Figures 5 and6, the optimal solutions and the number of iterations obtained with Q(δ, β) for K = 10 and K = 14 users while generating different initial input values for the parameters δ and β. The idea is to provide some empirical evidence indicating that the solutions obtained with Algorithm 2 are near the global optimal ones. In particular, in Figures 5 and6 we randomly generate 100 and 20 initial sample values for the parameters δ and β, respectively. All these values are generated while satisfying the conditions imposed in problem ( 28)-( 29). In Figure 6, we only generate 20 samples for δ and β since the CPU times required to solve Q(δ, β) with Algorithm 2 become larger for K = 14 users. From Figures 5 and6, we see that it is not always possible to obtain the same optimal solution when using different initial samples for the parameters δ and β. This confirms that Algorithm 2 does not guarantee finding the global optimal solution to the problem. It only allows finding locally optimal solutions to the problem. However, we also notice that a recurrent maximum value is attained for many samples which might suggest that it is indeed the global optimal solution. Regarding the number of iterations, we do not observe a clear pattern as these values oscillate within a large range of possible values. Finally, to report some empirical evidence that we observed when solving our instances regarding the slow convergence of Algorithm 2, in Figure 7 we plot 10 curves with the optimal solutions obtained in each iteration during 10 runs of Algorithm 2 while using different initial values for the parameters δ and β. From Figure 7, we mainly observe that only after 100 iterations do we start finding better objective function values for the problem. Also, we notice that only for 3 out of 10 runs do we obtain the larger objective function values which are above the value of 10. Most of the remaining runs converge in less than 250 iterations and with a lower objective function value. This explains the CPU time values reported in Table 6 for Algorithm 2. To conclude, we cannot ensure that global optimal solutions are obtained with Algorithm 2. However, for a few small-size instances, we provide some empirical evidence indicating that the solutions obtained are near-optimal ones.

Conclusions

In this paper, we consider the problem of power allocation for 5G wireless networks using massive multiple input multiple output technologies. We propose two non-linear optimization models where the aims are to maximize the worst user signal-to-interference noise ratio and the total capacity of the network subject to power constraints. In particular, our first model is transformed into an equivalent geometric programming problem. Whereas the second one is formulated as an equivalent signomial programming problem. Since the latter problem is non-convex, we propose an approximated geometric model to compute lower bounds. This is achieved by using the arithmetic-geometric mean inequality. Subsequently, we derive stochastic formulations for each geometric programming model to deal with the uncertainty of wireless channels. More precisely, we include individual and joint probabilistic constraints on each GP model. Finally, we compute lower and upper bounds by using sequential convex and piece-wise linear methods. Substantial numerical experiments are conducted to compare all the proposed models and algorithms. Our numerical experiments indicate the following relevant conclusions.

1. First, we conclude that the use of joint chance constraints allows obtaining more conservative solutions for the power allocation problems than using individual chance constraints or using directly the deterministic models. 2. Next, we observe that the sequential approximation Algorithm 1 together with the piece-wise GP models allows obtaining tight bounds while using joint chance constraints. The latter is a crucial fact as it shows that we can approximate the globally optimal solutions to the problems in an efficient manner. Moreover, we observe that this sequential method requires only a few iterations to converge for most of the tested instances. 3. For the signomial problem, we mainly observe that the condensed sequential approximation method requires a significantly larger number of iterations to converge than the other sequential method. Although, as shown in the numerical results section, we can still use it with different initial values for the required parameters to improve the solutions obtained with the deterministic model. 4. We observe that independently of the parameter values obtained when solving the deterministic signomial problem, we can still obtain tight bounds for the stochastic models derived from it. 5. Finally, as part of future research, we remark that all the modeling approaches and solution methods studied in this paper can be further investigated and applied to many other relevant problems related to resource allocation in wireless networks. Ultimately, we believe that novel algorithmic approaches including meta-heuristics should also be part of future research, especially for solving the signomial and stochastic programming models.

10 : 1 : 2 : 3 : 5 :

 101235 Set n = n + 1. 11: end while 12: return the best solution found (p n , t n , y n )Algorithm 2 Single condensed sequential approximation methodRequire: An instance of problem Q(δ, β) Ensure: A lower bound for Q(δ, β) Initialization step: Randomly generate positive values for the parameters δ and β satisfying the conditions[START_REF] Prékopa | Stochastic Programming[END_REF]. Solve the GP problem Q(δ, β). 4: Let (p, t) ∈ R 2K and V Q denote the optimal solution and objective function value obtained. Iterative step: 6: while (True) do

Fig. 1 1 Fig. 2 1 Fig. 3

 11213 Fig. 1 Number of violated constraints for each randomly generated instance using K = 30, α = 0.25 and σ = 0.1

Fig. 4

 4 Fig.4The total amount of constraint violation for each randomly generated instance using K = 30, α = 0.25 and σ = 1

Fig. 5

 5 Fig. 5 Optimal solutions and the number of iterations obtained with Algorithm 2 for 100 initial samples of parameters δ and β using K = 10.

Fig. 6

 6 Fig.6Optimal solutions and the number of iterations obtained with Algorithm 2 for 20 initial samples of parameters δ and β using K = 14.

Fig. 7

 7 Fig. 7 Slow convergence of the condensed sequential approximation method reporting optimal solutions and the number of iterations obtained for an instance of K = 16 users while using different initial values of parameters δ and β.

  or number of Iterations ≥ M axIt) then

	16:	Break
	17:	else
	18:

Acknowledgments

The authors acknowledge the financial support from Project: FONDECYT No. 11180107.

Declarations

Data availability

The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

Conflict of interest:

The authors have no conficts of interest to declare that are relevant to the content of this article.


Clément Iphar
email: clement.iphar@univ-brest.fr

Coraline Soul

Iwan Le Berre

Manuel Sahuquet

Éric Foulquier

Aldo Napoli

Unveiling of Maritime Networks from the Differentiated Analysis of Vessel Traffic in the Caribbean

Introduction

The Caribbean region is a critical area of the global maritime traffic, acting as an interface between the Atlantic and the Pacific through the Panama Canal, for various types of vessels ranging from cargo ships to oil tankers and cruise liners. Local ports serve not only as infrastructure nodes for goods and passenger transit but also as complex hubs that see intricate networks of maritime routes, both global and local. These networks are useful for understanding global and local logistics, and environmental impact.

The presence of different types of vessels, each having distinct operational behaviours and routing protocols, suggests that multiple complex networks coexist within this maritime space, though data on port activity, and all the more so on maritime routes activity, lack, or are incomplete and hard to reach.

In this paper, we use a dataset of port calls in the Caribbean that we previously developed in order to provide a comprehensive analysis of the networks formed by different types of vessels in the Caribbean. Investigating how these distinct networks interact and coexist in the same geographical space can open up interesting perspectives for policy development in fields of maritime safety, environmental protection or regional logistics.

A visualisation dashboard has been developped to facilitate the geovisualisation of these ports calls as a proxy for port activity [START_REF] Bai | Quantifying the impact of pandemic lockdown policies on global port calls[END_REF] and the resulting traffic networks, as well as personalised indicators on maritime traffic in the area.

Data

We have extracted a set of port calls with a method presented in [START_REF] Iphar | Escales portuaires dans les Petites Antilles : une approche désagrégée[END_REF], using three main data sources, AIS vessel position data, a set of ports and a vessel registry, described in the next three paragraphs. All the data will be catalogued and made available shortly, along with an information document.

The AIS dataset has been purchased from the firm exactEarth. Its timespan covers the whole of 2019, and its spatial extent spans from 3.86 degrees North to 34.05 degrees North in terms of latitude, and from 98.02 degrees West to 51.17 degrees West in terms of longitude. It covers the entire Caribbean, including the Gulf of Mexico. The whole dataset encompasses a grand total of 641,709,724 data contacts, which represents a mean value of 1.76 million AIS messages per day. Figure 1 represents all data contacts for the sole month of January 2019 (46 million messages), that clearly shows the bounding box of our dataset. In order to have a cartographic reference frame for calculating port calls, a dataset describing the infrastructure (quays and wharves) of all ports of interest in the Caribbean has been produced by the authors using GIS photo-interpretation. Digitalisation was carried out at a scale of 1:5000, using Google Earth and Bing Map satellite imagery, and the OpenStreetMap cartographic repository. A total of 528 ports, spanning across 41 territories (sovereign nations or dependencies) have been identified and labelled.

As we want to have a differentiated approach of the networks on the basis of the vessel type, it is important to label vessels. For this purpose, the IHS Markit database classification, otherwise purchased, has been used. Within all categories of vessels, the types we investigate are bulk carriers, container ships, cruise vessels, general cargo vessels, inter-island ships and tankers.

Once the port calls are computed using state-of-the-art parameters [START_REF] Li | Semantic Recognition of Ship Motion Patterns Entering and Leaving Port Based on Topic Model[END_REF], the semantic trajectories of the vessels going from port to port are reconstructed, and networks of ports are built by putting together all trajectories. Therefore, a total of 298,804 directional journeys from one port to another have been computed, out of 9,969 individual pair of ports.

A dashboard for network analysis

In order to display information on port calls and local networks, a geovisualisation dashboard has been developed. This dashboard enables the visualisation of the number of port calls at each port, characteristics such as the mean age and tonnage of vessels calling, and the exploration of data in an interactive way. More particularly, it is possible, as shown in Figure 2 to focus on a particular type of vessel, only showing edges with a cardinality within a Fig. 2. Local networks for all vessels (top), tankers (middle) and general cargo vessels (bottom). Only edges with more than 50 journeys for all vessels, 5 journeys for tankers and 4 journeys for general cargo vessels are shown in this figure .   certain range, or visualise the ports from which originate the vessels calling at a given port, and inversely the destination ports of the vessels when leaving a given port, with the corresponding cardinalities of each link.

Perspectives

This dashboard enables to dig into data and perspectives are manifold. One promising direction involves understanding the spatiotemporal organisation of maritime traffic, taking into account the variability introduced by different types of traffic, particularly in an archipelago which undergoes a two-tier maritimisation system, propelled by domestic needs on the one hand and the needs of a globalised trade on the other hand [START_REF] Foulquier | Maritimization in The Lesser Antilles: Between Global and Local Stakes[END_REF]. This approach could offer critical insights for the effective deployment of maritime surveillance resources and logistical planning.

Another concern is the evaluation of the environmental pressures generated by this maritime traffic. Potential impacts on marine ecosystems and vulnerable coastal populations could be modeled by linking vessel types to pollutant emissions, which would be particularly useful for environmental regulation and policy implementation.

The identification of maritime traffic hotspots in terms of biodiversity also constitutes a research priority, particularly in the Caribbean where the proliferation of protected marine areas testify the will of neighbouring states to protect the environment. Overlaying shipping routes with ecologically significant areas could help identify highrisk regions where maritime traffic could have adverse effects on local biodiversity.

Fig. 1 .

 1 Fig. 1. AIS contacts of exactEarth dataset, for January 2019

Acknowledgments

This work has benefited from several state fundings managed by the Agence Nationale de la Recherche under the "Investissements d'avenir" program: the Labex DRIIHM ANR-11-LABX-0010, and the ISblue project "Interdisciplinary graduate school for the blue planet" ANR-17-EURE-0015. It has also received funding from the European Union's Horizon 2020 research and programme under the Marie Sklodowska-Curie grant agreement No 899546.