MOSEK ApS

Optimization for multiuser communications

We consider a multiple-input multiple-output (MIMO) communication system with $n$ transmitters and $n$ receivers. Each transmitter transmits with power $p_j$ and the gain from transmitter $j$ to receiver $i$ is $G_{ij}$. The signal power from transmitter $i$ to receiver $i$ is then

$$ S_i = G_{ii} p_i $$

and the interference is

$$ I_i = \sum_{j\neq i} G_{ij} p_j + \sigma_i $$

where $\sigma_i$ is an additive noise component. In this notebook we consider different strategies for optimizing the signal-to-inference-plus-noise ratio (SINR)

$$ s_i = \frac{G_{ii} p_i}{\sum_{j\neq i} G_{ij} p_j + \sigma_i} $$

with a bound on the total transmitted power $ \sum_i p_i \leq P $.

Minimizing total power for given SINRs

Suppose we are given lower bounds $s_i \geq \gamma_i$. We can then minimize the required power $$ \begin{array}{ll} \text{minimize} & \sum_i p_i \\ \text{subject to} & s_i \geq \gamma_i \\ & \sum_i p_i \leq P, \end{array} $$ which is equivalent to a linear optimization problem $$ \begin{array}{ll} \text{minimize} & \sum_i p_i \\ \text{subject to} & G_{ii} p_i \geq \gamma_i\left ( \sum_{j\neq i} G_{ij} p_j + \sigma_i \right ) \\ & \sum_i p_i \leq P. \end{array} $$

Maximizing the worst SINR

Alternatively we can maximize the smallest $s_i$,

$$ \begin{array}{ll} \text{maximize} & t \\ \text{subject to} & s_i \geq t \\ & \sum_i p_i \leq P. \end{array} $$

Equivalently we can minimize the inverse,

$$ \begin{array}{ll} \text{minimize} & t^{-1} \\ \text{subject to} & t \left ( \sum_{j\neq i} G_{ij} p_j + \sigma_i \right ) G_{ii}^{-1} p_i^{-1} \leq 1 \\ & \sum_i p_i \leq P, \end{array} $$

which can be rewritten as a geometric programming problem

$$ \begin{array}{ll} \text{minimize} & -z\\ \text{subject to} & \log \left ( \sum_{j\neq i}e^{z + q_j - q_i + \log(G_{ij}/G_{ii})} + e^{z - q_i + \log(\sigma_i/G_{ii})} \right ) \leq 0\\ & \log \left ( \sum_i e^{q_i-\log P}\right) \leq 0 \end{array} $$

with $p_i := e^{q_i}$ and $t := e^z$. To rewrite the geometric program into conic form, we note that

$$ \log \left( \sum_{i=1}^n e^{a_i^T x + b_i}\right) \leq 0 \qquad \Longleftrightarrow \qquad \sum_i u_i\leq 1, \quad (u_i, 1, a_i^Tx + b_i)\in K_\text{exp}, \: i=1,\dots n. $$

In [1]:
import sys
import numpy as np
from mosek.fusion import *
import matplotlib.pyplot as plt
In [2]:
def logsumexp(M, A, x, b):    
    u = M.variable(A.shape[0])
    M.constraint( Expr.sum(u), Domain.lessThan(1.0))
    M.constraint( Expr.hstack(u,
                              Expr.constTerm(A.shape[0], 1.0),
                              Expr.add(Expr.mul(A, x), b)), Domain.inPExpCone())
In [3]:
def max_worst_sinr(G, P, sigma):
    n = G.shape[0]
    with Model('MAX_WORST_SINR') as M:
        qz = M.variable('qz', n+1) # concatenation of q and z
        M.objective('Objective',ObjectiveSense.Minimize,Expr.neg(qz.index(n)))
        for i in range(n):
            A = np.zeros((n,n+1))
            b = np.zeros(n)
            for j in [k for k in range(n) if k!=i]:
                A[j,[i,j,n]] = [-1, 1, 1]
                b[j] = G[i,j]/G[i,i]
            A[i, [i, n]] = [-1, 1]
            b[i] = sigma[i]/G[i,i]
            # If any Gij == 0, then we filter out row j
            idx = np.nonzero(b)[0] 
            logsumexp(M, A[idx,:], qz, np.log(b[idx]))
    
        logsumexp(M, np.eye(n), qz.slice(0,n), -np.log(P)*np.ones(n))
        M.setLogHandler(sys.stdout)

        M.solve()
        pt = np.exp(qz.level())
        return (pt[0:n], pt[n])
In [4]:
P = 0.5

G = np.array([[1.0,0.1,0.2,0.1,0.0],
              [0.1,1.0,0.1,0.1,0.0],
              [0.2,0.1,2.0,0.2,0.2],
              [0.1,0.1,0.2,1.0,0.1],
              [0.0,0.0,0.2,0.1,1.0]])

sigma = 0.01*np.ones(G.shape[0])
In [5]:
p1, t1 = max_worst_sinr(G, P, sigma)
Problem
  Name                   : MAX_WORST_SINR  
  Objective sense        : min             
  Type                   : CONIC (conic optimization problem)
  Constraints            : 84              
  Cones                  : 26              
  Scalar variables       : 110             
  Matrix variables       : 0               
  Integer variables      : 0               

Optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator started.
Freed constraints in eliminator : 1
Eliminator terminated.
Eliminator - tries                  : 1                 time                   : 0.00            
Lin. dep.  - tries                  : 1                 time                   : 0.00            
Lin. dep.  - number                 : 0               
Presolve terminated. Time: 0.01    
Problem
  Name                   : MAX_WORST_SINR  
  Objective sense        : min             
  Type                   : CONIC (conic optimization problem)
  Constraints            : 84              
  Cones                  : 26              
  Scalar variables       : 110             
  Matrix variables       : 0               
  Integer variables      : 0               

Optimizer  - threads                : 20              
Optimizer  - solved problem         : the primal      
Optimizer  - Constraints            : 26
Optimizer  - Cones                  : 26
Optimizer  - Scalar variables       : 84                conic                  : 78              
Optimizer  - Semi-definite variables: 0                 scalarized             : 0               
Factor     - setup time             : 0.00              dense det. time        : 0.00            
Factor     - ML order time          : 0.00              GP order time          : 0.00            
Factor     - nonzeros before factor : 273               after factor           : 274             
Factor     - dense dim.             : 0                 flops                  : 6.44e+03        
ITE PFEAS    DFEAS    GFEAS    PRSTATUS   POBJ              DOBJ              MU       TIME  
0   8.3e+00  2.0e+00  3.6e+01  0.00e+00   -2.256346207e+00  -2.484467505e+01  2.3e+00  0.03  
1   9.8e-01  2.3e-01  2.5e+00  2.84e-01   -1.209531123e-01  -3.870013357e+00  2.7e-01  0.05  
2   2.0e-01  4.8e-02  9.5e-02  1.53e+00   -1.990579522e-01  -7.785188954e-01  5.6e-02  0.05  
3   4.9e-02  1.2e-02  1.1e-02  1.43e+00   -6.685338326e-01  -7.871072260e-01  1.4e-02  0.05  
4   9.1e-03  2.1e-03  9.3e-04  1.10e+00   -7.570036183e-01  -7.776695214e-01  2.5e-03  0.05  
5   6.7e-04  1.6e-04  1.9e-05  1.02e+00   -7.751417794e-01  -7.766579413e-01  1.9e-04  0.06  
6   1.5e-05  3.6e-06  6.4e-08  1.00e+00   -7.765076649e-01  -7.765421638e-01  4.3e-06  0.06  
7   2.3e-06  5.3e-07  3.7e-09  1.00e+00   -7.765307842e-01  -7.765358682e-01  6.3e-07  0.06  
8   1.6e-07  3.8e-08  7.0e-11  1.00e+00   -7.765341634e-01  -7.765345269e-01  4.5e-08  0.06  
9   2.1e-08  5.0e-09  3.3e-12  1.00e+00   -7.765343712e-01  -7.765344185e-01  5.8e-09  0.06  
10  1.8e-08  4.1e-09  2.5e-12  1.00e+00   -7.765343736e-01  -7.765344130e-01  4.7e-09  0.06  
Optimizer terminated. Time: 0.07    


Interior-point solution summary
  Problem status  : PRIMAL_AND_DUAL_FEASIBLE
  Solution status : OPTIMAL
  Primal.  obj: -7.7653437363e-01   nrm: 5e+00    Viol.  con: 1e-08    var: 0e+00    cones: 0e+00  
  Dual.    obj: -7.7653441298e-01   nrm: 3e-01    Viol.  con: 3e-17    var: 1e-09    cones: 0e+00  
In [6]:
p1, t1
Out[6]:
(array([0.10756629, 0.09148653, 0.09004626, 0.12322311, 0.0876778 ]),
 2.173925182931828)
