On Improving the Capacity of Solving Large-scale Wireless Network Design Problems by Genetic Algorithms

TL;DR

This paper introduces a novel genetic algorithm designed to efficiently solve large-scale wireless network design problems, specifically focusing on power, frequency, and modulation assignment, demonstrating improved exploration and numerical stability.

Contribution

First genetic algorithm tailored for large-scale wireless network design problems, enhancing solution exploration and numerical stability over previous methods.

Findings

Effective in solving large realistic WiMAX network instances

Outperforms previous approaches in solution quality and stability

Enables wider exploration of power solutions

Abstract

Over the last decade, wireless networks have experienced an impressive growth and now play a main role in many telecommunications systems. As a consequence, scarce radio resources, such as frequencies, became congested and the need for effective and efficient assignment methods arose. In this work, we present a Genetic Algorithm for solving large instances of the Power, Frequency and Modulation Assignment Problem, arising in the design of wireless networks. To our best knowledge, this is the first Genetic Algorithm that is proposed for such problem. Compared to previous works, our approach allows a wider exploration of the set of power solutions, while eliminating sources of numerical problems. The performance of the algorithm is assessed by tests over a set of large realistic instances of a Fixed WiMAX Network.

Tables

Table 1. Table 1: Comparisons between (BM) and WPLAN formulations

							$\left|\text{T*}\right|$
ID	$\left|\text{T}\right|$	$\left|\text{B}\right|$	$\left|\text{F}\right|$	$\left|\text{H}\right|$		(BM-PFMAP)	WPLAN [7]	(GA-PFMAP)
S1	100	12	1	1		63 (78)	74	70
S2	169	12	1	1		99 (100)	107	103
S3	196	12	1	1		108	113	113
S4	225	12	1	1		93	111	115
S5	289	12	1	1		77	86	88
S6	361	12	1	1		154	170	175
S7	400	18	1	1		259 (266)	341	319
R1	400	18	3	4		370	400	400
R2	441	18	3	4		302 (303)	441	441
R3	484	27	3	4		99 (99)	427	434
R4	529	27	3	4		283 (286)	529	462
Q1	400	36	1	4		0	67	72
Q2	441	36	1	4		191	211	222
Q3	484	36	1	4		226	463	466
Q4	529	36	1	4		145 (147)	491	491

Equations

$$x_{tb}^{fh}$$

$$\frac{a_{t\beta }\cdot p_{\beta }^f}{N+{\sum }_{b\in B\setminus \{\beta \}}{a}_{tb}\cdot p_b^f}\ge {\delta }_h.$$

$$a_{t\beta }\cdot p_{\beta }^f-{\delta }_h\sum _{b\in B\setminus \{\beta \}}{a}_{tb}\cdot p_b^f\ge {\delta }_h\cdot N.$$

$$\underset{\forall (\beta ,f,h):\beta \in B,f\in F,h\in H}{\bigvee }\left(a_{t\beta }\cdot p_{\beta }^f-{\delta }_h\sum _{b\in B\setminus \{\beta \}}{a}_{tb}\cdot p_b^f\ge {\delta }_h\cdot N\right).$$

$$a_{t\beta }\cdot p_{\beta }^f-{\delta }_h\sum _{b\in B(t)\setminus \{\beta \}}{a}_{tb}\cdot p_b^f+M\cdot (1-x_{t\beta }^{fh})\ge {\delta }_h\cdot N.$$

$$\max$$

$$a_{t\beta }\cdot p_{\beta }^f-{\delta }_h\sum _{b\in B\setminus \{\beta \}}{a}_{tb}\cdot p_b^f+M\cdot (1-x_{t\beta }^{fh})\ge {\delta }_h\cdot N$$

$$t\in T,b\in B,f\in F,h\in H$$

$$\sum _{b\in B}\sum _{f\in F}\sum _{h\in H}{x}_{tb}^{fh}\le 1$$

$$\sum _{t\in T}\sum _{h\in H}{d}_t\cdot \frac{1}{s_h}\cdot x_{t\beta }^{fh}\le D$$

$$p_b^f\in [0,P_{max}]$$

$$x_{tb}^{fh}\in \{0,1\}$$

$$p{3}_{\beta }^{\varphi }$$

$$p{4}_{\beta }^{\varphi }$$

Full text

1 Introduction

