fmincon is not converge

Hi every one. I have a minimization problem with linear inequality and nonlinear inequality. when I run the code. it does not converge. you can consider H1_n_all, H12_n_all, H2_n_all and H21_n all as any random 3D matrix.

I was wondering what is the problem with my code.

thank you in advance for your help.

%% Define parameters
Mp = 32;  % number of antennas at PBS 
Ms = 20;  % number of antennas at SBS
Np = 4;   % number of primary users
Ns = 4;   % number of secondary users
Nsc = 4;  % number of subcarriers
dp = 1000;  % distance parameter in meters
ds = 500;   % distance parameter in meters
tilde_P2_tot_dB = 20; % transmit SN power in dBW
% Define tilde_P1_tot_dB as a row vector
tilde_P1_tot_dB_vec = [35, 40, 45, 50, 55]; % Example row vector of values for tilde_P1_tot_dB

eta = 0.5;  % parameter eta for tightened constraints
rho = 0.9;  % parameter rho for PN rate adjustment
Nt = 12;    % number of training symbols
Nup = 200;  % number of uplink symbols
Ndp = 200;  % number of downlink symbols
sigma2_dBm = -104;  % noise power in dBm
a_dB = -137;   % path loss at reference distance in dB
path_loss_exp = 3.5;  % path loss exponent
d0 = 1000;  % reference distance in meter

% Convert dBW and dBm to linear scale
tilde_P2_tot = 10^(tilde_P2_tot_dB / 10);
sigma2 = 10^(sigma2_dBm / 10 - 3);  % dBm to mW then to W
a = 10 ^ (a_dB / 10);

% Power parameters
alpha = tilde_P2_tot / (Ns * Nsc); % training power of each SUs' sub-carrier uniformly allocated
%% Preallocate cell arrays to store parameters
H1_n_all = cell(Nsc, 1);   
H12_n_all = cell(Nsc, 1);  
H2_n_all = cell(Nsc, 1);   
H21_n_all = cell(Nsc, 1);  
mu_l_n_all = zeros(Nsc, Np); 
phi_l_n_all = zeros(Nsc, Np); 
D1_n_all = cell(Nsc, 1); 
W1_n_all = cell(Nsc, 1); 
W2_n_all = cell(Nsc, 1); 
p1_l_n_all = zeros(Nsc, Np);
r1_l_max_all = zeros(1, Np);
pi_l_all = zeros(1, Np);  
r1_l_all = zeros(1, Np);

% Preallocate arrays to store results
num_tests = length(tilde_P1_tot_dB_vec); % Number of different tilde_P1_tot_dB values
P2_opt_all = cell(num_tests, 1);  % To store optimized P2 for each value
W2_opt_all = cell(num_tests, 1);  % To store optimized W2 for each value
R2_opt = zeros(num_tests, 1); %to store optimum sum_rate for Su

%% Load channel matrices
for n = 1:Nsc
    file_name = sprintf('fixed_channel_n_%d.mat', n);
    [H1_n_all{n}, H12_n_all{n}, H2_n_all{n}, H21_n_all{n}] = load_channel(file_name);
end

% Convert cell arrays to 3D matrices
H1_n_matrix = cat(3, H1_n_all{:});
H12_n_matrix = cat(3, H12_n_all{:});
H2_n_matrix = cat(3, H2_n_all{:});
H21_n_matrix = cat(3, H21_n_all{:});

