# Why does this transfer function estimation not work? System identification

Goal:

I have an unknow dynmical system $$G(s)$$ and I want to find it from measurement data, output $$y(t)$$ and input $$u(t)$$. The data is frequency responses.

Method:

I begun first with creating the data.

$$u(t) = A sin(2\pi \omega (t) t)$$

Where $$\omega(t)$$ is frequency in Hz over time and $$A$$ is fixed amplitude. Let's say that we know our model, just to make our data inside the computer.

t = linspace(0.0, 50, 2800);
w = linspace(0, 100, 2800);
u = 10*sin(2*pi*w.*t);
G = tf([3], [1 5 30]);
y = lsim(G, u, t);


Now when we have our data $$u(t)$$ and $$y(t)$$ and also $$\omega(t)$$. We can use Fast Fourier Transform to estimate the model.

First we find the complex ratio between $$u(t)$$ and $$y(t)$$ in frequency domain.

$$G(z) = \frac{FFT(y(t))}{FFT(u(t))}$$

  % Get the size of u or y or w
r = size(u, 1);
m = size(y, 1);
n = size(w, 2);
l = n/2;

% Do Fast Fourier Transform for every input signal
G = zeros(m, l*m); % Multivariable transfer function of magnitudes
for i = 1:m
% Do FFT
fy = fft(y(i, 1:n));
fu = fft(u(i, 1:n));

% Create the complex ratios between u and y and cut it to half
G(i, i:m:l*m) = (fy./fu)(1:l); % This makes so G(m,m) looks like an long idenity matrix
end

% Cut the frequency into half too and multiply it with 4
w_half = w(1:l)*4;


Wee need to divide it into half due to frequencies have mirrors.

Now when we got our complex ratios. We need to create a discrete transfer function on this form:

$$G(z^{-1}) = \frac{B(z^{-1})}{A(z^{-1})}$$

$$A(z^{-1}) = 1 + A_1 z^{-1} + A_2 z^{-2} + A_3 z^{-3} + \dots + A_p z^{-p}$$ $$B(z^{-1}) = B_0 + B_1 z^{-1} + B_2 z^{-2} + B_3 z^{-3} + \dots + B_p z^{-p}$$

Where $$p$$ is the model order.

Now we are going to solve this as least squares.

$$A(z^{-1})G(z^{-1}) = B(z^{-1})$$

$$G(z^{-1}) = -A_1G(z^{-1})z^{-1} - \dots -A_pG(z^{-1})z^{-p} + B_0 + B_1 z^{-1} + \dots + B_p z^{-p}$$

Like this: $$\begin{bmatrix} G(z_1^{-1})z_1^{-1} & \dots & G(z_1^{-1})z_1^{-p} & 1 & z_1^{-1} & \dots & z_1^{-p} \\ G(z_2^{-1})z_2^{-1} & \dots & G(z_2^{-1})z_2^{-p} & 1 & z_2^{-1} & \dots & z_2^{-p} \\ G(z_3^{-1})z_3^{-1} & \dots & G(z_3^{-1})z_3^{-p} & 1 & z_3^{-1} & \dots & z_3^{-p} \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ G(z_l^{-1})z_l^{-1} & \dots & G(z_l^{-1})z_l^{-p} & 1 & z_l^{-1} & \dots & z_l^{-p} \end{bmatrix}$$

$$\begin{bmatrix} -A_1\\ \vdots \\ -A_p\\ B_0\\ B_1\\ \vdots \\ B_p \end{bmatrix}$$

$$= \begin{bmatrix} G(z_1^{-1})\\ G(z_2^{-1})\\ G(z_3^{-1})\\ \vdots \\ G(z_l^{-1}) \end{bmatrix}$$

Where $$z_i = e^{j\omega_i T}$$ where $$T$$ is the sample ratio of measurement.

Let's call this equation above for $$Ax=B$$

