Parameter Estimation for Oscillatory Systems

PROCEEDINGS OF THE YEREVAN STATE UNIVERSITY Physical and Mathematical Sciences 2021, 55(2), p. 131–140 Mechanics PARAMETER ESTIMATION FOR OSCILLATORY SYSTEMS A. A. MATEVOSYAN 1∗ , A. G. MATEVOSYAN 2∗∗ 1 Department of Physics, University of Cambridge, UK 2 Chair of Mechanics, YSU, Armenia Simple harmonic motion was investigated of a rotational oscillating system. The effect of dumping and forcing on motion of the system was examined and measurements were taken. Resonance in a oscillating system was investigated and quality factor of the dumping system was measured at different damping forces using three different methods. Resonance curves were constructed at two different damping forces. A probabilistic model was built and system parameters were estimated from the resonance curves using Stan sampling platform. The quality factor of the oscillating system when the additional dumping was turned off was estimated to be Q = 71 ± 1 and natural frequency ω0 = 3.105 ± 0.008 s−1 . https://doi.org/10.46991/PYSU:A/2021.55.2.131 Keywords: statistical models, resonance, Bayesian methods, sampling, Stan. Introduction. The aim of this experiment is to estimate the parameters of the mechanical oscillating system by taking measurements. The experiment used a torsion pendulum, a variable frequency motor for driving force, and an "eddy brake" [1] for dumping. The pendulum mplitude was measures versus the number of complete oscillations, as well as the damping constant. Resonance curves were found at two different damping forces. The quality factor was measures in three ways: by the decay constant, by resonance amplitude and by the resonance curve. The Stan package [2] was used to analyze the resonance curve. This paper is organised as follows. The necessary theory is summarized in the Section 1. Section 2 describes apparatus used in the experiment, illustrates the method used and the results of the experiment. Section 3 is dedicated to the derivation of parameters from the resonance curve. Section 4 discusses the results, and the general conclusions are presented at the end of the main body of the article. 1. Theoretical Background. The Equation of Motion. The equation of motion for an undamped torsional oscillator has the form I θ̈ + τθ = 0 , (1) ∗ ∗∗ E-mail: [email protected] E-mail: [email protected] 132 A. A. MATEVOSYAN, A. G. MATEVOSYAN where I is the moment of inertial of the oscillator, θ is the angular displacement from the equilibrium, and τ is the torsion constant. The general solution to this equation is θ = θ0 cos(ω0t + φ ) , (2) s τ is the natural frequency of the oscillator and φ is the initial phase. I The system is damped if there is an opposing force proportional to θ̇ . Now the equation of motion takes the following form where ω0 = I θ̈ + bθ̇ + τθ = 0 , (3) where b is a constant. Finally, when a sinusoidal driving force with angular frequency ω is applied, we get I θ̈ + bθ̇ + τθ = A cos(ωt) . (4) In other notation, Eqs. 3 and 4 can be written in the following form θ̈ + 2γ θ̇ + ω02 θ = 0 , θ̈ + 2γ θ̇ + ω02 θ = f cos(ωt) , (5) (6) s τ b is the natural frequency of undamped oscillations, 2γ = is a I I A damping measure and f = is a measure of driving amplitude. I Solution of Equation of Motion. The general solution [3] to Eq. 6 is the sum of the solution of Eq. 5 (called transient response) and a particular solution to Eq. 6 (called steady state or forced response). The Eq. 5 has three different solutions depending on the magnitude of damping. These different cases are called underdamped, overdamped and critically damped. The last one is the special case: a critically damped system does not oscillate, but it returns to equilibrium faster than an overdamped system. In the underdamped case, the transient response has the form where ω0 = θ = θ0 e−γt cos(ω1t) , (7) q where ω1 = ω02 − γ 2 . γ is called the decay constant. Quality factor is a dimensionless parameter that describes how underdamped the oscillator is, how well it oscillates [3]. The quality factor Q is defined as ω0 Q= . (8) 2γ For a good oscillators Q  1. The forced response has the form θ (ω) = X(ω) cos(ωt − φ (ω)) , where ω is the angular frequency of the driving force, and f X(ω) = q , 2 2 2 2 (ω0 − ω ) + (2γω) (9) (10) PARAMETER ESTIMATION FOR OSCILLATORY SYSTEMS 133 π 1 γ = 0.2 γ = 0.3 γ = 0.5 φ (ω) X(ω) γ = 0.1 γ = 0.5 γ = 0.1 0 0 ω0 0 frequency ω frequency ω Figure 1: Dependence of amplitude (on the left) and phase difference (on the right) on the frequency of the driving force. tan φ (ω) = 2γω ω02 − ω 2 . (11) X(ω) is the amplitude and φ (ω) is the phase difference between the driving force and oscillatory response. The graphs of these functions are shown in Fig. 1. X(ω) gets its maximal value f Xmax ≡ X(ωmax ) = q 2γ ω02 − γ 2 (12) q at ωmax = ω02 − 2γ 2 . Note that when γ  ω0 we can write f f ω0 Xmax ≈ = · = X(0) · Q . 2γω0 ω02 2γ (13) 2. Methods and Results. Apparatus. The experiments were carried out on an oscillatory system consisting of: • a bronze disc, which undergoes rotational motion. The restoring couple provided by a coiled spring is linearly proportional to the angular displacement of the pendulum, which is measured by the scale around the disc. • a motor, with frequency control. The pendulum can be driven sinusoidally by an arm connected to the motor. The frequency of the motor is controlled by the voltage across it. • two coils of wire, which form an “eddy brake”, as a damping couple. Increasing the current in the eddy brake increases the damping of the system. A detailed diagram of the oscillating system is shown in Fig. 2. 134 A. A. MATEVOSYAN, A. G. MATEVOSYAN Figure 2: Diagram of the oscillatory system used in the experiments. Measurement of Period of the Pendulum. Four time measurements were made for ten complete periods of the pendulum. The results are shown in Tabs. 1 and 2. Time for 10 oscillations 20.35 s 20.22 s 20.10 s 20.27 s Table 1: Measurements of 10 oscillations. mean σ σm Period Frequency ω1 20.24 s 0.10 s 0.05 s 2.024 ± 0.005 s 3.105 ± 0.008 s−1 Table 2: Period and angular frequency of the pendulum. "Eddy brake" is off. Even though "eddy brake" is off, i.e. current through it is zero, the damping force is still present, such as friction in the axle. But it is very small (γ  ω0 ) so we can assume ω0 ≈ ω1 (see Eq. 7). Measurement of Amplitude. According to the Eq. 7, n-th successive oscillation has an amplitude an = a0 e−nγT , (14) 2π is the period and a0 is the initial amplitude. Taking the logarithm on ω1 both sides, we get a linear dependence on n: ln(an ) = ln(a0 ) − nγT . (15) The experiment was carried out under two conditions: with no current through the eddy brake and with a braking current Ib = 0.6 A, without driving force. In both cases the pendulum was released from rest at 10 on the scale. For more accurate measurements, the process was recorded by a video camera, then the amplitude was determined by playback in slow motion. The results are shown in Fig. 3. were T = 135 PARAMETER ESTIMATION FOR OSCILLATORY SYSTEMS ln(amplitude an / arbitrary units) 2.5 braking current = 0.0 A braking current = 0.6 A 2 −4.39 · 10−2 · x + 2.31 1.5 1 0.5 0 −0.5 −0.805 · x + 2.38 −1 −1.5 0 2 4 8 10 6 Number of oscillations, n 12 14 16 Figure 3: Dependence of ln(amplitude) on the oscillation number n. The relationship is linear as we would expect from the Eq. 15. First Estimation of Quality Factor. From Eq. 15 a graph of ln(an ) against n should be a straight line with 2π slope of the line = −γ T = −γ . (16) ω1 π For low damping ω1 ≈ ω0 , so the slope is approximately − . By statistical analysis Q of the data, we got slope [Ib = 0.0A] slope [Ib = 0.6A] and by assumption slope ≈ − −0.044 ± 0.001 −0.80 ± 0.03 π , Q Q [Ib = 0.0A] Q [Ib = 0.6A] 71 ± 1 3.9 ± 0.1 Forced Oscillations. Now a sinusoidal driving force is applied to the pendulum. The general solution to the equation of motion 6 is a superposition of the transient response (7) and the forced response (4). As there is a damping force, the transient response vanishes over time, and only the forced response remains. 136 A. A. MATEVOSYAN, A. G. MATEVOSYAN The experiment was carried out with a damping current of 0.3A and 0.6A. The frequency of the driving force was changed by the supply voltage and was determined by taking time of 10 complete rotations of the drive wheel (see the white mark on the wheel of the motor, Fig. 2). Resonance Curve. The measurements were made by varying the voltage of the driving motor, thus changing the frequency of the driving force. After some time, only the steady state solution (forced response) was seen. The amplitude was measured. At first measurements were taken over the whole frequency range. Then more points were taken close to the resonance frequency to define the shape of the resonance curve more accurately. The results are shown in Fig. 4 and Tab. 3. Ib = 0.3A Ib = 0.6A Xmax ωmax 11.5 ± 0.1 3.3 ± 0.1 3.11 ± 0.02 s−1 3.05 ± 0.05 s−1 Table 3: Resonance frequency and amplitude at resonance. Amplitude / arbitrary units 12 braking current = 0.3 A braking current = 0.6 A 10 8 6 4 2 0 1.5 2 2.5 3 3.5 4 Angular frequency / 4.5 5 5.5 s−1 Figure 4: Amplitude of forced oscillations as a function of angular frequency. Second Estimation of Quality Factor. The quality factor can be estimated using Eq. 13. But first it needs to estimate X(0). To measure the amplitude at ω = 0, the drive wheel was slowly rotated by hand and the maximum displacement of the pendulum was read on both sides of zero. PARAMETER ESTIMATION FOR OSCILLATORY SYSTEMS 137 The amplitude at ω = 0 was estimated as X(0) = 0.7 ± 0.1. Thus Q2nd [Ib = 0.3A] Q2nd [Ib = 0.6A] 16 ± 2 4.7 ± 0.7 3. Statistical Analysis of Resonance Curve. Model. From our experiment, we have (ωi , Xi ) pairs for i = 1, ..., n, where n is the number of measurements taken. It can be assumed that our measurements come from the distribution Xi ∼ P( · | ωi , ω0 , X0 , γ) ≡ N(X(ωi ; ω0 , X0 , γ), σ 2 ) for i = 1, ..., n , (17) 2 X0 ω0 X(ω; ω0 , X0 , γ) = q , (18) (ω02 − ω 2 )2 + (2γω)2 where N(µ, σ 2 ) is the Normal distribution and X(ω; ω0 , X0 , γ) (reparameterized version of Eq. 10) depends only on the parameters ω0 , γ and X0 . The latter is the amplitude at zero frequency. Since our measurements are accurate to 0.1, we take σ = 0.1. Our task is to estimate the parameters ω0 , γ and X0 from the resonance curve (ωi , Xi )1...n . More formally, find the posterior distribution ω0 , γ, X0 ∼ P( · | ω1 ...ωn , X1 ...Xn ) . (19) Unfortunately, the posterior distribution 19 is intractable. Therefore, we do not have its analytical form for calculating the moments of the parameters. However, without having explicit form of the posterior distribution, we can sample from that distribution. Using Stan Sampling Platform. Stan is a state-of-the-art platform for statistical modeling and high-performance statistical computation [2]. The Stan language is used to specify a (Bayesian) statistical model with an imperative program calculating the log probability density function. More about Stan at mc-stan.org. From the Bayes theorem P(X1 ...Xn | ω0 , γ, X0 , ω1 ...ωn ) · P(ω0 , γ, X0 ) P(ω0 , γ, X0 | X1 ...Xn , ω1 ...ωn ) = margnal n (20) ∏ P(Xi | ω0 , γ, X0 , ωi ) · P(ω0 , γ, X0 ) = i=1 , margnal where the first term on the right is defined by our model (Eq. 17) and the second term is the prior distribution – the prior confidence in the parameters. We used the Gamma distribution as a prior. Therefore, all terms in the numerator of Eq. 20 are defined, so we can do sampling. Stan uses MC-MC algorithm to draw samples from the distribution 20. Sampling Results. The mean and standard deviation of the samples were calculated and the values obtained are shown in Tab. 4. The shape of the curve corresponding to inferred parameters is shown in Fig. 5 to compare with experimental data. 138 A. A. MATEVOSYAN, A. G. MATEVOSYAN ω0 / s−1 γ / s−1 X0 ≡ X(0) Q3rd / s−1 Ib = 0.3 A Ib = 0.6 A 3.105 84 ± 0.000 12 0.090 59 ± 0.000 16 0.6675 ± 0.0010 17.14 ± 0.03 3.1025 ± 0.0013 0.3363 ± 0.0018 0.708 ± 0.003 4.61 ± 0.02 Table 4: Resonance frequency and amplitude at resonance. Amplitude / arbitrary units 12 Ib = 0.3 A Ib = 0.6 A 3 10 8 2 6 4 1 2 0 0 2 4 0 6 Angular frequency / s−1 2 4 Figure 5: Resonance curve shape using parameters found by inference. 4. Discussion. Consistency with Theory. From Fig. 3 we see that the points lie quite well on the line predicted by Eq. 15. Also, when we increase the damping force, we get a lower quality factor, as expected from Eq. 8. In the case of forced oscillations, our resonance curves shown in Fig. 4 have the same shape as in Fig. 1 derived from theory. In the case of lower damping, we have a higher resonance amplitude, which is also consistent with Eq. 10 (also visualized in Fig. 1). However, at lower damping, we should have a higher resonance frequency, which follows from Eq. 12, but from the results of Tab. 3 this is not obvious. Three Estimations of Quality Factor. We have estimated the quality factor in three different ways at same damping: • Q = 3.9 ± 0.1 – transient response method, • Q = 4.7 ± 0.7 – forced response method, • Q = 4.61 ± 0.02 – resonance curve method. PARAMETER ESTIMATION FOR OSCILLATORY SYSTEMS 139 The second and third values are consistent with each other. Note that the large error in the second estimate comes from the measurement of X(0) with high uncertainty. However, the first and third values are very different, there could have been systematic error in picking the braking current. Improvements to the Experiment. At high damping, e.g. Ib = 0.6A, approximaq tion ω1 = ω02 − γ 2 ≈ ω0 is not accurate anymore. On the other hand, we can easily measure the frequency ω1 . This way we will get the best Q-factor estimate. We can construct a better resonance curve (Fig. 4) simply by taking more measurements. We can also calculate the quality factor from the bandwidth of the resonance curve [3] and compare it with previous methods. Conclusions. We have shown that our measurements are consistent with theory, e.g. exponentially decaying amplitude, shape of the resonance curve. Different methods have been provided for analyzing measurements and deriving system parameters such as natural frequency, decaying factor, etc. The probabilistic method for inferring system parameters described in Sec. 4 gave estimates with high precision. Thus, we can conclude that 10 points of the resonance curve contain all information about the system. Received 15.01.2021 Reviewed 29.03.2021 Accepted 12.04.2021 REFERENCES 1. Cadwell L.H. Magnetic Damping: Analysis of an Eddy Current Brake Using an Airtrack. Am. J. Phys. 64 (1996), 917–923. http://doi.org/10.1119/1.18122 2. Carpenter B., Gelman A., et al. Stan: A probabilistic programming language. J. Stat. Softw. 76(1) (2017), 1–32. http://doi.org/10.18637/jss.v076.i01 3. Feynman R.P., Leighton R.B., Sands M.L. The Feynman lectures on physics. Reading, Mass: Addison-Wesley Pub. Co. (1963) 140 A. A. MATEVOSYAN, A. G. MATEVOSYAN A. A. MAEVOSYAN, A. G. MAEVOSYAN TATANOAKAN HAMAKARGERI PARAMETRERI GNAHATOWM Owsowmnasirvel  pttakan tatanoakan hamakargi parz nerdanak arowm: Hetazotvel en maro harkadrakan owerii azdecowyown hamakargi vra, inpes na katarvel en aowmner: Ereq tarber meodnerov hamakarg ditarkvel  ezonansayin eimowm maro hamakargi tatanowmneri orakakan gorakic avel  tarber maro oweri depqowm: Erkow tarber maro oweri depqowm kaowcvel en ezonansayin korer: gtagorelov Stan nmowaman mijavayr` ezonansayin koreri himan vra kaowcvel  havanakanayin model gnahatvel en hamakargi parametrer: Lracowci maro oweri bacakayowyan depqowm hamakargi tatanowmneri orakakan gorakic gnahatvel  orpes Q = 71 ± 1, isk seakan haaxowyown orpes ω0 = 3.105 ± 0.008 s−1 : А. А. МАТЕВОСЯН, А. Г. МАТЕВОСЯН ОЦЕНКА ПАРАМЕТРОВ КОЛЕБАТЕЛЬНЫХ СИСТЕМ Было исследовано простое гармоническое движение вращательноколебательной системе. Также исследовались затухание и вынужденные колебания системы, и были произведены измерения. Тремя разными методами исследовался резонанс в колебательной системе и измерялся коэффициент качества (добротность) затухающей системы при различных параметрах затухания. При двух различных демпфирующих силах были построены резонансные кривые. Построена вероятностная модель и оценены параметры системы по резонансным кривым с использованием платформы для отбора Stan. При отсутствии дополнительной силы сопротивления, коэффициент качества (добротность) колебательной системы оценивалась в Q = 71 ± 1, а собственная частота как ω0 = 3.105 ± 0.008 s−1 .