The CVX optimization rate is negative

1

Hello, everyone! I am looking forward to your help! When I using CVX to solve my optimization problem, the optimal rate value of its output is negative. And I don’t know why this is happening.

I also have another question. The order of magnitude between the channel and the solved precoded matrix varies greatly. And I can’t handle them by scaling. This problem causes me to get a value of 0 when I multiply the channel by the precoding matrix and then find the trace, which I suspect may be the reason for the negative rate of solving.

The values solved are as follows:
Uc:
val(:,:,1,1) =
1.0e-16 *
0.2628 -0.2126
-0.2126 0.2628
val(:,:,1,2) =
0 0
0 0
val(:,:,1,3) =
1.0e-16 *
0.1718 -0.1404
-0.1404 0.1718
Un:
val(:,:,1,1,1) =
1.0e-16 *
0.8039 -0.6748
-0.6748 0.8039
val(:,:,1,2,1) =
0 0
0 0
val(:,:,1,3,1) =
0 0
0 0
val(:,:,1,1,2) =
1.0e-16 *
0.5550 -0.4624
-0.4624 0.5550
val(:,:,1,2,2) =
0 0
0 0
val(:,:,1,3,2) =
0 0
0 0
val(:,:,1,1,3) =
0 0
0 0
val(:,:,1,2,3) =
0 0
0 0
val(:,:,1,3,3) =
0 0
0 0
val(:,:,1,1,4) =
0 0
0 0
val(:,:,1,2,4) =
0 0
0 0
val(:,:,1,3,4) =
0 0
0 0
val(:,:,1,1,5) =
0 0
0 0
val(:,:,1,2,5) =
0 0
0 0
val(:,:,1,3,5) =
1.0e-16 *
0.2755 -0.2288
-0.2288 0.2755
val(:,:,1,1,6) =
0 0
0 0
val(:,:,1,2,6) =
0 0
0 0
val(:,:,1,3,6) =
1.0e-16 *
0.6211 -0.5183
-0.5183 0.6211
Rc:
-2.34635110647975 0 -0.694971923611425
Rn:
val(:,:,1) =
-0.6010 0 0
val(:,:,2) =
-0.6851 0 0
val(:,:,3) =
0 0 0
val(:,:,4) =
0 0 0
val(:,:,5) =
0 0 -0.3022
val(:,:,6) =
0 0 -1.9534

The output of MOSEK solver is as follows:
Calling Mosek 9.1.9: 248 variables, 34 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 : 34
Cones : 29
Scalar variables : 248
Matrix variables : 0
Integer variables : 0

Optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator - tries : 2 time : 0.00
Lin. dep. - tries : 1 time : 0.00
Lin. dep. - number : 0
Presolve terminated. Time: 0.00
Problem
Name :
Objective sense : min
Type : CONIC (conic optimization problem)
Constraints : 34
Cones : 29
Scalar variables : 248
Matrix variables : 0
Integer variables : 0

Optimizer - threads : 14
Optimizer - solved problem : the primal
Optimizer - Constraints : 20
Optimizer - Cones : 15
Optimizer - Scalar variables : 67 conic : 45
Optimizer - Semi-definite variables: 0 scalarized : 0
Factor - setup time : 0.00 dense det. time : 0.00
Factor - ML order time : 0.00 GP order time : 0.00
Factor - nonzeros before factor : 210 after factor : 210
Factor - dense dim. : 0 flops : 1.18e+04
ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME
0 1.7e+01 1.5e+08 7.0e+08 0.00e+00 7.037003310e+08 4.953821135e-01 1.0e+00 0.00
1 3.2e+00 2.8e+07 3.0e+08 -1.00e+00 7.037001042e+08 -1.852622882e+02 1.9e-01 0.00
2 6.0e-02 5.3e+05 4.2e+07 -1.00e+00 7.036849851e+08 -1.211829212e+04 3.5e-03 0.00
3 6.7e-05 5.9e+02 1.4e+06 -1.00e+00 6.899292444e+08 -1.084602622e+07 3.9e-06 0.00
4 6.7e-06 5.9e+01 4.3e+05 -9.92e-01 5.760392211e+08 -9.738836061e+07 4.0e-07 0.00
5 1.1e-06 9.6e+00 1.4e+05 -8.87e-01 3.391866105e+08 -1.465338130e+08 6.5e-08 0.00
6 2.7e-07 2.4e+00 3.7e+04 -3.16e-01 2.119657195e+08 -2.526949919e+07 1.6e-08 0.00
7 4.1e-08 3.6e-01 2.7e+03 3.96e-01 3.845493026e+07 -7.154369251e+06 2.4e-09 0.00
8 1.1e-09 1.0e-02 9.5e+00 9.88e-01 1.132051551e+06 -1.532086876e+05 6.7e-11 0.00
9 2.3e-12 2.0e-05 8.6e-04 1.00e+00 2.296979606e+03 -3.109707989e+02 1.4e-13 0.00
10 6.7e-14 5.9e-07 4.3e-06 1.00e+00 6.723862119e+01 -8.837303163e+00 4.0e-15 0.00
11 1.1e-14 9.8e-08 3.2e-07 9.71e-01 1.419426904e+01 1.301369885e+00 6.5e-16 0.00
12 1.8e-15 1.6e-08 2.2e-08 9.15e-01 6.834785659e+00 4.654774796e+00 1.1e-16 0.00
13 5.9e-16 8.4e-10 3.7e-10 8.54e-01 6.144493468e+00 6.018156975e+00 9.7e-18 0.00
14 4.2e-16 3.1e-11 3.2e-12 1.05e+00 6.110433088e+00 6.105980852e+00 1.3e-18 0.00
15 4.2e-16 2.6e-11 2.5e-12 1.00e+00 6.120404420e+00 6.116803851e+00 1.1e-18 0.00
16 4.2e-16 2.6e-11 2.5e-12 3.22e-02 6.120404420e+00 6.116803851e+00 1.1e-18 0.00
17 5.0e-16 9.0e-12 6.5e-13 1.00e+00 6.119778342e+00 6.118597148e+00 3.7e-19 0.00
18 5.4e-16 8.8e-13 2.0e-14 1.00e+00 6.120836490e+00 6.120721363e+00 3.6e-20 0.00
19 6.7e-16 1.1e-13 8.8e-16 1.00e+00 6.121491716e+00 6.121477561e+00 4.5e-21 0.00
20 7.3e-16 6.7e-15 6.0e-18 1.00e+00 6.122087580e+00 6.122087083e+00 1.6e-22 0.00
21 1.0e-15 2.7e-15 2.5e-22 1.00e+00 6.122694289e+00 6.122694288e+00 2.9e-25 0.00
Optimizer terminated. Time: 0.00

