Is there anything wrong with my code or my design?

1 
clear; clc;
% ========== System Parameters Initialization ==========
I = 3; % Number of UAVs
J = 2; % Number of RSUs
M = 2; % Number of vehicles
T = 5; % Time steps
D = randi([4e6,18e6], T, I); % Data volume
S = randi([4e2,1e3], T, I); % Task size
sigma2 = 1e-9; % Noise power
B_mj = 1e6 * ones(J, M); % Bandwidth
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 = 0.8; % Maximum transmit power
eta_max = 2e10;% Resource allocation upper limit
max_iter = 20; % Maximum iterations
tol = 1e-3; % Convergence tolerance

delta = 1e9* (1 + 0.2 * rand(I, 1)); % UAV local computing resources (CPU Hz)
P = rand(T,I,J);% Line-of-sight probability
r_ij = 0.06e6 + (2.34e6 - 0.06e6) * rand(T,I,J);% Transmission power between UAV i and RSU j
r_im = 0.06e6 + (2.34e6 - 0.06e6) * rand(T,I,J,M);% Transmission power between UAV i and vehicle m
g_mj = 0.02 * ones(J,M);% Channel gain between vehicle m and RSU j
tau_i = zeros(T,I);
for t = 1:T
for i = 1: I
tau_i(t,i) = D(t,i) * S(t,i)/delta(i); % Task completion time threshold
end
end
v = randi([0,1], T,I); % Task priority (0: low, 1: high)

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


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

% ========== Initialize Variables ==========
theta1_prev = 0.5 * ones(T,I); % Initial task offloading ratio
for t = 1:T
for i = 1:I
for j =1:J
for m =1 : M
z_mj1_prev(t,i,j,m) = (log(1 + (p_max * g_mj(j,m) / sigma2)) ) / log(2); % Initial auxiliary variable value
end
end

end
end
T_1_prev = tau_i;
eta1_prev = 0.5 * eta_max *ones(T,I,J);
for t = 1:T
for i = 1:I
for j = 1:J
lambda1(t,i,j) = (theta1_prev(t,i) * S(t,i) * D(t,i)) / eta1_prev(t,i,j);
end
end
end
% ========== SCA Iteration ==========
for k = 1:max_iter

% --- Build convex optimization problem for current iteration ---
cvx_begin quiet
variable theta1(T,I) % Optimization variable: task offloading ratio
variable z_mj1(T,I,J,M) % Auxiliary variable
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);% Big-M auxiliary variable
variable u1(T, I);

variable T_off1(T, I);
% --- 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
% Calculate Taylor-expanded C term
C = (theta1_prev(t,i)*D(t,i))/(B_mj(j, m)*z_mj1_prev(t,i,j,m)) + ...
(theta1(t,i) * D(t,i))/(B_mj(j, m)*z_mj1_prev(t,i,j,m)) - ...
(theta1_prev(t,i)*D(t,i))/(B_mj(j, m)*(z_mj1_prev(t,i,j,m))^2) * z_mj1(t,i,j,m);

G = P(t,i,j) * ((theta1(t,i) * D(t,i))*inv_pos(r_ij(t,i,j))) + (1 - P(t,i,j)) * (((theta1(t,i) * D(t,i))*inv_pos(r_im(t,i,j,m))) + C);

T_off1(t, i) <= G + theta1(t,i)*S(t,i)*D(t,i) - lambda1(t,i,j)*eta1(t,i,j);% Minimum offloading time

end
end
T_1(t,i) >= T_off1(t,i);

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

