3. Identification — Five Methods

☀️ Light

The goal of identification

Identification is an optimisation problem. Given a parametric model \(\mathcal{M}(\mathbf{p})\) producing a predicted trajectory \(\hat\omega(t; \mathbf{p})\) for a given input \(u(t)\), find parameters \(\mathbf{p}\) that minimise the error against the measurement:

$$ \mathbf{p}^{\star} = \arg\min_{\mathbf{p}}; \frac{1}{N}\sum_{n=1}^{N}\left(\hat\omega(t_n; \mathbf{p}) - \omega_{\text{meas}}(t_n)\right)^{2} $$

That’s it. Identification is curve fitting in the state-space sense. Everything beyond that is about how you pick the structure of \(\mathcal{M}\) and what solver you use for the minimisation.

Black box vs. grey box

Two schools of thought:

For the Picooz we go grey box — we know the physics of a DC motor and a rotor with aerodynamic drag, we just don’t know tau, k, k2. Grey-box needs far less data, has parameters with physical meaning (so a change in rotor -> you know which value to re-fit), and yields a compact model that embeds well on a low-cost MCU.

The Picooz grey-box model

A brushed DC motor driving a rotor with non-linear aerodynamic drag obeys:

$$ \frac{d\omega}{dt} = -\frac{1}{\tau}, e^{k_{2}\omega} + k, u \qquad (\omega > 0) $$

with \(\omega = 0\) when at rest. Three unknown parameters:

• \(\tau\) — mechanical time constant [s], dominated by rotor inertia + nominal damping.
• \(k_2\) — exponential growth coefficient of the aerodynamic drag with speed.
• \(k\) — input gain [rotor-speed units per PWM duty].

The exponential drag term is what makes this non-linear: a linear first-order model would mis-predict the behaviour at high speed (where the drag grows faster than linearly).

The MATLAB implementation is a three-line ODE function shipped with the tutorial:

function [dx, y] = picooz_motorBlade_m(t, x, u, tau, tau2, k, k2, k3, varargin)
  x(1) = max(x(1), 0);            % Speed always positive
  y    = x(1);                    % Output = speed
  if x(1) > 0
    dx(1) = -(1/tau) * exp(k2*x(1)) + k*u(1);
  else
    dx(1) =  k*u(1);
  end
end

Five ways to fit the same model

We compare five methods in order of increasing sophistication.

Key takeaway upfront: L2, L3 and L4 all produce essentially the same parameter values on this dataset. The question is therefore only one of convenience, speed, and which toolboxes you happen to own — not which method is “more accurate”.

Level 0 — Paper / pencil

Before opening MATLAB, do a dimensional analysis. A brushed-DC-motor-plus-rotor settles to steady-state speed at some duty \(u_0\). From a step response:

• The time to reach 63 % of the steady-state speed gives a rough estimate of \(\tau\).
• The steady-state speed divided by the duty gives a rough estimate of \(k\).

This is not optimisation — it’s the sanity check that tells you whether \(\tau \sim 0.1,\text{s}\) or \(\tau \sim 10,\text{s}\) before the solver wastes iterations going to the wrong basin.

You will always do this first on paper. Every time.

Level 1 — Linear regression with \ (base MATLAB)

Simplification: ignore the exponential drag. Model the rotor as a linear first-order system around its nominal operating point:

$$ \frac{d\omega}{dt} \approx -\frac{\omega}{\tau_{\text{lin}}} + k,u $$

Discretised by the backward Euler rule at sample period \(T_s\):

$$ \omega[n] = a,\omega[n-1] + b,u[n-1] \quad\text{with}\quad a = 1 - T_s/\tau_{\text{lin}},; b = k T_s $$

This is linear in \((a, b)\). Build the regression matrix from the chirp data and solve in closed form:

$$ \begin{bmatrix}\omega[1]\\omega[2]\\vdots\\omega[N]\end{bmatrix} ;=; \underbrace{\begin{bmatrix}\omega[0] & u[0]\ \omega[1] & u[1]\ \vdots & \vdots \ \omega[N-1] & u[N-1]\end{bmatrix}}_{A} \begin{bmatrix}a\b\end{bmatrix} \quad\Rightarrow\quad \begin{bmatrix}a\b\end{bmatrix} = A \backslash Y $$

load tutorial_picooz_chirp.mat          % provides t, u, omega_meas
Ts = mean(diff(t));
Y  = omega_meas(2:end);
A  = [omega_meas(1:end-1)  u_abs(1:end-1)];
p  = A \ Y;
tau = Ts / (1 - p(1));
k   = p(2) / Ts;

10 lines of code, ~30 ms runtime, no non-linear solver.

Pedagogic point. On the Picooz this linear model is a poor fit — try it. The exponential drag dominates at high speed and the best linear tau is physically meaningless. That is exactly why we need the next levels. Showing a method that fails is as pedagogically valuable as one that works.

Level 2 — Non-linear fminsearch (base MATLAB)

The full non-linear model. No small-signal approximation. Still no add-on toolbox — fminsearch (Nelder-Mead simplex) ships with base MATLAB.

Recipe:

1. Wrap the grey-box ODE picooz_motorBlade_m.m in a simulator that returns \(\hat\omega(t)\) for a given parameter vector.
2. Define the cost as the mean-squared error against omega_meas.
3. Call fminsearch.

load tutorial_picooz_chirp.mat

% Cost = MSE between simulation and measurement
cost = @(p) sum((simulate_picooz(p, t, u_abs) - omega_meas).^2) / numel(t);

% Initial guess from L0 pencil work
p0 = [0.10;   % tau
      0.01;   % k2
      1.00];  % k

opts  = optimset('Display','iter','TolFun',1e-8,'MaxIter',400);
p_fit = fminsearch(cost, p0, opts);

