Размещено на http://www.allbest.ru/

Содержание

• Введение
• 1. Выбор и исследование моделей систем для идентификации
• 2. Метод идентификации с проектированием вход-выходных последовательностей в пространство состояний (subspace identification)
• 3. Исследование метода идентификации с проектированием вход-выходных последовательностей в пространство состояний при наличии шума измерений
• 4. Сравнение метода идентификации с проектированием вход-выходных последовательностей в пространство состояний с методом наименьших квадратов
• 4.1 Сравнение метода идентификации с проектированием вход-выходных последовательностей в пространство состояний с методом наименьших квадратов при наличии шума измерений
• Заключение
• Список литературы
• Приложение А
• Приложение Б

Введение

Проблема идентификации является в настоящее время одной из основных проблем теории и практики управления.

Идентификация - это определение математической модели объекта, или, точнее - определение оптимальной, в известном смысле, оценки истинного оператора реального объекта из заданного класса операторов с помощью анализа входных и выходных сигналов этого объекта. Без знания оператора объекта управления, ставящего в соответствие входные и выходные сигналы, нельзя эффективно решать задачи расчета и проектирования систем управления. Задачей идентификации является построение математической модели объекта, под которой понимается оператор, определяющий поведение объекта и описывающий все его информационные свойства.

Исходными данными для расчета и проектирования САУ является математическая модель объекта управления, которая обычно известна лишь частично (например, известна структура, но неизвестны параметры). Поэтому прибегают к экспериментальным исследованиям для определения динамики объектов управления. В некоторых случаях вообще отсутствует какая-либо априорная информация об объекте, поэтому как его структура, так и параметры определяются из эксперимента.

1. Выбор и исследование моделей систем для идентификации

Для исследования метода идентификации параметров системы зададим несколько произвольных передаточных функций высоко колебательных моделей систем (1), (2).

(1)

Данная система будет иметь следующие нули (Рисунок 1):

И полюса (Рисунок 1):

Рисунок 1 - Нули и полюса ПФ(1)

Для системы (1) построим логарифмические характеристики (Рисунок 2) и отклик системы на единичное ступенчатое воздействие, передаточную функцию (Рисунок 3).

Рисунок 2 - ЛЧХ и ЛФЧХ модели ПФ(1)

Рисунок 3 - Передаточная характеристика ПФ (1)

Далее рассмотрим вторую передаточную функцию для исследования (2).

(2)

Данная система будет иметь следующие нули (Рисунок 4):

И полюса (Рисунок 4):

Рисунок 4 - Нули и полюса ПФ(2)

Для системы (2) построим логарифмические характеристики (Рисунок 5) и отклик системы на единичное ступенчатое воздействие, передаточную функцию (Рисунок 6).

Рисунок 5 - ЛЧХ и ЛФЧХ модели ПФ(2)

Рисунок 6 - Передаточная характеристика ПФ (2)

Для получения выходных значений систем подадим на вход кусочно-постоянный сигнал с амплитудой равной 1 со случайным значением ширины (Рисунок 7).

Рисунок 7 - Входное воздействие

Рисунок 8 - Выход системы (1)

В результате подачи на вход системы (1) кусочно-постоянного сигнала (Рисунок 7) получим следующий выход (Рисунок 8).

В результате подачи на вход системы (2) кусочно-постоянного сигнала (Рисунок 7) получим следующий выход (Рисунок 9).

Рисунок 9 - Выход системы (2)

идентификация математическая последовательность управление

2. Метод идентификации с проектированием вход-выходных последовательностей в пространство состояний (subspace identification)

Метод идентификации линейных систем проектированием на подпространства (ИПП) был предложен бельгийскими учёными Бартом де Муром (Bart de Moor) и Питером Ван Оверши (Peter van Overschee) в 90-х годах прошлого века и в оригинале получил название "subspace identification" (букв. "подпространственная идентификация"). [1]

В основе данного метода лежат:

- геометрический подход, представленный ортогональным, либо косым проектированием строкового подпространства одной матрицы на строковое подпространство другой матрицы, а также определением углов между подпространствами.

- сингулярное разложение матриц для нахождения порядка идентифицируемой системы.

Ортогональное проектирование. Обозначим - оператор проектирования строкового пространства некоторой матрицы на строковое пространство матрицы:

:

, (3)

где обозначает "равенство по определению", - псевдообращение Мура-Пенроуза матрицы . Далее, обозначим ортогональную проекцию строкового пространства матрицы на строковое пространство матрицы .

Именно,

. (4)

На рисунке 10 представлена ортогональная проекция для , там же изображён оператор проектирования строкового пространства матрицы на ортогональное дополнение строкового пространства матрицы , т.е.

, (5)

, (6)

- единичная матрица размерности .

Рисунок 10 - Интерпретация ортогональной проекции на пространство

Косое проектирование. Определим косое проектирование строкового пространства матрицы на строковое пространство матрицы вдоль (параллельно) строкового пространства матрицы :

. (7)

Геометрическая интерпретация с "расшифровкой" выражения (7) представлена на рисунке 11.

Рисунок 11 - Интерпретация косой проекции на пространство

Если на рисунке 11 отмеченные точки считать границами соответствующих отрезков, то из подобия треугольников имеем:

. (8)

Выражая из последнего равенства , мы получим формулу (7) уже в матричном виде с учётом строкового базиса матрицы .

Идентификация линейных дискретных систем: детерминированный случай [1]

Постановка задачи: Для заданных измерений входного сигнала:

и выходного сигнала:

линейной системы:

(9)

определить:

· порядок системы ;

· системные матрицы:

с точностью до преобразования подобия.

На рисунках 7,8,9 изображены выборки входных и выходных сигналов систем, используемых для идентификации.

Решение. Введём блочные матрицы Ганкеля по входу (выходу), играющие важнейшую роль в методе ИПП. Данные матрицы формируется только из измеренных входных (и выходных) данных. Именно,

, (10)

, (11)

где:

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

- Число столбцов типично равно , что говорит о том, что все измерений используются для идентификации.

- Нижние индексы обозначают первый и последний элементы первого столбца матрицы Ганкеля. Индекс "p" (past, прошлое) определяет прошлые измерения, а "f"(future, будущее) для будущих измерений. Матрицы (прошлые входы) и (будущие входы) определяются разделением матрицы на две равные части по блочных строк. Матрицы определяются сдвигом вниз на одну блочную строку прошлых измерений и удалением верхней блочной строки у будущих измерений. [1]

Блочные ганкелевы матрицы выхода определяются аналогично блочным матрицам Ганкеля входа. Далее, определим блочные ганкелевы матрицы, состоящие из данных входа и выхода. Именно,

. (12)

Добавляя по одной блочной строке к соответствующим матрицам входа и выхода "прошлого":

, (13)

определим составную матрицу:

. (14)

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

Введём матрицу, состоящую из некоторой последовательности векторов состояния:

, (15)

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

. (16)

Матрицы, связанные со свойствами идентифицируемой системы. Метод идентификации проектированием на подпространства (ИПП) активно использует матрицы управляемости и наблюдаемости, но в расширенной форме. Расширенная ,порядок системы, матрица наблюдаемости , где индекс снизу обозначает число блочных строк, определяется следующим образом:

. (17)

Мы предполагаем, что пара наблюдаема, т.е. ранг матрицы равен . Вторая матрица связана с матрицей управляемости, но столбцы расширенной матрицы управляемости расположены в обратном (реверсном) порядке:

. (18)

Предполагаем, что ранг матрицы равен . И, наконец, сформируем нижнетреугольную блочную матрицу Тёплица:

. (19)

Матричные уравнения, связывающие входные и выходные данные. Представленные ниже матричные уравнения, получены простой подстановкой в дискретной системе (9) прошлых состояний для получения будущих состояний. Например, для двух шагов получаем:

(20)

Используя данную процедуру, имеем следующую систему матричных уравнений, связывающих полученные выше матрицы:

