Cvx problem optimization

Hi all. I have a problem with CVX. How to scale the parameter values of the optimization problem specially d_RE و d_SR و d_SD و d_RD و d_ER , Ne , Nd , sap in MATLAB?
When my program is run, it sometimes returns NaN for one of the values.
How can I solve this problem?
clc
clear all

M = 4 ;
K = 1:M;
%n =rand(M)
%n = rand(M,M)+1i*randM,M)

%n = randi(5,M,M)
%c1 = diag(n)
%theta_r = diag(c1)
theta_t=exp(1j.*rand(M,M).*2.*pi);
theta_r=exp(1j.*rand(M,M).*2.*pi);
%c2 = diag(n)
%theta_t = diag(c2)

NR = 4;
NT = 4;
%w_t = rand(NT,1);
%p_E = 10;
%Q_E = 10.^(p_E./10);
Pe_vdb = 0:2:4
;%dB
Pe_v = 10.^(Pe_vdb./10);%watt
Ne = 10.^(-80./10);
Nd = 10.^(-80./10);
Ps = 10.^(20./10);
sap = 10.^(-60./10);
%w_r = rand(NR,1);
itration = 20;

d_RE = sqrt(409); %distance of R & E (m)
d_SR = sqrt(59); %distance of S & R (m)
d_SD = 10; %distance of S & D (m)
d_RD = sqrt(59); %distance of R & D (m)
d_ER = sqrt(409); %distance of E & R (m)
Ex = zeros(1,length(Pe_vdb));
Ex1 = zeros(1,length(Pe_vdb));
Ex2 = zeros(1,length(Pe_vdb));
Ex3 = zeros(1,length(Pe_vdb));
t = 3;
for ii = 1:length(Pe_vdb)

ii
Q_E = Pe_v(ii);
c = 0;    
j  = 0; 

c1 = 0;
c2 = 0;
c3 = 0;

for i = 1:itration

h_RD=(1/sqrt(2))*(randn(1,M)+sqrt(-1)*randn(1,M));

h_SR=(1/sqrt(2))*(randn(M,1)+sqrt(-1)*randn(M,1));

h_SD=(1/sqrt(2))*(randn(1,1)+sqrt(-1)*randn(1,1));

H_ER=(sqrt(1)/sqrt(2))*(randn(M,NT)+sqrt(-1)*randn(M,NT));

H_RE=(sqrt(1)/sqrt(2))*(randn(NR,M)+sqrt(-1)*randn(NR,M));

H_EE=(sqrt(sap)/sqrt(2))*(randn(NR,NT)+sqrt(-1)*randn(NR,NT));

u = H_EE + H_REtheta_rH_ER;

%SINR_d = Ps*(abs(h_SD + h_RDtheta_rh_SR)^2)/(Q_E*(abs( h_RDtheta_tH_ER*w_t)^2) + Nd)

