CVX-Invalid quadratic form(s): not a square - CVX Forum: a community-driven support forum

L

Luchu Li

Sep 11, 2023 #1

How to fix it?? I’m confused

CODE:

X1 = zeros(K,K);
X2 = zeros(K,K);

% Beta0 = zeros(K,1); %beta’
% gamma = 0;
cvx_begin
variables Beta0(K) g
variable X1(K,K) complex
variable X2(K,K) complex
maximize ( g )
subject to

        for k=1:1:K
            Beta0(k)>=0;
        end
        sum(Beta0)<=1;
        
        for k=1:1:K
            for i=1:1:K
                X1(k,i) == Rho_d*Beta0(i)*trace(C_ek_ek(:,:,k)*C_vkk_vkk(:,:,i))/trace(C_vkk_vkk(:,:,i));
                if k~=i
                    X2(k,i) == Rho_d*Beta0(i)*I(k,i)/trace(C_vkk_vkk(:,:,i));
                else
                    X2(k,i) == 0;
                end
            end
        end
        
        for k=1:1:K
            Rho_d*Beta0(k)*trace(C_vkk_vkk(:,:,k))>=(Rho_d*Beta0(k)*(I(k,k)- (trace(C_vkk_vkk(:,:,k)))^2)/trace(C_vkk_vkk(:,:,k))...
                +sum(X1(k,:))+sum(X2(k,:))+1)*g;
        end 

cvx_end

for k=1:1:K
    Beta(k) = Beta0(k)/trace(C_vkk_vkk(:,:,k));
end

error .* (line 262)
Disciplined convex programming error:
Invalid quadratic form(s): not a square.

error * (line 36)
z = feval( oper, x, y );

