# System identification of a resonant system

I want to write a matlab script that would identify a system from it's inputs and outputs. I have so far had good results with simple systems, but with this slightly more complex one I'm not able to arrive at the original equation. To do identification, I use a chirp signal and collect output data that covers a range of input frequencies. I then use arx() to estimate a transfer function from the input data.

My octave script is as follows:

pkg load signal
pkg load control

s = tf('s');

L = 0.1;
R = 80;
C = 1e-4;

% Alt2: LC in series, resistor in parallel with output (LRC circuit)
%G = (s * L + 1 / (s * C)) / (R + s * L + 1 / (s * C));

% Alt1: simple RC circuit
G = 1 / (1 + R * C * s);

% simplify G
G = minreal(G)

% number of samples
N = 10000;

% frequency range in hz
Fstart = 15; Fend = 200;

t = linspace(0, 1, N);
sample_time = (t(end) - t(1)) / N

% generate stimulus signal and generate the output data from "real" system
stimulus = chirp(t, Fstart, t(end), Fend)';
[response, T, X] = lsim(G, stimulus, t);

%plot(t, stimulus, t, response, 'r');

% now try to go back to the original transfer function based on the collected data
[Gest, x0] = arx(iddata(response, stimulus, sample_time), 2);
Gest = minreal(d2c(Gest)); [b, a] = ss2tf(Gest); Gest = tf(b, a)

% compare the transfer functions
bode(G, Gest)

input("..");


If I run the script using simple transfer function (Alt1), I'm able to identify the original system precisely using arx (original transfer function G and estimated transfer function Gest are identical): Short of an error from arx about matrix A being unstable the approximation is 100% accurate (error is due to the fact that I'm using second order approximation on a system Alt1 which is really a first order system so this is irrelevant. If I pass 1 as last argument to arx this error disappears).

However, if I try to do the same experiment with Alt2, which has a resonance at around 25Hz, the method fails: The estimated transfer function is way off the target.

What is going on? Why am I not able to identify the second system using this method? Is it because I need more data around the resonance? How can I get an accurate estimate of the second system?

• I tested this in matlab and I think your problem is that you did not increase the model order for Alt2. – Kwin van der Veen May 3 '18 at 17:23
• @KwinvanderVeen please test your suggestions before you post them. The original equation is second order so second order identification should suffice. I think the problem was that arx algorithm is simply not suitable for identifying a discontinuous model. I got the perfect result when I used moen4 instead. – Martin May 4 '18 at 6:57
• ARX tries to fit $A(q)\,y(t)=B(q)\,u(t)+e(t)$. I am not sure what octave does when you only specify n, but when you use arx(dat, 'na', na, 'nb', nb) instead, with na the order of the A(q) polynomial and nb the order of the B(q) polynomial plus one with the constraint that nb<=na. For Alt2 the order of B(q) is two, so nb=3, which would also imply na=3 (but at least in matlab I was also able to use na=2 as well). So I would expect that arx(iddata(response, stimulus, sample_time), 'na', 2, 'nb', 3) would give the correct results. – Kwin van der Veen May 4 '18 at 8:02
• Ah I see 'na', 2, 'nb', 3 worked. na = 3 did not work. The result is best so far with arx, but not a good fit as moen4 estimate. However, I'm confused. It says in docs for octave that following must be true: nb <= na yet I had to use nb = 3 and na = 2? It does make sense to me that the transfer function must be proper, however that would imply nb >= na instead. – Martin May 4 '18 at 9:06
• I also do not know why. Maybe it meant to say the order of B(q)<=order of A(q), because otherwise the model would not be causal. But since you have all the data it is still possible to fit a noncausal filter. – Kwin van der Veen May 4 '18 at 9:11

## 1 Answer

I think the problem was that arx identification algorithm is just not suitable for solving discontinuous models. I changed algorithm to moen4 with same order and got the perfect solution.

[Gest, x0] = moen4(iddata(response, stimulus, sample_time), 2); 