智能网络计算 实验3-物理层信道容量分析实验

一、实验目的与要求

1、解什么是凸优化问题；

2、学会使用 Matlab CVX 工具箱解决最优功率分配问题，使得信道容量最大化；

3、了解注水算法。

二、题目要求

1、最优功率分配问题

       考虑T = 10个时隙。在每个时隙 i, 发射机的发射功率为 ，发射机到接收 机的信道状态与接收机的背景噪声的比值为。假设单位带宽，则收发机之间T个时隙的总信道容量可表示为

                                                        

        考虑发射机的发射功率之和不能超过 。发射机的最优功率分配问题可表示成如下的凸优化问题：

                                        ​​​​​​​        

                                                subject to 

        使用 CVX 找出最优的功率分配，其中每一个时隙的ai可用命令a=rand(T,1)产生。给出代码与运行结果，以及最优功率分配。

2、注水算法

        步骤 1 的最优功率分配问题除了可以用 CVX 找出数值解，还可以应用凸优化数学分析的方法找出解析解。其解析解可以表示为：

        ​​​​​​​        ​​​​​​​        ​​​​​​​        ​​​​​​​        ​​​​​​​        ​​​​​​​        ​​​​​​​    ​​​​​​​

        其中，表示max(0, x)，λ是一个保证成立的常数。称为注水线。

       基于此解析解，可得出实现功率分配的注水算法

（1）思考为什么称作注水算法，注水线，信道状态 a，与功率分配的关系如何？

（2）通过程序找出λ

（3）利用 Matlab 画出结果图

         其中，横轴是时隙，纵轴是 1/a，虚线是注水线。由图可以直观地看出应该如何根据信道状态“注水”，实现最优功率分配。

三、实验过程

1、下载并安装cvx

首先需要去官网下载CVX软件包http://cvxr.com/cvx/download/，下载可能需要科学上网，否则经常下载中断

http://cvxr.com/cvx/download/

下载后将软件包解压至特定路径，在MATLAB中将当前文件位置改到这个路径，如我自己是解压到D:\file\MATLAB_CVX\cvx，

在matlab 命令行中输入命令cvx_setup

cvx_setup

回车后如下所示

输入cvx_version发现还未注册，不能使用cvx

进入cvx网址http://cvxr.com/cvx/licensing/进行申请许可证，选择学术用户进行申请

http://cvxr.com/cvx/licensing/

申请需要获取电脑主机用户名和mac地址，在cmd命令行窗口输入“set username”可以获取主机用户名，输入“

”可以获取主机地址，一般取第一个

set username

getmac

填写完毕后点击下方按钮提交

提交后可能会失败，因为邮箱不允许是私人邮箱，需要是大学提供的官方地址，要用它来验证学术资格，通过学校内部网链接查看自己学校提供的学术邮箱

重新申请后会出现如下文字，意思是文件已经发送到目的邮箱，需要登录学校提供的学生邮箱进行查收。

邮件内容如下，点击附件进行下载保存。

将下载的文件移动至先前下载cvx同一个主目录下

然后在MATLAB命令行输入cvx_setup D:\file\MATLAB_CVX/cvx_license.dat

结果如下则注册成功，在接下来的程序中可以使用cvx解决凸优化问题：

2、求解最优功率分配问题——利用cvx求解

该方法总体公式如下图所示，由于matlab中每次生成或定义的数据会被保留，也就是说在随机生成信道状态a之后，把生成a的语句注释掉，继续运行会发现a矩阵数据仍然存在不会丢失，为了验证后续的注水算法，同时为了让注水算法绘制的图比较直观，所生成的a尽量不要存在过小的数。

先将后续部分代码注释掉，对信道状态a进行生成：

点击运行后可以在右边工作区查看随机生成的10个信道状态a的情况

如下生成的两次a对比，右边明显比左边更均匀，可以重复生成直到生成较合适的a

                                    

完成生成后如下，将a生成部分的代码注释掉，保持a的一致，方便与后续的注水算法进行检验比较结果，打开cvx部分的代码注释，运行。

运行结果如下，得到的信道容量最大是1.08189

下面是各个时隙的功率分配

3、求解最优功率分配问题——注水算法

（1）注水算法1：

①总代码如下图所示

②算法思想：

