Is there a way to determine an equivalent circuit for Nyquist diagrams (I have the data set) using Matlab (or Simulink)?

29 次查看（过去 30 天）
Arthur FdA
Arthur FdA 2024-10-4
评论： Star Strider 2024-10-15
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?

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

采纳的回答

Star Strider
Star Strider 2024-10-4
Ran in:
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)
Component Values: 2*s s/20 (40*s)/9 (9*s)/140 (70*s)/9
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 个评论
显示 3更早的评论
Arthur FdA
Arthur FdA 2024-10-15
编辑：Arthur FdA 2024-10-15
Wow! @Star Strider Thank you so much! I didn't reply to you before because I've just seen your answer. Now I will study your code and try to apply to other cases I encountered.
Eventhough it is not possible to estimate an equivalent circuit directly, it gives me a way to work on that.
Star Strider
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 个）

Sam Chak
Sam Chak 2024-10-4
Ran in:
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
William Rose 2024-10-4
@Arthur FdA,
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.
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 个评论
William Rose
William Rose 2024-10-7
Ran in:
@Arthur FdA,
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
Local minimum possible. Constraints satisfied. fmincon stopped because the size of the current step is less than the value of the step size tolerance and constraints are satisfied to within the value of the constraint tolerance. R=100.0 ohm, L=1.992e-01 H, C=1.424e-05 F, RMSE=0.131.

请先登录，再进行评论。

Arthur FdA
Arthur FdA 2024-10-7
Dear @Sam Chak, @Star Strider, and @William Rose,
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
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 FdA
Arthur FdA 2024-10-8
Oh, okay! All I have is the data I attached before as a .txt document. I'm also attaching it as .xlsx. One of the seven columns is time(s); i'm just not sure if it helps it work.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Linear Model Identification 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by