Find real valued parameters of a complex equation/model by using optimization or curve fitting.

Question

Ryban1 2025-1-8

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2172750-find-real-valued-parameters-of-a-complex-equation-model-by-using-optimization-or-curve-fitting

编辑： Torsten 2025-1-10

Curve_Fit_complex_model.m

I am trying to replicate finding the optimized parameters of a Lorentz model defined in the paper titled, "EXTRACTION OF EFFECTIVE METAMATERIAL PARAMETERS BY PARAMETER FITTING OF DISPERSIVE MODEL" (linked here). I've provided the equation and table of paramters below:

I've tried the curve fitting tool and get different results especially when I try to change the bounds. Sometimes I am able to get paramters that match the table above but my confidence in repeatbility isn't high. I've also tried scripting code to use the "fittype", "lsqnonlin", or "fminsearch" but also get different answers as well. In the code below, I'm comparing the fit to the real and imaginary parts of the Lorentz model. The fit to the real data is way off and even though the fit to the imaginary values looks qualitatively good. For the fit to the imaginary part, I get different values for the estimated parameters compared to what's reported in the paper.

%% Parameters in Table
e_inf = 1.62;
wp = 2*pi*14.63*1e9; % GHz
vc = 30.69*1e6; %MHz
mu_s = 1.26;
mu_inf = 1.12;
wo = 2*pi*9.67*1e9; % GHz
delta = 1.24*1e9; % GHz
%%
fo = [7:0.01:12].*1e9; % GHz
w = 2*pi*fo; 
mu_eff = mu_inf + ((mu_s - mu_inf)*wo.^2)./(wo.^2 + 1i*w*delta - w.^2); % equation in paper / Lorentz model
real_mu_eff = real(mu_eff);
imag_mu_eff = imag(mu_eff);
x = fo;
y = real_mu_eff; 
x2 = fo;
y2 = imag_mu_eff; 
myfittype = fittype("real(a+((b-a)*(2*pi*1e9*c).^2)./((2*pi*1e9*c).^2+1i*2*pi*x*d*1e9-x.^2))",...
    dependent="y",independent="x",...
    coefficients=["a" "b" "c" "d"])
myfit = fit(x',y',myfittype)
figure
plot(myfit,x,y)
myfittype2 = fittype("imag(a2+((b2-a2)*(2*pi*1e9*c2).^2)./((2*pi*1e9*c2).^2+1i*2*pi*x2*d2*1e9-x2.^2))",...
    dependent="y2",independent="x2",...
    coefficients=["a2" "b2" "c2" "d2"])
myfit2 = fit(x2',y2',myfittype2)
figure
plot(myfit2,x2,y2)

The results I was expecting were: a=1.12, b=1.26, c=9.67, and d=1.24. I 've atttached my attempts with using "lsqnonlin", "lsqcurvefit", and "fminsearch"; however I didn't have success in using those methods either. What is best way to find the parameters of the Lorentz model to get the values in the table above?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Torsten 2025-1-8

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2172750-find-real-valued-parameters-of-a-complex-equation-model-by-using-optimization-or-curve-fitting#answer_1557089

编辑：Torsten 2025-1-8

I don't know if this helps to get better results, but of course you have to fit real and imaginary part of mu_eff simultaneously, not with two separate calls to "fit".

And I think it's impossible to get a good result for the 4 parameters without supplying bounds and good initial guesses. Both can be supplied via the "fitoptions" ("Lower", "Upper","StartPoint").

5 个评论
显示 3更早的评论隐藏 3更早的评论

Torsten 2025-1-10

编辑：Torsten 2025-1-10

在 MATLAB Online 中打开

I cancelled a factor of 1e18 in numerator and denominator.

For lsqnonlin, the region of attraction to the "correct" parameters is not that large. Maybe "fit" or "fitnlm" behave better.