, (21)

, (22)

. (23)

Уравнению (22) можно дать следующую геометрическую интерпретацию косого проектирования (рисунок 11)

Рисунок 11 - Векторы строкового пространства блочной ганкелевой матрицы получаются как сумма линейной комбинации векторов строкового пространства и линейной комбинации векторов строкового пространства блочной ганкелевой матрицы

Покажем, что векторы строкового пространства можно получить непосредственно из вход-выходных данных. Имеем:

(24)

. (25)

Тогда формула (22) примет вид:

. (26)

Спроектируем левую и правую части (26) на ортогональное дополнение блочной матрицы "будущего " входа:

По определению:

. (27)

Косая проекция будущих выходов на строковое блочное пространство прошлых входа и выхода вдоль будущих входов для (27) с использованием формулы (7) может быть найдена как:

, (28)

где обозначено:

. (29)

Таким образом, получаем следующее матричное соотношение:

, (30)

которое показывает, что матрица равна произведению расширенной матрицы наблюдаемости и последовательности будущих векторов состояния. А формула (29) говорит о том, что данная матрица полностью определятся косым проектированием вход-выходных данных.

Расширенную матрицу наблюдаемости определим следующим образом. После получения из соотношения (29) матрицы проведём её сингулярное разложение:

. (31)

Порядок системы определяется числом сингулярных чисел отличных от нуля (представлены матрицей ). С учётом выражение (30) представим (31) в виде:

.

Из чего следует, что:

(32)

Найденные выражения и соотношения позволяют сформировать следующий алгоритм идентификации с использованием метода ИПП.

Алгоритм ИПП. 1. Вычислить косые проекции:

,

.

2. Вычислить сингулярное разложение матрицы :

.

3. Определить порядок системы просмотром сингулярных чисел матрицы и получить разложение (31).

4. Найти матрицы и из соотношений:

,

где обозначает матрицу без последних строк.

5. Определить как:

.

6. Решить систему линейных матричных уравнений относительно искомых матриц системы:

.

Для модели (1) матрицы A, B, C, D будут иметь вид:

,

,

Для модели (2) матрицы A, B, C, D будут иметь вид:

,

,

,

Отметим, что здесь:

.

Алгоритм реализован с помощью пакета прикладных программ Matlab 7.9.0 (R2012b), код алгоритма в приложении А.

В результате идентификации модели системы (1) получили дискретную передаточную функцию: сначала в обратной разности:

.

В прямой разности:

Осуществим переход дискретной передаточной функции к непрерывной.

Непрерывная ПФ модели:

(33)

Данная идентифицированная система будет иметь следующие нули (Рисунок 12):

и полюса (Рисунок 12):



Рисунок 12 - Нули и полюса идентифицированной ПФ (33)

Логарифмические частотные характеристики полученной непрерывной модели представлены на рисунке 13.

Рисунок 13 - ЛАЧХ и ЛФЧХ идентифицированной модели (33)

Для идентифицированной модели получена реакция на ступенчатое входное воздействие, сравним ее с реакцией ступенчатое входное воздействие для исходной модели (1).

Рисунок 14 - ПХ идентифицированной (33) и исходной системы (1)

Далее таким же образом проведем идентификацию для модели системы (2).

В результате идентификации модели системы (2) получили дискретную передаточную функцию: сначала в обратной разности:

.

В прямой разности:

Осуществим переход дискретной передаточной функции к непрерывной.

Непрерывная ПФ модели:

(34)

Данная идентифицированная система будет иметь следующие нули:

(Рисунок 15):

И полюса (Рисунок 15):

.

Рисунок 15 - Нули и полюса идентифицированной ПФ (34)

Логарифмические частотные характеристики полученной непрерывной модели представлены на рисунке 16.

Рисунок 16 - ЛАЧХ и ЛФЧХ идентифицированной модели (34)

Для идентифицированной модели получена реакция на ступенчатое входное воздействие, сравним ее с реакцией ступенчатое входное воздействие для исходной модели (2).