%SINR_e = Ps*(abs(w_r’H_REtheta_th_SR)^2)/(Q_E(abs(w_r’H_EEw_t + w_r’H_REtheta_rH_ERw_t )^2) + Ne)
%---------------------------- BF Design 2_TZF ---------------------------------------------------------------------------------------------------------------------------

    Wr_TZF = (H_RE*theta_t*h_SR)/norm((H_RE*theta_t*h_SR));
   
    H_EE_hat = (H_EE + H_RE*theta_r*H_ER);
    
    delta_tz = eye(NT)-(  ( (H_EE_hat')*Wr_TZF*(Wr_TZF')*(H_EE_hat) )/ ( ( norm (H_EE_hat*Wr_TZF ) ).^2 ) );
    
    Wt_TZF = ( delta_tz*( (h_RD*theta_t*H_ER)' ) )/norm(( delta_tz*( (h_RD*theta_t*H_ER)' ) ));
    
   SINR_d_TZF = ( ( Ps/(d_SD.^t) )*( ( abs( h_SD + h_RD*theta_r*h_SR ) ).^2 ) )/ ( ( Q_E/( (d_RD.^t)*(d_ER.^t) ) )*( ( abs( h_RD*theta_t*H_ER*Wt_TZF) ).^2 ) + Nd  );

   
   %SINR_e_TZF = ( ( Ps/( (d_SR.^t)*(d_RE.^t) ) )*( ( norm( H_RE*theta_t*h_SR ) ).^2 ) )/ ( Ne );
   SINR_e_TZF = ( ( Ps/( (d_SR.^t)*(d_RE.^t) ) )*( ( abs( Wr_TZF'*H_RE*theta_t*h_SR ) ).^2 ) )/ ( Q_E*( ( abs( Wr_TZF'*H_EE*Wt_TZF + Wr_TZF'*H_RE*theta_r*H_ER*Wt_TZF ) ).^2 )+ Ne );
   
   
  %---------------------------------BF Design 3_MRC_MRT---------------------------------------------------------------------------------------------------------------------
  
  
  %% ................... BF Design 1_RZF .........................
   
      Wt_RZF = (h_RD*theta_t*H_ER)'/norm(( h_RD*theta_t*H_ER)');
     
      H_EE_hat = (H_EE + H_RE*theta_r*H_ER);
     
      delta_rz = eye(NR)-( (H_EE_hat*Wt_RZF*(Wt_RZF')*(H_EE_hat'))/( (  norm ( H_EE_hat*Wt_RZF ) ).^2 ) ) ;
      
      %Wr_RZF = (H_RE*theta_t*h_SR)/norm(( H_RE*theta_t*h_SR ));
 
       Wr_RZF = ( delta_rz*( (H_RE*theta_t*h_SR) ) )/norm(( delta_rz*( (H_RE*theta_t*h_SR) ) ));
     
      %Wt_ZF = ( delta2.*(h_RD*theta_t*H_ER) )./( norm(( delta2.*(h_RD*theta_t*H_ER) ))) ; 

% SINR_d_RZF = ( ( Ps/(d_SD.^t) )( ( abs( h_SD + h_RDtheta_rh_SR ) ).^2 ) )/ ( ( Q_E/( (d_RD.^t)(d_ER.^t) ) )( ( norm( h_RDtheta_t*H_ER) ).^2 ) + Nd );

SINR_d_RZF = ( ( Ps/(d_SD.^t) )( ( abs( h_SD + h_RDtheta_rh_SR ) ).^2 ) )/ ( ( Q_E/( (d_RD.^t)(d_ER.^t) ) )( ( abs( h_RDtheta_tH_ER Wt_RZF) ).^2 ) + Nd );

      SINR_e_RZF = ( ( Ps/( (d_SR.^t)*(d_RE.^t) ) )*( ( abs( Wr_RZF'*H_RE*theta_t*h_SR ) ).^2 ) )/ ( Q_E*( ( abs( Wr_RZF'*H_EE*Wt_RZF + Wr_RZF'*H_RE*theta_r*H_ER*Wt_RZF) ).^2 )+ Ne );
  
     
    
     %Wr_ZF = (delta*((h_RD*theta_t*H_ER)'))/(norm( delta*((h_RD*theta_t*H_ER)') ));
     %hse_hat1 = abs(Wr_ZF'*H_RE*theta_t*h_SR)^2;
     %SINRd_MRT = (Ps .*(abs(h_SD + h_RD*theta_r*h_SR)^2))./( (Q_E*(norm( h_RD*).^2) + Nd) );
     %SINRd_MRT =(Ps*(abs(h_SD + h_RD*theta_r*h_SR).^2).*(d4.^t).*(d5.^t) )/( ( (d3.^t).*Q_E*(norm( h_RD*theta_t*H_ER).^2) ) + Nd )
     %SINRe_MRT = (Ps .* hse_hat1) ./ ( Ne);
     %SINRd_MRT = (Ps*(abs(h_SD + h_RD*theta_r*h_SR).^2).*(d4.^t).*(d5.^t) )/( (d3.^t).*(  (Q_E)*( norm( h_RD*theta_t*H_ER).^2 ) + Nd ) );
     %SINRe_MRT = ( Ps*( abs(Wr_ZF'*H_RE*theta_t*h_SR).^2) ) /(  (d1.^t).*(d2.^t).*( Ne ) );
     %SINR_d_tzf = (Ps*(abs(h_SD + h_RD*theta_r*h_SR).^2).*(d4.^t).*(d5.^t) )/( (d3.^t).*(  Q_E*( abs( h_RD*theta_t*H_ER*Wt_MRT).^2 ) + Nd ) );
     %SINR_e_tzf = (Ps*(norm(H_RE*theta_t*h_SR).^2)) /(  (d1.^t).*(d2.^t).*( Ne ) );

% Wt_MRT = rand(NT,1)
%
% delta2 = eye( NR )-( (H_EEWt_MRTWt_MRT’H_EE’) /(norm(H_EEWt_MRT)^2 ) );
%
% %Wr_ZF = (delta2*(h_RDtheta_tH_ER)‘)/norm(delta2*(h_RDtheta_tH_ER)’);
%
% Wr_ZF = rand(NR,1)
% %hse_hat1 = abs(Wr_ZF’H_REtheta_th_SR)^2;
%
% H_RE_hat = (H_RE’)diag(NR)H_RE
%
% hse_hat1 = norm(H_RE_hat
theta_t
h_SR)^2;
%
% %( Q_E .
abs(h_SD+h_RDtheta_rh_SR).^2 ) ./ ( ( (Q_E) .* norm(h_RDtheta_tH_ER).^2 + (Nd) ));
% SINRd_MRT = Ps*(abs(h_SD + h_RDtheta_rh_SR)^2)/(Q_E*(abs( h_RDtheta_tH_ER* Wt_MRT)^2) + Nd)
% SINRe_MRT = Ps*(abs( Wr_ZF’H_REtheta_th_SR)^2)/(Q_E(abs( Wr_ZF’H_EEWt_MRT + Wr_ZF’H_REtheta_rH_ERWt_MRT )^2) + Ne)
% %%

y_max = 1 + (Q_E/Nd)( ( norm( h_RDtheta_t*H_ER) ).^2 );

for z = 1:20;
y = 1 + (y_max./20).*(z-1);

   %............. Optimization Problem ..............
      
      cvx_begin  sdp quiet
      variable Z(NT,NT) hermitian;
      variable s;
      variable f_y;
         
      dual variables p1 p2 p3 p4 p5 

      minimize f_y; 
      
      subject to
      
      p1: s == hermitian_semidefinite(1);
      p2: trace(Z) == s;
      p3: s + (Q_E/Ne)*trace(Z*( u'  )*(  u ) ) == 1
      p4: (Q_E/Ne)*trace(Z*H_ER'*theta_t'*h_RD'*h_RD*theta_t*H_ER) == s .* (y-1)
      p5: Z == hermitian_semidefinite(NT)
          f_y -  (Q_E/Ne)*trace(Z*(u'*H_RE*theta_t*h_SR*h_SR'*theta_t'*H_RE'*u))== hermitian_semidefinite(1)
      cvx_end

 if( strcmp( cvx_status, 'Infeasible' ) == 1 );
         
         Exy(z) = 0;
     else
         Exy(z) = f_y + ( ((Ne/Nd).*(abs(h_SD + h_RD*theta_r*h_SR)^2))./y) ; 
     end       

end

E_xy = nonzeros(Exy);
[l,e] = min(E_xy) ;

       if mean(Exy) ~= 0
           
           j = find(Exy==l);
       end
       
      y = 1 + (y_max-1)/20*(j-1);

cvx_begin sdp quiet
variable Z(NT,NT) hermitian;
variable s;
variable v;

      dual variables p1 p2 p3 p4 p5 

      minimize v ; 
      
      subject to
      
      p1: s == hermitian_semidefinite(1);
      p2: trace(Z) == s;
      p3: s + (Q_E/Ne)*trace(Z*( u'  )*(  u ) ) == 1
      p4: (Q_E/Ne)*trace(Z*H_ER'*theta_t'*h_RD'*h_RD*theta_t*H_ER) == s .* (y-1)
      p5: Z == hermitian_semidefinite(NT)
          v -  (Q_E/Ne)*trace(Z*(u'*H_RE*theta_t*h_SR*h_SR'*theta_t'*H_RE'*u))== hermitian_semidefinite(1)
      cvx_end
      %% .......................... W and w_t .............................
     if length(cvx_status)==6 && sum(cvx_status == 'Solved')==6;
         W     = Z/s ;
         h_t   = real(cvx_optval);
         [V D] = eigs(W);
         d     = real(diag(D));
         [d I] = sort(d,'descend');
         w_opt = V(:,I(1));
         W_opt = W;
     else
         w_opt = zeros(NT,1);
         W_opt = zeros(NT);
         h_t   = 0;
     end
       
      %B = (H_EE + H_RE*theta_r*H_ER)*w_opt*(w_opt')*((H_EE + H_RE*theta_r*H_ER)')
     
       %Wr = ((inv( (Q_E/Nd).*(B) + eye(NR) )) *H_RE*theta_t*h_SR )/norm( ((inv( (Q_E/Nd).*(B) + eye(NR) )) *H_RE*theta_t*h_SR ))
      
     A = Q_E*( H_EE*w_opt + H_RE*theta_r*H_ER*w_opt )*( (H_EE*w_opt + H_RE*theta_r*H_ER*w_opt)')+Ne;
       
     Wr = inv(A)*H_RE*theta_t*h_SR/norm(inv(A)*H_RE*theta_t*h_SR);

% SINR_d = Ps*(abs(h_SD + h_RDtheta_rh_SR)^2)/(Q_E*(abs( h_RDtheta_tH_ERw_opt)^2) + Nd);
% SINR_e = Ps
(abs(Wr’H_REtheta_th_SR)^2)/(Q_E(abs(Wr’H_EEw_opt + Wr’H_REtheta_rH_ERw_opt )^2) + Ne);
%
SINR_d = ( ( Ps/(d_SD.^t) )( ( abs( h_SD + h_RDtheta_rh_SR ) ).^2 ) )/ ( ( Q_E/( (d_RD.^t)(d_ER.^t) ) )( ( abs( h_RDtheta_tH_ERw_opt) ).^2 ) + Nd )

      SINR_e = ( ( Ps/( (d_SR.^t)*(d_RE.^t) ) )*( ( abs( Wr'*H_RE*theta_t*h_SR ) ).^2 ) )/ ( Q_E*( ( abs( Wr'*H_EE*w_opt + Wr'*H_RE*theta_r*H_ER*w_opt) ).^2 )+ Ne )
  
       
       
         if SINR_e >= SINR_d;
        
        c = c + 1;
    end
       
      if SINR_e_TZF >= SINR_d_TZF;
        
        c1 = c1 + 1;
       end 
       
     if SINR_e_RZF >= SINR_d_RZF;
        
        c2 = c2 + 1;
     end 

% if SINRe_MRT >= SINRd_MRT;
%
% c3 = c3 + 1;
% end

end

Ex(ii) = c./itration;
 Ex1(ii) = c1./itration;
 Ex2(ii) = c2./itration;

% Ex3(ii) = c3./itration;

end

xq2 = 0:0.01:4;
p1 = pchip(Pe_vdb,Ex,xq2);
plot(xq2,p1,‘b’)
hold on
% plot(Pe_vdb,Ex)
hold on
xq3 = 0:0.01:4;
p2 = pchip(Pe_vdb,Ex1,xq3);
plot(xq3,p2,‘r’)
hold on
% plot(Pe_vdb,Ex1,‘r’)
% hold on

xq4 = 0:0.01:4;
p3 = pchip(Pe_vdb,Ex2,xq4);
plot(xq4,p3,‘g’)
hold on
% plot(Pe_vdb,Ex2,‘g’)
% hold on

% xq5 = 0:0.01:4;
% p4 = pchip(Pe_vdb,Ex3,xq5);
% plot(xq5,p4,‘k’)
% hold on
% plot(Pe_vdb,Ex3,‘k’)
% hold on
grid on
box on


command window
SINR_d =

1.0913e+05

SINR_e =

1.9248e+04

SINR_d =

2.2522e+08

SINR_e =

NaN

If possible, choose units so that all the non-zero input data in the optimization problem is within a few orders of magnitude of 1.