EM 406 Vibrations Laboratory 2 - System Identification of a Mass-Spring-Damper System We will investigate the effects of varying the parameters of a physical spring mass damper system, and see how its behavior is different from the lumped parameter model. Objectives: The objectives of this lab are to: --Estimate mass, stiffness, damping, and static gain of a system --Compare experimental data of system impulse response to a theoretical simulation. --Explain some of the limitations of lumped parameter models. Analysis: For the data from Test 2 determine ζ and ωn using log decrement 1) Measure the period, T, of the experimental response and determine the damped natural frequency (ωd) knowing that the period is given by Eq. 1. dTωπ2= (1) 2) Using the log decrement determine the damping ratio n0xxlnn1=δ (2) Knowing δ you can determine the damping ratio from Ea. 3. 224πδδζ+= (3) 3) Determine the natural frequency knowing 21ζωω−=nd (4) After completing these calculations for Test2a, Test2b, and Test2c fill out Table 1 in the Worksheet for Lab #2. Determine the equivalent mass of the cart 1) Use the results in Table 1 estimate the actual cart and damper mass. You can use a procedure as follows: a) Generate three equations by applying the known values of damped natural frequency and damping ratio to Eq. 121di iikmωζ=− (1) b) These equations can be manipulated to the following form which is linear in the unknowns (k and meq,cart) as shown in Eqns. 2-5. ()()012121=−−ζωkmdcart,eq (2) ()()0112222=−−+ζωkmdcart,eq (3) ()()01512323=−−+ζωk.mdcart,eq (4) Equations 2-4 can be written in matrix form as shown in Eq. 5. ⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧−−=⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡−−−22232222212323215110111ddcart,eqddd.kmωωζωζωζω (5) The estimates of mass, and stiffness can be obtained using Eq. 5 above. The way to use Matlab to solve an overdetermined set of equation is shown at the end of this document. The damping constant c can be found by matching coefficients between the expected system model as shown in Eq. 6 and standard form of the model as shown in Eq. 7. 1()mcxxx ftkk k++=&& & (6) 212()nnxxxKftζωω++=&& & (7) This comparison results in the following relationship between m, k, and ωn: mkn=ω (8) Solving for the damping constant gives iniikcωζ2= (9)We will now use several techniques to do a full system identification for Test2b. We will only perform a complete system identification for the run “Test 2b”. Assume the equivalent mass of the cart is the value determine above. Method 1 – Log decrement method Using the meq,cart determined earlier and shown in Table 2 determine the spring stiffness, k, and viscous damping coefficient, c, for Test2b using the frequency and damping ratio determined using the log decrement. Method 2 – Using a performance index in Excel Generally, analytical solutions and experimental measurements differ. The theoretical curve can be made to more closely (though usually not exactly) predict the experimental measurements by adjusting the model parameters (in this case, the mass, stiffness and damping). A measure of experimental/theoretical closeness is the performance index J given by ()21testmodel1∑=−=niiixxnJ (10) where: modelix = cart position predicted by the model at the ith point in time dataix = cart position experimentally determined at the ith point in time n = number of data points used in comparing experiment to theory A second order system forced with an impulse is shown in Eq. 11 ()tkFxxxnnδωζω=++&&&22 (11) With zero initial conditions Eq. 11 has a solution given by Eq. 12 (see page 314 in the Text) () ()tmFetxddeqtnωωζωsin−= (12) Your task is to find the value of F, ζ and ωn that minimizes J. The easiest way to do this is in Excel, although Matlab can be also be used. Let’s first look at how to do this in Excel. Finding parameters using Excel Format a spreadsheet as shown in Figure 4. Columns A and B contain your experimental data. Be sure your experimental data is in meters and not centimeters. Column C contains the theoretical prediction (your analytical solution of the ODE) which is affected by your choice of parameters. The formula entered in Column C should reference the cells containing omegan, zeta and omegad which in turn reference F, k, m and c. Column D is the square of the differencebetween theoretical and experimental. Cell E6 contains the sum of Column D. Cells J2, J3, and J4 contains your initial guesses of c, k and F respectively (use the values you found using log decrement and a least squares solution as a first guess). Figure 4: Spreadsheet example for tuning the model. Use Solver to minimize the sum of the square error by varying the value of your parameters. If you do not have Solver under your “tools” menu go to “Add-ins” and add it. To use Solver set the “target cell” as the cell in which the cost function is computed and set the “by changing cell” to the cell in which the parameters k, c and mag are located. To start the optimization, click on the Solve key. A sample Solver window is shown in Figure 5.Figure 5: Spreadsheet example for tuning the model. Method 3 – Using a performance index in Matlab (viscous damping) Matlab can also be used to minimize a performance index. One advantage of using Matlab is that it allows us to use Simulink to generate our theoretical model so we can use a more complex model that may not have a closed form solution. Since you know the input was actually a pulse of width 0.05 s rather than an impulse you can use this in your Simulink model. You will need three files: File 1: You need a Simulink model as shown in Figure 6. This model has been developed so that it can have viscous damping, Coulomb damping, or a combination of the two. a. Use a fixed times step of magnitude 0.01 or smaller. You can do this in the “Simulation” menu. The parameters, “mag”, “fd”, and “c” will be assigned in an m-file discussed next. Insert these into the step and Coulomb/viscous block. $F$4$G$2File 2: A Matlab m-file to set initial parameters and call “fminsearch” as shown below. % m-file to define initial parameters and call fminsearch Clear load data2b % this loads the data file with variables time and x1 m = 1.7;