@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
