Dual-Based Bounds for Resource Allocation in Zero-forcing Beamforming OFDMA-SDMA Systems

We consider multi-antenna base stations using orthogonal frequency division multiple access and space division multiple access techniques to serve single-antenna users. Some users, called real-time users, have minimum rate requirements and must be served in the current time slot while others, called non real-time users, do not have strict timing constraints and are served on a best-effort basis. The resource allocation problem is to find the assignment of users to subcarriers and the transmit beamforming vectors that maximize the total user rates subject to power and minimum rate constraints. In general, this is a nonlinear and non-convex program and the zero-forcing technique used here makes it integer as well, exact optimal solutions cannot be computed in reasonable time for realistic cases. For this reason, we present a technique to compute both upper and lower bounds and show that these are quite close for some realistic cases. First, we formulate the dual problem whose optimum provides an upper bound to all feasible solutions. We then use a simple method to get a primal-feasible point starting from the dual optimal solution, which is a lower bound on the primal optimal solution. Numerical results for several cases show that the two bounds are close so that the dual method can be used to benchmark any heuristic used to solve this problem. As an example, we provide numerical results showing the performance gap of the well-known weight adjustment method and show that there is considerable room for improvement.


Introduction
With the ubiquitous use of smart phones, tablets, laptops and other devices, traffic demand on wireless access networks is increasing geometrically. Multi-antenna base stations using orthogonal frequency-division multiple access (OFDMA) and space division multiple access (SDMA) can transmit at the same time to different sets of users on multiple subcarriers. This diversity increases the system throughput by assigning transmitting resources to users with good channel conditions. High data rates are thus made possible by exploiting the degrees of freedom of the system in time, frequency and space dimensions. OFDMA-SDMA is also supported by WiMAX and LTE-Advanced systems, which are the technologies currently being set up to implement the fourth generation (4G) cellular networks [1,2].
Due to these increased degrees of freedom, it is critical to use a dynamic and efficient scheduling and resource allocation (RA) mechanism that takes full advantage of all OFDMA-SDMA transmitting resources [3]. In this mechanism, the scheduler selects the users that are served at each frame and the RA algorithm allocates to these users the transmission resources required to meet the quality of service (QoS) requested from the upper layers.
In this paper, we focus on the RA part of the problem for an OFDMA-SDMA system supporting real time traffic with minimum rate requirements. We propose a dual-based method to get an upper bound to all feasible solutions. Because the dual solution is not necessarily primal feasible, we also present an algorithm to get a feasible point from the dual solution, which is a lower bound to the optimal primal solution. The importance of these bounds is that they give us limits on the duality gap and we can use them to estimate how far the solution given by any heuristic method is from the optimum.

State of the Art
The general resource allocation is a nonlinear, non-convex program so that it is almost impossible to solve it directly for any realistic number of subcarriers, users and antennas. For this reason, most research work focuses on developing heuristic algorithms. In this context, an important question is always how accurate are the results. For the resource allocation problem with rate constraints, it turns out that there are very few results of that kind, as we shall see. Traffic in the system can be divided into two main groups: delay-sensitive real time (RT) services and delay-insensitive non real time (nRT) services. Early work on OFDMA systems focused on solving the RA problem for nRT services only, where the objective is to maximize the total throughput with only power constraints and possibly minimum BER constraints. In [4], the complete optimization problem is divided into two sub-problems: Selection of users for each carrier and then power allocation to these users, which are both solved by a heuristic. A similar approach using zero-forcing (ZF) beamforming is reported in [5]. The work of [4,5] does not solve the complete optimization problem; instead, it separates it into uncoupled subproblems that provide suboptimal solutions.
There is a definite need to benchmark the performance of these heuristic algorithms. For the RA problem with nRT traffic only, several methods have been proposed to compute near-optimal solutions. For example, genetic algorithms are proposed in [6] while [7,8,9,10] provide methods to compute a near-optimal solution based on dual decomposition. In addition to providing a benchmark, near-optimal algorithms can also lead to the design of efficient RA methods as shown in [10], where heuristic algorithms derived from the dual decomposition method are proposed.
Several methods have been used to solve the RA problem for OFDMA-SDMA systems with both RT and nRT traffic. In [11], the objective is to maximize the sum of the user rates subject to per-user minimum rate constraints that model the priority assigned to each user at each frame. The optimization problem is solved approximately for each frame by minimizing a cost function representing the increase in power needed when increasing the number of users or the modulation order. The advantages of this approach are that it handles user scheduling and resource allocation together and supports Rt and nRT traffic. Its weaknesses are that no comparison is made against a near-optimal solution and the method used to determine user priorities at every frame is very complex.
In [12], both RT and nRT traffic are supported. Priorities are set according to the remaining deadline time for RT users and to the difference between the achieved rate and the desired rate required for nRT users. The user with the highest priority is paired with the subchannel with the highest vector norm and semi-orthogonal users are multiplexed on the same subchannel. Comparisons against the algorithm of [11] show that the packet drop rate for RT users is significantly reduced. The algorithm's complexity is also reduced because of the semi-orthogonality criteria used to add users. However, as in [11], a performance comparison with a near-optimal solution is not provided.
In [13], the objective is to maximize a utility function without any hard minimum rate constraints for the RT users. The channel quality information is added to the utility function to favor users with good channel conditions and priorities are set by increasing user weights in the utility function. The advantage is that the per-frame optimization problem has only a power constraint and no rate constraints, which makes it simpler to solve. The disadvantage with this reactive approach is that RT users with poor channel conditions are backlogged until their delay is close to the deadline, increasing the average delay and jitter.
In [14], a heuristic algorithm is proposed for the sum rate maximization problem with proportional rates among the user data rates, i.e., the ratio among allocated user rates is predetermined. The criteria to form user groups includes semi-orthogonality as in [12], but also fairness through proportional rate constraints. This method is extended to include hard minimum rates in [15] which attempts to solve exactly the same problem we formulate in section 2. Again, there is not reported method to evaluate the accuracy of these heuristics, except by comparing them to each other.
Another approach is to maximize the sum rate subject to constraints on the average rate delivered to a user [16]. However, unlike the work presented in [17] where an optimal solution is provided for the single antenna OFDMA RA problem with average rate constraints, the algorithm presented in [16] is a heuristic approximation. Note also that with average rate constraints, RT users tend to be served when they have good channel conditions which can create unwanted delay violations and jitter.

