我觉得你问错了问题。lmfit在这里没有任何优势，所以简单的决定是废弃它，直接使用scipy;但是-你现在的代码被诅咒了，所以应该修复。向量化，分解出常见的表达式，并减少到一个内部的lstsq调用。下面显示旧方法和新方法在数值上等效于~13位有效数字，但新方法将更快，更易于调试。在许多其他原因中，您的旧代码经历了许多内部线性拟合的所有麻烦，但随后总是丢弃结果。

import numpy as np
import scipy.optimize as so

data_1 = np.array([
    [117.417, 117.423, 117.438, 117.501],
    [124.160, 124.231, 124.089, 124.100],
    [115.632, 115.645, 115.828, 115.947],
    [118.314, 118.317, 118.287, 118.228],
    [108.407, 108.419, 108.396, 108.564],
])
data_2 = np.array([
    [9.050, 9.044, 9.057, 9.079],
    [9.178, 9.167, 9.160, 9.176],
    [7.888, 7.893, 7.911, 7.895],
    [7.198, 7.202, 7.197, 7.213],
    [7.983, 7.976, 7.979, 8.020],
])

def setup_for_chi(initial: np.ndarray, io: float) -> tuple[
    np.ndarray,
    np.ndarray,
    np.ndarray,
]:
    k, k1, x, y = initial
    xy1 = np.array((1, x, y, x*y))
    k4 = k*xy1
    k14 = k1*xy1

    partial_free_concentration = (
        np.sqrt(
            (k4*k14)**2
            + 8*io*k4*k14*(1 + k14)
        ) - k4*k14
    )/(4*(1 + k14))

    partial_closed_concentration = (
        4*io + 4*io*k14
    )/(
        4*(1 + k14)**2
    ) - partial_free_concentration/(1 + k14)

    partial_open_concentration = k14 * partial_closed_concentration

    populations = np.stack((
        partial_free_concentration,
        partial_open_concentration,
        partial_closed_concentration,
    ), axis=1)

    return (
        np.broadcast_to(populations/io, (len(data_1), *populations.shape)),
        data_1*80,
        data_2*800,
    )

def fit_populations(initial: np.ndarray, io: float) -> tuple[np.ndarray, float]:
    populations_io, experimental_shifts_n, experimental_shifts_h = setup_for_chi(initial, io)

    # 4x3
    A = populations_io[0, ...]

    # 4x10
    b = np.vstack((
        experimental_shifts_n,
        experimental_shifts_h,
    )).T

    x, chi2, rank, s = np.linalg.lstsq(a=A, b=b, rcond=None)
    return x, chi2.sum()

def get_populations(initial: np.ndarray, io: float) -> float:
    x, chi2 = fit_populations(initial, io)
    return chi2

