How to extracting parameters in iterations in genetic algorithms?

Heiio everyone. I am trying extracting some parameters which is calculated by the input parameters in iterations in genetic algorithms. I trying to display them in objective function (It's name is objectiveFunction) or in another ouput function (It's name is myOutputFunctionWithPower) during every iteration. The parameters I want to extract are AF_variance_per_gen, AF_mean_per_gen, and penalty_per_gen ,which I set them as global parameters. But I find that the progamm can't show them, so as global parameter: gen_count. I don't know why. Please help me, Thanks a lot!

The part of result I try to print in live script is as following. 'Test Generation is printed in "objectiveFunction" and Standard,Mean, and Penalty are printed in "myOutputFunctionWithPower"

The whole code is shown bellow:

clear; clc; close all;
global AF_variance_per_gen penalty_per_gen AF_mean_per_gen gen_count;
AF_variance_per_gen = [];
penalty_per_gen = [];
AF_mean_per_gen = [];
gen_count = 1;
% Basic parameters
lambda = 0.0107142857142857;  % Wavelength (28GHz)
Frequency = 28e9;             % Frequency
Lightspeed = 3e8;             % Speed of light
k = 2 * pi / lambda;          % Wavenumber

% Antenna array configuration
M = 10; % MxN antenna array
N = 10;

% Target angle ranges (in radians)
theta_target_min = deg2rad(1);       % Minimum elevation angle
theta_target_max = deg2rad(80);      % Maximum elevation angle
phi_target_min = deg2rad(-180);      % Minimum azimuth angle
phi_target_max = deg2rad(180);       % Maximum azimuth angle

% Number of antennas
num_antennas = M * N;

% Define the 361 possible phase values (0 to 360 degrees)
possiblePhases_deg = 0:1:188;  % 361 possible phase values (in degrees)
possiblePhases_rad = deg2rad(possiblePhases_deg);  % Convert to radians

% ------------------ Antenna Layout Range ------------------
% Define the range of the area
x_range = [-15 * lambda, 15 * lambda];   % x-coordinate range
y_range = [-15 * lambda, 15 * lambda];   % y-coordinate range
% ----------------------------------------------------------

% Generate initial antenna positions that satisfy the distance constraints
[initial_x, initial_y] = generateInitialPositions(num_antennas, lambda, x_range, y_range);

% Initial random phase indices between 1 and 188.62 
initialPhaseIndices = randi(length(possiblePhases_deg), 1, num_antennas);
%這樣寫會讓項位破千，爆掉，結果會很差。
% initialPhaseIndices = randi([0 188],1,100); 

% Optimization variables: phase_indices (num_antennas), x_positions (num_antennas), y_positions (num_antennas)
initialParams = [initialPhaseIndices, initial_x, initial_y];
zerophase = zeros(1,100);
% Define lower and upper bounds for phase indices and positions
lb = [ones(1, num_antennas), x_range(1) * ones(1, num_antennas), y_range(1) * ones(1, num_antennas)];
ub = [length(possiblePhases_deg) * ones(1, num_antennas), x_range(2) * ones(1, num_antennas), y_range(2) * ones(1, num_antennas)];

% Define integer constraint variable indices (phase indices are integers)
intcon = 1:num_antennas;

% Define GA options
ga_options = optimoptions('ga', 'Display', 'iter', ...
                          'MaxGenerations', 100, 'PopulationSize', 800, ...
                          'FunctionTolerance', 1e-7, 'ConstraintTolerance', 1e-7, ...
                          'UseParallel', true, 'MaxStallGenerations', 200, ...
                          'InitialPopulationMatrix', generateInitialPopulation(num_antennas, x_range, y_range, lambda, length(possiblePhases_deg)),...
                          'CreationFcn','gacreationuniformint',...
                          'SelectionFcn','selectiontournament',...
                          'MutationFcn','mutationgaussian',...
                          'CrossoverFcn','crossoverscattered',...
                          'OutputFcn', @myOutputFunctionWithPower,...
                          'PlotFcn',{'gaplotscorediversity','gaplotscores'});

% Run the genetic algorithm optimization
[optimalParams_ga, ~] = ga(@(params) objectiveFunction(params, num_antennas, k, ...
                                                      theta_target_min, theta_target_max, ...
                                                      phi_target_min, phi_target_max, lambda, possiblePhases_rad), ...
                           length(initialParams), [], [], [], [], lb, ub, ...
                           @(params) distanceConstraint(params, num_antennas, lambda, x_range, y_range), ...
                           intcon, ga_options);

