why function giving error?
显示 更早的评论
%% Clear previous data and close all figures
clear all
close all
clc
%% Phase 1: Input parameters
FoodSource = 100; % Food Sources Bees
D = 3; % Dimension of the problem (K, z, p)
lb = [0, 0, 0]; % Lower bound of the variables (K, z, p)
ub = [100, 10, 10]; % Upper bound of the variables (K, z, p)
max_iter = 100; % Maximum number of iterations/max no of cycle
N = FoodSource/2; % Population Size
limit = (N * D); % Used for Scout Phase
trial = zeros(N, 1); % Initialize to zero
% Simulink model name (custom Simulink model with the lead compensator)
simulink_model = 'lead_compensator_model';
%% Phase 2: Load or define the fitness function
% Note: The fitness function is already defined in the 'fitness_function.m' file.
%% Phase 3: Generation of initial population
pos = zeros(N, D);
for i = 1:N
for j = 1:D
pos(i, j) = lb(j) + rand() * (ub(j) - lb(j));
end
end
%% Phase 4: Calculate the fitness value for each solution in the initial population
fx = zeros(N, 1); % Initialize fitness values
for i = 1:N
% Evaluate fitness for each parameter set (row in "pos")
fx(i) = fitness_function(pos(i, :), simulink_model);
end
%% Main ABC loop
maxCycle = max_iter / N; % Maximum number of cycles
cycle = 0;
while cycle < maxCycle
% Employed bee phase
for i = 1:N
phi = -1 + 2 * rand(1, D); % Random parameter for mutation
vi = pos(i, :) + phi .* (pos(i, :) - pos(mod(i-1, N) + 1, :)); % Mutation
vi = max(vi, lb); % Ensure that the new solution is within bounds
vi = min(vi, ub);
fi = fitness_function(vi, simulink_model); % Calculate fitness for the new solution
if fi < fx(i) % Greedy selection
pos(i, :) = vi;
fx(i) = fi;
end
end
% Onlooker bee phase
prob = fx ./ sum(fx); % Probability for selecting each solution
onlookers = rand(N, 1); % Random numbers for selection
for i = 1:N
selected = find(onlookers <= cumsum(prob), 1, 'first'); % Roulette wheel selection
phi = -1 + 2 * rand(1, D); % Random parameter for mutation
vi = pos(selected, :) + phi .* (pos(selected, :) - pos(mod(selected-1, N) + 1, :)); % Mutation
vi = max(vi, lb); % Ensure that the new solution is within bounds
vi = min(vi, ub);
fi = fitness_function(vi, simulink_model); % Calculate fitness for the new solution
if fi < fx(selected) % Greedy selection
pos(selected, :) = vi;
fx(selected) = fi;
end
end
% Scout bee phase
[~, minIndex] = min(fx); % Find the index of the worst solution
if trial(minIndex) >= limit
% If the solution has not improved for "limit" cycles, create a new random solution
for j = 1:D
pos(minIndex, j) = lb(j) + rand() * (ub(j) - lb(j));
end
fx(minIndex) = fitness_function(pos(minIndex, :), simulink_model);
trial(minIndex) = 0; % Reset trial counter
else
trial(minIndex) = trial(minIndex) + 1; % Increment trial counter
end
% Update cycle count
cycle = cycle + 1;
end
%% Phase 5: Plot the result
plot(fx, 'r', 'LineWidth', 2);
xlabel('Iteration');
ylabel('TIAE Fitness Value');
grid on;
title('Fitness Value Convergence');
% Retrieve the best lead compensator parameters from the last iteration
best_fxval
3 个评论
Sam Chak
2023-8-1
Hi @ANDREW
Might be a bit of trouble to simulate the Simulink model lead_compensator_model.slx on the forum.
To test the ABC optimizer, please provide the Plant (transfer function or state-space also acceptable) that you are trying to control with the Lead Compensator.
Also specify the Performance Cost function, either a custom one, or the standard ones like ISE, ITAE, etc.
Rik
2023-8-7
I recovered the removed content from the Google cache (something which anyone can do). Editing away your question is very rude. Someone spent time reading your question, understanding your issue, figuring out the solution, and writing an answer. Now you repay that kindness by ensuring that the next person with a similar question can't benefit from this answer.
ANDREW
2023-8-19
采纳的回答
Hi @ANDREW
Your incomplete Artificial Bee Colony (ABC) algorithm is largely unmodified. I just added Phase #6 (to find the best pole for the Compensator based on ITAE) and a Fitness Function (based on ITAE) to supersede the unsupplied Simulink file. The Plant 
is exactly the same as shown in your block diagram.
is exactly the same as shown in your block diagram.%% Phase 1: Input parameters
FoodSource = 100; % Food Sources Bees
D = 1; % Dimension of the problem (K, z, p)
lb = 1.00; % Lower bound of the variables (K, z, p)
ub = 2.00; % Upper bound of the variables (K, z, p)
max_iter = 100; % Maximum number of iterations/max no of cycle
N = FoodSource/2; % Population Size
limit = N*D; % Used for Scout Phase
trial = zeros(N, 1); % Initialize to zero
% Simulink model name (custom Simulink model with the lead compensator)
% simulink_model = 'lead_compensator_model';
%% Phase 2: Load or define the fitness function
% Note: The fitness function is already defined in the 'fitness_function.m' file.
%% Phase 3: Generation of initial population
pos = zeros(N, D);
for i = 1:N
for j = 1:D
pos(i, j) = lb(j) + rand()*(ub(j) - lb(j));
end
end
%% Phase 4: Calculate the fitness value for each solution in the initial population
fx = zeros(N, 1); % Initialize fitness values
for i = 1:N
% Evaluate fitness for each parameter set (row in "pos")
fx(i) = fitfun(pos(i));
end
%% Main ABC loop
maxCycle = max_iter / N; % Maximum number of cycles
cycle = 0;
while cycle < maxCycle
% Employed bee phase
for i = 1:N
phi = -1 + 2*rand(1, D); % Random parameter for mutation
vi = pos(i, :) + phi.*(pos(i, :) - pos(mod(i-1, N) + 1, :)); % Mutation
vi = max(vi, lb); % Ensure that the new solution is within bounds
vi = min(vi, ub);
fi = fitfun(pos(i)); % Calculate fitness for the new solution
if fi < fx(i) % Greedy selection
pos(i, :) = vi;
fx(i) = fi;
end
end
% Onlooker bee phase
prob = fx./sum(fx); % Probability for selecting each solution
onlookers = rand(N, 1); % Random numbers for selection
for i = 1:N
selected = find(onlookers <= cumsum(prob), 1, 'first'); % Roulette wheel selection
phi = -1 + 2 * rand(1, D); % Random parameter for mutation
vi = pos(selected, :) + phi .* (pos(selected, :) - pos(mod(selected-1, N) + 1, :)); % Mutation
vi = max(vi, lb); % Ensure that the new solution is within bounds
vi = min(vi, ub);
fi = fitfun(pos(i)); % Calculate fitness for the new solution
if fi < fx(selected) % Greedy selection
pos(selected, :) = vi;
fx(selected) = fi;
end
end
% Scout bee phase
[~, minIndex] = min(fx); % Find the index of the worst solution
if trial(minIndex) >= limit
% If the solution has not improved for "limit" cycles, create a new random solution
for j = 1:D
pos(minIndex, j) = lb(j) + rand()*(ub(j) - lb(j));
end
fx(minIndex) = fitfun(pos(minIndex, :));
trial(minIndex) = 0; % Reset trial counter
else
trial(minIndex) = trial(minIndex) + 1; % Increment trial counter
end
% Update cycle count
cycle = cycle + 1;
end
%% Phase 5: Plot the result
plot(fx, 'r.', 'LineWidth', 2);
xlabel('Iteration');
ylabel('TIAE Fitness Value');
grid on;
title('Fitness Value Convergence');
% Phase 6
idx1 = find(fx == min(fx));
bestp = pos(idx1);
for k = 1:length(bestp)
J(k) = fitfun(bestp(k));
end
idx2 = find(J == min(J));
Best_Pole = - bestp(idx2)
Gp = tf(1, [1 0]) % Plant
Gc = zpk([], Best_Pole, 1) % Compensator
Gcl = tf(feedback(Gc*Gp, 1)) % Closed-loop System
step(Gcl, 20), grid on
% Fitness Function
function J = fitfun(p)
Gp = tf(1, [1 0]); % Plant
Gc = zpk([], -p, 1); % Compensator
Gcl = minreal(feedback(Gc*Gp, 1)); % Closed-loop TF
[y, t] = step(Gcl, 20);
err = y(end) - y;
ifcn = t.*abs(err);
J = trapz(t, ifcn);
end
8 个评论
ANDREW
2023-8-6
Sam Chak
2023-8-6
You are welcome, @ANDREW. If you find the solution helpful, please consider accepting ✔ and voting 👍 on the answer. Thanks a bunch!
By the way, I originally used the 3-parameter config of your compensator, 
, but it didn't go well. After reviewing your plant, 
is an integrator, I determined that a zero is unnecessary, and 
. If K is a free parameter, then the minimum cost occurs at 
, which is impractical. To improve the convergence time, you can select a higher value of 
, then let the Bee algorithm finds the best pole, p.
, but it didn't go well. After reviewing your plant, is an integrator, I determined that a zero is unnecessary, and . If K is a free parameter, then the minimum cost occurs at , which is impractical. To improve the convergence time, you can select a higher value of , then let the Bee algorithm finds the best pole, p.
ANDREW
2023-8-6
Hi @ANDREW, Are you trying to fit some real values to the parameters 
in the Compensator 
to control the Plant 
based on the Mathematical Control Theory?
in the Compensator to control the Plant based on the Mathematical Control Theory? The Artificial Bee Colony algorithm only searches, but it does not know whether it is good to use the proposed Compensator or not because its programmed algorithm has no knowledge about Control Theory.
If there are no constraints, the Bee algorithm might just find 
so that the zero and pole cancel each other out, resulting in a Proportional controller, 
.
so that the zero and pole cancel each other out, resulting in a Proportional controller, .
ANDREW
2023-8-6
ANDREW
2023-8-6
Sam Chak
2023-8-6
Hi @ANDREW, you can try to fit them using ABC algorithm. After all, exploring options is a crucial approach in research to understand the behavior of a system, as well as the limitations of ABC. You won't know if you didn't try it out.
Could you please try modifying the code for your proposed Compensator and then insert it here 
? If there is an error, we can troubleshoot it here.
? If there is an error, we can troubleshoot it here.
ANDREW
2023-8-7
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Physics 的更多信息
Translated by 