Interior-point solution summary
Problem status : PRIMAL_AND_DUAL_FEASIBLE
Solution status : OPTIMAL
Primal. obj: 6.1226942889e+00 nrm: 4e+01 Viol. con: 2e+01 var: 1e+01 cones: 0e+00
Dual. obj: 6.1226942883e+00 nrm: 4e+07 Viol. con: 0e+00 var: 1e-07 cones: 0e+00
Optimizer summary
Optimizer - time: 0.00
Interior-point - iterations : 21 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: Solved
Optimal value (cvx_optval): -6.58304

And my code is as follows:

cvx_clear;
rng(42);

S = 1;
C = 3;  % 点波束2
U = 6;
T = 1;  
% 天线数量
Ntx = 1; 
Nty = 2;
Nt = Ntx * Nty;

S_height = 2000;  % km   卫星 x y z
s1_pos = [106, 95]; 
offset = 4; 
S1_pos = zeros(2, T);  % 卫星1每个时隙的位置
S1_pos(:, 1) = s1_pos';
% for t = 2:T
%     S1_pos(:, t) = S1_pos(:, t-1) + offset;
% end

% 簇内用户位置  是否合适  -----------------------------------
C1_user_pos = [
    173.67173899833338, 108.37931961155354;
    174.57454953552102, 110.32493726720683
    ];
C2_user_pos = [
    136.32860551422957, 42.16289711432755;
    137.82744015490138, 49.40391946858323
    ];
C3_user_pos = [
    44.000564580961424, 129.4559770874314;
    43.8824784449596, 125.8989690461711
    ];

% 将所有簇的用户位置合并为一个cell array
user_positions = {
    C1_user_pos, C2_user_pos, C3_user_pos
};
all_user_positions = vertcat(user_positions{:});  % 合并所有簇内用户位置
user_counts_per_cluster = [2, 2, 2];  % 各簇内用户数量
% 用户索引   cluster_user_indices
cluster_user_indices = cell(1, length(user_counts_per_cluster));  % cluster_user_indices{c}用它取出每组用户的索引
% 起始用户索引
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:U  
        for t = 1:T
            % 获取卫星的位置，S1_pos或S2_pos，根据卫星编号
            if m == 1
                satellite_pos = S1_pos(:, t);  % 卫星1的位置
            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 * 1e3;  % km转换为m
        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:U
        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))^2);
        end
    end
end

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

doppler_shift = 100;  
time_delay = 1e-6;    
% 定义阵列响应角度 (假设值)
theta_x = 0.1;         % x方向角度
theta_y = 0.2;         % y方向角度
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);
% 计算 UPA 响应向量 v_k (使用克罗内克积)
v_k = kron(v_x, v_y);  % v_k的维度为 Nt^2 x 1

for m = 1:S
    for u = 1:U
        for t = 1:T
            % fprintf('m = %d, i = %d, n = %d, t = %d\n', m, i, n, t);
            g_m_i_n_t = gain(m, u, t);
            % 计算相位因子 e^{j 2 pi (doppler_shift * t - f_c * time_delay)}
            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;
        end
    end
end

h1 = h(:, 1, :, :);   % (Nt, U, T) 
H1 = zeros(Nt, Nt, T);  % 后面求tr用 卫星1的
for t = 1:T
    h1_t = h1(:, :, t);         % 维度 (Nt, U)
    H1(:, :, t) = h1_t * h1_t';  %  Nt*Nt
end

% 信道归一化             % 为什么预编码和信道的数量级总是相差悬殊 导致求迹为0
% magnitude_H1 = abs(H1);
% max_magnitude = max(magnitude_H1(:));
% if max_magnitude > 0
%     H1_normalized = H1 / max_magnitude;
% else
%     H1_normalized = H1; 
% end

H1 = H1 / norm(H1, 'fro');
H1 = H1 * 1e8;  % 可乘1e7到1e11

% 点波束调度
spotBeamCover = zeros(S, C, T);
spotBeamCover(:,:,1) = [
    1 0 1
];

% 点波束覆盖的用户总数 total_users_covered_by_spotbeam(s, t)
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;      % 带宽（MHz）
F = 10;        % 噪声系数
% 计算噪声功率（以瓦特为单位）
P_noise_dbw = k * noise_T * B * F;  
P_noise = 10^(P_noise_dbw / 10);  % 噪声功率P_noise
P_beam = 40 * ones(Nt, 1);  % 总功率

Bc2_prev = rand(S, C, U, T);
Cn2_prev = rand(S, C, U, T);
TAOc_prev = rand(S, C, U, T);
TAOn_prev = rand(S, C, U, T);

v_max_Uc = ones(Nt, S, C, T) / sqrt(Nt);    
v_max_Un = ones(Nt, S, C, U, T) / sqrt(Nt);    

