G =
I haven't defined syms x but solve function gives me values depend on x
Mostrar comentarios más antiguos
I am designing a PID controller for which I set the denominator of my closed-loop transfer function equal to the product of the dominant poles and the residue polynomial.
My goal is to get values that depend on Ki, but even though I don't use the syms x command, my values depend on x and Ki.
clear all;
clc;
syms zeta positive
syms s;
syms Kp Ki Kd;
% Given values
ts = 4; % Settling time
Mp = 0.05; % Maximum overshoot
% Define the overshoot equation
eq1 = exp(-zeta*pi / sqrt(1 - zeta^2)) == Mp;
% Solve for zeta numerically (only real positive solution makes sense)
zeta = double(vpasolve(eq1, zeta, [0 1]));
% Now calculate omega_n using Ts = 4 / (zeta * omega_n)
w_n = 4 / (zeta * ts);
G= (s+1)/(s^3 +6*s^2 +10*s -15)
H= 10/(s+10)
F= (Kd*s^2 + Kp*s + Ki)/s
Tcl= (F*G)/(1 + F*G*H)
[num, den] = numden(simplifyFraction(Tcl));
num = expand(num)
den = expand(den)
num_s = coeffs(num, s, 'All')
den_s = coeffs(den, s, 'All')
Pd=s^2 +2*zeta*w_n*s+(w_n)^2
syms e0 e1 e2 positive;
Pr= (s+e0)*(s+e1)*(s+e2)
Pe=expand(Pd*Pr)
Pe_s= coeffs(Pe,s,"All")
eq2 = Pe_s == den_s
sol2 = solve(eq2, [Kp, Kd, e0, e1, e2], 'Real', true);
% den_s=subs(den_s,[Kd,Kp],[sol2.Kd,sol2.Kp])
% F = subs(F,[Kd,Kp],[sol2.Kd,sol2.Kp])
Respuesta aceptada
Walter Roberson
el 3 de Mayo de 2025
The solutions are most naturally parameterized by e2, but you have specified that you want to solve for e2. So solve() has to invent a new variable 'x' and parameterize with respect to 'x' and then say "oh, and e2 is x"
sol2 = solve(eq2, [Kp, Kd, e0, e1, e2], 'Real', true, 'ReturnConditions', true);
subs([Kp, Kd, e0, e1, e2], subs(sol2, sol2.parameters, e2))
gives
5 comentarios
Cahit Semih
el 3 de Mayo de 2025
Remove the restrictions that e0 e1 e2 are positive, and remove the 'real', true restriction
When you plot the results, you will see that in each case, one of e0, e1, or e2 is never positive at the beginning. In the first four plots, the never-positive variable shows up again after a gap (the gap represents sections where the variable is complex-valued) -- however, its companion variable has vanished. In the last two plots, the never-positive variable does not show up again after a gap.
So... it appears that there are no combinations of values that make e0, e1, e2 simultaneously positive.
syms s;
syms Kp Ki Kd;
% Given values
ts = 4; % Settling time
Mp = 0.05; % Maximum overshoot
syms zeta
% Define the overshoot equation
eq1 = exp(-zeta*pi / sqrt(1 - zeta^2)) == Mp;
% Solve for zeta numerically (only real positive solution makes sense)
zeta = double(vpasolve(eq1, zeta, [0 1]));
% Now calculate omega_n using Ts = 4 / (zeta * omega_n)
w_n = 4 / (zeta * ts);
G= (s+1)/(s^3 +6*s^2 +10*s -15);
H= 10/(s+10);
F= (Kd*s^2 + Kp*s + Ki)/s;
Tcl= (F*G)/(1 + F*G*H);
[num, den] = numden(simplifyFraction(Tcl));
num = expand(num);
den = expand(den);
num_s = coeffs(num, s, 'All');
den_s = coeffs(den, s, 'All');
Pd=s^2 +2*zeta*w_n*s+(w_n)^2;
syms e0 e1 e2
Pr= (s+e0)*(s+e1)*(s+e2);
Pe=expand(Pd*Pr);
Pe_s= coeffs(Pe,s,"All");
eq2 = Pe_s == den_s
sol2 = solve(eq2, [Kp, Kd, e0, e1, e2], 'real', false, 'maxdegree', 3)
result = subs([Kp, Kd, e0, e1, e2], sol2)
The challenge after this is to isolate the ranges of Ki that make everything come out real-valued and make e0 e1 e2 come out positive
format long g
Kilin = linspace(-20,100,501).';
for K = 1 : 6
figure(K)
%y = double(subs(result(K,:),Ki,Kilin.'));
y = double(subs(result(K,:), Ki, Kilin));
mask = abs(imag(y)) < 1e-12;
y(mask) = real(y(mask));
y(~mask) = nan;
plot(Kilin, y);
xline(0, 'k:')
yline(0, 'k:')
legend({'Kp', 'Kd', 'e0', 'e1', 'e2'}, 'location', 'northwest')
end
I actually do not understand @Cahit Semih's approach. From his presentation, this math problem appears to be specifically about solving simultaneous equations. In Part 1, when the three PID variables are symbolically inserted into the closed-loop formula, the denominator is a polynomial of degree 5. It is clear that the three PID variables can only affect four terms: {s³, s², s¹, s⁰}, as the coefficient of s⁴ is a constant (16). Moreover, the three PID variables can only satisfy the conditions for both the s³ and s⁰ terms, and either one of the two terms, s² or s¹.
Next, in Part 2, a dominant polynomial of degree 2 is created after the desired zeta and 
are determined. Then, a residue polynomial of degree 3 is introduced with three symbolic roots {-e₁, -e₂, -e₃}, which is multiplied by the dominant polynomial, resulting in an expanded polynomial of degree 5. In Part 3, all six terms in the expanded polynomial are then compared with the denominator from Part 1, term by term. Regardless of the values assigned to the roots, the three PID variables still cannot change the coefficient of s⁴.
are determined. Then, a residue polynomial of degree 3 is introduced with three symbolic roots {-e₁, -e₂, -e₃}, which is multiplied by the dominant polynomial, resulting in an expanded polynomial of degree 5. In Part 3, all six terms in the expanded polynomial are then compared with the denominator from Part 1, term by term. Regardless of the values assigned to the roots, the three PID variables still cannot change the coefficient of s⁴.In Part 4, I attempt to determine the P, I, and D values from the Hurwitz matrix (Pascal Triangle), which is known to produce negative real parts of the roots for the polynomial (or positive real values for e₁, e₂, and e₃). Only the settling time (under 4 seconds) is satisfied.
%% PART 1: Find the denominator polynomial of a Closed-loop Transfer function
% --------------------------------
syms s
syms K [1 3]
% G rational function
G = (s + 1)/(s^3 + 6*s^2 + 10*s - 15);
% H rational function
H = 10/(s + 10);
% PID control equation
C = K1/s + K2 + K3*s; % K1 = Ki, K2 = Kp, K3 = Kd
% Closed-loop Transfer function
Tcl = C*G/(1 + C*G*H);
Tcl = simplifyFraction(Tcl)
[~, den] = numden(Tcl);
[den, term] = coeffs(den, s, 'All') % polynomial degree 5
%% PART 2: Convolute two polynomials such that it has the highest degree of the Closed-loop Transfer function
% --------------------------------
syms zeta
syms e [1 3]
sympref('FloatingPointOutput', true);
% Desired performance requirements
ts = 4; % Settling time
Mp = 0.05; % Maximum overshoot% Desired performance requirements
% Find dominant zeta and wn
eq1 = exp(-zeta*pi / sqrt(1 - zeta^2)) == Mp; % overshoot equation
zeta = double(vpasolve(eq1, zeta, [0 1])); % 0.6901
wn = 4/(zeta*ts); % 1.4491
% Dominant polynomial (quadratic polynomial)
Pd = s^2 + 2*zeta*wn*s + wn^2;
% Residue polynomial (cubic polynomial)
Pr = (s + e1)*(s + e2)*(s + e3);
% Expanded polynomial (quintic polynomial)
Pe = expand(Pd*Pr)
[Pe, term] = coeffs(Pe, s, "All")
%% PART 3: Identify simultaneous equations before solving them
% --------------------------------
eq2 = Pe == den;
Eq1 = eq2(1)
Eq2 = eq2(2)
Eq3 = eq2(3)
Eq4 = eq2(4)
Eq5 = eq2(5)
Eq6 = eq2(6)
%% PART 4: Determine P, I, D values from Hurwitz Matrix (not from e1, e2, e3)
% --------------------------------
m = 5;
HurwitzMatrix = sym(zeros(m + 1, m + 1));
for n = 0:m
for k = 0:n
HurwitzMatrix(n + 1, k + 1) = nchoosek(n, k);
end
end
disp('Hurwitz Matrix:'); disp(HurwitzMatrix);
% Hurwitz characteristic polynomial for degree 5
p5 = double([HurwitzMatrix(6, :)])
w = 16/p5(2); % s⁴ term, wn can be determined
Ki = p5(6)*w^5/10; % s⁰ term, Ki can be determined
Kd = (p5(3)*w^2 - 70)/10; % s³ term, Kd can be determined
Kp = (p5(5)*w^4 + 150 - 10*Ki)/10; % s¹ term, option 1 to get Kp
% Kp = (p5(4)*w^3 - 85 - 10*Kd)/10; % s² term, option 2 (different Kp value)
C = pid(Kp, Ki, Kd)
s = tf('s');
G = (s + 1)/(s^3 + 6*s^2 + 10*s - 15);
H = 10/(s + 10);
Tcl = feedback(series(C, G), H)
pole(Tcl)
figure
step(Tcl)
grid on
stepinfo(Tcl)
Cahit Semih
el 5 de Mayo de 2025
Hi @Cahit Semih
The step response of the closed-loop 5th-order TF will resemble that of the dominant TF only if the roots -e₀, -e₁, and -e₂ are freely selected, allowing the residue TF to achieve a very fast settling time. At the same time, the selection of the roots is constrained by your PID-induced condition:
It is important to note that the 2nd-order TF is legitimately referred to as 'dominant' if and only if the closed-loop TF behaves accordingly.
Is there a rationale for making the gains 
and 
dependent on another gain 
? Since the parameters 
and are arbitrarily selected under certain conditions, I am curious to know how you would apply @Walter Roberson's solution from the Answer to determine the gains 
and 
in your upcoming sequel. At least, from the last term of the polynomial, 
, we know that 
.
and dependent on another gain ? Since the parameters and are arbitrarily selected under certain conditions, I am curious to know how you would apply @Walter Roberson's solution from the Answer to determine the gains and in your upcoming sequel. At least, from the last term of the polynomial, , we know that .In the end, you will need to determine whether the proposed pure PID controller is algebraically sufficient to stabilize the system and satisfy the performance criteria (
and 
). Please keep us posted.
and ). Please keep us posted.%% Dominant TF (2nd-order)
zeta= 0.690106730559822;
wn = 1.49884262536285;
G2 = tf(wn^2, [1, 2*zeta*wn, wn^2])
[~, dG2] = tfdata(G2, 'v')
S2 = stepinfo(G2);
figure
step(G2, 7), grid on
title('Performance Criteria')
xline(S2.SettlingTime, '--', sprintf('Settling Time: %.2f sec', S2.SettlingTime), 'color', '#7F7F7F', 'LabelVerticalAlignment', 'bottom');
yline(S2.Peak, '--', sprintf('Max Overshoot: %.2f %%', S2.Overshoot), 'color', '#7F7F7F', 'LabelHorizontalAlignment', 'left');
%% Residue TF (3rd-order)
e2 = (16 - dG2(2))/3; % selected to satisfy the condition
e1 = e2;
e0 = e2;
G3 = tf(zpk([], [-e0, -e1, -e2], e0*e1*e2))
%% Expanded TF (5th-order Closed-loop)
G5 = series(G3, G2)
S5 = stepinfo(G5);
figure
step(G5, 7), grid on
title('Desired Closed-loop 5th-order transfer function')
xline(S5.SettlingTime, '--', sprintf('Settling Time: %.2f sec', S5.SettlingTime), 'color', '#7F7F7F', 'LabelVerticalAlignment', 'bottom');
yline(S5.Peak, '--', sprintf('Max Overshoot: %.2f %%', S5.Overshoot), 'color', '#7F7F7F', 'LabelHorizontalAlignment', 'left');
Más respuestas (1)
Hi @Cahit Semih
Do you have any updates regarding your optimization efforts?
The following solution is derived from your eigenvalue injection approach. However, instead of introducing three symbolic eigenvalues (resulting in more unknowns than the available useful equations), only one parameter is introduced to represent the three eigenvalues. This slight modification is not my typical approach for control design, but it helps you in illustrating a mathematical method for solving problems.
%% PART 1: Find the denominator polynomial of a Closed-loop Transfer function
% --------------------------------
syms s
syms K [1 3]
% G rational function
G = (s + 1)/(s^3 + 6*s^2 + 10*s - 15);
% H rational function
H = 10/(s + 10);
% PID control equation
C = K1/s + K2 + K3*s; % K1 = Ki, K2 = Kp, K3 = Kd
% Closed-loop Transfer function
Tcl = C*G/(1 + C*G*H);
Tcl = simplifyFraction(Tcl)
[~, den] = numden(Tcl);
[den, term] = coeffs(den, s, 'All') % polynomial degree 5
%% PART 2: Convolute two polynomials such that it has the highest degree of the Closed-loop Transfer function
% --------------------------------
syms zeta
syms e
sympref('FloatingPointOutput', true);
% Desired performance requirements
ts = 4; % Settling time
Mp = 0.05; % Maximum overshoot% Desired performance requirements
% Find dominant zeta and wn
eq1 = exp(-zeta*pi / sqrt(1 - zeta^2)) == Mp; % overshoot equation
zeta = double(vpasolve(eq1, zeta, [0 1])); % 0.6901
wn = 4/(zeta*ts); % 1.4491
% Dominant polynomial (quadratic polynomial)
Pd = s^2 + 2*zeta*wn*s + wn^2;
%% Critical part in the entire control design problem
% Residue polynomial (cubic polynomial)
% Pr = (s + e)^3; % It works, but no design parameter for user to manipulate
img = 15.9063; % introduce a design parameter (imaginary coefficient)
Pr = (s + e)*(s + e + img*1i)*(s + e - img*1i);
% Expanded polynomial (quintic polynomial)
Pe = expand(Pd*Pr);
[Pe, term] = coeffs(Pe, s, "All")
%% PART 3: Identify 4 useful equations to simultaneously solve them for 4 unknowns {Kp, Ki, Kd, e}
% --------------------------------
eq2 = Pe == den;
Eq1 = eq2(1) % not useful
Eq2 = eq2(2) % constraint, but not useful, so we choose not to comply
Eq3 = eq2(3) % 1st equation
Eq4 = eq2(4) % 2nd equation
Eq5 = eq2(5) % 3rd equation
Eq6 = eq2(6) % 4th equation
Eqs = [Eq3, Eq4, Eq5, Eq6];
% Solve 4 equations for 4 unknowns
sol = solve(Eqs, [K1, K2, K3, e], 'Real', true, 'ReturnConditions', true)
%% PART 4: Obtain the Closed-loop TF and clean up if necessary
% --------------------------------
% the value e = 1.6131 is not used, it is a by-product.
% no matter what, the coefficient of s⁴ is always 16.
% the pure PID controller
Kp = double(sol.K2);
Ki = double(sol.K1);
Kd = double(sol.K3);
C = pid(Kp, Ki, Kd)
s = tf('s');
% the Plant
G = (s + 1)/(s^3 + 6*s^2 + 10*s - 15)
% the Sensor
H = 10/(s + 10)
% Closed-loop TF (too many zeros affect the overshoot, see if they are cancellable)
Tcl = feedback(series(C, G), H)
S1 = stepinfo(Tcl);
if zero(Tcl) < 0
[ncl, dcl] = tfdata(Tcl, 'v');
end
% if the zeros are cancellable, then design a Pre-filter
Gf = tf(dcl(end), ncl)
% Filtered closed-loop (no more disturbing zeros, very 'CLEAN' now)
Fcl = minreal(series(Gf, Tcl))
S2 = stepinfo(Fcl)
%% PART 5: Simulation
% --------------------------------
figure
hold on
step(Tcl, 2*round(S2.SettlingTime))
step(Fcl, 2*round(S2.SettlingTime))
hold off
legend('Tcl', 'Fcl', 'location', 'east')
grid on
8 comentarios
Cahit Semih
el 9 de Mayo de 2025
Editada: Sam Chak
el 9 de Mayo de 2025
Obtain the PID controller that maximizes the dominance of the residual poles. Give the parametric expression of the controller. Give all the necessary steps.The form you need to use:
m_pole = |(min|Re(pole_residue)|)/sigma|
sigma = real part of the dominant pole.
My friend suggested this solution for this question, what do you think?
clear all; clc; close all;
syms s kp ki kd a b c real
syms zeta
ts = 4;
Mp = 0.05;
% Zeta & wn
zeta = double(vpasolve(exp(-zeta*pi/sqrt(1 - zeta^2)) == Mp, zeta, [0 1]));
wn = 4/(zeta*ts);
Pd = s^2 + 2*zeta*wn*s + wn^2;
% system
Gs = tf([1 1],[1 6 10 -15]);
Gss = collect(expand((s+1)/(s^3 +6*s^2 +10*s -15)));
Hss = 10/(s+10);
% PID
Fss = (kd*s^2 + kp*s + ki)/s;
%close loop
Tss = (Gss*Fss)/(1 + Gss*Fss*Hss);
[~, den] = numden(Tss);
% characteristic eq
Pr = s^3 + a*s^2 + b*s + c;
target_poly = expand(Pd * Pr);
eqs = coeffs(den, s, 'All') == coeffs(target_poly, s, 'All');
sol = solve(eqs, [a, b, c, kp, kd]);
disp('Kp(Ki):'); vpa(simplify(sol.kp))
disp('Kd(Ki):'); vpa(simplify(sol.kd))
%% Q2 – Increase the dominance of residual poles
fun_kp = matlabFunction(sol.kp);
fun_kd = matlabFunction(sol.kd);
fun_b = matlabFunction(sol.b);
zwn = zeta*wn;
ki_vec = logspace(1.5, 3, 100); %Ki on stable interval
eta_vals = zeros(size(ki_vec));
for i = 1:length(ki_vec)
ki_val = ki_vec(i);
kp_val = fun_kp(ki_val);
kd_val = fun_kd(ki_val);
b_val = fun_b(ki_val);
% residue characteristic equation
res_poly = double(coeffs(s^3 + 14*s^2 + b_val*s + ki_val, s, 'All'));
res_roots = roots(res_poly);
eta_vals(i) = min(abs(real(res_roots))) / zwn; % m_pole
end
% Best Ki for the most dominant residue pole
[max_eta, idx] = min(eta_vals);
best_ki = ki_vec(idx);
best_kp = fun_kp(best_ki);
best_kd = fun_kd(best_ki);
fprintf('En iyi Ki: %.4f, Kp: %.4f, Kd: %.4f, m_kutup: %.4f\n', best_ki, best_kp, best_kd, max_eta);
% Transfer function
Fs_best = tf([best_kd best_kp best_ki], [1 0]);
Ts_best = minreal(feedback(Gs * Fs_best, tf(10, [1 10])));
figure; step(Ts_best); title('Step Response ');
figure; pzmap(Fs_best*Gs); title('Poles and Zeros – PID');
%
figure;
semilogx(ki_vec, eta_vals, 'LineWidth', 2);
hold on; plot(best_ki, max_eta, 'ro', 'MarkerSize', 8);
xlabel('Ki'); ylabel('m_{pole}'); grid on;
title('Dominance coefficient of residual poles');
Hi @Cahit Semih
I have reviewed your code, and the approach may be acceptable for undergraduates. However, in the code, you neglected to use the Routh–Hurwitz technique to determine the minimum 
gain that stabilizes the system. It appears that you have simply calculated the gain from the solution of b (which you thought to be the minimum), resulting in 
or 
. However, this gain does not stabilize the system.
gain that stabilizes the system. It appears that you have simply calculated the gain from the solution of b (which you thought to be the minimum), resulting in or . However, this gain does not stabilize the system.Additionally, you should observe the changes in the residue poles as the value of is increased. Without the prefilter, and at high values, the imaginary part becomes significant, causing the transient performance to be oscillatory. As I mentioned previously, the sum of the real parts of the residue poles are constrained to 16 and cannot be increased arbitrarily. To satisfy the equations, the imaginary parts are increased instead, leading to unwanted oscillations.
is increased. Without the prefilter, and at high values, the imaginary part becomes significant, causing the transient performance to be oscillatory. As I mentioned previously, the sum of the real parts of the residue poles are constrained to 16 and cannot be increased arbitrarily. To satisfy the equations, the imaginary parts are increased instead, leading to unwanted oscillations.If you find the explanation and the added code helpful in this comment, please consider voting 👍 for the solution presented in my answer.
b: (22517998136852480*ki)/4728215395125409 - 373130563904622547/2476415581440161
clear all; clc; close all;
syms s kp ki kd a b c real
syms zeta
ts = 4;
Mp = 0.05;
% Zeta & wn
zeta = double(vpasolve(exp(-zeta*pi/sqrt(1 - zeta^2)) == Mp, zeta, [0 1]));
wn = 4/(zeta*ts);
Pd = s^2 + 2*zeta*wn*s + wn^2
% system
Gs = tf([1 1],[1 6 10 -15]);
Gss = collect(expand((s+1)/(s^3 +6*s^2 +10*s -15)));
Hss = 10/(s+10);
% PID
Fss = (kd*s^2 + kp*s + ki)/s;
%close loop
Tss = (Gss*Fss)/(1 + Gss*Fss*Hss);
[~, den] = numden(Tss)
% characteristic eq
Pr = s^3 + a*s^2 + b*s + c;
target_poly = expand(Pd * Pr);
eqs = coeffs(den, s, 'All') == coeffs(target_poly, s, 'All');
sol = solve(eqs, [a, b, c, kp, kd])
disp('Kp(Ki):'); vpa(simplify(sol.kp))
disp('Kd(Ki):'); vpa(simplify(sol.kd))
%% Not advisable to use this Routh–Hurwitz technique if not taught by your Professor
B = sol.b/sol.c^(2/3);
A = sol.a/sol.c^(1/3);
sol2 = solve(A*B - 1 == 0, ki, 'ReturnConditions', true);
min_Ki = double(sol2.ki)
initKi = log10(min_Ki)
%% Q2 – Increase the dominance of residual poles
fun_kp = matlabFunction(sol.kp);
fun_kd = matlabFunction(sol.kd);
fun_b = matlabFunction(sol.b);
zwn = zeta*wn;
% ki_vec = logspace(1.5, 3, 100) % Ki on stable interval
ki_vec = logspace(initKi, 4, 10); % true Ki on stable interval (10 points)
eta_vals = zeros(size(ki_vec));
for i = 1:length(ki_vec)
ki_val = ki_vec(i);
kp_val = fun_kp(ki_val);
kd_val = fun_kd(ki_val);
b_val = fun_b(ki_val);
% residue characteristic equation
res_poly = double(coeffs(s^3 + 14*s^2 + b_val*s + ki_val, s, 'All'));
res_roots = roots(res_poly)
eta_vals(i) = min(abs(real(res_roots))) / zwn; % m_pole
end
% Best Ki for the most dominant residue pole
[max_eta, idx] = min(eta_vals);
best_ki = ki_vec(idx);
best_kp = fun_kp(best_ki);
best_kd = fun_kd(best_ki);
fprintf('En iyi Ki: %.4f, Kp: %.4f, Kd: %.4f, m_kutup: %.4f\n', best_ki, best_kp, best_kd, max_eta);
% Transfer function
Fs_best = tf([best_kd best_kp best_ki], [1 0]);
Ts_best = minreal(feedback(Gs * Fs_best, tf(10, [1 10])))
% check the zeros
zero(Ts_best)
figure; step(Ts_best); title('Step Response ');
figure; pzmap(Fs_best*Gs); title('Poles and Zeros – PID');
figure;
semilogx(ki_vec, eta_vals, 'LineWidth', 2);
hold on; plot(best_ki, max_eta, 'ro', 'MarkerSize', 8);
xlabel('Ki'); ylabel('m_{pole}'); grid on;
title('Dominance coefficient of residual poles');
%% Should add a Pre-filter since all zeros have negative real parts
[ncl, dcl] = tfdata(Ts_best, 'v');
% if the zeros are cancellable, then design a Pre-filter
Gf = tf(dcl(end), ncl)
% Filtered closed-loop (no more disturbing zeros, very 'CLEAN' now)
Fcl = minreal(series(Gf, Ts_best))
S2 = stepinfo(Fcl)
figure
subplot(211)
[y1, t1] = step(Ts_best, 2*round(S2.SettlingTime));
plot(t1, y1)
title('Closed-loop TF')
grid on
subplot(212)
[y2, t2] = step(Fcl, 2*round(S2.SettlingTime));
plot(t2, y2)
title('Filtered closed-loop TF')
grid on
Cahit Semih
el 9 de Mayo de 2025
Sam Chak
el 9 de Mayo de 2025
The codes for the roots and the coefficients of the residue polynomial are all correct. However, you neglected to determine the minimum value of the 
gain using the Routh–Hurwitz technique. Instead, you directly used 
(without verifying its stability) in the creation of the 
vector for the for-loop. This value of 
leads to instability. Consequently, I have added my computation of the lowest value of 
gain that stabilizes the system in the code, using the lesser-known Routh–Hurwitz technique.
gain using the Routh–Hurwitz technique. Instead, you directly used (without verifying its stability) in the creation of the vector for the for-loop. This value of leads to instability. Consequently, I have added my computation of the lowest value of gain that stabilizes the system in the code, using the lesser-known Routh–Hurwitz technique.In summary, use
ki_vec = logspace(initKi, 3, 100); % true Ki on stable interval
instead of
ki_vec = logspace(1.5, 3, 100) % Ki on stable interval
Cahit Semih
el 9 de Mayo de 2025
Torsten
el 9 de Mayo de 2025
Make a google search:
matlab & Routh Hurwitz stability
There are several codes you could try from the File Exchange.
Walter Roberson
el 9 de Mayo de 2025
Is there a function in matlab that examines the stability of polynomials containing symbolic expressions, especially for Routh-Hurwitz stability criteria?
No, there is not. The Control Systems Toolbox is purely numeric and does not interact with the Symbolic Toolbox at all, and the Symbolic Toolbox does not have an control systems functionality built in.
Hi @Cahit Semih
Typically, we do not use the Routh–Hurwitz technique to determine the stability of a SISO linear system, as directly finding the eigenvalues is more computationally efficient. When there is only one or two control parameters, the Routh–Hurwitz technique can be employed to find the range of values for the stability region. This allows the designer to select values and manually tune the control gains until satisfactory performance is achieved.
If the control problem is well-conditioned, it is entirely possible to determine the unique control gains using the algebraic approach. However, in your case, the challenge of controlling a 5-parameter closed-loop transfer function of degree 5 is significantly constrained by the 3-parameter PID controller. The dominant and residue poles approach does not apply here because the selection of the 3 residue poles is limited by the requirement that 
. In order to ensure that the dominant poles to significantly influence the system's response, the residue poles must be at least 100 times faster than the real parts of the dominant poles.
. In order to ensure that the dominant poles to significantly influence the system's response, the residue poles must be at least 100 times faster than the real parts of the dominant poles.Your problem-solving approach is similar to a 'brute-force search,' in which you conduct an exhaustive search of values for the 
gain. However, no values of gain will yield performance that is close to Ts = 4 and OS = 5%, unless you relax the requirements to Ts = 4 and OS < 5%.
gain. However, no values of gain will yield performance that is close to Ts = 4 and OS = 5%, unless you relax the requirements to Ts = 4 and OS < 5%.A true dominant and residue poles approach would resemble the following. This method is effective because the 8-parameter pole/zero compensator can uniquely shape the response of the 8-parameter closed-loop transfer function.
%% Transfer function with desired Dominant poles (Ts = 4, OS = 5%)
zeta = 0.690106730559822;
wn = 1.49884262536285;
Gdom = tf(wn^2, [1, 2*zeta*wn, wn^2])
disp('Dominant poles'), pdom = pole(Gdom)
disp('Check Settling Time and Overshoot of Gdom'), Str = stepinfo(Gdom)
[ydom, tdom] = step(Gdom, linspace(0, 12, 61));
%% Unstable 3rd-order Plant
s = tf('s');
Gp = (s + 1)/(s^3 + 6*s^2 + 10*s - 15)
%% Stable 1st-order Sensor
H = 10/(s + 10)
%% Compensator, Gc
cz = [-9.67806208684678 + 8.60645555609334i % compensator zeros
-9.67806208684678 - 8.60645555609334i
-1.08237879358106 + 0i];
cp = [-460.919397379058 + 0i % compensator poles
-149.13100642372 + 325.266512508622i
-149.13100642372 - 325.266512508622i
143.547728200419 + 0i];
ck = 999283439.026759; % compensator gain
Gc = tf(zpk(cz, cp, ck))
%% Pre-filter, Gf (for cleaning up the unwanted zeros in Gcl)
fz = []; % No pre-filter zeros
fp = [-9.67806208684677 + 8.60645555609335i % Pre-filter poles
-9.67806208684677 - 8.60645555609335i
-10 + 0i
-1.08237879358104 + 0i
-1.00000000000002 + 0i];
fk = 3087.1779210446; % Pre-filter gain
Gf = tf(zpk(fz, fp, fk))
%% Closed-loop TF (with unwanted zeros)
Gcl = feedback(Gc*Gp, H)
%% Filtered Closed-loop
Fcl = minreal(series(Gf, Gcl))
disp('Check if Residue poles are approximately 100 times faster than Dominant poles'), pcl = pole(Fcl)
disp('Check Settling Time and Overshoot of Fcl'), stepinfo(Fcl)
%% Plot results
figure
hold on
plot(tdom, ydom, 'ro')
step(Fcl, 12), grid on
hold off
legend('Dominant TF', 'Filtered CLTF', 'location', 'east')
Categorías
Más información sobre Stability Analysis en Centro de ayuda y File Exchange.
el 3 de Mayo de 2025
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
