DiffEqFlux.jl实现模型预测控制

DiffEqFlux.jl简介

DiffEqFlux.jl 是 SciML 生态系统的参数估计系统。它是一个高级接口，将所有工具与结合在一起。神经网络可以是模型的全部或一部分。它们可以在微分方程周围、成本函数中或微分方程内部。

控制问题

控制方程:
$$x^{\prime \prime }=u^3(t)$$

通过优化控制变量: $u(t)$ 使得损失函数最小:

$$Loss(\theta )=\sum _i||4-x(t_i)||+2||x^{\prime }(t_i)||+||u(t_i)||$$

在求解时，把微分方程都将为一阶：
$$x^{\prime }=v$$

$$v^{\prime }=u^3(t)$$

损失函数变为：

$$Loss(\theta )=\sum _i||4-x(t_i)||+2||v(t_i)||+||u(t_i)||$$

基于DifferentialEquations.jl建模方式实现控制

代码分析：

使用相关的函数包：

using DiffEqFlux, DifferentialEquations, Plots, Statistics

定义与初始化神经网络ann
ann为一个3层神经网络，第一层1输入32输出，第二次32输入32输出，第三层32输入1输出。

可将其看作一个黑盒函数，一个输入对应一个输出。在控制问题中，即输入时间域上的某个时间t，神经网络通过计算可以输出在该时间点最优的一个控制值$u$。经过训练的神经网络ann可以代表经过计算后的一段时间域上的控制曲线$u(t)$

ann = FastChain(FastDense(1,32,tanh), FastDense(32,32,tanh), FastDense(32,1))
θ = initial_params(ann) #初始化参数

定义微分方程，以DifferentialEquations方式定义
ann在问题中相当于一个函数u(t)，p是神经网络的参数，t是自变量。$u(t)=ann(t,p)$

function dxdt_(dx,x,p,t)
    x1, x2 = x
    dx[1] = x[2]
    dx[2] = ann([t],p)[1]^3
end

定义相关参数与构造问题：

tspan = (0.0f0,8.0f0)
x0 = [-4f0,0f0]
prob = ODEProblem(dxdt_,x0,tspan,θ)

尝试求解：

solve(prob,Vern9(),abstol=1e-10,reltol=1e-10)

定义预测函数与损失函数：
predict_adjoint函数计算在当前神经网络参数下u，v的值，即求解微分方程
loss_adjoint函数计算损失函数的值

ts = Float32.(collect(0.0:0.01:tspan[2]))
function predict_adjoint(θ)
  Array(solve(prob,Vern9(),p=θ,saveat=ts,sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true))))
end
function loss_adjoint(θ)
  x = predict_adjoint(θ)
  mean(abs2,4.0 .- x[1,:]) + 2mean(abs2,x[2,:]) + mean(abs2,[first(ann([t],θ)) for t in ts])/10
end

尝试计算在初始参数下，loss的大小：

l = loss_adjoint(θ)

设置callback
该函数可以显示优化过程中，loss的变化情况以及函数图像（画图可以删除省略）

cb = function (θ,l)
  println(l)
  p = plot(solve(remake(prob,p=θ),Tsit5(),saveat=0.01),ylim=(-6,6),lw=3)
  plot!(p,ts,[first(ann([t],θ)) for t in ts],label="u(t)",lw=3)
  display(p)
  return false
end

最后，训练优化:
res1,res2,res3进行了三轮优化

res1 = DiffEqFlux.sciml_train(loss_adjoint, θ, ADAM(0.005), cb = cb,maxiters=100)
res2 = DiffEqFlux.sciml_train(loss_adjoint, res1.u,
                              BFGS(initial_stepnorm=0.01), cb = cb,maxiters=100,
                              allow_f_increases = false)
function loss_adjoint(θ)
  x = predict_adjoint(θ)
  mean(abs2,4.0 .- x[1,:]) + 2mean(abs2,x[2,:]) + mean(abs2,[first(ann([t],θ)) for t in ts])
end

res3 = DiffEqFlux.sciml_train(loss_adjoint, res2.u,
                              BFGS(initial_stepnorm=0.01), cb = cb,maxiters=100,
                              allow_f_increases = false)

全部代码（训练完成之后会有画图查看结果等操作）：

using DiffEqFlux, DifferentialEquations, Plots, Statistics
tspan = (0.0f0,8.0f0)
ann = FastChain(FastDense(1,32,tanh), FastDense(32,32,tanh), FastDense(32,1))
θ = initial_params(ann)
function dxdt_(dx,x,p,t)
    x1, x2 = x
    dx[1] = x[2]
    dx[2] = ann([t],p)[1]^3
