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!