During the last years, wireless communications have experienced an explosive growth thus rapidly leading to a dramatic congestion of radio resources. In such complex scenario, the trial-and-error approach commonly adopted by practitioners to design networks has clearly shown its limitations. Telecommunications companies and authorities are thus searching for more effective and efficient design approaches, also looking on Optimization (as shown by the recent call for tenders for developing a Digital Video Broadcasting simulator by the Italian Communications Regulatory Authority [2]). Many models and solution methods have been proposed for solving the problem of designing a wireless network. However, solving to optimality the overall problem is still a big challenge in the case of large instances. In this paper we present a Genetic Algorithm for solving the Power, Frequency and Modulation Assignment Problem, a relevant problem that arises in the design of wireless networks and captures specific features of Next Generation Networks.

2 The Wireless Network Design Problem

For modeling purposes, a wireless network can be described as a set of transmitters $B$ that provide for a telecommunication service to a set of receivers $T$. A receiver $t\in T$ is said to be covered or served if it receives the service within a minimum level of quality. Transmitters and receivers are characterized by a location and a number of radio-electrical parameters (e.g., power emission and frequency). The Wireless Network Design Problem (WND) consists in establishing the location and suitable values for the parameters of the transmitters with the goal of maximizing the number of covered receivers or a revenue associated with coverage.

In this work we consider a generalization of the so-called Power and Frequency Assignment Problem (PFAP), a version of the WND that is known to be NP-hard [12]. In addition to power emission and frequency, we also consider the transmission scheme (burst profile) as parameter to be established, modeling a feature of Next Generation Networks, such as WiMAX [3, 6]. We indicate this generalization of the PFAP by the name Power, Frequency and Modulation Assignment Problem (PFMAP).

In the PFMAP two decisions must be taken: (1) establishing the power emission of each transmitter on each available frequency and (2) assigning served receivers to activated transmitters specifying the frequency and the burst profile used to transmit. To model these decisions, as first step we introduce the set $F$ of available frequencies and the set $H$ of available burst profiles. Each frequency $f\in F$ has a constant bandwidth $D$ and each burst profile $h\in H$ is associated with a spectral efficiency $s_{h}$, which is the bandwidth required to satisfy one unit of traffic demand generated by a receiver.

Then we can introduce two typologies of decision variables, namely:

• •

a continuous power variable $p_{b}^{f}\in[0,P_{\max}]\hskip 5.69046pt\forall\hskip 2.84544ptb\in B,f\in F$ representing the power emission of a transmitter $b$ on frequency $f$;

• •

a binary service assignment variable $x_{tb}^{fh}\in\{0,1\}\hskip 5.69046pt\forall\hskip 2.84544ptt\in T,b\in B,f\in F,h\in H$ defined in the following way:

[TABLE]

Given a frequency $f\in F$, every receiver $t\in T$ picks out signals from every transmitter $b\in B$ and the power $P_{b}^{f}(t)$ that $t$ receives from $b$ on $f$ is proportional to the emitted power $p_{b}^{f}$ by a factor $a_{tb}\in[0,1]$, i.e. $P_{b}^{f}(t)=a_{tb}\cdot p_{b}^{f}$. The factor $a_{tb}$ is called fading coefficient and summarizes the reduction in power that a signal experiences while propagating from $b$ to $t$ [14].

Among the signals received from transmitters in $B$, a receiver $t$ can select a reference signal (or server), which is the one carrying the service. All the other signals are interfering. We remark that, since each transmitter in $B$ is associated with a unique signal, in what follows we will also refer to $B$ as the set of signals received by $t$.

A receiver $t$ is regarded as served by the network, specifically by server $\beta\in B$ on frequency $f\in F$ through burst profile $h\in H$, if the ratio of the serving power to the sum of the interfering powers (signal-to-interference ratio or SIR) is above a threshold $\delta_{h}$ (SIR threshold) whose value depends on the used burst profile $h$ [14]:

[TABLE]

Note that in the denominator we highlight the presence of the system noise $N>0$. By simple algebra operations, inequality (2) can be transformed into the following linear inequality, commonly called SIR inequality:

[TABLE]

