# using Pkg; Pkg.add("ControlSystemIdentification"); Pkg.add("Optim"); Pkg.add("Plots")
In this notebook, we will explore system identification of an LTI state-space model using ControlSystemIdentification.jl Before any identification can begin, we need to load all the relevant packages. The package Optim does the heavy lifting under the hood, and we load it to be able to s
using ControlSystemIdentification, ControlSystems, Optim, Plots
using Random, LinearAlgebra
using ControlSystemIdentification: newpem
default(size=(500,280))
We start by creating a system to use as the subject of identification
function generate_system(nx,ny,nu)
A = randn(nx,nx)
A = A - A'
A -= 0.1I
A = exp(A)
B = randn(nx,nu)
C = randn(ny,nx)
sys = ss(A,B,C,0,1)
end
Random.seed!(1)
T = 300 # Number of time steps
nx = 3 # Number of poles in the true system
nu = 1 # Number of control inputs
ny = 1 # Number of outputs
x0 = randn(nx) # Initial state
sim(sys,u,x0=x0) = lsim(sys, u, 1:T, x0=x0)[1] # Helper function
sys = generate_system(nx,nu,ny)
u = randn(nu,T) # Generate random input
y = sim(sys, u, x0) # Simulate system
d = iddata(y,u,1)
sysh,opt = newpem(d, nx) # Estimate model
yh = predict(sysh, d) # Predict using estimated model
plot([y; yh]', lab=["y" "ŷ"]) # Plot prediction and true output, see also predplot
Iter Function value Gradient norm
0 3.752219e-01 2.846173e+00
* time: 7.390975952148438e-5
We could also use the function predplot to plot the predicted output together with the measurements.
In the example above, there was no noise on the measurement signal. This is of course not a realistic scenario, and we can increase the realism by simulating some measurement noise
# Both noises
σu = 0.2
σy = 0.2
sysn = generate_system(nx,ny,ny)
u = randn(nu,T)
un = u + sim(sysn, σu*randn(size(u)),0*x0)
y = sim(sys, un, x0)
yn = y + sim(sysn, σy*randn(size(u)),0*x0)
dn = iddata(yn, un, 1)
sysh,opt = newpem(dn,nx, iterations=400)
f1 = bodeplot([sys,sysh.sys], lab=["True system" "estimated system"])
f2 = predplot(sysh,dn)
f3 = simplot(sysh,dn)
plot(f1,f2,f3, layout=(1,3), size=(900,300))
Iter Function value Gradient norm
0 5.206939e+01 9.417315e+01
* time: 7.295608520507812e-5
50 2.823161e+01 6.061666e+02
* time: 0.042353153228759766
We can also optimize other metrics, such as the sum of absolute errors. One can also add a regularization function:
u = randn(nu,T)
un = u + sim(sysn, σu*randn(size(u)),0*x0)
y = sim(sys, un, x0)
yn = y + sim(sysn, σy*randn(size(u)),0*x0)
d = iddata(y, u, 1)
dn = iddata(yn, un, 1)
sysh,opt = newpem(dn,nx, metric=abs, regularizer=(p, P)->(0.1/T)*norm(p), iterations=3000)
plot(predplot(sysh,d),simplot(sysh,d))
Iter Function value Gradient norm
0 1.041061e+02 2.935336e+02
* time: 6.103515625e-5
50 7.822966e+01 9.149926e+01
* time: 0.2074291706085205
100 7.786358e+01 9.284712e+01
* time: 0.23303818702697754
We now create both test data and validation data and fit 4 different models of increasing order.
σy = 0.5
sim(sys,u) = lsim(sys, u, 1:T)[1]
sys = tf(1,[1,2*0.1,0.1])
sysn = tf(σy,[1,2*0.1,0.3])
# Training data
u = randn(nu,T)
un = u + 0.1randn(size(u))
y = sim(sys, u)
yn = y + sim(sysn, σy*randn(size(u)))
dn = iddata(yn, un, 1)
# Validation data
uv = randn(nu,T)
yv = sim(sys, uv)
ynv = yv + sim(sysn, σy*randn(size(uv)))
dnv = iddata(ynv, uv, 1)
InputOutput data of length 300 with 1 outputs and 1 inputs
# fit 4 different models
res = [newpem(dn,nx, iterations=50) for nx = 1:4]
Iter Function value Gradient norm
0 3.252749e+03 9.862239e+03
* time: 6.890296936035156e-5
50 2.166017e+02 8.595861e+02
* time: 0.03091287612915039
Iter Function value Gradient norm
0 5.947661e+01 5.473508e+02
* time: 4.410743713378906e-5
50 5.892303e+00 3.814770e+02
* time: 0.020408153533935547
Iter Function value Gradient norm
0 2.451672e+02 2.157385e+04
* time: 4.100799560546875e-5
50 1.634334e+01 8.497676e+03
* time: 0.0279541015625
Iter Function value Gradient norm
0 1.477232e+01 7.723547e+01
* time: 4.38690185546875e-5
50 4.040811e+00 3.030091e+01
* time: 0.031090974807739258
4-element Vector{Tuple{PredictionStateSpace{Discrete{Float64}, StateSpace{Discrete{Float64}, Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}, Vector{Float64}, Optim.MultivariateOptimizationResults{BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Nothing, Flat}, Vector{Float64}, Float64, Float64, Vector{OptimizationState{Float64, BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Nothing, Flat}}}, Bool, NamedTuple{(:f_limit_reached, :g_limit_reached, :h_limit_reached, :time_limit, :callback, :f_increased), NTuple{6, Bool}}}}}:
(PredictionStateSpace{Discrete{Float64}, StateSpace{Discrete{Float64}, Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}
A =
0.9670657208474821
B =
-0.0004164374412114029
C =
-2.4167593262889273
D =
0.0
Sample Time: 1.0 (seconds)
Discrete-time state-space model, [-0.040944514749169485], * Status: failure (reached maximum number of iterations)
* Candidate solution
Final objective value: 2.166017e+02
* Found with
Algorithm: BFGS
* Convergence measures
|x - x'| = 3.80e-11 ≰ 0.0e+00
|x - x'|/|x'| = 1.57e-11 ≰ 0.0e+00
|f(x) - f(x')| = 2.21e-09 ≰ 1.0e-16
|f(x) - f(x')|/|f(x')| = 1.02e-11 ≰ 0.0e+00
|g(x)| = 8.60e+02 ≰ 1.0e-12
* Work counters
Seconds run: 0 (vs limit 100)
Iterations: 50
f(x) calls: 506
∇f(x) calls: 51
)
(PredictionStateSpace{Discrete{Float64}, StateSpace{Discrete{Float64}, Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}
A =
0.7963016096125531 0.3269158686742356
-0.2385574675969522 0.8873845317749124
B =
-0.36974987482103927
0.2505451148539292
C =
2.4489005282409986 5.351530284538445
D =
0.0
Sample Time: 1.0 (seconds)
Discrete-time state-space model, [0.18289280793433368, -0.12418559151566994], * Status: failure (reached maximum number of iterations)
* Candidate solution
Final objective value: 5.892303e+00
* Found with
Algorithm: BFGS
* Convergence measures
|x - x'| = 3.40e-11 ≰ 0.0e+00
|x - x'|/|x'| = 6.35e-12 ≰ 0.0e+00
|f(x) - f(x')| = 2.12e-10 ≰ 1.0e-16
|f(x) - f(x')|/|f(x')| = 3.61e-11 ≰ 0.0e+00
|g(x)| = 3.81e+02 ≰ 1.0e-12
* Work counters
Seconds run: 0 (vs limit 100)
Iterations: 50
f(x) calls: 521
∇f(x) calls: 51
)
(PredictionStateSpace{Discrete{Float64}, StateSpace{Discrete{Float64}, Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}
A =
1.0099478050135495 -0.7043340672475571 0.0
0.018366358326455938 0.8106598838062368 0.2503963593824119
0.0 -0.24482325662901439 0.9416679797386908
B =
0.28064744382478163
-0.2583276995996772
0.2366985016459889
C =
0.29195903655208594 3.521459951378126 5.49123028244722
D =
0.0
Sample Time: 1.0 (seconds)
Discrete-time state-space model, [-63.34348312735423, -1.0998183695345702, 3.902723664042919], * Status: failure (reached maximum number of iterations)
* Candidate solution
Final objective value: 1.634334e+01
* Found with
Algorithm: BFGS
* Convergence measures
|x - x'| = 1.69e-10 ≰ 0.0e+00
|x - x'|/|x'| = 2.66e-12 ≰ 0.0e+00
|f(x) - f(x')| = 2.01e-09 ≰ 1.0e-16
|f(x) - f(x')|/|f(x')| = 1.23e-10 ≰ 0.0e+00
|g(x)| = 8.50e+03 ≰ 1.0e-12
* Work counters
Seconds run: 0 (vs limit 100)
Iterations: 50
f(x) calls: 574
∇f(x) calls: 51
)
(PredictionStateSpace{Discrete{Float64}, StateSpace{Discrete{Float64}, Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}
A =
0.8387483525696412 0.2566460269271575 0.0 0.0
-0.25736006451912846 0.8880914118295825 -0.09243079834418748 0.0
0.0 -0.0037598421121620143 0.35341217779684697 0.45792253124694243
0.0 0.0 -0.9377734800011042 1.285978706698473
B =
-0.4015425182052098
0.3526630684900554
0.042167901810490245
0.09026722979057122
C =
2.721582746521805 4.83653285714682 0.7933047899721857 -2.1992727898394184
D =
0.0
Sample Time: 1.0 (seconds)
Discrete-time state-space model, [0.19117193256134293, 0.6371869296145248, 0.3900848445398745, 1.896900818429575], * Status: failure (reached maximum number of iterations)
* Candidate solution
Final objective value: 4.040811e+00
* Found with
Algorithm: BFGS
* Convergence measures
|x - x'| = 3.37e-09 ≰ 0.0e+00
|x - x'|/|x'| = 6.96e-10 ≰ 0.0e+00
|f(x) - f(x')| = 7.58e-09 ≰ 1.0e-16
|f(x) - f(x')|/|f(x')| = 1.88e-09 ≰ 0.0e+00
|g(x)| = 3.03e+01 ≰ 1.0e-12
* Work counters
Seconds run: 0 (vs limit 100)
Iterations: 50
f(x) calls: 484
∇f(x) calls: 51
)
We can inspect the estimated noise model and compare it to the true noise system
ω = exp10.(range(-2, stop=log10(pi), length=150))
fig = plot(layout=4, size=(1000,600))
for i in eachindex(res)
(sysh,opt) = res[i]
predplot!(sysh,dnv; subplot=1, ploty=i==1)
simplot!(sysh,dnv; subplot=2, ploty=i==1)
end
bodeplot!(getindex.(res,1), ω, plotphase=false, subplot=3, title="Process", linewidth=2*[4 3 2 1])
bodeplot!([noise_model(r[1]) for r in res], ω, plotphase=false, subplot=4, linewidth=2*[4 3 2 1])
bodeplot!(sys, ω, plotphase=false, subplot=3, lab="True", linecolor=:blue, l=:dash, legend = :bottomleft, title="System model")
bodeplot!(ControlSystems.innovation_form(ss(sys),syse=σy^2*ss(sysn), R2=I), ω, plotphase=false, subplot=4, lab="True", linecolor=:blue, l=:dash, ylims=(0.1, 100), legend = :bottomleft, title="Noise model")
┌ Warning: This call is deprecated due to ambiguity, use covar(ss(D), R) or covar(ss(D, Ts), R) instead └ @ ControlSystems /home/fredrikb/.julia/dev/ControlSystems/src/ControlSystems.jl:191
We see that adding more than two poles to the model does little to improve the prediction and simulation performance, it is thus a good idea to stick with a second or third-order model.