Is there a way to determine an equivalent circuit for Nyquist diagrams (I have the data set) using Matlab (or Simulink)?
이전 댓글 표시
I am working on gathering some Nyquist diagrams experimentally and facing some difficulties on simulating equivalent circuits for my data. What I have done so far is simply trying to come up with an idea of an equivalent circuit and seen if it fits the case. Well this is obscure and it is not always representative of what has been going on. So is it possible to use Matlab to construct/simulate equivalent circuits from experimental data?
채택된 답변
I am not certain what you want witth respect to an equivalent circuit. I believe you intend circuit synthesis or circuit realisation. If so, this may do what you want
First, using your data, identify the system using the System Identification Toolbox, specifically using the iddata and tfest functions to return a transdfer function.
When you have the transfer funciton that gives an appropriate fit to your data (check this using the using the compare function) you can then use circuit synthesis to identify the component values.
Illustrating that with a Cauer expansion (you will probably need to extract the transfer function numerator and denominator polynomials from the identified transfer function and then use the poly2sym function on each of them first) —
syms s
H = (2*s^5+40*s^3+128*s)/(s^4+10*s^2+9)
[Hn1,Hd1] = numden(H)
T1 = quorem(Hn1,Hd1)
H2 = Hn1/Hd1 - T1
[Hn2,Hd2] = numden(H2)
T2 = quorem(Hd2,Hn2)
H3 = Hd2/Hn2 - T2
[Hn3,Hd3] = numden(H3)
T3 = quorem(Hd3,Hn3)
H4 = Hd3/Hn3 - T3
[Hn4,Hd4] = numden(H4)
T4 = quorem(Hd4,Hn4)
H5 = Hd4/Hn4 - T4
[Hn5,Hd5] = numden(H5)
T5 = quorem(Hd5,Hn5)
That can be coded as a loop —
denpoly = [1 1]; % Define First (Used To Terminate The ‘while’ Loop)
syms s
H = (2*s^5+40*s^3+128*s)/(s^4+10*s^2+9)
k = 0;
H = 1/H;
while numel(denpoly) > 1 % Loop Version (Can Be A Function)
k = k + 1;
% disp("k = "+k)
[Hn,Hd] = numden(1/H);
denpoly = sym2poly(Hd);
Comp(k) = quorem(Hn,Hd);
H = Hn/Hd - Comp(k);
if k > 50 % Fail-Safe To Prevent Infinite Recursion
disp("Exiting loop (k = "+k+" )")
break
end
end
fprintf(['Component Values: \n',repmat('\t\t\t%s\n',1, numel(Comp)), '\n'], Comp)
These would be configured as:
>---- 2 H ----- 40/9 H ---- 70/9 H ---->
| |
1/20 F 9/40 F
| |
>-------------------------------------->
This is a textbook example, and gives the same result as in W-K Chen, Passive and Active Filters Wiley 1986 (ISBN 0-471-82352-X), P. 113. There are other approaches that may also be appropriate, such as the Foster expansion.
This may not always work as desired in ‘real world’ transfer functions, in that the component (‘Comp’ vector) values may not be simple functions of s as they are here, and may instead be vectors. Those results would not be valid.
(There is some sort of bug in Answers such that the results of symbolic calculations do not appear in the posted results, although they appear in the results while being run or edited. This has been reported, although apparently not yet fixed.)
.
댓글 수: 5
I did not see the .txt file previously.
T1 = readtable('EIS_Data.txt')
figure
semilogy(T1{:,1}, T1{:,2:end})
grid
It would help significantly to know what the variable names are.
Arthur
2024년 10월 9일
I doubt that it is possible to do what you want.
The System Identification Toolbox will not permit impropeer transfer functions (more zeros than poles), however the invfreqz funciton will, and that allows the use of the example loop in my original answer. Plotting the imaginary part of thee response function, there are two obvious polees and two obvious zeros, however eveene allowing for three of each, the result is not good, and it is likely not poossible to estimate the equivalent circuit compnent values.
T1 = readtable('EIS_Data.txt');
T1.Properties.VariableNames = {'Index','Frequency (Hz)','Z’','-Z’’','Z','Phase (°)','Time(s)'};
T1 = flipud(T1)
Zz = T1.Z .* exp(1j*deg2rad(T1.('Phase (°)')))
[pks,plocs] = findpeaks(imag(Zz))
[vys,vlocs] = findpeaks(-imag(Zz))
figure
semilogx(T1.('Frequency (Hz)'), imag(Zz))
xline(T1{vlocs,2}, "--", "Zero at "+T1{vlocs,2}+" Hz")
xline(T1{end,2}, "--", "Zero at \infty")
xline(T1{1,2}, "--", "Pole at "+0+" Hz")
xline(T1{plocs,2}, "--", "Pole at "+T1{plocs,2}+" Hz")
grid
Ts = 1/(2*T1{end,2})
FRD = idfrd(Zz, T1{:,2}, 0, FrequencyUnit="Hz")
H = tfest(FRD, 3, 3)
figure
compare(FRD, H)
[n,d] = invfreqz(T1.Z, pi*T1{:,2}/max(T1{:,2}), 3, 2)
syms s
H = poly2sym(n, s) / poly2sym(d, s);
H = vpa(H, 5)
denpoly = [1 1]; % Define First (Used To Terminate The ‘while’ Loop)
% H = (2*s^5+40*s^3+128*s)/(s^4+10*s^2+9)
k = 0;
H = 1/H;
while numel(denpoly) > 1 % Loop Version (Can Be A Function)
k = k + 1;
% disp("k = "+k)
[Hn,Hd] = numden(1/H);
denpoly = sym2poly(Hd);
Comp(k,:) = quorem(Hn,Hd);
H = Hn/Hd - Comp(k);
if k > numel(n) % Fail-Safe To Prevent Infinite Recursion
disp("Exiting loop (k = "+k+" )")
break
end
end
Comp
fprintf(['Component Values: \n',repmat('\t\t\t%s\n',1, numel(Comp)), '\n'], Comp)
[n,d] = invfreqz(T1.Z, pi*T1{:,2}/max(T1{:,2}), 4, 3)
syms s
H = poly2sym(n, s) / poly2sym(d, s);
H = vpa(H, 5)
denpoly = [1 1]; % Define First (Used To Terminate The ‘while’ Loop)
% H = (2*s^5+40*s^3+128*s)/(s^4+10*s^2+9)
k = 0;
H = 1/H;
while numel(denpoly) > 1 % Loop Version (Can Be A Function)
k = k + 1;
% disp("k = "+k)
[Hn,Hd] = numden(1/H);
denpoly = sym2poly(Hd);
Comp(k,:) = quorem(Hn,Hd);
H = Hn/Hd - Comp(k);
if k > numel(n) % Fail-Safe To Prevent Infinite Recursion
disp("Exiting loop (k = "+k+" )")
break
end
end
Comp
fprintf(['Component Values: \n',repmat('\t\t\t%s\n',1, numel(Comp)), '\n'], Comp)
This is the best I can do with these data.
(I initially began this comment yesterday, however Windows 11 crashed — again — so I lost everytthing and had to start over. My apologies for the delay.)
.
Star Strider
2024년 10월 15일
As always, my pleasure!
This is not exactly a straightforward problem. I am not certain how best to solve it, however hypothesizing a two-pole and two-zero equivalent ciircuit of some sort (there are likely several ways to configure it) and then doing nonlinear parameter estimation to identify the component values would likely be the way to go. The pole at the origin is likely close to it (so the negative real component is likely very small). however the zero at 0.0032 Hz and the pole at 0.015 Hz have much larger (greater magnitude) negative real parts, since the zero value does not go to zero and the pole does not go to infinity, so they are relatively far away from the imaginary axis.
The fun lies in estimating them! Fortunately, MATLAB has a number of nonlinear parameter estimation and optimisation algorithms (including Global Optimization Toolbox functions such as ga) that make that task much easier.
No worries about the time — we all have lives away from MATLAB Answers!
.
추가 답변 (3개)
Hi @Arthur FdA
I'm unfamiliar with your circuit system. However, if you have the data in the form of frequency-response, then you can apply the frd() and nyquistplot() commands
data = [
0.0010 1.562e+01-1.9904i
0.0018 1.560e+01-2.0947i
0.0034 1.513e+01-3.3670i
0.0062 1.373e+01-5.4306i
0.0113 1.047e+01-7.5227i
0.0207 5.829e+00-7.6529i
0.0379 2.340e+00-5.6271i
0.0695 7.765e-01-3.4188i
0.1274 2.394e-01-1.9295i
0.2336 7.216e-02-1.0648i
0.4281 2.157e-02-0.5834i
0.7848 6.433e-03-0.3188i
1.4384 1.916e-03-0.1740i
2.6367 5.705e-04-0.0950i
4.8329 1.698e-04-0.0518i
8.8587 5.055e-05-0.0283i
16.2378 1.505e-05-0.0154i
29.7635 4.478e-06-0.0084i
54.5559 1.333e-06-0.0046i
100.0000 3.967e-07-0.0025i
];
freq = data(:,1);
resp = data(:,2);
sys = frd(resp, freq);
figure
np = nyquistplot(sys); grid on
figure
bd = bodeplot(sys); grid on, grid minor
William Rose
2024년 10월 4일
1 개 추천
Yes you can. The experimental data needed for Nyquist and Bode plots is the same: a set of magnitudes and phases, or (equivalently) real and imaginary parts, at a set of frequencies, f. The Bode plot shows the frequencies on the plot. The frequencies are not explicit on the Nyquist plot, but I assume you have the frequencies in your experimental data.
You make a guess, based on experience and the appearance of th Bode plot or Nyquist plot, of the circuit type (high pass, low pass, etc) and number of components. Then you compute the expected theoretical frequency response 
, which will be a complex function involving the values of R. C, L, etc. in your circuit. Then you compare it to the measured frequency response, 
, which is also complex, and which you have measured in your experiment at a set of frequencies. Then you adjust the values of R, C, L, etc., in order to minimize the error between 
and 
. The error is the sum of the squared distances in the complex plane between 
and 
at the measured frequencies. You can use Matlab's fmincon() to find the values of the unknown parameters (R, L, C,...) that minimize the error.
, which will be a complex function involving the values of R. C, L, etc. in your circuit. Then you compare it to the measured frequency response, , which is also complex, and which you have measured in your experiment at a set of frequencies. Then you adjust the values of R, C, L, etc., in order to minimize the error between and . The error is the sum of the squared distances in the complex plane between and at the measured frequencies. You can use Matlab's fmincon() to find the values of the unknown parameters (R, L, C,...) that minimize the error.If the system is a filter, where the input and output have the same units (typically volts), then you have no informatopn about the overal circuit impedance because you don;t know what current is flowing. Then there may be many equivalent solutions that give the same frequency response. If this is true, then there is no unique solution forr the component values. So you fix one of the component values (for example, you specify that you want impedance Z=100 Ω at DC, which fixes one of the R values in the circuit). Then you solve for the other unknown parameters.
If your best-fit model response does not look like the experimental data, then try a different theoretical circuit.
댓글 수: 1
Here is an example of what I meant in my answer above.
Files simulatedFilterResponse{0,1,2}.txt contains Nyquist plot data from three simulated experiments. Each file has 3 columns: frequency (rad/s), real part of Gmeas, imaginary part of Gmeas, where Gmeas is the measured frequency response. The measurement noise has relative amplitude 0, 1, 2 on files 0, 1, 2 respectively.
The main script, fitFreqResponse.m, minimizes the error between a model and the data, by adjusting the values of R, L, and C for a second order filter circuit. After finding the best-fit values, the script displays the best-fit values, and plots the measured and best-fit transfer functions. The script calls external function sumsqerr(p,w,Gmeas), where p=[R;L;C] are the adjustable parameters, w is the vector of frequencies, and Gmeas is the experimental data: a complex vector of measured frequency responses. Function sumsqerr() calculates the model transfer function and the sum squared error betwen the model and the measured data.
If your experimental data does not look like a second order lowpass filter, then adjust the line Gmodel=..., in function file sumsqerr1.m. Also change the line Gfit=... in the main script, so that Gfit and Gmodel match. (Gfit is used to plot the best-fit transfer function after the best fit circuit component values have been found.)
The simulated Nyquist plot data is from the circuit below, with R=100 Ω, L=0.1125 H, C=22.51 
.
.
When I estimate R, L, C from the noise-free Nyquist plot data (simulatedFilterResponse0.txt), I get exactly correct estimates for R, L, C. When I edit the main script so that it reads a file of noisy Nyquist plot data (simulatedFilterResponse1.txt), the estimated R, L, C values are not quite the same. See below.
fitFreqResponse
Arthur
2024년 10월 7일
0 개 추천
Thank you for your responses. This is new to me, my they are widening my mind. Some of you said you are unfamiliar with my data. So here it goes as attachment. What I am trying to do is to propose a circuit that represents it. This kind of data is from electrochemical phenomena and I usually try to guess the circuit and evaluate some reports about errors to decide if I can accept my guess or if I should reject it and propose another circuit. This is my first attempt at doing it numerically so I would really appreciate it if you could help me out by using my experimental data. The Nyquist Diagram and the Bode Plot are as follows.
Nyquist Diagram
Bode Plot
댓글 수: 2
Star Strider
2024년 10월 7일
It would be helpful to have your actual data. The best option is time-domain data, however if all you have is frequency-domain data, that could also work.
With time-domain data we can simulate time-domain and frequency-domain results, however with only frequency-domain data, we can only simulate frequency-domain results.
Arthur
2024년 10월 8일
카테고리
도움말 센터 및 File Exchange에서 Simscape Electrical에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
