Fitting Impedance data file to parallel RLC circuit model using fmincon

Hello

I am trying to fit an impedance data file from simuation to a parallel equivalent RLC circuit model. attached is files needed to run the fmincon optimaztion command.

Please any help will be appreciated.

[Star Strider]

It would be helpful to know the circuit to be modeled, since it is difficult for me to understand it from the code.

Also, it may not be possible to identify the model because there is simply too much missing information. The frequency vector should be complete from D-C to the Nyquist frequency, and while the maximum of the frequency vector may be the Nyquist frequency (I don’t know that for sure), the information from 0 to 5.6 Hz (assuming unstated frequency units) is missing.

I need complete data and a graphic or functional description of the circuit to be estimated.

I cannot go further without all those.

[Abdullah Qaroot]

1) The circuit I am trying to model is a simple Parallel RLC circuit. Please check the link below.

https://www.electronics-tutorials.ws/accircuits/parallel-resonance.html

2) regarding the frequency vector. I am interested in 5.6 GHz to 6 GHz frequency range. the circuit should resonate at 5.85 GHz. So no need to have the entire frequency range starting from DC.

3) if you run the testplot.m file attached above. you will have the complete plots for the impedance data. Namly figure 3, 4, and 5 represent the real part, Imaginary part, and magnitude of Impedance data. similarly figure 1, 2, and 6 represent the real part, Imaginary part, and magnitude of the fited circuit model.

also please note that x(1) represent the resistance

x(2) represent the capacitance

x(3) represent the inductance

I hope these information are enough, please let me know if you have any other questions.

[Star Strider]

‘So no need to have the entire frequency range starting from DC.’

I disagree. The system cannot be modeled if the data are incomplete.

[Abdullah Qaroot]

The data are complete. you have the admittance data vs. frequency over the band from 5.6 GHz upto 6 GHz.

you can check the attached Y_adm.mat file. these are the measured/sim data. we don't have the measurments over a band statring from DC!!

again we are interested in finding the equivalent circuit model over the frequency band 5.6 GHz to 6 GHz.

[Star Strider]

The parameter estimation or optimisation will not work with incomplete data.

I will leave you to explore this at your leisure.

Good luck, and have fun!

Answer 1

[Mathieu NOE]

hello

