在本例中，我将使用curve_fit或lmfit;我很快就把它展示给第一个。

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

def func(x, a, b, c, d):
  return a + b * x + c * x ** 2 + d * x ** 3

x = np.array([0.0, 1.0, 2.0, 3.0,  4.0,  5.0])
y = np.array([0.0, 0.8, 0.9, 0.1, -0.8, -1.0])

print(np.polyfit(x, y, 3))

popt, _ = curve_fit(func, x, y)
print(popt)

popt_cons, _ = curve_fit(func, x, y, bounds=([-np.inf, 2, -np.inf, -np.inf], [np.inf, 2.001, np.inf, np.inf]))
print(popt_cons)

xnew = np.linspace(x[0], x[-1], 1000)

plt.plot(x, y, 'bo')
plt.plot(xnew, func(xnew, *popt), 'k-')
plt.plot(xnew, func(xnew, *popt_cons), 'r-')
plt.show()

这将打印：

[ 0.08703704 -0.81349206  1.69312169 -0.03968254]
[-0.03968254  1.69312169 -0.81349206  0.08703704]
[-0.14331349  2.         -0.95913556  0.10494372]

因此，在无约束的情况下，polyfit和curve_fit给予相同的结果（只是顺序不同），在约束的情况下，固定参数为2，如所需。
该图看起来如下：

在lmfit中，您还可以选择是否应该拟合参数，因此您也可以将其设置为所需的值（选中this answer）。

为了完整起见，使用lmfit的解决方案如下所示：

import numpy as np
import matplotlib.pyplot as plt
from lmfit import Model

def func(x, a, b, c, d):
    return a + b * x + c * x ** 2 + d * x ** 3

x = np.array([0.0, 1.0, 2.0, 3.0,  4.0,  5.0])
y = np.array([0.0, 0.8, 0.9, 0.1, -0.8, -1.0])

pmodel = Model(func)
params = pmodel.make_params(a=1, b=2, c=1, d=1)

params['b'].vary = False 

result = pmodel.fit(y, params, x=x)

print(result.fit_report())

xnew = np.linspace(x[0], x[-1], 1000)
ynew = result.eval(x=xnew)

plt.plot(x, y, 'bo')
plt.plot(x, result.best_fit, 'k-')
plt.plot(xnew, ynew, 'r-')
plt.show()

它将打印一份综合报告，包括不确定性、相关性和拟合统计数据，如：

[[Model]]
    Model(func)
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 10
    # data points      = 6
    # variables        = 3
    chi-square         = 0.066
    reduced chi-square = 0.022
    Akaike info crit   = -21.089
    Bayesian info crit = -21.714
[[Variables]]
    a:  -0.14331348 +/- 0.109441 (76.37%) (init= 1)
    b:   2 (fixed)
    c:  -0.95913555 +/- 0.041516 (4.33%) (init= 1)
    d:   0.10494371 +/- 0.008231 (7.84%) (init= 1)
[[Correlations]] (unreported correlations are <  0.100)
    C(c, d)                      = -0.987
    C(a, c)                      = -0.695
    C(a, d)                      =  0.610

并生成

图
请注意，lmfit.Model比curve_fit有许多改进，包括基于函数参数自动命名参数，允许任何参数有边界或简单地固定，而不需要像具有几乎相等的上限和下限这样的废话。关键是lmfit使用具有属性的Parameter对象，而不是拟合变量的普通数组。lmfit还支持数学约束，复合模型（例如，添加或乘以模型），并具有上级的报告。

抱歉复活

但我觉得这个答案不见了

为了拟合多项式，我们求解以下方程组：

a0*x0^n + a1*x0^(n-1) .. + an*x0^0 = y0
a0*x1^n + a1*x1^(n-1) .. + an*x1^0 = y1
                 ...
a0*xm^n + a1*xm^(n-1) .. + an*xm^0 = ym

这是一个形式为V @ a = y的问题
其中“V”是范德蒙矩阵：

