### arxar: Complete example with simulation and estimation

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/api

Demonstrates ARXAR model estimation workflow including transfer function definition, system simulation with process and noise models, and data generation. Shows setup of A, B, G, D, and H polynomials for a discrete-time system with sample time 1 second.

```julia
N = 500

sim(G, u) = lsim(G, u, 1:N)[1][:]

A = tf([1, -0.8], [1, 0], 1)

B = tf([0, 1], [1, 0], 1)

G = minreal(B / A)

D = tf([1, 0.7], [1, 0], 1)

H = 1 / D

u, e = randn(1, N), randn(1, N)

y, v = sim(G, u), sim(H * (1/A), e) # simulate process
```

--------------------------------

### Install ControlSystemIdentification Package

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/index

This code snippet shows how to install the ControlSystemIdentification package using the Julia Package Manager (Pkg). It's the first step to using the system identification functionalities.

```julia
using Pkg
Pkg.add("ControlSystemIdentification")
```

--------------------------------

### Install Optional ControlSystems.jl and Plots Packages

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/index

This snippet details the installation of optional packages, ControlSystemsBase.jl and Plots, which are recommended for working with linear systems and visualization, respectively. These packages build upon the core functionality.

```julia
Pkg.add(["ControlSystemIdentification", "ControlSystemsBase", "Plots"])
```

--------------------------------

### Load and Prepare Glass Furnace Data

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/glass_furnace

Loads furnace data from a URL, downloads and unzips it, then reads the data into a format suitable for system identification. The data is structured into input (u) and output (y) matrices, and then converted into an iddata object.

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

--------------------------------

### Load and Prepare Robot Arm Data

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/flexible_robot

Downloads robot arm data from a URL, decompresses it, and loads it into a format suitable for system identification. It then separates the input (torque) and output (acceleration) and creates an 'iddata' object with an estimated sample time.

```julia
using DelimitedFiles, Plots
using ControlSystemIdentification, ControlSystemsBase

url = "https://ftp.esat.kuleuven.be/pub/SISTA/data/mechanical/robot_arm.dat.gz"zipfilename = "/tmp/flex.dat.gz"
path = Base.download(url, zipfilename)
run(`gunzip -f $path`)
data = readdlm(path[1:end-3])
u = data[:, 1]' # torque
y = data[:, 2]' # acceleration
d = iddata(y, u, 0.01) # sample time not specified for data, 0.01 is a guess
```

--------------------------------

### Visualize Input-Output Data

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/glass_furnace

Generates a plot to visualize the input and output data of the identified system. This helps in understanding the data's characteristics before model estimation.

```julia
plot(d, layout=9)
```

--------------------------------

### Analyze Input-Output Spectra

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/flexible_robot

Generates Welch plots for the input and output data to visualize their power spectra. This helps understand the frequency content of the excitation signal and the system's response, correlating spectral characteristics with coherence observations.

```julia
welchplot(d, yticks=(xt,xt))
```

--------------------------------

### Data Scaling Example for NewPEM

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/api

Demonstrates how to scale identification data before using the `newpem` method to improve estimation accuracy. It involves normalizing the variance of the output and input signals by pre- and post-multiplying with diagonal matrices.

```julia
# Normalize variance
Dy = Diagonal(1 ./ vec(std(d.y, dims=2)))
Du = Diagonal(1 ./ vec(std(d.u, dims=2)))
d̃ = Dy * d * Du
```

--------------------------------

### Julia iddata practical examples and usage

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/iddata

Demonstrates practical usage of iddata function with random data examples showing SISO and MIMO configurations, data manipulation, indexing operations, and accessing object properties like number of inputs/outputs and time vectors.

```julia
julia> iddata(randn(10))
Output data of length 10 with 1 outputs, Ts = nothing

julia> iddata(randn(10), randn(10), 1)
InputOutput data of length 10, 1 outputs, 1 inputs, Ts = 1

julia> d = iddata(randn(2, 10), randn(3, 10), 0.1)
InputOutput data of length 10, 2 outputs, 3 inputs, Ts = 0.1

julia> [d d] # Concatenate along time
InputOutput data of length 20, 2 outputs, 3 inputs, Ts = 0.1

julia> d[1:3]
InputOutput data of length 3, 2 outputs, 3 inputs, Ts = 0.1

julia> d.nu
3

julia> d.t # access time vector
0.0:0.1:0.9
```

--------------------------------

### Analyze Model in Frequency Domain

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/glass_furnace

Visualizes the estimated system model in the frequency domain using a Bode plot. This helps in understanding the system's frequency response characteristics, such as gain and phase across different frequencies.

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

--------------------------------

### Scale input-output data for estimation

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/ss

Example showing how to normalize variance of input and output signals before estimation. The scaling matrices Dy and Du are used to pre and post-multiply the data, with inverse scaling applied to the resulting system.

