Unable to find transfer functon using model linearizer

Hello,

I tried to find the transfer function of a system, using the Model Linearizer, by setting the input pertuberation and Output Measurement Points.

After running the model linearizer for a sinestream frequency of 10 to 10000. Using the tfest function the transfer function is estimated.

The obtained transfer function, is wrong and it highly inaccurate.

Any help would be appreciated.

Thank You

Answer 1

[William Rose]

@S J Lohit Prakash,

[edit: Fix a typo in the y-axis label of two plots.]

[edit: Fix mistake in equation - thank you @Paul.]

You and @Paul are skilled in using Mat;lab's tools for modeling and system identificaiton. I am impressed.

@Paul gave a very nice compact demonstration that tfest() was not able to obtain a reasonable estimate of the system transfer function, when using time domain data from an unstable system, but tfest() was able to obtain a very good estimate of the transfer function when using ideal frequency domain data.

@Paul also showed that the frequency domain estimate of the transfer function, generated by tfestimate() using the time domain data, was very inaccurate.

Therefore we all wonder what could be done to get better esitmates.

@Paul generated a chirp signal which varied from 0.016 to 16 Hz over 10 seconds. Since the system is unstable, I wonder if the output does anything really wierd, so let's inspect it:

h = tf([0.5 3],[1 2.5 -10]);
f0 = 0.1/2/pi;
f1 = 100/2/pi;
t = 0:.001:10;
u1 = chirp(t,f0,t(end),f1);
y1 = lsim(h,u1,t);
figure
subplot(311), plot(t,u1,'-r'); grid on; ylabel('u1(t)')
subplot(312), plot(t,y1,'-r'); grid on; ylabel('y1(t)')
subplot(313), semilogy(t,y1,'-r'); xlabel('Time (s)'); ylabel('y1(t)'); grid on

The plots above show that the output, y1(t), never oscillates, despite the oscillatory input, and it grows exponentially (linear, on a semi-log plot) from t=1 to the end. It is therefore not suprising that we cannot use this time domain data to get a good estimate of the transfer function.

Would we do any better with a white noise input? I doubt it, but let's try.

u2 = rand(1,length(t))-0.5;  % uniform white noise, range [-0.5, +0.5]
y2 = lsim(h,u2,t);
figure
subplot(311), plot(t,u2,'-r'); grid on; ylabel('u2(t)')
subplot(312), plot(t,y2,'-r'); grid on; ylabel('y2(t)')
subplot(313), semilogy(t,y2,'-r'); grid on; ylabel('y2(t)'); xlabel('Time (s)')

The output is not much different with white noise input than it was with a cirp input. With white noise, the output takes a little longer to start its exponential growth phase.

By the way, here is a pole-zero plot:

figure; pzplot(h)

The denominator polynomial factors to , which means the unstable pole is associated with the frequency 2.15 rad/s = 0.34 Hz. What if we keep the input below this frequency?

f0=0.2;                 % frequency (Hz)
u3 = sin(2*pi*f0*t);    % input=sine wave at f0
y3 = lsim(h,u3,t);      % output
figure
subplot(311), plot(t,u3,'-r'); grid on; ylabel('u3(t)')
subplot(312), plot(t,y3,'-r'); grid on; ylabel('y3(t)')
subplot(313), semilogy(t,y3,'-r'); grid on; ylabel('y3(t)'); xlabel('Time (s)')

Warning: Negative data ignored

Even in this case, the output grows exponentially. One explanation for this result is that any abrupt transient contains high frequencies. In the case above, there is a sharp transient from zero activity before t=0 to a sine wave after t=0. The high frequencies in that transient are enough to excite the unstable mode of the system.

[Paul]

@William Rose

Looks like a typo: s^2 + 2.5*s - 10 = (s + 4.65)*(s - 2.15).

The pole at s = 2.15 rad/sec will result in a K*exp(2.15*t) term in the output in response to any input (except an input that has a zero in its Laplace transfrom at s = 2.15). The output should also contain an oscillating component due to the input as well, even if it's hard to see.