Paper Contribution and Organization
None of the previous work provide a near-optimal solution to the OFDMA-SDMA RA problem with minimum rate requirements. This is important not only as a benchmark for any heuristic algorithms, but also to get a better insight into the problem and to help devise efficient heuristics.
The main contribution of this paper is a method that provides an upper bound to the following OFDMA-SDMA RA problem for mixed RT and nRT traffic: For a given time slot, find the user selection and beamforming vectors that maximize a linear function of the users rates, given a total transmit power constraint and minimum rate constraints for RT users. The user weights in the linear utility function are arbitrary and can be the result of a prioritization or fairness policy by the scheduler. We solve the RA problem by Lagrange decomposition and we show, for small cases where we can find the optimal solution, that the duality gap is small.
A second contribution is a simple off-line heuristic algorithm to compute a feasible point based on the dual solution. This point is a lower bound for the optimal solution. This lower bound, in conjunction with the upper bound from the dual, can be used to limit the optimality gap in larger cases where an optimal solution is not available.
We then study several cases where we compare the performance of the upper and lower bounds. The results show that the two bounds are tight when the number of RT users is small and that their difference increases for larger values but that it stays quite small. Thus, the dual method provides a good approximation to the optimal solution. We also compare the performance of the weight adjustment algorithm versus the solution provided by the two bounds. The results indicate that adjusting the user weights to prioritize RT users can lead to significantly sub-optimal solutions. We describe the system and formulate the optimization problem we want to solve in Section 2 where we also briefly discuss a direct enumeration method to find the optimal solution for small problems. We present in Section 3 the dual method and two algorithms: One that finds the dual solution (the upper bound) and the other that finds a dual-based primal feasible solution (the lower bound). In Section 4, we present numerical results showing the accuracy of the upper and lower bounds and of the weight adjustment algorithm for different scenarios. Finally, we present our conclusions in Section 5 .

System Description and Problem Formulation
We consider the resource allocation problem for the downlink transmission in a MISO system with a single base station. There are K users, some of which have RT traffic with minimum rate requirements while the others have nRT traffic that can be served on a best-effort basis. The BS is equipped with M transmit antennas and each user has one receive antenna. In this configuration, the BS can transmit in the downlink data to different users on each subcarrier by performing linear transmit beamforming precoding. At each OFDM symbol, the BS can change the beamforming vector for each user on each subcarrier to maximize some performance function. In this paper, we assume that we use a channel coding that reaches the channel capacity. The data rate are in units of bits per OFDM symbol, or equivalently bits per second per Hertz (bps/Hz).

