2.004 Problem Set #9 Solutions

MASSACHUSETTS INSTITUTE OF TECHNOLOGY Department of Mechanical Engineering 2.004 Dynamics and Control II Problem Set #9 Solution 1. The 2.004 Tower system. The system parameters are m1 = 5.11 kg, b1 = 0.767 N · sec/m, k1 = 2024 N/m; m2 = 0.87 kg, b1 = 8.9 N · sec/m, k1 = 185 N/m. The system model schematic is shown again below for your convenience. Eigenvalues and Eigenvectors in state space representation The eigenvalues of A matrix are poles. Hence,  the imaginary parts of the poles are the damped oscillation frequencies (ωd = ωn 1 − ζ 2 ), and the real parts are ζωn . The eigenvectors of A matrix represent the modes of the system. Each element of a certain eigenvector indicates the way that state variables change when the system is excited with the corresponding damped natural frequency. For example, if an eigenvector is [−2 1]T , then the two state variables move in opposite way and the ratio of their magnitudes is 2. The real and imaginary parts of the eigenvalues represent decay rate (or settling time) and natural frequency. For the complex eigenvector, you may interpret that the real and imaginary part correspond to cos and sin components of the response. Since cos θ = sin(θ + π/2) and exp {jπ/2} = j, the real part represents the magnitude of cos(ωt) motion, and the imaginary part represents the magnitude of sin(ωt) motion. To give you an easy example, let’s consider the 2.004 damping. Then the system matrix is ⎡ 0 1 0 0 ⎢ −432.2896 0 36.2035 0 A=⎢ ⎣ 0 0 0 1 212.6437 0 −212.6437 0 Tower system without ⎤ ⎥ ⎥, ⎦ the eigenvectors are ⎡ ⎤ ⎡ ⎤ ⎡ j0.0354 −j0.0354 −j0.0106 ⎢ −0.7614 ⎥ ⎢ −0.7614 ⎥ ⎢ 0.1427 ⎥ ⎢ ⎥ ⎢ v1 = ⎢ ⎣ −j0.03 ⎦ , v2 = ⎣ j0.03 ⎦ , v3 = ⎣ −j0.0732 0.6466 0.6466 0.9870 1 ⎤ ⎡ ⎤ j0.0106 ⎥ ⎢ ⎥ ⎥ , v4 = ⎢ 0.1427 ⎥ . ⎦ ⎣ j0.0732 ⎦ 0.9870 and the eigenvalues are {±j21.5183, ±j13.4870}. Note that the eigenvalues are purely imaginary because there is no damping. Let’s look at v3 corresponding to the lower natural frequency 13.4870 (rad/s). The first (−j0.0106) and third (−j0.0732) elements are the displacement of the tower and slider, respectively. So they move along the same direction in sin fashion (because they are imaginary) without phase delay. The second (0.1427) and fourth (0.9870) elements are velocity of the tower and the slider, respectively. They change with the same sign in cos fashion (because they are real) without phase delay. Since velocity is a time derivative of displacement, it totally makes sense that the first and third elements are purely imaginary while the second and fourth ones d are pure real, because (sin ωt) ∼ cos(ωt). dt Therefore at the lower natural frequency, the tower and slider move in–phase with each other and they oscillate in sin(ωt) fashion. Their velocities oscillate also in–phase with each other but in cos(ωt) fashion. Therefore the velocity has phase delay of π/2 with respect to the displacement. (You would get the same conclusion from v4 , which corresponds to −13.4870 (rad/s).) If you look at v1 and v2 , then you would notice that the slider and tower dis­ placements and velocities move in opposite direction. But you still have the same phase delay between the displacement and the velocity (one goes like sin, the other like cos.) We are now ready to attack the given problem, i.e. the tower and slider with damping. a) Eigenvectors and eigenvalues of matrix A. Answer: From A matrix given by ⎛ 0 1 0 0 ⎜ −(k1 + k2 )/m1 −(b1 + b2 )/m1 k2 /m1 b /m 2 1 A = ⎜ ⎝ 0 0 0 1 b2 /m2 −k2 /m2 −b2 /m2 k2 /m2 ⎛ ⎞ 0 1 0 0 ⎜ −432.2896 −1.8918 36.2035 1.7417 ⎟ ⎟, = ⎜ ⎝ ⎠ 0 0 0 1 212.6437 10.2299 −212.6437 −10.2299 the eigenvectors and eigenvalues are computed by Matlab. 2 ⎞ ⎟ ⎟ ⎠ The eigenvectors are ⎡ 0.0319 + j0.0037 ⎢ −0.1410 + j0.6315 v1 = ⎢ ⎣ −0.0039 − j0.0376 0.7608 ⎡ −0.0067 + j0.0074 ⎢ −0.0757 − j0.1214 v3 = ⎢ ⎣ 0.0188 + j0.0659 −0.9873 ⎤ ⎥ ⎥ , v2 ⎦ ⎤ ⎥ ⎥ ⎦ , v4 ⎡ ⎤ 0.0319 − j0.0037 ⎢ −0.1410 − j0.6315 ⎥ ⎢ ⎥, =⎣ −0.0039 + j0.0376 ⎦ 0.7608 ⎡ ⎤ −0.0067 − j0.0074 ⎢ −0.0757 + j0.1214 ⎥ ⎥ =⎢ ⎣ 0.0188 − j0.0659 ⎦ , −0.9873 and the eigenvalues are {−2.1014 ± j20.0261, −3.9595 ± j13.8582}. b) Justify the following statement: “The 2.004 Tower has two modes of oscilla­ tion, one slow and one fast.” Compute the damped and natural frequencies of oscillation of the two modes. Answer: We have 4 poles at s = −2.1014 ± j20.0261, −3.9595 ± j13.8582. Thus we have two damped natural frequencies (13.8582, 20.0261 (rad/s)) and two modes associated with these natural frequencies. Keep in mind that the imaginary part of the pole is ωd = ωn 1 − ζ 2 and the real part is ζωn . (Or you can use the geometrical relation that cos θ = ζ, where θ is the angle to the pole from the negative real axis.) • Mode 1: ωd = 13.8582 (rad/s), ζ = 0.2747, and ωn = 14.4128 (rad/s) • Mode 2: ωd = 20.0261 (rad/s), ζ = 0.1044, and ωn = 20.1360 (rad/s). c) Justify the following statement: “In the slow mode the tower and slider mass oscillate in phase, while in the fast mode the tower and slider mass oscillate out of phase.” Answer: The eigenvectors v3 and v4 of the slow mode are given by ⎡ ⎤ ⎡ ⎤ −0.0067 + j0.0074 −0.0067 − j0.0074 ⎢ −0.0757 − j0.1214 ⎥ ⎢ −0.0757 + j0.1214 ⎥ ⎥ ⎢ ⎥ v3 = ⎢ ⎣ −0.0188 + j0.0659 ⎦ , v4 = ⎣ −0.0188 − j0.0659 ⎦ . −0.9873 −0.9873 Let’s look at the position of the tower(−0.0067 ± j0.0074) and the slider (−0.0188 ± j0.0659). The real parts have the same sign and the imaginary part does too. Also note that the imaginary part dominates because it is bigger then the real part. Hence, they move a kind of in–phase. (Not exactly in–phase because the real part is not zero.) 3 At the fast mode, the eigenvectors are, ⎡ ⎤ ⎡ −0.0319 + j0.0037 −0.0319 − j0.0037 ⎢ −0.1410 + j0.6315 ⎥ ⎢ −0.1410 − j0.6315 ⎥ ⎢ v1 = ⎢ ⎣ −0.0039 − j0.0376 ⎦ , v2 = ⎣ −0.0039 + j0.0376 0.7608 0.7608 ⎤ ⎥ ⎥. ⎦ Let’s look at the first (−0.0319 + j0.0037) and third (−0.0039 − j0.0376) elements in the same context. The first element is dominated by the real part, and the third one is dominated by the imaginary part. Thus we expect that the displacement of the tower and the slider have a phase delay about π/2, which means that they are out of phase (or not in–phase). (You would get the same conclusion from v2 , whose damped frequency is negative of the damped frequency of v1 .) Therefore, the tower and slider oscillate in–phase in the slow mode, while out of phase in the fast mode. d) Which mode has been excited by the impulse input? Compare the damped frequency of oscillation of the modethat you think has been excited with the frequency of oscillation that you measure from the impulse response simulation. Answer: −3 10 Impulse Response x 10 tower displacement 8 6 Amplitude 4 2 0 −2 −4 System: tower1 Time (sec): 1.19 Amplitude: −0.000891 System: tower1 Time (sec): 0.874 Amplitude: −0.00163 −6 −8 0 0.5 1 1.5 Time (sec) 2 2.5 3 2π = 19.8835 (rad/s). It (1.19 − 0.874) is very close to the higher mode frequency. We conclude that the second mode is predominantly excited by the impulse. The measured damped frequency ωd = e) Is the phase relationship between the tower1 and tower2 responses consis­ tent with the oscillation mode that the system is in? Answer: 4 Impulse Response 0.01 tower displacement slider displacement 0.008 0.006 Amplitude 0.004 0.002 0 −0.002 System: tower2 Time (sec): 0.967 Amplitude: −0.00132 −0.004 System: tower2 Time (sec): 0.655 Amplitude: −0.00372 −0.006 −0.008 −0.01 0 0.5 1 1.5 Time (sec) 2 2.5 3 2π = 20.1384 (rad/s). (0.967 − 0.655) It is very close to the higher mode frequency and the second (faster) mode is excited. The displacement of the slider is delayed with respect to the displacement of the tower. Note that the position of the tower (blue) and the slider (green) are out of phase by about π/2 as we expected in (c). The measured damped frequency: ωd = 2. Nise chapter 12, problem 4 (page 777). Answer: G(s) = 20 20 = 3 . (s + 1)(s + 3)(s + 7) s + 11s2 + 31s + 21 Please refer to the notes (Lecture 18) or chapter 3.5 in the Nise textbook about converting a transfer function to state space representation. ⎡ ⎤ ⎡ ⎤ 0 1 0 0 0 1 ⎦q + ⎣ 0 ⎦u q̇ = ⎣ 0 −21 −31 −11 20   y= 1 0 0 q With a state–variable feedback loop, ⎡ ⎤ 0 1 0 ⎦ 0 0 1 A − BK = ⎣ −(k1 + 21) −(k2 + 31) −(k3 + 11) det(sI − (A − BK)) = s3 + (k3 + 11)s2 + (k2 + 31)s + k3 + 21 = 0. The desired %OS is 15% and a settling time is 0.75 (s). From these desired values, ζ= − ln(%OS/100) π2 2 + ln (%OS/100) 5 = 0.5169, and 4 = 10.3173. (0.75)ζ The desired characteristic equation is ωn = (s + p)(s2 + 2ζωn s + ωn2 ) = (s + p)(s2 + 10.666s + 106.4467) The third pole should be 10 times as far from the imaginary axis as the dominant pole pair s = −5.333 ± j8.8321. Thus we choose p = 53.3. The characteristic equation should be identical to the desired one, s3 + (10.666 + p)s2 + (106.4467 + 10.666p)s + 106.4467p = s3 + (k3 + 11)s2 + (k2 + 31)s + k1 + 21. Thus k3 + 11 = 10.666 + p, k2 + 31 = 106.4467 + 10.666p, k1 + 21 = 106.4467p. Using p = 53.3, we find k1 = 5625.6, k2 = 643.9, and k3 = 53.0. The beauty of the state space representation and state feedback is that you can observe and control individual state variables (e.g. position, velocity, acceleration, etc.) of systems. In complex systems, the state space representation may give you easier understanding of the system dynamics and better. 3. Nise chapter 9, problem 1 (page 674). Answer: a. G(s) = 1 s(s + 2)(s + 4)   1 1   = =   3 2 s(s + 2)(s + 4) s=jω s + 6s + 8s s=jω 1 1 −6ω 2 − j(8ω − ω 3 ) = = . −jω 3 − 6ω 2 + 8jω −6ω 2 + j(−ω 3 + 8ω) 36ω 4 + (8ω − ω 3 )2 G(jω) = Hence, 1 M (ω) =  , 4 36ω + (8ω − ω 3 )2  and −1 φ(ω) = tan b. G(s) = ω2 − 8 −6ω  . s+5 s(s + 2)(s + 4) −6ω 2 − j(8ω − ω 3 ) 5 + jω (5 + jω) = G(jω) = −6ω 2 + j(−ω 3 + 8ω) 36ω 4 + (8ω − ω 3 )2 {−30ω 2 + ω 2 (8 − ω 2 )} + j {−6ω 3 − 5(8ω − ω 3 )} = 36ω 4 + (8ω − ω 3 )2 −ω 2 (ω 2 + 22) + jω(−ω 2 − 40) = . 36ω 4 + (8ω − ω 3 )2 6 Hence,  M (ω) = {−ω 2 (ω + 22)}2 + {ω(−ω 2 − 40)}2 36ω 4 + (8ω − ω 3 )2  and −1 φ(ω) = tan c. G(s) = −ω 2 − 40 −ω(ω + 22) ,  . (s + 3)(s + 5) s2 + 8s + 15 = 3 s(s + 2)(s + 4) s + 6s2 + 8s  (15 − ω 2 ) + j8ω −6ω 2 − j(8ω − ω 3 )  2 G(jω) = (15 − ω ) + j8ω = −6ω 2 + j(−ω 3 + 8ω) 36ω 4 + (8ω − ω 3 )2 {−6ω 2 (15 − ω 2 ) + 8ω(8ω − ω 3 )} + j {−48ω 3 − (15 − ω 2 )(8ω − ω 3 )} = 36ω 4 + (8ω − ω 3 )2 ω 2 (−2ω 2 − 26) + jω(−ω 4 − 25ω 2 − 120) = 36ω 4 + (8ω − ω 3 )2 Hence,  M (ω) = {ω 2 (−2ω 2 − 26)}2 + {ω(−ω 4 − 25ω 2 − 120)}2 36ω 4 + (8ω − ω 3 )2  and −1 φ = tan −ω 4 − 25ω 2 − 120 ω(−2ω 2 − 26) ,  . 4. Nise chapter 9, problem 2 (page 674). Use Matlab to generate the plots. Answer: From the analytical expression in question 3, the magnitude and phase of the system are plotted here without using bode command in Matlab. Note that the magnitude is not plotted in (dB). 1 s(s + 2)(s + 4) 0 Magnitude 10 −5 10 −1 10 0 1 10 10 2 10 Frequency (rad/s) −100 Phase(degree) a. G(s) = −150 −200 −250 −1 10 0 1 10 10 Frequency (rad/s) 7 2 10 b. G(s) = s+5 s(s + 2)(s + 4) 2 Magnitude 10 0 10 −2 10 −4 10 −1 10 0 1 10 10 2 10 Phase(degree) Frequency (rad/s) −100 −120 −140 −160 −180 −1 10 0 1 10 10 2 10 Frequency (rad/s) c. G(s) = s+5 s(s + 2)(s + 4) 2 Magnitude 10 0 10 −2 10 −1 10 0 1 10 10 2 10 Frequency (rad/s) Phase(degree) −90 −95 −100 −105 −110 −1 10 0 1 10 10 2 10 Frequency (rad/s) See the following results generated by built–in bode command of Matlab. Note that the magnitude are plotted in (dB). a. G(s) = 1 s(s + 2)(s + 4) 8 Bode Diagram 50 Magnitude (dB) 0 −50 −100 −150 −90 Phase (deg) −135 −180 −225 −270 −1 10 0 1 10 10 2 10 Frequency (rad/sec) b. G(s) = s+5 s(s + 2)(s + 4) Bode Diagram 20 Magnitude (dB) 0 −20 −40 −60 Phase (deg) −80 −90 −135 −180 −1 10 0 1 10 10 2 10 Frequency (rad/sec) s+5 s(s + 2)(s + 4) Bode Diagram 40 Magnitude (dB) 20 0 −20 −40 −90 −95 Phase (deg) c. G(s) = −100 −105 −110 −1 10 0 1 10 10 Frequency (rad/sec) 9 2 10 5. Nise chapter 9, problem 4 (page 674). a. G(s) = 1 s(s + 2)(s + 4) Answer: (Note that this is a color plot) b. G(s) = s+5 s(s + 2)(s + 4) Answer: (Note that this is a color plot) c. G(s) = s+5 s(s + 2)(s + 4) Answer: (Note that this is a color plot) 10 6. For the following transfer function G(s) = (s + 10) (s + 1)(s2 + 40s + 104 ) sketch the Bode asymptotic magnitude and phase plots, using appropriate cor­ rections, and compare with the exact result using Matlab. Answer: s  + 1 10 (s + 10)  102  = G(s) = s 40s (s + 1)(s2 + 40s + 104 ) 104 (s + 1) + +1 104 104 The system has one zero at s = −10, one real pole at s = −1, and one pole pair whose natural frequency ωn = 100 (rad/s) and damping ratio ζ = 40/2/ωn = 0.2. Hence, the system has three break frequencies: 1, 10, 100 (rad/s). The magnitude at the low frequency asymptote is 20 log10 10−3 = −60 (dB). The correction −20 log 2ζ = 7.96 (dB). The Bode asymptotic plot: 11 The exact Bode plot by Matlab: Bode Diagram −60 Magnitude (dB) −70 −80 −90 −100 −110 −120 0 Phase (deg) −45 −90 −135 −180 −2 10 −1 10 0 1 10 10 Frequency (rad/sec) 12 2 10 3 10