Unit 5.2: More on the Qualitative and Quantitative Response of First- and Second-Order Poles#
Follow along at cpjobling.github.io/eg-150-textbook/poles_and_zeros/second_order_properties
Acknowledgements#
The notes for this unit have been influenced by the MATLAB LiveScript PoleZeroAnalysis.mlx from the MathWorks curriculum module Transfer Function Analysis of Dynamic Systems [Allie, 2024]. You will find some useful interactive applications with which you can explore the concepts presented here. Some of these will be deomostrated in the lecture.
The examples have been inspired by the recording of the second half of the video recording of Stephen Boyd’s 1999 lecture 10 (delivered at Stanford University) which itself is based on Qualitative properties of signals & Laplace transforms [Boyd, 1993].
You can watch the video on YouTube here:
Agenda#
Quantitative behaviour of real poles#
In Qualitative properties of terms we stated that a response term that corresponds to a real pole \(s = \sigma\) will result in a growing expontial if \(\sigma > 0\) and a decaying exponential if \(\sigma < 0\).
We didn’t quantify these effects, so it is useful to do so here.
Rate of exponential growth#
Given that the rational polynomial \(F(s)\) representing the Laplace transfer of a signal or system has a pole \(s = \sigma\), one of the terms in the partial fraction expansion will be
where \(r_\sigma\) is the residue of the pole.
The corresponding term in \(f(t)\) will be \(r_\sigma e^{\sigma t}\).
This will be a growing exponential for which the doubling time is a useful measure of the growth rate.
Given that at time \(t=0\), \(r_\sigma e^{0} = r_\sigma\), the doubling time \(T\) is the time for which
Example 7#
Plot the response of the pole
and confirm that the doubling time is around 7 seconds.
Solution to example 7#
Done in MATLAB
format compact
clear variables
t = linspace(0,15,100);
sigma = 0.1;
% Doubling time
T = log(2)/sigma;
% Plot
plot(t,10*exp(sigma*t)),ylim([0,50]),grid
% Plot and label lines: first doubling
line([T,T],[0,20],'Color','r','LineStyle','--')
line([0,T],[20,20],'Color','r','LineStyle','--')
hold on
plot(T,20,'o')
text(3,22,'First doubling')
% Plot and label lines: second doubling
plot(2*T,40,'o')
line([2*T,2*T],[0,40],'Color','g','LineStyle','--')
line([0,2*T],[40,40],'Color','g','LineStyle','--')
text(10,42,'Second doubling')
% Label graph
title('Exponential growth'),
ylabel('f(t) = 10*exp(-0.1*t)'),xlabel('Time t [s]')
hold off
Given that \(\sigma = 0.1\), the doubling time \(T \approx 0.7/\sigma = 7\) s. The initial value is \(10\) at \(t=0\) s. It has doubled to \(20\) at \(t\approx 7\) s, and has doubled again to \(40\) at \(t\approx 14\).
We have added the exact solutions to the plot.
The MATLAB code to reproduce this result is given in example7.mlx
Rate of exponential decay#
If \(F(s)\) has a pole \(s = -\sigma\), one of the terms in the partial fraction expansion will be
where \(r_\sigma\) is the residue of the pole.
The corresponding term in \(f(t)\) will be \(r_\sigma e^{-\sigma t}\).
This will be a decaying exponential for which the time constant \(\tau = 1/|\sigma|\) is a useful measure of the decay rate.
Putting the time constant into the response equation
Thus the response reaches around \(37\%\) of it’s initial value in \(\tau = 1/|\sigma|\) s.
Another measure of decay that is sometimes used is the value of \(t\) for which the response reaches \(1\%\) of its initial value.
Example 8#
The current in the RC circuit shown in Fig. 41 has the transform
If \(R = 1\) M\(\Omega\) and \(C = 10\) \(\mu\)F, determine:
(a) The initial current \(i(0)\).
(b) The time constant \(\tau\).
(c) The time at which the current decays to 1% of \(i(0)\).
Solution to example 8#
(a) The initial current is given by1
Given the component values of \(R\) and \(C\), \(i(0) = 1/RC = 100\) mA.
(b) The time constant \(\tau = RC = 10\) s.
(c) The time at which the current decays to 1% of \(i(0)\) is 46 s.
MATLAB confirmation
t = linspace(0,50,100);
R = 1e6; C = 10e-6;
sigma = -1/(R*C);
% Time constant and initial current
tau = 1/abs(sigma); i0 = 1/(R*C);
% Plot
plot(t,i0*exp(sigma*t)),grid
% Plot and label lines: time constant
line([tau,tau],[0,i0*exp(-1)],'Color','r','LineStyle','--')
line([0,tau],[i0*exp(-1),i0*exp(-1)],'Color','r','LineStyle','--')
text(tau+1,i0*exp(-1),'Time constant: 37% of i0 in 10 seconds')
hold on
plot(tau,i0*exp(-1),'o')
% Plot and label lines: 1% point
line([4.6*tau,4.6*tau],[0,0.01*i0],'Color','g','LineStyle','--')
line([0,4.6*tau],[0.01*i0,0.01*i0],'Color','g','LineStyle','--')
plot(4.6*tau,0.01*i0,'o')
text(32,0.01*i0+0.005,'1% of initial current in 46 seconds')
% Label graph
title('Exponential decay of current in an RC circuit'),
ylabel('Current i(t) [A]'),xlabel('Time t [s]')
hold off
The MATLAB code to reproduce this result is given in example8.mlx
Natural frequency#
The locations of the poles in the \(s\)-plane determine the natural oscillation frequencies present in the system.
Consider a pole in the complex plane illustrated in Fig. 79
Fig. 79 Phase \(\phi\) and magnitude \(r\) of pole \(p\) shown in the complex plane#
The natural (undamped) frequency (in rad/s) contributed by a pole \(p\) is the pole magnitude, \(r=\left|p\right|\).
The phase, \(\phi=\arg p\), determines the relative strength of the oscillatory component.
Analysis of a second-order system#
Second-order systems with a pair of conjugate poles and no zeros are useful for illustrating the natural frequency response analytically. It is helpful to write this type of system in the form:
where
\(K\) is the gain coefficient
\(\zeta\) is the damping ratio already introduced in Complex poles: Damping ratio \zeta and quality factor Q
and \(\omega_n\) is the (undamped) natural frequency.
These parameters are helpful for drawing a connection between the natural response of the system and the properties of the poles. In this section, the poles of the second-order system will be analytically computed and related to the natural response of the system.
Pole analysis#
The poles of the system occur when the denominator is zero:
For simplicity, consider \(\zeta \in [0,1)\). In this case, the poles are
As stated above, the natural frequency is defined as the magnitude of the pole. The magnitude is computed as :
Using the definitions illustrated in Fig. 76, \(\mathbf{Re} = -\sigma = -\zeta\omega_n\) and \(\mathbf{Im} = \omega = \omega_n \sqrt{1 - \zeta^2}\), so :
This illustrates why this form of the equation is used: the natural frequency is the parameter \(\omega_n\).
The damping ratio \(\zeta\) determines the relative strength of the exponential part of the response.
As \(\zeta\to 1\), the complex part of the pole tends to zero, implying less oscillatory and stronger exponential behavior.
For stable systems, that implies greater damping.
The gain coefficient \(K\) affects the magnitude of the response, not the time-dependent behavior.
Natural response#
The impulse response of a system \(H(s)\) is also called the natural response. The natural response \(y(t)\) is computed by taking the inverse Laplace transform of
because the impulse \(x(t)=\delta(t)\) has a Laplace transform \(1\). You can find the response’s analytic form by referring to a table or using the MATLAB ilaplace function.
syms s t zeta omega_n K Y(s) y(t)
assume(t > 0)
Y(s) = K/(s^2 + 2*zeta*omega_n*s + omega_n^2)
Y(s) =
K/(omega_n^2 + 2*zeta*omega_n*s + s^2)
y(t) = ilaplace(Y); % The impulse response in the time domain
Gives the result
Notice that the impulse response is a product of an exponential and sine function. From these functions observe that:
The sine function frequency (in rads/s) is \(\omega_n \,\sqrt{1-\zeta^2 }\). If the damping \(\zeta=0\), then the frequency is \(\omega_n\). This is the reason \(\omega_n\) is often referred to as the undamped natural frequency.
The rate of decay of the exponential damping function is \(\zeta\omega_n\).
The gain parameter \(K\) only contributes as a constant multiplier to the magnitude of the response.
Step-response of a second-order system#
The step response of a system \(H(s)\) is often computed. The step response \(y_s(t)\) is determined by taking the inverse Laplace transform of
because the step function \(x(t)=u_0(t)\) has a Laplace transform \(1/s\). You can find the response’s analytic form by taking partial fraction expansion of \(Y(s)\) and then referring to a table of Laplace transforms.
Here we will state, without proof, that the step response will be
where :
\(\omega = \omega_n\sqrt{1 - \zeta^2}\) is the imaginary part of the system’s complex pole pair, sometimes called the damped natural frequency
\(\sigma = \zeta\omega_n\) is the magnitide of the real part of the system’s complex pole pair
\(\phi = \tan^{-1} \zeta/(1 - \zeta^2)\) is the phase shift.
Note
The proof requires us to complete the square in the denominator of the term with the complex roots
Then take the PFE of the terms assuming that the complex poles will yield terms:
Making the substitutions \(\sigma = \zeta\omega_n\) and \(\omega = \omega_n \sqrt{1 - \zeta^2}\) we can simplify (62) to
The step-response terms will then be
Having computed the residues \(r_1\) and \(r_2\) the sine and cosine terms can be combined to give the final \(e^{-\sigma t}\cos(\omega t + \phi)\) expression.
Example 9#
Determine the step response for a second order system with \(\omega_n = 10\) rad/s and \(\zeta = 0.5\). Confirm your result using (61) and the tf and step functions.
Solution to example 9#
Analytical solution using the symbolic math toolbox
syms Y_s(s) y_s(t)
zeta = 0.5; omega_n = 10;
X(s) = 1/s;
H(s) = omega_n^2/(s^2 + 2*zeta*omega_n*s + omega_n^2);
Y_s(s) = H(s)*X(s);
% step response
y_s(t) = ilaplace(Y_s(s));
Gives the result
Plot the solution
fplot(y_s(t)*heaviside(t),[0,1.2]),ylim([0,1.2]),grid
xlabel('Time (seconds)'),ylabel('y_s(t)'),...
title('Step Response: Determined analytically')
Confirming the result with (61) we get:
sigma = zeta*omega_n;
omega = omega_n*sqrt(1 - zeta^2);
% Compute response using the formula
t = linspace(0,1.2,100);
yst = (1 - exp(-sigma*t).*(cos(omega*t)+(sigma/omega)*sin(omega*t)));
Plot the result
plot(t,yst),grid,...
xlabel('Time (seconds)'),ylabel('y_s(t)'),...
title('Step Response: Computed with formula')
Using the tf and step functions we get:
num = omega_n^2; den = [1 2*zeta*omega_n,omega_n^2];
Hs = tf(num,den)
Hs =
100
----------------
s^2 + 10 s + 100
Continuous-time transfer function.
step(Hs),title('Step Response: Computed with tf function')
The MATLAB code to reproduce this result is given in example9.mlx
Useful quantifiers that result from step response#
The relationship between the pole locations and the step response are summarized in Fig. 80 2
Fig. 80 Relationship between poles and step response#
Knowledge of the location of the poles and their damping ratio and natural frequency are useful for qualtitaive anaysis of the response of a complex pole pair. There are also a number of parameters that quantively define the features of the step response which we may be interested in. These are shown in the right-hand picture of Fig. 81. Knowledge of the poles allows us to predict the step response. But also, usefully for many design problems, knowledge of the step response allows us to predict the location of the poles.
We will explore these ideas in the following sections.
Poles#
The location of the system poles is important and can be obtained by factorizing the denominator of \(H(s)\) either symbolically or numerically. There is also a handy function pole which will take these values from a transfer function.
We will illustrate these with the example used in Example 9.
First we set up the system function \(H(s)\)
syms s t
zeta = 0.5; omega_n = 10;
H = omega_n^2/(s^2 + 2*zeta*omega_n*s + omega_n^2)
H =
100/(s^2 + 10*s + 100)
Find the poles symbolically
[num,den] = numden(H);
sym_poles = factor(den,'FactorMode','full') % FactorMode needed to reduce quadratic
sym_poles =
[s + 5 + 3^(1/2)*5i, s + 5 - 3^(1/2)*5i]
Factors interpreted as
Now find the poles numerically
% Convert symbolic polynomials to numeric polynomials
n = sym2poly(num); d = sym2poly(den);
% Compute the roots
den_poles = roots(d)
den_poles =
-5.0000 + 8.6603i
-5.0000 - 8.6603i
FInally from the transfer function model
Hs = tf(n,d);
poles = pole(Hs)
poles =
-5.0000 + 8.6603i
-5.0000 - 8.6603i
Damping ratio and natural frequencies#
These are most conveniently obtained from the transfer function using the damp function.
[wn,z] = damp(Hs)
wn =
10.0000
10.0000
z =
0.5000
0.5000
The function damp can also return the poles
[wn,z,p] = damp(Hs)
wn =
10.0000
10.0000
z =
0.5000
0.5000
p =
-5.0000 + 8.6603i
-5.0000 - 8.6603i
Rise-tme#
The rise time \(T_r\) is a measure of the speed of response of a system. It is usually taken to be the time taken to transition from 10% to 90% of the final value in the initial rise of the response. It is shown in Fig. 81.
The rise-time depends on \(\omega_n\) but its actual value is also dependent on the damping ratio \(\zeta\) so we rely on a calibration curve such as that shown in Fig. 81.
Fig. 81 Damping ratio as a function of normalised rise-time for a second-order underdamped response#
For the problem being considered \(\zeta = 0.5\), so \(\omega_n T_r \approx 1.65\) giving
step(Hs),line([0,0.1],[0.1,0.1]),line([0,0.25],[0.9,0.9])
Settling time#
The settling time \(T_s\) is defined as the time taken for the peaks of the oscillations in the step response to be bounded by some arbitrary limit. In Fig. 81 the bounds have been set to 2%.
The actual setting time is related to the real part of the poles and, for this case,
For our example
Ts = 4/(zeta*omega_n) % seconds
Ts =
0.8000
step(Hs),line([0,1.2],[1.02,1.02]),line([0,1.2],[0.98,0.98])
The value looks about right!
Peak overshoot#
The peak overshoot is a measure of damping in a system and is the height of the first peak (\(C_\max\) in Fig. 81). It is usually quoted as a percentage of the final value.
If we know the peak overshoot, we can calculate the damping ratio \(\zeta\) using the formula:
For our example
POS = exp(-zeta*pi/sqrt(1 - zeta^2))*100
POS =
16.3034
step(Hs),line([0,1.2],[1.16304,1.16304])
Reverse formula check
z = -log(POS/100)/sqrt(pi^2 + log(POS/100)^2)
z =
0.5000
Peak time#
We sometimes compute the peak time \(T_p\), which is the time at which the first peak occurs. It depends on the damped natural frequency:
For our example
Tp = pi/(omega_n*sqrt(1 - zeta^2))
Tp =
0.3628
step(Hs),line([Tp,Tp],[0,1.16304])
List all properties#
MATLAB provides a useful function stepinfo that computes a step response for a system and takes measurements of the response to summarize the useful quantitative data defined above as well as a few more quantities.
stepinfo(Hs)
ans =
struct with fields:
RiseTime: 0.1639
TransientTime: 0.8076
SettlingTime: 0.8076
SettlingMin: 0.9315
SettlingMax: 1.1629
Overshoot: 16.2929
Undershoot: 0
Peak: 1.1629
PeakTime: 0.3592
You should compare the values computed with the approximations given in this section.
Exercises 16#
Exercise 16.1: Qualitative analysis#
A signal or system response \(f(t)\) contains the terms
(a) How many poles are there in \(F(s)\)?
(b) What is the nature of the response due to the terms given in (66)?
(c) Which are the dominant poles?
Exercise 16.2: Qualitative and quantitative analysis#
A system has poles \(s = 0,\,1\pm j,\,-1,\, -2\pm j2,\,-3\) and zeros at \(s = -1.5,\,-3\pm j3\).
(a) Give the transfer function \(F(s)=b(s)/a(a)\)
(b) Plot the poles on a pole zero map
(c) Give the time constant \(\tau\), doubling time \(T\), damping ratio \(\zeta\), angle \(\theta\), natural frequency \(\omega_n\), and quality factor \(Q\) as appropriate for each pole or pole pair.
(d) How do the zeros affect the response \(f(t)\)?
(e) Which is the dominant pole?
(f) Will the system be stable or unstable?
Partial solution to Exercise 16.2#
The poles and zeros can be entered into MATLAB
z = [-1.5; -3 - 3j; -3 + 3j];
p = [0; 1 + j; 1 - j; -1; -2 + 2j; -2 - 2j; -3];
% Define a transfer function in zero-pole-gain form
Fs = zpk(z,p,1)
Fs =
(s+1.5) (s^2 + 6s + 18)
-------------------------------------------
s (s+1) (s+3) (s^2 - 2s + 2) (s^2 + 4s + 8)
Continuous-time zero/pole/gain model.
(a) The transfer function is
tf(Fs)
ans =
s^3 + 7.5 s^2 + 27 s + 27
-----------------------------------------------------
s^7 + 6 s^6 + 13 s^5 + 6 s^4 - 10 s^3 + 40 s^2 + 48 s
Continuous-time transfer function.
(b) The pole-zero map
pzmap(Fs),xlim([-4,1]),ylim([-4,4]),title('Pole zero map for Example 16.2')
Bonus: the step response is
step(Fs)
Exercise 16.3: Spring-Mass-Damper System#
You can analyze a mass-spring-damper system (Fig. 82) by looking at the poles of its transfer function.
Fig. 82 Diagram of a spring-mass-damper system#
Consider a mass-spring-damper that
is dynamically forced by an arbitrary function
has zero initial conditions: \(x(0)=0\) and \(x'(0)=0\)
The position transfer function is
(a) Write the mass-spring-damper transfer function in the form
and solve for the expressions of the gain \(K\), damping ratio \(\zeta\), and natural frequency \(\omega_n\) in terms of the mass-spring-damper parameters. Write your answers using the symbolic variables \(c\), \(m\), and \(k\).
% Use these symbolic variables
syms m c k
% Replace the NaNs with your expressions
K = NaN;
omega_n = NaN;
zeta = NaN;
(b) Solve for the symbolic expressions of the poles of \(G\) in terms of the mass-spring-damper parameters \(m\), \(c\), and \(k\). Store the expressions below in pplus and pminus where pplus stores the positive root.
(c) Plot the step response of the system \(G\) starting with \(k = m = 1\) and \(c=0\). Note the values of the poles, damping ratio, and natural frequency obtained. Observe the step-response parameters obtained using sysinfo. Adjust the values of \(k\), \(m\) and \(c\) and comment on the effects on the step response observed.
Summary#
Unit 5.2: Take aways#
Real pole \(s=\sigma\):#
Growth rate: \(\sigma > 0\), the exponential signal doubles every \(T\approx 0.7/\sigma\) s.
Time constant: \(\sigma < 0\), the exponential signal decays to \(37\%\) of it’s original value in \(\tau = 1/\sigma\) s. The exponential signal reaches \(1\%\) of it’s original value in \(t\approx 4.6\tau\) s.
Complex pole pair \(s = -\sigma \pm j\omega\):#
Representation of denominator: \((s + \sigma + j\omega_n)(s - \sigma + j\omega) = s^2 + 2\sigma s + \left(\sigma^2 + \omega_n^2\right)\)
Natural frequency: \(\omega_n = \sqrt{\sigma^2 + \omega^2}\)
Damping ratio: \(\zeta = -\sigma/\sqrt{\sigma^2 + \omega^2}\)
Standard second-order system transfer function#
Transfer function:
Poles: \(p_{1,2} = -\zeta\omega_n \pm j\omega_n\sqrt{1 - \zeta^2}\)
Performance parameters for second-order system#
Rise time: \(T_r\) is estimated from graph of normaized rise-time \(\omega_n T_r\) v \(\zeta\) shown in Fig. 81.
2% Setting time: \(T_s \approx 4/\sigma\)
Percentage overshoot: \(\%\mathrm{OS} = \exp\left(-\zeta\pi/\sqrt{1 - \zeta^2}\right) \times 100\)
Peak time: \(T_p = \pi/\omega\)
Useful MATLAB commands#
tf: defines a system as a transfer function \(n(s)/d(s)\)step: plots the step response of a systempole: lists the poles of a systemdamp: lists the natural frequency and damping ratios of the poles of a systemstepinfo: lists the properties of the step response of a system
Unit 5.2: Further exploration#
If you have access to MATLAB (desktop of online), you can install the MATLAB curriculum module Transfer Function Analysis of Dynamic Systems [Allie, 2024] and further explore the concepts given here. The MATLAB Live Script PoleZeroAnalysis.mlx provides some interactive tools that enable you to interactively explore second order system response and includes the impact of zeros which we have not covered here.
References#
- All24(1,2)
Chad Allie. Transfer function analysis of dynamic systems. 2024. Retrieved April 3, 2024. URL: MathWorks-Teaching-Resources/Transfer-Function-Analysis-of-Dynamic-Systems.
- Boy93
Stephen Boyd. Ee-102: introduction to signals and systems. 1993. Retrieved April 3, 2024. URL: https://web.stanford.edu/~boyd/ee102/.
Footnotes#
- 1
This is easily determined from \(\left.f(t)\right|_{t\to 0} = \left.\frac{1}{RC} e^{-t/RC}\right|_{t\to 0} = 1/RC\).
- 2
In the diagram we use a slightly changed notation: \(\omega_d = \omega\) is the damped natural frequency given by the imaginary part of the pole and \(\sigma_d = \sigma\) is the real part of the pole.