1
$\begingroup$

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.

enter image description here

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

$\endgroup$
  • $\begingroup$ 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 $\endgroup$ – Dsp guy sam Apr 18 at 13:42
  • $\begingroup$ @Dspguysam Yes. It's uncosntrained. As the report says. $\endgroup$ – Daniel Mårtensson Apr 18 at 13:42
  • $\begingroup$ Added my thoughts in the answer $\endgroup$ – Dsp guy sam Apr 18 at 13:48
0
$\begingroup$

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

| improve this answer | |
$\endgroup$
  • $\begingroup$ Can you show me? I must be a frequency response :) Thank you for your reply. $\endgroup$ – Daniel Mårtensson Apr 18 at 13:21
  • $\begingroup$ You are looking for MATAB code?, Plot? $\endgroup$ – Dsp guy sam Apr 18 at 13:22
  • $\begingroup$ There is no problem with FFT here. I have showed that. $\endgroup$ – Daniel Mårtensson Apr 18 at 13:22
  • $\begingroup$ 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}$ ? $\endgroup$ – Daniel Mårtensson Apr 18 at 13:23
  • 1
    $\begingroup$ By the way! Assume that we have a MIMO TF. Can we convert it to MIMO state space model then? $\endgroup$ – Daniel Mårtensson Apr 18 at 13:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.