Prof Thomas Bauschert

Prof Akli Fundo

Dr Jean-Benoist Leger

Prof Andreas Kassler

Dr Liyang Xiao

Superfluid F 5g

less Networks. 5G SF Network 5G Superfluid wireless Networks

Keywords: BLP Binary Linear Programming GA Genetic Algorithm. GMM Gaussian Mixture Model. HD High Definition ILP Integer Linear Programming. InfraP Infrastructure Provider. IP Internet Protocol. ISPs Internet Service Providers. LP Linear Programming. LTE Long Term Evolution. MC RFB Reusable Functional Block. RINS Relaxation Induced Neighborhood Search. RO Robust Optimization RVSM Reduced Vertical Separation Minima. SN Substrate Network VNF Virtual Network Function. VNFC Virtual Network Function Component. VNFP Virtual Network Function Placement. VNO Virtual Network Operator Optimization réseaux, Optimisation robuste, Chance-constrained Programming, Approches exactes, Metaheuristiques

Fabio D'Andreagiovanni for his valuable advice, constructive assistance, and precious help. I would like to thank him for his optimism, positivity, and enthusiasm. He always encouraged me whenever I felt nervous somehow. I am sincerely grateful for his company of the interview of P.h.D candidate, laboratory days of Heudiasyc, preparation and setdown in the

List of Tables

3.8.10Continuous table of Table

 3.8.9 . . . . . . . . . . . . . . . . . . . . 4.5.1 Input parameters in line with these presented in [START_REF] Chiaraviglio | Optimal design of 5g superfluid networks: Problem formulation and solutions[END_REF]

Introduction

In our everyday life, infrastructures using networks have deeply pervaded our activities by putting at our disposal a wide range of services of very different nature that have become essential for us. If we want to give only a few major examples, we can think of telecommunications networks, both wireless and wired, by which we communicate via smartphones and internet; transportation networks, which allow us to move by air, ground and sea; and power networks, which transport and distribute the electrical energy powering a high number of personal and family devices that we use everyday. All these infrastructures are connected through a network lying at their core, which can be generically defined as a set of entities in which couples of them are connected according to some relation. When dealing with networks, it comes natural to model them by graphs made up of a set of nodes and a set of edges expressing the connection between pairs of nodes. Network optimization problems are thus naturally associated with problems related to graphs and typically model events in the network under the form of flows moving across the network. For an exhaustive introduction and discussion about network optimization problems, we refer to the famous book [START_REF] Ahuja | Network Flows: Theory, Algorithms, and Applications[END_REF].

Given the interest of this Ph.D. Thesis in studying and proposing new results about optimization methods for network optimization problems arising in realworld applications, an important aspect that we had to take into account in our research has been the presence of data uncertainty, namely the fact that (part of) the input data of the problem at hand are typically not exactly known in value when the problem is solved.

To tackle a real-world problem, we commonly establish a deterministic mathematical model by assuming that the data inputs are precisely known and try to solve such model in an exact or approximate way in order to obtain a qualified so-called "optimal solution" in each context. However the impact of data uncer-

Chenghao WANG

Contribution to robust network optimization tainty in these constructed models on the a-posteriori optimality and feasibility of the obtained solutions are ignored. Thereby, the solution generated from reference data that differ from the actual values encountered in reality may lead to identify solutions that result far away from actual optimal solutions and may even violate feasibility constraints. [START_REF] Ben | Robust solutions of linear programming problems contaminated with uncertain data[END_REF] have shown that even small (0.01%) perturbations in the value of input data can result in highly infeasible solutions for some benchmark optimization problems. Consequently, the need for methodologies capable of generating a robust solution, namely a solution immune to data uncertainty, naturally arises. The basic idea of robust optimization is to seek a solution which remains feasible and optimal even when deviations in the input data of an optimization problem occur.

To deal with relevant sources of data uncertainty arising in real-world network optimization problems related to transportation and communications networks, we have obtained a number of original modeling, algorithmic and computational results that are presented in the reminder of this Ph.D. Thesis, for which we provide an overview of the structure and of the main original contributions in the next sections.

Thesis organization

The thesis is organized in 6 chapters, including the introduction and the conclusion, as follows:

• Chapter 2: Optimization Methods under Data Uncertainty. This chapter is devoted to provide a concise introduction to optimization under uncertainty, highlighting the issues of dealing with uncertain data in optimization and offering an overview of major methodologies proposed over time.

• Chapter 3 : Robust Flight Level Assignment problem. This chapter is devoted to the flight level assignment problem. The problem arises in the Air Traffic Management (ATM) context where several flights compete to share the airspace resources, that is the flight levels. An appropriate assignment will lead to less conflicts and reduce the delays. We study in detail the Flight Level Assignment (FLA) problem and its robust variant. In practice, we model and solve a Chance-Constrained Programming (CCP) linear problem where the coefficients are uncertain but follow a known Gaussian Mixture distribution.

Our solution methodology stands in solving the CCP problem through an iterative heuristic approach. In each iteration a solution is provided and its feasibility probability evaluated. Then, an important element of the approach is the feasibility probability estimation where both analytical and sampling methods are proposed and experimented.

• Chapter 4: Resource Allocation in 5G Superfluid Wireless Networks. This chapter is devoted to address an optimization problem related to the design of 5 th generation of wireless networks (5G). Specifically, we have considered the 5G network architecture proposed and studied in the European Horizon 2020 project "Superfluidity" , detailed in [START_REF] Bianchi | Superfluidity: a flexible functional architecture for 5g networks[END_REF]. The 5G Superfluid architecture is based on atomic virtual entities called Reusable Functional Block (RFB)s, which are able to support the high level of flexibility, agility, portability and high performance required by 5G. We investigate the problem of minimizing the total installation costs of a Superfluid network composed of virtual and realized over a physical network, while guaranteeing constraint on user coverage, downlink traffic performance and technical constraints establishing relations between RFBs of different nature. We propose a new mathematical formulation which can enhance the possibility of optimally solving realistic networks instances and a Benders-like decomposition approach for accelerating the solution process.

• Chapter 5: Green and Robust 5G Virtual Network Function Placement Problem. We investigate the problem of optimally placing virtual network functions in 5G-based virtualized infrastructures according to a green paradigm that pursues energy-efficiency. This optimization problem can be modeled as an articulated Mixed Integer Linear Programming problem with a multicommodity flow model at its core. To model the data uncertainty that naturally affects the volume and features of the traffic associated with the requests of establishing virtual networks generated by users, we rely on adopting a robust optimization approach according to the Γ-robustness paradigm.

Since the resulting robust counterpart may easily become challenging to solve even for instances of moderate size for state-of-the-art solvers, we propose a new matheuristic for its solution. The matheuristic is based on combining a genetic algorithm with an exact large neighborhood search. Computational tests on realistic instances returned good results, showing that our algorithm can find better solutions in sensibly less time than a state-of-the-art solver.

• Chapter 6: Conclusions. This chapter concludes the Thesis, proposing a number of final considerations and possible directions for future work.

Contributions of this Ph.D. Thesis

The contribution of this thesis is threefold. We first investigate the Flight Level Assignment problem and its Chance-Constrained variant, proposing a compact deterministic mathematical formulation for the problem and a Chance-Constrained Programming counterpart. For its solution, we propose an iterative approach and another central original contribution is represented by the specific constraint generation approach solving the CCP problem associated with a given flight level.

Using the results of [START_REF] Klopfenstein | Tractable algorithms for chance-constrained combinatorial problems[END_REF], we reformulate the CCP as an ILP, in which the probability constraints are replaced taking into account the worst case, according to the procedure proposed in [START_REF] Allen | Convex programming with set-inclusive constraints and applications to inexact linear programming[END_REF]. These constraints are dynamically included when needed, through a constraint generation approach, until the desired feasibility probability is reached. All this gives a practical approach which may find application in a class of CCP problems. Another novelty of our approach is represented by the method used to check the feasibility probability of constraints. We study the case with uncertain parameters following a Gaussian Mixture Model (GMM) and propose an approximated method to estimate the feasibility probability of the solution. This may be of interest since the GMM is a powerful tool to capture characteristics distribution of a large number of real situations.

Concerning the design of 5G Networks based on the Superfluid architecture, the major contributions are constituted by proposing an alternative formulation for the optimization model by [START_REF] Bianchi | Superfluidity: a flexible functional architecture for 5g networks[END_REF] that we have taken as reference. Specifically, we were able to propose alternative feasibility constraints and characterize valid inequalities that express in a simpler way the technological constraints on the installation of the basic virtual entities, the so-called Reusable Functional Blocks, in distinct 5G network nodes. This simpler constraints lead to a new mathematical model that can be better handled by solvers. Moreover, to accelerate and improve the capacity of solving realistic instances, we have also proposed a new solution approach based on Benders decomposition that breaks the complete model into a master and a slave problems decoupling the complicated relations linking the installation of reusable functional blocks of different nature. Results of computational tests show the advantages of this new modeling and solution approach that we proposed.

Finally, the last contribution is related to the design of virtual networks according to a green network paradigm that pursues energy minimization. In this context, we have taken as reference state-of-the-art works and proposed a new effective and efficient matheuristic for solving the robust counterpart of the problem, exploiting the integration of a genetic algorithm with exact neighborhood searches, which formulate the exploration of (very) large neighborhood as mathematical programming problems solved to optimality by means of state-of-the-art solvers. The rationale at the basis of such matheuristic is that, while a state-of-the-art solver may not be able to solve the complete problem, it can instead efficiently solve to optimality suitable subproblems. Computational tests over realistic virtual network instances confirm the advantages of adopting such new matheuristic integration. The original results outlined above have been presented in the publications: and we are currently working to complete and finalize the journal versions of the conference papers.

Chenghao WANG Contribution to robust network optimization 2. Optimization Methods under Data Uncertainty

An introduction to data uncertainty

A fundamental assumption that is made in classical optimization is that all the data input of a problem are known exactly and precisely when the problem is solved. However, when dealing with optimization problems arising in the real world, one can often directly experience that such assumption does not hold and the input data are subject to uncertainty at some level, meaning that the value of (a subset of) coefficients appearing in the problem is not exactly known. The hypothesis of data certainty is typically adopted since this commonly leads to more tractable and less complex problems. However, obtained solutions under this artificial data certainty assumption may prove to be not useful in practice. As a consequence, adopting suitable modeling and solution techniques is imperative.

For an exhaustive introduction to the challenges and issues associated with optimization under data uncertainty, we refer the reader to the book [START_REF] Ben-Tal | Robust optimization[END_REF] and to the survey [Bertsimas et al., 2011a], which, though not very recent, still constitute major references for fundamentals of optimization under uncertainty. For more up-to-date surveys, we also refer to the work, [START_REF] Chenghao | Contribution to robust network optimization[END_REF], [START_REF] Yanıkoğlu | A survey of adjustable robust optimization[END_REF] and [START_REF] Leyffer | A survey of nonlinear robust optimization[END_REF]. As discussed in [START_REF] Ben-Tal | Robust optimization[END_REF], the presence of data uncertainty may be attributed to many causes, among which the most remarkable are:

• Errors due to prediction: in this case, the data are not known since they are related to events that will take place in the future and can be only guessed or forecast on the basis of historical data (when available). This is the case, for

Chenghao WANG Contribution to robust network optimization example, of traffic conditions in telecommunications networks which depend on the future behavior of customers and users and requires to be forecast (we note that in this case, practitioners tend to provide conservative predictions). Moreover, it can also be cited the case of resolution of flight conflicts in Air Traffic Management, which depends directly or indirectly on the prediction of weather conditions and flight trajectories.

• Errors due to measurements: it may happen that the data required in the optimization problem cannot be precisely measured and are naturally subject to measurement errors that must be taken into account when solving the problem. Referring to the case of telecommunications networks, such as the 5G virtualized networks that we consider, this could be represented by realtime measurements of the delay within the network.

• Errors due finite precision numerical representation: another source of uncertainty may be simply represented by the fact that optimization problems are commonly modeled and solved with by means of ad-hoc computer software, which rely on a finite precision representation of numbers and on finite precision arithmetic. In most applications, such error may be neglected. However, in some other, like wireless network design including signal-to-interference constraints, this error cannot be neglected, since it would lead to non accurate solutions, and must be taken into account (see [START_REF] Fabio | GUB covers and power-indexed formulations for wireless network design[END_REF]).

A very important observation that should be made at this point is that, over the years, many methodologies have been proposed in order to deal with data uncertainty in optimization problems and there is no just one "major and more effective" method that could adopted in all contexts. In what follows, we thus try to provide a concise overview of the main and most used methodologies, focusing more on those that are adopted as basis for our original developments.

Stochastic Programming

The first methodology that has considered the issue of data uncertainty in an optimization problem can be considered Stochastic Programming (SP), which can

Chenghao WANG

Contribution to robust network optimization be traced back to the seminal work of Dantzig presented in [START_REF] Dantzig | Linear Programming under Uncertainty[END_REF]. The main assumption at the basis of SP is that the stochastic distribution of data subject to uncertainty is known. The general form of a SP problem can be written as:

min E [c(x, R)] s.t. f i (x, R) ≥ b ∀i ∈ I x ≥ 0
In this problem, which aims at minimizing the expected value of the cost function, x ∈ R n is the vector of decision variables, c : R n → R is the cost objective function, f i : R n → R is the constraint function of the i th constraint. Also, R ∈ R |I| is a vector of random variables affecting the input coefficient and for which we know the corresponding distributions. For an exhaustive introduction to modeling and solution principles of Stochastic Programming, we refer the reader to [START_REF] Shapiro | Lectures on Stochastic Programming: Modeling and Theory, Second Edition[END_REF].

For many years, SP has represented the main methodology for dealing with data uncertainty in optimization and has been heavily studied and improved. However, a major limitation of SP that has been identified over the years is that it requires to know the probability distribution followed by the uncertain data and, as it is know, in many real-world application such distributions are not known. As a consequence, as also discussed in [START_REF] Ben-Tal | Robust optimization[END_REF], the application of SP is not so straightforward and accurate in a consistent number of relevant applications. Furthermore, it is also known that solving SP problems may be computationally challenging and expensive in practice, since typical solution approaches need to consider a high number of scenarios that are representative of the probability distributions of the uncertain data. This leads to very large problems that commonly require to be solved by suitable decomposition methods (see [START_REF] Shapiro | Lectures on Stochastic Programming: Modeling and Theory, Second Edition[END_REF]).

Robust Optimization

With the aim of overcoming some of the drawbacks of Stochastic Programming, the paradigm of Robust Optimization has been introduced at the beginning of the

Chenghao WANG

Contribution to robust network optimization new millennium. Specifically, a cornerstone of RO is to not deal with probability distributions, but to assume that the possible realizations of the uncertain data are completely specified through a so-called uncertainty set U . Then RO has the objective of finding the optimal solution of a robust counterpart of the original problem, identifying the solution that is optimal under all possible realizations of the data specified by the uncertainty set U . More formally, let us consider a general optimization problem:

min c(x) s.t. f i (x) ≥ b ∀i ∈ I x ≥ 0,
in which we aim at minimizing the cost function c : R n → R involving the vector of decision variables x ∈ R n , and in which f i : R n → R is the constraint function of the i th constraint. For a given uncertainty set U , the general form of the robust counterpart of the previous problem can be written as:

min max u∈U c(x) s.t. f i (x, u) ≥ b ∀i ∈ I, u ∈ U x ≥ 0
Feasible solutions of the previous problem must be feasible for all the realizations of the uncertain data specified by the uncertainty set U (indeed, we can notice that each feasibility constraint must be satisfied for each realization u ∈ U ). Such feasible solutions are called robust feasible solutions. Moreover, we can notice that the objective function takes the form of a min-max problem and identifies as robust optimal that solution granting the best cost value under the worst realization of the uncertain data of U . We remark that, without loss of generality, we can focus on the case of problems in which the uncertainty is just present in the input data appearing in the constraints. Indeed, if uncertain data are present in the objective or in the right-hand-sides of constraints, such uncertainty can be easily reformulated as uncertainty affecting only the constraints, as detailed in [Ben-Tal et al., 2009, Bertsimas and[START_REF] Bertsimas | The price of robustness[END_REF].

Chenghao WANG Contribution to robust network optimization

A critical task of RO is to "rationally " define the uncertainty set U . Indeed, since a solution is feasible for the robust counterpart only if it maintains its feasibility for all the realizations of the data specified by U , taking into account extreme and unlikely realizations of the data and including them in U may lead to very conservative solutions that are not useful in practice. As a consequence, a very substantial part of the literature of RO has been devoted to identifying those rational uncertainty sets that could provide satisfying level of robustness against data uncertainty without leading to over-conservative robust solutions. For an overview of the various kind of alternative uncertainty sets that have been proposed over time, we refer the reader to the book [START_REF] Ben-Tal | Robust optimization[END_REF] and the more recent survey [START_REF] Gabrel | Recent advances in robust optimization: An overview[END_REF]. In particular, these works proposing uncertainty sets have tried to overcome the limits of what can be considered the first (very conservative) example of robust optimization, proposed in 1973 [START_REF] Allen | Convex programming with set-inclusive constraints and applications to inexact linear programming[END_REF]. Here, we focus on so-called cardinality-constrained uncertainty sets that have constituted the most used sets and to which belongs the famous Γ-Robustness model by [START_REF] Bertsimas | The price of robustness[END_REF]] that we have taken as reference in our studies. For defining a cardinality-constrained set in a more formal way, let us introduce two vectors ū, û ∈ R p and let us impose that all the realizations u of the uncertain data must belong to the set [ū -û, ū + û]. The central value ū is commonly called nominal value, whereas û is commonly called highest or worst deviation. Furthermore, let us introduce an integer value β, called robustness budget, which indicate the highest number of input coefficients that may deviate from their nominal value simultaneously, then the cardinality-constrained uncertainty set may be written as:

U =            u ∈ R p : ūk -δ k ûk ≤ u k ≤ ūk + δ k ûk , p k=1 δ k ≤ β, δ k ∈ {0, 1}           
According to this formalization, at most β coefficients may deviate up to their worst deviation. We note also that such set could be easily defined also for fractional value of the coefficients δ (see e.g., [START_REF] Bertsimas | The price of robustness[END_REF]). After having provided such general definition of cardinality-constrained uncertainty set, we proceed to concentrate our attention on the famous Γ-robustness model.

Chenghao WANG Contribution to robust network optimization

As we discussed in the previous section, Γ-Robustness, proposed originally in the papers [START_REF] Bertsimasand | Robust discrete optimization and network flows[END_REF] and [START_REF] Bertsimas | The price of robustness[END_REF] belongs to the family of cardinality-constrained uncertainty sets. Formally, it can defined as follows. First of all, the assumptions at the basis of this model of uncertainty are:

1. Each entry u k of a vector u of the uncertainty set U represent a random variable following un unknown bounded and symmetrical probability distribution. These variables are assumed to be independent.

2. For each entry u k , the distribution is symmetrical with respect to a nominal value ūk and, given a maximum deviation ûk , is defined over the interval

[ū -û, ū + û].
3. Each realization of the data uncertainty u may have at most Γ ≥ 0 elements that deviate from the nominal value ūk .

Under these assumptions and exploiting the theoretical results presented in [START_REF] Bertsimasand | Robust discrete optimization and network flows[END_REF] and [START_REF] Bertsimas | The price of robustness[END_REF] , let us derive the robust of the following generic Linear Programming:

max c x s.t. a i x ≤ b i ∀i ∈ I x ≥ 0
Its generic robust counterpart can be written as:

max c x s.t. a i x ≤ b i ∀i ∈ I, ∀a i ∈ U i x ≥ 0,
in which we consider all the possible realizations of the row vector of constraint i that respect the assumptions on the uncertainty set. As next step, a critical observation that can be made is that the previous program may have exponentially

a i ∈U i a i x ≤ b i ∀i ∈ I x ≥ 0,
which, however, is non-linear due to the presence of the max term in the constraints. A major result by Bertsimas and Sim has been to elegantly prove that such non-linearity may be tackled by noticing that, for a constraint i and for a fixed vector x, computing the worst deviation expressed by the max term can be formulated as the following Binary Linear Programming problem:

max a i ∈U i a i x =                    max n j=1 āij xj + max n j=1 âij xj y j s.t. n j=1 y j ≤ Γ y j ∈ {0, 1} ∀j = 1, . . . , n                   
As well-noted by Bertsimas and Sim, an integral optimal solution of the previous problem can be obtained by solving its linear relaxation, which is:

max n j=1 āij xj + max n j=1 âij xj y j s.t. [z i ≥ 0] n j=1 y j ≤ Γ [v ij ≥ 0] 0 ≤ y j ≤ 1 ∀j = 1, . . . , n
Since the previous problem is linear, bounded and admits an optimal solution, we can define its dual problem which is also bounded and admits an optimal solution of identical value:

min n j=1 āij xj + Γz i + n j=1 v ij s.t. z i + v ij ≥ âij xj z i ≥ 0 Chenghao WANG Contribution to robust network optimization v ij ≥ 0 ∀j = 1, . . . , n
Finally, we can then reinsert the minimization problem in our original nonlinear problem to substitute the max term, obtaining the following compact and linear robust counterpart:

max n j=1 c j x j s.t. n j=1 āij x j + Γz i + n j=1 v ij ≤ b i ∀i ∈ I z i + v ij ≥ âij xj ∀i ∈ I, ∀j = 1, . . . , n x j ≥ 0 ∀j = 1, . . . , n z i ≥ 0 ∀i ∈ I v ij ≥ 0 ∀i ∈ I, ∀j = 1, . . . , n,
which has the advantage of being easy to derive and to provide to a state-of-the-art optimization solver like CPLEX or GUROBI for being solved. For an exhaustive formal description of all the passages, we refer the reader to the papers [START_REF] Bertsimasand | Robust discrete optimization and network flows[END_REF] and [START_REF] Bertsimas | The price of robustness[END_REF] which have formally presented the Γ-robustness approach. To conclude this subsection, we highlight that the results illustrated above are at the basis of the original robust optimization results that we present in the chapters that follow.

Chance-constrained Optimization

We finally introduce some fundamentals of Chance-constrained Optimization, referring the reader to the next chapter for a deeper coverage, taking also into account the specific features of the real-world application that was considered. A natural way for including stochastic data in a mathematical optimization problem is represented by the definition of probabilistic constraints, imposing that, for a given value ∈ (0.1), a constraint including stochastic data should be satisfied with a probability at least equal to 1 -. More in detail, if, for example, we refer to a generic linear constraint a x ≤ b, where a is a vector of stochastic data for which we know the corresponding probability distribution, then we would consider the following chance-constrained version of the linear constraint:

P [a x ≤ b] ≥ 1 -
expressing that we identify as feasible solutions those x that satisfy the linear constraint a x ≤ b with probability at least 1 -. The concept of chance-constrained mathematical programming problems is commonly traced back to the seminal work [START_REF] Charnes | Chance-constrained programming[END_REF], which has dealt with the definition of an optimization approach to manage heating oil production while taking into account weather and demand uncertainty. While Chance-constrained optimization can be recognized as a natural way of including stochastic data in mathematical programming, at the same time, it is known to pose a number of computational challenges. Specifically, as discussed in [START_REF] Ben-Tal | Robust optimization[END_REF], it typically leads to problems that are computationally intractable for two main reasons:

1. it may result hard to evaluate with accuracy the probability of the stochastic data appearing in the constraints, even when the probability distribution is "simple"; 2. the feasible set associated with chance-constrained models is typically nonconvex, thus leading to problems that are hard to solve.

To tackle such computational intractability, a way could be constituted by trying to define suitable convex approximations of the problem, as done for example in [START_REF] Klopfenstein | Tractable algorithms for chance-constrained combinatorial problems[END_REF]. However, we remark that these challenges constitute a stimulus to develop new more effective and efficient solutions approaches, as we do in the next chapter.

Robust Flight Level Assignment Problem

Introduction