```julia
Dy = Diagonal(1 ./ vec(std(d.y, dims=2))) # Normalize variance
Du = Diagonal(1 ./ vec(std(d.u, dims=2))) # Normalize variance
d̃ = Dy * d * Du
```

--------------------------------

### Estimate and Validate Glass Furnace Model

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/glass_furnace

Splits the data into training and validation sets, estimates a system model using the subspaceid function, and then visualizes the model's predictions on the validation data for different prediction horizons.

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

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

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

--------------------------------

### Estimate and Compare System Models

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/flexible_robot

Splits the data into training and validation sets. It then estimates two models: one using subspace identification (`subspaceid`) and another using the prediction-error method (`newpem`) with the Nelder-Mead optimizer. The results are visualized using prediction plots.

```julia
dtrain = d[1:end÷2]
dval = d[end÷2:end]

# A model of order 4 is reasonable, a double-mass model. We estimate two models, one using subspace-based identification and one using the prediction-error method
using Optim
model_ss = subspaceid(dtrain, 4, focus=:prediction)
model_pem, x0_pem = newpem(dtrain, 4; focus=:prediction, optimizer=NelderMead(), iterations=100000, show_every=50000)

predplot(model_ss, dval, h=1)
predplot!(model_ss, dval, h=5, ploty=false)
simplot!(model_ss, dval, ploty=false)
```

--------------------------------

### Linear system identification example with noisy data

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/api

Demonstrates system identification using newpem with simulated data from a linear system. Generates random input, simulates system response, adds measurement noise, and estimates a model. Returns a PredictionStateSpace model with system matrices and Kalman filter parameters.

```julia
using ControlSystemIdentification, ControlSystemsBase Plots
G = DemoSystems.doylesat()
T = 1000  # Number of time steps
Ts = 0.01 # Sample time
sys = c2d(G, Ts)
nx = sys.nx
nu = sys.nu
ny = sys.ny
x0 = zeros(nx) # actual initial state
sim(sys, u, x0 = x0) = lsim(sys, u; x0)[1]

σy = 1e-1 # Noise covariance

u  = randn(nu, T)
y  = sim(sys, u, x0)
yn = y .+ σy .* randn.() # Add measurement noise
d  = iddata(yn, u, Ts)

sysh, x0h, opt = ControlSystemIdentification.newpem(d, nx, show_every=10)

plot(
    bodeplot([sys, sysh]),
    predplot(sysh, d, x0h), # Include the estimated initial state in the prediction
)
```

--------------------------------

### Inspect Model's Direct Feedthrough Matrix (D)

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/glass_furnace

Extracts and displays the DDD matrix (direct feedthrough matrix) from the estimated system model. This matrix is important for understanding direct influence of inputs on outputs, particularly relevant when `zeroD=false`.

```julia
model.D
```

--------------------------------

### Estimate impulse response using impulseest and okid in Julia

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/impulse

Demonstrates how to estimate and plot impulse responses for SISO and MIMO systems using the ControlSystemIdentification package. The examples generate random input data, simulate system responses, and visualize both the true and estimated impulse responses. Requires Julia, ControlSystemIdentification, ControlSystemsBase, and Plots packages.

```julia
T = 200
h = 1
t = 0:h:T-h
sys = c2d(tf(1,[1,2*0.1,0.1]),h)

u  = randn(1, length(t))
res = lsim(sys, u, t)
d  = iddata(res)

impulseestplot(d,50, lab=\"Estimate\", seriestype=:steppost)
plot!(impulse(sys,50), lab=\"True system\", l=:dash)
```

```julia
using ControlSystemIdentification, ControlSystemsBase, Plots
T = 200
h = 1
t = 0:h:T-h
sys = ssrand(2,2,4, proper=true, Ts=h)

u  = randn(sys.nu, length(t))
res = lsim(sys, u, t)
d  = iddata(res)

H = okid(d, sys.nx)
plot(impulse(sys,50), lab=\"True system\", layout=sys.ny+sys.nu, sp=(1:4)')
plot!(reshape(H, sys.nu+sys.ny, :)', lab=\"OKID Estiamte\", seriestype=:steppre, l=:dash)
```

--------------------------------

### Visualize Estimated Models in Frequency Domain

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/flexible_robot

Compares the frequency response of the estimated models (PEM and Subspace) with a nonparametric estimate (tfest). Bode plots are used to visualize gain and phase characteristics across a range of frequencies, highlighting how well the models capture the system dynamics.

```julia
w = exp10.(LinRange(-1, log10(pi/d.Ts), 200))
bodeplot(model_pem.sys, w, lab="PEM", plotphase=false, hz=true)
bodeplot!(model_ss.sys, w, lab="Subspace", plotphase=false, hz=true)
plot!(tfest(d), legend=:bottomleft, hz=true, xticks=(xt,xt))
```

