Java程序辅导

ME 349
Dynamic Systems Laboratory
Manual
Boyd D. Schimel
Jow-Lian Ding
Michael J. Anderson
and
Walter J. Grantham
School of Mechanical and Materials Engineering
Washington State University
April 1997
Copyright °c 1997 by B. D. Schimel, J. L. Ding, M. J. Anderson and W. J. Grantham.
All Rights Reserved.
0-2 Schimel, Ding, Anderson and Grantham
Acknowledgment
This manual updates and upgrades a previous ME 349 Lab Manual by J. L. Ding and M. J. Anderson.
Some sections in this manual draw equations and text from the previous manual. The authors would like to
express their thanks for these contributions.
Contents
1 Free Vibration 1-7
1.1 Undamped Free Vibrations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1-7
1.2 Damped Free Vibrations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1-9
1.3 Apparatus, Measurements, and Calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1-12
Measurement of Damped Natural Frequency . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1-12
Linear Regression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1-14
Logarithmic Decrement Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1-15
Measurement of Viscous Damping Coefficient . . . . . . . . . . . . . . . . . . . . . . . . . . . 1-15
1.4 Laboratory #1 Procedures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1-16
Spring Constant . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1-16
Damped Natural Frequency and Damping Ratio . . . . . . . . . . . . . . . . . . . . . . . . . 1-17
Viscous Damping Coefficient . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1-18
2 Forced Vibration 2-20
2.1 Equations of Motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-20
Vibration Generated by a Force Input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-20
Vibration Generated by a Displacement Input . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-20
2.2 Laplace Transforms and the Transfer Function . . . . . . . . . . . . . . . . . . . . . . . . . . 2-22
2.3 Solutions to the Equations of Motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-23
Force Input System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-23
Displacement Input System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-24
2.4 The Frequency Response Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-25
Generalized System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-25
Force Input System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-26
Displacement Input System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-27
2.5 Apparatus, Measurements and Calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-28
Measurement of the Response Amplitude and Frequency . . . . . . . . . . . . . . . . . . . . . 2-28
Measurement of the Input Amplitude . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-29
Force Input System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-29
Displacement Input System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-31
Generation of the FRF Magnitude Plot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-31
Determination of Damping Ratio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-32
Determination of Natural Frequency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-33
2.6 Laboratory #2 Procedures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-33
Lab #2A - Vibrations Generated by Force Input . . . . . . . . . . . . . . . . . . . . . . . . . 2-33
Lab #2B - Vibrations Generated by Displacement Input . . . . . . . . . . . . . . . . . . . . . 2-34
3 Free Vibration of a Two-Degree-of-Freedom System 3-36
3.1 Equations of Motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-36
3.2 Solutions to the Equations of Motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-36
3.3 Analysis of the Free-Vibration Response . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-41
3.4 The Beating Phenomenon . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-43
0-3
0-4 Schimel, Ding, Anderson and Grantham
3.5 Apparatus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-43
Equations of Motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-45
Linearization of the Equations of Motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-45
Solution to the Linearized Equations of Motion . . . . . . . . . . . . . . . . . . . . . . . . . . 3-46
The Beating Phenomenon in the Double-Pendulum System . . . . . . . . . . . . . . . . . . . 3-47
3.6 Measurements and Calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-47
Measurement of the Torsional Spring Constant . . . . . . . . . . . . . . . . . . . . . . . . . . 3-47
Estimation of the Natural Frequencies and Beating Phenomenon Frequencies . . . . . . . . . 3-47
Measurement of the Natural Frequencies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-48
Measurement of the Beating Phenomenon Frequencies . . . . . . . . . . . . . . . . . . . . . . 3-48
3.7 Laboratory #3 Procedures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-49
Physical Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-49
Torsional Spring Constant . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-50
Natural Frequencies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-50
Beating Phenomenon Frequencies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-51
4 Multi-Degree-of-Freedom Vibrations and Spectral Analysis 4-52
4.1 Multi-Degree-of-Freedom Vibrations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-52
Free Vibrations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-52
Equation of Motion, Eigenvalues, and Eigenvectors . . . . . . . . . . . . . . . . . . . . 4-52
The Eigenproblem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-54
Solution to the Equation of Motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-55
Three-Degree-of-Freedom-Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-55
Forced Vibrations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-57
Equation of Motion and the Transfer Function . . . . . . . . . . . . . . . . . . . . . . 4-57
The Forced System Frequency Response Function . . . . . . . . . . . . . . . . . . . . . . . . . 4-58
4.2 Spectral Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-59
Discrete-Time-Measurements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-60
Transducers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-60
Discrete-Time Data Sampling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-61
The Discrete Fourier Transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-62
How the DFT Works . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-62
The Frequency Spectrum and Spectral Leakage . . . . . . . . . . . . . . . . . . . . . . 4-66
Frequency Resolution and the Nyquist Frequency . . . . . . . . . . . . . . . . . . . . . 4-67
Computing the Frequency Response Function from the DFT . . . . . . . . . . . . . . 4-68
The Fast Fourier Transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-70
4.3 Apparatus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-70
The Physical System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-70
Theoretical Structural Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-71
Free Vibrations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-71
Forced Vibrations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-73
4.4 Laboratory #4 Procedures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-75
Free-Vibration Testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-75
Forced-Vibration Testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-75
5 Linear System Identification 5-78
5.1 System Identification and the Transfer Function . . . . . . . . . . . . . . . . . . . . . . . . . . 5-78
5.2 The Frequency Response Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5-78
5.3 The Time Constant . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5-79
First-Order Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5-79
Higher-Order Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5-81
5.4 The Bode Plot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5-82
First-Order Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5-83
Higher-Order Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5-85
ME 349 Dynamic Systems Lab Manual 0-5
Identification of Systems with Closely Spaced Time Constants . . . . . . . . . . . . . . . . . . 5-89
5.5 Apparatus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5-89
Governing Differential Equations and Transfer Functions . . . . . . . . . . . . . . . . . . . . . 5-91
DC Motor and Pump . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5-91
DC Generator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5-91
Voltage to Current Converter and Power Supply . . . . . . . . . . . . . . . . . . . . . 5-92
The System Transfer Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5-93
5.6 Measurements and Calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5-93
Theoretical Transfer Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5-93
Time Constant of the Motor-Pump-Generator system . . . . . . . . . . . . . . . . . . 5-93
Determination of the Soft-Start Time Constant . . . . . . . . . . . . . . . . . . . . . . 5-94
Determination of the Transfer Function Leading Constant . . . . . . . . . . . . . . . . 5-95
5.7 Laboratory #5 Procedures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5-95
Time Constant of the DC Motor, Pump and Generator System . . . . . . . . . . . . . . . . . 5-95
Power Supply Soft-Start Time Constant . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5-96
Transfer Function Leading Coefficient . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5-96
Frequency Response Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5-96
6 Control of Dynamic Systems 6-98
6.1 The Transfer Function and Block Diagram Analysis . . . . . . . . . . . . . . . . . . . . . . . 6-98
6.2 Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6-100
6.3 The Initial and Final Value Theorems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6-101
6.4 Linear Control Laws . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6-101
Proportional Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6-102
Integral Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6-103
Proportional-Plus-Integral Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6-104
Derivative Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6-105
Proportional-Plus-Derivative Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6-105
6.5 PID-Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6-107
6.6 Transient Response Performance Criteria . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6-107
6.7 Discrete-Time Control Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6-108
6.8 Linearization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6-110
6.9 Apparatus, Measurements, and Calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6-112
Governing Differential Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6-112
Linearization of the Governing Differential Equation . . . . . . . . . . . . . . . . . . . . . . . 6-114
Proportional Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6-115
PI-control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6-116
6.10 Laboratory #6 Procedures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6-116
Proportional Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6-116
PI-control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6-117
7 State Space Modeling 7-118
7.1 State-Space Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7-118
Model Format . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7-118
State Transition Matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7-119
Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7-120
Converting State-Space Models to Transfer Function Models . . . . . . . . . . . . . . . . . . . 7-121
Nonlinear State-Space Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7-121
Equilibrium Point Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7-122
Linearization of Nonlinear State-Space Models . . . . . . . . . . . . . . . . . . . . . . . . . . 7-122
7.2 Applications of State-Space Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7-123
Linearization of an Automobile Cruise Control Model . . . . . . . . . . . . . . . . . . . . . . 7-123
A State-Space Model of an Armature Controlled DC Motor . . . . . . . . . . . . . . . . . . . 7-124
7.3 Apparatus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7-125
0-6 Schimel, Ding, Anderson and Grantham
Motor-Pendulum Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7-125
Proportional Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7-126
PD-Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7-128
PID-Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7-129
7.4 Defining Terms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7-130
7.5 Further Information . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7-130
Chapter 1
Free Vibration
1.1 Undamped Free Vibrations
Consider the theoretical system shown in Figure 1-1a. The system is composed of a mass m which is
supported by a spring that has negligible inertia and a stiffness of k. The system is assumed to have a single
degree of freedom of motion along a vertical axis. Let x denote displacements of m, with the downward
direction positive, and let x˙ and x¨ denote velocity and acceleration ofm respectively. Suppose the mass is set
into motion by an impulse or by suddenly releasing it from a point of nonzero displacement. The resulting
motion, in the absence of any additional external inputs, is referred to as free vibration. If there is no source
of energy dissipation for this theoretical system, then the motion is referred to as undamped free vibration;
and once set into motion it will vibrate indefinitely.
Observe the free-body diagram of the system in Figure 1-1b and note that the direction of the spring
force kx is opposite in direction to the displacement x. Application of Newton’s second law (force = mx¨)
gives the following equation of motion:
−kx = mx¨, (1-1)
or after rearrangement
mx¨+ kx = 0. (1-2)
Note that any deflection of the spring due to gravitational force on the mass is not included in the above
equations. Gravitational force has no effect on the dynamic motion of the system; it results only in the
presence of a static displacement, mgk , in the solution x (t) . For the present analysis, the point x = 0 is
assumed to coincide with the point of rest for the system when gravitational displacement is included and
therefore, this displacement does not appear explicitly in Eqs. (1-1) and (1-2).
Equation (1-2) is a linear, homogenous, second-order ordinary differential equation. The solution to this
equation is of the form
x(t) = Ceλt, (1-3)
where C and λ are constants and t is time. Substitution of Eq. (1-3) and its second time-derivative into Eq.
(1-2) gives
(mλ2 + k)Ceλt = 0. (1-4)
The case of C = 0 is a trivial solution. Nontrivial solutions are obtained by solving the equation
mλ2 + k = 0 (1-5)
for λ. Application of the quadratic formula gives
λ = ±i
r
k
m
, (1-6)
where i is equal to
√
−1, and hence, λ is purely imaginary for this system. For the sake of convenience let
ωn =
r
k
m
. (1-7)
1-7
1-8 Schimel, Ding, Anderson and Grantham
Figure 1-1: Theoretical single-degree-of-freedom mass-spring system.
Therefore,
λ = ±iωn. (1-8)
Equation (1-5) is known as the characteristic equation for the system depicted in Figure 1-1a, and the two
possible values of λ in Eq. (1-8) are known as the eigenvalues of the system. As will be seen in the following
development, the presence of purely imaginary eigenvalues will result in a purely oscillatory solution x(t).
Substitution of either the positive or negative value of λ into Eq. (1-3) will result in an equation for
x (t) which satisfies Eq. (1-2). In fact, since Eq. (1-2) is a linear ordinary differential equation, any linear
combination of the of the form
x(t) = C1e
iωnt + C2e
−iωnt (1-9)
will satisfy Eq. (1-2).
Observe that x(t) is a real-valued function which describes the motion of the mass m as a function of
time. Therefore, it is useful to convert the complex components of Eq. (1-9) into a real representation for
x(t). Substitution of Euler’s formula,
e±iθ = cos θ ± i sin θ, (1-10)
into Eq. (1-9) gives
x(t) = (C1 + C2) cosωnt+ i(C1 − C2) sinωnt. (1-11)
For a physically meaningful x(t) Eq. (1-11) must be real-valued. Observe that this equation will be real-
valued only if the constant coefficients (C1 + C2) and i(C1 − C2) are both real-valued. Suitable values for
constants C1 and C2 are therefore given by the complex conjugate pair
C1 = A−iB2
C2 = A+iB2 ,
(1-12)
where A and B are real constants. Substitution of Eqs. (1-12) into (1-11) yields the following equation:
x(t) = A cosωnt+B sinωnt. (1-13)
ME 349 Dynamic Systems Lab Manual 1-9
An alternate form of Eq. (1-13) is
x(t) = ρ cos(ωnt+ φ) (1-14)
where
ρ =
p
A2 +B2 (1-15)
and
φ = − tan−1
µ
B
A
¶
. (1-16)
Note that Eqs. (1-13) and (1-14) real-valued representations of the function x(t). Also observe that these
equations satisfy the equation of motion, Eq. (1-2), in general and without regard to the value or physical
meaning of the constants A and B. This observation can be verified by simple substitution of Eqs. (1-13) or
(1-14) into Eq. (1-2). The physical of meaning these constants, as they apply to the system of Figure 1-1a,
becomes clear when initial conditions are considered. Suppose that the mass is pulled down to the point xo,
and then is suddenly released at time t = 0. The initial conditions for this case are displacement x(0) = xo
and velocity x˙(0) = 0. Substitution of these two initial conditions into Eq. (1-13) at t = 0 gives: A = xo
and B = 0. For this example, the motion of the system is governed by
x(t) = xo cosωnt. (1-17)
As an alternative, suppose that the mass is initially undisplaced, and that an impulse imparts initial velocity
vo to the mass. In this case x(0) = 0 and x˙(0) = vo, constants A and B become 0 and voωn respectively, and
the motion of the system is governed by
x(t) =
vo
ωn
sinωnt. (1-18)
Observe that combinations of the initial conditions from the preceding two examples can also be described
by Eq. (1-13).
Equation (1-13) and the specific cases of Eqs. (1-17) and (1-18) are all sinusoidal in form. The frequency
of the motion defined by these equations is determined by the constant ωn, which is a function of the system
mass and stiffness as defined in Eq.(1-7). The constant ωn has special significance in the study of vibrations,
and is defined as a system’s natural frequency.
Natural frequency has dimensions of angle/time, and units of radians per second are conventionally used.
Note that there is a peculiarity involved in the use of the unit radians. When stating a frequency the units
are stated as radians per second. However, when a frequency is used in a calculation, such as voωn in Eq.
(1-18), the unit radians is dropped or ignored.
Recall that another common unit for frequency is the unit Hertz, and that one Hertz is equivalent one
cycle per second. Hertz and radians per second are related by the following:
1Hz = 2π rad/sec (1-19)
Beware that a common error is made by using frequencies in units of Hz in equations such as (1-13) through
(1-18), which results in incorrect computations. If a frequency is given or measured in Hz, it must be to
converted to radians per second before the frequency is used in calculations involving the above equations.
Another common mistake is made by assuming that computations made with Eq. (1-7) result in units of
Hz. A correct interpretation of the units for this computation given in the following example: If stiffness k
is given in units of N/m, and mass m is given in units of kg, the units on the right-hand-side of Eq. (1-7)
reduce to sec−1, while the resulting natural frequency is read in units of radians per second. Here, the unit
radians is inserted.
1.2 Damped Free Vibrations
Damping is a physical property which accounts for energy dissipation in a system. Viscous damping is a
common model in which the dissipative force is linearly related to velocity:
F = −cx˙. (1-20)
1-10 Schimel, Ding, Anderson and Grantham
  .  
Figure 1-2: Theoretical single-degree-of-freedom mass-spring-dashpot system.
Figure 1-2a shows the theoretical system of Figure 1-1a, with the addition of a viscous damping element, c.
Observation of the free-body diagram in Figure 1-2b and application Newtons’ second law gives the equation
of motion
−cx˙− kx = mx¨. (1-21)
After rearrangement the equation of motion becomes
mx¨+ cx˙+ kx = 0. (1-22)
As in the undamped case, it is assumed that solution x(t) is of the form of Eq. (1-3). For the viscously
damped case the characteristic equation is obtained by substituting Eq. (1-3) into Eq. (1-22), which gives
mλ2 + cλ+ k = 0. (1-23)
The eigenvalues of Eq. (1-23) are
λ =
−c±√c2 − 4mk
2m
. (1-24)
By defining the two eigenvalues in Eq. (1-24) as λ1 and λ2 the general solution can, in most cases, be written
as
x(t) = C1e
λ1t + C2e
λ2t, (1-25)
where C1 and C2 are constants. Several possibilities exist for the behavior of the system based on the
eigenvalues, or alternately, the values of viscous damping constant c, the mass m, and the stiffness k. When
c = 0 the system behavior reverts to the undamped case. If c2 > 4mk then both eigenvalues are real and
less than zero, and the solution given by Eq. (1-25) decays to zero in a nonoscillatory manner. This type of
behavior is referred to as overdamped motion. If c2 < 4mk then the eigenvalues are complex conjugates with
negative real parts, and the solution x (t) decays to zero in an oscillatory manner, This behavior is referred
to as underdamped motion.
ME 349 Dynamic Systems Lab Manual 1-11
The case of c2 = 4mk results in equal eigenvalues,
λ1 = λ2 = −
c
2m
, (1-26)
and this case has a somewhat different solution
x(t) = C1eλt + C2teλt. (1-27)
Since the eigenvalues in this case are both real and negative, the solution x(t) decays to zero with time as
in the overdamped case. The behavior in this case is known as critically damped motion. It possesses the
minimum value of damping coefficient, c, for which the motion remains nonoscillatory. It is also the case
which decays the most rapidly. The value of c associated with critically damped motion,
c = 2
√
km, (1-28)
is known as the critical damping coefficient.
A variation of Eq. (1-22) can be formed by first dividing by m to obtain
x¨+
c
m
x˙+
k
m
x = 0. (1-29)
Defining the damping ratio, ζ, as the ratio of damping coefficient, c, over critical damping coefficient, Eq.
(1-28), gives
ζ =
c
2
√
km
. (1-30)
By observing that cm = 2ζωn and that ω
2
n =
k
m , Eq. (1-29) can be rewritten as
x¨+ 2ζωnx˙+ ω2nx = 0. (1-31)
Equation (1-31) is referred to as the standard form of the equation of motion for the system defined by Eq.
(1-22). The characteristic equation for the standard form is
λ2 + 2ζωnλ+ ω2n = 0, (1-32)
and the resulting eigenvalues are
λ = −ζωn ± ωn
q
ζ2 − 1 (1-33)
By substituting the eigenvalues for the undamped case, Eq. (1-33) with (ζ < 1), into Eq. (1-25) the
standard form of the solution becomes
x(t) = e−ζωnt
³
C1eiωn
√
1−ζ2t + C2e−iωn
√
1−ζ2t
´
. (1-34)
By applying an argument similar to that which was applied in the development of Eqs. (1-9) through (1-13),
the solution for the underdamped case can be written as
x(t) = e−ζωnt
∙
A cos
µ
ωn
q
1− ζ2 t
¶
+B sin
µ
ωn
q
1− ζ2 t
¶¸
, (1-35)
where A and B are real-valued constants. Two primary differences exist between the solution for the
undamped case, Eq. (1-13), and the underdamped case, Eq. (1-35): First, the exponential term in Eq.
(1-35) causes the amplitude of the underdamped solution to decay with time, whereas the amplitude of the
undamped case is constant. Secondly, the frequency of oscillation has changed from ωn in the undamped
case to
ωd = ωn
q
1− ζ2 (1-36)
in the damped case. The reader should be aware that, for the underdamped case, ωn still refers to the
natural frequency as defined in Eq. (1-7). The frequency, ωd, in Eq. (1-36) is the frequency at which the
1-12 Schimel, Ding, Anderson and Grantham
x(0)
> 1
< 1
= 1
t
x(t)
Figure 1-3: Displacement vs. time traces for various damping ratios.
underdamped system vibrates and it is known as the damped natural frequency. In cases where the damping
is light (i.e. ζ < 0.10) the natural frequency and damped natural frequency are approximately equal.
The standard form for the overdamped case can be written as
x(t) = Ae
−ωn

ζ−
√
ζ2−1

t
+Be
−ωn

ζ+
√
ζ2−1

t
, (1-37)
and for the critically damped case
x(t) = (A+Bt)e−ζωnt, (1-38)
where A and B are again real-valued constants that are determined from initial conditions.
A plot of the solutions x (t) for the underdamped, critically damped, and overdamped cases appears in
Figure 1-3.
1.3 Apparatus, Measurements, and Calculations
Measurement of Damped Natural Frequency
For this laboratory, a mechanical strip-chart recorder is to be used to capture the free vibration motion of a
simple single-degree-of-freedom system. The strip-chart recorder consists of a marking pen affixed to a mass
carriage, which is movable in the vertical direction, and a continuous strip of paper which, when driven by
a rotating drum, moves in a horizontal direction past the pen. A drawing of the single-degree-of-freedom
system, with the strip-chart recorder, is shown in Figure 1-4, and an example free vibration record is shown
in Figure 1-5.
The damped natural frequency of the vibrating system is easily computed from the strip-chart recording,
if the velocity of the paper as it moves horizontally past the pen is known. The paper velocity, vp, is
determined by measuring the length of paper, lp, that moves past the pen in a set amount of time, tp. The
paper velocity is simply
vp =
lp
tp
. (1-39)
ME 349 Dynamic Systems Lab Manual 1-13
spring
mass
carriage
masses
dashpot
controller
pen
assemblystrip
chart
Figure 1-4: Apparatus for testing the single-degree-of-freedom system.
5 1/2  cycles
x1
xn
Figure 1-5: Example of a strip-chart recording from a free vibration test.
1-14 Schimel, Ding, Anderson and Grantham
The damped natural frequency is determined from the free vibration record by measuring the length, losc,
over which a given number of oscillations occur. The damped natural frequency in radians per second is
given by the equation
ωd = 2πvp
µ
osc. count
losc
¶
. (1-40)
The accuracy of this computation is improved by including all measurable cycles from the strip-chart record-
ing in the oscillation count, including half cycles.
Linear Regression
Linear regression methods will be utilized in some of the calculations for this laboratory. Therefore, a review
of these methods is presented in this section.
Consider a set of data, (xj , yj), j = 1, . . . , q, for which it is desirable to compute a best fit linear
relationship between the ordered pairs. Recall that the slope-intercept equation for a line is given by
y = ax+ b, (1-41)
where y is the dependent variable, x is the independent variable, a is the slope, and b is the intercept. Assume
that it is the objective to determine values of a and b which will minimize the error between the right- and
left-hand sides of Eq. (1-41) for all pairs (xj , yj) by the least squares criterion. That is, the objective is to
minimize the function
E =
qX
j=1
(yj − axj − b)2. (1-42)
Values of a and b which minimize E occur where the partial derivatives of E with respect to a and b are
equal to zero:
∂E
∂a
= 2
qX
j=1
(yj − axj − b)(−xj) = 0, (1-43)
and
∂E
∂b
= 2
qX
j=1
(yj − axj − b)(−1) = 0. (1-44)
Solving for a gives
a =
q
Pq
j=1 xjyj −
Pq
j=1 xj
Pq
j=1 yj
q
Pq
j=1 x
2
j −
³Pq
j=1 xj
´2 , (1-45)
or, alternately,
a =
Pq
j=1(xj − x¯)(yj − y¯)Pq
j=1(xj − x¯)2
, (1-46)
where x¯ and y¯ denote mean values for the data xj and yj , respectively. The variable a is known as the
regression coefficient, and is equivalent the slope of the regression line through the data points. Finally,
substituting x¯, y¯ and a into Eq. (1-41) and solving for b gives the y-intercept of the regression line
b = y¯ − ax¯. (1-47)
As an example, consider the system in Figure 1-1a. The spring constant of this type of system is easily
determined by a static displacement test, where the applied force is gravitational. Analysis of the free-body
diagram in Figure 1-1b, with gravity included, gives the equation
mjg = kxj . (1-48)
The subscripting is used to denote the jth displacement xj associated with the applied gravitational force
mjg. A sequence of displacement-force pairs (xj ,mjg), j = 1, . . . , q, will theoretically fall on a straight line.
Linear regression analysis can be applied directly to the data set. The gravitational force mjg acts as the
jth dependent variable, and xj as the jth independent variable. Computation of the regression coefficient
gives the slope of the line, which in this example, is a least squares estimate of the stiffness k.
ME 349 Dynamic Systems Lab Manual 1-15
Logarithmic Decrement Method
The logarithmic decrement method is a useful technique for estimating viscous damping in a single-degree-
of-freedom system by using the free vibration trace. A modified version of this method, which incorporates a
linear regression, will be presented here. Consider the solution to the equation of motion for the underdamped
case, when x(0) = xo and x˙(0) = 0 :
x(t) = xoe
−ζωnt cos (ωdt) . (1-49)
A positive peak occurs at t = 0, and then at intervals of t = 2π/ωd. Negative peaks occur starting at
t = π/ωd, and then proceed at the same interval. Consider a ratio of the height of the first negative peak to
the height of a some arbitrary negative peak occurring later in the trace (at t = (2n+ 1)π/ωd):
x1
xn
=
x(π/ωd)
x ((2n+ 1)π/ωd)
=
xoe−ζωnπ/ωd cos (π)
xoe−ζωn(2n+1)π/ωd cos ((2n+ 1)π)
. (1-50)
By canceling like terms the ratio reduces to
x1
xn
= e2πnζωn/ωd (1-51)
Taking the natural logarithm of both sides, and applying the definition of damped natural frequency, gives
the nth logarithmic decrement
δn = ln
µ
x1
xn
¶
=
2πnζp
1− ζ2
. (1-52)
For small damping, Eq. (1-52) can be approximated by
δn ∼= 2πnζ. (1-53)
An equivalent relationship can be obtained for the positive peaks.
Observe that for Eq. (1-53) a sequence of pairs, (δn, n), n = 1, . . . , r, taken from the same free vibration
trace, fall approximately on a straight line. Therefore, linear regression analysis can be applied to the data
set. The nth logarithmic decrement, δn, acts as the dependent variable and 2πn acts as the independent
variable. Computation of the regression coefficient gives the slope of the line, which in this case is a least
squares estimate of the damping ratio, ζ.
Measurement of Viscous Damping Coefficient
It is possible to compute the viscous damping coefficient, c, for a dashpot from a constant force vs. velocity
record. For this laboratory the dashpot of the single-degree-of-freedom system, shown in Figure 1-4, will be
subject to a constant gravitational load in a drop test. The gravitational load is applied by the mass carriage,
with the spring detached. The test is conducted by raising the mass carriage to the top of its travel, with
the dashpot attached, and then releasing the carriage. After a constant velocity is achieved, the equation of
motion for the mass-dashpot system reduces to
mg = cx˙. (1-54)
Solving for viscous damping coefficient gives
c =
mg
x˙
. (1-55)
The drop velocity, x˙, is determined from a strip-chart recording of a drop test. An example recording is
shown in Figure 1-6. If the length of the drop, ld, and the length of the paper travel, lp, are measured from
the drop test record, and the paper velocity is determined as in Section 1.3, then the drop velocity can be
computed from the equation
x˙ = vp
ld
lp
. (1-56)
1-16 Schimel, Ding, Anderson and Grantham
ld
lp
Figure 1-6: Example of a strip-chart recording from a drop test.
1.4 Laboratory #1 Procedures
This laboratory will consist of a series of experiments directed at determining physical properties of a single-
degree-of-freedom mass-spring-damper system. Two types of tests will be conducted: tests in which the
spring or dashpot are isolated so that coefficients can be computed independently, and tests in which mea-
surements from a free vibration test are used to determine the coefficients through specialized computations.
The measurement methods and results will be compared and contrasted for a more in-depth understanding
of simple vibratory systems and measurement techniques.
Spring Constant
1. Install a pen in the pen holder, but do not allow the pen tip to rest on the paper.
2. Remove any of the slotted masses from the mass carriage.
3. By moving the mass carriage up and down, ensure that the dashpot does not interfere with the motion
of the mass carriage. Check the dashpot valve, which is located on the piston, to ensure that it is open
sufficiently, and check that the piston does not drag on the dashpot cylinder wall. If necessary, open
the valve or loosen the thumb-nut at the dashpot base slightly.
4. By means of the thumb nuts on the spring support, adjust the height of the mass carriage so that the
pen will make contact near the top of, but not above, the graph paper.
5. Position the pen so that it contacts the graph paper.
6. Switch on the graph-paper feed, and allow the pen to draw a baseline approximately 5-cm. in length.
7. Switch the paper feed off.
8. Attach one slotted 1-kg mass to the mass carriage, and allow the carriage to stabilize.
ME 349 Dynamic Systems Lab Manual 1-17
9. Record the displacement of the mass carriage by advancing the graph-paper feed approximately 2-cm.
10. Add the remaining weights, one at a time, and record each displacement each time. Note: Only two
to three weights can be supported by the light spring.
11. Remove the mass-displacement record from the strip-chart recorder, and use a ruler to measure the
total vertical displacement from the baseline for each added mass.
12. Record the measured displacements and corresponding masses in Tables 1-1 through 1-3.
13. Repeat steps 3 through 13 for the two other springs.
14. Remove the pen from contact with the paper.
15. Use a linear regression to compute the stiffness of each spring.
Table 1-1. Light Spring Data
Added Mass (kg) Total Mass (kg) Force (N) Deflection (m)
Table 1-2. Medium Spring Data
Added Mass (kg) Total Mass (kg) Force (N) Deflection (m)
Table 1-3. Heavy Spring Data
Added Mass (kg) Total Mass (kg) Force (N) Deflection (m)
Damped Natural Frequency and Damping Ratio
1. Select a mass-spring-dashpot combination from the tables below. Perform the cases which include the
use of dashpot oil last, because the same setup will be used later.
2. By means of the thumb nuts on the spring support, adjust the height of the mass carriage so that the
pen will make contact near the center of the graph paper.
3. Ensure that the dashpot does not interfere with the motion of the mass carriage by moving the mass
carriage up and down as before.
4. Move the pen into contact with the graph paper
5. Switch on the graph paper feed and allow the pen to draw a centerline of approximately 5-cm.
6. Switch the unit off.
7. Depress the mass carriage to a point near but not below the bottom of the graph paper and hold it
stationery.
1-18 Schimel, Ding, Anderson and Grantham
8. Switch the graph paper feed on and then release the mass carriage. Do not interfere with the system
in any way while it vibrates.
9. Allow the graph paper trace to continue approximately 2-cm beyond the point where the carriage stops
vibrating.
10. Switch the graph-paper feed off, and remove the pen from contact with the graph paper.
11. Perform damped natural frequency and log decrement measurements and record them in the tables
below.
12. Repeat steps 2 through 10 for the mass-spring-dashpot combinations remaining in the tables.
Table 1-4. Damped Natural Frequency Data
Case Light Spring Medium Spring Heavy Spring Medium Spring
no Dashpot Oil no Dashpot Oil no Dashpot Oil w/ Dashpot Oil
Mass 1.9kg 2.9kg 3.9kg 5.9kg 3.9kg 6.9kg 3.9kg 5.9 kg
Cycle
Count
Length,
losc, (m)
Table 1-5. Free Vibration Peak Amplitudes (m)
Light Spring Medium Spring Heavy Spring Medium Spring
Peak no Dashpot Oil no Dashpot Oil no Dashpot Oil w/ Dashpot Oil
Number 1.9kg 2.9kg 3.9kg 5.9kg 3.9kg 6.9kg 3.9kg 5.9 kg
1
2
3
4
5
6
7
8
9
10
11
12
Viscous Damping Coefficient
1. Do not change the dashpot valve setting from its setting in the previous experiment.
2. Remove any of the slotted 1-kg masses form the mass carriage.
3. Disconnect and remove the spring from the apparatus.
4. Raise the mass carriage so that the pen will make contact near the top of, but not above, the graph
paper.
5. Put the pen into contact with the graph paper and start the graph paper feed.
6. Release the mass carriage suddenly and do not disturb it during its downward travel.
ME 349 Dynamic Systems Lab Manual 1-19
7. Switch off the graph paper feed and remove the pen from contact with the paper.
8. Record the length of the drop, ld, and the length of the paper travel, lp.
ld = __________
lp = __________
Chapter 2
Forced Vibration
2.1 Equations of Motion
Vibration Generated by a Force Input
Recall the single-degree-of-freedom mass-spring-dashpot system from Chapter 1. Consider the effects of
a time-varying force, f(t), applied to a mass, starting from rest as shown in Figure 2-1a. Application of
Newton’s second law to the free-body diagram shown in Figure 2-1b gives the equation
−kx− cx˙+ f(t) = mx¨, (2-1)
or, after rearrangement,
mx¨+ cx˙+ kx = f(t). (2-2)
Equation (2-2) is the equation of motion for the system shown in Figure 2-1. In standard form the equation
becomes
x¨+ 2ζωnx˙+ ω
2
nx =
f(t)
m
. (2-3)
Observe that Eq. (2-2) has the same form as the equation of motion for the freely vibrating mass-spring-
dashpot system in Eq. (1-22), with the exception of the force input on the right-hand side.
Vibration Generated by a Displacement Input
Consider a single-degree-of-freedom mass-spring-dashpot system with a time varying displacement, y(t),
which is input to one end of the spring as shown in Figure 2-2a. The spring force, fs, imparted to the mass
is a function of the relative displacement between ends of the spring, and is given by
fs = −k(x− y(t)). (2-4)
Application of Newton’s second law to the free-body diagram shown in Figure 2-2b gives the equation
−k(x− y(t))− cx˙ = mx¨, (2-5)
or, after rearrangement,
mx¨+ cx˙+ kx = ky(t). (2-6)
Equation (2-6) is the equation of motion for the single-degree-of-freedom system with displacement input,
as shown in Figure 2-2. In standard form Eq. (2-6) becomes
x¨+ 2ζωnx˙+ ω2nx = ω
2
ny(t) (2-7)
Observe that the spring force, ky (t), generated by the displacement input in Eq. (2-6), has the same effect
as the force input, f(t), in the equation of motion given by Eq. (2-2).
2-20
ME 349 Dynamic Systems Lab Manual 2-21
  .  
