curve fitting of a nyquist plot

Question

Ruben Messner 2023-8-29

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2014266-curve-fitting-of-a-nyquist-plot

评论： Star Strider 2023-9-6

采纳的回答： Star Strider

hi,

I have impedance data points in a nyquist plot:

Anybody an idea how to fit a curve to that data? The data is stored in a array with complex impedance values. I tried a few things, so far without success...

Highly appreciate your support!

2 个评论
显示无隐藏无

Walter Roberson 2023-8-29

When I follow the plot by eye, I get the impression that for any given Re value that the Im value is either positive or negative, but not both. I do not notice any locations in which the same Re has both a positive and negative Im reading.

Am I mistaken about that, or am I correct but it is an artifact, or is it an actual property of the system?

The question becomes whether there is a continuous curve that is vaguely oval-ish, or if instead it is rather discontinuous and we have to find something in the single parameter Re that jumps back and forth between + and -

Ruben Messner 2023-8-29

Thx for your comment! The inductive half circle is an artifact of the measurement...

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2023-8-29

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2014266-curve-fitting-of-a-nyquist-plot#answer_1295621

I doubt that there is any direct way to fit it, for example to a nyquist plot. You would have to estimate the system first. The data in the plot are apparently from a Fourier transform, so you have frequency response data. If frequency response data are all you have (rather than the original time-domain data), then use the System Identification Toolbox to identify it, beginning with the idfrd function and going from there, probably ssest since it is the most robust in my experience, although tfest is also an option. Use the compare function to understand graphically how well the estimated system fits the data. If you have the original time-domain data, start with the iddata function.

15 个评论
显示 13更早的评论隐藏 13更早的评论

Star Strider 2023-9-1

在 MATLAB Online 中打开

pemfcDataSet.mat

I am not happy with these results (I expected better), however they are the best I can come up with. This approach makes the code relatively straightforward to experiment with. The one improvement I added was to add the sampling interval to the estimate functions. This speeds up the code considerably, and provides much better results than considering it as a continuous model.

I included a Fourier transform of the signals, demonstrating band-limited noise. A lowpass filter with a cutoff of 2 Hz should elliminate most of it. Teh best way to do that would be:

uy = lowpass(TT{:,[1 2]}, 2, 1/Ts, 'ImpulseResponse','iir');

u1 = uy(:,1);

y1 = uy(:,2);

You have the relevant signals, so you can create a timetable from the data to replace ‘TT’ here, since the iddata2timetable function was introduced in R2023a. (I needed to do that to see what the signals were, and to get the sampling interval, ‘Ts’.) I printed a segment of it here so you can create one to match it. You can either create it directly, or first create a table and then use the table2timetable function.

These appear to be the only models (I considered several and read the documentation on them) that are appropriate for your data —

LD = load('pemfcDataSet.mat')

LD = struct with fields:

estimateData: [64516×1×1 iddata]

estimateData = LD.estimateData;

TT = iddata2timetable(estimateData); % New In R2023

TT(1:5,:)

ans = 5×2 timetable

Time u1 y1 _____________ _________ _______ 7.75e-05 sec 0.12462 0.45093 0.000155 sec 0.11658 0.46701 0.0002325 sec 0.11658 0.45897 0.00031 sec 0.11859 0.45897 0.0003875 sec 0.0080402 0.42681

% TT.Properties

% format long

% t_stats = [mean(diff(TT.Time)); mode(diff(TT.Time)); std(diff(TT.Time)); min(diff(TT.Time)); max(diff(TT.Time))]

% format short

[h,m,s] = hms(TT.Time);

Ts = mean(diff(s))

Ts = 7.7500e-05

[FTs1(:,1),Fv] = FFT1(TT{:,1}, s);

[FTs1(:,2),Fv] = FFT1(TT{:,2}, s);

figure

semilogx(Fv, mag2db(abs(FTs1)*2))

grid

xlabel('Frequency (Hz)')

ylabel('Power (dB)')

legend('u_1','y_1', 'Location','best')