% T_1(t,i)> tau_i
f2 = -v(t,i)*Gamma_H + (1-v(t,i))*Gamma_L * exp(-c*(T_1_prev(t,i) - tau_i(t,i))) ...
- c*(1-v(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
0 <= eta1 <= eta_max;% Constraint C3
% Constraint C4
for t = 1:T
for j = 1:J
sum(eta1(t,:,j)) <= eta_max;
end
end
0 <= p_mi1 <= p_max;% Constraint C5
% Constraint C6
for t = 1:T
for m = 1:M
sum(p_mi1(t,:,m)) <= p_max;
end
end

% Auxiliary variable constraints
for t = 1:T
for i = 1:I
for j =1:J
for m =1 : M
z_mj1(t,i,j,m) <= (log(1 + (p_mi1(t,i,m) * g_mj(j,m) / sigma2)) ) / log(2);
end
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);
% Relaxation
0<= z1(t,i) <= 1;
end
end
cvx_end

% --- Check convergence ---
if norm(theta1 - theta1_prev, 'fro') < tol
break;
end
% --- Save previous iteration variable values ---
theta1_prev = theta1;
z_mj1_prev = z_mj1;
T_1_prev = T_1;
eta1_prev = eta1;


% --- Update Dinkelbach parameters ---
for t = 1:T
for i = 1:I
for j = 1:J
lambda1(t,i,j) = (theta1(t,i) * S(t,i) * D(t,i)) / eta1(t,i,j);
end
end
end
end

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


% ========== Traditional method ==========

% ========== Initialize Variables ==========
theta2_prev = 0.5 * ones(T,I); % Initial task offloading ratio
for t = 1:T
for i = 1:I
for j =1:J
for m =1 : M
z_mj2_prev(t,i,j,m) = (log(1 + (p_max * g_mj(j,m) / sigma2)) ) / log(2); % Initial auxiliary variable value
end
end
end
end

T_2_prev = tau_i;
eta2_prev = 0.5 * eta_max *ones(T,I,J);
for t = 1:T
for i = 1:I
for j = 1:J
lambda2(t,i,j) = (theta2_prev(t,i) * S(t,i) * D(t,i)) / eta2_prev(t,i,j);
end
end
end
% ========== SCA Iteration ==========
for k = 1:max_iter

% --- Build convex optimization problem for current iteration ---
cvx_begin quiet
variable theta2(T,I) % Optimization variable: task offloading ratio
variable z_mj2(T,I,J,M) % Auxiliary variable
variable eta2(T,I,J) % RSU resource allocation
variable p_mi2(T,I,M) % Vehicle transmit power
variable T_2(T,I);% Time required for task ki(t)
variable z2(T, I);% Big-M auxiliary variable
variable u2(T, I);

variable T_off2(T, I);
% --- Objective function definition ---
total_utility2 = 0;
for t = 1:T
for i = 1:I
term_local = (1 - theta2(t,i)) * S(t,i) * D(t,i) / delta(i);
T_2(t,i) >= term_local;
for j = 1:J
for m = 1:M
% Calculate Taylor-expanded C term
C = (theta2_prev(t,i)*D(t,i))/(B_mj(j, m)*z_mj2_prev(t,i,j,m)) + ...
(theta2(t,i) * D(t,i))/(B_mj(j, m)*z_mj2_prev(t,i,j,m)) - ...
(theta2_prev(t,i)*D(t,i))/(B_mj(j, m)*(z_mj2_prev(t,i,j,m))^2) * z_mj2(t,i,j,m);

G = P(t,i,j) * ((theta2(t,i) * D(t,i))*inv_pos(r_ij(t,i,j))) + (1 - P(t,i,j)) * (((theta2(t,i) * D(t,i))*inv_pos(r_im(t,i,j,m))) + C);

T_off2(t, i) <= G + theta2(t,i)*S(t,i)*D(t,i) - lambda2(t,i,j)*eta2(t,i,j);% Minimum offloading time

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

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

u2(t,i) <= f1 + beta * (1-z2(t,i));

u2(t,i) <= f2 + beta * z2(t,i);

total_utility2 = total_utility2 + u2(t,i);
end
end
maximize(total_utility2)

% --- Constraints ---
subject to

% Constraint C2
0 <= theta2 <= 1;
% Constraint C3
0 <= eta2 <= eta_max;
% Constraint C4
for t = 1:T
for j = 1:J
sum(eta2(t,:,j)) <= eta_max;
end
end
% Constraint C5
0 <= p_mi2 <= p_max;
% Constraint C6
for t = 1:T
for m = 1:M
sum(p_mi2(t,:,m)) <= p_max;
end
end

% Auxiliary variable constraints
for t = 1:T
for i = 1:I
for j =1:J
for m =1 : M
z_mj2(t,i,j,m) <= (log(1 + (p_mi2(t,i,m) * g_mj(j,m) / sigma2)) ) / log(2);
end
end

end
end

% Big-M method
for t = 1:T
for i = 1:I
T_2(t,i) <= tau_i(t,i) + beta * (1 - z2(t,i));
T_2(t,i) >= tau_i(t,i) - beta * z2(t,i);
% Relaxation
0<= z2(t,i) <= 1;
end
end
cvx_end

% --- Check convergence ---
if norm(theta2 - theta2_prev, 'fro') < tol
break;
end
% --- Save previous iteration variable values ---
theta2_prev = theta2;
z_mj2_prev = z_mj2;
T_2_prev = T_2;
eta2_prev = eta2;


% --- Update Dinkelbach parameters ---
for t = 1:T
for i = 1:I
for j = 1:J
lambda2(t,i,j) = (theta2(t,i) * S(t,i) * D(t,i)) / eta2(t,i,j);
end
end
end
end
% ========== Output Results ==========
fprintf('Traditional method total utility: %f\n', cvx_optval);


% Generate bar chart
e1 = 0;
e2 = 0;
for t = 1:T
    for i = 1:I
        if v(t,i) == 1
            e1 = e1 + eta1(t,i);
            e2 = e2 + eta2(t,i);
        end
    end
end

pic = bar(1,[e1;e2]);
legend("Propose","Traditional",'Location','northwest');
title('Time comparasion')
xlabel('High priority task')
ylabel('Task execution time(s)')
xtips1 = pic(1).XEndPoints;
ytips1 = pic(1).YEndPoints;
labels1 = string(pic(1).YData);
a=text(xtips1,ytips1,labels1,'HorizontalAlignment','center',...
'VerticalAlignment','bottom','FontSize',10)
xtips2 = pic(2).XEndPoints;
ytips2 = pic(2).YEndPoints;
labels2 = string(pic(2).YData);
c=text(xtips2,ytips2,labels2,'HorizontalAlignment','center',...
'VerticalAlignment','bottom','FontSize',10)
ax = gca;
ax.FontSize = 10;

% Add axis labels
fontname("Times New Roman");

I am doing researches about task offloading and resource allocation and I wrote this code to compare resources allocated to high-priority tasks using two different ways to calculate the utility of completing tasks on time(see variable “f1”).

From “f1”, you can see that compared with traditional method,

log(1 + tau_i(t,i) - T_2(t,i))/log(2)

I add the utility of completing a low-priority task on time(Gamma_L, a constant) to the utility of completing a high-priority task on time in the proposed method to make sure the benefit of completing a high-priority task is always higher than completing a low-priority task (which cannot be guaranteed in the traditional method).

log(1 + tau_i(t,i) - T_1(t,i))/log(2) + Gamma_L

In my opinion, this could make sure the system allocate more resources to the high-priority tasks than low-priority tasks compared with the traditional method. However, the experimental results did not align with my expection. Is there anything wrong with my code or my design?

2

I have not carefully studied your logic and implementation. However, given that you are applying an iterative method (and with some relaxations?), perhaps it can’t be assured that you are finding the global optimum of your 2 different optimization models. Even if both model solutions converge, they might not converge to the global optima; therefore the ordering of optimal objective values from your iterations might not the same as the ordering for the global optima. Perhaps you could try runs in which the solution from one of the formulations is fed in as the starting value for the other formulation, and see what happens.

3

I will try it, thank you.

4

Hello Mark, I followed your suggestion, but the result still differs from my expectation: the method allocating more resources to high-priority tasks seems to be random, regardless of which method’s solution I use as the initial input for the other method. Current comparison strategy doesn’t demonstrate the superiority of my proposed method. Could I ask for additional guidance?

Another question, I observe that my relaxation of z1/z2 (they acctually are binary, but I set them in [0,1]) causes beta, the big-M constant, influence the value of T1/T2 and u1/u2. Could I fix it without using Gurobi or other solver (I encountered issues with their installation)?

5

As I wrote before, you don’t know that you are finding the global optimum of either model, let alone both. Between your relaxations and iteration, do you really understand what any solution of your iterative process means? Also, do not use the quiet option, so that you can see the CVX and solver output, which might be helpful.

Gurobi and Mosek are the only solvers which can be called by CVX which can handle binary or integer variables (and for Mosek, only if there are no SDP constraints). I don’t know what you think you are accomplishing applying Big M to a relaxed variable, rather than to a binary variable; but you’re probably not accomplishing anything good.

6

Thank you for your guidance.
I downloaded Gurobi 12.0.2 and activated it with license code, I could run the example offered by Gurobi with matlab (but not with cvx) properly.

function mip1()
% Copyright 2025, Gurobi Optimization, LLC
% This example formulates and solves the following simple MIP model:
%  maximize
%        x +   y + 2 z
%  subject to
%        x + 2 y + 3 z <= 4
%        x +   y       >= 1
%        x, y, z binary

names = {'x'; 'y'; 'z'};

model.A = sparse([1 2 3; 1 1 0]);
model.obj = [1 1 2];
model.rhs = [4; 1];
model.sense = '<>';
model.vtype = 'B';
model.modelsense = 'max';
model.varnames = names;

gurobi_write(model, 'mip1.lp');

params.outputflag = 0;

result = gurobi(model, params);

disp(result);

for v=1:length(names)
    fprintf('%s %d\n', names{v}, result.x(v));
end

fprintf('Obj: %e\n', result.objval);
end

But when I ran cvx_grbgetkey, it showed :

---------------------------------------------------------------------------
CVX/Gurobi license key installer
---------------------------------------------------------------------------
Contacting the Gurobi Optimization license server...done.
The attempt to retrieve the license key failed with the following error:

    info  : grbgetkey version 9.0.0, build v9.0.0rc2
info  : Contacting Gurobi key server...
info  : Key for license ID 2672699 was successfully retrieved
info  : License expires at the end of the day on 2026-06-09
info  : Saving license key...
info  : License 2672699 written to file C:\Users\lenovo\AppData\Local\Temp\tp899b6f22_0eb7_4521_9b7b_2ec2a4c4f9d9\gurobi.lic

For information about this error, please consult the Gurobi documentation

    http://www.gurobi.com/documentation/5.5/quick-start-guide/node5

Once you have rectified the error, please try again.
---------------------------------------------------------------------------
警告: 转义字符 '\U' 无效。有关支持的特殊字符，请参阅 'doc sprintf'。 
> 位置：cvx_grbgetkey (第 338 行) 
info  : grbgetkey version 9.0.0, build v9.0.0rc2
info  : Contacting Gurobi key server...
info  : Key for license ID 2672699 was successfully retrieved
info  : License expires at the end of the day on 2026-06-09
info  : Saving license key...
info  : License 2672699 written to file C:---------------------------------------------------------------------------
Now that the license has been retrieved, please run CVX_SETUP to configure
CVX to use the Gurobi solver with the new license.
---------------------------------------------------------------------------

According to Problem with grbgetkey - #5 by krishnasandeep, I tried to run Gurobi outside the CVX and ran cvx_setup, but the Gurobi solver initialized is still 9.0.0 :

>> cvx_setup

---------------------------------------------------------------------------
CVX: Software for Disciplined Convex Programming       (c)2014 CVX Research
Version 2.2, Build 1148 (62bfcca)                  Tue Jan 28 00:51:35 2020
---------------------------------------------------------------------------
Installation info:
    Path: C:\Users\lenovo\Documents\MATLAB\cvx
    MATLAB version: 9.14 (R2023a)
    OS: Windows 10 amd64 version 10.0
    Java version: 1.8.0_202
Verfying CVX directory contents:
    WARNING: The following files/directories are missing:
        C:\Users\lenovo\Documents\MATLAB\cvx\sedumi\.travis.yml
    These omissions may prevent CVX from operating properly.
Preferences: 
    Path: C:\Users\lenovo\AppData\Roaming\MathWorks\MATLAB\cvx_prefs.mat
License host:
    Username: lenovo
    Host ID: 745d226b91db (eth0)
Installed license:
    No license installed.
No valid licenses found.
    Click here to fill out an academic license request
    for the username and first hostid listed above.
---------------------------------------------------------------------------
Setting CVX paths...already set!
Searching for solvers...5 shims found.
3 solvers initialized (* = default):
    Gurobi     9.00       {cvx}\gurobi\w64
 *  SDPT3      4.0        {cvx}\sdpt3
    SeDuMi     1.3.4      {cvx}\sedumi
4 solvers skipped:
    GLPK                  
        Could not find a GLPK installation.
    Gurobi_2   unknown    C:\gurobi1202\win64\matlab
        result = gurobi(model, params)
    Gurobi_3   unknown    C:\gurobi1202\win64\matlab
        result = gurobi(model, params)
    Mosek      9.1.9      {cvx}\mosek\w64
        Using MOSEK with CVX requires an academic MOSEK license
        or a MOSEK-enabled CVX Professional license.
1 solver issued warnings:
    Gurobi     9.00       {cvx}\gurobi\w64
        Using CVX with a commercial Gurobi license also requires
        a Gurobi-enabled CVX Professional license. Commercial
        users with Gurobi licenses should contact sales@cvxr.com.
Saving updated preferences...done.
Testing with a simple model...done!
---------------------------------------------------------------------------
To change the default solver, type "cvx_solver <solver_name>".
To save this change for future sessions, type "cvx_save_prefs".
Please consult the users' guide for more information.
---------------------------------------------------------------------------

and when I ran a test code:

clear all;
cvx_clear

cvx_solver Gurobi
cvx_begin
    variable z(3) binary
    minimize( sum(z) )
    subject to
        z(1) + z(2) >= 1
        z(3) - z(1) <= 0
cvx_end

disp(z);

an error occurred :

Calling Gurobi 9.00: 5 variables, 2 equality constraints

------------------------------------------------------------

------------------------------------------------------------

Status: Error

Optimal value (cvx_optval): NaN

错误使用 cvx_end

model.quadcon must be a struct array with fields q, and rhs

出错 test1 (第 11 行)

cvx_end

Now what could I do to fix it ?

7

Try using CVX 2.2.2 .

8

Thank you for your help.

10

I tried Mosek and runned a test code properly

clear all;
cvx_clear

cvx_solver Mosek
cvx_begin
    variable z binary;
    variable b;
    minimize( z + b )
    subject to
        z + log(1) > 0
        log(b) >= 1
cvx_end

disp(z);

Now I correct the big_M method auxiliary variable ‘z1’ to binary.

clear; clc;
% ========== System Parameters Initialization ==========
I = 3; % Number of UAVs
J = 2; % Number of RSUs
M = 2; % Number of vehicles
T = 5; % Time steps
D = randi([4e6,18e6], T, I); % Data volume
S = randi([4e2,1e3], T, I); % Task size
sigma2 = 1e-9; % Noise power
B_mj = 1e6 * ones(J, M); % Bandwidth
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 = 0.8; % Maximum transmit power
eta_max = 2e10;% Resource allocation upper limit
max_iter = 20; % Maximum iterations
tol = 1e-3; % Convergence tolerance

delta = 1e9* (1 + 0.2 * rand(I, 1)); % UAV local computing resources (CPU Hz)
P = rand(T,I,J);% Line-of-sight probability
r_ij = 0.06e6 + (2.34e6 - 0.06e6) * rand(T,I,J);% Transmission power between UAV i and RSU j
r_im = 0.06e6 + (2.34e6 - 0.06e6) * rand(T,I,J,M);% Transmission power between UAV i and vehicle m
g_mj = 0.02 * ones(J,M);% Channel gain between vehicle m and RSU j
tau_i = zeros(T,I);
for t = 1:T
    for i = 1: I
        tau_i(t,i) = D(t,i) * S(t,i)/delta(i); % Task completion time threshold
    end
end
v = randi([0,1], T,I); % Task priority (0: low, 1: high)

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


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

% ========== Initialize Variables ==========
theta1_prev = 0.5 * ones(T,I); % Initial task offloading ratio
for t = 1:T
    for i = 1:I
        for j =1:J
            for m =1 : M
                z_mj1_prev(t,i,j,m) = (log(1 + (p_max * g_mj(j,m) / sigma2)) ) / log(2); % Initial auxiliary variable value
            end
        end

    end
end
T_1_prev = tau_i;
eta1_prev = 0.5 * eta_max *ones(T,I,J);
for t = 1:T
    for i = 1:I
        for j = 1:J
            lambda1(t,i,j) = (theta1_prev(t,i) * S(t,i) * D(t,i)) / eta1_prev(t,i,j);
        end
    end
end
% ========== SCA Iteration ==========
for k = 1:max_iter

    % --- Build convex optimization problem for current iteration ---
    cvx_begin
    variable theta1(T,I) % Optimization variable: task offloading ratio
    variable z_mj1(T,I,J,M) % Auxiliary variable
    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 T_off1(T, I);
    % --- 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
                    % Calculate Taylor-expanded C term
                    C = (theta1_prev(t,i)*D(t,i))/(B_mj(j, m)*z_mj1_prev(t,i,j,m)) + ...
                        (theta1(t,i) * D(t,i))/(B_mj(j, m)*z_mj1_prev(t,i,j,m)) - ...
                        (theta1_prev(t,i)*D(t,i))/(B_mj(j, m)*(z_mj1_prev(t,i,j,m))^2) * z_mj1(t,i,j,m);

                    G = P(t,i,j) * ((theta1(t,i) * D(t,i))*inv_pos(r_ij(t,i,j))) + (1 - P(t,i,j)) * (((theta1(t,i) * D(t,i))*inv_pos(r_im(t,i,j,m))) + C);

                    T_off1(t, i) <= G + theta1(t,i)*S(t,i)*D(t,i) - lambda1(t,i,j)*eta1(t,i,j);% Minimum offloading time

                end
            end
            T_1(t,i) >= T_off1(t,i);

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

            % T_1(t,i)> tau_i
            f2 = -v(t,i)*Gamma_H + (1-v(t,i))*Gamma_L * exp(-c*(T_1_prev(t,i) - tau_i(t,i))) ...
                - c*(1-v(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
    0 <= eta1 <= eta_max;% Constraint C3
    % Constraint C4
    for t = 1:T
        for j = 1:J
            sum(eta1(t,:,j)) <= eta_max;
        end
    end
    0 <= p_mi1 <= p_max;% Constraint C5
    % Constraint C6
    for t = 1:T
        for m = 1:M
            sum(p_mi1(t,:,m)) <= p_max;
        end
    end

    % Auxiliary variable constraints
    for t = 1:T
        for i = 1:I
            for j =1:J
                for m =1 : M
                    z_mj1(t,i,j,m) <= (log(1 + (p_mi1(t,i,m) * g_mj(j,m) / sigma2)) ) / log(2);
                end
            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
    cvx_end

    % --- Check convergence ---
    if norm(theta1 - theta1_prev, 'fro') < tol
        break;
    end
    % --- Save previous iteration variable values ---
    theta1_prev = theta1;
    z_mj1_prev = z_mj1;
    T_1_prev = T_1;
    eta1_prev = eta1;


    % --- Update Dinkelbach parameters ---
    for t = 1:T
        for i = 1:I
            for j = 1:J
                lambda1(t,i,j) = (theta1(t,i) * S(t,i) * D(t,i)) / eta1(t,i,j);
            end
        end
    end
end

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

But it shows

警告: A non-empty cvx problem already exists in this scope.
   It is being overwritten. 
> 位置：cvxprob (第 28 行)
位置: cvx_begin (第 41 行)
位置: task_offloading0521 (第 61 行) 
 
Calling Mosek 11.0.22: 575 variables, 372 equality constraints
------------------------------------------------------------

MOSEK Version 11.0.22 (Build date: 2025-6-4 10:41:04)
Copyright (c) MOSEK ApS, Denmark WWW: mosek.com
Platform: Windows/64-X86

MOSEK warning 62 (MSK_RES_WRN_LARGE_AIJ): The A matrix contains a large value of 1.1e+10 in constraint '' (54) at variable '' (181).
MOSEK warning 62 (MSK_RES_WRN_LARGE_AIJ): The A matrix contains a large value of 1.1e+10 in constraint '' (55) at variable '' (181).
MOSEK warning 62 (MSK_RES_WRN_LARGE_AIJ): The A matrix contains a large value of 1.1e+10 in constraint '' (56) at variable '' (181).
MOSEK warning 62 (MSK_RES_WRN_LARGE_AIJ): The A matrix contains a large value of 1.1e+10 in constraint '' (57) at variable '' (181).
MOSEK warning 62 (MSK_RES_WRN_LARGE_AIJ): The A matrix contains a large value of 1.2e+10 in constraint '' (38) at variable '' (185).
MOSEK warning 62 (MSK_RES_WRN_LARGE_AIJ): The A matrix contains a large value of 1.2e+10 in constraint '' (39) at variable '' (185).
MOSEK warning 62 (MSK_RES_WRN_LARGE_AIJ): The A matrix contains a large value of 1.2e+10 in constraint '' (40) at variable '' (185).
MOSEK warning 62 (MSK_RES_WRN_LARGE_AIJ): The A matrix contains a large value of 1.2e+10 in constraint '' (41) at variable '' (185).
MOSEK warning 62 (MSK_RES_WRN_LARGE_AIJ): The A matrix contains a large value of 1.0e+10 in constraint '' (61) at variable '' (186).
MOSEK warning 62 (MSK_RES_WRN_LARGE_AIJ): The A matrix contains a large value of 1.0e+10 in constraint '' (62) at variable '' (186).
Warning number 62 is disabled.
MOSEK warning 54 (MSK_RES_WRN_LARGE_CON_FX): The equality constraint '' (128) is fixed to numerically large value -2.0e+10.
MOSEK warning 54 (MSK_RES_WRN_LARGE_CON_FX): The equality constraint '' (129) is fixed to numerically large value -2.0e+10.
MOSEK warning 54 (MSK_RES_WRN_LARGE_CON_FX): The equality constraint '' (130) is fixed to numerically large value -2.0e+10.
MOSEK warning 54 (MSK_RES_WRN_LARGE_CON_FX): The equality constraint '' (131) is fixed to numerically large value -2.0e+10.
MOSEK warning 54 (MSK_RES_WRN_LARGE_CON_FX): The equality constraint '' (132) is fixed to numerically large value -2.0e+10.
MOSEK warning 54 (MSK_RES_WRN_LARGE_CON_FX): The equality constraint '' (133) is fixed to numerically large value -2.0e+10.
MOSEK warning 54 (MSK_RES_WRN_LARGE_CON_FX): The equality constraint '' (134) is fixed to numerically large value -2.0e+10.
MOSEK warning 54 (MSK_RES_WRN_LARGE_CON_FX): The equality constraint '' (135) is fixed to numerically large value -2.0e+10.
MOSEK warning 54 (MSK_RES_WRN_LARGE_CON_FX): The equality constraint '' (136) is fixed to numerically large value -2.0e+10.
MOSEK warning 54 (MSK_RES_WRN_LARGE_CON_FX): The equality constraint '' (137) is fixed to numerically large value -2.0e+10.
Warning number 54 is disabled.
Problem
  Name                   :                 
  Objective sense        : minimize        
  Type                   : CONIC (conic optimization problem)
  Constraints            : 372             
  Affine conic cons.     : 0               
  Disjunctive cons.      : 0               
  Cones                  : 67              
  Scalar variables       : 575             
  Matrix variables       : 0               
  Integer variables      : 15              

Optimizer started.
Mixed integer optimizer started.
Threads used: 24
Presolve started.
Presolve terminated. Time = 0.00, probing time =  0.00
Presolved problem: 247 variables, 152 constraints, 322 non-zeros
Presolved problem: 0 general integer, 15 binary, 232 continuous
Presolved problem: 67 cones
Clique table size: 0
BRANCHES RELAXS   ACT_NDS  DEPTH    BEST_INT_OBJ         BEST_RELAX_OBJ       REL_GAP(%)  TIME  
0        0        1        0        -2.8866096215e+01    NA                   NA          0.0   
0        1        1        0        -2.8866097731e+01    -2.8866097731e+01    1.48e-13    0.0   
An optimal solution satisfying the relative gap tolerance of 1.00e-02(%) has been located.
The relative gap is 1.48e-13(%).

Objective of best integer solution : -2.886609773097e+01     
Best objective bound               : -2.886609773097e+01     
Initial feasible solution objective: Undefined
Construct solution objective       : Not employed
User objective cut value           : Not employed
Number of cuts generated           : 0
Number of branches                 : 0
Number of relaxations solved       : 1
Number of interior point iterations: 10
Number of simplex iterations       : 0
Time spend presolving the root     : 0.00
Time spend optimizing the root     : 0.00
Mixed integer optimizer terminated. Time: 0.01

Optimizer terminated. Time: 0.02    


Integer solution solution summary
  Problem status  : PRIMAL_FEASIBLE
  Solution status : INTEGER_OPTIMAL
  Primal.  obj: -2.8866097731e+01   nrm: 2e+10    Viol.  con: 9e-07    var: 5e-15    cones: 0e+00    itg: 9e-08  
Optimizer summary
  Optimizer                 -                        time: 0.02    
    Interior-point          - iterations : 10        time: 0.00    
    Simplex                 - iterations : 0         time: 0.00    
    Mixed integer           - relaxations: 1         time: 0.01    

------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +28.8661
 
错误使用  .* 
Disciplined convex programming error:
    Cannot perform the operation: {real affine} ./ {zero}

出错  ./  (第 19 行)
z = times( x, y, './' );

出错  *  (第 36 行)
    z = feval( oper, x, y );

出错  /  (第 15 行)
z = mtimes( x, y, 'rdivide' );

出错 task_offloading0521 (第 82 行)
                        (theta(t,i) * D(t,i))/(B_mj(j, m)*z_mj_prev(t,i,j,m)) - ...

I am trying to fix it. Could I ask for more suggestions? I’m still not familiar with using Mosek.

11

I suppose you take the solution from one iteration Mosek solved and use it in the next, but you use it so that division by zero appears. You should probably check that the data from the previous iteration satisfies whatever conditions you need for the next iteration, maybe inspect it and see what is not as you expected. The huge bound warnings don’t help, perhaps Mosek giving you a solution you didn’t expect for numerical reasons.

12

I will check it. Thank you.

14

Well, I checked the data from the first iteration and found that eta1 and p_mi1 are two zero matrices. So that, z_mj1 is also a zero matrix and lambda1 is a infinite matrix.(z_mj1 is an auxiliary variable, and its product with B_mj represents the transmission rate between vehicle m and RSU j.)

I wanna to know why I got this result? Are eta1 and p_mi1 gotten optimized in the first iteration?