With the highly increasing demand for commercial flights each year, the air traffic has been heavily increased by around 14% in Europe in 2019 compared to 2014, with a total of over 11.1 million flights [EUROCONTROL, 2019b]. Leaving out the catastrophic 2020, the average annual growth was forecast at 2.0% per year for the next five years [EUROCONTROL, 2019c]. Although the current situation of traffic airspace is greatly underloaded due to the world pandemic situation, one may expect that there will still be a high level of congestion in airspace in a few years, leading to important delays. Among all-cause delays for airlines, the enroute Air Traffic Flow Management (ATFM) delay is a significant cause of delay to airlines. It is still far from the reference values (0.5 minute per flight), though it decreased to 1.6 minutes per flight in 2019, where the total ATFM delays (airport, en-route, and weather delay) reported by airlines decreased to 2.7 minutes. Moreover, the level of delay was the third-worst in the last 10 years, behind 2010 and 2018, with en-route ATFM delays during summer season remaining a problem for airlines [EUROCONTROL, 2019a]. Several solutions have already been proposed to deal with en-route congestion such as reducing the size of control sectors or the distance of separation, while the current Air Traffic Management (ATM) system seems to have reached the structural limits of the system. Apart from the aforementioned approaches, several degrees of freedom on the trajectories can be exploited to regulate the traffic in order to reduce the potential conflicts and hereafter to improve the Air Traffic Control (ATC) capacity, such as re-routing and

Chenghao WANG Contribution to robust network optimization flight level allocation or assignment. Our focus is on flight level assignment, that is assigning each flight to an appropriate level (cruising altitude level), such that the total en-route ATFM delays are reduced. The uncertainty is an important factor to be taken into account when dealing with air traffic issues -there were over 12.5% of flights delayed by an ATFM regulation in 2019 where the average flight departure time delay is 13.1 minutes per flight and about 40% of these flights were delayed by more than 15 minutes [EUROCONTROL, 2019b]. All this pleads for careful modeling of the problem including uncertainty, and the need for robust optimization.

Related works

Optimization problems in ATM in relation to en-route congestion have been widely studied these last decades. Let us cite first the fundamental work of Bertsimas and Patterson, where the Traffic Flow Management Re-routing Problem (TFMRP) [START_REF] Bertsimas | The traffic flow management rerouting problem in air traffic control: A dynamic network flow approach[END_REF]] is formulated and in-depth investigated. The authors show how to optimally control aircraft by re-routing, delaying, or adjusting the speed of the aircraft in the ATC system to avoid airspace regions that have reduced capacities, primarily due to dynamically changing weather conditions. This work has been extended in [START_REF] Bertsimas | The air traffic flow management problem: An integer optimization approach[END_REF]Bertsimas et al. [ , 2011b] ] where different delay causes with respect to all flight phases of a flight have been included in a single optimization problem. Agustín et al. [2012a,b] time savings per flight [START_REF] Li | Evaluating of the benefits of allowing flight level and mach number adjustment for fuel efficiency for flight operations in oceanic airspace[END_REF]. [START_REF] Nace | A linear programming approach for route and level flight assignment[END_REF] proposed a linear programming approach for route and flight-level assignment in a trajectory-based ATM environment. Barnier and Brisset [2004] investigated the problem of flight level allocation considering direct routes only and vertically separate intersecting ones by allocating distinct flight levels. [START_REF] Abad | Using tactical flight level allocation to alleviate airspace corridor congestion[END_REF] proposed en-route flight level allocation for aircraft to mitigate air traffic congestion and airline operating costs in airspace corridors. [START_REF] Chenghao | Applying genetic techniques to the tactical flight level assignment[END_REF] of departures has been included in the solution method to improve the accuracy of the solution.

Contribution and organization

In this paper, we study the robust Flight Level Assignment problem dealing with uncertainty in flight departure time. The contribution of this work is twofold. We investigate first the FLA problem and its Chance-Constrained variant. A compact deterministic mathematical formulation and the Chance-Constrained Programming (CCP) counterpart problem are given, and a heuristic approach is provided.

In essence, the approach stands in separating the problem by level altitudes and solve each of them consecutively in a certain order. A second contribution that can be drawn is the specific approach to solving the CCP FLA problem associated with each flight level. Using the results of Klopfenstein, we reformulate the CCP as an Integer Linear Programming (ILP) one where the probability constraints are replaced by the worst case like Soyster model ones [START_REF] Allen | Convex programming with set-inclusive constraints and applications to inexact linear programming[END_REF]. These constraints are added as they are needed through a constraint generation approach until the desired feasibility probability is reached. All this gives a practical approach that may find application in a class of CCP problems. Another novelty of the approach is the method used to check the feasibility probability of constraints. We consider the case with uncertain parameters following a truncated Gaussian Mixture Model (GMM) distribution and propose an approximated method to estimate the feasibility probability of the solution. This may be of high interest in practice since the GMM is a powerful tool to capture the characteristics distribution of a large number of real situations. The remainder is organized as follows: Section 3.4 reports a brief description of airspace configuration. Section 3.5 presents the compact deterministic mathematical formulation of the FLA problem and its N P-hardness is present. The robust counterpart of the FLA problem taking account of the uncertainty of the flight departure time delay is established in Section 3.6. We report in Section 3.7 the general procedure to tackle the robust FLA problem. This procedure includes a heuristic estimation method is proposed to estimate the feasibility probability of each obtained solution from the corresponding robust subproblem. Section 3.8 reports the computational results used to validate the proposed approaches. Some discussion on the generality of the proposed method is reported together with concluding remarks in Section 3.9.

Airspace configuration

We present here a brief description of airspace configuration including flight mission, separation minima for security, potential conflict, and the computation of induced delay in line with a resolution of pairwise conflict. We assume known or given what follows.

• a set of flights with their traffic trajectory;

• a set of waypoints;

• a set of feasible flight levels for each concerned flight;

• an acceptable upper bound of en-route ATFM delay for each flight;

• the trajectory between two waypoints is straight;

• two aircraft are assumed having a potential conflict at a crossing waypoint if and only if there is a non null probability that the minimum separation distance between them is less than minimum separation (denoted with S);

• a potential conflict between two flights is occurred during the "Cruise Phase" for flights flying at the same level;

• pairwise potential conflicts are assumed independent events;

• flight departure delays are propagated constantly through the flight. Separation Minima Separation minima defines the minimum separation distance between two aircraft for a safety consideration, including vertical separation, lateral separation and longitudinal separation. In this study, for the vertical separation we use the Reduced Vertical Separation Minima (RVSM) which is established between Flight Level (FL)2901 and FL410 in order to increase the airspace capacity. In RVSM airspace, the minimum vertical separation is 10 FL in contrast to 20 FL in non-RVSM airspace. The lateral separation minima describe the minimum separation between aircraft in a horizontal plane such that the spacing between aircraft is never less than a specified amount where lateral separation is applied for aircraft following different tracks while the longitudinal separation is applied for aircraft following the same, converging or diverging tracks. When surveillance systems are used, the minimum separation prescribed in [ICAO, 2016] is 5 nm 2 . With respect to the flight level assignment, the so-called "Semicircular/hemispheric" rule is also applied. The eastbound (respectively, southbound) flights (Magnetic Track 0 to π) use odd FL (e.g., FL290, FL310) and westbound (respectively, northbound) flights (Magnetic Track π to 2π) use even FL (e.g., FL300, FL320) for the airspace of Europe. However, 20 FL intervals are resumed to separate same-direction aircraft and only odd FLs are assigned at FL410 or above, depending on the direction of flight: Magnetic Track 0 to π uses FL410, FL450, etc; Magnetic Track π to 2π uses FL430, FL470, etc.

Potential conflict

We identify four types of pairwise air conflicts due to a loss of separation at cruise phase (see Figure 3.4.3).

• Crossing conflict-may occur if the two aircraft cross at some point o and diverge afterwards.

• Trailing conflict-may occur if two aircraft follow the same route, as it is often the case on airways.

• Converging conflict-may occur if the two aircraft join at some point and remain the same afterwards, at least for a portion of the flight. Two aircraft are involved simultaneously in a crossing and trailing conflict around the waypoint o.

• Diverging conflict-may occur if the two aircraft share the same track and diverge afterwards. Two aircraft are involved simultaneously in a crossing and trailing conflict around the waypoint o 2 .

In practice, all conflicts are solved through predetermined conflict solving procedures. We have analyzed these situations and computed the delay associated with such procedure as described in the following:

En-route ATFM delay ω ij of potential conflict for two aircraft. Figure 3.4.3 shows us geometrically the conflict situation during cruise phase assuming the trajectory between two waypoints is straight. We assume that the potential conflict occurred at the first time for two aircraft is the only one taken into account if there are more than one potential conflict occur between the two aircraft. (indeed having more than one is practically improbable). Let assume that some potential conflict has encountered at waypoint o, and θ is the crossing angle. 3 -En-route potential conflict between two aircraft cruising at same level for sake of simplicity, the delay at other phases of a flight due to the departure delay is omitted, we consider only the en-route ATFM delay which is caused by the flight departure delay of flights. Note t o i , t o j is the time that aircraft i and j passes the conflict point o, respectively. Let t msd ij specify the minimum separation time instead of minimum separation distance for two aircraft to pass safely the conflict point, which is computed exhaustively in Appendix I. The induced en-route ATFM delay ω ij of resolution for pairwise conflict is then formulated by:

ω ij = (t msd ij -t o i + t o j )1(t msd ij -t o i + t o j ) [0,∞) 1(t o i -t o j ) [0,∞) ∈ [0, t msd ij ], (3.1)
where 1(x) A is an indicator function: 1(x) A = {1, x ∈ A; 0, x / ∈ A}, the second term ) specifies that the induced en-route ATFM delay of associated flight due to resolution should be positive which in turn means that there exists a potential conflict, and 1(t o i -t o j ) [0,∞) denotes whether the aircraft i arrives latter than j at the potential conflict point o. The above formula can be expressed as a function of departure times and we assume that the flight departure time delay follows a GMM, then ω ij follows a truncated GMM distribution.

1(t msd ij -t o i + t o j ) [0,∞

The flight level assignment problem

In this section, we show first the N P-hardness of the FLA problem via a 3-Graph (vertex) Coloring (3-GC) problem. A corresponding compact formulation is then described mathematically.

Complexity issues

The complexity of the FLA problem is N P-hard in the strong sense. The proof is based on the 3-GC known to be N P-Complete in the strong sense even for a planar graph [START_REF] Wegener | Complexity theory: exploring the limits of efficient algorithms[END_REF]. Before the details of how 3-GC can be polynomially reduced to the FLA decision problem, we first define formally these problems.

The 3-GC decision problem Given an undirected graph G = (V, E) with V denoting the set of vertices and E specifying the set of edges, a set of 3 colors denoted C and the coloring function c : V → C, is there any assignment of colors from C such that we have c(u) = c(v) for any two adjacent vertices u, v ∈ V ?

The FLA decision problem. Given a set of flights denoted F, a set of eligible flight levels for each flight denoted L i (usually, |L i | = 3, one is the most preferred level, and the other two feasible candidate levels) included in a set L (clearly, |L| ≥ 3) , and a function P l : F × F → {0, 1} (i.e., for each pair (i, j) of flights flying at the same level l ∈ L, P l (i, j) takes value 0 if there is no en-route conflict between them and 1 if there is a potential one), is there any assignment of flights to their eligible levels such that we have P(f, f ) = 0 for any pair (f, f ) of flights the same level? More in specifically, the airspace is defined by a set of waypoints (that is a reference point in the airspace used for purposes of navigation), and we assume that trajectory between two waypoints is straight. Moreover, two flights passing through the same waypoint within a short interval of time are assumed in potential conflict.

Proposition 3.5.1 The FLA problem is N P-complete in the strong sense.

Proof We will show that for any instance of the 3-GC decision problem, we can construct in polynomial time an instance of the FLA decision problem accepting a solution if and only if the 3-GC accepts a solution. Let first show that the FLA decision problem is in N P, that is we can check in polynomial time with respect to the size of the problem instance that a given solution accepts a yes answer to the question [START_REF] Wegener | Complexity theory: exploring the limits of efficient algorithms[END_REF]. It is obvious that a verification of existence of potential conflict for a given assignment of flights to levels can be done in polynomial time. We just need to check if for each pair of flights assigned at the same level there is some potential conflict, which gives at most n(n -1)/2 verifications (n = |F|).

Let us consider in the following the equivalence of both decision problems. For any instance of 3-GC we can construct an instance of FLA as follows: Let be given a graph G = (V, E) with respect to an instance of the 3-GC problem. Let us now define a bijection A between vertices in V and flights in the instance of FLA that we are supposed to construct, that is |F| = |V|. Further, we have three levels in the FLA instance that we will build, which corresponds to the number of colors needed for the graph. Let us note these colors r, g, b and build a correspondence with levels in the FLA problem. Let us use now the vertices and edges in G to build an airspace network with respect to the FLA instance. The airspace network will be composed of |V| flights, |E| waypoints and only 3 flight levels (|L| = 3, for sake of simplicity). For each flight f we have a specific origin and destination airport, denoted respectively f o and f d . The waypoints are identified by a pair of flights corresponding to the extremities of edges in G, for instance some waypoint in level l corresponding to edge e = (v, v ) ∈ E will be denoted with (A(v), A(v )) l = (f, f ) l . We suppose that at each waypoint there is a potential en-route conflict between the corresponding flights when they are at the same level. Then, with each flight f = A(v) at level l we associate a route from f o to f d passing through waypoints corresponding the adjacent links to node v in the initial graph G. The waypoints are traversed in increasing order of the corresponding adjacent nodes in G (e.g., < 1, 2 >). Given an instance of the 3-GC decision problem (see Fig. 3.5.4), the route of all flights involved in corresponding instance of the FLA decision problem are summarized as follows:

• Flight 1:

1 o →< 1, 2 >→< 1, 3 >→< 1, 4 >→ 1 d • Flight 2: 2 o →< 1, 2 >→ 2 d • Flight 3: 3 o →< 1, 3 >→< 3, 4 >→ 3 d • Flight 4: 4 o →< 1, 4 >→< 3, 4 >→ 4 d

Chenghao WANG

Contribution to robust network optimization

2 o 1 o 3 o 4 o 1 d 2 d 3 d 4 d 1,4/1 1,4/2 1,4/3 1,3/1 1,3/2 1,3/3 1,2/1 1,2/2 1,2/3 3,4/1 3,4/2 3,4/3 (b) The corresponding FLA instance Figure 3.5.4 -FLA instance construction from 3-GC
It is straightforward that this construction can be done in polynomial time. Let us show now that solving the 3-GC problem provides a valid assignment, i.e., avoiding conflicts, for the corresponding instance of FLA. Indeed, let c be a valid color assignment to vertices in G. Then, if we choose for each flight f in the airspace network the level l which corresponds to color l used to color vertex A -1 (f ) ∈ V , we will never find two flights in the same level passing through the same waypoint, and hence there are no conflicts. Vice-versa, if there is a solution avoiding conflicts in our airspace network, let associate with flights on level l the color l. By using this coloring for the corresponding vertices in G we obtain a solution for this instance of 3-GC. As the 3-GC problem is N P-complete in the strong sense, this proves that the FLA is as well, concluding thus the proof.

Compact mathematical formulation

Sets and indices: η i : The penalizing cost for a flight cancellation (it happens when the total induced en-route ATFM delays exceeding this maximum value W i ). This is called cancellation in the sense that this flight will not be counted as realized and its potential conflicts with other flights not accounted.

α i : The average unit cost per minute and per flight of en-route ATFM delay.

b l i : The estimated profit corresponding to assignment of flight i at level l. ω ij : The induced en-route ATFM delay of flight i when resolving a potential conflict with j flying at the same level l.

M l

i : A sufficient large number, e.g.

M l i = Σ j∈S l i ω ij . Variables:
x l i : A binary variable taking value 1 if the flight i flies on level l, 0 otherwise. y i : A positive continuous variable indicating the cumulative induced en-route ATFM delay for an assigned flight.

Using the above-mentioned notation, mathematical model associated with the deterministic compact FLA problem denoted by CP is described as below:

max i∈F l∈L i (b l i -α i y i )x l i - i∈F η i (1 - l x l i ) (3.2) s.t. j∈S l i ω ij x l j + M l i (x l i -1) ≤ y i ∀l ∈ L, ∀i ∈ F l (3.3) l∈L i x l i ≤ 1 ∀i ∈ F (3.4) 0 ≤ y i ≤ W i ∀i ∈ F (3.5) x l i ∈ {0, 1} ∀l ∈ L, ∀i ∈ F l (3.6)
The objective function maximizes the assignment profit as well as minimizes the penalizing cost due to cancellation and en-route ATFM delays. Clearly the above objective function can be written as i∈F l∈L i (b l i +η i -αy i )x l i . Constraints (3.3) calculate the cumulative induced delay of each flight i assigned at level l and specify that the flight is canceled if the induced en-route ATFM delay exceeds the maximum acceptable upper bound W i . Constraints (3.4) specify that each flight

Chenghao WANG

Contribution to robust network optimization is assigned at most to one of its eligible levels. Constraints (3.5) and (3.6) define the feasible domain of decision variables. The above model is nonlinear due to the bilinear terms y i x l i , hence a linearization is applied thanks to McCormick Envelopes [START_REF] Garth | Computability of global solutions to factorable nonconvex programs: Part i-convex underestimating problems[END_REF]:

z l i ≥ 0; z l i ≤ W i x l i ; z l i ≤ y i ; z l i ≥ W i (x l i -1) + y i (3.7)
where z l i is the product of y i and x l i .

3.6

The Robust counterpart of FLA problem

Modeling the Robust FLA problem

In a real ATM environment, flight departure delays are usually uncertain due to various sources of nature, such as airport weather conditions, airport capacity, airport disruption, and airport staffing. The Gaussian Mixture Model is a powerful tool to capture main characteristics of departure delay distribution, as investigated by [START_REF] Tu | Estimating flight departure delay distributions-a statistical approach with long-term trend and short-term pattern[END_REF]. We apply the coefficients of aforementioned GMM (present in Table 3.6.1) for the distribution of flight departure time delay in this study, where components C 1 and C 2 capture mainly the negative value of departure delays where C 2 has a higher peak to shape the skewness, component C 3 specifies some medium delays and C 4 accounts for the very large delays. All this makes gives the maximum induced en-route ATFM delay occurring in case of conflict between flights i and j. Furthermore, ω ij follows a truncated GMM distribution since the flight departure time delay is GMM distributed (see Equation (3.1) and Appendix I).

Assuming separate probability conditions and taking account of above uncertainty of flight departure time delay, the compact mathematical formulation of the robust FLA problem denoted by RP can be expressed via the Chance-Constrained Programming as follows:

max i∈F l∈L i (b l i + η i -α i y i )x l i - i∈F η i (3.8) s.t. P   j∈S l i ω ij x l j + M l i (x l i -1) ≤ y i   ≥ 1 - ∀l ∈ L, ∀i ∈ F l (3.9) l∈L i x l i ≤ 1 ∀i ∈ F (3.10) 0 ≤ y i ≤ W i ∀i ∈ F (3.11) x l i ∈ {0, 1} ∀l ∈ L, ∀i ∈ F l (3.12)
where probabilistic constraints (3.9) ensure for each flight that the sum of induced en-route ATFM delays will not exceed the given upper bound of acceptable delay with a probability at least 1-( gives the desired infeasibility (violation) tolerance of constraint (3.3)). For sake of simplicity, we will allow ourselves to use the same notation for ω ij , M l i as in CP, but here ω ij is a random value (truncated GMM distributed) bounded in [0, ω ij ] due to the uncertainty of flight departure delay and

M l i = Σ j∈S l i ω ij .
Clearly the problem is difficult and becomes intractable even for moderate instances. Our solution approach stands in two paradigms: first, we proceed to a heuristic decomposition approach separating the problem per flight level, and then the specific (CCP) FLA problem for each single level is solved through Robust Optimization (RO) methods.

The key issue of above solution approach is how to formulate the associated robust subproblem, called RP l , into an ILP problem and solve it efficiently. This is in the focus of the next section.

Subproblem associated with a single flight level (RP l )

Before detailing the mathematical formulation of this subproblem, let give some precision on the notation. As there is no need to distinguish flight levels, the binary variable x l i is now replaced by x i , and as before it takes value 1 when the flight i flies on level l and 0 otherwise. For sake of simplicity, we will allow ourselves to use the same notation for F l as in CP, but here it groups only flights 1) being eligible to fly at this level and not yet assigned to other flight level, or 2) whose most preferred flight level is the processing one. Notice that y i are now parameters of cumulative induced en-route ATFM delay for the corresponding flight i in RP l , (denoted as vector y F l ). The mathematical formulation associated with the CCP FLA restricted to level l denoted by RP l (F l , y F l ) then follows:

max i∈F l (b l i + η i -α i y i )x i (3.13) s.t. P   j∈S l i ω ij x j + M l i (x i -1) ≤ y i   ≥ 1 - ∀i ∈ F l (3.14) x i ∈ {0, 1} ∀i ∈ F l (3.15)
The above RP l (F l , y F l ) is still a very difficult CCP one. Recall first that CCP is a specific model of stochastic optimization looking to optimize the objective, given an infeasibility probability tolerance. It is a very hard optimization area and one way to tackle it, largely studied these last two decades, is to go through Robust Optimization. In RO, the Soyster model gives a very conservative model preserving the feasibility at all scenarios. The Γ-robustness [START_REF] Bertsimas | The price of robustness[END_REF] model looks for a better trade off between the feasibility and optimality. Concerning the probabilistic constraint (3.14), we have = exp(-Γ 2 i /(2|S l i |)) from Γ-robustness model if the ω ij is independently symmetrically distributed. However, this model is not directly usable because of asymmetric uncertain interval of ω ij . In practice, given the desired value of for each concerned flight i, we can not easily compute the corresponding coefficient Γ i for the involved constraint of each flight i in each associated subproblem in our study. On the contrary, the feasibility probability 1 -for the involved probabilistic constraint can be posteriorly estimated with a realized value of the corresponding Γ i by sampling the possible value of ω ij . A robust solution can be obtained through Γ-robustness, but it may be a time-costly solution due to the high combination of Γ i values for the involved constraints of associated flights and the probability estimation sampling procedure. Hence, to solve the above problem we have opted to use the model introduced in Klopfenstein [2009], which proposes a simple heuristic algorithm to find good solutions to general chance-constrained integer linear problems. In this model, the author has introduced a parameter vector γ ∈ [0, 1] |F l | which allows tuning the robustness of the solution in a convenient way. Applying this idea, we obtain the following

Chenghao WANG

Contribution to robust network optimization alternative model:

max i∈F l (b l i + η i -α i y i )x i (3.16) s.t. min    j∈S l i ω ij x j , γ i j∈S l i ω ij    + M l i (x i -1) ≤ y i ∀i ∈ F l (3.17) x i ∈ {0, 1} ∀i ∈ F l . (3.18)
Let us focus on the robust constraint associated with flight i in (3.17). Either we consider the worst case (maximum induced en-route ATFM delay, ω ij , for all flights j ∈ S l i in conflict with flight i), or we have a constraint:

γ i j∈S l i ω ij +M l i (x i -1
) ≤ y i . In this latter case, two sub-cases occur: when γ i j∈S l i ω ij > y i , then x i = 0; when γ i j∈S l i ω ij ≤ y i , we have a dummy constraint which can be ignored. These three cases are in fact summarized in the two following ones:

• either flight i has cumulative induced en-route ATFM delays less than the given en-route ATFM delay y i and no constraint is necessary to model this situation;

• or flight i is associated with maximum induced en-route ATFM delay with flights j ∈ S l i , i.e., j∈S l i

ω ij x j + M l i (x i -1) ≤ y i .
Hence, the analysis of the above robust model leads to a new one, which is very simple. Indeed, for a given value of γ i we know in advance if the constraint corresponding to flight i is necessary to be put in the model or not. Let denote with I c ⊆ F l a subset of concerned flights with respect to a given vector γ. In this way, instead of vector γ we use the subset I c as a parameter enabling to tune robustness. The corresponding problem denoted by RP l (I c , F l , y F l ) is then formulated as follows:

max i∈F l (b l i + η i -α i y i )x i (3.19) s.t. j∈S l i ω ij x j + M l i (x i -1) ≤ y i ∀i ∈ I c (3.20) x i ∈ {0, 1} ∀i ∈ F l (3.21)
With respect to set I c considered, the size of the above ILP varies between few constraints (initially I c is empty) and all constraints (i.e., I c = F l ). The set I c is said valid if the constraints (3.22) are satisfied for the obtained solution x * of RP l (I c , F l , y F l ).

P   j∈S l i ω ij x * j + M l i (x * i -1) ≤ y i   ≥ 1 - ∀i ∈ F l \ I c (3.22)
A natural approach to solve the subproblem RP l is a constraint generation approach following the strategy proposed in [START_REF] Klopfenstein | Tractable algorithms for chance-constrained combinatorial problems[END_REF]. The idea is to start with an empty set of constraints (I c = ∅), check the obtained and add the most violated one into the subset I c until the feasibility probability of the obtained solution, the constraint (3.22) for each concerned flight, is satisfied.