Is tfestimate, which does essentially the same as your cpsd approach that you showed in another question, inherently inapporiate for estimating the response an unstable system? tfestimate is based on Fourier transform concepts, but an unstable system doesn't have a Fourier transform.

[William Rose]

@Paul, thanks for that correction.

I think tfestimate() and the cpsd() approach are not going to work with time-domain data from an unstable system. I am not good enough at math to explain exactly why not. One explanation is that the exp(2.45t) term in the output swamps everything else. But that is not a very mathematical argument.

You said "an unstable system doesn't have a Fourier transform". I didn't know that. Its transfer function is well-defined at all frequencies, as shown in the Bode plot below. In what sense does it not have a Fourier transform?

h = tf([0.5 3],[1 2.5 -10]);
bode(h)

[Paul]

I should have been more precise and said the impulse response of an unstable system doesn't have a Fourier transform.

A continuous-time signal has a Fourier transform (in the sense of the Fourier transform integral converging) if the imaginary axis is contained in the region of convergence of its Laplace transform. If the signal is the impulse response of an unstable system, the ROC of its Laplace transform is somwhere to the right of the imaginary axis, hence the signal doesn't have a Fourier transform.

For example, suppose we have

syms t alpha real
x(t) = exp(-alpha*t)*heaviside(t);

x(t) has a Laplace transform regardless of the value alpha.

laplace(x(t))

ans = 

Although Matlab doesn't show it, the ROC is Re(s) > -alpha. So if alpha is positive

assume(alpha > 0)

the imaginary axis is contained within the ROC and x(t) has a Fourier transform

fourier(x(t))

ans = 

But if alpha is negative

assume(alpha < 0)

the imaginary axis is not contained within the ROC and the Fourier transform is not defined.

fourier(x(t))

ans = 

Note that entry 205 in the Fourier transform table only applies if alpha > 0.

For a discrete-time signal, the analogous result is that the unit circle has to be constained within the ROC of its z-transform for the signal to have a convergent Discrete Time Fourier Transform.

It may be interesting to define a system with one unstable pole with a very small residue so that the unstable component of the output due to the impulse response doesn't swamp out the componen of the output due to the input and see if tfestimate can be made to work.

[William Rose]

@Paul, thank you very much for that explanation.

[William Rose]

Mistake in my earlier comment fixed now. Thanks @Paul.

Answer 2

[Paul]

Hi S L,

Can you provide the parameters for your sinestream input signal? I couldn't get the Model Linearizer to run to completion for what I tried, probably because the system is unstable. Even if it did run to completion I suspect that there might be problems with the estimation because the system is unstable.

Having said that, I don't see anything in the documentation of tfest that says it assumes the system to be identified is stable (though the result you obtained is stable), but I'd be concerned about how the Model Linearizer is developing estsys in the first place.

Along these lines we can experiment with tfest as follows.

First, a stable system (same tf as yours but with the sign reversed on the third coefficient in the denominator)

h = tf([0.5 3],[1 2.5 10]);

Using time domain data

f0 = 0.1/2/pi;
f1 = 100/2/pi;
t = 0:.001:10;
u = chirp(t,f0,t(end),f1);
y = lsim(h,u,t);
sys = tfest(iddata(y(:),u(:),0.001),2)

sys =
  From input "u1" to output "y1":
     0.501 s + 2.983
  ---------------------
  s^2 + 2.497 s + 9.994
 
Continuous-time identified transfer function.

Parameterization:
   Number of poles: 2   Number of zeros: 1
   Number of free coefficients: 4
   Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                   
Estimated using TFEST on time domain data.
Fit to estimation data: 99.71%            
FPE: 2.643e-08, MSE: 2.64e-08             
 

Or using perfect frd data

w = linspace(0.1,100,1000);
hresp = frd(freqresp(h,w),w);
sys = tfest(hresp,2)

