求解Beamforming-SOCP（CVX求解）

news2024/11/18 3:20:59

时间：2023年11月23日14:00:16：
直接上代码（辛苦两天才改出来的）


clear all;
K = 4; %user number
N=4; %base station number
var=1e-9;
H = []; %initialize H matrix
for i=1:K
    h = 1/sqrt(2*K)*mvnrnd(zeros(N,1),eye(N),1)'+1i/sqrt(2*K)*mvnrnd(zeros(N,1),eye(N),1)';
    H = [H h];
end
H = H';
gamma_dB = 20; %SINR / dB
gamma = db2mag(2*gamma_dB);
gammavar=gamma*var;
POWER=40;
Kr = size (H,1); %Number of users
% N = size (H,2); %Number of transmit antennas ( in total )
D = repmat ( eye (N),[1 1 Kr ]);


cvx_begin
cvx_solver mosek
cvx_quiet ( true ); 
variable W(N, Kr ) complex ;  %多天线的话：variable W(L*N1,K) complex;
% variable POWER
minimize 0
subject to
%SINR constraints ( Kr constraints )
    for k =1: Kr
    %Channels of the signal intended for user i when it reaches user k
        hkD = zeros ( Kr ,N);
        for i =1: Kr
            hkD ( i ,:)= H( k ,:) * D (:,:, i );
        end
        imag ( hkD ( k ,:) *W(:, k ))==0; %Useful link is assumed to be real-valued
        %SOCP formulation for the SINR constraint of user k
        real ( hkD ( k ,:) *W(:, k ))>= sqrt ( gammavar ) *norm ([1 hkD( k ,:)*W(:,[1: k-1 k+1: Kr ])/sqrt(var)]);
%         real ( hkD ( k ,:) *W(:, k ))>= sqrt ( gamma) *norm ([1 hkD( k ,:)*W(:,[1: k-1 k+1: Kr ])]);
%         norm ([hkD( k ,:)*W(:,[1: k-1 k+1: Kr ])])<=sqrt(1+1/gamma)*real ( hkD ( k ,:) *W(:, k ));
    end
    for a = 1:N
        norm(W(a,:),'fro') <= POWER;
    end
%Power constraints (L constraints) scaled by the variable betavar
% norm(W,'fro') <= POWER;
% POWER >= 0; %Power constraints must be positive

cvx_end

%Analyze result and prepare the output variables.
if isempty(strfind(cvx_status,'Solved')) %Both power minimization problem and feasibility problem are infeasible.
    feasible = false;
    Wsolution = [];
else %Both power minimization problem and feasibility problem are feasible.
    feasible = true;
    Wsolution = W;
    p=norm(W,'fro');
    disp(['Power:',num2str(p)]);
    sinr=zeros(K,1);
    for k =1:K
        noise=(H( k ,:) *W(:,[1: k-1 k+1: Kr ]));
%         sinr(k)= abs(H ( k ,:) *W(:, k ))^2/(sum(abs(noise))+var);
        sinr(k)= abs(H ( k ,:) *W(:, k ))^2/((norm(noise)^2)+var);
    end
    disp("SINR of Users:");
    disp(sinr);
end

disp(["feasible:",num2str(feasible)]);