# Solving ODE parameters using experimental data with control inputs

**URL:** <https://discourse.julialang.org/t/solving-ode-parameters-using-experimental-data-with-control-inputs/66614>\
**Category:** Numerics\
**Tags:** control, diffeq, sciml, experimentation\
**Created:** [August 18, 2021, 5:25pm UTC](https://discourse.julialang.org/t/solving-ode-parameters-using-experimental-data-with-control-inputs/66614 "2021-08-18T17:25:58Z")\
**Posts on this page:** 1\
**Page:** 2

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [August 21, 2021, 6:44am UTC](https://discourse.julialang.org/t/solving-ode-parameters-using-experimental-data-with-control-inputs/66614/21 "2021-08-21T06:44:28Z")

</div>

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

```julia
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]))

```

 ![tempplot](https://global.discourse-cdn.com/julialang/original/3X/4/6/46a267f1cd00b3f5d5d0a4396c63992009da596c.jpeg)

With printout

```julia
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.

[Previous page](https://discourse.julialang.org/t/solving-ode-parameters-using-experimental-data-with-control-inputs/66614.md?page=1)