cvx_solver mosek
cvx_begin sdp
    % -----先去除超公有流----------
    % variable Usc(Nt, Nt, S, C, U, T) semidefinite;  % C21
    variable Uc(Nt, Nt, S, C, T) semidefinite; 
    variable Un(Nt, Nt, S, C, U, T) semidefinite;

    % variable Qsc(S, C, U, T);
    % variable Rsc(S, C, U, T);  % R 为相关的 Q 取最小  不能按用户数均分(用这个)  目标函数里再均分  或者  均分 之后约束里乘用户数
    variable Qc(S, C, U, T);  % 存计算出的用户速率
    variable Rc(S, C, T);   % Qc 属于一个簇里的用户中的  存取最小值的速率  未均分的
    variable Rn(S, C, U, T);

    % variable Asc1(S, C, U, T); 
    % variable Asc2(S, C, U, T);
    variable Bc1(S, C, U, T);
    variable Bc2(S, C, U, T);
    variable Cn1(S, C, U, T);
    variable Cn2(S, C, U, T);

    % variable TAOsc(S, C, U, T);
    variable TAOc(S, C, U, T);
    variable TAOn(S, C, U, T);

    expression rateSum;
    expression f;  % 惩罚函数
    expression P_tot(Nt, S, T);

    % expression A2LHS(S, C, U, T);  % C24用
    expression B2LHS(S, C, U, T);  % C25用
    expression B2RHS(S, C, U, T);  % C29用
    expression C2LHS(S, C, U, T);  % C26用
    expression C2RHS(S, C, U, T);  % C31用

    expression auxiliary_var_37;
    expression auxiliary_var_38;

    expression common_fenmu(S, C, U, T);
    expression private_fenmu(S, C, U, T);


    for s = 1:S  % 优化目标和速率
        for t = 1:T
            covered_clusters = find(spotBeamCover(s, :, t) == 1);   
            for c = covered_clusters
                users = cluster_user_indices{c};
                for u = users
                    % fprintf('c = %d, u = %d\n', c, u);
                    rateSum = rateSum + Rc(s, c, t) / user_counts_per_cluster(c) + Rn(s, c, u, t);  % cvx中不能用+=
                end
            end
        end
    end

    for s = 1:S  % 迭代惩罚函数
        for t = 1:T
            covered_clusters = find(spotBeamCover(s, :, t) == 1);   
            for c = covered_clusters
                f = f + trace(Uc(:, :, s, c, t)) - v_max_Uc(:, s, c, t)' * Uc(:, :, s, c, t) * v_max_Uc(:, s, c, t);
                users = cluster_user_indices{c};
                for u = users
                    % fprintf('c = %d, u = %d\n', c, u);
                    f = f + trace(Un(:, :, s, c, u, t)) - v_max_Un(:, s, c, u, t)' * Un(:, :, s, c, u, t) * v_max_Un(:, s, c, u, t);
                end
            end
        end
    end

    for s = 1:S  % C23左侧 功率计算
        for t = 1:T
            covered_clusters = find(spotBeamCover(s, :, t) == 1);   
            for c = covered_clusters
                P_tot(:, s, t) = P_tot(:, s, t) + diag(Uc(:, :, s, c, t));
                users = cluster_user_indices{c};
                for u = users
                    % fprintf('c = %d, u = %d\n', c, u);
                    P_tot(:, s, t) = P_tot(:, s, t) + diag(Un(:, :, s, c, u, t));
                end
            end
        end
    end

    maximize(rateSum - 0.2 * f);

    subject to
        
        for s = 1:S  % C23 功率约束
            for t = 1:T
                P_tot(:, s, t) <= P_beam;
            end
        end
        
        % Bc1(S, C, U, T);  Bc2(S, C, U, T);   Rc(S, C, T);  % 已经取最小了  ----------------
        for s = 1:S  % C13
            for t = 1:T
                covered_clusters = find(spotBeamCover(s, :, t) == 1);   
                for c = covered_clusters
                    users = cluster_user_indices{c};
                    for u = users
                        % fprintf('c = %d, u = %d\n', c, u);
                        Bc1(s, c, u, t) - Bc2(s, c, u, t) >= Rc(s, c, t);
                    end
                end
            end
        end

        for s = 1:S  % C14
            for t = 1:T
                covered_clusters = find(spotBeamCover(s, :, t) == 1);   
                for c = covered_clusters
                    users = cluster_user_indices{c};
                    for u = users
                        % fprintf('c = %d, u = %d\n', c, u);
                        Cn1(s, c, u, t) - Cn2(s, c, u, t) >= Rn(s, c, u, t);
                    end
                end
            end
        end

        for s = 1:S  % C25左侧
            for t = 1:T  
                covered_clusters = find(spotBeamCover(s, :, t) == 1); 
                all_users = [];  % 用于存储调度的所有用户
                user_to_cluster_map = [];  % 用于记录用户与所属簇的映射
        
                % 收集调度的所有簇中的用户
                for c = covered_clusters
                    users = cluster_user_indices{c};  
                    all_users = [all_users, users];  
                    user_to_cluster_map = [user_to_cluster_map; [users(:), repmat(c, length(users), 1)]];
                end
        
                % 遍历每个簇
                for c = covered_clusters
                    users = cluster_user_indices{c};  % 当前簇的用户列表
                    other_clusters = setdiff(covered_clusters, c);  % 找到除当前簇以外的其他簇
                    
                    for u = users
                        for other_c = other_clusters
                            B2LHS(s, c, u, t) = B2LHS(s, c, u, t) + trace(H1 * Uc(:, :, s, other_c, t));
                            % B2LHS(s, c, u, t) = B2LHS(s, c, u, t) + 2;
                        end
                    end
                end

                for c = covered_clusters
                    users = cluster_user_indices{c};  % 当前簇的用户列表
                    for u = users
                        for j = covered_clusters
                            users = cluster_user_indices{j};  % 当前簇的用户列表
                            for k = users
                                B2LHS(s, c, u, t) = B2LHS(s, c, u, t) + trace(H1 * Un(:, :, s, j, k, t));
                            end
                        end
                        B2LHS(s, c, u, t) = B2LHS(s, c, u, t) + P_noise;
                    end
                end
            end
        end
        
        for s = 1:S  % C25
            for t = 1:T
                covered_clusters = find(spotBeamCover(s, :, t) == 1);   
                for c = covered_clusters
                    users = cluster_user_indices{c};
                    for u = users
                        % fprintf('c = %d, u = %d\n', c, u);
                        B2LHS(s, c, u, t) <= exp(Bc2_prev(s, c, u, t)) * (Bc2(s, c, u, t) - Bc2_prev(s, c, u, t) + 1);
                    end
                end
            end
        end

        for s = 1:S  % C26左侧
            for t = 1:T  
                covered_clusters = find(spotBeamCover(s, :, t) == 1); 
                all_users = [];  % 用于存储调度的所有用户
                user_to_cluster_map = [];  % 用于记录用户与所属簇的映射
        
                % 收集调度的所有簇中的用户
                for c = covered_clusters
                    users = cluster_user_indices{c};  
                    all_users = [all_users, users];  
                    user_to_cluster_map = [user_to_cluster_map; [users(:), repmat(c, length(users), 1)]];
                end
        
                % 遍历每个簇
                for c = covered_clusters
                    users = cluster_user_indices{c};  % 当前簇的用户列表
                    
                    for u = users
                        % 找到除用户 u 以外的其他用户
                        other_users_in_all_clusters = setdiff(all_users, u);
                        for other_user = other_users_in_all_clusters
                            % 利用 user_to_cluster_map 找到 other_user 所属的簇
                            other_user_cluster = user_to_cluster_map(user_to_cluster_map(:, 1) == other_user, 2);
                            C2LHS(s, c, u, t) = C2LHS(s, c, u, t) + trace(H1 * Un(:, :, s, other_user_cluster, other_user, t));
                            % C2LHS(s, c, u, t) = C2LHS(s, c, u, t) + 2;  %
                            % 没有累加成功是因为求trace这项一直为0 循环应该没问题 好像是预编码在1e-8数量级
                        end
                        C2LHS(s, c, u, t) = C2LHS(s, c, u, t) + P_noise;
                    end
                end
            end
        end
        
        for s = 1:S  % C26
            for t = 1:T
                covered_clusters = find(spotBeamCover(s, :, t) == 1);   
                for c = covered_clusters
                    users = cluster_user_indices{c};
                    for u = users
                        % fprintf('c = %d, u = %d\n', c, u);
                        C2LHS(s, c, u, t) <= exp(Cn2_prev(s, c, u, t)) * (Cn2(s, c, u, t) - Cn2_prev(s, c, u, t) + 1);
                    end
                end
            end
        end
        
        % C29 C31
        for s = 1:S  % C29右侧
            for t = 1:T  
                covered_clusters = find(spotBeamCover(s, :, t) == 1); 
                all_users = [];  % 用于存储调度的所有用户
                user_to_cluster_map = [];  % 用于记录用户与所属簇的映射
        
                % 收集调度的所有簇中的用户
                for c = covered_clusters
                    users = cluster_user_indices{c};  
                    all_users = [all_users, users];  
                    user_to_cluster_map = [user_to_cluster_map; [users(:), repmat(c, length(users), 1)]];
                end

                for c = covered_clusters
                    users = cluster_user_indices{c};  % 当前簇的用户列表
                    for u = users
                        for j = covered_clusters
                            B2RHS(s, c, u, t) = B2RHS(s, c, u, t) + trace(H1 * Uc(:, :, s, j, t));
                            users = cluster_user_indices{j};  % 当前簇的用户列表
                            for k = users
                                B2RHS(s, c, u, t) = B2RHS(s, c, u, t) + trace(H1 * Un(:, :, s, j, k, t));
                            end
                        end
                        B2RHS(s, c, u, t) = B2RHS(s, c, u, t) + P_noise;
                    end
                end
            end
        end
        
        for s = 1:S  % C29
            for t = 1:T
                covered_clusters = find(spotBeamCover(s, :, t) == 1);   
                for c = covered_clusters
                    users = cluster_user_indices{c};
                    for u = users
                        % fprintf('c = %d, u = %d\n', c, u);
                        TAOc(s, c, u, t) <= B2RHS(s, c, u, t);
                    end
                end
            end
        end

        for s = 1:S  % C31右侧
            for t = 1:T  
                covered_clusters = find(spotBeamCover(s, :, t) == 1); 
                all_users = [];  % 用于存储调度的所有用户
                user_to_cluster_map = [];  % 用于记录用户与所属簇的映射
        
                % 收集调度的所有簇中的用户
                for c = covered_clusters
                    users = cluster_user_indices{c};  
                    all_users = [all_users, users];  
                    user_to_cluster_map = [user_to_cluster_map; [users(:), repmat(c, length(users), 1)]];
                end

                for c = covered_clusters
                    users = cluster_user_indices{c};  % 当前簇的用户列表
                    for u = users
                        for j = covered_clusters
                            users = cluster_user_indices{j};  % 当前簇的用户列表
                            for k = users
                                C2RHS(s, c, u, t) = C2RHS(s, c, u, t) + trace(H1 * Un(:, :, s, j, k, t));
                            end
                        end
                        C2RHS(s, c, u, t) = C2RHS(s, c, u, t) + P_noise;
                    end
                end
            end
        end

        for s = 1:S  % C31
            for t = 1:T
                covered_clusters = find(spotBeamCover(s, :, t) == 1);   
                for c = covered_clusters
                    users = cluster_user_indices{c};
                    for u = users
                        % fprintf('c = %d, u = %d\n', c, u);
                        TAOn(s, c, u, t) <= C2RHS(s, c, u, t);
                    end
                end
            end
        end

        for s = 1:S  % C37-38
            for t = 1:T
                covered_clusters = find(spotBeamCover(s, :, t) == 1);   
                for c = covered_clusters
                    users = cluster_user_indices{c};
                    for u = users
                        % fprintf('c = %d, u = %d\n', c, u);
        
                        auxiliary_var_37 = log(TAOc_prev(s, c, u, t)) + 1;
                        norm([TAOc(s, c, u, t) + Bc1(s, c, u, t) - auxiliary_var_37, 2 * sqrt(TAOc_prev(s, c, u, t))], 2) ...
                            <= TAOc(s, c, u, t) - Bc1(s, c, u, t) + auxiliary_var_37;
        
                        auxiliary_var_38 = log(TAOn_prev(s, c, u, t)) + 1;
                        norm([TAOn(s, c, u, t) + Cn1(s, c, u, t) - auxiliary_var_38, 2 * sqrt(TAOn_prev(s, c, u, t))], 2) ...
                            <= TAOn(s, c, u, t) - Cn1(s, c, u, t) + auxiliary_var_38;
                    end
                end
            end
        end
        
        % 公有流SINR分母计算 无噪声
        for s = 1:S  % 
            for t = 1:T  
                covered_clusters = find(spotBeamCover(s, :, t) == 1); 
                all_users = [];  % 用于存储调度的所有用户
                user_to_cluster_map = [];  % 用于记录用户与所属簇的映射
        
                % 收集调度的所有簇中的用户
                for c = covered_clusters
                    users = cluster_user_indices{c};  
                    all_users = [all_users, users];  
                    user_to_cluster_map = [user_to_cluster_map; [users(:), repmat(c, length(users), 1)]];
                end
        
                % 遍历每个簇
                for c = covered_clusters
                    users = cluster_user_indices{c};  % 当前簇的用户列表
                    other_clusters = setdiff(covered_clusters, c);  % 找到除当前簇以外的其他簇
                    
                    for u = users
                        for other_c = other_clusters
                            common_fenmu(s, c, u, t) = common_fenmu(s, c, u, t) + trace(H1 * Uc(:, :, s, other_c, t));
                        end
                    end
                end

                for c = covered_clusters
                    users = cluster_user_indices{c};  % 当前簇的用户列表
                    for u = users
                        for j = covered_clusters
                            users = cluster_user_indices{j};  % 当前簇的用户列表
                            for k = users
                                common_fenmu(s, c, u, t) = common_fenmu(s, c, u, t) + trace(H1 * Un(:, :, s, j, k, t));
                            end
                        end
                    end
                end
            end
        end

        for s = 1:S  % 公有流SINR约束
            for t = 1:T
                covered_clusters = find(spotBeamCover(s, :, t) == 1);   
                for c = covered_clusters
                    users = cluster_user_indices{c};
                    for u = users
                        % fprintf('c = %d, u = %d\n', c, u);
                        trace(H1 * Uc(:, :, s, c, t)) - 10 * common_fenmu(s, c, u, t) >= 10 * P_noise;  % 100对应SINR
                    end
                end
            end
        end
        
        % 私有流SINR分母计算 无噪声
        for s = 1:S  % 
            for t = 1:T  
                covered_clusters = find(spotBeamCover(s, :, t) == 1); 
                all_users = [];  % 用于存储调度的所有用户
                user_to_cluster_map = [];  % 用于记录用户与所属簇的映射
        
                % 收集调度的所有簇中的用户
                for c = covered_clusters
                    users = cluster_user_indices{c};  
                    all_users = [all_users, users];  
                    user_to_cluster_map = [user_to_cluster_map; [users(:), repmat(c, length(users), 1)]];
                end
        
                % 遍历每个簇
                for c = covered_clusters
                    users = cluster_user_indices{c};  % 当前簇的用户列表
                    
                    for u = users
                        % 找到除用户 u 以外的其他用户
                        other_users_in_all_clusters = setdiff(all_users, u);
                        for other_user = other_users_in_all_clusters
                            % 利用 user_to_cluster_map 找到 other_user 所属的簇
                            other_user_cluster = user_to_cluster_map(user_to_cluster_map(:, 1) == other_user, 2);
                            private_fenmu(s, c, u, t) = private_fenmu(s, c, u, t) + trace(H1 * Un(:, :, s, other_user_cluster, other_user, t));
                            % C2LHS(s, c, u, t) = C2LHS(s, c, u, t) + 2;  %
                            % 没有累加成功是因为求trace这项一直为0 循环应该没问题 好像是预编码在1e-8数量级
                        end
                    end
                end
            end
        end
        
        for s = 1:S  % 私有流SINR约束
            for t = 1:T
                covered_clusters = find(spotBeamCover(s, :, t) == 1);   
                for c = covered_clusters
                    users = cluster_user_indices{c};
                    for u = users
                        % fprintf('c = %d, u = %d\n', c, u);
                        trace(H1 * Un(:, :, s, c, u, t)) - 20 * private_fenmu(s, c, u, t) >= 20 * P_noise;
                    end
                end
            end
        end
