SDPT3: linsysolve: solution contains NaN or inf stop: difficulty in computing corrector directions; Mosek: Status: Infeasible Optimal value (cvx_optval): -Inf

Hello, everyone! I am looking forward to your help! When I using CVX to solve my optimization problem, I’m having the following problem.

The output of SDPT3 solver is as follows:


version predcorr gam expon scale_data
HKM 1 0.000 1 0
it pstep dstep pinfeas dinfeas gap prim-obj dual-obj cputime

0|0.000|0.000|1.7e+02|2.0e+02|6.6e+07|-5.547336e+03 0.000000e+00| 0:0:00| splu 737 1
linsysolve: solution contains NaN or inf
stop: difficulty in computing corrector directions

number of iterations = 1
primal objective value = -5.54733568e+03
dual objective value = 0.00000000e+00
gap := trace(XZ) = 6.59e+07
relative gap = 1.19e+04
actual relative gap = -1.00e+00
rel. primal infeas (scaled problem) = 1.67e+02
rel. dual " " " = 1.95e+02
rel. primal infeas (unscaled problem) = 0.00e+00
rel. dual " " " = 0.00e+00
norm(X), norm(y), norm(Z) = 5.3e+04, 0.0e+00, 1.2e+03
norm(A), norm(b), norm(C) = 7.3e+01, 5.9e+02, 6.2e+00
Total CPU time (secs) = 0.97
CPU time per iteration = 0.97
termination code = 0
DIMACS: 1.1e+03 0.0e+00 4.6e+02 0.0e+00 -1.0e+00 1.2e+04


Status: Solved
Optimal value (cvx_optval): +5547.34

The output of Mosek solver is as follows:
Calling Mosek 9.1.9: 6275 variables, 938 equality constraints

MOSEK Version 9.1.9 (Build date: 2019-11-21 11:34:40)
Copyright (c) MOSEK ApS, Denmark. WWW: mosek.com
Platform: Windows/64-X86

Problem
Name :
Objective sense : min
Type : CONIC (conic optimization problem)
Constraints : 938
Cones : 174
Scalar variables : 2375
Matrix variables : 390
Integer variables : 0

Optimizer started.
Presolve started.
Eliminator - tries : 0 time : 0.00
Lin. dep. - tries : 0 time : 0.00
Lin. dep. - number : 0
Presolve terminated. Time: 0.00
Optimizer terminated. Time: 0.00

Interior-point solution summary
Problem status : PRIMAL_INFEASIBLE
Solution status : PRIMAL_INFEASIBLE_CER
Dual. obj: 1.0000000000e-05 nrm: 1e+00 Viol. con: 0e+00 var: 0e+00 barvar: 0e+00 cones: 0e+00
Optimizer summary
Optimizer - time: 0.00
Interior-point - iterations : 0 time: 0.00
Basis identification - time: 0.00
Primal - iterations : 0 time: 0.00
Dual - iterations : 0 time: 0.00
Clean primal - iterations : 0 time: 0.00
Clean dual - iterations : 0 time: 0.00
Simplex - time: 0.00
Primal simplex - iterations : 0 time: 0.00
Dual simplex - iterations : 0 time: 0.00
Mixed integer - relaxations: 0 time: 0.00


Status: Infeasible
Optimal value (cvx_optval): -Inf

My code is as follows:

clear all;

S = 2;
C = 8;
U = 30;
T = 5; 
Ntx = 2; 
Nty = 2;
Nt = Ntx * Nty;  

C_pos = [175.63, 106.87; 136.30, 43.52; 45.78, 130.01; 112.98, 160.91; 203.56, 61.04; 211.69, 102.28; 199.17, 156.30; 250.84, 82.47];

S_height = 2000;  % km

s1_pos = [106, 95]; 
offset = 2; 
S1_pos = zeros(2, T);  
S1_pos(:, 1) = s1_pos';
for t = 2:T
    S1_pos(:, t) = S1_pos(:, t-1) + offset;
end

s2_pos = [106, 95];  
S2_pos = zeros(2, T); 
S2_pos(:, 1) = s2_pos';
for t = 2:T
    S2_pos(:, t) = S2_pos(:, t-1) + offset;
end


C1_user_pos = [
    
    173.67173899833338, 108.37931961155354;
    174.57454953552102, 110.32493726720683
    ];
C2_user_pos = [
    139.40651298082275, 38.67803930763517;
   
    137.82744015490138, 49.40391946858323
    ];
C3_user_pos = [
    50.220941592900196, 127.20160520535724;
   
    43.8824784449596, 125.8989690461711
    ];
C4_user_pos = [
    113.08499873947316, 158.80371594827875;
   
    116.41936814395419, 163.36229468704227;
    116.19780560454484, 160.09214237182525
    ];
C5_user_pos = [
    205.220426703268, 65.84596511048743;
    202.8614646798177, 60.19082268959938;
   
    203.66097558905423, 60.33373793784346;
    ];
C6_user_pos = [
   
    212.47281462868258, 101.69097593870022
    ];
C7_user_pos = [
    198.13023757449002, 156.0800008435594;
   
    ];
C8_user_pos = [
   
    255.8605921815085, 81.64237655577355;
    249.118031669222, 82.20905242308542
    ];

user_positions = {
    C1_user_pos, C2_user_pos, C3_user_pos, C4_user_pos, ...
    C5_user_pos, C6_user_pos, C7_user_pos, C8_user_pos
};
all_user_positions = vertcat(user_positions{:}); 
% user_counts_per_cluster = [4, 3, 4, 4, 4, 3, 4, 4];  
user_counts_per_cluster = [2, 2, 2, 3, 3, 1, 1, 2];

cluster_user_indices = cell(1, length(user_counts_per_cluster));  


start_index = 1;


for c = 1:length(user_counts_per_cluster)

    end_index = start_index + user_counts_per_cluster(c) - 1;
   
    cluster_user_indices{c} = start_index:end_index;

    start_index = end_index + 1;
end


d = zeros(S, U, T); 
for m = 1:S
    for u = 1:size(all_user_positions, 1)  
        for t = 1:T
        
            if m == 1
                satellite_pos = S1_pos(:, t); 
            else
                satellite_pos = S2_pos(:, t);  
            end
            

            user_pos = all_user_positions(u, :);  
            
     
            distance = sqrt((satellite_pos(1) - user_pos(1))^2 + ...
                            (satellite_pos(2) - user_pos(2))^2 + ...
                            (S_height - 0)^2);  
            d(m, u, t) = distance;
        end
    end
end


c = 3e8;            
f_c = 2e9;       
g_t = 10;          
g_r = 5;           

gain = zeros(S, U, T);  
for m = 1:S
    for u = 1:size(all_user_positions, 1) 
        for t = 1:T
            delta = abs(randn); 
            dist = d(m, u, t); 
            gain(m, u, t) = delta * sqrt(g_t * g_r * (c / (4 * pi * f_c * dist * 1e3))^2);
        end
    end
end
disp("gain")
disp(gain)  % 1e-5

magnitude_gain = abs(gain);
max_magnitude = max(magnitude_gain(:));
if max_magnitude > 0
    gain_normalized = gain / max_magnitude;
else
    gain_normalized = gain; 
end


% h = zeros(S, Nt, U, T); 
h = zeros(Nt, S, U, T);

doppler_shift = 100;   
time_delay = 1e-6;    


theta_x = 0.1;         
theta_y = 0.2;        


array_response = @(N, theta) 1/sqrt(N) * exp(-1j * pi * theta * (0:N-1)).';


v_x = array_response(Ntx, theta_x);
v_y = array_response(Nty, theta_y);


v_k = kron(v_x, v_y); 

for m = 1:S
    for u = 1:size(all_user_positions, 1)
        for t = 1:T
            % fprintf('m = %d, i = %d, n = %d, t = %d\n', m, i, n, t);

            g_m_i_n_t = gain_normalized(m, u, t);
       phase_factor = exp(1j * 2 * pi * (doppler_shift * t - f_c * time_delay));

            h(:, m, u, t) = g_m_i_n_t * phase_factor * v_k * 1e7;
        end
    end
end

h1 = h(:, 1, :, :);  
h2 = h(:, 2, :, :);
H1 = zeros(Nt, Nt, T);  
H2 = zeros(Nt, Nt, T);  
for t = 1:T
    h1_t = h1(:, :, t);      
    h2_t = h2(:, :, t);        
    
    H1(:, :, t) = h1_t * h1_t';  
    H2(:, :, t) = h2_t * h2_t';  
end