% Extract optimal phase indices and x, y positions
optimalPhaseIndices = round(optimalParams_ga(1:num_antennas));
optimalPhases_ga = possiblePhases_rad(optimalPhaseIndices);  % Map phase indices to phase values
optimal_x = optimalParams_ga(num_antennas+1:2*num_antennas);
optimal_y = optimalParams_ga(2*num_antennas+1:end);

% Generate complex phase weights
w_ga = exp(1i * optimalPhases_ga);  % Optimal antenna phases from GA
w_original = exp(1i * deg2rad(initialPhaseIndices));
% Display each antenna's phase angle and optimal x, y positions
disp('Phase angles and optimal x, y positions for each antenna (in degrees):');
for idx = 1:num_antennas
    phaseAngle = mod(rad2deg(optimalPhases_ga(idx)), 360);    % Convert to degrees
    fprintf('Antenna %d: Phase Angle: %.2f degrees, x: %.5f meters, y: %.5f meters\n', idx, phaseAngle, optimal_x(idx), optimal_y(idx));
end

% Save antenna phases and positions to a CSV file
combined_csv_filename = 'C:\Users\Kiven Yllow\Documents\新研究Reconfigurable Intelligent Surface\matlab最佳化\兆富的程式\GA_10by10_File\ChangePhasePos_188degree_10by10_GA.csv';

% Convert phases to degrees
phase_vector_deg = mod(rad2deg(optimalPhases_ga), 360);  % Convert phases to degrees
phase_vector_deg_original = mod(initialPhaseIndices, 360);