Figure 2-1: The single-degree-of-freedom mass-spring-dashpot system with force input.
  .  
Figure 2-2: The single-degree-of-freedom mass-spring-dashpot system with displacement input.
2-22 Schimel, Ding, Anderson and Grantham
2.2 Laplace Transforms and the Transfer Function
Equations (2-2) and (2-6) are linear, inhomogeneous, second-order differential equations in input-output
form. The general form for an equation of this type is written as
an
dnx
dtn
+ an−1
dn−1x
dtn−1
+ · · ·+ a1 dxdt + a0x = bm
dmf
dtm
+ bm−1
dm−1f
dtm−1
+ · · ·+ b1 dfdt + b0f, (2-8)
where d
j
dtj denotes the j
th-order time derivative with aj and bj as a constant coefficients. The solution, x(t),
for an equation in the form of Eq. (2-8) is readily obtained by application of the Laplace transform method.
The Laplace transform of a function x (t) is defined as
L (x(t)) = X(s) ≡
Z ∞
0
x(t)e−stdt. (2-9)
The Laplace transform transforms a function by changing its dependence on the time-variable t to depen-
dence on the variable s. Once transformed into the s-domain certain operations become easier to perform.
Application of the Laplace transform method to an equation in the form of Eq. (2-8) will generate an
equation that can be solved for X(s) algebraically. Application of the inverse Laplace transform,
L−1 (X(s)) = x(t) ≡
Z ∞
0
X(s)estds, (2-10)
results in the solution to the equation of motion.
An important property of the Laplace transform is linearity:
L (au(t) + bv(t)) = aL (u(t)) + bL (v(t)) = aU(s) + bV (s). (2-11)
Also, the Laplace transform of the nth-order time derivative of a function with zero initial conditions is
L
µ
dnu
dtn
¶
= snU(s). (2-12)
As will be seen later, a useful inverse Laplace transform is the following one, which gives the exponential
function:
L−1
µ
1
s− a
¶
= eat. (2-13)
Consistent with the Laplace transform and integration in general, the inverse Laplace transform also possesses
the linearity property
L−1 (aU(s) + bV (s)) = aL−1 (U(s)) + bL−1 (V (s)) = au(t) + bv(t). (2-14)
Any of the other rules of integration apply to the Laplace transform and inverse Laplace transform as well.
With zero initial conditions, application of Eqs. (2-11) and (2-12) to Eq. (2-8) gives¡
ansn + an−1sn−1 + · · ·+ a1s+ a0
¢
X(s) =
¡
bmsm + bm−1sm−1 + · · ·+ b1s+ b0
¢
F (s), (2-15)
or, after rearrangement,
G(s) =
X(s)
F (s)
=
bmsm + bm−1sm−1 + · · ·+ b1s+ b0
ansn + an−1sn−1 + · · ·+ a1s+ a0 . (2-16)
Equation (2-16) is known as the transfer function for the system with equation of motion given by Eq. (2-8).
The transfer function, G(s), is the ratio of the Laplace transform of the output, X(s), over the input, F (s).
ME 349 Dynamic Systems Lab Manual 2-23
2.3 Solutions to the Equations of Motion
Force Input System
The Laplace transform of the equation of motion for the system with force input, Eq. (2-3), is
(s2 + 2ζωns+ ω
2
n)X(s) =
1
m
F (s) (2-17)
and the resulting transfer function is
Gf (s) =
X(s)
F (s)
=
1/m
s2 + 2ζωns+ ω2n
. (2-18)
A specific equation for the forcing function must be selected before a particular solution to the equation of
motion can be obtained. For this laboratory it will be assumed that the forcing function is the sinusoid
f(t) = Fo sinΩt, (2-19)
where Fo is the force amplitude, and Ω the driving frequency. The Laplace transform of Eq. (2-19) is
F (s) =
FoΩ
s2 +Ω2 . (2-20)
Solving Eq. (2-18) for X(s), and substituting the right-hand-side of Eq. (2-20) gives the equation
X(s) = Gf (s)F (s) =
1/m
s2 + 2ζωns+ ω2n
FoΩ
s2+Ω2 . (2-21)
Application of the method of partial fraction expansion to Eq. (2-21) enables X(s) to be expressed in
the following format, for which the underdamped case (ζ < 1) has been assumed:
X(s) =
C1
s−
³
−ζωn + iωn
p
1− ζ2
´ + C2
s−
³
−ζωn − iωn
p
1− ζ2
´ + C3
s− iΩ +
C4
s+ iΩ . (2-22)
The inverse Laplace transform results in the equation
x(t) = C1e

−ζωn+iωn
√
1−ζ2

t
+ C2e

−ζωn−iωn
√
1−ζ2

t
+ C3eiΩt + C4e−iΩt. (2-23)
Observe that Eq. (2-23) can be rewritten as
x(t) = e−ζωn
³
C1e
iωn
√
1−ζ2 + C2e
−iωn
√
1−ζ2
´
+ C3e
iΩt + C4e
−iΩt. (2-24)
By applying Euler’s formula, Eq. (1-10), Eq. (2-24) becomes
x(t) = e−ζωnt
∙
A cos
µ
ωn
q
1− ζ2t
¶
+B sin
µ
ωn
q
1− ζ2t
¶¸
+D1 sinΩt+D2 cosΩt, (2-25)
where A, B, D1, and D2 are constants. Equation (2-25) can be broken into two parts: a time-decaying
or transient part, xT (t), and a constant amplitude or steady-state part, xSS(t). Hence, Eq. (2-25) can be
written as
x(t) = xT (t) + xSS(t), (2-26)
where, in terms of solutions to differential equations, xT (t) forms the homogeneous solution, and xSS(t) forms
the particular solution to the equation of motion, Eq. (2-3). The sum of the homogenous and particular
solutions form the general solution, Eq. (2-26).
As t approaches infinity, the contribution from the transient part of the solution decays to zero, leaving
the steady-state part of the solution to satisfy any effects from the force input. The last two terms of Eq.
(2-25) are constant amplitude sinusoids, which form the steady-state part of the solution
xSS(t) = D1 sinΩt+D2 cosΩt. (2-27)
2-24 Schimel, Ding, Anderson and Grantham
Since the steady-state part of the solution must account for the effects due to the force input, the constantsD1
and D2 are determined by substituting Eqs. (2-27) and (2-19) into Eq. (2-3 ). Performing this substitution
and separating the sine and cosine terms (which are independent of each other) gives the equations
¡
−D1Ω2 −D2Ω(2ζωn) +D1ω2n
¢
sinΩt = Fo
m
sinΩt (2-28)
and ¡
−D2Ω2 +D1Ω(2ζωn) +D2ω2n
¢
cosΩt = 0. (2-29)
Solving for D1 and D2 results in
D1 =
(Fo/k)(1− α2)
(1− α2)2 + (2αζ)2 (2-30)
and
D2 =
−(Fo/k)(2αζ)
(1− α2)2 + (2αζ)2 , (2-31)
where
α =
Ω
ωn
(2-32)
defines the dimensionless ratio of input frequency to natural frequency. Substituting D1 and D2 into Eq.
(2-27) gives the steady-state part of the solution
xSS (t) =
Fo/k
(1− α2)2 + (2αζ)2
¡
(1− α2) sinΩt− 2αζ cosΩt
¢
. (2-33)
Observe that the steady-state response is composed of both sine and cosine terms and, therefore, differs in
phase from the input function.
The first term of Eq. (2-25) forms the transient part of the solution,
xT (t) = e−ζωnt
∙
A cos
µ
ωn
q
1− ζ2 t
¶
+B sin
µ
ωn
q
1− ζ2 t
¶¸
, (2-34)
which in this case is equivalent to the solution for underdamped free vibration, Eq. (1-35). The constants
A and B are determined from initial conditions in the general solution
x(t) = e−ζωnt
h
A cos
³
ωn
p
1− ζ2 t
´
+B sin
³
ωn
p
1− ζ2 t
´i
+ Fo/k(1−α2)2+(2αζ)2
¡
(1− α2) sinΩt− 2αζ cosΩt
¢ . (2-35)
Displacement Input System
The Laplace transform of the equation of motion for the system with displacement input, Eq. (2-7), is
(s2 + 2ζωns+ ω
2
n)X(s) = ω
2
nY (s) (2-36)
and the resulting transfer function is
Gd(s) =
X(s)
Y (s)
=
ω2n
s2 + 2ζωns+ ω2n
. (2-37)
It will be assumed that the displacement input is the sinusoid
y(t) = Yo sinΩt, (2-38)
which has the Laplace transform
Y (s) =
YoΩ
s2 +Ω2 . (2-39)
ME 349 Dynamic Systems Lab Manual 2-25
t
t
t
t
f(t)
xSS(t)
xT(t)
x(t)
Figure 2-3: Input and response for the single degree of freedom system with force input.
Solving Eq. (2-37) for X(s), and substituting the right-hand-side of Eq. (2-39) gives the equation
X(s) = Gd(s)Y (s) =
ω2n
s2 + 2ζωns+ ω2n
YoΩ
s2 +Ω2 . (2-40)
Following a procedure similar to that applied to Eqs. (2-22) through (2-35) yields the general solution for
the displacement input system
x(t) = e−ζωnt
h
A cos
³
ωn
p
1− ζ2t
´
+B sin
³
ωn
p
1− ζ2t
´i
+ Yo(1−α2)2+(2αζ)2
¡
(1− α2) sinΩt− 2αζ cosΩt
¢
.
(2-41)
A plot of the force input and response for an underdamped system, started from a state of rest, appears
in Figure 2-3. The response is typical of both the force input and displacement input systems.
2.4 The Frequency Response Function
Generalized System
The frequency response function is closely related to the transfer function of a system. For the frequency
response function it is assumed that the input to a system is a sinusoidal function, and that the system
output is a steady state sinusoidal response at the frequency of the input. By making the assumption that
the input and response are equivalent in frequency and at steady-state conditions, the resulting solution
ignores any possible transient behavior. In forced vibration analysis the transient behavior often dies out
quickly, and it is the steady-state response that is of primary interest. For such cases the frequency response
function or FRF gives a compact representation of the magnitude and phase angle relationships between
the input and output.
Consider the generalized equation of motion given by Eq. (2-8) which results in the transfer function
given by Eq. (2-16). Substituting s = iΩ into Eq. (2-16) gives
G (iΩ) = X (iΩ)
F (iΩ) =
bm (iΩ)m + bm−1 (iΩ)m−1 + · · ·+ b1 (iΩ) + b0
an (iΩ)n + an−1 (iΩ)n−1 + · · ·+ a1 (iΩ) + a0
. (2-42)
Equation (2-42) defines the frequency response function, G (iΩ), for the generalized system defined by Eq.
(2-8). The frequency response function is, in general, a complex-valued function, which represents the
2-26 Schimel, Ding, Anderson and Grantham
magnitude and phase angle relationship between the input and output. The complex-valued FRF can be
split into the real and imaginary parts, Re {G (iΩ)} and Im {G (iΩ)} respectively. It can be shown that the
magnitude of the FRF ,
|G (iΩ)| =
q
Re {G (iΩ)}2 + Im {G (iΩ)}2, (2-43)
defines the ratio of the response amplitude over the input amplitude at a particular frequency Ω. The
function,
∠G (iΩ) = tan−1 Im {G (iΩ)}
Re {G (iΩ)} , (2-44)
defines the phase angle by which the response lags behind the input at the same frequency.
Force Input System
Replacing every occurrence of s in Eq. (2-18) with iΩ gives the frequency response function for the force
input system:
Gf (iΩ) =
X(iΩ)
F (iΩ) =
1/m
−Ω2 + 2iζωnΩ+ ω2n
=
1/m
(ω2n − Ω2) + 2iζωnΩ
. (2-45)
Dividing numerator and denominator by ω2n and recalling the definition for the frequency ratio, α =
Ω
ωn
,
gives
Gf (iΩ) =
1/k
(1− α2) + 2iζα (2-46)
The real and imaginary parts of Gf (iΩ) can be obtained by multiplying the numerator and denominator
by the complex conjugate of the denominator of Eq. (2-46), yielding
Gf (iΩ) =
1/k
(1− α2) + 2iζα
µ
(1− α2)− 2iζα
(1− α2)− 2iζα
¶
=
(1/k)
¡
(1− α2)− 2iζα
¢
(1− α2)2 + (2ζα)2 . (2-47)
The real and imaginary parts are
Re {Gf (iΩ)} = (1/k)(1− α
2)
(1− α2)2 + (2ζα)2 (2-48)
and
Im {Gf (iΩ)} = − (1/k)(2ζα)
(1− α2)2 + (2ζα)2 . (2-49)
Application of Eq. (2-43) gives the magnitude
|Gf (iΩ)| =
s
(1/k)2 ((1− α2)2 + (2ζα)2)
((1− α2)2 + (2ζα)2)2
=
1/kp
(1− α2)2 + (2ζα)2
. (2-50)
A more direct approach to obtaining the magnitude is to take the square root of the quantity Gf (iΩ)
multiplied by its own complex conjugate G∗f (iΩ):
|Gf (iΩ)| =
q
Gf (iΩ)G∗f (iΩ) =
sµ
1/k
(1− α2) + 2iζα
¶µ
1/k
(1− α2)− 2iζα
¶
=
1/kp
(1− α2)2 + (2ζα)2
(2-51)
Observe that the magnitude of the FRF, |Gf (iΩ)|, is highly dependent on the frequency ratio, α = Ωωn .
Values of α near one result in large-amplitude steady-state responses, while values of α away from one result
in reduced amplitudes. The value of the frequency ratio, α, equal to one has a special significance; it is
termed resonance: the condition where the driving frequency, Ω, is equal to the system’s natural frequency,
ωn, and it is marked by large response amplitudes. An engineer may, in some cases, want to design for
resonant behavior, such as in industrial shaker design, but in most cases an engineer would want to design
away from resonance, such as in automotive suspension or machine tool designs. Figure 2-4 shows example
ME 349 Dynamic Systems Lab Manual 2-27
0 0.5 1 1.5 2 2.5
0
1
2
3
4
5
6
|G(i   )|
  n 
 = 
= 0.1
= 0.2
= 0.4
= 0.6
= 0.8
Figure 2-4: Frequency response function magnitude vs. frequency ratio for the single-degree-of-freedom
mass-spring-dashpot system.
FRF magnitude traces as a function of frequency ratio for various levels of damping. Note that the peaks
tend toward infinity as the damping ratio approaches zero. Application of Eq. (2-44) gives the phase angle
relationship
∠Gf (iΩ) = tan−1
µ
−2ζα
1− α2
¶
. (2-52)
Figure 2-5 shows phase angle traces as a function of the frequency ratio, α, for various levels of damping.
Observe that at resonance, the phase angle is equal to π/2 (90◦). That is, when the steady-state response is
at a maximum, the force input is at a minimum, and vice-versa.
Displacement Input System
A development similar to that in the previous section gives the FRF for the displacement input system
Gd(iΩ) =
ω2n
(ω2n − Ω2) + 2iζωnΩ
; (2-53)
or after division of the numerator and denominator by ω2n
Gd(iΩ) =
1
(1− α2) + 2iζα, (2-54)
where α = Ωωn is the frequency ratio. The FRF magnitude and phase angle relationships for the displacement
input system are given by
|Gd(iΩ)| = 1p
(1− α2)2 + (2ζα)2
, (2-55)
and
∠Gd (iΩ) = tan−1
µ
−2ζα
1− α2
¶
. (2-56)
2-28 Schimel, Ding, Anderson and Grantham
 = 
  n  
G(i   )
Figure 2-5: Frequency response function phase angle vs. frequency ratio for the single-degree-of-freedom
mass-spring-dashpot system.
Plots of the magnitude and phase angle as a function of frequency ratio for the displacement input FRF are
similar to those shown in Figures 2-4 and 2-5.
2.5 Apparatus, Measurements and Calculations
Some of the same techniques for frequency and damping ratio measurement which were used in Laboratory
#1 will be used in this laboratory. Any new techniques are described below. The primary objective of the
measurement and calculation methods presented below are to generate a FRF magnitude plot, from which
damped natural frequency, ωd, and damping ratio, ζ, estimates will be made.
Measurement of the Response Amplitude and Frequency
This laboratory will use a single-degree-of-freedom testing apparatus, like the one used in Laboratory #1,
but equipped for either force or displacement input tests. The testing machines are described in more detail
below. The frequency and amplitude of the steady-state response, for a particular input frequency, are
measured from a strip-chart record, an example of which is shown in Figure 2-6. The record is generated
by applying a input at a selected frequency and, after waiting for any transients to die away, recording the
steady-state motion of the mass carriage on the strip-chart recorder. The response amplitude, Xo, is given
by the equation
Xo =
Xp-p
2
, (2-57)
where Xp-p is the peak-to-peak height of the response. The response frequency (which is equivalent to the
input frequency), is calculated as in Eqs. (1-39) and (1-40) of Chapter 1. However, the driving frequency, Ω,
replaces the damped natural frequency, ωd, on the left-hand-side of Eq. (1-40). The frequency measurement
baseline (the horizontal line which crosses the sinusoidal response trace in Figure 2-6) needs to be parallel
to, but not necessarily equidistant from, the two amplitude measurement lines in the figure.
ME 349 Dynamic Systems Lab Manual 2-29
xp-p
Figure 2-6: An example strip-chart recording of a forced vibration test.
Measurement of the Input Amplitude
Force Input System
For the theoretical system depicted in Figure 2-1, it has been assumed that the force input, f(t), is a sinusoid.
For the single-degree-of-freedom testing machine used in the laboratory, Figure 2-7, the sinusoidal force, f(t),
is imparted onto the mass carriage through a set of counter-rotating masses. The masses are equal in size
and rotate at the same rate; but in opposite directions, as shown in Figure 2-8. Each mass is oriented so
that, as it rotates, it mirrors the other mass in position on a vertically symmetric axes. At some instant in
time the mass centers point in exactly opposite directions. Assume this instant corresponds to t = 0. At
this instant the centrifugal force generated by the rotating masses are in the opposite directions, and the net
force is given by
f(0) =
mr
2
Ω2r − mr
2
Ω2r = 0, (2-58)
where mr2 is the mass of a single counter rotating weight, Ω is the rotation frequency, and r is the distance
from rotation center to mass center of a rotating weight. At some arbitrary time later, one mass has rotated
through angle Ωt while the other has rotated through angle −Ωt. At this instant the vector diagram shown
in Figure 2-9 reveals that the net vertical force is given by
f(t) =
mr
2
Ω2 sinΩt− mr
2
Ω2r sin(−Ωt) = mrΩ2r sinΩt, (2-59)
and that the net horizontal force cancels:
mr
2
Ω2r cosΩt− mr
2
Ω2r cos(−Ωt) = 0. (2-60)
Therefore, the counter rotating masses impart to the mass carriage a sinusoidally varying force in the vertical
direction and zero force in the horizontal direction.
The equation of motion for this specific case can be written as
mx¨+ cx˙+ kx = mrΩ2r cosΩt, (2-61)
2-30 Schimel, Ding, Anderson and Grantham
spring
mass
carriage
masses
dashpot
controller
pen
assemblystrip
chart
counter rotating
      masses
Figure 2-7: The single degree of freedom testing machine configured for force input tests.
r r
Figure 2-8: Arrangement of the counter rotating masses.
ME 349 Dynamic Systems Lab Manual 2-31
spring
mass
carriage
masses
dashpot
controller
pen
assemblystrip
chart crank
arm
connecting
rod
slider
Figure 2-9: The single degree of freedom apparatus configured for displacement input tests.
where, in this case, the mass, m, includes the mass of the carriage, the mass of any of the slotted disks, and
the mass of the rotating masses, mr. As in Laboratory #1 the spring mass is considered to be negligible.
Determination of the amplitude of the sinusoidal force, Fo, applied to the mass carriage is achieved by
first measuring the frequency of the steady-state response. The response frequency is measured from a strip-
chart record of the response, Figure 2-6, and is theoretically the same as the force input frequency, Ω. From
the input frequency, Ω, and the constants that appear in Eq. (2-58) the sinusoidal force amplitude, Fo, can
be calculated.
Displacement Input System
The displacement input for the theoretical system depicted in Figure (2-2) occurs at the top of the spring.
For the single-degree-of-freedom testing machine, this input is achieved through the vertically oriented slider-
crank mechanism as shown in Figure (2-9). The spring has been moved in Figure (2-9) for clarity. However,
it retains its normal position during testing. For this apparatus the displacement input amplitude, Yo, is
equal to the crank arm radius, which is a constant.
Generation of the FRF Magnitude Plot
A magnitude plot of the frequency response function, for either the force or displacement input systems, can
be generated from a series of steady-state tests in which the frequency is incrementally increased. The ratio
of response amplitude to input amplitude gives the FRF magnitude at a particular frequency ratio. For the
force input test the FRF magnitude is given by¯¯
G¯f (iΩ)
¯¯
=
Xo
Fo
, (2-62)
and for the displacement input test ¯¯
G¯d(iΩ)
¯¯
=
Xo
Yo
. (2-63)
2-32 Schimel, Ding, Anderson and Grantham
0 0.5 1 1.5
0
1
2
3
4
5
6
7
8
2
1
  n 
  =  
Figure 2-10: An example FRF magnitude plot.
Keep in mind that Eqs. (2-62) and (2-63) represent experimental values of the FRF at specific driving
frequencies, Ω. The experimental representations are used to predict the values of the theoretical coefficients
in Eqs. (2-50) and (2-55), respectively. An example FRF magnitude plot appears in Figure 2-10 The
frequency ratio, α, in the plot can be computed from the driving frequency, Ω, and the natural frequency,
ωn. The natural frequency is computed from a free-vibration test. Note that it is useful to collect as much
frequency data near the resonant frequency as is possible, since the accuracy of the calculations described
will depend on accurate determination of the FRF peak location.
Determination of Damping Ratio
Consider the FRF magnitude for the force input case, Eq.(2-50). When the frequency ratio, α, is equal to
one Eq.(2-50) reduces to
|Gf (iΩ)|α=1 =
1/k
2ζ
. (2-64)
Solving for the damping ratio, ζ, gives
ζ =
1
2
1
k |Gf (iΩ)|α=1
. (2-65)
Therefore, if the stiffness k is known or can be computed, the damping ratio can be readily computed from
the FRF magnitude at α = 1. A more simple result is obtained for the displacement input case:
ζ =
1
2
1
|Gd(iΩ)|α=1
(2-66)
The damping ratio measurement is depicted in Figure 2-10 for the displacement input case,
¯¯
G¯d(iΩ)
¯¯
. As
will be seen below the value of the frequency ratio α = 1 corresponds approximately to the resonance peak
maximum for lightly damped systems.
ME 349 Dynamic Systems Lab Manual 2-33
Determination of Natural Frequency
The frequency ratio that corresponds to the maximum of the theoretical FRF magnitude plot, Figure 2-10
can be determined by the first derivative test. Consider the system with force input. Computing the first
derivative of Eq. (2-50) with respect to frequency ratio, α, yields
∂ |Gf (iΩ)|
∂α
=
−(1/2k)
¡
2(1− α2)(−2α) + 2(2ζα)(2ζ)
¢
((1− α2)2 + (2ζα)2)−
3
2
. (2-67)
Application of the first derivative test gives the equation
−4α(1− α2) + 8ζ2α = 0, (2-68)
which has a root at
α =
q
1− 2ζ2. (2-69)
This is the value of α at which the magnitude of the FRF is a maximum. Applying the definition of
frequency ratio yields
ωn =
Ωp
1− 2ζ2
. (2-70)
For light amounts of damping, the denominator is approximately equal to one. Equation (2-70) also holds
for the displacement input case. Therefore, in the case of light damping, it can be assumed that the natural
frequency corresponds to the frequency at the maximum of
¯¯
G¯f (iΩ)
¯¯
or,
¯¯
G¯d(iΩ)
¯¯
as shown in Figure 2-10.
2.6 Laboratory #2 Procedures
Lab #2A - Vibrations Generated by Force Input
1. Check to see that the apparatus is setup for force input testing, and that it is configured for one of the
spring-mass-dashpot combinations listed in the tables below.
2. Set the vertical position of the strip-chart pen so that it will contact the graph paper near the center.
3. Conduct a free vibration test, as described in Chapter 1, and record the results in Tables 2-1 and 2-2.
4. Turn the motor speed adjustment knob to its lowest setting.
5. Switch on the motor drive and increase the speed until a measurable level of steady-state motion is
observed on the strip-chart recorder.
6. Once the motion has reached steady-state, switch on the recorder drive.
7. For accuracy, record at least two complete cycles and record over at least of 20-cm of strip-chart paper.
8. Switch off the recorder drive.
9. Increase the frequency by a small amount.
10. Wait for the system to reach steady-state.
11. Repeat steps 6 through 10 until the drive frequency is well beyond the natural frequency of the appa-
ratus.
12. Remove the pen from contact with the strip-chart recorder.
13. Perform frequency and amplitude measurements for each drive frequency and record the results in
table 2-3.
14. Repeat steps 2 through 13 for the remaining spring-mass-dashpot combinations.
2-34 Schimel, Ding, Anderson and Grantham
Lab #2B - Vibrations Generated by Displacement Input
1. Check to see that the apparatus is set up for displacement input testing. Before switching on the motor
drive, make sure that the slider mechanism is free to move in the vertical direction and has not been
locked in position.
2. Configure the apparatus for one of the spring-mass-dashpot combinations listed in the tables below.
3. Conduct a free vibration test, as described in Chapter 1, and record the results in Tables 2-1 and 2-2.
4. Turn the motor speed adjustment knob to its lowest setting.
5. Switch on the motor drive and advance the motor speed knob until the crankshaft begins to rotate.
6. When the slider mechanism reaches the middle of its travel, return the motor speed knob to the low
position and switch off the motor drive.
7. Set the vertical position of the strip-chart pen so that it will contact the graph-paper near the center,
when the slider mechanism is near the middle of its travel.
8. Switch on the motor drive and increase the speed until the mass carriage is moving in a smooth
sinusoidal motion.
9. Once the motion has reached steady-state, switch on the strip-chart recorder drive.
10. For accuracy, record at least two complete cycles and record over at least of 20-cm or strip-chart paper,
whichever is greater.
11. Switch off the recorder drive.
12. Increase the frequency by a small amount.
13. Wait for the system to reach steady-state.
14. Repeat steps 10 through 14 until the drive frequency is well beyond the natural frequency of the
apparatus.
15. Remove the pen from contact with the strip-chart recorder.
16. Perform frequency and amplitude measurements for each driving frequency and record the results in
the Table 2-3.
17. Repeat steps 2 through 17 for the remaining spring-mass-dashpot combinations.
ME 349 Dynamic Systems Lab Manual 2-35
Table 2-1. Damped Natural Frequency Data
Case Medium Spring Medium Spring
no Dashpot Oil w/ Dashpot Oil
Mass 3.9kg 4.9kg 3.9kg 4.9 kg
Cycle
Count
Length,
losc, (m)
Table 2-2. Free Vibration Peak Amplitudes (m)
Medium Medium Spring
Peak no Dashpot Oil w/ Dashpot Oil
Number 3.9kg 4.9kg 3.9kg 4.9 kg
1
2
3
4
5
6
7
8
9
10
11
12
Table 2-3. Forced Vibration Data (m)
Medium Spring Medium Spring
no Dashpot Oil w/ Dashpot Oil
Number 3.9kg 4.9kg 3.9kg 4.9 kg
# osc. losc Xp-p # osc. losc Xp-p # osc. losc Xp-p # osc. losc Xp-p
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
Chapter 3
Free Vibration of a
Two-Degree-of-Freedom System
3.1 Equations of Motion
Consider the undamped dual-mass-and-spring system shown in Fig 3-1a. Since either of the masses can be
independently positioned along the horizontal axis, the system has two degrees of freedom, which will be
denoted as x1 and x2. Application of Newton’s second law to the free body diagram for the mass associated
with variable x1, Figure 3-1 part (b), gives the equation of motion
−kx1 − k (x1 − x2) = mx¨1. (3-1)
Similarly, application of Newton’s second law to the free-body diagram for the mass associated with variable
x2, Figure 3-1c, gives
−k(x2 − x1)− kx2 = mx¨2 (3-2)
After rearrangement the two equations of motion become
mx¨1 + 2kx1 − kx2 = 0, (3-3)
and
mx¨2 − kx1 + 2kx2 = 0. (3-4)
Equations (3-3) and (3-4) can be represented in matrix form as
∙
m 0
0 m
¸½
x¨1
x¨2
¾
+
∙
2k −k
−k 2k
¸½
x1
x2
¾
=
½
0
0
¾
. (3-5)
The matrix containing the mass constants in Eq. (3-5) is known as the mass matrix, and the matrix
containing the spring constants is known as the stiffness matrix.
3.2 Solutions to the Equations of Motion
Equation (3-5) forms a coupled pair of linear, homogenous, second-order ordinary differential equations. The
solution to Eq. (3-5) is of the form ½
x1
x2
¾
=
½
A1
A2
¾
eλt. (3-6)
Substituting Eq. (3-6), and its second-order time derivative, into Eq. (3-5) gives
∙
m 0
0 m
¸½
A1
A2
¾
λ2eλt +
∙
2k −k
−k 2k
¸½
A1
A2
¾
eλt =
½
0
0
¾
. (3-7)
3-36
ME 349 Dynamic Systems Lab Manual 3-37
Figure 3-1: The two degree of freedom mass-spring system.
Combining terms and simplifying gives
∙
mλ2 + 2k −k
−k mλ2 + 2k
¸½
A1
A2
¾
eλt =
½
0
0
¾
. (3-8)
Since eλt is a nonzero number, the remaining part of Eq. (3-8) must equal zero. Therefore,
∙
mλ2 + 2k −k
−k mλ2 + 2k
¸½
A1
A2
¾
=
½
0
0
¾
. (3-9)
Cases of Eq. (3-9) in which A1 and A2 are not both zero are referred to as nontrivial solutions. Nontrivial
solutions to Eq. (3-9) are obtained by setting the determinant of the two-by two matrix in Eq. (3-9) equal
to zero, and solving for λ. The equation for the determinant is¯¯¯¯
mλ2 + 2k −k
−k mλ2 + 2k
¯¯¯¯
= 0, (3-10)
and computing the determinant gives
m2λ4 + 4mkλ2 + 3k2 = 0. (3-11)
Equation (3-11) is quadratic in λ2. Application of the quadratic formula gives
λ2 =
−4mk ±√16m2k2 − 12m2k2
2m2
, (3-12)
which yields
λ2 = − k
m
, (3-13)
or
λ2 = −3k
m
. (3-14)
3-38 Schimel, Ding, Anderson and Grantham
Therefore, there are four possible values of λ that will yield solutions to Eq. (3-5):
λ = ±i
r
k
m
, (3-15)
and
λ = ±i
r
3k
m
. (3-16)
These four values of λ are the eigenvalues of the system defined by Eq. (3-5). For brevity, define
ω1 =
r
k
m
, (3-17)
and
ω2 =
r
3k
m
. (3-18)
Therefore, the eigenvalues become
λ = ±iω1, (3-19)
and
λ = ±iω2. (3-20)
Note that Eqs. (3-17) and (3-18) define the natural frequencies of the two degree-of-freedom system.
Relationships between the constants A1 and A2 can be derived from the eigenvalues and Eq. (3-9).
Performing the vector-matrix multiplication in Eq. (3-9) yields two equations:¡
mλ2 + 2k
¢
A1 − kA2 = 0, (3-21)
and
−kA1 +
¡
mλ2 + 2k
¢
A2 = 0. (3-22)
Substitution of either the positive or negative value of λ from Eq. (3-19) into Eqs. (3-21) and (3-22) gives
the equations
kA1 − kA2 = 0, (3-23)
and
−kA1 + kA2 = 0. (3-24)
Both Eq. (3-23) and (3-24) yield A1 = A2, or in vector form½
A1
A2
¾
= A1
½
1
1
¾
, (3-25)
where A1 is a constant. Similarly, substitution of either the positive or negative value of λ from Eq. (3-20)
gives ½
A1
A2
¾
= A1
½
1
−1
¾
. (3-26)
The vector part of Eqs. (3-25) and (3-26),
{Ψ1} =
½
1
1
¾
, (3-27)
and
{Ψ2} =
½
1
−1
¾
, (3-28)
are the eigenvectors of the system defined by Eq. (3-5). For the above solution approach, each eigenvector is
associated with two eigenvalues, since the eigenvectors were obtained by substituting the eigenvalues of Eqs.
(3-19) and (3-20) into Eqs. (3-21) and (3-22). The magnitude of an eigenvector is arbitrary; for instance,
ME 349 Dynamic Systems Lab Manual 3-39
the multiplier A1 could have been included. However, eigenvectors are often normalized so that their first
component is equal to one.
Observe that the eigenvectors are orthogonal. That is, the eigenvectors possess the properties
{Ψq}T {Ψr} = 0 if q 6= r, (3-29)
and
{Ψq}T {Ψr} 6= 0 if q = r. (3-30)
The eigenvectors can also be represented in matrix form as
[Ψ] = [{Ψq} {Ψr}] =
∙
1 1
1 −1
¸
. (3-31)
The matrix of eigenvectors, [Ψ] , can be used to diagonalize the mass matrix (which is already diagonal) and
the stiffness matrix of Eq. (3-5):
[Ψ]
T
∙
m 0
0 m
¸
[Ψ] =
∙
1 1
1 −1
¸ ∙
m 0
0 m
¸ ∙
1 1
1 −1
¸
=
∙
2m 0
0 2m
¸
, (3-32)
and
[Ψ]T
∙
2k −k
−k 2k
¸
[Ψ] =
∙
1 1
1 −1
¸ ∙
2k −k
−k 2k
¸ ∙
1 1
1 −1
¸
=
∙
2k 0
0 6k
¸
. (3-33)
Also, note that the components of the diagonalized mass and stiffness matrices on the right-hand-side of
Eqs. (3-32) and (3-33) can be used to obtain the natural frequencies:
ω1 =
r
k11
m11
=
r
2k
2m
=
r
k
m
, (3-34)
and
ω2 =
r
k22
m22
=
r
6k
2m
=
r
3k
m
, (3-35)
which are equivalent to Eqs. (3-17) and (3-18).
The above facts concerning the eigenvectors will not be of immediate use in this chapter. However, such
properties are of use when analyzing multi-degree-of-freedom systems, which will be dealt with later in the
text.
Returning to the pursuit of a solution to the equations of motion, observe that Eq. (3-5) is a linear
ordinary differential equation. Therefore, the general solution,½
x1
x2
¾
= C1
½
1
1
¾
ei
√
k
m
t + C2
½
1
1
¾
e−i
√
k
m
t +D1
½
1
−1
¾
ei
√
3k
m
t +D2
½
1
−1
¾
e−i
√
3k
m
t, (3-36)
is a superposition of all four of the possible solutions. Each component of Eq. (3-36) is in the form of
Eq. (3-6), and C1, C2,D1, and D2 are constants which represent the magnitude of each component. By
substituting the natural frequencies as defined in Eqs. (3-17) and (3-18); and collecting the terms with like
eigenvectors, Eq. (3-36) becomes½
x1
x2
¾
=
½
1
1
¾¡
C1eiω1t + C2e−iω1t
¢
+
½
1
−1
¾¡
D1eiω2t +D2e−iω2t
¢
. (3-37)
By the same arguments as applied in Eqs. (1-9) through (1-16), Eq. (3-37) can be rewritten as½
x1
x2
¾
=
½
1
1
¾
C cos(ω1t+ φ1) +
½
1
−1
¾
D cos(ω2t+ φ2), (3-38)
where C, D, φ1, and φ2 are constants.
3-40 Schimel, Ding, Anderson and Grantham
The four constants in Eq. (3-38) can be re-expressed in terms of the initial displacements and velocities,
x1 (0) , x2 (0) , x˙1 (0) , and x˙2 (0) . Writing Eq. (3-38) as two separate displacement equations gives
x1 = C cos(ω1t+ φ1) +D cos(ω2t+ φ2), (3-39)
and
x2 = C cos(ω1t+ φ1)−D cos(ω2t+ φ2). (3-40)
Differentiating Eqs. (3-39) and (3-40) gives the velocities
x˙1 = −Cω1 sin(ω1t+ φ1)−Dω2 sin(ω2t+ φ2), (3-41)
and
x˙2 = −Cω1 sin(ω1t+ φ1) +Dω2 sin(ω2t+ φ2). (3-42)
At t = 0 Eqs. (3-39) through (3-42) become:
x1(0) = C cosφ1 +D cosφ2, (3-43)
x2(0) = C cosφ1 −D cosφ2, (3-44)
x˙1 (0) = −Cω1 sinφ1 −Dω2 sinφ2, (3-45)
and
x˙2 (0) = −Cω1 sinφ1 +Dω2 sinφ2; (3-46)
or, by solving for the functions containing φ1 and φ2,
C cosφ1 =
x1(0) + x2(0)
2
, (3-47)
D cosφ2 =
x1(0)− x2(0)
2
, (3-48)
C sinφ1 = −
x˙1(0) + x˙2(0)
2ω1
, (3-49)
and
D sinφ2 = −
x˙1(0)− x˙2(0)
2ω2
. (3-50)
Applying the trigonometric identity
cos (u+ v) = cosu cos v − sinu sin v (3-51)
to Eqs. (3-39) and (3-40) gives
x1 = C cosφ1 cosω1t− C sinφ1 sinω1t+D cosφ2 cosω2t−D sinφ2 sinω2t, (3-52)
and
x2 = C cosφ1 cosω1t− C sinφ1 sinω1t−D cosφ2 cosω2t+D sinφ2 sinω2t. (3-53)
Finally, substituting Eqs. (3-47) through (3-50) into Eqs. (3-52) and (3-53) gives the solutions to the
equations of motion:
x1 =
³
x˙1(0)+x˙2(0)
2ω1
´
sinω1t+
³
x1(0)+x2(0)
2
´
cosω1t
+
³
x˙1(0)−x˙2(0)
2ω2
´
sinω2t+
³
x1(0)−x2(0)
2
´
cosω2t,
(3-54)
and
x2 =
³
x˙1(0)+x˙2(0)
2ω1
´
sinω1t+
³
x1(0)+x2(0)
2
´
cosω1t
−
³
x˙1(0)−x˙2(0)
2ω2
´
sinω2t−
³
x1(0)−x2(0)
2
´
cosω2t.
(3-55)
ME 349 Dynamic Systems Lab Manual 3-41
3.3 Analysis of the Free-Vibration Response
Equation (3-38) or Eqs. (3-54) and (3-66) show that the displacements, x1 and x2, are the superposition of
two simple harmonic motions, one with frequency ω1, and one with frequency ω2. The frequencies ω1 and
ω2 are the natural frequencies of the system defined by Eq. (3-5). The natural frequencies are equivalent to
the magnitude of the imaginary part of the eigenvalues in Eqs. (3-19) and (3-20) (Note that in the absence
of damping there is no real part). In contrast to the single-degree-of-freedom system, which has only one
natural frequency, the two-degree-of-freedom system has two natural frequencies. Depending on the initial
conditions, the system can vibrate freely with either one or both of the natural frequencies active.
Each natural frequency for a two-degree-of-freedom system has a particular mode of motion associated
with it known as a normal mode of vibration. The normal mode is often represented by a mode shape . The
mode shapes for this system are equivalent to the eigenvectors given in Eqs. (3-27) and (3-28). Consider the
case where the initial displacements, x1(0) and x2(0), are equal to an arbitrary nonzero constant A, and the
initial velocities, x˙1(0) and x˙2(0), are both zero. In this case, Eqs. (3-54) and (3-55) yield
x1 = A cosω1t, (3-56)
and
x2 = A cosω1t; (3-57)
or in matrix form ½
x1
x2
¾
= A
½
1
1
¾
cosω1t. (3-58)
Equation (3-58) shows that for the given set of initial conditions the masses vibrate in unison, that is, at
the same frequency, ω1, and in phase with each other. This is known as the first normal mode of vibration.
Figure 3-2a shows the displacement vs. time traces of the two degree-of-freedom system for this normal
mode.
If the initial displacements are such that x1(0) = A and x2(0) = −A, and the initial velocities remain
zero, Eqs. (3-54) and (3-55) yield
x1 = A cosω2t, (3-59)
and
x2 = −A cosω2t; (3-60)
or in matrix form ½
x1
x2
¾
= A
½
1
−1
¾
cosω2t. (3-61)
In this case the masses vibrate at the second natural frequency, ω2, and out of phase with each other. Figure
3-2b shows the displacement vs. time traces for the second normal mode of vibration. Other cases of initial
conditions, which result in normal mode vibration at frequencies of ω1 or ω2, are also possible.
Figure 3-2c and d show a common method of representing the mode shapes. Figure 3-2c represents the
mode shape associated with ω1, which is often referred to in vibration analysis as the first vibration mode,
or simply as mode one. The nodes in Figure 3-2c are both on the positive side of the axis, and they are
consistent with the two positive components of the eigenvector in Eq. (3-27). Figure 3-2d represents the
mode shape associated with ω2, commonly referred to as mode two. The nodes are on opposite sides of the
axis, consistent with the two opposite-in-sign components of the eigenvector in Eq. (3-28).
Consider a case where the initial conditions which resulted in the above two normal modes of vibration
are combined additively. This gives x1(0) = 2A, x2(0) = 0, x˙1(0) = 0, and x˙2(0) = 0, and the solution
becomes ½
x1
x2
¾
= A
½
1
1
¾
cosω1t+A
½
1
−1
¾
cosω2t. (3-62)
Equation (3-62) is a superposition of the previous two normal mode solutions given in Eqs. (3-58) and
(3-61). Equation (3-62) is also a particular case of general motion for the two-degree-of-freedom system,
where general motion refers to motion in which both normal modes of vibration are active. Many other sets
of initial conditions will also produce general motion. Figure 3-3 shows the displacement vs. time traces of
the general motion for the case given in Eq. (3-62).
3-42 Schimel, Ding, Anderson and Grantham
x1(t)
x1(t)
x2(t)
x2(t)
  t  
  t  
  t  
  t  
