SciPy 的 scipy.optimize.minimize 不遵守约束

Question

我试图根据一些约束，使用SciPy的optimize.minimize来拟合一些分散的数据的抛物线：曲线下的面积以及曲线穿过数据的初始点和终点。这里

curve_fit

工作得很好，但我也希望在其他一些情况下使用带有约束的

minimize

。约束版本似乎没有成功（可能是因为

x

或

y

...？的值范围）。知道如何帮助它收敛吗？

import scipy
import numpy as np
import matplotlib.pyplot as plt

x = np.array([258.6669469, 258.831, 259.129, 259.428 , 259.726, 260.025, 260.324, 260.622,261.15238824])
y = np.array([0, -0.062, -0.139, -0.181, -0.197, -0.193, -0.17 , -0.129, 0])

def Func(x,a,b,c):
    return a*x**2 + b*x + c

def FuncCons(x,params):
    return params[2]*x**2 + params[1]*x + params[0]

def ConstraintIntegral(params):
    integral = scipy.integrate.quad(FuncCons, x[0], x[-1], args=(params,))[0]
    return integral- np.trapz(y,x)

def ConstraintBegin(params):
    return y[0] - FuncCons(x[0],params)

def ConstraintEnd(params):
    return y[-1] - FuncCons(x[-1],params)

def Objective(params,x,y):
    y_pred = FuncCons(x,params)
    return np.sum((y_pred - y) ** 2) # least squares

cons = [{'type':'eq', 'fun': ConstraintIntegral},
        {'type':'eq', 'fun': ConstraintBegin},
        {'type':'eq', 'fun': ConstraintEnd}]

sigma = np.ones(len(x))
sigma[[0, -1]] = 0.0001 # first and last points

popt1, _ = scipy.optimize.curve_fit(Func, x, y, sigma=sigma)

new = scipy.optimize.minimize(Objective, x0=popt1, args=(x,y), constraints=cons)
popt3 = new.x

y_fit1 = Func(x, *popt1)  
y_fit3 = FuncCons(x, popt3) 

fig, ax = plt.subplots(1)
ax.scatter(x,y)
ax.plot(x,y_fit1, color='g', alpha=0.75, label='curve_fit')
ax.plot(x,y_fit3, color='r', alpha=0.75, label='generally constrained')
plt.legend()

Answer 1

这不应该使用迭代拟合；问题是精确确定的，应该进行分析计算：

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit


def quadratic(x: np.ndarray, a: float, b: float, c: float) -> np.ndarray:
    return x*(a*x + b) + c


def estimate(x: np.ndarray, y: np.ndarray, p: float, q: float) -> np.ndarray:
    iv = y.argmin()
    vx = x[iv]
    vy = y[iv]

    return np.linalg.solve(
        a=[
            [  p*p,  p, 1],  # p root
            [  q*q,  q, 1],  # q root
            [vx*vx, vx, 1],  # vertex
        ],
        b=[0, 0, vy],
    )


def analytic_integral(a: float, b: float, c: float, p: float, q: float) -> float:
    return a*(q**3 - p**3)/3 + b*(q**2 - p**2)/2 + c*(q - p)


def solve(p: float, q: float, integ: float) -> tuple[float, float, float]:
    b = integ/(
        -1/(q + p)*(q**3 - p**3)/3
        + (q**2 - p**2)/2
        - p*q/(q + p)*(q - p)
    )
    a = -b/(p + q)
    c = -b*p*q/(p + q)
    return a, b, c


def main() -> None:
    x = np.array((258.6669469, 258.831, 259.129, 259.428, 259.726,
                  260.025, 260.324, 260.622, 261.15238824))
    y = np.array((0, -0.062, -0.139, -0.181, -0.197, -0.193, -0.17, -0.129, 0))
    integ = np.trapz(y, x)
    p = x[0]
    q = x[-1]

    est = estimate(x, y, p, q)
    popt_exact = solve(p, q, integ)
    popt_approx, _ = curve_fit(f=quadratic, xdata=x, ydata=y, p0=est)

    print('     Estimate:', est)
    print('Unconstrained:', popt_approx)
    print('  Constrained:', np.array(popt_exact))
    print('Empirical integral:', integ)
    print(' Analytic integral:', analytic_integral(*popt_exact, p, q))
    print('Root 1:', quadratic(p, *popt_exact))
    print('Root 2:', quadratic(q, *popt_exact))

    xlong = np.linspace(p, q, 1000)

    fig, ax = plt.subplots()
    ax.scatter(x, y, label='Empirical data')
    ax.plot(xlong, quadratic(xlong, *est), label='Estimate')
    ax.plot(xlong, quadratic(xlong, *popt_exact), label='Constrained fit')
    ax.plot(xlong, quadratic(xlong, *popt_approx), label='Unconstrained fit')
    ax.legend()
    plt.show()


if __name__ == '__main__':
    main()

     Estimate: [ 1.30409954e-01 -6.77896154e+01  8.80938681e+03]
Unconstrained: [ 1.28055584e-01 -6.65592985e+01  8.64866088e+03]
  Constrained: [ 1.29167467e-01 -6.71437466e+01  8.72545495e+03]
Empirical integral: -0.3305311875799978
 Analytic integral: -0.330531187588349
Root 1: 0.0
Root 2: 0.0