跳至主要內容

组织方程的艺术——组件化（入门）

Ai4Energy大约 9 分钟

组织方程的艺术——组件化（入门）

!!! tip
Contents：建模思想

Contributor: YJY

Email:522432938@qq.com

如有错误，请批评指正。

前言

Ai4Energy开发者总是张口闭口组件化建模，那么组件化建模的内涵到底是什么？为什么要组件化建模？这篇我们从头至尾介绍做一个组件化思想的入门介绍。

建模用ModelingToolkit.jlopen in new window；求解用DifferentialEquations.jlopen in new window。

电路问题

以电路系统切入，深入探讨组件化建模的内涵是最合适不过的了。

首先准备2个电路元件的方程。

  • 欧姆定律：电阻两端的电压与流过电阻电流成正比。

U=IRU=IR

  • 电容器：电容两端电压的变化率与电流成正比。

dVdt=IC\frac{dV}{dt} = \frac{I}{C}

!!! tip
C是电容量，单位为法拉（F）。

我们再准备一个电路系统，如下图:

图 5
图 5

需要求解的问题是：电容初始电压为0，当开关闭合后，电容电压电流的变化情况是？

手搓求解

手搓的意思是电路中学的传统解法。开关闭合后，视为导线。采用节点电压法求解。标记三个节点1，2，3。节点电压分别为V1,V2,V3V_1,V_2,V_3，如下图。整个电路只有一条主路，所以电流为II。

图 6
图 6

所以电路的方程为：

V1−V3=1(V)V1−V2=1(Ω)∗Id(V2−V3)dt=I1(F)=IV3=0(1)V_1-V_3=1 (V) \\ V_1-V_2=1 (\Omega) * I\\ \frac{d(V_2-V_3)}{dt} = \frac{I}{1 (F)} = I\\ V_3 = 0 \tag{1}

手动化简之后可以得到2种情况：

  • 保留II，化简掉V1,V3,V2V_1,V_3,V_2。开关刚闭合，有电容视为导线，则0时刻，I=1I=1

1−Idt=−dIdt=II0=1(A)(2)\frac{1-I}{dt} = -\frac{dI}{dt}= I\\ I_0 = 1(A) \tag{2}

  • 保留V2V_2，化简掉V1,V3,IV_1,V_3,I

dV2dt=1−V2V2,0=0(V)(3)\frac{dV_2}{dt} = 1-V_2\\ V_{2,0}=0(V) \tag{3}

以上微分方程式2或3，求解任何一个，把结果代回方程组1，整个电路系统的情况就可以得到了。这个问题就解决了。

下面看一看组件化怎么做！

组件化求解

在之前的基础之上，我们只做一件事情——引入额外的变量。

为每个器件添加两个端口变量（地只加一个端口）。正极叫p，负极叫n。同时为器件取名，电源叫做source，电阻叫做resistor，电容叫做capacitor，地叫做ground。同时规定，器件索引其子要素的方式为“.”。例如，电阻正极电压的变量全名为resistor.p.Vresistor.p.V。如下图所示。

图 2
图 2

那么我们新引入了14个变量：

  • 电源的4个

source.n.Vsource.n.isource.p.Vsource.p.isource.n.V\\ source.n.i\\ source.p.V\\ source.p.i\\

  • 电阻的4个

resistor.n.Vresistor.n.iresistor.p.Vresistor.p.iresistor.n.V\\ resistor.n.i\\ resistor.p.V\\ resistor.p.i\\

  • 电容的4个

capacitor.n.Vcapacitor.n.icapacitor.p.Vcapacitor.p.icapacitor.n.V\\ capacitor.n.i\\ capacitor.p.V\\ capacitor.p.i\\

  • 地的两个

ground.port.Vground.port.iground.port.V\\ ground.port.i\\

再研究一下变量关系，我们可以发现。所有端口电流（地端口除外）的数值都等于II，端口电压值都等于对应的电压节点的电压。

这，貌似就是给V1,V3,V2,IV_1,V_3,V_2,I这几个变量加了很多别名？一顿操作猛如虎，细想发现全白干，形式主义我最行？

下面才是主菜！

Step1: 把V1,V3,V2,IV_1,V_3,V_2,I全踢掉，如下图。

图 3
图 3

Step2. 分别写组件的方程，规定电流入流为正，流出为负：

  • 电源：

source.p.V−source.n.V=1(V)source.p.I+source.n.I=0source.p.V - source.n.V = 1 (V)\\ source.p.I + source.n.I = 0

  • 电阻：