(a)
(b)
  1  
(c)
  1  
 -1  
(d)
Figure 3-2: Normal mode response and mode shapes for the two degree of freedom system.
x1(t)
x2(t)
  t  
  t  
Figure 3-3: General motion for the two degree of freedom system.
ME 349 Dynamic Systems Lab Manual 3-43
3.4 The Beating Phenomenon
Consider a two-degree-of-freedom system different from the one defined by Eq. (3-5), in which the two
natural frequencies are very close to each other. For this system, the general solution given in Eqs. (3-54)
and (3-55) is still valid, but the natural frequencies are not the same as defined in Eqs. (3-34) and (3-35).
When the natural frequencies of a two-degree-of-freedom system are very close, and that system is subject
to particular sets of initial conditions, a general motion results which is marked by the masses alternating
gradually between relatively large amplitude motion and a near state of rest. This type of motion is known
as the beating phenomenon, and the theory behind its existence will be shown in the following development.
Consider Eqs. (3-54) and (3-55) with the initial conditions x1 (0) = 12A, x2(0) = 0, x˙1 (0) = 0, and
x˙2 (0) = 0. For these initial conditions the equations become:
x1 =
A
2
cosω1t+
A
2
cosω2t (3-63)
and
x2 =
A
2
cosω1t−
A
2
cosω2t. (3-64)
Applying the trigonometric identities
cosu+ cos v = 2 cos
u+ v
2
cos
u− v
2
(3-65)
and
sinu− sin v = −2 sin u+ v
2
sin
v − u
2
(3-66)
to Eqs. (3-63) and (3-64) gives
x1 = A cos
µ
ω1 + ω2
2
t
¶
cos
µ
ω2 − ω1
2
t
¶
(3-67)
and
x2 = −A sin
µ
ω1 + ω2
2
t
¶
sin
µ
ω2 − ω1
2
t
¶
. (3-68)
Observe that Eqs. (3-67) and (3-68) both contain a sinusoidal function of the frequency sum, ω1+ω22 ,
multiplied by a sinusoidal function of the frequency difference, ω2−ω12 . Since the frequencies ω1 and ω2 are
very close, the frequency sum ω1+ω22 is relatively large when compared to the frequency difference
ω2−ω1
2 .
Therefore, the sinusoidal terms which contain the frequency difference act as slowly varying functions, and
the sinusoidal terms which contain the frequency sum act as quickly varying functions. In general, when a
quickly varying function is multiplied by a slowly varying function, the local behavior looks primarily like
the quickly varying function. Globally, the slowly varying function acts as an envelope to the quickly varying
function. For the case of the beating phenomenon, the slowly varying sinusoidal functions in Eqs. (3-67)
and (3-68) modulate the amplitude of the quickly varying sinusoidal functions.
An example of the beating phenomenon response is shown in Figure 3-4. The high frequency trace is the
response, which is at ω1+ω22 . The curve which envelopes the response is at the low frequency,
ω2−ω1
2 . As can
be seen in the figure, the low frequency controls, or modulates, the amplitude of the high frequency.
3.5 Apparatus
Figure 3-5 shows the two-degree-of-freedom double-pendulum apparatus which will be used in this laboratory.
Each pendulum consists of a slender steel rod-arm connected on one end to a brass bob. The opposite end
of each of the rod-arms connect to separate pivots. The pivots constrain the pendulums to swing on parallel
paths. A single steel band attaches to each pendulum at its pivot and acts as a torsional spring between the
two.
3-44 Schimel, Ding, Anderson and Grantham
x1(t)
x2(t)
  t  
  t  
envelope
Figure 3-4: Beating phenomenon response for the two degree of freedom system.
rod-arm
torsional
 spring pivot
bob
Figure 3-5: The double pendulum apparatus.
ME 349 Dynamic Systems Lab Manual 3-45
mg sin
mg sin
kT(     -     ) kT(     -     )
mass
center
lcm
lcm
(a) (b) (c)
Figure 3-6: Free-body diagram of the double-pendulum apparatus.
Equations of Motion
Figure 3-6a shows an end view of the undamped two-degree-of-freedom pendulum system. Denote that
the angles of rotation associated with the two pendulums as θ1 and θ2. The moment generated by a linear
torsional spring by the first pendulum is given by
M1 = −kT (θ1 − θ2) , (3-69)
where θ1− θ2 is the net angle of rotation in the torsional spring and kT is the torsional spring constant. The
negative sign indicates that the moment, M1, resists any applied rotation.
The free body diagram for the first pendulum is shown in Figure 3-6b. Application of Newton’s second
law gives
−mglcm sin θ1 − kT (θ1 − θ2) = Ipθ¨1, (3-70)
where m is the total mass of the pendulum, lcm is the distance from the pivot to the pendulum’s center of
mass, Ip is the moment of inertia of the pendulum, and θ¨1 is the first pendulum’s angular acceleration. After
rearrangement, Eq. (3-70) becomes
Ipθ¨1 + kT (θ1 − θ2) +mglcm sin θ1 = 0. (3-71)
A similar derivation for the free body diagram shown in Figure 3-6c yields the equation of motion for the
second pendulum:
Ipθ¨2 + kT (θ2 − θ1) +mglcm sin θ2 = 0. (3-72)
Linearization of the Equations of Motion
Equations (3-71) and (3-72) both contain sine terms, and therefore these equations are not linear ordinary
differential equations. If the solution techniques developed in the previous sections are to be used, these
equations must be linearized.
3-46 Schimel, Ding, Anderson and Grantham
To perform the linearization, consider the case when θ1 and θ2 are small angles. For small angles the
following approximation can be used for the sine function:
sinα ≈ α. (3-73)
Application of this approximation to Eqs. (3-71) and (3-72) yields
Ipθ¨1 + kT (θ1 − θ2) +mglcmθ1 = 0 (3-74)
and
Ipθ¨2 + kT (θ2 − θ1) +mglcmθ2 = 0. (3-75)
In matrix form, Eqs. (3-74) and (3-75) are
∙
Ip 0
0 Ip
¸½
θ¨1
θ¨2
¾
+
∙
kT +mglcm −kT
−kT kT +mglcm
¸½
θ1
θ2
¾
=
½
0
0
¾
. (3-76)
Solution to the Linearized Equations of Motion
The form of Eq. (3-76) is equivalent to that of Eq. (3-5). Therefore, the general solution given by Eqs. (3-54)
and (3-55) applies. For the case of the linearized double-pendulum model, θ1 and θ2 replace the variables
x1 and x2 respectively. Applying Eqs. (3-10) through (3-16) gives the natural frequencies
ω1 =
s
mglcm
Ip
, (3-77)
and
ω2 =
s
mglcm + 2kT
Ip
. (3-78)
The general solution to the equation of motion for the linearized double-pendulum model is
θ1 =
Ã
θ˙1(0) + θ˙2(0)
2ω1
!
sinω1t+
µ
θ1(0) + θ2(0)
2
¶
cosω1t
+
Ã
θ˙1(0)− θ˙2(0)
2ω2
!
sinω2t+
µ
θ1(0)− θ2(0)
2
¶
cosω2t,
(3-79)
and
θ2 =
Ã
θ˙1(0) + θ˙2(0)
2ω1
!
sinω1t+
µ
θ1(0) + θ2(0)
2
¶
cosω1t
−
Ã
θ˙1(0)− θ˙2(0)
2ω2
!
sinω2t−
µ
θ1(0)− θ2(0)
2
¶
cosω2t.
(3-80)
The first normal mode of vibration is obtained by applying θ1(0) = Θ◦, θ2(0) = Θ◦, θ˙1(0) = 0, and
θ˙2(0) = 0 as initial conditions. In matrix form, these initial conditions result in½
θ1
θ2
¾
= Θ◦
½
1
1
¾
cosω1t. (3-81)
Similarly, the second normal mode of vibration can be obtained by applying θ1(0) = Θ◦, θ2(0) = −Θ◦,
θ˙1(0) = 0, and θ˙2(0) = 0 as initial conditions which gives½
θ1
θ2
¾
= Θ◦
½
1
−1
¾
cosω2t. (3-82)
ME 349 Dynamic Systems Lab Manual 3-47
The Beating Phenomenon in the Double-Pendulum System
For Eqs. (3-77) and (3-78), consider the case in which kT is small in comparison to mgl . This results in
relatively close values for the natural frequencies ω1 and ω2. For the initial conditions θ1(0) = Θo, θ2(0) = 0,
θ˙1(0) = 0, and θ˙2(0) = 0, Eqs. (3-77) and (3-78) reduce to
θ1 =
Θ◦
2
cosω1t+
Θ◦
2
cosω2t (3-83)
and
θ2 =
Θ◦
2
cosω1t−
Θ◦
2
cosω2t. (3-84)
In matrix form Eqs. (3-83) and (3-84) become½
θ1
θ2
¾
=
Θ◦
2
½
1
1
¾
cosω1t+
Θ◦
2
½
1
−1
¾
cosω2t. (3-85)
Observe that Eqs. (3-83) and (3-84) are equivalent in form to Eqs. (3-63) and (3-64). Therefore, the linearized
equations for the double-pendulum system predict the presence of the beating phenomenon. However,
since the true double-pendulum model is actually nonlinear, the above linear analysis does not assure that
the beating phenomenon will be present. The presence of the beating phenomenon might possibly be
detected through numerical simulation of the nonlinear model, or, as will be done in this laboratory, through
experiment.
3.6 Measurements and Calculations
Measurement of the Torsional Spring Constant
The torsional spring constant, kT , will be used in the following section for computing the theoretical natural
and beating phenomenon frequencies. The torsional spring constant can be easily determined by a static
displacement test. Consider the system as shown in Figure 3-6, for the case when the second pendulum is
suspended at an angle of θ2 = π2 rad. (90
◦). For this case, the spring torsion on the first pendulum is given
by
M1 = −kT
³
θ1 −
π
2
´
. (3-86)
The spring torsion is cancelled by a gravitational moment. The gravitational moment is induced by the
gravitational force located at the pendulums center of mass:
lcm =
mrlr/2 +mb(lr + rb)
m
, (3-87)
where mr is the mass of the rod-arm, lr is the length of the rod-arm, mb is the mass of the bob, rb is the
radius of the bob, and m is the total mass of the pendulum. The gravitational moment is given by
Mg = −mglcm sin θ1. (3-88)
Adding the two moments in Eqs. (3-86) and (3-88), setting the sum equal to zero for the static case, and
then solving for the torsional spring constant, kT , gives
kT =
mglcm sin θ1¡
π
2 − θ1
¢ . (3-89)
Estimation of the Natural Frequencies and Beating Phenomenon Frequencies
The natural frequencies, ω1 and ω2, of the double-pendulum system can be computed from Eqs. (3-77) and
(3-78). With exception to the moment if inertia, Ip, all the variables in these equations have been defined.
3-48 Schimel, Ding, Anderson and Grantham
The moment of inertia is obtained by adding the moment of inertia for a slender rod of length lr rotating
about one end,
Ir =
1
3
mrl
2
r , (3-90)
the moment of inertia for a cylindrical bob of length a and radius rb rotating end-over-end,
Ib =
1
12
mb
¡
3r2b + a
2
¢
, (3-91)
and the parallel axes theorem contribution of the cylindrical bob
Ia = mb (lr + rb)
2 . (3-92)
Therefore, the moment of inertia for a pendulum is given by
Ip =
1
3
mrl
2
r +
1
12
mb
¡
3r2b + a
2
¢
+mb (lr + rb)
2 , (3-93)
where mr and mb are the rod and bob masses respectively.
The beating phenomenon frequencies, which will be denoted as ωL for the low frequency and ωH for the
high frequency are given by
ωL =
ω2 − ω1
2
, (3-94)
and
ωH =
ω1 + ω2
2
. (3-95)
Measurement of the Natural Frequencies
The approximate natural frequencies for the double-pendulum system are readily determined from experi-
ment, by applying a knowledge of the normal mode shapes. The first mode natural frequency, ω1, is obtained
by applying initial rotations in proportion to the first normal mode shape given in Eq. (3-81). This requires
giving both pendulums the same initial displacement and then releasing them at the same instant. By mea-
suring the amount of time it takes for one of the pendulums to complete a set number for cycles, the first
mode natural frequency can be computed from
frequency = 2π
cycle count
elapsed time
rad/sec (3-96)
Similarly, the second mode natural frequency, ω2, is determined by applying the mode shape given in Eq.
(3-82). This requires deflecting the pendulums in opposite in directions, but by the same magnitude, and
then releasing both pendulums at the same instant. Application of the frequency measurement technique
described above, along with Eq. (3-96), yields ω2.
Measurement of the Beating Phenomenon Frequencies
The beating phenomenon frequencies are determined by applying a superposition of the mode shapes given
in Eq. (3-85). This superposition of mode shapes results in one pendulum being displaced by an arbitrary
angle, while the other pendulum is held in a position of zero angular displacement. Releasing both pendulums
at the same instant initiates the beating phenomenon motion.
Both the low frequency and the high frequency associated with the beating phenomenon are measured
during the same test. Both these frequencies are computed as in Eq. (3-96). Measurement of the cycle
counts associated with these frequencies is described in the following paragraph.
The frequency ωL is determined as in the case of the natural frequencies, that is, by measuring the
elapsed time required to complete a set number of cycles. Observe, as shown in Figure 3-7, that it takes two
transitions from low to high amplitude and then back to low, to constitute one cycle at the low frequency,
ωL.
ME 349 Dynamic Systems Lab Manual 3-49
Figure 3-7: Cycle count measurements for the beating phenomenon.
The cycle count for the high frequency, ωH , is determined by counting the total number of cycles which
occur during the elapsed time described in the previous paragraph, that is, the elapsed time measured for
the low frequency cycle count.
Note that the cycles associated with the high frequency may not be clearly visible during transitions
through a minimum amplitude. Also, it is important to observe, as shown in Figure 3-8, that there is a
phase shift as a pendulum moves through the smallest amplitude. Cycle a in the figure is π radians (180◦)
out of phase with cycle c. An explanation for the phase shift is that the low frequency envelope function
changes sign at the minimum amplitude, cycle b in Figure 3-8, which reverses the sign of the peaks of the
high frequency sinusoid.
3.7 Laboratory #3 Procedures
Physical Properties
Put the physical properties of the dual pendulum system in Table 3-1.
Table 3-1. Physical Properties
lr (m)
mr (kg)
rb (m)
a (m)
mb (kg)
lcm (m)
Ir (kgm2)
Ib (kgm2)
Ia (kgm2)
Ip (kgm2)
kT (Nm/rad)
3-50 Schimel, Ding, Anderson and Grantham
x2(t)
  t  
