First Order System: Identifying System Parameters from Scattered Response

Problem Background
Determining System Parameters

The idea behind these tutorials is that you can view them in one window while running Simulink in another window.  Do not confuse the windows, icons, and menus in the tutorials for your actual Simulink windows.  Most images in these tutorials are not live - they simply display what you should see in your own Simulink windows.  All Simulink operations should be done in your Simulink windows.

Problem Background

Recall that we are investigating the characteristics of a propeller-shaft system to be used in a boat.  The damping moment on the system was assumed to be proportional to its angular velocity, which gave us the following free body diagram of the moments:

In this diagram:

Looking at the total moment present on the system resulted in the following system equation:

In our previous analysis (First Order System: System Identification), we applied a step input of T=800 N-m to the propeller-shaft system (by using the Simulink block we provided to approximate the actual system), and we were able to determine estimates for b and I by finding the steady-state value and time constant of the response.  Now, let's take a look at what happens if we wanted to estimate b and I using the much smaller input torque of T=100 N-m (for instance, we would do this if we were using a more complicated system model and wanted to linearize it about this operating point).  If we set up a Simulink model with a step input of 100 N-m into the Propeller-Shaft system block and monitor the output of this block with a scope, we see a response that looks like:

Determining System Parameters

Note that in the plot above, the response's fluctuations are a substantial proportion of its steady-state value, and this makes it hard to make an accurate estimation of the steady-state angular velocity and the time constant.  So, instead of using this plot to determine our system parameters, let's try an alternate method.  Notice that by the time the simulation has run for 3 seconds, the system appears to have reached its steady-state (the fluctuations in the response we see there are caused by disturbance inputs and noise in the tachometer signal).  Thus, if we average the angular velocity from 3 seconds forward, we can get a pretty good estimate of the steady-state response.  To do this, we will use the "mean" command in Matlab.  Modify your Simulink model of the system so that the angular velocity and the system time are output to the MATLAB workspace.  This is accomplished by using "To Workspace" blocks from the Sinks subfolder.  For the angular velocity, you can branch off the signal travelling to the Scope and bring it to a To Workspace block.   To have the system time get to its To Workspace block, you will need to use a clock (as below).  Open the Block Parameters windows for each of the To Workspace blocks, and change their Variable Names to "time" and "RPM", as appropriate.   Also change their Save Format to "Matrix"  Your model should now look like:

When we run the simulation, the shaft RPM will now be saved in a single-column MATLAB matrix named "RPM", and the system time will be stored in a matrix called "time".  Verify that the simulation Stop Time is set to 6.0 seconds by opening the Simulation Parameters window (click on Simulation, then Parameters).  Also make sure that the ode5 (Dormand-Prince) fixed-step solver is selected, and that the step size is 0.05 seconds.  Run the simulation, then enter the following in the MATLAB Command Window:

This displays all the measured values of the shaft RPM on the screen as a single column matrix.  Now, type the command:

MATLAB responds by telling you the dimensions of the matrix "RPM". The first value is its number of rows, and the second value indicates the number of columns. Thus, the matrix "RPM" is 121 rows by 1 column (which makes sense, since sampling every 0.05 seconds on the interval from 0 to 6 seconds gives 121 measurements).   We would like to average the values taken from t = 3 seconds and on in the simulation.  So, we enter the following command in the MATLAB window:

MATLAB responds by telling us the steady-state angular velocity is about 314 RPM (32.8 rad/sec).

Next, we must determine a way to find the time constant of the system.  To do this, let's consider the solution of the first-order system equation we are approximating the system with using a step input (T) applied when the system is initially at rest.  Using Laplace transforms, we convert the system equation to the Laplace domain and solve for the angular velocity as a function of time:

If we take a closer look at this equation, we see that it can be used to determine an expression for the steady-state angular velocity of the system:

Also notice that the exponential term in the solution gives us an expression for the time constant:

We can use these expressions to rewrite the closed-form solution in terms of the steady-state angular velocity and the time constant:

Now, let's see what happens if we move the steady-state term to the other side of the equation and take the natural logarithm of the result:

Notice that this equation has the form of a straight line, y(x) = mx + b, with time as the independent variable (x) and the natural logarithm of the difference between the steady-state response and the actual response as the dependent variable (y).  Thus, if we plot (ln [Omegass - Omega(t)]) vs. t (or equivalently, y vs. x), we should be able to estimate the time constant by the slope of the graph, since:

To get the plot of (ln [Omegass - Omega(t)]), create a MATLAB M-file with the following commands:

Run this M-file (if you are unfamiliar with M-files or plotting in MATLAB, click here), and MATLAB should display the following figure:

Notice that at small time values, the plot appears to approximate a straight line, but it becomes increasingly scattered and even has discontinuities as t gets larger.   This sporadic behavior can be explained if we take a closer look at what we're plotting in this analysis.  As time passes, and Omega(t) approaches Omegass, the fluctuations in the angular velocity readings become larger in comparison to the difference between the steady-state and actual speeds, and this makes the natural log of this value more unpredictable.  In cases where the fluctuations in Omega(t) are large enough to make Omega(t) larger than Omegass, the model attempts to take the natural log of a negative number, and thus we have a discontinuity in our plot.   Thus, to get the best estimate for the slope of the graph, we will focus on the beginning of the plot.  To do so, modify your M-file so it looks like:

When you run this file, you should get a plot that looks like:

To find the slope of this plot, enter the following in the MATLAB command window:

This command uses a least-squares method to fit a first order polynomial (i.e. a straight line) to the plot of 'graph' vs. 't'.  MATLAB responds to this command with:

This gives us the coefficients in front of descending powers of t for the best-fit line.  In our case, we asked for the linear best-fit line, so the first number is the slope (-1.7476) and the second is the y-intercept (3.5323) of the line.  You can display the least-squares linear best-fit line on the figure by entering the following in the MATLAB command window:

As we pointed out previously, knowing the slope of this line helps us determine the system's time constant:

Using the relationships developed earlier, we can use the steady-state response and the time constant to estimate b and I for the system:

Notice that these values for the parameters are very similar to those we found before using the much larger step input of 800 N-m.  The fact that we ended up with consistent estimates for the system parameters using a wide range of step inputs suggests that our first order system model approximates this particular system pretty well for the type of inputs likely to be applied to it.

 

Author:  RDM
Updated:  6/6/00