L. -N. Tran, M. Juntti, M. Bengtsson and B. Ottersten, “Weighted Sum Rate Maximization for MIMO Broadcast Channels Using Dirty Paper Coding and Zero-forcing Methods,” in IEEE Transactions on Communications, vol. 61, no. 6, pp. 2362-2373, June 2013,
我们考虑为最大化连续零强迫脏纸编码(SZF-DPC)的加权和速率(WSR)设计预编码器。对于这个问题,现有的预编码器设计通常假设有总功率约束(SPC),并依赖于奇异值分解(SVD)。基于 SVD 的设计被认为是最优的,但需要高复杂度。我们首先提出了一种低复杂度的最优预编码器设计,用于在 SPC 下的 SZF-DPC,采用 QR 分解。然后,我们提出了一种高效的数值算法,以找到满足每天线功率约束(PAPCs)的最优预编码器。为此,PAPCs 的预编码器设计被表述为一个具有协方差矩阵秩约束的优化问题。解决这个问题的一个著名方法是放宽秩约束并解决放宽后的问题。有趣的是,对于 SZF-DPC,我们能够证明秩放宽是紧的。因此,PAPCs 的最优预编码器设计是通过解决放宽后的问题来计算的,我们提出了一种定制的内点法,具有超线性收敛速率。还提出了两种次优预编码器设计,并与最优设计进行了比较。 我们还表明,所提出的数值方法适用于寻找块对角化方案的最优前置编码器。
代码解释
这段代码实现了一个障碍方法(Barrier Method)优化算法,用于求解具有总功率约束和逐个天线功率约束(Per-Antenna Power Constraint, PAPC)的MIMO(多输入多输出)系统中的预编码矩阵优化问题。其目的是最大化信道中的和速率(Sum Rate)。
1. 基本变量初始化
channel
: 输入的信道矩阵。Po
: 各个天线的功率约束。nTx
: 发射天线的数量。nUsers
: 用户的数量。effectivechannel
: 每个用户的有效信道矩阵,经过处理得到的信道表示。B
: 存储用于正交化信道的矩阵。A
: 用于表示每个天线与用户相关联的矩阵,计算基于特定用户的复合信道。
2. 阶段I和阶段II
-
阶段I:
- 初始化变量,并通过计算牛顿步(Newton Step)来更新原始变量
Omega
和对偶变量mymu
。 - 通过回溯线搜索确定适当的步长,确保更新的变量满足约束。
- 继续迭代直到达到容忍的误差阈值。
- 初始化变量,并通过计算牛顿步(Newton Step)来更新原始变量
-
阶段II:
- 在阶段I的基础上,通过进一步增加障碍函数的权重
t
,继续进行优化。 - 每次迭代都会更新
Omega
和mymu
,直到达到设定的最大值或达到收敛标准。
- 在阶段I的基础上,通过进一步增加障碍函数的权重
3. 辅助函数
computeresidualnorm
: 计算残差的范数,用于判断算法的收敛性。computeinfeasiblestartNewtonstep
: 计算用于更新原始和对偶变量的牛顿步。computeobjval
: 计算目标函数的值,即和速率的负值。computemodifiedobjval
: 计算带有障碍项的目标函数的值,同时计算对应的梯度。
总结
这段代码实现了一个基于障碍方法的复杂优化算法,主要用于MIMO系统的预编码矩阵设计。在求解过程中,算法需要解决具有逐个天线功率约束的优化问题,通过逐步更新原始和对偶变量,最大化系统的和速率。
function [SumRate,Precoders,gap] = szfdpcpapcbarriermethod(channel,Po)
global A effectivechannel
nTx = size(channel{1},2);
nUsers = length(channel);
effectivechannel = cell(nUsers,1); % 用户的有效信道
B = cell(nUsers,1);
B{1} = eye(nTx);
A = cell(nUsers,nTx); % (33)中的矩阵A
h_bar=[];
Nr=zeros(1,nUsers);
for iUser=1:nUsers
h_bar = [h_bar;channel{iUser}];
if( iUser < nUsers)
B{iUser+1} = null(h_bar); % B{iUser} 是从用户1到(iUser-1)的复合信道的正交基
end
effectivechannel{iUser} = channel{iUser}*B{iUser};
for iTx=1:nTx
A{iUser,iTx} = B{iUser}'*diag(1:nTx==iTx)*B{iUser}; % 计算 (33) 中的矩阵 A
end
end
initialOmega=cell(nUsers,1); % (31)中 \Omega 的初始值
for iUser=1:nUsers
initialOmega{iUser} = eye(size(B{iUser},2)); % 设置为单位矩阵(也可以是其他值)
end
MAXITERS = 100;
ALPHA = 0.01; % 回溯线搜索参数
BETA = 0.5; % 回溯线搜索参数
NTTOL = 1e-5;
myphi = 50;
TOL = 1e-5;
%% 初始化
Omega=cell(nUsers,1);
Omega_next=cell(nUsers,1);
DeltaOmega=cell(nUsers,1);
t=1; % 算法2中的初始 t0 值
for iUser = 1:nUsers
Omega{iUser} = initialOmega{iUser};
DeltaOmega{iUser} = zeros(size(Omega{iUser}));
end
mymu = zeros(nTx,1);
mymunext = zeros(nTx,1);
iCenteringstep = 1;
iTotaliteration = 0;
sumratebest = [];
%% 阶段I
for iIteration = 1:MAXITERS
iTotaliteration = iTotaliteration+1;
modifiedobjval = computemodifiedobjval(t,Omega);
sumrate(iTotaliteration) = computeobjval(Omega);
sumratebest = [sumratebest min(sumrate)];
sumratebarrier(iTotaliteration) = modifiedobjval;
[DeltaOmega,Deltamu] = computeinfeasiblestartNewtonstep(t,Omega,mymu,Po);
residualnorm = computeresidualnorm(t,Omega,mymu,effectivechannel,A,Po);
if (residualnorm < NTTOL),
break; % 如果 residualnorm 很小,则退出 for 循环
end;
s=1;
while(1) % 回溯线搜索
feasibility_flag = 0;
for iUser = 1:nUsers
feasibility_flag = feasibility_flag+(min(real(eig(Omega{iUser}+s*DeltaOmega{iUser}))) <= 0);
Omega_next{iUser} = Omega{iUser}+s*DeltaOmega{iUser};
end
mymunext = mymu+s*Deltamu;
r_norm_next = computeresidualnorm(t,Omega_next,mymunext,effectivechannel,A,Po);
minimum_flag = (r_norm_next > (1-ALPHA*s)*residualnorm);
if(feasibility_flag==0) && (minimum_flag==0)
break
else
s=s*BETA;
end
end
% 更新原始变量和对偶变量
mymu = mymunext;
for iUser=1:nUsers
Omega{iUser}=Omega_next{iUser};
end
end
% 第一步居中步骤结束
iCenteringstep = iCenteringstep + 1;
t=t*myphi;
%% 阶段II
while(t<1e8)
% 内循环开始
for iIteration=1:MAXITERS
iTotaliteration=iTotaliteration+1;
[modifiedobjval,Grad_barrier] = computemodifiedobjval(t,Omega);
[DeltaOmega,Deltamu] = computeinfeasiblestartNewtonstep(t,Omega,mymu,Po);
sumrate(iTotaliteration) = computeobjval(Omega);
sumratebest=[sumratebest min(sumrate)];
sumratebarrier(iTotaliteration) = modifiedobjval;
% 计算牛顿减量
Newtondecrement = 0;
for iUser=1:nUsers
Newtondecrement = Newtondecrement + ...
real(trace(Grad_barrier{iUser}' * DeltaOmega{iUser}));
end
% 回溯线搜索
if (abs(Newtondecrement) < NTTOL),
% SumRate_barrier_centering(iCenteringstep) = modifiedobjval;
break;
end;
s=1;
while(1)
feasibility_flag=0;
for iUser=1:nUsers
feasibility_flag=feasibility_flag+(min(real(eig(Omega{iUser}+s*DeltaOmega{iUser}))) <= 0);
end
for iUser=1:nUsers
Omega_next{iUser}=Omega{iUser}+s*DeltaOmega{iUser};
end
mymunext=mymu+s*Deltamu;
val_next = computemodifiedobjval(t,Omega_next);
minimum_flag =(val_next > modifiedobjval + s*ALPHA*Newtondecrement);
if(feasibility_flag==0) && (minimum_flag==0)
break
else
s=s*BETA;
end
end
% 更新原始变量和对偶变量
mymu = mymunext;
for iUser=1:nUsers
Omega{iUser} = Omega_next{iUser};
end
end
% 内循环结束
iCenteringstep = iCenteringstep+1;
t = t*myphi;
end
gap = sumratebest - sumratebest(end);
SumRate = -sumrate(end);
Precoders = Omega;
function residualnorm = computeresidualnorm(t,Omega,Deltamu,effectivechannel, A, Po)
nUsers=length(Omega);
Grad_orginal = cell(nUsers,1);
Grad_barrier = cell(nUsers,1);
rdual = cell(nUsers,1);
for iUser = 1:nUsers
[Nr]=size(Omega{iUser},1);
Grad_orginal{iUser} = -effectivechannel{iUser}'*effectivechannel{iUser}...
*inv(eye(Nr)+Omega{iUser}*effectivechannel{iUser}'*effectivechannel{iUser});
Grad_barrier{iUser} = t*Grad_orginal{iUser}-inv(Omega{iUser});
rdual{iUser} = Grad_barrier{iUser};
end
nTx=size(Omega{1},2);
v=zeros(nTx,1);
for i = 1:nTx
v(i) = Po(i);
for iUser = 1:nUsers
rdual{iUser} = rdual{iUser}+Deltamu(i)*A{iUser,i};
v(i) = v(i)-real(trace(A{iUser,i}*Omega{iUser}));
end
end
r = zeros(nUsers,1);
for iUser = 1:nUsers
r(iUser) = r(iUser)+norm(rdual{iUser},'fro');
end
residualnorm = sum(r) + norm(v);
function [DeltaOmega,Deltamu] = computeinfeasiblestartNewtonstep(t,Omega,mymu,Po)
global effectivechannel A
nUsers=length(Omega);
Grad_orginal=cell(nUsers,1);
nTx=size(Omega{1},1);
myXi=cell(nUsers,nTx+1);
DeltaOmega=cell(nUsers,1);
for iUser=1:nUsers
[Nr]=size(Omega{iUser},1);
Grad_orginal{iUser}= -effectivechannel{iUser}'*effectivechannel{iUser}*inv(eye(Nr)+Omega{iUser}*effectivechannel{iUser}'*effectivechannel{iUser});
P = t*Omega{iUser}*(-Grad_orginal{iUser});
C = t*Omega{iUser}*(-Grad_orginal{iUser})*Omega{iUser}+Omega{iUser};
for i=1:nTx
C = C - mymu(i)*Omega{iUser}*A{iUser,i}*Omega{iUser};
end
myXi{iUser,1}=dlyap(P,-P'/t,C); % 求解 (42) 的第一个方程
for i=1:nTx
P=t*Omega{iUser}*(-Grad_orginal{iUser});
C=-Omega{iUser}*A{iUser,i}*Omega{iUser};
myXi{iUser,i+1}=dlyap(P,-P'/t,C); % 求解 (42) 的其余方程
end
end
myPsi=zeros(nTx,nTx);
mypsi=zeros(nTx,1);
for i=1:nTx
mypsi(i)=Po(i);
for iUser=1=nUsers
mypsi(i) = mypsi(i)-trace(A{iUser,i} * myXi{iUser,1})-trace(A{iUser,i}*Omega{iUser});
for j=1:nTx
myPsi(i,j)=myPsi(i,j)+trace(A{iUser,i} * myXi{iUser,j+1});
end
end
end
% 对偶变量的牛顿步
Deltamu=real(myPsi\mypsi);
% 原始变量的牛顿步
for iUser=1=nUsers
DeltaOmega{iUser}=myXi{iUser,1};
for i=1:nTx
DeltaOmega{iUser}=DeltaOmega{iUser}+Deltamu(i)*myXi{iUser,i+1};
DeltaOmega{iUser} = (DeltaOmega{iUser}+DeltaOmega{iUser}')/2; % 消除数值误差
end
end
%%
function [val]=computeobjval(Omega)
global effectivechannel
nUsers = length(Omega);
val = 0;
for iUser = 1:nUsers
Nr = size(effectivechannel{iUser},1);
val = val+log(det(eye(Nr) + effectivechannel{iUser} ...
* Omega{iUser} * effectivechannel{iUser}'));
end
val = -real(val);
%%
function [ modifiedobjval ,Grad_barrier] = computemodifiedobjval(t,Omega)
global effectivechannel
nUsers = length(Omega);
Grad_orginal = cell(nUsers,1);
Grad_barrier = cell(nUsers,1);
modifiedobjval = 0;
for iUser = 1:nUsers
[Nr,Nr_bar] = size(effectivechannel{iUser});
modifiedobjval = modifiedobjval+t*log(det(eye(Nr) + effectivechannel{iUser} * Omega{iUser} *...
effectivechannel{iUser}'))+log(det(Omega{iUser}));
Grad_orginal{iUser} = -effectivechannel{iUser}'*effectivechannel{iUser}*...
inv(eye(Nr_bar)+Omega{iUser}*effectivechannel{iUser}'*effectivechannel{iUser});
Grad_barrier{iUser} = t*Grad_orginal{iUser}-inv(Omega{iUser});
end
modifiedobjval = -real(modifiedobjval);
clear all
clc
rng(1)
nTx = 6; % number of users
nRxarray = [2 2 2]; % number of receive antennas for each user, i.e.,
% % nRxarray(i) is the number of rx. antennas at the ith user
% nRxarray = [2 2 ]; % number of receive antennas for each user, i.e.,
nUsers = length(nRxarray); % number of users
P = 10; % sum power, dB scale
P = 10.^(P./10);
Po = P/nTx*ones(nTx,1); % maximum power per antenna
channel = cell(nUsers,1); % channel matrices of users
tic
for iUser=1:nUsers
channel{iUser} = sqrt(1/2)*(randn(nRxarray(iUser),nTx) + 1i*randn(nRxarray(iUser),nTx));
end
[sumrate,precoder,gap] = szfdpcpapcbarriermethod(channel,Po);
toc
semilogy(gap)
xlabel('Iteration')
ylabel('Residual error')
saveas(gcf,'../results/convergence.png')
%% YALMIP code for solving (31)
effectivechannel = cell(nUsers,1);
effectivechannel{1} = channel{1};
B = cell(nUsers,1);
B{1} = eye(nTx);
h_bar = [];
A = cell(nUsers,nTx); % matrix A in (33)
for iUser=1:nUsers
h_bar = [h_bar;channel{iUser}];
if( iUser < nUsers)
B{iUser+1} = null(h_bar);
end
effectivechannel{iUser} = channel{iUser}*B{iUser};
for iTx=1:nTx
A{iUser,iTx}=B{iUser}'*diag(1:nTx==iTx)*B{iUser}; % compute matrix A in (33)
end
end
X=cell(nUsers,1);
F=[];
for iUser = 1:nUsers
X{iUser} = sdpvar(nTx-sum(nRxarray(1:iUser-1)),nTx-sum(nRxarray(1:iUser-1)),'hermitian','complex');
F=[F,X{iUser}>=0];
end
obj=0;
for iUser=1:nUsers
obj=obj+logdet(eye(nRxarray(iUser)) + effectivechannel{iUser}*X{iUser}*effectivechannel{iUser}');
end
for iTx=1:nTx
y = 0;
for iUser=1:nUsers
y = y + real(trace(A{iUser,iTx}*X{iUser}));
end
F=[F,y <= Po(iTx)];
end
ops = sdpsettings('solver','sdpt3','verbose',0);
tic
sol = optimize(F,-obj,ops);
sumrateyalmip = real(double(obj));
precodersoftware = cell(nUsers,1);
for iUser = 1:nUsers
precodersoftware{iUser} = B{iUser}*double(X{iUser})*B{iUser}';
end
toc
norm(sumrate-sumrateyalmip)