envelope
Figure 3-8: The beating phenomenon response at near a minimum amplitude transition.
Torsional Spring Constant
1. Displace one pendulum π2 radians (90
◦).
2. Measure the angular displacement of the other pendulum and record the angle in Table 3-2.
3. Swap the tasks performed on each pendulum and repeat steps one and two.
Table 3-2. Static Spring Deflections
Trial θ (degrees) θ (radians)
1
2
Natural Frequencies
1. Displace both pendulums by π6 radians (30
◦) in the same direction and simultaneously release.
2. Measure the elapsed time for ten oscillations and record the results in table Table 3-3.
3. Perform five trials of steps one and two.
4. Displace both pendulums by π6 radians (30
◦) but in opposite directions and simultaneously release.
5. Measure the elapsed time for ten oscillations and record the results in table Table 3-4.
6. Perform five trials of steps four and five.
ME 349 Dynamic Systems Lab Manual 3-51
Table 3-3. Period Measurement for ω1
Trial Time for 10 cycles (sec.)
1
2
3
4
5
Ave.
Table 3-4. Period Measurement for ω2
Trial Time for 10 cycles (sec.)
1
2
3
4
5
Ave.
Beating Phenomenon Frequencies
1. Displace one pendulum by π6 radians (30
◦) and hold the other at 0, then release both at the same
instant.
2. Measure the time required for five low frequency cycles and count the total number of high frequency
cycles which occur during this time. Record the results in Table 3-5. Be sure to account for the phase
shift in the high frequency cycles while counting.
3. Perform five trials of steps one and two.
Table 3-5. Beating Phenomenon Measurements
Trial Time for 5 Low Number of
# Freq. Cycles, (sec) Fast Cycles
1
2
3
4
5
Ave.
Chapter 4
Multi-Degree-of-Freedom Vibrations
and Spectral Analysis
4.1 Multi-Degree-of-Freedom Vibrations
Free Vibrations
Equation of Motion, Eigenvalues, and Eigenvectors
Consider a linear N -dimensional mass-spring system. The equation of motion for such a system may be
written as
⎡
⎢⎢⎢⎢⎣
m1 0 · · · 0
0 m2
. . .
...
...
. . .
. . . 0
0 · · · 0 mN
⎤
⎥⎥⎥⎥⎦
⎧
⎪⎪⎨
⎪⎪⎩
x¨1
x¨2
...
x¨N
⎫
⎪⎪⎬
⎪⎪⎭
+
⎡
⎢⎢⎢⎢⎣
k11 k12 · · · k1N
k21 k22
...
...
. . .
...
kN1 · · · · · · kNN
⎤
⎥⎥⎥⎥⎦
⎧
⎪⎪⎨
⎪⎪⎩
x1
x2
...
xN
⎫
⎪⎪⎬
⎪⎪⎭
=
⎧
⎪⎪⎨
⎪⎪⎩
0
0
...
0
⎫
⎪⎪⎬
⎪⎪⎭
, (4-1)
where the mq are masses and kqr, q, r = 1, . . . , N , are stiffnesses, respectively. Equation (4-1) can be written
in a more compact form as
[m] {x¨}+ [k] {x} = {0} , (4-2)
where [m], [c] and [k] are N -by-N matrices, and {x¨}, {x}, and {0} are N -dimensional vectors. Some
example multi-degree-of-freedom systems are the multiple spring and mass system shown in Fig 4-1a, the
lumped-parameter flexible beam model shown in Fig 4-1b, and the two-degree-of-freedom system shown in
Figure 3-1.
Assume that the solution to Eq. (4-1) is of the form
⎧
⎪⎪⎪⎨
⎪⎪⎪⎩
x1(t)
x2(t)
...
xN (t)
⎫
⎪⎪⎪⎬
⎪⎪⎪⎭
=
⎧
⎪⎪⎪⎨
⎪⎪⎪⎩
ψ1
ψ2
...
ψN
⎫
⎪⎪⎪⎬
⎪⎪⎪⎭
eλt; (4-3)
which can also be written as
{x(t)} = {Ψ} eλt, (4-4)
where {x(t)} is an N -dimensional vector, and {Ψ} is an N -dimensional vector of constants. Substitution of
Eq. (4-4), and its second time derivative, into Eq. (4-2) yields¡
[m]λ2 + [k]
¢ {Ψ}X◦eλt = 0. (4-5)
Since there are no values of λ for which eλt = 0 for all values of t, it must be that¡
[m]λ2 + [k]
¢ {Ψ} = 0 (4-6)
4-52
ME 349 Dynamic Systems Lab Manual 4-53
E,I,L,m
Figure 4-1: Example multi-degree of freedom systems.
Nontrivial solutions to Eq. (4-6) require that the following determinant is equal to zero:¯¯
[m]λ2 + [k]
¯¯
= 0, (4-7)
which yields the polynomial
d2Nλ
2N + d2(N−1)λ
2(N−1) + · · ·+ d2λ2 + d0 = 0, (4-8)
with d2r, r = 1, . . . ,N as constant coefficients. Solving Eq. (4-8) for λ results in N complex conjugate pairs,
each pair in the form
λr = iωr, r = 1, . . . , N, (4-9)
and
λ∗r = −iωr, r = 1, . . . , N. (4-10)
Each complex conjugate pair of eigenvalues constitute two out of 2N possible eigenvalues, all of which satisfy
Eq. (4-5).
For the multi-degree-of-freedom mass-spring model in Eq. (4-1) there exist N independent eigenvectors
of the form {Ψr}. Each eigenvector, {Ψr}, is obtained by substituting ether λr or λ∗r into Eq. (4-6). The
resulting system of N equations and N unknowns is then solved for the components of the rth eigenvector.
Each eigenvector, {Ψr} , is real valued, because λ2r is real valued in Eq. (4-6).
The eigenvectors are typically presented in the following N -by N matrix form:
[Ψ] = [{Ψ1} {Ψ2} . . . {ΨN}] , (4-11)
Correspondingly, the eigenvalues can also be represented in vector form as
{Λ} =
⎧
⎪⎪⎪⎨
⎪⎪⎩
λ1
λ2
...
λN
⎫
⎪⎪⎪⎬
⎪⎪⎭
(4-12)
4-54 Schimel, Ding, Anderson and Grantham
and
{Λ∗} =
⎧
⎪⎪⎨
⎪⎪⎩
λ∗1
λ∗2
...
λ∗N
⎫
⎪⎪⎬
⎪⎪⎭
. (4-13)
The Eigenproblem
Equation (4-6) is known as an eigenproblem, and any eigenvalue-eigenvector pair, (λr, {Ψr}) or (λ∗r , {Ψr}) ,
which satisfies Eq. (4-6) is known as an eigensolution. For N -dimensional equations in the form of Eq. (4-6)
there are 2N possible eigensolutions: N solutions of the form (λr, {Ψr}) , and another N solutions of the
form (λ∗r , {Ψr}).
The eigensolution and the eigenproblem have some important properties. Consider Eq. (4-6) and sub-
stitute the eigensolution (λr, {Ψr}), which gives¡
[m]λ2r + [k]
¢ {Ψr} = {0} . (4-14)
Premultiplying both sides of Eq. (4-14) by the transpose of the eigenvector {Ψr}T gives
{Ψr}T
¡
[m]λ2r + [k]
¢ {Ψr} = {Ψr}T {0} . (4-15)
Although λr is complex valued, it is still a scalar quantity. Therefore, Eq. (4-15) can also be written as
{Ψr}T [m] {Ψr}λ2r + {Ψr}T [k] {Ψr} = 0 (4-16)
A vector-matrix-vector multiplication, such as {Ψr}T [m] {Ψr}, always results in a scalar quantity. There-
fore, assign
Mr = {Ψr}T [m] {Ψr} (4-17)
and
Kr = {Ψr}T [k] {Ψr} . (4-18)
The scalar constantsMr, andKr are associated with the rth mode of vibration of the multi-degree-of-freedom
system. Therefore, Mr is referred to as the modal mass constant, and Kr is referred to as the modal stiffness
constant for the rth mode. The modal mass and stiffness constants can also be represented in matrix form,
by assigning
[M] = [Ψ]
T
[m] [Ψ] , (4-19)
and
[K] = [Ψ]
T
[k] [Ψ] , (4-20)
where [M], and [K] are diagonal matrices consisting of values of Mr and Kr respectively.
Substitution of the rth modal mass and stiffness constants into Eq. (4-16) gives
Mr(λr)
2 +Kr = 0. (4-21)
Equivalent, the substitution of the eigensolution (λ∗r , {Ψr}) instead of λr into Eq. (4-6) would have yielded
Mr(λ
∗
r)
2 +Kr = 0 (4-22)
Equations (4-21) or (4-22) can be used to obtain an expression for the rth natural frequency, ωr. Substituting
the eigenvalue λr = iωr into Eq. (4-21) yields
−Mrω2r +Kr = 0. (4-23)
Then, solving for ωr gives
ωr =
r
Kr
Mr
. (4-24)
Observe that the form of the equation for the rth mode natural frequency, Eq. (4-24), is equivalent to the
form of its single-degree-of-freedom counterpart given in Eq. (1-7).
ME 349 Dynamic Systems Lab Manual 4-55
Solution to the Equation of Motion
Since Eq. (4-6) is a system of linear ordinary differential equations, the general solution can be written as the
superposition of all possible solutions. Each solution is in the form of Eq. (4-6), and there is one solution for
each eigenvector-eigenvalue pair. Therefore, the general solution is a linear combination of the 2N possible
solutions, containing the (λr, {Ψr}) and (λ∗r , {Ψr}) . The general solution can be written as
{x(t)} = A1 {Ψ1} eλ1t +B1 {Ψ1} eλ∗1t + . . .+AN {ΨN} eλN t +BN {ΨN} eλ∗N t, (4-25)
or, by substituting the definitions of λr and λ
∗
r and collecting like terms,
{x(t)} =
NX
r=1
{Ψr}
¡
Areiωrt +Bre−iωrt
¢
; (4-26)
where Ar and Br are constants. By applying the same argument as in Eqs. (1-9) through (1-16), Eq. (4-26)
can be written as
{x(t)} =
NX
r=1
Dr {Ψr} cos (ωrt+ φr) . (4-27)
Observe that Eq. (4-27) is simply an N -degree-of-freedom extension of the solution for free vibration of the
two-degree-of-freedom system, Eq. (3-38).
Three-Degree-of-Freedom-Example
Consider the three-degree-of-freedom system shown in Figure 4-2a. Analysis of the free-body diagrams in
Figure 4-2b gives the following equations of motion:
−kx1 − k (x1 − x2) = mx¨1, (4-28)
−k (x2−x1)− k (x2 − x3) = mx¨2, (4-29)
and
−k (x3 − x2)− kx3 = mx¨3. (4-30)
Following rearrangement and conversion to matrix form, Eqs. (4-28) through (4-30) become
⎡
⎣
m 0 0
0 m 0
0 0 m
⎤
⎦
⎧
⎨
⎩
x¨1
x¨2
x¨3
⎫
⎬
⎭+
⎡
⎣
2k −k 0
−k 2k −k
0 −k 2k
⎤
⎦
⎧
⎨
⎩
x1
x2
x3
⎫
⎬
⎭ =
⎧
⎨
⎩
0
0
0
⎫
⎬
⎭ (4-31)
Consider the case when m = 1kg, and k = 1N/m. The eigenproblem to solve for this system of equations is
in the form of Eq. (4-6), and is given by
⎡
⎣
⎡
⎣
1 0 0
0 1 0
0 0 1
⎤
⎦λ2 +
⎡
⎣
2 −1 0
−1 2 −1
0 −1 2
⎤
⎦
⎤
⎦
⎧
⎨
⎩
ψ1
ψ2
ψ3
⎫
⎬
⎭ =
⎧
⎨
⎩
0
0
0
⎫
⎬
⎭ (4-32)
Setting the determinant equal to zero, as in Eq. (4-7), and solving the resulting sixth-degree polynomial for
roots λ gives:
{Λ} = − {Λ∗} =
⎧
⎪⎨
⎪⎩
i
p
2−
√
2
i
√
2
i
p
2 +
√
2
⎫
⎪⎬
⎪⎭
. (4-33)
Therefore, the natural frequencies are
⎧
⎨
⎩
ω1
ω2
ω3
⎫
⎬
⎭ =
⎧
⎪⎨
⎪⎩
p
2−
√
2√
2p
2 +
√
2
⎫
⎪⎬
⎪⎭
rad/sec. (4-34)
4-56 Schimel, Ding, Anderson and Grantham
Figure 4-2: A three-degree-of-freedom system.
Inserting values of λr or λ
∗
r for each r = 1, . . . , N into Eq. (4-32) and solving for the respective eigenvector
gives
{Ψ1} =
⎧
⎨
⎩
1√
2
1
⎫
⎬
⎭ {Ψ2} =
⎧
⎨
⎩
1
0
−1
⎫
⎬
⎭ {Ψ3} =
⎧
⎨
⎩
1
−
√
2
1
⎫
⎬
⎭ . (4-35)
In matrix form Eq. (4-35) becomes
[Ψ] =
⎡
⎣
1 1 1√
2 0 −
√
2
1 −1 1
⎤
⎦ . (4-36)
The modal mass matrix is
[M] = [Ψ]
T
[m] [Ψ] =
⎡
⎣
1
√
2 1
1 0 −1
1 −
√
2 1
⎤
⎦
⎡
⎣
1 0 0
0 1 0
0 0 1
⎤
⎦
⎡
⎣
1 1 1√
2 0 −
√
2
1 −1 1
⎤
⎦ =
⎡
⎣
4 0 0
0 2 0
0 0 4
⎤
⎦ kg, (4-37)
and similarly, the modal stiffness matrix is
[K] = [Ψ]
T
[k] [Ψ] =
⎡
⎣
8− 4
√
2 0 0
0 4 0
0 0 8 + 4
√
2
⎤
⎦N/m. (4-38)
Calculating the natural frequencies as in Eq. (4-24) gives
⎧
⎨
⎩
ω1
ω2
ω3
⎫
⎬
⎭ =
⎧
⎪⎨
⎪⎩
p
2−
√
2√
2p
2 +
√
2
⎫
⎪⎬
⎪⎭
rad/sec =
⎧
⎨
⎩
0.765
1.414
1.848
⎫
⎬
⎭ rad/sec, (4-39)
ME 349 Dynamic Systems Lab Manual 4-57
which is the same as in Eqs. (4-34). Note that obtaining the natural frequencies by Eqs. (4-37) through
(4-39) serves only as a check on calculations, since the natural frequencies are also available in Eq. (4-34)
However, if viscous damping had been included in the model, the eigenvalues would have been in the form
λr = −ζrωr + iωr
q
1− ζ2r, (4-40)
and the most direct way to obtain the natural frequencies would have been through Eqs. (4-24).
Finally, the solution to the equation of motion for the freely vibrating system can be written as
{x(t)} =
⎧
⎨
⎩
x1(t)
x2(t)
x3(t)
⎫
⎬
⎭ = D1
⎧
⎨
⎩
1√
2
1
⎫
⎬
⎭ cos
³³p
2−
√
2
´
t− φ1
´
+ D2
⎧
⎨
⎩
1
0
1
⎫
⎬
⎭ cos
¡√
2t− φ2
¢
+ D3
⎧
⎨
⎩
1
−
√
2
1
⎫
⎬
⎭ cos
³³p
2 +
√
2
´
t− φ3
´
,
(4-41)
where D1 through D3 and φ1 through φ3 are determined through initial conditions.
Forced Vibrations
Equation of Motion and the Transfer Function
If a general set of force inputs are applied to the multi-degree-of-freedom system, the equation of motion
becomes
⎡
⎢⎢⎢⎢⎣
m1 0 · · · 0
0 m2
. . .
...
...
. . .
. . . 0
0 · · · 0 mN
⎤
⎥⎥⎥⎥⎦
⎧
⎪⎪⎨
⎪⎪⎪⎩
x¨1
x¨2
...
x¨N
⎫
⎪⎪⎬
⎪⎪⎪⎭
+
⎡
⎢⎢⎢⎢⎣
k11 k12 · · · k1N
k21 k22
...
...
. . .
...
kN1 · · · · · · kNN
⎤
⎥⎥⎥⎥⎦
⎧
⎪⎪⎨
⎪⎪⎪⎩
x1
x2
...
xN
⎫
⎪⎪⎬
⎪⎪⎪⎭
=
⎧
⎪⎪⎨
⎪⎪⎪⎩
f1
f2
...
fN
⎫
⎪⎪⎬
⎪⎪⎪⎭
; (4-42)
or in compact form
[m] {x¨}+ [k] {x} = {f} , (4-43)
where {f} is an N -dimensional vector of forcing functions, and the other components are as defined for Eq.
(4-2). Taking the Laplace transform of Eq. (4-43) yields¡
[m] s2+ [k]
¢ {X (s)} = {F (s)} , (4-44)
and solving Eq. (4-44) for {X (s)} gives
{X (s)} = ¡[m] s2+ [k]¢−1 {F (s)} . (4-45)
The N -by-N
[G(s)] =
¡
[k] s2 + [m]
¢−1
(4-46)
is known as the transfer matrix of the multi-degree-of-freedom system. Element Gpq(s) of [G(s)] is the
transfer function between the pth output, xp, and qth input, fq.
4-58 Schimel, Ding, Anderson and Grantham
The Forced System Frequency Response Function
Recall, from Chapter 2, that the frequency response function is a special case of the transfer function that
represents the relationship between a steady-state sinusoidal input and response at an arbitrary forcing
frequency Ω. As was noted for Eq. (2-42), the frequency response function can be obtained from the transfer
function by substituting iΩ for each occurrence of s. The matrix form of the frequency response functions
is obtained by substituting iΩ for each occurrence of s in Eq. (4-46). Performing the substitution gives
[G(iΩ)] =
¡
[k]− Ω2 [m]
¢−1
. (4-47)
This form of the frequency response function is not particularly useful; it requires the inversion of the right-
hand-side of the equation for each frequency of interest, and it offers little insight into the effects that the
value of the frequency has on the response. A more useful equation for the frequency response function will
be derived below.
Consider the inverse of Eq. (4-47),
[G(iΩ)]−1 =
¡
[k]− Ω2 [m]
¢
, (4-48)
and a scaled version of the eigenvectors known as mass normalized eigenvectors,
[Φ] = [{Φ1} , {Φ2} , . . . , {ΦN}] =
∙
1√
M1
{Ψ1} 1√
M2
{Ψ2} , . . . , 1√
MN
{ΨN}
¸
. (4-49)
Observe that
{Φr}T [m] {Φr} = 1√
Mr
{Ψr}T [m] 1√
Mr
{Ψr} = MrMr = 1 (4-50)
and
{Φr}T [k] {Φr} = 1√
Mr
{Ψr}T [k] 1√
Mr
{Ψr} = KrMr = ω
2
r. (4-51)
Therefore,
[Φ]
T
[m] [Φ] = [I] (4-52)
and
[Φ]
T
[k] [Φ] = [ω]2 =
⎡
⎢⎢⎢⎢⎣
(ω1)
2 0 · · · 0
0 (ω2)
2 . . .
...
...
. . .
. . . 0
0 · · · 0 (ωN )2
⎤
⎥⎥⎥⎥⎦
. (4-53)
Postmultiplying both sides of Eq. (4-48) by [Φ] and premultiplying by [Φ]T gives
[Φ]
T
[G(iΩ)]−1 [Φ] = [Φ]T
¡
[k]− Ω2 [m]
¢
[Φ] =
³
[ω]2 − Ω2 [I]
´
. (4-54)
Assigning
[Ω]2 =
⎡
⎢⎢⎢⎢⎣
Ω2 0 · · · 0
0 Ω2 . . .
...
...
. . .
. . . 0
0 · · · 0 Ω2
⎤
⎥⎥⎥⎥⎦
(4-55)
gives
[Φ]T [G(iΩ)]−1 [Φ] =
³
[ω]2 − [Ω]2
´
. (4-56)
Premultiplying both sides of Eq. (4-56) by
³
[ω]2 − [Ω]2
´−1
and postmultiplying both sides by [Φ]−1 results
in ³
[ω]2 − [Ω]2
´−1
[Φ]T [G(iΩ)]−1 = [Φ]−1 . (4-57)
ME 349 Dynamic Systems Lab Manual 4-59
Premultiplying both sides of Eq. (4-57) by [Φ] , postmultiplying by [G(iΩ)] , and exchanging right- and
left-hand sides yields the frequency response function,
[G(iΩ)] = [Φ]
³
[ω]2 − [Ω]2
´−1
[Φ]T . (4-58)
Since
³
[ω]2 − [Ω]2
´
is diagonal,
³
[ω]2 − [Ω]2
´−1
=
⎡
⎢⎢⎢⎢⎢⎣
1
ω 21 −Ω2
0 · · · 0
0 1ω 22 −Ω2
. . .
...
...
. . .
. . . 0
0 · · · 0 1ω2
N
−Ω2
⎤
⎥⎥⎥⎥⎥⎦
, (4-59)
and therefore, an individual element of [G(iΩ)] can be written as
Gpq(iΩ) =
NX
r=1
1
Mr
ψprψqr
ω 2r − Ω2
, (4-60)
where ψpr is the p
th and ψqr is the q
th element of the rth eigenvector. The element Gpq(iΩ) is the frequency
response function relating the qth sinusoidal input, Fq (iΩ) , to the pth sinusoidal response Xp (iΩ) . It is
equivalent to a ratio of the response amplitude over the input amplitude. Observe that when the driving
frequency is near a given natural frequency, ωr, the magnitude of the frequency response for that mode,
1
Mr
ψprψqr
ω 2r −Ω2
, will dominate the overall magnitude of the frequency response function, Gpq(iΩ).
The frequency response function of Eq. (4-60) is for an undamped system. If viscous damping is included
in the model, the frequency response function is of the form
Gpq(iΩ) =
NX
r=1
Apqr
(ω 2r − Ω2) + 2iΩωrζr
, (4-61)
where Apqr is a complex-valued quantity. In the presence of viscous damping, the magnitude of Gqr(iΩ) is
of the form
|Gqr(iΩ)| =
NX
r=1
Bpqrq
(ω 2r − Ω2)
2
+ (2Ωωrζr)
2
. (4-62)
Observe that the denominator of Eq. (4-62) is equivalent in form to its single-degree-of-freedom counterpart
in Eq. (2-50). Figure 4-3 shows an example frequency response function magnitude plot for a three-degree-
of-freedom system with viscous damping. The figure plots the frequency response function magnitude as a
function of driving frequency, Ω. Each natural frequency appears as a peak on the magnitude plot. Note
that if a natural frequency is well separated from other modes the response in the neighborhood of that
frequency approximates a single-degree-of-freedom system.
4.2 Spectral Analysis
The general motion of a linear multi-degree-of-freedom system can be much more complex than that of a
one- or two-degree-of-freedom system. The presence of multiple natural frequencies in a system can preclude
the application of any simple analysis technique which is based on the time-domain response traces. In cases
where multi-degree-of-freedom behavior is present, spectral analysis techniques can be applied to recover the
response at individual frequencies.
Spectral analysis consists of converting time-domain data into frequency-domain data. For instance,
time-domain vibration data in the form of displacement, velocity, or acceleration verses time, etc., can
converted into the frequency-domain, which consists of frequency, amplitude, and phase angle information.
The transformation into the frequency domain is known as a Fourier transform.
4-60 Schimel, Ding, Anderson and Grantham
Driving Frequency,
|G qr(    )|
Figure 4-3: Example frequency response function magnitude plot for a three-degree-of-freedom system with
viscous damping.
A modern approach to spectral analysis utilizes computers to perform the data collection and subsequent
transformations into the frequency domain. The data collection of time-domain signals, which are produced
by transducers, is usually done by sampling the data at discrete time intervals. The transformation into
the frequency domain is done with software that is designed to perform a Fourier transform on discretely
sampled data; it is known as a discrete Fourier transform or DFT .
Transducers and discrete data sampling, as well as the DFT and its applications to multi-degree-of-
freedom vibrations, will be described in more detail in the following sections.
Discrete-Time-Measurements
A computerized data sampling system, commonly known as a data acquisition system, consists of several
primary devices. A transducer, which interfaces with the test system for making physical measurements, a
device that converts the output of the transducer into numeric data which can be stored on a computer, and
software that controls collection and storage of the data.
Transducers
Physical measurements on a system can be made with transducers. A type of transducer that is commonly
used in vibration sensing is the accelerometer, which as the name suggests is used to measure acceleration.
Other types of transducers commonly used in vibration studies are load cells, used for measuring forces, and
linear variable differential transducers, used for measuring displacements. Accelerometers will be used in
this laboratory, and therefore they will be described in more detail here.
A schematic of an accelerometer is shown in Figure 4-4. The accelerometer consists of two main parts:
a piezoelectric crystal structure and a small mass. Piezoelectric crystals generate a charge proportional to
the force applied to them. One side of the crystal structure is mounted to a rigid base in the body of the
transducer. The mass is attached to, and is completely supported by, the other side of the crystal structure.
The base of the transducer is attached to the test system. Piezoelectric crystals are very stiff. Therefore,
ME 349 Dynamic Systems Lab Manual 4-61
housing
mass
piezoelectric
crystalsbase
leads
apparatus
Figure 4-4: Schematic of an accelerometer.
the acceleration experienced by the mass attached to the crystals is essentially the same as the acceleration
experienced by the base of the transducer. Acceleration in the direction of the axis on which the transducer
base, piezoelectric crystals, and mass lie will generate an inertia force in the mass which is countered by
tension or compression of the crystal structure. By Newton’s second law, the inertia force is proportional to
the acceleration experienced by the mass and, hence, the charge produced by the crystals is proportional to
the acceleration. The charge is measured, amplified, and converted to a voltage proportional to the charge
by a device which is generally known as a signal conditioner. The output of the signal conditioner is suitable
as input to the data acquisition system.
Discrete-Time Data Sampling
Consider the theoretical output of a transducer shown in Figure 4-5a. The output is in the form of a voltage,
such as the output of a signal conditioned accelerometer, which varies with time in a continuous manner. A
typical discrete time data sampling system reads the voltage output of the transducer at discrete intervals
of time as shown in Fig 4-5b. The most common case is to use time intervals of constant length. The time
interval, which is commonly called the sampling interval, will be denoted here as τ s.When stating the speed
of the data collection, it is common to state the inverse of the sampling interval as a frequency,
fs =
1
τs
, (4-63)
usually in cycles per second or Hertz. The total number of samples in the data set is of importance when
computing the discrete Fourier transform or DFT . The number of samples will be denoted as Ns.
At each sampling interval the transducer output voltage is converted to numeric data by a device known
as an analog-to-digital converter or A/D converter. The numeric data is in the form of integers, and the value
of the integer stored is proportional to the transducer voltage. Ten-bit converters, which have a numeric
range of zero to 210− 1 = 4095, are common. If, for example, the transducer output ranges from zero to ten
volts, then five volts would be stored numerically as approximately 2048. For calculation or other purposes,
the numeric data can be easily converted back into a voltage or physical quantity by a simple conversion
factor.
4-62 Schimel, Ding, Anderson and Grantham
x(n)
n
t
x(t)
(a)
(b)
Figure 4-5: Example continious-time and sampled data.
The resolution of the numeric data is dependent on the bit-resolution of the A/D converter. If the
transducer output ranges from zero to ten volts, for example, then for a 10-bit converter the voltage resolution
is 10 volts/4096 ≈ 0.0025 volts. Higher resolution converters are available. However, resolution on this order
poses few if any problems when computing the DFT .
The Discrete Fourier Transform
How the DFT Works
The primary objective of the discrete Fourier transform, or DFT, is to extract amplitude, frequency, and
phase angle information from a sampled data set which is periodic. To understand how the DFT works,
some simple theoretical examples will be presented in this section. The first example will consist of applying
the DFT to a simple sine function, and it will be shown that the output of the DFT generates the amplitude
and frequency of that function.
Consider the output, x (t) , of a theoretical transducer which is in the form of the sine function: x(t) =
A sin (ωt) . Assume that x(t) is sampled at a sampling interval of length τs. The value of x(t) at start of the
nth sampling interval is given by
x (n) = A sin (ωnτs) , (4-64)
where the integer n ranges from 0 to Ns − 1.
It will be useful to convert the right-hand side of Eq. (4-64) into a form that offers insight into the
DFT . To achieve this, assume that the length of the test, Nsτs, is an integer multiple of the period of x(t).
Therefore,
l =
Nsτ sω
2π
, (4-65)
where l is an integer and the period of x(t) is given by 2πω . The assumption of the integer relationship
between period and test length is useful for the theoretical analysis. It may be relaxed for the analysis of
ME 349 Dynamic Systems Lab Manual 4-63
actual data sets. Solving Eq. (4-65) for ω gives
ω = ωl =
2πl
Nsτ s
. (4-66)
Substituting this representation of ωl into Eq. (4-64) yields
x(n) = A sin
µ
2πln
Ns
¶
. (4-67)
Before presenting the discrete Fourier transform, some simple relationships for the theoretical sampled
sine function in Eq. (4-67) will be presented. First, observe that
1
Ns
Ns−1X
n=0
A sin
µ
2πln
Ns
¶
= 0. (4-68)
The sum is equal to zero due to the previous assumption that the test length is an integer multiple of the
period of x (n) . This means that the sum of the positive points is canceled by the sum of the negative data
points. The example in Figure 4-6 shows that point x (1) is exactly cancelled by point x (Ns − 1), or in
general point x (n) is canceled by point x (Ns − n) . The same conclusion can be reached for higher multiples
of the same frequency. That is,
2
Ns
Ns−1X
n=0
A sin
µ
2qπln
Ns
¶
= 0, (4-69)
where q is an integer and |q| ≥ 1. By a similar approach, it can be shown that
2
Ns
Ns−1X
n=0
A cos
µ
2qπln
Ns
¶
= 0. (4-70)
The usefulness of Eqs. (4-69) and (4-70) will be seen in the following development.
Consider the multiplication of an element of Eq. (4-68) by the element of another sine function with
frequency
ωj =
2πj
Nsτ s
. (4-71)
An element of this multiplication is given by
2A
Ns
sin
µ
2πln
Ns
¶
sin
µ
2πjn
Ns
¶
, (4-72)
and the summation over all the elements is
2
Ns
Ns−1X
n=0
A sin
µ
2πln
Ns
¶
sin
µ
2πjn
Ns
¶
. (4-73)
It will be shown next that the sum in Eq. (4-73) is equal to the amplitude A. Applying the trigonometric
relationship
sinu sin v =
1
2
(cos (u− v)− cos (u+ v)) (4-74)
to Eq. (4-73) gives
A
Ns
Ns−1X
n=0
µ
cos
µ
2π (l − j)n
Ns
¶
− cos
µ
2π (l + j)n
Ns
¶¶
. (4-75)
If j 6= l then the frequency of the two sine functions in Eq. (4-73) are unequal, and Eq. (4-75) can be
rewritten as
A
Ns
Ns−1X
n=0
µ
cos
µ
2πpn
Ns
¶
− cos
µ
2πqn
Ns
¶¶
= 0, (4-76)
4-64 Schimel, Ding, Anderson and Grantham
  1    2    3    0 
Ns  4    5    ...    ...  
 
n 
x(n)
Figure 4-6: An example revealing the zero sum of the elements of x(n).
where, in this equation, p = l − j and q = l + j are both nonzero. The sum is equal to zero by Eq. (4-70)
since both p and q are integers. If j = l then the frequency of the two sine functions in Eq. (4-73) are equal
and Eq. (4-75) becomes
A
Ns
Ns−1X
n=0
µ
cos (0)− cos
µ
4πln
Ns
¶¶
, (4-77)
which, by Eq. (4-70), reduces to
A
Ns
Ns−1X
n=0
cos (0) =
ANs
Ns
= A. (4-78)
Thus, if the summation of the discrete sine function x (n) is multiplied on an element-by-element basis by a
sine function of equivalent frequency and amplitude one, the amplitude of x(t) can be recovered. As will be
seen, the approach of multiplying a discrete function x (n) by a discrete sinusoidal function in this fashion is
central to the DFT method. Some other considerations will be made first however.
Consider multiplying each element of Eq. (4-68) by the elements of a cosine function. For this case the
sum becomes
2
Ns
Ns−1X
n=0
A sin
µ
2πln
Ns
¶
cos
µ
2πjn
Ns
¶
. (4-79)
By applying the trigonometric relationship
sinu cos v =
1
2
(sin (u+ v) + sin (u− v)) , (4-80)
Eq. (4-79) may be written as
A
Ns
Ns−1X
n=0
µ
sin
µ
2π (l − j)n
Ns
¶
− sin
µ
2π (l + j)n
Ns
¶¶
. (4-81)
ME 349 Dynamic Systems Lab Manual 4-65
Application of Eq. (4-69) would reveal that Eq. (4-81) is equal to zero for any value of the integer j˙.
Therefore, the summation of the element wise multiplication of a sine function by the cosine function is
always equal to zero under the above assumptions.
Similarly, if x (t) is a cosine function
x(t) = A cos (ωt) , (4-82)
its discretely sampled counterpart can be written as
x(n) = A cos
µ
2πln
Ns
¶
. (4-83)
By a development similar to that in Eqs. (4-73) through (4-81) it can be shown that the summation
2
Ns
Ns−1X
n=0
A cos
µ
2πln
Ns
¶
cos
µ
2πjn
Ns
¶
(4-84)
is equal to A if j = l and is equal to 0 if j 6= l. Thus, if the summation of the discrete cosine function x (n) is
multiplied on an element by element basis by a cosine function of equivalent frequency and amplitude one,
the amplitude of x(t) can be recovered.
Consider a function x (t) that is the sum of a cosine and a sine function of equal frequency:
x (t) = A cos (ωt) +B sin (ωt) (4-85)
As in Eq. (4-67) the discrete sampled version of x (t) can be written as
x(n) = A cos
µ
2πln
Ns
¶
+B sin
µ
2πln
Ns
¶
. (4-86)
Now consider the sum
2
Ns
Ns−1X
n=0
µ
A cos
µ
2πln
Ns
¶
+B sin
µ
2πln
Ns
¶¶µ
cos
µ
2πjn
Ns
¶
− i sin
µ
2πjn
Ns
¶¶
. (4-87)
Denote this sum as X¯ (iωj) where the frequency ωj is equal to
2πj
Nsτs
. Expanding terms gives
X¯ (iωj) = 2Ns
PNs−1
n=0
³
A cos
³
2πln
Ns
´
cos
³
2πjn
Ns
´
+B sin
³
2πln
Ns
´
cos
³
2πjn
Ns
´´
− 2iNs
PNs−1
n=0
³
A cos
³
2πln
Ns
´
sin
³
2πjn
Ns
´
+B sin
³
2πln
Ns
´
sin
³
2πjn
Ns
´´
.
(4-88)
Application of the conclusions drawn for Eq. (4-73), Eq. (4-79) and Eqs. (4-84) to (4-88) yields, for j = l,
X¯ (iωl) = A− iB (4-89)
and X¯ (iωj) = 0 when j 6= l . Observe that the magnitude of X¯ (iωl) is given by¯¯
X¯ (iωl)
¯¯
=
p
A2 +B2, (4-90)
and the phase angle by
∠ X¯ (iωl) = − tan−1 BA. (4-91)
Equations (4-90) and (4-91) give the magnitude and phase of x (t) , as the function is defined in Eq. (4-85).
Therefore, the transformation in Eq. (4-87) preserves the magnitude, phase angle, and frequency of the
time-domain function x (t) .
Now imagine a general function that is given by the sum of several sinusoids:
x (t) =
LX
l=0
(Al cos (ωlt) +Bl sin (ωlt)) . (4-92)
4-66 Schimel, Ding, Anderson and Grantham
There is an upper bound of Ns2 on the summation limit L due to the frequency dependence on l. The reason
for this will be dealt with in a later section under the topic of aliasing. For now just assume that L < Ns2 .
The sampled counterpart to Eq. (4-92) is given by
x(n) =
LX
l=0
µ
Al cos
µ
2πln
Ns
¶
+Bl sin
µ
2πln
Ns
¶¶
. (4-93)
The summation
X¯ (iωq) =
2
Ns
Ns−1X
n=0
LX
l=0
µ
Al cos
µ
2πln
Ns
¶
+Bl sin
µ
2πln
Ns
¶¶µ
cos
µ
2πqn
Ns
¶
− i sin
µ
2πqn
Ns
¶¶
(4-94)
results in
X¯ (iωq) = Aq − iBq, (4-95)
when q = l. The coefficients Aq and Bq in Eq. (4-95) are equivalent to the coefficients Al and Bl in Eq.
(4-93) when q = l. The coefficients are equivalent since the contribution of any terms where l 6= q to X¯ (iωq)
is zero by the argument developed in Eqs. (4-85) to (4-89).
The function defined by Eq. (4-94) can be used to determine the magnitude and phase of any particular
frequency component of a function that is composed of the summation of multiple sinusoids. This fact makes
Eq. (4-94) particularly useful in the analysis of multi-degree-of-freedom systems, because the response of
these systems is often the summation of multiple vibration modes.
It is often the case that the frequency components of the sampled function x (n) are not explicitly known.
In such cases Eq. (4-94), which is a form of the discrete Fourier transform, is used to compute the magnitude,
phase, and frequency of the components. Therefore, the discrete Fourier transform of a discretely sampled
general periodic function x (n) can be defined as:
X¯ (iωq) =
2
Ns
Ns−1X
n=1
x(n)
µ
cos
µ
2πqn
Ns
¶
− i sin
µ
2πqn
Ns
¶¶
, q = 0, . . . ,
Ns
2
. (4-96)
The range of the index q will be discussed below. Various forms of Eq. (4-96) appear in the literature. The
leading coefficient 2Ns is sometimes changed or dropped altogether. The value
2
Ns
was used here because it
properly scales the amplitude components in Eq. (4-95). Also, the sign of the sine term in Eq. (4-96) is
sometimes switched. These variations tend to make little difference in the practical use of the DFT , since
the resulting frequencies, phase differences between signals, and ratios of their amplitudes are of primary
importance and the computation of these values tends to cancel the variations in form. The DFT also often
appears in complex form
X¯ (iωq) =
2
Ns
Ns−1X
n=1
x(n)e−i(
2πqn
Ns
), q = 0, . . . ,
Ns
2
, (4-97)
which is equivalent to Eq. (4-96).
The Frequency Spectrum and Spectral Leakage
The output of the discrete Fourier transform, X¯ (iωq), at a particular frequency is referred to as a spectral
estimate, and the entire set of X¯ (iωq) , q = 0, . . . , Ns2 , is referred to as the spectrum or frequency spectrum
of the sampled function x (n) . Magnitude plots of some typical spectrums, along with their time-domain
functions, appear in Figure 4-7. Data at a single frequency, ωq, will produce a magnitude plot with a
single spike at ωq, as in Figure 4-7a. Multiple frequency data will produce a magnitude plot with spikes
corresponding to each frequency, as in Figure 4-7b.
If the data is in the form of a pure sinusoid with frequency ω, but ω is intermediate to one of the
frequencies given by
ωq =
2πq
Nsτs
, (4-98)
ME 349 Dynamic Systems Lab Manual 4-67
(a)
(b)
(c)
(d)
Figure 4-7: Continious time signals and their spectral estimates. (a) x(t) = sinωt. (b) x(t) = sinωt+sin 2ωt.
(c) x (t) = e−ζωt cosωdt. (d) x (t) = sinωt; ω 6= 2πjNsτ , (leakage).
then a phenomenon known as leakage occurs. The spectral estimate of such a sinusoid becomes spread over
the entire spectrum, but the estimates become very small at spectral frequencies far from ω, as in Figure
4-7d. The maximum spectral estimate still occurs at the frequency ωq that is closest to ω. Therefore, an
estimate of the frequency ω can still be obtained. However, the magnitude and phase angle information are
not directly available from the spectrum.
Note that in cases where the time-domain signal is not a pure sinusoid or a sum of pure sinusoids, for
example the function
x (t) = Ae−ζωnt cos (ωdt) , (4-99)
the spectral estimate again becomes spread over the entire spectrum, as in Figure 4-7c. This function is
not truly periodic, but spectral analysis can still be used to estimate the frequency ωd. As with leakage, an
estimate of the frequency ωd is obtained by using the largest spectral estimate.
Frequency Resolution and the Nyquist Frequency
The DFT produces amplitude and phase information at frequencies that are limited in resolution and
maximum value. The frequency resolution is a function of the sampling rate, given by fs as defined in Eq.
(4-63), and the sample count, Ns. The resolution is given by
∆ω = 2πfs
Ns
. (4-100)
The frequency maximum, known as the Nyquist frequency is given by
ωmax =
2πfs
2
= πfs. (4-101)
The Nyquist frequency limits the number of spectral estimates computed by Eq. (4-96) to Ns2 . The reason for
the presence of the Nyquist frequency can be seen in Figure 4-8. Figure 4-8a shows sampling of a sinusoidal
4-68 Schimel, Ding, Anderson and Grantham
  t  
  t  
