使用lmfit向我的拟合模型添加约束

Question

我正在尝试使用lmfit.minimize来拟合复杂的电导率模型(德鲁德-史密斯-安德森模型)。在拟合中，我希望对参数c和c1进行约束，使得0

#reference: Juluri B.K. "Fitting Complex Metal Dielectric Functions with Differential Evolution Method". http://juluribk.com/?p=1597.
#reference: https://lmfit.github.io/lmfit-py/fitting.html
#import libraries (numdifftools needs to be installed but doesn't need to be imported)
import matplotlib.pyplot as plt
import numpy as np
import lmfit as lmf
import math as mt
#define the complex conductivity model
def model(params,w):
sigma0 = params["sigma0"].value
tau = params["tau"].value
c = params["c"].value
d = params["d"].value
c1 = params["c1"].value
druidanderson = (sigma0/(1-1j*2*mt.pi*w*tau))*(1 + c1/(1-1j*2*mt.pi*w*tau)) - sigma0*c/(1-1j*2*mt.pi*w*d*tau)
return druidanderson
#defining the complex residues (chi squared is sum of squares of residues)
def complex_residuals(params,w,exp_data):
delta = model(params,w)
residual = (abs((delta.real - exp_data.real) / exp_data.real) + abs(
(delta.imag - exp_data.imag) / exp_data.imag))
return residual
# importing data from CSV file
importpath = input("Path of CSV file: ") #Asking the location of where your data file is kept (give input in form of pathname.csv)
frequency = np.genfromtxt(rf"{importpath}",delimiter=",", usecols=(0)) #path to be changed to the file from which data is taken
conductivity = np.genfromtxt(rf"{importpath}",delimiter=",", usecols=(1)) + 1j*np.genfromtxt(rf"{importpath}",delimiter=",", usecols=(2)) #path to be changed to the file from which data is taken
frequency = frequency[np.logical_not(np.isnan(frequency))]
conductivity = conductivity[np.logical_not(np.isnan(conductivity))]
w_for_fit = frequency
eps_for_fit = conductivity
#defining the bounds and initial guesses for the fitting parameters
params = lmf.Parameters()
params.add("sigma0", value = float(input("Guess for u03C3u2080: ")), min =10 , max = 5000) #bounds have to be changed manually
params.add("tau", value = float(input("Guess for u03C4: ")), min = 0.0001, max =10) #bounds have to be changed manually
params.add("c1", value = float(input("Guess for c1: ")), min = -1 , max = 0) #bounds have to be changed manually
params.add("constraint", value = float(input("Guess for constraint: ")), min = 0, max=1)
params.add("c", expr="1+c1-constraint", min = 0, max = 1) #bounds have to be changed manually
params.add("d", value = float(input("Guess for u03C4_1/u03C4: ")),min = 100, max = 100000) #bounds have to be changed manually

# minimizing the chi square
minimizer_results = lmf.minimize(complex_residuals, params, args=(w_for_fit, eps_for_fit), method = 'differential_evolution', strategy='best1bin',
popsize=50, tol=0.01, mutation=(0, 1), recombination=0.9, seed=None, callback=None, disp=True, polish=True, init='latinhypercube')
lmf.printfuncs.report_fit(minimizer_results, show_correl=False)

作为匹配的结果，我得到以下输出:

sigma0:      3489.38961 (init = 1000)
tau:         1.2456e-04 (init = 0.01)
c1:         -0.99816132 (init = -1)
constraint:  0.98138820 (init = 1)
c:           0.00000000 == '1+c1-constraint'
d:           7333.82306 (init = 1000)

这些值没有任何意义，因为1+c1-c = -0.97954952不是0，因此无效。如何解决这个问题?

Answer 1

您的代码不可运行。input()的使用有点令人震惊-请不要这样做。编写易于阅读的代码，并将i/o与逻辑分离。

使用complex_array.view(float)

从复杂数组中提取浮点余数猜测任何参数值处于或非常接近其极限(这里是c)是一个非常糟糕的主意，可能使拟合更加困难。

对于您的问题，您将c定义为"评估1+c1-constant，然后应用min=0, max=1"的边界。这就是你的

params.add("c", expr="1+c1-constraint", min = 0, max = 1)

表示:将c计算为1+c1-constraint，然后应用边界[0,1]。代码完全按照你的要求去做。

除非你知道你在做什么(我怀疑可能不是;))，我强烈建议在尝试使用differential_evolution之前使用默认的leastsq方法。事实证明，differential_evolution不是一个很好的全局拟合方法(shgo通常更好，尽管没有"全局")。求解器应该被认为是非常可靠的)。但是，除非您知道您需要这样的方法，否则您可能不会这样做。

我还强烈建议你绘制你的数据和一些模型，用你认为合理的参数进行评估。