As we do not know a priori which transmitter $b\in B$ will be the server of a receiver $t\in T$ and which frequency $f\in F$ and burst profile $h\in H$ will be used, given a receiver $t\in T$ we have one SIR inequality (3) for each potential server $\beta\in B$ and potentially usable frequency $f$ and burst profile $h$. To ensure that $t$ is covered, at least one of such inequalities must be satisfied. This requirement can be equivalently expressed through the following disjunctive constraint:

[TABLE]

Adopting a standard approach used in Mixed-Integer Programming (see [13]), the above disjunction can be represented by a family of linear constraints in the $p$ variables by introducing a large positive constant $M$, the so-called big-M coefficient. Specifically, given a receiver $t\in T$ we use the assignment variable $x_{t\beta}^{fh}$ to introduce the following constraint for each potential 3-ple $(\beta,f,h)$:

[TABLE]

It is easy to check that if $x_{t\beta}^{fh}=1$ then (5) reduces to the simple SIR constraint (3). If instead $x_{t\beta}^{fh}=0$ and $M$ is sufficiently large111For example, we can set $M=\delta_{h}\cdot N+\delta_{h}\sum_{b\in B\setminus\{\beta\}}a_{tb}\cdot P_{max}$ , then (5) is satisfied by any feasible power vector $p$ and becomes redundant.

By using constraints (5) and by introducing a parameter $r_{t}$ to denote revenue associated with receiver $t\in T$ (e.g., population, number of customers) , we can define the following natural formulation (BM-PFMAP) for the PFMAP [6, 7]:

[TABLE]

The objective function is to maximize the total revenue obtained by serving receivers and constraint (7) ensures that each receiver is served at most once. Each receiver generates a traffic demand $d_{t}$ and each frequency has a bandwidth equal to $D$. Constraint (8) ensures that the sum of traffic demands (re-sized by the spectral efficiency $s_{h}$ of the used burst profile) generated by the receivers served by a transmitter does not exceed the bandwidth of the frequency. Finally, (9) and (10) define the decision variables of the problem.

Drawbacks of the natural formulation. The natural formulation (BM-PFMAP) expands a basic model that is widely used for the WND in different application contexts, such as DVB, (e.g., [12]), UMTS (e.g., [1, 11]) and WiMAX ([6, 7]). In principle, such basic model and (BM-PFMAP) can be solved by commercial solvers such as IBM ILOG CPLEX [10]. However, it is well-known (see [7]) that: (i) the fading coefficients may vary in a wide range leading to (very) ill-conditioned coefficient matrices that make the solution process numerically unstable; (ii) the big-M coefficients generate poor quality bounds that dramatically reduce the effectiveness of standard solution approach [13]; (iii) the resulting coverage plans are often unreliable and may contain errors (e.g., [7, 11]). In practice, the basic model and (BM-PFMAP) can be solved to optimality only when used for small-sized instances. In the case of large real-life instances, even finding feasible solutions can represent a difficult task, also for state-of-the-art commercial solvers like CPLEX. Though these drawbacks are well-known, it is interesting to note that just a relatively small part of the wide literature devoted to WND has tried to overcome them. We refer the reader to [6] for a review of works that have tried to tackle these drawbacks.

2.1 Contribution of this work and review of related literature

In this paper, we develop our original contribution by starting from a recent work, [7], that proposes a family of strong valid inequalities for tackling the drawbacks of (BM-PFMAP) that we have pointed out. The idea at the basis of [7] is to quit modeling emission power as a continuous variable $p_{b}$ and to use instead a set of discrete power levels ${\cal P}=\{P_{1},\dots,P_{|{\cal P}|}\}$, with $P_{1}=0$ (switched-off value), $P_{|{\cal P}|}=P_{max}$ and $P_{i}>P_{i-1}$, for $i=2,\dots,|{\cal P}|$. This basic operation allows the authors to define a family of lifted GUB cover inequalities that are used in a solution algorithm that drastically enhances the quality of solutions found.

The solution algorithm proposed in [7] is motivated by a trade-off that arises from discretization: larger sets of discrete levels lead in principle to better solutions, but on the other hand the corresponding 0-1 Linear Program gets larger and harder to solve. The computational experience shows that very good solutions can be found by considering small sets with well-spaced power values, but that no improvement is obtained within the time limit when a number of levels higher than six is considered.