[[x0^n  x0^(n-1)  1],
 [x1^n  x1^(n-1)  1],
        ...
 [xm^n  xm^(n-1)  1]]

“y”是保持y值的列向量：

[[y0],
 [y1],
  ...
 [ym]]

..并且“a”是我们正在求解的系数的列向量：

[[a0],
 [a1],
  ...
 [an]]

这个问题可以用线性最小二乘法解决，如下所示：

import numpy as np

x = np.array([0.0, 1.0, 2.0, 3.0,  4.0,  5.0])
y = np.array([0.0, 0.8, 0.9, 0.1, -0.8, -1.0])

deg = 3
V = np.vander(x, deg + 1)
z, *_ = np.linalg.lstsq(V, y, rcond=None)

print(z)
# [ 0.08703704 -0.81349206  1.69312169 -0.03968254]

..产生与polyfit方法相同的解决方案：

z = np.polyfit(x, y, deg)

print(z)
# [ 0.08703704 -0.81349206  1.69312169 -0.03968254]

相反，我们需要一个解决方案，其中a2 = 1

将a2 = 1从答案的开始代入方程组，然后将相应的项从lhs移动到rhs，我们得到：

a0*x0^n + a1*x0^(n-1) + 1*x0^(n-2) .. + an*x0^0 = y0
a0*x1^n + a1*x1^(n-1) + 1*x0^(n-2) .. + an*x1^0 = y1
                 ...
a0*xm^n + a1*xm^(n-1) + 1*x0^(n-2) .. + an*xm^0 = ym

=>

a0*x0^n + a1*x0^(n-1) .. + an*x0^0 = y0 - 1*x0^(n-2)
a0*x1^n + a1*x1^(n-1) .. + an*x1^0 = y1 - 1*x0^(n-2)
                 ...
a0*xm^n + a1*xm^(n-1) .. + an*xm^0 = ym - 1*x0^(n-2)

这对应于如下从范德蒙矩阵中移除列2并从y向量中减去它：

y_ = y - V[:, 2]
V_ = np.delete(V, 2, axis=1)
z_, *_ = np.linalg.lstsq(V_, y_, rcond=None)
z_ = np.insert(z_, 2, 1)

print(z_)
# [ 0.04659264 -0.48453866  1.          0.19438046]

请注意，我在解决线性最小二乘问题后在系数向量中插入了1，我们不再求解a2，因为我们将其设置为1并将其从问题中删除。
为了完整起见，这是绘制时的解决方案：

我使用的完整代码：

import numpy as np

x = np.array([0.0, 1.0, 2.0, 3.0,  4.0,  5.0])
y = np.array([0.0, 0.8, 0.9, 0.1, -0.8, -1.0])

deg = 3
V = np.vander(x, deg + 1)
z, *_ = np.linalg.lstsq(V, y, rcond=None)

print(z)
# [ 0.08703704 -0.81349206  1.69312169 -0.03968254]

z = np.polyfit(x, y, deg)

print(z)
# [ 0.08703704 -0.81349206  1.69312169 -0.03968254]

y_ = y - V[:, 2]
V_ = np.delete(V, 2, axis=1)
z_, *_ = np.linalg.lstsq(V_, y_, rcond=None)
z_ = np.insert(z_, 2, 1)

print(z_)
# [ 0.04659264 -0.48453866  1.          0.19438046]

from matplotlib import pyplot as plt

plt.plot(x, y, 'o', label='data')
plt.plot(x, V @ z, label='polyfit')
plt.plot(x, V @ z_, label='constrained (a2 = 0)')

plt.legend()

plt.show()

下面是一个使用scipy.optimize.curve_fit的方法：
首先，让我们重新创建您的示例（作为健全性检查）：

import numpy as np
from scipy.optimize import curve_fit
​
def f(x, x3, x2, x1, x0):
    """this is the polynomial function"""
    return x0 + x1*x + x2*(x*x) + x3*(x*x*x)
