CVX says my problem is solved but returns infeasible solution

I am trying to solve the following problem.

clear all
clc

cvx_solver mosek

M = 4;
K_max = 2;
K = 2;
C = 1e3;
k = 1e-27;
scale = 1e8; % scale for k
scaled_k = k * scale^3;
p_0_r = 0.002;
p_r = 0.001;
p_t = 0.002;
p_l = 0.001;
O = 0.1;
B_set = 2e5:2e5:10e5;
beta = 0.05;
k_M = 1e-30;
f_M = 5e9;
T = 1;
E_vec = 0.1;
sigma_square = 1e-12;
P_dBm = 40;
P_i_dBm = 23;
P = 10^((P_dBm - 30)/10);
P_i = 10^((P_i_dBm - 30)/10);
I_ins_dBm = 0;
I_avg_dBm = 0;
I_ins = 10^((I_ins_dBm - 30)/10);
I_avg = 10^((I_avg_dBm - 30)/10);
R_min = 0.5e5;
H = 1e-3 * …
[0.0301 + 0.1077i, -0.1332 - 0.1977i;
-0.0158 - 0.0268i, -0.0936 - 0.0374i;
0.0591 - 0.1034i, 0.0918 + 0.0425i;
0.0087 + 0.0435i, 0.1574 + 0.0118i];
H_A_to_p = 1e-7 * …
[0.0108 + 0.0000i -0.0023 + 0.0253i -0.0194 + 0.0041i -0.0003 - 0.0330i;
-0.0023 - 0.0253i 0.0594 + 0.0000i 0.0137 + 0.0445i -0.0770 + 0.0077i;
-0.0194 - 0.0041i 0.0137 - 0.0445i 0.0365 + 0.0000i -0.0120 + 0.0594i;
-0.0003 + 0.0330i -0.0770 - 0.0077i -0.0120 - 0.0594i 0.1008 + 0.0000i];
h_i_to_p_sq = 1e-8 * [0.4700; 0.3187];
G = H;
Z = (G’ * G) \ G’; % [K X M]
ch_gamma_sdma = 1./( sigma_square * sum(abs(Z).^2, 2) ); % [K X 1]
H_i = zeros(M, M, K);
for i = 1:K
H_i(:,:,i) = H(:,i) * H(:,i)';
end
B = 8e5;
zeta = k_M * C * B * (f_M^2);
tau_D_L = 5.3835e-07;
Gamma_I_L = [0; 0];
q = 0;

cvx_begin
variable tau_D nonnegative;
variable tau_U nonnegative;
variable e_U(K) nonnegative;
variable b(K) nonnegative;
variable Gamma(K) nonnegative;
variable Gamma_I(K) nonnegative;
variable S(M, M, K) complex semidefinite;
MEC_comp_time = (C * B * b)/f_M;
expression tr_S(K);
expression tr_HS(K);
expression tr_HA_S(K);
for i = 1:K
tr_S(i) = trace(S(:,:,i));
tr_HS(i) = real( trace(H_i(:,:,i) * S(:,:,i)) );
tr_HA_S(i) = real( trace(H_A_to_p * S(:,:,i)) );
end

sum_S = sum(S, 3);  
expression not_i_sum_S(M, M, K);
expression tr_intf_HS(K)
for i = 1:K
    not_i_sum_S(:,:,i) = sum_S - S(:,:,i);            
    tr_intf_HS(i) = real( trace( H_i(:,:,i)*not_i_sum_S(:,:,i) ) );
end

%% grad_Q_D
%[dQ_1/dtau_D  ...  dQ_K/dtau_D  ]
%[dQ_1/dGamma_1 ... dQ_K/dGamma_1]
%[...                   ...      ]
%[dQ_1/dGamma_K ... dQ_K/dGamma_K]
grad_Q_mat = zeros(K+1, K);   
grad_Q_mat(1, :) = log2(Gamma_I_L./tau_D_L + sigma_square) - 1/log(2) * (Gamma_I_L./tau_D_L)./(Gamma_I_L./tau_D_L + sigma_square); % grad w.r.t. tau
for i = 1:K
    grad_Q_mat(i+1, i) = 1/log(2) * 1./(Gamma_I_L(i)./tau_D_L + sigma_square); % grad w.r.t. Gamma
end

%% tilde_Q_D (hat_Q_D)

Q_l = zeros(K, 1);
for i = 1:K
    Q_l(i) = tau_D_L .* log2(Gamma_I_L(i)./tau_D_L + sigma_square);
end
d_i = [tau_D; Gamma_I];  d_i_l = [tau_D_L; Gamma_I_L];

expression hat_Q(K);    
for i = 1:K 
    hat_Q(i) = Q_l(i) + grad_Q_mat(:, i).'*(d_i - d_i_l); 
end

%% Problem
num = sum( (B/scale) * b );
den = sum( e_U + p_t * tau_U + p_r * tau_D + beta * zeta * b + beta * tr_S + beta * p_0_r * tau_U );
maximize (num - q * den);
subject to
e_U + p_t * tau_U + p_r * tau_D - E_vec <= 0 % K constraints
R_min./scale - ( (B/scale)b ) <= 0 % K constraints
O * b + rel_entr(tau_D, tau_D * sigma_square + Gamma + Gamma_I) / log(2) + hat_Q <= 0
b + rel_entr(tau_U, tau_U + ch_gamma_sdma .
e_U) / log(2) <= 0 % K constraints
sum(tr_S) - tau_D * P <= 0
e_U - tau_U * P_i <= 0
tau_U + sum(MEC_comp_time) + tau_D - T <= 0
sum(e_U .* h_i_to_p_sq) - tau_U * I_ins <= 0
sum(tr_HA_S) - tau_D * I_ins <= 0
sum(e_U .* h_i_to_p_sq) + sum(tr_HA_S) - T * I_avg <= 0
Gamma - tr_HS <= 0
tr_intf_HS - Gamma_I <= 0
cvx_end

CVX says the problem is solved, but the values of the third constraints are shown as follows:

O * b + rel_entr(tau_D, tau_D * sigma_square + Gamma + Gamma_I) / log(2) + hat_Q

ans =
0.3083
0.3041

which violate the constraints. A small value can be considered as a reasonable numerical error, but this seems to be an infeasible solution. How can I get a feasible solution?

The input data is horribly scaled. That could cause all sorts of havoc.

Thank you for your comments!
Properly scaling the problem, I think I can get a feasible solution.