Рисунок 17 - ПХ идентифицированной (34) и исходной системы (2)

3. Исследование метода идентификации с проектированием вход-выходных последовательностей в пространство состояний при наличии шума измерений

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

(35)

Рисунок 18 - Шум измерений

Идентифицируя систему (1) при уровне шума n=3 и n=5 получаем следующие передаточные характеристики (Рисунок 19,20)



Рисунок 19 - ПХ идентифицированной (33) и исходной системы (1) при наличии шума измерений (n=3)

Рисунок 20 - ПХ идентифицированной (33) и исходной системы (1) при наличии шума измерений (n=5)

При идентификации системы (2) при уровне шума n=7 и n=3 в (35) получаем следующие передаточные характеристики (Рисунок 21, 22).

Рисунок 21 - ПХ идентифицированной (34) и исходной системы (2) при наличии шума измерений (n=7)

Рисунок 22 - ПХ идентифицированной (34) и исходной системы (2) при наличии шума измерений (n=3)

В результате идентификации можно отметить, что при относительно среднем уровне шума система идентифицируется довольно хорошо.

4. Сравнение метода идентификации с проектированием вход-выходных последовательностей в пространство состояний с методом наименьших квадратов

Метод наименьших квадратов (МНК) один из методов идентификации параметров систем по сигналам входа и выхода. Суть метода заключается на вычисление матрицы регрессоров (36) с использованием входного и выходного сигналов, на основе, которой вычисляется вектор неизвестных параметров системы (37) по формуле (38).

(36)

(37)

(38)

Проведем идентификацию методом МНК для моделей систем (1) и (2).

Алгоритм реализован с помощью пакета прикладных программ Matlab 7.9.0 (R2012b), код алгоритма в приложении Б.

Для системы (1) получаем следующую непрерывную ПФ модели (39).

(39)

Для системы (2) получаем следующую непрерывную ПФ модели (40).

(40)

Передаточные характеристики идентифицированных систем представлены на рисунках 23, 24.

Рисунок 23 - ПХ идентифицированной методом МНК (39) рассматриваемым методом и исходной системы (1)

Рисунок 24 - ПХ идентифицированной методом МНК (40), рассматриваемым методом и исходной системы (2)

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

4.1 Сравнение метода идентификации с проектированием вход-выходных последовательностей в пространство состояний с методом наименьших квадратов при наличии шума измерений

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

Добавим шум измерений к одной из заданных моделей (1) для исследования.

Рисунок 25 - ПХ идентифицированной методом МНК, рассматриваемым методом и исходной системы (1)

Заключение

В научной работе был подробно рассмотрен метод идентификации с проектированием вход - выходных последовательностей в пространство состояний.

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

В результате исследования было установлено, что метод идентификации с проектированием вход - выходных последовательностей в пространство состояний наиболее точно идентифицирует параметры моделей систем по сравнению с методом наименьших квадратов. Также можно сказать, что метод идентификации с проектированием вход - выходных последовательностей в пространство состояний хорошо восприимчив к шумам измерений.

Список литературы

1. De Moor B., Van Overschee P. Subspace identification for linear systems: Theory. Implementation. Applications. - Boston/London/Dordrecht: Kluwer Academic Publishers, 1996. - 248 p.

2. Кожинский О. С, Грабарь Л.П. Об одном методе идентификации линейных объектов управления // Идентификация. Докл. II Всесоюзн. Совещ. По статистическим проблемам теории управления. - М.: Наука, 1970. - С. 39-50.

3. Льюнг Л. Идентификация систем. Теория для пользователей. - М.: Наука, 1991.

4. Методы классической и современной теории автоматического управления: Учебник в 3-х т. Т.З: Методы современной теории автоматического управления / Под ред. Н.Д. Егупова. - М.: Изд-во МГТУ им. Н.Э. Баумана, 2000. - 748 с, ил.

5. Методы классической и современной теории автоматического управления: Учебник в 3-х т. Т. 1: Анализ и статистическая динамика систем автоматического управления / Под ред. Н.Д. Егупова. - М: Изд-во МГТУ им. Н.Э. Баумана, 2000, - 748 с, ил.

