У меня есть функция 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 со случайными значениями от 0 до 1 для x и использовать их в функции ограничения def der().Идея заключалась в том, что каждое значение T и p имеет диапазон x от 0 до 1. К сожалению, я получаю следующую ошибку:
ValueError: operands could not be broadcast together with shapes (6,5) (6,)
Я знаю, что есть много других вопросов с этимошибка, но я не могу перенести ее в мою реальную проблему, поэтому любая помощь будет полезна.