“scipy.optimize.minimize”如何强制系数不为零

Question

我正在尝试拟合一个函数 p，它取决于两个变量 x,T。 p,T,x 的数据通过带有pandas 的 Excel 工作表提供。以下代码运行良好。

import pandas as pd
import os
from scipy.optimize import minimize
import numpy as np

df = pd.read_excel(os.path.join(os.path.dirname(__file__), "./DataTest.xlsx"))
df = df.sort_values('x')

T = np.array(df['T'], dtype=float)
x = np.array(df['x'], dtype=float)
p = np.array(df['p'], dtype=float)
p0 = 67.17

def cav2(pars, T, x): # function p(T,x)
    a,b,c,d,e,f = pars
    return x * p0 + x * (1 - x) * (a + b * T + c * T ** 2 + d * x + e * x * T + f * x * T ** 2) * p0

def resid(pars, T, x):
    return ((p - cav2(pars, T, x)) ** 2).sum()

def constr(pars):
    return np.gradient(cav2(pars, T, x))

con1 = {'type': 'ineq', 'fun': constr}
pars0 = np.array([0.1, 0.1, 0.1, 0.1, 0.1, 0.1], dtype=float)
res = minimize(resid, pars0, args=(T, x), method='cobyla', options={'maxiter': 50000}, 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]))

最后一个print给了我函数的系数：

a = 2.891584 , b = 0.000000, c = -0.000000, d = 0.792256, e = -0.000000, f = 0.000000

这让我想到了我的实际问题。因为一些系数变为零，它使函数p(T,x) 独立于T，这是我不想要的。需要明确的是，目前cav2(res.x, 300, 0.1) 给出的结果与例如cav2(res.x, 500, 0.1) 相同。

scipy.optimize.minimize 中是否有一种（简单的）方法可以强制所有系数接受大于零的值？

谢谢

Answer 1

一些优化器支持对系数的约束（例如 L-BFGS-B）。

import pandas as pd
import os
from scipy.optimize import minimize
import numpy as np

T = np.random.normal(10)
x = np.random.normal(10)

p0 = 67.17

# Fake true parameters
a, b, c, d, e, f = np.random.uniform(-1, 1, size=6)

# targets
p = x * p0 + x * (1 - x) * (a + b * T + c * T ** 2 + d * x + e * x * T + f * x * T ** 2) * p0


def cav2(pars, T, x): # function p(T,x)
    a, b, c, d, e, f = pars
    return x * p0 + x * (1 - x) * (a + b * T + c * T ** 2 + d * x + e * x * T + f * x * T ** 2) * p0


def resid(pars, T, x):
    return ((p - cav2(pars, T, x)) ** 2).sum()


def constr(pars):
    return np.gradient(cav2(pars, T, x))

# this will force all parameters to be positive
bounds = [(0, None), (0, None), (0, None), (0, None), (0, None), (0, None)]
pars0 = np.array([0.1, 0.1, 0.1, 0.1, 0.1, 0.1], dtype=float)

res = minimize(resid, pars0, args=(T, x), method='L-BFGS-B', options={'maxiter': 50000}, bounds=bounds)

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]))

边界的工作方式是(lower, upper)，设置None 意味着没有应用边界。因此，如果您不想在第一个参数上设置边界，例如，您可以将边界替换为：

[(None, None), (0, None), (0, None), (0, None), (0, None), (0, None)]