我想通过 4 个参数
dchi(a, b, c, d)abcdres = sp.optimize.minimize(dchi, first_guess)import functools
import scipy as sp
def CH(chi: float, T: float, a: float, b: float, c: float, d: float, H: float, TC: float, C:float) -> float:
return (a*(T - TC)*chi**(1/c)+ b * chi**(1/c) * chi**(1/d) * H**(1/d)- 1)+C/(T-TC)
exp=[[10.0052, 87.9716], [11.0433, 86.3866], [12.359, 84.8015], [13.7127, 83.2164], [14.9939, 82.0276], [16.3367, 80.8388], [17.6282, 80.0463], [18.9464, 78.8575], [18.9934, 78.8575], [18.9969, 78.8575], [18.9962, 78.8575], [19.4221, 78.4612], [20.8034, 77.6687], [22.1306, 76.4799], [23.4352, 75.6873], [24.7503, 75.291], [26.0295, 74.4985], [27.3692, 73.706], [30.0264, 72.9134], [31.3262, 72.5172], [32.6313, 72.1209], [33.9404, 71.7246], [35.2494, 71.3284], [36.551, 70.9321], [37.8445, 70.5358], [40.5813, 69.7433], [41.7696, 69.347], [43.0624, 69.347], [44.3609, 68.9507], [45.6536, 68.5545], [46.9616, 68.5545], [48.254, 68.1582], [49.5428, 67.7619], [52.1298, 67.3657], [53.4147, 67.3657], [54.7091, 66.9694], [56.0001, 66.9694], [57.2883, 66.5731], [58.5865, 66.5731], [59.8911, 66.1769], [62.4856, 66.1769], [63.7665, 65.7806], [65.0689, 65.7806], [66.3588, 65.3843], [67.6581, 65.3843], [68.9546, 65.3843], [70.2572, 64.9881], [72.8482, 64.9881], [74.1637, 64.9881], [75.4546, 64.5918], [76.7495, 64.5918], [78.0347, 64.5918], [79.3376, 64.1955], [80.6437, 64.1955], [83.5861, 64.1955], [84.8957, 63.7993], [86.2123, 63.7993], [87.5317, 63.7993], [88.8468, 63.7993], [90.1674, 63.7993], [91.4775, 63.403], [94.2769, 63.403], [95.5863, 63.403], [96.9007, 63.403], [98.2013, 63.403], [99.523, 63.0067], [100.839, 63.0067], [102.147, 63.0067], [104.812, 63.0067], [106.131, 63.0067], [107.448, 62.6104], [108.789, 62.6104], [110.099, 62.6104], [111.414, 62.6104], [112.756, 62.6104], [115.404, 62.6104], [116.737, 62.6104], [118.059, 62.6104], [119.376, 62.2142], [120.7, 62.2142], [122.025, 62.2142], [123.367, 62.2142], [126.006, 62.2142], [127.352, 62.2142], [128.674, 62.2142], [129.999, 62.2142], [131.333, 62.2142], [132.652, 61.8179], [133.987, 61.8179], [136.652, 61.8179], [137.987, 61.8179], [139.324, 61.8179], [140.645, 61.8179], [141.989, 61.8179], [143.316, 61.8179], [144.646, 61.8179], [147.323, 61.8179], [148.651, 61.4216], [149.978, 61.4216], [151.303, 61.4216], [152.637, 61.4216], [153.949, 61.4216], [155.268, 61.4216], [157.929, 61.4216], [159.252, 61.4216], [160.581, 61.0254], [161.889, 61.0254], [163.219, 61.0254], [164.546, 61.0254], [165.873, 61.0254], [168.524, 61.0254], [169.852, 60.6291], [171.173, 60.6291], [172.494, 60.6291], [173.83, 60.6291], [175.16, 60.6291], [176.473, 60.6291], [179.122, 60.2328], [180.459, 60.2328], [181.78, 60.2328], [183.106, 59.8366], [184.439, 59.8366], [185.757, 59.8366], [187.079, 59.8366], [189.73, 59.4403], [191.06, 59.4403], [192.386, 59.044], [193.717, 59.044], [195.037, 59.044], [196.372, 58.6478], [197.693, 58.6478], [200.361, 58.2515], [201.689, 57.8552], [203.008, 57.8552], [204.335, 57.459], [205.651, 57.459], [206.996, 57.0627], [208.324, 57.0627], [210.977, 56.6664], [212.292, 56.2702], [213.602, 55.8739], [214.947, 55.8739], [216.259, 55.4776], [217.583, 55.0813], [218.895, 55.0813], [221.56, 54.2888], [222.865, 53.8925], [224.195, 53.4963], [225.505, 53.1], [226.824, 53.1], [228.192, 52.7037], [229.508, 52.3075], [232.191, 51.5149], [233.507, 51.1187], [234.836, 50.7224], [236.161, 50.3261], [237.479, 49.9299], [238.807, 49.5336], [240.13, 48.741], [242.774, 47.9485], [244.096, 47.5522], [245.417, 47.156], [246.74, 46.7597], [248.057, 46.3634], [249.382, 45.5709], [250.457, 46.2508], [251.117, 46.0365], [252.897, 45.4166], [253.112, 45.3033], [254.115, 44.9356], [255.112, 44.5707], [256.107, 44.2179], [257.936, 43.554], [259.077, 43.1208], [260.395, 42.633], [260.555, 42.5688], [261.881, 42.0749], [263.515, 41.4803], [264.389, 41.108], [264.559, 41.047], [264.724, 40.9939], [264.891, 40.9305], [265.058, 40.8548], [265.219, 40.807], [265.384, 40.7401], [265.549, 40.6799], [265.707, 40.6116], [265.872, 40.5381], [266.041, 40.4919], [266.218, 40.4131], [266.379, 40.365], [266.532, 40.2826], [266.705, 40.2386], [266.876, 40.1762], [267.038, 40.1109], [267.201, 40.0532], [267.373, 39.9972], [267.543, 39.9247], [267.71, 39.8555], [267.874, 39.7915], [268.035, 39.7336], [268.194, 39.6811], [268.367, 39.6056], [268.542, 39.5387], [268.705, 39.4762], [268.866, 39.3997], [269.029, 39.3423], [269.201, 39.2861], [269.374, 39.22], [269.539, 39.1538], [269.698, 39.0823], [269.87, 39.0231], [270.036, 38.9709], [270.194, 38.8925], [270.371, 38.8394], [270.536, 38.7638], [270.697, 38.7051], [270.868, 38.6425], [271.038, 38.574], [271.208, 38.5115], [271.367, 38.4455], [271.54, 38.3782], [271.711, 38.3156], [271.868, 38.2524], [272.03, 38.1917], [272.2, 38.1262], [272.361, 38.0706], [272.525, 38.0021], [274.156, 37.3821], [274.229, 37.3188], [274.362, 37.27], [274.524, 37.2048], [274.705, 37.1411], [274.877, 37.0827], [275.029, 37.0104], [275.196, 36.9607], [275.37, 36.8851], [275.526, 36.8172], [275.683, 36.7629], [275.855, 36.6982], [276.025, 36.6346], [276.185, 36.5695], [276.345, 36.5063], [276.506, 36.4347], [276.68, 36.3789], [276.845, 36.315], [277.012, 36.2632], [277.184, 36.1689], [277.35, 36.1295], [277.517, 36.0462], [277.676, 35.9975], [277.839, 35.9304], [278.009, 35.8579], [278.18, 35.8111], [278.341, 35.7453], [278.5, 35.6707], [278.67, 35.6269], [278.836, 35.5542], [279.005, 35.486], [279.176, 35.4266], [279.333, 35.3622], [279.502, 35.2938], [279.671, 35.2304], [279.833, 35.1688], [279.997, 35.1014], [280.169, 35.0437], [280.337, 34.9787], [280.498, 34.9181], [280.67, 34.8453], [280.837, 34.7763], [280.995, 34.7195], [281.164, 34.6492], [281.332, 34.5889], [281.5, 34.523], [281.664, 34.4628], [281.83, 34.3953], [282.004, 34.3438], [282.164, 34.2835], [282.324, 34.214], [282.493, 34.1503], [282.67, 34.087], [282.833, 34.0348], [282.996, 33.9549], [283.164, 33.9051], [284.809, 33.291], [284.864, 33.2491], [284.993, 33.1839], [285.16, 33.111], [285.323, 33.0541], [285.487, 32.9876], [285.655, 32.9227], [285.821, 32.8639], [285.991, 32.8064], [286.164, 32.7342], [286.329, 32.6806], [286.494, 32.606], [286.66, 32.562], [286.824, 32.4928], [286.985, 32.4283], [287.149, 32.3634], [287.321, 32.2873], [287.487, 32.2322], [287.644, 32.1765], [287.817, 32.1208], [287.988, 32.0552], [288.146, 31.9866], [288.306, 31.9345], [288.476, 31.8621], [288.646, 31.8045], [288.815, 31.7536], [288.985, 31.6971], [289.14, 31.6201], [289.307, 31.551], [289.481, 31.4893], [289.645, 31.444], [289.81, 31.3815], [289.975, 31.311], [290.141, 31.252], [290.306, 31.1829], [290.466, 31.1326], [290.636, 31.0702], [290.807, 31.011], [290.967, 30.9446], [291.139, 30.8823], [291.304, 30.8283], [291.47, 30.7563], [291.646, 30.706], [291.809, 30.6516], [291.974, 30.578], [292.141, 30.5155], [292.306, 30.4571], [292.47, 30.3986], [292.634, 30.3283], [292.805, 30.2715], [292.971, 30.212], [293.127, 30.1523], [293.299, 30.0891], [293.475, 30.0263], [293.634, 29.9634], [293.794, 29.9016], [295.41, 29.3469], [295.467, 29.2859], [295.612, 29.2319], [295.787, 29.1886], [295.954, 29.1198], [296.116, 29.0857], [296.277, 29.0178], [296.445, 28.9598], [296.617, 28.9062], [296.78, 28.8415], [296.946, 28.7719], [297.111, 28.7037], [297.277, 28.6442], [297.444, 28.5661], [297.603, 28.5065], [297.772, 28.438], [297.945, 28.3781], [298.11, 28.3183], [298.27, 28.2553], [298.436, 28.1991], [298.608, 28.1362], [298.775, 28.089], [298.934, 28.0215], [299.089, 27.9536], [299.26, 27.8898], [299.431, 27.8379], [299.596, 27.7736], [299.77, 27.7224], [299.939, 27.6707], [300.1, 27.6075], [300.268, 27.5584], [300.436, 27.4834], [300.595, 27.445], [300.761, 27.3805], [300.928, 27.3191], [301.085, 27.2599], [301.248, 27.2049], [301.416, 27.1372], [301.582, 27.0904], [301.742, 27.0381], [301.902, 26.9595], [302.077, 26.9182], [302.254, 26.8639], [302.409, 26.7965], [302.573, 26.751], [302.744, 26.6902], [302.904, 26.6289], [303.075, 26.5635], [303.24, 26.5116], [303.398, 26.4599], [303.571, 26.3952], [303.75, 26.3447], [303.915, 26.2879], [304.07, 26.2416], [304.236, 26.1701], [304.396, 26.1088]]
#print(exp[0][0])
def dchi(a, b, c, d):
sd=0
for i in range(len(exp)):
CH_bound = functools.partial(CH, T=exp[i][0], a=a, b=b, c=c, d=d, H=315, TC=1000, C=2)
ch1 = sp.optimize.root_scalar(f=CH_bound, x0=0)
sd=sd+(ch1.root - exp[i][1]) ** 2
return sd
#print(dchi(1,1,1,1))
first_guess = [1,1,1,1]
res = sp.optimize.minimize(dchi, first_guess)
#print(res)
错误：
Traceback (most recent call last):
File "C:\Users\...\PycharmProjects\pythonProject6\main.py", line 23, in <module>
res = sp.optimize.minimize(dchi, first_guess)
File "C:\Users\...\PycharmProjects\pythonProject6\venv\lib\site-packages\scipy\optimize\_minimize.py", line 705, in minimize
res = _minimize_bfgs(fun, x0, args, jac, callback, **options)
File "C:\Users\...\PycharmProjects\pythonProject6\venv\lib\site-packages\scipy\optimize\_optimize.py", line 1419, in _minimize_bfgs
sf = _prepare_scalar_function(fun, x0, jac, args=args, epsilon=eps,
File "C:\Users\...\PycharmProjects\pythonProject6\venv\lib\site-packages\scipy\optimize\_optimize.py", line 383, in _prepare_scalar_function
sf = ScalarFunction(fun, x0, args, grad, hess,
File "C:\Users\...\PycharmProjects\pythonProject6\venv\lib\site-packages\scipy\optimize\_differentiable_functions.py", line 158, in __init__
self._update_fun()
File "C:\Users\...\PycharmProjects\pythonProject6\venv\lib\site-packages\scipy\optimize\_differentiable_functions.py", line 251, in _update_fun
self._update_fun_impl()
File "C:\Users\...\PycharmProjects\pythonProject6\venv\lib\site-packages\scipy\optimize\_differentiable_functions.py", line 155, in update_fun
self.f = fun_wrapped(self.x)
File "C:\Users\...\PycharmProjects\pythonProject6\venv\lib\site-packages\scipy\optimize\_differentiable_functions.py", line 137, in fun_wrapped
fx = fun(np.copy(x), *args)
TypeError: dchi() missing 3 required positional arguments: 'b', 'c', and 'd'
Process finished with exit code 1
要解决至少几个问题，请使用复杂的、正确解压参数，并检查根收敛性：
import functools
import numpy as np
import scipy as sp
def CH(
chi: float, T: float, a: float, b: float, c: float, d: float, H: float, TC: float, C: float,
) -> float:
chi = chi + 0j
err = (
a*(T - TC)*chi**(1/c)
+ b * chi**(1/c) * chi**(1/d) * H**(1/d)
+ C/(T-TC)
- 1
)
return np.abs(err)
exp = (
(10.0052, 87.9716), (11.0433, 86.3866), (12.359, 84.8015),
(169.852, 60.6291), (171.173, 60.6291), (172.494, 60.6291), (173.83, 60.6291),
)
def dchi(params: np.ndarray) -> float:
a, b, c, d = params
sd = 0
x0 = 0
for exp_a, exp_b in exp:
CH_bound = functools.partial(CH, T=exp_a, a=a, b=b, c=c, d=d, H=315, TC=1000, C=2)
ch1 = sp.optimize.root_scalar(
f=CH_bound, x0=x0,
)
assert ch1.converged
sd += (ch1.root - exp_b)**2
x0 = ch1.root
return sd
first_guess = (1,1,1,1)
res = sp.optimize.minimize(
fun=dchi,
x0=first_guess,
)
assert res.success, res.message
print(res.x)
这确实更进一步，但在上下文中你的函数是有问题的并且最终不会收敛。您显示的任何变量是否有界限？你确实需要找出这些应该是什么。
无论如何，运行内部寻根并不是一个好主意。您应该运行一个带有矢量化根约束的上层
minimize您的气功能为
import numpy as np
from scipy.optimize import NonlinearConstraint, minimize
def CH(
chi: np.ndarray, T: float, H: float, TC: float, C: float, a: float, b: float, c: float, d: float,
) -> np.ndarray:
err = (
a*(T - TC)*chi**(1/c)
+ b * chi**(1/c) * chi**(1/d) * H**(1/d)
+ C/(T-TC)
- 1
)
return err
然后（缓慢但）成功的优化看起来像这样
def dchi(params: np.ndarray) -> float:
chi = params[4:]
chi_err = chi - exp_chi
err = chi_err.dot(chi_err)
return err
def dchi_jac(params: np.ndarray) -> np.ndarray:
chi = params[4:]
jac = np.zeros_like(params)
jac[4:] = 2*(chi - exp_chi)
return jac
def dchi_root(params: np.ndarray) -> np.ndarray:
# for each value of exp_T, there is some unknown optimal value of chi such that
# CH(chi, T) == 0
abcd, chi = np.split(params, (4,))
a, b, c, d = abcd
root_err = CH(
chi=chi, T=exp_T, a=a, b=b, c=c, d=d, H=315, TC=1000, C=2,
)
return root_err
exp_T, exp_chi = np.array(exp).T
first_guess = np.concatenate((
(0.79, 0.00893, 0.428, 1.38),
exp_chi,
))
result = minimize(
fun=dchi, jac=dchi_jac,
x0=first_guess,
constraints=NonlinearConstraint(
fun=dchi_root,
lb=0, ub=0,
),
options={'maxiter': 10, 'disp': True},
)
print(result.message)
print(f'Chi error: {result.fun:.3f}')
print(f'Root error: {np.abs(dchi_root(result.x)).sum():.3f}')
print(f'Parameters: {result.x[:4]}')
Iteration limit reached (Exit mode 9)
Current function value: 1543.01943809294
Iterations: 10
Function evaluations: 16
Gradient evaluations: 10
Iteration limit reached
Chi error: 1543.019
Root error: 57.920
Parameters: [-1.61941337e-05 -7.62007792e-06 6.02274220e-01 1.32563818e+00]
如果你给它更多的迭代次数来运行，效果会更好，特别是如果你为根编写第二个雅可比行列式。