How to minimize a function by several parameters? (part 1)

99 Views Asked by At

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
1

There are 1 best solutions below

7
Reinderien On BEST ANSWER

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.