In the present work, we investigate the possibility of using a Genetic Algorithm (GA) [8] as a fast heuristic to widen the exploration of the discrete power solution space: the aim is to exploit the entire set of discrete power levels and thus to evaluate power configurations with levels not included in the best solutions found in [7]. In particular, our aim is to improve the capacity of solving large realistic instances by finding higher value solutions. We thus design a GA that takes into account the specific features of the PFMAP and we test its performance on the same set of realistic WiMAX instances used in [7].

Heuristics have been extensively used to tackle large instances of different versions of the WND problem. Two relevant cases are provided by [1], where a two-stage Tabu Search algorithm is proposed to solve the base station location and power assignment problem in UMTS networks, and by [12], where a GRASP algorithm is proposed to solve the PFAP arising in the planning of the Italian National DVB network. The use of GA to solve versions of the WND is not a novelty as well and many works can be found in literature. However, to our best knowledge, no GA has been yet developed to solve the PFMAP and the algorithm that we propose is the first for solving this level of generalization of the WND. Until now, GAs were indeed developed to solve just single aspects of the PFMAP: (i) the transmitter location problem (e.g., [4]); (ii) the service assigment problem (e.g., [9]); (iii) the frequency assignment problem (e.g., [5]); (iv) the power assignment problem (e.g., [15]). Moreover, we remark that our algorithm is the first to be designed with the specific aim of improving the capacity of solving instances, while tackling the numerical problems pointed out in Section 2. We now proceed to present our original contributions for the WND.

3 A Genetic Algorithm for the PFMAP

A Genetic Algorithm (GA) is a heuristic method for solving optimization problems that resembles the evolution process of a population of individuals (for a comprehensive introduction to the topic we refer the reader to [8]). At any iteration, a GA maintains a population whose individuals represent feasible solutions to the problem. The solution is encoded in a chromosome associated with each individual. The genetic strength of an individual is evaluated by a fitness function that establishes the quality of the corresponding solution to the problem. A GA starts by defining an initial population, that iteration after iteration changes by crossover, mutation and death of individuals, according to a natural selection Darwinistic mechanism.

We develop a GA for the PFMAP that presents the following general structure:

Creation of the initial population 2. 2.

UNTIL the arrest condition is not satisfied DO

1. (a)

Selection of individuals who generate the offspring 2. (b)

Generation of the offspring by crossover 3. (c)

Mutation of part of the population 4. (d)

Death of part of the population

We now characterize the elements and the phases presented above for the algorithm (GA-PFMAP) that we propose to solve the PFMAP.

3.1 Characterization of the population

Individual representation. As we have explained in Section 2.1, our aim is to conduct a wider exploration of the power solution space, trying to obtain solutions with higher value. To this end, we establish that the chromosome of an individual corresponds to a power vector $p$ of size $|B|\cdot|F|$. Specifically, the chromosome presents one locus for each transmitter $b\in B$ and frequency $f\in F$ and each locus stores the power $p_{b}^{f}$ emitted by $b$ on $f$, namely $p=(p_{1}^{1},\ldots,p_{1}^{|F|},p_{2}^{1},\ldots,p_{2}^{|F|},\ldots,p_{|B|}^{|F|})$. Such power belongs to the set of discrete power levels $\cal{P}$, i.e. $p_{b}^{f}\in$ $\cal{P}$ $=\{P_{1},\ldots,P_{|\cal{P}|}\}$.

We remark that establishing the power emission $p_{b}^{f}\hskip 5.69046pt\forall\hskip 2.84544ptb\in B,\hskip 2.84544ptf\in F$ does not completely characterize a solution of the PFMAP. We indeed have to fix the value of the assignment variables $x_{tb}^{fh}$ and thus we need to set some assignment rule. First, note that given a power vector $p=(p_{1}^{1},p_{1}^{2},\ldots,p_{|B|}^{|F|})$ and a receiver $t\in T$ we can compute the power $P_{b}^{f}(t)$ that $t$ receives from $b$ on $f$, $\hskip 2.84544pt\forall b\in B,f\in F$. Through $P_{b}^{f}(t)$, if we fix the server $\beta\in B$ of $t$, we can check if there exists a SIR inequality (3) that is satisfied for some frequency $f\in F$ and burst profile $h\in H$. We establish the following assignment rule: as server of $t$ we choose the transmitter $b$ that ensures the highest received power $P_{b}^{f}(t)$ on some $f$. This in fact ensures the highest serving power. Once that the server $\beta$ is chosen, we can identify the SIR inequalities (3) that are satisfied by $p$ when $t$ is served by $\beta$ for some $f\in F$ and $h\in H$. If the SIR inequality is satisfied for a multiplicity of frequencies and/or burst profiles, we first choose as serving frequency $\hat{f}$ the one that ensures the highest value for the left-hand-side of (3) and then we choose as burst profile $\hat{h}$ the one that ensures the highest spectral efficiency (see Section 2). Thus for $t$ served by $\beta$ we have $x_{t\beta}^{\hat{f}\hat{h}}=1$ and $x_{t\beta}^{fh}=0$ $\forall f\in F\setminus\{\hat{f}\},h\in H\setminus\{\hat{h}\}$.

