Лекции по «Моделирование и идентификация объектов систем автоматики»

20. Из соответствующих уравнений определяются значения α₁ и b.
21. Для кривых типа II и III определяются значения δ.
22. В зависимости от типа кривой по соответствующим формулам определяют коэффициент а.
23. Определяются коэффициенты дифференциального уравнения а₁, а₂, а₃ из соотношений (5). Для кривых типа II и III определяют значение b₁.
24. Определяются постоянные интегрирования в зависимости от типа.
25. Производится проверка с помощью графического построения, проверяется решение дифференциального уравнения.

 

 

6. ИДЕНТИФИКАЦИЯ  ЛИНЕЙНЫХ ДИНАМИЧЕСКИХ СИСТЕМ

 

 1. Условия идентифицируемости   линейных динамических систем

 

 

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

                                                            .

 

Состояние системы задается вектором x(t)   n- мерного эвклидового пространства и известны путем измерений следующие векторы

                                                                                       (1)

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

                                                                                                    (2)

 

Состояние системы, заданное известными векторами (1), назовем идентифицируемым, если для него существует матрица  А, удовлетворяющая равенствам (2).

Систему уравнений (2) для каждой строки матрицы А     

сложно записать систему уравнений:

                        для j =1,2,..., n                                     (3)

Условие существования решения задачи идентификации, т.е. систем (3), можно записать в виде

                                                                                        (4)

Из соотношений (3) параметры модели можно определить по формулам

 

                               для j =1,2,..., n                                  (5)

Подставим соотношения (2) в условие (4), получим неравенство

        det (x, Ax, ... , A^n-1x) ¹0                                                                             (6)

Укажем на связь между  задачами идентификации и управляемости. Эта связь  формулируется так:

для существования  решения  задачи идентификации в виде математической модели

при наблюдении вектора  состояния x(t)  достаточно, чтобы матрица А и вектор x(t) удовлетворяли условия вполне управляемости

                                      det (b, Ab, ... , A^n-1b) ¹0

где  b = x(t).

Поскольку матрица А наперед  неизвестна, то на практике условия  идентификации проверяют на основании условия (4).

Для дискретных линейных систем управления определение идентифицируемости формулируется  следующим образом:  если по известным  значениям векторов  x(k), x(k+1), ... , x(k+n) состояние дискретной линейной стационарной системы  n-го порядка

x(k+1) = Ax(k)

можно восстановить матрицу А, то систему  называют идентифицируемой.

Теорема. Решение задачи идентификации  дискретной линейной системы имеет  место, если

det (x(k), x(k+1), ... , x(k+n-1)) ¹0

 

 

2. Определение весовой функции из уравнения свертки

 

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

 

 

 

 

 

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

                                                                    (1)

предполагается, что      при , .

Введем теперь аппроксимацию входной  функции времени w(t) кусочно-постоянной функцией в N точках с шагом , причем 

  при                                                 (2)

h(t) примем постоянной между точками разбиения.

 

   при  nD<t<(n+1)D                                       (3)

В терминах ступенчатой аппроксимации  и интеграл (1) при

t = nD приближенно запишется в виде

                                            (4)

Обозначим вектор наблюдения выхода (размерности N) через

Y^T(T) = [y(D) y(2D) ... y(ND)]

и значений весовой функции

.

 

Перепишем уравнение (4) в векторно- матричном  виде

Y(T) = DWh(T)                                                                            (5)

Матрица W  определяется равенством

 

Отметим, что W –   левая треугольная матрица, w(0) на диагонали.

Теперь задача сведена к определению  из уравнения (5) вектора h значений весовой функции в точках фиксации. Т.к.   w –невырожденная  и .

Поэтому формально решение уравнения (5) можно записать в виде

                                                                                    (6)

Благодаря левой треугольной форме  W выражение для h можно переписать в рекуррентном виде

                                                              (7)

        где 

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

Если входной сигнал является функцией единичного скачка, алгоритм (7) заметно упрощается. В  этом случае w(iD) = 1 для всех i, (7) примет вид

.

Определив величину                     ,

получим                       при этом  H_n = H_n-1 +h_n-1.

 

 

3. Оценивание весовой функции по методу наименьших квадратов

 

Рассмотрим систему  с одним входом и одним выходом. Объект предполагается линейным и стационарным.

Выход системы запишем  в виде:

                                                                (1)

 h(t) - весовая функция (импульсная  переходная функция)

 x(t-t)- вход,  n(t) - невязка (иногда называют шумом),

 

TS - время установления, определяется как min интервал времени, измеренный от момента подачи импульсного сигнала до момента, когда реакция системы составит 5% пикового значения (рис. ). Входная и выходная переменные представлены в  формуле (1)  в виде отклонений от  своих математических ожиданий, т.е.

 

Представим уравнение (1)  в дискретном виде

 или                                                            (2)