%% Parameters in Table
mu_inf = 1.12;
mu_s = 1.26;
wo = 2*pi*9.67; % GHz
delta = 1.24; % GHz
%%
fo = 7:0.01:12; % GHz
w = 2*pi*fo; 
mu_eff = mu_inf + ((mu_s - mu_inf)*wo^2)./(wo^2 + 1i*w*delta - w.^2); % equation in paper / Lorentz model
real_mu_eff = real(mu_eff);
imag_mu_eff = imag(mu_eff);
x = fo;
y = mu_eff; 
%p = lsqnonlin(@(p)fun(p,x,y),[1.1, 1.4, 9, 2],[],[],[],[],[],[],[],optimset('MaxFunEvals',10000,'MaxIter',10000))
p = lsqnonlin(@(p)fun(p,x,y),[1.1, 1.4, 9, 2])
Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.
p = 1×4
    1.1200    1.2600    9.6700    1.2400
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
function res = fun(p,x,y)
  mu_inf = p(1);
  mu_s = p(2);
  wo = 2*pi*p(3);
  delta = p(4);
  w = 2*pi*x;
  mu_eff = mu_inf + ((mu_s - mu_inf)*wo^2)./(wo^2 + 1i*w*delta - w.^2);
  res = [real(mu_eff-y),imag(mu_eff-y)];
end

Gayathri 2025-1-10

编辑：Gayathri 2025-1-10

在 MATLAB Online 中打开

@Ryban1, @Torsten "fminsearch" also gives the same output as "lsqonlin" function.

% Parameters in Table
mu_inf = 1.12;
mu_s = 1.26;
wo = 2*pi*9.67; % GHz
delta = 1.24; % GHz
% Frequency range
fo = 7:0.01:12; % GHz
w = 2*pi*fo;
% Calculate mu_eff using the given parameters
mu_eff = mu_inf + ((mu_s - mu_inf)*wo^2)./(wo^2 + 1i*w*delta - w.^2);
% Lorentz model function
lorentzModel = @(p, w) p(1) + ((p(2) - p(1)) * (2*pi*p(3))^2) ./ ...
    ((2*pi*p(3))^2 + 1i * w * p(4) - w.^2);
% Error function for fminsearch
errorFunction = @(p) sum(abs(lorentzModel(p, w) - mu_eff).^2);
% Initial guesses for the parameters
initialGuesses = [1.1, 1.4, 9, 2];
% Perform optimization using fminsearch
fittedParams = fminsearch(errorFunction, initialGuesses);
% Display the fitted parameters
disp('Fitted parameters using fminsearch:');
Fitted parameters using fminsearch:
disp(fittedParams);
    1.1200    1.2600    9.6700    1.2401

I had tried with the "GlobalSearch" function as well

% Parameters from the paper
mu_inf = 1.12;
mu_s = 1.26;
wo = 2*pi*9.67; % GHz
delta = 1.24; % GHz
% Frequency range
fo = 7:0.01:12; % GHz
w = 2*pi*fo;
% Calculate mu_eff using the given parameters
mu_eff = mu_inf + ((mu_s - mu_inf)*wo^2)./(wo^2 + 1i*w*delta - w.^2);
% Lorentz model function
lorentzModel = @(p, w) p(1) + ((p(2) - p(1)) * (2*pi*p(3))^2) ./ ...
    ((2*pi*p(3))^2 + 1i * w * p(4) - w.^2);
% Error function for optimization
errorFunction = @(p) sum(abs(lorentzModel(p, w) - mu_eff).^2);
% Initial guesses for the parameters
initialGuesses = [1.1, 1.4, 9, 2];
% Lower and upper bounds for the parameters
lb = [0, 0, 0, 0];
ub = [Inf, Inf, Inf, Inf];
% Set up the optimization problem
problem = createOptimProblem('fmincon',...
    'objective', errorFunction,...
    'x0', initialGuesses,...
    'lb', lb,...
    'ub', ub,...
    'options', optimoptions('fmincon', 'Algorithm', 'sqp'));