x(t)
x(t)
(a)
(b)
Figure 4-8: Effects of sampling signals with frequencies at and above the Nyquist frequency.
function that is at the Nyquist frequency; there are two points per cycle. If the frequency is greater than
ωmax, such as in Figure 4-8b, where the frequency is equal to 3ωmax, then many of the oscillations are missed
and the resulting sampled data set looks as if it is at a frequency of ωmax. This problem is known as aliasing,
and it can result in false data appearing at frequencies anywhere in the DFT spectrum. The location in
the spectrum of the aliased frequency data is dependent upon the actual frequency of the time-domain data.
Aliased frequencies can be eliminated by either increasing the sampling frequency, or by pre-filtering the
continuous time signal to remove frequencies above ωmax.
Computing the Frequency Response Function from the DFT
The frequency response function can be looked at as a ratio of the input over the response, where the input
and response are complex numbers representing the magnitude and phase angle of sinusoidal functions. Since
the DFT can be used to compute a complex representation of the magnitude and phase angle of sinusoidal
functions, the DFT 0s of the input and response can be used to compute a value of the frequency response
function at a given frequency. Consider the qth a sinusoidal input to a multi-degree-of-freedom system,
fq (t) = C cosΩt+D sinΩt. (4-102)
For a general linear system, the steady-state time-domain response of the pth mass has the form
xp (t) = A cosΩt+B sinΩt. (4-103)
Equations (4-102) and (4-103) are typical of what would be observed in an experiment on a linear system.
The coefficients A, B, C, and D are constants for a particular value of Ω but their relative values may vary
as Ω varies. In an experimental setting, the input coefficients C and D might be held constant throughout
a series of tests over a frequency range. The constants A and B would be observed to vary as Ω is adjusted.
Consider the Laplace transform of Eqs. (4-102) and (4-103). Since the Laplace transforms
L (cosΩt) = sΩ2 +Ω2 (4-104)
ME 349 Dynamic Systems Lab Manual 4-69
and
L (sinΩt) = Ω
s2 +Ω2 , (4-105)
the Laplace transform of fp (t) is given by
Fp (s) =
Cs+DΩ
s2 +Ω2 . (4-106)
Likewise, for xq (t) ,
Xq (s) =
As+BΩ
s2 +Ω2 . (4-107)
The transfer function for this experimentally observed system is given by
G¯pq (s) =
Xq (s)
Fp (s)
=
As+BΩ
Cs+DΩ . (4-108)
Keep in mind that the experimentally observed transfer function given in Eq. (4-108) differs from the
theoretical transfer function Gpq (s) given in Eq. (4-61). The transfer function G¯pq (s) defines the value that
the transfer function, Gpq (s), of a general linear system takes on if the forcing function and response are
the steady-state sinusoids defined in Eqs. (4-102) and (4-103), respectively. Thus, G¯pq (s) is a particular
instance of Gpq (s), with coefficients A, B, C, and D that may vary as the forcing frequency is varied.
The relationship between input and response is more easily seen if the experimentally observed transfer
function G¯pq (s) is converted into the frequency response function G¯pq (iΩ). Performing the substitution of
s by iΩ in Eq. (4-108) gives
G¯ (iΩ) = X (iΩ)
F (iΩ) =
AiΩ+BΩ
CiΩ+DΩ =
iA+B
iC +D
. (4-109)
Multiplying both the numerator and denominator of Eq. (4-108) by the imaginary number −i gives a more
useful form:
G (iΩ) = A− iB
C − iD . (4-110)
Equation (4-110) is equivalent to the ratio of the discrete Fourier transforms of the sampled versions of the
input, f (t) , and the output, x (t). That is, computing the discrete Fourier transform of the sampled signals
x (n) and f (n) and computing the ratio gives
G¯ (iΩ) = X¯ (iΩ)
F¯ (iΩ)
=
A− iB
C − iD . (4-111)
Therefore, the discrete Fourier transform can be used to obtain an experimental value of the frequency
response function at a particular driving frequency, Ω. Complex algebra can be used to obtain the magnitude
and phase angle, ¯¯
G¯ (iΩ)
¯¯
=
r
A2 +B2
C2 +D2
, (4-112)
and
∠G¯ (iΩ) = tan−1 B
A
− tan−1 D
C
, (4-113)
respectively. As can be seen in Eq. (4-112), the magnitude of G¯ (iΩ) gives the ratio of the amplitude of
the steady-state response function, x (t), over the amplitude of the sinusoidal forcing function, f (t) . The
argument of G¯ (iΩ), Eq. (4-113), is equal to the phase difference between the two sinusoids. Thus, the
frequency response function given in Eq. (4-111) is a compact representation of the magnitude and phase
angle relationship between the input and output of a linear system at a particular forcing frequency.
If a sequence of tests are conducted over a range of discrete forcing frequencies, Ω, then a discrete
version of the frequency response function can be obtained. The process of measuring input and steady-
state response as the frequency is incremented through some desired range is used in practice to obtain the
frequency response function experimentally. The method is known as stepped sine testing.
4-70 Schimel, Ding, Anderson and Grantham
Computer &
  Plotting
   Device 
Spectrum
Analyzer
     Data
Acquistion
   System
Vibration
 System
   Signal
Conditioner
Acceler-
ometers
Vibration
 Exciter
  Power
Amplifier
Sinusoidal
 Function
 Generator
Figure 4-9: Bolck diagram of the apparatus for multi-degree of freedom vibration testing.
Note that if a driving frequency, Ω, is not equal to one of the frequencies output by the DFT , then
leakage will be observed in the DFT ’s of both the input and response. However, this does not generate
problems with the computation of the discrete frequency response function, because the effects of leakage
tend cancel when the ratio is computed.
The Fast Fourier Transform
In practice, the calculation of the DFT is seldom performed by computing Eq. (4-96) directly. Instead a
variant, which produces exactly the same spectral estimates, X¯ (ωj), is used because it is much more efficient
computationally. The advantage becomes more apparent as the number of samples, Ns, becomes large. The
variant is known as the fast Fourier transform, or FFT. The only limitation that the FFT has beyond those
of the DFT is that the sample count, Ns, must be a power of two. This is limitation is not particularly
restrictive for most applications. Therefore, in situations where the entire discrete frequency spectrum is
needed, the FFT is used almost exclusively.
4.3 Apparatus
The Physical System
This experiment will consist of conducting measurements on a multi-degree-of-freedom apparatus using a
computerized data acquisition and spectral analysis system. A block diagram of the entire system is shown
in Figure 4-9.
The vibrating portion of the apparatus will consist of the four-degree-of-freedom model structural system
shown in Figure 4-10. The system consists of four steel masses connected by flexural aluminum members.
The masses model floor and/or roof masses in a four-story structure, and the aluminum members model
support columns and between levels. The base of the structure rests on a slider connected to bearings which
allow horizontal travel. For forced vibration testing, the slider can be connected through a linkage to a
ME 349 Dynamic Systems Lab Manual 4-71
structural
 masses
columns
slider
linkage
vibration
 exciter
model structure
Figure 4-10: The four-degree of freedom model structure and vibration input apparatus.
vibration exciter. The slider can be clamped in position for free-vibration testing. Accelerometers are used
to measure both input and response.
The computerized portion of the apparatus consists of a data acquisition system, a spectrum analyzer,
and a personal computer. The computer is used for initiating data collection and outputting the spectral
analysis results. An example magnitude plot of a spectrum from a free vibration test is shown in Figure 4-11.
Note that the spectral estimates are plotted as a continuum along the horizontal axis instead of as spikes
at each discrete frequency. This is the typical approach, although it is understood that the DFT produces
spectral estimates at only discrete frequencies.
Theoretical Structural Model
Free Vibrations
A small-deflection, free vibration model of the four-degree-of-freedom structural system can be developed as
in Figure 4-12. Figure 4-12a shows the model structure. The coordinates, x1 through x4, are affixed to each
mass, m1 through m4. The mass of the material contained within the dotted rectangles in Figure 4-12a is
assumed to define each mass. Therefore, each mass includes the mass of one steel floor and/or roof member,
plus the mass of half the length of each aluminum column in contact with the steel.
Each of the aluminum columns, which span between two steel masses or a steel mass and the base, is
assigned a stiffness of k2 . Since there are two columns per level, the overall stiffness between each level is
equal to k. Figure 4-12b shows the stiffness model for each column. Each column is assumed equivalent to
a vertically oriented beam in which transverse deflection is allowed, but end rotation is constrained to zero.
Figure 4-12c shows a four-degree-of-freedom spring-mass model equivalent to the structure in Figure
4-12a with the same underlying mass and stiffness assumptions. The direction of motion of the mass-spring
system is along a different axis. However, the models are equivalent. A free-body-diagram analysis part from
4-72 Schimel, Ding, Anderson and Grantham
Figure 4-11: Example output from the spectrum analyzer.
 x4 
k/2=
12E I
(b)
 x2 
 x1 
(a) (c)
Figure 4-12: Free vibration model for the four-degree-of-freedom structural system: (a) structural model,
(b) column stiffness, (c) equivqlent mass-spring model.
ME 349 Dynamic Systems Lab Manual 4-73
Figure 4-12c yields the following equations for motion for the four-degree-of-freedom mass-spring system:
−kx1 − k (x1 − x2) = m1x¨1
−k (x2 − x1)− k (x2 − x3) = m2x¨2
−k (x3 − x2)− k (x3 − x4) = m3x¨3
−k (x4 − x3) = m4x¨4.
(4-114)
Rearrangement of Eq. (4-114) and subsequent conversion to matrix form gives
⎡
⎢⎢⎣
m1 0 0 0
0 m2 0 0
0 0 m3 0
0 0 0 m4
⎤
⎥⎥⎦
⎧
⎪⎨
⎪⎪⎩
x¨1
x¨2
x¨3
x¨4
⎫
⎪⎬
⎪⎪⎭
+
⎡
⎢⎢⎣
2k −k 0 0
−k 2k −k 0
0 −k 2k −k
0 0 −k k
⎤
⎥⎥⎦
⎧
⎪⎨
⎪⎪⎩
x1
x2
x3
x4
⎫
⎪⎬
⎪⎪⎭
=
⎧
⎪⎨
⎪⎪⎩
0
0
0
0
⎫
⎪⎬
⎪⎪⎭
. (4-115)
An experimental and theoretical investigation will yield the following approximate values for the mass
and spring constants: m1,m2,m3 = 1.458 kg; m4 = 1.322 kg; k = 8313. N/m. Substituting these values
into Eq. (4-115) and solving the Eigenproblem as defined in Eqs. (4-4) through (4-11) gives the natural
frequencies ⎧
⎪⎪⎨
⎪⎪⎩
ω1
ω2
ω3
ω4
⎫
⎪⎪⎬
⎪⎪⎭
=
⎧
⎪⎪⎨
⎪⎪⎩
26.77
76.75
116.8
142.3
⎫
⎪⎪⎬
⎪⎪⎭
rad/sec =
⎧
⎪⎪⎨
⎪⎪⎩
4.260
12.22
18.59
22.65
⎫
⎪⎪⎬
⎪⎪⎭
Hz, (4-116)
and their respective mode shapes
{Ψ1} =
⎧
⎪⎪⎨
⎪⎪⎩
1.000
1.874
2.513
2.836
⎫
⎪⎪⎬
⎪⎪⎭
, {Ψ2} =
⎧
⎪⎪⎨
⎪⎪⎩
1.000
0.967
−0.065
−1.030
⎫
⎪⎪⎬
⎪⎪⎭
, {Ψ3} =
⎧
⎪⎪⎨
⎪⎪⎩
1.000
−0.392
−0.846
0.724
⎫
⎪⎪⎬
⎪⎪⎭
, {Ψ4} =
⎧
⎪⎪⎨
⎪⎪⎩
1.000
−1.552
1.408
−0.634
⎫
⎪⎪⎬
⎪⎪⎭
, .
(4-117)
Forced Vibrations
For this laboratory a vibration exciter will be used to drive the four-degree-of-freedom model structure along
a horizontal axis at its base. This input type is referred to by structural engineers as base excitation and
the input location is typical of that seen when buildings are subjected to an earthquake. The base input for
a model structure can be in the form of a force, displacement, velocity, or acceleration. In this laboratory,
accelerometers will be used to measure the input (as well as the response). Therefore, the input will be
modeled as an acceleration, denoted as y¨. Figure 4-13a shows the structural model and Figure 4-13b shows
the equivalent spring-mass system with base input y¨.
The input, y¨, enters the spring-mass system through the spring located at the base of the structure.
Input through a spring is entered into the equations of motion as a displacement, much like input through
a viscous damper is entered into the equations of motion as a velocity. The displacement is obtained by
integrating y¨ twice with respect to time.
A free-body-diagram analysis of the four-degree-of-freedom system gives the following equations of mo-
tion:
−k
¡
x1 −
R R
y¨dt2
¢
− k (x1 − x2) = m1x¨1
−k (x2 − x1)− k (x2 − x3) = m2x¨2
−k (x3 − x2)− k (x3 − x4) = m3x¨3
−k (x4 − x3) = m4x¨4.
(4-118)
Rearrangement of Eq. (4-118) followed by conversion to matrix form gives
⎡
⎢⎢⎣
m1 0 0 0
0 m2 0 0
0 0 m3 0
0 0 0 m4
⎤
⎥⎥⎦
⎧
⎪⎪⎨
⎪⎪⎩
x¨1
x¨2
x¨3
x¨4
⎫
⎪⎪⎬
⎪⎪⎭
+
⎡
⎢⎢⎣
2k −k 0 0
−k 2k −k 0
0 −k 2k −k
0 0 −k k
⎤
⎥⎥⎦
⎧
⎪⎪⎨
⎪⎪⎩
x1
x2
x3
x4
⎫
⎪⎪⎬
⎪⎪⎭
=
⎧
⎪⎪⎨
⎪⎪⎩
k
R R
y¨dt2
0
0
0
⎫
⎪⎪⎬
⎪⎪⎭
. (4-119)
4-74 Schimel, Ding, Anderson and Grantham
 x4 
 x2 
 x1 
 y 
 .. 
(a) (b)
Figure 4-13: The four-degree-of-freedom system subject to base excitation: (a) structural model,
(b)equivalent mass-spring system.
The next objective is to obtain a transfer function for this system. Since accelerometers will be used to
measure both the input and response, the transfer function used will relate input acceleration to response
acceleration
Gq (s) =
X¨q(s)
Y¨ (s)
. (4-120)
To obtain this transfer function first note that the Laplace transform identity
L
µZ
udt
¶
=
U (s)
s
(4-121)
gives
L
µZ Z
y¨dt2
¶
=
Y¨ (s)
s2
. (4-122)
Also note that, for zero initial conditions,
X¨q (s) = s2Xq (s) , (4-123)
which can be rearranged to give
Xq (s) =
X¨q (s)
s2
. (4-124)
Taking the Laplace transform of Eq. (4-119), with zero initial conditions, and substituting from Eq. 4-114,
yields
⎡
⎢⎢⎣
m1 0 0 0
0 m2 0 0
0 0 m3 0
0 0 0 m4
⎤
⎥⎥⎦
⎧
⎪⎨
⎪⎪⎩
X¨1 (s)
X¨2 (s)
X¨3 (s)
X¨4 (s)
⎫
⎪⎬
⎪⎪⎭
+
1
s2
⎡
⎢⎢⎣
2k −k 0 0
−k 2k −k 0
0 −k 2k −k
0 0 −k k
⎤
⎥⎥⎦
⎧
⎪⎨
⎪⎪⎩
X¨1 (s)
X¨2 (s)
X¨3 (s)
X¨4 (s)
⎫
⎪⎬
⎪⎪⎭
=
1
s2
⎧
⎪⎨
⎪⎪⎩
kY¨ (s)
0
0
0
⎫
⎪⎬
⎪⎪⎭
,
(4-125)
ME 349 Dynamic Systems Lab Manual 4-75
or after rearrangement
s2
⎡
⎢⎢⎣
m1 0 0 0
0 m2 0 0
0 0 m3 0
0 0 0 m4
⎤
⎥⎥⎦
⎧
⎪⎪⎨
⎪⎪⎩
X¨1 (s)
X¨2 (s)
X¨3 (s)
X¨4 (s)
⎫
⎪⎪⎬
⎪⎪⎭
+
⎡
⎢⎢⎣
2k −k 0 0
−k 2k −k 0
0 −k 2k −k
0 0 −k k
⎤
⎥⎥⎦
⎧
⎪⎪⎨
⎪⎪⎩
X¨1 (s)
X¨2 (s)
X¨3 (s)
X¨4 (s)
⎫
⎪⎪⎬
⎪⎪⎭
=
⎧
⎪⎪⎨
⎪⎪⎩
kY¨ (s)
0
0
0
⎫
⎪⎪⎬
⎪⎪⎭
. (4-126)
In compact form Eq. (4-126) can be written as¡
[m] s2 + [k]
¢n
X¨ (s)
o
= k
n
Y¨ (s)
o
. (4-127)
Note that the form of Eq. (4-127) is equivalent to Eq. (4-44). Therefore, components of the frequency
response function are equivalent in form to Eq. (4-60). The frequency response function for this system is
given by
Gp(iΩ) =
X¨p (iΩ)
Y¨ (iΩ)
=
NX
r=1
kψprψ1r
ω 2r − Ω2
. (4-128)
It is important to note that the natural frequencies are the same for the four-degree-of-freedom system in
both the free and forced vibration cases.
4.4 Laboratory #4 Procedures
Free-Vibration Testing
This test will involve conducting approximately five free-vibration tests, each initiated by the impact of a
hammer. The computer system computes the spectrum for each test and then presents the average spectrum
on the screen. Averaging several spectrums tends to smooth the overall spectrum, making it more readable.
1. Prepare the computerized data acquisition system for a free vibration tests, as instructed in the labo-
ratory.
2. Lightly strike the four-degree-of-freedom structure with the supplied impulse hammer. You should
hear a beep from the computer system if data collection was initiated by the strike. If the message
OVERLOAD appears on the computer screen, the structure was struck too hard with the hammer. If
there was an overload, reject the test at the prompt.
3. If there was no overload, accept the test at the prompt.
4. After five tests have been accepted, record the observed natural frequencies in table 4-1 below. You
will be instructed on how to zoom in on each displayed peak in the laboratory. Do this to improve the
accuracy of the readings.
Table 4-1. Natural Frequencies
Mode Frequency, (hz.)
1
2
3
4
Forced-Vibration Testing
Stepped sine testing will be used to obtain the frequency response function and natural frequencies from the
forced tests.
1. Prepare the computerized data acquisition system for forced vibration testing as instructed in the
laboratory.
4-76 Schimel, Ding, Anderson and Grantham
2. Set the sinusoidal output to a frequency well below the first natural frequency. Set the amplitude as
instructed.
3. Initiate the output and allow the system to reach steady state.
4. Initiate the data collection.
5. In Table 4-2 record the frequency response function magnitude at the driving frequency.
6. Increment the frequency. When the driving frequency is close to a natural frequency use the smallest
frequency increments possible. When the driving frequency is far away from natural frequency, use
increments of 0.25hz to 0.5hz.
7. Repeat steps three through six until all four natural frequencies have been scanned.
ME 349 Dynamic Systems Lab Manual 4-77
Table 4-2. Frequency response function.
Freq.(Hz) Magnitude Freq.(Hz) Magnitude Freq.(Hz) Magnitude
Chapter 5
Linear System Identification
5.1 System Identification and the Transfer Function
System identification is a practice concerned with determining the dynamic model of a physical system.
There are many procedures utilized, since systems and their models vary widely. As an example, each of
the procedures used in the previous chapters for determining natural frequency and damping ratio can be
considered to be part of the system identification process. Another part of the linear system identification
process is concerned with determining the order of a system. For a linear system model, the order is a function
of the order and number of governing differential equations it possesses. For example, the single-degree-of-
freedom vibrating system of Chapter 1 is of order two since it has one governing differential equation (or
equation of motion) that contains a second time-derivative. By an equivalent argument, the N -degree-of-
freedom systems discussed in Chapter 4 are of order 2N .
The system identification approach used in this laboratory will utilize the concepts of the transfer function
and the frequency response function. Therefore, before proceeding it will be very useful to re-read Section
2.2 on the Laplace transform and transfer functions, and Section 2.4 on the frequency response function
presented in Chapter 2.
Recall the transfer function for the general linear differential equation, Eq. (2-8), from Chapter 2:
G(s) =
X(s)
F (s)
=
bmsm + bm−1sm−1 + · · ·+ b1s+ b0
ansn + an−1sn−1 + · · ·+ a1s+ a0 . (5-1)
The objective of linear system identification stated in general terms, is to determine the transfer function
constants a0, . . . , an and b0, . . . , bm that will best predict a given output X (s) , resulting from a given input
F (s) .
5.2 The Frequency Response Function
An experimental approach to determining values of the constants in Eq. (5-1) as well as a more detailed
understanding of the transfer function can be gained by studying the frequency response function. As
presented in Chapter 2, the frequency response function for the general system defined by Eq. (2-8) is given
by
G (iΩ) = X (iΩ)
F (iΩ) =
bm (iΩ)m + bm−1 (iΩ)m−1 + · · ·+ b1 (iΩ) + b0
an (iΩ)n + an−1 (iΩ)n−1 + · · ·+ a1 (iΩ) + a0
. (5-2)
Values of the frequency response function at particular forcing frequencies, Ω, can be determined by experi-
ment. As was presented in Chapter 4, for a sinusoidal input
f (t) = A cosΩt+B sinΩt (5-3)
and steady-state response
x (t) = C cosΩt+D sinΩt,
5-78
ME 349 Dynamic Systems Lab Manual 5-79
the value of the experimentally determined frequency response function at the forcing frequency is given by
G¯ (iΩ) = A− iB
C − iD . (5-4)
The magnitude and phase angle are given by
¯¯
G¯ (iΩ)
¯¯
=
r
A2 +B2
C2 +D2
, (5-5)
and
∠Ge (iΩ) = tan−1 BA − tan
−1 D
C
, (5-6)
respectively. As can be seen in Eq. (5-5) the magnitude of G¯ (iΩ) gives the ratio of the amplitude of
the steady-state response function, x (t), over the amplitude of the sinusoidal forcing function, f (t) . The
argument of G¯ (iΩ), Eq. (5-6), is equal to the phase difference between the two sinusoids. Thus, the
frequency response function given in Eq. (5-4) is a compact representation of the magnitude and phase angle
relationship existing between the input and response functions.
A useful approach for experimentally computing the frequency response function G¯ (iΩ) is given in Section
4.2 on spectral analysis in Chapter 4. It will be useful to review topics on the discrete Fourier transform and
frequency response function from that section at this point.
The experimentally observed frequency response function in Eq. (5-4) is a particular instance of the
frequency response function for the general system, as given in Eq. (5-2). For the process of system identifi-
cation, the coefficients a0, . . . , an and b0, . . . , bm in Eq. (5-2) constitute m+n unknowns. To determine these
coefficients it is conceivable that the experimental frequency response function, G¯ (iΩ), could be measured
at m + n different driving frequencies. Then, the system of m + n equations, given by G (iΩ) = G¯ (iΩ) at
the each of the driving frequencies, would be solved for the coefficients. However, as will be seen in a later
section, it is more practical to generate a frequency response plot consisting of a series of values of G¯ (iΩ),
and then deduce the values of the coefficients from the plot.
5.3 The Time Constant
First-Order Systems
It is useful to rearrange the coefficients a0, . . . , an and b0, . . . , bm of the transfer function, Eq. (5-1), into a
more physically meaningful form. To develop an understanding of this process consider the first-order linear
differential equation
a1x˙+ a0x = b0f (t) . (5-7)
Equation (5-7) is equivalent in form to the general differential linear equation given by Eq. (2-8). Therefore,
the transfer function is in the form of Eq. (5-1), and is given by
G (s) =
b0
a1s+ a0
=
b0
a1
1
s+ a0a1
. (5-8)
Defining
x0 =
b0
a1
(5-9)
and
τ =
a1
ao
, (5-10)
gives
G (s) =
x0
s+ 1τ
. (5-11)
5-80 Schimel, Ding, Anderson and Grantham
(a)
(b)
Figure 5-1: First orser system: (a) impulse input; (b) response.
Now consider the case when the input, f (t), is the unit impulse function as shown in Figure 4-1a. The
impulse function can be used to model short duration impact forces, such as those induced by the hammer
tap in the free vibration tests of Chapter 4. The Laplace transform of the impulse is given by
L (f (t)) = 1. (5-12)
For this case the Laplace transform of the input, x (t), is given by
X (s) = G (s)F (s) =
x0
s+ 1τ
. (5-13)
Application of the inverse Laplace transform
L−1
µ
1
s+ α
¶
= e−αt (5-14)
to X (s) gives the time domain solution
x (t) = x0e
−
t
τ , (5-15)
which is shown in Figure 4-1b. The constant τ defined in Eq. (5-10) is known as the time constant for the
first-order linear system defined by Eq. (5-7). Note that the relative magnitude of the time constant defines
how fast the solution given by Eq. (5-15) decays to zero. Relatively large values of τ cause the solution to
decay to zero slowly. Conversely, relatively small values of τ cause the solution to decay to zero quickly.
Note that the solution given in Eq. (5-15) has the form
x (t) = x0eλt, (5-16)
where λ is the eigenvalue of this system. Therefore, the eigenvalue and time constant are related in this
system by the equation
λ = −1
τ
. (5-17)
ME 349 Dynamic Systems Lab Manual 5-81
Note also, that the eigenvalue, or alternatively − 1τ , is a root to the denominator of the transfer function in
Eq. (5-13). For this form of the transfer function, the denominator is known as the characteristic equation
or characteristic polynomial of the system.
Higher-Order Systems
The first-order concepts developed above can be extended to the transfer function of the general linear
system, Eq. (5-1). Equation (5-1) can be rewritten as
G(s) =
X(s)
F (s)
=
dmsm + dm−1sm−1 + · · ·+ d1s+ d0
sn + cn−1sn−1 + · · ·+ c1s+ c0 , (5-18)
where
cr =
ar
an
(5-19)
and
dr =
bq
an
. (5-20)
Factoring the characteristic polynomial in the denominator of Eq. (5-18) gives
G(s) =
dmsm + dm−1sm−1 + · · ·+ d1s+ d0µ
s+
1
τ1
¶µ
s+
1
τ2
¶
· · ·
µ
s+
1
τn
¶ , (5-21)
and, if the time constants take on distinct values, application of the method of partial fraction expansion
allows Eq. (5-21) to be rewritten as
G (s) =
D1µ
s+
1
τ1
¶ + D2µ
s+
1
τ2
¶ + . . .+ Dnµ
s+
1
τn
¶ = nX
r=1
Drµ
s+
1
τr
¶ . (5-22)
If two or more of the time constants are equal then the Eq. (5-22) takes on a different form. However, this
case will not be considered here. For an impulse input, the solution for systems which can be written in the
form of Eq. (5-22) may be written as
x (t) =
nX
r=1
Dre
−
t
τr , (5-23)
or alternatively as
x (t) =
nX
r=1
Dreλrt. (5-24)
Thus, the general solution can be written as the superposition of exponentials, with time constants τ r.
Note that the generalization presented in Eqs. (5-18) through (5-24) does not necessarily assume that
roots of the characteristic polynomial, the eigenvalues, are real valued. It may be that some of the roots
appear as complex conjugate pairs. This case has been already dealt with in previous chapters, where the
eigenvalues are written as
λr = −ζrωr ± ωr
q
1− ζ2r. (5-25)
Thus, in the case of complex eigenvalues the concept of the time constant is replaced by damping ratio and
natural frequency, and part of the general solution is oscillatory. System models that possess both real and
complex eigenvalues are used in practice. However, for the remainder of this chapter only systems with real
valued eigenvalues will be considered.
5-82 Schimel, Ding, Anderson and Grantham
An alternate form of the transfer function, given by Eq. (5-21), will prove useful in system identification.
By factoring the numerator in Eq. (5-21), the transfer function may be rewritten as
G (s) =
K
µ
s+
1
τn+1
¶µ
s+
1
τn+2
¶
· · ·
µ
s+
1
τn+m
¶
µ
s+
1
τ1
¶µ
s+
1
τ2
¶
· · ·
µ
s+
1
τn
¶ , (5-26)
or, after rearrangement,
G (s) = K
Ã
1
s+ 1τ1
!
· · ·
Ã
1
s+ 1τn
!µ
s+
1
τn+1
¶
· · ·
µ
s+
1
τn+m
¶
. (5-27)
Eq. (5-27) may be written more compactly as
G (s) = KG1 (s)G2 (s) . . . Gn+m (s) , (5-28)
where
Gq (s) =
1
s+ 1τq
(5-29)
for q = 1, . . . , n, and
Gq (s) = s+
1
τq
(5-30)
for q = n+ 1, . . . , n+m.
The linear system identification objective may now be restated for the transfer function defined by Eq.
(5-27). The objective is to determine the constants τ1, . . . , τn+m and the coefficient K, which will best
predict the response, X (s) , resulting from a given input, F (s). The frequency response function approach
to determining the transfer function coefficients is also applicable for this formulation. The procedure for
determining the constants is given in the next section.
5.4 The Bode Plot
The Bode plot is a special form of a frequency response plot in which the response data is presented in a
logarithmic format. A complete Bode plot actually contains two plots, one for the magnitude and one for
the phase angle of the frequency response function. An example Bode plot for a first-order system appears in
Figure 5-2. As will be seen, the logarithmic presentation of the data will enable the graphical determination
of system time constants for specific forms of Eq. (5-28).
The Bode plot utilizes a specialized unit known as the decibel, which is denoted as db. Only the magnitude
portion of the frequency response function is converted to decibels on a Bode plot. The frequency response
function magnitude is converted to decibels by the following equation:
g (Ω) = 20 log10 |G (iΩ)| db. (5-31)
If the frequency response function, G (iΩ) , carries its own units, then the units on g (Ω) are often presented
as db-units. For example, if G (iΩ) caries the units of Volts then the units of g (Ω) are presented as dbV.
Note that once a data set is converted to decibels by Eq. (5-31), the results are plotted on a linear scale,
not on a logarithmic scale.
The frequency axis of the Bode plot is logarithmic for both the magnitude and phase portions, but the
frequencies are not converted to units of decibels. The frequencies are simply plotted on a logarithmic axes.
The units for the frequency on a Bode plot must be radians per second to produce correct results from the
system identification procedures described below. The phase angle is plotted on a linear axes, usually in
units of degrees.
A decade is terminology that is often used to describe factor of ten changes along the frequency axis.
For instance, the phrase, ”The frequency response asymptote decreases by 20db per decade.”, means that
the frequency response function magnitude asymptotically approaches a line which decreases by 20db every
time the frequency is increased by a factor of ten.
ME 349 Dynamic Systems Lab Manual 5-83
G
(i 
   
) ,
 d
eg
G
(i 
   
)
, d
b
, rad/sec
Figure 5-2: The Bode plot of a first-order system.
First-Order Systems
Recall the first-order system model given in Eq. (5-7). Since a Bode plot presents the magnitude and phase of
the frequency response function, it will be necessary to obtain these relationships from the transfer function
given by Eq. (5-11). The frequency response function for this system is
G (iΩ) = x0
iΩ+ 1τ
. (5-32)
The frequency response function is converted to a slightly different form for Bode plot analysis:
G (iΩ) = τx0
1 + iΩτ . (5-33)
The magnitude and phase of Eq. (5-33) are given by
|G (iΩ)| = τx0q
1 + (Ωτ)2
, (5-34)
and
∠G (iΩ) = − tan−1 (Ωτ) , (5-35)
respectively. Converting |G (iΩ)| in Eq. (5-34) to decibels gives
g (Ω) = 20 log10 (|G (iΩ)|) = 20 log10 (x0) + 20 log10 (τ)− 10 log10
³
1 + (Ωτ)2
´
, (5-36)
where the following logarithmic identities have been applied:
log
u
v
= log u− log v, (5-37)
log uv = log u+ log v, (5-38)
5-84 Schimel, Ding, Anderson and Grantham
G
(i 
   
)
, d
b
, rad/sec
0 db/decade-20 db/decade
  corner