Signal Model
First, we describe the model used to compute the bit rate received by each user. Definẽ s k,n the symbol transmitted by the BS to user k on subcarrier n. We assume that thes k,n are independently identically distributed (i.i.d) random variables withs k,n ∼ CN (0, 1). w k,n an M -component column vector representing the beamforming vectors for user k on subcarrier n. Unless otherwise noted, we denote w the vector made up by the column stacking of the vectors w k,n . x n an M -component column vector representing the signal sent by the array of M antennas at the BS for each subcarrier n. h k,n an M -component row vector representing the channel between the M antennas at the BS and the receive antenna at user k for each subcarrier n. z k,n the additive white gaussian noise at the receiver for user k on subcarrier n. The z k,n are i.i.d. Gaussian random variables and, without loss of generality, we assume that they have unit variance, that is z k,n ∼ CN (0, 1). y k,n the signal received by user k on subcarrier n. r 0 k,n the rate of user k on subcarrier n in bps/Hz.
The signal vector x n is built by a linear precoding scheme which is a linear transformation of the information symbolss k,n : The signal received by user k on subcarrier n is then given by y k,n = h k,n x n + z k,n = h k,n w k,nsk,n + j =k h k,n w j,nsj,n + z k,n .
The second and third terms in the right-hand side of (2) correspond to the interference and noise terms. Since the signals and noise are Gaussian, their sum is also Gaussian and the data rate of user k for subcarrier n is given by the Shannon channel capacity for an additive white Gaussian noise channel:

Rate Maximization Problem
The general rate maximization problem corresponding to the OFDMA-SDMA RA problem with mixed RT and nRT traffic is to find a set of beamforming vectors w k,n that will maximize the weighted sum of user rates. This is limited by the total power available for the transmission at the base station and some users with real time QoS requirements must receive a minimum rate. More precisely, we assume that we know K Number of users in the cell. K Set of users in the cell: {1, . . . , K}. D A subset of K containing the users that have minimum rate requirements. We define D = |D|. d k Minimum rate requirement for user k. M Number of antennas at the BS. N Number of subcarriers available. P Total power available at the base station for transmitting over all channels. c k Weight of the user rates in the objective function. These could be computed by the scheduler to implement prioritization or fairness.
We then want to solve the following optimization problem to obtain the resource allocation The total transmit power is represented by the sum of the L 2 norms of the beamforming vectors in constraint (5). The achievable rate over all subcarriers should be higher or equal than the required minimum rate per user as in (6). Problem (4-6) is a non-convex, non-linear optimization problem. Using an exact algorithm to find a global optimal solution is very hard considering the size of a typical problem where there can be up to a hundred users and hundreds of sub-channels. Another option is to use a standard non-linear program (NLP) solver to compute a local optimal solution and use different starting points in the hope of finding a good global solution. The problem with this approach is that 1) we don't know how close we are to the true optimum and 2) the technique is quite time-consuming.
Nevertheless, we tried this approach for some small cases and observed that most users end up with a zero beamforming vector and only a small subset of users actually get some rate, often no more than M . Furthermore, in accordance to what was reported for the SDMA problem in [18], we observed that at high SNR, a so-called zero-forcing (ZF) solution is very close to the local optimum. This ZF solution can be easily computed by channel diagonalization and water-filling power allocation and is near-optimal compared to the general beamforming solution in the high SNR regime. Moreover, due to multi-user diversity, there is a good chance of finding users with high SNR channels when the total number of users increase. Therefore, in a multi-user Rayleigh fading environment, the ZF beamforming technique can provide results close to the general beamforming solution, even in the moderate SNR regime. For these reasons, we now turn to the so-called Zero-Forcing beamforming strategy.

Zero-Forcing Beamforming
In general, user k is subject to the interference from other users which reduces its bit rate, as indicated by the denominator in (3). Zero-forcing beamforming is a strategy that completely eliminates interference from other users. For each subcarrier n, we choose a set s of g ≤ M users which are allowed to transmit. This is called an SDMA set. We then impose the condition that for each user k in this set, the beamforming vector of user k must be orthogonal to the channel vectors of all the other users of the set. This amounts to adding the orthogonality constraints h k,n w j,n = 0 j = k, j, k ∈ s (7) and the user k data rate for subcarrier n then simplifies to: With ZF beamforming, the problem is now made up of two parts. We need to select a SDMA set s(n) for each subcarrier n and, for each selected SDMA set, we must compute the beamforming vectors in such a way that the total rate received by all users is maximized. Because of this, we need to add another set of decision variables α k,n a binary variable that is 1 if we allow user k to transmit on subcarrier n and zero otherwise. We denote the collection of α k,n by the vector α .
This results in the following ZF problem where A and B are some large positive constants. Constraint (12) guarantees that we do not choose more than M users for each subcarrier, constraint (13) guarantees that if we have chosen two users k and j, they meet the ZF constraints and is redundant for other users, and constraint (14) guarantees that the beamforming vector is zero for users that are not chosen. Problem (9)(10)(11)(12)(13)(14) is a non-linear mixed integer program (MIP) and these are known to be very hard to solve exactly. In this paper, whenever we need to get an exact solution, we use a complete enumeration of the binary variables α satisfying the constraints (12). For each α, we then find the beamforming vectors maximizing the utility function given the power, minimum rate and ZF constraints. The optimal solution is the vector α that maximizes the utility function. The problem of finding the optimal w variables for a particular α is a relatively small non-convex problem. Using the pseudo-inverse approach to satisfy the ZF constraint (see Section 3.3), it can also be transformed into an approximate convex problem and solved using a dual algorithm similar to the one described in Section 3.4.
It would seem that the zero-forcing model is not improving things much: We have gone from a non-convex nonlinear program to a non-convex mixed nonlinear program that can be solved by picking the global solution from a large collection of small convex problems. However, as we explain in section 3, this allows us to design an efficient and accurate algorithm.

