Higher order systems tend to be difficult to model because, in general, there are no analytical formulas to compute the parameters of a model directly from the experimental data. However it is often possible to break down a modeling problem into the problems of deriving the parameters of first and second order systems.

In this example, the Bode plot for the frequency response of a given set of experimental data  will be analyzed similar to the previous examples but with unknown transfer function and steady state parameters.  Download the data from the link provided:

                 higherdata.m

This file contains only the data from the experiment.  This data in this file must be modified to be stored in a MATLAB variable.

Open the higherdata.m file place the following around the data:

data = [ (data) ];

Be sure to put the close bracket and semi-colon at the end of the data, then save the file and run it.

The data must now be turned and stored as a Frequency Response Data (FRD) model.  To store the data as a FRD model, the frequency and response must be stored separately with the following commands"

freq = data(:,1);
resp = data(:,2);

The colon used above is a wildcard term.  Here it represents all rows.  So, data(:,1), means all rows in column one. 

Now that the frequency and corresponding response are separate, they can be stored using the FRD command.  Remember when using the FRD command to always put the response first and the frequency after.

 frd_model = frd(resp, freq);

With this model we will use the Bode plot of the FRD model to create and compare a second order system and its corresponding Bode plot.  After adding the Bode command, the end of your M-file should look something like the following.

The following will be the resulting Bode plot.

Notice on the above magnitude plot the fairly level graph from 0.1 to 3 followed by a peak and a drop.  At high frequencies, the magnitude graph decreases at about 40 dB/dec.  Also, notice that the phase plot decreases from 0 to -180.  These are all characteristics of a second order system.    So to start the modeling process, construct a new second order transfer function to capture these features of the Bode plot. 

To begin we will use the following second order transfer function model:

Using the graph we can locate the frequency, w_p, and magnitude, Mp, of the peak and the DC gain, k_dc.  By dragging the cursor along the plot of the graph, one can find the different values and frequencies.  Find the magnitude at the lowest frequency of the graph which will be a good approximation of the DC gain.  Also find the frequency and magnitude of the peak as shown.

 The relationship between features of the Bode plot and the parameters of the second order system are as follows:

Where w_n is the natural frequency where the graph begins to turn and ζ is the damping ratio.

From the graph the  k_dc = 1 (i.e. 0 dB), w_p = 2.54 rad/sec, and Mp = 1.33 (i.e. 2.45 dB).

Adding the following commands into the M-file will define the parameters and create a second order transfer function model.  It is helpful to leave out the semi-colon at the end of the lines defining w_n and ζ because those values, which will appear in the command window, may be needed later

        kdc = 10^(0/20);
        wp = 2.54;
        Mp = 10^(2.45/20);
        zeta = sqrt(0.5 - 0.5 * sqrt(1 - (1/Mp^2)))
        wn = wp / sqrt(1 - 2 * zeta^2)
        s = tf('s');
        second_tf = kdc /((s/wn)^2 + 2 * zeta * (s/wn) + 1);
        figure(1);
        bode(frd_model, second_tf);

The following Bode plot will appear:

The second order model clearly captures some of the features of the FRD model, specifically the DC gain, the magnitude peak, and the slope of the high frequency asymptote of the magnitude plot.  The second order model also matches the phase plot at very high and low frequencies. 

Between 10 and 1000 rad/sec there is a large difference between the second order plot and the frequency response data model. In that region the magnitude plot of the FRD model decreases at a rate of 20 dB/dec rather then the expected 40 dB/dec of the second order plot.  Also in that range there is a large phase discrepancy.  The phase of the FRD model is significantly greater then the second order model.

These differences suggest adding a so called lead compensator transfer function to the second order model. A lead compensator has the following form:

An example of a lead compensator can be seen below:

The lead compensator has a zero at -z and a pole at -p, which are defined by:

      
            

Here, f is the amount of phase at the phase peak and w_l is the frequency of the phase peak.

The next step of the system identification process is to find the frequency at which the phase difference between the FRD model and the second order model is greatest.  This frequency will be the w_l of the lead compensator, and the difference in phase will be the f of the lead compensator.

A good estimate for these values would be: w_l = 30 rad/sec, j = (-91) - (-175) = 84 degrees.  To create the transfer function model of the lead compensator, add the following commands to the M-file to define the appropriate parameters.  Remember that MATLAB does trigonometry in radians, and so degrees must be multiplies with pi/180.  Then graph the Bode plot the of the second order transfer function multiplied with the lead compensator as shown.

        wlead = 30;
        philead = 84 * (pi / 180);
        a = (1 - sin(philead))/(1+sin(philead));
        zlead = wlead * sqrt(a);
        plead = wlead / sqrt(a);
        leadcomp = (s/zlead + 1)/(s/plead + 1);
        figure(1);
        bode(frd_model, second_tf * leadcomp);
        title('highermodel');

Run the M-file and the following Bode plot will result:

 

Notice that the peak is higher than expected.  This is because of the influence of the lead compensator.  Use just the resulting values of ζ =0.4127 and w_n=3.1282 without the equations as a inital values for further iterations. 

A good starting point for a iteration would be to find the difference between peaks of the FRD model and higher order model and multiply that value to ζ.  Remember to convert from dB.  Here the difference it found to be 6.11 dB (i.e. 2.02).  Multiplying zeta by this value gives the following plots.

Above, the peaks are nearly coincident and may be a sufficient model for many cases.  Often times, matching the lower frequency values where the magnitude is highest is more important then matching the high ones where the magnitude drops off rapidly.  Proceeding two iteration as follows results in an even better graph.

% Iteration 1
zeta = zeta * 10^((8.56-2.45)/20)
second_tf = 1/((s/wn)^2 + 2*zeta*(s/wn) + 1);

% Iteration 2
wlead = 35;
philead = 81 * (pi / 180);
a = (1 - sin(philead))/(1+sin(philead));
zlead = wlead * sqrt(a);
plead = wlead / sqrt(a);
leadcomp = (s/zlead + 1)/(s/plead + 1);

% ...Iteration X
zeta = zeta / 10^((3.0)/20)
second_tf = 1/((s/wn)^2 + 2*zeta*(s/wn) + 1);
figure(1);
bode(frd_model, second_tf * leadcomp,{1, 10});
title('highermodel');

Further adjustments may be necessary depending on problem requirements.  If a specific area needs to be very exact, one can zoom in on the graph with the addition to the Bode command shown below, and iterations can continue.

        bode(frd_model, second_tf * leadcomp,{1, 10});

This would result in a closer view of desired points such as the one below:

Carnegie Mellon University | University of Michigan

Mechanical Engineering Department