In [7]:
SINR1 = (np.diagonal(G)*p1)/(np.dot(G,p1) - np.diagonal(G)*p1 + sigma)
SINR1
Out[7]:
array([2.17392523, 2.17392526, 2.17392523, 2.17392522, 2.17392524])

Maximizing the best SINR

The solution to $$ \begin{array}{ll} \text{maximize} & t_i \\ \text{subject to} & s_i \leq t_i \\ & \sum_i p_i \leq P \end{array} $$ is trivial; we choose the index $k$ maximizing $P_{ii}/\sigma_i$ and take $p_k=P$ and $p_j=0,\: j\neq k$.

In [8]:
def max_best_SINR(G,P,sigma):
    GSD = [G[i][i]/sigma[i] for i in range(G.shape[0])]
    P_max = max(GSD)
    #Thus, maximum of the best SINR is equal to...
    return(P_max*P)
In [9]:
max_best_SINR(G,P,sigma)
Out[9]:
100.0

Maximizing average SINR

We can maximize the average SINR as $$ \begin{array}{ll} \text{maximize} & \sum_i t_i \\ \text{subject to} & s_i \geq t_i \\ & 0 \leq p_i \leq P_i \\ & \sum_i p_i \leq P, \end{array} $$ which corresponds to an intractable non-convex bilinear optimization problem. However, in the low-SINR regime, we can approximate the above problem by maximizing $\sum_i \log t_i$, or equivalently minimizing $\prod_i t_i^{-1}$: $$ \begin{array}{ll} \text{minimize} & \prod_i t_i^{-1} \\ \text{subject to} & t_i \left ( \sum_{j\neq i} G_{ij} p_j + \sigma_i \right ) G_{ii}^{-1} p_i^{-1} \leq 1 \\ & 0 \leq p_i \leq P_i \\ & \sum_i p_i \leq P, \end{array} $$ which again corresponds to a geometric programming problem.

In [10]:
def min_Geo_mean(G,P,sigma):
    n = G.shape[0]
    with Model('MIN_GEO_MEAN') as M:
        t = M.variable('t',n)
        x = M.variable('x',n)
        q = M.variable('q',n)
            
        logsumexp(M,np.eye(n),q,-np.log(P)*np.ones(n))
        
        M.constraint(Expr.hstack(x,Expr.constTerm(n, 1.0),Expr.neg(t)),Domain.inPExpCone())
        M.objective('Objective',ObjectiveSense.Minimize,Expr.sum(x))
        
        for i in range(n):
            A = np.zeros((n,n+1))
            b = np.zeros(n)
            for j in [k for k in range(n) if k!=i]:
                A[j,[i,j,n]] = [-1,1,1]
                b[j] = G[i,j]/G[i,i]
            A[i,[i,n]] = [-1,1]
            b[i] = sigma[i]/G[i,i]
            idx = np.nonzero(b)[0]
            logsumexp(M,A[idx,:],Expr.vstack(q,t.index(i)),np.log(b[idx]))
        
        M.setLogHandler(sys.stdout)

        M.solve()
        T = t.level()
        p = np.exp(q.level())
        return(T,p)
In [11]:
t2,p2 = min_Geo_mean(G, P, sigma)
Problem
  Name                   : MIN_GEO_MEAN    
  Objective sense        : min             
  Type                   : CONIC (conic optimization problem)
  Constraints            : 99              
  Cones                  : 31              
  Scalar variables       : 134             
  Matrix variables       : 0               
  Integer variables      : 0               

Optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator - tries                  : 1                 time                   : 0.00            
Lin. dep.  - tries                  : 1                 time                   : 0.00            
Lin. dep.  - number                 : 0               
Presolve terminated. Time: 0.00    
Problem
  Name                   : MIN_GEO_MEAN    
  Objective sense        : min             
  Type                   : CONIC (conic optimization problem)
  Constraints            : 99              
  Cones                  : 31              
  Scalar variables       : 134             
  Matrix variables       : 0               
  Integer variables      : 0               