​
popt, pcov = curve_fit(f, x, y)

print(popt)
#array([ 0.08703704, -0.81349206,  1.69312169, -0.03968254])

它与从np.polyfit()得到的值匹配。
现在添加x1的约束：

popt, pcov = curve_fit(
    f, 
    x,
    y,
    bounds = ([-np.inf, -np.inf, .999999999, -np.inf], [np.inf, np.inf, 1.0, np.inf])
)
print(popt)
#array([ 0.04659264, -0.48453866,  1.        ,  0.19438046])

我必须使用.999999999，因为下限必须严格小于上限。
或者，您可以将约束系数定义为常数，并获得其他3个的值：

def f_new(x, x3, x2, x0):
    x1 = 1
    return x0 + x1*x + x2*(x*x) + x3*(x*x*x)
popt, pcov = curve_fit(f_new, x, y)
print(popt)
#array([ 0.04659264, -0.48453866,  0.19438046])

这里有一个使用scipy.optimize.curve_fit的通用方法，目的是固定所需的多项式系数。

import numpy as np
from scipy.optimize import curve_fit

def polyfit(x, y, deg, which=-1, to=0):
    """
    An extension of ``np.polyfit`` to fix values of the vector 
    of polynomial coefficients. By default, the last coefficient 
    (i.e., the constant term) is kept at zero.

    Parameters
    ----------
    x : array_like
        x-coordinates of the sample points.
    y : array_like
        y-coordinates of the sample points.
    deg : int
        Degree of the fitting polynomial.
    which : int or array_like, optional
        Indexes of the coefficients to remain fixed. By default, -1.
    to : float or array_like, optional
        Values of the fixed coefficients. By default, 0.

    Returns
    -------
    np.ndarray
        (deg + 1) polynomial coefficients.
    """
    p0 = np.polyfit(x, y, deg)

    # if which == None it is reduced to np.polyfit
    if which is None: 
        return p0
    
    # indexes of the coeffs being fitted
    which_not = np.delete(np.arange(deg + 1), which)
    
    # create the array of coeffs
    def _fill_p(p):
        p_ = np.empty(deg + 1)  # empty array
        p_[which] = to          # fill with custom coeffs
        p_[which_not] = p       # fill with `p`
        return p_

    # callback function for fitting
    def _polyfit(x, *p):
        p_ = _fill_p(p)
        return np.polyval(p_, x)

    # get the array of coeffs
    p0 = np.delete(p0, which)                # use `p0` as initial condition
    p, _ = curve_fit(_polyfit, x, y, p0=p0)  # fitting
    p = _fill_p(p)                           # merge fixed and no-fixed coeffs
    return p

两个简单的例子来说明如何使用上面的函数：

import matplotlib.pyplot as plt

# just create some fake data (a parabola)
np.random.seed(0)                               # seed to reproduce the example
deg = 2                                         # second order polynomial
p = np.random.randint(1, 5, size=deg+1)         # random vector of coefficients
x = np.linspace(0, 10, num=20)                  # fake data: x-array
y = np.polyval(p, x) + 1.5*np.random.randn(20)  # fake data: y-array
print(p)                                        # output:[1, 4, 2]

# fitting
p1 = polyfit(x, y, deg, which=2, to=p[2])        # case 1: last coeff is fixed
p2 = polyfit(x, y, deg, which=[1,2], to=p[1:3])  # case 2: last two coeffs are fixed
y1 = np.polyval(p1, x)                           # y-array for case 1
y2 = np.polyval(p2, x)                           # y-array for case 2
print(p1)                                        # output: [1.05, 3.67, 2.]
print(p2)                                        # output: [1.08, 4.,   2.]

# plotting
plt.plot(x, y, '.', label='fake data: y = p[0]*x**2 + p[1]*x + p[2]')
plt.plot(x, y1, label='p[2] fixed at 2')
plt.plot(x, y2, label='p[2] and p[1] fixed at [4, 2]')
plt.legend()
plt.show()