3.3. Идентификация тоу с помощью пакета sit

1.Загрузка в рабочую область MATLAB исходных данных эксперимента занесенных в файл DATA с помощью команды:

>> load datta

2. Вводим интервал дискретизации (интервал через который производились измерения входной и выходной величин):

>> ts=7

3. Формируем файл исходных данных для пакета SIT:

>> dan13 = iddata (y(1401:1500), u(1401:1500), ts)

4. Вводим названия входной и выходной переменных (u, y), при помощи:

>> dan13.inputn = 'Rashod para'

>> dan13.outputn = 'Temperatura'

5. Вводим единицы измерения переменных (u, y):

>> dan13.inputuntil = 'm^3/hour'

>> dan13.outputuntil = 'grad'

6. Посмотрим полную информацию о файле:

>> get(dan13)

ans =

Domain: 'Time'

Name: []

OutputData: [100x1 double]

y: 'Same as OutputData'

OutputName: {'Temperatura'}

OutputUnit: {'grad'}

InputData: [100x1 double]

u: 'Same as InputData'

InputName: {'Rashod para'}

InputUnit: {'m^3/hour'}

Period: Inf

InterSample: 'zoh'

Ts: 7

Tstart: []

SamplingInstants: [100x0 double]

TimeUnit: ''

ExperimentName: 'Exp1'

Notes: []

UserData: []

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

>> plot (dan13)

8. При использовании исходных данных произведём предварительную обработку этих данных с целью удаления тренда из набора, и если необходимо, отфильтруем данные с помощью средств имеющихся в пакете “SIT”, разделив данные на две половины:

- Dan13v – используется “Matlab” для построения модели объекта

- Dan13e – используется для проверки адекватности полученной модели

Открываем графический интерфейс пользователя пакета SIT:

>> ident

9. Проводим предварительную обработку файла исходных данных

Operations → Quick start:

dand – исходные данные без тренда

dande – половина файла с которой работает пакет SIT (с этой сравниваем)

dandv – половина файла которая оставлена для сравнения с результатами регулирования (с этой работаем)

10. Приступаем к параметрической идентификации (получение различных видов мат. модели объектов). Нахождение оценок различных моделей производится нажатием на клавишу Estimate:

- Parametric models – параметрическая идентификация

- Spectral model – производится оценивание частотных характеристик модели

- Correlation model – оценивание модели импульсной характеристики

Для нахождения параметров уравнений, описывающих ТОУ, необходимо провести параметрическую индикацию:

А) ARX – модель авторегрессии с дополнительным входом:

Б) ARMAX – модель регрессии скользящего среднего:

В) ОЕ – модель выход-ошибка:

Г) BJ – модель Бокса-Дженкинса:

Д) State space – модель для переменных состояний :

11. После того, как получили необходимое количество моделей, определим наиболее близкую модель путем сравнения коэффициентов адекватности (Model Output):

Выбираем модель с хорошими показателями адекватности (в моем случае это модель “n4s3”) т.к. ее показатели адекватности составляют 83%.

Преобразование модели

Проведем преобразование модели. Все полученные модели в рабочей области Matlab представлены в θ-формате (внутренний вид матричных моделей) и являются дискретными. Для преобразования моделей из тетта-формата в вид, удобный для дальнейшего использования при анализе и синтезе систем автоматики в пакте SIT есть специальные функции:

12. Преобразуем модель тета-формата в вектор передаточных функций:

>> [num,den]=th2tf(n4s3)

num =

0 0.0081 0.0022 0.0600

den =

1.0000 -1.9609 1.4657 -0.4401

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

>> zn4s3=tf(num,den,7)

Transfer function:

0.008076z^2 + 0.002226 z + 0.05999

z^3 - 1.961 z^2 + 1.466 z - 0.4401

Sampling time: 7 (Осуществление выборки времени)

14. Преобразуем дискретную передаточную функцию в непрерывную:

>> Ws=d2c(zn4s3)

Transfer function:

0.00608 s^2 - 0.001792 s + 0.000318

s^3 + 0.1173 s^2 + 0.01401 s + 0.0002927

15. Преобразуем дискретную модель тета-формата в непрерывную модель:

>> sn4s=thd2thc(n4s3)

State-space model: dx/dt = A x(t) + B u(t) + K e(t)

y(t) = C x(t) + D u(t) + e(t)

A =

x1 x2 x3

x1 -0.018423 0.0062662 -0.019034

x2 0.033099 -0.011022 -0.074865

x3 -0.080802 0.17317 -0.087814

B =

Rashod para

x1 -0.0062025

x2 0.019561

x3 -0.037115

C =

x1 x2 x3

Temperatura -2.8015 -0.33278 0.12898

D =

Rashod para

Temperatura 0

K =

Temperatura

x1 -0.0052945

x2 -0.035445

x3 0.16992

x(0) =

x1 0

x2 0

x3 0

Estimated using N4SID from data set dan13de

Loss function 0.00347817 and FPE 0.00500517

16. Получим числитель и знаменатель непрерывной передаточной функции:

>> [num1,den1]=th2tf(n4s3)

num1 =

0 0.0081 0.0022 0.0600

den1 =

1.0000 -1.9609 1.4657 -0.4401

17. Получаем передаточную функцию непрерывной системы:

>> sysn4=tf(num1,den1)

Transfer function:

0.008076s^2 + 0.002226 s + 0.05999

