My simulation results are inconsistent with my expectation

[CVX_user12]

clear all; clc;
% ========== System Parameters Initialization ==========
I = 3; % Number of UAVs
J = 2; % Number of RSUs
M = 2; % Number of vehicles
T = 5; % Time steps

D = 4e4 * ones(T, I); % Data volume(3e4-5e4)4e4
S = 1e3 * ones(T, I); % Task size(1000)

sigma2 = 1e-14; % Noise power(W)(-110dBm)
B_mj = 4e5 * ones(J, M); % Bandwidth(2MHz)
Gamma_L = 0.5; % Utility parameter of completing low-priority task on time
Gamma_H = 1; % Utility parameter of completing low-priority task past the time
c = 0.1; % Attenuation coefficient
p_max = 1; % Maximum transmit power（1W）
eta_max = 1e7;% Resource allocation upper limit(100-150MHz)
max_iter = 20; % Maximum iterations
tol = 1e-3; % Convergence tolerance

delta = 7e6 * ones(I, 1); % UAV local computing resources (CPU Hz)(6-8MHz)

P = 0.5 * ones(T,I,J);% Line-of-sight probability

ratio = 0.42;%Constant

r_ij = 3.226e6 * ones(T,I,J);% Transmission power between UAV i and RSU j
r_im = 9.678e6 * ones(T,I,J,M);% Transmission power between UAV i and vehicle m
g_mj = 7.14e-9 * ones(J,M);% Channel gain between vehicle m and RSU j(1/(1.4e8))

tau_i = zeros(T,I);
for t = 1:T
    for i = 1: I
        tau_i(t,i) =ratio * D(t,i) * S(t,i)/delta(i); % Task completion time threshold
    end
end
% pri = randi([0,1], T,I); % Task priority (0: low, 1: high)
pri = [0 1 1
    1 0 0
    1 0 1
    0 1 0
    1 1 0];%Priorities of tasks


beta = 1e3;%max(tau_i(:)) * 2;%1e3; % Big-M constant

% ========== Proposed method ==========

% ========== Initialize Variables ==========
theta1_prev = [0 0.8 0.8
    0.8 0 0
    0.8 0 0.8
    0 0.8 0
    0.8 0.8 0]; % Initial task offloading ratio

T_1_prev = tau_i/ratio;
% eta1_prev = eta_max/I *ones(T,I,J);
eta1_prev(:,:,1)= [1e-6 1 1e-6
    1 1e-6 1e-6
    1 1e-6 1e-6
    1e-6 1 1e-6
    1 1e-6 1e-6]* eta_max;
eta1_prev(:,:,2)=[
    1e-6 1e-6 1
    1e-6 1e-6 1e-6
    1e-6 1e-6 1
    1e-6 1e-6 1e-6
    1e-6 1 1e-6] * eta_max;
% p_mi1_prev = p_max/I * ones(T,I,M);
p_mi1_prev(:,:,1)= [1e-6 1 1e-6
    1 1e-6 1e-6
    1 1e-6 1e-6
    1e-6 1 1e-6
    1 1e-6 1e-6]*p_max;
p_mi1_prev(:,:,2)=[1e-6 1e-6 1
    1e-6 1e-6 1e-6
    1e-6 1e-6 1
    1e-6 1e-6 1e-6
    1e-6 1 1e-6] *p_max;

%Each UAV corresponds to one vehicle and one RSU.
% R1_prev = zeros(T,I,J);
% v1_prev = zeros(T,I,M);
% for t = 1:T
%     for i = 1:I
%         rand_j = randi(J);
%         rand_m = randi(M);
%         R1_prev(t,i,rand_j) = 1;
%         v1_prev(t,i,rand_m) = 1;
%     end
% end
R1_prev(:,:,1) = [0 1 0
    1 0 0
    1 0 0
    0 1 0
    1 0 0];
R1_prev(:,:,2) =[
    0 0 1
    0 0 0
    0 0 1
    0 0 0
    0 1 0];
v1_prev(:,:,1) =[0 1 0
    1 0 0
    1 0 0
    0 1 0
    1 0 0];
v1_prev(:,:,2) =[0 0 1
    0 0 0
    0 0 1
    0 0 0
    0 1 0];

%Record the target value for each iteration
iter_utility = zeros(max_iter, 1);
converged = false;