function omega_sim = simulate_picooz(p, t, u)
    tau = p(1); k2 = p(2); k = p(3);
    uinterp = @(tt) interp1(t, u, tt, 'linear', 'extrap');
    function dx = ode(tt, x)
        x = max(x,0);
        if x > 0
            dx = -(1/tau) * exp(k2*x) + k * uinterp(tt);
        else
            dx = k * uinterp(tt);
        end
    end
    [~, X] = ode45(@ode, t, 0);
    omega_sim = X(:,1);
end

~40 lines of code, runs in 1-2 s on a decent laptop, no add-on toolbox.

• fminsearch is a derivative-free solver (Nelder-Mead). Slow on big problems but always works as long as your initial guess \(p_0\) is in the right basin.
• Use an absolute value of the PWM duty: the rotor spins the same way for \(|u|\) and \(-|u|\) (single-quadrant drive on the Picooz). This is why the dataset ships u_abs.
• The ode45 call is what dominates the runtime. Decimating the dataset by 5-10× for identification is a standard speed-up; validate the fit on the full-resolution data afterwards.

This is the headline level. L2 proves the point: you don’t need a paid toolbox. You are doing exactly what the toolbox does — simulate the model, compute the error, minimise — just without the GUI.

Level 3 — lsqnonlin (Optimization Toolbox)

Same model, same cost — but a gradient-based solver with numerical finite-difference Jacobians. Typically 5-10× faster than fminsearch because it converges in a handful of iterations rather than hundreds.

lsqnonlin expects the residual vector (length \(N\)) rather than a scalar cost.

residual = @(p) simulate_picooz(p, t, u_abs) - omega_meas;

opts = optimoptions('lsqnonlin', ...
    'Algorithm','levenberg-marquardt', ...
    'Display','iter','TolFun',1e-10,'MaxIterations',200);

lb = [1e-3; 0;  0];   ub = [10; 1; 10];
p_fit = lsqnonlin(residual, p0, lb, ub, opts);

• Same script, same simulator as L2 — only the solver call changes.
• Accepts bounds (lb, ub), so you can constrain parameters to physical ranges (\(\tau > 0\), etc.).
• Requires Optimization Toolbox (widely available in academic MATLAB licences).

Level 4 — Parameter Estimator GUI (System Identification Toolbox)

The System Identification Toolbox (paid add-on) plus Simulink Design Optimization offer the visual Parameter Estimator app:

1. Drop the grey-box ODE into a Simulink model with tunable parameters.
2. Open Apps → Parameter Estimator.
3. Import (t, u, omega_meas).
4. Select the parameters to estimate, choose an algorithm (Levenberg-Marquardt, nonlinear least-squares, trust-region-reflective, pattern search).
5. Press Estimate — live-updating convergence plots.

What you gain: visual feedback, easy parameter bounds and initial-guess sweeps, one-click validation on held-out data. What you pay: a licence for the toolbox, the Simulink wrap around your ODE, and some vendor-lock — the script is hidden behind the GUI.

Every result is reproducible on L2 or L3 with a dozen lines of code. The GUI is a quality-of-life win, not a theoretical one.

Comparison on the real dataset

Running all three script methods (L1, L2, L3) on the chirp dataset in this tutorial:

Exact numbers depend on initial guess and solver tolerances — re-run the scripts on your own machine to see the values converge.

L1 fails — the “tau” it returns is an artefact of fitting an inadequate linear model to non-linear data. This is the pedagogic red flag: always plot the fit residual. If it is structured (not noise-like), your model structure is wrong.

L2, L3, L4 agree on the same non-linear parameters. Pick whichever tool you own.

Validation on held-out data — don’t overfit

Split the 30 s chirp into two halves, fit on the first 15 s, simulate on the second 15 s with the fitted parameters, and compare. If the second-half fit looks similar to the first-half fit, your model generalises. If it does not, you overfit: either reduce the number of free parameters, or enrich the model structure.

This is always the final sanity check. Three lines of MATLAB:

t_fit = t(t <= 15);    u_fit = u_abs(t <= 15);    w_fit = omega_meas(t <= 15);
t_val = t(t >  15);    u_val = u_abs(t >  15);    w_val = omega_meas(t >  15);
p_fit = lsqnonlin(@(p) simulate_picooz(p, t_fit, u_fit) - w_fit, p0, lb, ub);
w_pred = simulate_picooz(p_fit, t_val, u_val);
rmse_val = sqrt(mean((w_pred - w_val).^2));   % should be close to rmse_fit

Which method should you use?

• Just starting out? L2 — base MATLAB, 40 lines, works on anything.
• Automated CI / batch? L2 or L3 — scripts, no GUI dependency.
• Largest problem / fastest solver? L3 — lsqnonlin with bounds.
• Understanding the math? L0 → L1 → L2 — build up in layers.
• Have the toolbox? L4 — not more accurate but the visual feedback is nice.

Common pitfalls

• Poor initial guess. fminsearch can land in a local minimum. Do L0 first.
• Unit mismatch. BEMF voltage is not rad/s directly. The scaling is absorbed into k — don’t try to physically interpret k unless you also calibrated the BEMF → speed constant.
• Over-fitting. Adding parameters always improves the fit on the training set. Validate on held-out data.
• Numerical stiffness. The exponential drag term can blow up for large k2. Bound k2 ≤ 1 in lb/ub with L3.
• Sign conventions. Rotor speed is always positive. The Picooz tutorial uses u_abs = abs(u) to avoid confusion.

What’s next

You now have an identified model. In the next step we design a controller that keeps the rotor at a commanded speed setpoint, rejecting disturbances (wind, battery voltage drift, rotor wear).

Next → 4. Controller Design

Back → 2. Acquisition with External Mode