Dual-Based Solution Method
Obviously, we cannot solve problem (9-15) fast enough to use it in real time. Not only is it NP-complete [4] but the actual computation time becomes quickly prohibitive for realistic sizes, even for off-line computations. Still, we need to compute solutions so that we can use them as benchmarks to evaluate the quality of real time heuristic approximations. We now present two off-line solution techniques that are tractable for problems of realistic size based on the Lagrange relaxation of the primal. Because problem (9-15) is not convex, there will often be a strictly positive duality gap at the solution of the dual. However, if it is small enough, we can use the solution provided by the dual method as a useful benchmark to check the accuracy of heuristic methods. Results in section 4 show that in many cases, the duality gap is in fact less than a few percent.
Solving the ZF problem will require some form of search over the α k,n variables. Note that this ranges over all subsets with a number of users smaller than or equal to M , so that the search space is going to be fairly large. Our first transformation is thus to separate the problem into single-subcarrier subproblems. For this, we dualize the constraints (10) and (11) since they are the ones that couple the subcarriers. Define the dual variables λ Lagrange multiplier for power constraint (10). µ k Lagrange multipliers for minimum rate constraint (11) of user k. The collection of µ k is denoted µ.
In order to simplify the derivation, we also define the dual variables µ k for all users k ∈ K. For users with no minimum rate requirements (k / ∈ D), we have µ k = 0 andď k = 0. In what follows, we use the standard form of Lagrangian duality which is expressed in terms of minimization with inequality constraints of the form ≤. Under these conditions, the multipliers λ, µ ≥ 0. We get the partial Lagrangian with constraints (12)(13)(14)(15). The value of the dual function Θ at some point (λ, µ) is obtained by minimizing the Lagrange function over the primal variables and the dual problem is which we can solve by the well known subgradient algorithm [19]. From now on, we concentrate on the calculation of the subproblem (17).

Subchannel Subproblem
Because of the relaxation of the carriers coupling constraints (10-11), the subproblems in (17) decouple by subcarrier since the objective (16) is separable in n and so are constraints (12)(13)(14). Computing the dual function then requires the solution of N independent subproblems. For each subcarrier n, this has the form min wn,αn where w n is the vector made up by the column stacking of the vectors w k,n for subcarriers n and α n denote the collection of α k,n for subcarrier n. Problem (20-23) is still a mixed NLP, albeit of a smaller size.

SDMA Subproblem
A simple solution procedure is to enumerate all possible choices for α k,n that meet constraint (12). This is called the extensive formulation of the problem. Each choice defines a SDMA set s and κ = |s| and here, we present a solution technique for the subproblems with given s. For each s, the problem separates into κ independent problems to compute the optimal beamforming vector w k,n,s for each user k ∈ s. For each user k ∈ s we know the set of channel vectors for the other members of s and we collect these vectors in the where f * k,n,s is given by the solution of the optimization H k,n,s w k,n,s = 0 where c ′ k = c k + µ k , w k,n,s is the beamforming vector for user k on subcarrier n for SDMA set s, and w n,s is the vector made up by the column stacking of the vectors w k,n,s for the κ users in s. Note that constraint (21) is automatically satisfied by the construction of s, constraint (23) simply drops out since w k,n,s = 0 for k ∈ s and constraint (22) remains only for k ∈ s, but we write it as (27) because we are considering only users that belong to SDMA set s. This is certainly not a feasible real-time algorithm, but for realistic values of K and M the number of SDMA sets is still manageable and the optimization sub-problem (24-27) is a small nonlinear program with M variables and κ − 1 linear constraints. It can thus be solved quickly by a number of techniques. Still, the overall computation load can be quite large. There will be κ such problems to solve for each SDMA set, and there are S = M i=1 K i such sets for each of the N subcarriers so that we have to solve the problem κ × S × N times, and this for each iteration of the subgradient algorithm. Clearly, any simplification of the beamforming subproblem can reduce the overall computation time significantly.

