curve fitting of a nyquist plot

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!

[Walter Roberson]

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]

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

Answer 1

[Star Strider]

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.

[Ruben Messner]

Thx! I have the time-domain data as well, that's what I measured...

So iddata instead of idfrd, the rest is the same?

I also tried this script before, but it didn't worked.

[Ruben Messner]

I tried following code:

estimateData = iddata(transpose(voltageFilteredAC(2,:)),transpose(currentFilteredAC(2,:)),Ts(2));
sys = ssest(estimateData, 10);
compare(estimateData, sys);
nyquist(sys);

comparison looks like this:

I used different sys orders, but this is the best. Not good enough I guess.

my manually calculated nyquist:

via nyquist function:

Any idea how to improve the model?

[Star Strider]

First, using time-domain data is preferable, especially if you have the input signal as well as the output. Having the input data would be important in this instance because of the initial transient.

I do not have your data, so I cannot experiment with it. One option would be to increase the system order until you get a decent fit on the compare plot. (I also do not know what your particular system is, and while ssest is usually the best option, experiment with tfest if you cannot get a good fit with ssest and an acceptable order.)

Slightly editing the data to begin wiith the start of the initial square pulse rather than the initial transient may also improve the estimate. (That would have the start time at about 0.5 seconds.) Just send that part of the data (input and output) to iddata first (as a separate iddata structure), estimate that, then include the initial transient later (with the iddata structure you are using here) if you want to.

[Ruben Messner]

Thx!

I attached my "estimateData" (iddata entitiy) in case you have time to have a look...

It contains synchronous measured current and voltage data of a small PEM fuel cell stack while being excited with a square signal. Goal is a EIS nyquist plot. So far I calculated it manually via FFT and multiple datasets with different excitation frequencies.

[Star Strider]

Thank you!

I need to run this elsewhere (currently MATLAB Online), since it takes more than the allotted 55 seconds to run here. I will keep working on this.

Also, although I assume the data are sampled, the system Identification Toolbox is running this as a continuous-time model. I am not certain if specifying a sampling interval (assumed to be uniform over the entire data set) would be helpful, however it might be worth experimenting with a sampled data set as well. If the data are not uniformly sampled, specify a sampling frequency and then use the Signal Processing Toolbox resample function to do that. The syntax I usually use is:

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

[sr,tr] = resample(s,t,Fs);

where ‘t’ is the original time vector, ‘s’ is the original signal vector, and ‘tr’ and ‘sr’ are their uniformly-resampled versions.

Meanwhile, I will keep working with the provided data. (It could take a while.)

[Les Beckham]

@Star Strider I'm thinking that you meant

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

Instead of

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

An example of why I think this

t = 0:0.0015:100;
Fs = 1/round(mean(diff(t)))

Fs = Inf

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

Fs = 667

[Star Strider]

Yes!

[Ruben Messner]

Thx a bunch!

Yes the data is uniformly sampled. And as already mentioned I have multiple datasets with different excitation frequencies. Maybe its possible to estimate the model using all 6 datasets for it?

[Star Strider]

I discovered iddata2timetable (new in R2023a) and was able to determine the sampling interval from the available times. Adding that as an additional name-value pair argument to the estimation functions helped considerably with the accuracy of the estimation, and significantly reduced the time for the estimations to converge (well within the 55-second limit to run them here, although running them in MATLAB Online for convenience). I am still not getting what I consider to be good fits to the data, so I am going to see if additional model options will provide better fits (n4sid, polynomial and ARIMA models, among others).

I worked on this for several hours yesterday. This is not a trivial problem. Since you accepted my Answer, I intend to keep working on it until I have a workable solution.

[Ruben Messner]

Thx again! That would be great! Although I don't want to use too much of your time...

Back to the question: Is there a way to train a model with multiple datasets. I got six, excited at different frequencies.

[Star Strider]

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.

.

[Ruben Messner]

Thx, this is great :)

I'll have a closer look at it soon...

Like already asked, is it possible to estimate a model based on several dataSets each containing different excitation frequencies?

[Star Strider]

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]

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]

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).