ARMAX System Identification

System identification for two physical plants — a tank filling valve (SISO) and a KUKA youBot arm (MIMO) — using hand-implemented Gauss-Newton and Gradient Descent optimizers, benchmarked against MATLAB's built-in armax().

---

What it does

Given input/output measurements from a real or simulated system, the goal is to find a discrete-time model that explains how the plant responds. The model structure chosen here is ARMAX — it captures system dynamics, accounts for external excitation, and explicitly models colored output noise.

Two plants are studied:

Plant	Type	Input	Output	Sampling
Tank filling valve	SISO	Valve voltage (V)	Water level (cm)	1 s
KUKA youBot arm (5 DOF)	MIMO	Joint commands (rad)	Joint angles (rad)	100 ms

Data for the tank was collected from a Factory IO + TIA Portal setup over Modbus TCP. Data for the KUKA was collected from a Webots simulation using PRBS excitation via webot_sim/youbot.c.

---

The Model — ARMAX(na, nb, nc)

The ARMAX model in the backward-shift operator q^-1:

A(q) * y(k) = B(q) * u(k) + C(q) * e(k)

where the polynomials are:

A(q) = 1 + a1*q^-1 + a2*q^-2         (autoregressive part)
B(q) =     b1*q^-1 + b2*q^-2         (exogenous input part)
C(q) = 1 + c1*q^-1                   (moving average / noise part)

Expanded into a one-step-ahead predictor:

y_hat(k) = -a1*y(k-1) - a2*y(k-2)
          + b1*u(k-1) + b2*u(k-2)
          + c1*e(k-1)

This is linear in the parameters, so it can be written as:

y_hat(k) = phi(k)^T * theta

where the regressor vector is:

phi(k) = [-y(k-1), -y(k-2), u(k-1), u(k-2), e(k-1)]^T

and the parameter vector is:

theta = [a1, a2, b1, b2, c1]^T

The prediction error at time k:

e(k) = y(k) - phi(k)^T * theta

MIMO Extension (KUKA)

For the 5-joint arm, y(k) and u(k) are both vectors in R^5. The scalar polynomials become matrix polynomials:

A(q) = I + A1*q^-1 + A2*q^-2      (5x5 matrices)
B(q) =     B1*q^-1 + B2*q^-2      (5x5 matrices)
C(q) = I + C1*q^-1                (5x5 matrix)

Each output j has its own parameter column. The regressor for all outputs at time k:

phi(k) = [-y(k-1); -y(k-2); u(k-1); u(k-2); e(k-1)]   in R^25

Parameters are stored in Theta (25x5), so y_hat(k)^T = phi(k)^T * Theta. Total: 125 parameters.

---

Cost Function

All three methods minimize the sum of squared prediction errors over the training set:

J(theta) = sum_k e(k)^2 = sum_k [y(k) - phi(k)^T * theta]^2

In matrix form, stacking N samples into Phi (NxNp) and E (Nx1):

J(theta) = ||E||^2 = ||Y - Phi*theta||^2

The gradient is:

grad J(theta) = -2 * Phi^T * E

The Gauss-Newton approximation to the Hessian is H ~= 2 * Phi^T * Phi, which avoids computing second derivatives.

---

Optimization Methods

1. Gauss-Newton + Armijo Line Search

The Gauss-Newton step solves the normal equations at each iteration:

Phi^T * Phi * delta = Phi^T * E

The step delta = (Phi^T * Phi)^-1 * Phi^T * E is the Newton direction. To guarantee cost reduction, a backtracking line search enforces the Armijo (sufficient decrease) condition:

J(theta + alpha*delta) <= J(theta) - c * alpha * ||delta||^2

Starting from alpha = 1, the step is halved (alpha <- rho*alpha, rho = 0.5) until the condition holds. Convergence is declared when alpha * ||delta|| < tol.

For the tank script, Gauss-Newton uses Levenberg-Marquardt regularization instead of pure Armijo — the normal equations become:

(Phi^T * Phi + lambda*I) * delta = Phi^T * E

lambda is initialized as lambda0 = 1e-3 * trace(Phi^T*Phi) / Np and adapts: if the step reduces cost, lambda <- lambda/5; otherwise lambda <- 5*lambda. This blends between Gauss-Newton (small lambda) and gradient descent (large lambda), giving robust convergence from a cold start.

