# How to optimize a parameter, influencing the dynamic response of a control system

4 views (last 30 days)
javier rivas on 23 Feb 2024
Commented: javier rivas on 26 Feb 2024
In the following program, I need to optimize the Kc value, so the variable ratio has a value of 4.
Kc=1;
s = tf('s');
G = 1/(s^2 + s + 1) * exp(-0.1 * s); %plant
GCL = feedback(Kc * G, 1); %closed loop response
t = linspace(0, 15);
u = ones(1, length(t));
out = lsim(GCL, u, t); %system response to step
peak = findpeaks(out); %peaks
ratio = (peak(1) - dcgain(GCL)) / (peak(2) - dcgain(GCL)); %Ratio of two first peaks

Aquatris on 23 Feb 2024
Here is a brute force way:
desired_ratio = 4;
Kc_vec = 0.1:0.01:10; % Kc search range
ratio = NaN(size(Kc_vec)); % initialize ratio vector
% loop through Kc_vec to see what the resulting ratio is for each Kc
for i = 1:length(Kc_vec)
Kc = Kc_vec(i);
s = tf('s');
G = 1/(s^2 + s + 1) * exp(-0.1 * s); %plant
GCL = feedback(Kc * G, 1); %closed loop response
t = linspace(0, 15);
u = ones(1, length(t));
out = lsim(GCL, u, t); %system response to step
peak = findpeaks(out); %peaks
ratio(i) = (peak(1) - dcgain(GCL)) / (peak(2) - dcgain(GCL)); %Ratio of two first peaks
end
tmp = abs(ratio-desired_ratio); % create a 'cost function'
idx = find(tmp==min(tmp));% find index where ratio is 4
figure(1)
plot(Kc_vec,ratio,Kc_vec(idx),ratio(idx),'x')
text(Kc_vec(idx),ratio(idx)+1,sprintf('Kc = %.2f ratio = %.2f',Kc_vec(idx),ratio(idx)));
xlabel('Kc'),ylabel('ratio')
javier rivas on 26 Feb 2024
Thank you very much for your help, I went through this different approach avoiding the use of a loop:
% Initial value of Kc
Kc= 1;
% Define the plant transfer function
s = tf('s');
G = 1/(s^2 + s + 1) * exp(-0.1 * s);
% Create the closed-loop function
GCL = feedback(Kc * G, 1);
% Define the objective function
objective_function = @(Kc) abs(ratio_objective(Kc, G));
% Choose the optimization algorithm
optimizer = @fminsearch; % or @fminunc
% Initial value of Kc for optimization
Kc_initial = Kc;
% Find the optimal Kc
[Kc_optimal, ~] = optimizer(objective_function, Kc_initial);
% Simulate the step response with the optimal Kc
t = linspace(0, 15);
u = ones(1, length(t));
out_optimal = lsim(feedback(Kc_optimal * G, 1), u, t);
% Calculate the peak ratio with the optimal Kc
peak_optimal = findpeaks(out_optimal);
ratio_optimal = (peak_optimal(1) - dcgain(feedback(Kc_optimal * G, 1))) / (peak_optimal(2) - dcgain(feedback(Kc_optimal * G, 1)));
% Display the results
fprintf('Optimal Kc: %2.3f ', Kc_optimal);
fprintf('\n Peak ratio with optimal Kc: %2.3f', ratio_optimal);
% Plot the step response
figure;
plot(t, out_optimal, 'b', 'LineWidth', 2);
grid on;
xlabel('Time (s)');
ylabel('Response');
title('Step Response with Optimal Kc');
% Define the objective function for the peak ratio
function value = ratio_objective(Kc, G)
GCL = feedback(Kc * G, 1);
t = linspace(0, 15);
u = ones(1, length(t));
out = lsim(GCL, u, t);
peak = findpeaks(out);
ratio = (peak(1) - dcgain(GCL)) / (peak(2) - dcgain(GCL));
value = abs(ratio - 4);
end

Sam Chak on 23 Feb 2024
I just wanted to point out that if you use the step() command instead of the lsim() command, the resulting output will be slightly different. However, I cannot provide a definitive answer regarding which of the two commands, step() or lsim(), is more accurate in generating the step response for time-delayed systems at this moment.
On another note, it's beneficial to familiarize yourself with @Aquatris's brute force searching approach as it offers valuable learning opportunities. Moreover, you can apply the methodology in Python as well.
klb = 1; % lower bound
kub = 3; % upper bound
[Kc, fval] = fminbnd(@cost, klb, kub)
Kc = 2.2308
fval = 8.7254e-12
s = tf('s');
G = 1/(s^2 + s + 1)*exp(-0.1*s);
GCL = feedback(Kc*G, 1);
t = linspace(0, 15, 1501);
out = step(GCL, t);
u = ones(1, length(t));
out2 = lsim(GCL, u, t);
TL = tiledlayout(2, 2, 'TileSpacing', 'Compact');
nexttile([1 2]);
plot(t, out), hold on
plot(t, out2), grid on, hold off
legend('step', 'lsim', 'location', 'SE')
nexttile([1 1]);
plot(t, out), hold on
plot(t, out2), grid on, hold off
xlim([1.8 2.0]), title('Peak 1')
nexttile([1 1]);
plot(t, out), hold on
plot(t, out2), grid on, hold off
xlim([5.3 5.8]), title('Peak 2')
title(TL, 'Step Responses')
xlabel(TL, 'Time (sec)')
peak = findpeaks(out); % peaks
ratio = (peak(1) - dcgain(GCL))/(peak(2) - dcgain(GCL)) % Ratio of two first peaks
ratio = 4.0000
function J = cost(Kc)
s = tf('s');
G = 1/(s^2 + s + 1)*exp(-0.1*s);
GCL = feedback(Kc*G, 1);
t = linspace(0, 15, 1501);
out = step(GCL, t);
peak = findpeaks(out);
ratio = (peak(1) - dcgain(GCL))/(peak(2) - dcgain(GCL));
J = (ratio - 4)^2;
end
javier rivas on 26 Feb 2024
Thank you very much for your suggestions, I finally considered a different approach creating a function and using the fminsearch command.
Regards