frequency,   1  
x0 , db
Figure 5-3: Bode magnitude plot with high-frequency and low-frequency asymptotes.
and
log uα = α log u. (5-39)
Some highly important facts can be noted about Eq. (5-36): First, the limit
lim
Ω→0
20 log10 (|G (iΩ)|) = 20 log10 (x0) + 20 log10 (τ) , (5-40)
shows that as the frequency Ω approaches zero on a Bode plot, g (Ω) approaches a value of 20 log10 (τx0) .
Another way of looking at the conclusion drawn from Eq. (5-40) is that as Ω becomes very small, the
function g (Ω) asymptotically approaches a horizontal line emanating from the value 20 log10 (τx0) on the
db axes. This line is known as the low-frequency asymptote. Secondly, the limit
lim
Ω→∞
20 log10 (|G (iΩ)|) = 20 log10 (x0) + 20 log10 (τ)− 20 log10 (Ωτ) (5-41)
shows that as the frequency Ω approaches infinity, g (Ω) approaches another asymptotic line that decreases
by 20db for each decade that Ω increases. This line is known as the high-frequency asymptote. Figure 5-3
shows a Bode plot of the frequency response function magnitude with the two asymptotes drawn in.
The asymptotic lines can be used to determine the value of the constants τ and x0 for the first-order
system. To see this, consider the frequency at which the high and low frequency asymptotes intersect. The
point of intersection is obtained by equating Eqs. (5-40) and (5-41), which gives:
20 log10 (x0) + 20 log10 (τ) = 20 log10 (x0) + 20 log10 (τ)− 20 log10 (Ωτ) (5-42)
or
log10 (Ωτ) = 0. (5-43)
Applying the logarithmic identities in Eqs. (5-37) and (5-38) gives
log10 (Ω) = log10
µ
1
τ
¶
. (5-44)
ME 349 Dynamic Systems Lab Manual 5-85
Taking the inverse-logarithm of both sides of Eq. (5-44) yields
Ω = 1
τ
(5-45)
at the point where the high- and low-frequency asymptotes intersect. The frequency at the point of inter-
section is known as the corner frequency, and it is obviously useful in obtaining the time constant τ . Also,
since τ is now known, as is the value of the low-frequency asymptote at the db axis intercept Eq. (5-40) can
be used to compute the value of the constant x0.
The system identification procedure for a first-order linear system may be outlined as follows:
1. Conduct a sequence of steady-state frequency response tests (i.e. stepped sine tests) and collect sinu-
soidal input and response data, f (n) and x (n) respectively. The range of frequencies must be wide
enough to determine the position of the high and low frequency asymptotes on the Bode Plot.
2. For each frequency, compute the discrete Fourier transform of x (n) and f (n) as defined in Chapter 4,
and use the results to compute the value of the frequency response function at each driving frequency,
Ω.
3. Compute the magnitude and phase of the frequency response function at each driving frequency.
4. Construct a Bode plot from the experimental data.
5. On the magnitude plot, draw the zero and −20db per decade asymptotes.
6. Determine the corner frequency at the intercept of the two asymptotes and compute the time constant,
τ .
7. Determine g (Ω) at the point where the low frequency asymptote intercepts the db axis and compute
the value of x0 from Eq. (5-40).
8. Insert the values of x0 and τ into the transfer function given in Eq. (5-7).
Higher-Order Systems
As is the case for the first-order system, Bode plots of higher-order systems are constructed from the frequency
response function. The transfer function for the general system is given in Eq. (5-26), substituting iΩ for s
in this equation gives the frequency response function:
G (iΩ) =
K
µ
iΩ+ 1
τn+1
¶µ
iΩ+ 1
τn+2
¶
· · ·
µ
iΩ+ 1
τn+m
¶
µ
iΩ+ 1
τ1
¶µ
iΩ+ 1
τ2
¶
· · ·
µ
iΩ+ 1
τn
¶ , (5-46)
or, after rearrangement,
G (iΩ) = K τ1
(1 + iΩτ1)
· · · τn
(1 + iΩτn)
· 1
τn+1
(1 + iΩτn+1) · · · 1τn+1 (1 + iΩτn+1) . (5-47)
The frequency response function of Eq. (5-47) can be written more compactly as
G (iΩ) = KG1 (iΩ)G2 (iΩ) . . . Gn+m (iΩ) , (5-48)
where
Gq (iΩ) =
τ q
(1 + iΩτq)
, (5-49)
for q = 1, . . . , n and
Gq (iΩ) =
1
τq
(1 + iΩτq) , (5-50)
5-86 Schimel, Ding, Anderson and Grantham
for q = n+ 1, . . . , n+m.
The magnitude of G (iΩ) is given by
|G (iΩ)| = K |G1 (iΩ)| |G2 (iΩ)| . . . |Gn+m (iΩ)| , (5-51)
where
|Gq (iΩ)| = τ qq
1 + (Ωτq)2
, (5-52)
for q = 1, . . . , n and
|Gq (iΩ)| = 1τ q
q
1 + (Ωτ q)2, (5-53)
for q = n+1, . . . , n+m, and the values τ q are assumed positive. Similarly, the phase, or argument, of G (iΩ)
is given by
∠G (iΩ) = ∠G1 (iΩ) +∠G2 (iΩ)+, . . . ,+∠Gn+m (iΩ) , (5-54)
where
∠Gq (iΩ) = − tan−1 (Ωτ q) (5-55)
for q = 1, . . . , n and
∠Gq (iΩ) = tan−1 (Ωτq) (5-56)
q = n+ 1, . . . , n+m.
Converting |G (iΩ)| to decibels results in
g (Ω) = 20 log10K +
n+mX
q=1
20 log10 |Gq (iΩ)| (5-57)
or
g (Ω) = 20 log10K +
n+mX
q=1
gq (Ω) , (5-58)
where
gq (Ω) = 20 log10 (τ q)− 10 log10
³
1 + (Ωτq)2
´
(5-59)
for q = 1, . . . , n, and
gq (Ω) = −20 log10 (τq) + 10 log10
³
1 + (Ωτ q)2
´
(5-60)
for q = n+ 1, . . . , n+m.
Observe the similarities between the frequency response magnitude in decibels of the higher-order system,
Eq. (5-57), and that of the first-order system Eq. (5-36). Converting to decibels allows the frequency response
magnitude of a higher-order system to be looked at as the summation of a series of first-order frequency
response magnitudes. Each first-order component of the summation will show behavior similar to that of
the first-order system presented earlier.
Consider a first-order component of a higher-order system in the limit as driving frequency, Ω, approaches
infinity:
lim
Ω→∞
gq (Ω) = 20 log10 (τ q)− 20 log10 (Ωτ q) (5-61)
for q = 1, . . . , n, or
lim
Ω→∞
gq (Ω) = −20 log10 (τq) + 20 log10 (Ωτq) (5-62)
for q = n+ 1, . . . , n+m.
Recall that the gq (Ω) in Eq. (5-61) are related to terms in the denominator of the frequency response
function of Eq. (5-46). Observe that if q = 1, . . . , n as in Eq. (5-61) then the magnitude of gq (Ω) approaches
a −20db per decade asymptote as Ω approaches infinity.
Similarly, recall that the gq (Ω) in Eq. (5-62) are related to terms in the numerator of the frequency
response function of Eq. (5-46). When the frequency response function or the transfer function has non
ME 349 Dynamic Systems Lab Manual 5-87
(a)
(b)
db
db
rad/sec
rad/sec
Figure 5-4: Third-order system behavior: (a) Bode magnitude plot of the first-order components; (b) Bode
magnitude plot of the third-order system.
constant terms in the numerator it is said to posses numerator dynamics. Numerator dynamics occur
when one or more time-derivatives of the input function appear in the input equation. Observe that if
q = n+ 1, . . . , n +m as in Eq. (5-62) then the magnitude of gq (Ω) approaches a positive 20db per decade
asymptote as Ω approaches infinity. Systems with numerator dynamics will not be tested in this laboratory,
so only gq (Ω) in the form of Eq. (5-61) will be given further consideration.
Figure 5-4a shows the components gq (Ω) as in Eq. (5-61) for the case of a third-order system. Each
component possesses a −20db per decade asymptote as discussed above. Each component also has a zero
slope asymptote given by
lim
Ω→0
gq (Ω) = 20 log10 (τq) . (5-63)
As is the case for the first-order system, the intercept between the zero slope asymptote and the −20db per
decade asymptote gives the corner frequency for each component. By the same argument as in Eqs. (5-43)
to (5-45), the corner frequency of each component is equivalent to the inverse of the time constant for that
component.
The solution to the system identification problem is not as simple as determining the τ q from each of
the independent components of Figure 5-4. In general the transfer function data is given as a summation of
all the components, as in Figure 5-4b. In practice it would be difficult to separate the components of Figure
5-4b into those appearing in Figure 5-4a, so an alternate approach in used.
The Bode plot of the third-order system in Figure 5-4b contains four asymptotes. A zero-slope asymptote,
and three others with slopes of −20db per decade, −40db per decade, and −60db per decade respectively.
The slope of the asymptotes come a result of the summation of the first-order components of the system.
By looking at Figure 5-4a it can be seen that the frequency region up to 1τ1 contains no −20db per decade
asymptotes. Hence, the slope of the frequency response function in Figure 5-4b is asymptotic with a zero
slope line in this reigon. In the frequency region between 1τ1 and
1
τ2
in Figure 5-4a, a single −20db per decade
asymptote is present in one of the components. Therefore, the summation of the components leads to an
asymptotic slope of −20db per decade in Figure 5-4b. Similarly between 1τ2 and
1
τ3
there are two −20db per
decade asymptotes present, and the summation of components leads to a −40db per decade asymptote in
5-88 Schimel, Ding, Anderson and Grantham
Figure 5-4b. Finally, beyond 1τ3 there are three −20db per decade asymptotes in Figure 5-4a, which leads to
a −60db per decade asymptote in Figure 5-4b. Obviously, an n-order system, in the absence of numerator
dynamics, will have asymptotes ranging from zero to −20ndb per decade slope.
The asymptotes can be used to determine the time constants and the constant K in Eq. (5-46). Consider
the zero-slope asymptote of the third-order system in Figure 5-4b, and denote this asymptote as L0. The
asymptote is obtained from Eq. (5-57) and is given by
L0 = limΩ→0
g (Ω) = 20 log10K +
3X
q=1
20 log (τ q) . (5-64)
The −20db per decade asymptote is given by
L1 = 20 log10K +
3X
q=1
20 log (τq)− 20 log10 (Ωτ1) . (5-65)
The corner frequency between these two asymptotes occurs at the point L0 = L1, or where
20 log10 (Ωτ1) = 0. (5-66)
Solving for the frequency gives
Ω = 1
τ1
. (5-67)
Therefore, τ1 can be determined from the corner frequency given by the intersection of the zero and the
−20db per decade asymptotes, or at the intersection of L0 and L1.
The −40db per decade asymptote is given by
L2 = 20 log10K +
3X
q=1
20 log (τq)−
2X
q=1
20 log10 (Ωτ q) . (5-68)
The corner frequency at the intersection of the −20db and −40db per decade asymptotes occurs at the point
L1 = L2, or where
20 log10 (Ωτ2) = 0. (5-69)
Solving for frequency gives
Ω = 1
τ2
, (5-70)
and therefore, τ2 can be determined from the corner frequency given by the intersection of the −20db and
the −40db per decade asymptotes.
In general, for an n-order system without numerator dynamics, the asymptotes are given by
Lr = 20 log10K +
nX
q=1
20 log (τ q)−
rX
q=1
20 log10 (Ωτq) , (5-71)
and the rth corner frequency occurs at the point Lr−1 = Lr. The inverse of the rth corner frequency gives
the time constant τr. Once the time constants are determined, the value of the coefficient K in Eq. (5-46)
can be determined from intercept of the zero-slope asymptote on the bode plot and Eq. (5-71) evaluated at
r = 0.
The system identification procedure for higher-order systems is much the same as that given above for
the first-order system. The only difference involves construction of additional asymptote lines to determine
the additional time constants.
The order of a linear system may also be estimated from the slope of the steepest asymptote. Dividing
the steepest slope by the quotient (−20db per decade) will yield an estimate of the system order. However,
keep in mind that the experimental data used to construct the Bode plot of the frequency response function
magnitude has a limited frequency range. It may be that higher-order behavior exists at forcing frequencies
beyond those investigated. The exact order of a linear system may be determined only if the governing
differential equations are known.
ME 349 Dynamic Systems Lab Manual 5-89
rad/sec
rad/sec
db
db
(a)
(b)
Figure 5-5: Behavior of a system with close time constants: (a) Bode magnitude plot of the first-order
components; (b) Bode magnitude plot of the second-order components.
Identification of Systems with Closely Spaced Time Constants
If the time constants of a higher-order system are closely spaced, as for the second-order system in Figure
5-5a, then the location of the asymptotes will obscured. For the second-order system in Figure 5-5b, the
location of the −20db per decade asymptote is not clearly defined by a flat segment on the frequency response
function magnitude plot. In cases such as this, an asymptote is positioned tangent to the frequency response
function while maintaining an appropriate slope (in this case −20db per decade). The time constants are
then determined from the corner frequencies, as usual. This approach works, even in cases where multiple
time constants are closely spaced. Simply draw each asymptote, at the proper slope and tangent to the
response function, and then compute the time constants as usual.
5.5 Apparatus
The apparatus used in this laboratory will consist of DC voltage controlled variable speed pump system,
and the computerized data acquisition and spectral analysis system used in Chapter 4. The pump system
in shown in Figure 5-6. A block diagram of the overall apparatus is shown in Figure 5-7, with the system
to be identified is enclosed within the dashed lines. It consists of a voltage to current converter, a power
supply, a DC motor, a water pump, and a small DC generator which functions as a tachometer. The voltage
to current converter accepts a control input in the form of a voltage, and generates a current proportional
to this input. The power supply utilizes the current as a control signal and produces a drive current for the
motor. The DC motor drive shaft turns the water pump, and the DC generator. The output of the DC
generator and the input to the voltage to current converter are measured by the data acquisition system.
The objective of this laboratory will be to identify the transfer function between the input to the voltage
to current converter, Vin (s) , and the output of the DC generator, Vout (s). Thus, the transfer function to
be identified is
G (s) =
Vout (s)
Vin (s)
= Gp (s)Gm (s)Gg (s) , (5-72)
5-90 Schimel, Ding, Anderson and Grantham
switch
pump
     DC
generator
tubing
   speed
adjustment
Figure 5-6: The pumping system.
Sinusoidal
  Function
 Generator
Voltage-to-
   Current
  Converter
DC Power
  Supply
DC Motor
  & Pump
     DC
Generator
      Data
Acquisition
    System
Spectrum
Analyzer
Computer &
  Plotting
    Device
Figure 5-7: The system identification apparatus.
ME 349 Dynamic Systems Lab Manual 5-91
where Gp (s) is the transfer function for the voltage to current converter and power supply, Gm (s) denotes
the transfer function for the DC motor and pump, and Gg (s) denotes the transfer function of the DC
generator. .
The Bode plot and frequency response techniques describe above will be employed to determine the
constants associated with the overall transfer function.
Governing Differential Equations and Transfer Functions
DC Motor and Pump
The torque output, Γ, of an armature controlled DC motor is a linear function of the current, Ia, which
moves through the armature windings, that is,
Γ = kΓIa, (5-73)
where kΓ is constant. The DC motor drives the pump, which for this laboratory will be modeled as a linear
system. The pump places an load on the motor, primarily due to the inertia of the pump rotor and the work
expended in pushing water through the system. The armature of the DC motor, and the armature of the
DC generator also place an inertia load on the motor. Since it is not the purpose of this lab to determine
the value each of these quantities, the inertia from all these sources will be lumped into a single rotational
moment of inertia, J . Likewise, there are sources of energy loss in the system. For example, flow resistance
from pumping water through the lines and bearing drag. The losses in the system will be assumed to result
from a linear function of the shaft speed, and will be lumped into a single coefficient, c. The resistive torque
generated by both loads
Γ = Jθ¨ + cθ˙, (5-74)
where θ˙ is the motor shaft rotation speed, and θ¨ is the shaft acceleration. Equating Eqs. (5-73) and (5-74)
yields the relationship between the input current and the shaft motion
kΓIa = Jθ¨ + cθ˙. (5-75)
Finally, the linear transfer function between input armature current and motor shaft speed is given by
Gm (s) =
Θ˙ (s)
Ia (s)
=
kΓ
Js+ c
. (5-76)
DC Generator
The DC generator is in fact a DC armature controlled motor, that is being used as a generator. The model
for this type of motor is depicted in Figure 5-8. The generator shaft, which spins at θ˙, generates the back
emf voltage VB proportional to the speed. Thus,
VB = kB θ˙ (5-77)
The voltage available to the measurement system is affected by the generator’s armature inductance, Lg,
the armature resistance, Rg, and a noise filtering capacitance, Cg. If the current in the circuit in Figure 5-8
is denoted by Ig, then the voltage across the inductor is
VL = Lg
dIg
dt
, (5-78)
the voltage across the resistor is
VR = RgIg, (5-79)
and the voltage across the capacitor is
VC = Cg
Z
Igdt. (5-80)
5-92 Schimel, Ding, Anderson and Grantham
Application of Kirkoff’s voltage law yields the following equation
Cg
Z
Igdt+RgIg + Lg
dIg
dt
= VB. (5-81)
Since
Vout = Cg
Z
Igdt, (5-82)
differentiating and rearrangement yields
Ig =
1
Cg
dVout
dt
. (5-83)
By substituting Eqs. (5-77), (5-82), and (5-83) into Eq. (5-81), it can be rewritten as
Vout +
Rg
Cg
dVout
dt
+
Lg
Cg
d2Vout
dt2
= kB θ˙. (5-84)
The transfer function between the shaft speed and Vout can be written as
Gg (s) =
Vout (s)
Θ˙ (s)
=
kB
Lg
Cg
s2 + RgCg s+ 1
. (5-85)
The roots of the denominator (and therefore, the time constants) can be found by applying the quadratic
formula. However, the generator used in this laboratory has very small time constants. This is equivalent to
having very large corner frequencies. The frequency range used in the experiment will be well below these
corner frequencies. Therefore, the response of the generator within the tested frequency range is essentially
a flat line, and the transfer function for the generator can be modeled as a simple constant:
Gg (s) =
Vout (s)
Θ˙ (s)
≈ kB (5-86)
Voltage to Current Converter and Power Supply
The voltage to current converter is an electronic device that converts the input voltage into a current suitable
for as an input to the motor power supply. The converter lacks any significant dynamics within the frequency
range that will be tested. Therefore, it will be modeled as a single constant gain for this laboratory. Denote
Iv as the current out of the converter, and kv as the conversion constant. The relationship between input
voltage and output current is
Iv = kvVin, (5-87)
and the transfer function is
Iv (s)
Vin (s)
= kv. (5-88)
The motor power supply accepts the output of the voltage to current converter as an input and produces
the DC motor armature current as an output. The power supply possesses a soft-start feature which limits
sudden changes in the armature current, and therefore output torque, of the DC motor. The slow start
function is performed by a low pass filter built into the electronics of the power supply. A low pass filter can
be modeled as a first-order linear device. Therefore, the transfer function for the voltage to current converter
and power supply can be written as
Gp (s) =
Iv (s)
Vin (s)
Ia (s)
Iv (s)
=
kvkp
s+ 1RfCf
, (5-89)
where kp is a conversion constant between output and input power supply currents, and Ia is the motor
armature current. The time constant of a first-order low pass filter can be determined by the filters equivalent
resistance, Rf , and capacitance, Cf , as shown in Eq. (5-89).
ME 349 Dynamic Systems Lab Manual 5-93
Cg Ig
Lg
VBVout
Rg
Figure 5-8: The DC generator model.
The System Transfer Function
The overall transfer function for the system can be obtained by combining the component transfer functions
in Eqs. (5-76), (5-86), and (5-88). That is,
G (s) =
Vout (s)
Vin(s)
=
Ia (s)
Vin (s)
Θ˙ (s)
Ia (s)
Vout (s)
Θ˙ (s)
= kB
Ã
kvkp
s+ 1RfCf
!µ
kΓ
Js+ c
¶
, (5-90)
or
G (s) =
kpkΓkB
J
Ã
1
s+ 1RfCf
!µ
1
s+ cJ
¶
. (5-91)
Equation (5-91) has the same form as Eq. (5-27). Therefore, the Bode plot analysis techniques will apply.
5.6 Measurements and Calculations
There are two primary objectives for this laboratory. The first objective will be to obtain a theoretical
estimate for G (s) in Eq. (5-91). The second objective will be to obtain G (s) experimentally using the
stepped-sine testing procedure presented in Chapter 4.
Theoretical Transfer Function
Time Constant of the Motor-Pump-Generator system
The time constant of the motor and pump system, Jc , will be determined by the first-order equivalent of a
free vibration test. Consider the equation of motion for the motor-pump system when no armature current
is applied
Jθ¨ + cθ˙ = 0. (5-92)
5-94 Schimel, Ding, Anderson and Grantham
The solution to this equation is given by
θ˙ (t) = θ˙ (0) e−
c
J
t. (5-93)
At steady state, and at low frequencies, the rotational speed of the generator shaft and the output voltage
are related by
θ˙ =
Vout
kB
. (5-94)
Substitution of Eq. (5-94) into the left-hand side of Eq. (5-93) followed by rearrangement gives
Vout (t) = θ˙ (0) kBe−
c
J
t. (5-95)
Equation (5-95) shows that if the motor shaft is rotating at initial speed, θ˙ (0) , with no armature current
applied, then the output of the generator will show an exponential decay. The time constant associated with
the decay, Jc , is same as the time constant that appears in the transfer function for the overall system, Eq.
(5-91).
Equation (5-95) can be used to develop experimental method for estimating the time constant Jc . The
initial condition can be achieved by applying power to the motor until it reaches a desired speed, and then
switching the motor power off. The motor will rotate freely with the speed, and the voltage out of the
generator, decaying to zero. Observation of Figure 5-1 shows that after a period of four time constants, the
response of a theoretical system will decay to less than 2% of its initial value. Therefore, measuring the
time, td, that it takes for the output of the generator to decay down to 2% of the initial value will yield an
estimate of time constant:
J
c
=
td
4
. (5-96)
Determination of the Soft-Start Time Constant
The power supply for the DC motor contains a low pass filter which performs the soft-start function described
above. The filter contains a known capacitance, Cf . However, the resistance Rf is a combination of several
resistors and other electronic devices built into the power supply circuitry. The analysis of this circuit will
not be performed here. Instead, the time constant RfCf will be determined experimentally by applying a
step voltage at the input.
Consider the response of the entire system, as defined in Eq. (5-91), when it is subject to a step input
voltage, V◦, at time t = 0. The Laplace transform of the step input is given by
Vin (s) =
V◦
s
, (5-97)
and the output is given by
Vout (s) = G (s)Vin (s) =
kpkΓkB
J
Ã
1
s+ 1RfCf
!µ
1
s+ cJ
¶
V◦
s
. (5-98)
Application of the method of partial fraction expansion results in an equation of the form
Vout (s) =
K1
s
− K2
s+ 1RfCf
− K3
s+ cJ
, (5-99)
where the Ki are constants. Computing the inverse Laplace transform of Eq. (5-99) gives
Vout (t) = K1 −K2e
− t
RfCf −K3e−
c
J
t. (5-100)
For the system being tested in this laboratory, the time constant cJ is several times larger than the time
constant tRfCf . Therefore, the exponential containing the time constant
c
J in Eq. (5-100) decays to zero
ME 349 Dynamic Systems Lab Manual 5-95
much more quickly than the exponential containing 1RfCf . This means that for t > 4
c
J Eq. (5-100) can be
approximated by
Vout (t) ≈ K1 −K2e
− t
RfCf . (5-101)
By using an argument similar to the one that lead to Eq. (5-96), it can be seen that
RfCf =
td
4
. (5-102)
The time td in Eq. (5-102) is that required for the system modeled by Eq. (5-101) to reach approximately
98% of its final value.
Determination of the Transfer Function Leading Constant
Application of the method of partial fraction expansion to Eq. (5-98) will show that
K1 =
kvkpkΓkBRfCfV◦
c
, (5-103)
where K1 is as defined in Eq. (5-100). Equation (5-103) gives the value that the system modeled by Eq.
(5-100) reaches at steady state conditions. Dividing K1 in Eq. (5-103) by the two system time constants
and the step input voltage gives
K1¡
V◦RfCf Jc
¢ = kvkpkΓkB
J
. (5-104)
This is the leading coefficient to the transfer function G (s) in Eq. (5-91). Therefore, by applying a step
input to the overall system and measuring the steady-state response, it is possible to calculate the leading
coefficient to G (s) from known quantities.
5.7 Laboratory #5 Procedures
Time Constant of the DC Motor, Pump and Generator System
1. Connect the output of the DC Generator to the oscilloscope, as instructed in the laboratory.
2. Press the acquire button located on the upper right of the oscilloscopes front panel.
3. Set the toggle switch on the back side of the DC motor housing to the Int. position.
4. Set the speed knob on the DC motor housing front panel to a value of 10.
5. Set the power switch on the DC motor housing front panel to the forward position.
6. Wait approximately 10-15 seconds for the motor speed to reach steady-state.
7. Once the motor speed has reached steady-state, set the power switch to the off position. Note: Be
careful not to advance the power switch beyond the off position and into the reverse position.
8. Once the decay portion of the generator output trace appears on the oscilloscope screen, press the save
button at the upper right of the oscilloscope’s front panel.
9. Determine the time it takes for the trace to reach approximately 2% of its initial value. Note that
there is some noise in the DC generator output. However, using the mean value of the signal will yield
reasonably accurate results.
5-96 Schimel, Ding, Anderson and Grantham
Power Supply Soft-Start Time Constant
1. Leave the generator output connected to the oscilloscope.
2. Change the oscilloscopes time base as instructed.
3. Press the oscilloscopes acquire button.
4. Set the toggle switch on the back side of the DC motor housing to the Ext. position.
5. Turn the pump motor switch to the forward position.
6. Initiate a voltage step output from the computer system as instructed.
7. Wait until the output trace from the generator reaches an approximate steady-state condition on the
oscilloscope screen.
8. Press the save button on the oscilloscope.
9. Cancel the step output from the computer system.
10. Determine the time it takes for the trace to reach 98% of its approximate steady-state value.
Transfer Function Leading Coefficient
1. Initiate a voltage step output from the computer system for each of the voltage levels in table 5-1
below.
2. Wait for the system to achieve steady-state conditions in each case.
3. Record the steady-state voltage output of DC Generator from the oscilloscope in Table 5-1.
Table 5-1: Steady-State Step Response Voltages
Input Voltage (V) Response Voltage (V)
2.5
5.0
7.5
10.0
Frequency Response Function
Use the stepped-sine testing procedure from Chapter 4 to determine the frequency transfer function for this
system. For this experiment, the motor shaft speed, θ˙, will be varied sinusoidally. The mean level of the
input sinusoid will be offset from zero volts, so that the minimum motor speed will be approximately zero.
The discrete Fourier transform is able to separate the offset from the remainder of the input at the driving
frequency so no effects should be observed in the transfer function.
Record results in Table 5-2 for the test frequency values given during the laboratory.
ME 349 Dynamic Systems Lab Manual 5-97
Table 5-2. Frequency response function.
Freq.(hz) Magnitude Freq.(hz) Magnitude Freq.(hz) Magnitude
Chapter 6
Control of Dynamic Systems
In previous chapters the properties of dynamic systems were determined by applying some type of input, such
as an impulse from a hammer tap or harmonic excitation, and then measuring and analyzing the response.
However, there are times when it is the objective to obtain a specific response from a dynamic system.
In such cases a specific input or control signal may be applied with the intent of producing this response.
However, for a given dynamic system it may turn out that the process of directly applying a control signal
leads to an inappropriate response. This may be due to the time constants or stability characteristics of the
original system. In these cases the behavior of the system may be altered by incorporating a device known
as a controller. The controller may in itself be looked at as a dynamic system that, when combined with the
original dynamic system, gives the desired response to an input function. This chapter will deal primarily
with an introduction to linear controllers and the control of linear or linearized dynamic systems.
As a controls example, consider an automobile. The driver can control the speed of the automobile by
positioning the accelerator pedal; or, if present, a cruise control can be used to automatically control the
speed. A very simple cruise-control system can be devised by determining a set of throttle positions that
give specific speeds and then incorporating some mechanism into the automobile that sets the appropriate
throttle position when a certain speed is desired. This control scheme, known as open-loop control, lacks
robustness however. If a hill or winds are encountered then the speed of the automobile may change and a
new throttle setting would be needed to return the vehicle to the desired speed. A more robust approach
to controlling the vehicle is to make active measurements of the speed, and then automatically adjust the
throttle setting based on what speed is. If the measured speed is too low then the throttle may be opened by
some amount, or, if the speed is too high then it may be closed by some amount. The process of continuously
measuring the speed and automatically adjusting the throttle based on the speed is an example of closed-loop
control. Figure 6-1a shows a basic closed-loop control system. For the automobile, the desired speed is the
input, d (t), and the measured speed, x (t), is the feedback signal. Subtracting the measured speed from
the desired speed generates the error signal, e (t). The controller processes current and past error signals
to compute a control signal, c (t), to the which is sent to actuator. The actuator produces the output f (t)
which drives the plant. For the automobile, the control signal is the throttle position, the actuator is the
engine and drivetrain, and the plant is a model which includes the vehicle mass, wind drag, and gravitational
forces. Specific linear control schemes will be discussed in a later section.
6.1 The Transfer Function and Block Diagram Analysis
The feedback control system shown in Figure 6-1a can be looked at as an ordinary dynamic system with
general input d (t) and response x (t). Assume that the system is linear, and recall that the transfer function
is defined by a ratio of the Laplace transform the input, D (s) , over the Laplace the output, X (s),
G (s) =
X (s)
D (s)
. (6-1)
6-98
ME 349 Dynamic Systems Lab Manual 6-99
Controller Actuator Plant
e(t)d(t) c(t) f(t) x(t)
 - 
 Gc(s) Ga(s)  Gp(s)
E(s)D(s) C(s) F(s) X(s)
 - 
 Gc(s)
 H(s)
Ga(s)  Gp(s)
E(s)D(s) C(s) F(s) X(s)
 - 
