optimize.basinhoppingの[setting an array element with a sequence]エラー

前提・実現したいこと

逆問題の求解をしています。
未知のパラメーターを含む連立微分方程式が、実測で得られているデータと最もFitするように、
scipy.optimize.basinhoppingで最適パラメーターを見つけたいです。
リファレンス：
https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.basinhopping.html#scipy.optimize.basinhopping

プログラムを実行すると、for ループが2回終了したところで、
「setting an array element with a sequence.」と表示されてしまいます。

途中からエラーを起こす理由がよくわかりません。
修正点を教えていただけないでしょうか？

コード

Python3
1import numpy as np
2from scipy.integrate import odeint
3import pandas as pd
4from scipy import optimize
5from scipy.optimize import basinhopping
6
7n = 3 #データ数
8date = "1812" #シート名
9
10D = pd.read_excel('前データ.xlsx', sheet_name=date,dtype='float64')
11DD = np.array(D)
12c = DD[5:12,0]
13
14def f(y, t, z): 
15#単位式
16    R1 = c[3] * y[0] * y[1]      
17    R2 = c[4] * c[0] * y[1]  
18    R3 = c[5] * c[1] * y[1]    
19    R4 = c[6] * y[4] * y[1]       
20    R5 = z[0] * y[0] * y[2]          
21    R6 = z[1] * y[0] * y[3]          
22    R7 = z[2] * y[1] * y[2]       
23    R8 = z[3] * y[1] * y[3]      
24#連立微分方程式
25    e0   = -R1 -c[2]*R5 -c[2]*R6
26    e1   = -R1 -R2 -R3 -R4 +z[4]*c[2]*R5 +z[5]*c[2]*R6 -c[2]*R7 -c[2]*R8
27    e2   = -z[6]*R5 -z[7]*R7
28    e3   = -z[8]*R6 -z[9]*R8
29    e4   = -R4
30    return [e0, e1, e2, e3, e4]
31
32def ODE(x, teta, i): 　　　　　　　　#数値積分の実行
33    f2 = lambda y,t: f(y, t, teta)
34    y0 = DD[0:5, i]
35    r = odeint(f2,y0,x)
36    r_O = scipy.stats.zscore(r[:,0])
37    r_P = scipy.stats.zscore(r[:,4])
38    return np.concatenate([r_O, r_P])
39
40
41def ODE_resid(p):　　　　　　　　　　#実測値と計算値の比較
42    for i in range(n):
43        x_data =  np.array(DD[22:31, i])
44        data_Oa = scipy.stats.zscore(DD[33:42,i])
45        data_Pa = scipy.stats.zscore(DD[44:53,i])
46        y_data =  np.concatenate([data_Oa, data_Pa])
47        if i == 0:
48            print(i) #何回ぐらい計算されているのか把握するために書いています。
49            Res0 = y_data - ODE(x_data,p,i))
50        if i == 1:
51            print(i)
52            Res1 = y_data - ODE(x_data,p,i))
53        if i == 2:
54            Res2 = y_data - ODE(x_data,p,i))
55            print(i)
56            Q = np.concatenate([Res0,Res1,Res2])
57            return Q
58    
59g = DD[12:22,0] #推定するパラメーター
60j = basinhopping(ODE_resid,g)

###エラー文

ValueError                                Traceback (most recent call last)
<ipython-input-32-a4a40533ee77> in <module>
     68           (g[6],g[6]*100),(g[7],g[7]*100),(g[8],g[8]*100),(g[9],g[9]*100)]
     69 
---> 70 j = basinhopping(ODE_resid,g)
     71 # jj = j.x
     72 # df = pd.DataFrame(jj)

~\Anaconda3.7\lib\site-packages\scipy\optimize\_basinhopping.py in basinhopping(func, x0, niter, T, stepsize, minimizer_kwargs, take_step, accept_test, callback, interval, disp, niter_success, seed)
    667 
    668     bh = BasinHoppingRunner(x0, wrapped_minimizer, take_step_wrapped,
--> 669                             accept_tests, disp=disp)
    670 
    671     # start main iteration loop

~\Anaconda3.7\lib\site-packages\scipy\optimize\_basinhopping.py in __init__(self, x0, minimizer, step_taking, accept_tests, disp)
     72 
     73         # do initial minimization
---> 74         minres = minimizer(self.x)
     75         if not minres.success:
     76             self.res.minimization_failures += 1

~\Anaconda3.7\lib\site-packages\scipy\optimize\_basinhopping.py in __call__(self, x0)
    284             return self.minimizer(x0, **self.kwargs)
    285         else:
--> 286             return self.minimizer(self.func, x0, **self.kwargs)
    287 
    288 

~\Anaconda3.7\lib\site-packages\scipy\optimize\_minimize.py in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
    593         return _minimize_cg(fun, x0, args, jac, callback, **options)
    594     elif meth == 'bfgs':
--> 595         return _minimize_bfgs(fun, x0, args, jac, callback, **options)
    596     elif meth == 'newton-cg':
    597         return _minimize_newtoncg(fun, x0, args, jac, hess, hessp, callback,

~\Anaconda3.7\lib\site-packages\scipy\optimize\optimize.py in _minimize_bfgs(fun, x0, args, jac, callback, gtol, norm, eps, maxiter, disp, return_all, **unknown_options)
    968     else:
    969         grad_calls, myfprime = wrap_function(fprime, args)
--> 970     gfk = myfprime(x0)
    971     k = 0
    972     N = len(x0)

~\Anaconda3.7\lib\site-packages\scipy\optimize\optimize.py in function_wrapper(*wrapper_args)
    298     def function_wrapper(*wrapper_args):
    299         ncalls[0] += 1
--> 300         return function(*(wrapper_args + args))
    301 
    302     return ncalls, function_wrapper

~\Anaconda3.7\lib\site-packages\scipy\optimize\optimize.py in approx_fprime(xk, f, epsilon, *args)
    728 
    729     """
--> 730     return _approx_fprime_helper(xk, f, epsilon, args=args)
    731 
    732 

~\Anaconda3.7\lib\site-packages\scipy\optimize\optimize.py in _approx_fprime_helper(xk, f, epsilon, args, f0)
    668         ei[k] = 1.0
    669         d = epsilon * ei
--> 670         grad[k] = (f(*((xk + d,) + args)) - f0) / d[k]
    671         ei[k] = 0.0
    672     return grad

ValueError: setting an array element with a sequence.