by Robert Gowers, Roger Hill, Sami Al-Izzi, Timothy Pollington and Keith Briggs
from Boyd and Vandenberghe, Convex Optimization, exercise 4.62 page 210
Consider a system in which a central node transmits messages to $n$ receivers. Each receiver channel $i \in \{1,...,n\}$ has a transmit power $P_i$ and bandwidth $W_i$. A fraction of the total power and bandwidth is allocated to each channel, such that $\sum_{i=1}^{n}P_i = P_{tot}$ and $\sum_{i=1}^{n}W_i = W_{tot}$. Given some utility function of the bit rate of each channel, $u_i(R_i)$, the objective is to maximise the total utility $U = \sum_{i=1}^{n}u_i(R_i)$.
Assuming that each channel is corrupted by Gaussian white noise, the signal to noise ratio is given by $\beta_i P_i/W_i$. This means that the bit rate is given by:
$R_i = \alpha_i W_i \log_2(1+\beta_iP_i/W_i)$
where $\alpha_i$ and $\beta_i$ are known positive constants.
One of the simplest utility functions is the data rate itself, which also gives a convex objective function.
The optimisation problem can be thus be formulated as:
minimise $\sum_{i=1}^{n}-\alpha_i W_i \log_2(1+\beta_iP_i/W_i)$
subject to $\sum_{i=1}^{n}P_i = P_{tot} \quad \sum_{i=1}^{n}W_i = W_{tot} \quad P \succeq 0 \quad W \succeq 0$
Although this is a convex optimisation problem, it must be rewritten in DCP form since $P_i$ and $W_i$ are variables and DCP prohibits dividing one variable by another directly. In order to rewrite the problem in DCP format, we utilise the $\texttt{kl_div}$ function in CVXPY, which calculates the Kullback-Leibler divergence.
$\text{kl_div}(x,y) = x\log(x/y)-x+y$
$-R_i = \text{kl_div}(\alpha_i W_i, \alpha_i(W_i+\beta_iP_i)) - \alpha_i\beta_iP_i$
Now that the objective function is in DCP form, the problem can be solved using CVXPY.
#!/usr/bin/env python3
# @author: R. Gowers, S. Al-Izzi, T. Pollington, R. Hill & K. Briggs
import numpy as np
import cvxpy as cvx
def optimal_power(n, a_val, b_val, P_tot=1.0, W_tot=1.0):
# Input parameters: α and β are constants from R_i equation
n=len(a_val)
if n!=len(b_val):
print('alpha and beta vectors must have same length!')
return 'failed',np.nan,np.nan,np.nan
P=cvx.Variable(n)
W=cvx.Variable(n)
alpha=cvx.Parameter(n)
beta =cvx.Parameter(n)
alpha.value=np.array(a_val)
beta.value =np.array(b_val)
# This function will be used as the objective so must be DCP;
# i.e. elementwise multiplication must occur inside kl_div, not outside otherwise the solver does not know if it is DCP...
R=cvx.kl_div(cvx.mul_elemwise(alpha, W),
cvx.mul_elemwise(alpha, W + cvx.mul_elemwise(beta, P))) - \
cvx.mul_elemwise(alpha, cvx.mul_elemwise(beta, P))
objective=cvx.Minimize(cvx.sum_entries(R))
constraints=[P>=0.0,
W>=0.0,
cvx.sum_entries(P)-P_tot==0.0,
cvx.sum_entries(W)-W_tot==0.0]
prob=cvx.Problem(objective, constraints)
prob.solve()
return prob.status,-prob.value,P.value,W.value
Consider the case where there are 5 channels, $n=5$, $\alpha = \beta = (2.0,2.2,2.4,2.6,2.8)$, $P_{\text{tot}} = 0.5$ and $W_{\text{tot}}=1$.
np.set_printoptions(precision=3)
n=5 # number of receivers in the system
a_val=np.arange(10,n+10)/(1.0*n) # α
b_val=np.arange(10,n+10)/(1.0*n) # β
P_tot=0.5
W_tot=1.0
status,utility,power,bandwidth=optimal_power(n,a_val,b_val,P_tot,W_tot)
print('Status: ',status)
print('Optimal utility value = %.4g '%utility)
print('Optimal power level:\n', power)
print('Optimal bandwidth:\n', bandwidth)