resistor.p.V−resistor.n.V=resistor.p.I∗1(Ω)resistor.p.I+resistor.n.I=0resistor.p.V - resistor.n.V = resistor.p.I * 1 (\Omega)\\ resistor.p.I + resistor.n.I = 0

  • 电容：

d(capacitor.p.V−capacitor.n.V)dt=capacitor.p.I1(F)capacitor.p.I+capacitor.n.I=0\frac{d(capacitor.p.V - capacitor.n.V)}{dt} = \frac{capacitor.p.I}{1 (F)} \\ capacitor.p.I + capacitor.n.I = 0

  • 地（规定电势为0的点，电流不做要求）：

ground.port.V=0ground.port.I=0ground.port.V = 0\\ ground.port.I = 0

Step3: 构建连接关系

  • 电源和电阻连（等价于节点1，节点电压相等、电流净流量为0）

source.n.V=resistor.p.Vsource.n.I+resistor.p.I=0source.n.V = resistor.p.V\\ source.n.I + resistor.p.I = 0

  • 电阻和电容连（等价于节点2，节点电压相等、电流净流量为0）

resistor.n.V=capacitor.p.Vresistor.n.I+capacitor.p.I=0resistor.n.V = capacitor.p.V\\ resistor.n.I + capacitor.p.I = 0

  • 电容、电源和地连（等价于节点3，节点电压相等、电流净流量为0）

capacitor.n.V=sourc.p.V=ground.port.Vcapacitor.n.I+sourc.p.I+ground.port.I=0capacitor.n.V = sourc.p.V = ground.port.V\\ capacitor.n.I + sourc.p.I + ground.port.I= 0\\

稍微思考一下可以发现，上述10多个方程组成的方程大组，其数学本质与方程小组1完全等价！唯一的区别是前者多引入了很多变量，并且这些变量完全可以消去！

但是，前者完成了组件化的封装，由此连接关系与组件内部之间成了相互独立的关系。

什么意思呢？以V1V_1为例。V1V_1既是电压源的出口电压又是电阻的正极电压。同时V1V_1作为节点电压，已经消去了电流的概念。这就意味每节点是不同的。而组件化的内涵就是把这一点拆解成了两部分：连接与组件内部。以电阻的正负极端口为例，对内建立组件内部的方程：

resistor.p.V−resistor.n.V=resistor.p.I∗1resistor.p.V - resistor.n.V = resistor.p.I * 1

其本质就是U=IRU=IR；对外建立与其他组件连接：

source.p.V=resistor.p.Vsource.n.I+resistor.p.I=0resistor.n.V=capacitor.p.Vresistor.n.I+capacitor.p.I=0source.p.V = resistor.p.V\\ source.n.I + resistor.p.I = 0\\ resistor.n.V = capacitor.p.V\\ resistor.n.I + capacitor.p.I = 0

其本质就是节点电压1和2.

这么多方程很复杂但逻辑清晰。其实它很简单，因为化简计算机会去做。

求对比解

手搓版

把节点电压法的分析过程，写成ModelingToolkit的代码。

!!! tip
* 在VScode中，输入入V_3就能输入V₃
图 1
* 在VScode中，按下Ctrl+Enter可以执行光标所在行的一行代码。以下的代码可以一行一行执行，看中间结果。

using ModelingToolkit, DifferentialEquations

# 定义变量
@variables t
@variables V₁(t) V₂(t) V₃(t) I(t)

# 定义微分
D = Differential(t)

# 设置参数
R = 1.0
C = 1.0
V = 1.0

# 输入方程
rc_eqs = [
    V₁ - V₃ ~ V
    V₁ - V₂ ~ I * R
    D(V₂) ~ I / C
    V₃ ~ 0
]

# 构建系统
@named rc_model = ODESystem(rc_eqs, t)

# 系统化简
sys = structural_simplify(rc_model)

# 设置初值
u0 = [
    V₂ => 0.0
]
# 求解时间范围
tspan = (0.0, 10.0)

# 构建问题并求解
prob = ODAEProblem(sys, u0, tspan)
sol = solve(prob, Tsit5())

# 分别查看 V₁, V₂, V₃ 的变化
sol[V₂]
sol[V₁]
sol[V₃]

再查看一下化简前的方程：

equations(rc_model)

再查看一下化简后的方程：

equations(sys)

方程中V1=1V_1=1，和推导的结果是一致的。

组件版

把组件化的过程，写成代码。

using ModelingToolkit, DifferentialEquations

# 定义独立时间变量
@variables t