Approximate Solution to the Beamforming Problem
This can be done by the following construction. Instead of searching in the whole orthogonal subspace of H k,n,s as defined by (27), we pick a direction vector in that subspace and search only on its support. This will give a good approximation to the extent that the direction vector is close to the optimal vector. The choice of direction is motivated by the fact that the objective function depends only on the product h k w k,n,s . We then introduce a new independent variable q k,n,s = h k w k,n,s and because this variable is not independent of w k,n,s , we add (28) as a constraint. We then get the equivalent problem max w k,n,s ,q k,n,s h k w k,n,s = q k,n,s H k,n,s w k,n,s = 0.
Constraints (30) and (31) can then be rewritten in the standard form G k,n,s w k,n = b k,n,s where the G k,n,s matrix is the concatenation of h k and H k,n,s and b k,n,s = [q k,n,s , 0, 0 . . . 0] T . Since we are proposing to transform the constrained optimization over the κ variables into an unconstrained optimization over q k,n,s only, we must be able to express w k,n,s as a function of q k,n,s . The linear system being under-determined, this is obviously not unique. We propose to use G + k,n,s , the pseudo-inverse of G k,n,s , for the back-transformation w k,n,s = G + k,n,s b k,n,s . The pseudo-inverse picks the vector of minimum norm compatible with the linear system. In other words, choosing this transformation will minimize w k,n,s so that it is minimizing the power term in the objective function in (29). Because λ ≥ 0, this has the effect of contributing to the maximization of f * k,n,s . Note that this technique provides only an approximate solution of the beamforming problem; we cannot invoke Theorem 1 from [20] which shows that in certain cases, the pseudo-inverse transformation is optimal. A strong assumption for the theorem is that the objective function depends only on the q k,n,s variable, which is not the case here since (29) also depends on w k,n,s 2 . However, we observed from numerical results that the difference between the pseudo-inverse solution and the optimal solution is small. With this approximation we fix the direction of the beamforming vectors to w k,n,s =G + k,n,s b k,n,s =q k,n,s [G + k,n,s ] 1 where [G + k,n,s ] 1 denotes the first column of G + k,n,s . Now, we can obtain a problem formulation in terms of the user powers only by replacing the following expression in (29): w k,n,s 2 = γ 2 k,n,s p k,n,s , where γ k,n,s = [G + k,n,s ] 1 and p k,n,s = |q k,n,s | 2 . Adding the constraint p k,n,s ≥ 0, we get the equivalent problem max p k,n,s c ′ k log(1 + p k,n,s ) − λγ 2 k,n,s p k,n,s so that the computation time is basically that for computing G + k,n,s . Also note that using G + k,n,s we can also find the optimal beamforming vectors for all users in s, the only difference being that γ k,n,s is computed using the column of G + k,n,s corresponding to the channel vector of this user.