--------------------------------

### PEM: System Identification with Prediction Error Method

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/ss

Demonstrates using the `newpem` function from ControlSystemIdentification.jl to estimate a discrete-time LTI system model using the Prediction-Error Method. This example simulates a system, adds noise, and then estimates and predicts using the identified model.

```julia
using ControlSystemIdentification, ControlSystemsBase, Random, LinearAlgebra
using ControlSystemIdentification: newpem
sys = c2d(tf(1, [1, 0.5, 1]) * tf(1, [1, 1]), 0.1)

Random.seed!(1)
T   = 1000                      # Number of time steps
nx  = 3                         # Number of poles in the true system
nu  = 1                         # Number of inputs
x0  = randn(nx)                 # Initial state
sim(sys,u,x0=x0) = lsim(ss(sys), u, x0=x0).y # Helper function
u   = randn(nu,T)               # Generate random input
y   = sim(sys, u, x0)           # Simulate system
y .+= 0.01 .* randn.()          # Add some measurement noise
d   = iddata(y,u,0.1)

sysh,opt = newpem(d, nx, focus=:prediction) # Estimate model

yh = predict(sysh, d) # Predict using estimated model
predplot(sysh, d)     # Plot prediction and true output
```

--------------------------------

### Estimate Model Spectrum with model_spectrum (Julia)

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/api

The `model_spectrum` function estimates the model spectrum using a specified estimation function and sample time. It provides an example using `ar` and `spectrogram`.

```Julia
using ControlSystemIdentification, DSP

T = 1000
s = sin.((1:T) .* 2pi/10)
S1 = spectrogram(s,window=hanning)
estimator = model_spectrum(ar,1,2)
S2 = spectrogram(s,estimator,window=rect)
plot(plot(S1),plot(S2))
```

--------------------------------

### Visualize and Analyze Input-Output Data

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/flexible_robot

Plots the input-output data and its coherence function. The coherence plot helps identify frequency ranges where the output is well-correlated with the input, indicating potential signal-to-noise ratio issues or nonlinearities at low coherence frequencies.

```julia
xt = [2,5,10,20,50]
plot(
    plot(d),
    coherenceplot(d, xticks=(xt,xt), hz=true),
)
```

--------------------------------

### Setup True System and Generate Identification Data in Julia

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/temp

This snippet defines a continuous-time first-order transfer function, discretizes it, simulates output data using a chirp-like input signal, adds measurement noise, and creates an identification data object for model estimation. It requires ControlSystemsBase, ControlSystemIdentification, and Plots packages. Inputs include frequency range and time vector; outputs are the data object ready for identification. Limitation: Assumes linear time-invariant system and may require tuning for different noise levels.

```julia
using ControlSystemsBase, ControlSystemIdentification, Plots
w = 2pi .* exp10.(LinRange(-3, log10(0.5), 500))
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(G0)

u = sign.(sin.((0:0.01:20) .^ 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
plot(data)
```

--------------------------------

### Perform subspace identification with frequency-response data (Julia)

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/freq

Shows how to convert continuous-time frequency data to discrete-time using a bilinear transform and then fit a model with `subspaceid`. The second example generates simulated frequency-response data, fits a model, and visualizes results. Requires ControlSystemIdentification, ControlSystemsBase, and Plots packages.

```Julia
Ts    = 0.01 # Sample time
frd_d = c2d(frd_c::FRD, Ts) # Perform a bilinear transformation to discrete-time frequency vector
Ph, _ = subspaceid(frd_d, Ts, nx)
```

```Julia
using ControlSystemIdentification, ControlSystemsBase, Plots
ny,nu,nx = 2,3,5                        # number of outputs, inputs and states
Ts = 1                                  # Sample time
G = ssrand(ny,nu,nx; Ts, proper=true)   # Generate a random system

N = 200             # Number of frequency points
w = freqvec(Ts, N)
frd = FRD(w, G)     # Build a frequency-response data object (this should in practice come from an experiment)

Gh, x0 = subspaceid(frd, G.Ts, nx; zeroD=true) # Fit frequency response
sigmaplot([G, Gh], w[2:end], lab=[\"True system\" \"Estimated model\"])
```

--------------------------------

### Create ARX model regressor matrix

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/api

Creates a regressor matrix for fitting ARX models of the form A(z)y = B(z)u. Returns a shortened output signal and regressor matrix for least-squares estimation. Supports input delays which modify the effective model order. Includes example usage with signal filtering.

```julia
getARXregressor(y::AbstractVector,u::AbstractVecOrMat, na, nb; inputdelay = ones(Int, size(nb)))
```

--------------------------------

### ERA/OKID: Basic Usage and Multiple Datasets

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/ss

