Skip to content

Support online parameter estimation with an API (a.k.a. adaptive state estimation) #432

Description

@franckgaga

Supporting adaptive state estimation is quite straightforward on any state estimator that works with NonLinModel: you just augment the state vector with the estimated parameters $p_1$, $p_2$, etc. and augment the $\mathbf{f}$ equation with a time-invariant assumption $dp_1/dt = 0$, $dp_2/dt = 0$, etc.

Because of this, I first thought that adding an explicit API that would do this automatically is a very low priority feature. But, while writing this simple example of online parameter estimation with the MHE on the inverted pendulum:

using ModelPredictiveControl, Plots, JuMP
theme(:dark); default(fontfamily="Computer Modern"); scalefontsizes(1.1)
function f!(ẋ, x, u, _ , p)
    g, L, m = p                 # [m/s²], [m], [kg]
    θ, ω, K = x[1], x[2], x[3]  # [rad], [rad/s], [kg/s]
    τ  = u[1]                   # [Nm]
    ẋ[1] = ω
    ẋ[2] = -g/L*sin(θ) - K/m*ω + τ/m/L^2
    ẋ[3] = 0                    # dK/dt = 0 (the estimated parameter)
end
h!(y, x, _ , _ ) = (y[1] = 180/π*x[1]; nothing) # [°]
p = [9.8, 0.4, 0.3]
nu, nx, ny, Ts = 1, 3, 1, 0.1
vu, vx, vy = ["\$τ\$ (Nm)"], ["\$θ\$ (rad)", "\$ω\$ (rad/s)", "\$K\$ (kg/s)"], ["\$θ\$ (°)"]
model = setname!(NonLinModel(f!, h!, Ts, nu, nx, ny; p=p); u=vu, x=vx, y=vy)
σQ=[0.1, 1.0, 1.0]; σR=[5.0]; nint_ym=0; nint_u=0; 
transcription = OrthogonalCollocation(); hessian = true;
mhe = MovingHorizonEstimator(model; σQ, σR, nint_ym, nint_u, He=10, transcription, hessian)
setconstraint!(mhe, x̂min=[-Inf, -Inf, 1.0], x̂max=[+Inf, +Inf, 2.0], v̂min=[-2.5], v̂max=[2.5])
mhe |> display
N = 30; u = [1.10];
unset_time_limit_sec(mhe.optim)
res = sim!(mhe, N, u, x_0=[0, 0, 1.5], x̂_0=[0, 0, 1.2], y_noise=[0.5])
T = @elapsed sim!(mhe, N, u, x_0=[0, 0, 1.5], x̂_0=[0, 0, 1.2], y_noise=[0.5])
println("Sampling time: $Ts s, Average compute time: $(T/N) s")
plot(res, plotu=false, plotxwithx̂=true) |> display

I noticed that the states associated to the estimated parameter are in fact LTI models (the 0 in the RHS, that is, an integrator). It means that we could specialized the equality constraint as linear one on these states, for non-SingleShooting transcription. This is very similar to what is already done with the generic integrating states from nint_ym and nint_u options. The goal is to reduces the dimensions of the Jacobian thus more efficient AD. So this feature is more interesting to implement than I initially thought.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions