'ValueError: could not broadcast input array from shape (5,5) into shape (5)' for scipy.optimize-minimize

'ValueError: could not broadcast input array from shape (5,5) into shape (5)' for scipy.optimize-minimize

我有一个函数y(T,x,p)。我有 Tpxy 的数据。有了这些数据,我想知道系数,这样我就可以使用该函数来获得我想要的任何 y 。到目前为止,我使用 scipy.optimize.minimizemethod='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 值,它不符合。 所以我想我定义了多维数组 T2x2p2,其中 x 的随机值介于 0 和 1 之间,并在约束函数 def der() 中使用它们。这个想法是每个 Tp 值都有一个从 0 到 1 的 x 范围。 不幸的是,我收到以下错误:

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

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

出现错误是因为求解器试图将 der(pars) 函数的输出与全为零的向量相匹配。本质上,您必须构建导数,以便 der(pars) 函数的 return 的形状为 (6,1)。您已经正确计算了函数的雅可比矩阵。要展平雅可比矩阵,您可以使用以下方法:

由于您只对受约束的不等式感兴趣,而不是对雅可比矩阵的每个值感兴趣,因此您可以 return 每行的最小值。如果最小值大于零,则整个 jabobian 也是如此。 为您的函数尝试此代码 der(pars)。函数 .min(axis=1) returns 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