Приложение А

Код программы идентификации с проектированием вход-выходных последовательностей в пространство состояний

clear all

close all

clc

h=0.001;

t=0:h:3;

u0(1)=1;

A0=1*u0(1);

A=A0*sign(randn);

u(1)=1;

for i=2:numel(t)

sw=rand;

if sw>=0.5 & sw<=0.6

A=A0*sign(randn);

end

u0(i)=0;

u(i)=u0(i)+A;

end

%

% figure

% plot(t,u)

% title('Вход системы');

% xlabel('u(t)')

% ylabel('t')

%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

%ПФ ОУ

% num=10*conv([1 3], [1-1]);

% den=conv([1 1 100.25], [1 10 925]);

%

num=55*conv([1 2 5], [1 1]);

den=conv([conv([1 4], [1 0.4 25.04])], [1 8 160]);

W_start=tf(num,den)

zerosW=zero(W_start)

polesW=pole(W_start)

%Модель ОУ в пр-ве состояний

Wo_ss=ss(W_start);

Ao=Wo_ss.a;

bo=Wo_ss.b;

co=Wo_ss.c;

% X(:,1)= [0 0 0 0]';

X(:,1)= [0 0 0 0 0]';

% ШУМ

Disp=1e-5;

v=randn(1,numel(u))*sqrt(Disp);

for i=2:numel(t)

X(:,i)=X(:,i-1)+h*(Ao*X(:,i-1)+bo*u(i));

y(i)=co*X(:,i-1)+v(1,i);

end

figure

plot(t,v)

grid on

% grid on

title('Шум измерений');

xlabel('v')

ylabel('t')

% %

figure

step(W_start)

grid on

title('ПХ исходной системы');

% % ЛЧХ

% figure

% bode(W_start)

% title('ЛФЧХ ');

grid on

% % нули, полюса

% figure

% pzmap(W_start)

% title('Нули и полюса системы');

%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

m=1; %dim u

l=1; %dim y

la=length(y);

mi=1;

Nf=fix(la);

%!!!!!!!!!!!!!!

i_past=370;

n_model=5;

% n_model=4;

%!!!!!!!!!!!!

N=fix(Nf/1);

M=fix(N/2);

U_0_2i_end=zeros(2*i_past*m,M);

Y_0_2i_end=zeros(2*i_past*l,M);

k=0;

for ii=1:m:2*i_past*m;

k=k+1;

U_0_2i_end(ii:ii+m-1,:)=u(:,k:M+k-1);

end

k=0;

for ii=1:l:2*i_past*l;

k=k+1;

Y_0_2i_end(ii:ii+l-1,:)=y(:,k:M+k-1);

end

Up=U_0_2i_end(1:m*i_past,:);

Uf=U_0_2i_end(m*i_past+1:end,:);

[U_f,S_f,V_f]=svd(Uf);

Yp=Y_0_2i_end(1:l*i_past,:);

Yf=Y_0_2i_end(l*i_past+1:end,:);

Wp= [Up;Yp];

Uf_minus=Uf(m+1:end,:);

Yf_minus=Yf(l+1:end,:);

Up_plus= [Up;Uf(1:m,:)];

Yp_plus= [Yp;Yf(1:l,:)];

Wp_plus= [Up_plus;Yp_plus];

Oi=oblique_projection(Yf,Uf,Wp);

Oi_1=oblique_projection(Yf_minus,Uf_minus,Wp_plus);

[U,S,V]=svd(Oi);

% S(1,1)

% S(2,2)

% S(3,3)

% S(4,4)

% S(5,5)

% S(6,6)

% S(7,7)

U1=U(:,1:n_model);

S1=S(1:n_model,1:n_model);

V1=V(1:n_model,:)';

Gi=U1*sqrtm(S1);

Gi_1=Gi(1:end-l,:);

Xi=pinv(Gi)*Oi; %Xi

Xi_plus=pinv(Gi_1)*Oi_1; %Xi_plus

