2026a
# dae_index_lowering
自动降低 dae 指标
函数库: TySymbolicMath
# 语法
dae_index_lowering(sys::ODESystem; kwargs...) -> ODESystem
# 说明
自动降低 dae 指标，执行 Pantelides 算法将更高索引的 DAE 转换为索引 1 DAE。kwargs 被转发到 pantelides！鼓励最终用户改为调用 structure_simplify，它在内部调用此函数。示例
# 示例
自动减低 DAE 指标
在许多情况下，人们可能会不小心写下一个不容易通过数值方法求解的 DAE。
先将求解完整代码进行展示，后面讲逐步进行分析。
using TySymbolicMath
using DifferentialEquations
using LinearAlgebra
using TyPlot
using ModelingToolkit
function pendulum!(du, u, p, t)
x, dx, y, dy, T = u
g, L = p
du[1] = dx
du[2] = T*x
du[3] = dy
du[4] = T*y - g
du[5] = x^2 + y^2 - L^2
return nothing
end
pendulum_fun! = ODEFunction(pendulum!, mass_matrix=Diagonal([1,1,1,1,0]))
u0 = [1.0, 0, 0, 0, 0]
p = [9.8, 1]
tspan = (0, 10.0)
pendulum_prob = ODEProblem(pendulum_fun!, u0, tspan, p)
traced_sys = modelingtoolkitize(pendulum_prob)
pendulum_sys = structural_simplify(dae_index_lowering(traced_sys))
prob = ODAEProblem(pendulum_sys, [], tspan)
solver = pkgversion(ModelingToolkit) >= v"9.0.0" ? Rodas5P() : Tsit5()
sol = solve(prob, solver,abstol=1e-8,reltol=1e-8)
第一步构建摆系统
function pendulum!(du, u, p, t)
x, dx, y, dy, T = u
g, L = p
du[1] = dx; du[2] = T*x
du[3] = dy; du[4] = T*y - g
du[5] = x^2 + y^2 - L^2
end
pendulum_fun! = ODEFunction(pendulum!, mass_matrix=Diagonal([1,1,1,1,0]))
u0 = [1.0, 0, 0, 0, 0]; p = [9.8, 1]; tspan = (0, 10.0)
pendulum_prob = ODEProblem(pendulum_fun!, u0, tspan, p)
solve(pendulum_prob,Rodas4())
会出现如下 warning ，求解效果也并不好。
┌ Warning: dt(1.7763568394002505e-15) <= dtmin(1.7763568394002505e-15) at t=4.549693054277249e-5, and step error estimate = 225.75513606217686. Aborting. There is either an error in your model specification or the true solution is unstable.
事实证明，高阶 DAE 可以转换为低阶 DAE。
请注意，这在数学上等同于我们之前的方程，但 Jacobian w.r.t. 由于代数方程的 T 不再为零。 这意味着如果你写下这个版本的模型，它将是 index-1 并正确求解！ 事实上，DAE 指数通常是这样定义的：将 DAE 转换为 ODE 所需的微分次数，其中 ODE 是指数为 0 的 DAE，通过替换所有代数关系。
然而，要求用户坐在那里并在可能数以百万计的方程上完成这个过程是一种深不可测的精神负担。 但是，我们可以通过使用像 Pantelides 算法这样的方法来自动执行到索引 1 的缩减。虽然这需要 ModelingToolkit 符号形式，但我们使用modelingtoolkitize 将数字代码转换为符号代码，运行 dae_index_lowering 降低，然后转换回 使用 ODEProblem 计算数值代码，并使用数值求解器求解。
自动减少索引
traced_sys = modelingtoolkitize(pendulum_prob)
pendulum_sys = structural_simplify(dae_index_lowering(traced_sys))
prob = ODEProblem(pendulum_sys, Pair[], tspan)
sol = solve(prob, Rodas4())
plot(sol.t,sol[1,:],sol.t,sol[2,:],sol.t,sol[3,:],sol.t,sol[4,:],sol.t,sol[5,:])
/dae_index_lowering/dae1.svg)
# 输入参数
sys - 常微分方程系统
sys 为第一个参数，ode_order_lowering 将包含高指标的常微分方程系统降低为低指标常微分方程系统。
# 提示
如果想要了解 DAE 索引问题的详细情况，请参见这篇文章 (opens new window)。