Algorithm 1: Heuristic method for solving RP l Input: F l : A set of flights eligible for processing level l, y i : The given cumulative induced delay for associated flights f ∈ F l . Output: x * : An optimal solution of flight level assignment In other words, by refining the set I c such that constraint (3.22) is satisfied for obtained solution x * from RP l (I c , F l , y F l ), we then have a certified feasible solution for the associated robust subproblem. As highlighted in Algorithm 1, a solution is obtained initially by an empty set I c . For each obtained solution, the concerned feasibility probability is estimated by the methods detailed in Section 3.6.3. Set I c is updated iteratively by inserting the flight i that the associated constraint (3.22) is the most violated and the associated problem is resolved until the constraints (3.22) for all concerned flights are satisfied. Finally, a robust solution is obtained for the corresponding subproblem RP l , where the robustness (feasibility probability) of the solution is guaranteed by the minimum value of the posteriorly estimated feasibility probability over all associated constraints.

1 Set I c ← ∅; 2 Unsolved ← True; 3 do 4 Solve RP l (I c , F l , y F l );

Estimation of feasibility probability of solution of RP l

Given a solution of an instance of RP l (I c , F l , y F l ), only constraints (3.22) may be violated due to uncertainty of ω ij . Observe first that the constraint is dummy if x i takes value of zero, which is always feasible regardless the value of ω ij . Therefore, we restrict ourselves in ensuring that the feasibility probability of the concerned constraint is at least 1 -(i.e., P( j∈S l i ω ij x j ≤ y i ) ≥ 1 -) for all x i = 1 and i ∈ F l \I c . The flowchart to estimate the feasibility probability of an obtained solution of RP l by different methods is highlighted in Figure 3.6.5, Figure 3.6.6, and Figure 3.6.8. There are three methods used to estimate the feasibility probability which are: the first one is using Hoeffding's Inequality, the second is based on Monte-Carlo simulation and the last is a heuristic estimation method evaluating the sum of random variables following a truncated GMM distribution.

Conservative robust method For a comparison purpose, we consider first the Soyster model [START_REF] Allen | Convex programming with set-inclusive constraints and applications to inexact linear programming[END_REF], which solves a specific deterministic variant instead of seeking the best solution remaining feasible over all possible scenarios. We refer it as Soyster feasibility probability estimation method (Soyster method). The obtained optimal solution is surely feasible over all scenarios as all data takes their worst-case value (i.e., ω ij ). The biggest drawback is that it may lead to a costly solution that strays far away from the optimal one of a given scenario due to the over-conservatism.

Hoeffding's Inequality method Hoeffding's Inequality [START_REF] Hoeffding | Probability inequalities for sums of bounded random variables[END_REF] is a result in probability theory that gives an upper bound on the probability for

Chenghao WANG

Contribution to robust network optimization the sum of identical independent random variables to deviate from its expected value. We hereby refer it as Hoeffding's inequality feasibility probability estimation method (Hoeffding method). Let us recall the Hoeffding's Inequality: Let X 1 , X 2 , • • • X n be the identical independent random variables. Assume that X i are almost surely bounded; that is , for

1 ≤ i ≤ n we have P(X i ∈ [a i , b i ]) = 1
for some finite a i , b i . Let be S = n i=1 X i and E[S] its expected value. Then we have the inequality:

P(S -E[S] ≥ nt) ≤ exp -2n 2 t 2 n i=1 (b i -a i ) 2 ∀t > 0 (3.23)
We apply the Hoeffding's Inequality. Noting that,

P   j∈S l i ω ij x j ≥ y i   = P   j∈S l i ω ij x j - j∈S l i E[ω ij ]x j ≥ y i - j∈S l i E[ω ij ]x j   (3.24)
we obtain:

P   j∈S l i ω ij x j ≥ y i   ≤ exp    - 2 y i -j∈S l i E[ω ij ]x j 2 j∈S l i ω 2 ij x 2 j    = (3.25) When t = y i -j∈S l i E[ω ij ]
x j ≤ 0, the probability P( j∈S l i ω ij x j ≤ y i ) is set to zero. Whereas, in case that y i is bigger than the sum of all upper bounds of random variables, then the probability is surely one. Thus, we obtain a piece-wise probability function as follows: Monte-Carlo simulation method The Monte-Carlo Simulation method is frequently used in mathematical problems such as optimization, generating draws from a probability distribution. We refer it as Monte-Carlo simulation feasibility probability estimation method (MC method). The main idea behind this method is randomly generating a sufficiently large number of scenarios to obtain numerical results. In this study, we randomly generate the departure time delay following a Gaussian Mixture Model (see Table .3.6.1) for the concerned flights in each randomly generated scenario. We have assumed that for each flight, data uncertainty coming from departure time delay will be propagated constantly through the flight trajectory.

P   j∈S l i ω ij x j ≤ y i   =          0, if y i ≤ j∈S l i E[ω ij ]x j 1, if j∈S l i ω ij x j ≤ y i 1 -, otherwise
Note p the probability P j∈S l i ω ij x j ≤ y i , N mc is the number of scenarios for simulation, 1 -α is the confidence level of p, u α is the corresponding quantile such that Φ(u α ) = 1 -α and p is the frequency of event j∈S l i ω ij x j ≤ y i of N mc scenarios. Applying the central limit theorem [START_REF] Rosenblatt | A central limit theorem and a strong mixing condition[END_REF] and Slutsky's theorem [START_REF] Delbaen | A remark on slutsky's theorem[END_REF], we have:

p -p p(1-p) Nmc -→ N (0, 1) (3.27)
Then we deduce that:

Chenghao WANG Contribution to robust network optimization

P   p -p p(1-p) Nmc ≤ u α   ≥ 1 -α ⇔ P   p ≥ p -u α p(1 -p) N mc   ≥ 1 -α (3.28) Note that t = Nmc(1-)+0.5u 2 α Nmc+u 2 α and t = Nmc(1-) 2 Nmc(1-)+0.5u 2 α
, we obtain:

p -u α p(1 -p) N mc ≥ 1 -⇔ p ≥ t + √ t 2 -t t (3.29)
Thus the concerned constraint is said feasible with probability 1 -with a confidence level at least 1-a when testing it for N mc scenarios, and getting constraints (3.22) satisfied for at least N mc p scenarios (where p is given in Table 3.6.2).

* i = 1, i ∈ F l .
The induced ATFM en-route delay is then calculated in each generated scenario, and the number of scenario such that

j∈S l i ω ij x j ≤ y i , x * i = 1, i ∈ F l \ I c is countered.
The feasibility probability of constraint (3.22) for corresponding x * i is hereby obtained.

Heuristic estimation method-Computing an estimation of sum of random values following truncated GMM distribution The Hoeffding's inequality gives us an upper bound of the probability of feasibility of obtained solution while it may be quite weak and lead to a costly solution. The Monte-Carlo simulation provides us a good robust solution than the former method, while the computation time may be extremely higher than the others. Therefore, it strives us to calculate an estimation of the sum of random variables ω ij , based on a data driven approach taking advantage of their known distribution of each random variable ω ij (i.e., following a truncated GMM distribution) while the corresponding computation remains tractable. The distribution for the sum of identical independent random variables is usually unknown or quite difficult to be characterized by a closed-form expression, even if the distribution of each random variable is known. However the sum of identical independent Gaussian Mixture Model distributed variables follows also a GMM, so as the Gaussian distribution, Poisson distribution and Gamma distribution do. Saying that, our intention is to approximate the distribution of uncertain data, i.e., the induced en-route ATFM delay ω ij of flight i for a resolution of potential conflict with flight j at the same level, as a GMM variable which will give a GMM as well for their summation. We refer it as heuristic feasibility probability estimation method (Heuristic method). To do this we need to go through two steps (the algorithm's general scheme is described in Algorithm 2) :

• First, we approximate each ω ij (which are truncated GMM and bounded in [0, ω ij ]), into a GMM one. We do this by applying an approximation operator implemented through a modified-Expectation-Maximization (modified-EM) algorithm (see Algorithm 3).

• For the second step, we need to approximate the sum of GMM, i.e., convolution of GMM. Remind here that the convolution of any two GMMs with K 1 and K 2 components will produce a GMM of K 1 * K 2 components, which leads

Chenghao WANG Contribution to robust network optimization to an extremely large number of components for a sum of n GMM distributed variables. Then, to make the computation tractable, we propose to approximate the sum of two GMMs with K components (which gives a GMM with K 2 components due to convolution) to another GMM with only K components. This comes to an approximation of GMM with K 2 components into an another one with K components. We denote this operation with ⊕. This operation consists in merging components of the convolution according to the Hellinger distance (look at [START_REF] Cutler | Minimum hellinger distance estimation for finite mixture models[END_REF]). Applying this to all sum will give a GMM of K components (see Algorithm 4).

Algorithm 2: The general scheme of approximation of ω ij into a Gaussian Mixture Model of K components Input: The distribution of n random variables ω ij denoted ω p . Each of them corresponds to a truncated GMM of K components.

Output: A Gaussian Mixture Model of K components 1 Compute ωp for all 1 ≤ p ≤ n (ω p is an approximation of ω p following a GMM distribution of K components); 2 Compute ω * = ω1 ⊕ ω2 ⊕ • • • ⊕ω n (⊕
is a merging operation using the Hellinger distance, the operation priority is illustrated in Figure 3.6.7); 3 return ω * , a GMM of K components.

Approximate the distribution ω p by the modified-EM algorithm. The modified-EM algorithm is detailed in Algorithm 3. To start with, the coefficients µ k , σ k , c k for the GMM of K components of given distribution ω p are initialized, and the involved distribution is sampled by a set of nodes X i (i.e., X i ∈ [µ-6σ, µ+6σ], where µ, σ 2 are mean and variance of involved distribution, respectively). We then calculate and normalize the conditional probability P k (X = X i ) of each sampling point under each component k of targeted GMM. Furthermore, the coefficients µ k , σ k and c k of each component k of targeted GMM are updated iteratively according to the conditional probability and sampling points by the equations (3.31). The approximated GMM of K components is obtained by N em iterations.

Approximate the convolution of two GMMs by the K-means algorithm. The convolution of any two GMMs of K 1 and K 2 components will produce a GMM of K 1 * K 2 components, which leads to an extremely large number

p ik = c k P k (X = X i ); p ik = p ik Nsamp i=1 K k=1 p ik (3.30) 4
Update coefficient of GMM: µ k , σ k , and c k :

µ k = Nsamp i=1 p ik X i ; σ 2 k = Nsamp i=1 p ik X 2 i -µ 2 k ; c k = Nsamp i=1 p ik (3.31)
5 while number of iterations < N em ; 6 return an approximated distribution ωp (a GMM of K components).

of components (some of them contribute few for the distribution) for a sum of n GMM distributed variables. Therefore, a merging operator is introduced to reduce the big number of components for a GMM due to convolution, hereby an approximation of GMM with K 2 components into K components such that the distribution of sum of finite GMM distributed variables follows a GMM of K components. The procedure of such approximation is given in Algorithm 4. Initially, we randomly choose K components (i.e., a Gaussian distribution N ) of a GMM with K 2 components as cluster, then for each component of involved GMM (i.e., a Gaussian distribution N ), the Hellinger distance [START_REF] Cutler | Minimum hellinger distance estimation for finite mixture models[END_REF] of each pair of (N , N ) between each component of involved GMM and each cluster in targeted GMM is calculated. Find the nearest cluster N for each N .

For each cluster N , merge all weighted Gaussian distributions N and update the cluster N by equations (3.33). The approximated GMM of K components is then obtained by N km iterations.

An example illustrated the computation for probability of sum of induced delay for an assigned flight is given in Figure 3.6.7. The computation of P(ω 1 + ω 2 + ω 3 + Algorithm 4: K-means algorithm to approximate a Gaussian Mixture Model of K 2 components into the one of K components Input: G:

A GMM of K 2 components Output: G : A GMM of K components 1 Randomly choose K components from G as clusters for G ; 2 do 3 foreach Gaussian distribution N i ∈ G do 4
Calculate the Hellinger distance between N i and cluster N j :

H 2 (N i , N j ) = 1 2 dN i -dN j 2 (3.32) 5
Find the nearest cluster N j for N i ;

6 end 7 foreach cluster N j ∈ G do 8
Let S be Gaussian distributions N i whose nearest cluster is N j ;

9

Update coefficients of the cluster N j : µ j , σ j , and c j :

µ j = N i ∈S c i µ i ; σ 2 j = N i ∈S c i σ 2 i + µ 2 i ; c j = N i ∈S c i (3.33)
10 end 11 while number of iterations < N km ; 12 return the new GMM G .

ω 4 ) is transformed into calculation of P((A(ω 1 ) ⊕ A(ω 2 )) ⊕ (A(ω 3 ) ⊕ A(ω 4 )))
, applying an approximation operator for each ω i , i = 1, 2, 3, 4 and a merging operator for each pair of composed summation.

Therefore the computation of P( j∈S l i ω ij x j ≤ y i ) for corresponding flight i is then tractable and easily calculated by applying the approximation operator (Algorithm 3) and merging operator (Algorithm 4).

The flow chart of estimating the feasibility probability of obtained solution by means of heuristic estimation method is present in Figure 3.6.8. To estimate the feasibility probability of the involved constraints (3.22), we approximate the distribution of summation of independent GMM-approximated random variables by proposed K-Means algorithm. Especially, these GMM-approximated random variables are approximated by proposed modified-EM algorithm from the original random variables, a GMM approximation pool is applied to avoid a redundant approximation for the same variable which may be occurred in the constraint (3.22). The feasibility probability of constraint (3.22) for corresponding x * i = 1 is hereby tractable and obtained as the final evaluated term in constraint (3.22) approximately follows a GMM distribution.

Probability: P((A(ω1) ⊕ A(ω2)) ⊕ (A(ω3) ⊕ A(ω4)) ≤ yi) Approximate distribution ((A(ω1) ⊕ A(ω2)) ⊕ (A(ω3) ⊕ A(ω4))) Approximate distribution (A(ω1) ⊕ A(ω2)) Approximate distribution A(ω1) Distribution ω1 A() Approximate distribution A(ω2) Distribution ω2 A()

Putting all pieces together

A general approximation framework to tackle the involved robust FLA problem is reported in Figure 3.7.9. Clearly the Algorithm 1 is the key brick of the solution approach.

The order examination of flight levels is done starting from the maximum loaded one (with the largest number of flights whose the most preferred level is the processing one) and so on. Then y i is initialized as 0 for each concerned

j∈S l i ω ij x j > y i ; 2.
j∈S l i ω ij x j ≤ y i , but there may be some flights k ∈ S l i such that x k = 1, j ∈S l k \{i} ω kj x j ≤ y k and j ∈S l k \{i} ω kj x j + ω ki > y k .

For the first case, we increment y i for the unassigned flight i. For the second case, we increment y k for the associated assigned flight k rather than the unassigned flight i.

The modified heuristic feasibility probability estimation method

One may argue that the aforementioned robust FLA procedure involving the proposed heuristic feasibility probability estimation method can not guarantee that the posterior solution feasibility probability is larger than the desired one (1-).

To assure that the obtained solution is robust with a probability at least 1-, we repeat the aforementioned robust FLA procedure with a higher probability target. More precisely, if at k th iteration, the solution feasibility probability (denoted by P k ) has not the desired probability, we require a higher probability feasibility to

Chenghao WANG Contribution to robust network optimization the Heuristic method for the concerned constraints in the next iteration of the robust FLA procedure. We refer to this method as the modified heuristic feasibility probability estimation method (HeuMod method).

Computational results

Numerical experiments are reported in this section. All algorithms were implemented in C++. Experiments were carried out on a server with Intel Xeon Gold 6138 2.0-GHz CPU and 125 GB of RAM under Linux. In this study, the same aircraft of Airbus A320 is considered for all flights. The number of available flight seats is 180, and the Revenue per Available Seat Kilometer (RASK) is 0.681 euros [START_REF] Airfrance | Airfrance operating result -full year[END_REF] 3 . The en-route ATFM delay cost per minute per flight (α i ) is set to 70 euros if the en-route ATFM delay (y i ) is less than 30 minutes, and the cost is doubled for each supplementary 30 minutes (i.e., α i = 140 for 30 < y i ≤ 60, and so on). The cancellation cost η i for a flight i (when y i > W i , the flight i is then considered as canceled) is set to 10 times of b 0 i (the cancellation cost here is much larger than the assignment profit at their most preferred level as we want to minimize first the number of canceled flight, and the number of flight changing their level in order to maximize the total (robust) revenue), where b 0 i indicates the estimated profit at their most preferred level. Posterior feasibility probability 39.94% 38.98% 38.44% 38.95%

Test instances

3 In our calculations we have not involved different costs of running the flight as those connected to personnel, airports, airspace use, kerosene, etc.

4 MPM: Maximum Permitted Miles Table 3.8.4 indicates the posterior feasibility probability of all flights assigned at their most preferred flight level under the uncertainty of departure delays. Such assignment leads to a large number of flights (more than 60%) that become infeasible under real conditions involving uncertainty of flight departure time, hereby a costly maintainable solution may be induced to deal with data uncertainties. Table 3.8.6 presents the average percentage of flights changing level (CHA%), the average percentage of canceled flights (CAN%), the average induced cumulative en-route delay per flight (ỹ i ) and the average percentage of gain of total revenue (GanRev%=(total revenue A-total revenue B)/total revenue B*100%, where total revenue A is obtained by Robust FLA with a corresponding configuration and total revenue B is obtained by assigning all flights to their most preferred level with zero tolerance of infeasibility) when the robust flight level assignment is taken into account. The average delay cost per flight and the average total induced cost per flight including en-route ATFM delay, level changing and cancellation are also reported. The tolerance of solution infeasibility are parameterized as 1%, 2%, 3%, 4%, 5%, 10%, 15%, 20% and 25%. For each subproblem involved at a processing flight level, the Soyster model is used to deal with the flight departure delay uncertainty where all random values take their worst-case value. For all test instances, we have that: a) the value ỹi are significantly reduced (ỹ i is up to 0.29 minute per flight by applying Robust FLA within Soyster model) for different configuration of maximum acceptable en-route ATFM delay per flight (W i ), compared to ỹi in the obtained solution in Table 3.8.5. Moreover, the associated average delay

Effectiveness of robust flight level assignment

Chenghao WANG

Contribution to robust network optimization cost per flight is drastically reduced and the average total induced cost per flight is reduced from 3 to 7 times; b) there may be until 23.97% (CHA%=number of flights changing level/total number of flights*100%) flights changing their flight level and may be until 0.49% flight canceled (CAN%=number of canceled flights/total number of flights*100%) for their induced en-route ATFM delay exceeding the maximum acceptable en-route ATFM delay (i.e., y i > W i , W i = 5 minutes), 0.12% and 0.16% for W i = 10 and 15 minutes, respectively; d) there may be until 46.33%, 51.43% and 50.39% gain of total revenue for W i = 5, 10 and 15 minutes, respectively. Therefore, the robust flight level assignment can significantly reduce the en-route ATFM delay also its total induced cost per flight. The gain of total revenue by Robust FLA can vary between 20% and 50%.

Validation of the heuristic feasibility probability estimation method

Table 3.8.7 presents the comparison results in terms of the gain of total revenue (GainRev%) of Robust FLA, average percentage of flights changing level (CHA%), average percentage of canceled flights (CAN%), average en-route ATFM delay per flight (ỹ i ), average delay cost per flight and average total induced cost, CPU time, minimum percentage of difference between posterior probability and desired feasibility probability (Min Diff.Proba.% = Minimum of Solution feasiblity probability-(1-)

1- * 100%) and the maximum percentage for all three W i (5, 10 and 15) and all infeasibility tolerance (1%, 2%, 3%, 4%, 5%, 10%, 15%, 20%, 25%) configurations, where Hoeffding method, MC method (N mc = 10000), proposed Heuristic method (N em = 10, N km = 10), and the modified heuristic feasibility probability estimation method are applied. We find that the robust FLA with Soyster model performs not well in terms of gain of total revenue and gives the largest induced cost per flight due to the over-conservatism of uncertainty on ω ij , which requires more flights canceled or their flight level changed to satisfy the constraints (3.22). The robust FLA using Hoeffding method to estimate feasibility probability is slightly better than the robust FLA with Soyster method, but spends roughly 3 times more computation time. The robust FLA with MC method performs significantly better than the two former, however enormous computation time is required to obtain such solution with certified feasibility probability. The robust FLA with heuristic feasibility probability estimation method performs much better in terms of gain of total revenue, the en-route ATFM delay, its induced cost and computation time. Still, we notice that the solution feasibility probability has an error of at most 1.83% for their lower bound value compared to the desired feasibility probability. For example, the posterior solution feasibility probability is at least 99%*(1-1.83%) = 97.19%, where the infeasibility tolerance ( ) is 1%. More specifically, this approximation error has its roots in the approximation operator and the merging operator. A trade-off between the robustness of solution and the total revenue (the operational cost could be level change, flight cancellation and cumulative en-route delay) is applied for the robust flight level assignment problem by the proposed modified heuristic estimation method in case when the posterior solution feasibility probability must be at least 1-for each concerned probabilistic constraint. This method can have a quasi-certified solution compared to the Monte-Carlo simulation method but requires much less computation time.

Comparison results of different configurations of W i