error main (line 317)
Rho_dBeta0(k)trace(C_vkk_vkk(:,:,k))>=(Rho_dBeta0(k)(I(k,k)-(trace(C_vkk_vkk(:,:,k)))^2)/trace(C_vkk_vkk(:,:,k))…

I don’t know why I can’t build it successfully. Hope somebody can help me. It’s painful…

L

Luchu Li

Sep 11, 2023 #2

the ordinary question:

M

Mark L. Stone

Sep 11, 2023 #3

I’m having trouble figuring out which parentheses go with which. Are you multiplying variables? That is not allowed in CVX unless gp mode is specified, using cvx_begin gp. That plus all the rules of Geometric Programming are discussed in Geometric programming mode — CVX Users' Guide .

Did you notice the mention of geometric programming in the image you posted? That should give you a good hint that you need to specify gp mode, and follow gp mode’s rules.

L

Luchu Li

Sep 12, 2023 #4

Replying to #3

I add the cvx_begin gp , but it still have some problems…

In main (line 293)
error cvxprob/newobj (line 57)
Disciplined convex programming error:
Cannot maximize a(n) convex expression.

error maximize (line 21)
newobj( prob, ‘maximize’, x );

error main (line 297)
maximize ( g )

M

Mark L. Stone

Sep 12, 2023 #5

Did you read the link?

L

Luchu Li

Sep 12, 2023 #6

why there is no warn in the original one? The question need “maximize (g)”, and I think g is affine which can be maximized

M

Mark L. Stone

Sep 12, 2023 #7

Please sow a complete reproducible problem, with all input data. Then show all the CVX output.

L

Luchu Li

Sep 12, 2023 #8

The figure above presents the complete problem:

The variable g is an auxiliary variable that represents the minimum user rate.
The variables X1 and X2 are intermediate variables that do not have practical meaning and are used for auxiliary calculations.
The variables I(k,i), C_vkk_vkk(:,:,k), and C_ek_ek(:,:,k) can be treated as constant complex matrices.
Rho_d is a constant parameter.
Beta0(k) is a crucial variable that represents the user rate allocation coefficients.

CVX1559×621 41.9 KB

M

Mark L. Stone

Sep 12, 2023 #9

That is not code which generated an error message/.

L

Luchu Li

Sep 13, 2023 #10

clc
clear

%% parameter define
global M M1 M2 N N1 N2 K B noise NF Tau_c Tau_u lambda

P=1;
N=1600;N1=40;N2=40;
M=900;M1=30;M2=30;
K=5;

noise=3.98e-14; %-104dbm
NF=5.012;
B=10e+6;
P_noise=10log10(noiseB*NF);

Tau_c=196; %相干间隔
Tau_u=K; %导频长度
d_irs_user=[70, 50, 50, 80, 60]; %RIS-user distance
d_BS_user=[80, 80, 60, 60, 70]; %BS-user distance
d_BS_ris=50; %BS-RIS间距

lambda = 0.3; %假设电磁波频率为1GHz
d_az = lambda;
d_el = lambda; %假设天线间距都等于一个波长
d = lambda/4;
d_H = d;d_V = d;A = d*d; %RIS元件的间距/面积

zeta_0=1.58e+4; %路径损耗 zeta_x = zeta_0 + 10k0log (d);
k0=3.5;
zeta_b = zeta_0 + 10k0log (d_BS_ris); %BS-user路径损耗
zeta_f = zeros(1,K); %RIS-user路径损耗
zeta_h = zeros(1,K); %BS-user路径损耗
for k=1:1:K
zeta_f(k) = zeta_0 + 10k0log (d_irs_user(k));
zeta_h(k) = zeta_0 + 10k0log (d_BS_user(k));
end

%% variable Setup
variable1 = 1:1:70;
x_axis = variable1; %横坐标——物理意义为rho_d(对数)
rate = zeros(1,length(variable1)); %纵坐标——物理意义为最小用户速率

for mm=1:1:length(variable1) %对于每个确定的横坐标值variable(mm)，去确定纵坐标的值
%% Update the variable %%更新横坐标值

Rho_d=10^(variable1(mm)/10);

%% phase shifts/transmit power allocation – definition，Random
%phi = exp(2pilirand(1,N)); %变量-相位向量
Beta = zeros(K,1); %beta(列向量)
% phi = 2pi*rand(1,N);
% for k=1:1:K
% Beta(k) = 1/K;
% end

%% 计算信道协方差矩阵

% R_r=zeros(N,N);
% R_b=zeros(M,M);
R_tl_az=zeros(M2,M2);
R_tl_el=zeros(M1,M1);

%BS部分——
psi_l = 2*pi/8;
omega_l = pi/3;
alpha_l = pi/12;
nu_l = pi/36;

for u=1:1:M2
    for v=1:1:M2
        delta1 = 2*pi*d_az*(v-u)*sin(psi_l);
        delta2 = 2*pi*d_az*(v-u)*nu_l*cos(psi_l);
        delta3 = (delta2*alpha_l*sin(omega_l))^2+1;
        R_tl_az(u,v) = ((delta3)^(-0.5))*exp(-(delta2*cos(omega_l))^2/(2*delta3)-((delta1*alpha_l*sin(omega_l))^2)/(2*delta3))...
            *exp(1i*delta1*cos(omega_l)/delta3);
    end
end
for x=1:1:M1
    for y=1:1:M1
        R_tl_el(x,y) = exp(1i*2*pi*d_el*(y-x)*cos(psi_l))*exp(-0.5*((2*pi*nu_l*d_el)^2)*((y-x)^2)*((sin(psi_l))^2));
    end
end
R_tl = kron(R_tl_az,R_tl_el);

%RIS部分——
R = zeros(N,N);
u = zeros(3,N);
for n=1:1:N
u(1,n) = 0;u(2,n) = mod(n-1,N2)d_H;u(3,n) = d_Vfloor((n-1)/N2);
end
for p=1:1:N
for q=1:1:N
R(p,q) = sinc(2*norm((u(:,p)-u(:,q)),2)/lambda);
end
end

%计算R_b:
R_b = zeta_b*R_tl;
    
%计算R_r:
R_r = A*R;

%计算R_h:BS的VR区域——1:(1-10)*30,2:(6-15)*30,3:(11-20)*30,4:(16-25)*30,5:(21-30)*30
R_h=zeros(M,M,K);


dx = zeros(M2,K);
dy = ones(M1,K);
Dx = zeros(M2,M2,K);
Dy = zeros(M1,M1,K);
for n=1:1:M2
    if (1<=n)&&(n<=10)
        dx(n,1)=1;
    end
    if (6<=n)&&(n<=15)
        dx(n,2)=1;
    end
    if (11<=n)&&(n<=20)
        dx(n,3)=1;
    end
    if (16<=n)&&(n<=25)
        dx(n,4)=1;
    end
    if (21<=n)&&(n<=30)
        dx(n,5)=1;
    end
end
for k=1:1:K
    Dx(:,:,k) = diag(dx(:,k));
    Dy(:,:,k) = diag(dy(:,k));
    
    R_h(:,:,k) = kron(((Dx(:,:,k))^0.5)*R_tl_az*((Dx(:,:,k))^0.5),((Dy(:,:,k))^0.5)*R_tl_el*((Dy(:,:,k))^0.5));
end


%计算R_f:RIS的VR区域——1:(1-8)*40,2:(9-16)*40,3:(17-24)*40,4:(25-32)*40,5:(33-40)*40
R_f=zeros(N,N,K);
R_f_0=zeros(N,N,K);
for k=1:1:K
    R_f_0(:,:,k) = zeta_f(k)*A*R;
end

for p=1:1:N
    for q=1:1:N
        if (1<=p)&&(p<=320)&&(1<=q)&&(q<=320)
            R_f(p,q,1)=R_f_0(p,q,1);
        end
        if (321<=p)&&(p<=640)&&(321<=q)&&(q<=640)
            R_f(p,q,2)=R_f_0(p,q,2);
        end
        if (641<=p)&&(p<=960)&&(641<=q)&&(q<=960)
            R_f(p,q,3)=R_f_0(p,q,3);
        end
        if (961<=p)&&(p<=1280)&&(961<=q)&&(q<=1280)
            R_f(p,q,4)=R_f_0(p,q,4);
        end
        if (1281<=p)&&(p<=1600)&&(1281<=q)&&(q<=1600)
            R_f(p,q,5)=R_f_0(p,q,5);
        end
    end
end

%% CVX_Phase shift optimization

% phi_t = exp(2pi1irand(1,N)); % ?θt
% phi_t_1 = exp(2pi1irand(1,N)); %?θ?t
% %phi_best = zeros(1,N);%最佳相位
%
% t=0;
% Nu_t = rand;
% Nu = rand;
% %Ak = zeros(N,N,K);
%
% while 1
% cvx_begin
% variable phi(1,N) complex % ?θt
% variable Nu
% maximize ( Nu )
% subject to
% for n=1:1:N
% square_abs(phi(n))<=1;
% end
% for k=1:1:K
% %Ak(:,:,k)=R_r.R_f(:,:,k);
% real(trace(R_h(:,:,k))+trace(R_b)(phi_t(R_r.R_f(:,:,k))(phi_t’)+2real((phi-phi_t)(R_r.R_f(:,:,k))(phi_t’))))>=Nu;
% end
% cvx_end
%
%
% for n=1:1:N
% if phi(n)==0
% phi_t_1(n) = 0;
% else
% phi_t_1(n) = phi(n)/(abs(phi(n)));
% end
% end
% phi_t = phi_t_1;
%
% if abs(Nu-Nu_t)/abs(Nu_t)<0.005 %如果收敛，则跳出
% break
% end
% Nu_t = Nu;
% t = t + 1;
% end
%
% phi_best = phi_t;
% Phi = diag(phi_best); %使用最佳相位优化的结果

% phiAk(phi’)>=(phi_tAk(phi_t’)+2*real((phi-phi_t)Ak(phi_t’))))

