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

luchu_li · September 11, 2023, 9:37am

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…

luchu_li · September 11, 2023, 9:43am

the ordinary question:

Mark_L_Stone · September 11, 2023, 11:03am

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.

luchu_li · September 12, 2023, 1:54am

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 )

Mark_L_Stone · September 12, 2023, 1:58am

Did you read the link?

luchu_li · September 12, 2023, 2:21am

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

Mark_L_Stone · September 12, 2023, 2:32am

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

luchu_li · September 12, 2023, 2:44am

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

Mark_L_Stone · September 12, 2023, 12:05pm

That is not code which generated an error message/.

luchu_li · September 13, 2023, 2:00am

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

Mark_L_Stone · September 13, 2023, 2:05am

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.

luchu_li · September 13, 2023, 2:34am