I haven't defined syms x but solve function gives me values depend on x

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)

G = 

H= 10/(s+10)

H = 

F= (Kd*s^2 + Kp*s + Ki)/s

F = 

Tcl= (F*G)/(1 + F*G*H)

Tcl = 

[num, den] = numden(simplifyFraction(Tcl));
num = expand(num)

num = 

den = expand(den)

den = 

num_s = coeffs(num, s, 'All')

num_s = 

den_s = coeffs(den, s, 'All')

den_s = 

Pd=s^2 +2*zeta*w_n*s+(w_n)^2

Pd = 

syms e0 e1 e2 positive;

Pr= (s+e0)*(s+e1)*(s+e2)

Pr = 

Pe=expand(Pd*Pr)

Pe = 

Pe_s= coeffs(Pe,s,"All")

Pe_s = 

eq2 = Pe_s == den_s

eq2 = 

sol2 = solve(eq2, [Kp, Kd, e0, e1, e2], 'Real', true);

Warning: Solutions are parameterized by the symbols: x. To include parameters and conditions in the solution, specify the 'ReturnConditions' value as 'true'.

Warning: Solutions are only valid under certain conditions. To include parameters and conditions in the solution, specify the 'ReturnConditions' value as 'true'.

% den_s=subs(den_s,[Kd,Kp],[sol2.Kd,sol2.Kp])
% F = subs(F,[Kd,Kp],[sol2.Kd,sol2.Kp])

Answer 1

[Walter Roberson]

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

[Cahit Semih]

As far as I understand, here we accept e2 as a free variable, but I want to write all variables only depending on Ki.

[Walter Roberson]

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

eq2 = 

sol2 = solve(eq2, [Kp, Kd, e0, e1, e2], 'real', false, 'maxdegree', 3)

sol2 = struct with fields:
    Kp: [6x1 sym]
    Kd: [6x1 sym]
    e0: [6x1 sym]
    e1: [6x1 sym]
    e2: [6x1 sym]

result = subs([Kp, Kd, e0, e1, e2], sol2)

result = 

●

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

[Sam Chak]

Hi @Walter Roberson

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⁴.

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)

Tcl = 

[~, den]    = numden(Tcl);
[den, term] = coeffs(den, s, 'All')     % polynomial degree 5

den = 

term = 

%% 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 = 

[Pe, term] = coeffs(Pe, s, "All")

Pe = 

term = 

%% PART 3: Identify simultaneous equations before solving them
% --------------------------------
eq2 = Pe == den;

Eq1 = eq2(1)

Eq1 = 

Eq2 = eq2(2)

Eq2 = 

Eq3 = eq2(3)

Eq3 = 

Eq4 = eq2(4)

Eq4 = 

Eq5 = eq2(5)

Eq5 = 

Eq6 = eq2(6)

Eq6 = 

%% 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 Matrix:

% Hurwitz characteristic polynomial for degree 5
p5  = double([HurwitzMatrix(6, :)])

p5 = 1×6

1 5 10 10 5 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

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)

C =
 
             1          
  Kp + Ki * --- + Kd * s
             s          

  with Kp = 33.9, Ki = 33.6, Kd = 3.24
 
Continuous-time PID controller in parallel form.

s   = tf('s');
G   = (s + 1)/(s^3 + 6*s^2 + 10*s - 15);
H   = 10/(s + 10);
Tcl = feedback(series(C, G), H)

Tcl =
 
    3.24 s^4 + 69.51 s^3 + 438.6 s^2 + 707.8 s + 335.5
  ------------------------------------------------------
  s^5 + 16 s^4 + 102.4 s^3 + 456.1 s^2 + 524.3 s + 335.5
 
Continuous-time transfer function.

pole(Tcl)

ans = 
  -9.7731 + 0.0000i
  -2.4495 + 5.2943i
  -2.4495 - 5.2943i
  -0.6640 + 0.7537i
  -0.6640 - 0.7537i

figure
step(Tcl)
grid on

stepinfo(Tcl)

ans = struct with fields:
         RiseTime: 0.2011
    TransientTime: 3.9072
     SettlingTime: 3.9072
      SettlingMin: 0.9818
      SettlingMax: 1.4203
        Overshoot: 42.0275
       Undershoot: 0
             Peak: 1.4203
         PeakTime: 0.5264

[Cahit Semih]

