基于组件模型的参数辨识

DiffEqParamEstim简介

DiffEqParamEstim.jl 是DifferentialEquations.jl 生态中的一个组件包。它拥有参数估计工具。

下面以RC电路模型为例，进行电路中正弦电压源的电压值估计。

RC组件建模

基于ModelingToolkit，建立RC电路模型，模型如下图所示。这是一个震荡电路

代码

首先完成电路模型的组件编写，一共有5个组件，分别是电源、电阻、电容、电感、接地

using ModelingToolkit, DifferentialEquations
# Basic electric components
@parameters t

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

function ModelingToolkit.connect(::Type{Pin}, ps...)
    eqs = [
           0 ~ sum(p->p.i, ps) # KCL
          ]
    # KVL
    for i in 1:length(ps)-1
        push!(eqs, ps[i].v ~ ps[i+1].v)
    end

    return eqs
end

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

function Resistor(;name, R = 1.0)
    val = R
    @named p = Pin()
    @named n = Pin()
    @variables v(t)
    @parameters R
    eqs = [
           v ~ p.v - n.v
           0 ~ p.i + n.i
           v ~ p.i * R 
          ]
    ODESystem(eqs, t, [v], [R], systems=[p, n], defaults=Dict(R => val), name=name)
end
function Capacitor(;name, C = 1.0)
    val = C
    @named p = Pin()
    @named n = Pin()
    @variables v(t)
    @parameters C
    D = Differential(t)
    eqs = [
           v ~ p.v - n.v
           0 ~ p.i + n.i
           D(v) ~ p.i / C
          ]
    ODESystem(eqs, t, [v], [C], systems=[p, n], defaults=Dict(C => val), name=name)
end

function Inductor(; name, L = 1.0)
    val = L
    @named p = Pin()
    @named n = Pin()
    @variables v(t) i(t)
    @parameters L
    D = Differential(t)
    eqs = [
           v ~ p.v - n.v
           0 ~ p.i + n.i
           i ~ p.i
           D(i) ~ v / L
          ]
    ODESystem(eqs, t, [v, i], [L], systems=[p, n], defaults=Dict(L => val), name=name)
end

function ChangeableVoltage(;name,vol=16.0)
    @named p = Pin()
    @named n = Pin()
    @variables v(t)
    @parameters Vol
    eqs = [
           v ~ Vol*sin(2π*t)
           v ~ p.v - n.v
           0 ~ p.i + n.i
          ]
    ODESystem(eqs, t, [v], [Vol], systems=[p, n], defaults=Dict(Vol => vol),name=name)
end

可以看到，系统一共定义有4个参数，分别是电阻值，电容值，电感值以及电压的最大值

---

然后，通过连接函数，组建系统

@named resistor = Resistor(R=3.0)
@named capacitor = Capacitor(C=1.0/24)
@named source = ChangeableVoltage(vol=10.0)
@named inductor = Inductor(L=0.1)
@named ground = Ground()


rc_eqs = [
          connect(source.p, capacitor.p)
          connect(capacitor.n, inductor.p)
          connect(inductor.n, resistor.p)
          connect(source.n,resistor.n,ground.g)
         ]

sys = [resistor, inductor, capacitor,source,ground]
@named rc_model = ODESystem(rc_eqs, t, systems=sys)

sys = structural_simplify(rc_model)

using Plots
u0 = [
    capacitor.v => 0.0
    capacitor.p.i => 0.0
    inductor.i => 0
    inductor.v => 0
     ]
P=[3.0,1.0/24,0.1,10.0]  #系统的参数
prob = ODAEProblem(sys, u0, (0, 10.0),P)
sol = solve(prob, Tsit5())

p1 = plot(sol,vars=[capacitor.v capacitor.p.i],xlims = (0,10),ylim = (-10,15))
p2 = plot(sol,vars=[inductor.v inductor.i],xlims = (0,10),ylim = (-5,5))
plot(p1,p2,layout=(2,1))

运行代码可以得到结果图：

电容的电压电流、电感的电压电流都是震荡的。

初值向量P分别对应4个参数的参数值，电压值为最后一个，设定为10.0

在不知道参数的顺序时，可以使用parameters函数查看系统参数。

parameters(sys)

建立问题

为了能够模拟参数辨识，给其中一个变量加上扰动。可以使用states函数查看系统变量（在组件设计过程中定义的）。

states(sys)

这里对两个变量都添加扰动。

using RecursiveArrayTools
t = collect(range(0,stop=10,length=1000)) # 建立时间向量
randomized = VectorOfArray([(sol(t[i]) + .5randn(2)) for i in 1:length(t)])
data = convert(Array,randomized)  

对sol的结果，选取了1000个点作为样本点加入扰动。

加入的扰动的方式是：通过生成绝对值小于0.5的随机数加入到从sol中选取出来的样本点中去，将其作为需要参数辨识的样本数据。

看一看扰动生成的结果：

p1=plot(t,data[1,:],ylim=(-8,8))
p1=plot!(t,data[2,:],ylim=(-8,8))
p2=plot(sol,vars=[inductor.i,capacitor.v],ylim=(-8,8))
plot(p1,p2,layout=(2,1))

产生的结果图为：

辨识过程

参数辨识使用DiffEqParamEstim，

using DiffEqParamEstim
cost_function = build_loss_objective(prob,Tsit5(),
L2Loss(t,data),maxiters=10000,verbose=false)

通过build_loss_objective建立了辨识模型，即通过该函数可以求处理后的样本点与模型计算结果的方差和（L2Loss），当然方差和只是作为拟合的评价指标之一，可以选取不同的函数或者构建不同的指标。

接下来，看一看不同的电压值，带来的方差和的变化。

vals = 0:0.1:20.0
plot(vals,[cost_function([3.0,1.0/24,0.1,i]) for i in vals],yscale=:log10,
     xaxis = "Parameter", yaxis = "Cost", title = "1-Parameter Cost Function",
     lw = 3)

val为从0开始，到20，步长为0.1。对这些电压值，去匹配模型，算出来的与样本值的方差和会不同。

结果为：

可以看到，在10.0的地方，方差和最小（因为设定的值就是10.0，扰动是在此基础之上叠加的）。说明辨识结果有效。

以上是通过“肉眼”观察的出来的结果。科学的方法是求一组参数，使得方差和最小。这是一个最优化问题。使用Optim

using Optim
result = optimize(cost_function, [3.0,1.0/24,0.1,15])
result.minimizer

将电压的初始值设为15去寻优，得到的结果为：

4-element Vector{Float64}:
 2.975801390547115
 0.04266590691854741
 0.09999565290932136
 9.923607580088433

和真实值是非常接近的！寻优有效！

也可以将4个初始值都改变去寻优。

result = optimize(cost_function, [2.0,0.1,0.5,16])
result.minimizer

得到的结果也很好

4-element Vector{Float64}:
 2.972037379476059
 0.04300492580277582
 0.09998659716890046
 9.911215185371102