Table 3.8.8 presents the computational result for robust FLA with different configurations of W i (i.e., W i =5, 10, and 15 minutes) and for different probability estimation method for probabilistic constraints (3.22). For all probability estimation methods, the solution obtained from Robust FLA have higher gain in terms of objective value and lower induced cost per flight (payed to flight level changing,

Chenghao WANG Contribution to robust network optimization the robust FLA with the proposed heuristic probability estimation method always performs the best (in terms of gain of total revenue, number of flight changing flight level, number of canceled flights, the en-route ATFM delay, its induced cost and computation time) within an error of solution feasibility probability less than 2%.

Comparison results of different configurations of when W i = 5 minutes Table 3.8.9-3.8.10 show the computational result for Robust FLA of different configurations of solution infeasibility tolerance ( ) (i.e., =1%, 2%, 3%, 4%, 5%, 10%, 15%, 20% and 25%) when W i = 5 minutes for different probability estimation method for probabilistic constraints (3.22). The Robust FLA has higher gain in terms of total revenue and lower induced cost per flight when the desired solution feasibility probability is lower for all probability estimation methods. Focusing on = 25%, when MC method, Heuristic method or HeuMod method is applied, we can observe there is no flights cancelled, the robust FLA achieve more than 40% gain in terms of total assignment revenue by a few flight changing their level, compare to the solution where all flights are assigned to their most preferred level with zero tolerance of infeasibility (i.e., =0%). When we allow some flight to change their level, and apply a large tolerance of infeasibility such as 25%, the total assignment revenue can be significantly improved. When focusing on the proposed heuristic feasibility probability estimation method, the robust FLA can assign the flight into an appropriate level to make all flight conflict-free (i.e., ỹi =0) by changing at most 4.15% flights without cancellation for ≥ 10%. However, when looking at the Robust FLA with the MC method we notice that it requires more flights changing and cancellations to make all flights conflict-free for the solution infeasibility tolerance not less than 15%, while the computation time is much higher compared to the proposed Heuristic method. Concerning the error of posterior solution feasibility probability (i.e., less than 2%), it is found that the largest error comes from the configuration of at 4%. This error of posterior solution feasibility probability is actually less than 1% for all other configurations of excepting for =4%. Comparison results of lower bound of solution feasibility probability for different methods when W i = 5 minutes Figure 3.8.10 presents the lower bound of posterior solution feasibility probability for different methods when W i = 5 minutes. The curve of posterior solution feasibility probability for the proposed heuristic feasibility probability estimation method is below the curve of desired solution feasibility probability (orange curve) excepting for 25%. Moreover, when = 4%, the error between the posterior probability for Heuristic method and the target one reaches the highest one (roughly 2%) as already discussed in the previous section. With respect to the Hoeffding's inequality feasibility probability estimation method, we can see that the curve of posterior feasibility probability (red curve) is always close to 99%-100% regardless of the desired feasibility probability. Concerning the Monte-Carlo simulation feasibility probability estimation method, the curve of posterior feasibility probability (blue curve) is always above the target curve which guarantees the robustness of solution. One issue may be argued, the proposed heuristic probability estimation can not guarantee the lower bound of the posterior solution feasibility probability as there is an error less than 2%. The obtained solution is not as robust as what we expect. Therefore, a modified variant which penalizes the estimated probability once the solution is not sufficiently robust and we repeat the Robust FLA until the solution becomes totally robust. The modified variant of the proposed Heuristic method now can guarantee the lower bound of the posterior solution feasibility probability (i.e., the green curve is always above the target curve) and has a better solution (in terms of total assignment revenue, induced cost per flight and computation time) than the MC method.

Conclusion and discussion

As conclusion, we may say that the robust FLA problem expressed as CCP, is a difficult one. All this justifies the use of heuristics as developed in the paper. The proposed approach, inspired from the work of Klopfenstein, yields a general CCP approach for a specific class of ILP problems. Suppose it is given an ILP (or LP) with uncertain coefficients. According to Soyster model, one can easily write the corresponding robust optimization problem by putting each constraint in its "worst" scenario version. On way to remedy the very conservative nature of the Soyster model, is not to put all constraints in the beginning but to start with a few of them and add consecutively the others as they are needed, while keeping an eye at the feasibility. Obviously this leads to a constraint generation approach working as below:

1. Start with the master problem P. Solve P and denote with x * the obtained solution.

2. Check the feasibility probability of x * with respect to the desired feasibility probability p (p is the smallest one of the posterior probability for all associated probabilistic constraints).

Chenghao WANG Contribution to robust network optimization 3. If p < 1 -, and c is the constraint achieving it, then add c to P. Repeat above procedure until p ≥ 1 -.

Clearly the above procedure is extremely simple. Still, it is naturally deduced as well from Klopfenstein as it is for our application. It can be especially usable for robust binary linear programming with if-else constraints. Such constraints raise generally in N P hard problems and are formally represented by a constraint involving binary variables and big M. All this yields difficult optimization problems which often become intractable when robustness has to be considered. In our focus are conditions of type: «if cond then thenexp else elseexp endif » where «cond» is a linear formulation involving uncertain coefficients while « thenexp » and « elseexp » are related to a decision binary variable, for instance thenexp ∼ x ≤ 1 and elseexp ∼ x = 0. Then we will need to add M (x -1) in the left hand of the condition. For instance if we have:

A 11 x 1 + A 12 x 2 ≤ B =⇒ x ≤ 1 else x = 0.
This may be expressed as: A 11 x 1 + A 12 x 2 + M (x -1) ≤ B, which are similar to (3.20) and where the same arguments may be used to strive a simplified version of Klopfenstein's model. We intend to investigate this issue in detail in the future.

Concerning the feasibility probability estimation method, the proposed heuristic approach pursues approximating the distribution of a summation for the independent asymmetrically distributed random variables into a GMM distribution by the proposed approximation operator and merging operator. The modified variant with penalization (for the estimated probability) of the proposed heuristic approach has a guarantee of the posterior solution feasibility probability, which means the obtained solution is robust as expected. This method is more efficient to get a good robust solution compared to the Hoeffding method while it remains much faster than the MC method.

Resource Allocation in 5G Superfluid Wireless Networks

Introduction

With the exponential rise of mobile users, overall mobile data traffic is expected to grow to 49 exabytes 1 per month by 2021, a seven-fold increase over 2016 [Cisco Mobile, 2017]. This is the natural result of the evolution and incredible expansion that mobile networks have known since the beginning of the new millennium: from the 1980s, new generations of mobile network technologies have been proposed to continuously offer higher speed, greater capacity and new innovative sets of services as presented in Figure 4.1.1. The 1 st generation of wireless networks (1G) has granted the first (expensive) mobile services, offering voice-only calls by means of bulky cellular phones which, however, were considered a status symbol of successful business people. Then, the 2G appeared in 1991 and allowed a more democratic and cheaper access to cellular phones and services, better supporting texting. Since 1998, the 3G, besides voice calls (wireless connection can be achieved at high speed from 20 Kbps to 42.2Mbps), introduced mobile access to the internet and, during its technological life cycle, smartphones appeared. The 4G technologies starting in the 2008, offered enhanced performance especially for smartphones with a peak data rate of 100 Mbps (4G LTE), 150 Mbps (4G LTE Cat.4), and 1000 Mbps (4G LTE Advanced). During the 4G era, an impressive surge in the number of mobile users, asking for higher performance, has taken place. This has led to the concrete need for a 5G that could offer unprecedented performance and features. "Ubiquitous connectivity", "Zero latency"and "High-speed Gigabyte connection", [START_REF] Gohar | The role of 5g technologies in a smart city: The case for intelligent transportation system[END_REF] introduced in [START_REF] Panwar | A survey on 5g: The next generation of mobile communication[END_REF], are considered to be essential features that a 5G technology should support. Ubiquitous connectivity can be essentially described as the capacity of granting connectivity to every device, everywhere, at every time; we remark that this question may result very challenging in 5G networks, which should effectively support the Internet-of-Things thus requiring to grant simultaneous wireless connectivity to an extremely high number of devices. Zero latency refers to the capacity of a 5G network of supporting null or negligible time lag during data transfer; this is considered particularly critical in applications such as telemedicine (e.g., telesurgery), which requires to remotely intervene on patients in real-time without suffering from negative effects of latency that could threat the life of a patient. High-speed Gigabyte refers to the capacity of a 5G network to establish wireless connection characterized by extremely high speed that is currently far beyond the possibility of 4G networks. To achieve this 5G network, lots of novel architectures (e.g., Superfluidity [START_REF] Bianchi | Superfluidity: a flexible functional architecture for 5g networks[END_REF], SELFNET [START_REF] Jiang | Intelligent network management for 5g systems: The selfnet approach[END_REF], Flexible Functional Split [START_REF] Harutyunyan | Flexible functional split in 5g networks[END_REF]) and new technologies (e.g., mMIMO (massive Multiple-Input Multiple-Output) [START_REF] Erik G Larsson | Massive mimo for next generation wireless systems[END_REF], mmWave (millimeter Wave) and beam-forming [START_REF] Roh | Millimeterwave beamforming as an enabling technology for 5g cellular communications: Theoretical feasibility and prototype results[END_REF]) are proposed by active researchers.

Among the different investigated 5G architectures, it is worthwhile to mention 5G Superfluid [START_REF] Bianchi | Superfluidity: a flexible functional architecture for 5g networks[END_REF] (SF) architecture. The goal of 5G Superfluid architecture is to design a new 5G network architecture, which ensures the required levels of flexibility, agility, portability and high performance. In a nutshell, 5G Superfluid architecture aims to achieve a superfluid state of the network, which is the ability to instantiate services on-the-fly, run them anywhere in the network (core, aggregation, edge) and shift them transparently to different locations. The key brick of the 5G Superfluid architecture is the definition of the concept of Reusable Functional Block (RFB), which is a virtualized entity, used to decompose network functions and services, and it is deployed on top of a physical node. The deployment of the 5G network through the RFBs have notable features, including: i) the possibility to build chain of RFBs, in order to implement more complex functionalities and to provide different services to users; ii) the independence of the RFBs from a specific platform, i.e., RFBs can be realized via software functions, and they can be run on several hardware architectures; and, iii) the introduction of high levels of flexibility and performance, thanks to the fact that the RFBs can be deployed where and when they are really needed. Furthermore, the RFB concept is a generalization of the Virtual Network Function (VNF) concept proposed by [ETSI, 2014]. In particular, RFBs can be arbitrarily decomposed in other RFBs, while VNFs in the ETSI model cannot be composed or decomposed in other VNFs in a flexible way (see Figure 4.1.2 for their different class diagram). Moreover, the RFBs can be mapped to different software and hardware execution environments (see [START_REF] Bianchi | Superfluidity: a flexible functional architecture for 5g networks[END_REF]), while the ETSI model focuses on mapping VNFs to Virtual Machines (or Containers) in traditional cloud infrastructures. [START_REF] Chiaraviglio | Optimal design of 5g superfluid networks: Problem formulation and solutions[END_REF]. Since 5G-RFB-RA is N P-hard, heuristic algorithms were proposed to solve it efficiently in [START_REF] Chiaraviglio | Algorithms for the design of 5g networks with vnf-based reusable functional blocks[END_REF]. The problem was also tackled by a Particle Swarm Optimization algorithm in [START_REF] Shojafar | P5g: A bio-inspired algorithm for the superfluid management of 5g networks[END_REF]. However, one drawback was represented by the fact that some users could obtain limited downlink bandwidth. In order to improve the capacity of solving instances of larger size associated with 5G Superfluid network design problems, we analyzed the polyhedral structure of the original model and propose new more efficient solution methods. Specifically, our original contributions can be summarized as follows:

• we propose an alternative formulation for the problem of minimizing the installation costs of a 5G SF network, taking as reference the model proposed in [START_REF] Chiaraviglio | Optimal design of 5g superfluid networks: Problem formulation and solutions[END_REF]];

• we strengthen the above proposed model, identifying new valid inequalities;

• we propose a hybrid Benders decomposition approach to tackle the simplified model.

Results of computational tests show that our new improved modeling and solution approach offers a higher performance both in computational time and used memory.

The remainder of this Chapter is organized as follows: in Section 4.2, we report a short discussion on 5G Superfluid architecture. The reference model and its improved version are proposed in Section 4.3. We then introduce a Benderslike decomposition approach to solve the problem more efficiently in Section 4.4. Finally, in Section 4.5, we report and discuss computational results. For an exhaustive introduction to Superfluid networks and Reusable Functional Blocks, we refer the reader to [START_REF] Bianchi | Superfluidity: a flexible functional architecture for 5g networks[END_REF]]. Here we recall some main facts about them. A Reusable Functional Block in 5G SF architecture is a logical entity that performs a set of functionalities and has a set of logical input/output ports. In general, a RFB can hold state information and can be combined with other RFBs to form other more complex and performing RFBs (see An RFB results analogous to a traditional VNF or VNFC, implemented as a fully-fledged VM running on a hypervisor or in a container. An RFB can correspond to a small footprint Unikernel VM running in a specialized hypervisor. RFBs can also be modules or components of special purpose execution environments, like extended finite state machines based on OpenFlow for packet processing, software routers, or radio signal processing chains. Hence the RFB concept can be applied to different heterogeneous environments, according to their specification as reported in [START_REF] Bianchi | Superfluidity: a flexible functional architecture for 5g networks[END_REF]. In this study, we consider the RFBs that perform specific tasks in the network architecture, such as processing the video to users, or per-

Chenghao WANG

Contribution to robust network optimization [START_REF] Bianchi | Superfluidity: a flexible functional architecture for 5g networks[END_REF] forming networking and physical layer tasks [START_REF] Chiaraviglio | Optimal superfluid management of 5g networks[END_REF]. Focusing on the tasks realized by these RFBs, the following RFB are identified:

• Resource Radio Head RFB (RRH-RFB): it is in charge of providing the physical signal to the users. Specifically, it handles a set of Radio Frequency (RF) channels established with users and the corresponding baseband channels with the BBU-RFBs;

• Mobile/Multi-access Edge Computing RFB (MEC-RFB): it is in charge of managing an amount of traffic, such as the provisioning of a HD video service to users.

• Base Band Unit RFB (BBU-RFB): it acts as a middle interface between the RRH-RFBs and the MEC-RFBs. Specifically, the BBU-RFB exchanges an amount of IP traffic with the MEC-RFBs, and a baseband signal with the RRH-RFBs.

Using the above-mentioned notation, assuming that each RRH-RFB is compatible with the same type BBU-RFB and the same type of MEC-RFB, and a user is compatible with all types of RRH-RFB, the existing model presented in [START_REF] Chiaraviglio | Optimal design of 5g superfluid networks: Problem formulation and solutions[END_REF] (denoted by Chiaraviglio's Model) can be rewritten as follows:

min n∈N q∈Q c RRH q y RRH nq + n 1 ∈N n 2 ∈N q∈Q c BBU q v BBU n 1 n 2 q + c MEC q v MEC n 1 n 2 q (4.1) s.t. n∈N x un ≤ 1 ∀u ∈ U (4.2) u∈U n∈N x un ≥ δ|U | (4.3) q∈Q y RRH nq ≤ 1 ∀n ∈ N (4.4) x un ≤ q∈Q y RRH nq ∀u ∈ U, ∀n ∈ N (4.5) u∈U x un ≤ q∈Q U max q y RRH nq ∀n ∈ N (4.6) n∈N y RRH nq ≤ A RRH q ∀q ∈ Q (4.7) n 1 ∈N n 2 ∈N v BBU n 1 n 2 q ≤ A BBU q ∀q ∈ Q (4.8) n 1 ∈N n 2 ∈N v MEC n 1 n 2 q ≤ A MEC q ∀q ∈ Q (4.9) y RRH n 2 q ≤ n 1 ∈N v BBU n 1 n 2 q ∀n 2 ∈ N, ∀q ∈ Q (4.10) v BBU n 1 n 2 q ≤ y RRH n 2 q ∀n 1 , n 2 ∈ N, ∀q ∈ Q (4.11) y RRH n 2 q ≤ n 1 ∈N v MEC n 1 n 2 q ∀n 2 ∈ N, ∀q ∈ Q (4.12) v MEC n 1 n 2 q ≤ y RRH n 2 q ∀n 1 , n 2 ∈ N, ∀q ∈ Q (4.
φ un ≤ t u ∀u ∈ U, ∀n ∈ N (4.16) φ un ≥ t u -(1 -x un )CAP max u ∀u ∈ U, ∀n ∈ N (4.17) CAP max u = max q∈Q,n∈N {CAP unq } ∀u ∈ U (4.18) u∈U CAP unq θ unq ≤ CAP RRH q ∀n ∈ N, ∀q ∈ Q (4.19) θ unq ≤ x un ∀u ∈ U, ∀n ∈ N, ∀q ∈ Q (4.20) θ unq ≤ y RRH nq ∀u ∈ U, ∀n ∈ N, ∀q ∈ Q (4.21) θ unq ≥ x un + y RRH nq -1 ∀u ∈ U, ∀n ∈ N, ∀q ∈ Q (4.22) u∈U n 2 ∈N ϕ un 1 n 2 q ≤ CAP MEC q n 2 ∈N v MEC n 1 n 2 q ∀n 1 ∈ N, ∀q ∈ Q (4.23) ϕ un 1 n 2 q ≤ CAP max u v MEC n 1 n 2 q ∀u ∈ U, ∀n 1 , n 2 ∈ N, ∀q ∈ Q (4.24) ϕ un 1 n 2 q ≤ φ un 2 ∀u ∈ U, ∀n 1 , n 2 ∈ N, ∀q ∈ Q (4.25) ϕ un 1 n 2 q ≥ φ un 2 -(1 -v MEC n 1 n 2 q )CAP max u ∀u ∈ U, ∀n 1 , n 2 ∈ N (4.26) y RRH n 1 q + y RRH n 2 q ≤ 1 ∀q ∈ Q, (n 1 , n 2 ) ∈ CONF q (4.27) q∈Q n 2 ∈N v BBU n 1 n 2 q ≤ q∈Q y RRH n 1 q ∀n 1 ∈ N (4.28)
Chenghao WANG Contribution to robust network optimization

q∈Q n 2 ∈N v MEC n 1 n 2 q ≤ q∈Q y RRH n 1 q ∀n 1 ∈ N (4.29) t u ≥ t MIN x un ∀u ∈ U, ∀n ∈ N (4.30)
x un ∈ {0, 1} ∀u ∈ U, ∀n ∈ N (4.31)

y RRH nq ∈ {0, 1} ∀n ∈ N, ∀q ∈ Q (4.32) v BBU n 1 n 2 q ∈ {0, 1} ∀n 1 , n 2 ∈ N, ∀q ∈ Q (4.33) v MEC n 1 n 2 q ∈ {0, 1} ∀n 1 , n 2 ∈ N, ∀q ∈ Q (4.34) t u ≥ 0 ∀u ∈ U (4.35) φ un ≥ 0 ∀u ∈ U, ∀n ∈ N (4.36) θ unq ≥ 0 ∀u ∈ U, ∀n ∈ N, ∀q ∈ Q (4.37) ϕ un 1 n 2 q ≥ 0 ∀u ∈ U, ∀n 1 , n 2 ∈ N, ∀q ∈ Q (4.38)
The objective function minimizes the total installation cost of a 5G SF Network. Constraint (4.2) specifies that each user is served by at most one node. Constraint (4.3) assures that he minimum number of users has to be served. Constraint (4.4) indicates that at most one type of RRH-RFB can be installed in each node. Constraint (4.5) denotes that if the node is serving a user, a RRH-RFB then has to be installed on the node. Constraint (4.6) guarantees that the number of users served by each RRH-RFB is then bounded by the maximum number of users that can be supported by that RRH-RFB. Constraint (4.7)-(4.9) specify that the number of installed RFB module of type q ∈ Q is then bounded by its maximum available number. Constraints indicate that for each node, the total capacity provided to the connected users has to be lower than the maximum total capacity managed by a RRH-RFB of type q ∈ Q. (4.23)-( 4.26) denote that the total traffic from users connected to the RRH-RFB placed at node n 2 ∈ N has to be lower than the maximum total capacity managed by a MEC-RFB of type q ∈ Q in the chain. Constraint (4.27) specifies that if a pair (n 1 , n 2 ) ∈ CONF q , then at most one RRH-RFB of type q can be installed either in n [START_REF] Chiaraviglio | Optimal design of 5g superfluid networks: Problem formulation and solutions[END_REF] The problem solved in [START_REF] Chiaraviglio | Optimal design of 5g superfluid networks: Problem formulation and solutions[END_REF] can be seen as composed of three parts (see Figure 4.3.5) : A) User-Node Assignment problem that optimally decides which user is served by which 5G node; B) RRH-RFB Allocation problem that optimally allocate the RRH-RFB at 5G node in order to serve the connected users; C) RFB Chain Construction problem that optimally places the BBU-RFBs and MEC-RFB to provide a complete RFB chain to satisfy the traffic demanded by users. However, there is an extremely high number of combinations of RFBs in the solution space, in particular concerning the RRH-RFB allocation problem and RFB chain construction problem. In order to tackle this, we propose to represent in an alternative way the definition of RFB chains as presented in the next section.

A New Mathematical Model

New variables:

x k nq : A binary variable taking value 1, if an RFB module k is placed on node n serving a RFB chain in type q; 0, otherwise.

y RRH

unq : A binary variable taking value 1, if an RRH-RFB is installed on node n serving an user u with a RFB chain in type q; 0, otherwise.

y MEC

unq : A binary variable taking value 1, if a MEC-RFB is installed on node n serving an user u with a RFB chain in type q; 0, otherwise.

t MEC

unq : A continuous variable indicating the amount of downlink traffic served to user u at node n with a MEC-RFB of type q, defined as t MEC unq = t u y MEC unq . Using the same notation, but the different decision variables, the compact model of RFB-based Resource Allocation problem under 5G Superfluid wireless Networks denoted by 5G-RFB-RA-CP1 can be then described as below:

min n∈N k∈K q∈Q c k q x k nq (4.39) s.t. q∈Q x k nq ≤ 1 ∀n ∈ N, ∀k ∈ K (4.40) n∈N x k nq ≤ A qk ∀k ∈ K, ∀q ∈ Q (4.41) x RRH n 1 q + x RRH n 2 q ≤ 1 ∀q ∈ Q, ∀(n 1 , n 2 ) ∈ CONF q , n 1 = n 2 (4.42) y RRH unq t MIN ≤ x RRH nq CAP unq ∀u ∈ U, ∀n ∈ N, ∀q ∈ Q (4.
t MEC unq ≤ t u ∀u ∈ U, ∀n ∈ N, ∀q ∈ Q (4.57) t MEC unq ≤ y MEC unq CAP max u ∀u ∈ U, ∀n ∈ N, ∀q ∈ Q (4.58) t MEC unq ≥ t u + (y MEC unq -1)CAP max u ∀u ∈ U, ∀n ∈ N, ∀q ∈ Q (4.59) u∈U q∈Q t MEC unq ≤ q∈Q CAP MEC q x MEC
x k nq ∈ {0, 1} ∀n ∈ N, ∀k ∈ K, ∀q ∈ Q (4.63) y RRH unq ∈ {0, 1} ∀u ∈ U, ∀n ∈ N, ∀q ∈ Q (4.64) y MEC unq ∈ {0, 1} ∀u ∈ U, ∀n ∈ N, ∀q ∈ Q (4.65) t MEC unq ≥ 0 ∀u ∈ U, ∀n ∈ N, ∀q ∈ Q (4.66) t u ≥ 0 ∀u ∈ U (4.67)
The objective function minimizes the total installation cost of a 5G SF Network. Constraints (4.40)-(4.41) specify each RFB module can be installed on a 5G node at most one of its type, and have an upper bound for available number of over the network in each type q. Constraint (4.42) denotes the conflict limitation of RRH-RFB installation in each type q due to radio interference. Constraints (4.43)-(4.45) assure that a user is served by at most one activated RRH-RFB placed on node n in type q, and only if the provided radio link capacity (CAP unq ) by this activated RRH-RFB to user is larger than the minimum required traffic demand (t M IN ). Constraint (4.46) specifies that activated RRH-RFB on node n has an upper bound for number of served users (U max q ). Constraint (4.47) guarantees a minimum required coverage of users by activated RRH-RFBs in the network. Constraint (4.48) specifies that total radio link capacity of served users by an activated RRH-RFB on node n is then bounded by the maximal managed radio link capacity of this RRH-RFB. Constraints

Simplification of the proposed Model

In an optimal solution of the model 5G-RFB-RA-CP, the traffic of each user is set to t MIN , as the objective function aims at minimizing minimize the total construction cost of a 5G SF network. Hence, constraints (4.56)-(4.60) can be rewritten as:

u∈U q∈Q y MEC unq t MIN ≤ q∈Q CAP MEC q x MEC nq ∀n ∈ N (4.68)
Furthermore, given decision value of x MEC nq , then x MEC nq = 1, ∀n ∈ N , we have:

u∈U y MEC unq ≥ 1 ∀n ∈ N, ∀q ∈ Q (4.69) n∈N y MEC unq ≤ 1 ∀u ∈ U, ∀q ∈ Q (4.70) u∈U y MEC unq t MIN ≤ CAP MEC q ∀n ∈ N, ∀q ∈ Q (4.71)
Then the maximum number of served users by the activated MEC-RFBs in

Chenghao WANG Contribution to robust network optimization type q is an optimal solution of the following problem:

max u∈U n∈N y MEC unq (4.72) s.t. n∈N y MEC unq ≤ 1 ∀u ∈ U, ∀q ∈ Q (4.73) u∈U y MEC unq t MIN ≤ CAP MEC q ∀n ∈ N, ∀q ∈ Q (4.74) y MEC unq ∈ {0, 1} ∀u ∈ U, n ∈ N, q ∈ Q, (4.75)
where constraint (4.73) is obviously implied by the objective function, is a so-called uniform 0-1 knapsack problem. Since all the weights of the items in the knapsack are identical, the problem is easy to solve and its optimal value is then defined by

CAP MEC q t MIN .
Therefore, the constraints (4.53)-(4.56), (4.62) and (4.65) can be replaced by the following constraints:

n∈N q∈Q x MEC nq CAP MEC q t MIN ≥ δ|U | (4.76)
The overall mathematical model denoted by 5G-RFB-RA-CP2 can be then rewritten as:

min n∈N k∈K q∈Q c k q x k nq (4.77) s.t. q∈Q x k nq ≤ 1 ∀n ∈ N, ∀k ∈ K (4.78) n∈N x k nq ≤ A qk ∀k ∈ K, ∀q ∈ Q (4.79) x RRH n 1 q + x RRH n 2 q ≤ 1 ∀q ∈ Q, ∀(n 1 , n 2 ) ∈ CONF q , n 1 = n 2 (4.80) y RRH unq t MIN ≤ x RRH nq CAP unq ∀u ∈ U, ∀n ∈ N, ∀q ∈ Q (4.81)
y RRH unq CAP unq ≤ q∈Q CAP RRH q x RRH nq ∀n ∈ N (4.86) n∈N x BBU nq ≥ n∈N x RRH nq ∀q ∈ Q (4.87) q∈Q x BBU nq ≤ q∈Q x RRH nq ∀n ∈ N (4.88) n∈N x MEC nq ≥ n∈N x RRH nq ∀q ∈ Q (4.89) q∈Q x MEC nq ≤ q∈Q x RRH nq ∀n ∈ N (4.90) n∈N q∈Q x MEC nq CAP MEC q t MIN ≥ δ|U | (4.91) x k nq ∈ {0, 1} ∀n ∈ N, ∀k ∈ K, ∀q ∈ Q (4.92) y RRH unq ∈ {0, 1} ∀u ∈ U, ∀n ∈ N, ∀q ∈ Q (4.93)

Deriving new valid inequalities

Considering the RRH-RFBs allocation, the explicit upper bound of maximum number of served users of an activated RRH-RFB in type of q is equal to U max q . The minimum number of activated RRH-RFBs is then defined by:

n∈N q∈Q x RRH nq U max q ≥ δ|U | (4.94)
Naturally, the maximum number of served users of an activated RRH-RFB in

Chenghao WANG Contribution to robust network optimization type q is actually limited by the minimal radio link capacity of eligible served users (where CAP unq ≥ t MIN ), hence UB q = CAP RRH q min{CAPunq: CAPunq≥t MIN ,∀u∈U } . A better upper bound is then U max q = min U max q , UB q , therefore, a stronger constraint and valid inequality are obtained compared to constraint 4.84 and (4.94) by writing:

u∈U q∈Q y RRH unq ≤ q∈Q x RRH nq min U max q UB q ∀n ∈ N (4.95) n∈N q∈Q x RRH nq min U max q , UB q ≥ δ|U | (4.96)
Similarly to MEC-RFB allocation, for a given valorization of x RRH nq , then x RRH nq = 1, ∀n ∈ N and y RRH unq for CAP unq ≥ t MIN imply:

u∈U y RRH unq ≥ 1 ∀n ∈ N, ∀q ∈ Q (4.97) n∈N y RRH unq ≤ 1 ∀u ∈ U, ∀q ∈ Q (4.98) u∈U y RRH unq CAP unq ≤ CAP RRH q ∀n ∈ N, ∀q ∈ Q (4.99)
Then, the maximum number of served users by the activated RRH-RFBs in type q (UB nq ) is an optimal solution of the problem below:

max u∈U n∈N y RRH unq (4.100) s.t. n∈N y RRH unq ≤ 1 ∀u ∈ U, ∀q ∈ Q (4.101) u∈U y RRH unq CAP unq ≤ CAP RRH q ∀n ∈ N, ∀q ∈ Q (4.102) y RRH unq ∈ {0, 1} ∀u ∈ U, ∀n ∈ N, ∀q ∈ Q (4.103)
which is a uniform knapsack problem [START_REF] Martello | Knapsack problems: algorithms and computer implementations[END_REF] that can be solved polynomially via sorting increasingly the selected users (CAP unq > t MIN ) by their size. So a tighter upper bound and valid inequality can be defined as: Thus, the overall mathematical model denoted by 5G-RFB-RA-CP3 can be rewritten as:

min n∈N k∈K q∈Q c k q x k nq (4.106) s.t. q∈Q x k nq ≤ 1 ∀n ∈ N, ∀k ∈ K (4.107) n∈N x k nq ≤ A qk ∀k ∈ K, ∀q ∈ Q (4.108) x RRH n 1 q + x RRH n 2 q ≤ 1 ∀q ∈ Q, ∀(n 1 , n 2 ) ∈ CONF q , n 1 = n 2 (4.109) y RRH unq t MIN ≤ x RRH nq CAP unq ∀u ∈ U, ∀n ∈ N, ∀q ∈ Q (4.
y RRH unq ≤ q∈Q x RRH nq min U max q , UB nq ∀n ∈ N (4.120) n∈N q∈Q x RRH nq min U max q , UB nq ≥ δ|U | (4.121) x k nq ∈ {0, 1} ∀n ∈ N, ∀k ∈ K, ∀q ∈ Q (4.122) y RRH unq ∈ {0, 1} ∀u ∈ U, ∀n ∈ N, ∀q ∈ Q (4.123)

A Benders decomposition Approach

Benders Decomposition is a major solution method used for optimization problems proposed by Benders, which has been intensively investigated over the five decades and used in many different application contexts, such as power system and network design ( [START_REF] Alguacil | Multiperiod optimal power flow using benders decomposition[END_REF], [START_REF] Binato | A new benders decomposition approach to solve power transmission network design problems[END_REF], [START_REF] Alysson | A survey on benders decomposition applied to fixed-charge network design problems[END_REF], [START_REF] Shahidehopour | Benders decomposition: applying benders decomposition to power systems[END_REF], [START_REF] Fortz | An improved benders decomposition applied to a multi-layer network design problem[END_REF], [START_REF] Rahmaniani | Accelerating the benders decomposition method: Application to stochastic network design problems[END_REF]), planning and scheduling ( [START_REF] Hooker | Planning and scheduling by logic-based benders decomposition[END_REF]), routing and scheduling ( [START_REF] Cordeau | Benders decomposition for simultaneous aircraft routing and crew scheduling[END_REF], [START_REF] Mercier | A computational study of benders decomposition for the integrated aircraft routing and crew scheduling problem[END_REF], [START_REF] Xin Cao | The integrated yard truck and yard crane scheduling problem: Benders' decomposition-based methods[END_REF]). For an exhaustive introduction to it, we refer the reader to the survey [START_REF] Rahmaniani | The benders decomposition algorithm: A literature review[END_REF]. Here, we proceed to recall some fundamentals.

Let us consider the following general LP:

min f T y + c T x (4.124) s.t. By + Dx = d (4.125) Ay = b (4.126)
x ≥ 0 (4.127) y ≥ 0 and integer, (4.128)

Chenghao WANG Contribution to robust network optimization where x ∈ R + is a real valued variable and y ∈ Z + is a non-negative integer complicating variable whose domain is defined by polyhedron Y : {Ay = b} with a known matrix A and a given vector b. A linking constraint By + Dx = d between

x and y must be satisfied with known matrix B and D for the associated x and y variables, and a given vector d. The objective function here is to minimize the total cost with the given cost vector f and c associated with x and y variables, respectively.

The Benders decomposition method partitions the problem in two: a master problem containing the y variables and a subproblem containing the x variables. With q(y) as the incumbent value for the x part, thus, we can define a LP using only variable y:

min f T y + q(y) (4.129) s.t. Ay = b (4.130)
y ≥ 0 and integer (4.131)

Then, we have the subproblem in terms of x. Note that if the subproblem is unbounded, the original problem is unbounded as well. If the problem is bounded, we calculate the value of q(y) by solving the following LP:

min c T x (4.132) s.t. Dx = d -By (4.133) x ≥ 0 (4.134)
Considering the dual variable π, associated with the subproblem, the dual form of q(y) can be rewritten as:

max π T (d -By) (4.135) s.t. D T π ≤ c (4.136) π ∈ R (4.137)
When the solution space is not empty, we can enumerate all extreme rays

Chenghao WANG Contribution to robust network optimization (denoted by a vector ρ) and extreme points (denoted by a vector π) of the feasible region D T π ≤ c. Given a solution vector y * , we can solve the dual problem by checking if we can find that:

1. a ρ i , i ∈ I such that ρ i (d -By * ) > 0, in which case the dual problem of q(y * ) is unbounded;

2. a π j , j ∈ J maximizing π j (d -By * ), in which case both the primal and dual of q(y * ) reach their optimality at this extreme point.

In the former case, there is a direction of boundlessness ρ i , this should be avoided due to the infeasibility of y in the master problem. A feasible cut ρ i (d -By * ) ≤ 0 is then added to the master problem to restrict the movement in this direction. In the latter case, the solution of q(y * ) (denoted by q) is one of the extreme points π j , j ∈ J. With the above notation of ρ an π, the master problem can be rewritten as:

min f T y + q (4.138) s.t. Ay = b (4.139) ρ i (d -By) > 0 ∀i ∈ I (4.140) π j (d -By) ≤ q ∀j ∈ J (4.141) q ≥ 0 ∀q ∈ R (4.142)
y ≥ 0 and integer (4.143)

Since there is an exponential number of extreme points and extreme rays, and because their enumeration is N P-hard, Benders decomposition starts without any cut and solves a relaxed master problem which gives an eligible candidate solution (y * , q * ). Taking this solution, it solves the subproblem to obtain an optimal value q(y * ). If q * is equal to q(y * ) then the candidate solution is optimal for the original problem. Otherwise, two cases may occur:

1. the dual is unbounded, then we select an extreme ray to generate a feasibility cut (i.e., (4.140));

Chenghao WANG Contribution to robust network optimization 2. the subproblem is solved optimally given a solution of y, then we select an extreme point to generate an optimality cut (i.e., (4.141)), which gives a lower bound of the incumbent q.

Master problem

Benders Decomposition

Dual subproblems information (solutions) feedbacks (feasibility cuts and/or optimality cuts) A general scheme of Benders decomposition is presented in Figure 4.4.7. The algorithm is based on an iterative constraint generation approach that runs until a convergence condition is met (e.g., q * = q(y * )), during which the objective value of the master problem provides an upper bound on the global solution while the value of the subproblem gives a lower bound at each iteration.

The key to generate a so-called Benders cut is to construct the inference dual (which is the problem of inferring a strongest possible bound from the constraint set, first proposed by Hooker and Yan) of the associated subproblems. In the case of BLPsubproblem, a logic-based Benders decomposition was proposed in [START_REF] John | Logic-based benders decomposition[END_REF] to obtain an inference dual β (e.g., obtained by a branch-andbound algorithm with a given solution y * ) such that f T y * +c T x ≥ β provides a tighter lower bound of the master problem. However, when the integer or boolean decision variables do not appear in the objective function of the original problem, the logic-based Benders decomposition could be useless because only feasibility cut will be generated, the inference dual cannot be obtained via the subproblems as no need of optimality cuts. The Combinatorial Benders decomposition has been proposed in [START_REF] Codato | Combinatorial benders' cuts for mixedinteger linear programming[END_REF] to tackle a set of conditional linear constraints (i.e., if x i = 1, then A i y = b i ) between decision variables x and y occurred in the original problem. Thanks to the minimal (or irreducible) infeasible subsystem of the associated constraints of variables y, if the involved subproblem is infeasible, then we can observe that at least one binary variable x i has to be changed to break the infeasibility, leading to the definition of a Combinatorial Benders' cut:

i:x * i =0 x i + i:x * i =1 (1 -x i ) ≥ 1.
Trying to apply these considerations to our Superfluid resource allocation problem, we can first note that our problem can be seen as made up of two main problems. One is the RFB allocation problem (solution space X), the other is the User-RRH-RFB assignment problem (solution space Y ) for a given solution of RRH-RFBs. Both the problems are BLP problems. Therefore, we propose a hybrid Benders decomposition to first project solution space Ω r (X, Y ) (where decision variable in Y is relaxed) into Π(X). For an optimal solution s * from Π(X) is feasible for Y, then we check whether the obtained decision values from Y are integral or not; otherwise we add a feasible cut into Π(X) found from above infeasibility. If the integrality requirement is met, then s * is also an optimal solution for the whole problem; otherwise, we add a Combinatorial Benders cut to Π(X), indicating at most one changeable decision variables that should be considered in the master problem.

More formally, we identify the RFB allocation problem as our master problem (denoted by P M ) in Benders decomposition and the User-RRH-RFB assignments problem as the slave problem (or subproblem) (denoted by P S (x RRH nq * )). The master problem can be written as follows:

min n∈N k∈K q∈Q c k q x k nq (4.144) s.t. q∈Q x k nq ≤ 1 ∀n ∈ N, ∀k ∈ K (4.145) n∈N x k nq ≤ A qk ∀k ∈ K, ∀q ∈ Q (4.146) x RRH n 1 q + x RRH n 2 q ≤ 1 ∀q ∈ Q, ∀(n 1 , n 2 ) ∈ CONF q , n 1 = n 2 (4.
CAP RRH q x RRH nq π * n + n∈N q∈Q min U max q , U nq x RRH nq η * n ≥ δ|U | (4.176)
With the notation introduced above, the proposed Benders decomposition can be described as is Algorithm 22. The main idea is to solve first the RFB allocation problems with the known explicit upper bound of maximum served users by an activated RRH-RFB. Using the solution of the RFB allocation problem, especially, we try to solve the linear relaxation variant of the involved problem. If the user coverage constraint is violated, then a Benders feasibility cut (4.176) is added into the master problem to improve the lower bound of total number of served users by a realized configuration of RRH-RFB allocation. Otherwise, we solve the associate subproblem with a given solution of RRH-RFB allocation. If the user coverage constraint is violated, then a Benders feasibility cut (4.162) is added to the master problem to make at least one change of RRH-RFB allocation to eliminate the incumbent infeasible solution. Otherwise, the optimal solution is obtained for the original problem.

Algorithm 5: Hybrid Benders Decomposition 4.5 Computational Results

5G test scenarios

In order to evaluate the performance of the proposed decomposition approach, we considered a 5G scenario including a hexagonal cellular geometry as in [START_REF] Thomas | Noncooperative cellular wireless with unlimited numbers of base station antennas[END_REF]. In this scenario, 9x9 candidate 5G nodes were considered, distributed according to a regular grid, and the users are quasi-uniformly located around the center of this network as shown in Fig. 4.5.8. In the test instances that we considered, the number of users varies from 25 to 100 with a step of 25, from 100 to 300 with a step of 50, and from 500 to 1000 with a step of 250. In line with this, 11 5G instances have been derived, in which we set up the minimum user downlink traffic t MIN varying from 10 to 50 Mbps with a step of 20, and the minimum required ratio of coverage α ∈ [0.1, 1.0] with a step of 0.1. Furthermore, concerning the computational settings of the machine Table 4.5.1 -Input parameters in line with these presented in [START_REF] Chiaraviglio | Optimal design of 5g superfluid networks: Problem formulation and solutions[END_REF] Parameters value Micro Macro

U MAX q 42 126 A RRH q 81 5 A BBU q 81 5 A MEC q 81 5 c RRH q 53951[e]
133951 that we used, we imposed a limit of 8Gb of maximum virtual memory and a time limit limitation of 900 seconds for experimentation (based on discussions with experts of the Superfluid architecture taking into account realistic real-world requirements). The value of other coefficients and parameters appearing in the problem are detailed in Table 4.5.1, where M denotes the mathematical model involved in [START_REF] Chiaraviglio | Optimal design of 5g superfluid networks: Problem formulation and solutions[END_REF], AM specifies the proposed model (5G-RFB-RA-CP1), and AMS is the simplified version (i.e., 5G-RFB-RA-CP3), respectively.

N opt calculates the number of optimal solutions, N inf eas counts the number of infeasible solutions, N timeout indicates the instances exceeding time limitation, and N memout specifies the number of instances that cannot be solved within 8G virtual memory. Moreover, T opt avg and T inf eas avg indicate the total averaged elapsed time on seconds for an optimal solution and the infeasible one.

The reference formulation proposed in [START_REF] Chiaraviglio | Optimal design of 5g superfluid networks: Problem formulation and solutions[END_REF] may present a performance that can be very costly in terms of computational time and RAM memory. The new formulation that we have proposed has instead the advantage of being able to identify an optimal solution for a much larger number of instances as it can be seen from the table of results. This improved performance can be attributed to the new valid inequalities that we have characterized and that allow to more effectively catch and express the Combinatorial relations that link the RFB). However, it should also be noted that some hard instances cannot still be solved within the time limit, thus encouraging us to refine and deepen the analysis of the new proposed approach as future work. Furthermore, we intend also to extend the tests to instances based on different 5G scenarios, thus obtaining the possibility of obtaining further insights about the behavior of the algorithm.

Green and Robust 5G Virtual Network Function Placement Problem

Introduction

The Fifth Generation of wireless telecommunications systems, widely known as 5G, has attracted a lot of attention in recent times, since it is largely considered as a crucial element for a full realization of a digital society and a critical technology to support the deployment of smart cities (see, for example, the work of the European Commission about this topic, e.g. [European-5G-Observatory]). 5G is going to offer enhanced service performances unknown to previous wireless technologies, such as data rates of at least 40 Mbps for tens of thousands of users, data rates of 100 Mbps for metropolitan areas, enhanced spectral efficiency and a dramatic reduction of latency (see e.g, [START_REF] Larsson | 5G Networks -Planning, Design and Optimization[END_REF][START_REF] Dahlman | 5G NR: The next generation wireless access technology[END_REF]).

5G will be strongly based on Network Function Virtualization (NFV), according to which network functions run on a set of Virtual Machines that are hosted in cheap commodity hardware servers [START_REF] Sherif Abdelwahab | Network function virtualization in 5g[END_REF]. This will sensibly reduce the cost of network infrastructures, decreasing the need for expensive dedicated hardware. For a very effective and accessible introduction to the main concepts, principles and features of network virtualization, we refer the reader to [START_REF] Mosharaf | Network virtualization: state of the art and research challenges[END_REF] and [START_REF] Schaffrath | Network virtualization architecture: Proposal and initial prototype[END_REF].

A central entity of network virtualization is represented by the Virtual Network (VN), which can be defined as a combination of network elements (network nodes and network links) realized over a Substrate Network (SN). Virtual nodes are inter-connected through virtual links, giving raise to a virtual topology. A determinant advantage associated with virtualization of elements like nodes and links of a SN is that a (high) number of distinct virtual network topologies characterized by very different performance characteristics may be defined using the same physical hardware as basis. Furthermore, another major advantage is represented by the fact that the characteristics and performance of these virtual topologies can be very easily changed by network operators, flexibly allocating or de-allocating resources from the physical hardware, allowing to vary them in very fast ways, giving a high degree of flexibility to adapt dynamically to changing requirements of users and traffic conditions.

A further critical advantage that is commonly attributed to Network Virtualization is represented by its potential of strongly supporting the concept of Infrastructure as a Service (IaaS) (see e.g., [START_REF] Bhardwaj | Cloud computing: A study of infrastructure as a service (iaas)[END_REF]), which is considered a very desirable property for next generation of internet architectures, in which the current figure of the Internet Service Providers (ISPs) should be split into two separate figures: 1) an Infrastructure Provider (InfraP) that has the task of constructing and maintaining the network equipment and 2) a Stochastic Programming (SP) that has the task of offering and managing end-to-end services. We note that this separation between the figure that maintains the infrastructure and the figure that offers services is spreading in many different engineering areas (e.g., public transportation and energy systems). As pointed out in works such as those cited above, the advent of network virtualization brings toward the identification of three distinct principal players that substitute the unique traditional service provider: a Virtual Network Provider (VNP) that has the task of arranging the virtual resources from a number of InfraPs, a Virtual Network Operator (VNO) that makes the virtual networks available according to the requirements of the SP, and an SP that can instead focus just on the offer of tailored virtual networks to the customers of the service.

The optimization of virtual network placement essentially requires to decide how to map a number of virtual network requests, corresponding to requested virtual topologies, to the available substrate network, while satisfying the demanded network resources with the available substrate network resources. In a more formal way, we can essentially and generally describe the placement problem as follows: we are given a substrate network SN(N, L), in which N is the set of substrate S N

In fr a P 1

In fr a P 2 su bs tr at e no de :

su bs tr at e lin k:

V N 1 vi rt ua l no de : vi rt ua l lin k:

V N 2 S P 2 vi rt ua l no de : vi rt ua l lin k: . So this actually requires to define a mapping of virtual nodes to substrate nodes and of virtual links to paths in the substrate network, while taking into account the constraints of satisfying the demand of resources by the virtual elements without exceeding the capacity of the substrate elements. Such problem of virtual network placement is known to be N P- Hard [Kolliopoulos and Stein, 1997]. Even remarkable subproblems of it, such as finding a virtual link mapping for a given node mapping, are known to be N P- Hard [Kolliopoulos and Stein, 1997].

Concerning the definition of a taxonomy of Network Virtualization problems in terms of their objectives, we can identify the following major classes of problems:

• maximization of the quality-of-service compliance of requests, in which the requests must be satisfied so that their features result as close as possible to the specifications fixed by the customers in terms of measures such as bandwidth, delay and CPU requirements (e.g., [START_REF] Inführ | Introducing the virtual network mapping problem with delay, routing and location constraints[END_REF]);

• maximization of the profit of the infrastructure provider, in which the provider naturally attempts at obtaining the highest economical return and must carefully consider how to accept and manage user requests over a time horizon, in order to maintain sufficient spare resources for dealing with (more profitable) unexpected requests and future unknown requests (see e.g, [START_REF] Chowdhury | Vineyard: Virtual network embedding algorithms with coordinated node and link mapping[END_REF]);

• maximization of the survivability of the user requests, which requires to setup specific backup resources in case of possible failures of elements of the substrate network; similarly to the previous class of problems, also in this case the provider must carefully choose how to reserve backup resources: indeed, the reserved resources cannot be used to satisfy other requests and thus limit the possibility of accepting future requests and may lead to sensible reduction in profit if not carefully dimensioned (see e.g, [START_REF] Shahriar | Joint backup capacity allocation and embedding for survivable virtual networks[END_REF][START_REF] Li | Survivable services oriented protection level-aware virtual network embedding[END_REF]);

• minimizing the total energy consumption, adopting a green network paradigm -this class of problems has emerged as one of the most important in computer networks in general, since the increase of such networks in size and complexity has lead to the consumption of huge amounts of energy. Such energy consumption not only represents a major cost for providers, but is also not sustainable from an environmental point of view (see e.g., [START_REF] Bilal | A taxonomy and survey on green data center networks[END_REF][START_REF] Rosario | Robust green wireless local area networks: A matheuristic approach[END_REF][START_REF] Wang ; Khalid | Towards green communication in 5g systems: Survey on beamforming concept[END_REF]);

As we started to discuss above, the problem of optimally designing virtual networks, allocating Virtual Network Function Component (VNFC)s to physical servers and managing the data flows between servers has received great attention

Chenghao WANG

Contribution to robust network optimization in recent times, in particular focusing on adopting a green networking perspective aiming at minimizing the overall power consumption (see e.g., [START_REF] Caggiani Luizelli | Piecing together the nfv provisioning puzzle: Efficient placement and chaining of virtual network functions[END_REF][START_REF] Juliver | Resource allocation in nfv: A comprehensive survey[END_REF][START_REF] Mechtri | A scalable algorithm for the placement of service function chains[END_REF], Marotta et al., 2017a[START_REF] Baumgartner | Towards robust network slice design under correlated demand uncertainties[END_REF]). However, while purely heuristic solution approaches for virtual network design have been quite widely investigated, the development of hybrid exact-heuristic algorithms exploiting the potentialities of mathematical programming (so-called matheuristic) has received limited attention, as it may also found in a recent survey such as [START_REF] Yang | Recent advances of resource allocation in network function virtualization[END_REF] and in other works that have focused on directly using state-of-the-art solvers to solve instances of small size or adopt simple ad-hoc heuristics (e.g., [START_REF] Zhang | Network service mapping and scheduling under uncertain processing time[END_REF]).

By our work presented here, we aim to start to fill this gap by proposing a new matheuristic for the green placement of Virtual Network Function in 5G, while taking into account the uncertainty of function requests which has been identified as a critical source of uncertainty (see e.g., [Marotta et al., 2017a[START_REF] Baumgartner | Towards robust network slice design under correlated demand uncertainties[END_REF]).

The remainder of this chapter is organized as follows: in Section 5.2, we describe a Binary Linear Programming (BLP) model for modeling the green and robust placement of VNFCs; in Section 5.3, we present a new matheuristic to fast solve the placement problem: finally, in Section 5.4, we report preliminary computational results and derive some conclusions.

A Binary Linear Programming model for VNFC Placement

From a modeling point of view, we can essentially describe the topology of the 5G network that we consider through a graph G(N, L), where N is the node set and L is the link set. Each link ∈ L corresponds to a pair (i, j), where i, j ∈ N are the nodes it connects. Each link is associated with a bandwidth b . The network interconnects a set of servers S and the node to which a server s is connected is denoted by n(s) ∈ N . Each server offers an amount of computational resources (e.g., CPU and RAM): denoting by R the set of resource types, the amount of resources available for each type r ∈ R at a server s ∈ S is denoted by a sr . The set of VNFCs is denoted by V and the set of service chains offered in the network ) and individual CPU, memory and disk capacities. In the example given in [Marotta et al., 2017a], server s i has installed a 1i CPU, a 2i RAM and a 3i disk. Each server is connected to a unique router (for example, s 1 is connected to n 1 ). Each link has a dedicated capacity and latency (for example, the latency for the link between n 1 and n 2 is denoted as l 12we omit bandwidth from Figure 5.2.2 to maintain readability). The servers, their capacities, together with the network nodes and links with their capacities form the

• variables f (v 1 ,v 2 ) ij ∈ {0, 1}, ∀(i, j) ∈ L, (v 1 , v 2 ) ∈ C∈C representing that link (i, j
) is used for data exchange between v 1 and v 2 belonging to some C ∈ C.

These variables are employed in the following Binary Linear Programming, denoted by BLP-VP, modeling the VNFC optimal placement problem:

(n,i)∈L b v 2 v 1 f v 1 ,v 2 ni - (i,n)∈L b v 2 v 1 f v 1 ,v 2 in = s∈S:n(s)=n b v 2 v 1 (x v 1 s -x v 2 s ) ∀n ∈ N, ∀(v 1 , v 2 ) ∈ C∈{C} C (5.6) (v 1 ,v 2 )∈ C∈{C} C b v 2 v 1 f v 1 ,v 2 ij ≤ B ij w ij ∀(i, j) ∈ L (5.7) w ij ≤ z i ∀(i, j) ∈ L (5.8) w ij ≤ z j ∀(i, j) ∈ L (5.9) f v 1 ,v 2 ij ≤ z i ∀(i, j) ∈ L (5.10) f v 1 ,v 2 ij ≤ z j ∀(i, j) ∈ L (5.11)
y s ∈ {0, 1} ∀s ∈ S (5.12)

x vs ∈ {0, 1} ∀s ∈ S, ∀v ∈ V (5.13)

w ij ∈ {0, 1} ∀(i, j) ∈ L (5.15) f v 1 ,v 2 ij ∈ {0, 1} ∀(i, j) ∈ L, ∀(v 1 , v 2 ) ∈ C∈{C} C (5.16)
The objective function (5.1) pursues the minimization of the total power consumption that is expressed as the sum of 1) a fixed power component, which represents the power consumed by activated servers, 2) a variable power component, which takes into account the amount of resources of servers that are used by the VNFCs assigned to the server, 3) the power consumed by the nodes of the substrate network, 4) the power consumed by the arcs of the substrate network.

