In this example, we will estimate a model for a glass furnace.

We will get the data from STADIUS's Identification Database

using DelimitedFiles, Plots
using ControlSystemIdentification, ControlSystemsBase

url = "https://ftp.esat.kuleuven.be/pub/SISTA/data/process_industry/glassfurnace.dat.gz"
zipfilename = "/tmp/furnace.dat.gz"
path = Base.download(url, zipfilename)
run(`gunzip -f $path`)
data = readdlm(path[1:end-3])
u = data[:, 2:4]'
y = data[:, 5:10]'
d = iddata(y, u, 1)

InputOutput data of length 1247, 6 outputs, 3 inputs, Ts = 1

The input consists of two heating inputs and one cooling input, while there are 6 outputs from temperature sensors in a cross section of the furnace.

Before we estimate any model, we inspect the data

plot(d, layout=9)

We split the data in two, and use the first part for estimation and the second for validation. This system requires zeroD=false to be able to capture a direct feedthrough to output 4, otherwise the fit for output 4 will always be rather poor.

dtrain = d[1:2end÷3]
dval = d[2end÷3:end]

model = subspaceid(dtrain, 7, zeroD=false)

We can have a look at the $D$ matrix in the estimated model

model.D

6×3 Matrix{Float64}:
 -0.0162229   -0.0143999   0.214359
 -0.00356019  -0.00725287  0.0991567
 -0.00343267  -0.00755522  0.0986383
 -0.0712456   -0.0359864   0.665711
 -0.011395    -0.0131599   0.127809
 -0.0103743   -0.00876221  0.133655

indeed, the (4,3) element is rather large.

We validate the model by prediction on the validation data:

predplot(model, dval, h=1, layout=6)
predplot!(model, dval, h=10, ploty=false)

The figures above show the result of predicting $h={1, 10}$ steps into the future.

We can visualize the estimated model in the frequency domain as well.

w = exp10.(LinRange(-3, log10(pi/d.Ts), 200))
sigmaplot(model.sys, w, lab="MOESP")