Solving the Dual Problem
To summarize, the dual function Θ(λ, µ) is obtained for the current values of the multipliers by finding for each subcarrier n = 1, . . . , N the optimal SDMA set s * (n) to the minimization problem in (24), where and p k,n,s is given by (35). Substituting back in (17), the dual function is with f n,s given by (36). For any value of the dual variables (λ, µ) we can determine the optimal primal variables (α, w); α is obtained by the optimal subcarrier assignment vector s(n) after performing the minimization over s in (37), and the optimal beamforming vectors w * n,k for the users k ∈ s * (n) are given by w * n,k = G + k,n,s * (n) [p The largest part of the computation to evaluate the dual function is the calculation of G + k,n,s which has to be done for each subchannel and each possible SDMA set. The number of evaluations can become quite large but the size of each matrix is relatively small so that the calculation remains feasible for medium-size networks. Furthermore, while solving the dual problem requires multiple subgradient iterations, the calculation of the pseudo-inverses is independent of the value of the multipliers. This means that the calculation of G + k,n,s can be done only once in the initialization step of the subgradient procedure.
Note that if the algorithm finds an optimal solution, the corresponding primal computed from the optimization of the Lagrange function will be feasible since the subgradients g (k) µ and g λ are small . We use the (somewhat inappropriate) expression "dual feasible" to denote such a solution. If, on the other hand, the algorithm stops before reaching the optimum, the primal will generally not be feasible.

Algorithm 1 Calculation of the Dual Problem
Construct the set S of all subsets of users of size 1 ≤ κ ≤ M for all n = 1 . . . N do for all s ∈ S do for all k ∈ s do Compute the pseudo-inverse G + k,n,s and γ k,n,s end for end for end for Choose an initial value λ 0 and µ 0 Subgradient iterations. We set a limit of Im on the number of iterations for all i = 1 . . . Im do Solve the N subproblems (24) to compute the dual function Compute the subgradients: g (k) µ =ď k − n r k,n g λ = n k∈s * (n) w * k,n 2 −P if g µ ≤ ǫ and g λ ≤ ǫ then Exit. A dual optimal solution has been found. else Update the multipliers Exit. A dual optimal solution was NOT found.
Algorithm 1 finds the optimal dual variables (λ * , µ * ) that solve the dual problem (18) using the subgradient method [19] with a fixed step size δ and therefore yields an upper bound to problem (9)(10)(11)(12)(13)(14)(15). In algorithm 1, the pseudoinverse matrices G + k,n,s can be computed by any number of well known algorithms. Here, we have used the Matlab command pinv [21]. Note that this algorithm can be used to solve the beamforming problem or equivalently the power allocation problem for a fixed SDMA set assignment. The only difference is that (24) becomes trivial since there is only one possible SDMA set per subcarrier.
It is not possible to evaluate the overall complexity of algorithm 1 because we don't have a bound on the number of dual iterations. We can nevertheless give a bound on the complexity of one iteration since it is dominated by the pseudo-inverse matrices computation. We have to compute N S ∼ N K M matrix inversions each one with complexity O(M 3 ), giving a total computational complexity O(N K M M 3 ). Computing the subgradient vectors and updating the dual variables have much lower complexity than the pseudo-inverse matrices computation.

Analysis of the Dual Function
Lets denote Θ * the maximum of the dual function Θ(λ, µ) over (λ, µ) ≥ 0. If U 1 is the weighted sum-rate objective function achieved by any feasible The value −Θ * , or any feasible approximation −Θ(λ, µ) to it, is thus an upper bound to the optimum value of the primal problem U * . Thus, −Θ * can be used to benchmark any other solution method. Figure 1 shows a contour plot of the dual function Θ(λ, µ) for the case of one RT user. The diamond marker shows the maximum.
We can get some insight on the shape of the dual function from figure 2 where we plotted the function with respect to µ for a fixed value of λ. The solid line curve corresponds to the same dual function as in figure 1, where the rate constraint is active,ď 1 = 50 bps/Hz. We see that the dual function goes through a maximum at µ = 0.24. We have also shown the case where we increase the minimum rate constraint so much that the problem becomes infeasible, e.g., we makeď 1 = 100 bps/Hz. As expected from duality theory, the dual function has no maximum since lim µ→∞ Θ(λ, µ) = ∞ as shown by the dash-dotted curve. Finally, the dashed curve at the bottom corresponds toď 1 = 0 bps/Hz such that the constraint is inactive and the solution where the maximum occurs is located at µ 1 = 0.
For completeness we show in Figure 3 the dual function as a function of λ when the rate constraints are feasible. The dual function increases rapidly and reaches a maximum at λ = 1.83.

Dual-Based Primal Feasible Method
The SDMA set selection and beamforming vectors found by Algorithm 1 do not always provide a primal feasible solution. The rate or power constraints might be violated whenever the algorithm stops because the number of iterations has been reached before the convergence rule is met. In this subsection, we propose a simple procedure to obtain a feasible point to problem (9-15) from the dual solution found with Algorithm 1. This point is not optimal but because we start from the dual optimal solution, we expect that it will be close to the optimal solution. Obviously this will give us a lower bound to the optimal primal solution.
Algorithm 2 summarizes this method. The algorithm begins by solving the dual problem (18) using Algorithm 1. If the solution is not feasible either directly or by recomputing the power allocation for the SDMA set assignment found in the dual problem, the algorithm performs a search by increasing the dual variables associated to the users whose QoS constraints are not met until a new SDMA set assignment is found. It then solves the power allocation problem for this new SDMA set assignment and checks if the minimum rate constraints are met. The search for new SDMA sets continues using this method until a feasible SDMA set assignment is found or a maximum number of iteration is reached. In contrast to the enumeration method described in Section 2.3 which performs an enumeration of all possible SDMA set assignments, the dual-based Algorithm 2 is a method that finds new candidate SDMA set assignments close to the dual optimal and then uses them to solve a simple power allocation problem until the rate and power constraints are met. This makes the search for a near-optimal feasible point much faster than finding the exact solution.

Algorithm 2 Calculating a Feasible Point from the Dual Solution
Solve the dual problem (18) using algorithm 1. This yields the optimal dual variables λ * , µ * k and a SDMA set assignment vector s * (n) for each subchannel n. Set s o 0 (n) = s * (n) Evaluate total power and user rate constraints (10)(11) if All constraints are met then Exit. A feasible solution has been found. end if Compute power allocation problem for s o 0 (n) and evaluate total power and user rate constraints (10)(11) if All constraints are met then Exit. A feasible solution has been found. end if Compute the multipliers µ k for users k such that r k <ď k for j = 1 to J do µ k = µ k + δ Find s o j = arg mins{fn,s} where fn,s is given by (25) for the current dual variables λ, µ Let s o j (n) be the SDMA assignment found if s o j (n) = s o j−1 (n) then We have found a new SDMA assignment Compute power allocation problem for s o j (n) and evaluate total power and user rate constraints (10)(11) if All constraints are met then Exit. A feasible solution has been found. end if end if end for Exit. A feasible solution was not found.

Performance Analysis
In this section, we present some numerical results on the performance of the dual-based algorithm and the accuracy of the upper and lower bounds. To show how they can be used to evaluate heuristic algorithms, we also compare them with the solution provided by a weight adjustment method described in Section 4.2.

Convergence of the Dual Algorithm
First we present in Figure 4 the value of the dual function and Lagrange multipliers as a function of the number of iterations for a given channel realization. The corresponding transmit power and the rate received by the RT user are shown in Figure 5. The parameters used for the calculation are listed in the figure titles. We see that the algorithm converges very quickly to a solution that is both close to the minimum value and feasible. This is typical of several other configurations, except that the number of iteration increases with the number of RT users.

Weight Adjustment Heuristic
As discussed in Section 1.1, several RA algorithms provide support for users with RT traffic by increasing the user weights in the utility function until they receive enough transmission resources. In this section, we describe a generic weight adjustment method which will be used to show that this technique leaves much room for improvement.

Algorithm 3 Weight Adjustment Algorithm
Solve RA problem (9-15) without minimum rate constraints constraints (11) c ′ ← c Let r k be the achieved rate for user k at every iteration iteration ← 1 while (r k <ď k for one or more users k ∈ D) AND (iteration ≤Ī) do Increase user weight using c ′ k = c ′ k + ǫ ď k − r k for users in need, where 0 < ǫ ≤ 1 Solve RA problem (9-15) without minimum rate constraints (11) using user weights c ′ iteration ← iteration +1 end while We want to find a set of weights in the utility function (9) such that the rate requirements of the RT users are met when we solve problem (9-15) without the rate constraints (11). Also, the set of weights must not be very different from one user to the other in order to maximize the multiuser diversity gain. Algorithm 3 implements a generic method that can do this. It increases the user weights for RT users until enough resources are allocated to meet the minimum rate requirements. The parameter ǫ controls how much the weights are increased with respect to the rate bounds. The rates achieved by Algorithms 3 and 2 are different since they solve different problems. Algorithm 3 can be seen as solving problem (9-15) by using a linear penalty method for constraints (11) of the form P k = min 0, r k −ď k .
The modified objective function is then At each iteration of the penalty method, whenever rate constraints are active, the solution of (40) cannot be smaller than that of (9-15) since it is a relaxation. Notice that problem (40) is quite simple since it has a single constraint for the transmit power but it has to be solved many times to adjust the weights of the real time users. In weight adjustment algorithms such as [13], the user weights are increased at each time slot using an increasing function of the packets delay, so the computation task is distributed over time. However, this distributed approach does not guarantee that the rate requirements are met in a given time slot which can lead to delay violations and jitter.

Parameter Setup and Methodology
We now present the method and parameter values used to compare the performance of the different methods to solve problem (9). We used a Rayleigh fading model to generate the user channels such that each component of the channel vectors h k,n are i.i.d. random variables distributed as CN (0, 1). We also assumed independent fading between users, antennas and subcarriers. Unless otherwise noted, we used a configuration with M = 3 antennas, K = 16 users and N = 16 subcarriers. We have only one RT user when we examine the effect of various parameters and we also look at the impact of increasing the number of RT users. The minimum rate constraint was set at 40 bps/Hz unless otherwise stated. We also fixed the power constraint toP = 20 dBm and used a large-scale attenuation of 0 dB for all users. The user weights in (9) were set to c k = 1 for all users. The results are the average over the feasible cases from 100 independent channel realizations. We compared the performance of the different methods for various scenarios where we increased the resource requirements for the RT users until the minimum rate requirements can no longer be met for all RT users. For each scenario and channel realization, the upper bound was computed from the dual solution using Algorithm 1 described in Section 3.4. For small systems, we also found the exact solution using the primal enumeration method given in Section 2.3. We also computed the lower bound given by dual-based primal feasible Algorithm 2 and the heuristic solution provided by the weight adjustment Algorithm 3 described in Section 3.6 and 4.2. We use the upper bound given by the dual optimal solution as the reference point when computing the gap when the exact solution is not available.

Single User, Increasing Minimum Rate
In this first scenario, we have a single RT user and we increase its minimum rateď 1 . First we consider a small system with K = 4 users and N = 2  Fig. 6 Average total rate as a function of the minimum rate requirements subcarriers where it is possible to compute an exact solution. We present in table 1 the average gap in percent between the three methods used to find feasible solutions against the dual upper bound for a small system. As the required minimum rate increases from 13.33 to 20 bps/Hz, the upper bound decreases as more resources need to be assigned to the RT user until the problem is no longer feasible. For this small configuration, we see that all methods give excellent results and the duality gap is very small.
In the remaining results, we use a larger system with K = 16 users and N = 16 subcarriers. With these values, it is no longer possible to use the primal enumeration method and we compute the gap relative to the dual upper bound. We present in Table 2 the difference in percentage between the bound and the solutions of the dual feasible and the weight adjustment algorithms. The dual feasible algorithm gives a lower bound within 0.25% of the dual upper bound while the weight adjustment solutions difference can be almost 10%. As discussed in Section 4.2, this is due to the fact that the weight adjustment algorithm stops as soon as it finds a feasible solution and does not have the option of finding a better assignment. As a result, the objective does not change much when the minimum rate is increased. This can be seen from Figure 6 which shows the sum rate achieved by the dual feasible algorithm and the weight modification method against the minimum rate requirement.

Single User, Increasing Attenuation
As the user moves away from the BS and the channel attenuation increases, the RA algorithm dedicates more resources to the RT user until the problem is unfeasible. Figure 7 shows the average total rate when the large-scale channel attenuation of the RT user varies from 0 to 15 dB. The results show that for all SNR, the dual lower bound provides a tight solution with the upper bound while the weight adjustment method shows a large performance gap. Table 3 shows the error in percentage between the objective and the upper bound. For an attenuation of 15 dB, neither method is able to find a feasible solution; the problem is feasible because the dual upper bound is around 140, but the algorithms cannot find a solution. That is, algorithm 1 for the dual upper bound finds a solution as long as the primal problem is feasible. On the other hand, algorithm 2 cannot find a solution when the solution point is close to the infeasibility region. In that case, the algorithm exits declaring that a feasible solution could not be found. Also, note that the proposed dual-based algorithms provides a near-optimal to the ZF problem (9)(10)(11)(12)(13)(14) for all SNRs where a feasible solution can be found by the algorithm. Figure 7 shows average sum rates. The average is performed over the channel realizations that provide feasible points. When the attenuation of the RT user is low, most of the channel realizations produce feasible points. When it is high, some of the channel realizations do not produce feasible points and are discarded. For example, in figure 7 when the attenuation is 15 dBs, none or very few of the channel realizations produce feasible points. Therefore, figure  7 does not show a sum rate for that point.

Increasing Number of RT Users
Finally, figure 8 shows the upper dual bound, the lower bound and the solution given by weight adjustment methods as a function of the number of RT users.  For a single RT user, we have seen in tables 2 and 3 that the difference between the upper and lower solution is small. In figure 8, we see that this difference increases for three or more RT users. Still, this growth is not large and we can consider that a 3.52% is an acceptable error tolerance. Based on this, we can claim that it is possible to find a quasi-optimal solution to problem (9-15) with the proposed method, albeit with an off-line algorithm.
Furthermore, the results show that the weight modification method has a large performance gap which becomes more significant as the number of RT users increase. Also, the dual method can find feasible solutions for cases where the weight adjustment method cannot. This shows that the weight modification method should be used carefully for RA in OFDMA-SDMA systems with RT users and that more efficient heuristics should be developed to approach the performance of the dual-based feasible solution.

Conclusion
In this work, we proposed a method to compute the beamforming vectors and the user selection in an OFDMA-SDMA MISO system with minimum rate requirements for some RT users. We gave a rigorous mathematical formulation of the Zero-forcing model and then used a Lagrangian relaxation of the power and rate constraints to solve the dual problem using a subgradient algorithm. The Lagrange decomposition yields sub-problems separated per subcarrier, SDMA sets and users which substantially reduces the computational complexity. We obtained a simple expression of the dual function for the beamforming problem for a given SDMA set based on a pseudo-inverse condition on the beamforming vectors. The dual optimum can then be used as a benchmark to compare against any other solution methods and heuristics. The dual function also gives us a better understanding of the problem. Its shape is related to the rate constraint activation and problem feasibility, and it also justifies the splitting of the subcarrier assignment and power allocation processes used in several heuristic methods. We then proposed an algorithm which finds a feasible point by starting from the dual-based optimum solution and searching among the dual variables of the rate constraints. Numerical results indicate that the two bounds are tight so that the feasible point is near-optimal. As a point of comparison, we also evaluated the performance of a weight adjustment method which uses weight adjustments in the objective function to achieve the required rates. Our results show that the performance gap of this approach is large and grows when the SNR of a single RT user increases or when the number of RT users increases.
In addition, the weight adjustment method requires many time slots to adjust the weights and schedule real time users. The dual-based method explicitly includes the minimum rate constraints which allows RT users to be scheduled in the current slot, which decreases the overall packet delay and jitter.
The significant gap between the weight adjustment algorithm and the optimal RA solution, suggests that there is a need to find better heuristics. The dual approach looks promising to guide the design of efficient novel heuristics. To implement the RA algorithm in real time, we also need to design fast methods to reduce the number of SDMA sets to be searched. The design of these heuristic algorithms is outside the scope of this paper and is part of our current efforts. Finally, the upper and lower bounds also provide a very useful benchmark to compare the performance of any heuristic method.

Competing interests
The authors declare that they have no competing interests.