在 python 中求解没有初始解的非线性方程（matlab 有效，python 无效）

Posted 2023-02-18

【中文标题】在 python 中求解没有初始解的非线性方程（matlab 有效，python 无效）

【英文标题】：Solve nonlinear equations without initial solution in python (matlab works, python does not)

【发布时间】：2021-11-05 11:10:23

【问题描述】：

目前我尝试计算以下形式的指数函数： f(x) = aexp(-bx)+c

对于上下文：这应该是二氧化碳排放量的减少函数。

条件如下：

f(x=0) = 700 Mt（德国当前的二氧化碳排放量） f(x=30) = 0 Mt（30 年内无二氧化碳排放） 积分 0->30 f(x) = 预算（剩余的二氧化碳预算 - 曲线下面积）

我可以使用 matlab (solve-function) 解决这个问题

但我想用 python 来做，但我尝试的所有求解器都失败了。当我使用 matlab 脚本的解决方案作为初始猜测时，一切正常。但我不想那样做。如果我不知道一个好的初始猜测，我想使用一个 python 求解器。

我应该使用哪一个？到目前为止，我从 scipy 和 gekko 尝试了 fsolve。

这是我的壁虎尝试的代码：

def solve(self, end_year, budget):
        self.end_year = end_year
        self.budget = budget
        d_year = self.end_year - self.start_year
        
        m = GEKKO()
        a,b,c = [m.Var(1) for i in range(3)]
        
        eq1 = a + c - self.co2_start
        eq2 = a*m.exp(-b*d_year) + c
        eq3 = a/b*(1-m.exp(-b*d_year))+c*d_year-self.budget
        
        m.Equation([eq1==0, eq2==0, eq3==0])
        m.solve(disp=False)
        self.a, self.b, self.c = a.value, b.value, c.value
´´´

【问题讨论】：

在我看来，分析（精确）解决方案应该是可能的。 a = -c 因为exp(0) = 1 不管b 的值如何，都可以从那里拿走。

@Thomas 实际上是 a+c = co2_start。我认为没有分析解决方案。 matlab 告诉我找不到解析解

嗯，好的，再想一想，由于积分，我不太确定是否存在解析解。但无论如何，您都可以通过消除其中一个未知数来简化求解器的工作。也许这允许你使用一个简单的初始猜测，比如将剩余的两个未知数设置为 0。

@Thomas 哦，太好了，谢谢，这确实有效！我设法消除了 a 和 c，最后只得到了 b 的方程式

【参考方案1】：

我没有认真考虑过这个问题的形式，但这里有一段使用scipy 的代码对我有用，不需要分析解决方案或简化：

from math import exp
import numpy as np
import scipy.integrate
import scipy.optimize
import matplotlib.pyplot as plt

co2_start = 700
co2_end = 0
year = 30
budget = 5000 # Not mentioned in the question

def func(x, a, b, c):
    return a * exp(-b * x) + c

def func_integral(x, a, b, c):
    # In this case it's easy to work out the integral analytically, but maybe
    # you want a more sophisticated `func` later.
    return scipy.integrate.quad(func, 0, x, (a, b, c))[0]

def err_func(x):
    try:
        a, b, c = x
        err_start = func(0, a, b, c) - co2_start
        err_end = func(year, a, b, c) - co2_end
        err_area = func_integral(year, a, b, c) - budget
        return [err_start**2, err_end**2, err_area**2]
    except OverflowError:
        return [np.inf, np.inf, np.inf]

solution = scipy.optimize.fsolve(err_func, [0, 1, 0])
print(solution)

a, b, c = solution
xs = np.arange(0, year)
ys = [func(x, a, b, c) for x in xs]
plt.plot(xs, ys)
plt.show()

我设置 a, b, c = 0, 1, 0 是因为我们知道 b 必须是正数，而 b = 0 区域对求解器来说是不愉快的。

它从quad 函数中吐出一个IntegrationWarning，但我希望这仅适用于远离实际解决方案的极端情况。

来自err_func，我在这里返回一个误差向量，其中包含每个条件的平方误差。这是一种常见的方法，但在这种情况下，它似乎也可以在没有平方的情况下工作。

另一个因素是处理OverflowError，因为求解器会尝试b 的值，其顺序为1e6，对于exp 函数来说太大了。我们通过返回无穷大来引导求解器远离该区域。