How to extracting parameters in iterations in genetic algorithms?
11 ビュー (過去 30 日間)
古いコメントを表示
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
0 件のコメント
回答 (1 件)
Swastik Sarkar
2024 年 10 月 21 日
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:
参考
カテゴリ
Help Center および File Exchange で Waveform Design and Signal Synthesis についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!