The approach I am aiming for is this: Obtaining other variables dependent on the coefficient of Ki and observing how changes in Ki affect them. I did not specify it when I asked the question, but in the sequel I used the Routh-Hurwitz stability criteria to determine the appropriate range for Ki depending on Ki. I then designed the ideal PID controller that would dominate the desired design criteria by minimizing the roots of the residue polynomial that depends on the residue Ki.

[Sam Chak]

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 .

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.

%% Dominant TF (2nd-order)
zeta= 0.690106730559822;
wn  = 1.49884262536285;
G2  = tf(wn^2, [1, 2*zeta*wn, wn^2])

G2 =
 
          2.247
  ---------------------
  s^2 + 2.069 s + 2.247
 
Continuous-time transfer function.

[~, dG2] = tfdata(G2, 'v')

dG2 = 1×3

1.0000 2.0687 2.2465

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

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))

G3 =
 
                100.1
  ---------------------------------
  s^3 + 13.93 s^2 + 64.69 s + 100.1
 
Continuous-time transfer function.

%% Expanded TF (5th-order Closed-loop)
G5  = series(G3, G2)

G5 =
 
                          225
  ----------------------------------------------------
  s^5 + 16 s^4 + 95.76 s^3 + 265.3 s^2 + 352.5 s + 225
 
Continuous-time transfer function.

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');

Answer 2

[Sam Chak]

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)

Tcl = 

[~, den]    = numden(Tcl);
[den, term] = coeffs(den, s, 'All')     % polynomial degree 5

den = 

term = 

%% 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")

Pe = 

term = 

%% 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

Eq1 = 

Eq2 = eq2(2)        % constraint, but not useful, so we choose not to comply

Eq2 = 

Eq3 = eq2(3)        % 1st equation

Eq3 = 

Eq4 = eq2(4)        % 2nd equation

Eq4 = 

Eq5 = eq2(5)        % 3rd equation

Eq5 = 

Eq6 = eq2(6)        % 4th equation

Eq6 = 

Eqs = [Eq3, Eq4, Eq5, Eq6];

% Solve 4 equations for 4 unknowns
sol = solve(Eqs, [K1, K2, K3, e], 'Real', true, 'ReturnConditions', true)

sol = struct with fields:
            K1: 86.5771
            K2: 65.6520
            K3: 20.2595
             e: 1.6131
    parameters: [1x0 sym]
    conditions: symtrue

%% 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)

C =
 
             1          
  Kp + Ki * --- + Kd * s
             s          

  with Kp = 65.7, Ki = 86.6, Kd = 20.3
 
Continuous-time PID controller in parallel form.

s   = tf('s');

% the Plant
G   = (s + 1)/(s^3 + 6*s^2 + 10*s - 15)

G =
 
           s + 1
  -----------------------
  s^3 + 6 s^2 + 10 s - 15
 
Continuous-time transfer function.

% the Sensor
H   = 10/(s + 10)

H =
 
    10
  ------
  s + 10
 
Continuous-time transfer function.

% Closed-loop TF (too many zeros affect the overshoot, see if they are cancellable)
Tcl = feedback(series(C, G), H)

Tcl =
 
    20.26 s^4 + 288.5 s^3 + 1011 s^2 + 1609 s + 865.8
  -----------------------------------------------------
  s^5 + 16 s^4 + 272.6 s^3 + 944.1 s^2 + 1372 s + 865.8
 
Continuous-time transfer function.

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)

Gf =
 
                        865.8
  -------------------------------------------------
  20.26 s^4 + 288.5 s^3 + 1011 s^2 + 1609 s + 865.8
 
Continuous-time transfer function.

% Filtered closed-loop (no more disturbing zeros, very 'CLEAN' now)
Fcl = minreal(series(Gf, Tcl))

Fcl =
 
                          865.8
  -----------------------------------------------------
  s^5 + 16 s^4 + 272.6 s^3 + 944.1 s^2 + 1372 s + 865.8
 
Continuous-time transfer function.

S2  = stepinfo(Fcl)

S2 = struct with fields:
         RiseTime: 1.8913
    TransientTime: 4.0858
     SettlingTime: 4.0858
      SettlingMin: 0.9077
      SettlingMax: 1.0200
        Overshoot: 2.0000
       Undershoot: 0
             Peak: 1.0200
         PeakTime: 4.0855

%% 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

[Cahit Semih]

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))

Kp(Ki):

ans = 

