A Matheuristic for Joint Optimal Power and Scheduling Assignment in DVB-T2 Networks †

: Because of the introduction and spread of the second generation of the Digital Video Broadcasting—Terrestrial standard (DVB-T2), already active television broadcasters and new broadcasters that have entered in the market will be required to (re)design their networks. This is generating a new interest for effective and efﬁcient DVB optimization software tools. In this work, we propose a strengthened binary linear programming model for representing the optimal DVB design problem, including power and scheduling conﬁguration, and propose a new matheuristic for its solution. The matheuristic combines a genetic algorithm, adopted to efﬁciently explore the solution space of power emissions of DVB stations, with relaxation-guided variable ﬁxing and exact large neighborhood searches formulated as integer linear programming (ILP) problems solved exactly. Computational tests on realistic instances show that the new matheuristic performs much better than a state-of-the-art optimization solver, identifying solutions associated with much higher user coverage.


Introduction
Television represents one of the most impacting inventions of contemporary times: developed between the end of the 19th century and the beginning of the 20th century, in the postwar period of World War II, it had rapidly spread, first in the US and then continued its diffusion in the rest of the world, bringing programs and shows into private households that have deeply influenced people's lives and habits.With the passing of the years, the number of television broadcasting companies has (greatly) increased, leading to a high demand and competition for obtaining scarce frequency resources that are fundamental for the effective operation of broadcasting.Due to such spectrum scarcity and to intrinsic technological features of analogue television broadcasting, which is known to need a considerable amount of frequency spectrum for broadcasting the content of a television channel with sufficient quality, it was considered very important to pass to a new generation of (terrestrial) television broadcasting characterized by more spectral-efficient digital transmissions.Thus, the study and introduction of digital terrestrial television (DTT) was strongly supported and the switch from analogue to digital television around the world has represented a fundamental passage to pursue a more efficient use of the spectrum and the support of a new range of higher quality television services [1,2].
It is interesting to note that many different DTT standards have been proposed and adopted in different countries.However, among such standards, Digital Video Broadcasting-Terrestrial (DVB-T) [3] is the most adopted around the world, being used in Europe, Asia, Africa and Oceania.The DVB-T standard was published in 1997 and is strongly based on the use of orthogonal frequency-division multiplexing (OFDM).OFDM is a modulation technique that, essentially speaking, subdivides the digital data stream into multiple streams that are transmitted on multiple frequency channels (subcarriers) that are close in the spectrum and are characterized by orthogonality, i.e. they do not interfere with one another despite the absence of guard bands between them (we refer the reader to [4] for an exhaustive introduction to OFDM).Thanks to its attractive features, DVB-T has rapidly spread, but after a few years the need for an updated standard able to follow the evolution of next generation telecommunications services, in particular, by offering improved performance in terms of data rates and spectral efficiency, arose.As a consequence, in 2009 the second generation of DVB-T, namely DVB-T2, was released and published [5], offering improved features, such as more sophisticated modulation and coding schemes that can support up to 50% higher data rates and more flexible bandwidth choices (see e.g., [6]).
Passing from the first to the second generation of DVB-T has represented an important technological evolution that, in Europe, has been highly pushed by the European Commission.Indeed, in the European Union, the switch from DVB-T to DVB-T2 is expected to grant a capacity expansion of the digital television system of not lower than 25%, thus allowing a higher number of broadcasters to be active and compete in the market, while also enhancing the pluralism of information [7].The switch to DVB-T2 will require the television broadcasting companies that are already active in the market to reconfigure their networks, since they will have to embrace the innovative technological features of the new standard.Moreover, there will be the need to configure from scratch the new networks installed by new broadcasting companies entering the market.
The technological context of the switch to DVB-T2 is attracting new interest in the development of DVB-T network design tools, especially those including mathematical optimization features, that can improve the outcome of the design.In this work, we thus address the question of developing an innovative optimization algorithm to tackle the complex mathematical models associated with the design of DVB-T2 networks.More in detail, the contributions of our work are the following:

•
As first step, we show how a mathematical optimization model for DVB-T design can be derived.We address the two crucial decisions of establishing which DVB broadcasting station must serve a specific portion of territory, through binary decision variables, and the power setting of the stations, expressed as continuous decision variables.Furthermore, we include the design of the horizontal antenna diagram of the transmitters.The resulting optimization model is a Mixed Integer Linear Programming (MILP) problem, which includes canonical signal-to-interference ratios to assess service coverage of the territory.

•
We provide a discussion about the mathematical weakness of the MILP model, due to a typically ill-conditioned coefficient matrix and the presence of the notorious big-M coefficients, and we discuss how to strengthen it through the discretization of the power emissions of the stations, deriving a purely binary power-indexed model, and how to include further strengthening valid inequalities expressing conflicts between station power configurations and territory coverage.

