'ValueError：无法将形状（5,5）的输入数组广播为形状（5）'，用于scipy.optimize-minimize

Question

我有一个功能y(T,x,p)。我有T，p，x，y的数据。有了这些数据，我想知道系数，所以我可以使用该函数来获得我想要的任何y。到目前为止，我使用scipy.optimize.minimize和method='cobyla'：

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

T = np.array([262,257,253,261,260,243], dtype=float)
p = np.array([25,22,19,24,24,14], dtype=float)
x = np.array([0.1,0.1,0.2,0.2,0.3,0.3], dtype=float)
y = np.array([10,9,13,16,20,12], dtype=float)

T2 = np.array([[262,262,262,262,262],[257,257,257,257,257],[253,253,253,253,253],[261,261,261,261,261],[260,260,260,260,260],[243,243,243,243,243]])
p2 = np.array([[25,25,25,25,25],[22,22,22,22,22],[19,19,19,19,19],[24,24,24,24,24],[24,24,24,24,24],[14,14,14,14,14]])
x2 = np.array([[0,0.25,0.5,0.75,1],[0,0.25,0.5,0.75,1],[0,0.25,0.5,0.75,1],[0,0.25,0.5,0.75,1],[0,0.25,0.5,0.75,1],[0,0.25,0.5,0.75,1]])

def func(pars, T, x, p): #my actual function
    a,b,c,d,e,f = pars
    return  x * p + x * (1 - x) * (a + b * T + c * T ** 2 + d * x + e * x * T + f * x * T ** 2) * p


def resid(pars): #residual function
    return ((func(pars, T, x, p) - y) ** 2).sum()

def der(pars): #constraint function: Derivation of func() after x positive everywhere
    a,b,c,d,e,f = pars
    return p2+p2*(2*x2*a+2*x2*b*T2+2*x2*c*T2**2+3*x2**2*d+3*x2**2*e*T2+3*x2**2*f*T2**2)+p2*(a+b*T2+c*T2**2+2*x2*d+2*e*x2*T2+2*f*x2*T2**2)


con1 = (dict(type='ineq', fun=der))
pars0 = np.array([0,0,0,0,0,0])
res = minimize(resid, pars0, method='cobyla',options={'maxiter': 500000}, constraints=con1)

print("a = %f , b = %f, c = %f, d = %f, e = %f, f = %f" % (res.x[0], res.x[1], res.x[2], res.x[3], res.x[4], res.x[5]))

T0 = 262.741   # plot an example graph y(x) for a certain T and p
x0 = np.linspace(0, 1, 100)
p0 = 26
fig, ax = plt.subplots()
fig.dpi = 80
ax.plot(x0, func(res.x, T0, x0, p0), '-')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

因为我对x的数据只达到0.3，所以约束（x之后的推导是正向的）只符合这个区域。对于更高的x值，它不符合。所以我认为我为T2定义了多维数组x2，p2，x，其随机值介于0和1之间，并在约束函数def der()中使用它们。想法是每个T和pvalue有一个x范围从0到1.不幸的是我收到以下错误：

ValueError: operands could not be broadcast together with shapes (6,5) (6,)

我知道这个错误还有很多其他问题，但我不能真正将它转移到我的实际问题，所以任何帮助都会很好。

Answer 1

发生错误是因为求解器尝试将der(pars)函数的输出与满为零的向量匹配。基本上你必须构造你的导数，所以der(pars)函数的返回是形状的（6,1）。你已经正确计算了函数的雅可比。要拼合雅可比，您可以使用以下方法：

由于您只对受限制的不等式感兴趣而不是对jacobian的每个值感兴趣，因此您可以返回每行的最小值。如果最小值大于零，那么整个jabobian就是这样。尝试使用此函数der(pars)的代码。函数.min(axis=1)返回jacobian中每行的最小值：

def der(pars): #constraint function: Derivation of func() after x positive everywhere
    a,b,c,d,e,f = pars
    jacobian = p2+p2*(2*x2*a+2*x2*b*T2+2*x2*c*T2**2+3*x2**2*d+3*x2**2*e*T2+3*x2**2*f*T2**2)+p2*(a+b*T2+c*T2**2+2*x2*d+2*e*x2*T2+2*f*x2*T2**2)
    return jacobian.min(axis=1)

这将产生以下结果，使用剩余的代码：

a = 1.312794 , b = -0.000001, c = -0.000084, d = 1.121216, e = -0.003784, f = 0.000129