% Create a GlobalSearch object
gs = GlobalSearch('Display', 'iter');
% Perform the global search
[globalParams, fval] = run(gs, problem);
 Num Pts                 Best       Current    Threshold        Local        Local                 
Analyzed  F-count        f(x)       Penalty      Penalty         f(x)     exitflag        Procedure
       0       74        1423                                    1423            2    Initial Point
     200     1461   3.677e-09                               3.677e-09            2    Stage 1 Local
     300     1563   3.677e-09          2224        1.074                              Stage 2 Search
     400     1663   3.677e-09     2.178e+09         4.16                              Stage 2 Search
     500     1763   3.677e-09          1455        11.84                              Stage 2 Search
     600     1863   3.677e-09         416.9        30.95                              Stage 2 Search
     700     1963   3.677e-09         959.8        65.25                              Stage 2 Search
     800     2063   3.677e-09     4.668e+10        163.8                              Stage 2 Search
     900     2163   3.677e-09           145        135.2                              Stage 2 Search
     954     2374   3.677e-09         133.5        158.6     3.57e-09            2    Stage 2 Local
     962     2539   3.677e-09         132.8        133.5    3.593e-09            2    Stage 2 Local
     965     2703   3.677e-09         131.6        132.8    3.689e-09            2    Stage 2 Local
     972     2862   3.677e-09         131.3        131.6    3.275e-09            2    Stage 2 Local
     977     3023   3.677e-09         131.3        131.3    3.592e-09            2    Stage 2 Local
     979     3186   3.677e-09         131.1        131.3    3.598e-09            2    Stage 2 Local
     980     3340   3.677e-09         128.5        131.1    3.018e-09            2    Stage 2 Local
     986     3506   3.677e-09         125.8        128.5    3.527e-09            2    Stage 2 Local
    1000     3520   3.677e-09         131.2        125.8                              Stage 2 Search

GlobalSearch stopped because it analyzed all the trial points.

All 10 local solver runs converged with a positive local solver exit flag.
% Display the fitted parameters
disp('Fitted parameters using GlobalSearch:');
Fitted parameters using GlobalSearch:
disp(globalParams);
    1.1200    1.2600    9.6700    1.2400

All of the methods greatly depend on the starting point chosen.

Ryban1 2025-1-10

Thanks you both @Torsten and @Gayathri I appreciate the help.

One last question to @Torsten on your comment about "For lsqnonlin, the region of attraction to the "correct" parameters is not that large. ", how are you defining region of attraction? Is that determined from the bounds of the problem, cancelling out the 1e18 factor, or something else?

Torsten 2025-1-10

编辑：Torsten 2025-1-10

how are you defining region of attraction?

I just tested out a bit on how far the initial guess values for the parameter can be away from the solution for that the method still converges to the correct parameter set.

请先登录，再进行评论。

Answer 2

Gayathri 2025-1-8

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2172750-find-real-valued-parameters-of-a-complex-equation-model-by-using-optimization-or-curve-fitting#answer_1557097

编辑：Torsten 2025-1-8

在 MATLAB Online 中打开

Hi @Ryban1,

If the output of your fitting process is inconsistent, the fitting algorithm might be sensitive to the initial guess of parameters. We can also experiment with different fitting algorithms that suits the data to be fitted. We can use the "fitoptions" function to specify the "StartPoint", "Method", Tolerance etc.

Please refer to the below code which implements the same.

%% Parameters in Table

e_inf = 1.62;

wp = 2*pi*14.63*1e9; % GHz

vc = 30.69*1e6; %MHz

mu_s = 1.26;

mu_inf = 1.12;

wo = 2*pi*9.67*1e9; % GHz

delta = 1.24*1e9; % GHz

%%

fo = [7:0.01:12].*1e9; % GHz

w = 2*pi*fo;