cvx_end

Thanks for all suggestions!

2

Looking at the MOSEK solution summary then I can see MOSEK reports an invalid solution.

Either you have a nasty problem or you hit a bug in MOSEK. Therefore, I suggest you upgrade to the latest MOSEK version 10.2.

You problem looks extremely badly scaled e.g.

1.0e-16 *
0.2628 -0.2126
-0.2126 0.2628

It makes no sense to work with such numbers when solving the problem using floating numbers IMO.

3

Thank you. And I have another question. That is … I’ve installed version 10.2 and the path to matlab points to version 10.2, but version 9.1.9 is still automatically called when solving. I don’t know why this happened.

4

Re-run cvx_setup
The output of cvx_solver will show all the versions of Mosek. Choose the version, such as Mosek_2, which corresponds to the Mosek version you want, by using
cvx_solver Mosek_2. You can do cvx_save_prefs to make this the default for future MATLAB sessions.

Or rename (or delete) the mosek directory (which is Mosek 9.1.9) under the CVX installation, then re-run cvx_setup. Then the Mosek 10.2 should be the only Mosek version, and you can just specify cvx_solver mosek`.

Or install CVX 2.2.2.

5

Thank you for your suggestions.

When I looked at my optimization problem again, I realized that there might be a problem with the notation of the constraints related to the exponential function, so I modified it.

If the symbol I modified is correct, the corresponding code is as follows.

cvx_clear;
rng(42);

S = 1;
C = 3; % 点波束2
U = 6;
T = 1;
% 天线数量
Ntx = 1;
Nty = 2;
Nt = Ntx * Nty;

S_height = 2000; % km 卫星 x y z
s1_pos = [106, 95];
offset = 4;
S1_pos = zeros(2, T); % 卫星1每个时隙的位置
S1_pos(:, 1) = s1_pos’;
% for t = 2:T
% S1_pos(:, t) = S1_pos(:, t-1) + offset;
% end

% 簇内用户位置 是否合适 -----------------------------------
C1_user_pos = [
173.67173899833338, 108.37931961155354;
174.57454953552102, 110.32493726720683
];
C2_user_pos = [
136.32860551422957, 42.16289711432755;
137.82744015490138, 49.40391946858323
];
C3_user_pos = [
% 44.000564580961424, 129.4559770874314;
% 43.8824784449596, 125.8989690461711
198.13023757449002, 156.0800008435594;
192.6303832555794, 157.87173209488702;
];

% 将所有簇的用户位置合并为一个cell array
user_positions = {
C1_user_pos, C2_user_pos, C3_user_pos
};
all_user_positions = vertcat(user_positions{:}); % 合并所有簇内用户位置
user_counts_per_cluster = [2, 2, 2]; % 各簇内用户数量
% 用户索引 cluster_user_indices
cluster_user_indices = cell(1, length(user_counts_per_cluster)); % cluster_user_indices{c}用它取出每组用户的索引
% 起始用户索引
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:U
for t = 1:T
% 获取卫星的位置，S1_pos或S2_pos，根据卫星编号
if m == 1
satellite_pos = S1_pos(:, t); % 卫星1的位置
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:U
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
% 较1220_02新增
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(Nt, S, U, T);

% 替换到138行
v = 7500;
lambda = c / f_c;
doppler_shift = v / lambda; % 多普勒频移 (假设值，根据具体情况设定)
time_delay = S_height / c; % 时延 (秒)，假设值

% 定义阵列响应角度 (假设值)
theta_x = 0.1; % x方向角度
theta_y = 0.2; % y方向角度

% 定义阵列响应向量函数
% v_d_k为一个复数向量，使用定义式 (4)
array_response = @(N, theta) 1/sqrt(N) * exp(-1j * pi * theta * (0:N-1)).';

% 计算 x 和 y 方向的阵列响应向量
v_x = array_response(Ntx, theta_x);
v_y = array_response(Nty, theta_y);

% 计算 UPA 响应向量 v_k (使用克罗内克积)
v_k = kron(v_x, v_y); % v_k的维度为 Nt^2 x 1

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)
g_m_i_n_t = gain_normalized(m, u, t);
% 计算相位因子 e^{j 2 pi (doppler_shift * t - f_c * time_delay)}
phase_factor = exp(1j * 2 * pi * (doppler_shift * t - f_c * time_delay));
% h(m, :, u, t) 代表 h 中第 m 个索引、第 u 个用户和第 t 个时间片的切片，提取了该维度中的所有天线元素。
% h(m, :, u, t) 的大小为 1 x Nt，也就是一个包含 Nt 个元素的行向量。
h(:, m, u, t) = g_m_i_n_t * phase_factor * v_k * 1e7;
end
end
end

h1 = h(:, 1, :, :); % (Nt, U, T) ------------------- 数量级1e-8
H1 = zeros(Nt, Nt, T); % 后面求tr用 卫星1的
for t = 1:T
h1_t = h1(:, :, t); % 维度 (Nt, U)
H1(:, :, t) = h1_t * h1_t’; % Nt*Nt
end

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

% H1 = H1 / norm(H1, ‘fro’);
% H1 = H1 * 1e8; % 可乘1e7到1e11

% 点波束调度
spotBeamCover = zeros(S, C, T);
spotBeamCover(:,:,1) = [
1 0 1
];

% 点波束覆盖的用户总数 total_users_covered_by_spotbeam(s, t)
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; % 带宽（MHz）
F = 10; % 噪声系数
% 计算噪声功率（以瓦特为单位）
P_noise_dbw = k * noise_T * B * F;
P_noise = 10^(P_noise_dbw / 10); % 噪声功率P_noise
P_beam = 30 * ones(Nt, 1); % 总功率
% P_beam = 30; % 总功率

% Bc2_prev = rand(S, C, U, T);
% Cn2_prev = rand(S, C, U, T);
% TAOc_prev = rand(S, C, U, T);
% TAOn_prev = rand(S, C, U, T);

Bc1_prev = rand(S, C, U, T);
Cn1_prev = rand(S, C, U, T);

TAOc_prev = 40 + (20 - 10) * rand(S, C, U, T);
TAOn_prev = 40 + (20 - 10) * rand(S, C, U, T);

v_max_Uc = ones(Nt, S, C, T) / sqrt(Nt);
v_max_Un = ones(Nt, S, C, U, T) / sqrt(Nt);

% mosek_path = ‘C:/Program Files/Mosek/10.2/toolbox/r2017a’;
% addpath(mosek_path);
% [res, rinfo] = mosekopt(‘version’);
% disp(rinfo);

cvx_solver mosek

cvx_begin sdp
% -----先去除超公有流----------
% variable Usc(Nt, Nt, S, C, U, T) semidefinite; % C21
variable Uc(Nt, Nt, S, C, T) semidefinite;
variable Un(Nt, Nt, S, C, U, T) semidefinite;

% variable Qsc(S, C, U, T);
% variable Rsc(S, C, U, T);  % R 为相关的 Q 取最小  不能按用户数均分(用这个)  目标函数里再均分  或者  均分 之后约束里乘用户数
variable Qc(S, C, U, T);  % 存计算出的用户速率
variable Rc(S, C, T);   % Qc 属于一个簇里的用户中的  存取最小值的速率  未均分的
variable Rn(S, C, U, T);

% variable Asc1(S, C, U, T); 
% variable Asc2(S, C, U, T);
variable Bc1(S, C, U, T);
variable Bc2(S, C, U, T);
variable Cn1(S, C, U, T);
variable Cn2(S, C, U, T);

% variable TAOsc(S, C, U, T);
variable TAOc(S, C, U, T);
variable TAOn(S, C, U, T);

expression rateSum;
expression f;  % 惩罚函数
expression P_tot(Nt, S, T);

% expression A2LHS(S, C, U, T);  % C24用
expression B2LHS(S, C, U, T);  % C25用
expression B2RHS(S, C, U, T);  % C29用
expression C2LHS(S, C, U, T);  % C26用
expression C2RHS(S, C, U, T);  % C31用

expression auxiliary_var_37;
expression auxiliary_var_38;

expression common_fenmu(S, C, U, T);
expression private_fenmu(S, C, U, T);


for s = 1:S  % 优化目标和速率
    for t = 1:T
        covered_clusters = find(spotBeamCover(s, :, t) == 1);   
        for c = covered_clusters
            users = cluster_user_indices{c};
            for u = users
                % fprintf('c = %d, u = %d\n', c, u);
                rateSum = rateSum + Rc(s, c, t) / user_counts_per_cluster(c) + Rn(s, c, u, t);  % cvx中不能用+=
            end
        end
    end
end

for s = 1:S  % 迭代惩罚函数
    for t = 1:T
        covered_clusters = find(spotBeamCover(s, :, t) == 1);   
        for c = covered_clusters
            f = f + trace(Uc(:, :, s, c, t)) - v_max_Uc(:, s, c, t)' * Uc(:, :, s, c, t) * v_max_Uc(:, s, c, t);
            users = cluster_user_indices{c};
            for u = users
                % fprintf('c = %d, u = %d\n', c, u);
                f = f + trace(Un(:, :, s, c, u, t)) - v_max_Un(:, s, c, u, t)' * Un(:, :, s, c, u, t) * v_max_Un(:, s, c, u, t);
            end
        end
    end
end

for s = 1:S  % C23左侧 功率计算
    for t = 1:T
        covered_clusters = find(spotBeamCover(s, :, t) == 1);   
        for c = covered_clusters
            P_tot(:, s, t) = P_tot(:, s, t) + diag(Uc(:, :, s, c, t));
            users = cluster_user_indices{c};
            for u = users
                % fprintf('c = %d, u = %d\n', c, u);
                P_tot(:, s, t) = P_tot(:, s, t) + diag(Un(:, :, s, c, u, t));
            end
        end
    end
end

% maximize(rateSum - 0.2 * f);

subject to  % 改公式 e的多少次方
    
    for s = 1:S  % C23 功率约束
        for t = 1:T
            P_tot(:, s, t) <= P_beam;
        end
    end
    
    % Bc1(S, C, U, T);  Bc2(S, C, U, T);   Rc(S, C, T);  % 已经取最小了  ----------------
    for s = 1:S  % C13
        for t = 1:T
            covered_clusters = find(spotBeamCover(s, :, t) == 1);   
            for c = covered_clusters
                users = cluster_user_indices{c};
                for u = users
                    % fprintf('c = %d, u = %d\n', c, u);
                    Bc1(s, c, u, t) - Bc2(s, c, u, t) >= Rc(s, c, t);
                end
            end
        end
    end

    for s = 1:S  % C14
        for t = 1:T
            covered_clusters = find(spotBeamCover(s, :, t) == 1);   
            for c = covered_clusters
                users = cluster_user_indices{c};
                for u = users
                    % fprintf('c = %d, u = %d\n', c, u);
                    Cn1(s, c, u, t) - Cn2(s, c, u, t) >= Rn(s, c, u, t);
                end
            end
        end
    end

    for s = 1:S  % C29右侧
        for t = 1:T  
            covered_clusters = find(spotBeamCover(s, :, t) == 1); 
            all_users = [];  % 用于存储调度的所有用户
            user_to_cluster_map = [];  % 用于记录用户与所属簇的映射
    
            % 收集调度的所有簇中的用户
            for c = covered_clusters
                users = cluster_user_indices{c};  
                all_users = [all_users, users];  
                user_to_cluster_map = [user_to_cluster_map; [users(:), repmat(c, length(users), 1)]];
            end
    
            % 遍历每个簇
            for c = covered_clusters
                users = cluster_user_indices{c};  % 当前簇的用户列表
                other_clusters = setdiff(covered_clusters, c);  % 找到除当前簇以外的其他簇
                
                for u = users
                    for other_c = other_clusters
                        B2RHS(s, c, u, t) = B2RHS(s, c, u, t) + trace(H1_normalized * Uc(:, :, s, other_c, t));
                        % B2LHS(s, c, u, t) = B2LHS(s, c, u, t) + 2;
                    end
                end
            end

            for c = covered_clusters
                users = cluster_user_indices{c};  % 当前簇的用户列表
                for u = users
                    for j = covered_clusters
                        users = cluster_user_indices{j};  % 当前簇的用户列表
                        for k = users
                            B2RHS(s, c, u, t) = B2RHS(s, c, u, t) + trace(H1_normalized * Un(:, :, s, j, k, t));
                        end
                    end
                    B2RHS(s, c, u, t) = B2RHS(s, c, u, t) + P_noise;
                end
            end
        end
    end
    
    for s = 1:S  % C25左侧
        for t = 1:T  
            covered_clusters = find(spotBeamCover(s, :, t) == 1); 
            all_users = [];  % 用于存储调度的所有用户
            user_to_cluster_map = [];  % 用于记录用户与所属簇的映射
    
            % 收集调度的所有簇中的用户
            for c = covered_clusters
                users = cluster_user_indices{c};  
                all_users = [all_users, users];  
                user_to_cluster_map = [user_to_cluster_map; [users(:), repmat(c, length(users), 1)]];
            end

            for c = covered_clusters
                users = cluster_user_indices{c};  % 当前簇的用户列表
                for u = users
                    for j = covered_clusters
                        B2LHS(s, c, u, t) = B2LHS(s, c, u, t) + trace(H1_normalized * Uc(:, :, s, j, t));
                        users = cluster_user_indices{j};  % 当前簇的用户列表
                        for k = users
                            B2LHS(s, c, u, t) = B2LHS(s, c, u, t) + trace(H1_normalized * Un(:, :, s, j, k, t));
                        end
                    end
                    B2LHS(s, c, u, t) = B2LHS(s, c, u, t) + P_noise;
                end
            end
        end
    end

    for s = 1:S  % C25
        for t = 1:T
            covered_clusters = find(spotBeamCover(s, :, t) == 1);   
            for c = covered_clusters
                users = cluster_user_indices{c};
                for u = users
                    % fprintf('c = %d, u = %d\n', c, u);
                    B2LHS(s, c, u, t) <= exp(Bc1_prev(s, c, u, t)) * (Bc1(s, c, u, t) - Bc1_prev(s, c, u, t) + 1);
                end
            end
        end
    end

    for s = 1:S  % C31右侧
        for t = 1:T  
            covered_clusters = find(spotBeamCover(s, :, t) == 1); 
            all_users = [];  % 用于存储调度的所有用户
            user_to_cluster_map = [];  % 用于记录用户与所属簇的映射
    
            % 收集调度的所有簇中的用户
            for c = covered_clusters
                users = cluster_user_indices{c};  
                all_users = [all_users, users];  
                user_to_cluster_map = [user_to_cluster_map; [users(:), repmat(c, length(users), 1)]];
            end
    
            % 遍历每个簇
            for c = covered_clusters
                users = cluster_user_indices{c};  % 当前簇的用户列表
                
                for u = users
                    % 找到除用户 u 以外的其他用户
                    other_users_in_all_clusters = setdiff(all_users, u);
                    for other_user = other_users_in_all_clusters
                        % 利用 user_to_cluster_map 找到 other_user 所属的簇
                        other_user_cluster = user_to_cluster_map(user_to_cluster_map(:, 1) == other_user, 2);
                        C2LHS(s, c, u, t) = C2LHS(s, c, u, t) + trace(H1_normalized * Un(:, :, s, other_user_cluster, other_user, t));
                        % C2LHS(s, c, u, t) = C2LHS(s, c, u, t) + 2;  %
                        % 没有累加成功是因为求trace这项一直为0 循环应该没问题 好像是预编码在1e-8数量级
                    end
                    C2LHS(s, c, u, t) = C2LHS(s, c, u, t) + P_noise;
                end
            end
        end
    end
    
    for s = 1:S  % C26左侧
        for t = 1:T  
            covered_clusters = find(spotBeamCover(s, :, t) == 1); 
            all_users = [];  % 用于存储调度的所有用户
            user_to_cluster_map = [];  % 用于记录用户与所属簇的映射
    
            % 收集调度的所有簇中的用户
            for c = covered_clusters
                users = cluster_user_indices{c};  
                all_users = [all_users, users];  
                user_to_cluster_map = [user_to_cluster_map; [users(:), repmat(c, length(users), 1)]];
            end

            for c = covered_clusters
                users = cluster_user_indices{c};  % 当前簇的用户列表
                for u = users
                    for j = covered_clusters
                        users = cluster_user_indices{j};  % 当前簇的用户列表
                        for k = users
                            C2RHS(s, c, u, t) = C2RHS(s, c, u, t) + trace(H1_normalized * Un(:, :, s, j, k, t));
                        end
                    end
                    C2RHS(s, c, u, t) = C2RHS(s, c, u, t) + P_noise;
                end
            end
        end
    end

    for s = 1:S  % C26
        for t = 1:T
            covered_clusters = find(spotBeamCover(s, :, t) == 1);   
            for c = covered_clusters
                users = cluster_user_indices{c};
                for u = users
                    % fprintf('c = %d, u = %d\n', c, u);
                    C2RHS(s, c, u, t) <= exp(Cn1_prev(s, c, u, t)) * (Cn1(s, c, u, t) - Cn1_prev(s, c, u, t) + 1);
                end
            end
        end
    end
    
    % C29 C31
    
    
    for s = 1:S  % C29
        for t = 1:T
            covered_clusters = find(spotBeamCover(s, :, t) == 1);   
            for c = covered_clusters
                users = cluster_user_indices{c};
                for u = users
                    % fprintf('c = %d, u = %d\n', c, u);
                    TAOc(s, c, u, t) <= B2RHS(s, c, u, t);
                end
            end
        end
    end

    

    for s = 1:S  % C31
        for t = 1:T
            covered_clusters = find(spotBeamCover(s, :, t) == 1);   
            for c = covered_clusters
                users = cluster_user_indices{c};
                for u = users
                    % fprintf('c = %d, u = %d\n', c, u);
                    TAOn(s, c, u, t) <= C2LHS(s, c, u, t);
                end
            end
        end
    end

    for s = 1:S  % C37-38
        for t = 1:T
            covered_clusters = find(spotBeamCover(s, :, t) == 1);   
            for c = covered_clusters
                users = cluster_user_indices{c};
                for u = users
                    % fprintf('c = %d, u = %d\n', c, u);
    
                    auxiliary_var_37 = log(TAOc_prev(s, c, u, t)) + 1;
                    norm([TAOc(s, c, u, t) + Bc2(s, c, u, t) - auxiliary_var_37, 2 * sqrt(TAOc_prev(s, c, u, t))], 2) ...
                        <= TAOc(s, c, u, t) - Bc2(s, c, u, t) + auxiliary_var_37;
    
                    auxiliary_var_38 = log(TAOn_prev(s, c, u, t)) + 1;
                    norm([TAOn(s, c, u, t) + Cn2(s, c, u, t) - auxiliary_var_38, 2 * sqrt(TAOn_prev(s, c, u, t))], 2) ...
                        <= TAOn(s, c, u, t) - Cn2(s, c, u, t) + auxiliary_var_38;
                end
            end
        end
    end
    
    % 公有流SINR分母计算 无噪声
    for s = 1:S  % 
        for t = 1:T  
            covered_clusters = find(spotBeamCover(s, :, t) == 1); 
            all_users = [];  % 用于存储调度的所有用户
            user_to_cluster_map = [];  % 用于记录用户与所属簇的映射
    
            % 收集调度的所有簇中的用户
            for c = covered_clusters
                users = cluster_user_indices{c};  
                all_users = [all_users, users];  
                user_to_cluster_map = [user_to_cluster_map; [users(:), repmat(c, length(users), 1)]];
            end
    
            % 遍历每个簇
            for c = covered_clusters
                users = cluster_user_indices{c};  % 当前簇的用户列表
                other_clusters = setdiff(covered_clusters, c);  % 找到除当前簇以外的其他簇
                
                for u = users
                    for other_c = other_clusters
                        common_fenmu(s, c, u, t) = common_fenmu(s, c, u, t) + trace(H1_normalized * Uc(:, :, s, other_c, t));
                    end
                end
            end

            for c = covered_clusters
                users = cluster_user_indices{c};  % 当前簇的用户列表
                for u = users
                    for j = covered_clusters
                        users = cluster_user_indices{j};  % 当前簇的用户列表
                        for k = users
                            common_fenmu(s, c, u, t) = common_fenmu(s, c, u, t) + trace(H1_normalized * Un(:, :, s, j, k, t));
                        end
                    end
                end
            end
        end
    end
   
    for s = 1:S  % 公有流SINR约束
        for t = 1:T
            covered_clusters = find(spotBeamCover(s, :, t) == 1);   
            for c = covered_clusters
                users = cluster_user_indices{c};
                for u = users
                    % fprintf('c = %d, u = %d\n', c, u);
                    trace(H1_normalized * Uc(:, :, s, c, t)) - 0.9 * common_fenmu(s, c, u, t) >= 0.9 * P_noise;  % 100对应SINR
                end
            end
        end
    end
    
    % 私有流SINR分母计算 无噪声
    for s = 1:S  % 
        for t = 1:T  
            covered_clusters = find(spotBeamCover(s, :, t) == 1); 
            all_users = [];  % 用于存储调度的所有用户
            user_to_cluster_map = [];  % 用于记录用户与所属簇的映射
    
            % 收集调度的所有簇中的用户
            for c = covered_clusters
                users = cluster_user_indices{c};  
                all_users = [all_users, users];  
                user_to_cluster_map = [user_to_cluster_map; [users(:), repmat(c, length(users), 1)]];
            end
    
            % 遍历每个簇
            for c = covered_clusters
                users = cluster_user_indices{c};  % 当前簇的用户列表
                
                for u = users
                    % 找到除用户 u 以外的其他用户
                    other_users_in_all_clusters = setdiff(all_users, u);
                    for other_user = other_users_in_all_clusters
                        % 利用 user_to_cluster_map 找到 other_user 所属的簇
                        other_user_cluster = user_to_cluster_map(user_to_cluster_map(:, 1) == other_user, 2);
                        private_fenmu(s, c, u, t) = private_fenmu(s, c, u, t) + trace(H1_normalized * Un(:, :, s, other_user_cluster, other_user, t));
                        % C2LHS(s, c, u, t) = C2LHS(s, c, u, t) + 2;  %
                        % 没有累加成功是因为求trace这项一直为0 循环应该没问题 好像是预编码在1e-8数量级
                    end
                end
            end
        end
    end
    
    for s = 1:S  % 私有流SINR约束
        for t = 1:T
            covered_clusters = find(spotBeamCover(s, :, t) == 1);   
            for c = covered_clusters
                users = cluster_user_indices{c};
                for u = users
                    % fprintf('c = %d, u = %d\n', c, u);
                    trace(H1_normalized * Un(:, :, s, c, u, t)) - 0.2 * private_fenmu(s, c, u, t) >= 0.2 * P_noise;
                end
            end
        end
    end

cvx_end

I have a problem that after I’ve modified it: with an optimization objective, the solution state is unbounded, positive infinity, while it can be solved without the objective function, but the optimal value is 0, so I want to know how I can obtain a reasonable rate value instead of positive infinity.

Thank you very much.

6