How to extracting parameters in iterations in genetic algorithms?

16 次查看（过去 30 天）
Kiven
Kiven 2024-10-21
评论： Kiven 2024-10-22，7:47
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

请先登录，再回答此问题。

回答（1 个）

Swastik Sarkar
Swastik Sarkar 2024-10-21，9:16
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

类别

在 Help Center 和 File Exchange 中查找有关 Waveform Design and Signal Synthesis 的更多信息

标签

产品


版本

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by