Цель:

У меня есть неизвестная динамическая система $G(s)$, и я хочу найти ее по данным измерений, вывести $y( t)$ и введите $u(t)$. Данные представляют собой частотные характеристики.

Метод:

Сначала я начал с создания данных.

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

Где $\omega(t)$ — частота в Гц во времени, а $A$ — фиксированная амплитуда. Скажем, что мы знаем нашу модель, просто чтобы сделать наши данные внутри компьютера.

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);

Теперь, когда у нас есть данные $u(t)$ и $y(t)$, а также < span class="math-container">$\omega(t)$. Мы можем использовать быстрое преобразование Фурье для оценки модели.

Сначала мы находим сложное соотношение между $u(t)$ и $y(t)$ по частоте домен.

$$G(z) = \frac{БПФ(y(t))}{БПФ(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;

Нам нужно разделить его пополам из-за того, что частоты имеют зеркала.

Теперь, когда мы получили наши сложные отношения. Нам нужно создать дискретную передаточную функцию на этой форме:

$$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}$$

Где $p$ — порядок модели.

Теперь решим это методом наименьших квадратов.

$$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}$$

Нравится: $$ \begin{bmatrix} G(z_1^{-1})z_1^{-1} & \dots & G(z_1^{-1})z_1^{-p} & 1 & z_1^{-1} & \dots & z_1^{ -п} \\ G(z_2^{-1})z_2^{-1} & \dots & G(z_2^{-1})z_2^{-p} & 1 & z_2^{-1} & \dots & z_2^{ -п} \\ G(z_3^{-1})z_3^{-1} & \dots & G(z_3^{-1})z_3^{-p} & 1 & z_3^{-1} & \dots & z_3^{ -п} \\ \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^{ -п} \end{bmatrix}$$

$$ \begin{bmatrix} -А_1\\ \vdots\\ -А_п\\ 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}$$

Где $z_i = e^{j\omega_i T}$, где $T$ — отношение выборки измерение.

Назовем это уравнение выше для $Ax=B$

Код MATLAB/Octave для этого:

  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];

Теперь я собираюсь решить это уравнение. Мы должны принять во внимание, что здесь есть только комплексные значения. Итак, мы решим это как:

$$\begin{bmatrix} реальный(А)\\ изображение(А) \end{bmatrix}x = \begin{bmatrix} реальный(В)\\ изображение (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

И числитель и знаменатель от $x$ равны

  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

И вот в чем проблема.

У переменной $den$ полюса больше 1 в единичной окружности. Это означает, что модель нестабильна.

Вопрос:

Что я пропустил? Что нужно сделать?

Я предполагаю, что метод наименьших квадратов был сделан неправильно. Правильно?

Что я проверил:

Я проверил правильность этого кода:

  % 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;

  % 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

Предположим, что наша модель

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

Поэтому наша диаграмма тела из данных будет выглядеть так. На левом рисунке показана диаграмма тела данных, а на правом рисунке показана диаграмма тела из модели передаточной функции.

enter image description here

Вы можете следовать математической логике в уравнении 14 здесь: https://ntrs.nasa.gov/archive/ nasa/casi.ntrs.nasa.gov/19920023413.pdf

2
MrYui 18 Апр 2020 в 16:21
Как я вижу, ваше ограничение оптимизации не учитывает причинно-следственную связь и стабильность, поэтому вы обязательно увидите результаты, которые имеют полюса за пределами единичного круга. Не так ли? Это простая неограниченная задача наименьших квадратов.
 – 
Dsp guy sam
18 Апр 2020 в 16:42
Да. Это неограничено. Как говорится в отчете.
 – 
MrYui
18 Апр 2020 в 16:42
Добавил свои мысли в ответ
 – 
Dsp guy sam
18 Апр 2020 в 16:48

1 ответ

Ясно, это простая подгонка кривой линии, вам нужно будет ограничить полюса внутри единичного круга (это можно превратить в выпуклое ограничение), целью наименьших квадратов является $ l_2$ минимизация нормы (которая также является выпуклой), поэтому вам нужно будет настроить задачу выпуклой оптимизации, чтобы обеспечить стабильность и полюса внутри единичного круга.

Одним из более простых подходов будет следующий:

Формулировка выпуклой задачи может быть не такой тривиальной, особенно если не с оптимизацией, поэтому я предлагаю вам

продолжайте решать эту неограниченную проблему, если вы получите полюс за пределами единичного круга в плоскости z, держите полюс на той же частоте и масштабируйте величину полюса так, чтобы он находился только в пределах единичного круга, это должно дать вам очень приличное приближение частотной характеристики .

Вообще в стороне:

Поскольку вы упомянули, что системная функция связана с вводом и выводом следующим образом, в значительной степени описывая систему LTI как $$G(z) = \frac{БПФ(y(t))}{БПФ(u(t))}$$

Затем я бы предложил следующее: вместо того, чтобы брать синусоиду в качестве входных данных, возьмите белый гауссовский шум, предположим, что $u(t)$ является гауссовым процессом, который является IID для разного времени экземпляров, то его преобразование Фурье просто $\frac{N_o}{2}$ для всех частот. Это означает преобразование Фурье, если выходные данные $y(t)$ просто $\frac{N_o}{2}G( f)$, так что просто беря БПФ выхода системы, когда через нее проходит белый гауссовский шум, мы напрямую получаем передаточную функцию системы.

Я думаю, что это очень прямолинейный и простой подход. Можно легко смоделировать в MATLAB. Обязательно запустите симуляцию Монте-Карло с шумом

0
Dsp guy sam 18 Апр 2020 в 17:53
Можешь ли ты показать мне? У меня должна быть частотная характеристика :) Спасибо за ответ.
 – 
MrYui
18 Апр 2020 в 16:21
Вы ищете код MATAB?, Сюжет?
 – 
Dsp guy sam
18 Апр 2020 в 16:22
Здесь нет проблем с БПФ. Я показал это.
 – 
MrYui
18 Апр 2020 в 16:22
Да. Я ищу код MATLAB. Не сюжет. Я знаю, что мои $G(z)$ верны. Это больше то, как я решаю эту систему. Правильно ли это $z_i^{-p} = (e^{j\omega_i T})^{-p}$ ?
 – 
MrYui
18 Апр 2020 в 16:23
1
Кстати! Предположим, что у нас есть MIMO TF. Можем ли мы преобразовать его в модель пространства состояний MIMO?
 – 
MrYui
18 Апр 2020 в 16:56