disp('Kd(Ki):'); vpa(simplify(sol.kd))

Kd(Ki):

ans = 

%% 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);

En iyi Ki: 31.6228, Kp: 13.4828, Kd: -3.9971, m_kutup: 0.0813

% 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');

[Sam Chak]

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.

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.

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

Pd = 

% 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)

den = 

% 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])

sol = struct with fields:
     a: 14
     b: (22517998136852480*ki)/4728215395125409 - 373130563904622547/2476415581440161
     c: (22517998136852480*ki)/4728215395125409
    kp: (4503599627370496*ki)/4728215395125409 - 927782854911531624052135807257523/55763921448941996348406004449280
    kd: (2251799813685248*ki)/4728215395125409 - 1062714778089363148174872296622783/55763921448941996348406004449280

disp('Kp(Ki):'); vpa(simplify(sol.kp))

Kp(Ki):

ans = 

disp('Kd(Ki):'); vpa(simplify(sol.kd))

Kd(Ki):

ans = 

%% 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)

min_Ki = 34.0714

initKi  = log10(min_Ki)

initKi = 1.5324

Use the Table instead.

%% 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

res_roots = 
 -13.3220 + 0.0000i
  -0.3390 + 1.5629i
  -0.3390 - 1.5629i

res_roots = 
  -6.7844 +10.1251i
  -6.7844 -10.1251i
  -0.4312 + 0.0000i

res_roots = 
  -6.8563 +19.2857i
  -6.8563 -19.2857i
  -0.2875 + 0.0000i

res_roots = 
  -6.8775 +29.6143i
  -6.8775 -29.6143i
  -0.2450 + 0.0000i

res_roots = 
  -6.8864 +42.7337i
  -6.8864 -42.7337i
  -0.2272 + 0.0000i

res_roots = 
  -6.8906 +60.0875i
  -6.8906 -60.0875i
  -0.2188 + 0.0000i

res_roots = 
  -6.8927 +83.4576i
  -6.8927 -83.4576i
  -0.2146 + 0.0000i

res_roots = 
   1.0e+02 *

  -0.0689 + 1.1521i
  -0.0689 - 1.1521i
  -0.0021 + 0.0000i

res_roots = 
   1.0e+02 *

  -0.0689 + 1.5852i
  -0.0689 - 1.5852i
  -0.0021 + 0.0000i

res_roots = 
   1.0e+02 *

  -0.0689 + 2.1777i
  -0.0689 - 2.1777i
  -0.0021 + 0.0000i

% 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);

En iyi Ki: 10000.0000, Kp: 9508.3083, Kd: 4743.4156, m_kutup: 0.2107

% Transfer function
Fs_best = tf([best_kd best_kp best_ki], [1 0]);
Ts_best = minreal(feedback(Gs * Fs_best, tf(10, [1 10])))

Ts_best =
 
    4743 s^4 + 6.169e04 s^3 + 1.62e05 s^2 + 2.051e05 s + 1e05
  -------------------------------------------------------------
  s^5 + 16 s^4 + 4.75e04 s^3 + 1.426e05 s^2 + 1.949e05 s + 1e05
 
Continuous-time transfer function.

% check the zeros
zero(Ts_best)

ans = 
 -10.0000 + 0.0000i
  -1.0023 + 1.0505i
  -1.0023 - 1.0505i
  -1.0000 + 0.0000i

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)

Gf =
 
                            1e05
  ---------------------------------------------------------
  4743 s^4 + 6.169e04 s^3 + 1.62e05 s^2 + 2.051e05 s + 1e05
 
Continuous-time transfer function.

% Filtered closed-loop (no more disturbing zeros, very 'CLEAN' now)
Fcl = minreal(series(Gf, Ts_best))

Fcl =
 
                              1e05
  -------------------------------------------------------------
  s^5 + 16 s^4 + 4.75e04 s^3 + 1.426e05 s^2 + 1.949e05 s + 1e05
 
Continuous-time transfer function.

S2  = stepinfo(Fcl)

S2 = struct with fields:
         RiseTime: 2.4211
    TransientTime: 4.1649
     SettlingTime: 4.1649
      SettlingMin: 0.9054
      SettlingMax: 0.9978
        Overshoot: 0
       Undershoot: 0
             Peak: 0.9978
         PeakTime: 6.8617

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]

