In this example, we will estimate a model for a four-stage evaporator to reduce the water content of a product, for example milk. The 3 inputs are feed flow, vapor flow to the first evaporator stage and cooling water flow. The three outputs are the dry matter content, the flow and the temperature of the outcoming product.

The example comes from STADIUS's Identification Database

Zhu Y., Van Overschee P., De Moor B., Ljung L., Comparison of three classes of identification methods. Proc. of SYSID '94,

using DelimitedFiles, Plots
using ControlSystemIdentification, ControlSystemsBase

url = "https://ftp.esat.kuleuven.be/pub/SISTA/data/process_industry/evaporator.dat.gz"
zipfilename = "/tmp/evaporator.dat.gz"
path = Base.download(url, zipfilename)
run(`gunzip -f $path`)
data = readdlm(path[1:end-3])
# Inputs:
# 	u1: feed flow to the first evaporator stage
# 	u2: vapor flow to the first evaporator stage
# 	u3: cooling water flow
# Outputs:
# 	y1: dry matter content
# 	y2: flow of the outcoming product
# 	y3: temperature of the outcoming product
u = data[:, 1:3]'
y = data[:, 4:6]'
d = iddata(y, u, 1)

InputOutput data of length 6305, 3 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=6)

We split the data in two, and use the first part for estimation and the second for validation. A model of order around 8 is reasonable (the paper uses 6-13). This system requires the option zeroD=false to be able to capture a direct feedthrough, otherwise the fit will always be rather poor.

dtrain = d[1:3300] # first experiment ends after 3300 seconds
dval = d[3301:end]

model,_ = newpem(dtrain, 8, zeroD=false)

┌ Warning: x_tol is deprecated. Use x_abstol or x_reltol instead. The provided value (0) will be used as x_abstol.
└ @ Optim ~/.julia/packages/Optim/8dE7C/src/types.jl:110
Iter     Function value   Gradient norm
     0     7.796309e+02     2.166337e+02
 * time: 4.601478576660156e-5
    50     7.515894e+02     7.099459e+02
 * time: 4.6249308586120605
   100     7.370387e+02     9.003382e+01
 * time: 8.230781078338623
   150     7.319665e+02     6.576515e+01
 * time: 11.751981973648071
   200     7.300869e+02     1.292878e+01
 * time: 15.277271032333374
   250     7.297392e+02     1.681885e+01
 * time: 18.81157398223877
   300     7.296118e+02     1.322206e+01
 * time: 22.335543870925903
   350     7.295775e+02     5.851955e+00
 * time: 25.87279200553894
   400     7.295727e+02     2.074114e+00
 * time: 29.388813972473145
   450     7.295717e+02     1.132825e+00
 * time: 32.91208791732788

predplot(model, dval, h=1, layout=d.ny)
predplot!(model, dval, h=5, ploty=false)

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

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

w = exp10.(LinRange(-2, log10(pi/d.Ts), 200))
sigmaplot(model.sys, w, lab="PEM", plotphase=false)

Let's compare prediction performance to the paper

ys = predict(model, dval, h=5)
ControlSystemIdentification.mse(dval.y-ys)

3×1 Matrix{Float64}:
 0.05761868864741539
 0.1563550174672708
 0.019269556712616372

The authors got the following errors: [0.24, 0.39, 0.14]