Здесь - время установления,  - время измерения выхода,  N_i - содержит не только невязку в моменты времени n(iD), но и ошибки аппроксимации функции x(t-t).

В результате аппроксимации задача оценивания непрерывной функции h(t)  заменяется  (параметризуется) оцениванием конечного множества параметров  h₀ ,..., h_Ns-1 (называется дискретной импульсной переходной функцией)

Для упрощения представления запишем  уравнения (2) в матричном виде:

                                 (3)

Запишем матричное представление  в символическом виде

                                                                                        (4)

Задача сводится к определению вектора параметров b при заданной матрице А и вектора измерений z.

Критерием при оценивании вектора  параметров  b является выбор таких b,    которые минимизируют  сумму квадратов невязок на интервале измерений.

Положим        в матричном виде   J = n^T n    (5)

подставляя (4) в (5) получим  J=(z-Ab)^T (z-Ab)

Необходимо определить b^*, удовлетворяющие условию

                          .

Необходимым условием вычисления  J является выполнение условия экстремума                    .

 

Запишем уравнения (3) в виде сумм

                             (6)

 

Продифференцируем (6) J по компонентам вектора b:

=

= -2     при m = 0,1,..., N_S –1                    (7)

Формулу (7) представим в  матричном виде:

                                                     (8)

(8) является необходимым  условием экстремума J.

Достаточным условием при  расчете min J  является положительная определенность квадратной матрицы

Если формулу (8) продифференцируем  еще раз по b, то получим

                                                                  (9)

Если матрица  A^TA  - неособенная, и также правая часть не зависит от b в (9), то условие экстремума (8) является необходимым и достаточным условием минимума.

Перепишем (8) в виде:

A^TA b* = A^Tz  ,    отсюда   b* = (A^TA)^-1 A^Tz                                 (10)

Напомним, что

. 

В непрерывной форме уравнение (10) принимает вид

                     (11)

Уравнение (11) – уравнение Винера-Хопфа  и может быть переписано в виде

где R_xx(t) -  автокорреляционная функция и R_xz(t) - взаимная корреляционная функция.

 

Автокорреляционной функцией случайного процесса X(t) называется неслучайная функция двух аргументов K_xx(t,t‘), которая при каждой паре значений t,t‘ равна корреляционному моменту  соответствующему случайному сечению функции X(t).

K_xx(t,t‘) = .

 

7. Регрессионный метод идентификации линейных систем  (Метод наименьших квадратов)

Рассмотрим уравнение:

 

                 ,                         (2.1)

 

В уравнении заменим  производную конечной разностью:

 

                               ,                                        (2.2)

 

                               ,                                                 (2.3)

 

                                    ,                                                      (2.4)

Подставим выражения (2.2)-(2.4) в уравнение (2.1):

,               (2.5)

Приведя подобные, получим:

,        (2.6)

Представим (2.6) в виде:

 

                 ,                                         (2.7)

где                                               

      .

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

Рассматриваем систему, которая задана уравнением:

                                                       

Минимизируемая функция  имеет вид:

    

                                 (2.8)

Запишем систему уравнений  для нахождения , для этого найдем частные производные:                            (2.9)             

                   (2.10)

Запишем систему в матричной  форме:

          (2.11)

                .                                           (2.12)

 

По полученным найдем коэффициенты а_i

.

 

1. Градиентные методы идентификации  нелинейных систем

 

Для общей задачи минимизации функционала

                                                                             

 

при ограничениях

где f – нелинейная вектор - функция

Случай  p_i  - const  будем задаваться начальным значением pⁱ , i- номер итерации, и решив систему дифференциальных уравнений оценим величину функции штрафа Jⁱ. Слегка изменяя pⁱ ,  для нового значения найдем штраф

Jⁱ +h_j,    j = 1,n,   n - число неизвестных коэффициентов

j-ю компоненту вектора-градиента функции штрафа можно приближенно оценить как

                                                                      (3)

Повторяя эту процедуру для  возмущений различных компонент  вектора параметров, определим приближенное значение вектора-градиента dJ/dpⁱ  Первое приращение вектора параметров в направлении наискорейшего спуска к минимуму функции штрафа составит

                                                (4)

и  Kⁱ – выбирается из условия

                                         (5)

 

а новое приближение  для вектора параметров определится  как

pⁱ⁺¹ = pi +Dpⁱ.

Простота приближенного  метода позволила положить его в  основу нескольких итерационных схем, однако трудно оценить ошибку, связанную  с приближенным вычислением dJ/dpⁱ  (процедура точного вычисления свелась бы к уже известным алгоритмам решения динамических задач). Приближенная процедура приводит  к существенным ошибкам, особенно тогда, когда функция штрафа не очень чувствительна к изменению вектора параметров. Последнее, к сожалению, довольно часто имеет место, если измерения или наблюдения искажены помехой и имеются неизвестные входные сигналы.