In the first question, when we found roots for the residue polynomial, the roots contained dependent variables, but in the last code, the coefficients of the residue polynomial were completely dependent on Ki. Was there a mistake in the syntax?

[Sam Chak]

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.

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]

Is there a function in matlab that examines the stability of polynomials containing symbolic expressions, especially for Routh-Hurwitz stability criteria? I obtained the transfer function in terms of ki and calculated it using paper and pencil, it was really troublesome.

[Torsten]

Make a google search:

matlab & Routh Hurwitz stability

There are several codes you could try from the File Exchange.

[Walter Roberson]

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.

[Sam Chak]

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.

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%.

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])

Gdom =
 
          2.247
  ---------------------
  s^2 + 2.069 s + 2.247
 
Continuous-time transfer function.

disp('Dominant poles'), pdom = pole(Gdom)

Dominant poles

pdom = 
  -1.0344 + 1.0847i
  -1.0344 - 1.0847i

disp('Check Settling Time and Overshoot of Gdom'), Str = stepinfo(Gdom)

Check Settling Time and Overshoot of Gdom

Str = struct with fields:
         RiseTime: 1.3992
    TransientTime: 4.0000
     SettlingTime: 4.0000
      SettlingMin: 0.9005
      SettlingMax: 1.0500
        Overshoot: 5.0000
       Undershoot: 0
             Peak: 1.0500
         PeakTime: 2.8939

[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)

Gp =
 
           s + 1
  -----------------------
  s^3 + 6 s^2 + 10 s - 15
 
Continuous-time transfer function.

%% Stable 1st-order Sensor
H   = 10/(s + 10)

H =
 
    10
  ------
  s + 10
 
Continuous-time transfer function.

%% 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))

Gc =
 
   9.993e08 s^3 + 2.042e10 s^2 + 1.886e11 s + 1.814e11
  -----------------------------------------------------
  s^4 + 615.6 s^3 + 1.565e05 s^2 + 2.09e07 s - 8.472e09
 
Continuous-time transfer function.

%% 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))

Gf =
 
                           3087
  ------------------------------------------------------
  s^5 + 31.44 s^4 + 423.5 s^3 + 2461 s^2 + 3884 s + 1816
 
Continuous-time transfer function.

%% Closed-loop TF (with unwanted zeros)
Gcl = feedback(Gc*Gp, H)

Gcl =
 
                   9.993e08 s^5 + 3.142e10 s^4 + 4.232e11 s^3 + 2.46e12 s^2 + 3.881e12 s + 1.814e12
  ------------------------------------------------------------------------------------------------------------------
  s^8 + 631.6 s^7 + 1.665e05 s^6 + 2.345e07 s^5 + 1.867e09 s^4 + 8.016e10 s^3 + 1.498e12 s^2 + 2.977e12 s + 3.085e12
 
Continuous-time transfer function.

%% Filtered Closed-loop
Fcl = minreal(series(Gf, Gcl))

Fcl =
 
                                                       3.085e12
  ------------------------------------------------------------------------------------------------------------------
  s^8 + 631.6 s^7 + 1.665e05 s^6 + 2.345e07 s^5 + 1.867e09 s^4 + 8.016e10 s^3 + 1.498e12 s^2 + 2.977e12 s + 3.085e12
 
Continuous-time transfer function.

disp('Check if Residue poles are approximately 100 times faster than Dominant poles'), pcl = pole(Fcl)

Check if Residue poles are approximately 100 times faster than Dominant poles

pcl = 
   1.0e+02 *

  -1.0569 + 0.0000i
  -1.0531 + 0.0067i
  -1.0531 - 0.0067i
  -1.0454 + 0.0066i
  -1.0454 - 0.0066i
  -1.0415 + 0.0000i
  -0.0105 + 0.0110i
  -0.0105 - 0.0110i

disp('Check Settling Time and Overshoot of Fcl'), stepinfo(Fcl)

Check Settling Time and Overshoot of Fcl

ans = struct with fields:
         RiseTime: 1.3795
    TransientTime: 4.0000
     SettlingTime: 4.0000
      SettlingMin: 0.9002
      SettlingMax: 1.0500
        Overshoot: 4.9968
       Undershoot: 0
             Peak: 1.0500
         PeakTime: 2.9123

%% Plot results
figure
hold on
plot(tdom, ydom, 'ro')
step(Fcl, 12), grid on
hold off
legend('Dominant TF', 'Filtered CLTF', 'location', 'east')