disp('H1:');
disp(H1) 
disp('H2:');
disp(H2)

magnitude_H1 = abs(H1);
max_magnitude = max(magnitude_H1(:));
if max_magnitude > 0
    H1_normalized = H1 / max_magnitude;
else
    H1_normalized = H1; 
end

magnitude_H2 = abs(H2);
max_magnitude = max(magnitude_H2(:));
if max_magnitude > 0
    H2_normalized = H2 / max_magnitude;
else
    H2_normalized = H2; 
end


wide_beam_cover_matrix = [
    [1, 1, 1, 1, 0, 0, 0, 0];
    [1, 0, 0, 0, 1, 1, 1, 1]
]; 
wideBeamCover = zeros(S, C, T);

for t = 1:T
    wideBeamCover(:, :, t) = wide_beam_cover_matrix;
end

spotBeamCover = zeros(S, C, T);

for t = 1:T
    for s = 1:S
       
        available_clusters = find(wideBeamCover(s, :, t) == 1);
        
       
        if length(available_clusters) <= 3
            selected_clusters = available_clusters;
        else
            selected_clusters = randsample(available_clusters, 3);
        end
        
       
        spotBeamCover(s, selected_clusters, t) = 1;
    end
end


total_users_covered_by_spotbeam = zeros(S, T);

for s = 1:S
    for t = 1:T
 
        covered_clusters = find(spotBeamCover(s, :, t) == 1);        

        total_users_covered = sum(user_counts_per_cluster(covered_clusters));        

        total_users_covered_by_spotbeam(s, t) = total_users_covered; 
    end
end


k = 1.38e-23; 
noise_T = 290;    
B = 20;    
F = 10;     

P_noise_dbw = k * noise_T * B * F;  
P_noise = 10^(P_noise_dbw / 10);

P_beam = 90 * ones(Nt, 1);


Asc2_prev = rand(U, T);  
Bc2_prev = rand(U, T);
Cn2_prev = rand(U, T);
TAOsc_prev_external = 10 + (20 - 10) * rand(U, T);
% TAOc_prev_external = rand(U, T);
TAOc_prev = 10 + (20 - 10) * rand(U, T);
TAOn_prev = 10 + (20 - 10) * rand(U, T);


total_users_per_timeslot = zeros(1, T);

for t = 1:T
    scheduled_clusters = false(1, C);
    for s = 1:S
        scheduled_clusters = scheduled_clusters | (spotBeamCover(s, :, t) > 0);  % 或 true
    end

    unique_users = [];
  
    for c = 1:C
        if scheduled_clusters(c) 
           
            cluster_users = cluster_user_indices{c};
            
            
            unique_users = union(unique_users, cluster_users);
        end
    end

    total_users_per_timeslot(t) = length(unique_users);
end 

lambda_max_Usc = rand(S, T);      
v_max_Usc = ones(Nt, S, T) / sqrt(Nt);        

lambda_max_Uc = rand(S, C, T);   
v_max_Uc = ones(Nt, S, C, T) / sqrt(Nt);     

lambda_max_Un = rand(S, U, T);  
v_max_Un = ones(Nt, S, U, T) / sqrt(Nt);     