Demonstrates the basic usage of ERA and OKID for system identification. It also shows how to utilize multiple datasets to improve estimation accuracy by averaging impulse responses and then estimating a model, or by passing a vector of datasets directly to ERA.

```julia
using ControlSystemIdentification, ControlSystemsBase, Plots, Statistics
Ts = 0.1
G = c2d(tf(1, [1,1,1]), Ts) # True system

# Create several "experiments"
ds = map(1:5) do i
    u = randn(1, 1000)
    y, t, x = lsim(G, u)
    yn = y + 0.2randn(size(y))
    iddata(yn, u, Ts)
end

Ys = okid.(ds, 2, round(Int, 10/Ts), smooth=true, λ=1)    # Estimate impulse response for each experiment
Y = mean(Ys)                            # Average all impulse responses

imp = impulse(G, 10)
f1 = plot(imp, lab="True", l=5)
plot!(imp.t, vec.(Ys), lab="Individual estimates", title="Impulse-response estimates")
plot!(imp.t, vec(Y), l=(3, :black), lab="Mean")

models = era.(Ys, Ts, 2, 50, 50)    # estimate models based on individual experiments
meanmodel = era(Y, Ts, 2, 50, 50)   # estimate model based on mean impulse response

f2 = bodeplot([G, meanmodel], lab=["True" "" "Combined estimate" ""], l=2)
bodeplot!(models, lab="Individual estimates", c=:black, alpha=0.5, legend=:bottomleft)

plot(f1, f2)
```

```julia
era(ds, 2, 50, 50, round(Int, 10/Ts), p=1, λ=1, smooth=true) # Should be identical to meanmodel above
```

--------------------------------

### Load and Prepare Evaporator Data in Julia

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/evaporator

This snippet downloads and prepares the evaporator dataset from a remote source. It reads the data, defines inputs and outputs, and creates an iddata object for system identification.

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

--------------------------------

### One-Step Ahead Prediction Example in Julia

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/validation

This example demonstrates one-step ahead prediction using the predict function on a transfer function and input-output data. The output is a vector of predicted values, shorter by the maximum lag due to initial conditions. Limited to ARX models; requires iddata for input data structure.

```Julia
predict(tf(1, [1, -1], 1), iddata(1:10, 1:10))
```

--------------------------------

### Kautz Basis Function

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/freq

This snippet appears to be the start of a function definition for generating Kautz basis functions, which are often used in system identification and signal processing for function approximation.

```julia
function kautz(a::AbstractVector)

```

--------------------------------

### Simulate Models and Compare (Julia)

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/delayest

This snippet simulates the estimated models and compares their performance with the original data.

```Julia
plot(res, lab="Data")

```

```Julia
simplot!(model_arx, d, ploty=false, sysname="ARX model")

```

```Julia
simplot!(model_arx_2, d, ploty=false, sysname="ARX model with delay")

```

--------------------------------

### Create PI controller and closed-loop system, plot Nyquist and step response

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/delayest

Defines a PI controller using pid, forms the loop with plant P, obtains the closed-loop feedback system, and visualizes the Nyquist plot of the loop and the step response of the closed-loop for 50 samples. Requires ControlSystems.jl and Plots.

```Julia
C = pid(0.1, 0.5; Ts)
L = P*C
G = feedback(L)
plot(nyquistplot(L), plot(step(G, 50)))
```

--------------------------------

### Control System Identification Initialization with Options

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/ss

Initializes the `newpem` function for system identification with various keyword arguments to control optimization, focus, regularization, and model stability. Dependencies include Optim.jl for optimization and automatic differentiation for gradient calculation. The default solver is BFGS().

```julia
newpem(..., sys0=initial_guess, focus=:prediction, regularizer=reg_val, stable=true, zeroD=true, h=horizon_value)
```

--------------------------------

### Estimate Model Parameters with nonlinear_pem

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/modelingtoolkit

This snippet shows how to prepare input/output data, define initial guesses for parameters and states, set up noise covariance matrices (R1 for process noise, R2 for measurement noise), and finally call `nonlinear_pem` to estimate the model parameters. The `iddata` object bundles the data, and `Diagonal` creates covariance matrices.

```julia
d = iddata(Y, U, Ts)
x0_guess = [2.5, 1.5, 1, 2]     # Guess for the initial condition (initial state)
p_guess  = [1.4, 1.4, 5.1, 0.25] # Initial guess for the parameters

R1 = Diagonal([0.1, 0.1, 0.1, 0.1]) # This controls how much we trust the model (covariance of the process noise)
R2 = 100*Diagonal(0.5^2 * ones(ny)) # This controls how much we trust the measurements (covariance of the measurement noise)

model = ControlSystemIdentification.nonlinear_pem(d, discrete_dynamics, measurement, p_guess, x0_guess, R1, R2, nu)
```

--------------------------------

### Obtain Dynamics Functions with MTK

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/modelingtoolkit

