pyomo的IPOPT求解微分方程（基于scipy的示例修改）

网路末端遗传因子

已于 2023-11-01 13:51:27 修改

阅读量932

点赞数 7

文章标签： scipy

于 2023-11-01 12:54:39 首次发布

本文链接：https://blog.csdn.net/qq_54394719/article/details/134158983

版权

示例源于scipy的微分方程组求解包solve_ivp使用详解_网路末端遗传因子的博客-CSDN博客

在scipy下使用的代码为：

from scipy.integrate import solve_ivp
import numpy as np

new_v = 50
def upward_cannon(t, y): return [y[1], -0.5]
def hit_ground(t, y): 
    if (y[0] - 2000)>=0:
        y[1] = new_v
    return y[0] - 2000
hit_ground.terminal = True
hit_ground.direction = -1
sol = solve_ivp(upward_cannon, [0, 100], [0, 10], events=hit_ground)

在pyomo下改写为：

from pyomo.environ import *
from pyomo.dae import *

m = ConcreteModel()

m.t = ContinuousSet(bounds=(0, 100)) # 时间范围
m.y = Var(m.t, domain=Reals) # y的定义
m.v = Var(m.t, domain=Reals) # v的定义
m.dydt = DerivativeVar(m.y, wrt=m.t) # y对t的导数
m.dvdt = DerivativeVar(m.v, wrt=m.t) # v对t的导数

# 初始条件
m.y[0] = 0
m.v[0] = 10

# 微分方程
def _ode1(m, t):
    return m.dydt[t] == m.v[t]
m.ode1 = Constraint(m.t, rule=_ode1)

def _ode2(m, t):
    return m.dvdt[t] == -0.5
m.ode2 = Constraint(m.t, rule=_ode2)

# 终止条件
# 当炮弹击到达指定高度时，速度会突然改变。
# 这是一个不连续的行为，不能直接用连续的微分方程来描述。
# 因此这里忽略速度变更。
def _terminal(m, t):
    return m.y[t] <= 2000
m.terminal = Constraint(m.t, rule=_terminal)

# 目标函数
m.obj = Objective(expr=m.t[-1], sense=minimize)

# 求解
# 创建一个离散化器对象，用于将连续的微分方程转化为离散的差分方程。*1
discretizer = TransformationFactory('dae.finite_difference')
# 将离散化器应用于模型的时间变量，使用前向差分法进行离散化。*2*5
discretizer.apply_to(m, nfe=50, method='forward')

# 创建一个求解器对象，使用的求解器是IPOPT（Interior Point OPTimizer）。*3 *4
solver = SolverFactory('ipopt')
results = solver.solve(m)

# 使用 m.load(results) 将求解结果加载到模型中
m.load(results)

# 查看变量 y 的结果
for i in m.y:
    print(f"y[{i}] = {m.y[i].value}")

# 查看变量 t 的结果
for i in m.t:
    print(f"t[{i}] = {m.t[i].value}")

在数学建模中，我们经常遇到连续的微分方程，但是在求解这些方程时，通常需要将其转化为离散的差分方程。离散化的过程将连续的时间或空间分割成一系列离散的点，然后通过近似方法来计算这些点上的函数值。

在这个例子中，dae.finite_difference 是一种常用的离散化方法，它使用有限差分法来近似微分方程。通过使用 TransformationFactory 类创建 discretizer 对象，并指定离散化方法为 dae.finite_difference，我们可以将原始的连续微分方程转化为离散的差分方程，以便进行数值求解。

在离散化过程中，我们需要指定离散化的参数，如离散化的时间步数、离散化的变量以及使用的差分方案。在这个例子中，model 是待离散化的模型对象，nfe=100 指定了离散化的时间步数为100，wrt=model.t 指定了离散化的变量为模型的时间变量，scheme='BACKWARD' 指定了使用后向差分方案进行离散化。

通过调用离散化器对象的 apply_to 方法，并传入相应的参数，我们可以将模型中的时间变量离散化为一系列离散的时间点，并使用后向差分方案来近似微分方程。这样，我们就可以在离散的时间点上进行数值求解，以获得微分方程的近似解。

在数学建模和优化问题中，选择合适的求解器是非常重要的。IPOPT 是一种开源的非线性优化求解器，被广泛应用于各种数学建模和优化问题中。它采用了内点法（Interior Point Method）来求解非线性优化问题，具有较高的求解效率和稳定性。

在这个例子中，我们使用 IPOPT 求解器来求解离散化后的微分方程模型。IPOPT 可以通过迭代的方式，寻找使得目标函数最小化或约束条件满足的最优解。通过将离散化后的微分方程模型转化为一个非线性优化问题，我们可以使用 IPOPT 求解器来寻找模型的最优解，即在离散的时间点上求解微分方程的近似解。

选择 IPOPT 求解器的原因是它在处理非线性优化问题方面表现良好，并且具有广泛的应用和支持。然而，根据具体的问题和求解要求，也可以选择其他的求解器来求解离散化后的微分方程模型。

solver = SolverFactory(‘ipopt’)这一步除了使用ipopt还能使用什么求解器？

除了 IPOPT 求解器，还可以使用其他求解器来求解优化问题。在 Pyomo 中，可以使用以下求解器：

1. `glpk`：GNU Linear Programming Kit，用于线性规划问题。
2. `gurobi`：Gurobi Optimizer，用于线性规划、整数规划和混合整数规划问题。
3. `cplex`：IBM ILOG CPLEX Optimization Studio，用于线性规划、整数规划和混合整数规划问题。
4. `cbc`：COIN-OR Branch and Cut，用于混合整数规划问题。
5. `baron`：BARON，用于全局优化问题。
6. `scip`：SCIP Optimization Suite，用于整数规划和混合整数规划问题。
7. `knitro`：Artelys Knitro，用于非线性规划问题。
8. `bonmin`：COIN-OR Bonmin，用于非线性整数规划问题。
9. `ipopt`：Interior Point OPTimizer，用于非线性规划问题。

这些求解器具有不同的特点和适用范围，可以根据具体的问题需求选择合适的求解器。在使用这些求解器之前，需要确保已经安装了相应的求解器软件，并将其配置到 Pyomo 中。

method 参数是用来选择离散化方法的。Pyomo提供了多种离散化方法，包括：

‘forward’：前向差分法
‘backward’：后向差分法
‘central’：中心差分法
‘lagrange’：拉格朗日插值法
‘pseudo_spectral’：伪谱法
这些方法各有优缺点，选择哪一种取决于你的具体问题和需求。例如，‘forward’ 和 ‘backward’ 方法比较简单，但可能不够精确；‘central’ 方法精度较高，但计算量也较大；‘lagrange’ 和 ‘pseudo_spectral’ 方法可以达到很高的精度，但计算量更大，且可能需要更复杂的算法支持。