IRS中波束赋形设计源代码之AO算法学习(持续更新，多多交流)

论文：Weighted Sum-Rate Maximization for Reconfigurable Intelligent Surface Aided Wireless Networks
源代码

2020/7/12

完成：
1、搞清楚了论文中向量，矩阵的设置在代码中的对应

2、信道参数设置

function [K,N,M,chanRealNum,pathLossDire,pathLossIndire,...pd,pid,UserArriAngle,ApTranAngle,IrsArriAngle,...Hd_w,theta_init,RicFac,RicFac1,RicFac2,G_sig,Hr_sig] = parameters()
%重要参数生成函数K = 4; %用户数目，单天线
N = 100; %IRS元件个数
M = 4; %发射端发射天线数目
chanRealNum = 100; %信道实现数目
%损耗
pathLossDire = [0.9, 0.95, 0.99, 1.1];
%path_d = 10^((-noise-Ld)/10) 为什么这样算路径功率损耗？
pathLossIndire = [0.00012,0.00014,0.00016,0.00018];
pd = sqrt(pathLossDire);%直接链路损耗(AP-user)
pid = sqrt(pathLossIndire);%间接链路损耗(AP-IRS-user)
pd = repmat(pd.',1,M);%repmat(A,m,n) 将矩阵A复制m x n 块，并且平铺格式是m x n
pid = repmat(pid.',1,N);
%出发角与到达角
%注：因为考虑的是MISO系统，所以用户接收端(单天线)，不考虑到达角
UserArriAngle = rand(1,K); %Hd信道出发角
ApTranAngle= rand(1,1); %G信道出发角
IrsArriAngle = rand(1,1);%G信道到达角%信道Hd(瑞利信道)初始化
Hd_w = zeros(K,M,chanRealNum); %hdk转置
theta_init = zeros(N,1,chanRealNum);
for i = 1:chanRealNumHd = sqrt(1/2).*(randn(K,M) + 1j.*randn(K,M));theta = exp(1j.* rand(N,1).*2.*pi); %随机产生theta%初始化每个信道实现下的Hd，thetaHd_w(:,:,i) =  Hd; theta_init(:,:,i) = theta;
end
%NLOS下的G,hrk
RicFac = 10; %Rician factor
RicFac1 = sqrt(RicFac/(1+RicFac));
RicFac2 = sqrt(1/(1+RicFac));
G_sig = zeros(N,M,chanRealNum);
for i = 1: chanRealNumG_sig(:,:,i) =  sqrt(1/2).*(randn(N,M) + 1j.*randn(N,M)); %NLOS组件(CN(0,1))
end
Hr_sig = zeros(K,N,chanRealNum);% hrk转置
for i = 1: chanRealNumHr_sig(:,:,i) = sqrt(1/2).*(randn(K,N) + 1j.*randn(K,N));
end
end

搞懂了前面的矩阵维度，这里参数设置基本都没什么问题。

3、毫米波信道方向矢量的实现

function [h] = ULA_fun(phi,N)
%几何信道(一般在毫米波通信中)方向矢量
%关于方向矢量见有道云笔记中毫米波信道的记录h = exp(1j*pi*sin(phi).*(0:N-1)');
end

4、间接链路与直接链路初始化实现

    G = RicFac1.*ULA_fun(IrsArriAngle, N).*(ULA_fun(ApTranAngle, M))' + RicFac2.*G_sig;Hr = pid.*getHr(UserArriAngle,Hr_sig(:,:,t),RicFac1,RicFac2,K,N);

function [Hr] = getHr(UserArriAngle,Hr_sig,RicFac1,RicFac2,K,N)
%获取hr信道
%循环每个到达用户的到达角方向矢量Hr = zeros(K,N);for k = 1:Kuser_arri_k = UserArriAngle(k); %取到达角hr_sig = Hr_sig(k,:); %取NLOS组件hr = RicFac1.*(ULA_fun(user_arri_k, N))+RicFac2.*hr_sig; %代入公式Hr(k,:) = hr; %赋值end
end

关于毫米波通信中信道建模，可以参考这篇博客。

2020/7/13

完成：
1、学习了预编码器W的初始化部分内容
(1). γ的计算：

function [ga,gat,FA] = update_SINR(H,W,K,omega)
%计算SINR
%公式(2)ga = zeros(1,K); %存储K个用户的γ(gamma)值He = abs(H*W).^2; %He = |(hdk + hrk·Θ·G)W|^2 → (K,K)for k = 1:Ktmp = He(k,:);  %取出第k行，tmp(k)第k行第k列ga(k) = tmp(k)/(sum(tmp)-tmp(k)+1);endFA = sum(omega.*log(1+ga)); %目标函数gat = omega.*(1+ga); %更新β需要
end

其中，下面是我对He的理解：

通过以上的分析，可以很理解代码中tmp的含义：正确信号的幅度；ga(k)是γ_k；FA是f_A(W,θ)。gat是为了求解β设置的参数。

(2). β的计算：

function [beta] = update_beta(H,W,K,gat)Hc = H*W;He = abs(Hc).^2;beta = zeros(1,K);for k = 1:Ktmp = Hc(k,:);tmpe = He(k,:);beta(k) = sqrt(gat(k)) * tmp(k) / (sum(tmpe(k)) + 1); %公式(12)%注：观察公式9，αk=γk，所以公式(12)β中根号下的就是gat(k)end
end

Hc, He的理解同(1)中的图，根据公式(9)：

将α_k = γ_k代入是符合的，所以gat(k)为

2020/7/26

完成：
W的迭代更新函数(update_beam_v2)
W的求解函数(downbeam_lambda)
λ的更新函数(update_lambda_v2)
W的初始化函数(init_W)

downbeam_lambda：

function [W] = downbeam_lambda(A,H,K,M,gat,beta,lambda)
%求解WW = zeros(M,K);A = A + lambda.*eye(M); %λIA = A^(-1);for k = 1:K%wk = √(ωk(1+αk)·βk)……(在C. Prox-linear Update for W章节中)W(:,k) = sqrt(gat(k)) * beta(k) * A * H(k,:)';end
end

W的求解公式：

从公式里就能看出代码中A的含义。

update_lambda_v2

function [lambda,lambda_min] = update_lambda_v2(lambda_min,lambda_max,A,M)
%更新λ
%tmp = (lambda_min + lambda_max)/2; %初始中间值while(1)tA = A + tmp.*eye(M);if rank(tA) < Mtmp = (tmp + lambda_max)/2;  %更新中间值lambda_min = tmp; %更新最小值elsebreakendendlambda = tmp;
end

论文中关于λ的介绍是：λ > 0 is the optimal dual variable for the transmit power constraint. 所以我认为应该从这句话中得到获得λ 的方法，然鹅在论文中并不能找到如何设计λ。所以就对着源代码重新敲一遍。其中有以下不理解的地方：
①. tA含义是什么？
②. 为什么要判断rank(tA)与M的大小？

update_beam_v2

function [Wopt] = update_beam_v2(H,K,M,gat,Pt,beta,omega)
% 更新W%C.Prox-linear Update，wk = √(ωk(1+αk)·βk) 其中的求和项A = zeros(M,M);for k = 1:Ktmp = abs(beta(k))^2 * H(k,:)' * H(K,:); A = A + tmp; end% 部分参数lambda_min = 0;lambda_max = real(sum(sum(A))); %获取矩阵A所有复数和的实部flag = 0;%使优化的W功率小于或等于Ptwhile(1)%C.Prox-linear Update，wk = √(ωk(1+αk)·βk)[Wn] = downbeam_lambda(A,H,K,M,gat,beta,lambda_max); % 这一段判断是干什么用的？ 降低预编码器损耗功率？power = real(trace(Wn'*Wn));if power > Ptlambda_min = lambda_max;lambda_max = lambda_max * 2;elseif power == Ptflag = 1;break;else break;endend%W功率小于Ptif power < Ptrho = Pt/power; %ρWopt = Wn.*sqrt(rho);[~,~,f0] = update_SINR(H,Wopt,K,omega); %获取目标函数while(1)%更新λ(用的啥方法没看懂)[lambda, lambda_min] = update_lambda_v2(lambda_min, lambda_max,A,M);  %更新W[Wn] = downbeam_lambda(A,H,K,M,gat,beta,lambda);power = real(trace(Wn'*Wn));rho = Pt/power;Wn = Wn.*sqrt(rho); %优化的W[~,~,f1] = update_SINR(H,Wn,K,omega); %获取目标函数值%根据功率设置λif power > Ptlambda_min = lambda;else lambda_max = lambda;endf0 = f1;%迭代更新结束判断if abs(lambda_max-lambda_min)<1e-3 && abs(power-Pt)<1e-5 && abs(f1-f0)<1e-5breakelseif abs(lambda_max-lambda_min)<1e-7 && lambda_min == 0breakendendend
end

update_beam_v2函数是用来更新W获取初始优化值Wopt，其中不明白的一下几点：
①. 为什么要进行power与Pt的比较，判断完后对lambda设置的依据是什么。
②. λ的更新方法不清楚。
③. ρ(rho)的含义。

init_W：

function [W,gat,fint] = init_W(H,M,K,Pt,omega)
%初始化W
W = H'/(H*H'); %这是什么操作？为什么？？
W = W./sqrt(real(trace(W'* W))); %对W进行归一化for k = 1:KW(:,k) = W(:,k).*sqrt(Pt);end[~,gat,f0] = update_SINR(H,W,K,omega);  %获取γ(ga)，非log γ(gat)，目标函数f0while(1)%%更新w，beta[beta] = update_beta(H,W,K,gat);[W] = update_beam_v2(H,K,M,gat,Pt,beta,omega);[~,gat,fint] = update_SINR(H,W,K,omega);if abs(f0 - fint) < 1e-3breakendf0 = fint;end
end

其中不明白的地方：
①. W = H’ / (H * H’)的原因是什么。

IRS中波束赋形设计源代码之AO算法学习(持续更新，多多交流)相关推荐

工作中php遇到的问题以及常用函数整理(持续更新)
说明以下整理的文档是本人2017年从事php开发到目前遇到的问题的部分整理,因为上家公司有改错本这个东西,偶然间翻开,整理了一部分,后续遇到问题会持续更新,最新更新的内容会放到最前面. php开启错 ...
【中创】壹起共享“免费”网络资源库-持续更新中
想找资源又不知道在哪下载?今天整理分享16个可以免费下载的资源库, 从学习资源,电影,动漫,实用工具,大学考研,软件下载,素材资源,大学资源,网盘资源,网易,极客等...用到的各个类目的课件都能找到, ...
实习中遇到值得积累下来的编程习惯（持续更新已结束-实习结束进入正式工作）
出来实习后遇到很多在编程上的问题.有很多是在编程上的小细节,这些细节有些很关键,有些无关紧要,但是均是值得积累下来并不断坚持下去的好习惯.好方法.所以在此积累下来,并持续更新~ magic numbe ...
C++设计原则——开闭原则（持续更新中）
系列文章目录 C++开闭原则 C++迪米特法则文章目录系列文章目录前言一.开闭原则到底有什么用? 二.开闭原则有什么特点前言在一个新人刚入编程这一行的时候,可能在正式参加到一个项目之前跟着 ...
以下构成python循环结构的方法中_超星尔雅初级英语口语（持续更新中）选修课答案...
套期具有"对冲""互抵"的关系答:正确差异化战略核心是企业在市场营销.研究和开发.产品技术和工艺设计以及服务等方面具有强大的实力. 答:正确以下构成Pyt ...
Android开发学习持续更新中
Android开发单个Activity界面内的操作控件1TextView控件使用控件2Button控件使用 1首先对于android的按键格式 2对按键监听事件进行绑定控件3EditText文 ...
jtabel 遍历_Swing中经常会遇到的若干问题——JTable（持续更新） | 学步园
(1)让组件在屏幕中央显示 public static void setContainerCenter(Container container) { Dimension screenSize = To ...
Java spark中的各种范型接口Function的区别(持续更新中)
表格来自[2] Class Function Type Function<T, R> T => R DoubleFunction<T> T => Double Pa ...
javascript算法汇总(持续更新中)
1. 线性查找 1 <!doctype html> 2 <html lang="en"> 3 <head> 4 <meta charset ...

IRS中波束赋形设计源代码之AO算法学习(持续更新，多多交流)

IRS中波束赋形设计源代码之AO算法学习(持续更新，多多交流)

IRS中波束赋形设计源代码之AO算法学习(持续更新，多多交流)相关推荐

最新文章

热门文章