I would like to minimize a function dchi(a, b, c, d) by 4 parameters a,b,c and d.
It looks like there is an mistake in this line res = sp.optimize.minimize(dchi, first_guess), could you please tell me what exactly is written incorrectly and how to fix it?
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)
the errors:
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
To solve at least a couple of your problems, use a complex, properly unpack your parameters, and check root convergence:
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)
This does get further, but in context your function is problematic and eventually doesn't converge. Are there bounds on any of the variables you've shown? You really do need to find out what these should be.
Anyway, it's a bad idea to run an inner root-find. You should be running a single upper-level minimize with a vectorized root constraint; and don't loop.
With your chi function as
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
then a (slow, but) successful optimization can look like
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]
This will do better if you give it more iterations to run, and especially if you write a second Jacobian for the root.