How to solve orthogonal channel sum-rate maximization with CVXPY?

The optimization tools and numerical computation tools for Python have come a long way. These days you can formulate and solve complex problems with only a couple of lines of code.

In this post, I am going to show you how to formulate and solve the sum rate (throughput) maximization problem over inter-channel interference (ICI), i.e. orthogonal, channels with sum power constraint. That is, the channels are assumed to be in a technical sense orthogonal. This can mean OFDMA, TDMA, ZF-MIMO, or any other means of removing the ICI by design.

Let's first define what we mean by channel capacity in this context. We can write the throughput for \(N\) orthogonal channels by the sum of their individual rates. That is,

\[ \sum_{n=1}^N \text{log}_2 (1 + \text{SNR}_n) \]

To get to the capacity of the said system, we need to set some constraints on the SNR. We will use the sum power budget (\(P\)), which limits the available power across the channels as

\[ \sum_{n=1}^N p_n = P \]

Further, we must define the \(\text{SNR}_n\) in terms of the \(p_n\) as

\[ \text{SNR}(p_n) = {g_n p_n \over N_0} \]

where \(g_n\) denotes the channel gain for the \(n^\text{th}\) channel, \(p_n\) is the corresponding power allocation and \(N_0\) is a fixed noise-level determined by the given SNR.

Now, we can write the capacity as

\[ C = \sum_{n=1}^N \text{log}_2 (1 + \text{SNR}(p^*_n)) \]

where \(p^*_n\) is the sum-rate maximizing power allocation for channel \(n\).

Problem formulation

Next, we formulate the capacity-achieving power allocation as a throughput maximization problem, where we try to maximize the sum throughput over \(N\) channels with varying channel gains while being limited by a power budget. This problem can be mathematically formulated as

\[ \begin{array}{rl} \underset{p_n}{\text{max.}} & \displaystyle \sum_{n=1}^N \text{log}_2 (1 + {g_n p_n \over N_0}) \\ \text{s.t.} & \displaystyle \sum_{n=1}^N p_n = P \\ & \displaystyle p_n \geq 0 \end{array} \]

Formulation in python

Now, we have the problem formulation down, and we can proceed to write the code to actually solve the problem. We start with the boilerplate, where we import the necessary package, namely, NumPy for numerical operations and CVXPY for solving the convex optimization problem.

import cvxpy as cp
import numpy as np

Next, we need to set up the simulation environment. We gave the number of channels (N) and SNR in dB as input parameters. Then, we derive the power budget (P) and by fixing the unit noise floor. There would be other means to derive the relationship between the SNR, noise, and power. It's all pretty much a matter of definition. Finally, we draw the channel gains from the positive side of the Gaussian distribution with zero mean and unit variance.

# Number of channels 
N = 10
N0 = 1 # Normalized noise level
SNR_dB = 10 # The signal to noise ratio in dB
P = 10**(SNR_dB/10) # Sum power budget defined via the SNR

# The channel specific gains drawn from Gaussian distribution
g = np.abs(np.random.randn(N, 1)) 
G = np.diag(g[0:,0]) # Make gains a diagonal matrix

We reformulated the channel gains to a diagonal matrix with the channel gains on the diagonal. This will be helpful in the next step, where we write the optimization objective using CVXPY.

The optimization problem is formulated using the CVXPY package as

# Problem construction

# The individual power allocations as variables (N x 1 real vector)
p = cp.Variable(N)

objective = cp.Maximize(cp.sum(cp.log(1 + (G @ p) / N0)))
constraints = []
constraints.append(p >= 0) # Positive powers constraint
constraints.append(cp.sum(p) <= P) # Total power budget

problem = cp.Problem(objective, constraints)

CVXPY uses @ operator for matrix-matrix and matrix-vector multiplication. This is something that you need to be careful with. Otherwise, we will get a warning for using deprecated notation. I've replaced log2 it with the natural logarithm as it is really just scaling and does not affect the solution.

Finally, we are ready to solve the problem with

# Solve the problem
result = problem.solve()

When the problem has been solved, we can print out the solution

# Print the optimal sum rate 
# This has now natural log basis making the unit nats
print('Sum-rate (capacity) in nats', result)

# Show the power allocation
print('Power allocation', p.value)

The complete code can be found here.

Conclusion

We were able to simulate the sum-rate maximization problem for orthogonal channels by using a couple of lines of Python, where half of the code was for setting up the time environment and printing the results.