# CVX-Invalid quadratic form(s): not a square

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…

the ordinary question:

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.

## 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 )

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

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

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.

That is not code which generated an error message/.

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 = 2
pi*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(2
pi1irand(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_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_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) %画出理论速率

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.

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 = 2
pi*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) %画出理论速率

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