Generates continuous and discrete dynamics functions from an MTK model. It handles input/output collection, parameter tuning, and creates a wrapper for efficient computation. Dependencies include ModelingToolkit and SeeToDee.

```julia
mtkmodel = complete(mtkmodel)
inputs = [collect(mtkmodel.u);]
outputs = [collect(mtkmodel.h[1:2]);]
tunable_p = [mtkmodel.k1, mtkmodel.k2, mtkmodel.A, mtkmodel.γ] # Provided in the same order as p_guess

function get_mtk_dynamics(mtkmodel, inputs, outputs, tunable_p) # A wrapper function to avoid using global variables

    mtkmodel = mtkcompile(mtkmodel; inputs, outputs, split=false)
    (f_oop, f_ip), statevars, p, io_sys = ModelingToolkit.generate_control_function(mtkmodel)

    continuous_dynamics = f_oop # This is ẋ = f(x, u, p, t)
    inner_discrete_dynamics = SeeToDee.Rk4(continuous_dynamics, Ts::Float64) # x⁺ = f(x, u, p, t)
    tunable_indices = [findfirst(isequal(pi), p) for pi in tunable_p] # Figure out what indices of the parameter array correspond to our tunable parameters
    p0 = [ModelingToolkit.defaults(io_sys)[pi] for pi in p]
    full_p = deepcopy(p0)
    output_indices = [findfirst(isequal(yi), statevars) for yi in outputs] # Figure out what indices of the state array correspond to our outputs

    function discrete_dynamics_wrapper(x, u, p, t)
        opt_p = similar(p, length(full_p)) # Create an array of correct length and element type to host the full parameter vector
        opt_p .= full_p # Write all the initial parameters into this new array
        opt_p[tunable_indices] .= p # Overwrite the tunable parameters with the optimization variable
        inner_discrete_dynamics(x, u, opt_p, t) # Call the discretized dynamics function from MTK
    end

    mtk_measurement(x,u,p,t) = x[output_indices]

    discrete_dynamics_wrapper, mtk_measurement
end

Ts = 1.0    # Sample interval used for discretization
discrete_dynamics, measurement = get_mtk_dynamics(mtkmodel, inputs, outputs, tunable_p)

```

--------------------------------

### Generate Dataset for Identification (Julia)

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/delayest

This snippet generates a dataset using `lsim` to simulate the discrete-time system with a delay.

```Julia
u = sign.(repeat(randn(100), inner=5))

```

```Julia
res = lsim(Pd, u')

```

```Julia
plot(res, plotu=true)

```

--------------------------------

### Validate Evaporator Model in Julia

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/evaporator

Generates prediction plots for 1-step and 5-step ahead predictions on the validation data, and compares the model's performance with results from the original paper.

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

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

--------------------------------

### Compare Estimated Parameters

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/nonlinear

Compares the true parameters, initial guess parameters, and the estimated parameters obtained from the nonlinear PEM. This visualizes the accuracy of the parameter estimation.

```julia
[p_true p_guess model.p]
```

--------------------------------

### Model Simulation and Comparison

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/nonlinear

Simulates the estimated model and compares its output with the original data and a model based on an initial guess. This helps in evaluating the performance of the estimated model.

```julia
simplot(model, d, layout=2)

x_guess = LowLevelParticleFilters.rollout(discrete_dynamics, x0_guess, u, p_guess)[1:end-1]
y_guess = measurement.(x_guess, u, 0, 0)
plot!(reduce(hcat, y_guess)', lab="Initial guess")
```

--------------------------------

### Estimate statespace model using n4sid in Julia

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/ss

Shows alternative subspace-based identification using the n4sid function on the same dataset. This older implementation produces comparable results to subspaceid. The example demonstrates basic usage and plots the estimated model for comparison with the true system and subspaceid results.

```julia
sys2 = n4sid(d, :auto; verbose=false, zeroD=true)
bodeplot!(sys2.sys, lab=["n4sid" ""])
```

--------------------------------

### Load and Prepare Data for Belt-Drive System Identification

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/hammerstein_wiener

Loads data from a specified URL, unzips it, and prepares it into a format suitable for system identification. It creates `iddata` objects for input-output pairs and plots them. Dependencies include `DelimitedFiles`, `Plots`, `ControlSystemIdentification`, and `ControlSystemsBase`.

```julia
using DelimitedFiles, Plots
using ControlSystemIdentification, ControlSystemsBase

url = "https://drive.google.com/uc?export=download&id=10Z83_2qHTSJ_F9usqTsPohjZ1xtXCtE3"
zipfilename = "/tmp/bd.zip"
cd("/tmp")
path = Base.download(url, zipfilename)
run(`unzip -o $path`)
data = readdlm("/tmp/DATAUNIF.csv", ',')[2:end, 1:4]
iddatas = map(0:1) do ind
    u = data[:, 1 + ind]' .|> Float64 # input
    y = data[:, 3 + ind]' .|> Float64 # output
    iddata(y, u, 1/50)
end

plot(plot.(iddatas)...)

d = iddatas[1] # We use one dataset for estimation
coherenceplot(d)
```

