doc/source/examples/applications/maximise_minimum_SINR_BV4.20.rst
by Robert Gowers, Roger Hill, Sami Al-Izzi, Timothy Pollington and Keith Briggs
from Boyd and Vandenberghe, Convex Optimization, exercise 4.20 page 196
Convex optimization can be used to maximise the minimum signal to
inteference plus noise ratio (SINR) of a wireless communication system.
Consider a system with :math:n transmitters, each with power
:math:p_j \geq 0, transmitting to :math:n receivers. Let
:math:G_{ij} \geq 0 denote the path gain from transmitter :math:j to
receiver :math:i. These path gains form the matrix
:math:G \in \mathbb{R}^{n \times n}.
Each receiver is assigned to a transmitter such that the signal power at
receiver :math:i, :math:S_i = G_{ii}p_i and the interefence power at
receiver :math:i is :math:I_i = \sum_{k\neq i} G_{ik}p_k. Given a
noise power :math:\sigma_i at each receiver, the SINR at receiver
:math:i, :math:\gamma_i = \frac{S_i}{I_i + \sigma_i}.
The objective is to maximise the minimum SINR of the system under certain power constraints. These constraints are:
i - Each transmitter power :math:p_j \leq P_j^{\text{max}}
ii - If the transmitters are partitioned into :math:m nonoverlapping
groups, :math:K_1, ..., K_m, which share a common power supply with
total power :math:P_l^{\text{gp}}:
:math:\sum_{k\in K_l}p_k \leq P_l^{\text{gp}}.
iii - There is a maximum power that each receiver can receive
:math:P_i^{\text{rc}},
:math:\sum_{k=1}^{n}G_{ik}p_k \leq P_i^{\text{rc}}.
The objective function can be rewritten as:
minimise :math:\max_{i=1,...,n}\frac{I_i + \sigma_i}{S_i}
However, since this is a quasiconvex objective function we cannot solve
it directly using CVXPY. Instead we must use a bisection method. First
we take the step of rewriting the objective,
:math:\alpha = \gamma^{-1} \geq 0, as a constraint:
:math:I_i+\sigma_i \leq S_i\alpha
Then we choose initial lower and upper bounds :math:L_0 and
:math:U_0 for :math:\alpha, which should be chosen such that
:math:L < \alpha^* < U, where :math:\alpha^* is the optimal value of
:math:\alpha. Starting with an initial value
:math:\alpha_0 = \frac{1}{2}(L_0+U_0), feasibility is checked for
:math:\alpha_0 by using an arbitrary objective function. The new upper
and lower bounds are determined from the feasibility:
If :math:\alpha_0 is feasible then :math:L_1 = L_0,
:math:U_1 = \alpha_0 and :math:\alpha_1 = \frac{1}{2}(L_1+U_1).
If :math:\alpha_0 is infeasible then :math:L_1 = \alpha_1,
:math:U_1 = U_0 and :math:\alpha_1 = \frac{1}{2}(L_1+U_1).
This bisection process is repeated until :math:U_N - L_N < \epsilon,
where :math:\epsilon is the desired tolerance.
.. code:: python
#!/usr/bin/env python3
# @author: R. Gowers, S. Al-Izzi, T. Pollington, R. Hill & K. Briggs
import cvxpy as cp
import numpy as np
.. code:: python
def maxmin_sinr(G, P_max, P_received, sigma, Group, Group_max, epsilon = 0.001):
# find n and m from the size of the path gain matrix
n, m = np.shape(G)
# Checks sizes of inputs
if m != np.size(P_max):
print('Error: P_max dimensions do not match gain matrix dimensions\n')
return 'Error: P_max dimensions do not match gain matrix dimensions\n', np.nan, np.nan, np.nan
if n != np.size(P_received):
print('Error: P_received dimensions do not match gain matrix dimensions\n')
return 'Error: P_received dimensions do not match gain matrix dimensions', np.nan, np.nan, np.nan
if n != np.size(sigma):
print('Error: σ dimensions do not match gain matrix dimensions\n')
return 'Error: σ dimensions do not match gain matrix dimensions', np.nan, np.nan, np.nan
#I = np.zeros((n,m))
#S = np.zeros((n,m))
delta = np.identity(n)
S = G*delta # signal power matrix
I = G-S # interference power matrix
# group matrix: number of groups by number of transmitters
num_groups = int(np.size(Group,0))
if num_groups != np.size(Group_max):
print('Error: Number of groups from Group matrix does not match dimensions of Group_max\n')
return ('Error: Number of groups from Group matrix does not match dimensions of Group_max',
np.nan, np.nan, np.nan, np.nan)
# normalising the max power of a group so it is in the range [0,1]
Group_norm = Group/np.sum(Group,axis=1).reshape((num_groups,1))
# create scalar optimisation variable p: the power of the n transmitters
p = cp.Variable(shape=n)
best = np.zeros(n)
# set upper and lower bounds for sub-level set
u = 1e4
l = 0
# alpha defines the sub-level sets of the generalised linear fractional problem
# in this case α is the reciprocal of the minimum SINR
alpha = cp.Parameter(shape=1)
# set up the constraints for the bisection feasibility test
constraints = [I*p + sigma <= alpha*S*p, p <= P_max, p >= 0, G*p <= P_received, Group_norm*p <= Group_max]
# define objective function, in our case it's constant as only want to test the solution's feasibility
obj = cp.Minimize(alpha)
# now check whether the solution lies between u and l
alpha.value = [u]
prob = cp.Problem(obj, constraints)
prob.solve()
if prob.status != 'optimal':
# in this case the level set u is below the solution
print('No optimal solution within bounds\n')
return 'Error: no optimal solution within bounds', np.nan, np.nan, np.nan
alpha.value = [l]
prob = cp.Problem(obj, constraints)
prob.solve()
if prob.status == 'optimal':
# in this case the level set l is below the solution
print('No optimal solution within bounds\n')
return 'Error: no optimal solution within bounds', np.nan, np.nan, np.nan
# Bisection algortithm starts
maxLoop = int(1e7)
for i in range(1,maxLoop):
# First check that u is in the feasible domain and l is not, loop finishes here if this is not the case
# set α as the midpoint of the interval
alpha.value = np.atleast_1d((u + l)/2.0)
# test the size of the interval against the specified tolerance
if u-l <= epsilon:
break
# form and solve problem
prob = cp.Problem(obj, constraints)
prob.solve()
# If the problem is feasible u -> α, if not l -> α, best takes the last feasible value as the optimal one as
# when the tolerance is reached the new α may be out of bounds
if prob.status == 'optimal':
u = alpha.value
best = p.value
else:
l = alpha.value
# final condition to check that the interval has converged to order ε, i.e. the range of the optimal sublevel set is <=ε
if u - l > epsilon and i == (maxLoop-1):
print("Solution not converged to order epsilon")
return l, u, float(alpha.value), best
As a simple example, we will consider a case with :math:n=5, where
:math:G_{ij} = 0.6 if :math:i=j and :math:0.1 otherwise.
:math:P_j^{\text{max}} = 1 for all transmitters and the transmitters
are split into two groups, each with :math:P_l^{\text{gp}} = 1.8. The
first group contains transmitters 1 & 2, while the second group contains
3,4 & 5.
For all receivers :math:P_i^{\text{rc}} = 4 and
:math:\sigma_i = 0.1.
.. code:: python
np.set_printoptions(precision=3)
# in this case we will use a gain matrix with a signal weight of 0.6 and interference weight of 0.1
G = np.array([[0.6,0.1,0.1,0.1,0.1],
[0.1,0.6,0.1,0.1,0.1],
[0.1,0.1,0.6,0.1,0.1],
[0.1,0.1,0.1,0.6,0.1],
[0.1,0.1,0.1,0.1,0.6]])
# in this case m=n, but this generalises if we want n receivers and m transmitters
n, m = np.shape(G)
# set maximum power of each transmitter and receiver saturation level
P_max = np.array([1.]*n)
# normalised received power, total possible would be all power from all transmitters so 1/n
P_received = np.array([4.,4.,4.,4.,4.])/n
# set noise level
sigma = np.array([0.1,0.1,0.1,0.1,0.1])
# group matrix: number of groups by number of transmitters
Group = np.array([[1.,1.,0,0,0],[0,0,1.,1.,1.]])
# max normalised power for groups, number of groups by 1
Group_max = np.array([1.8,1.8])
# now run the optimisation problem
l, u, alpha, best = maxmin_sinr(G, P_max, P_received, sigma, Group, Group_max)
print('Minimum SINR={:.4g}'.format(1/alpha))
print('Power={}'.format(best))
.. parsed-literal::
Minimum SINR=1.148
Power=[0.8 0.8 0.8 0.8 0.8]