cvx_solver mosek
% cvx_solver_settings('verbose', 3); 
% for iter = 1:5
    % TAOsc_prev_external = max(TAOsc_prev_external, 1e-6);
    cvx_begin sdp
        variable Usc(Nt, Nt, S, T) semidefinite;  
        variable Uc(Nt, Nt, S, C, T) semidefinite;
        variable Un(Nt, Nt, S, U, T) semidefinite;
        % variable Qsc(U, T);  
        variable Rsc(T);  
        % variable Qc(U, T);  
        variable Rc(C, T); 
        variable Qn(U, T);
        variable Asc1(U, T);  
        variable Asc2(U, T);
        variable Bc1(U, T);
        variable Bc2(U, T);
        variable Cn1(U, T);
        variable Cn2(U, T);
        variable TAOsc(U, T);
        variable TAOc(U, T);
        variable TAOn(U, T);
        expression O; 
        % expression H(Nt, Nt);
        expression Qsc(U, T);
        expression Qc(U, T);
        expression P_total(Nt, S, T);
        expression A2LHS(U, T);
        expression A1RHS(U, T);
        expression B2LHS(U, T);
        expression B1RHS(U, T);
        expression C2LHS(U, T);
        expression C1RHS(U, T);
        expression TAOsc_prev(U, T);
        expression value;
        expression auxiliary_var_25;
        expression auxiliary_var_26;
        expression auxiliary_var_27;

        TAOsc_prev = TAOsc_prev_external;
        
        
        for t = 1:T  % O
            O = 0;
            for u = 1:U
                for c = 1:C
                    if ismember(u, cluster_user_indices{c})  
                        O = O + Rsc(t) / total_users_per_timeslot(t) + Rc(c, t) / user_counts_per_cluster(c) + Qn(u, t);
                    end
                end
            end
        end

        for t = 1:T  % C12
            for s = 1:S 
                P_total(:, s, t) = 0;
                P_total(:, s, t) = P_total(:, s, t) + diag(Usc(:,:,s,t));
                for c = 1:C  
                    if spotBeamCover(s, c, t) == 1
                        P_total(:, s, t) = P_total(:, s, t) + diag(Uc(:,:,s,c,t));
                        for u = 1:U
                            if ismember(u, cluster_user_indices{c})
                                P_total(:, s, t) = P_total(:, s, t) + diag(Un(:,:,s,u,t));
                            end
                        end
                    end
                end
            end
        end

        f = 0;
        for s = 1:S
            for t = 1:T
                f = f + trace(Usc(:, :, s, t)) - v_max_Usc(:, s, t)' * Usc(:, :, s, t) * v_max_Usc(:, s, t);
                for c = 1:C
                    f = f + trace(Uc(:, :, s, c, t)) - v_max_Uc(:, s, c, t)' * Uc(:, :, s, c, t) * v_max_Uc(:, s, c, t);
                end
                for u = 1:U
                    f = f + trace(Un(:, :, s, u, t)) - v_max_Un(:, s, u, t)' * Un(:, :, s, u, t) * v_max_Un(:, s, u, t);
                end
            end
        end
        
        
        maximize(O - 1e-2 * f);

        
    
        subject to

            for t = 1:T
                for s = 1:S
                    for c = 1:C
                        if spotBeamCover(s, c, t) == 1
                            user_indices = cluster_user_indices{c};
                            for u = user_indices
                                Rsc(t) <= Qsc(u, t);
                            end
                        end
                    end
                end
            end
    
            for t = 1:T
                for c = 1:C
                    user_indices = cluster_user_indices{c};
                    for u = user_indices
                        Rc(c, t) <= Qc(u, t);
                    end
                end
            end
    


            for t = 1:T
                for c = 1:C
                    user_indices = cluster_user_indices{c};
                    for u = user_indices
                        Rsc(t) >= 1e-5;
                    end
                end
            end
    


              for t = 1:T
                for c = 1:C
                    user_indices = cluster_user_indices{c};
                    for u = user_indices
                        Rc(c, t) >= 1e-5;
                    end
                end
            end
    

            
            for t = 1:T
                for c = 1:C
                    user_indices = cluster_user_indices{c};
                    for u = user_indices
                         Qn(u, t) >= 1e-5;
                    end
                end
            end
    

            for t = 1:T  % C12
                for s = 1:S
                    P_total(:, s, t) <= P_beam;
                end
            end



            for t = 1:T  % C1
                for u = 1:U
                    Asc1(u, t) - Asc2(u, t) >= Qsc(u, t); 
                end
            end

        
            
            for t = 1:T  % C2
                for u = 1:U
                    Bc1(u, t) - Bc2(u, t) >= Qc(u, t); 
                end
            end

          

            for t = 1:T  % C3
                for u = 1:U
                    Cn1(u, t) - Cn2(u, t) >= Qn(u, t); 
                end
            end



            for t = 1:T  % C13左侧
                for s = 1:S
                    if s == 1
                        H = H1_normalized(:, :, t);
                    else
                        H = H2_normalized(:, :, t);
                    end
   
                    for c = 1:C
                        if spotBeamCover(s, c, t) == 1
                            A2LHS(u, t) = A2LHS(u, t) + trace(H * Uc(:, :, s, c, t));
                            for u = 1:U
                                if ismember(u, cluster_user_indices{c})
                                    A2LHS(u, t) = A2LHS(u, t) + trace(H * Un(:, :, s, u, t));
                                end
                            end
                        end
                    end
                    A2LHS(u, t) = A2LHS(u, t) + P_noise;
                    % ------
                end
            end

            for t = 1:T  % C16右侧
                for s = 1:S
                    if s == 1
                        H = H1_normalized(:, :, t);
                    else
                        H = H2_normalized(:, :, t);
                    end

                    A1RHS(u, t) = A2LHS(u, t) + trace(H * Usc(:, :, s, t));
                end
            end



            for t = 1:T  % C13
                for s = 1:S
                    for c = 1:C
                        if spotBeamCover(s, c, t) == 1
                            for u = 1:U
                                if ismember(u,cluster_user_indices{c})
                                    value = exp(Asc2_prev(u, t));
                                    % Asc2_prev(u, t) + rel_entr(1, value) <= 0;
                                    A2LHS(u, t) <= value * (Asc2(u, t) - Asc2_prev(u, t) + 1);
                                end
                            end
                        end
                    end
                end
            end


            for t = 1:T  % C14左侧
                for s = 1:S
                    scheduled_clusters = [];
                    if s == 1
                        H = H1_normalized(:, :, t);
                    else
                        H = H2_normalized(:, :, t);
                    end

                    for c = 1:C
                        if spotBeamCover(s, c, t) == 1
                            scheduled_clusters = [scheduled_clusters, c];
                        end
                    end

                    for c = 1:C
                        if ismember(c, scheduled_clusters)
                            for j = 1:C
                                if ismember(j, scheduled_clusters) && j ~= c
                                    B2LHS(u, t) = B2LHS(u, t) + trace(H * Uc(:, :, s, j, t));
                                end
                            end

                            for u = 1:U
                                if ismember(u, cluster_user_indices{c})
                                    B2LHS(u, t) = B2LHS(u, t) + trace(H * Un(:, :, s, u, t));
                                end
                            end
                        end
                    end
                    B2LHS(u, t) = B2LHS(u, t) + P_noise;
                end
            end

            for t = 1:T  % C18右侧 --
                for s = 1:S
                    scheduled_clusters = [];
                    if s == 1
                        H = H1_normalized(:, :, t);
                    else
                        H = H2_normalized(:, :, t);
                    end
                    % ---外层求和内
                    for c = 1:C
                        if spotBeamCover(s, c, t) == 1
                            scheduled_clusters = [scheduled_clusters, c];
                            user_indices = cluster_user_indices{c};  % 属于簇 c 的用户索引列表
                            for u = user_indices
                                B1RHS(u, t) = B1RHS(u, t) + trace(H * Uc(:, :, s, c, t));
                            end
                        end
                    end
                    for c = 1:C
                        if ismember(c, scheduled_clusters)
                            for j = 1:C
                                if ismember(j, scheduled_clusters) && j ~= c
                                    B1RHS(u, t) = B1RHS(u, t) + trace(H * Uc(:, :, s, j, t));
                                end
                            end

                            for u = 1:U
                                if ismember(u, cluster_user_indices{c})
                                    B1RHS(u, t) = B1RHS(u, t) + trace(H * Un(:, :, s, u, t));
                                end
                            end
                        end
                    end
                    B1RHS(u, t) = B1RHS(u, t) + P_noise;
                    % ------
                end
            end

            for t = 1:T  % C14
                for s = 1:S
                    for c = 1:C
                        if spotBeamCover(s, c, t) == 1
                            for u = 1:U
                                if ismember(u,cluster_user_indices{c})
                                    value = exp(Bc2_prev(u, t));
                                    % Bc2_prev(u, t) + rel_entr(1, value) <= 0;
                                    B2LHS(u, t) <= value * (Bc2(u, t) - Bc2_prev(u, t) + 1);
                                end
                            end
                        end
                    end
                end
            end


            for t = 1:T  % C15左侧
                for s = 1:S
                    if s == 1
                        H = H1_normalized(:, :, t);
                    else
                        H = H2_normalized(:, :, t);
                    end

                    for u = 1:U
                        for c = 1:C
                            if ismember(u, cluster_user_indices{c}) 
                                cluster_users = cluster_user_indices{c};
                                other_users = cluster_users(cluster_users ~= u);  % 排除用户 u 本身
                                for v = other_users
                                    C2LHS(u, t) = C2LHS(u, t) + trace(H * Un(:, :, s, v, t));
                                end
                                break;  % 用户只属于其中一个簇
                            end
                        end
                    end
                end
            end

            for t = 1:T  % C20右侧
                for s = 1:S
                    if s == 1
                        H = H1_normalized(:, :, t);
                    else
                        H = H2_normalized(:, :, t);
                    end

                    for c = 1:C
                        if spotBeamCover(s, c, t) == 1
                            for u = 1:U
                                if ismember(u, cluster_user_indices{c})
                                    C1RHS(u, t) = C1RHS(u, t) + trace(H * Un(:, :, s, u, t));
                                end
                            end
                        end
                    end
                    C1RHS(u, t) = C1RHS(u, t) + P_noise;
                end
            end

            for t = 1:T  % C15
                for s = 1:S
                    for c = 1:C
                        if spotBeamCover(s, c, t) == 1
                            for u = 1:U
                                if ismember(u,cluster_user_indices{c})
                                    value = exp(Cn2_prev(u, t));
                                    % Cn2_prev(u, t) + rel_entr(1, value) <= 0;
                                    C2LHS(u, t) <= value * (Cn2(u, t) - Cn2_prev(u, t) + 1);
                                end
                            end
                        end
                    end
                end
            end

            for t = 1:T  % C16
                for s = 1:S
                    for c = 1:C
                        if spotBeamCover(s, c, t) == 1
                            for u = 1:U
                                if ismember(u,cluster_user_indices{c})
                                    TAOsc(u, t) <= A1RHS(u, t);
                                end
                            end
                        end
                    end
                end
            end

            for t = 1:T  % C18
                for s = 1:S
                    for c = 1:C
                        if spotBeamCover(s, c, t) == 1
                            for u = 1:U
                                if ismember(u,cluster_user_indices{c})
                                    TAOc(u, t) <= B1RHS(u, t);
                                end
                            end
                        end
                    end
                end
            end

            for t = 1:T  % C20
                for s = 1:S
                    for c = 1:C
                        if spotBeamCover(s, c, t) == 1
                            for u = 1:U
                                if ismember(u,cluster_user_indices{c})
                                    TAOn(u, t) <= C1RHS(u, t);
                                end
                            end
                        end
                    end
                end
            end

            for t = 1:T  % C25
                for s = 1:S
                    for c = 1:C
                        if spotBeamCover(s, c, t) == 1
                            for u = 1:U
                                if ismember(u,cluster_user_indices{c})
                                    auxiliary_var_25 = log(TAOsc_prev(u, t)) + 1;
                                    norm([TAOsc(u, t) + Asc1(u, t) - auxiliary_var_25, 2 * sqrt(TAOsc_prev(u, t))], 2) ...
                                        <= TAOsc(u, t) - Asc1(u, t) + auxiliary_var_25;  % log(TAOsc_prev(u, t))
                                end
                            end
                        end
                    end
                end
            end

            for t = 1:T  % C26
                for s = 1:S
                    for c = 1:C
                        if spotBeamCover(s, c, t) == 1
                            for u = 1:U
                                if ismember(u,cluster_user_indices{c})
                                    auxiliary_var_26 = log(TAOc_prev(u, t)) + 1;
                                    norm([TAOc(u, t) + Bc1(u, t) - auxiliary_var_26, 2 * sqrt(TAOc_prev(u, t))], 2) ...
                                        <= TAOc(u, t) - Bc1(u, t) + auxiliary_var_26;
                                end
                            end
                        end
                    end
                end
            end

            for t = 1:T  % C27
                for s = 1:S
                    for c = 1:C
                        if spotBeamCover(s, c, t) == 1
                            for u = 1:U
                                if ismember(u,cluster_user_indices{c})
                                    auxiliary_var_27 = log(TAOn_prev(u, t)) + 1;
                                    norm([TAOn(u, t) + Cn1(u, t) - auxiliary_var_27, 2 * sqrt(TAOn_prev(u, t))], 2) ...
                                        <= TAOn(u, t) - Cn1(u, t) + auxiliary_var_27;
                                        % <= TAOn(u, t) - Cn1(u, t) - rel_entr(1, TAOn_prev(u, t)) + 1;
                                end
                            end
                        end
                    end
                end
            end

    cvx_end

    % Asc2_prev = Asc2 / 1e9;
    % Bc2_prev = Bc2 / 1e9;
    % Cn2_prev = Cn2 / 1e13; 

    % 计算 Asc2 的数量级并归一化
    if max(abs(Asc2(:))) > 0
        scale_Asc2 = 10^floor(log10(max(abs(Asc2(:)))));
        Asc2_prev = Asc2 / scale_Asc2;
    else
        Asc2_prev = Asc2;
    end
    
    % 计算 Bc2 的数量级并归一化
    if max(abs(Bc2(:))) > 0
        scale_Bc2 = 10^floor(log10(max(abs(Bc2(:)))));
        Bc2_prev = Bc2 / scale_Bc2;
    else
        Bc2_prev = Bc2;
    end
    
    % 计算 Cn2 的数量级并归一化
    if max(abs(Cn2(:))) > 0
        scale_Cn2 = 10^floor(log10(max(abs(Cn2(:)))));
        Cn2_prev = Cn2 / scale_Cn2;
    else
        Cn2_prev = Cn2;
    end

    TAOsc_prev_external = abs(TAOsc);
    % TAOc_prev_external = abs(TAOc);
    TAOc_prev = abs(TAOc);
    TAOn_prev = abs(TAOn);

    for s = 1:S
        for t = 1:T
            % 计算 Usc 的数量级(以 10 为底的对数)
            scale_Usc = floor(log10(max(abs(Usc(:,:,s,t)), [], 'all')));

            % 如果数量级不为零,进行数量级移除
            if scale_Usc ~= -Inf
                Usc(:,:,s,t) = Usc(:,:,s,t) / 10^scale_Usc;
            end
        end
    end
    
    % 调整 Uc
    for s = 1:S
        for c = 1:C
            for t = 1:T
                % 计算 Uc 的数量级(以 10 为底的对数)
                scale_Uc = floor(log10(max(abs(Uc(:,:,s,c,t)), [], 'all')));

                % 如果数量级不为零,进行数量级移除
                if scale_Uc ~= -Inf
                    Uc(:,:,s,c,t) = Uc(:,:,s,c,t) / 10^scale_Uc;
                end
            end
        end
    end
    
    % 调整 Un
    for s = 1:S
        for u = 1:U
            for t = 1:T
                % 计算 Un 的数量级(以 10 为底的对数)
                scale_Un = floor(log10(max(abs(Un(:,:,s,u,t)), [], 'all')));

                % 如果数量级不为零,进行数量级移除
                if scale_Un ~= -Inf
                    Un(:,:,s,u,t) = Un(:,:,s,u,t) / 10^scale_Un;
                end
            end
        end
    end

    for s = 1:S
        for t = 1:T
            % Usc 的处理
            P = Usc(:, :, s, t);
            % disp("P")
            % disp(P)
            [V, D] = eig(double(P));
            [lambda_max, idx] = max(diag(D));
            v_max = V(:, idx);
            v_max = v_max / norm(v_max); % 归一化
            lambda_max_Usc(s, t) = lambda_max;
            v_max_Usc(:, s, t) = v_max;
        end
    end

    for s = 1:S
        for c = 1:C
            for t = 1:T
                P = Uc(:, :, s, c, t);
                [V, D] = eig(double(P));
                [lambda_max, idx] = max(diag(D));
                v_max = V(:, idx);
                v_max = v_max / norm(v_max);
                lambda_max_Uc(s, c, t) = lambda_max;
                v_max_Uc(:, s, c, t) = v_max;
            end
        end
    end

    for s = 1:S
        for u = 1:U
            for t = 1:T
                P = Un(:, :, s, u, t);
                [V, D] = eig(double(P));
                [lambda_max, idx] = max(diag(D));
                v_max = V(:, idx);
                v_max = v_max / norm(v_max);
                lambda_max_Un(s, u, t) = lambda_max;
                v_max_Un(:, s, u, t) = v_max;
            end
        end
    end
% end

Does the infeasible status of Mosek solver means that there is a problem with my problem modeling or the process of converting a math problem into code? Or is there still a problem with my numeric scaling? What can I do?

Thanks for all suggestions!

Did you follow the advice in the link I previously provided Debugging infeasible models - YALMIP .?

Your model is complicated. You presumably understand what your model is doing, whereas forum readers don’t. So you need to figure out how to simplify things to figure out what’s going on and diagnose the infeasibility reported by Mosek.