matlab注水算法的三种实现&与waterfill() 函数的对比

参考以下博客进行适当修改，添加了正确的容量计算，以及二分法和梯度法求解，并与matlab自带waterfill算法的对比

http://t.csdn.cn/Bmz8s

``
clear
clc
channel_n=10;
M=[1 5 10 20 50 100];%不同的总功率
N0=0.5;
h_1=random('rayleigh',1,1,channel_n);
h_2=h_1.^2;
h_2_sorted=sort(h_2);
h=h_2_sorted/N0; 
h_3=h_2/N0;
capacity=zeros(1,length(M));
delta=1e-5;
syms mu %符号变量
for m=1:length(M)
    mu2Max=1e4;
    mu2Min=1e-4;
    fprintf('》》》》》》total transmitted power : %d watt\n',M(m));
    capa=0;
    %======================使用排序+临界信道法求解=======================
    for k=1:channel_n %从差到好遍历信道起点，找到目标临界注水点
        p=zeros(1,channel_n);
        sum1=0;
        for i=k:channel_n  %假设从第k信道开始注水，每个信道按照（1/mu-1/h(i)）注水，得到sum of best powers
            sum1=sum1+(1/mu-1/h(i));% \sum(P_n^*)
        end
        f=sum1-M(m)==0; %由于该式相对1/mu线性递增，因此存在唯一解
        x=solve(f,mu); %求解方程f(mu)=0,得到当前对偶变量，后面检验功率非负条件是否满足
        if k==1 
            left=0; %从第一信道开始注水
        else
            left=h(k-1);%从第k（~=1）开始注水
        end
        if vpa(x)>left && vpa(x)<h(k) % h是经过排序的,找到注水的临界信道
                for i=k:channel_n %只有该临界以上的信道才会被分配功率
                    p(i)=1/x-1/h(i);
                    capa=capa+log2(1+p(i)*h(i));
                end
                p(i)=vpa(p(i),3);
                capa=vpa(capa,3);
                break
        end
    end
    %=====================二分法求解==========================
    while mu2Max-mu2Min>delta
         pw_bisec=zeros(1,channel_n);
         mu2=(mu2Max+mu2Min)/2;
         sum2=0;
         for k=1:channel_n
             if mu2<h(k) %保证功率非负约束
                 pw_bisec(k)=1/mu2-1/h(k);
                 sum2=sum2+pw_bisec(k);%log2(1+pw_bisec(k)*h(k));
             end
         end
         if abs(sum2-M(m))<delta
             break
         end
         if sum2>M(m)
             mu2Min=mu2;
         else 
             mu2Max=mu2;
         end
    end
%=====================梯度法求解====================
    maxiter=100000;
    mu3=2; %初始对偶变量
    ita=0.01; %固定步长
    for i=1:maxiter %总循环
        pw_grad=zeros(1,channel_n);
        sum3=0;
        for k=1:channel_n
            if mu3<h(k)
                pw_grad(k)=1/mu3-1/h(k);
                sum3=sum3+pw_grad(k);
            end
        end
        diff=M(m)-sum3;
        mu3=mu3-ita*diff;%max(,0);
        if abs(diff)<delta
            disp(['current step',num2str(i),' current diff',num2str(diff)]);
            break
        end
    end
    
    %打印功率分配结果并与matlab自带注水算法对比
    disp(['Power allocation in ',num2str(M(m)),' Watt:']);
    p
    disp('Power allocation using waterfill function inside MATLAB:');
    pw=waterfill(M(m),1./h) 
    disp('Power allocation using bisection method:');
    pw_bisec
    disp('Power allocation using gradient method:');
    pw_grad
    capacity(m)=capa;
end

上述给出了三种注水算法的实现方式：

1. 排序+找临界信道的注水求解
2. 二分法
3. 梯度法

其原理都是根据KKT条件求最佳功率值$p^{\ast }$和对偶变量$v^{\ast }$，可以参考原博文和

S.Boyd，Convex Optimization，Cambridge Press