% Combine into a matrix: columns are phase (degrees), x position, y position
combined_matrix = [phase_vector_deg', optimal_x', optimal_y'];

% Save to CSV file
writematrix(combined_matrix, combined_csv_filename);

fprintf('Antenna phases and positions have been saved to:\n%s\n', combined_csv_filename);

% ------------------ Plot Antenna Layout with Phases ------------------
figure;
%將帶有初始的相位跟位置的點波源畫出來
scatter(initial_x, initial_y, 50, zerophase, 'filled');
xlabel('X Position (meters)');
ylabel('Y Position (meters)');
title('Antenna 2D Layout with Excitation Phases (Original)');
colorbar;
colormap('jet');
grid on;
axis equal;
hold on;
% Annotate each antenna with its phase value (degrees)
for idx = 1:num_antennas
    phase_deg = zerophase(idx);
        text(initial_x(idx), initial_y(idx), sprintf('%.0f°', phase_deg), ...
        'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', 'FontSize', 6);
end
hold off;


figure;
%將帶有初始的相位跟位置的點波源畫出來(此初始值沒有介入計算)
scatter(initial_x, initial_y, 50, phase_vector_deg_original, 'filled');
xlabel('X Position (meters)');
ylabel('Y Position (meters)');
title('Antenna 2D Layout with Excitation Phases (Original)');
colorbar;
colormap('jet');
grid on;
axis equal;
hold on;
% Annotate each antenna with its phase value (degrees)
for idx = 1:num_antennas
    phase_deg = phase_vector_deg_original(idx); % Get phase value (degrees)
    text(initial_x(idx), initial_y(idx), sprintf('%.0f°', phase_deg), ...
        'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', 'FontSize', 6);
end
hold off;

figure;
scatter(optimal_x, optimal_y, 50, phase_vector_deg, 'filled');
xlabel('X Position (meters)');
ylabel('Y Position (meters)');
title('Antenna 2D Layout with Excitation Phases (GA)');
colorbar;
colormap('jet');
grid on;
axis equal;
hold on;

% Annotate each antenna with its phase value (degrees)
for idx = 1:num_antennas
    phase_deg = phase_vector_deg(idx); % Get phase value (degrees)
    text(optimal_x(idx), optimal_y(idx), sprintf('%.0f°', phase_deg), ...
        'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', 'FontSize', 7);
end
hold off;

% Compute distance matrix between optimized antenna positions
positions_opt = [optimal_x', optimal_y'];
distance_matrix_opt = squareform(pdist(positions_opt));

% Check nearest neighbor distances for each antenna
for i = 1:num_antennas
    distances_i = distance_matrix_opt(i, :);
    distances_i(i) = inf; % Exclude self-distance
    min_distance = min(distances_i);
    if min_distance <= (0.004) || min_distance >= (2 * lambda)
        fprintf('Antenna %d has a nearest neighbor distance of %.6f meters, which does not satisfy the constraint.\n', i, min_distance);
    end
end

%-------------新增---------畫出初始解和最佳化後解的power比較圖-------------
    % 计算对应的 theta 和 phi（整个球面）
    u = 1:-0.01:-1;
    v = -1:0.01:1;
    theta_ob = acos(u).';    % θ ∈ [0, π]  變成行向量，同一行不同列的數字代表固定水平角的鉛直角度數
    phi_ob = pi * v ;                              %列向量，同一列不同行的數字代表固定鉛質角的水平角度數
    % 計算初始相位與位置的power，還有全0相位與初始位置的power
    AF_values = zeros(length(u),length(v)); %測試先註解
    original_AF_values = zeros(length(u),length(v));
    zerophase_AF_values = zeros(length(u),length(v));
    for idx = 1:num_antennas
        %psi是用所有觀察角代表的相位二維陣列
        psi = k * ( optimal_x(idx) * sin(theta_ob) * cos(phi_ob) + optimal_y(idx) * sin(theta_ob) * sin(phi_ob) );
        original_psi = k * ( initial_x(idx) * sin(theta_ob) * cos(phi_ob) + initial_y(idx) * sin(theta_ob) * sin(phi_ob) );
        %w_ga(idx)是第idx個點波源的自帶相位(落後第1個點波源某某度)，AF_values是跟psi一樣是二維陣列
        AF_values = AF_values + w_ga(idx) * exp(1i * psi); %測試先註解
        original_AF_values = original_AF_values + w_original(idx) * exp(1i * original_psi);
        zerophase_AF_values = zerophase_AF_values + 1 * exp(1i * original_psi);
    end
    AF_values = abs(AF_values);            original_AF_values = abs(original_AF_values);           zerophase_AF_values = abs(zerophase_AF_values);
    AF_square = AF_values.^2;              original_AF_square = original_AF_values.^2;             zero_AF_square = zerophase_AF_values.^2;
    AF_square_dB = 10*log10(AF_square);    original_AF_square_dB = 10*log10(original_AF_square);   zero_AF_square_dB = 10*log10(zero_AF_square);
    POWER_values = AF_values.^2 / 377;     original_POWER_values = original_AF_values.^2 / 377;    zero_POWER_values = zerophase_AF_values.^2 / 377;
    POWER_dB = 10*log10(POWER_values);     original_POWER_dB = 10*log10(original_POWER_values);    zero_POWER_dB = 10*log10(zero_POWER_values);
figure;
hold on
plot(theta_ob.'/pi*180,POWER_dB(:,101),'color',"#77AC30",'LineWidth',2.0) 
plot(theta_ob.'/pi*180,original_POWER_dB(:,101),'b','LineWidth',2.0)
plot(theta_ob.'/pi*180,zero_POWER_dB(:,101),'m','LineWidth',2.0)
hold off
title('比較(沒有正歸化 X方向)')
legend('調整過的','原本的')
grid on
xlim([0 180])

xlabel('theta,°')
ylabel('power_d_B')


figure;
hold on
plot(theta_ob.'/pi*180,POWER_dB(:,151),'color',"#77AC30",'LineWidth',2.0) 
plot(theta_ob.'/pi*180,original_POWER_dB(:,151),'b','LineWidth',2.0)
plot(theta_ob.'/pi*180,zero_POWER_dB(:,151),'m','LineWidth',2.0)
hold off
title('比較(沒有正歸化 Y方向)')
legend('調整過的','原本的')
grid on
xlim([0 180])

xlabel('theta,°')
ylabel('power_d_B')

figure;
patternCustom(zero_POWER_dB.',theta_ob.'/pi*180,phi_ob/pi*180, 'CoordinateSystem', 'polar');
    title('3D Array Factor Pattern, zero phase and random position, same with origin');
    caxis([-60,15]);

figure;
patternCustom(original_POWER_dB.',theta_ob.'/pi*180,phi_ob/pi*180, 'CoordinateSystem', 'polar');
    title('3D Array Factor Pattern, original');
    caxis([-60,15]);
figure;
patternCustom(POWER_dB.', theta_ob.'/pi*180, phi_ob/pi*180, 'CoordinateSystem', 'polar');
    title('3D Array Factor Pattern, optimization');
    caxis([-60,15]);
%--------------繪製最佳化前後的pdf------------
PAT_optimize_all = POWER_values;  PAT_origin_all = original_POWER_values; PAT_zero_all = zero_POWER_values;
PATAFSQUARE_optimize_all = AF_square;   PATAFSQUARE_origin_all = original_AF_square;  PATAFSQUARE_zero_all = zero_AF_square;
for i = 2:201  %把南極點北極點重複的數值去掉
    PAT_zero_all(1,i) = -10; PAT_origin_all(1,i) = -10;   PAT_optimize_all(1,i) = -10;
    PAT_zero_all(201,i) = -10;   PAT_origin_all(201,i) = -10;   PAT_optimize_all(201,i) = -10;
    PATAFSQUARE_zero_all(1,i) = -10; PATAFSQUARE_origin_all(1,i) = -10;   PATAFSQUARE_optimize_all(1,i) = -10;
    PATAFSQUARE_zero_all(201,i) = -10;   PATAFSQUARE_origin_all(201,i) = -10;   PATAFSQUARE_optimize_all(201,i) = -10;
end
point_all_optimize = zeros(1,2653);  point_all_original = zeros(1,2653); point_all_zerophase = zeros(1,2653);
point_all_optimize_AFSQUARE = zeros(1,100000);  point_all_original_AFSQUARE = zeros(1,100000); point_all_zerophase_AFSQUARE = zeros(1,100000);
interval_width = 0.01;
for j = 1:length(u)
    for i = 1:length(v)
        value_all = PAT_zero_all(i,j); 
        if value_all >= 0 
            % 計算該值應該放入的區間
            interval_index = ceil(value_all / interval_width);
            % 更新對應區間的計數
            point_all_zerophase(interval_index) = point_all_zerophase(interval_index) + 1;
        end
        value_all = PAT_origin_all(i,j); 
        if value_all >= 0 
            % 計算該值應該放入的區間
            interval_index = ceil(value_all / interval_width);
            % 更新對應區間的計數
            point_all_original(interval_index) = point_all_original(interval_index) + 1;
        end
        value_all = PAT_optimize_all(i,j);
        if value_all >= 0 
            % 計算該值應該放入的區間
            interval_index = ceil(value_all / interval_width);
            % 更新對應區間的計數
            point_all_optimize(interval_index) = point_all_optimize(interval_index) + 1;
        end
    end
end
interval_width = 1;
for j = 1:length(u)
    for i = 1:length(v)
        value_all = PATAFSQUARE_zero_all(i,j); 
        if value_all >= 0 
            % 計算該值應該放入的區間
            interval_index = ceil(value_all / interval_width);
            % 更新對應區間的計數
            point_all_zerophase_AFSQUARE(interval_index) = point_all_zerophase_AFSQUARE(interval_index) + 1;
        end
        value_all = PATAFSQUARE_origin_all(i,j); 
        if value_all >= 0 
            % 計算該值應該放入的區間
            interval_index = ceil(value_all / interval_width);
            % 更新對應區間的計數
            point_all_original_AFSQUARE(interval_index) = point_all_original_AFSQUARE(interval_index) + 1;
        end
        value_all = PATAFSQUARE_optimize_all(i,j);
        if value_all >= 0 
            % 計算該值應該放入的區間
            interval_index = ceil(value_all / interval_width);
            % 更新對應區間的計數
            point_all_optimize_AFSQUARE(interval_index) = point_all_optimize_AFSQUARE(interval_index) + 1;
        end
    end
end

x = 0:1:2652; 
figure;
hold on
plot(x,point_all_optimize,'color',"#77AC30",'Marker','o','linewidth',1.2,'MarkerSize',1.6);%20231023註解：x,point_70_30degree,'r-o'表示一組曲線
plot(x,point_all_original,'b-o','linewidth',1.2,'MarkerSize',1.6)
plot(x,point_all_zerophase,'m-o','linewidth',1.2,'MarkerSize',1.6)
hold off
title('PDF of power are sampled with uv transform');
xlim([0 200])
%ylim([0 400])
xlabel('power(mag)');
ylabel('個數');
legend('最佳化後','最佳化前','最佳化(全0相位)');
grid on;

x_AF = 0:1:99999; 
figure;  %沒有除以377的power,0.01為一個間隔太密！
hold on
plot(x_AF,point_all_optimize_AFSQUARE,'color',"#00AC30",'Marker','o','linewidth',1.2,'MarkerSize',1.6);%20231023註解：x,point_70_30degree,'r-o'表示一組曲線
plot(x_AF,point_all_original_AFSQUARE,'b-o','linewidth',1.2,'MarkerSize',1.6)
plot(x_AF,point_all_zerophase_AFSQUARE,'r-','linewidth',1.2,'MarkerSize',1.6)
hold off
title('PDF of E-field square are sample with uv transform');
%xlim([0 200])
%ylim([0 400])
xlabel('power(mag)');
ylabel('個數');
legend('最佳化後','最佳化前','最佳化(全0相位)');
grid on;

PATdB_optimize_all = POWER_dB;  PATdB_origin_all = original_POWER_dB; PATdB_zero_all = zero_POWER_dB;
PATAFSQUAREdB_optimize_all = AF_square_dB;  PATAFSQUAREdB_origin_all = original_AF_square_dB; PATAFSQUAREdB_zero_all = zero_AF_square_dB;
for i = 2:201  %把南極點北極點重複的數值去掉
    PATdB_zero_all(1,i) = -60; PATdB_origin_all(1,i) = -100;   PATdB_optimize_all(1,i) = -100;
    PATdB_zero_all(201,i) = -100;   PATdB_origin_all(201,i) = -100;   PATdB_optimize_all(201,i) = -100;
    PATAFSQUAREdB_zero_all(1,i) = -100; PATAFSQUAREdB_origin_all(1,i) = -100;   PATAFSQUAREdB_optimize_all(1,i) = -100;
    PATAFSQUAREdB_zero_all(201,i) = -100;   PATAFSQUAREdB_origin_all(201,i) = -100;   PATAFSQUAREdB_optimize_all(201,i) = -100;
end
pointdB_all_optimize = zeros(1,100);  pointdB_all_original = zeros(1,100); pointdB_all_zerophase = zeros(1,100);
pointAFSQUAREdB_all_optimize = zeros(1,100);  pointAFSQUAREdB_all_original = zeros(1,100); pointAFSQUAREdB_all_zerophase = zeros(1,100);
value_all = 0;
interval_width = 1;
for j = 1:length(u)
    for i = 1:length(v)
        value_all = PATdB_zero_all(i,j); 
        if value_all >= -50
            % 計算該值應該放入的區間
            interval_index = floor(value_all / interval_width);
            interval_index = interval_index + 51;
            % 更新對應區間的計數
            pointdB_all_zerophase(interval_index) = pointdB_all_zerophase(interval_index) + 1;
        end
        value_all = PATdB_origin_all(i,j); 
        if value_all >= -50 
            % 計算該值應該放入的區間
            interval_index = floor(value_all / interval_width);
            interval_index = interval_index + 51;
            % 更新對應區間的計數
            pointdB_all_original(interval_index) = pointdB_all_original(interval_index) + 1;
        end
        value_all = PATdB_optimize_all(i,j);
        if value_all >= -50 
            % 計算該值應該放入的區間
            interval_index = floor(value_all / interval_width);
            interval_index = interval_index + 51;
            % 更新對應區間的計數
            pointdB_all_optimize(interval_index) = pointdB_all_optimize(interval_index) + 1;
        end
    end

    for i = 1:length(v)
        value_all = PATAFSQUAREdB_zero_all(i,j); 
        if value_all >= -50
            % 計算該值應該放入的區間
            interval_index = floor(value_all / interval_width);
            interval_index = interval_index + 51;
            % 更新對應區間的計數
            pointAFSQUAREdB_all_zerophase(interval_index) = pointAFSQUAREdB_all_zerophase(interval_index) + 1;
        end
        value_all = PATAFSQUAREdB_origin_all(i,j); 
        if value_all >= -50 
            % 計算該值應該放入的區間
            interval_index = floor(value_all / interval_width);
            interval_index = interval_index + 51;
            % 更新對應區間的計數
            pointAFSQUAREdB_all_original(interval_index) = pointAFSQUAREdB_all_original(interval_index) + 1;
        end
        value_all = PATAFSQUAREdB_optimize_all(i,j);
        if value_all >= -50 
            % 計算該值應該放入的區間
            interval_index = floor(value_all / interval_width);
            interval_index = interval_index + 51;
            % 更新對應區間的計數
            pointAFSQUAREdB_all_optimize(interval_index) = pointAFSQUAREdB_all_optimize(interval_index) + 1;
        end
    end
end
xdB = -50:1:49; 
figure;
hold on
plot(xdB,pointdB_all_optimize,'color',"#77AC30",'Marker','o','linewidth',1.2,'MarkerSize',1.6);%20231023註解：x,pointdB_70_30degree,'r-o'表示一組曲線
plot(xdB,pointdB_all_original,'b-o','linewidth',1.2,'MarkerSize',1.6)
plot(xdB,pointdB_all_zerophase,'m-o','linewidth',1.2,'MarkerSize',1.6)
hold off
title('PDF of power are sampled with uv transform');
%xlim([0 200])
%ylim([0 400])
xlabel('power(dB)');
ylabel('個數');
legend('最佳化後','最佳化前');
grid on;

figure;
hold on
plot(xdB,pointAFSQUAREdB_all_optimize,'color',"#77AC30",'Marker','o','linewidth',1.2,'MarkerSize',1.6);%20231023註解：x,pointdB_70_30degree,'r-o'表示一組曲線
plot(xdB,pointAFSQUAREdB_all_original,'y-o','linewidth',1.2,'MarkerSize',1.6)
plot(xdB,pointAFSQUAREdB_all_zerophase,'r-o','linewidth',1.2,'MarkerSize',1.6)
hold off
title('PDF of E-field square are sample with uv transform');
%xlim([0 200])
%ylim([0 400])
xlabel('power(dB)');
ylabel('個數');
legend('最佳化後','最佳化前','最佳化前(全0相位)');
grid on;

figure;
yyaxis left
plot(AF_variance_per_gen, 'LineWidth', 2);
ylabel('AF Variance Target');
xlabel('Generation');
yyaxis right
plot(AF_mean_per_gen, 'LineWidth', 2);
ylabel('mean');
title('AF Variance Target and mean over Generations');
legend('AF Variance Target', 'Mean');
grid on;

figure;
plot(penalty_per_gen, 'LineWidth', 2);
title('penalty over generations');
grid on;

figure;
hold on
plot(AF_variance_per_gen, 'LineWidth', 2);
plot(AF_mean_per_gen, 'LineWidth', 2);
plot(penalty_per_gen, 'LineWidth', 2);
plot(AF_variance_per_gen + 1./AF_mean_per_gen + penalty_per_gen,'LineWidth', 1.8)
hold off
title('AF Variance Target, mean, penalty, and cost over generations');
legend('AF Variance Target', 'Mean', 'penalty','cost');
grid on;
%------------------ Generate Initial Antenna Positions Function ------------------


function [initial_x, initial_y] = generateInitialPositions(num_antennas, lambda, x_range, y_range)
    %把初始群種的其中一個個體寫來這裡，位置都是排列整齊的。每個個體的間隔是以等差數列擴大，最密集的個體是間格0.4波長，公差為0.0098個波長，引數為i。最稀疏的間隔是佈滿30波長乘30波長的正方形。
    %一個個體中，波源以等差數列排列，引數為j。
    %這邊的個體會參與最佳化。
    initial_x = zeros(1,100); initial_y = zeros(1,100);
    for j = 1:1:100
        initial_x(j) = (-1.8*lambda) + 0.4*lambda*mod(j-1,10);
        initial_y(j) = (-1.8*lambda) + 0.4*lambda*floor((j-1)/10);   
    end
end

% ------------------ Generate Initial Population for GA ------------------

function initialPopulation = generateInitialPopulation(num_antennas, x_range, y_range, lambda, num_phases)
      popSize = 300; % Population size 初始的種群有300個
      initialPopulation = zeros(popSize, num_antennas * 3); % Each individual has phases and positions
      lamda = 0.01071412;
    %初始群種設定：
    for i = 1:popSize
        % 全部都是同樣的相位
        phaseIndices = zeros(1,100);
        x_positions = zeros(1,100); y_positions = zeros(1,100);
        for j = 1:1:100
            x_positions(j) = (-1.8*lamda-0.0441*(i-1)*lamda) + (0.4+0.0098*(i-1))*lamda*mod(j-1,10);
            y_positions(j) = (-1.8*lamda-0.0441*(i-1)*lamda) + (0.4+0.0098*(i-1))*lamda*floor((j-1)/10);   
        end

        % Combine into one individual
        initialPopulation(i, :) = [phaseIndices, x_positions, y_positions];
    end
end

% ------------------ Objective Function 成本函式------------------

function cost = objectiveFunction(params, num_antennas, k, theta_target_min, theta_target_max, phi_target_min, phi_target_max, lambda, possiblePhases_rad)
    global AF_variance_per_gen AF_mean_per_gen penalty_per_gen gen_count;    
     % Extract phase indices and x, y positions
    phase_indices = params(1:num_antennas);
    x_positions = params(num_antennas+1:2*num_antennas);
    y_positions = params(2*num_antennas+1:end);

    % Map phase indices to phase values in radians
    phases = possiblePhases_rad(round(phase_indices));

    % Compute complex phases
    w = exp(1i * phases);

    % Number of sampling points
    num_samples = 6000;

    % Generate random sampling points over the sphere
    u = rand(num_samples, 1);
    v = rand(num_samples, 1);

    % Compute corresponding theta and phi
    theta = acos(1 - 2 * u);    % θ ∈ [0, π]
    phi = 2 * pi * v - pi;      % φ ∈ [-π, π]

    % Compute array factor
    AF_values = zeros(num_samples, 1);
    for idx = 1:num_antennas
        psi = k * ( x_positions(idx) * sin(theta) .* cos(phi) + y_positions(idx) * sin(theta) .* sin(phi) );
        AF_values = AF_values + w(idx) * exp(1i * psi);
    end
    AF_values = abs(AF_values);
    AF_mean = mean(AF_values,"all");
    % Compute variance of AF in target angle range
    target_indices = (theta >= theta_target_min) & (theta <= theta_target_max);
    AF_values_target = AF_values(target_indices);

    %if isempty(AF_values_target) %判斷AF_values_target是不是空陣列，是的話AF的變異數就是0，不是的話就取變異數。
    %    AF_variance_target = 0;
    %else
    AF_variance_target = var(AF_values_target);
    %end

    %-------------新增權重，避免相位只有兩個----------------%
    % 設置目標theta範圍
    theta_min = deg2rad(30);  % 30 degrees
    theta_max = deg2rad(60);  % 50 degrees
    sigma = deg2rad(5);       % 控制範圍外的權重下降（可調整）

    % 使用權重函數來計算在目標theta範圍附近的權重
    %weights = some_weight_function(theta(target_indices), theta_min, theta_max, sigma);

    % 將權重應用到AF_values_target
    %AF_weighted_target = AF_values_target .* weights;

    % 計算加權後的方向圖變異性
    %AF_variance_weighted = var(AF_weighted_target);

    % Penalty for violating nearest neighbor constraints 
    penalty = 0;
    [c, ~] = distanceConstraint(params, num_antennas, lambda, [], []);
    violation = sum(c(c > 0)); % Sum of constraint violations
    penalty = penalty + 10 * violation; % Large penalty

    % Store values for the current generation
    AF_variance_per_gen(gen_count) = AF_variance_target;
    AF_mean_per_gen(gen_count) = AF_mean;
    penalty_per_gen(gen_count) = penalty;
    % Final cost function 基因演算法會保留適應值較大的個體，所以要把變異數和懲罰數倒數再相加。
    cost = AF_variance_target + 1/AF_mean + penalty;%到底要不要倒數?
    disp(['Test Generation：',num2str(gen_count)]);
end

% ------------------ Distance and Position Constraint Function ------------------

function [c, ceq] = distanceConstraint(params, num_antennas, lambda, x_range, y_range)
    x_positions = params(num_antennas+1:2*num_antennas);
    y_positions = params(2*num_antennas+1:end);

    c = [];
    ceq = [];

    % Get all antenna positions
    positions = [x_positions', y_positions'];

    % Compute distance matrix between antennas
    distance_matrix = squareform(pdist(positions));

    % For each antenna, find the nearest neighbor distance
    for i = 1:num_antennas
        distances_i = distance_matrix(i, :);
        distances_i(i) = inf; % Exclude self-distance
        min_distance = min(distances_i);

     
        c = [c; 0.004*sqrt(2) - min_distance]; 
        %c = [c; min_distance - (2 * lambda)]; 
    end

    % Add position constraints if x_range and y_range are provided
    if ~isempty(x_range)
        c = [c; x_positions' - x_range(2)];  % x <= x_max
        c = [c; x_range(1) - x_positions'];  % x >= x_min
    end
    if ~isempty(y_range)
        c = [c; y_positions' - y_range(2)];  % y <= y_max
        c = [c; y_range(1) - y_positions'];  % y >= y_min
    end
end

function [state, options, optchanged] = myOutputFunctionWithPower(options, state, flag)
    optchanged = false; % 不改變遺傳演算法的運行過程
    global AF_variance_target AF_mean_per_gen penalty_per_gen gen_count
    disp(['Standard : ', num2str(AF_variance_target)]);
    disp(['Mean: ', num2str(AF_mean_per_gen)]);
    disp(['Penalty: ', num2str(penalty_per_gen)]);
    disp(['cost:',num2str(AF_variance_target + 1./AF_mean_per_gen + penalty_per_gen)]);
    disp(['1/cost :',num2str(1./(AF_variance_target + 1./AF_mean_per_gen + penalty_per_gen))]);
    
    switch flag
        case 'iter'
            % 提取當前群體的適應值和群體
            scores = state.Score;
            population = state.Population;
            %diff = abs(state.Population - state.PreviousPopulation);
            % 找到適應值最佳的個體
            [bestScore, bestIdx] = min(scores);
            bestIndividual = population(bestIdx, :);

            % 顯示當代最佳個體及其適應值
            disp(['Generation: ', num2str(state.Generation)]);
            disp(['Best individual: ', num2str(bestIndividual)]);
            %disp(['Scores of this generation:',num2str(scores)]); %陣列維度不符合，需要再除錯
            disp(['Best fitness score: ', num2str(bestScore)]);
            disp('-----------------------------------');
            %disp(['delta with parents: ',diff]);
            % 根據最佳個體的資料計算AF平方
%             calculatePower(bestIndividual);  % 假設有一個計算功率的函數 (程式會卡死)
        case 'done'
            disp('GA process finished.');
    end
    if strcmp(flag, 'iter')  % Each generation iteration
        gen_count = gen_count + 1;
    end
end

Answer 1

[Swastik Sarkar]

Hi @Kiven,

The Genetic Algorithm appears to be running in parallel, as indicated by the following configuration:

ga_options = optimoptions('ga', 'Display', 'iter', ...
                          'MaxGenerations', 100, 'PopulationSize', 800, ...
                          'FunctionTolerance', 1e-7, 'ConstraintTolerance', 1e-7, ...
                          'UseParallel', true, 'MaxStallGenerations', 200, ...
                          'InitialPopulationMatrix', generateInitialPopulation(num_antennas, x_range, y_range, lambda, length(possiblePhases_deg)),...
                          'CreationFcn','gacreationuniformint',...
                          'SelectionFcn','selectiontournament',...
                          'MutationFcn','mutationgaussian',...
                          'CrossoverFcn','crossoverscattered',...
                          'OutputFcn', @myOutputFunctionWithPower,...
                          'PlotFcn',{'gaplotscorediversity','gaplotscores'});

When running in parallel, global variables are not synchronized across workers, leading to potential inconsistencies. To address this, consider disabling parallel execution by modifying the code as follows:

ga_options = optimoptions('ga', 'Display', 'iter', ...
                          'MaxGenerations', 100, 'PopulationSize', 800, ...
                          'FunctionTolerance', 1e-7, 'ConstraintTolerance', 1e-7, ...
                          'UseParallel', false, 'MaxStallGenerations', 200, ...
                          'InitialPopulationMatrix', generateInitialPopulation(num_antennas, x_range, y_range, lambda, length(possiblePhases_deg)),...
                          'CreationFcn','gacreationuniformint',...
                          'SelectionFcn','selectiontournament',...
                          'MutationFcn','mutationgaussian',...
                          'CrossoverFcn','crossoverscattered',...
                          'OutputFcn', @myOutputFunctionWithPower,...
                          'PlotFcn',{'gaplotscorediversity','gaplotscores'});

This change should ensure that the value of gen_count is displayed correctly. For further insights, refer to the following MATLAB Answer:

https://www.mathworks.com/matlabcentral/answers/215763-does-using-global-variables-in-objective-function-risk-error-when-using-ga-in-parallel

[Kiven]

Oh! Thank you for your reply and solution! I will try that, Thank you very much!