使用lmfit为Drude-Smith-Anderson拟合模型添加约束遇问题求助

我尝试用lmfit.minimize拟合Drude-Smith-Anderson复电导率模型，需要给参数c和c1施加约束：0<c<1、-1<c1<0且0<1+c1-c<1，编写了如下代码：

#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 path\name.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 σ₀: ")), min =10 , max = 5000) #bounds have to be changed manually
params.add("tau", value = float(input("Guess for τ: ")), 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 τ₁/τ: ")),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到1之间，该如何解决？

问题根源 #

你当前通过expr="1+c1-constraint"定义c，同时给c加了min=0, max=1，但这个设置并没有正确实现0<1+c1-c<1的约束。实际拟合中，当c被压到下限0时，1+c1-c = 1+c1，而c1接近-1，导致结果为负，直接违反约束。

正确的约束实现方式 #

方式1：重新定义变量，自动满足所有约束 #

把0<1+c1-c<1变形为c1 < c < 1+c1，结合已有约束0<c<1、-1<c1<0，可以换一种变量定义逻辑：

• 保留c1的约束：-1 < c1 < 0
• 定义新变量k，范围0 < k < 1，令c = c1 + k

这样既自动满足c1 < c < 1+c1（因为k∈(0,1)），同时结合c1∈(-1,0)，c的范围会自动落在(0,1)，无需额外设置c的上下限。

修改后的参数定义代码：

params = lmf.Parameters()
params.add("sigma0", value = float(input("Guess for σ₀: ")), min =10 , max = 5000)
params.add("tau", value = float(input("Guess for τ: ")), min = 0.0001, max =10)
params.add("c1", value = float(input("Guess for c1: ")), min = -1 , max = 0)
# 定义k∈(0,1)，通过k实现c = c1 + k，自动满足所有约束
params.add("k", value = float(input("Guess for k (0<k<1): ")), min=0, max=1)
params.add("c", expr="c1 + k")
params.add("d", value = float(input("Guess for τ₁/τ: ")),min = 100, max = 100000)

方式2：给残差函数添加约束惩罚（适用于复杂约束场景） #

如果需要保留原有变量逻辑，可以在残差函数中加入约束惩罚项，当约束被违反时，给残差添加一个大的惩罚值，迫使拟合过程遵守约束：

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))
    
    # 添加约束惩罚：当1+c1-c不在(0,1)时，追加惩罚值
    c = params["c"].value
    c1 = params["c1"].value
    constraint_val = 1 + c1 - c
    if constraint_val <=0 or constraint_val >=1:
        # 惩罚系数可根据数据规模调整，确保足够约束拟合方向
        residual += 1000 * abs(constraint_val - 0.5)
    
    return residual

方式3：调整微分进化算法参数 #

你使用的differential_evolution方法中，tol=0.01可能过于宽松，导致算法提前停止在不满足约束的点。可以尝试：

• 减小tol值（比如tol=1e-4）
• 增大popsize（比如popsize=100）
• 关闭polish选项（polish=False），避免抛光阶段违反边界约束

约束验证 #

拟合完成后，可手动计算约束值验证是否符合要求：

c = minimizer_results.params['c'].value
c1 = minimizer_results.params['c1'].value
print(f"1+c1-c = {1 + c1 - c}")

内容的提问来源于stack exchange，提问作者Garvit Bansal