def get_populations_old(initial,io, debug=False):
    from scipy.sparse.linalg import lsmr
    k,k1,x,y=initial[0],initial[1],initial[2],initial[3]
    kx,k1x=k*x,k1*x
    ky,k1y=k*y,k1*y
    kxy,k1xy=k*x*y,k1*x*y
    partial_free_concentration_WT=(np.sqrt((k*k1)**2+(8*io*k*k1)+(8*io*k*k1**2))-(k*k1))/(4*(1+k1))
    partial_closed_concentration_WT=(((4*io)+(4*io*k1))/(4*(1+k1)**2))-(partial_free_concentration_WT/(1+k1))
    partial_open_concentration_WT=k1*partial_closed_concentration_WT
    partial_free_concentration_L273A=(np.sqrt((kx*k1x)**2+(8*io*kx*k1x)+(8*io*kx*k1x**2))-(kx*k1x))/(4*(1+k1x))
    partial_closed_concentration_L273A=(((4*io)+(4*io*k1x))/(4*(1+k1x)**2))-(partial_free_concentration_L273A/(1+k1x))
    partial_open_concentration_L273A=k1x*partial_closed_concentration_L273A
    partial_free_concentration_I272A=(np.sqrt((ky*k1y)**2+(8*io*ky*k1y)+(8*io*ky*k1y**2))-(ky*k1y))/(4*(1+k1y))
    partial_closed_concentration_I272A=(((4*io)+(4*io*k1y))/(4*(1+k1y)**2))-(partial_free_concentration_I272A/(1+k1y))
    partial_open_concentration_I272A=k1y*partial_closed_concentration_I272A
    partial_free_concentration_ILAA=(np.sqrt((kxy*k1xy)**2+(8*io*kxy*k1xy)+(8*io*kxy*k1xy**2))-(kxy*k1xy))/(4*(1+k1xy))
    partial_closed_concentration_ILAA=(((4*io)+(4*io*k1xy))/(4*(1+k1xy)**2))-(partial_free_concentration_ILAA/(1+k1xy))
    partial_open_concentration_ILAA=k1xy*partial_closed_concentration_ILAA
    local_chi2=0
    for experimental_shifts_n,experimental_shifts_h in zip(data_1,data_2):
        populations=np.array([[partial_free_concentration_WT,partial_open_concentration_WT,partial_closed_concentration_WT],[partial_free_concentration_L273A,partial_open_concentration_L273A,partial_closed_concentration_L273A],[partial_free_concentration_I272A,partial_open_concentration_I272A,partial_closed_concentration_I272A],[partial_free_concentration_ILAA,partial_open_concentration_ILAA,partial_closed_concentration_ILAA]])
        experimental_shifts_n=(np.array([experimental_shifts_n])/10)*800
        experimental_shifts_h=(np.array([experimental_shifts_h]))*800
        least_squared_fit_n=lsmr(populations/io,experimental_shifts_n,maxiter=10)
        least_squared_fit_h=lsmr(populations/io,experimental_shifts_h,maxiter=10)
        if debug:
            print(least_squared_fit_n[0])
            print(least_squared_fit_h[0])
        local_chi2+=((least_squared_fit_n[3])**2+least_squared_fit_h[3]**2)
    return local_chi2

io = 270_000

def new_optimize() -> None:
    print('New:')
    solution = so.minimize(
        fun=get_populations,
        args=(io,),
        x0=np.array([500, 0.02, 7, 30]),
        bounds=((0, np.inf),)*4,
        options={'maxiter': 1_000},
    )
    assert solution.success, solution.message
    print('Initial:', solution.x)
    x, chi2 = fit_populations(solution.x, io)
    print('Residual:', chi2)

def regression_test() -> None:
    popio, exp_n, exp_h = setup_for_chi(initial=np.array([500, 0.02, 7, 30]), io=io)

    assert len(popio) == 5

    for pop in popio:
        assert np.allclose(
            pop,
            np.array([
                [0.00425185, 0.01952447, 0.97622368],
                [0.02781779, 0.11939080, 0.85279142],
                [0.09698655, 0.33863005, 0.56438341],
                [0.32547629, 0.54480761, 0.1297161]],
            ),
        )

    assert np.allclose(
        np.vstack(exp_n),
        np.array([
            [9393.36, 9393.84, 9395.04, 9400.08],
            [9932.8 , 9938.48, 9927.12, 9928.  ],
            [9250.56, 9251.6 , 9266.24, 9275.76],
            [9465.12, 9465.36, 9462.96, 9458.24],
            [8672.56, 8673.52, 8671.68, 8685.12],
        ]),
    )
    assert np.allclose(
        np.vstack(exp_h),
        np.array([
            [7240. , 7235.2, 7245.6, 7263.2],
            [7342.4, 7333.6, 7328. , 7340.8],
            [6310.4, 6314.4, 6328.8, 6316. ],
            [5758.4, 5761.6, 5757.6, 5770.4],
            [6386.4, 6380.8, 6383.2, 6416. ],
        ]),
    )

def old_optimize() -> None:
    scipy_solution = so.minimize(get_populations_old, args=(io,), x0=np.array([500, 0.02, 7, 30]),
                                 bounds=((0, np.inf),) * 4, method='Nelder-Mead',
                                 options={'maxiter': 1000})
    assert scipy_solution.success, scipy_solution.message
    print('OLD -')
    print('Initial:', scipy_solution.x)
    print('Residual:', scipy_solution.fun)
    print('Inner factor:')
    get_populations_old(scipy_solution.x, io, debug=True)
    x, chi2 = fit_populations(scipy_solution.x, io)
    print('chi2 from new, to compare:', chi2)
    print('Inner factor from new:')
    print(x)
    print()

if __name__ == '__main__':
    regression_test()
    old_optimize()
    new_optimize()

个字符