Concerning the feasibility constraints, (5.2) impose that each VNFC must be allocated on exactly one server. Furthermore, constraints (5.3) are adopted to logically link the server activation and VNFC allocation decision variables: the activation of a server imposes that at least one VNFC is assigned to it. The constraints (5.4) are instead logical constraints imposing that, if the server is not activated, then no binary variable allocating a VNFC to the server can be activated. The constraints (5.5) model the resource capacity for each server and resource type, imposing that the summation of the amount of a specific resource consumed by VNFCs assigned to a server must not exceed the capacity of the server for that kind of resource. The constraints (5.6) are flow conservation constraints that regulate how links are used for transferring data between VNFCs, depending on which server the VNFCs are allocated to: the value of the right-hand-side depends upon the value of the difference of variables reflecting to which server the two involved VNFC of a chain are allocated. The set of constraints (5.7) is needed to represent the bandwidth capacity of each link, by imposing that it cannot be exceeded by the summation of bandwidth consumed by communications between VNFCs that are connected by means of link (i, j). The constraints (5.8) and (5.9) impose that using a link (i, j) requires to activate the two end nodes of the link. A similar role is played by the constraints (5.10) and (5.11), which impose to activate the end nodes of a link when traffic is sent over the link for supporting the communication between two VNFCs. Finally the constraints (5.12)-(5.16) define the feasible domain of all the decision variables involved in the problem.

Guaranteeing protection against resource uncertainty

Traffic routed through telecommunication networks is normally subject to uncertainty since the behavior of the users of the network is typically not exactly known in advance. Concerning this point, we refer the reader to [START_REF] Bauschert | Network planning under demand uncertainty with robust optimization[END_REF] for an accessible discussion about the motivations behind traffic uncertainty in telecommunications networks. Concerning the design of virtual network functions, the major uncertainty question is represented by the resources that are necessary to the functions generated by users and we reflect this in the robust optimization model that is derived here. Specifically, referring to the notation that we have introduced, we assume that the coefficients a vr representing the amount of resource r requested by a VNFC v is subject to uncertainty and its exact value is not known when the design problem is solved. In order to tackle such data uncertainty, we rely on the well-known Γ-Robustness model that was initially proposed by [START_REF] Bertsimas | The price of robustness[END_REF]. Since this is an interval uncertainty model, we must introduce a reference value āvr (nominal value) for the unknown value and a value ∆a vr representing the highest deviation that may occur from the nominal value. As a consequence, we assume that the actual value of the uncertain coefficient a vr lies in the symmetric interval: a vr ∈ [ā vr -∆a vr , āvr + ∆a vr ]

(5.17)

Once this deviation interval is defined, we proceed to show how it may be used as basis to derive robust counterpart of the uncertain resource constraints. We now proceed to focus attention on the constraints that we need to robustify, namely: v∈V āvr x vs ≤ a rs y s (5.18) that are the capacity constraints fixing the capacity of a server s ∈ S for each type of resource r ∈ R and containing the uncertain resource coefficients a vr . We note that in these constraints we have included the nominal values āvr .

As first step to derive the robust counterpart of these constraints, we can rewrite the constraints adding one term DEV rs (Γ, x) that indicates the worst deviation in value that the coefficients may attain for a given solution x and for Γ coefficients The worst deviation value DEV rs (Γ, x) can be computed as the optimal value of the following Combinatorial optimization problem in which we remark that the worst deviation of coefficient is represented by an increase in the resource that is requested by a VNFC v, since it tend to lead to a violation of the installed capacity:

max v∈V (∆a vr x vs ) φ rsv (5.20) s.t. v∈V φ rsv ≤ Γ (5.21) φ rsv ∈ {0, 1} v ∈ V (5.22)
In the problem above, binary variables φ rsv are introduced to represent whether a resource coefficient a vr is subject to its worst deviation ∆a vr and the unique cardinality constraint imposes that at most Γ coefficients may deviate. This value Γ may range from 0 (no protection against deviation and lowest price of robustness) to |V | (i.e., full protection against all VNFC deviating and highest price of robustness. By highlighting the fact that DEV rs (Γ, x) corresponds with an optimization problem, it becomes evident that (5.19) is actually non-linear. However, as elegantly proved by Bertsimas and Sim, it is possible to rely on duality theory to produce a linear robust constraint as follows. In primary, we may note that for given x, the value DEV rs (Γ, x) equals the optimal value of its linear relaxation, namely:

max v∈V (∆a vr x vs ) φ rsv (5.23) s.t. [δ rs ≥ 0] v∈V φ rsv ≤ Γ (5.24) [η rsv ≥ 0] 0 ≤ φ rsv ≤ 1 v ∈ V (5.25)
The dual variables δ rs , η rsv ∀v ∈ V are introduced for the corresponding constraints (5.24) and (5.25), respectively. The dual of DEV rs (Γ, x) is then formulated as follows:

w ij ∈ {0, 1} ∀(i, j) ∈ L (5.52) f v 1 ,v 2 ij ∈ {0, 1} ∀(i, j) ∈ L, ∀(v 1 , v 2 ) ∈ C∈{C} C
(5.53)

A New Matheuristic for Green Robust Virtual Network Function Placement

We present here a new matheuristic for optimal VNFC placement that is based on the integration of a Genetic Algorithm (GA) with an exact large neighborhood search, namely a search formulated as an optimization problem solved by a stateof-the-art solver such as CPLEX [IBM-ILOG-CPLEX]. The solver is also used for completing partial solutions of (ROB-BLP-VP) in an optimal way: for a fixed value configuration of a subset of decision variables, we employ the solver to find a feasible valorization of all the remaining variables while optimizing the objective function. At the basis of this matheuristic there is the consideration that, while a state-of-the-art solver may find difficulties in identifying good quality solutions for ROB-BLP-VP, it is instead able to efficiently identify good quality solutions for appropriate subproblems of ROB-BLP-VP, obtained by fixing the value of a consistent number of decision variables.

GAs are widely known meta-heuristics that draw inspiration from the evolution of a population (see [START_REF] David | Genetic algorithms and machine learning[END_REF] for an exhaustive introduction to the topic). The individuals of the population represent solutions of the optimization problem and the chromosome of an individual corresponds to a valorization of decision variables of a solution. The quality of an individual/solution is assessed through a fitness function. The GA begins with the definition of an initial population that then changes through evolutionary mechanisms like crossover and mutation of individuals, until some stopping criterion is met.

The general structure of the GA that we take as reference and adapt that is presented in Algorithm 6. We now proceed to detail how the elements and the phases that have been presented above for the generic hybrid genetic algorithm are adapted to be applied to the problem ROB-BLP-VP.

obtain the optimal solution of (ROB-BLP-VP) for a fixed (ȳ, x). We denote the set of individuals constituting the population at a generic iteration of the algorithm by POP.

Evolution of the population

Selection. The individuals chosen for being combined and generating the new individuals are chosen according to a tournament selection principle: we first create a number β of (small cardinality) groups of individuals by randomly selecting them from POP. Then the γ individuals in each group associated with the best fitness value are combined through crossover.

Crossover. We form the couples that generate the offspring according to the following procedure. From the tournament selection, we obtain βγ individuals that are randomly paired in couples, each generating one offspring. Assuming that the two parents are associated with chromosomes/partial solutions (y 1 , x 1 ) and (y 2 , x 2 ), the chromosome of the offspring (y off , x off ) is defined according to two rules:

1. if the parents have the same binary value in a position j, then the offspring inherits such value in its position j (i.e., if (y 1 , x 1 ) j = (y 2 , x 2 ) j then (y off , x off ) j = (y 1 , x 1 ) j );

2. if the parents have distinct binary values in a position j, then the offspring inherits a null value (i.e., if (y 1 , x 1 ) j = (y 2 , x 2 ) j then (y off , x off ) j = 0).

Possible violations in the constraints (5.2) and (5.5) associated with (y off , x off ) are then repaired. The main rationale at the basis of this procedure is assuming that two solutions having the same valorization of a variable is a good indication that such valorization should be maintained also in the offspring.

Exact Improvement Search

We attempt at improving the best solution found by the GA through an exact large neighborhood search, namely a search that is formulated as a suitable Binary Linear Programming problem solved by a state-of-the-art optimization solver

Chenghao WANG

Contribution to robust network optimization

Computational results

We assessed the performance of the proposed matheuristic by considering 20 instances that refer to a network made up of 10 nodes to which 50 servers are connected and that are defined for different VNFC features, defined referring to the works [Marotta et al., 2017a,b]. To execute the tests, we employed a Windows machine with 2.70 GHz professor and 8 GB of RAM. As optimization solver, we relied on IBM ILOG CPLEX 12.5, which is interfaced through Concert Technology with a C/C++ code. The global time limit imposed to CPLEX to solve (ROB-BLP-VP) is set to 3600 seconds. The same time limit is set for the matheuristic (denoted here by MatHeu), assigning 3000 seconds to the GA phase and 600 to the improvement phase based on RINS (in which we set = 0.1). The initial population includes 100 individuals/solutions and, at each iteration, we consider β = 10 groups from each of which γ = 2 individuals are chosen.

The results of the computational tests are presented in Table 1, where: ID identifies the instance; T * (CPLEX) and T * (MatHeu) are the time (in seconds) that CPLEX and MatHeu needs to find the best solution within the time limit, respectively, whereas ∆T * % is the percentage reduction in time that MatHeu grants to find a solution that is at least as good as the best solution found by CPLEX. Finally, ∆P * % is the reduction in power consumption that the best solution found by MatHeu grants with respect to the best solution found by MatHeu within the time limit.

We have paid particular attention to the computational time aspect, since, according to discussions that we had with professionals of the sector, identifying high quality solutions within limited amount of time is considered a particularly important objective when establishing sets of virtual networks in a business context.

Concerning the difficulty of solving (ROB-BLP-VP), as highlighted in several works such as [START_REF] Caggiani Luizelli | Piecing together the nfv provisioning puzzle: Efficient placement and chaining of virtual network functions[END_REF], Marotta et al., 2017a], even simplified deterministic versions of (ROB-BLP-VP) may prove difficult to solve for state-of-the-art optimization solvers also in the case of instances. We confirm such behavior in the case of our instances, which highlights the need for fast (heuristic) solution algorithms. On the basis of the results, we can say that M atHeu, for all the instances, is able to return a solution that is at least as good as the best solution found by CPLEX within the time limit in 24% less time, on average. Concerning the reduction in consumed power, we can instead notice that M atHeu allows to find better quality solution than CPLEX within the time limit, with a reduction in power consumption that can reach 24% and on average is equal to about 13%. The better performance of the matheuristic results particularly evident for the second half of instances.

These results, which have been presented in our publication [START_REF] Bauschert | A matheuristic for green and robust 5g virtual network function placement[END_REF], have resulted remarkable and, as future work, have encouraged us to attempt at refining the solution construction mechanism, better exploiting the specific features of the mathematical model of (ROB-BLP-VP) to define the rules adopted to generate the initial population and to generate the offspring solutions by crossover.

As an alternative to the exact search based on RINS, we have also evaluated the possibility of adopting another exact neighborhood search based on using hamming distance constraints. Specifically, given a feasible solution associated with a fixing (ȳ, x) of the binary variables x and y, we considered the neighborhood of feasible solutions that can be obtained by modifying at most a number D of valorization of binary variables. Such neighborhood can be formally defined by adding the following constraint: j:(ȳ,x) j =0 (y, x) j + j:(ȳ,x) j =1

(1 -(y, x) j ) ≤ D, which counts the number of binary variables that have switched their value from 0 to 1 and from 0 to 1, imposing that such number must not exceed the value D. Similarly to the RINS approach, the resulting neighborhood is explored by means of a state-of-the-art solver. 16.8 9.9

In Table 5.4.2, we report comparisons of the performance granted by the two alternative exact searches, in terms of power reduction that the best solution found within the time limit by each search is able to grant with respect to the best solution found by CPLEX within the time limit. Specifically, for each instance, we report the percentage power consumption reduction P * (MatHeu-RINS) granted by the GA with RINS with respect to CPLEX, whereas P * (MatHeu-HD) reports the percentage power consumption reduction P * (MatHeu-HD) granted by the GA with the hamming distance constraint set to D = 2 with respect to CPLEX.

Evaluating the results, the adoption of RINS guarantees a higher reduction in power for all but four instances, indicating that neighborhood search defined with D = 2, while reducing the consumption with respect to the best solution found by CPLEX for instances, results not competitive with respect to RINS. However, the four cases in which the hamming distance-based constraint neighborhood performs better may suggest that tuning the value D and refining the definition of the neighborhood, involving a different set of variables, could possibly allow to define an alternative effective exact search. We consider this a possible subject of future research.

Conclusions and Future Work

Infrastructures based on networks nowadays constitute a fundamental component of our everyday life and continue to grow in size and complexity. Designing and managing the networks at their basis have become a very complex task and the adoption of mathematical optimization approaches has clearly shown to grant big advantages for identifying high quality design and management solutions. This Ph.D. Thesis has been focused on proposing new optimization modeling and algorithmic approaches for dealing with real-world network optimization problems arising in the transportation and telecommunications field. Since the focus has been on real-world applications, a relevant aspect that we have taken into account has been represented by data uncertainty, i.e. the fact that the value of a subset of input data of the problem is not exactly known when the problem is solved. In order to deal with such data uncertainty, we have also investigated the development of new modeling and algorithmic robust optimization approaches, which aim at identifying solutions that maintain their feasibility and optimality even when input data are subject to deviations in value.

More precisely, in the context of transportation problems, we have considered the flight level assignment problem, which arises in air traffic management and consists of establishing the flight levels of a set of aircraft in order to improve the total avenue from assignment, reduce the total number of potential en-route Air Traffic Flow Management (ATFM) conflicts and also the total ATFM delay. In this context, we proposed a new chance-constrained optimization problem and iterative solution heuristic which is based on both analytical and sampling methods. Besides transportation problems, this Thesis has also focused on the optimal design of 5 th generation of wireless networks considering Superfluid and virtual architectures. Specifically, the 5G Superfluid architecture is based on atomic virtual entities called Reusable Functional Block (RFB)s and we investigated the prob-

Chenghao WANG

Contribution to robust network optimization lem of minimizing the total installation costs of a Superfluid network composed of virtual entities and realized over a physical network, while guaranteeing constraint on user coverage, downlink traffic performance and technical constraints on RFBs of different nature. To solve this hard problem, we proposed a Benders decomposition approach. Concerning instead the design of general virtual networks, we adopted a green paradigm that pursues energy-efficiency and tackled a state-ofthe-art robust mixed integer linear programming formulation of the problem, by means of a new matheuristic based on combining a genetic algorithm with exact large neighborhood searches.

Results of computational tests executed considering realistic problem instances have shown the validity of all the new optimization modeling and algorithmic approaches proposed in this Thesis for the transportation and telecommunications problems sketched above.

As ongoing work, we are studying the generation of the application used to solve the CCP FLA problem. Also, we are extending the computational tests to larger set of instances and attempting at strengthening the performance of the algorithm, investigating tuning strategies for setting the parameters at the basis of the various solution approaches. The aim is to include such new results in the journal versions of the conference papers that have been published during the Ph.D. Furthermore, as future work, we intend to investigate the identification of other class of valid inequalities in the context of 5G SF networks, better catching resource interactions between RFBs of different nature, with the aim of improving the convergence of our Benders decomposition solution approach. Moreover, concerning the matheuristic, we believe that a significant improvement in performance could be obtained by better integrating tight formulations in the solution process: stronger formulations could be used as basis for defining an initial population characterized by individuals with higher fitness, exploiting the valuable information coming from stronger linear relaxations, thus providing stronger individuals from the first iterations. (6.6) where v i , v j are the velocity of aircraft i and j, respectively, ρ = v j /v i , λ = sin θ/ ρ 2 -2ρ cos θ + 1.

If sin 2 θ = 0, then cos θ = 1 under the "Semicircular/hemispheric" rule, and the trailing conflict may occur for t ∈ [max(t o i , t o j ), min(t o i + d s /v i , t o j + d s /v j )] where d s specifies the distance between two waypoints when a trailing conflict is encountered. We then have: Therefore, the minimum distance for two aircraft of crossing conflict and trailing conflict is formulated: To cope with the converging conflict and diverging conflict, we put together the four types of conflict in a unified scenario in Figure 6.1.2. Then, the minimum

d t ij 2 = v 2 i (1 -ρ) 2 t

3. 6 . 1

 61 Coefficients of GMM for flight departure time delays . . . . . . . . 3.6.2 Corresponding observed p value when α = 0.05, N mc = 10000 . . . . 3.8.3 Characteristics of test instances . . . . . . . . . . . . . . . . . . . . 3.8.4 Posterior feasibility probability for all flights assigned at their most preferred flight level . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8.5 Average en-route ATFM delay per flight and induced cost per flight when all ω ij takes their worst-case value ω ij and all flights are assigned at their most preferred level . . . . . . . . . . . . . . . . . . 3.8.6 Computation results for robust FLA using the Soyster method . . . 3.8.7 Computational result for robust FLA by different feasibility estimation methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8.8 Computational result for Robust FLA by different feasibility estimation methods with different configurations of W i . . . . . . . . . 3.8.9 Computational result for Robust FLA by different estimation methods with different configurations of when W i = 5 minutes . . . . .

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5.2 Comparison on solved instances for each method . . . . . . . . . . . . . 5.4.1 Computational results . . . . . . . . . . . . . . . . . . . . . . . . . 5.4.2 Comparison of alternative exact search -power reduction . . . . . . List of Figures 3.4.1 The mission of a flight . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.2 Air transport delay . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.3 En-route potential conflict between two aircraft cruising at same level 3.5.4 FLA instance construction from 3-GC . . . . . . . . . . . . . . . . . 3.6.5 Estimation of feasibility probability of obtained solution by the Hoeffding's Inequality . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6.6 Estimation of feasibility probability of obtained solution by Monte-Carlo simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6.7 An example illustrating the heuristic estimation of probability of the expression ω 1 + ω 2 + ω 3 + ω 4 . . . . . . . . . . . . . . . . . . . 3.6.8 Estimation of feasibility probability of obtained solution by heuristic estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7.9 General approximation framework tackling the Robust FLA problem 3.8.10Lower bound of posterior solution feasibility probability for the robust FLA by different estimation methods with different configurations of when W i = 5 minutes . . . . . . . . . . . . . . . . . . . . 4.1.1 Evolution of wireless network from 1G to 5G, source from [Gohar and Nencioni, 2021] . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.2 Class diagram for the ETSI VNF (left side) and the 5G Superfluid RFB (right side), source from [Bianchi et al., 2016] . . . . . . . . . 4.2.3 5G SF Architecture, source from [Bianchi et al., 2016] . . . . . . . . 4.2.4 A complete RFB chain and the exchanged information among them 4.3.5 Problem structure solved in [Chiaraviglio et al., 2018] . . . . . . . . 4.3.6 Problem structure solved in the new proposed model . . . . . . . . 4.4.7 General schema of Benders decomposition method . . . . . . . . . . Chenghao WANG Contribution to robust network optimization VIII 4.5.8 A 5G network instance with 81 candidate nodes and 500 users . . . . . . 5.1.1 Mapping of Virtual Networks to a Substrate Network . . . . . . . . 5.2.2 an example of the VNF placement and mapping problem in [Marotta et al., 2017a] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1.1 En-route potential conflict between two aircraft cruising at same level124 6.1.2 A unified en-route potential conflict between two aircraft cruising at same level . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1.3 GMM distributions associated with flight departure delays . . . . . 6.1.4 truncated GMM distribution of ω ij . . . . . . . . . . . . . . . . . .

  Chenghao WANGContribution to robust network optimization 1.2 -Contributions of this Ph.D. Thesis

1 .

 1 Akli Fundo, Dritan Nace, and Chenghao Wang. A heuristic approach for the robust flight level assignment problem. In International Conference on Belief Functions, pages 86-94. Springer, 2018 2. Thomas Bauschert, Fabio D'andreagiovanni, Andreas Kassler, and Chenghao Wang. A matheuristic for green and robust 5g virtual network function placement. In International Conference on the Applications of Evolutionary Computation (Part of EvoStar), pages 430-438. Springer, 2019 3. Chenghao Wang, Fabio D'Andreagiovanni, and Dritan Nace. Solving a resource allocation problem in rfb-based 5g wireless networks. In Third International Balkan Conference on Communications and Networking (BalkanCom 2019), 2019 4. Akli Fundo, Jean-Benoist Leger, Dritan Nace, and Chenghao Wang. Dealing with uncertainty in atm-the flight level assignment problem. In 21e congrès annuel de la Société Française de Recherche Opérationnelle et d'Aide à la Décision (ROADEF 2020), 2020 5. Liyang Xiao, Zhengpei Wang, Zheyi Tan, and Chenghao Wang. A solution method for the maritime pilot scheduling problem with working hour regulations. Asia-Pacific Journal of Operational Research, 38(03):2040015, 2021

  studied the deterministic and stochastic TFMRP, where flight cancellation and re-routing, and arrival and departure capacity at the airport are taken into account in order to reduce the ground holding cost and air delays imposed on flights. Furthermore, the stochastic counterpart is reformulated by a deterministic equivalent model (a medium-scale mixed 0-1 model) considering the uncertainties of departure and arrival capacity of the airports, the capacity of sectors, and flight demands. More recently, Chen et al. [2017] proposed a polynomial approximation-based chance-constrained optimization method to address the uncertainty of capacity of sectors for the TFMRP. Sandamali et al. [2017] introduced a flight routing and scheduling model taking into account the uncertainty due to aircraft departures. Flight Level Allocation or Flight Level Assignment (FLA) is another important approach to reduce the average fuel consumption per flight and to improve total Chenghao WANG Contribution to robust network optimization

  applied a genetic technique to the tactical Flight Level Assignment. Vela et al. [2009] proposed a compact formulation of a complete optimization model for speed control and flight level assignment to reduce fuel burn over time horizons between 15-45 minutes. The combining of flight level allocation and ground holding Barnier and Allignol [2011] was investigated to reduce the 4 Dimensions (4D) trajectorybased conflict. Allignol et al. [2012] proposed a flight level allocation schema to avoid the conflicts occurring during the cruise phase of intersecting flights. Soler-Arnedo et al. [2013] studied the contrail sensitive 4D trajectory planning problem performing the permitted step climbs to change flight level in order to minimize the fuel consumption and Carbon dioxide (CO2) emissions. Moreover, a two-step hybrid metaheuristic is proposed to solve the flight level allocation problem in order to avoid most losses of separation occurring between cruising flights before running the automated conflict resolution. More recently, Gimenez-Guzman et al. [2020] study the joint optimization of fuel consumption and Flight Level Assignment using graph coloring. However, most of them reformulated the flight level allocation/assignment problem considering only the crossing conflict (Figure 3.4.3a) at the same level (whereas trailing conflict (Figure 3.4.3b), converging conflict (Figure 3.4.3c) and hybrid conflict (Figure 3.4.3d) may encounter), without consideration of data uncertainties. To tackle the uncertainty of conflict around crossing waypoints,[START_REF] Constans | Optimal flight level assignment: introducing uncertainty[END_REF] studied the optimal Flight Level Assignment taking account of uncertainty determination of flight crossing time at a given point.[START_REF] Klopfenstein | The robust flight level assignment problem[END_REF] introduced a mathematical model for the robust FLA via a chance-constrained optimization approach and proposed a fast approximation framework to solve the robust FLA efficiently, assuming the induced cost due to resolution of potential conflict is uncertain and bounded.[START_REF] Fundo | A heuristic approach for the robust flight level assignment problem[END_REF] have investigated the robust FLA problem under the angle of uncertainty due to flight departure time delays. Specifically, the statistical delay model3.3 -Contribution and organization