# 器件端口作为连接点
@connector function Pin(; name)
    sts = @variables v(t) = 1.0 i(t) = 1.0 [connect = Flow]
    return ODESystem(Equation[], t, sts, []; name=name)
end

# 地，地的端口名字改为g，不是port
function Ground(; name)
    @named g = Pin()
    eqs = [g.v ~ 0; g.i ~ 0]
    return compose(ODESystem(eqs, t, [], []; name=name), g)
end
# 电阻元件
function Resistor(; name, R=1.0)
    @named p = Pin()
    @named n = Pin()
    ps = @parameters R = R
    eqs = [
        p.v - n.v ~ p.i * R
        0 ~ p.i + n.i
    ]
    return compose(ODESystem(eqs, t, [], ps; name=name), p, n)
end
# 电容元件，因为不能对表达式进行微分，所以再引入一个电容的电压差变量v
function Capacitor(; name, C=1.0)
    @named p = Pin()
    @named n = Pin()
    ps = @parameters C = C
    sts = @variables v(t) = 1.0
    D = Differential(t)
    eqs = [
        v ~ p.v - n.v
        D(v) ~ p.i / C
        0 ~ p.i + n.i
    ]
    return compose(ODESystem(eqs, t, sts, ps; name=name), p, n)
end
# 电压源
function ConstantVoltage(; name, V=1.0)
    @named p = Pin()
    @named n = Pin()
    ps = @parameters V = V
    eqs = [
        V ~ p.v - n.v
        0 ~ p.i + n.i
    ]
    return compose(ODESystem(eqs, t, [], ps; name=name), p, n)
end

# 定义组件
R = 1.0
C = 1.0
V = 1.0
@named resistor = Resistor(R=R) # @name可以视为给定义的组件起了个名字，Resistor(R=R)返回的组件名字就叫resistor
@named capacitor = Capacitor(C=C)
@named source = ConstantVoltage(V=V)
@named ground = Ground()

# 构建连接关系
rc_eqs = [
    connect(source.p, resistor.p)
    connect(resistor.n, capacitor.p)
    connect(capacitor.n, source.n, ground.g)
]
@named _rc_model = ODESystem(rc_eqs, t) #连接关系也需要放到ODESystem中。

# 组件与组件连接关系一起构建系统
@named rc_model = compose(_rc_model,[resistor, capacitor, source, ground])
equations(rc_model) # 查看方程

# 系统化简
sys = structural_simplify(rc_model)
equations(sys) # 查看方程

# 定义初值
u0 = [
    capacitor.v => 0.0
]

# 求解
prob = ODAEProblem(sys, u0, (0, 10.0))
sol = solve(prob, Tsit5())

# 查看
sol[capacitor.v]
sol[resistor.n.v]
sol[capacitor.p.v]

小结

引入的组件端口变量作为内部方程与外部连接关系的桥梁，使得组件完整又独立，可以与任意外部组件相连接。连接就去干连接的活，组件内部就管组件内部的方程。

从另一个角度说，节点电压法需要针对每一个系统做比较复杂的系统分析，因为其包含了某种化简的内涵；而组件化模式则是批量定义组件和批量处理连接关系，最后统一化简。后者更加易于理解与设计计算算法。

组件化的意义就在于高效复用以提高系统构建效率。因为组件都是由开发者写好的，用户只需要构建连接关系，系统即可搭建好。这也是为什么说组件化是关于方程组织的艺术了，本质还是那个本质，只是换一种变量与方程的组织模式。电路系统的组件化是组件化思想的体现之一。

最后，思考如何用节点电压法解下图的电路，以及如何用组件搭建。

图 2
图 2

其实思考之后就会发现组件化构建的优点：

  • 高复用性：节点电压法构建需要自己写组件的特性方程；组件化不需要，不管它内部是什么牛鬼蛇神，组件写好了就能拿来用。之后就一个事——连接。定义组件的过程就是对组件内部方程的重用。
  • 高易用性：一旦组件设计好了，构建连接关系不需要再动脑子了，照着系统图哪里有节点连哪里不就完事了？剩下精髓，是如何去设计一个系统，不再需要考虑组件底层。

!!! tip
提问：我们上面手敲代码定义组件和构建连接的过程，与在Simulink中通过图形化界面拖一个组件连一个线的过程，有本质的区别吗？

如果你觉得没有区别，恭喜，你已经大致掌握流程模拟仿真软件底层架构的逻辑了。

练习

尝试用上面写好的组件构建上图系统并求解。同样的，都假设电容电压初值为0。

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