--------------------------------

### Simulate Subspace Identification Models (Julia)

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/delayest

This snippet simulates the subspace identification models and compares their performance with the original data.

```Julia
plot(res, lab="Data")

```

```Julia
simplot!(model_sim, d, ploty=false, sysname="Simulation model")

```

```Julia
simplot!(model_pred, d, ploty=false, sysname="Prediction model")

```

--------------------------------

### Evaluate Model Prediction and Simulation Performance

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/flexible_robot

Quantifies the prediction and simulation errors of the estimated models for various prediction horizons. It calculates the Root Mean Square (RMS) error for multi-step predictions and open-loop simulation, comparing the performance of the PEM and subspace-identified models.

```julia
using Statistics
hs = [1:40; 45:5:80]
perrs_pem = map(hs) do h
    yh = predict(model_pem, d, x0_pem; h)
    ControlSystemIdentification.rms(d.y - yh) |> mean
end
perrs_ss = map(hs) do h
    yh = predict(model_ss, d; h)
    ControlSystemIdentification.rms(d.y - yh) |> mean
end

serr_pem = ControlSystemIdentification.rms(d.y - simulate(model_pem, d)) |> mean
serr_ss = ControlSystemIdentification.rms(d.y - simulate(model_ss, d)) |> mean

plot(hs, perrs_pem, lab="Prediction errors PEM", xlabel="Prediction Horizon", ylabel="RMS error")
plot!(hs, perrs_ss, lab="Prediction errors Subspace")
hline!([serr_pem], lab="Simulation error PEM", l=:dash, c=1, ylims=(0, Inf))
hline!([serr_ss], lab="Simulation error Subspace", l=:dash, c=2, legend=:bottomright, ylims=(0, Inf))
```

--------------------------------

### Create Discrete-Time System with Delay (Julia)

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/delayest

This snippet converts the continuous-time system to a discrete-time equivalent with a one-sample delay using `c2d` and visualizes the result.

```Julia
Ts = 0.01

```

```Julia
Pd = c2d(P*delay(Ts), Ts)

```

```Julia
bodeplot!(bp, Pd, lab="Discrete-time system with delay")

```

--------------------------------

### Download and preprocess ball and beam data in Julia

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/ballandbeam

Downloads the ball and beam dataset from STADIUS's database, decompresses it, and loads it into an iddata object for system identification. The input is the beam angle and the output is the ball position.

```julia
using DelimitedFiles, Plots
using ControlSystemIdentification, ControlSystemsBase

url = "https://ftp.esat.kuleuven.be/pub/SISTA/data/mechanical/ballbeam.dat.gz"
zipfilename = "/tmp/bb.dat.gz"
path = Base.download(url, zipfilename)
run(`gunzip -f $path`)
data = readdlm(path[1:end-3])
u = data[:, 1]' # beam angle
y = data[:, 2]' # ball position
d = iddata(y, u, 0.1)
```

--------------------------------

### Nonlinear Prediction-Error Method parameter estimation

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/nonlinear

Performs parameter estimation using the prediction-error method with Unscented Kalman Filtering. This function minimizes one-step prediction errors to estimate model parameters and initial conditions. Requires LeastSquaresOptim.jl to be manually installed and loaded. The method works well with poor initial guesses, unstable systems, and chaotic systems.

```julia
nonlinear_pem(
    d::IdData,
    discrete_dynamics,
    measurement,
    p0,
    x0,
    R1,
    R2,
    nu;
    optimizer = LevenbergMarquardt(),
    λ = 1.0,
    optimize_x0 = true,
    kwargs...,
)
```

--------------------------------

### estimate_x0

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/api

Estimates the initial state of a system given input-output data.

```APIDOC
## estimate_x0

### Description
Estimate the initial state of the system.

### Method
Function (specific HTTP method not applicable for library functions)

### Endpoint
N/A (Function within ControlSystemIdentification class)

### Parameters
#### Path Parameters
None

#### Query Parameters
None

#### Request Body
None

### Arguments
  * `sys` : System object
  * `d` : `iddata` object containing input-output data.
  * `n`: Number of samples to use for estimation.
  * `fixed`: Optional vector to specify fixed initial state values.

### Request Example
```julia
# Assuming sys is a system object and d is iddata
sys   = ssrand(2,3,4, Ts=1)
x0    = randn(sys.nx)
u     = randn(sys.nu, 100)
y,t,x = lsim(sys, u; x0)
d     = iddata(y, u, 1)
x0h   = estimate_x0(sys, d, 8, fixed=[Inf, x0[2], Inf, Inf])
```

### Response
#### Success Response (200)
- **x0h** (Array) - Estimated initial state vector.

#### Response Example
```json
{
  "x0h": [0.5, -1.2, 0.8, ...]
}
```
```