Flight

  phases A flight is subdivided into different phases as illustrated in Figure 3.4.1. The conflicts concerned in this work occur during the "Cruise Phase" for flights flying at same level. An induced en-route ATFM delay may occur due to the flight departure delay as illustrated in Figure 3.4.2. A set of feasible cruise altitudes exist depending on the aircraft type of which the operating costs are minimized at the optimum cruise altitude.

  Figure 3.4.1 -The mission of a flight

  Figure 3.4.3 -En-route potential conflict between two aircraft cruising at same level

F:

  The set of flights, indexed by i, j. L: The set of flight levels, indexed by l. F l : The subset of flights allowed flying at the flight level l ∈ L. L i : The subset of eligible flight levels for each flight i ∈ F, |L i | = 3. For instance for an aircraft of type Airbus A320, the three preferred level flights are FL390, FL410 and FL370. Chenghao WANG Contribution to robust network optimization S l i : The subset of flights having a potential conflict with flight i at level l. Parameters:W i : The maximal acceptable en-route ATFM delay for a given flight i.

5Fix

  Let x * be the optimal solution found; 6 Estimate the concerned feasibility probability of obtained solution x * ; 7 if Constraints (3.22) for all concerned flights in F l \ I c are respected then 8 i as corresponding to the highest violation of constraints (3.22); 12 I c ← I c ∪ {i}; 13 end 14 while Unsolved ; 15 return Last obtained solution x * .

  Figure 3.6.5 -Estimation of feasibility probability of obtained solution by the Hoeffding's Inequality

Figure 3

 3 Figure 3.6.6 -Estimation of feasibility probability of obtained solution by Monte-Carlo simulation

  Chenghao WANGContribution to robust network optimization Algorithm 3: A modified-EM algorithm to approximate the distribution ω p into a Gaussian Mixture Model of K components Input: The distribution ω p Output: A Gaussian Mixture Model of K components 1 Initialize the coefficients µ k , σ k , c k for the GMM of K components; 2 do 3 Calculate the conditional probability of sampling points (X i ) under each component k:

Figure 3

 3 Figure 3.6.7 -An example illustrating the heuristic estimation of probability of the expression ω 1 + ω 2 + ω 3 + ω 4

Figure 3 Figure 3

 33 Figure 3.6.8 -Estimation of feasibility probability of obtained solution by heuristic estimation

  3 and FR 4 by randomly generating 15%, 30%, and 50% supplementary traffics, respectively. Each instance is characterized by the number of flights (|F|), the number of flights having at least one conflict with others at their most preferred level (|F c |), the maximum number of potential conflicts per flight at their most preferred level (|S i | max ), total number of potential conflicts at their most preferred level ( i∈F |S i |).

b 0 i

 0 is calculated as: b 0 i = 180 * 0.681 * MPM i 4 . The level changing cost for a flight from the most preferred level to the other feasible alternative levels is considered as 10% of b 0 i . For example, a flight i flies from Bastia-Poretta to Lyon-Saint-Exupéry at level FL180 (the most preferred level, the two feasible alternative levels are FL160 and FL200), the corresponding flight miles is 588 Kilometers, it has a cumulative en-route ATFM delay of 140 (45 and 15, respectively) minutes with other flights assigned at the same level FL290 (FL270 and FL310, respectively). If we assign this flight i at level FL290, then we have the assignment profit: b 0 i = 180 * 0.681 * 588 = 72077.04 euros, the en-route ATFM delay cost: y i * α i = 30 * 70 + 30 * 140 + 30 * 280 + 30 * 560 + 20 * 1120 = 53900 euros, hereby the revenue of this assigned flight i at level FL180 is 72077.04 -53900 = 18177.04 euros. The revenue is reduced heavily due to the en-route ATFM delay. If this flight is canceled, then we considered a penalization as 72077.04 * 10 = 720770.4 euros. Furthermore, if we assign this flight at level FL160 (FL200, respectively), we have a revenue 72077.04 * 0.9 -30 * 70 -15 * 140 = 67877.04 (72077.04 * 0.9 -15 * 70 = 71027.04, respectively) euros. Table3.8.4 -Posterior feasibility probability for all flights assigned at their most preferred

  8.5 specifies the average en-route ATFM delay per flight (ỹ i ) and the corresponding induce cost per flight (the sum of the level changing cost, the enroute ATFM delay cost, and the cancellation cost per flight) when a Soyster model (all ω ij takes their worst-case value ω ij ) is used to deal with the flight departure delay uncertainty for all flights assigned at their most preferred flight level, which means there is neither level changing their flight level nor flight cancellation. For these solutions of all different instances, there are in average 18852.02 euros of total induced cost per flight (equal to delay cost per flight as there is neither flight changing their flight level nor any flight canceled) due to en-route ATFM delay, where the average en-route ATFM delay in such solution is much far from the reference value (0.5 minute per flight).

  flight en-route ATFM delay, and the flight cancellation) when the upper bound of cumulative en-route ATFM delay (y i ) per flight (W i ) is larger, since increasing W i hwill result in possible longer ATFM delays. In the other hand, for the configuration of larger value of W i (e.g., 10 minutes), they requires more computation time because the associated y i is incremented step by step at end of each iteration when the corresponding constraint (3.22) is violated. For each configuration of W i ,Chenghao WANGContribution to robust network optimization

  Figure 3.8.10 -Lower bound of posterior solution feasibility probability for the robust FLA by different estimation methods with different configurations of when W i = 5 minutes

1 1

 1 exabyte = 10 6 Gigabytes Chenghao WANG Contribution to robust network optimization

Figure 4 . 1 . 1 -

 411 Figure 4.1.1 -Evolution of wireless network from 1G to 5G, source from[START_REF] Gohar | The role of 5g technologies in a smart city: The case for intelligent transportation system[END_REF]

  Figure 4.1.2 -Class diagram for the ETSI VNF (left side) and the 5G Superfluid RFB (right side), source from[START_REF] Bianchi | Superfluidity: a flexible functional architecture for 5g networks[END_REF]

  Figure 4.1.2). A RFB Description and Composition Language (RDCL) is introduced for characterizing and describing each RFB in a formal manner at the platform-agnostic node-level and network-level. Moreover, a RFB Execution Environments (REE) is specified to support the execution and deployment of the RDCL scripts and the relevant coordination of the signal/radio/packet/flow/network processing primitives as shown in Figure 4.2.3. The figure visualizes that the Reusable Functional Block (RFB)s are orchestrated recursively at each level (i.e., network level and node level). At each level, two main actors are identified: REE User and REE Manager. The REE User requests the deployment/execution of a service/service component described using a RDCL script to the REE manager, and the REE Manager is in charge of deploying/executing the RDCL script using the resources within its REE. Two main interfaces (denoted by UM API and MR API) are characterized to support the interaction between REE User and REE Manager, so that the REE User can deploy a service or a component into an REE, and between REE Manager and REE Resource Entity, so that the REE Manager can interact with the resources in its REE.

Figure 4

 4 Figure 4.2.3 -5G SF Architecture, source from[START_REF] Bianchi | Superfluidity: a flexible functional architecture for 5g networks[END_REF]

  13) Chenghao WANG Contribution to robust network optimization φ un ≤ q∈Q CAP unq y RRH nq ∀u ∈ U, ∀n ∈ N (4.14) φ un ≤ CAP max u x un ∀u ∈ U, ∀n ∈ N (4.15)

  (4.10) and (4.11) indicate that a BBU-RFB of type q ∈ Q installed in node n 1 ∈ N can be part of the chain serving the RRH-RFB of the same type installed in node n 2 ∈ N . Similarly, constraints (4.12) and (4.13) indicate that a MEC-RFB of type q ∈ Q installed in node n 1 ∈ N can be part of the chain serving the RRH-RFB of the same type installed in node n 2 ∈ N . Constraint (4.14)-(4.18) specify that for each user, the amount of downlink traffic is then limited by the maximum radio link capacity on node n ∈ N .

Figure 4

 4 Figure 4.3.5 -Problem structure solved in[START_REF] Chiaraviglio | Optimal design of 5g superfluid networks: Problem formulation and solutions[END_REF]

  Figure 4.3.6 -Problem structure solved in the new proposed model

,

  UB nq ≥ δ|U | (4.105)

Figure 4

 4 Figure 4.4.7 -General schema of Benders decomposition method

Figure 4

 4 Figure 4.5.8 -A 5G network instance with 81 candidate nodes and 500 users

Figure 5 . 1 . 1 -

 511 Figure 5.1.1 -Mapping of Virtual Networks to a Substrate Network

  Chenghao WANGContribution to robust network optimization is denoted by C (a service chain corresponds with a subset of VNFC that must be executed to provide a requested service to a user). When executed, a VNFC v ∈ V requires an amount a vr of each resource type r ∈ R. Each chain C ∈ C corresponds to a subset of pairs (v 1 , v 2 ) belonging to V × V . The exchange of data between v 1 and v 2 in a pair (v 1 , v 2 ) requires an amount of bandwidth b v 2 v 1 in each traversed network link. Concerning power consumption, every node n ∈ N and link ∈ L consumes P n and P when used, respectively. Each server s ∈ S has a consumption that is a linear function in the range [P min s

Figure 5

 5 Figure 5.2.2 -an example of the VNF placement and mapping problem in[Marotta et al., 2017a]

  Chenghao WANGContribution to robust network optimization allowed to differ from their nominal values:v∈V āvr x vs + DEV rs (Γ, x) ≤ a rs y s ∀s ∈ S, ∀r ∈ R (5.19)

Figure 6 . 1 . 1 -

 611 Figure 6.1.1 -En-route potential conflict between two aircraft cruising at same level

  o i -t o j |, crossing conflict: sin θ > 0 (v j |t o i -t o j | -d s |1 -ρ|)/ max(ρ,[START_REF]European-5G-Observatory. European 5g observatory[END_REF], trailing conflict: sin θ = 0 (6.10)

  Figure 6.1.3 -GMM distributions associated with flight departure delays

Figure 6

 6 Figure 6.1.4 -truncated GMM distribution of ω ij

Table 3

 3

	Component	C 1	C 2	C 3	C 4
	c i	0.37	0.40	0.15	0.08
	µ i	-17.15	-7.31	19.57	69.13
	σ 2 i	88.20	89.33	1007.73	3926.00
	ω

.6.1 -Coefficients of GMM for flight departure time delays ij parameters to take random values in given interval. So, from now ω ij express the random parameter varying in [0, ω ij ], where ω ij = t msd ij

Table 3 .

 3 6.2 -Corresponding observed p value when α = 0.05, N mc = 10000 Calculte induced delay ω ij for each pair of flights having potential conflict Calculate the probability of j∈S l i ω ij x j ≤ y i for all assigned flights i such that i ∈ F l \ I c Let feasibility probability of x * be the minimum obtained probability

		0.02	0.03	0.05	0.10	0.15	0.20	0.25
	p(%)	98.22	97.27	95.35	90.48	85.58	80.65	75.71
				Estimate robustness-Feasibility Probability
					Sampling Method
			)
			, y F l
			(I c , F l
		Start	= Solution of RP l					End
			x *

Generate randomly N mc (e.g., 10000) scenarios for each assigned flight i such that i ∈ F l

Table 3

 3

			.8.3 -Characteristics of test instances
	Instance			Characteristics
		|F|	|F c |	|S i | max	i∈F |S i |
	FR_1	1641	687	9	1612
	FR_2	1887	810	9	1888
	FR_3	2133	931	9	2090
	FR_4	2461	1150	10	2782

Table 3 .

 3 8.3 presents the characteristics of test instances corresponding to dailyFrench air traffic with 134 airports and 715 waypoints. The instance FR 1 corresponds to French air traffic of August 12th, 1999. To accommodate current air traffic, we generate the instances FR 2 , FR

Table 3 .

 3 8.5 -Average en-route ATFM delay per flight and induced cost per flight when all ω ij takes their worst-case value ω ij and all flights are assigned at their most preferred level

	Instance	FR_1 FR_2 FR_3 FR_4 Avg.
	ỹi [minutes]	2.42	2.45	2.58	2.95	2.63
	Induced cost per flight [e]	15079.34 13934.76 23282.89 21221.09 18852.02

Table 3 .

 3

Table 3 .

 3 8.6 -Computation results for robust FLA using the Soyster method

		W i /Instance	FR_1 FR_2 FR_3 FR_4 Avg.
		GainRev%	20.85% 18.45% 46.33% 33.58% 29.80%
	i =5 W	CHA% CAN% ỹi [minutes]	20.66% 21.63% 22.32% 23.81% 22.28% 0.24% 0.21% 0.28% 0.49% 0.32% 0.25 0.23 0.19 0.23 0.23
		Delay cost per flight [e]	17.72	16.27	13.42	16.03	15.76
		Induced cost per flight [e] 4852.70 4730.81 5191.72 7102.88 5596.29
	i =10 W	GainRev% CHA% CAN% ỹi [minutes]	24.29% 21.38% 51.43% 41.36% 34.62% 20.54% 21.52% 22.36% 23.97% 22.29% 0.06% 0.05% 0.05% 0.12% 0.07% 0.29 0.27 0.23 0.29 0.27
		Delay cost per flight [e]	20.43	18.73	16.27	20.00	18.81
		Induced cost per flight [e] 3344.73 3412.44 3404.44 4211.46 3639.04
	i =15 W	GainRev% CHA% CAN% ỹi [minutes]	23.22% 20.46% 50.39% 40.52% 33.65% 20.48% 21.47% 22.32% 23.93% 22.24% 0.12% 0.11% 0.09% 0.16% 0.12% 0.28 0.26 0.22 0.28 0.26
		Delay cost per flight [e]	19.63	18.02	15.65	19.47	18.16
		Induced cost per flight [e] 3815.45 3823.97 3766.58 4525.33 4019.93
		GainRev%	22.79% 20.10% 49.38% 38.49% 32.69%
	Avg.	CHA% CAN%	20.56% 21.54% 22.33% 23.91% 22.27% 0.14% 0.12% 0.14% 0.26% 0.17%
		ỹi [minutes]	0.28	0.25	0.22	0.26	0.25
		Delay cost per flight [e]	19.26	17.70	15.11	18.50	17.58
		Induced cost per flight [e] 4004.29 3989.07 4120.91 5279.89 4418.42

Table 3 .

 3 8.7 -Computational result for robust FLA by different feasibility estimation methods

	Method	Soyster	Hoeffding	MC	Heuristic	HeuMod
	GainRev%	32.69% 33.73% 38.22% 39.04% 38.62%
	CHA%	22.27% 22.26% 11.08% 8.68%	9.42%
	CAN%	0.17%	0.12%	0.06%	0.05%	0.06%
	ỹi [minutes]	0.25	0.36	0.10	0.09	0.10
	Delay cost per flight [e]	17.58	25.33	6.81	6.28	7.36
	Induced cost per flight [e] 4418.42 4007.27 2067.06 1701.57 1871.62
	CPU Time [seconds]	286.46	668.78	25481.05 472.50	1776.32
	Min Diff.Proba.%	-	1.01%	0.00%	-1.83%	0.00%
	Max Diff.Proba.%	-	33.00% 2.20%	2.39%	6.74%

Table 3 .

 3

	Soyster	Hoeffding MC	Heuristic	HeuMod

8.8 -Computational result for Robust FLA by different feasibility estimation methods with different configurations of W i W i /Method

Table 3 .

 3 8.9 -Computational result for Robust FLA by different estimation methods with different configurations of when W i = 5 minutes

		Method/	1.00% 2.00% 3.00% 4.00% 5.00%
		GainRev%	31.32% 31.32% 31.31% 31.32% 31.33%
		CHA%	22.26% 22.26% 22.28% 22.25% 22.23%
	Hoeffding	CAN% ỹi [minutes] Delay cost per flight [e] Induced cost per flight [e] 4992.03 4991.94 4993.15 4989.47 4985.86 0.25% 0.25% 0.25% 0.25% 0.25% 0.34 0.34 0.33 0.33 0.33 23.74 23.65 23.42 22.94 22.86 CPU Time [seconds] 702.69 685.54 293.78 671.20 285.47
		Min Diff.Proba.%	1.01%	2.04%	3.09%	4.17%	5.25%
		Max Diff.Proba.%	1.01%	2.04%	3.09%	4.17%	5.26%
		GainRev%	31.62% 32.99% 33.64% 34.15% 34.73%
		CHA%	21.86% 19.05% 16.75% 14.93% 12.88%
		CAN%	0.23%	0.21%	0.21%	0.21%	0.21%
	MC	ỹi [minutes] Delay cost per flight [e] Induced cost per flight [e] 4867.50 4276.24 3991.47 3762.65 3506.85 0.22 0.17 0.14 0.10 0.08 15.21 12.00 9.58 7.29 5.42
		CPU Time [seconds]	59440.88 40787.45 32542.84 30854.94 19281.12
		Min Diff.Proba.%	0.03%	0.15%	0.03%	0.21%	0.18%
		Max Diff.Proba.%	0.20%	0.21%	0.14%	0.33%	0.43%
		GainRev%	30.77% 34.14% 34.32% 36.40% 37.70%
		CHA%	20.43% 15.85% 13.23% 10.98% 8.99%
		CAN%	0.30%	0.18%	0.21%	0.14%	0.10%
	Heuristic	ỹi [minutes] Delay cost per flight [e] Induced cost per flight [e] 5188.70 3801.41 3674.57 2775.55 2219.17 0.27 0.20 0.13 0.07 0.02 19.17 13.95 9.26 5.16 1.39 CPU Time [seconds] 818.12 806.13 366.96 844.04 265.08
		Min Diff.Proba.%	-0.28% -1.72% -0.39% -0.26% -0.34%
		Max Diff.Proba.%	-0.09% -0.15% -0.24% -0.07% -0.09%
		GainRev%	30.85% 32.98% 34.43% 34.45% 37.52%
		CHA%	21.03% 17.90% 14.07% 11.44% 9.47%
		CAN%	0.28%	0.20%	0.18%	0.21%	0.10%
	HeuMod	ỹi [minutes] Delay cost per flight [e] Induced cost per flight [e] 5173.48 4201.16 3642.31 3574.66 2295.92 0.30 0.24 0.16 0.08 0.02 20.81 16.81 11.34 5.87 1.69
		CPU Time [seconds]	1813.64 1713.43 747.28	4877.05 1336.60
		Min Diff.Proba.%	0.01%	0.08%	0.06%	0.04%	0.00%
		Max Diff.Proba.%	0.22%	1.70%	0.19%	0.23%	0.25%

Table 3 .

 3 8.10 -Continuous table of Table 3.8.9

		Method/	5.00% 10.00% 15.00% 20.00% 25.00%
		GainRev%	31.33% 31.32% 31.53% 31.92% 31.62%
		CHA%	22.23% 22.26% 22.28% 22.45% 22.62%
	Hoeffding	CAN% ỹi [minutes] Delay cost per flight [e] Induced cost per flight [e] 4985.86 4990.66 4895.68 4738.19 4874.92 0.25% 0.25% 0.23% 0.21% 0.22% 0.33 0.32 0.32 0.31 0.28 22.86 22.44 22.27 21.49 19.91 CPU Time [seconds] 285.47 525.61 279.67 478.21 261.64
		Min Diff.Proba.%	5.25%	11.10% 17.62% 24.69% 32.91%
		Max Diff.Proba.%	5.26%	11.11% 17.65% 24.79% 33.00%
		GainRev%	34.73% 36.69% 38.87% 39.01% 42.59%
		CHA%	12.88% 7.52%	3.08%	1.84%	1.31%
		CAN%	0.21%	0.20%	0.16%	0.17%	0.00%
	MC	ỹi [minutes] Delay cost per flight [e] Induced cost per flight [e] 3506.85 2637.60 1678.99 1620.62 176.50 0.08 0.02 0.00 0.00 0.00 5.42 1.06 0.16 0.09 0.21
		CPU Time [seconds]	19281.12 13750.77 5265.16 5436.28 2291.50
		Min Diff.Proba.%	0.18%	0.08%	1.01%	0.55%	1.53%
		Max Diff.Proba.%	0.43%	0.81%	1.51%	0.90%	2.20%
		GainRev%	37.70% 41.63% 42.33% 42.57% 42.70%
		CHA%	8.99%	4.15%	2.10%	1.36%	0.91%
		CAN%	0.10%	0.00%	0.00%	0.00%	0.00%
	Heuristic	ỹi [minutes] Delay cost per flight [e] Induced cost per flight [e] 2219.17 596.21 0.02 0.00 1.39 0.14 CPU Time [seconds] 265.08 230.04	0.00 0.00 290.07 68.21	0.00 0.00 187.07 115.86	0.00 0.00 127.77 63.38
		Min Diff.Proba.%	-0.34% -0.18% -0.39% -0.64% 1.28%
		Max Diff.Proba.%	-0.09% 0.08%	0.51%	0.49%	1.92%
		GainRev%	37.52% 40.74% 42.32% 42.55% 42.70%
		CHA%	9.47%	5.23%	2.12%	1.39%	0.91%
		CAN%	0.10%	0.02%	0.00%	0.00%	0.00%
	HeuMod	ỹi [minutes] Delay cost per flight [e] Induced cost per flight [e] 2295.92 984.17 0.02 0.01 1.69 0.35	0.00 0.00 294.29	0.00 0.00 193.89	0.00 0.00 127.77
		CPU Time [seconds]	1336.60 439.21	110.58	266.65	87.37
		Min Diff.Proba.%	0.00%	0.01%	0.08%	0.05%	1.28%
		Max Diff.Proba.%	0.25%	4.91%	0.51%	0.69%	1.92%

  1 or in n 2 . Constraint (4.28)-(4.29) make sure that the BBU-RFBs and MEC-RFBs can be installed only in nodes already storing RRH-RFB. Constraint (4.30) assures that the traffic assigned to a user has to be higher than a minimum value. Constraints (4.31)-(4.38) define the feasible domain of decision variables.

	A. User-Node Assignment
	Problem: Constraints
	(4.2), (4.3), and (4.30)
	Users	x un , t u	5G Nodes
				C. RFB Chain Construction
	B. RRH-RFB Allocation Problem: Constraints (4.4)-(4.7), (4.14)-(4.22), and (4.27)	y RRH nq	Problem: Constraints (4.8)-(4.13), (4.23)-(4.26), and (4.28)-(4.29)
		v BBU n1n2q	v MEC n1n2q
	BBU-RFBs		RRH-RFBs	MEC-RFBs

  43)

			y RRH unq ≥ δ|U |	(4.47)
	u∈U n∈N q∈Q
	y RRH unq CAP unq ≤	CAP RRH q	x RRH nq	∀n ∈ N	(4.48)
	u∈U q∈Q				q∈Q
	x BBU nq	≥	x RRH nq	∀q ∈ Q	(4.49)
	n∈N		n∈N
	x BBU nq		≤	x RRH nq	∀n ∈ N	(4.50)
	q∈Q		q∈Q
	x MEC nq	≥	x RRH nq	∀q ∈ Q	(4.51)
	n∈N		n∈N
	x MEC nq	≤	x RRH nq	∀n ∈ N	(4.52)
	q∈Q		q∈Q
	y MEC unq ≤ x MEC nq			∀u ∈ U, ∀n ∈ N, ∀q ∈ Q	(4.53)
	y MEC unq ≥		x MEC nq	∀n ∈ N	(4.54)
	u∈U q∈Q			q∈Q
	y MEC unq ≤ 1		∀u ∈ U	(4.55)
	n∈N q∈Q
	t MIN y RRH unq ≤ t u ≤	CAP unq y RRH unq	∀u ∈ U	(4.56)
	n∈N q∈Q				n∈N q∈Q
	y RRH unq ≥		x RRH nq		∀n ∈ N	(4.44)
	u∈U q∈Q	q∈Q
	y RRH unq ≤ 1			∀u ∈ U	(4.45)
	n∈N q∈Q
	y RRH unq ≤		U max q	x RRH nq	∀n ∈ N	(4.46)
	u∈U q∈Q	q∈Q
	Chenghao WANG	Contribution to robust network optimization

  110)

	n∈N q∈Q	x MEC nq	CAP MEC q t MIN	≥ δ|U |	(4.119)
	u∈U q∈Q
			y RRH unq ≥		x RRH nq	∀n ∈ N	(4.111)
	u∈U q∈Q			q∈Q
			y RRH unq ≤ 1	∀u ∈ U	(4.112)
	n∈N q∈Q
				y RRH unq ≥ δ|U |	(4.113)
	u∈U n∈N q∈Q
		x BBU nq	≥	x RRH nq	∀q ∈ Q	(4.114)
	n∈N		n∈N
			y RRH unq CAP unq ≤	CAP RRH q	x RRH nq	∀n ∈ N	(4.115)
	u∈U q∈Q				q∈Q
		x BBU nq	≤	x RRH nq	∀n ∈ N	(4.116)
		q∈Q		q∈Q
		x MEC nq	≥	x RRH nq	∀q ∈ Q	(4.117)
	n∈N		n∈N
		x MEC nq	≤	x RRH nq	∀n ∈ N	(4.118)
		q∈Q		q∈Q
	Chenghao WANG		Contribution to robust network optimization

  147) in which the objective function maximizes the number of users served by the activated RRH-RFBs. Constraint (4.156) indicates the minimum traffic demand required by users. Constraint (4.157) specifies that one user can be served at most by one 5G node. Constraint (4.158) gives an upper bound of served number of users by the activated RRH-RFBs. Constraint (4.159) denotes that the total served traffic at a 5G node by an RRH-RFB should not violate its capacity. Finally, constraint (4.160) defines the domain of decision variables y.< δ|U | , then P M is infeasible; otherwise P M is optimal. For the former case, Benders cuts, where U * nq is obtained from P M . Furthermore, consider the linear relaxation of P S (x RRH nq * ) (denoted by P SC (x RRH nq * )), if the maximum number of served users by the given activated RRH-RFBs is less than the desired covered number of users, then the master problem is infeasible as the solution cost of an LP problem is always less than or equal to the cost of its linear relaxation variant in sense of maximization optimization. In this case, the Benders feasibility cut (i.e., constraints (4.161) and (4.162)) can be added into the master problem such that at least one RRH-RFB assignment changed, and the user coverage is guaranteed by this new change of x RRH nq . Let us then consider the linear relaxation of P S (x RRH nq * ), denoted by P SC (x RRH nq * ), which can be written as below:Let λ un ≥ 0, ξ u ≥ 0, η n ≥ 0, π n ≥ 0, ∀u ∈ U, ∀n ∈ N be the dual variables associated with constraints (4.156), (4.157), (4.158), and (4.159), respectively. Then the dual form of P SC (x RRH nq * ), denoted by P DSC (x RRH nq * ), can be written as:

	x BBU nq [π n ≥ 0]	≤	x RRH nq y un x RRH nq	* CAP unq ≤	CAP RRH q	x RRH nq	*	∀n ∈ N ∀n ∈ N	(4.150) (4.167)
	q∈Q	q∈Q u∈U q∈Q				q∈Q
		x MEC nq	≥ y un ≥ 0 x RRH nq				∀q ∈ Q ∀u ∈ U, ∀n ∈ N	(4.151) (4.168)
	n∈N		n∈N
		x MEC nq	≤		x RRH nq				∀n ∈ N	(4.152)
	q∈Q n∈N q∈Q x k n∈N q∈Q x MEC q∈Q nq CAP MEC q t MIN nq ∈ {0, 1} If the solution cost of P S (x RRH nq * ), denoted by Z * ≥ δ|U | n∈N q∈Q P S (x RRH ∀n ∈ N, ∀k ∈ K, ∀q ∈ Q, nq * ), satisfies Z * P S (x RRH nq x RRH nq * =1:n∈N, q∈Q (1 -x RRH nq ) + x RRH nq * =0:n∈N, q∈Q min u∈U n∈N q∈Q λ un x RRH nq * CAP unq + ξ u t MIN u∈U x RRH nq ≥ 1 (4.161) (4.153) (4.154) * ) (4.169) + CAP RRH q x RRH nq * π n + min U max q , U nq x RRH nq * η n (4.170)
	where the objective function minimizes the total installation cost. Constraint (4.145) indicates that one 5G node can hold at most one RRH-RFB, one BBU-x RRH nq * =1:n∈N, q∈Q x RRH nq U * nq + x RRH nq x RRH nq min U max q , U nq ≥ δ|U | (4.162) s.t. π n CAP unq x RRH nq * + η n + ξ u + λ un ≥ 1 ∀u ∈ U, ∀n ∈ N (4.171) * =0:n∈N, q∈Q q∈Q RFB, and one MEC-RFB. Constraint (4.146) specifies the available number of each RFB module. Constraint (4.147) denotes the interference constraint among RRH-RFB of each type placed on the node. Constraints (4.148) and (4.153) specify the minimum number of activated RRH-RFBs and BBU-RFBs. Constraints (4.149)-(4.152) indicate the activation of BBU-RFBs and MEC-RFBs. Finally, constraint (4.154) defines the domain of decision variables x. The slave problem is instead: max u∈U n∈N y un (4.155) s.t. y un t MIN ≤ q∈Q x RRH nq * CAP unq ∀u ∈ U, n ∈ N (4.156) n∈N y un ≤ 1 ∀u ∈ U (4.157) min U max q , U nq x RRH nq * ∀n ∈ N (4.158) λ un ≥ 0 ∀u ∈ U, ∀n ∈ N (4.172) η n ≥ 0 ∀n ∈ N (4.173) π n ≥ 0 ∀n ∈ N (4.174) ξ u ≥ 0 ∀u ∈ U (4.175) Focus on the solution cost, denoted by Z * P DSC (x RRH nq * ), obtained from the dual form of relaxation variant of the subproblem, if we have Z * P DSC (x RRH nq * ) < δ|U | are the feasibility max y un holding, then P M is clearly infeasible. Then a Benders feasibility cut is generated (4.163) as below: u∈U n∈N s.t. [λ un ≥ 0] y un ≤ q∈Q x RRH nq * CAP unq t MIN ∀u ∈ U, ∀n ∈ N (4.164) u∈U n∈N q∈Q λ * un x RRH nq CAP unq + ξ * u + t MIN u∈U n∈N q∈Q
				q∈Q
	n∈N q∈Q [ξ u ≥ 0]	x RRH nq n∈N	min U max q y un x RRH nq * CAP unq ≤ , U nq ≥ δ|U | y un ≤ 1	CAP RRH q	x RRH nq	*	(4.148) ∀n ∈ N (4.159) ∀u ∈ U (4.165)
	u∈U q∈Q					q∈Q
	n∈N [η n ≥ 0] x BBU nq y un ∈ {0, 1} ≥ n∈N x RRH nq u∈U q∈Q y un ≤	min U max q	, U nq x RRH nq	*	∀q ∈ Q ∀u ∈ U, ∀n ∈ N, (4.149) (4.160) ∀n ∈ N (4.166)
	Chenghao WANG Chenghao WANG	Contribution to robust network optimization Contribution to robust network optimization

