跳至主要內容

MTK稳态系统仿真实例

Ai4Energy大约 4 分钟

MTK稳态系统仿真实例

!!! tip
Contents:MTK、组件化建模

Contributor: YJY

Email:522432938@qq.com

如有错误,请批评指正。

!!! note

MTK = ModelingToolkit.jl

[ModelingToolkit.jl Acausal Component-Based Modeling the RC Circuit](https://mtk.sciml.ai/stable/tutorials/acausal_components/)介绍了组件化建模的例子

动态与瞬态

从内涵来说,稳态问题是动态问题的子集。

对于稳态动态问题,当微分代数方程中对时间项的微分为0时,则由动态问题过渡为稳态问题。也就是说,稳态问题是动态问题的特殊形式,动态问题是问题的一般形式。微分在实际中代表着“变化”,当该“变化”项为0时,意味着“变化”停止,此时则进入了稳态阶段。

例如,对于这样一个问题:

dxdt=1xx0=0 \frac{dx}{dt} = 1-x\\x_0=0

从物理演化的角度看,初始时刻,x的变化率为1,x将会不断增大。系统发展到最后,当x增大到1的时候。dxdt=0\frac{dx}{dt}=0,此时此刻变化就停止了。所以x将进入稳态,即x=1x=1

这是因为,从物理的角度理解,微分代表一种变化趋势!上述方程的走势如下图,与分析的情况一致。

图 1
图 1

MTK的稳态求解技巧

求解以下稳态电路问题:

图 4
图 4

该问题的电路状态是确定值,不随时间而改变。而DE和MTK处理的是关于时间的微分方程,单纯的代数方程无法求解。对于这样稳态问题,应该怎样求解呢?

这里主要运用一个技巧——设定一个变化为0的量加入到方程中

什么意思呢?问题必须含有关于时间微分的微分方程,在这一基础之上,让dxdt=0\frac{dx}{dt}=0不就实现了一种稳态么?xx可以代表任何变量,电阻值、电源电压值、电流源值等等。于是,微分变为一种形式,只是让问题能够被求解器求解,对于问题的内涵并没有任何改变。

@MTK版本:8.5.5

using ModelingToolkit, Plots, DifferentialEquations

@variables t
@connector function Pin(; name)
    sts = @variables v(t) = 1.0 i(t) = 1.0 [connect = Flow]
    ODESystem(Equation[], t, sts, []; name=name)
end

function Ground(; name)
    @named g = Pin()
    eqs = [g.v ~ 0]
    compose(ODESystem(eqs, t, [], []; name=name), g)
end

function GroundDIFF(; name)
    @named g = Pin()
    sts = @variables v(t)
    D = Differential(t)
    eqs = [
        g.v ~ 0,
        D(v) ~ 0
    ]
    compose(ODESystem(eqs, t, sts, []; name=name), g)
end

function OnePort(; name)
    @named p = Pin()
    @named n = Pin()
    sts = @variables v(t) = 1.0 i(t) = 1.0
    eqs = [
        v ~ p.v - n.v
        0 ~ p.i + n.i
        i ~ p.i
    ]
    compose(ODESystem(eqs, t, sts, []; name=name), p, n)
end

function Resistor(; name, R=1.0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    ps = @parameters R = R
    eqs = [
        v ~ i * R
    ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function ConstantVoltage(; name, V=1.0)
    @named oneport = OnePort()
    @unpack v = oneport
    ps = @parameters V = V
    eqs = [
        V ~ v
    ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function ConstantCurrent(; name, I=1.0)
    @named oneport = OnePort()
    @unpack i = oneport
    ps = @parameters I = I
    # D = Differential(t)
    eqs = [
        i ~ I
    ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end


@named resistor1 = Resistor(R=5.0)
@named resistor2 = Resistor(R=5.0)
@named resistor3 = Resistor(R=30.0)
@named resistor4 = Resistor(R=20.0)
@named Isource = ConstantCurrent()
@named Vsource1 = ConstantVoltage(V=30.0)
@named Vsource2 = ConstantVoltage(V=5.0)
@named ground = GroundDIFF()


rc_eqs = [
    connect(Vsource1.p, resistor1.p)
    connect(resistor1.n, resistor2.p, Isource.p)
    connect(resistor2.n, resistor3.p, resistor4.p)
    connect(resistor4.n, Isource.n, Vsource2.p)
    connect(Vsource1.n, Vsource2.n, resistor3.n, ground.g)
]

@named _rc_model = ODESystem(rc_eqs, t)
@named rc_model = compose(_rc_model,
    [resistor1, resistor2, resistor3, resistor4,
        Isource, Vsource1, Vsource2, ground])

sys = structural_simplify(rc_model)

using Plots

u0 = [
    ground.v => 0.0
]
prob = ODAEProblem(sys, u0, (0, 10.0))
sol = solve(prob, Tsit5())
observed(sys)
plot(sol[resistor4.n.i])

代码中的关键在于函数GroundDIFF(),这个函数中添加的微分项为电势变化为0,即

dVdt=0 \frac{dV}{dt}=0

并且这个V没有任何意义,既不是接地的电势,也没有和电路中的任何变量产生关联。它只是为了让问题能够被求解而存在。


最后求解的结果为:

图 1
图 1

可以看到,求解的结果在时间轴上是不变的。随然问题看似是一个动态问题,本质上其实是稳态。因为其没有任何改变,时间轴上的任何状态与其前后状态都一样。

思想的核心其实就是——稳态是动态的一个子集,当动态过程中状态不再改变则进入稳态。这里不过是将问题设置为状态不发生变化的动态,其本质就是稳态

!!! tip
稳态电路本质是代数方程组,若加上了微分方程,则变成了微分代数方程组。代码中定义的问题是ODAEProblem而不是ODEProblem。DE只有求解微分代数方程的接口而没有求解代数方程的接口。所以需要构造成微分代数问题。

上次编辑于:
贡献者: Mingtao Li