•
In order to solve the mathematically stronger yet still complex power-indexed model, we propose a new matheuristic, i.e., a heuristic algorithm based on combining mathematical programming techniques and metaheuristics (see e.g., [8][9][10]).Specifically, our new matheuristic is based on combining a Genetic Algorithm (GA) with variable fixing heuristics exploiting a suitable (tight) linear relaxation of the problem and an Integer Linear Programming (ILP) improvement heuristic.The purpose of the GA is providing an efficient exploration of the solution space associated with the discrete power emissions of the stations.The relaxation-based fixing heuristic and the ILP heuristic have the task to pursue the improvement of the best solution found by the GA, in particular, by operating a so-called exact large neighborhood search: In an exact search, the search is modeled as an ILP problem which is tackled by exploiting the computational power of a modern state-of-the-art solver [11].

•
We report results that we obtained from computational tests considering realistic DVB-T data.
The results clearly indicated that our new matheuristic offers a superior performance with respect to a state-of-the-art ILP solver, returning solutions associated with much better values of the objective function expressing service coverage of the users.
We remark that this paper is an extended journal version of [12].The remainder of the paper is organized as follows: in Section 2, we introduce an optimization model for DVB-T2 network design, while in Section 3, we discuss how to mathematically strengthen it; in Sections 4 and 5, we introduce the new matheuristic and report the computational tests, respectively.At last, we discuss conclusions and future directions of research in Section 6.

Related Works
The use of mathematical optimization techniques in the context of wireless network design is not new and can be traced back to at least twenty years ago, as highlighted and surveyed for quite a wide range of problems in [13] and in [14].Also, the adoption of bio-inspired and genetic heuristics is not new: if we focus on genetic and evolutionary algorithms, we can report the remarkable cases of: [15], which addresses the problem of the optimal positioning of base stations in a mobile network, encoding the location of base stations in the chromosomes of the genetic algorithm; [16], which addresses the decision problem of how establishing the optimal assignment of users to deployed transmitters, in particular in the context of WiMAX networks, encoding the assignment in the chromosomes; [17], which focuses on the frequency assignment problem (FAP), proposing a permutation-based genetic algorithm to solve minimum span and fixed spectrum variants of the FAP; [18], which focuses on the problem of setting the power emissions of base stations, proposing two distributed power control algorithms that are based on evolutionary computation techniques to fast solve the linear equation systems associated with power updates of the stations; [19], which proposes a genetic algorithm for addressing the joint problem of power, frequency and modulation scheme assignment in fixed networks based on the WiMAX technology.Moreover, it is interesting to note that other works have also tried to tackle sources of data uncertainty within wireless network design problems, such as [20,21], which adopt a stochastic optimization approach to find a robust location plan of base stations to deal with fluctuations in traffic demand, [22], which deals with the stochastic scheduling of 5G multimedia services and [23][24][25], which propose stochastic programming and robust optimization approaches to deal with signal propagation uncertainty of wireless technologies in real-world environments.
After having discussed the previous works, it is very interesting to note that, to our best knowledge, in literature there are no matheuristics proposed for DVB network design problems.More specifically, it is also very important to remark that there are no cases of matheuristic expressly designed for tackling the numerical difficulties (ill-conditioned matrices and big-M coefficients) associated with classical optimization models for wireless network design problems (we refer the reader to Section 3 for a discussion about such numerical difficulties).Our new matheuristic, which combines genetic algorithms and ILP heuristics, is thus the first of its genre in the wireless network design context that we study.
Within the DVB network design context, we can position our matheuristic between exact optimization algorithms, guaranteeing, at least theoretically, convergence to an optimal solution, and (bio-inspired) heuristics, which instead produce solutions without quality guarantees.Concerning exact algorithms, we can cite the case of [26], which proposes an innovative cutting-plane method based on generalized-upper bound cover inequalities that eliminates sources of numerical instabilities, and [27] which proposes a packing reformulation of the single frequency network design problem arising in DVB-T.Concerning heuristics, besides the already discussed case of [19], we can cite the cases of: [28], which proposes particle swarm optimization and adaptive simulated annealing algorithms for optimizing the coverage of a DVB-T/H network; [29], which proposes a genetic algorithm to solve a service coverage problem of single frequency digital radio and television networks, considering also the energy efficiency and carbon footprint; [30], which attempts to reduce interference in single frequency networks, by optimizing the orientation of transmitters exploiting a simulated annealing algorithm.
In contrast to exact and heuristic approaches, our matheuristic has the merit of running fast, while at the same time exploiting the precious information that a mathematically stronger model of the problem can provide.In particular, such information is obtained from suitably strengthened linear relaxations and can be used to operate a better setting of the values of the decision variables (fixing) within the genetic algorithm and the ILP improvement heuristics.Furthermore, by exploiting linear relaxations, we can also compute an optimality gap indicating how far the best feasible solution found is far from the optimum and thus providing an indication of how good are produced feasible solutions.
To our best knowledge, there is just one work that tried to combine a metaheuristic with an exact improvement heuristic for the design of DVB networks, namely [31], which proposes an improved ant colony algorithm.However, in contrast to our present work, it is important to highlight that [31] considers an approximated model where service coverage is provided by one single transmitter (differently from the present work, where we consider the combined effect of multiple serving transmitters as detailed in the DVB-T technical specifications) and was not developed with the aim of sorting out numerical instabilities that characterize wireless network design.