u∈U y un ≤

1

  Solve P M ; 2 if P M infeasible or unbounded then

	3	Stop;
	4 else
	5	isSolved ← False;
	6	repeat
	7 8 9	Obtain x RRH nq Solve P SC (x RRH * from solved P M ; * ); nq if Z * P SC (x RRH nq * ) < δ|U | then
	10	Add a feasible benders cut (4.176) into P M ;
	11	Solve P M ;
	12	else
	13 14	Solve P S (x RRH nq if Z * P S (x RRH nq * ) < δ|U | then * );
	15	Add the feasible benders cuts (4.161) and (4.162) into P M ;
	16	Solve P M ;
	17	else
	18	isSolved ← True;
	19	end
	20	end
	21	until isSolved ;

Table 4 .

 4 5.2 -Comparison on solved instances for each method

			M	AM	AMS
	N opt		6	142	212
	N inf eas	-	6	21
	N timeout	114	122	37
	N memout	150	-	-
	T opt avg (s)	737.98	126.18	45.99
	T inf eas avg	(s)	-	292.84	104.44
	different types of Reusable Functional Blocks (RRH-RFB, BBU-RFB and MEC-

  min Γδ rs + Noticing then that DEV rs (Γ, x)-primal is a feasible and bounded problem, strong duality can be exploited and states that DEV rs (Γ, x)-dual is also feasible and bounded and the optimal values of the two coincide problems. Following the Bertsimas and Sim procedure, the non-linear constraint (5.19) with the following robust version and with the related dual variables, getting the following compact model: rs ] v∈V a vr x vs ≤ ξ s , which indicates the power consumption for charging the different resource at each server s. The complete model that we consider is thus the following, which we denote by ROB-BLP-VP:

			η rsv	(5.26)
			v∈V
			v ∈ V	(5.27)
		δ rs ≥ 0		(5.28)
		η rsv ≥ 0		(5.29)
		δ rs ≥ 0	∀s ∈ S, ∀r ∈ R	(5.32)
		η rsv ≥ 0	∀s ∈ S, ∀r ∈ R, ∀v ∈ V	(5.33)
	Similarly to uncertain a vr in the objective function, we introduce an artificial
	decision variable ξ s such that [(P max s )/a min -P min s (i,j)∈L s∈S r=CPU P min s y s + ξ s + P ij w ij	(5.34)
	s.t.	P max s a rs -P min s	Γδ rs +

s.t. δ rs + η rsv ≥ ∆a vr x vs v∈V āvr x vs + Γδ rs + v∈V η rsv ≤ a rs y s ∀s ∈ S, ∀r ∈ R (5.30) δ rs + η rsv ≥ ∆a vr x vs ∀s ∈ S, ∀r ∈ R, ∀v ∈ V (5.31) n∈N P n z n + v∈V (ā vr x vs η rsv ) ≤ ξ s ∀s ∈ S (5.35) s∈S x vs = 1 ∀v ∈ V (5.36) Chenghao WANG Contribution to robust network optimization

Table 5 .

 5 4.1 -Computational resultsIDT

	I1	3322	2580	22.3	5.4
	I2	3194	2742	14.1	6.8
	I3	3157	2335	26.0	6.2
	I4	3552	2905	18.2	10.2
	I5	3513	2536	27.8	6.9
	I6	3402	2892	14.9	5.8
	I7	3475	2642	23.9	8.6
	I8	3362	3041	9.5	9.3
	I9	3595	2587	28.0	7.6
	I10	3488	2769	20.6	5.5
	I11	3302	2281	30.9	19.4
	I12	3260	2105	35.4	24.8
	I13	3512	2447	30.3	13.5
	I14	3396	2572	24.2	21.7
	I15	3471	2076	40.1	19.4
	I16	3395	2405	29.1	14.2
	I17	3185	2953	7.2	22.1
	I18	3338	2633	21.1	16.9
	I19	3198	2194	31.3	21.6
	I20	3056	2310	24.4	16.8

* (CPLEX) T * (MatHeu) ∆T * % ∆P * %

Table 5 .

 5 4.2 -Comparison of alternative exact search -power reduction

	ID	P * (MatHeu-RINS)	P * (MatHeu-HD)
	I1	5.4	2.1
	I2	6.8	2.6
	I3	6.2	3.8
	I4	10.2	6.3
	I5	6.9	8.5
	I6	5.8	3.9
	I7	8.6	10.2
	I8	9.3	3.3
	I9	7.6	4.8
	I10	5.5	1.9
	I11	19.4	13
	I12	24.8	11.5
	I13	13.5	16.3
	I14	21.7	8.4
	I15	19.4	23.2
	I16	14.2	5.9
	I17	22.1	13.2
	I18	16.9	10.7
	I19	21.6	11.6
	I20

  -t o i ) 2 + v 2 j (t -t o j ) 2 -2v i v j (t -t o i )(t -t o j ) cos θ (6.2) = v 2 i ([1 + ρ 2 -2ρ cos θ]t 2 -2[t o i + ρ 2 t o j -ρ cos θ(t o i + t o j

						)]t
						+ [t o i	2 + ρ 2 t o j	2 -2ρ cos θt o i t o j ])	(6.3)
	= v 2 i [1 + ρ 2 -2ρ cos θ] t -	1 + ρ 2 -2ρ cos θ i + ρ 2 t o t o j -ρ cos θ(t o i + t o j )	2
	+ v 2 i [t o i	2 + ρ 2 t o j	2 -2ρ cos θt o i t o j ] -v 2 i (t o i + ρ 2 t o j -ρ cos θ(t o i + t o j )) 2 (6.4)
	= v 2 i [1 + ρ 2 -2ρ cos θ] t -	1 + ρ 2 -2ρ cos θ i + ρ 2 t o t o j -ρ cos θ(t o i + t o j )	2
						+	v 2 j (t o j -t o i ) 2 sin 2 θ 1 + ρ 2 -2ρ cos θ	(6.5)
	≥ λ 2 v 2 j (t o j -t o
		t o i t i ) 2 + (d	t o j t j ) 2 -2(d	t o i t i )(d	t o j t j ) cos θ	(6.1)

Chenghao WANG

Contribution to robust network optimization

= v 2 i (t i ) 2 if sin 2 θ > 0,

  ) + d s (1 -ρ)) 2 , t = t i + d s /v i (v j (t o j -t o i ) + d s (1 -ρ)) 2 /ρ 2 , t = t j + d s /v j

		-	t o i -ρt o j 1 -ρ	2	(6.7)
	= min	  v 2 j (t o j -t o i ) 2 , t = t o i            v 2 i (t o j -t o i ) 2 , t = t o j (v j (t o j -t o i (6.8)
	= (v j (t o j -t o

i ) -d s |1 -ρ|) 2 / max(ρ 2 , 1) (6.9)

FL290 =

000 ft 2 1 nm = 1.852 km

3.6 -The Robust counterpart of FLA problem

The research work done in this thesis took place in Heudiasyc Laboratory, Université de Technologie de Compiègne from December 2017 to June 2021. It was funded by the French Ministry of High Education and Research.

At the end of this thesis, I would like to express my sincere gratitude to my supervisor Dr.

From a logical point of view, each 5G node is able to host one RRH-RFB, one MEC-RFB and one MEC-RFB. Moreover, it can pool also BBU-RFBs and MEC-RFBs from other nodes. On the other hand, the RFBs are organized in chains where each RRH-RFB is connected to a BBU-RFB, which is in turn linked to a MEC-RFB as shown in Figure 4.2.4. We then assume that the 5G node can provide a service to user if and only if there exists a complete RFB chain composed of one RRH-RFB, one BBU-RFB and one MEC-RFB, linked to each other in this order. The RFB chain is hence not constrained to be located on the same 5G node, but it can be realized across several nodes (e.g., each 5G node holds one RFB module from the complete RFB chain as shown in Figure 4.2.4). RFBs are characterized by their resource requirements (i.e., storage, processing). The requirements in terms of consumed resources by the RFBs are then used in this work to properly dimension the 5G nodes. Finally, we consider on additional classification of each RFB into Micro and Macro type, depending on the area size and user number that they may serve (this corresponds with a classification of 5G base stations into micro and macro).

Mathematical Formulation

The optimization model that we present is based on the following sets and indices: U : Set of users. N : Set of candidate 5G nodes. K: Set of RFB modules: RRH-RFB, BBU-RFB, and MEC-RFB. Q: Set of type for RFB chains/RFB modules: Micro, Macro. A k q : Available number of RFB module k in type q.

Chenghao WANG Contribution to robust network optimization U max q : Maximum number of users served by RRH-RFB in type q. CAP unq : Radio link capacity provided to user u at node n placed RRH-RFB in type q.

CAP k q : Maximum capacity managed by RFB module k in type q. CONF q : All the pairs of nodes that conflict for a RRH-RFB in type q.

t MIN : Minimum data traffic required by user. δ: Minimum fraction of users that has to be covered by the 5G service. c k q : Construction cost of each RFB module k in type q.

Reference Mathematical Model

In this section, we show the model presented in [START_REF] Chiaraviglio | Optimal design of 5g superfluid networks: Problem formulation and solutions[END_REF] that we have used as basis for our original developments.

Variables: t u : A continuous variable indicating the amount of downlink traffic served to user u ∈ U .

x un : A binary variable taking value 1, if the user u ∈ U is served by a RRH-RFB installed at node n ∈ N ; 0, otherwise.

y RRH

nq : A binary variable taking value 1, if the RRH-RFB of type q ∈ Q is installed at node n ∈ N ; 0, otherwise. v BBU n 1 n 2 q : A binary variable taking value 1, if an BBU-RFB of type q ∈ Q installed at node n 1 ∈ N serves the RFB chain originating from the RRH-RFB installed at node n 2 ∈ N ; 0, otherwise. v MEC n 1 n 2 q : A binary variable taking value 1, if a MEC-RFB of type q ∈ Q installed at node n 1 ∈ N serves the RFB chain originating from the RRH-RFB installed at node n 2 ∈ N ; 0, otherwise. φ un : A continuous variable indicating the amount of downlink traffic served to user u at node n, defined as φ un = t u x un . θ unq : An artificial binary variable indicating whether a user u at node n is served by the RRH-RFB of type q, defined as θ unq = x un y RRH nq . ϕ un 1 n 2 q : An artificial continuous variable indicating the amount of downlink traffic served to user u charged by a RFB chain which originates from the RRH-RFB installed at node n 2 and using as sink to the MEC-RFB installed at node n 1 , defined as

Chenghao WANG Contribution to robust network optimization NFV infrastructure in terms of computing power, storage and network topology.

In this example, we should map three service chains (denoted as c1, c2 and c3) into this NFV Infrastructure, each one with their own latency bounds. In total, we have three different VNFCs (v1, v2 and v3) and we assume that the traffic source for c 1 is the Sender S 1 , which is connected to router n 2 and injects a certain volume of traffic into the service chain towards v 1 . The VNFC v 1 processes the packets (for which it needs CPU, memory and disc) and forwards the processed traffic (which may have a different volume than the one injected) towards VNFC v 2 , which again processes it and forwards a certain volume to the Destination D 1 that is connected to router n 2 . Note that Figure 5.2.2 assumes additional source/sink nodes where traffic for a service chain is created/terminated, which are not explicitly mentioned in the model but they could be introduced by adding network nodes. The figure depicts an exemplary VNF placement and mapping into the physical Substrate Network. For example, the VNFC v 1 would be placed onto server s 3 , v 3 onto server s 4 and so on. Servers hosting no VNFC would be powered down (e.g., s 1 , s 2 or s 5 ) together with all the nodes not carrying any traffic (e.g., only n 1 in this case).

The optimization problem related to VNFC placement that we consider can be resumed as follows: given a 5G network interconnecting a set of servers, we want to decide how to establish a set of virtual chains in the network respecting the available resource budget of the servers and networks, while minimizing the overall power consumption.

In order to model the decision taken in the optimization problem, we introduce the following decision variables:

• variables y s ∈ {0, 1}, ∀s ∈ S representing the activation of a server (y s = 1 if s is turned on and 0 otherwise);

• variables x vs ∈ {0, 1}, ∀v ∈ V, s ∈ S representing the allocation of a VNFC v to server s (x vs = 1 if v is allocated to s and 0 otherwise);

• variables z n ∈ {0, 1}, ∀n ∈ N representing the activation of a node n (z n = 1 if n is turned on and 0 otherwise);

• variables w ij ∈ {0, 1}, ∀(i, j) ∈ L representing the activation of a link = (i, j) (g ij = 1 if = (i, j) is turned on and 0 otherwise);

Chenghao WANG Contribution to robust network optimization

Γδ rs + v∈V (ā vr x vs + η rsv ) ≤ a rs y s ∀s ∈ S, ∀r ∈ R (5.37)

x vs ≤ y s ∀s ∈ S, ∀v ∈ V (5.38)

v∈V a vr x vs ≤ a rs y s ∀s ∈ S, ∀r ∈ R (5.39)

(5.41)

(5.42) We decided that the chromosome of an individual corresponds with a valorization of the decision variables (y, x) (of ROB-BLP-VP): these variables are those used to represent whether a server is activated and whether a VNFC is allocated to a specific server. These two decisions are particularly critical for the problem and once their values have been fixed, we obtain an easier subproblem of (ROB-BLP-VP) . Specifically, once the value of the variables (y, x) is fixed, (ROB-BLP-VP) reduces to a kind of robust network flow problem and is easier to be solved by a state-of-the-art optimization solver and it is much faster to identify an optimal solution for the restricted (ROB-BLP-VP) problem with (y, x) fixed.

Fitness function. In order to assess the quality of an individual, a natural choice consists of adopting the objective function (5.1) of (ROB-BLP-VP) as fitness function. In this way, we can establish a very simple correspondence between the genetic algorithmic interpretation and the optimization model and it is immediate to evaluate how good is an individual.

Initial population.

The strategy that we explored to generate the initial group of individuals relies on the following principles: to generate an individual, we randomly activate a number σ < |S| of servers and then we randomly assign each VNFC in V to one single activated server, checking that the resource constraints (5.5) are not violated. In this way, we obtain a valorization (ȳ, x) of the server and allocation variables that we can then complete by solving the remaining subproblem of (ROB-BLP-VP) through a state-of-the-art solver. By this strategy, we can [START_REF] Blum | Hybrid metaheuristics in combinatorial optimization: A survey[END_REF]. The search is based on using the effective heuristic Relaxation Induced Neighborhood Search (RINS) (we refer the reader to [START_REF] Danna | Exploring relaxation induced neighborhoods to improve mip solutions[END_REF] for an exhaustive description of it). Specifically, given a partial solution (ȳ, x) of (ROB-BLP-VP) and (y TLR , x TLR ) an optimal solution of a tight linear relaxation (specifically, the optimal solution obtained by removing the integrality requirements on the binary variables and considering the basic model strengthened by the cuts identified by the state-of-the-art solver), we solve a subproblem of (ROB-BLP-VP) where the value of the j-th component of the vectors (y, x) is fixed according to the following two rules:

1. If (ȳ, x) j = 0 ∧ (y TLR , x TLR ) ≤ , then (y, x) j = 0;

The subproblem of (ROB-BLP-VP) subject to such variable fixing is then solved by the state-of-the-art solver, running with a time limit.

The overall pseudo-code of the hybrid genetic algorithm presented above is provided in Algorithm 7.

Algorithm 7: Adaption of the general GA algorithm 1 Generate a number n of individuals (each individual is defined by randomly activating σ < |S| servers and then randomly assign each VNFC in V to one single activated server, while checking that the capacity constraints are satisfied). These individuals constitute the starting population POP; 2 while a time limit is not reached do

, otherwise (6.12)

The two aircraft have a potential conflict at some intersecting waypoint o if and only if the minimum separation distance between them is less than minimum separation S. Then the probability of potential conflict is formulated as below:

The induced en-route ATFM delay ω ij of resolution for pairwise conflict is then formulated by: 6.14) where 1(x) A is a indicator function: 1(x) A = {1, x ∈ A; 0, x / ∈ A}, the second term denoted by 1(t msd ij -t o i + t o j ) [0,∞) specifies that the induced en-route ATFM delay of associated flight due to resolution should be positive meanwhile the potential conflict existed, 1(t o i -t o j ) [0,∞) denotes whether the aircraft i arrives latter than j at the potential conflict point o.

Assume that flight departure time delay follows a GMM (see Figure 6.1.3a),

Chenghao WANG Contribution to robust network optimization