issue with fmincon function.

I want to maximize rate using fmincon function. the problem is I have an error as "Arrays have incompatible sizes for this operation." and "Caused by: Failure in initial objective function evaluation. FMINCON cannot continue.". I provide a part of my code and the objective function here. I guess the problem is in the reshape part and it might be wrong. I would really appreciate it if someone could help me to write the objective function correctly.

when . I should mention that is a vector and is a matrix for .

here is my code.

%initial guess for P2 and W2
P2_0 = ones(Ns,1) * (tilde_P2 / Ns) / 2; % initial guess for P2
W2_0 = ones(Ms,Ns); %initial guess for W2

%combine initial guess in a single vector
x0 = [P2_0 ; W2_0(:)];

% Define the linear constraint (C2)
A = zeros(Np, Ns + Ms * Ns);
for l = 1:Np
    A(l, 1:Ns) = xi'; % Place xi' in the first Ns columns of each row
end

% Define the bounds(C1)
lb = [zeros(Ns, 1); -Inf(Ms * Ns, 1)]; % Lower bound (P2 >= 0, no lower bound for W2)
ub = [repmat(tilde_P2 / Ns, Ns, 1); Inf(Ms * Ns, 1)]; % Upper bound (P2 <= tilde_P2_max / Ns, no upper bound for W2)
%%%objective function
obj_fun = @(x) -((Nup - Nt) / Nup * sum(log(1 + x(1:Ns) ./ (sigma2 * reshape(x(Ns+1:end), Ms, Ns).^2))));

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

% Run the optimization
[x_opt, fval] = fmincon(obj_fun, x0, A, hat_zeta1, [], [], lb, ub, [], options);

% Extract optimized values for P2 and W2
P2_opt = x_opt(1:Ns);
W2_opt = reshape(x_opt(Ns+1:end), Ms, Ns);

[Nasrin]

@Torstensure. it is long so I write the objective function.

Here is the excutable code.

%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
dp = 1000;  %distance parameter in meters
ds = 500;   %distance parameter in meters
tilde_P2_dB = 20; %transmit power in dBW
tilde_P1_dB = 35; %example for tilde_P1 in dBW
p1_max_dB = tilde_P1_dB / Np; % example for p1_max in dBW
eta = 0.5;  %parameter eta
rho = 0.9;  %parameter rho
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 = 10^(tilde_P2_dB / 10);
tilde_P1 = 10 .^ (tilde_P1_dB / 10 );
sigma2 = 10^(sigma2_dBm / 10 -3 );  %dBm to mW then to W
a = 10 ^ (a_dB / 10);

%power parameters:
alpha = tilde_P2 / Ns; %parameter alpha
p2_max = tilde_P2 /Ns;
p1_max = tilde_P1 / Np; % example for p1_max


%%%%%%%%%%%%%%%%%%%%Channel modeling%%%%%%%%%%%%%%%%%%%%%%
% Geometry of the system
PBS_location = [0, 0]; %PBS location
SBS_location = [(-dp+ds)/2, (dp-ds)/2]; %SBS location

%generate distances
distance_pu_pbs = dp * rand(Np, 1);  %distances of PUs to PBS
distance_su_pbs = dp * rand(Ns, 1);  %distances of SUs to PBS
distance_su_sbs = ds * rand(Ns,1);  %distances of SUs to SBS
distance_pu_sbs = ds * rand(Ns, 1);  %distances of PUs to SBS

%calculate variance for the chanel model
variance_pu_pbs = a * (distance_pu_pbs/ d0).^(-path_loss_exp);
variance_su_pbs = a * (distance_su_pbs/ d0).^(-path_loss_exp);
variance_su_sbs = a * (distance_su_sbs/ d0).^(-path_loss_exp);
variance_pu_sbs = a * (distance_pu_sbs/ d0).^(-path_loss_exp);