Ui_i=Uf(1:m,:);

Yi_i=Yf(1:l,:);

% old var

% AA= [Xi;Ui_i];

% BB= [Xi_plus;Yi_i];

% XX=BB*pinv(AA);

% new variant

AA= [Xi;Ui_i]';

BB= [Xi_plus;Yi_i]';

XX=(pinv(AA)*BB)';

F_iden=XX(1:n_model,1:n_model);

G_iden=XX(1:n_model,n_model+1:end);

C_iden=XX(n_model+1:end,1:n_model);

D_iden=XX(n_model+1:end,n_model+1:end);

[num,den]=ss2tf(F_iden,G_iden,C_iden,D_iden,1);

W_sys_iden=tf(num,den,h)

pw_ssi=pole(W_sys_iden);

zw_ssi=zero(W_sys_iden);

W=(d2c(W_sys_iden)) %перевели в непререрывную

pw=pole(W)

zw=zero(W)

%

% % ЛЧХ

% figure

% bode(W)

% title('ЛФЧХ идентифицированной системы');

%

% grid on

%

% % нули, полюса

% figure

% pzmap(W)

% title('Нули и полюса идентифицированной системы');

% переходная характеристика

%

hold on

% figure

step(W)

legend('исходная ПХ системы','идентифицированная')

title('ПХ идентифицированной системы');

grid on

Приложение Б

Код программы идентификации методом наименьших квадратов

clear all

close all

h=0.01;

t=0:h:4;

u0(1)=1;

A0=1*u0(1);

A=A0*sign(randn);

u(1)=1;

for i=2:numel(t)

sw=rand;

if sw>=0.5 & sw<=0.6

A=A0*sign(randn);

end

u0(i)=0;

u(i)=u0(i)+A;

end

figure

plot(t,u)

title('Вход системы');

%ПФ ОУ

% num=10*conv([1 3], [1-1]);

% den=conv([1 1 100.25], [1 10 925]);

num=55*conv([1 2 5], [1 1]);

den=conv([conv([1 4], [1 0.4 25.04])], [1 8 160]);

W_start=tf(num,den)

zerosW=zero(W_start)

polesW=pole(W_start)

%Модель ОУ в пр-ве состояний

Wo_ss=ss(W_start);

Ao=Wo_ss.a;

bo=Wo_ss.b;

co=Wo_ss.c;

% X(:,1)= [0 0 0 0]';

X(:,1)= [0 0 0 0 0]';

% ШУМ

% Disp=1e-4;

% v=randn(1,numel(u))*sqrt(Disp);

for i=2:numel(t)

X(:,i)=X(:,i-1)+h*(Ao*X(:,i-1)+bo*u(i));

y(i)=co*X(:,i-1); %+v(1,i);

end

%

figure

plot(t,y)

grid on

title('Выход исходной системы');

figure

step(W_start)

grid on

title('ПХ исходной системы');

for i=6:numel(y);

F(i-5,:)= [y(i-1),y(i-2),y(i-3),y(i-4),y(i-5),u(i-1),u(i-2),u(i-3),u(i-4),u(i-5)];

y1(i-5,:)= [y(i)];

end

Q=((F'*F)^(-1))*F'*y1;

%

for k=1:5

a(k)=Q(k);

b(k)=Q(k+5);

end

% дискретная ПФ

p= [1 -a]; %знаменатель

n= [b]; %числитель

W=(d2c(tf(n,p,h)))

nuln=roots(n)

poln=roots(p)

% % ЛЧХ

% figure

% bode(W)

% title('ЛФЧХ идентифицированной МНК');

%

% grid on

%

% % нули, полюса

% figure

% pzmap(W)

% title('Нули и полюса идентифицированной МНК');

% переходная характеристика

hold on

step(W)

title('переходная характеристика идентифицированной МНК');

grid on

% % АФЧХ (nyquist)

% figure

% nyquist(W)

% title('АФЧХ идентифицированной МНК');

% grid on

Размещено на Allbest.ru