xline(2, '--k', '2 Hz')

figure

plot(TT.Time, TT{:,[1 2]})

grid

legend('Input','Output', 'Location','best')

% estimateData = iddata(transpose(voltageFilteredAC(2,:)),transpose(currentFilteredAC(2,:)),Ts(2));

% OPTIONS: 'ss', 'tf', 'armax'

Model = 'ss' % Select Model To Run

Model = 'ss'

switch Model

case 'ss'

tic

sys = ssest(estimateData, 6, 'Ts',Ts)

toc

% % % Order 6 works (sort of)

case 'tf'

tic

sys = tfest(estimateData, 6, 'Ts',Ts)

toc

% % % Order 4 just models the noise

case 'armax'

tic

sys = armax(TT.u1, TT.y1, [3 3 3 1], 'Ts',Ts)

toc

end

sys = Discrete-time identified state-space model: x(t+Ts) = A x(t) + B u(t) + K e(t) y(t) = C x(t) + D u(t) + e(t) A = x1 x2 x3 x4 x5 x6 x1 1.001 -0.0004582 -0.005969 0.0008462 -0.001282 -0.005756 x2 -0.005908 1.079 0.6373 0.1957 -0.3177 0.01806 x3 -0.009047 0.0004081 0.462 0.7997 -0.03163 -0.5867 x4 -0.005108 -0.01231 0.3179 -0.4289 0.8803 -0.02475 x5 -0.0007279 -0.03658 0.1868 -0.5846 -0.4707 -0.5666 x6 0.00373 0.06352 0.02189 -0.008845 0.2248 -0.3591 B = u1 x1 0.0006778 x2 -0.06637 x3 0.03677 x4 0.005505 x5 0.002452 x6 0.01331 C = x1 x2 x3 x4 x5 x6 y1 69.58 -0.6518 0.8493 -0.002054 -0.173 0.2586 D = u1 y1 0 K = y1 x1 0.00104 x2 0.03362 x3 -0.001961 x4 0.00103 x5 -0.002627 x6 0.001252 Sample time: 7.75e-05 seconds Parameterization: FREE form (all coefficients in A, B, C free). Feedthrough: none Disturbance component: estimate Number of free coefficients: 54 Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties. Status: Estimated using SSEST on time domain data "estimateData". Fit to estimation data: 98.38% (prediction focus) FPE: 2.556e-05, MSE: 2.554e-05

Elapsed time is 17.968366 seconds.

figure

compare(estimateData, sys)

figure

nyquist(sys)

function [FTs1,Fv] = FFT1(s,t)

s = s(:);

t = t(:);

L = numel(t);

Fs = 1/mean(diff(t));

Fn = Fs/2;

NFFT = 2^nextpow2(L);

FTs = fft((s - mean(s)).*hann(L), NFFT)/sum(hann(L));

Fv = linspace(0, 1, NFFT/2+1)*Fn;

Iv = 1:numel(Fv);

FTs1 = FTs(Iv);

end

I will be glad to contnue working on this as long as necessary to get a good result, however this is the best I can do for now. If you have any further ideas, please share them.

.

Star Strider 2023-9-1

在 MATLAB Online 中打开

Thank you!

Together, no, because the sampling frequency is an important parameter. Models with different sampling frequencies would have different parameters, however the model order can be the same.

An example with a completely synthetic signal demonstrates that the same essential signal (the lengths are different) with different sampling frequencies and the same model order will give models different parameters —

Fs1 = 1000;

L = 500;

t1 = linspace(0, Fs1*L, Fs1*L+1).'/Fs1;

u1 = [0; 1; zeros(L*Fs1-1,1)];

y1 = exp(-0.005*t1) .* sin(2*pi*t1*0.05);

data1 = iddata(y1, u1, 1/Fs1);

figure

plot(t1, y1)

grid

sys1 = ssest(data1,3)

