Octave中使用leasqr拟合复数EIS数据报错求助

问题背景 #

使用Windows 10下的Octave 7.2.0，借助Optim包的leasqr算法（Levenberg-Marquardt非线性回归）拟合复数形式的EIS数据（实部ReZ、虚部ImZ）时，直接传入复数ydata = ReZ + j*ImZ会触发错误：

error: weighted residuals are not real
error: called from
    __lm_svd__ at line 147 column 20
    leasqr at line 662 column 25
    Code_for_StackOverflow at line 47 column 73

分开拟合实部与虚部虽能成功，但会得到两组独立参数，不符合单组参数拟合的需求。

错误原因 #

1. leasqr的底层实现__lm_svd__要求加权残差必须为实数，直接传入复数数据时，残差的虚部无法被处理，导致报错。
2. 原代码中权重wt = abs(sqrt(ydata).^-1)虽用abs()取模，但ydata是复数，计算过程中仍可能引入不必要的复数处理逻辑，加重问题。

解决方案 #

将复数拟合转化为实数域的联合拟合：

• 把复数数据拆分为实部和虚部，拼接成一个一维实数向量作为新的ydata
• 修改拟合函数，使其输出对应长度的实数向量（先输出所有频率点的拟合实部，再输出所有频率点的拟合虚部）
• 使用实数权重，避免复数干扰

修改后完整代码 #

clear -a;
clf;
clc;
pkg load optim;
pkg load symbolic;

# 原始实验数据
Linear_freq = [1051.432, 394.2871, 112.6535, 39.42871, 11.59668, 3.458659, 1.065641, 0.3258571, 0.1000221];
ReZ = [84.10412, 102.0962, 178.8031, 283.0663, 366.7088, 431.3653, 514.4105, 650.5853, 895.9588];
MinusImZ = [27.84804, 59.56786, 116.5972, 123.2293, 102.6806, 117.4836, 178.1147, 306.256, 551.2337];
Z = [88.5946744, 118.2030626, 213.4606653, 308.7264008, 380.8131426, 447.0776424, 544.3739605, 719.0646495, 1051.950932];
MinusPhase = [18.32042302, 30.26135402, 33.1083029, 23.52528583, 15.64255593, 15.23515301, 19.09841797, 25.2082044, 31.60167787];
ImZ = -MinusImZ;
Angular_freq = 2*pi*Linear_freq;
xdata = Angular_freq;

# 重构ydata：先实部，后虚部，转为实数向量
ydata = [ReZ, ImZ];

# 修改拟合函数：输出实数向量，顺序为所有点的实部、所有点的虚部
Fitting_Function = @(xdata, p) ...
  (
    Z_fit = p(1) + ((p(2) + (1./(p(3)*(j*xdata).^0.5))).^-1 + (1./(p(4)*(j*xdata).^p(5))).^-1).^-1;
    [real(Z_fit), imag(Z_fit)]
  );

# 初始参数与设置
p = [80, 300, 6.63E-3, 5E-5, 0.8]; # 参考真实值：76, 283, 1.63E-3, 1.5E-5, 0.876
options.fract_prec = [0.0005, 0.0005, 0.0005, 0.0005, 0.0005].';
niter=400;
tol=1E-12;
dFdp="dfdp";
dp=1E-9*ones(size(p));
# 改用实数权重：基于阻抗模的倒数平方根
wt = 1./sqrt(Z);
# 扩展权重到和ydata同长度（实部和虚部用相同权重）
wt = [wt, wt];

# 执行联合拟合
[Fitted_Parameters pfit cvg iter corp covp covr stdresid z r2] = leasqr(xdata, ydata, p, Fitting_Function, tol, niter, wt, dp, dFdp, options);

#########################################################################
# 用拟合参数计算结果
#########################################################################
Z_fit_full = pfit(1) + ((pfit(2) + (1./(pfit(3)*(j*xdata).^0.5))).^-1 + (1./(pfit(4)*(j*xdata).^pfit(5))).^-1).^-1;
Fitted_Function_Real = real(Z_fit_full);
Fitted_Function_Imag = imag(Z_fit_full);
Fitted_Function_Mod = abs(Z_fit_full);
Fitted_Function_Phase = (-(angle(Z_fit_full))*(180./pi));

################################################################################
# 计算残差（参考IOP文献方法）
################################################################################
Residuals_Real = (ReZ-Fitted_Function_Real)./Fitted_Function_Mod;
Residuals_Imag = (ImZ-Fitted_Function_Imag)./Fitted_Function_Mod;