Optimizer  - threads                : 20              
Optimizer  - solved problem         : the primal      
Optimizer  - Constraints            : 27
Optimizer  - Cones                  : 31
Optimizer  - Scalar variables       : 99                conic                  : 93              
Optimizer  - Semi-definite variables: 0                 scalarized             : 0               
Factor     - setup time             : 0.00              dense det. time        : 0.00            
Factor     - ML order time          : 0.00              GP order time          : 0.00            
Factor     - nonzeros before factor : 185               after factor           : 201             
Factor     - dense dim.             : 0                 flops                  : 3.06e+03        
ITE PFEAS    DFEAS    GFEAS    PRSTATUS   POBJ              DOBJ              MU       TIME  
0   1.2e+01  2.8e+00  7.0e+01  0.00e+00   6.454638549e+00   -2.495816205e+01  4.6e+00  0.01  
1   1.3e+00  3.1e-01  3.5e+00  4.75e-01   4.463773399e+00   6.636911219e-02   5.1e-01  0.02  
2   2.6e-01  6.2e-02  2.3e-01  1.32e+00   2.810914833e+00   2.070129028e+00   1.0e-01  0.02  
3   4.0e-02  9.4e-03  1.3e-02  1.21e+00   2.324042616e+00   2.222420496e+00   1.6e-02  0.02  
4   4.6e-03  1.1e-03  4.8e-04  1.04e+00   2.242616061e+00   2.231211157e+00   1.8e-03  0.02  
5   2.0e-04  4.8e-05  4.5e-06  1.01e+00   2.233464306e+00   2.232960386e+00   8.0e-05  0.02  
6   1.5e-05  3.5e-06  8.7e-08  1.00e+00   2.233112252e+00   2.233075749e+00   5.8e-06  0.02  
7   1.8e-06  4.3e-07  3.7e-09  1.00e+00   2.233089609e+00   2.233085129e+00   7.1e-07  0.02  
8   1.5e-07  3.6e-08  9.3e-11  1.00e+00   2.233086826e+00   2.233086444e+00   6.0e-08  0.03  
9   4.1e-08  9.5e-09  1.3e-11  1.00e+00   2.233086648e+00   2.233086548e+00   1.5e-08  0.03  
10  4.3e-08  9.3e-09  1.2e-11  1.00e+00   2.233086650e+00   2.233086552e+00   1.5e-08  0.03  
11  3.6e-08  8.3e-09  1.0e-11  1.00e+00   2.233086644e+00   2.233086557e+00   1.3e-08  0.03  
12  3.6e-08  8.3e-09  1.0e-11  1.00e+00   2.233086644e+00   2.233086557e+00   1.3e-08  0.03  
13  3.6e-08  8.3e-09  1.0e-11  1.00e+00   2.233086644e+00   2.233086557e+00   1.3e-08  0.03  
14  7.6e-09  1.8e-09  1.0e-12  1.00e+00   2.233086591e+00   2.233086572e+00   2.8e-09  0.03  
Optimizer terminated. Time: 0.04    


Interior-point solution summary
  Problem status  : PRIMAL_AND_DUAL_FEASIBLE
  Solution status : OPTIMAL
  Primal.  obj: 2.2330865908e+00    nrm: 5e+00    Viol.  con: 3e-09    var: 0e+00    cones: 0e+00  
  Dual.    obj: 2.2330865723e+00    nrm: 1e+00    Viol.  con: 2e-16    var: 5e-10    cones: 0e+00  
In [12]:
p2,t2
Out[12]:
(array([0.10622214, 0.1062226 , 0.07511098, 0.10622196, 0.10622232]),
 array([0.83111115, 1.00826399, 0.57707338, 0.62443057, 1.09194245]))
In [13]:
SINR2 = (np.diagonal(G)*p2)/(np.dot(G,p2) - np.diagonal(G)*p2 + sigma)
SINR2
Out[13]:
array([2.29586837, 2.74083876, 1.78081902, 1.86718244, 2.98005707])

Comparing the SINR for the cases above...

In [14]:
fig,ax = plt.subplots(figsize = (13,10))
bar_width = 0.35
p_num = np.arange(1,G.shape[0]+1)

B1 = ax.bar(p_num,SINR1,bar_width,label = 'Max Worst SINR')
B2 = ax.bar(p_num+bar_width,SINR2,bar_width,label = 'Max Average SINR')

ax.set_ylabel('SINR')
ax.set_xticks(p_num + bar_width/2)
x_tiK = ['p{}'.format(i+1) for i in range(G.shape[0])]
ax.set_xticklabels(x_tiK)
ax.set_xlabel('Transmitter')
ax.legend()
plt.show()

Creative Commons License
This work is licensed under a Creative Commons Attribution 4.0 International License. The MOSEK logo and name are trademarks of Mosek ApS. The code is provided as-is. Compatibility with future release of MOSEK or the Fusion API are not guaranteed. For more information contact our support.