Note that this last rule may assign a receiver $t$ to a transmitter $\beta$ that violates the capacity constraint (8) of $\beta$ on the chosen frequency $\hat{f}$. If this is the case, we choose the second best couple of frequency and burst profile according to the rule. If this not possible, the third best and so on. In the end, if there is no capacity left for any valid couple $(f,h)$, $t$ is not considered covered by $\beta$.

Fitness function. As the aim of the WND is to maximize coverage, we adopt a fitness function that evaluates the coverage ensured by an individual. Specifically, the fitness $COV(p)$ of an individual is equal to the number of receivers that are covered when the power vector is $p$ and service assignment is done according to the previously introduced rules.

Initial population. Our aim is to consider all the feasible discrete power levels from the beginning. Therefore, the initial population is represented by the power vectors that are obtained by activating a single transmitter $b\in B$ on a single frequency $f\in F$ at each of the discrete power levels $P_{l}\in\cal{P}$. For every $b\in B,f\in F$, the corresponding individuals included in the initial population are thus: $(0,0,\ldots,p_{b}^{f}=P_{2},\ldots,0,0)\hskip 8.5359pt\cdots\hskip 8.5359pt(0,0,\ldots,p_{b}^{f}=P_{|\cal{P}|},\ldots,0,0)\hskip 5.69046pt$. Note that we exclude the individual corresponding to all transmitters turned off, i.e. $p_{b}^{f}=P_{1}=0\hskip 2.84544pt\forall b\in B,f\in F$. We thus have $|B|\cdot|F|\cdot|L-1|$ initial individuals. We denote the set of individuals representing the population at a generic iteration of the algorithm by $P$.

3.2 Evolution of the population

Selection. In order to select the individuals who give birth to the new generation, we adopt a tournament selection approach: given the set $P$ of individuals constituting the current population and a value $0<\alpha<1$, we first define a number $k\in\mathbb{Z}^{+}$ of groups including $\lfloor\alpha\cdot|P|\rfloor$ individuals who are randomly extracted from $P$. Then we extract $m<\lfloor\alpha\cdot|P|\rfloor$ individuals with the best fitness from every group. These are the individuals who generate offspring by crossover.

Crossover, mutation and death. The individuals selected for crossover are randomly paired up to constitute $\lfloor k\cdot m/2\rfloor$ couples. Each couple generates two offspring by mixing its chromosome. Given a couple of parents with power vectors $p1,p2$, the crossover operation consists in mixing power levels that are in the same position of $p1$ and $p2$ to generate two offspring with (possibly) higher fitness power vectors $p3,p4$.

Before presenting the crossover procedure, we need to define a measure that evaluates the impact of crossover on coverage. To this end, let $\Delta\mbox{COV}(p,p_{b}^{f}=P_{l})\in\mathbb{Z}$ denote the variation in the number of covered receivers caused by changing the power value $p_{b}^{f}$ in position $(b,f)$ of vector $p$ to the value $P_{l}$, while maintaining all the other power values unchanged. We can then propose the following crossover procedure, that concentrates the effort of creating a higher fitness individual on $p3$. At the beginning of the crossover, $p3$ and $p4$ have all elements equal to 0. Then, by following this ordering of indices $(b,f):b\in,f\in F$: $(1,1)\mbox{ }(1,2)\mbox{ }\ldots\mbox{ }(1,|F|)\mbox{ }(2,1)\mbox{ }\ldots\mbox{ }(2,|F|)\mbox{ }\ldots\mbox{ }(|B|,1)\mbox{ }\ldots\mbox{ }(|B|,|F|)$, each null value inherits the power value in the same position of one of the two parents.

