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