2. Gradient Descent + Armijo (Normalized Direction)

The descent direction is the normalized negative gradient:

d = -grad J / ||grad J||

Normalizing removes the scale sensitivity of raw gradient descent. The Armijo condition becomes:

J(theta + alpha*d) <= J(theta) + c * alpha * grad J^T * d

Since d is normalized, grad J^T * d = -||grad J||, which is always negative — the condition guarantees descent. Convergence is much slower than Gauss-Newton (no curvature information), but it serves as a useful baseline.

3. MATLAB built-in armax()

MATLAB's armax() uses a Prediction Error Method (PEM) internally. It is run per-joint in SISO mode for the KUKA arm and used purely as a reference benchmark.

---

Extended Least Squares (ELS) Regressor

The noise term e(k-1) in the regressor phi(k) creates a dependency problem: computing phi(k) requires the prediction error at k-1, which depends on theta. This is resolved by maintaining a persistent error buffer updated sample-by-sample as the forward pass proceeds — a scheme known as Extended Least Squares (ELS):

for k = n0 to N:
    phi(k) = [-y(k-1), -y(k-2), u(k-1), u(k-2), e_hat(k-1)]
    e_hat(k) = y(k) - phi(k)^T * theta

The error buffer e_hat is carried across iterations rather than re-zeroed, keeping the MA part properly informed throughout optimization.

---

Transfer Function and Stability

After identification, the Gauss-Newton parameters build the discrete transfer function:

          b1*z + b2
G(z) =  ─────────────────
          z^2 + a1*z + a2

For KUKA (MIMO), this becomes a 5x5 transfer function matrix G(z) where each output row shares the same denominator polynomial across all input columns.

Stability is verified by checking all poles lie inside the unit circle:

|z_i| < 1   for all roots of A(z) = 0

The continuous-time model G(s) is recovered via zero-order hold (ZOH) using d2c. The dominant pole gives the time constant tau = -1/Re(s*) and settling time ~= 4*tau.

---

Fit Metric

Model quality is reported as the NRMSE fit percentage (same convention as MATLAB's System Identification Toolbox):

Fit% = 100 * (1 - ||y - y_hat|| / ||y - mean(y)||)

A fit of 100% means perfect prediction; 0% means the model is no better than a constant.

---

Project Structure

armax-identification/
├── mat_scripts/
│   ├── tank.m              — SISO identification: tank valve -> water level
│   ├── kuka.m              — MIMO identification: joint commands -> joint angles
│   └── results_kuka/       — saved figures and CSVs from kuka.m
│
├── sim_models/
│   ├── armax_system.slx    — Simulink model: real-time identifier in the loop
│   ├── setup_workspace.m   — loads CSV data into Simulink workspace
│   ├── identifier.m        — SISO Gauss-Newton identifier (Simulink MATLAB Function block)
│   └── identifier_mimo.m   — MIMO Gauss-Newton identifier (Simulink MATLAB Function block)
│
├── webot_sim/
│   ├── youbot.c            — Webots C controller: PRBS excitation + CSV logging
│   └── Makefile
│
├── data_kuka/              — joint{1..5}_data.csv  (k, u_cmd, y_meas)
└── data_tank/              — Tank_Data.csv, fill*.csv, online_collected.csv

---

Running

Offline identification

% From mat_scripts/
run('tank.m')    % SISO — tank valve
run('kuka.m')    % MIMO — KUKA arm

Both scripts produce figures and save results to results_tank/ or results_kuka/.

Simulink real-time identifier

% From sim_models/
mode = 'tank';   % or 'kuka'
run('setup_workspace.m')
sim('armax_system')

The Simulink model calls identifier.m (SISO) or identifier_mimo.m (MIMO) as MATLAB Function blocks, updating theta every simulation step.

Webots data collection

Build and run youbot.c inside Webots. In OFFLINE mode the controller applies PRBS excitation and writes joint{1..5}_data.csv directly. In ONLINE mode, data is buffered and sent over TCP to MATLAB on keypress S.

---

Dependencies

• MATLAB R2021a or later
• System Identification Toolbox — armax(), iddata(), predict()
• Control System Toolbox — tf(), step(), d2c(), pole()
• Webots R2023a or later (simulation only)
• Factory IO + TIA Portal + Modbus TCP (tank online mode only)