Unable to find transfer functon using model linearizer

3 visualizzazioni (ultimi 30 giorni)
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

Risposte (2)

Paul
Paul il 25 Nov 2023
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
  6 Commenti
Paul
Paul il 29 Nov 2023
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.

Accedi per commentare.


William Rose
William Rose il 26 Nov 2023
Modificato: William Rose il 26 Nov 2023
[edit: Fix a typo in the y-axis label of two plots.]
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.
  4 Commenti
Paul
Paul il 26 Nov 2023
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
William Rose il 26 Nov 2023
@Paul, thank you very much for that explanation.

Accedi per commentare.

Community

Più risposte nel  Power Electronics Control

Prodotti


Release

R2023b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by