s^3 - 1.961 s^2 + 1.466 s - 0.4401

18. Построим переходную характеристику ТОУ для дискретной и непрерывной моделей и определим основные показатели переходного процесса. Для этого можно воспользоваться командой step(n4s4,sn4s), либо командой plot(n4s4,sn4s). Различие заключается в том, что в последнем случае представляется возможность использовать все достоинства LTI view:

>> plot(n4s3,sn4s)

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

Для дискретной модели:

• время нарастания переходного процесса (Rise time) – 84с.

• время регулирования (Setting time) – 175с.

•установившееся значение выходной координаты (Final value) – 1,09

Для непрерывной модели:

• время нарастания переходного процесса (Rise time) – 83,6с.

• время регулирования (Setting time) – 168с.

• установившееся значение выходной координаты (Final value) – 1,09

19. Для построения импульсной характеристики непрерывной модели необходимо воспользоваться командой impulse(sn4s), либо, щелкнув правой кнопкой мыши в поле графика LTI view, выбрать опцию Plot Types → Impulse:

• пиковая амплитуда (Peak amplitude) составляет – 30,9.

• время регулирования составляет (Setting time) – 175 с.

20. Определим частотные характеристики моделей с помощью команды bode(sn4s) либо, щелкнув правой кнопкой мыши в поле графика LTI view, выбрать опцию Plot Types → Bode:

Значение запасов устойчивости можно определить с помощью команды

>> [Gm,Pm,Wcg,Wcp] = margin (n4s3,sn4s)

Gm = 2,5875 – запас устойчивости по амплитуде в натуральных величинах на частоте Wcg

Pm = 145,4173 – запас устойчивости по фазе на частоте Wcp

Wcg = 0,0851

Wcp = 0,0110

В логарифмическом масштабе

Gmlog=20*log10(Gm) = 8,2577

На графиках частотных характеристик указаны значения запасов устойчивости по амплитуде (Gain Margin), которая для непрерывной модели ΔL = 9,09 dB.

21. Для построения АФХ необходимо воспользоваться командой nyquist(sysn4s), либо, щелкнув правой кнопкой мыши в поле графика LTI view, выбрать опцию Plot Types → Nyquist.

22. Для наглядности построим график изменения Y и определим основные статистические характеристики помехи с помощью команды:

>> plot (y)

Для получения статистических характеристик необходимо в строке меню графика в позиции Tools выбрать опцию Data statistics. Результатом выполнения команды явится окно, в котором будут указаны основные статистические характеристики случайного процесса изменения во времениy(t):

Статистические характеристики

• min и max – минимальное и максимальное значения помехи. Для нашего случая – 303.2 и 306.3 соответственно;

• mean – арифметическое среднее значение (304.9);

• median – медиана процесса (304.9);

• std – среднеквадратическое отклонение (0,8296);

• range – диапазон изменения помехи от минимального до максимального значения (3.05).

Управляемость и наблюдаемость

Для решения задач анализа и синтеза системы управления важно знать, является ли объект управляемым и наблюдаемым. Управляемость объекта заключается в способности объекта после подачи управляющего воздействия перейти из заданной начальной точки в заданную конечную. Наблюдаемость объекта заключается в возможности выяснить состояние объекта (вектора фазовых координат) по измеренным значениям выходной переменной на некотором временном интервале.

Объект называется вполне управляемым, если выбором управляющего воздействия U(t) можно перевести его на интервале времени от t₀ до t_к из любого начального состояния y(t₀) в произвольное заданное конечное состояние y(t_к). Критерием управляемости линейных стационарных объектов является условие: для того, чтобы объект был вполне управляемым, необходимо и достаточно, чтобы ранг матрицы управляемости M_и был равен размерности вектора состояний n.

Критерием наблюдаемости линейных стационарных объектов является условие: для того, чтобы объект был вполне наблюдаемым, необходимо и достаточно, чтобы ранг матрицы наблюдаемости М_Y = (C^T A^TC^T (A^T)²C^T … (A^T)^n-1C) равнялся размерности вектора состояния n = rang M_Y.

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

23. Для нахождения матрицы управляемости M_u и матрицы наблюдаемости М_Y необходимо знать матрицы A, B, C, D управления переменных состояний. Воспользуемся матрицами модели в пространстве состояния:

>> [A,B,C,d]=ssdata(sn4s)

A =

-0.0184 0.0063 -0.0190

0.0331 -0.0110 -0.0749

-0.0808 0.1732 -0.0878

B =

-0.0062

0.0196

-0.037

C =

-2.8015 -0.3328 0.1290

d =

0

Вычислим матрицу управляемости:

>> Mu=ctrb(A,B)

Mu =

-0.0062 0.0009 -0.0001

0.0196 0.0024 -0.0005

-0.0371 0.0071 -0.0003

Определим ранг матрицы управляемости:

>> n=rank(Mu)

n =

3

В данном случае ранг матрицы управляемости равен 3 и размерность вектора состояния равна 3, т.е. объект вполне управляем.

Вычислим матрицу наблюдаемости:

>> My=obsv(A,C)

My =

-2.8015 -0.3328 0.1290

0.0302 0.0084 0.0669

-0.0057 0.0117 -0.0071

Вычислим ранг матрицы наблюдаемости:

>> n=rank(My)

n =

3

В данном случае ранг матрицы наблюдаемости равен 3 и размерность вектора состояния равна 3, т.е. объект вполне наблюдаем.