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. 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.
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. HurwitzMatrix = sym(zeros(m + 1, m + 1));
HurwitzMatrix(n + 1, k + 1) = nchoosek(n, k);
disp('Hurwitz Matrix:'); disp(HurwitzMatrix);
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, :)]);
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. k4 = p5(end-3)*wn^2 - 0.02
K = [k1, k2, k3, k4, 0]/2660.4;
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. 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
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. 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.
Gcl = feedback(1/sys.B(end)*Gp, Hc);
Gcl = k1*minreal(Gcl)
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.
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.
step(Gcl, 2*ceil(Ts1)), grid on
legend('Full-state feedback', 'Output feedback', 'location', 'east')