my attempt below - with just fminsearch (as I don't have the Optimisation Toolbox)

I refined a bit the IC by manual serach , not even sure it was needed. ...

clear all
clc

load('Y_adm.mat')

figure(1)
freq = freq*1e9; % must be GHz
Zp = 1./Yp1(:); % measured 
R0 = 100;
L0 = 35e-012;
C0 = 1/(L0*(2*pi*5.85e9)^2);

x = [R0 C0 L0]; %from ADS
Zth = 1./((1./x(1)) + (1i.*2.*pi.*freq.*x(2)) - (1i./(2.*pi.*freq.*x(3))));

semilogy(freq,abs(Zp),'b',freq,abs(Zth),'r')
grid on

% curve fit using fminsearch
x = freq;
y = abs(Zp);

f = @(a,b,c,x) abs(1./((1./a) + (1i.*2.*pi.*freq.*b) - (1i./(2.*pi.*freq.*c))));

obj_fun = @(params) norm(f(params(1), params(2), params(3),x)-y);

sol = fminsearch(obj_fun, [R0,C0,L0]);
a_sol = sol(1);
b_sol = sol(2);
c_sol = sol(3);
y_fit = f(a_sol, b_sol,c_sol, x);
Rsquared = my_Rsquared_coeff(y,y_fit); % correlation coefficient

figure(2)
plot(x, y_fit, '-',x,y, 'r .', 'MarkerSize', 20)
legend('fit','data');
title(['Gaussian Fit / R² = ' num2str(Rsquared) ], 'FontSize', 15)
ylabel('Z', 'FontSize', 14)
xlabel('freq (GHz)', 'FontSize', 14)

R = a_sol
C = b_sol
L = c_sol

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Rsquared = my_Rsquared_coeff(data,data_fit)
    % R2 correlation coefficient computation
    
    % The total sum of squares
    sum_of_squares = sum((data-mean(data)).^2);
    
    % The sum of squares of residuals, also called the residual sum of squares:
    sum_of_squares_of_residuals = sum((data-data_fit).^2);
    
    % definition of the coefficient of correlation is
    Rsquared = 1 - sum_of_squares_of_residuals/sum_of_squares;
    
    
end

[Mathieu NOE]

and optimal values are

R = 71.4678

C = 2.6527e-11

L = 2.7941e-11

[Abdullah Qaroot]

Thank you very much Mathieu.

[Mathieu NOE]

My pleasure !

[Abdullah Qaroot]

Mathieu,

Please I need your help one more time!

I am trying to find the amplitude coefficient of 4x4 antenna array factor. attached is the code I am using. the problem is I need the maximum to be at theta 30 degree. also I attached the second code to plot the results.

thanks

Abdullah

[Mathieu NOE]

hello Abdullah

well, I am not an expert in beam forming but I can try to help you

do you mean you need to find a 4x4 antenna array factor to get a single lobe with max amplitude at 30° ?

[Mathieu NOE]

I tried a few things in your main code , but I don't get whe changing the AB_phases values does not have any impact on the lobe positions

%3D Array Factor of a 4x4 planar array antenna 
%Elements on the x and y axis

clc
clearvars
close all
       
%constants
c=3e8;
f=5.8e9;
lambda=c/f;
dx=  44e-3; %lambda/2; %distance between elements (X)
dy=  35e-3; %lambda/2;%distance between elements (Y)
AB=[6.136359408801308  -1.566346488973549   0.658011497614044  -5.230748528495800
   3.322000897915952  -6.681141375408192   6.262766180230134  -2.893350643833831
  -6.027064578550604  -0.275987581318047  -3.988065810115025  10.240520530592832
   0.869802182657515   2.084910570393388   2.238944073677493  -5.203033164059653]; % Array amplitudes

% AB=[1 1 1 1;
%     1 1 1 1;
%     1 1 1 1;
%     1 1 1 1].*rand; % Array amplitudes

AB=ones(4,4); % Array amplitudes


% AB_phase=[0 0 0 0;
%           0 0 0 0;
%           0 0 0 0;
%           0 0 0 0]*pi/180; %Array phases
      
AB_phase=eye(4,4)*30*pi/180; % Array phases
      
AB_coe=AB.*(cos(AB_phase)+1i*sin(AB_phase));

total_power= ( sum(sum(AB.^2)) ); 

%AF calculation

theta0= [-pi:pi/100:pi]; % increased resoltion from pi/50 to pi/100;
phi0= 0 ;
[phi,theta]=meshgrid(phi0,theta0);
sinU=sin(theta).*cos(phi);
sinV=sin(theta).*sin(phi);

AF_Field=0;
for n=1:4
    for m=1:4
        AF_Field = AB_coe(n,m)*exp(1i*2*pi*(n-1)/lambda*dx*(sinU)).*exp(1i*2*pi*(m-1)/lambda*dy*(sinV)) + AF_Field; 
    end
end
AF_power=abs(AF_Field.^2); % mod : abs
AF_power_normalized=AF_power/total_power;
 
 fitness = AF_power_normalized;
 [row,col] = find(fitness == max(AF_power_normalized));


if theta(row,col) >= 25*pi/180 & theta(row,col) <= 35*pi/180
Cost= max(fitness) - 16;
else
    Cost= 1000;
end
 
%Plotting

 figure(1)
plot(theta*180/pi,AF_power_normalized,theta(row)*180/pi,AF_power_normalized(row),'dr')
xlim([-180 180])
grid on

xlabel('\theta')
ylabel('AF')

 figure(2)
 
 polarplot(theta,AF_power_normalized,theta(row),AF_power_normalized(row),'dr')

[Abdullah Qaroot]

Hi Mathiue

changing to the phses has an impact please check the attached code. the phase shift should be progressive in dx. in our example here since phi = 0 the phase shift in dy = 0.

[Mathieu NOE]

hello again

ok, I can see the result

but I am a bit lost by the fact that you have some many degrees of freedom (AB_amplitudes, AB_phase, dx, dy,...) and your only target is to get a directivity plot centered at 30°

seems like an overdetermined system to me...