%% CVX_Transmit power allocation optimization
Phi = diag(exp(2*pi*1i*rand(1,N)));


C_yk_vk = zeros(M,M,K);    
C_yk_yk = zeros(M,M,K);    
C_vkk_vkk = zeros(M,M,K);    
C_ek_ek = zeros(M,M,K);
for k=1:1:K 
    C_yk_vk(:,:,k) = R_h(:,:,k)+R_b*trace(R_r*Phi*R_f(:,:,k)*(Phi'));%M*M
    C_yk_yk(:,:,k) = C_yk_vk(:,:,k)+P_noise*eye(M);%M*M
    C_vkk_vkk(:,:,k) = C_yk_vk(:,:,k)*((C_yk_yk(:,:,k))^(-1))*C_yk_vk(:,:,k);%M*M
    C_ek_ek(:,:,k) = C_yk_vk(:,:,k)-C_vkk_vkk(:,:,k);%M*M
end

S = zeros(N,N,K);
for k=1:1:K
    S(:,:,k) = R_f(:,:,k)*(Phi')*R_r*Phi;
end

C = zeros(M,M,K,K);
T = zeros(M,M,K,K);
for k=1:1:K
    for i=1:1:K
        C(:,:,k,i) = ((C_yk_yk(:,:,k))^(-1))*C_yk_vk(:,:,k)*C_yk_vk(:,:,i)*((C_yk_yk(:,:,i))^(-1));
        T(:,:,k,i) = R_b*C(:,:,k,i);
    end
end


Xi = P_noise;%噪声方差

I = zeros(K,K);
T_D_Sk = zeros(1,K);
T_D_Uk = zeros(1,K);
T_IU_Iki = zeros(K,K);
for k=1:1:K
    for i=1:1:K
        if i==k
            Tk = T(:,:,k,i);
            Sk = S(:,:,k);
            Ck = C(:,:,k,i);
            R_hk = R_h(:,:,k);
            I(k,i) = 2*trace(Sk)*trace(Tk)*trace(R_hk*Ck)+(trace(Sk)*trace(Tk))^2+...
                       trace(Sk^2)*trace(Tk^2)+((trace(Sk))^2)*trace(Tk^2)+...
                       ((trace(Tk))^2)*trace(Sk^2)+2*Xi*trace(Sk)*trace(Tk)*trace(Ck)+...
                       2*trace(Sk)*trace(Tk*R_hk*Ck)+2*Xi*trace(Sk)*trace(Tk*Ck)+...
                       2*Xi*trace(R_hk*Ck)*trace(Ck)+2*Xi*trace(R_hk*(Ck^2))+(abs(trace(R_hk*Ck)))^2+...
                       trace((R_hk*Ck)^2)+Xi*Xi*((abs(trace(Ck)))^2)+Xi*Xi*trace(Ck^2);
            
            T_D_Sk(k) = Rho_d*Beta(k)*((trace(C_vkk_vkk(:,:,k)))^2);
            T_D_Uk(k) = Rho_d*Beta(k)*(trace(C_ek_ek(:,:,k)*C_vkk_vkk(:,:,k))-(trace(C_vkk_vkk(:,:,k)))^2+I(k,i));
         
        else
            Tki = T(:,:,k,i);
            Sk = S(:,:,k);
            Si = S(:,:,i);
            Cki = C(:,:,k,i);
            Qki = R_b*(Cki');
            R_hk = R_h(:,:,k);
            R_hi = R_h(:,:,i);
            I(k,i) = trace(Sk*Si)*trace(Tki)*trace(Qki)+trace(Sk)*trace(Tki*R_hi*(Cki'))+...
                        Xi*trace(Sk)*trace(Tki*(Cki'))+Xi*Xi*trace(Cki*(Cki'))+trace(R_hk*Cki*R_hi*(Cki'))+...
                        trace(Sk)*trace(Si)*trace(Tki*Qki)+Xi*trace(Cki*R_hi*(Cki'))+Xi*trace(R_hk*Cki*(Cki'))+...
                        trace(Si)*trace(R_hk*Cki*Qki)+Xi*trace(Cki*Qki)*trace(Si);
                    
            T_IU_Iki(k,i) = Rho_d*Beta(i)*(trace(C_ek_ek(:,:,k)*C_vkk_vkk(:,:,i))+I(k,i));
        end
    end
end


X1 = zeros(K,K);
X2 = zeros(K,K);

% Beta0 = zeros(K,1); %beta’
% gamma = 0;
% cvx_begin
% variables Beta0(K) g
% variable X1(K,K) complex
% variable X2(K,K) complex
% maximize ( g )
% subject to
%
% for k=1:1:K
% Beta0(k)>=0;
% end
% sum(Beta0)<=1;
%
% for k=1:1:K
% for i=1:1:K
% X1(k,i) == Rho_d*Beta0(i)trace(C_ek_ek(:,:,k)C_vkk_vkk(:,:,i))/trace(C_vkk_vkk(:,:,i));
% if k~=i
% X2(k,i) == Rho_dBeta0(i)I(k,i)/trace(C_vkk_vkk(:,:,i));
% else
% X2(k,i) == 0;
% end
% end
% end
%
% for k=1:1:K
% Rho_dBeta0(k)trace(C_vkk_vkk(:,:,k))>=(Rho_dBeta0(k)(I(k,k)-(trace(C_vkk_vkk(:,:,k)))^2)/trace(C_vkk_vkk(:,:,k))…
% +sum(X1(k,:))+sum(X2(k,:))+1)*g;
% end
%
% cvx_end
cvx_begin
variables Beta0(K) g
variables X1r(K,K) X1i(K,K)
variables X2r(K,K) X2i(K,K)
variables Z1(K,K) Z2(K,K)
maximize(g)
subject to
for k = 1:K
Beta0(k) >= 0;
end
sum(Beta0) <= 1;

    for k = 1:K
        for i = 1:K
            X1r(k,i) >= 0;
            X1i(k,i) >= 0;
            X1r(k,i) == real(Rho_d * Beta0(i) * trace(C_ek_ek(:,:,k) * C_vkk_vkk(:,:,i)) / trace(C_vkk_vkk(:,:,i)));
            X1i(k,i) == imag(Rho_d * Beta0(i) * trace(C_ek_ek(:,:,k) * C_vkk_vkk(:,:,i)) / trace(C_vkk_vkk(:,:,i)));
            if k ~= i
                X2r(k,i) >= 0;
                X2i(k,i) >= 0;
                X2r(k,i) == real(Rho_d * Beta0(i) * I(k,i) / trace(C_vkk_vkk(:,:,i)));
                X2i(k,i) == imag(Rho_d * Beta0(i) * I(k,i) / trace(C_vkk_vkk(:,:,i)));
            else
                X2r(k,i) == 0;
                X2i(k,i) == 0;
            end
        end
    end

    for k = 1:K
        Z1(k,k) >= 0;
        Z2(k,k) >= 0;
        norm(vec([X1r(k,:); X1i(k,:)]), 2) <= Z1(k,k);
        norm(vec([X2r(k,:); X2i(k,:)]), 2) <= Z2(k,k);
        Rho_d * Beta0(k) * trace(C_vkk_vkk(:,:,k)) >= quad_over_lin([Z1(k,k); Z2(k,k)], 1)...
            - Rho_d * Beta0(k) * (I(k,k) - square_abs(trace(C_vkk_vkk(:,:,k))) / trace(C_vkk_vkk(:,:,k))) + g;
    end
cvx_end

for k=1:1:K
    Beta(k) = Beta0(k)/trace(C_vkk_vkk(:,:,k));
end
rate(mm) = gamma;

%% 计算用户速率
% Rate = zeros(1,K);
% for k=1:1:K
% Rate(k) = ((Tau_c-Tau_u)/Tau_c)*log2(1+T_D_Sk(k)/(T_D_Uk(k)+sum(T_IU_Iki(k,:)+1)));
% end
% Rate_a = sum(Rate)/K;

end

figure(1) %画出图2
plot(x_axis,rate,‘-p’,‘LineWidth’,1.5,‘MarkerSize’,8) %画出理论速率

M

Mark L. Stone

Sep 13, 2023 #11

That;s not gp mode. You said you were using it.

Ir would also make things easier for forum readers if you eliminated the commented out code from your posted code.

L

Luchu Li

Sep 13, 2023 #12

clc
clear

%% parameter define
global M M1 M2 N N1 N2 K B noise NF Tau_c Tau_u lambda

P=1;
N=1600;N1=40;N2=40;
M=900;M1=30;M2=30;
K=5;

noise=3.98e-14; %-104dbm
NF=5.012;
B=10e+6;
P_noise=10log10(noiseB*NF);

Tau_c=196; %相干间隔
Tau_u=K; %导频长度
d_irs_user=[70, 50, 50, 80, 60]; %RIS-user distance
d_BS_user=[80, 80, 60, 60, 70]; %BS-user distance
d_BS_ris=50; %BS-RIS间距

lambda = 0.3; %假设电磁波频率为1GHz
d_az = lambda;
d_el = lambda; %假设天线间距都等于一个波长
d = lambda/4;
d_H = d;d_V = d;A = d*d; %RIS元件的间距/面积

zeta_0=1.58e+4; %路径损耗 zeta_x = zeta_0 + 10k0log (d);
k0=3.5;
zeta_b = zeta_0 + 10k0log (d_BS_ris); %BS-user路径损耗
zeta_f = zeros(1,K); %RIS-user路径损耗
zeta_h = zeros(1,K); %BS-user路径损耗
for k=1:1:K
zeta_f(k) = zeta_0 + 10k0log (d_irs_user(k));
zeta_h(k) = zeta_0 + 10k0log (d_BS_user(k));
end

%% variable Setup
variable1 = 1:1:70;
x_axis = variable1; %横坐标——物理意义为rho_d(对数)
rate = zeros(1,length(variable1)); %纵坐标——物理意义为最小用户速率

for mm=1:1:length(variable1) %对于每个确定的横坐标值variable(mm)，去确定纵坐标的值
%% Update the variable %%更新横坐标值

Rho_d=10^(variable1(mm)/10);

%% phase shifts/transmit power allocation – definition，Random
%phi = exp(2pilirand(1,N)); %变量-相位向量
Beta = zeros(K,1); %beta(列向量)
% phi = 2pi*rand(1,N);
% for k=1:1:K
% Beta(k) = 1/K;
% end

%% 计算信道协方差矩阵

% R_r=zeros(N,N);
% R_b=zeros(M,M);
R_tl_az=zeros(M2,M2);
R_tl_el=zeros(M1,M1);

%BS部分——
psi_l = 2*pi/8;
omega_l = pi/3;
alpha_l = pi/12;
nu_l = pi/36;

for u=1:1:M2
    for v=1:1:M2
        delta1 = 2*pi*d_az*(v-u)*sin(psi_l);
        delta2 = 2*pi*d_az*(v-u)*nu_l*cos(psi_l);
        delta3 = (delta2*alpha_l*sin(omega_l))^2+1;
        R_tl_az(u,v) = ((delta3)^(-0.5))*exp(-(delta2*cos(omega_l))^2/(2*delta3)-((delta1*alpha_l*sin(omega_l))^2)/(2*delta3))...
            *exp(1i*delta1*cos(omega_l)/delta3);
    end
end
for x=1:1:M1
    for y=1:1:M1
        R_tl_el(x,y) = exp(1i*2*pi*d_el*(y-x)*cos(psi_l))*exp(-0.5*((2*pi*nu_l*d_el)^2)*((y-x)^2)*((sin(psi_l))^2));
    end
end
R_tl = kron(R_tl_az,R_tl_el);

%RIS部分——
R = zeros(N,N);
u = zeros(3,N);
for n=1:1:N
u(1,n) = 0;u(2,n) = mod(n-1,N2)d_H;u(3,n) = d_Vfloor((n-1)/N2);
end
for p=1:1:N
for q=1:1:N
R(p,q) = sinc(2*norm((u(:,p)-u(:,q)),2)/lambda);
end
end

%计算R_b:
R_b = zeta_b*R_tl;
    
%计算R_r:
R_r = A*R;

%计算R_h:BS的VR区域——1:(1-10)*30,2:(6-15)*30,3:(11-20)*30,4:(16-25)*30,5:(21-30)*30
R_h=zeros(M,M,K);


dx = zeros(M2,K);
dy = ones(M1,K);
Dx = zeros(M2,M2,K);
Dy = zeros(M1,M1,K);
for n=1:1:M2
    if (1<=n)&&(n<=10)
        dx(n,1)=1;
    end
    if (6<=n)&&(n<=15)
        dx(n,2)=1;
    end
    if (11<=n)&&(n<=20)
        dx(n,3)=1;
    end
    if (16<=n)&&(n<=25)
        dx(n,4)=1;
    end
    if (21<=n)&&(n<=30)
        dx(n,5)=1;
    end
end
for k=1:1:K
    Dx(:,:,k) = diag(dx(:,k));
    Dy(:,:,k) = diag(dy(:,k));
    
    R_h(:,:,k) = kron(((Dx(:,:,k))^0.5)*R_tl_az*((Dx(:,:,k))^0.5),((Dy(:,:,k))^0.5)*R_tl_el*((Dy(:,:,k))^0.5));
end


%计算R_f:RIS的VR区域——1:(1-8)*40,2:(9-16)*40,3:(17-24)*40,4:(25-32)*40,5:(33-40)*40
R_f=zeros(N,N,K);
R_f_0=zeros(N,N,K);
for k=1:1:K
    R_f_0(:,:,k) = zeta_f(k)*A*R;
end

for p=1:1:N
    for q=1:1:N
        if (1<=p)&&(p<=320)&&(1<=q)&&(q<=320)
            R_f(p,q,1)=R_f_0(p,q,1);
        end
        if (321<=p)&&(p<=640)&&(321<=q)&&(q<=640)
            R_f(p,q,2)=R_f_0(p,q,2);
        end
        if (641<=p)&&(p<=960)&&(641<=q)&&(q<=960)
            R_f(p,q,3)=R_f_0(p,q,3);
        end
        if (961<=p)&&(p<=1280)&&(961<=q)&&(q<=1280)
            R_f(p,q,4)=R_f_0(p,q,4);
        end
        if (1281<=p)&&(p<=1600)&&(1281<=q)&&(q<=1600)
            R_f(p,q,5)=R_f_0(p,q,5);
        end
    end
end

%% CVX_Phase shift optimization

phi_t = exp(2*pi*1i*rand(1,N));   % ?θt
phi_t_1 = exp(2*pi*1i*rand(1,N));   %?θ?t
%phi_best = zeros(1,N);%最佳相位

t=0;
Nu_t = rand;
Nu = rand;
%Ak = zeros(N,N,K);

while 1   
    cvx_begin
        variable phi(1,N) complex  % ?θ*t
        variable Nu 
        maximize ( Nu )
        subject to
            for n=1:1:N
                square_abs(phi(n))<=1;
            end
            for k=1:1:K
                %Ak(:,:,k)=R_r.*R_f(:,:,k);
                real(trace(R_h(:,:,k))+trace(R_b)*(phi_t*(R_r.*R_f(:,:,k))*(phi_t')+2*real((phi-phi_t)*(R_r.*R_f(:,:,k))*(phi_t'))))>=Nu;
            end
    cvx_end
          
    
    for n=1:1:N
        if phi(n)==0
            phi_t_1(n) = 0;
        else
            phi_t_1(n) = phi(n)/(abs(phi(n)));
        end
    end
    phi_t = phi_t_1;
    
    if abs(Nu-Nu_t)/abs(Nu_t)<0.005   %如果收敛，则跳出
       break
    end
    Nu_t = Nu;
    t = t + 1;
end

phi_best = phi_t;
Phi = diag(phi_best);   %使用最佳相位优化的结果




%% CVX_Transmit power allocation optimization



C_yk_vk = zeros(M,M,K);    
C_yk_yk = zeros(M,M,K);    
C_vkk_vkk = zeros(M,M,K);    
C_ek_ek = zeros(M,M,K);
for k=1:1:K 
    C_yk_vk(:,:,k) = R_h(:,:,k)+R_b*trace(R_r*Phi*R_f(:,:,k)*(Phi'));%M*M
    C_yk_yk(:,:,k) = C_yk_vk(:,:,k)+P_noise*eye(M);%M*M
    C_vkk_vkk(:,:,k) = C_yk_vk(:,:,k)*((C_yk_yk(:,:,k))^(-1))*C_yk_vk(:,:,k);%M*M
    C_ek_ek(:,:,k) = C_yk_vk(:,:,k)-C_vkk_vkk(:,:,k);%M*M
end

S = zeros(N,N,K);
for k=1:1:K
    S(:,:,k) = R_f(:,:,k)*(Phi')*R_r*Phi;
end

C = zeros(M,M,K,K);
T = zeros(M,M,K,K);
for k=1:1:K
    for i=1:1:K
        C(:,:,k,i) = ((C_yk_yk(:,:,k))^(-1))*C_yk_vk(:,:,k)*C_yk_vk(:,:,i)*((C_yk_yk(:,:,i))^(-1));
        T(:,:,k,i) = R_b*C(:,:,k,i);
    end
end


Xi = P_noise;%噪声方差

I = zeros(K,K);
T_D_Sk = zeros(1,K);
T_D_Uk = zeros(1,K);
T_IU_Iki = zeros(K,K);
for k=1:1:K
    for i=1:1:K
        if i==k
            Tk = T(:,:,k,i);
            Sk = S(:,:,k);
            Ck = C(:,:,k,i);
            R_hk = R_h(:,:,k);
            I(k,i) = 2*trace(Sk)*trace(Tk)*trace(R_hk*Ck)+(trace(Sk)*trace(Tk))^2+...
                       trace(Sk^2)*trace(Tk^2)+((trace(Sk))^2)*trace(Tk^2)+...
                       ((trace(Tk))^2)*trace(Sk^2)+2*Xi*trace(Sk)*trace(Tk)*trace(Ck)+...
                       2*trace(Sk)*trace(Tk*R_hk*Ck)+2*Xi*trace(Sk)*trace(Tk*Ck)+...
                       2*Xi*trace(R_hk*Ck)*trace(Ck)+2*Xi*trace(R_hk*(Ck^2))+(abs(trace(R_hk*Ck)))^2+...
                       trace((R_hk*Ck)^2)+Xi*Xi*((abs(trace(Ck)))^2)+Xi*Xi*trace(Ck^2);
            
            T_D_Sk(k) = Rho_d*Beta(k)*((trace(C_vkk_vkk(:,:,k)))^2);
            T_D_Uk(k) = Rho_d*Beta(k)*(trace(C_ek_ek(:,:,k)*C_vkk_vkk(:,:,k))-(trace(C_vkk_vkk(:,:,k)))^2+I(k,i));
         
        else
            Tki = T(:,:,k,i);
            Sk = S(:,:,k);
            Si = S(:,:,i);
            Cki = C(:,:,k,i);
            Qki = R_b*(Cki');
            R_hk = R_h(:,:,k);
            R_hi = R_h(:,:,i);
            I(k,i) = trace(Sk*Si)*trace(Tki)*trace(Qki)+trace(Sk)*trace(Tki*R_hi*(Cki'))+...
                        Xi*trace(Sk)*trace(Tki*(Cki'))+Xi*Xi*trace(Cki*(Cki'))+trace(R_hk*Cki*R_hi*(Cki'))+...
                        trace(Sk)*trace(Si)*trace(Tki*Qki)+Xi*trace(Cki*R_hi*(Cki'))+Xi*trace(R_hk*Cki*(Cki'))+...
                        trace(Si)*trace(R_hk*Cki*Qki)+Xi*trace(Cki*Qki)*trace(Si);
                    
            T_IU_Iki(k,i) = Rho_d*Beta(i)*(trace(C_ek_ek(:,:,k)*C_vkk_vkk(:,:,i))+I(k,i));
        end
    end
end


X1 = zeros(K,K);
X2 = zeros(K,K);
Beta0 = zeros(K,1);   %beta'
gamma = 0;
cvx_begin 
    variables Beta0(K) g
    variable X1(K,K) complex
    variable X2(K,K) complex
    maximize ( g )
    subject to
         
        for k=1:1:K
            Beta0(k)>=0;
        end
        sum(Beta0)<=1;
        
        for k=1:1:K
            for i=1:1:K
                X1(k,i) == Rho_d*Beta0(i)*trace(C_ek_ek(:,:,k)*C_vkk_vkk(:,:,i))/trace(C_vkk_vkk(:,:,i));
                if k~=i
                    X2(k,i) == Rho_d*Beta0(i)*I(k,i)/trace(C_vkk_vkk(:,:,i));
                else
                    X2(k,i) == 0;
                end
            end
        end
        
        for k=1:1:K
            Rho_d*Beta0(k)*trace(C_vkk_vkk(:,:,k))>=(Rho_d*Beta0(k)*(I(k,k)-(trace(C_vkk_vkk(:,:,k)))^2)/trace(C_vkk_vkk(:,:,k))...
                +sum(X1(k,:))+sum(X2(k,:))+1)*g;
        end 

cvx_end


for k=1:1:K
    Beta(k) = Beta0(k)/trace(C_vkk_vkk(:,:,k));
end
rate(mm) = gamma;

%% 计算用户速率
% Rate = zeros(1,K);
% for k=1:1:K
% Rate(k) = ((Tau_c-Tau_u)/Tau_c)*log2(1+T_D_Sk(k)/(T_D_Uk(k)+sum(T_IU_Iki(k,:)+1)));
% end
% Rate_a = sum(Rate)/K;

end

figure(1) %画出图2
plot(x_axis,rate,‘-p’,‘LineWidth’,1.5,‘MarkerSize’,8) %画出理论速率

L

Luchu Li

Sep 13, 2023 #13

Replying to #11

I want to use any possible mode to solve my problem, so I’m still debuging

CODE:

I add the cvx_begin gp , but it still have some problems…

error main (line 297) maximize ( g )

error main (line 297)
maximize ( g )