直接把所有Pi求和，先不考虑其正负，求和后可以得出T×1/(ln2×λ)-∑1/ai=1，此时T还是初始值10，ai是初始的信噪比，此时遍历Pi，如果有不符合Pi≥0的，直接设对应的1/ai为0，T减1，重新按求和后的式子再求λ，直到没有Pi不符合Pi≥0为止。

相当于利用公式解出1/ln2×λ也就是注水线的值，但此时所计算的注水线是通过所有信道状态得到的，某些信道实际上“无法注水”，去掉这些无法注水的信道，通过公式关系可知注水线会随之下降，直到某一次迭代不再出现无法被注水的ai，则为最终的注水线。

这个思路的核心想法是：最后的注水线应该只和“有用”的信道有关系(ai够小)，在有可能有负数Pi的情况下，注水线会更高，如果在注水线较高的时候Pi都是负的，删掉对应的ai高的信道，注水线变低后，更不可能选择这一信道。最终会选择出合适的Pi。

③运行结果：

由图中可看出将功率分配于第2、5、7、8四个时隙，且总功率等于1.4853*4-1.20362-1.09028-1.32065-1.32674=0.99991≈1，因为有强制转换类型，会丢失一定精度，总体在误差范围内。

  

通过其与cvx的结果相比较，4个时隙cvx的功率分配分别为0.28171、0.39504、0.16467、0.15859，1.4853-1.20362=0.28168、1.4853-1.09028=0.39502、1.4853-1.32065=0.16465、1.4853-1.32674=0.15856，做表格进行更直观的比较如下：可以发现除去精度误差，两者的答案是相同的，可以得出注水算法1是成立的正确的结论。

	时隙2	时隙5	时隙7	时隙8
CVX	0.28171	0.39504	0.16467	0.15859
注水算法1	0.28168	0.39502	0.16465	0.15856

（2）注水算法2

①总代码如下图所示

②算法思想：

MATLAB有自带的函数可以求解，只要能够合适的表达出等式关系，就让MATLAB自己求解。限制Pi是非负数只需用(Pi+abs(Pi))/2。

代码51行中的func为计算总功率之和，因为功率可能会出现负数，但是负值是我们所不需要的，将10个功率以及各自绝对值相加，则负值的功率会被抵消为0，正值的功率被计算了两次，所以要除以2，并最后令func为临界值1，即被分配的总功率为1，则可以得到含x(也就是λ)的方程eqn，求解该方程可得λ的值，进而求得注水线v。

③运行结果：

可以发现与注水算法1的结果几乎一模一样，同样可以得出注水算法2是正确的结论。

  

（3）注水算法3：

①总代码如下图所示

②算法思想：

我将1/ai的信道状态当作是一个宽度一致的杯子不规则的底部，而水的量即总功率是固定的为1，十个1/ai的倒数按照从小到大的顺序排列，观察把总量为1的水倒入这个特殊的杯子能装到第几格，也就是求当水的总量为1时倒入这个杯子的高度，即注水线。划分多条水平线去计算几个矩形的面积，面积和为1。

把1/ai进行从小到大排序是为了使底部呈阶梯状，方便计算面积，每次循环以信道状态的倒数作为注水线，若此时水量小于1，说明注水线太低，需要提高；则以下一个信道状态的倒数作为注水线……，直到以第n个信道状态的倒数为注水线时，水量大于1，说明注水线太高，同时说明了注水线的范围在1/a(n-1)到1/a(n)之间。利用矩形面积去推导注水线，设注水线为x，则有如下公式，解得x即为注水线。

③运行结果：

可以发现与注水算法1和2的结果基本一致，可以得出注水算法3是正确的结论。

  

  

四、讨论与总结

1、如果生成较小的信道状态a，则可能在注水算法中出现如下的情况，a太小导致1/a太大，对比起其他信道，注水的性质难以观察清楚。

​​​​​​​

2、在写代码过程中，不知出于什么原因，每次打开matlab运行程序，涉及到CVX的调用只能运行一次，已经运行过再次运行会出现如下报错，导致我每次测试代码都要不断重复打开、关闭软件的操作，后面在测试注水算法时，为保证数据一致，将a的生成与CVX两部分代码进行注释，防止受到影响需要重新测试。

3、在编写CVX部分的代码时，需要有底数为2的对数表达式，在网络上搜索到matlab支持函数log2(x)，但是在自己实验过程中会出现如下报错。

查阅多次资料仍无法解决，最后借助换底公式完成公式的编写。

4、在编写注水算法3时，对target有错误的代入，导致计算结果时总功率不为1，最后通过matlab的步进功能，在工作区观察每次target值的变化，发现在跳出循环时，target会多执行一次加1操作，所以在判断注水线位于哪两个时隙之间时会往右边多偏移一个时隙，导致结果出现偏差，检查后将target在外面设置为-2即可正确运行出结果。

5、思考题：思考为什么称作注水算法，注水线，信道状态 a，与功率分配的关系如何？

观察绘制的结果图，可以发现本问题就像是在往一个单位底面高度（信道状态 a的倒数）凹凸不平的容器中注水一样，能够分配的总水量（总功率）是一定的，当把水倒完会有一条注水线（功率分配基准），被水覆盖到的底面所分配到的水量（功率分配）就是注水线的高度减去底面的高度，没有被水覆盖到的地面显然分配到的水量就是0。

注水算法3部分可以很直观说明注水算法的命名由来

五、总代码展示

T=10;               %10个时间间隙     
a=rand(T,1);        %每个时隙的信道状态

cvx_begin  
variable p(T);
maximize(sum(log(1+p.*a)/log(2))); %最大化目标函数，利用换底公式
subject to
    p>=0;
    sum(p)==1;
cvx_end
disp(['最优解：',num2str(p')]);


%注水算法1%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
syms x
N=T;
inverse_a=1./a;
iterlamda=N/(1+sum(inverse_a));         %设一下初值，iter_lamda就是不断在更新的lamda值
for j=1:10                              %把所有不符合条件的Pi都剔除
    for i=1:10                          %检测有没有Pi不符合条件，有的话删除掉再计算参数lamda
        if(inverse_a(i)~=0&&1/iterlamda<inverse_a(i))%1/ai不等于0，说明信道还没被删掉，
                                                     %而1/lamda小于1/ai，说明该信道不该被选
            inverse_a(i)=0;
            N=N-1;        %1/ai设为9，N--
            iterlamda=N/(1+sum(inverse_a));%再次计算
        end
    end
end

lamda=iterlamda/(log(2));
v=1/((log(2)/log(exp(1)))*lamda);%计算注水线值
z=[];
for i=0:T-1
    y = 1/a(i+1);
    z=[z;i y;i+1 y];
end
figure(1);
plot(z(:,1),z(:,2));
line([0 T],[v v],'linestyle',':')
xlabel('Channel index a(i)');
legend('1/a(i)','注水线');
set(gca,'xtick',[],'ytick',[]);
text(-1.2,v,num2str(v));


%注水算法2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
syms x
func=0;
for i=1:10
    func=func+1/2*((1/(log(2)*x)-1/a(i))+abs((1/(log(2)*x)-1/a(i))));
    %求和，限制pi是非负数
end
eqn=func==1;
lamda=double(vpasolve(eqn,x));

v=1/((log(2)/log(exp(1)))*lamda);%计算注水线
z=[];
for i=0:T-1
    y=1/a(i+1);
    z=[z;i y;i+1 y];
end

figure(2);
plot(z(:,1),z(:,2));
line([0 T],[v v],'linestyle',':');
xlabel('Channel index a(i)');
legend('1/a','注水线');
set(gca,'xtick',[],'ytick',[]);
text(-1.2,v,num2str(v));


%注水算法3%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%首先对a进行处理，得到取倒数后的b数组，即注水的底，并从小到大进行排序
b=1/a;
for i=1:T
    b(i)=1/a(i);
end
b=sort(b);
sum=0;
target=1;

%找到注水线所在的大致范围
while(sum<1)
    sum=0;
    t=b(target);
    for i=1:target
        sum=sum+t-b(i);
    end
    target=target+1;
end

target=target-2;
syms x
sumb=0;
for i=1:target
    sumb=sumb+b(i);
end
enq=x*target-sumb==1;
x=solve(enq);
ans=double(x);

z=[];
for i=0:T-1
    y=1/a(i+1);%y=b(i+1)
    z=[z;i y;i+1 y];
end

figure(3);
plot(z(:,1),z(:,2));
line([0 T],[ans ans],'linestyle',':');
xlabel('Channel index a(i)');
legend('1/a','注水线');
set(gca,'xtick',[],'ytick',[]);
text(-1.2,ans,num2str(ans));