mu_eff = mu_inf + ((mu_s - mu_inf)*wo.^2)./(wo.^2 + 1i*w*delta - w.^2); % equation in paper / Lorentz model

real_mu_eff = real(mu_eff);

imag_mu_eff = imag(mu_eff);

x = fo;

y = real_mu_eff;

x2 = fo;

y2 = imag_mu_eff;

opts = fitoptions('Method', 'NonlinearLeastSquares', ...

'StartPoint', [1.5, 2, -1, 1], ... % Example initial guesses

'MaxIter', 100, ...

'TolFun', 1e-8);

myfittype = fittype("real(a+(((b-a)*(2*pi*1e9*c).^2)./((2*pi*1e9*c).^2+1i*2*pi*x*d*1e9-x.^2)))",...

dependent="y",independent="x",...

coefficients=["a" "b" "c" "d"], ...

options=opts);

myfit = fit(x',y',myfittype)

myfit =

General model: myfit(x) = real(a+(((b-a)*(2*pi*1e9*c).^2)./((2*pi*1e9*c).^2+1i*2*pi*x*d*1e9- x.^2))) Coefficients (with 95% confidence bounds): a = 1.12 (1.12, 1.12) b = 1.26 (1.26, 1.26) c = -1.539 (-1.539, -1.539) d = 0.03141 (0.03141, 0.03141)

figure

plot(myfit,x,y)

opts1 = fitoptions('Method', 'NonlinearLeastSquares', ...

'StartPoint', [0, -1, -2, 1], ... % Example initial guesses

'MaxIter', 100, ...

'TolFun', 1e-8);

myfittype2 = fittype("imag(a2+(((b2-a2)*(2*pi*1e9*c2).^2)./((2*pi*1e9*c2).^2+1i*2*pi*x2*d2*1e9-x2.^2)))",...

dependent="y2",independent="x2",...

coefficients=["a2" "b2" "c2" "d2"], ...

options=opts1);

myfit2 = fit(x2',y2',myfittype2)

myfit2 =

General model: myfit2(x2) = imag(a2+(((b2-a2)*(2*pi*1e9*c2).^2)./((2*pi*1e9*c2).^2+1i*2*pi*x2*d2*1e9- x2.^2))) Coefficients (with 95% confidence bounds): a2 = 20.63 (20.61, 20.64) b2 = 20.77 (20.75, 20.78) c2 = -1.539 (-1.539, -1.539) d2 = 0.03141 (0.03141, 0.03141)

figure

plot(myfit2,x2,y2)

I am able to obtain a consistent output as shown below using this approach.

You can also try experimenting with other parameters within "fitoptions" function for better and accurate curve fitting.

For more information on "fitoptions" function, please refer to the below documentation link.

https://www.mathworks.com/help/releases/R2022b/curvefit/fitoptions.html

Hope you find this information helpful.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Ryban1 2025-1-9

Thank you for your response, I know that the parameters cannot be negative and will investigate upper bounds to further constrain the problem. I would like the start point to be [0,0,0,0] or random. I experimented with this in the curve fitting app and ran into issues and will report back what they were. With the change to generating the fits through scripting, I'm still working out which method is best between "lsqnonlin", "lsqcurvefit", and "fminsearch" as I would like to verify and have confidence in my ability to find the fit for data where the prior would come from a simulation of a different design.

Do you know how to apply the different methods fit to a complex funtion so that way there is a similtaneous fit to the equation "mu_eff"? I only showed the individual fits for demonstration purposes and to understand the discrepancies.

请先登录，再进行评论。

Find real valued parameters of a complex equation/model by using optimization or curve fitting.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

5 个评论
显示 3更早的评论隐藏 3更早的评论

更多回答（1 个）

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Find real valued parameters of a complex equation/model by using optimization or curve fitting.

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

5 个评论 显示 3更早的评论隐藏 3更早的评论

更多回答（1 个）

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

5 个评论
显示 3更早的评论隐藏 3更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论