We now present the crossover rule for a generic position $(\beta,\phi)$. For indices $(b,f):b<\beta,f<\phi$, the crossover was already executed and thus the offspring vectors $p3,p4$ present power levels inherited by the parents $p1,p2$. Power levels of $p3,p4$ in positions $(b,f):b\geq\beta,f\geq\phi$ are instead still equal to zero. The rule to establish the power value inherited by $p3,p4$ in $(\beta,\phi)$ is the following:

[TABLE]

This ensures that, at any step of the crossover procedure, offspring $p3$ inherits the power level of the parent that allows the most favourable variation $\Delta COV$ in coverage.

In addition to crossover, we also allow to alter the power vector of single individuals by mutation. This introduces new genetic information in the population and helps to widen the solution space exploration and to avoid entrapment in local optima. At any iteration, a number of individuals $\lfloor\gamma\cdot|P|\rfloor$ with $0<\gamma<1$ is randomly chosen. Then, still by random selection, $|F|$ power levels corresponding with different frequencies are reduced to the immediately lower power level allowed in $\cal{P}$. This mutation rule is set with the aim of defining new power vectors that have lower powers but ensure the same coverage. The reduction in power is generally desirable as a signal that is useful for a receiver may be interfering for a different receiver.

Finally, after crossover and mutation, the weakest individuals die and are removed from $P$. Specifically, we choose to select and remove the $2\cdot\lfloor k\cdot m/2\rfloor$ individuals with the worst fitness function. The size of P is thus maintained constant over all the iterations.

4 Computational experience

We test the performance of our GA on a set of 15 realistic instances, developed with the Technical Strategy & Innovations Unit of British Telecom Italia (BT Italia SpA). All the instances refers to a Fixed WiMAX Network [3], deployable in an urban area corresponding to a residential neighborhood of Rome (Italy). The instances consider various scenarios with up to $|T|=529$ receivers, $|B|=36$ transmitters, $|F|=3$ frequencies, $|H|=4$ burst profiles (see Table 1). This leads to large formulations (BM-PFMAP) that are very hard to solve. For a detailed description of the instances, we refer the reader to [7].

For each instance, we run the proposed algorithm (GA-PFMAP) 50 times with a time limit of 1 hour by using a machine with a 1.80 GHz Intel Core 2 Duo processor and 2 GB of RAM. Each tournament selection involves $k=20$ groups that include a fraction $\alpha=0.05$ of the population $P$. The best $m=8$ individuals of each group are selected for crossover and, after the generation of the new individuals, mutation affects a fraction $\gamma=0.1$ of the population.

In Table 1, we compare the value of the best solution obtained through the three approaches that we consider, namely the direct solution of (BM-PFMAP) by ILOG Cplex 10.1, the solution of the Power-Indexed formulation by the algorithm WPLAN [7] and the solution of (BM-PFMAP) by the proposed algorithm (GA-PFMAP). Results for (BM-PFMAP) and WPLAN are derived from [7]. The presence of two values in some lines of the column of (BM-PFMAP) indicates that the coverage plans returned by Cplex contain errors and some receivers are actually not covered. We remark that (GA-PFMAP) provides solutions that always ensure a higher coverage than (BM-PFMAP) and without coverage errors. Making a comparison with WPLAN, we instead note that (GA-PFMAP), though in some cases finds solutions with lower coverage, for most of the cases is able to find solutions that ensure an equal or higher number of covered receivers than WPLAN. This is particularly evident for instances that seems to be very hard to solve through (BM-PFMAP). The algorithm is thus effective and is worth of further investigations.

5 Conclusion and future work

We presented a Genetic Algorithm (GA) to tackle large realistic instances of a relevant problem arising in wireless network design. We showed that a GA helps to improve the value of solutions found through a wider exploration of the power space. A future research path could be represented by the integration of a refined GA into an exact solution method. It is indeed common experience that the combination of fast heuristics with Mixed-Integer Linear Programming leads to a great reduction in the running times w.r.t pure exact optimization methods. A sensitivity analysis of the GA parameters and a study on the impact of different starting conditions and selection strategies would also constitute important subjects of investigations.

Acknowledgments.

The author thanks Carlo Mannino and Antonella Nardin for fruitful discussions. Thanks go also to the three anonymous referees who provided valuable comments.