654
.pdfгде U – единичный столбец, введенный для идентификации парамет-
|
|
1 |
u (1) |
u |
|
(1) |
K u |
n |
(1) |
|
|
|
|
|
1 |
|
2 |
|
|
|
|
|
|
ра |
a , |
U = 1 |
u1(2) |
u2 (2) |
K un (2) |
; |
|||||
|
0 |
M |
M |
|
M |
M |
M |
|
|
||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
1 u1(k) u2 (k) |
K un (k) |
|
a0
Am – матрица коэффициентов модели, Am = a1 .
Man
Ошибка оценивания, или функция невязки, имеет вид
|
y(1) |
|
ym (1) |
|
|
|
|
|
|
|
|
|
|
E =Y −Y |
= y(2) |
− ym (2) . |
(1.48) |
|||
m |
M |
|
|
M |
|
|
|
|
|
|
|
|
|
|
y(k) |
ym (k) |
|
Тогда критерий оптимальности определяется по методу наименьших квадратов:
J = ET E. |
(1.49) |
Таким образом, наилучшая (в смысле наименьших квадратов) оценка матрицы Am определяется как решение следующего уравнения:
|
|
|
|
|
|
дJ |
= 0. |
|
|
|
|
|
(1.50) |
|
|
|
|
|
|
|
|
|
|
|
|
|
|||
|
|
|
|
|
|
дA |
|
|
|
|
|
|
||
|
|
|
|
|
|
|
m |
|
|
|
|
|
|
|
После подстановки (1.48) в (1.49) уравнение (1.49) имеет вид |
||||||||||||||
|
дJ |
|
д{ET E} |
|
|
д{(Y −Y )T (Y −Y )} |
|
|||||||
|
|
|
= |
|
= |
|
|
m |
|
m |
|
= |
||
|
дA |
дA |
|
|
дA |
|
||||||||
|
|
|
|
|
|
|
|
|
||||||
|
|
m |
|
m |
|
|
|
|
|
|
m |
|
|
(1.51) |
|
|
д{(Y −UA )T |
(Y −UA )} |
|
|
|
||||||||
|
|
|
|
|
|
|||||||||
= |
|
|
m |
|
|
|
m |
|
|
= 0. |
|
|
|
|
|
|
дAm |
|
|
|
|
|
|
|
|||||
|
|
|
|
|
|
|
|
|
|
|
|
|||
Уравнение (1.50) представляет дифференцирование скалярной ве- |
||||||||||||||
личины ET E по матрице A , |
используя матричную алгебру, преобра- |
|||||||||||||
|
|
|
|
m |
|
|
|
|
|
|
|
|
|
|
зуем выражение (1.50):
41
дJ |
|
д{(Y −UA )T (Y −UA )} |
|
дtr{(Y −UA )(Y −UA )T } |
|
||||
|
= |
m |
m |
|
= |
m |
m |
|
= |
дA |
дA |
|
дA |
|
|||||
|
|
|
|
|
|
|
|||
m |
|
m |
|
|
|
m |
|
|
|
= дtr(YYT +UAm (UAm )T −Y (UAm )T −UAmY )T =
дAm
= дtr(YYT +UAm AmTU T −YAmTU T −UAmY )T =
дAm
|
дβααT γ |
|
|
T |
|
T |
|
|
|
|||
|
|
|
|
|
= |
(β γ |
|
+ γβ)α |
|
|
||
|
∂α |
|
|
|
|
|
||||||
= |
дβαT γ |
= γβ |
|
|
|
= |
|
|||||
дα |
|
|
|
|
|
|||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
дβαγ |
|
|
|
T |
|
T |
|
|
|
|
|
|
дα |
=β γ |
|
|
|
|
|
|
||||
|
|
|
|
|
||||||||
= 0 +(U T (U T )T +U TU )A |
−U TY −(U )T (YT )T = |
|
||||||||||
|
|
|
|
|
|
|
m |
|
|
|
|
|
= (U TU +U TU )A −U TY −U TY = |
|
|||||||||||
|
|
|
|
|
|
|
m |
|
|
|
|
|
|
= 2U TUA −2U TY = 0, |
|
||||||||||
|
|
|
|
|
m |
|
|
|
|
|
|
|
где tr – след матрицы. |
|
|
|
|
|
|
|
|
|
|
|
|
Уравнение (1.50) после преобразований имеет вид |
|
|||||||||||
|
U TUA −U TY = 0. |
(1.52) |
||||||||||
|
|
|
m |
|
|
|
|
|
|
|
|
Таким образом, наилучшая в смысле наименьших квадратов оценка матрицы параметров Am удовлетворяет уравнению
A =[U TU ]−1U TY = 0. |
(1.53) |
m |
|
Данное выражение и является основой для идентификации на основе линейной регрессии и метода наименьших квадратов.
Для того чтобы построенная модель была адекватна объекту, количество измерений k должно удовлетворять следующему условию:
k ≥ n +1. |
(1.54) |
Одномерный линейный регрессионный анализ может быть применим и для нелинейных систем. Проиллюстрируем следующим примером.
Пусть система описана квадратичным уравнением
42
y = a0 + a1x + a2 x2.
Нетрудно заметить, что уравнение (1.54) применимо для уравнения (1.45) при условии, что u1 = x,u2 = x2. Тогда матрица U определяется как
|
x(1) |
x |
2 |
(1) |
|
|
1 |
|
|
||||
1 x(2) |
x2 |
(2) |
||||
U = |
|
|
|
|
|
. |
M |
M |
|
|
|
M |
|
|
x(k) |
x |
2 |
|
|
|
1 |
|
|
(k) |
|||
|
|
|
|
|
|
|
a0
Оценивание вектора параметров A = a1 осуществляется по
a2
формуле (1.52).
Рассмотрим линейную статическую систему, представленную на рис. 1.15, имеющую n входов U и m выходов Y.
Рис. 1.15. Схема многомерной системы
Процесс, протекающий в многомерных системах, имеет n входов и m выходов и по аналогии с одномерным процессом может быть описан следующей системой уравнений:
n |
|
y1 = a10 +∑a1iui |
|
i=1 |
|
M |
(1.55) |
n
ym = am0 +∑amiui . i=1
43
u1
Используя серию измерений величин U = u2 и Y
Mun
y1 |
|
|
|
|
|
= y2 |
|
в тече- |
M |
|
|
|
|
|
ym |
|
ние определенного отрезка времени t [t0 , tk ], можно сформировать матрицу измерений величин U и Y следующим образом:
1
U = u1(1)M
un (k)
y1(1)
Y = Mym(1)
1 u1(2)
un (k)
y1(2)
ym(2)
L |
|
1 |
L |
u |
(k) |
|
1 |
|
|
|
|
|
|
|
L |
un (k) . |
Lym(k)
M
Lym(k)
Уравнение (1.55) можно переписать в матричном виде:
Ym = AmU. |
(1.56) |
Вычислим оптимальную оценку матрицы A способом, аналогичным способу для одномерных систем:
дJ |
|
д(ET E) |
|
д(Y −Y )T (Y −Y ) |
|
|
|
= |
|
= |
M |
M |
= |
дA |
дA |
дA |
|
|||
|
|
|
|
|||
m |
|
m |
|
m |
|
|
= дtr[(Y − AmU )(Y − AmU )T ] =
дAm
= дtr(YYT + AmU (AU )T −Y (AmU )T − AmUYT ) =
дAm
= дtr(YYT + AmUU T AmT −YU T AmT − AmUYT ) =
дAm
44
дαβαдα T = α(β+βT )
= |
дβαT |
=β |
= |
|
дα |
||||
|
|
|
||
|
дβαγ |
T |
T |
|
|
дα |
=β γ |
|
=0 + Am (UU T +(UU T )T ) −YU T −(I )T (UYT )T =
=Am (UU T +UU T ) −YU T −YU T =
= 2AmUU T −2UX T = 0.
Таким образом, наилучшая в смысле наименьших квадратов оценка матрицы параметров Am удовлетворяет уравнению
A =YU T (UU T )−1. |
(1.57) |
m |
|
Требования к количеству измерений состоят в том, что для адек- |
|
ватности модели необходимо, чтобы: |
|
l ≥ (n +1) m. |
(1.58) |
Рассмотренные выше методы оценивания одномерной и многомерной линейной модели реализуются по явной схеме идентификации, достоинства и недостатки которой были рассмотрены выше.
На практике целесообразнее применять итерационные методы, т.е. схемы с настраиваемой модели.
Линейный регрессионный анализ может быть реализован в замкнутой схеме идентификации (схеме с настраиваемой моделью).
Итерационные методы идентификации используют методы последовательного уточнения матрицы коэффициентов.
Формула итерационного регрессионного анализа имеет вид
A (k +1) = A (k) +[Y (k +1) − A (k)U (k)][P(k +1)U (k)]T , (1.59) |
||
m |
m |
m |
где Am (k +1), Am (k) – оцениваемая матрица коэффициентов в предыдущий и в текущий моменты наблюдений; Y (k +1) – вектор выхода, измеренный в текущий(k +1) момент наблюдений; U (k) – вектор входа, измеренный в предыдущий k-й момент наблюдений; P(k +1) – вспомогательная матрица, определяемая по формуле
45
P(k +1) = P(k) − P(k) U (k)×
×[1+U T (k)P(k)U (k)]−1U T (k)P(k). (1.60)
Итерационные вычисления осуществляются до тех пор, пока уточнение коэффициентов матрицы Am (k +1) не достигнет заданной точности:
Am (k +1) − Am (k) |
|
≤ ε. |
(1.61) |
|
Начальные значения матриц Am (0), P(0) можно определить произвольно, например, так:
A(0) |
= I, |
|
|
|
P(0) |
= |
1 |
I, |
(1.62) |
|
||||
|
|
ε |
|
|
где I – единичная матрица; ε – заданная точность оценивания.
1.4.6.Идентификация динамических систем
Воснове идентификации динамических систем, как правило, лежит описание данных систем в виде передаточной функции [7, 13, 21].
Допустим, динамическая система описана передаточной функцией следующего вида:
W(p)= |
b pm +b pm−1 |
+...+b |
p +b |
|
||
0 |
1 |
m+1 |
m |
. |
(1.63) |
|
a pn + a pn−1 |
+...+ a |
|
||||
|
p + a |
|
||||
|
0 |
1 |
n−1 |
n |
|
На основе передаточной функции (1.63) получаем систему дифференциальных уравнений первого порядка из n-уравнений
dx |
= AX(t)+ BU(t), |
(1.64) |
dt |
|
Y(t)= CX(t),
где Y(t) – вектор выходных переменных; U(t) – вектор входных пере-
менных; X(t) – вектор состояния (внутренние динамические перемен-
ные).
От дифференциальных уравнений переходим к разностным уравнениям
46
V(k +1)= Ф(Т )V(k),
0Y(k)= СV(k),
U
где V – обобщенный вектор, V = ; Ф(bi , ai , T0 )
X
(1.65)
– матрица перехо-
да, параметры которой bi , ai подлежат оцениванию.
В соответствии с (1.57) и (1.65) матрица перехода Ф определяется
как
Ф(Т0)=V T (k +1)V T (k)V(k)V T (k) −1, |
(1.66) |
||||||||
где |
|
|
|
|
|
|
|
|
|
V(k)= |
v1(1) |
v1(2) |
L v1(l −1) |
|
|||||
|
M |
|
|
|
|
M |
, |
|
|
|
|
|
|
|
|
|
|
|
|
|
vn(1) vn |
(2) |
L vn(l −1) |
|
|||||
|
|
v |
(2) |
v |
(3) |
L |
v (l) |
|
|
V(k +1)= |
1M |
1 |
|
|
1M |
. |
|
||
|
|
|
|
|
|
|
|
|
|
|
|
vn |
(2) vn |
(3) |
L vn(l) |
|
Количество измерений l определяется по формуле (1.58).
По полученной матрице перехода можно вычислить матрицу коэффициентов:
∞ (AT )i |
|
AT |
|
Ф(T ) − I |
|
|||
Ф(T0 ) = ∑ |
|
0 |
≈ I + |
0 |
A ≈ |
0 |
. |
(1.67) |
|
i! |
1! |
T |
|||||
i=0 |
|
|
|
|
0 |
|
|
Пример.
Пусть объект управления задан структурной схемой (рис 1.16):
u(t) |
|
x2 (t) |
x1(t) |
K1 |
T |
K2 |
|
− |
|
y(t) |
|
|
|
||
|
|
1T |
|
Рис. 1.16. Структурная схема САУ |
|
47
Необходимо определить параметры системы K1, K2 , T. Определим вектор состояния системы:
r V = x1 .
x2
Динамика процессов описывается с помощью системы дифференциальных уравнений:
drdt = 0,
|
|
|
dx1 |
|
= K |
2 |
x , |
|
|
|
|
|
|
|
|
||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||
|
|
|
dt |
|
|
|
|
|
2 |
|
|
|
|
|
|
|
|
||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
|
|
|
dx2 |
|
= |
|
K1 |
r − |
K1 |
x − |
1 |
x , |
|||||||
|
|
|
dt |
|
|
|
|
|
|
||||||||||
|
|
|
|
|
|
T |
|
|
|
T |
1 |
T 2 |
|||||||
|
|
|
y = x1. |
|
|
|
|
|
|
|
|
|
|
|
|
||||
Основные матрицы модели: |
|
|
|
|
|
|
|
|
|
|
|||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
0 |
|
0 |
|
|
|
0 |
|
|
|
||||||||
A = |
|
0 |
|
0 |
|
|
K |
2 |
,C =[0 1 0]. |
||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
|
K |
|
|
|
K |
|
|
|
|
|
1 |
|
|
|
|||||
|
|
1 |
|
|
− |
|
|
1 |
|
− |
|
|
|
|
|
|
|||
|
|
|
|
|
T |
|
|
T |
|
|
|||||||||
|
T |
|
|
|
|
|
|
|
|
|
|
||||||||
Матрица коэффициентов |
|
A |
имеет размерность n×n = 3×3, по- |
этому для достоверных результатов идентификации (1.58) количество наблюдений l должно быть не менее 10.
В результате серии экспериментов определен вектор состояния:
r(1) |
K |
r(10) |
|
||
V = x |
(1) |
K |
x |
(10) |
, |
1 |
|
|
1 |
|
|
|
(1) |
K |
x2 |
|
|
x2 |
(10) |
|
а матрицы V (k), V (k +1), соответственно, имеют вид
48
r(1) |
K |
r(9) |
|
|
(1) |
|
|
V (k) = x1 |
K |
x1(9) , |
|
|
(1) |
K |
|
x2 |
x2 (9) |
||
|
r(2) |
K r(10) |
|
V (k +1) = |
|
(2) |
|
x1 |
K x1(10) . |
||
|
|
(2) |
|
|
x2 |
K x2 (10) |
Тогда можно определить матрицу перехода:
Ф(T0 ) =V (k +1) V T (k) V (k) V T (k) −1,
а на основе матрицы перехода – и матрицу коэффициентов:
A ≈ Ф(T0 ) − I. T0
Зная численные значения и структуру матрицы коэффициентов А, можно определить параметры системы K1, K2 , T.
1.4.7. Идентификация нелинейных систем
Достаточно развитая теория идентификации нелинейных систем предлагает большое количество различных методов их исследования. Но ни один из методов не является универсальным и характеризуется своей областью применения. Наиболее распространены следующие методы [3, 7, 21]:
1.Метод прямого поиска.
2.Аппроксимация нелинейности.
3.Модель Гаммерштейна.
4.Метод Винера.
5.Двухэтапная процедура. Рассмотрим каждый из методов.
1.Метод прямого поиска.
Нелинейную функцию f(x) преобразуют в линейную функцию fЛ(x) . Далее применяют любой метод идентификации линейныхсистем.
Допустим, что модель объекта имеет вид |
|
||||
y = α |
0 |
xβ1 xβ2 |
, |
(1.68) |
|
|
1 |
2 |
|
|
|
|
|
|
|
|
49 |
где x1,x2 – входные переменные; y – выходная переменная; α0 , β1, β2 –
оцениваемые параметры модели.
Нетрудно заметить, что логарифмирование нелинейного уравне-
ния (1.68) приводит систему к линейному виду: |
|
z = a0 +a1r1 + a2r2 , |
(1.69) |
где z – выходная переменная линейной модели, |
z = lny ; |
r1, r2 – входные переменные линейной модели, r1 = lnx1, r2 = lnx2 ; a0 , a1, a2 – параметры линейнойсистемы, a0 = lnα, a1 =β1, a2 =β2.
Оценивание параметров линейной модели осуществляется явными методами одномерной линейной регрессии (1.52) или итерационным методом (схема с настраиваемой моделью) (1.59).
2. Аппроксимация нелинейностей.
В соответствии с теоремой Вейерштрасса любая непрерывная нелинейная функция можетбыть представлена в виде полинома
y = a |
+ a x + a |
2 |
x2 |
+b х |
2 |
+b х2 |
+K . |
(1.70) |
||
0 |
1 |
1 |
1 |
1 |
2 |
2 |
|
|
Полином (1.70) легко приводится к линейной регрессии, алгоритм оценивания которой был рассмотрен выше.
Аппроксимация с помощью полинома удобна, если значения нелинейной функции получены экспериментально.
3. Модель Гаммерштейна.
В соответствии с моделью Гаммерштейна нелинейная система приводится к виду, представленному на рис. 1.17.
Рис. 1.17. Структурная схема модели Гаммерштейна
Алгоритм идентификации зависит от априорной информации о виде нелинейности F(u(t)) .
Если известна функциональная зависимость F(u(t)) , то при вводе переменной z(t) = F(u(t)) идентификация сводится к определению параметров линейной части модели Гаммерштейна:
50