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

Question

I have a function y(T,x,p). I have the data for T,p,x,y. With this data i want to know the coefficients so I can use the function to get any y I want to. So far I have this with using scipy.optimize.minimize and 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()

Because my data for x only reach to 0.3, the constraint (that the derivation after x is positive everwhere) only comply for this area. For higher x values it does not comply. So i thought I define multidimensional arrays T2,x2,p2 with random values between 0 and 1 for x and use them in the constraint function def der(). The Idea was that every T and pvalue has a x range from 0 to 1. Unfortunateley I am getting the following error:

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

I know there are many other questions with this error, but I cannot realy transfer it to my actual problem, so any help would be nice.

Answer 1

The error occurs because the solvers try's to match the output of the der(pars) function to a vector full of zeros. Essentially you have to construct your derivative so the return of the der(pars) function is of shape (6,1). You've correctly calculated the jacobian of the function. To flatten the jacobian, you could use the following approach:

As you're just interested in the inequality constrained and not in each value of the jacobian, you can return the minimum of each row. If the minimum is greater zero, than so is the whole jabobian. Try this code for your function der(pars). the function .min(axis=1) returns the minimal value of each row in the 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)

This yield the following result, using the remaining of your code:

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