Analyzing and developing PID control for a 5th order system

18 次查看(过去 30 天)
Hi all,
I would like to know if there is a way of simulating and designing a controller for a 5th order system in MATLAB?
I have derived the open loop transfer function of my system and would like to analyse the response and develop a controller, but I have not been successful using the PID Tuner or the Control System Designer, the step responses are also very unstable:
G = tf([2660.4], [1 20.75 0.02 0 0 0])
G =
2660
--------------------------
s^5 + 20.75 s^4 + 0.02 s^3
Is it perhaps because its a 5th order system that Im encountering challenges?
Thanks
Niel

采纳的回答

Sam Chak
Sam Chak 2025-5-31
The characteristic polynomial of the 5th-order plant has zero coefficients for the terms , , and . The configuration of the PID controller can only affect the proportional term () and the derivative term (). The coefficient of the accelerative term () remains zero regardless of how we tune the PID controller. According to the Routh–Hurwitz stability criterion, this implies that such a system cannot be stabilized by the PID controller.
To stabilize the 5th-order plant, we can utilize a 5th-order compensator. After analyzing the locations of the plant's poles in the root locus diagram, you can position the compensator's poles and zeros until a ten-legged camel spider-like locus emerges. Ideally, they should be placed near the desired damping ratio grids and natural frequency contours. If the closed-loop zeros have negative real parts, they can be canceled out using a pre-filter to improve the transient response performance.
%% Plant
Gp = tf([2660.4], [1 20.75 0.02 0 0 0])
Gp = 2660 -------------------------- s^5 + 20.75 s^4 + 0.02 s^3 Continuous-time transfer function.
pp = pole(Gp)
pp = 5×1
0 0 0 -20.7490 -0.0010
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
%% Compensator
cz = [-20.7 + 0i % Gc zeros
-1.14 + 2.13i
-1.14 - 2.13i
-1.88 + 0i];
cp = [-31.3 + 0i % Gc poles
-22.1 + 18.1i
-22.1 - 18.1i
-3.32 + 18.7i
-3.32 - 18.7i];
ck = 22307; % Gc gain
Gc = tf(zpk(cz, cp, ck))
Gc = 22307 s^4 + 5.546e05 s^3 + 2.147e06 s^2 + 4.919e06 s + 5.067e06 ---------------------------------------------------------------- s^5 + 82.14 s^4 + 3062 s^3 + 6.738e04 s^2 + 9.63e05 s + 9.213e06 Continuous-time transfer function.
figure
rlocusplot(Gc*Gp), grid on
axis([-40 10 -30 30])
%% Closed-loop System
Gcl = feedback(Gc*Gp, 1)
Gcl = 5.935e07 s^4 + 1.475e09 s^3 + 5.711e09 s^2 + 1.309e10 s + 1.348e10 -------------------------------------------------------------------------------------------------------------------------------------------- s^10 + 102.9 s^9 + 4766 s^8 + 1.309e05 s^7 + 2.361e06 s^6 + 2.92e07 s^5 + 2.505e08 s^4 + 1.476e09 s^3 + 5.711e09 s^2 + 1.309e10 s + 1.348e10 Continuous-time transfer function.
damp(Gcl)
Pole Damping Frequency Time Constant (rad/seconds) (seconds) -4.50e+00 + 1.49e+00i 9.49e-01 4.74e+00 2.22e-01 -4.50e+00 - 1.49e+00i 9.49e-01 4.74e+00 2.22e-01 -5.27e+00 + 5.54e+00i 6.90e-01 7.65e+00 1.90e-01 -5.27e+00 - 5.54e+00i 6.90e-01 7.65e+00 1.90e-01 -7.82e+00 1.00e+00 7.82e+00 1.28e-01 -1.08e+01 + 9.61e+00i 7.47e-01 1.44e+01 9.27e-02 -1.08e+01 - 9.61e+00i 7.47e-01 1.44e+01 9.27e-02 -1.72e+01 1.00e+00 1.72e+01 5.81e-02 -1.84e+01 + 5.32e+00i 9.61e-01 1.91e+01 5.44e-02 -1.84e+01 - 5.32e+00i 9.61e-01 1.91e+01 5.44e-02
%% Pre-filter
fz = []; % Gf zeros
fp = cz; % Gf poles
fk = real(prod(cz)); % Gf gain
Gf = tf(zpk(fz, fp, fk))
Gf = 227.1 --------------------------------------------- s^4 + 24.86 s^3 + 96.23 s^2 + 220.5 s + 227.1 Continuous-time transfer function.
%% Filtered Closed-loop
Fcl = minreal(series(Gf, Gcl))
Fcl = 1.348e10 -------------------------------------------------------------------------------------------------------------------------------------------- s^10 + 102.9 s^9 + 4766 s^8 + 1.309e05 s^7 + 2.361e06 s^6 + 2.92e07 s^5 + 2.505e08 s^4 + 1.476e09 s^3 + 5.711e09 s^2 + 1.309e10 s + 1.348e10 Continuous-time transfer function.
Scl = stepinfo(Fcl);
figure
step(Fcl, 2*ceil(Scl.SettlingTime)), grid on
  4 个评论