--------------------------------

### Plot Evaporator Data in Julia

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/evaporator

Visualizes the input-output data of the evaporator system using a 6-panel plot layout to inspect the data before model estimation.

```julia
plot(d, layout=6)
```

--------------------------------

### Define Quad‑tank Model using ModelingToolkit in Julia

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/modelingtoolkit

Creates a ModelingToolkit model with default parameters, variables, and differential equations for a quad‑tank system. Requires the ControlSystemIdentification, ModelingToolkit, and related packages. The model is instantiated as `mtkmodel` for later use in parameter estimation.

```julia
using ControlSystemIdentification, ModelingToolkit, LeastSquaresOptim, SeeToDee, LowLevelParticleFilters, LinearAlgebra, Random, Plots # Load the
using ModelingToolkit: D_nounits as D
t = ModelingToolkit.t_nounits
ssqrt(x) = √(max(x, zero(x)) + 1e-3) # For numerical robustness at x = 0
@register_symbolic ssqrt(x)

@mtkmodel QuadtankModel begin
    @parameters begin
        k1 = 1.4
        k2 = 1.4
        g = 9.81
        A = 5.1
        a = 0.03
        γ = 0.25
    end
    begin
        A1 = A2 = A3 = A4 = A
        a1 = a3 = a2 = a4 = a
        γ1 = γ2 = γ
    end
    @variables begin
        h(t)[1:4] = 0
        u(t)[1:2] = 0
    end
    @equations begin
        D(h[1]) ~ -a1/A1 * ssqrt(2g*h[1]) + a3/A1*ssqrt(2g*h[3]) +     γ1*k1/A1 * u[1]
        D(h[2]) ~ -a2/A2 * ssqrt(2g*h[2]) + a4/A2*ssqrt(2g*h[4]) +     γ2*k2/A2 * u[2]
        D(h[3]) ~ -a3/A3*ssqrt(2g*h[3])                          + (1-γ2)*k2/A3 * u[2]
        D(h[4]) ~ -a4/A4*ssqrt(2g*h[4])                          + (1-γ1)*k1/A4 * u[1]
    end
end

@named mtkmodel = QuadtankModel()
```

--------------------------------

### Frequency Domain Analysis in Julia

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/evaporator

Performs frequency domain analysis of the estimated model by plotting the singular values across a range of frequencies.

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

--------------------------------

### Benchmark nonlinear_pem Estimation Speed

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/modelingtoolkit

This snippet uses the BenchmarkTools.jl package to measure the execution time and memory allocations of the `nonlinear_pem` function. This is useful for understanding the performance characteristics of the estimation process.

```julia
using BenchmarkTools
@btime ControlSystemIdentification.nonlinear_pem(d, discrete_dynamics, measurement, p_guess, x0_guess, R1, R2, nu)
```

--------------------------------

### Create Continuous-Time System with Delay (Julia)

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/delayest

This snippet creates a continuous-time double-mass system and visualizes its frequency response using `bodeplot` from ControlSystems.jl.

```Julia
P = DemoSystems.double_mass_model(outputs=2) |> minreal

```

```Julia
bp = bodeplot(P, lab="Continuous-time system without delay", legend=:bottomleft)

```

--------------------------------

### Estimate multiple models for echo data and overlay simulation plots

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/delayest

Fits an ARX model with input delay, a high‑order subspace model, and a low‑order subspace model with extended prediction horizon, then overlays their simulation results for visual comparison. Demonstrates trade‑offs between model complexity and fit quality.

```Julia
model3 = arx(decho, 4, 4, inputdelay=τ_samples)
model4 = subspaceid(decho, 24)
model5 = subspaceid(decho, 4, r=50)
figsim = simplot(model3, decho, zeros(ss(model3).nx), sysname="ARX 4")
simplot!(model4, decho, zeros(model4.nx), ploty=false, plotu=false, sysname="Subspace 24")
simplot!(model5, decho, zeros(model5.nx), ploty=false, plotu=false, sysname="Subspace 4")
```

--------------------------------

### Visualize Model Performance with simplot

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/modelingtoolkit

This snippet demonstrates how to visualize the performance of the estimated model by simulating it with the `simplot` function and comparing the simulated output against the measured data `d`. The `layout` parameter controls the arrangement of plots.

```julia
simplot(model, d, layout=2)
```

--------------------------------

### Generate reference signal, simulate closed-loop, create iddata, plot signals and cross-correlation

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/delayest