sys =
     0.5 s + 3
  ----------------
  s^2 + 2.5 s + 10
 
Continuous-time identified transfer function.

Parameterization:
   Number of poles: 2   Number of zeros: 1
   Number of free coefficients: 4
   Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                                  
Estimated using TFEST on frequency response data "hresp".
Fit to estimation data: 100%                             
FPE: 6.984e-30, MSE: 6.928e-30                           
 

Either work.

But using the transfer function in the question

h = tf([0.5 3],[1 2.5 -10]);
y = lsim(h,u,t);
sys = tfest(iddata(y(:),u(:),0.001),2)

sys =
  From input "u1" to output "y1":
       64.49 s + 100.1
  -------------------------
  s^2 + 2.264 s + 0.0007649
 
Continuous-time identified transfer function.

Parameterization:
   Number of poles: 2   Number of zeros: 1
   Number of free coefficients: 4
   Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                   
Estimated using TFEST on time domain data.
Fit to estimation data: 0.4825%           
FPE: 2.01e+15, MSE: 2.007e+15             
 

The estimation based on th time domain data is miserable.

hresp = frd(freqresp(h,w),w);
sys = tfest(hresp,2)

sys =
     0.5 s + 3
  ----------------
  s^2 + 2.5 s - 10
 
Continuous-time identified transfer function.

Parameterization:
   Number of poles: 2   Number of zeros: 1
   Number of free coefficients: 4
   Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                                  
Estimated using TFEST on frequency response data "hresp".
Fit to estimation data: 100%                             
FPE: 1.53e-31, MSE: 1.517e-31                            
 

But, the estimation based on ideal frequency domain data works fine.

So the question is whether or not the frequency domain estimate of the transfer function can be obtained from the time domain data. tfestimate with default parameters doesn't work so well, though I do think that tfestimate has an implicit assumption that that transfer function to be estimated has no poles in the right half plane. Maybe @William Rose can comment on this aspect of tfestimate (or work some magic with the parameters to get a better result).

[txy,fxy] = tfestimate(u,y,[],[],[],1000);
bode(h,frd(txy,fxy*2*pi))

Of course, tfest won't work well either given that poor estimate of the transfer function frequency response.

tfest(frd(txy,fxy*2*pi),2)

ans =
     3.325e12 s + 1.81e13
  ---------------------------
  s^2 - 1.659e05 s + 1.632e07
 
Continuous-time identified transfer function.

Parameterization:
   Number of poles: 2   Number of zeros: 1
   Number of free coefficients: 4
   Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                          
Estimated using TFEST on frequency response data.
Fit to estimation data: 6.35%                    
FPE: 9.485e+13, MSE: 9.448e+13                   
 

[S J Lohit Prakash]

Hello,

I tried changing the transfer function to , a stable one, and from using the in-built Model linearizer

with the sinestream parameters of

The obtained results are still inaccurate.

Moreover, i am using this to find the transfer function of a cascaded power converter using the simulink model, instead of MATLAB Code. So, it would be helpful to have some insights on how to find the transfer function of a power converter using a Simulink model.

[Paul]

For this stable model, I used a sinestream with 40 frequencies, logarithmically spaced from 0.1 to 100 rad/sec (based on the Bode plot of the transfer function, starting at 10 Hz or ~60 rad/sec is not sufficient), with all other parameters default. The result of tfest was not good

>> tfest(estsys1,2)

ans =

From input "Constant" to output "Transfer Fcn":

-0.04411 s + 0.007366

-------------------------

s^2 + 0.01029 s + 0.01077

Continuous-time identified transfer function.

Parameterization:

Number of poles: 2 Number of zeros: 1

Number of free coefficients: 4

Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:

Estimated using TFEST on frequency response data "Frequency response estimation using the input signal "in_sine1"".

Fit to estimation data: 69.95%

FPE: 0.1523, MSE: 0.1246

