Solving ODE parameters using experimental data with control inputs

For linear systems like this one, there is a lot of well developed methods to estimate the system parameters, as well as covariance properties and noise models etc. that are useful in a practical control scenario. For a temperature system like this, classical system identification could look something like

using ControlSystems, ControlSystemIdentification
plotly(show = false)
w = 2pi .* exp10.(LinRange(-3, log10(0.5), 500)) # frequency vector
G0 = tf(1, [10, 1]) # The true system, 10ẋ = -x + u
G = c2d(G0, 1) # discretize with a sample time of 1s
println("True system")
display(ss(G0))

# Data generation
u = sign.(sin.((0:0.01:10) .^ 2)) # sample a control input for identification
y, t, x = lsim(ss(G), u) # Simulate the true system to get test data
yn = y .+ 0.2 .* randn.() # add measurement noise
data = iddata(yn, u, t[2] - t[1]) # create a data object

# Estimation
na, nb = 1, 1 # number of parameters in denominator and numerator
Gh = arx(data, na, nb, estimator = wtls_estimator(data.y, na, nb)) # estimate an arx model
Gh2, noise_model = plr(data, na, nb, 1) # try another identification method

# Plot results
println("Estimated system")
display(ss(d2c(Gh))) # Convert from discrete to continuous time

bp = bodeplot(G, w, lab = "G (true)", hz = true, l = 5)
bodeplot!(Gh, w, lab = "arx", hz = true)
bodeplot!(Gh2, w, lab = "plr", hz = true, ticks = :default)

sp = stepplot(G, 150, lab="G (true)")
stepplot!(Gh, 150, lab = "arx")
stepplot!(Gh2, 150, lab = "plr", ticks = :default)
hline!([1], primary = false, l = (:black, :dash))

lp = lsimplot(ss(G), u, lab="G (true)")
lsimplot!(ss(Gh), u, lab = "arx")
lsimplot!(ss(Gh2), u, lab = "plr", ticks = :default)
plot!(data.t, yn[:], lab = "Estimation data")

plot(bp, sp, lp, layout = @layout([[a b]; c]))

With printout

True system
StateSpace{Continuous, Float64}
A = 
 -0.1
B = 
 1.0
C = 
 0.1
D = 
 0.0

Continuous-time state-space model

Estimated system
StateSpace{Continuous, Float64}
A = 
 -0.10260326778775565
B = 
 1.0
C = 
 0.10056519755546571
D = 
 0.0

Continuous-time state-space model

For practical identification, additional centering of the data is required.

2 Likes