(a)
(b)
(c)
Figure 6-1: General block diagrams for control systems.
The general form of the transfer function for a linear system is
G (s) =
bmsm + bm−1sm−1 + · · ·+ b1s+ b0
ansn + an−1sn−1 + · · ·+ a1s+ a0 , n > m. (6-2)
The transfer function can be looked at as a compact mathematical representation of a linear system. If the
system has been identified as defined in Chapter 5, then the constants a0, . . ., an as well as b0, . . ., bm are
known values. The denominator in Eq. (6-2) gives the characteristic equation for the system. Solving
ansn + an−1sn−1 + · · ·+ a1s+ a0 = 0 (6-3)
for roots of s gives the eigenvalues, λ’s, of the linear system. As will be seen later, the eigenvalues are very
useful in analyzing the behavior of control systems.
The plant of a general linear system is be modeled by a transfer function that relates the input driving
signal from the actuator to the response:
Gp (s) =
X (s)
F (s)
(6-4)
Similarly, an actuator that is linear may be modeled with a transfer function that relates the input control
signal to the output that drives the plant
Ga (s) =
F (s)
C (s)
. (6-5)
If the controller is linear, then it may also be modeled by a transfer function relating the input error signal
to the output control signal
Gc (s) =
C (s)
E (s)
. (6-6)
Some example linear controllers will be discussed in a later section. For now, just assume that the controller
has a linear transfer function in the form of Eq. (6-2). The presence of the controller as well as the feedback
6-100 Schimel, Ding, Anderson and Grantham
loop will alter the dynamic behavior of the system. In some cases the order of the system is changed and in
some cases the eigenvalues of the system may change or both.
If the time-domain feedback control system shown in Figure 6-1a is converted into the Laplace domain
then the system can be represented as in Figure 6-2b. The Laplace-domain control system shown in Figure
6-2b is commonly called a block diagram. A more general block diagram can result if the transducer used
to make the feedback measurement is also a linear dynamic system. In such cases the transfer function of
the transducer, which will be denoted by H (s) , is placed directly in the feedback path as is shown in figure
6-2c.
The system given by the more general block diagram of Figure 6-2c can be changed to a system with a
single input and a single output with no feedback loop through a process known as block diagram reduction.
Block diagram reduction can be achieved through the application of simple algebraic rules to obtain a transfer
function that relates the input desired response, D (s), to the actual response, X (s). Consider the error
signal
E (s) = D (s)−H (s)X (s) , (6-7)
and the relationship between the output and the error signal
X (s)
E (s)
=
C (s)
E (s)
F (s)
C (s)
X (s)
F (s)
= Gc (s)Ga (s)Gp (s) . (6-8)
Solving Eq. (6-8) for E (s) and substituting into Eq. (6-7) gives
X (s)
Gc (s)Ga (s)Gp (s)
= D (s)−H (s)X (s) . (6-9)
Solving Eq. (6-9) for the transfer function between the output X (s) and the input D (s) gives the transfer
function of the control system after block diagram reduction:
G (s) =
X (s)
D (s)
=
1
1
Gc (s)Ga (s)Gp (s)
+H (s)
=
Gc (s)Ga (s)Gp (s)
1 +Gc (s)Ga (s)Gp (s)H (s)
. (6-10)
Keep in mind that if the controller, actuator, plant, and transducer of a system are linear, then it is possible
to write the linear transfer function defined by Eq. (6-10). Since Eq. (6-10) is linear it will reduce to the
form of Eq. (6-2), and the eigenvalues of the system can be computed as described from Eq. (6-3). An
analysis of the eigenvalues will lead to an understanding to the behavior of the overall system.
6.2 Stability
Stability is an important concept in control system design. For a linear control system to be successful, the
output of the plant must asymptotically approach some value in the neighborhood of desired response with
time. This type of response is known as asymptotically stable response. Asymptotic stability in a linear
control system helps ensure that if the system is disturbed some arbitrary distance away from the desired
response then it will return to a level at or near desired response with time. If the desired response is a
constant, then it will be seen that a given asymptotically stable linear system consisting of a controller,
actuator, plant, transducer, and feedback loop will either asymptotically approach the desired response or
some static level near the desired response depending on the type of controller selected.
For a linear system, the eigenvalues from the solution to Eq. (6-3) may be used to determine the stability.
That is, the roots of the transfer functions’ characteristic equation are used to determine the stability of a
system. A linear system of order n will have n eigenvalues associated with it, and the eigenvalues may be
complex valued, i.e.
λr = αr + iβr, r = 1, . . . , n. (6-11)
For asymptotic stability in a linear system, it is required that the real part of each eigenvalue satisfy the
inequality
Re (λr) = αr < 0. (6-12)
ME 349 Dynamic Systems Lab Manual 6-101
If the system is asymptotically stable, then the approach of the plants’ response to a neighborhood of the
desired response may be purely exponential, or it may be oscillatory. However, if the approach is oscillatory
the amplitude of oscillation about the desired response will decay exponentially. If the real part one or more
of the eigenvalues is positive then the linear system will be unstable. If the system is unstable then the
response of the linear system may diverge away from the desired response, or, oscillate about the desired
response with increasing amplitude. If the real part of one or more of the eigenvalues is equal to zero then
the system can either be neutrally stable (which is sometimes denoted as stable), or, in some particular cases
involving repeated eigenvalues with zero real parts, unstable. Only the asymptotically stable case will be
given further consideration.
6.3 The Initial and Final Value Theorems
The steady-state response of a system is defined as the output generated by a system in the limit as time
approaches infinity. The final value theorem can be used to quickly access the steady-state response of
a system from its Laplace transform, without solving the differential equation. Suppose that the Laplace
transform of the input, D (s), and the transfer function, G (s) , are known. Given that the Laplace transform
of x (t) ,
X (s) = G (s)D (s) , (6-13)
and the Laplace transform of dx/dt exist, and given that the characteristic equation of sX (s) has asymp-
totically stable eigenvalues (Re (λr) < 0), the steady-state response for x (t) is given by the final value
theorem:
lim
t→∞
x (t) = lim
s→0
[sX (s)] = lim
s→0
[sG (s)D (s)] . (6-14)
The final value theorem can also be applied to signals other than the response. For example, the Laplace
transform of the error signal in a control system is defined as
E (s) = D (s)−X (s) = D (s) (1−G (s)) , (6-15)
and hence, the steady-state error between the output and input is given by
lim
t→∞
e (t) = lim
s→0
[sE (s)] = lim
s→0
[sD (s) (1−G (s))] . (6-16)
The initial value of the response, x (0), can be determined from the Laplace transform, X (s), by applying
the initial value theorem,
x (0) = lim
s→∞
[sX (s)] . (6-17)
The theorem is valid only if the Laplace transforms of x (t) and dx/dt exist, and if the limit in Eq. (6-17)
exists. As with the final value theorem, the initial value theorem can be applied to other signals in the
system such as the error signal, e (t).
6.4 Linear Control Laws
A control law is defined by a specific function that produces a control signal as an output. The control laws
that will be considered in this chapter operate linearly on the error signal, e (t), to produce the control signal,
c (t). The linear operations that make up a control law can include multiplication by a constant, integration,
and/or differentiation of the error signal, e (t) . These three operations are referred to as proportional control,
derivative control, and integral control, respectively. A linear controller can contain one or a summation of
two or more of the operations.
The control laws will be applied to an example theoretical first-order plant with an input f (t) from the
actuator:
x˙− x = f (t) (6-18)
The objective in this example is to make the response, x (t), follow some desired response at the input.
Consider a special case of a system defined by the plant only, without a controller, actuator, or feedback
6-102 Schimel, Ding, Anderson and Grantham
loop. The input that the plant output is to follow is the signal f (t) in this special case. The system defined
by Eq. (6-18) has the transfer function
Gp (s) =
X (s)
F (s)
=
1
s− 1 . (6-19)
For this example, the desired-response input will be a unit-step function (f (t) = 1 for t ≥ 0) which has the
transfer function
F (s) =
1
s
. (6-20)
Therefore, the Laplace transform of the output is
X (s) =
1
s (s− 1) . (6-21)
Note that the characteristic equation, as defined in Eq. (6-3), of the quantity sX (s) yields the positive real
eigenvalue λ = 1. Therefore, as discussed above for Eq. (6-14), the final value theorem cannot be applied.
However, the inverse Laplace transform gives the solution
x (t) = et − 1. (6-22)
It can be seen in Eq. (6-22) that x (t) goes toward infinity as time increases. Since the desired response is
x (t) = f (t) = 1, (6-23)
the plant alone will not follow the input. Therefore, a control system must be employed to achieve the
desired response. A block diagram of the control system that will be used in the proceeding examples is
shown in Figure 6-1c.
Proportional Control
A Proportional-control signal is generated by multiplying the error signal by a constant, Kp, as in the
equation
c (t) = Kpe (t) . (6-24)
The transfer function of the proportional-control equation is
Gc (s) =
C (s)
E (s)
= Kp. (6-25)
If the proportional controller is used in the system depicted in Figure 6-1c, then application of Eq. (6-10)
for block diagram reduction, with H (s) and Ga (s) equal to one gives the overall transfer function
G (s) =
Gc (s)Gp (s)
1 +Gc (s)Gp (s)
=
Kp
1
s− 1
1 +Kp
1
s− 1
=
Kp
s+Kp − 1
. (6-26)
An eigenvalue analysis of Eq. (6-26) results in
λ = 1−Kp. (6-27)
Thus, the system is asymptotically stable for values of Kp greater then one, and, since the eigenvalue is
real-valued, the solution will be nonoscillatory.
Consider the case where the desired response, d (t) , is a unit-step function. The unit-step has the transfer
function
D (s) =
1
s
. (6-28)
ME 349 Dynamic Systems Lab Manual 6-103
The Laplace transform of the response is given by
X (s) = G (s)D (s) =
Kp
s (s+Kp − 1)
, (6-29)
and the Laplace transform of the error signal is
E (s) = D (s)−G (s)D (s) = 1
s
µ
1− Kp
s+Kp − 1
¶
. (6-30)
Application of the final-value theorem to Eqs. (6-29) and (6-30) gives the steady-state response
lim
t→∞
x (t) = lim
s→0
s
µ
Kp
s (s+Kp − 1)
¶
=
Kp
Kp − 1
(6-31)
and the steady-state error signal
lim
t→∞
e (t) = lim
s→0
s
s
µ
1− Kp
s+Kp − 1
¶
= 1− Kp
Kp − 1
= − 1
Kp − 1
. (6-32)
The steady-state error signal is useful measure for determining how close the actual response will approach
the desired response level. Therefore, it will be used to analyze the steady-state behavior of many of the
controller types below.
Observe that the proportional-control error signal is nonzero as t→∞. This is referred to as a nonzero
steady-state error. A nonzero steady-state error is typical for systems that are subject to purely proportional
control. Note that as the positive valued constant Kp becomes large the steady-state error approaches zero.
For the example system, Kp can also be assigned an arbitrarily large value without affecting the stability of
the response. The value of Kp will affect the speed of the response however. Large values of Kp give a small
time constant and therefore a fast response for the first-order example. In practice however, increasingly
large values of Kp may not result in faster plant response times due to power limitations of the actuator.
Integral Control
An integral control signal is generated by integrating the error signal and multiplying it by a constant, KI .
For a zero initial error, this operation gives control signal
c (t) = KI
Z t
0
e (t) dt, (6-33)
which has the transfer function
Gc (s) =
C (s)
E (s)
=
KI
s
. (6-34)
Application of Eq. (6-10) to the example first-order system for the case of integral control with Ga (s) and
H (s) equal to one gives
G (s) =
Gc (s)Gp (s)
1 +Gc (s)Gp (s)
=
KI
s
1
s− 1
1 +
KI
s
1
s− 1
=
KI
s2 − s+KI
. (6-35)
Application of the quadratic formula to the characteristic equation of Eq. (6-35) gives
λ1,2 =
1±√1− 4KI
2
. (6-36)
Observe that regardless of the value selected for KI , at least one of the roots will have a positive real part.
This means that the transfer function of the example system with integral control is unstable and pure
integral control in not suitable for this system.
6-104 Schimel, Ding, Anderson and Grantham
Proportional-Plus-Integral Control
Although pure integral control is not suitable for the first-order example system, a combination of propor-
tional and integral control can be useful. Proportional-plus-integral control is often referred to as PI -control.
For zero initial error, the signal generated by PI-control is given by
c (t) = Kpe (t) +KI
Z t
0
e (t) dt, (6-37)
and the transfer function is given by
Gc (s) =
C (s)
E (s)
=
Kps+KI
s
. (6-38)
For the example system with PI-control, application of Eq. (6-10) gives
G (s) =
Gc (s)Gp (s)
1 +Gc (s)Gp (s)
=
Kps+KI
s
1
s− 1
1 +
Kps+KI
s
1
s− 1
=
Kps+KI
s2 + (Kp − 1) s+KI
. (6-39)
Application of the quadratic formula to Eq. (6-39) gives the eigenvalues
λ1,2 =
1−Kp ±
q
(Kp − 1)2 − 4KI
2
. (6-40)
Both the eigenvalues will have a negative real part if Kp > 1 and KI ≥ 0. Therefore, the overall system will
be stable if the appropriate constants are selected.
Note that if 4KI > (Kp − 1)2 in Eq. (6-40), the eigenvalues will be complex valued and the response
will be oscillatory. Improper selection of the control constants in this case can lead to large amplitudes of
oscillation and require long times for the system to stabilize on the desired response level. This can be seen
more clearly if the standard form of the characteristic equation for a second-order system is considered:
s2 + 2ζωns+ ω
2
n = 0, (6-41)
were the constant ζ defines the damping ratio and ωn defines the natural frequency. Comparing this equation
to the characteristic equation for the PI-control system,
s2 + (Kp − 1) s+KI = 0, (6-42)
gives
2ζωn = Kp − 1 (6-43)
and
ω2n = KI . (6-44)
Solving for the natural frequency and damping ratio gives
ωn =
p
KI (6-45)
and
ζ =
Kp − 1√
KI
. (6-46)
Observation of Eq. (6-46) reveals that for values of Kp near one and large values of KI the damping ratio
is small, which can lead to the response problems discussed above.
Application of the final value theorem to the PI-control system defined by Eq. (6-39) with step input
gives the steady-state error signal
lim
t→∞
e (t) = lim
s→0
s
s
µ
1− Kps+KI
s2 + (Kp − 1) s+KI
¶
= 1− KI
KI
= 0. (6-47)
Equation (6-47) reveals advantage that PI control has over proportional control; for some systems it can
result in a zero steady-state error.
ME 349 Dynamic Systems Lab Manual 6-105
Derivative Control
The derivative control signal is generated by differentiating the error signal
c (t) = KD
de
dt
(6-48)
and for zero initial error the transfer function is
Gc (s) =
C (s)
E (s)
= KD s. (6-49)
Application of Eq. (6-10) to the example first-order system for the case of PI control with Ga (s) and H (s)
equal to one gives
G (s) =
Gc (s)Gp (s)
1 +Gc (s)Gp (s)
=
KD s
1
s− 1
1 +KD s
1
s− 1
=
KD s
(KD + 1) s− 1
(6-50)
The transfer function is stable only if KD < -1. Application of the final value theorem to the first-order
system with derivative control and step-input gives the steady-state error signal
lim
t→∞
e (t) = lim
s→0
s
s
µ
1− KD s
(KD + 1) s− 1
¶
= 1. (6-51)
The system response does not approach the desired response since the steady-state error is equal to a
constant. Therefore, derivative control cannot be used for this system. Purely derivative control is not used
in practice because the control signal is dependent only on the rate of change of the error signal.
Proportional-Plus-Derivative Control
Proportional-plus-derivative control, which is referred to as PD-control, can produce a useful control action.
The PD-control signal is given by
c (t) = Kpe (t) +KD
de
dt
(6-52)
and the transfer function is
Gc (s) =
C (s)
E (s)
= Kp +KD s. (6-53)
Application of Eq. (6-10) to the example first-order system for the case of PD-control with Ga (s) and H (s)
equal to one gives
G (s) =
Gc (s)Gp (s)
1 +Gc (s)Gp (s)
=
(Kp +KD s)
1
s− 1
1 + (Kp +KD s)
1
s− 1
=
KD s+Kp
(KD + 1) s+Kp − 1
, (6-54)
which is stable for Kp > 1 and KD > −1. Application of the final value theorem gives the steady-state error
signal
lim
t→∞
e (t) = lim
s→0
s
s
µ
1− KD s+Kp
(KD + 1) s+Kp − 1
¶
=
1
Kp − 1
. (6-55)
For large values of Kp the system response approaches the desired response level.
Application of the initial value theorem for the step-input case gives the initial response
x (0) = lim
s→∞
s
s
µ
KD s+Kp
(KD + 1) s+Kp − 1
¶
=
KD
KD + 1
. (6-56)
The initial response for the previous controlled systems was zero for each case other than derivative control.
However, for PD-control of a first-order system, there is an initial step in the response. This type of control
action may or may not be desirable, depending on the type physical system the first-order model represents.
6-106 Schimel, Ding, Anderson and Grantham
Proportional-plus-derivative control is often used in higher-order systems to increase the level of damping
in the control action. Consider the following second-order plant:
mx¨+ cx˙+ kx = f (t) (6-57)
This system has the transfer function
Gp (s) =
X (s)
F (s)
=
1
ms2 + cs+ k
(6-58)
Application of Eq. (6-10) to this system for the case of PD-control with Ga (s) and H (s) equal to one gives
G (s) =
Gc (s)Gp (s)
1 +Gc (s)Gp (s)
=
(Kp +KD s)
1
ms2 + cs+ k
1 + (Kp +KD s)
1
ms2 + cs+ k
=
Kp +KD s
ms2 + (c+KD) s+ (k +Kp)
. (6-59)
By applying the same logic as in Eqs. (6-41) through (6-46) and assuming the undamped case, the natural
frequency for this system is given by
ωn =
r
k +Kp
m
(6-60)
and the damping ratio is given by
ζ =
c+Kd
2
p
(k +Kp)m
. (6-61)
Thus, increasing Kp will increase the natural frequency and reduce the damping ratio, and increasing KD
will increase the damping ratio of the system. Note that for the case of a unit-step input, the final value
theorem gives
lim
t→∞
e (t) = lim
s→0
s
s
µ
1− Kp +KD s
ms2 + (c+KD) s+ (k +Kp)
¶
= 1− Kp
k +Kp
. (6-62)
Therefore, the value of Kp should be kept relatively large to maintain a small steady-state error for this case.
It is very important to realize that selecting large values for KD in most systems can lead to undesirable
effects. This is due to the fact that most control systems rely on electric signals. Electric signals are easily
corrupted by noise from the environment, which can often be characterized by sudden and random variations
in the signal level. Sudden changes in the signal will result in locally steep slopes. Thus, noise in a function,
such as in the error signal, can result in very large slope magnitudes. If the value of KD is large and the
signal is noisy, large impulses can be transferred to the plant which can easily disturb the desired response.
Therefore, the amount of noise in a system can limit the maximum useful value of the derivative control
constant.
Another problem, associated with amplification of signal noise, exists when derivative control is applied
to a continuous time system. The problem occurs even when derivative control is applied in combination
with proportional control or with PI-control. Consider the frequency response function magnitude obtained
from the derivative control transfer function in Eq. (6-49):
|Gc (iω)| = KD ω. (6-63)
Note that the magnitude of the transfer function, which is equivalent to the output of the controller, increases
with increasing frequency. Therefore, if high frequency noise is present, the output of the controller may be
very large and unpredictable. This may lead to undesirable and unpredictable motion of the plant.
One remedy to the derivative control noise problem is to use a control law known as pseudo-derivative
feedback control. The transfer function for this type of controller is
Gc (s) = KD
s
1 + τds
, (6-64)
and the frequency response function magnitude is
|Gc (iω)| = KD ωp
1 + τ2d ω
2
, (6-65)
ME 349 Dynamic Systems Lab Manual 6-107
where τd is a constant. Application of the system identification techniques of Chapter 5 will show that
the pseudo-derivative feedback controller has a corner frequency of 1τd . At low frequencies, the frequency
response function magnitude of the pseudo-derivative feedback controller is asymptotic with 6-63, and at
high frequencies it approaches the constant
|Gc (iω)| = KDτ s . (6-66)
Therefore, a pseudo-derivative feedback controller approximates the derivative controller at low frequencies
but filters out high-frequency inputs. The method applies to systems where the spectrum of the primary
control signal is confined to a relatively low frequency range when compared to the inherent signal noise.
6.5 PID-Control
Combining the control actions of PD-control and PI-control results in PID-control. The PID-control signal
is given by
c (t) = Kpe (t) +KI
Z t
0
e (t) dt+KD
de
dt
(6-67)
and the transfer function by
Gc (s) =
C (s)
E (s)
= Kp +
KI
s
+KD s (6-68)
Application of Eq. (6-10) to the example first-order system for the case of PID-control with Ga (s) and
H (s) equal to one gives
G (s) =
Gc (s)Gp (s)
1 +Gc (s)Gp (s)
=
µ
Kp +
KI
s
+KD s
¶
1
s− 1
1 +
µ
Kp +
KI
s
+KD s
¶
1
s− 1
=
KD s2 +Kps+KI
(1 +KD) s2 + (Kp − 1) s+KI
, (6-69)
For the case of a unit-step input, the final value theorem gives
lim
t→∞
e (t) = lim
s→0
s
s
µ
1− KD s
2 +Kps+KI
(1 +KD) s2 + (Kp − 1) s+KI
¶
= 1− KI
KI
= 0. (6-70)
A PID-controller is a generalization of the three basic types of linear controllers discussed previously. For
a general linear system PID-control can be used to change the response time constants and damping levels
plus yield a nonzero steady-state error.
6.6 Transient Response Performance Criteria
The performance of a linear dynamic control system can be characterized by a set of criteria that are
measured directly from the response to a step function input. The criteria are maximum overshoot, peak
time, rise time, delay time, settling time and cycle time. The measurement of each criteria is shown in figure
6-2. Although the measurements are shown for the case of a second order system, they may be applied
to a system of any order. The appropriate performance criteria may be related directly to the natural
frequency and damping ratio for an underdamped system, and the time constants of an overdamped system
for controller design purposes.
Maximum overshoot, xm, defines the maximum amount by which the system response overshoots the
desired response. It is a property of second- and higher-order systems that are underdamped. For first-order
and overdamped systems, the maximum overshoot is zero. Peak time, tp, is the elapsed time it takes to
reach the point of maximum overshoot from the initiation of the step input.
The rise time, tr, defines the amount of time required for the response to pass from some initial percentage
of the desired response level to some final percentage of the desired response level. There are no defined
6-108 Schimel, Ding, Anderson and Grantham
maximum
overshootpeak
time
100% rise
    time
delay
 time
settling
  time
 desired
response
  level
x(t)
  t  
Figure 6-2: Control system performance criteria.
initial and final values for the percentages. However, 10% to 90% and 0% to 100% rise times are commonly
used. The 0% to 100% rise time will be used in this laboratory for underdamped response cases. The delay
time, td, can be defined as a special case of the rise time. The delay time is the time it takes for the response
to pass from 0% to 50% of the desired response level.
The settling time, ts, is the time elapsed required for the response to arrive and stay within some
percentage of the final value starting from the initiation of the step input. There is no defined percentage
for the settling time; 2% and 5% are common. For a first order system, the 2% settling time is equivalent to
approximately four system time constants.
The cycle time, tc, is used only for underdamped response cases. It is simply the time required for the
underdamped response to complete one full cycle, and is equivalent to the response period.
6.7 Discrete-Time Control Systems
Microprocessor-based control has become the predominant control method in many application areas. Mi-
croprocessors use sampled data, as described in Chapter 4, to compute the control signal. If the desired
response, d (t) , and the actual response, x (t), are sampled, the respective discrete-time counterparts are
given by d (k) and x (k), k = 0, 1, 2 . . . . The values of d (k) and x (k) are defined at discrete points in time,
kτs, where τ s defines the sampling rate. The the discrete-time error signal is given by
e (k) = d (k)− x (k) . (6-71)
The discrete-time error signal may be incorporated into a discrete-time counterpart of the PID controller
given in Eq. 6-67. One possible form of a discrete-time PID controller is given by
c (k) = Kpe (k) +
kX
n=0
KI τs e (n) +Kd
e (k)− e (k − 1)
τs
. (6-72)
Other possibilities exist for the definitions of the summation which approximates the integral, and the finite
difference derivative in Eq. (6-72). For instance, it may be desirable to replace the integral-action term with
ME 349 Dynamic Systems Lab Manual 6-109
cz(t)
time 2  s  3  s  s 
 . . . 
Figure 6-3: Example output of a discrete-time controller with zero-order hold circuitry.
a trapezoidal rule,
1
2
kX
n=0
KI τs [e (n) + e (n− 1)] ; (6-73)
or, replace the derivative-action term with a higher-order approximation, such as
Kd
6τ s
[e (k) + 3e (k − 1)− 3e (k − 2)− e (k − 3)] , (6-74)
which would lead to some smoothing of the signal noise.
The output of the discrete-time controller, c (k), is a series of pulses defined at times kτ s. This type of
output is, in many instances, not sufficient to control the actuator. The series of pulses may not contain
sufficient energy to drive the actuator, or, the magnitude of the pulses required to effectively drive the
actuator may in fact damage some of its components. The typical solution to this limitation is to pass the
signal c (k) through a device known as a zero-order hold. The zero-order hold sets the signal level at c (k)
starting at the time kτ s and maintains it at that level up until time (k + 1) τ s, when the signal level is
changed to c (k + 1) . This process is repeated at each time step to give an output such as that denoted by
cz (t) in Figure 6-3.
The response of the actuator and plant, when they are subject to the output of the zero-order hold, will
in general be different from the response due to the continuous time signal, c (t) . A complete analysis of
the response is beyond the scope of this text. However, it can be shown that in some cases the response
is essentially the same. Consider a discrete-time control signal c (k) that approximates points along a sine
wave, y (t) = A sinΩt. Assume that the period of the sine wave, Ts = 2πΩ , is greater than the sampling
period, τs, by a large factor, i.e. Ts ≈ 100τs. Also, assume that the sampling (and control) frequency,
fs = 2πτs , is well above any natural frequencies of the plant or control system. If c (k) is passed through
a zero-order hold circuit to produce cz (t) , the output will appear as in Figure 6-4a. A discrete Fourier
transform of the zero-order hold output, sampled at a rate greater than 2fs, will give a spectrum such as the
one shown in Figure 6-4b. A predominant spike exists at the spectral estimate associated with the frequency
6-110 Schimel, Ding, Anderson and Grantham
 A 
frequency
 (b) 
 A 
time
 (a) 
 fs  2fs  3fs 
Figure 6-4: Effect of a zero-order hold on a sinusoidal signal, a) time-domain signal, b) frequency domain
signal.
of the sine-wave, 2πTs , and additional spikes associated with the sampling rate are present near frequencies of
fs, 2fs, 3fs, . . . , due to the step-like nature of the zero-order-hold output.
Consider the response if the zero-order hold output is input into a second-order linear system with a bode
diagram as shown in Figure 6-5. Since the system is linear, its response can be looked at as a superposition
of the responses at the frequencies Ω, fs, 2fs, 3fs, . . .. Reading from the bode diagram, the amplitude of the
response at the frequency Ω is given by
|X (iΩ)| = |G (iΩ)| |Y (iΩ)| ≈ (6db)A ≈ 2A (6-75)
Considering the relative amplitude of the largest spectral estimate associated with the sampling frequency
fs, the response due to this excitation is given by
|X (ifs)| = |G (ifs)| |Y (ifs)| ≈ (−62db) (0.032A) ≈ 2.5A× 10−5 (6-76)
The contribution to the response of the higher frequencies associated with the sampling rate, 2fs, 3fs, . . . ,
will be even smaller due to their smaller spectral magnitudes and the -40db per decade rolloff of the example
system. Since, for this particular system, the response due to the discrete nature of the zero-order hold is
small, the system response can be said to approximate that of system with a continuous time input.
The example above should not be viewed as a method for designing a control system. It serves only
to illustrate how the response of an appropriately designed discrete-time control system can approximate
that of a continuous-time control system. The design of a discrete-time control system often proceeds by
converting continuous time components of the system into discrete time counterparts, and then proceeding
with the analysis of the entire system in the discrete domain.
6.8 Linearization
Nonlinear differential equations may occur in models of plants, actuators or transducers. However, in order to
perform a linear controls analysis for these systems, it is a necessity that the differential equations used in the
ME 349 Dynamic Systems Lab Manual 6-111
frequency  fs  
|G
(i 
   
)|,
 d
b
Figure 6-5: Bode plot of a second order system showing the response magnitudes due to the driving frequency,
Ω, and the forcing frequency, fs.
analysis are linear. While it is sometimes possible to convert a nonlinear governing differential equation into
a linear one, in general it is necessary to form a linear approximation to the nonlinear differential equation
and then proceed with the controls analysis. The linearization process involves the use of a coordinate system
transformation and a truncated Taylor series approximation to the nonlinear differential equations. Consider
a general nonlinear differential equation
x˙ = f (x, u) (6-77)
where the variable u defines an input. An example of such an equation might be
x˙ = cx2 + dux (6-78)
where c and d are constants.
The linearization is made about some reference level of the state-variable, x, and the input, u. The
reference level often describes some setpoint that the system approaches at steady-state. For Eq. (6-77),
define the reference levels as x¯ and u¯ for the state-variable and input, respectively. The Taylor series
expansion of Eq. (6-77) about the reference points is given by
x˙ = f (x, u) = f (x¯, u¯) +
∂f(x¯, u¯)
x
(x− x¯) + ∂f(x¯, u¯)
u
(u− u¯) + · · · (6-79)
The higher-order terms (not shown) are truncated to form the linear approximation. A change of variables
is then made to transform the truncated Taylor series into the form of a linear differential equation. By
defining
y = x− x¯
v = u− u¯
a =
∂f(x¯, u¯)
x
b =
∂f(x¯, u¯)
u
y˙ = x˙− f (x¯, u¯)
, (6-80)
6-112 Schimel, Ding, Anderson and Grantham
and substituting these values into the truncated version of Eq. (6-79), the following linearized equation is
obtained:
y˙ = ay + bv (6-81)
This equation is suitable for control system analysis by the methods detailed above. Applying the lineariza-
tion approach to the example in Eq. (6-78) gives:
y = x− x¯
v = u− u¯
a = 2cx¯+ du¯
b = dx¯
y˙ = x˙− cx¯2 − dx¯u¯
, (6-82)
and linearized equation is exactly as in Eq. (6-81). Keep in mind that the linearization applies only in some
small neighborhood of the reference values x¯ and u¯, and thus the control constants selected may only give
the designed performance in this neighborhood. However, it is often possible to select control constants for
a system that result in appropriate response over a wide range of reference values.
6.9 Apparatus, Measurements, and Calculations
This laboratory will consist of controlling the fluid level in a small water tank with a computerized control
system. A drawing of the mechanical portion of the apparatus appears in Figure 6-6, and a schematic of
the overall system, including controller appears in Figure 6-7. The system consists of a water tank that is
supplied with water by the variable speed motor-pump system used in Laboratory #5. A pressure transducer
at the bottom of the tank generates an electric signal proportional to the fluid level in the tank. The signal
from the pressure transducer is amplified by a signal conditioner, read into the computer through an analog-
to-digital convert, and This is then converted into the measured fluid level by a simple conversion constant.
The measured fluid level is then subtracted from a user-input desired fluid level to generate the error signal
which is input into a digital controller that has user-specified control constants. The digital controller is
actually a computer program that calculates the control signal from the error signal using an equation similar
to Eq. (6-72). Once the control signal has been calculated, it is output to the variable speed motor and
pumping system described in Chapter 5 via a digital-to-analog converter. The digital-to-analog converter
performs the service of a zero-order hold circuit.
The magnitude of the control signal determines how fast the pump operates, and thus, how fast the tank
fills. The magnitude of the control signal is computed from the control equation, and therefore depends on
existing and previous fluid levels in the tank. As a rule of thumb, the pumping rate will be fast when the
fluid level is far below the desired level, and it will slow as the fluid level stabilizes at or near the desired
fluid level. However, stability is not guaranteed. A set of control constants that result in a stable set of
eigenvalues for the overall system must be selected for the fluid level to reach a stable value. Given a set
of control constants that result in a stable system, the fluid level can either reach a level equivalent to the
desired fluid level if a method such as PI-control is used, or there may be some steady-state error present if
only proportional control is used.
For the apparatus used in this laboratory, the pump cannot be operated in reverse, i.e. it can only fill the
tank. However, a 1/4-inch drain tube allows water to escape from the bottom of the tank. Thus, starting
from a tank-empty condition, the controller is employed to fill the tank to some constant level, and maintain
the level in the presence of fluid draining from the tank bottom. If the selected control action causes the
tank to overfill, then at some point the control signal may become negative (instructing the pump to operate
in reverse). However, since the pump will not operate in reverse, the rate at which the tank empties is
determined only by the drain rate.
Governing Differential Equations
Figure 6-8 shows a schematic of the water tank. The variable h (t) is the height of the water in the tank
and q (t) is the water flow rate into the tank. From Bernoulli’s equation, the velocity of the fluid exiting the
ME 349 Dynamic Systems Lab Manual 6-113
tank
drain
water
 tank
inlet
tube
 flow
meter
  pressure
transducer
supply
 tube
