Fitting data to an equation with complex part
显示 更早的评论
Hi,
I hope you are all well.
I have experimental data for a diaphragm displacement and want to fit it to the following equation. The experimental data is displacement (H in meters) and frequency (w in rad/s).
I have a working code but it is changes a lot with the initial guess. I need help to improve this.
x(1) is the alpha (numerator) and x(2) is the damping. omega_r is the resonant frequency and it is known. Gamma should be between 0.04 and 0.07. I am interested in gamma, the damping term in the denominator.
% non linear fitting
fun = @(x,omega)(x(1)./(omega_r^2 - omega.^2 + 1i*(2*x(2)*omega_r.*omega)));
x0 = [-35.9811 + 1i*23.8154,0.06];
opts = optimoptions(@lsqcurvefit,'Display','off','Algorithm','trust-region-reflective');
[vestimated,resnorm] = lsqcurvefit(fun,x0,omega,H,[],[],opts);
Looking forward to your suggestions.
Best regards,
Baris
Edit: x0 is changed.
13 个评论
Alex Sha
2020-4-10
Post your data please
Baris Gungordu
2020-4-10
darova
2020-4-10
Can you please collaborate?
Baris Gungordu
2020-4-10
J. Alex Lee
2020-4-10
Before thinking about the optimization code, do you understand your problem well? For example you can analyze how the norm of your model-to-data residual changes with your parameters - luckily you only have 2 parameters so you can visualize it very easily in a contour plot, c.f.
produced by code below. It doesn't appear to me that gamma should be bounded by 
based on this plot...
based on this plot...data = readmatrix(datapath);
omega = data(:,1);
H = data(:,2);
omega_r = 10920;
fun = @(x,omega)(x(1)./(omega_r^2 - omega.^2 + 1i*(2*x(2)*omega_r.*omega)));
x0 = [-1.5780e+04,0.05];
N = 100;
M = 100;
alph = x0(1) + 1e5*linspace(-1,1,N);
gamm = x0(2) + 1e-1*linspace(-1,1,M);
[X1,X2] = ndgrid(alph,gamm);
obj = nan(size(X1));
for k = 1:numel(X1)
obj(k) = localobj([X1(k),X2(k)] , fun, H,omega);
end
figure(1); cla; hold on;
contourf(X1,X2,obj,logspace(log10(min(obj(:))),log10(max(obj(:))),10))
plot(x0(1),x0(2),'o','MarkerFaceColor','w','LineWidth',2)
set(gca,'ColorScale','log')
colorbar
% opts = optimoptions(@lsqcurvefit,'Display','off','Algorithm','trust-region-reflective');
% [vestimated,resnorm] = lsqcurvefit(fun,x0,omega,H,[],[],opts);
function obj = localobj(x , fun, H,omega)
res = localres(x , fun, H,omega);
obj = norm(res);
end
function res = localres(x , fun, H,omega)
Hmdl = fun(x,omega);
res = H - Hmdl;
end
Baris Gungordu
2020-4-10
编辑：Baris Gungordu
2020-4-10
J. Alex Lee
2020-4-10
Hi Baris,
Sorry, I missed that alpha is complex, somehow i thought only the model was complex. You could try the above exercise with the real and imaginary parts of alpha, if you have reason to believe gamma can be fixed to some value, or do it for several values of gamma within your bounds
Alex Sha
2020-4-10
if parameter x1 and x2 are all complex-type values, then there will be reults below:
Root of Mean Square Error (RMSE): 1.17407781380902E-7
Sum of Squared Residual: 4.54891375249931E-12
Correlation Coef. (R): 0.999813998640161
R-Square: 0.999628031876828
Adjusted R-Square: 0.999625756842435
Determination Coef. (DC): 0.999625508902353
F-Statistic: 290981.048034672
Parameter Best Estimate
-------------------- -------------
x1.realpart -386.462851740156
x1.imagpart 416.525052408265
x2.realpart 0.031303394360693
x2.imagpart 0.0114859438407334
Baris Gungordu
2020-4-10
Alex Sha
2020-4-10
Hi, Baris, rather than Matlab, the results above are obtained by using a package called "1stOpt", the code looks like:
ComplexStr = i;
ComplexPar x1,x2;
Constant omega_r = 10920;
Variable omega, H[realPart];
Function H=(x1/(omega_r^2 - omega^2 + 1*i*(2*x2*omega_r*omega)));
Data;
10737 4.6598e-05
10738 4.6863e-05
10740 4.7127e-05
10741 4.739e-05
10743 4.7654e-05
10744 4.7916e-05
10746 4.8178e-05
....
J. Alex Lee
2020-4-10
I am no longer clear on the problem or objective. Baris, what do you mean gamma is a "positive complex number"? If both alpha and gamma are complex, are there any constraints that you need to ensure that H_model remains real? Or as Alex has asserted, are you chopping off the imaginary part of the H_model? Actually, even if gamma is strictly real, the second part of the question remains.
Is the issue that your solution is too sensitive to initial guess?
Baris Gungordu
2020-4-10
Baris Gungordu
2020-4-10
采纳的回答
Alex Sha
2020-4-11
0 个投票
lsqcurvefit is seneitive to initial start-values since it use local optimization algorithm, same as cftool, there is a GA toolbox with global optimization algorithm in Matlab, unfortunately, the effects of GA are not so good as expected.
2 个评论
Baris Gungordu
2020-4-12
J. Alex Lee
2020-4-13
I don't know that it's helpful to try to characterize the methods' sensitivities...it's really sensitivity, i.e., the nature of the optimization problem that is important. If your optimization landscape is bumpy, as Alex Sha says, methods like Levenberg and trust-region may only find local bumps, depending on how much control you have over the size of the search space. On the other hand, if you have a very smooth optimization landscape, you might be fine. On the other hand, if your optimization landscape is smooth but flat, most algorithms will have trouble because no set of parameters will look any better than the next.
This is why I recommended characterizing your specific problem to the extent possible by studying norm(res) (and its gradient) as a function of your fit parameters.
It looks like you have a 4 parameter fit (2 real parts and 2 complex parts). But depending on how the model works (since you only care about real part of H), it's possible that you're over-specifying, in which case your sensitivity might have to do with multiple combinations of your 4 parameters result in the same real part of the model.
更多回答（1 个）
Baris Gungordu
2020-4-17
0 个投票
1 个评论
Alex Sha
2020-4-17
Hi, if all are as previous except the data, then the result will be:
Root of Mean Square Error (RMSE): 0.478025304038979
Sum of Squared Residual: 133.220275528809
Correlation Coef. (R): 0.994221627528188
R-Square: 0.9884766446448
Parameter Best Estimate
-------------------- -------------
x1.realpart 828282884.340485
x1.imagpart -136314338.186314
x2.realpart 0.681537210289169
x2.imagpart 9.3260946979766
类别
在 帮助中心 和 File Exchange 中查找有关 Linear Least Squares 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