################################################################################
# 计算卡方值
################################################################################
chi_squared_ReZ = sum(((ReZ-Fitted_Function_Real).^2)./Z.^2);
chi_squared_ImZ = sum(((ImZ-Fitted_Function_Imag).^2)./Z.^2);
Pseudo_chi_squared = sum((((ReZ-Fitted_Function_Real).^2)+((ImZ-Fitted_Function_Imag).^2))./Z.^2);

# 输出结果
disp('拟合得到的参数：'), disp(pfit);
disp('实部卡方值：'), disp(chi_squared_ReZ);
disp('虚部卡方值：'), disp(chi_squared_ImZ);
disp('联合卡方值：'), disp(Pseudo_chi_squared);
disp('决定系数R^2：'), disp(r2);

###################################################
## 绘图部分（与原代码一致）
###################################################
#Set plot parameters
set(0, "defaultlinelinewidth", 1);
set(0, "defaulttextfontname", "Verdana");
set(0, "defaulttextfontsize", 20);
set(0, "DefaultAxesFontName", "Verdana");
set(0, 'DefaultAxesFontSize', 12);

figure(1);
## Nyquist plot (Argand diagram)
subplot(1,2,1, "align");
plot((ReZ), (MinusImZ), "o", "markersize", 2, (Fitted_Function_Real), -(Fitted_Function_Imag), "-k");
axis ("square");
grid on;
daspect([1 1 2]);
title ('Nyquist Plot - Argand Diagram');
xlabel ('Z'' / \Omega' , 'interpreter', 'tex');
ylabel ('-Z'''' / \Omega', 'interpreter', 'tex');

## Bode Modulus
subplot (2, 2, 2);
loglog((Linear_freq), (Z), "o", "markersize", 2, (Linear_freq), (Fitted_Function_Mod), "-k");
grid on;
title ('Bode Plot - Modulus');
xlabel ('\nu (Hz)' , 'interpreter', 'tex');
ylabel ('|Z| / \Omega', 'interpreter', 'tex');

## Bode Phase
subplot (2, 2, 4);
semilogx((Linear_freq), (MinusPhase), "o", "markersize", 2, (Linear_freq), (Fitted_Function_Phase), "-k");
set(gca,'YTick',0:10:90);
grid on;
title ('Bode Plot - Phase');
xlabel ('\nu (Hz)' , 'interpreter', 'tex');
ylabel ('-\theta (°)', 'interpreter', 'tex');

figure(2)
## Bode Z'
subplot (2, 1, 1);
semilogx((Linear_freq), (ReZ), "o", "markersize", 2, (Linear_freq), (Fitted_Function_Real), "-k");
grid on;
title ('Bode Plot Z''');
xlabel ('\nu (Hz)' , 'interpreter', 'tex');
ylabel ('Z'' / \Omega', 'interpreter', 'tex');

## Bode -Z''
subplot (2, 1, 2);
semilogx((Linear_freq), (MinusImZ), "o", "markersize", 2, (Linear_freq), -(Fitted_Function_Imag), "-k");
grid on;
title ('Bode Plot -Z''''');
xlabel ('\nu (Hz)' , 'interpreter', 'tex');
ylabel ('-Z'''' / \Omega', 'interpreter', 'tex');

figure(3)
## Residuals Real
subplot (2, 1, 1);
semilogx((Angular_freq), (Residuals_Real), "-o", "markersize", 2);
grid on;
title ('Residuals Real');
xlabel ('\omega (Hz)' , 'interpreter', 'tex');
ylabel ('\Delta_{re} / \Omega', 'interpreter', 'tex');

## Residuals Imaginary
subplot (2, 1, 2);
semilogx((Angular_freq), (Residuals_Imag), "-o", "markersize", 2);
grid on;
title ('Residuals Imaginary');
xlabel ('\omega (Hz)' , 'interpreter', 'tex');
ylabel ('\Delta_{im} / \Omega', 'interpreter', 'tex');

关键修改点说明 #

1. ydata重构：将ReZ和ImZ拼接成一个一维实数向量，让leasqr在实数域处理所有数据。
2. 拟合函数修改：输出向量包含所有频率点的拟合实部和虚部，保证与ydata维度匹配，用同一组参数同时拟合两类数据。
3. 权重调整：将权重改为实数（基于阻抗模的倒数平方根），并扩展到与ydata同长度，避免复数权重引发的错误。

内容的提问来源于stack exchange，提问作者Corrado Locati