DVB-T2 Network Design
In this section, we derive the mixed integer linear programming model that we use as basis for our new matheuristic algorithmic proposal.To this end, we first describe the essential elements constituting the DVB-T2 network that we consider: the network is made up of a set of DVB-T transmitting stations S that have the task of broadcasting digital television signals to users located in a region.Each station s ∈ S is characterized by a number of parameters: primarily, its location and then its radio-electrical parameters, such as the emitted power and the channel used to transmit.In accordance with recommendations and regulations released by telecommunication authorities, like the Italian Authority for Telecommunications (AGCOM) [32], the region to be covered with DVB-T services is partitioned into a raster of pixels corresponding to small squared portions of territories.The center of each pixel is called a testpoint (TP) and is taken as representative of all points in the pixel, in particular concerning propagation behavior of signals throughout the pixel.The set of all testpoints included in the region is denoted by T.
The Wireless Network Design problem (WND) is an optimization problem that generally consists of choosing the parameters of a set of wireless transmitters (e.g., their location and radio-electrical configuration) with the aim of providing service coverage to a set of users while maximizing a revenue function (see e.g., [13,14,26,33]).Notwithstanding the quite high number of parameters that characterize a transmitter and that can be configured, all works considering the WND focuses just on a subset of them.In particular, the vast majority of works include the setting of the power emission of transmitters and the assignment of served users to deployed transmitters as decision variables, since they constitute critical decision aspects in the design of a wireless network (see, for example, [7,26,27,34] related to DVB, [35,36] related to 5G, [23,37] related to FTTx, [38,39] related to mesh networks, [40,41] related to UMTS, [33,42] related to WiMAX, [43,44] related to WLANs and [13,14,[45][46][47][48] related to other kind of wireless network technologies).
If the optimization problem requires to decide the power emitted by transmitters and to choose the assignment of users to transmitters, it is obtained the so-called Scheduling and Power Assignment Problem (SPAP), which is NP-hard [27] and is considered a crucial WND problem in a hierarchy that has been developed in [27,33,49].
Before introducing the decision variables that are included in the optimization model, we make additional considerations about the power emission of DVB-T stations: such stations typically correspond with very powerful transmitters that have the task of covering vast portions of territory and host antennas whose horizontal radiation patterns can be customized (see e.g., [50]).Such customization of the horizontal antenna diagram typically considers 36 directions associated with angles from 0 to 350 degrees separated by a 10-degree interval and, from a modeling point of view, can be represented as the possibility of specifying the power emitted by a station s ∈ S in each direction d ∈ D = {1, 2, . . ., 36}.As noted in [50], the power emitted by the same station in different directions cannot exhibit too large differences.In particular, the difference between adjacent directions cannot exceed a threshold denote by ∆P.For each station s ∈ S and TP t ∈ T, we denote by δ(s, t) ∈ D the (unique) direction on which station s emits signals to t.
To model the decision of setting the power emission of the transmitters and establishing the assignment of users to transmitters in the context of DVB-T2 networks, we define two families of decision variables: • a continuous decision variable p sd ∈ [0, P max ] that models the power emission of each transmitter s ∈ S in each direction d ∈ D; • a binary decision variable x ts ∈ {0, 1} modeling the assignment of a TP to a station and such that: for every TP t ∈ T and every s ∈ S.
Concerning the evaluation of service quality, the first consideration to be made is that every TP t ∈ T is able to pick up signals from all stations s ∈ S and the power P t (s) that t obtains from s is equal to the power emission of s multiplied by a coefficient a ts ∈ [0, 1] (i.e., P t (s) = a ts • p s δ(s,t) -note that here, we consider the specific direction δ(s, t) on which s emits to t).The factor a ts is commonly called fading coefficient and summarizes the power reduction experienced by a signal that propagates from s to t [51].
We have written that a TP picks up signals from all stations.However, in the case of DVB-T networks, only a subset of signals is useful for guaranteeing the service, while all other signals are interfering and thus deteriorate the quality of service.Specifically, in DVB-T, a TP has to decide when to start a time window for picking up signals: all signals received within the time window can be summed up and increase the quality of service obtained, whereas all signals received outside the window are interfering and decrease service quality.We refer the reader to [27,52] for a detailed discussion of how the time window mechanism works.
We note that the start of the time window associated with a TP can be positioned in a very high number of points along the time axis, depending on the time sensitivity of the DVB system.However, it is common to let the window start when the signal of a station is received [27,32].In the case of our optimization model, we follow this practice and thus the number of possible starts of the time window of a TP t ∈ T equals |S|, i.e., the number of stations.We call server or serving station of t, the station s whose signal reception is used as starting time of the time window of t.
For a given TP t ∈ T and serving station s ∈ S of t, we denote by U(s, t) ⊆ S the subset of stations that emit useful signals for t and by I(s, t) ⊆ S the subset of interfering stations.These two subsets form a partition of the station set S, i.e., S = U(s, t) ∪ I(s, t) and U(s, t) ∩ I(s, t) = ∅.Once the subsets of useful and interfering signals are established, a TP t is considered served by station s if its Signal-to-Interference Ratio (SIR) is above a threshold φ > 0. The SIR is the ratio of the sum of useful powers to the sum of interfering powers [51,52]: In the previous inequality, N > 0 is the system noise.Every TP t that is covered with service generates a revenue r t > 0 (typically, the higher the number of users located in t the higher r t ).By means of simple linear algebra operations, we can reorganize Equation (1) into the inequality that follows, which is usually named SIR inequality: As one of the decisions to be considered in the problem is establishing which is the serving station of a TP, we must include one SIR inequality for every potential serving station s ∈ S of every TP t ∈ T. As a consequence, for every t, we face the following disjunctive constraint that contains one SIR inequality for each station s: The previous single disjunctive constraint can be reformulated as a set of linear constraints, according to a standard MILP approach (see, for example, [26,53]), which is based on defining a large constant M (commonly known as big-M coefficient) and exploiting the service variable x ts as a binary variable that activates or deactivates the constraint.More in detail, Equation ( 3) can be replaced with the following set of modified SIR inequalities (one for each potential server s ∈ S) that also include the big-M coefficient and the variables x ts : We call such modified inequality SIR constraint and we remark that it represents a fundamental component of every WND problem that considers coverage of users by means of SIR quantities.It is straightforward to verify that in Equation (4), when x ts = 1, we obtain the original SIR inequality without a big-M coefficient, and such an inequality must be fulfilled, thus it is required to serve TP t through station s as a server.Instead, when x ts = 0, the large value of the big-M coefficient is added to the left-hand-side of the SIR inequality, making it satisfied for any setting of power variables of the stations and thus making the constraint redundant.
We provide an overview of the system elements and notation introduced above in Table 1.We can then proceed to define the problem of designing a DVB network as follows.
Definition 1.The DVB Network Design Problem (DVB-ND): Given a set of stations S, a set of TPs T, the set of directions D and δ(σ, t) ∀t ∈ T, s ∈ T, the fading coefficients a ts ∀t ∈ T, s ∈ T, the testpoint population r t ∀t ∈ T, the maximum power emission P max of each station, the system noise N, the SIR threshold δ and the antenna adjacent direction power threshold ∆P, the DVB network design problem consists of establishing the power emission of each station s ∈ S and the serving station of each TP t ∈ T, so that the total population covered with service is maximized, while the SIR constraints of served TPS are satisfied, each TP is served by at most one station and the power limits of the stations and antenna technological constraints are respected.
The DVB-ND problem can be modeled as a mixed integer linear programming of which we detail below the feasibility constraints and objective function.Besides the SIR constraints in Equation ( 4), we must ensure that each TP t ∈ T is served by at most one station, namely: Also, we need to model the power adjacency constraints of antennas, ensuring that adjacent directions stay under the power difference threshold (note that here we assume that d = 36 + 1 = 1, since, in an arc of 360 degrees, by summing the 10 degrees of direction d = 1 to the 360 degrees of d = 36 we obtain 10 degrees and thus direction d = 1): Finally, the objective function expresses the maximization of total revenue obtained from serving testpoints, namely: Resuming all the elements, we obtain the following MILP model that we shortly indicate as DVB-MILP:

Strengthening the Initial MILP Formulation
DVB-MILP represents a straightforward model for DVB design that is based on the direct inclusion of the SIR constraints, as made in most other works addressing the WND (e.g., [13,14,26,33]).Nevertheless, such straightforward modeling presents remarkable drawbacks:

•
The fading coefficients typically lead to (very) ill-conditioned coefficient matrices, since they commonly vary in really wide ranges, and this leads to numerical instabilities even in state-of-the-art optimization software used to solve the problem (see [54] for a detailed discussion).

•
The big-Ms constitute notorious coefficients whose inclusion in a model is well-known to lead to weak bounds that badly affects the performance of optimization software (see e.g., [55]).

•
Due to the combined presence of the fading and big-M coefficients and their effects on the solution process, solutions returned by solvers typically contains errors under the form of SIR constraints that are recognized as satisfied when they actually are not; coverage plans can thus actually be wrong (see [26,46,47] and especially [54] for a more specialistic discussion, attributing errors to floating-point arithmetic used in solvers).
As a consequence, if we use DVB-MILP, we can hope to get reliable and optimal solutions only for small-sized instances, while, in the case of large realistic instances, even getting solutions of reasonable quality can represent a major challenge for state-of-the-art optimization software.It is relevant to remark that, though such issues are known, only few works have tried to tackle them (we refer the reader to [26,33] for a survey of such works).Our new matheuristic addresses the highlighted difficulties, in particular by also using a strengthened optimization model obtained according to the procedure presented below.
The first step to obtain a model mathematically stronger than DVB-MILP is to discretize the power emissions of stations: we define a set P = {P 1 , . . ., P m } of power values that represent feasible power emissions of any station and we introduce the set of corresponding power level indices L = {1, 2, . . ., m}.The second step consists of replacing each continuous power variable p sd with the linear combination of a set of new binary variables z sdl ∈ {0, 1} that multiplies the feasible power values of P: A variable z sdl is such that: 1 if s ∈ S emits with power value P l in direction d 0 otherwise.
Moreover, we must impose that every station s ∈ S emits by one single power value.Therefore, for every station s and direction d, it is necessary to add the following constraint ∑ l∈L z sdl ≤ 1 for each station s and direction d (note that having ∑ l∈L z sdl = 0 corresponds to a non-activated station emitting at zero power).Such constraints go under the name of Generalized Upper Bound (GUB) constraints.
By exploiting the new binary power emission variables, we obtain the following new purely binary version of SIR constraints: We denote by DVB-PI the Power-Indexed (PI) version of DVB-MILP that employs the binary variables z sdl and includes the modified SIR constraints (8) and the corresponding GUB constraints.
Exploiting the presence of the GUB constraints, following the approach initially proposed in [26,33], we can replace the SIR constraints in Equation ( 8) with a set of GUB cover inequalities.For an exhaustive introduction to simple and GUB cover inequalities and to their GUB version, we refer the reader to [53,56].Here, we just concisely remind major well-known results about cover inequalities: if we are given a knapsack constraint ∑ j∈J a j x j ≤ b with a j , b ∈ R + and x j ∈ {0, 1}, ∀j ∈ J, then we can replace it with its cover inequalities ∑ j∈C x j ≤ |C| − 1, where C is a cover.A cover represents a subset of indices C ⊆ J with the property that the sum of the corresponding coefficients a j , j ∈ C does not satisfy the knapsack constraint (i.e., ∑ j∈C a j > b).Each cover inequality associated with a cover C can therefore be interpreted as a combination of binary variables x j that cannot be set to 1 at the same time and thus we can activate at most |C| − 1 of them.The GUB cover inequalities represent a stronger version of the simple cover inequalities, which one can define exploiting the presence of GUB constraints.
Proceeding as in [26,33], we can define the general form of the GUB cover inequalities through which we can replace constraints in Equation ( 8), obtaining a stronger formulation: in which Σ ⊆ S is a subset of stations emitting useful signals towards TP t and Γ ⊆ S\Σ is a subset of interfering stations.Moreover, (q 1 , . . ., q |Γ| ) ∈ L I (t, Σ, λ, Γ), with L I (t, Σ, λ, Γ) ⊆ L |Γ| represent the subset of interfering levels of interfering stations in Γ that destroy the service of t provided by the subset of serving stations Σ, emitting with power levels corresponding to indices λ = (λ 1 , . . ., λ |Σ| ).
Intuitively, for a fixed TP, subset of serving stations Σ and subset of interfering stations Γ, a GUB cover inequality is built by fixing a power setting of the serving stations and defining a power setting of the interfering stations that compromise the coverage of the considered TP.We remark that in what follows we do not completely replace the binary SIR constraints in Equation ( 9), but we add the GUB constraints that consider one single useful station and one single interferer to DVB-PI.Indeed, as shown in [26,33], such a GUB provides a complete convex hull characterization for the single-server-single-interferer case and, even in the case of multiple serving stations and multiple interferers, still provide (tight) valid inequalities.
We further strengthen the mathematical formulation by including additional valid inequalities that are aimed at modeling the presence of couples of SIR constraints that consider only two stations and that cannot be satisfied at the same time by any power configuration.In a more formal way, let us consider two SIR constraints corresponding with TPs t 1 and t 2 served by two stations s 1 and s 2 : The previous constraints respectively represent TP t 1 served by s 1 and interfered by s 2 and TP t 2 served by s 2 and interfered by s 1 .If there exists no assignment of discrete power levels in P to s 1 , s 2 that can satisfy the two previous constraints, then coverage variables x t 1 s 1 and x t 2 s 2 cannot be activated at the same time and the following is a valid inequality: We can identify the inequalities in Equation ( 10) through pre-processing and include them in the model, obtaining a remarkable strengthening, as discussed in [57].
We denote by DVB-PI+, the power-indexed model DVB-PI strengthened by including the power-indexed GUB cover inequalities in Equation ( 9) and the conflict inequalities in Equation (10).

A Matheuristic for DVB Design
Though DVB-PI+ constitutes a binary linear program that can be easily passed to any optimization software to be solved, even a sophisticated last generation solver like CPLEX [58] may experience a very hard time just trying to identify solutions of reasonable quality.Furthermore, CPLEX also exhibits a very slow convergence to optimal solutions.In order to tackle such computational issues, we propose to solve DVB-PI+ through a new matheuristic algorithm that combines a Genetic Algorithm (GA) with relaxation-based variable fixing and exact large neighborhood searches.
GAs correspond to heuristics for optimization problems that draw inspiration from the evolutionary dynamics of populations and have known a wide success in many different application contexts.For a thorough introduction to GAs, we refer an interested reader to [59][60][61].Any GA is based on maintaining a population, in which every individual provides a feasible solution to the optimization problem.The feasible solution is encoded in the chromosome of the individual and a fitness function is used to evaluate how good the chromosome is, and thus, the corresponding feasible solution.During the execution of the GA, the population evolves, iteration after iteration, by means of procedures that resemble the process of natural selection, implementing the crossover of individuals' chromosomes, chromosome mutation and the death of individuals.
In our case, we strengthen the performance of a GA by combining it with relaxation-based variable fixing, in which the value of decision variables is set a-priori considering the value of variables in an optimal solution to a suitable (tight) linear relaxation of the problem, followed by exact large neighborhood searches, where the search for solutions of improved quality in a large neighborhood is expressed as an optimization problem solved exactly through a state-of-the-art solver.The rationale at the basis of the approach is that, while a state-of-the-art solver may find difficulties in solving a complete hard optimization problem, it may instead tackle effectively and efficiently suitable subproblems of it, obtained by fixing variables and exploring related large neighborhoods.We were attracted by the possibility of defining a matheuristic by combining a GA with tight linear relaxations and exact ILP search, since such an original combination allows the ability to explore the solution space of DVB-PI fast and effectively, exploiting the valuable tight linear relaxations associated with the strengthened formulation DVB-PI+.
The overall scheme of the GA that we consider is provided in Algorithm 1.
1: Creation of the initial population 2: while an arrest condition is not satisfied do Death of a portion part of the population 7: end while 8: Improvement by ILP Heuristic We now proceed to discuss in detail the features of all the steps of the algorithm.

Representing the Individuals
In our GA algorithm, a chromosome encodes the power setting of the DVB stations and therefore coincides with a power vector p of size |S| • |D| that sets the emission of every station in every direction.A generic position (s, d) in the chromosome specifies the discrete power level p l ∈ P that station s ∈ S emits in direction d ∈ D.
After having completely characterized a chromosome/power vector, it is still necessary to establish how to choose the serving station of a TP and thus set the values of the variables x ts .To this end, we compute the signal-to-interference ratio SIR ts (p) for the power vector p and use the distinction between useful and interfering signals that is originated by selecting s as serving station of t.We then introduce Σ(t, p) ⊆ S to denote the subset of serving stations that are able to provide service coverage to TP t for the power vector p by satisfying the corresponding SIR inequality, i.e., Σ(t, p) = {s ∈ S : SIR ts (p).≥ δ} When Σ(t, p) is not empty, we select as the serving station σ of t, the station σ ∈ Σ(t, p) that grants the highest SIR value SIR ts (p).In terms of the binary variables representing the server-TP assignment, we therefore impose x tσ = 1 and x ts = 0, ∀s ∈ S \ {σ}.

Fitness Function
As a fitness function, it is natural to use the coverage that is provided by an individual as a function of the power vector/chromosome p.We denote such function by COV(p) and it returns the sum of the population of TPs for which the SIR inequalities are satisfied and whose serving stations are identified following the steps illustrated earlier.

Initial Population
We generate the initial population through inclusion of power vectors in which one single station is activated in only one direction and, for every station, we define one power vector for every discrete power value excluding the zero value.In a more formal way, for every s ∈ S and d ∈ D, the initial population POP includes the vectors provided below: (0, 0, . . ., p sd = P 1 , . . ., 0, 0) (0, 0, . . ., p sd = P 2 , . . ., 0, 0) . . .

Such vectors provide |S| • |D| • |P | initial individuals.
Moreover, besides such elementary individuals, we enrich the population through more complex individuals with better fitness, identified by a variable-fixing-based procedure defined as follows.Let DVB-PI+ RLX be the linear relaxation of DVB-PI+ and (x, z) RLX its optimal solution.We exploit DVB-PI+ RLX as starting point to fix the value of a subset of variables in DVB-PI+ and obtain a smaller version of the problem that is easier to solve.We denote by DVB-PI+ FIX the problem where some variables are fixed and the strategy for fixing consists of setting to 1 the variables whose value in DVB-PI+ RLX is sufficiently close to 1.In a more formal way, we adopt this fixing rule: where > 0. The rationale at the basis of this approach is that, if the value of a variable coming from a (tight) continuous relaxation of the problem (i.e., a version of the problem where the integrality requirement on variables is dropped and the mathematical strength of the relaxation is improved by adding valid inequalities) is sufficiently close to an integer value, then we have a good indication that the value of such variable should be set to the close integral value in a solution of good quality.
For the proposed fixing, we focus just on the binary power emission variables since, due to the presence of the GUB constraints, they are particularly critical and, once that we set z sdl = 1 for some triple (s, d, l), we know that we can set z sdλ = 0 for all other triples (s, d, λ) with λ ∈ L\{l}.By adopting this fixing strategy, we obtain DVB-PI+ FIX , where we do not need to decide the power setting of some stations in some direction anymore and thus we face a smaller and, in general, easier to solve problem.DVB-PI+ FIX is solved through CPLEX, with an imposed time limit, and all the feasible solutions identified during the solution process are added to the initial population.

Selection
We select the individuals to be combined to give birth to the new generation by operating a tournament selection: given the current population POP, as first step we create k > 0 subgroups by random selection of α • |P| individuals from POP, with α ∈ (0, 1).As second step, we extract m < α • |P| individuals offering the highest fitness value in every group.These extracted individuals are then combined to create the new generation by means of crossover.

Crossover, Mutation and Death
The individuals selected in the previous step are randomly associated in couples, so as to form k • m/2 pairs.After pairing, each couple generates two offspring through chromosome crossover.More in detail, let us consider the crossover of a pair (the parents) corresponding with power vectors p 1 , p 2 ; the crossover combines the power levels in the same position of p 1 , p 2 for generating two offspring characterized by new power vectors p 3 , p 4 that hopefully are associated with higher fitness value.
With the aim of evaluating the effect of crossover, we define the measure ∆COV(p, p s = P l ) ∈ Z to represent the variation in the number of covered users that is caused by changing the power value p s in position s of a power vector p to a value P l ∈ P, while keeping all other power values unchanged.Exploiting measure ∆COV(p, p s = P l ), we operate the following crossover, which tries to make p 3 the best individual in the offspring.When the crossover starts, p 3 and p 4 have all elements equal to 0. Then, by following an increasing order of the index s from 1 to |S| and, for each s, of the index d from 1 to |D|, each value 0 inherits the power value of one of the two parents in the same position.More in detail, let us assume to focus on the crossover procedure for a position (σ, δ) with σ ∈ {1, . . ., |S|} and δ ∈ {1, . . ., |D|}: For values of s ∈ {1, . . ., |S|} such that s < σ and, for fixed σ, for values of d ∈ {1, . . ., |D|} such that d < δ, the crossover has been operated and thus the vectors p 3 , p 4 include power levels inherited by the parents p 1 , p 2 ; in contrast, positions corresponding to values s ≥ σ and, for fixed σ, to values d ≥ δ have not yet been processed and are thus still equal to 0.
The inheritance of power values of p 3 and p 4 from their parents p 1 and p 2 is set according to these rules: sd otherwise which lead p 3 to inherit the power values offering the best variation in coverage ∆COV.
In addition to crossover, we also provide for varying the values of power vectors through mutation, since it allows the ability to more effectively explore the set of solutions and helps to avoid the trap of local optimal solutions.To implement mutation, at every iteration, we select γ • |POP| with 0 < γ < 1 individuals in a random way.After this, we randomly select |L| • |D| power levels and we change them to the immediately lower power value included in P. The aim of this mutation procedure is to create individuals that grant the same coverage but using lower power emission: this is beneficial since it allows the ability to reduce interference among stations.When some crossover or mutation implies a violation of the constraint imposing a maximum difference between powers in adjacent directions of a station, we reduce the power to the closest value in P that allows the ability to satisfy the constraint.
As the last step, some individuals die and are excluded from POP.In particular, we remove the 2 • k • m/2 individuals with lowest fitness function values.

ILP Improvement Heuristic
Using the best solution found by the GA as basis, we attempt to find a better solution by running an ILP heuristic that corresponds to executing a very large neighborhood search in an exact way, i.e., modeling the search as an ILP problem solved through a state-of-the-art solver [11].This is based on the observation that, in contrast to the complete hard-to-solve problem, suitable subproblems of the complete problem, derived through fixing the value of a subset of variables, may be efficiently solved by a state-of-the-art solver.
In the ILP heuristic, we use a modified Relaxation Induced Neighborhood Search (RINS) (see [62] for a full description of the algorithm).Let ( x, z) be a feasible solution of the strengthened power-indexed problem DVB-PI+ and let (x TLR , z TLR ) be an optimal solution of the linear relaxation, strengthened by the cuts found by CPLEX in the root node of the branch-and-bound tree.Moreover, let ( x, z) j , (x TLR , z TLR ) j denote the j-th component of the vectors.
The modified RINS that we use and that we denote by mod-RINS corresponds to solving a subproblem where we a-priori fix the value of those decision variables whose value in ( x, z) and (x TLR , z TLR ) differs of at most ρ : 0 < ρ < 1, as specified by the rules that follow: The resulting problem is then passed to CPLEX, to which a solution time limit is imposed.

Computational Tests
For testing the performance of the new matheuristic, we considered 25 instances corresponding with regional DVB-T networks from Italy.Every instance corresponds to a single frequency network that is based on terrestrial DVB technology and is made up of a number of stations with a power emission belonging to the range [−40, 26] dBkW.The number of stations involved in the instances varies from 127 to more than 500 and the number of testpoints vary from about 2000 to more than 7000.We conducted the computational tests on a notebook equipped with the Windows operating system, a 2.70 GHz Intel i7 and 8 GB of RAM.In order to solve the MILP and ILP problems, we used the optimization software IBM ILOG CPLEX 12.5.We implemented the optimization model and solution algorithms by means of C/C++, using IBM Concert Technology for interfacing the code with CPLEX.
A time limit of 3600 s was set for CPLEX when solving the instances.Also the new matheuristic runs with a time limit of 3600 s.Such computational time budget is distributed in the following way among the two main phases: the GA phase is subject to a time limit of 3000 s, while the ILP-based improvement heuristic exploiting mod-RINS is subject to a time limit of 600 s.The parameters of the algorithm were set as follows, on the basis of preliminary computational experience: (1) In every tournament selection we use k = 10 groups; (2) each group includes a fraction α = 0.1 of the population POP; (3) the best m = 10 individuals of every group are chosen for crossover; (4) after having generated the new individuals, a fraction γ = 0.2 of the population experiences mutation; (5) in mod-RINS, we set ρ = 0.1.
In Table 2, we report the results of the computational tests.ID is the identifier of the instances.COV-CPLEX% is the percentage of population covered by the best solution found by the optimizer CPLEX within the time limit.COV-mathGA% and COV-mathGA+ILP% are instead the percentages of population covered by the matheuristic before and after the execution of the mod-RINS improvement heuristic.Finally, ∆COV-mathGA% and ∆COV-mathGA+ILP% denote the percentage coverage increase that the matheuristic grant with respect to CPLEX before and after mod-RINS.
Considering the values reported in Table 2, it is evident that even an effective and efficient solver like CPLEX finds difficulties in getting feasible solutions of high quality for the DVB model and that for the instances of larger size (with ID from I16 to I25) the percentage coverage is quite low.In contrast, the new matheuristic is able to identify solutions of much higher quality for all instances, both before and after the application of the mod-RINS improvement heuristic.Specifically, before the execution of mod-RINS, the matheuristic grants an average percentage improvement that is about 25% and that reaches almost 40% in the best case.After the execution of the improvement phase, the average percentage increase raises to 48% with peaks that exceed 60% improvement with respect to CPLEX.
As discussed in the subsection devoted to strengthening the model, this can be explained by the fact that CPLEX must face the combined effect of complicating signal-to-interference constraints with a large number of stations involved, ill-conditioned matrices and weak big-M formulations, which greatly slow the convergence to optimal solutions, whereas the matheuristic can much faster and more effectively explore the discrete space of power emissions.The improvements that we are able to obtain are very remarkable, especially considering that, when dealing with instances that involve large portions of territories like our regional instances, even small increases in coverage may entail that thousands of additional users are able to obtain the service.

Conclusions and Future Work
Optimally designing DVB-T networks by means of Mathematical Programming is a very challenging task and resulting models may prove very tough even for state-of-the-art optimization commercial software.With the aim of tackling the unsatisfying performance of commercial solvers, we presented a new matheuristic for DVB network design that combines mathematical optimization with exact and heuristic solution strategies.Specifically, we combined a genetic algorithm with a linear relaxation-based variable fixing strategy and an exact large neighborhood search formulated as an Integer Linear Programming problem.The obtained matheuristic refers to a Power-Indexed formulation, a strong and numerically stable pure Binary Linear Programming model for a DVB network design that is obtained through discretization of the power emissions of the DVB stations from a Mixed Integer Linear Programming model that explicitly include signal-to-interference ratios.Computational tests on realistic instances indicate that our matheuristic is able to identify solutions of much higher quality than those identified by a state-of-the-art solver, thanks to a more efficient exploration of the power emission solution space.
As future work, we plan to further increase the performance of the matheuristic, in particular through integration of it within a branch-and-cut approach.To this end, we would plan also to consider more refined lifting techniques of the considered GUB cover inequalities, to obtain stronger valid inequalities that could improve the upper bounds available on the problem.Moreover, we would aim also at improving the lower bounds, by studying refinements of the matheuristic, in particular by investigating alternative integrations of the GA-phase with suitable exact ILP searches.Last but not least, we also plan to tackle the uncertainty of signal propagation by means of Multiband Robust Optimization [63,64], a robust model that allows the ability to easily represent empirical discrete distributions commonly built in real-world applications on the basis of historical data.

Table 1 .
System elements and notation.
ts binary decision variable representing whether TP t is served by station s