%Record the optimal solution
best_utility = -inf;
best_theta1 = [];
best_T_1 = [];
best_eta1 = [];
best_p_mi1 = [];
best_R1 = [];
best_v1 = [];
% ========== SCA Iteration ==========
for k = 1:max_iter
    fprintf('Iteration %d\n', k);

    % --- Build convex optimization problem for current iteration ---
    cvx_solver Mosek
    cvx_begin
    variable theta1(T,I) % Optimization variable: task offloading ratio

    variable eta1(T,I,J) % RSU resource allocation
    variable p_mi1(T,I,M) % Vehicle transmit power
    variable T_1(T,I);% Time required for task ki(t)
    variable z1(T, I) binary;% Big-M auxiliary variable
    variable u1(T, I);

    variable R1(T,I,J) binary;% The allocation of RSUs
    variable v1(T,I,M) binary;% The allocation of vehicels
    % --- Objective function definition ---
    total_utility1 = 0;
    for t = 1:T
        for i = 1:I
            term_local = (1 - theta1(t,i)) * S(t,i) * D(t,i) / delta(i);
            T_1(t,i) >= term_local;

            for j = 1:J
                for m = 1:M
                    weight = R1_prev(t,i,j) * v1_prev(t,i,m);
                    % Transmission time of task from vehicle m to RSU j
                    C = weight *((theta1_prev(t,i)*D(t,i))/(B_mj(j, m) * (log(1 + p_mi1_prev(t,i,m) * g_mj(j,m) / sigma2)/log(2)))+ ...
                        ((theta1(t,i) - theta1_prev(t,i))*D(t,i))/(B_mj(j, m) * (log(1 + p_mi1_prev(t,i,m) * g_mj(j,m) / sigma2)/log(2)))- ...
                        ((theta1_prev(t,i)*D(t,i)* g_mj(j,m))/(log(2) * B_mj(j, m) * sigma2 * (1 + p_mi1_prev(t,i,m) * g_mj(j,m) / sigma2) * (log(1 + p_mi1_prev(t,i,m) * g_mj(j,m) / sigma2)/log(2))^2) ) * (p_mi1(t,i,m) - p_mi1_prev(t,i,m)));

                    % Mathematical expectation of transmission time from UAV i to RSU j for task
                    G = weight *(P(t,i,j) * ((theta1(t,i) * D(t,i))/r_ij(t,i,j)) + (1 - P(t,i,j)) * (((theta1(t,i) * D(t,i))/r_im(t,i,j,m)) + C));

                    % Computation time of task in RSU j
                    E = weight *((theta1_prev(t,i)*S(t,i)*D(t,i))/eta1_prev(t,i,j) + ...
                        (theta1(t,i) - theta1_prev(t,i))*S(t,i)*D(t,i)/eta1_prev(t,i,j) - ...
                        (theta1_prev(t,i)*S(t,i)*D(t,i))/(eta1_prev(t,i,j)^2)*(eta1(t,i,j) - eta1_prev(t,i,j)));

                    T_1(t,i) >= G + E;
                end
            end

            % Calculate task utility based on big-M method
            % T_1(t,i)<= tau_i
            f1 = pri(t,i)*(log(1 + tau_i(t,i) - T_1(t,i))/log(2) + Gamma_L) + (1 - pri(t,i))*Gamma_L;

            % T_1(t,i)> tau_i
            f2 = -pri(t,i)*Gamma_H + (1-pri(t,i))*Gamma_L * exp(-c*(T_1_prev(t,i) - tau_i(t,i))) ...
                - c*(1-pri(t,i))*Gamma_L * exp(-c*(T_1_prev(t,i) - tau_i(t,i))) * (T_1(t,i) - T_1_prev(t,i));

            u1(t,i) <= f1 + beta * (1-z1(t,i));

            u1(t,i) <= f2 + beta * z1(t,i);

            total_utility1 = total_utility1 + u1(t,i);
        end
    end
    maximize(total_utility1)

    % --- Constraints ---
    subject to
    % Variable range constraints
    0 <= theta1 <= 1;% Constraint C2

    1e-6 <= eta1 <= eta_max;% Constraint C3

    for t = 1:T
        for j = 1:J
            sum(eta1(t,:,j)) <= eta_max * sum(R1(t,:,j));% Constraint C4：only the RSUs that have been assigned tasks perform computation
            sum(eta1(t,:,j)) <= eta_max
        end
    end


    for t = 1:T
        for i = 1:I
            sum(R1(t,i,:)) <= 1; %Constraint C6：Each task is assigned to one RSU at most
            for j = 1:J
                eta1(t,i,j) <= eta_max * R1(t,i,j) + 1e-6; %Only the RSUs that have been assigned tasks perform computation
            end
        end
    end

    1e-6 <= p_mi1 <= p_max;% Constraint C7


    for t = 1:T
        for m = 1:M
            sum(p_mi1(t,:,m)) <= p_max * sum(v1(t,:,m));% Constraint C8:only the vehicles that have been assigned tasks perform transmission
            sum(p_mi1(t,:,m)) <= p_max
        end
    end

    for t = 1:T
        for i = 1:I
            sum(v1(t,i,:)) <= 1; %Constraint C10：Each task is assigned to one vehicle at most
            for m = 1:M
                p_mi1(t,i,m) <= p_max * v1(t,i,m)+1e-6; % Only the vehicles that have been assigned tasks perform transmission
            end
        end
    end

    % Big-M method
    for t = 1:T
        for i = 1:I
            T_1(t,i) <= tau_i(t,i) + beta * (1 - z1(t,i));
            T_1(t,i) >= tau_i(t,i) - beta * z1(t,i);

        end
    end

    % Unload constraints - ensure consistency between theta 1 and R1/v1
    for t = 1:T
        for i = 1:I
            % When theta 1 (t, i)=0, no RSUs or vehicles are assigned
            sum(R1(t,i,:)) <= beta * theta1(t,i);
            sum(v1(t,i,:)) <= beta * theta1(t,i);

            % When theta 1 (t, i)>0, an RSU and a vehicle must be assigned
            sum(R1(t,i,:)) >= theta1(t,i);
            sum(v1(t,i,:)) >= theta1(t,i);

        end
    end
    cvx_end

    % Save the target value of the current iteration
    iter_utility(k) = cvx_optval;

    % Check the optimization status
    if strcmp(cvx_status, 'Solved') || strcmp(cvx_status, 'Optimal')
        fprintf('Optimization status: %s\n', cvx_status);

        % Save the target value of the current iteration
        iter_utility(k) = cvx_optval;

        % Calculate the relative rate of change
        if k > 1 && ~isnan(iter_utility(k-1)) && iter_utility(k-1) ~= 0
            rel_change = abs(iter_utility(k) - iter_utility(k-1)) / abs(iter_utility(k-1));
            fprintf('Relative change in objective: %.6f\n', rel_change);

            % Convergence check
            if rel_change < tol
                converged = true;
                fprintf('Converged after %d iterations\n', k);
                break;
            end
        end

        % Update variable values only when the results are valid
        if ~(anynan(theta1(:)) || anynan(T_1(:)) || anynan(eta1(:)) || anynan(p_mi1(:)) || anynan(v1(:)) || anynan(R1(:)))
            if cvx_optval > best_utility
                best_utility = cvx_optval;
                best_theta1 = theta1;
                best_T_1 = T_1;
                best_eta1 = eta1;
                best_p_mi1 = p_mi1;
                best_R1 = R1;
                best_v1 = v1;
            end

            theta1_prev = theta1;
            T_1_prev = T_1;
            eta1_prev = eta1;
            p_mi1_prev = p_mi1;
            v1_prev = v1;
            R1_prev = R1;
        else
            fprintf('Warning: NaN values detected in optimization results, skipping update\n');
            % Return to the historical optimal solution
            if ~isempty(best_theta1)
                theta1_prev = best_theta1;
                T_1_prev = best_T_1;
                eta1_prev = best_eta1;
                p_mi1_prev = best_p_mi1;
                v1_prev = best_v1;
                R1_prev = best_R1;
            end
        end
    else
        fprintf('Optimization failed with status: %s\n', cvx_status);

    end

    fprintf('Current utility: %.6f\n\n', cvx_optval);
