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

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?

Answer 1

[Torsten]

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

[Ryban1]

Thank you for your reply, I agree that the fit needs to include both real and imaginary parts of which what I'm asking help with. I'm showing the separate fits as an example and questioning why the fits resulted in different aswers for which you eluded to not applying sufficient bounds and good enough initial guesses. Do you have a reccomendation on which method to use between "lsqnonlin", "lsqcurvefit", and "fminsearch"? It would be great to see if "fminsearch" could be validated against the fitting methods or vice versa to get the same parameters reported in the paper.

I saw that you gave a simialr response in the link below but didn't see a solution.

https://www.mathworks.com/matlabcentral/answers/346835-how-does-one-specify-a-real-datatype-when-attempting-to-fit-real-parameters-to-a-complex-function-wi

[Torsten]

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]

@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]

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]

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]

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.

[Ryban1]

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.