end
x0 = [-4f0,0f0]
ts = Float32.(collect(0.0:0.01:tspan[2]))
prob = ODEProblem(dxdt_,x0,tspan,θ)
solve(prob,Vern9(),abstol=1e-10,reltol=1e-10)
function predict_adjoint(θ)
  Array(solve(prob,Vern9(),p=θ,saveat=ts,sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true))))
end
function loss_adjoint(θ)
  x = predict_adjoint(θ)
  mean(abs2,4.0 .- x[1,:]) + 2mean(abs2,x[2,:]) + mean(abs2,[first(ann([t],θ)) for t in ts])/10
end
l = loss_adjoint(θ)
cb = function (θ,l)
  println(l)
  p = plot(solve(remake(prob,p=θ),Tsit5(),saveat=0.01),ylim=(-6,6),lw=3)
  plot!(p,ts,[first(ann([t],θ)) for t in ts],label="u(t)",lw=3)
  display(p)
  return false
end
# Display the ODE with the current parameter values.
cb(θ,l)
loss1 = loss_adjoint(θ)
res1 = DiffEqFlux.sciml_train(loss_adjoint, θ, ADAM(0.005), cb = cb,maxiters=100)
res2 = DiffEqFlux.sciml_train(loss_adjoint, res1.u,
                              BFGS(initial_stepnorm=0.01), cb = cb,maxiters=100,
                              allow_f_increases = false)
function loss_adjoint(θ)
  x = predict_adjoint(θ)
  mean(abs2,4.0 .- x[1,:]) + 2mean(abs2,x[2,:]) + mean(abs2,[first(ann([t],θ)) for t in ts])
end

res3 = DiffEqFlux.sciml_train(loss_adjoint, res2.u,
                              BFGS(initial_stepnorm=0.01), cb = cb,maxiters=100,
                              allow_f_increases = false)

l = loss_adjoint(res3.u)
cb(res3.u,l)
p = plot(solve(remake(prob,p=res3.u),Tsit5(),saveat=0.01),ylim=(-6,6),lw=3)
plot!(p,ts,[first(ann([t],res3.u)) for t in ts],label="u(t)",lw=3)

基于ModelingToolkit建模方式实现控制

基于modelingtoolkit方式的建模与differentialequations大同小异。

由于modelingtoolkit是符号建模体系，关键问题在于处理符号问题。上面神经网络的参数在modelingtoolkit中都只能当作parameters处理。同时神经网络的计算式都会被纳入符号体系中去计算。所以针对神经网络这种大量参数的运算，运算速度会下降很快。这是问题所在。

代码：

using DiffEqFlux, DifferentialEquations, Plots, Statistics,ModelingToolkit
tspan = (0.0f0,8.0f0)
ann = FastChain(FastDense(1,6,tanh), FastDense(6,6,tanh), FastDense(6,1))
θ = initial_params(ann)
N = length(θ)
function an(t,p)
    return ann(t,p)[1]^3
end
@register an(t,p)
@variables t,x(t),xx(t)
@parameters p[1:N]
D = Differential(t)
eqs=[ 
    D(x)~ xx
    D(xx) ~ an(t,p[1:N])
]
sys = ODESystem(eqs,t,[x,xx],p)
u0 =[
    x => -4.f0
    xx => 0.f0
]
paras = [p[i]=>θ[i] for i in 1:N]
prob = ODEProblem(structural_simplify(sys),u0,tspan,paras)
sol = solve(prob)
ts = Float32.(collect(0.0:0.01:tspan[2]))
function predict_adjoint(θ)
  Array(solve(prob,Vern9(),p=θ,saveat=ts,sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true))))
end
function loss_adjoint(θ)
  x = predict_adjoint(θ)
  mean(abs2,4.0 .- x[1,:]) + 2mean(abs2,x[2,:]) + mean(abs2,[first(ann([t],θ)) for t in ts])/10
end
loss_adjoint(θ)
l = loss_adjoint(θ)
cb = function (θ,l)
  println(l)
  return false
end
res1 = DiffEqFlux.sciml_train(loss_adjoint, θ, ADAM(0.005), cb = cb,maxiters=10)

神经网络的项数由32变为6，精度是大大下降的。这是因为32项要等太久了，等不到结果。这是使用modelingtoolkit的问题之一。计算量太大。

上面的参数用参数向量来定义，是通过在传递初始参数时，为每一个参数都进行了初始赋值。

如果打印sys也可以发现，神经网络的计算公式已经被嵌入到符号计算体系中去。