Creates a sign‑sine reference signal, simulates the closed-loop system G, converts the output to identification data, and plots both the time series and the cross‑correlation to assess delay effects. Uses ControlSystems.jl identification utilities.

```Julia
ref = sign.(sin.(0.02 .* (0:Ts:50).^2)) # An interesting reference signal
res = lsim(G, ref')
dG = iddata(res)
plot(plot(dG), crosscorplot(dG))
```

--------------------------------

### Generate delayed system and simulate response in Julia

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/delayest

Creates a discrete‑time second‑order system with an input delay, generates a sinusoidal input, simulates the system response, and converts the result into identification data for plotting. Requires ControlSystemIdentification, ControlSystemsBase, and Plots packages.

```Julia
using DelimitedFiles, Plots
using ControlSystemIdentification, ControlSystemsBase

τ  = 2.0 # Delay in seconds
Ts = 0.1 # Sampling time
P  = c2d(tf(1, [1, 0.5, 1])*delay(τ), Ts) # Dynamics is given by a simple second-order system with input delay

u   = sin.(0.1 .* (0:Ts:30).^2) # An interesting input signal
res = lsim(P, u', x0=[0.5; 0; zeros(20)])
d   = iddata(res)
plot(d)
```

--------------------------------

### Pendulum Dynamics Simulation in Julia

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/unstable_systems

This snippet defines a pendulum dynamics model and demonstrates its simulation in Julia using the `LowLevelParticleFilters` package. It simulates the pendulum's motion and generates data for parameter estimation.

```Julia
using Plots, Statistics, DataInterpolations, LowLevelParticleFilters

Ts = 0.01 # Sample time
tsteps = range(0, stop=20, step=Ts)

x0 = [0.0, 3.0] # Initial angle and angular velocity

function pendulum(x, u, p, t) # Pendulum dynamics
    g = 9.82 # Gravitational constant
    L = p isa Number ? p : p[1] # Length of the pendulum
    gL = g / L
    θ = x[1]
    dθ = x[2]
    [dθ
     -gL * sin(θ)]
end
```

```Julia
using LowLevelParticleFilters: rk4, rollout

discrete_pendulum = rk4(pendulum, Ts) # Discretize the continuous-time dynamics using RK4

function simulate(fun, p)
    x = rollout(fun, x0, tsteps, p; Ts)[1:end-1]
    y = first.(x) # This is the data we have available for parameter estimation (angle measurement)
    x, y
end

x, y = simulate(discrete_pendulum, 1.0) # Simulate with L = 1.0
plot(tsteps, y, title = "Pendulum simulation", label = "angle")
```

```Julia
function simloss(p)
    x,yh = simulate(discrete_pendulum, p)
    yh .= abs2.(y .- yh)
    return mean(yh)
end
```

```Julia
Ls = 0.01:0.01:2
simlosses = simloss.(Ls)
fig_loss = plot(Ls, simlosses,
    title  = "Loss landscape",
    xlabel = "Pendulum length",
    ylabel = "MSE loss",
    lab    = "Simulation loss"
)
```

--------------------------------

### Create echo chamber system, discretize, simulate, and plot various responses

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/delayest

Constructs an echo‑chamber system with a resonant channel and delayed feedback, discretizes it, simulates with a reference signal, creates identification data, and plots Bode, pole‑zero, impulse response, and raw data. Highlights challenges of internal delays.

```Julia
ref = sign.(sin.(0.02 .* (0:Ts:100).^2)) # An interesting reference signal
Pc = feedback(tf(100, [1, 10, 100]), -tf(60, [1, 10, 100])*delay(τ)) # Feed 60% of the output back at the input with a delay of 2 seconds (like an echo)
Pd = c2d(Pc, Ts)
res = lsim(Pd, ref')
decho = iddata(res)
plot(bodeplot(Pd, lab=""), pzmap(Pd), plot(impulse(Pd, 10), title="Impulse response"), plot(decho))
```

--------------------------------

### Compare PEM and Fourier-based models in Julia

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/examples/ballandbeam

Compares the PEM-estimated model with a nonparametric Fourier-based estimate in the frequency domain. Shows differences in low-frequency estimation due to system instability.

```julia
w = exp10.(LinRange(-1.5, log10(pi/d.Ts), 200))
bodeplot(model.sys, w, lab="PEM", plotphase=false)
plot!(tfest(d))
```

--------------------------------

### Simulate System with LowLevelParticleFilters.simulate

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/api

Simulates a system's response given the system model, input, and optionally an initial state.

```julia
simulate(sys, u, x0 = nothing)
simulate(sys, d, x0 = nothing)
```

--------------------------------

### Create noise model in Julia

Source: https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/ss

Returns a model of the noise driving the system in innovation form. The model neglects input u and represents the noise dynamics. Internally calls ControlSystemsBase.innovation_form.

```julia
noise_model(sys::AbstractPredictionStateSpace)
```