% Loop over each value of tilde_P1_tot_dB
for idx = 1:num_tests
    % Convert current tilde_P1_tot_dB to linear scale
    tilde_P1_tot = 10^(tilde_P1_tot_dB_vec(idx) / 10);
    
    % Compute mu, D1, W1 for each subcarrier
    for n = 1:Nsc
        H1_n = H1_n_matrix(:, :, n);
        product_matrix = H1_n * inv(H1_n' * H1_n);
        mu_l_n = 1 ./ vecnorm(product_matrix, 2, 1);
        mu_l_n_all(n, :) = mu_l_n; 
        D1_n = diag(mu_l_n); 
        D1_n_all{n} = D1_n;  
        W1_n = product_matrix * D1_n;  
        W1_n_all{n} = W1_n;
    end

    % Convert W1_n_all to 3D matrix
    W1_n_matrix = cat(3, W1_n_all{:});

    % Compute phi_l_n
    phi_l_n_all = mu_l_n_all.^2 / sigma2;

    % Water-filling over all subcarriers
    for n = 1:Nsc
        [p1_l_n_all(n, :), lambda_opt, r1_l_max_all(n)] = waterfilling_bisection_UL(phi_l_n_all(n, :), tilde_P1_tot / Np, 1e-6);
        r1_l_all(n) = rho * r1_l_max_all(n);
    end

    % Compute pi_l for each primary user
    for l = 1:Np
        sum_term = 0;
        for n = 1:Nsc
            W1_n = W1_n_all{n};
            H12_n = H12_n_matrix(:, :, n);  
            norm_term = norm(W1_n(:, l)' * H12_n, 'fro')^2;
            log_term = log2(1 + ((mu_l_n_all(n, l)^2) * p1_l_n_all(n, l)) / (alpha * norm_term + sigma2));
            sum_term = sum_term + log_term;
        end
        pi_l_all(l) = Nt * sum_term;
    end

    % Compute zeta1_l
    zeta1_l = (Nup - Nt) ./ ((r1_l_all) * Nup - pi_l_all);

    % Compute W2_n using ZF beamforming
    H_n_matrix = [H21_n_matrix, H2_n_matrix];
    J = [zeros(Np,Ns); eye(Ns)];
    for n = 1:Nsc
        H_n = H_n_matrix(:, :, n);
        W2_n = H_n * inv(H_n' * H_n) * J;
        W2_n_all{n} = W2_n; 
    end
    W2_n_matrix = cat(3, W2_n_all{:});

    %% Optimization using fmincon
    p2_k_max = tilde_P2_tot / Ns;

    % Initial guess for the beamforming vectors and power allocations
    %w2_init = randn(Ms, Ns, Nsc);  
    % W2_init = ones(Ms, Ns, Nsc);
    p2_init = p2_k_max/2 * ones(Nsc, Ns);  

    x0 = p2_init(:);
    
    % Define bounds for the optimization variables
    lb = zeros(size(x0));  % Lower bound: power and beamforming must be non-negative
    % ub = [repmat(p2_max, Nsc * Ns, 1); Inf * ones(numel(w2_init), 1)];  % Upper bound for powers, no bound for beamforming
    ub = [];

    % Linear constraint C1
    A = zeros(Ns, numel(x0));  
    b = p2_k_max * ones(Ns, 1);  
    Aeq = [];
    Beq = [];

    for k = 1:Ns
        for n = 1:Nsc
            A(k, (n-1)*Ns + k) = 1;
        end
    end

    % Objective function
    objective = @(x) -computeObjective(x, Nup, Nt, Nsc, Ns, sigma2, W2_n_matrix);

    % Nonlinear constraints
    nonlcon = @(x) computeNonlinearConstraint(x, H12_n_matrix, W1_n_matrix, mu_l_n_all, sigma2, zeta1_l, p1_l_n_all, W2_n_matrix);

    % Optimization options
    options = optimoptions('fmincon', 'Display', 'iter', 'Algorithm', 'interior-point');

    % Call fmincon
    [x_opt, fval, ~] = fmincon(objective, x0, A, b, Aeq, Beq, lb, ub, nonlcon, options);

    % Extract results
    P2_opt = reshape(x_opt, [Nsc, Ns]);

    % Store results
    P2_opt_all{idx} = P2_opt;
    R2_opt(idx) = fval;
    end

% % Display results
% disp('Optimal power allocation P2 for each value of tilde_P1_tot_dB:');
% disp(P2_opt_all);

disp("Optimal SNs' achievable rate for each value of tilde_P1_tot_dB:");
disp(R2_opt);


%% plot
% Plot R2_opt for different values of tilde_P1_tot_dB
figure;

% Plot for the first curve (e.g., perfect CSI)
plot(tilde_P1_tot_dB_vec, R2_opt, 'b-o', 'LineWidth', 2, 'MarkerSize', 8);
hold on;

xlabel('Total power of PN in downlink, $\tilde{P}_{1}$ (dBm)', 'Interpreter', 'latex', 'FontSize', 12);
ylabel('SNs UL achievable rate', 'Interpreter', 'latex', 'FontSize', 12);

% Customize axes
set(gca, 'FontSize', 12);  % Adjust axis font size
grid on;  % Enable grid

% Set y-axis limits for better readability (optional)
ylim([min(R2_opt)-0.5, max(R2_opt)+0.5]);

%% fmincon functions
%%% Objective Function
function obj = computeObjective(x, Nup, Nt, Nsc, Ns, sigma2, W2_n_matrix)
    % Extract P2 and W2 from the optimization variable
    p2 = reshape(x, [Nsc, Ns]);
    
    % Compute the objective: sum over subcarriers and secondary users
    sum_term = 0;
    for n = 1:Nsc
        for k = 1:Ns
            sum_term = sum_term + log2(1 + p2(n,k) / (sigma2 * vecnorm(W2_n_matrix(:,k,n), 2).^2));
        end
    end
    obj = -(Nup - Nt) / Nup * sum_term; 
end

%%% Nonlinear Constraints
function [c, ceq] = computeNonlinearConstraint(x, H12_n_matrix, W1_n_matrix, mu_l_n_all, sigma2, zeta1_l, p1_l_n_all, W2_n_matrix)
    [Ms, Ns, Nsc] = size(H12_n_matrix);  % Ms: # antennas at SBS, Ns: # secondary users, Nsc: # subcarriers
    [Mp, Np] = size(W1_n_matrix(:, :, 1)); % Mp: # antennas at PBS, Np: # primary users
    
    %Only optimize over p2, W2 is fixed
    p2 = reshape(x, [Nsc, Ns]);
    
    c = zeros(Np, 1);  % Inequality constraints (C2)
    for l = 1:Np
        sum_term = 0;
        for n = 1:Nsc
            H12_k_n = H12_n_matrix(:,:,n); %matrix for different n
            W1_l_n = W1_n_matrix(:,:,n); %matrix for different n
            norm_term = zeros(1,Ns); % to store norm values for each secondary user
            for k = 1:Ns
                norm_term(k) = norm(W1_l_n(:,l)'.* H12_k_n(:, k), 'fro') .^2;
            end
            interference = sum(p2(n,:) .* norm_term);
            % p1_l_n = 1;  % Power for primary user, assuming it's already defined
            log_term = log2(1 + (mu_l_n_all(n,l)^2 * p1_l_n_all(n,l)) / (interference + sigma2));
            sum_term = sum_term + log_term;
        end
        c(l) = 1/zeta1_l(l) - sum_term;  % C2 inequality
    end
    ceq = [];  % No equality constraints
end

%% functions
%function for loading the channel matrices%%%
function [H_1_n, H_12_n, H_2_n, H_21_n] = load_channel(file_name)
    % Load the channel matrices from the specified .mat file
    data = load(file_name);
    H_1_n = data.H_1;
    H_12_n = data.H_12;
    H_2_n = data.H_2;
    H_21_n = data.H_21;
    disp(['Channel matrices loaded from ', file_name]);
end

Answer 1

[John D'Errico]

It is very rarely easy to know why a solver fails to converge on a messy thing like this. But you should never let it come to that point.

BEFORE you ever use a solver on a problem like this, first test your objective function. Evaluate it at several points in the domain of your parameter space, points where you expect it should be well-behaved. Make several observations at this stage.

Does the objective function produce a REAL result, and is it SCALAR valued?

If you perturb the parameters by a small amount, does the objective change?

Are there any issues you observe that might be cause for a problem? That is, if the function is essentially constant, then you should not expect to see fmincon fail. If the function returns complex results, then again, fail. If the function returns a vector or array, then again, fail.

Look carefully at the variation you see in the objective function. Remember that a tolerance will be applied, and if the function is seen to be changing not at all, or less than the tolerance applied, then what do you expect to happen?

Think about your code. Is there any point where you form a discretized operation inside? That, do you ever use floor, fix, round, ceil, etc.? If so, then your function will not e differentiable, or even continuous. Again, FAIL.

Is there any random aspect of the code inide the objective? Fmincon assumes the function is absolutely deterministic. If not, then fail.

The above are some basic tests you should ALWAYS perform. (Given some time, I could probably think of a few more, but the above tests would obviate most of the desperate questions we see here, about why some optimzer has failed to produce a viable result.)

NEVER just throw together a pile of code, and then throw it at an optimizer. Instead, think about what it does, and if it makes sense to employ fmincon there. Test the objective, and verify it does what it should do. (Do I do all of the above tests when I use a solver? Well, usually not. But I've been writing code and using optimizers since probably before you were born. Until you get to that point, then follow my advice.)

When you do see a failure to converge, first, go back and perform EXACTLY those same tests I just told you to run.

Next, look at the return flag. What does the exitflag tell you? If you don't understand what the various returns are for exitflag, then tell us what you got, and ask for help on interpreting that result.

Finally, if everything above seems reasonable, look at the final point fmincon returns. What is the objective value at that location? Is it just that your function has no minimum? For example, if I told fmincon to minimize exp(x), then I would expect a divergent result at x==-inf, or at least as close as it could get to -inf before it gives up.