Sam Chak
Sam Chak 2025-6-2
Thank you, @William Rose. The intention was to encourage the OP to explore the Root Locus Editor. Since the performance requirements are not specified, the poles and zeros can be placed arbitrarily, provided that the compensator stabilizes the plant and produces a response without overshoot.
Niel
Niel 2025-6-2
Thanks Sam, I played around a bit with the electrical and mechanical design, adding damping and stiffness into the system, I then used the Root Locus editor to help develop the controller, I think it came out quite good, thank you for the assistance:

请先登录,再进行评论。

更多回答(1 个)

Sam Chak
Sam Chak 2025-6-1
I am not an expert in motor hardware, so I cannot explain how the 5h-order transfer function is derived. Most mechanical motions can be described by 2nd-order differential equations using the Newton–Euler equations. The high order likely arises from the couplings of mechanical elements. Furthermore, adding physical damping and stiffness to the system does not alter the fact that the acceleration term () is still absent.
Previously, I suggested using the root locus technique because some individuals are visual learners. They tend to develop an intuitive understanding by observing how the locations of the poles and zeros affect the system response as they are interactively adjusted on the Root Locus Editor plot.
For those who use logical reasoning to solve problems, understanding governing laws and rules is beneficial. In your example, the terms , , and are absent from the 5th-order transfer function. However, the coefficient of the jerk term is also very small compared to the coefficient of the term . In the sense of control theory, we need a Proportional–Derivative–Acceleration–Jerk controller to compensate for the inadequate terms. This approach is easier to demonstrate using a full-state feedback controller and therefore, the plant's transfer function must first be converted to the state-space model.
%% Plant transfer function
s = tf('s');
Gp = tf(2660.4, [1 20.75 0.02 0 0 0])
Gp = 2660 -------------------------- s^5 + 20.75 s^4 + 0.02 s^3 Continuous-time transfer function.
%% Equivalent State-space model in Canonical form
sys = compreal(Gp);
sys.A = round(sys.A, 4);
sys = ss(sys.A', sys.C', sys.B', sys.D')
sys = A = x1 x2 x3 x4 x5 x1 0 1 0 0 0 x2 0 0 1 0 0 x3 0 0 0 1 0 x4 0 0 0 0 1 x5 0 0 0 -0.02 -20.75 B = u1 x1 0 x2 0 x3 0 x4 0 x5 2660 C = x1 x2 x3 x4 x5 y1 1 0 0 0 0 D = u1 y1 0 Continuous-time state-space model.
Before we can determine the control gains, we need to identify the guiding coefficients of the desired characteristic equation. These coefficients are typically found in certain control papers (though rarely in textbooks) that present the ISE and ITAE tables. One can think of the coefficients as the baking molds that pastry chefs often use to shape cakes and cookies. These coefficients are employed to calculate the gains using specific formulas.
The coefficients in the ISE or ITAE table enable us to design a full-state feedback controller that yields an almost minimum value of the ISE or ITAE cost function for linear systems only. Why "almost"? Because these coefficients are truncated values. Host people find it challenging to memorize these coefficients. So, I'd rather derive them from the Hurwitz-stable matrix. Technically, this matrix relates to Pascal's triangle, or Yang Hui's triangle, where the coefficients in the rows correspond to the number of states in the 2nd column. The coefficients in the Hurwitz matrix ensure exponential stability with a critically damped response. Given that the matrix possesses palindromic properties (reading the same forwards and backward), it is relatively easy to remember these coefficients for low-order systems.
%% Hurwitz Matrix Generation
m = 5; % number of poles (for proper TFs) or the number of states (for state-space)
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:
disp(sprintf('Hurwitz-stable Characteristic Polynomial of degree %d', m)); disp(HurwitzMatrix(m + 1, :));
Hurwitz-stable Characteristic Polynomial of degree 5
p5 = double([HurwitzMatrix(m + 1, :)]); % store coefficients of polynomial as an array
%% Critically-damped Transfer Function with wn = 1
G5 = tf(1, p5)
G5 = 1 --------------------------------------- s^5 + 5 s^4 + 10 s^3 + 10 s^2 + 5 s + 1 Continuous-time transfer function.
Transfer function design template:
From the transfer function design template above, it is easy to memorize the formula for each term, as the sum of the natural frequency's exponent and the Laplace variable's exponent always equals the number of states. To ensure a critically damped response in the closed-loop system by controlling only four terms of a 5th-order system, we need to determine the natural frequency from the term that the controller cannot manipulate, specifically the term. Once ​ is determined, the four control gains can be easily calculated.
%% Full-state feedback control, u = –K·x (only works if all states are measurable)
wn = 20.75/p5(2) % Desired natural frequency
wn = 4.1500
k1 = p5(end-0)*wn^5 % Proportional gain, kp
k1 = 1.2310e+03
k2 = p5(end-1)*wn^4 % Derivative gain, kd
k2 = 1.4831e+03
k3 = p5(end-2)*wn^3 % Acceleration gain, ka
k3 = 714.7338
k4 = p5(end-3)*wn^2 - 0.02 % Jerk gain, kj
k4 = 172.2050
K = [k1, k2, k3, k4, 0]/2660.4; % Feedback gain matrix
By applying the feedback linear state-space control law via , the system's open-loop dynamics are effectively modified to our favor in the form of the closed-loop system dynamics.
%% Closed-loop system (via A – B·K)
Cls = tf(ss(sys.A - sys.B*K, sys.B*K(1), sys.C, sys.D))
Cls = 1231 ------------------------------------------------------- s^5 + 20.75 s^4 + 172.2 s^3 + 714.7 s^2 + 1483 s + 1231 Continuous-time transfer function.
damp(Cls)
Pole Damping Frequency Time Constant (rad/seconds) (seconds) -4.15e+00 + 2.76e-03i 1.00e+00 4.15e+00 2.41e-01 -4.15e+00 - 2.76e-03i 1.00e+00 4.15e+00 2.41e-01 -4.15e+00 + 4.47e-03i 1.00e+00 4.15e+00 2.41e-01 -4.15e+00 - 4.47e-03i 1.00e+00 4.15e+00 2.41e-01 -4.15e+00 1.00e+00 4.15e+00 2.41e-01
S1 = stepinfo(Cls);
Ts1 = S1.SettlingTime;
disp('Settling time of full-state feedback system:'); disp(Ts1);
Settling time of full-state feedback system: 2.5495
In real-world scenarios, while we can measure velocity and acceleration using expensive sensors, such as velocity transducers and accelerometers, it may be challenging to accurately measure jerk. In such situations, we often need to estimate the unknown-but-observable states. The method of estimating these states is known as observer design. This topic is adequately covered in all modern control textbooks.
What I would like to demo is that, instead of using non-realizable full-state feedback, we can employ output feedback control by designing a lowpass filter of the form ​, where the time constant ​ is small enough to produce the same frequency response as the derivative s at the low frequencies of interest. I would like to thank @Paul, from whom I learned this approach three years ago. However, this method suffers from numerical precision issues when dealing with high-order systems that require a fast response (). For a relatively low 5th-order system like yours, we can already observe that the magnitude of some coefficients is on the order of . Therefore, use this approach at your discretion.
Please note that the output feedback controller must be placed in the feedback path. The concept of feedback compensation is discussed in Chapter 9: "Design via Root Locus" in Control Systems Engineering by Norman S. Nise.
%% Output feedback control, Hc
% Hc = k1 + k2*s + k3*s^2 + k4*s^3; % full-state feedback but not realizable (improper form)
wcf = 15*wn; % ensure that cutoff freq >> natural freq
Tf = 1/wcf; % Filter time constant
Hc = minreal(k1 + k2*s/(Tf*s + 1) + k3*(s/(Tf*s + 1))^2 + k4*(s/(Tf*s + 1))^3)
Hc = 4.44e07 s^6 + 8.476e09 s^5 + 5.51e11 s^4 + 1.292e13 s^3 + 4.88e13 s^2 + 9.32e13 s + 7.163e13 -------------------------------------------------------------------------------------------- s^6 + 373.5 s^5 + 5.813e04 s^4 + 4.824e06 s^3 + 2.252e08 s^2 + 5.609e09 s + 5.819e10 Continuous-time transfer function.
%% Closed-loop system
Gcl = feedback(1/sys.B(end)*Gp, Hc); % <-- place controller in the feedback path
Gcl = k1*minreal(Gcl) % multiply with k1 to ensure a zero steady-state error
Gcl = 1231 s^6 + 4.598e05 s^5 + 7.155e07 s^4 + 5.939e09 s^3 + 2.773e11 s^2 + 6.904e12 s + 7.163e13 ---------------------------------------------------------------------------------------------------------------------------------------------------------- s^11 + 394.3 s^10 + 6.588e04 s^9 + 6.031e06 s^8 + 3.254e08 s^7 + 1.033e10 s^6 + 1.83e11 s^5 + 1.758e12 s^4 + 1.292e13 s^3 + 4.88e13 s^2 + 9.32e13 s + 7.163e13 Continuous-time transfer function.
S2 = stepinfo(Gcl);
Ts2 = S2.SettlingTime;
disp('Settling time of output feedback system:'); disp(Ts2);
Settling time of output feedback system: 2.6203
The settling time of the practical output feedback system is only a fraction of a second longer than the settling time of the ideal full-state feedback system. This discrepancy is expected, as we introduced the derivative filter to approximate the non-causal ideal derivative action. If the cutoff frequency is much closer to the natural frequency, you may observe a deterioration in the step response performance.
%% Simulations
figure
hold on
step(Cls, 2*ceil(Ts1))
step(Gcl, 2*ceil(Ts1)), grid on
hold off
ylim([-0.2, 1.2])
legend('Full-state feedback', 'Output feedback', 'location', 'east')
  1 个评论
Sam Chak
Sam Chak 2025-6-2
In Simulink, the control block diagram appears as follows. The controller is placed in the feedback path. Please note that the baseline coefficients {1, 5, 10, 10, 5, 1} influence the overshoot, while the chosen natural frequency affects the settling time.
You may also try using the ITAE coefficients for the 5th-order system, which can be found in Stanley M. Shinners' Modern Control System Theory and Design, 2e. It is important to note that these coefficients are truncated values. Because the ITAE coefficents typically yields faster responses, it is necessary to design a much higher cutoff frequency.
% p5 = [1, 2.8, 5.0, 5.5, 3.4, 1];

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Classical Control Design 的更多信息

产品


版本

R2024b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by