MATLAB / Octave code for that:

  Gz = repmat(G', 1, p);
Ir = repmat(eye(r), l, 1); % Just a I column for size r and length l
Irz = repmat(eye(r), l, p);
for n = 1:l
for j = 1:p
z = (exp(1i*w_half(n)*sampleTime)).^(-j); % Do z = (e^(j*w*T))^(-p)
sn = (n-1)*m + 1; % Start index for row
tn = (n-1)*m + m; % Stop index for row
sj = (j-1)*m + 1; % Start index for columns
tj = (j-1)*m + m; % Stop index for columns
Gz(sn:tn, sj:tj) = Gz(sn:tn, sj:tj)*z;    % G'(z^(-1))*z^(-1)
Irz(sn:tn, sj:tj) = Irz(sn:tn, sj:tj)*z;  % Ir*z^(-1)
end
end
% Join them all
A = [Gz Ir Irz];


Now I going to solve this equation. We need to take accound that there are only complex values here. So we will solve this as:

$$\begin{bmatrix} real(A)\\ imag(A) \end{bmatrix}x = \begin{bmatrix} real(B)\\ imag(B) \end{bmatrix}$$

  Ar = real(A);
Ai = imag(A);
Gr = real(G');
Gi = imag(G');
A = [Ar; Ai];
B = [Gr; Gi];
x = (inv(A'*A)*A'*B)'; % Ordinary least squares


And the numerator and denominator from $$x$$ is

  den = [1 (x(1, 1:p))] % -A_1, -A_2, -A_3, ... , -A_p
num = (x(1, (p+1):end)) % B_0, B_1, B_2, ... , B_p


And here is the problem.

The variable $$den$$ have poles that are larger than 1 in unit circle. That' means that the model is unstable.

Question:

What have I missed? What need to be done?

I assume that the least squares was not made correct. Right?

What I have checked:

I have checked that this code is correct:

  % Get the size of u or y or w
r = size(u, 1);
m = size(y, 1);
n = size(w, 2);
l = n/2;

% Do Fast Fourier Transform for every input signal
G = zeros(m, l*m); % Multivariable transfer function of magnitudes
for i = 1:m
% Do FFT
fy = fft(y(i, 1:n));
fu = fft(u(i, 1:n));

% Create the complex ratios between u and y and cut it to half
G(i, i:m:l*m) = (fy./fu)(1:l); % This makes so G(m,m) looks like an long idenity matrix
end


Because I can plot the bode diagram of the measurement data

  % Cut the frequency into half too and multiply it with 4
w_half = w(1:l)*4;

% Plot the bode diagram of measurement data - This is not necessary for identification
if(w_half(1) <= 0)
w_half(1) = w_half(2); % Prevent zeros on the first index. In case if you used w = linspace(0,...
end
semilogx(w_half, 20*log10(abs(G))); % This have the same magnitude and frequencies as a bode plot


Assume that our model is

$$G(s) = \frac{3}{s^2 + 5s + 30}$$

There fore our bode diagram from data is going to look like this. The left picture shows the data-bode diagram and the right picture shows the bode diagram from the transfer function model.

You can follow the math logic at equation 14 here: https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19920023413.pdf

• You optimization constraint, does not take into account causality and stability as I see, so you are bound to see results which have poles outside unit circle. Isn't it? It's a plain unconstrained least sqaures problem – Dsp guy sam Apr 18 at 13:42
• @Dspguysam Yes. It's uncosntrained. As the report says. – Daniel Mårtensson Apr 18 at 13:42
• Added my thoughts in the answer – Dsp guy sam Apr 18 at 13:48

I see, it's a simple line curve fitting, you would need to cosntraint poles to be inside unit circle( this can be turned into a convex constraint), the objective of least sqaures is an $$l_2$$ norm minimization (which is also convex), so you would need to setup a convex optimization problem to ensure stability and poles inside unit circle.

One easier approach would be the following:

formulating the convex problem might be not so trivial, especially if not with optimization background, so I suggest that you

go ahead with this unconstrained problem, if you get a pole outside unit circle in the z plane, keep the pole at same frequency and scale magnitude of pole to lie just within unit circle, that should give you a very decent approximation of the frequency response.

Aside in general:

Since you mention that the system function is related to input and output as the following, pretty much describing an LTI system as $$G(z) = \frac{FFT(y(t))}{FFT(u(t))}$$

Then I would suggest the following, instead of taking a sinusoid as input, take white gaussian noise, suppose $$u(t)$$ is gaussian proceed that is IID for different time instances, then it's Fourier transform is simply $$\frac{N_o}{2}$$ for all frequencies. That means the Fourier transform if output $$y(t)$$ is simply $$\frac{N_o}{2}G(f)$$, so simply taking the FFT of the output of the system when white gaussian noise is passed through it, directly providers the system transfer function.

I think this is a much starightforward and easy approach. Can be easily simulates in MATLAB. Make sure to run Monte Carlo simulation over noise

• Can you show me? I must be a frequency response :) Thank you for your reply. – Daniel Mårtensson Apr 18 at 13:21
• You are looking for MATAB code?, Plot? – Dsp guy sam Apr 18 at 13:22
• There is no problem with FFT here. I have showed that. – Daniel Mårtensson Apr 18 at 13:22
• Yes. I'm looking for MATLAB code. Not plot. I know that my $G(z)$ are correct. It's more how I solve this system. Is this correct $z_i^{-p} = (e^{j\omega_i T})^{-p}$ ? – Daniel Mårtensson Apr 18 at 13:23
• By the way! Assume that we have a MIMO TF. Can we convert it to MIMO state space model then? – Daniel Mårtensson Apr 18 at 13:56