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

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

Answer 1

[Aquatris]

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]

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

Answer 2

[Sam Chak]

Hi @javier rivas

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]

Thank you very much for your suggestions, I finally considered a different approach creating a function and using the fminsearch command.

Regards