%generate chanel matrices
H_1 = (randn(Mp, Np) +1i *randn(Mp, Np)) .*sqrt(variance_pu_pbs.') / sqrt(2); %matrix 32*4
H_12 = (randn(Mp, Ns) +1i *randn(Mp, Ns)) .*sqrt(variance_su_pbs.') / sqrt(2); %matrix 32*4
H_2 = (randn(Ms, Ns) +1i *randn(Ms, Ns)) .*sqrt(variance_su_sbs.') / sqrt(2); %matrix 20*4
H_21 = (randn(Ms, Np) +1i *randn(Ms, Np)) .*sqrt(variance_pu_sbs.') / sqrt(2); %matrix 20*4

%%%%%%%%%%%  mu_l %%%%%%%%%%%%%%%%%%%
ZFBF_mu_l = H_1' * H_1; % matrix 4*4
[ZFBF_mu_l_inv, ZFBF_mu_l_pinv] = compute_inverse_pseudoinverse(ZFBF_mu_l);
mu_l_matrix = H_1 * ZFBF_mu_l_inv; % matrix 32*4

mu_l = zeros(Np, 1); %matrix 4*1
for l = 1:Np
    column_l = mu_l_matrix(:, l);
    mu_l(l) = (norm(column_l))^(-1);
end

D1 = diag(mu_l); %matrix 4*4
W1 = mu_l_matrix * D1; % matrix 32*4
%%%%%%%%%%%%%%%%%%%%%% r1_l_max %%%%%%%%%%%%%%%%%
r1_max = zeros(Np,1); %matrix Np*1
for l = 1:Np
    r1_l_max = log(1 + ((mu_l(l)^2 * p1_max) / sigma2));
    r1_max(l) = r1_l_max;
end
r1 = rho * r1_max;

%%%%%%%%%%%%%%%% waterfilling for tilde_p1_l %%%%%%%%%%%%%%%
%computing phi
phi_l = mu_l .^2 / sigma2;

%%computing psi
% Objective function for psi
obj_fun = @(psi) sum(max(0, psi - 1 ./ phi_l)) - tilde_P1;

% Initial guess for psi
psi_init = tilde_P1 / Np;

% Options
options = optimoptions('fsolve', 'Display', 'off');

% Solve psi
psi = fsolve(obj_fun, psi_init, options);

%tilde_p1_l
tilde_p1_l = max(0, psi - 1./phi_l);

%%%%%%%%%%%%%%%%%%%%% xi  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xi = zeros(Ns,1); % matrix Ns*1 
for k = 1:Ns
    xi_k = 0;
    for l = 1:Np
        w1_l = W1(:, l);
        h12_k = H_12(:, k);
    
        % Compute the inner product and absolute value squared
        abs_val_squared = abs(w1_l' * h12_k)^2;
        xi_k = xi_k + tilde_p1_l(l) * abs_val_squared;
    end
    xi(k) = xi_k ./ sum(tilde_p1_l);% matrix 4*1
end
     
%%%%%%%%%%%%%% hat_pi %%%%%%%%%%%%%%%%%%%
hat_pi= zeros(Np,1); %matrix Np*1
for l = 1:Np
    hat_pi_l = Nt * log(1 + (mu_l(l)^2 * p1_max) / (alpha * sum(xi) + sigma2));
    hat_pi(l) = hat_pi_l;
end

%%%%%%%%%%%%%% hat_zeta1  %%%%%%%%%%%%%%%%%%%
hat_zeta1 = zeros(Np,1);% matrix Np*1
for l = 1:Np
    num = mu_l(l)^2 * p1_max;
    denom = 2 ^((Nup * r1(l) - hat_pi(l))/(Nup - Nt));
    hat_zeta1_l = (num / (denom -1)) - sigma2;
    hat_zeta1(l) = hat_zeta1_l;
end

%%%%%%%%%%%%%%% Optimization problem %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%initial guess for P2 and W2
P2_0 = ones(Ns,1) * (tilde_P2 / Ns) / 2; % initial guess for P2
W2_0 = ones(Ms,Ns); %initial guess for W2

%combine initial guess in a single vector
x0 = [P2_0 ; W2_0(:)];

% Define the linear constraint (C2)
A = zeros(Np, Ns + Ms * Ns);
for l = 1:Np
    A(l, 1:Ns) = xi'; % Place xi' in the first Ns columns of each row
end

% Define the bounds(C1)
lb = [zeros(Ns, 1); -Inf(Ms * Ns, 1)]; % Lower bound (P2 >= 0, no lower bound for W2)
ub = [repmat(tilde_P2 / Ns, Ns, 1); Inf(Ms * Ns, 1)]; % Upper bound (P2 <= tilde_P2_max / Ns, no upper bound for W2)
%%%objective function
obj_fun = @(x) -((Nup - Nt) / Nup * sum(log(1 + x(1:Ns) ./ (sigma2 * reshape(x(Ns+1:end), Ms, Ns).^2))));

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

% Run the optimization
[x_opt, fval] = fmincon(obj_fun, x0, A, hat_zeta1, [], [], lb, ub, [], options);

Arrays have incompatible sizes for this operation.

Error in solution>@(x)-((Nup-Nt)/Nup*sum(log(1+x(1:Ns)./(sigma2*reshape(x(Ns+1:end),Ms,Ns).^2)))) (line 148)
obj_fun = @(x) -((Nup - Nt) / Nup * sum(log(1 + x(1:Ns) ./ (sigma2 * reshape(x(Ns+1:end), Ms, Ns).^2))));

Error in objfunEvaluator (line 5)
fval = feval(Objfun, x, self.FunArgs.AdditionalParameters{:});

Error in OptimFunctions/objective (line 271)
[fval_, fgrad_, hess_] = self.ObjectiveFunAndGrad(self,self.FunFcn{3},...

Error in OptimFunctions/objective_first_eval (line 614)
[fval,self] = self.objective(X0);

Error in fmincon (line 500)
[initVals.f,initVals.g,HESSIAN,funObj] = funObj.objective_first_eval(X);

Caused by:
Failure in initial objective function evaluation. FMINCON cannot continue.

% Extract optimized values for P2 and W2
P2_opt = x_opt(1:Ns);
W2_opt = reshape(x_opt(Ns+1:end), Ms, Ns);

% Display the result
disp('Optimized values of P2:');
disp(P2_opt);
disp('Optimized values of W2:');
disp(W2_opt);
disp('Optimized objective value:');
disp(-fval);


% Compute the inverse and pseudoinverse of a matrix
function [matrix_inv, matrix_pinv] = compute_inverse_pseudoinverse(matrix)

    % Initialize the inverse as empty, to be used if matrix is not invertible
    matrix_inv = [];
   
    % Check if the matrix is square
    [rows, cols] = size(matrix);
    if rows == cols
        % Check if the matrix is invertible by determining its determinant
        if det(matrix) ~= 0
            % Compute the inverse of the matrix
            matrix_inv = inv(matrix);
        else
            warning('Matrix is singular and cannot be inverted. Only pseudoinverse will be returned.');
        end
    else
        warning('Matrix is not square. Only pseudoinverse will be returned.');
    end

    % Compute the pseudoinverse of the matrix
    matrix_pinv = pinv(matrix);
end

Answer 1

[Torsten]

I guess the objective function should read

obj_fun = @(x) -(Nup - Nt) / Nup * sum(log2(1 + x(1:Ns) ./ (sigma2 * (vecnorm(reshape(x(Ns+1:end), Ms, Ns)).^2).')));

If you still encounter errors, please give us executable code. Above, many variables are undefined.

[Nasrin]

@Torsten It works now. Thank you. you helped me alot.