But, when I plotted estsys1 I noticed the that first few points were obviously not good estimates of the frequency response. I deleted the first five points and got a much better result:

>> tfest(frd(squeeze(estsys1.ResponseData(:,:,5:end)),estsys1.Frequency(5:end)),2)

ans =

0.4966 s + 3.019

---------------------

s^2 + 2.509 s + 10.04

Continuous-time identified transfer function.

Parameterization:

Number of poles: 2 Number of zeros: 1

Number of free coefficients: 4

Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:

Estimated using TFEST on frequency response data.

Fit to estimation data: 98.86%

FPE: 5.933e-06, MSE: 4.746e-06

It's nearly perfect.

[William Rose]

@S J Lohit Prakash,

I am again impressed by how @Paul uses the system identification routines, and by his idea, and efficient implementation, of fitting the frequency response without the 5 lowest frequencies. It worked really well.

I was not familiar with sinestream until @S J Lohit Prakash and @Paul mentioned it.

By the way, the frequencies that @Paul excluded were:

input = frest.Sinestream('Frequency',logspace(-1,2,40));
fprintf('Freqs 1-5 = %.3f, %.3f, %.3f, %.3f, %.3f Hz.\n',...
    input.Frequency(1:5)/2/pi)

Freqs 1-5 = 0.016, 0.019, 0.023, 0.027, 0.032 Hz.

Plot the input signal to see its duration:

plot(input)

@Paul's chirp signal was only 10 s long, and he got a good estimate of the stable system transfer function using the 10 s chirp input.

With a signal duration of 1500 s, it is not surprising that we get a good estimate of this stable system. (@S J Lohit Prakash's and @Paul's sinestream signals could have been shorter or longer, if they used non-default options for sinestream().)

In a real experiment, there are often tight limits for the duration of data collection. (With human or other animal subjects, for example.) The chirp was much more efficient than the sinestream at identifying the system quickly. Random noise is also efficient.

I wonder why the transfer function estimate was bad at the low frequencies in the sinestream.

[Paul]

Upon further inspection, it looks like only the first three frequencies needed to be excluded. Here's the plot.

I poked around a bit, and I'm still not sure why the frd estimation was so poor at those first three frequencies.

[William Rose]

Thanks, @Paul.

[Paul]

Now I'm even more puzzled by what's going on.

Create the Sinestream object with specified frequency vector and defaults otherwise

w = logspace(-1,2,40);
u = frest.Sinestream('Frequency',w);

Create a time vector corresponding to the first NumPeriods for the first frequency point

t = (0:u.NumPeriods*u.SamplesPerPeriod-1)/(u.NumPeriods*u.SamplesPerPeriod)*u.NumPeriods*2*pi/w(1);
t = t.';

Create a sine wave input for that time vector

uu = 1e-5*sin(w(1)*t);

Verify it's the same as the first NumPeriods of u

figure
plot(u),grid
c = get(gca,'Children');
c.Marker = 'o';
hold on;plot(t,uu,'x'),xlim([0 270])

Define the system transfer function

h = tf([0.5 3],[1 2.5 10]);

Generate the system output

yy = lsim(h,uu,t);

Generate the frequency response

HH = fft(yy)./fft(uu);
ww = (0:numel(HH)-1)/numel(HH)*2*pi/t(2);

The fifth element is the one we want

ww(5)

ans = 0.1000

Generate and plot the frequency response and overlay with the point estimate

[mag,phase,wout] = bode(h); mag = squeeze(mag); phase = squeeze(phase);
figure
subplot(211)
semilogx(wout,db(mag),w(1),db(abs(HH(5))),'.','MarkerSize',20),grid
axis padded
subplot(212)
semilogx(wout,(phase),w(1),180/pi*(angle(HH(5))),'.','MarkerSize',20),grid
axis padded

I think frestimate does a bit more than a simple FFT ratio, but not sure why it has such a problem at that first frequency point. I haven't investigated the second or third frequency points.