fluid reservoir
Figure 6-6: Pump and water tank system.
Figure 6-7: A schematic of the liquid-level control apparatus.
6-114 Schimel, Ding, Anderson and Grantham
tank through the drain tube is given by
vo = (2gh)
1
2 . (6-83)
Assuming that the cross-sectional area of the drain tube is Ao, and assuming that the density is ρ, conser-
vation of mass gives
m˙ = ρq − ρAo (2gh)
1
2 . (6-84)
Since the tank has a constant cross-sectional area the mass-flow rate, m˙, is a function of the rate of change
in the water height:
m˙ = ρAh˙. (6-85)
Substituting Eq. (6-79) into Eq. (6-78) and rearranging gives
h˙+
Ao
A
(2gh)
1
2 =
1
A
q. (6-86)
For this system, the tank acts as the plant, and the pump and motor system act as the actuator. The
control signal into the actuator is the voltage output from the digital-to-analog converter, Vin, as defined for
Eq. (5-72). Under some circumstances it may be necessary to model the portion of the system between the
control system and the inlet to the tank with a differential equation. The derivation leading up to Eq. (5-91)
shows that the equation of motion for the pump and motor is at least a second-order model. Additionally,
a complete dynamic model of this system would include a model relating the voltage out of the generator,
Vout in Eq. (5-91), to the pump flow rate q (t) in Eq. (5-91). However, the time constants associated with
part of the system modeled in Chapter 5 are very small in comparison to the time constant associated with
filling the tank. Therefore, for this application, the actuator will be modeled with a simple constant, Ka,
and the relationship between the control voltage and the flow rate into the tank is given by
q = KaVin. (6-87)
Substituting this relationship into Eq. (6-80) gives
h˙+
Ao
A
(2gh)
1
2 =
Ka
A
Vin; (6-88)
the governing differential equation for the actuator and plant of the liquid-level control system. Since the
pump motor cannot operate in reverse, Eq. (6-88) holds only for cases when the control input Vin is greater
than zero. For the case when Vin < 0, the governing differential equation reverts to
h˙+
Ao
A
(2gh)
1
2 = 0, (6-89)
and the equilibrium point (h˙ = 0) for this case is at h = 0. Thus, while the control input is negative, the
tank simply drains up to a point where control input becomes positive again. For the case of proportional
control and the initial conditions used in this experiment, the control input is always positive. However, a
negative control input can be observed experimentally for PI and PID-control.
Linearization of the Governing Differential Equation
Equation (6-82) defines the relationship between the control input, Vin, and the fluid level in the tank, h (t).
The equation is nonlinear and, before proceeding with a linear control system analysis, it must be linearized.
By rewriting Eq. (6-88) as
h˙ = −Ao
A
(2gh)
1
2 +
Ka
A
Vin, (6-90)
the linearization can proceed as defined in Eqs. (6-79) through (6-81). Define the reference fluid level as
h¯. The reference fluid level is equivalent to the user-input desired fluid level. The reference input, V¯in, can
be taken as zero since the input is already linear in form. Performing the truncated Taylor series expansion
about h¯ gives
h˙ ≈ −Ao
A
¡
2gh¯
¢ 1
2 − Aog
A
¡
2gh¯
¢− 12 ¡h− h¯¢+ Ka
A
Vin. (6-91)
ME 349 Dynamic Systems Lab Manual 6-115
q(t)
h(t)
Vo
Figure 6-8: Schematic of the water tank.
Defining
y = h− h¯
v = Vin
a = −AogA
¡
2gh¯
¢− 12
b = KaA
y˙ = h˙+ AoA
¡
2gh¯
¢ 1
2
, (6-92)
gives the linearized version of the governing differential equation:
y˙ = ay + bVin. (6-93)
This equation is linear in form, and hence it is suitable for the analysis by the linear control methods
described above.
Equation (6-93) applies only in cases where the control input, Vin, is greater than zero. Since the pump
motor cannot operate in reverse, keep in mind that the system behavior will revert to Eq. (6-89) when the
control input becomes negative.
Proportional Control
For the linearized system, the transfer function between input and output is given by
Ga (s)Gp (s) =
Y (s)
Vin (s)
=
b
s− a. (6-94)
Application of Eq. (6-10) for the case of proportional control gives
G (s) =
Kp bs−a
1 +Kp bs−a
=
Kpb
s− a+Kpb
. (6-95)
6-116 Schimel, Ding, Anderson and Grantham
By Eq. (6-12), the eigenvalue of the linearized system is real and negative, and thus, stable for proportional
control constants that satisfy
Kp >
a
b
= −Aog
Ka
¡
2gh¯
¢− 12 . (6-96)
The final value theorem shows that the steady-state error for the linearized system subject to a step input
h¯ at time t = 0 is
lim
t→∞
e (t) =
µ
1− Kpb−a+Kpb
¶
h¯ =
⎛
⎝1− KpKa
Aog
¡
2gh¯
¢− 12 +KpKa
⎞
⎠ h¯. (6-97)
The linearized model predicts that the steady-state error approaches zero as Kp becomes large.
PI-control
Application of Eqs. (6-10), (6-38), and (6-94) yield the transfer function for the linearized PI-control system
G (s) =
Kps+KI
s
b
s−a
1 +
Kps+KI
s
b
s−a
=
(Kps+KI) b
s2 + (Kpb− a) s+KIb
. (6-98)
The final value theorem shows that the steady-state error is zero for the linearized PI-control system subject
to a step input h¯ at t = 0:
lim
t→∞
e (t) =
µ
1− KIb
KIb
¶
h¯ = 0.
The eigenvalues of the characteristic equation of Eq. (6-98) are
λ1,2 =
− (Kpb− a)±
q
(Kpb− a)2 − 4KIb
2
. (6-99)
The linearized PI-control system is stable for cases where Kp > ab and KI > 0. If
KI >
(Kpb− a)2
4b
=
³
KpKa +Aog
¡
2gh¯
¢− 12´2
4Ka
(6-100)
then the eigenvalues of the linearized system are complex valued, indicating that the behavior of the actual
system may be oscillatory.
6.10 Laboratory #6 Procedures
Proportional Control
Using the liquid-level control system and the proportional control constants given in class, measure the
time constants and steady-state error each proportional control constant as instructed. Tabulate the results.
Note that the time constant for a linear system exhibiting first-order behavior is equivalent to any of the
the following: the 37% settling time, 1/2 the 14% settling time, or 1/4 the 2% settling time for a system
exhibiting first-order behavior.
Table 6-1. Proportional Control Performance Criteria
KP final height (in.) τ (sec.) 2τ (sec.) 4τ (sec.) ess (in.)
ME 349 Dynamic Systems Lab Manual 6-117
PI-control
Using the liquid-level control system and the PI-control constants supplied in class, measure the 100% rise
time, cycle time, 10% settling time, and maximum overshoot as shown in Figure 6-2. Tabulate the results.
Table 6-2. PI-Control Performance Criteria
Kp KI tr (sec.) xm (in.) tc (sec.) ts (sec.)
Chapter 7
State Space Modeling
The derivation of a mathematical model for a physical system can often lead to one or more time-varying
equations, that is, a dynamic model. Accurately describing the behavior of a physical system as it evolves
through time may require a multi-variable dynamic model with multiple inputs and multiple outputs. Any
of the inputs or outputs to the physical system may be uncertain. State-space methods are an approach
to dynamic modeling of physical systems, and these methods are highly applicable to modeling the above
conditions. State-space methods are a matrix oriented approach that directly utilize a system’s dynamical
equations, as opposed to indirectly as in transfer function analysis. State-space methods also are particularly
useful in control system modeling and analysis.
The concept of a system’s state is central to state-space methods, and is defined as:
Definition 1 (State) a set of physical quantities which describe the dynamic behavior of a system and,
given any inputs to the system, allow the future behavior to be completely determined.
State-space methods have been in use in the US since the late 1950’s where they were introduced by
R.E. Kalman. The methods originated in Russia in the early 1950’s and were later translated to French and
English. State-space methods have since found applications in many technologically advanced applications
including the US Apollo program and later space programs.
7.1 State-Space Theory
Model Format
A state-space model defines a system model in terms of a set of first-order differential equations. In the
simplest form a state-space model contains a set of state variables and a coefficient matrix. Consider a
general N thx order differential equation, such as a linear inhomogeneous ordinary differential equation with
time derivatives d
iz
dti , coefficients ai, and input function u,
dNx z
dtNx
+ aNx−1
dNx−1z
dtNx−1
+ · · ·+ a2 d
2z
dt2
+ a1
dz
dt
+ aoz = bu. (7-1)
Equation (7-1) can be expressed in state space form by defining a set of first order differential equations,
x1 = z
x2 = x˙1 = dzdt
x3 = x˙2 = d
2z
dt2
...
...
...
xn = x˙n−1 = d
Nx−1z
dtNx−1
x˙n = x˙n = d
Nxz
dtNx = −a1x1 − a2x2 − · · ·− aNxxNx + bu.
(7-2)
7-118
ME 349 Dynamic Systems Lab Manual 7-119
If the differential equations are linear, the state-space equations can be represented in matrix form. For Eq.
(7-1) the matrix form of the state-space equations is written as
⎧
⎪⎪⎪⎪⎪⎪⎨
⎪⎪⎪⎪⎪⎪⎪⎩
x˙1
x˙2
·
·
·
·
x˙Nx
⎫
⎪⎪⎪⎪⎪⎪⎬
⎪⎪⎪⎪⎪⎪⎪⎭
=
⎡
⎢⎢⎢⎢⎢⎢⎢⎢⎣
0 1 0 · · · 0
0 0 1 0 · · 0
· · · · ·
· · · · ·
· · · 0
0 0 · · · 0 1
−a1 −a2 · · · · −aNx
⎤
⎥⎥⎥⎥⎥⎥⎥⎥⎦
⎧
⎪⎪⎪⎪⎪⎪⎨
⎪⎪⎪⎪⎪⎪⎪⎩
x1
x2
x3
·
·
·
xNx
⎫
⎪⎪⎪⎪⎪⎪⎬
⎪⎪⎪⎪⎪⎪⎪⎭
+ b
⎧
⎪⎪⎪⎪⎪⎪⎨
⎪⎪⎪⎪⎪⎪⎪⎩
0
·
·
·
·
0
u
⎫
⎪⎪⎪⎪⎪⎪⎬
⎪⎪⎪⎪⎪⎪⎪⎭
(7-3)
The matrix format can also be used to express general coupled linear differential equations with multiple
deterministic inputs, ui, and uncertain inputs, vj, for a time invariant system:
⎧
⎪⎪⎨
⎪⎪⎪⎩
x˙1
x˙2
...
x˙Nx
⎫
⎪⎪⎬
⎪⎪⎪⎭
=
⎡
⎢⎢⎢⎣
a11 a12 · · · a1Nx
a21 a22 · · · a2Nx
...
...
. . .
...
aNx1 aNx2 · · · aNxNx
⎤
⎥⎥⎥⎦
⎧
⎪⎪⎨
⎪⎪⎪⎩
x1
x2
...
xNx
⎫
⎪⎪⎬
⎪⎪⎪⎭
+
⎡
⎢⎢⎢⎣
b11 b12 · · · b1Nu
b21 b22 · · · b2Nu
...
...
. . .
...
bNx1 bNx2 · · · bNxNu
⎤
⎥⎥⎥⎦
⎧
⎪⎪⎨
⎪⎪⎪⎩
u1
u2
...
uNu
⎫
⎪⎪⎬
⎪⎪⎪⎭
+
⎡
⎢⎢⎢⎣
r11 r12 · · · r1Nv
r21 r22 · · · r2Nv
...
...
. . .
...
rNx1 rNx2 · · · rNxNv
⎤
⎥⎥⎥⎦
⎧
⎪⎪⎪⎨
⎪⎪⎪⎩
v1
v2
...
vNv
⎫
⎪⎪⎪⎬
⎪⎪⎪⎭
(7-4)
The concepts of deterministic and uncertain inputs will become more clear in an example presented below.
The linear state equation can also be represented in compact notation as:
x˙ = Ax+Bu+Rv (7-5)
where x is an Nx-dimensional state vector, A is an Nx-by Nx-dimensional state coefficient matrix; u is an
Nu-dimensional deterministic input vector, B is an Nx-by Nu-dimensional matrix; v is an Nv-dimensional
uncertain input vector, R is an Nx-by Nv-dimensional matrix; and x˙ is a vector of time derivatives of x.
The model interfaces with its environment through a set of input equations, or input vectors, as shown
in Eq. (7-5) and a set of output equations, or the output vector. The output vector is written in compact
form as:
y = Cx+Du+Ev+ Sw (7-6)
where y is an Ny-dimensional state vector; x, u and v are as before, and C, D and E are their respec-
tive Ny-by Nx, Ny-by Nu, and Ny-by Nv-dimensional coefficient matrices; and w is an Nw-dimensional
measurement uncertainty vector with Nw-by Ny-dimensional coefficient matrix.
Equations (7-5) and (7-6) constitute the input-output form of the linear state-space model of a system.
If the model is time-varying the coefficient matrices in Eqs. (7-5) and (7-6) become functions of time. Note
that the dimensions of the vectors: y, u, v and w; need not equal the model order, which is the dimension,
Nx, of the state vector. Also, a given set of state variables is not unique for a particular system. For
instance, an alternate set of state vectors can be obtained from a linear combination of the original set. The
application of state space techniques is not dependent on the model dimension. However, the order of the
model for a particular system is unique.
State Transition Matrix
A solution to the state equation can always be obtained for linear systems. For the homogenous, or unforced
case
x˙ = Ax, (7-7)
and the solution can be expressed as
x(t) = x◦ e
At (7-8)
7-120 Schimel, Ding, Anderson and Grantham
where x◦ is a vector specifying the state at t=0, and eAt is the matrix exponential function. The matrix
exponential is defined as
eAt ≡ I+At+A2 t
2
2!
+A3
t3
3!
+ · · · (7-9)
If the state of the linear system is known at some arbitrary time, τ , the state can be determined at time, t,
with the relationship
x(t) = x(τ) eA(t−τ) (7-10)
The matrix exponential with the argument (t−τ) in Eq. (7-10) has special significance, it is called the state
transition matrix and is defined for a linear system as
Φ(t, τ) ≡ eA(t−τ). (7-11)
The state transition matrix shows that in the absence of inputs the state at time, t, depends only on the
state at an arbitrary time, τ , or alternately, the current state of an unforced system depends only on the
previous state.
For a forced system, the solution to the linear state-space equations at time t given a knowledge of the
input vector from time τ to time t, and state vector at time τ , is
x(t) = x(τ) eA(t−τ) +
Z t
τ
eA(t−τ) [Bu(α) +Rv(α)] dα (7-12)
The output vector, y(t), is obtained by substituting Eqs. (7-10) or (7-12) into Eq. (7-6).
Stability
The stability of an unforced linear state space model can be determined through an eigenvalue analysis of
the coefficient matrix A. The eigenvalues are the zeros of the characteristic equation
det [A−λI] =
¯¯¯¯
¯¯¯¯
¯
a11 − λ a12 · · · a1Nx
a21 a22−λ · · · a2Nx
...
...
. . .
...
aNx1 aNx2 · · · aNxNx − λ
¯¯¯¯
¯¯¯¯
¯ = 0 (7-13)
or
det [A−λI] = λNx + hNx−1 λNx−1 + · · ·+ h1 λ+ h◦ = 0 (7-14)
The analysis produces Nx eigenvalues, λi, some or all of which can be complex valued. From a standpoint of
stability, the sign on the real part of each eigenvalue determines the behavior of the system. The conclusions
which can be drawn from an eigenvalue analysis of a linear state-space model can be summarized as below.
The indices i and j range from 1 to Nx and Re(λi) represents the real part of the ith eigenvalue.
• Re(λi) < 0 for all i ⇔ The model is asymptotically stable.
• Re(λi) > 0 for one or more i ⇔ The model is unstable.
• Re(λi) = 0 for one or more i, but if there is more than one Re(λi) = 0, all such λi are unequal; and any
remaining nonzero Re(λj) are such that Re(λj) < 0 ⇔ The model is stable but not asymptotically
stable.
• Re(λi) = Re(λj) = 0 for some i 6= j; and some of the eigenvalues with zero real parts are such that,
λi = λj ⇔ The model is unstable.
ME 349 Dynamic Systems Lab Manual 7-121
Converting State-Space Models to Transfer Function Models
The input-output form of the linear state-space model of Eqs. (7-5) and (7-6) is, in the absence of input and
measurement uncertainty, directly convertible to a transfer function model. The transfer function model is
obtained by taking the Laplace transform of Eqs. (7-5) and (7-6), setting the initial conditions to zero, and
solving for the input-output relationship:
Y(s) =G(s)U(s) (7-15)
where
G(s) = C [sI−A]−1B+D (7-16)
Equation (7-16) is the Ny-by Nu-dimensional transfer matrix.
Nonlinear State-Space Models
It is often the case that the behavior of a physical system cannot be expressed in terms of a linear model. A
state-space model can still be derived from the governing differential equations, however, it may no longer
be possible to represent the model in matrix format. A generalized format for nonlinear models does not
exist, they can contain a rich variety of terms or functions. As an example consider the following pair of
coupled nonlinear differential equations
d2w
dt2 + a3
dw
dt
dz
dt + a2(w − 1) z = b1u1
d2z
dt2 + a1
dz
dt
dw
dt + ao(z − 1)w = b2u2.
(7-17)
The state vector might contain
x1 = w, x2 =
dw
dt
, x3 = z, x4 =
dz
dt
. (7-18)
The state-space model would then be
x˙1 = x2
x˙2 = −a3x2x4 − a2(x1 − 1)x3 + b1u1
x˙3 = x4
x˙4 = −a1x4x2 − ao(x3 − 1)x1 + b2u2.
(7-19)
The input-output form of the linear and nonlinear state-space models can be represented in compact
notation as
x˙ = f(x,u) (7-20)
and
y = g(x,u) (7-21)
where f(x,u), and output g(x,u) represent the vector functions
f(x,u) =
⎧
⎪⎪⎪⎨
⎪⎪⎩
f1(x,u)
f2(x,u)
...
fNx(x,u)
⎫
⎪⎪⎪⎬
⎪⎪⎭
and g(x,u) =
⎧
⎪⎪⎪⎨
⎪⎪⎩
g1(x,u)
g2(x,u)
...
gNy(x,u)
⎫
⎪⎪⎪⎬
⎪⎪⎭
. (7-22)
each of which may depend on the state and input vectors.
7-122 Schimel, Ding, Anderson and Grantham
Equilibrium Point Analysis
Equilibrium points are useful when characterizing the stability of a general state-space model and are of
particular application to control system design. By selecting a suitable controller it is possible to transform
a unstable or neutrally stable system into an asymptotically stable system or create a designer selected
equilibrium operating point.
Equilibrium points can obtained by assuming that the model rests in a steady-state condition. Mathe-
matically this can be achieved by setting the left-hand-side of Eq. (7-5) and the input to zero. Equilibrium
points can then obtained by algebraically solving the resulting equations for any state vector satisfying the
zero left-hand-side condition. For the compact form represented by Eq. (7-20) this means solving
x˙ = 0 = f(x,0) (7-23)
for the vector x. Unforced linear systems for which all eigenvalues are nonzero have one and only one isolated
equilibrium point, at x = 0. Other possibilities include a continuous line of equilibria when one eigenvalue
of a linear system is zero, a continuous plane of equilibria when two eigenvalues of a linear system are zero,
and multiple isolated equilibrium points for some nonlinear systems.
As an example consider the nonlinear state-space model of Eq. (7-19) for an unforced case (u1 and
u2 = 0). Setting the left-hand-side to zero gives:
x˙1 = 0 = x2
x˙2 = 0 = −a3x2x4 − a2(x1 − 1)x3
x˙3 = 0 = x4
x˙4 = 0 = −a1x4x2 − ao(x3 − 1)x1
(7-24)
Solving these equations gives two equilibrium points:
x =
⎧
⎪⎪⎨
⎪⎩
0
0
0
0
⎫
⎪⎪⎬
⎪⎭
and, x =
⎧
⎪⎪⎨
⎪⎩
1
0
1
0
⎫
⎪⎪⎬
⎪⎭
(7-25)
Linearization of Nonlinear State-Space Models
It is often desirable to control or analyze a nonlinear system about some reference operating level. The terms
of the input-output form of the nonlinear state-space model, Eq. (7-20), can be modeled as the sum of a
reference operating function, and a small perturbation function. That is, for the inputs,
U(t) = U¯(t) + u(t) (7-26)
V(t) = V¯(t) + v(t), (7-27)
and measurement uncertainty
W(t) = W¯(t) +w(t) (7-28)
the solution to the state-space model is of the form
X(t) = X¯(t) + x(t) (7-29)
Y(t) = Y¯(t) + y(t) (7-30)
Where the terms with the bar over-strike denote reference functions, and the lower case terms denote small
perturbations away from reference. In control system analysis, it is often the case that the reference functions
are reduce to fixed points in state-space about which the system operates.
ME 349 Dynamic Systems Lab Manual 7-123
Linearization and subsequent linear analysis will often produce useful information. Linearization of Eq.
(7-20) for is achieved utilizing the first term of a Taylor series expansion along the reference functions
x˙ =
£
∂f
∂X
¤
X¯
¡
X− X¯
¢
+
£
∂f
∂U
¤
U¯
¡
U− U¯
¢
+
£
∂f
∂V
¤
V¯
¡
V− V¯
¢
y =
h
∂g
∂X
i
X¯
¡
X− X¯
¢
+
h
∂g
∂U
i
U¯
¡
U− U¯
¢
+
h
∂g
∂V
i
V¯
¡
V− V¯
¢
+
h
∂g
∂W
i
W¯
¡
W− W¯
¢ (7-31)
Or applying equations Eqs. (7-26) through (7-30)
x˙ =
£
∂f
∂X
¤
X¯
x+
£
∂f
∂U
¤
U¯
u+
£
∂f
∂V
¤
V¯
v
y =
h
∂g
∂X
i
X¯
x+
h
∂g
∂U
i
U¯
u+
h
∂g
∂V
i
V¯
v +
h
∂g
∂W
i
W¯
w
(7-32)
where the partial derivative of a general n-dimensional vector function,
Ψ(Z) =
⎡
⎢⎢⎢⎣
Ψ1(Z)
Ψ2(Z)
...
ΨNx(Z)
⎤
⎥⎥⎥⎦ , (7-33)
evaluated along a reference function Z¯, and is given by:
∂Ψ
∂Z
=
⎡
⎢⎢⎢⎣
∂Ψ1
∂Z1
∂Ψ1
∂Z2
· · · ∂Ψ1∂Zn
∂Ψ2
∂Z1
∂Ψ2
∂Z2
· · · ∂Ψ2∂Zn
...
...
. . .
...
∂Ψn
∂Z1
∂Ψn
∂Z2
· · · ∂Ψn∂Zn
⎤
⎥⎥⎥⎦
Z¯
. (7-34)
7.2 Applications of State-Space Models
Linearization of an Automobile Cruise Control Model
An automobile being operated under cruise control is an example of a very commonplace system that can
be readily modeled by state-space methods. Once the driver has set a desired speed the controller attempts
to maintain this speed in the presence of changing road conditions. Development of one possible state-space
model proceeds as follows. Consider an automobile depicted in Figure 7-1 with mass m, traveling at velocity
X, along a road grade with angle V, is subject to a resistive wind force that varies with the square of
vehicle velocity, kX2, rolling resistance, µmg cosV, and gravitational force, mg sinV. The terms k, g and
µ are coefficient of drag, gravitational acceleration and coefficient of resistance to rolling respectively. The
forces are countered, in the sense drawn, by a motor-generated force U, delivered at the drive wheels of
the automobile. It is assumed that the motor can generate positive forces through the application of engine
power as well as negative forces through engine compression. Application of Newton’s second law results in
the following model:
mX˙ = −kX2−µmg cosV−mg sinV+U. (7-35)
The model is nonlinear and consists of a single state variable X, its time derivative X˙, plus known input
U. Additionally, the road grade angle V, can be modeled as an undetermined input. Simple algebraic
manipulation gives the state-space form of the model
X˙ = − 1
m
kX2 +
1
m
U−µg cosV−g sinV. (7-36)
Linearization of this model can be performed about preselected constant reference values for vehicle velocity
X¯, and road grade angle V¯. Setting X˙ =0 and substituting for reference values in Eq. (7-36) gives the
reference motor-generated force or control force necessary to maintain constant velocity under reference
conditions,
U¯ =kX¯2+µmg cos V¯+mg sin V¯. (7-37)
7-124 Schimel, Ding, Anderson and Grantham
  V  
  X 
  U  
 mg 
  N 2  
  N1  
 kX2
  mg  sin   V  
  mg  cos  V  
 (a) 
 (b) 
Figure 7-1: Model automobile for the cruise control example, (a) general configuration, (b) free-body diagram.
A linearized model can be used to predict the system response to minor perturbations about the reference
conditions. Application of Eq. (7-32) to Eq. (7-36) gives the linearized state-space model
x˙ = −
µ
2k
m
X¯
¶
x+
1
m
u+ (µg sin V¯)v− (g cos V¯)v.
This equation matches exactly the linear form of the state-space equation given in Eq. (7-5), where the
coefficient matrices for this model are one-by-one and given by
A = −
³
2kX¯
m
´
B = 1m
R = (µg sin V¯)− (g cos V¯).
(7-38)
To complete the state-space model in input-output format, an output equation must be specified. In this
example, a likely output would be the vehicle velocity. The output equation is then given by:
y = x (7-39)
which matches the form of Eq. (7-6) with C = 1 and, D,E and S equal to zero.
A State-Space Model of an Armature Controlled DC Motor
Consider an idealized model of an armature controlled DC motor.. The torque, Γ, generated by the Motor
is modeled as proportional to the armature current, IA, or
Γ = kΓ IA (7-40)
For this laboratory, the armature current will supplied by a programmable power supply with current, IA,
proportional to an input control voltage, Vin. Therefore,
IA = KAVin, (7-41)
ME 349 Dynamic Systems Lab Manual 7-125
where KA is a constant. The equation of motion for the motor shaft is given by
J θ¨ = Γ− Γex − cθ˙ (7-42)
where J is the mass moment of inertia of the rotor plus any attachments, Γex is the sum of any external
torques applied to the shaft, and cθ˙ accounts for sources of energy dissipation such as friction in the armature
bearings. Substituting Eq. (7-40) into Eq. (7-42) gives the governing differential equation
J θ¨ = kΓKAVin − Γex − cθ˙ (7-43)
Assigning ½
x1
x2
¾
=
½
θ
θ˙
¾
and putting Eq. (7-43) in state-space form yields½
x˙1
x˙2
¾
=
½
x2
kΓKAVin − Γex − cx2
¾
. (7-44)
In this case the reference function defined in Eq. (7-29) reduces to zero, and hence the uppercase notation
has been dropped. The output equation vector corresponding to Eq. 7-44 might be modeled as
y =
∙
1 0
0 1
¸½
x1
x2
¾
. (7-45)
7.3 Apparatus
Motor-Pendulum Model
For this laboratory the control of an inverted pendulum as depicted in Figure 7-2 will be investigated. The
pendulum system consists of a slender rod of length 2l and mass m to the output shaft of an armature-
controlled DC motor. The voltage across a rotary potentiometer is used to measure the angular position
of the pendulum. A free-body diagram of this system is shown in Figure 7-3. The control objective is to
balance the rod in an inverted position by adjusting the current input to the motor. The total moment of
inertia of the system is given by
J = Jm +
4
3
ml2 (7-46)
where Jm is the moment of inertia of the motor shaft and rotor. The rod, which acts as a pendulum, will
exert an external torque related to its angular displacement θ from the upright position
Γext = −mgl sin θ = −mgl sinx1. (7-47)
Applying Eq. (7-44) gives the state-space equations of motion for the motor-pendulum system,½
x˙1
x˙2
¾
=
½
x2
kΓKAVin +mgl sinx1 − cx2
¾
. (7-48)
In the absence of the input Vin, Eq. (7-48) becomes½
x˙1
x˙2
¾
=
½
x2
mgl sinx1 − cx2
¾
. (7-49)
An equilibrium point analysis shows that the stationary inverted position,½
x1
x2
¾
=
½
0
0
¾
, (7-50)
7-126 Schimel, Ding, Anderson and Grantham
Pendulum Rod
DC Servomotor
     Rotary
Potentiometer
Figure 7-2: The inverted pendulum apparatus.p
is an equilibrium point. Linearization about this equilibrium point by Eq. (7-32) gives½
x˙1
x˙2
¾
=
∙
0 1
mgl
J −
c
J
¸½
x1
x2
¾
. (7-51)
The linearized system has the eigenvalues
λ1,2 =
−c±
p
c2 + 4mgl
2J
. (7-52)
The eigenvalues are real, and since one of the eigenvalues is greater than zero, the pendulum system in the
inverted position is unstable. Therefore, some form of control input must be used to stabilize the system.
Proportional Control
For this laboratory, the voltage Vin is the output of a linear controller. A schematic diagram of the inverted
pendulum and control system is shown in Figure 7-4. For the case of proportional control, the control signal
is given by
Vin = Kpe (t) , (7-53)
where Kp is a constant, and e (t) is an error signal. The error signal for this laboratory is defined as
e (t) = θd − θ = θd − x1 (7-54)
where the angle θd is the desired angular position of the pendulum. The desired angular position will be
taken to be θd = 0, or the inverted position. Therefore, for the case of proportional control, Eq. (7-48)
becomes ½
x˙1
x˙2
¾
=
½
x2
−KΓKAKp
J x1 +
mgl
J sinx1 −
c
J x2
¾
. (7-55)
ME 349 Dynamic Systems Lab Manual 7-127
 mg 
  l 
 
  l 
 
Figure 7-3: A free-body diagram of the inverted pendulum system.
Equilibrium points exist at ½
x1
x2
¾
=
(
mgl
KΓKAKp
sinx1
0
)
. (7-56)
Thus, ½
x1
x2
¾
=
½
0
0
¾
(7-57)
is an equilibrium point, and if the inequality
mgl
KΓKAKp
> 1 (7-58)
is satisfied then two additional equilibrium points exist at equal but opposite angles from the upright position.
Linearization of the system about the upright position gives½
x˙1
x˙2
¾
=
∙
0 1
mgl−KΓKAKp
J −
c
J
¸½
x1
x2
¾
. (7-59)
This system has the eigenvalues
λ1,2 =
−c±pc2 + 4 (mgl −KΓKAKp)
2J
, (7-60)
and is therefore stable for
Kp >
mgl
KΓKA
. (7-61)
If
Kp >
c2 + 4mgl
4KΓKA
(7-62)
then the eigenvalues are complex conjugates and the response is oscillatory.
7-128 Schimel, Ding, Anderson and Grantham
 Vin  
Figure 7-4: The inverted pendulum control apparatus.
PD-Control
For the case of PD-control,
Vin = Kpe (t) +KD
de
dt
= Kp x1 +KD x2, (7-63)
and state space model is½
x˙1
x˙2
¾
=
½
x2
−KΓKAKp
J x1 +
mgl
J sinx1 −
c+KΓKAKD
J x2
¾
. (7-64)
This system retains the same equilibrium points as for the case of proportional control. Linearizing about
the upright position gives ½
x˙1
x˙2
¾
=
∙
0 1
mgl−KΓKAKp
J −
c+KΓKAKD
J
¸½
x1
x2
¾
. (7-65)
The eigenvalues of the linearizes system with PD-control are
λ1,2 =
− (c+KΓKAKD)±
q
(c+KΓKAKD)
2 + 4 (mgl −KΓKAKp)
2J
. (7-66)
The PD-controlled system is stable for the same range of values of Kp as for the proportional-control case,
given that KD > 0. As can be seen in Eqs. (7-65) and (7-66), increasing the derivative control constant,
KD, is equivalent to increasing the damping constant, c. Thus, positive values of KD serve to increase the
total amount of damping in the system.
If electrical noise or other sources of noise are present in the error signal, then care must be taken not
to use values of KD that are too large. This requirement is due to the fact that noise may cause sudden
changes in the error signal level. The effects of noise are most pronounced when derivative control is applied
since derivative-control is a function of the time-derivative of the error signal. The sudden changes in the
ME 349 Dynamic Systems Lab Manual 7-129
error signal slope associated with the noise are generally much more pronounced than changes in the error
signal level alone. Therefore, if KD is unsuitably large, the motion of the system may become unpredictable
due to throughput of the error signal time derivative to the system.
The presence of high-frequency noise presents another problem when using continuous time control meth-
ods that incorporate derivative control. The problem, which is associated with noise amplification, may be
remedied by application of pseudo-derivative feedback control. See Eqs. (6-63) through (6-66) and the
surrounding discussion for information on this topic.
PID-Control
For the case PID-control
Vin = −Kpe (t)−KI
Z t
0
e (τ) dτ −KD
de
dt
(7-67)
Substituting Eqs. (7-47) and (7-55) into Eq. (7-43) gives the governing differential equation for the case of
PID-control
J θ¨ = −KΓKA
µ
KP θ +KI
Z t
0
θ (τ) dτ +KDθ˙
¶
+mgl sin θ − cθ˙. (7-68)
Assigning the state variables ⎧
⎨
⎩
x0
x1
x2
⎫
⎬
⎭ =
⎧
⎨
⎩
R t
0
θ (τ) dτ
θ
θ˙
⎫
⎬
⎭ (7-69)
and converting Eq. (7-68) into the general nonlinear state-space form of Eq. (7-20) gives
⎧
⎨
⎩
x˙0
x˙1
x˙2
⎫
⎬
⎭ =
⎧
⎨
⎩
x1
x2
−KΓKAKIJ x1 −
KΓKAKP
J x2 +
mgl
J sinx2 −
c+KΓKAKD
J x3
⎫
⎬
⎭ . (7-70)
An equilibrium point analysis for this system shows that an equilibrium point exists at
⎧
⎨
⎩
x0
x1
x2
⎫
⎬
⎭ =
⎧
⎨
⎩
0
0
0
⎫
⎬
⎭ , (7-71)
which corresponds to the pendulum being in a stationary, inverted position. Application of Eq.(7-32) yields
the following set of linearized state-space equations
⎧
⎨
⎩
x˙0
x˙1
x˙2
⎫
⎬
⎭ =
⎡
⎣
0 1 0
0 0 1
−KΓKAKIJ
mgl−KΓKAKP
J −
c+KΓKAKD
J
⎤
⎦
⎧
⎨
⎩
x0
x1
x2
⎫
⎬
⎭ . (7-72)
Including integral-control in the state-space model has increased the order of model by one. A purely
theoretical analysis of the third-order linearized model produces a set of equations that are not easily analyzed
by symbolic algebra. However, if numeric values are introduced for the constants in Eq. (7-72), various
behaviors, including both stable oscillatory and stable nonoscillatory response, can be observed. For example,
consider the case when the control constants are selected so that Eq. (7-72) becomes
⎧
⎨
⎩
x˙0
x˙1
x˙2
⎫
⎬
⎭ =
⎡
⎣
0 1 0
0 0 1
−6 −11 −6
⎤
⎦
⎧
⎨
⎩
x0
x1
x2
⎫
⎬
⎭ . (7-73)
The eigenvalues for this case are ⎧
⎨
⎩
−1
−2
−3
⎫
⎬
⎭ , (7-74)
and therefore the model is stable and nonoscillatory.
7-130 Schimel, Ding, Anderson and Grantham
There are many possibilities for a set of output equations to the inverted pendulum model. A simple
output vector for the linearized inverted pendulum model with PID-control is
y =
⎡
⎣
1 0 0
0 1 0
0 0 1
⎤
⎦
⎧
⎨
⎩
x0
x1
x2
⎫
⎬
⎭ . (7-75)
7.4 Defining Terms
asymptotically stable: An equilibrium point is asymptotically stable if in addition to being stable all
initial states starting sufficiently close to the equilibrium point approach the equilibrium point as t→∞
input vector: A one-dimensional column array consisting of either deterministic or uncertain inputs to
the state-space model.
measurement uncertainty vector: A one-dimensional column array whose terms model noise and
other sources of uncertainty in the model output.
output vector: A one-dimensional column array whose elements model the measurements on a the
physical system. Elements of output vector can consist of any combination of state, input and measurement
uncertainty vectors.
physical system: The real system, consisting of physical components, real inputs and response
stable: An equilibrium point is stable if every initial state which begins near the equilibrium point stays
near the equilibrium point indefinitely.
state equations: A set of first order differential equations which model the behavior of a physical
system.
state-space: A geometric space with dimension equal to the number of state variables. Any possible
state of a dynamic model can be represented as a point in state-space.
state variables: The smallest set of variables which can completely describe the state of a dynamic
system. It is possible to define more than one set of state variables for any particular model, however the
number of state variables is a unique quantity of the model.
state vector: A one-dimensional column array consisting of the state variables of a model.
stiff system: Systems of differential equations having two or more time scales.
unstable: If an equilibrium point does not meet the definition of stable then it is unstable.
7.5 Further Information
The following two texts present additional examples and in-depth information on state-space models:
B. Friedland, Control System Design an Introduction to State-Space Methods, New York: McGraw-Hill,
Inc, 1986.
W.J. Grantham and T.L. Vincent, Modern Control Systems Analysis and Design, New York: John Wiley
and Sons, Inc, 1993.
A complete presentation of state-space models for discrete-time systems appears in chapters five and six
of:
K. Ogata, Discrete-Time Control Systems, Englewood Cliffs, NJ: Prentice-Hall, Inc. 1987.