end

% ========== Output Results ==========
fprintf('Proposed method total utility: %f\n', cvx_optval);

In the simulation results, for a low-priority task (where pri(t,i) = 0), 58.32% is offloaded to a RSU (denoted by j), and the computing resources allocated to it (eta1(t,i,j) and eta1_prev(t,i,j)) are approximately 0. According to the code for “E”, its computation time at RSU j should be large. And given the constraints “T_1(t,i) >= G + E” and “T_1(t,i) >= term_local;”, the final execution time of the task, T_1(t,i), should also be large.

However, in reality, T_1(t,i) is approximately equal to term_local—a small constant compared to the “E” I computed. This suggests that the actual value of G + E is nearly 0? Could there be an issue with my code or calculations?

[Mark_L_Stone]

It’s your model. I don’t think you can expect CVX forum readers to understand how to model whatever it is you are modeling.

Moreover, between poorly scaled input data and use of SCA, we have no idea whether SCA is converging to anything lest alone a global or even local optimum of your original problem.

The comments I provided to your earlier questions also apply to this question.

Why these variables are optimized to zero?

I don’t sufficiently understand your model to have an opinion as to whether your opinion is correct. But if your opinion is correct, can you produce values of the variables which your opinion says should exist, and which produce a better objective value. If not, perhaps your model is not suited for your purpose, and you need to figure out how to change it to make it so.

Algorithm cannot converge

SCA is generally unreliable, and might not converge to anything, let alone a global or even local optimum of the original problem. You should try to improve numerical scaling of input data, for instance by changing units, so that all non-zero input data is within a small number of orders of magnitude of 1. That might help your overall SCA performance, but not necessarily. Different (better) initialization values for the “SCA” variables might help, but not necessarily.

[CVX_user12]

I understand, thank you for your guidance.