[This linear matrix inequality appears to be unsymmetric. This is

very likely an error that will produce unexpected results. Please check

the LMI; and, if necessary, re-enter the model.]

I would really appreciate if anybody has a solution for this.

The Matalb code is:

```
function [Sum_Rate, P, Theta ] = Opt_func( M,N,K,P_max,sigma,G,f,Arr_BS,Arr_LIS,noise1,IRS,AP,rf,rg,rh)
%% init
a = sqrt(10/(1+10));
b = sqrt(1/(1+10));
Phi = diag(exp(1j*(rand(M,1)*2*pi))); % phase shifters
Theta=diag(Phi);
[L_f,LG,Lh]=pathloss(IRS,AP,rf,rg,rh);
% ZF
H_all = G'*Phi*f;
P = H_all'*pinv(H_all*H_all');
P = P/norm(P,'fro');
SumRate = [];
iter = 0;
%% Iterations UPDATE
while(1)
iter = iter+1;
[Alpha,sum_rate]=Efun(M,K,sigma,G,f,noise1,IRS,AP,rf,rg,rh,P,Phi);
SumRate = [SumRate, sum_rate];
%%
lamuda1=zeros(K,1);
for k=1:K
lamuda1(k) = noise1^2*LG(k)*L_f;
end
%%
ypsl=zeros(K,1);
for k=1:K
ypsl(k)=lamuda1(k)*a^2*b^2*M+lamuda1(k)*b^4*M;
end
%%
deta=zeros(K,1);%得儿塔
for k=1:K
deta(k) = lamuda1(k)*a^4*abs(Theta'*diag(G(:,k))*Arr_LIS)^2+lamuda1(k)*a^2*b^2*M;
end
%% Gama
Gama = zeros(K,1);
for k=1:K
[Q1wk,Q2wk] = Qw(M,P(:,k),K,a,b,lamuda1,ypsl,Arr_BS);
S_gain3 =Theta'*diag(G(:,k))*Arr_LIS;
IN_gain3 = 0;
for j=1:K
[Q1wj,Q2wj] = Qw(M,P(:,j),K,a,b,lamuda1,ypsl,Arr_BS);
IN_temp3 = Q1wj(k)*abs(Theta'*diag(G(:,k))*Arr_LIS)^2+Q2wj(k);
IN_gain3 = IN_gain3 + IN_temp3;
end
Gama(k) = sqrt((1+Alpha(k))*Q1wk(k))*S_gain3/(IN_gain3+sigma);
end
%% Doreta
Doreta = zeros(K,1);
for k=1:K
[Q1wk,Q2wk] = Qw(M,P(:,k),K,a,b,lamuda1,ypsl,Arr_BS);
IN_gain4 = 0;
for j=1:K
[Q1wj,Q2wj] = Qw(M,P(:,j),K,a,b,lamuda1,ypsl,Arr_BS);
IN_temp1 = Q1wj(k)*abs(Theta'*diag(G(:,k))*Arr_LIS)^2+Q2wj(k);
IN_gain4 = IN_gain4 + IN_temp1;
end
Doreta(k) = sqrt((1+Alpha(k))*Q2wk(k))/(IN_gain4+sigma);
end
%%
V=zeros(M,K);
for k=1:K
V(:,k)=diag(G(:,k))*Arr_LIS;
end
A = zeros(M, M);
bb = zeros(M, 1);
cc=0;
%C = 0;
for k=1:K
[Q1wk,Q2wk] = Qw(M,P(:,k),K,a,b,lamuda1,ypsl,Arr_BS);
for j=1:K
[Q1wj,Q2wj] = Qw(M,P(:,j),K,a,b,lamuda1,ypsl,Arr_BS);
cc=cc+Q1wj(k)*V(:,k)*V(:,k)';
end
A = A + cc*(abs(Gama(k))^2+abs(Doreta(k))^2);
bb = bb + sqrt((1+Alpha(k))*Q1wk(k))*conj(Gama(k))*V(:,k);
%C = C + abs(rho(k))^2*sigma;
end
cvx_begin sdp quiet
variable f_SDP
variable zeta(M)
maximize (f_SDP-trace(diag(zeta)))
subject to
for m=1:M
zeta(m) >= 0;
end
[A + diag(zeta) bb; ...
bb' -f_SDP] >= 0;
cvx_end
Theta_opt = inv(A + diag(zeta))*bb;%（25）
Theta_opt = exp(1i*angle(Theta_opt));
%%
Theta=Theta_opt;
Phi = diag(Theta');
%% if converge
if iter >=3
diff = SumRate(iter) - SumRate(iter-1);
if abs(diff)/SumRate(iter) <= 1e-2
Sum_Rate = SumRate(end);
break;
end
end
end
end
```