Warning: For transient data (step or impulse experiment), make sure that the change in input signal does not happen too early relative to the order of the desired model. You can achieve this by prepending sufficient number of zeros (equilibrium values) to the input and output signals. For example, a step input must be represented as [zeros(nx,1); ones(N,1)] rather than ones(N,1), such that nx > model order.

sys1 = Continuous-time identified state-space model: dx/dt = A x(t) + B u(t) + K e(t) y(t) = C x(t) + D u(t) + e(t) A = x1 x2 x3 x1 -8.47e-08 -0.3143 31.39 x2 0.3141 -0.01 -0.06021 x3 0 0 -2000 B = u1 x1 9.987e-16 x2 1.266e-12 x3 1.721e-14 C = x1 x2 x3 y1 222.7 -0.1064 3.496 D = u1 y1 0 K = y1 x1 2.647 x2 -1550 x3 -8.715e-15 Parameterization: FREE form (all coefficients in A, B, C free). Feedthrough: none Disturbance component: estimate Number of free coefficients: 18 Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties. Status: Estimated using SSEST on time domain data "data1". Fit to estimation data: 100% (prediction focus) FPE: 1.089e-28, MSE: 1.089e-28

figure

compare(data1, sys1)

Fs2 = 500;

t2 = linspace(0, Fs2*L, Fs2*L+1).'/Fs2;

u2 = [0; 1; zeros(L*Fs2-1,1)];

y2 = exp(-0.005*t2) .* sin(2*pi*t2*0.05);

data2 = iddata(y2, u2, 1/Fs2);

figure

plot(t2, y2)

grid

sys2 = ssest(data2, 3)

Warning: For transient data (step or impulse experiment), make sure that the change in input signal does not happen too early relative to the order of the desired model. You can achieve this by prepending sufficient number of zeros (equilibrium values) to the input and output signals. For example, a step input must be represented as [zeros(nx,1); ones(N,1)] rather than ones(N,1), such that nx > model order.

sys2 = Continuous-time identified state-space model: dx/dt = A x(t) + B u(t) + K e(t) y(t) = C x(t) + D u(t) + e(t) A = x1 x2 x3 x1 -1.303e-06 -0.3173 78.2 x2 0.3112 -0.009999 -0.06288 x3 0 0 -1000 B = u1 x1 0.001796 x2 -2.002e-06 x3 -0.02296 C = x1 x2 x3 y1 156 -0.1211 12.2 D = u1 y1 0 K = y1 x1 2.499 x2 -709.1 x3 7.823e-17 Parameterization: FREE form (all coefficients in A, B, C free). Feedthrough: none Disturbance component: estimate Number of free coefficients: 18 Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties. Status: Estimated using SSEST on time domain data "data2". Fit to estimation data: 100% (prediction focus) FPE: 1.502e-23, MSE: 1.502e-23

figure

compare(data2, sys2)

The pole-zero placements and other characteristics of the two models will be different as werll.

So, with that in mind, you can do something like this with your different inputs and outputs. This uses ssest, however you will need to determine the correct procedure for your data.

.

Ruben Messner 2023-9-6

thx for that demonstration!

Yeah so merging multiple iddata works only if the sampling frequency is the same I guess...

I also experimented with the system estimation and found out, that removing the offset of the current signal helps a lot. Here is an example:

I calculated the nyquist plot for all dataSets and that's what they look like:

The numbers don't look too bad. Some of them are pretty similar to my other calculations. Maybe a dumb question: Is the nyquist always symmetric when I use nyquist(sys)?

Star Strider 2023-9-6

As always, my pleasure!

Removing the D-C offset eliminates a significant problem in the estimation. I usually do that by subtracting the mean of the signal. (I did not do that because they were your data and you need to make those decisions.) There is absolutely nothing wrong with doing that providing it does what you want.

I believe the Nyquist diagram should always be symmetric, at least with respect to the real axis. (It is of course not going to be symmetric with respect to the imaginary axis.)

I am happy that it is working and matching your other calculations! I was concerned that it might not work that well, considering that fitting them is the objective here (although that cannot be done directly, at lest in my experience).

请先登录，再进行评论。