Построение и моделирование линейных математических моделей систем управления
Выбери формат для чтения
Статья: Построение и моделирование линейных математических моделей систем управления
Загружаем конспект в формате pdf
Это займет всего пару минут! А пока ты можешь прочитать работу в формате Word 👇
1
Курс «Математическое моделирование технических объектов и систем управления».
Лекция 2. Построение и моделирование линейных математических моделей систем
управления
Линеаризованные математические модели первого приближения. Математический аппарат анализа линейных
моделей. Методы синтеза линейных моделей. Краткое описание функций пакета Control System Toolbox.
Методы интегрирования детерминированных и стохастических линейных систем. Наблюдатели полного и
сокращенного порядков. Методы фильтрации входных сигналов. Фильтр Калмана.
Линеаризованные математические модели первого приближения.
Линейные математические модели систем применяются в некоторой малой окрестности
положения равновесия, когда управляющие сигналы ограничены по величине. В этом случае,
исходная нелинейная модель может быть достаточно хорошо описана с помощью правильно
построенной линейной модели. Если линейная модель будет построена плохо, то и результаты
построения системы будут плохими, и использование такой системы управления с нелинейным
объектом приведет к плохому качеству управления. Таким образом, процесс линеаризации
нелинейных уравнений является очень важным при проектировании систем управления. Процедура
линеаризации базируется на представлении нелинейных функций в виде ряда Тейлора,
относительно некоторой точки равновесия. В дальнейшем, принимая малыми возможные
отклонения переменных, мы ограничиваемся только линейными членами в ряде Тейлора, что
позволяет построить линейную модель объекта.
Линейная аппроксимация нелинейных математических моделей.
Рассмотрим нелинейное уравнение вида:
x(t ) F ( x(t ), u(t )) ,
где функция F определяет отображение F : . Точка x называется
изолированной или гиперболической точкой равновесия, если существует такой постоянный
n
m
n
n
входной сигнал u , что F ( x , u ) 0 , и в ближайшей окрестности точки ( x , u ) нет других точек
m
F x x
|u u ) 0 . Если, начальное состояние системы равно x(t0 ) x
x
u , то, очевидно, что x(t ) x для t t0 .
равновесия, а
воздействии
det(
при входном
Определим вариации переменных относительно точки равновесия:
x(t ) x(t ) x , u(t ) u(t ) u ,
x(t ) решение уравнения x(t ) F ( x(t ), u(t )) . Тогда можно записать следующее уравнение:
x(t ) F ( x x(t ), u u(t )) . Применим разложение в ряд Тейлора для правой части уравнения
где
и, оставим только линейные его члены. Тогда можно записать приближенное равенство:
x(t ) F ( x ), u )
Так как F ( x , u ) 0 , то получим:
F x x
F x x
|u u x(t )
|u u u
x
u
F x x
F x x
|u u x(t )
|u u u .
x
u
F x x
B
|u u nm являются постоянными, вещественными
u
x(t )
Матрицы A
F x x
|u u nn
x
и
матрицами, которые также называются матрицами Якоби. Линейная стационарная система
x(t ) A x(t ) B u называется линейной системой первого приближения. Используя такую
2
x(t )
u (t ) относительно точки равновесия ( x , u ) .
модель можно спроектировать систему управления в области вариации
применяя малые управления
ее переменных и
Пример.
Пусть нелинейная система имеет уравнения (при u (t ) 0 ):
x x(6 2 x y)
y y(4 x y)
Здесь имеется две изолированные точки равновесия x1 y1 0 и x2 y2 2 .
В окрестности первой точки x1 y1 0 матрица Якоби равна:
x x 0 6 0
6 4x y
| y 0
.
J ( x1 , y1 )
4 x 2 y
y
0 4
Очевидно, что Якобиан: det(J ( x1 , y1 )) 0 . Тогда можно записать следующую систему первого
x 6 0 x
.
приближения:
y 0 4 y
В окрестности второй точки равновесия x2
y2 2 матрица Якоби равна:
6 4x y
x x2 4 2
J ( x2 , y2 )
|
; det( J ( x2 , y2 )) 0 .
y
4
x 2 y y 2 2 2
x 4 2 x
.
Тогда система уравнений первого приближения имеет вид:
y 2 2 y
В среде Matlab имеется символьная функция
jacobian (.) , которая позволяет вычислять матрицу
Якоби в символьном виде. Обычно, синтаксис команды вызова этой функции имеет вид:
jacobian ( F , v) , где F
- вектор – функция (в символьном виде),
которым вычисляется матрица Якоби. Элемент
(i, j )
Пример.
Вычислим матрицу Якоби векторной функции F ( x, y )
и
x2 y2 2 .
Файл E16_Jacob.
v
- вектор переменных по
результата определяется, как
x(6 2 x y)
y (4 x y )
F (i )
.
v( j )
в окрестности точек x1
y1 0
3
К сожалению, символьное вычисление матрицы Якоби в среде Matlab, не всегда возможно. В
этом случае, для ее вычисления необходимо использовать приближенные методы, где вычисление
производных заменяется их разностной аппроксимацией.
Довольно часто, при линеаризации, возникают более сложные случаи. Предположим, что
исследуемая система (объект) следует вдоль некой программной траектории
xˆ (t ), uˆ (t ) , где
xˆ (t ) F ( xˆ (t ), uˆ (t )) . Определим, как уже указывалось выше, x(t ) x(t ) xˆ (t ) ,
u(t ) u(t ) uˆ (t ) . Тогда мы можем записать уравнения:
F x xˆ (t )
F x xˆ (t )
xˆ (t ) x(t ) F ( xˆ (t ), uˆ (t ))
|u uˆ (t ) x(t )
|u uˆ (t ) u .
x
u
Так как xˆ (t ) F ( xˆ (t ), uˆ (t )) , то получим:
F x x
F x x
x(t )
|u u x(t )
|u u u .
x
u
F x xˆ (t )
F x xˆ(t )
В этом случае матрицы Якоби A(t )
|u uˆ (t ) nn и B(t )
|u uˆ(t ) nm зависят от
x
u
параметра t (или от времени). То есть линеаризованная модель первого приближения не является
LTI моделью. Такие модели называются линеаризованными моделями первого приближения
относительно программной траектории xˆ (t ), uˆ (t ) .
Основные модели представления детерминированных линейных систем.
Представление линейных систем в виде моделей состояния.
Представление детерминированной линейной системы в виде модели уравнений состояния
можно записать, в общем случае, в виде:
x A(t ) x B(t )u, y C (t ) x D(t )u ,
где x n , u m , y p , A(t ) nn , B(t ) nm , C (t ) pn , D pm
Модель состояния стационарной линейной системы имеет, соотвественно, форму
x Ax Bu , y Cx Du , где x n , u m
,
y p , A nn , B nm , C pn ,
D pm .
Операторное представление линейных стационарных систем.
Прямое преобразование Лапласа функции f (t ) определяется с помощью соотношения:
L{ f (t )} e st f (t )dt F ( s)
,
4
где обозначение
L{ f (t )} используется для обозначения интеграла Лапласа для краткости.
Функция: F ( s) является функцией комплексной переменной s i . Функция f (t )
называется оригиналом, а функция F ( s) - изображением.
Преобразование Лапласа определено для функции f (t ) , удовлетворяющей следующим
условиям:
1.
Для значений
2.
Для значений
t0
t0
3.
Функция
ограничена в росте, при
M
f (t )
функция имеет значение
функция
f (t )
f (t ) 0 ;
имеет конечное число точек разрыва непрерывности;
t , экспонентой | f (t ) | M et , где | | ,
- неотрицательная константа
Обратное преобразование Лапласа, при заданном изображении F ( s) определяется
следующим образом:
i
1
f (t ) L {F ( s)}
F ( s)e st ds
2 i i
больше, чем реальная часть сингулярности изображения F ( s) .
1
где
,
Рассмотрим LTI (Linear Time Invariance) модель системы, которая описывается линейным
дифференциальным уравнением n - го порядка с нулевыми начальными условиями:
x( n) an 1 x( n 1) ... a1 x a0 x bmu ( m) bm1u ( m1) ... b1u b0u ;
x(0) 0 , x(0) 0 , …, x( n 1) (0) 0 , u(0) 0 , u(0) 0 ,…, u ( m1) (0) 0 .
Применив преобразование Лапласа, получим, что:
bm s m bm 1s m 1 ... b0
X (s) n
U (s) .
s an 1s n 1 ... a0
Таким образом, каждой паре вход - выход {x(t ), u (t )} для LTI модели можно поставить в
X (s)
соответствие пару { X ( s),U ( s)} на s – плоскости. При этом соотношение
будет
U (s)
определяться дробно – рациональной функцией от параметра s в виде:
X ( s) bm s m bm 1s m 1 ... b0
n
U ( s)
s an 1s n 1 ... a0
,
которая называется передаточной функцией.
Определение.
Передаточной функцией системы управления называется отношение изображения X ( s)
выходного сигнала к изображению
U ( s)
входного сигнала при нулевом состоянии системы
(нулевых начальных условиях).
Для LTI модели системы передаточная функция представляет собой дробно – рациональную
функцию параметра s в виде:
bm s m bm 1s m 1 ... b0
W ( s) n
s an 1s n 1 ... a0
.
Такое представление системы в среде Matlab называется tf – формой.
5
Наивысший порядок n знаменателя передаточной функции определяет порядок системы. Для
физически реализуемых систем порядок числителя m не должен превышать порядок знаменателя,
то есть m n . В этом случае передаточная функция называется правильной. Если выполняется
соотношение m n , то такая функция называется строго правильной. Величина называется
относительной степенью передаточной функции.
Модель системы в виде передаточной функции (tf – модель) реализуется в среде Matlab с
помощью следующей команды:
W tf [num, den] ,
где num [bm , bm 1 ,..., b0 ] ,
den [1, an1 , an 2 ,..., a0 ] .
Пример.
Передаточную функцию (точнее передаточную матрицу) для систем с многими входами и
выходами можно получить из уравнений линеаризованной модели первого приближения системы
при нулевых начальных условиях. Пусть задано уравнения состояния системы вида:
x Ax Bu , y Cx Du , где x n , u m
,
y p .
Такая форма записи системы в среде Matlab называется ss – моделью. Применяя
преобразование Лапласа к левой и правой части уравнения найдем:
Y (s) (C (s) B D) U (s) ,
(s) ( sE A)1 . Отсюда получим, что передаточная матрица системы равна:
W (s) C (s) B D .
где:
Пример.
1 1 1
0
0 2
x
u , y
x .
2 3 0 2
1 0
Характеристический полином системы имеет вид: det(sE A) (s 1)(s 2) . Тогда преобразование Лапласа
s 3 1
1
. Так как D 0 , то
переходной матрицы системы можно записать в виде: ( s )
( s 1)(s 2) 2 s
Пусть задано уравнение состояния следующей LTI системы: x
4
( s 1)( s 2)
передаточная функция (матрица) равна: W ( s ) C ( s ) B
s3
( s 1)( s 2)
Вычислим передаточную функцию W (s) с помощью команд Matlab.
Файл E4_TrFunc
4
s 2 W11( s) W12 ( s)
.
1 W21( s) W22 ( s)
s 2
6
Анализ линейных стационарных систем при операторном представлении.
Частотные характеристики
Под термином «частотный выход системы» понимают установившееся значение выхода системы
при подаче на ее вход синусоидального сигнала с различной частотой.
Рассмотрим установившийся режим при синусоидальном входе для устойчивых, линейных,
стационарных систем, которые описываются с помощью передаточной функции. Напомним, что
устойчивые системы имеют все полюса, расположенные в левой части комплексной плоскости, то
есть вещественные части полюсов являются отрицательными. Обозначим передаточную функции
системы, как W ( s ) , а, соответственно изображение входа через L{x(t )} X ( s) и выхода, как
L{ y(t )} Y (s)
. Если вход
x(t )
является синусоидальным сигналом некоторой частоты
, то
выход системы в установившемся режиме будет также синусоидальным с той же частотой, но с
другой амплитудой и фазой сдвига. Пусть входной сигнал имеет вид:
x(t ) A sin( t ) ,
где частота
измеряется в rad/sec. Представим передаточную функцию
полюсами в виде:
W ( s)
B( s )
B( s )
A( s) ( s s1 )( s s2 ) ... ( s sn )
Изображение выхода системы тогда равно:
W ( s)
с простыми
7
Y ( s) W ( s) X ( s) .
Декомпозицию изображения выхода, в виде разложения на простые дроби, можно
представить как:
Y (s) W (s) X (s) W (s)
где
c, c
полюса
t
dn
d1
A
c
c
...
s 2 2 s i s i s s1
s sn
,
и
d i являются, в общем случае, комплексными константами. Для устойчивой системы
s1 , s2 ,..., sn имеют отрицательные вещественные части. Поэтому, при больших значениях
, составляющие выхода
es1t , es2t ,..., esnt
можно принять равными нулю. Таким образом,
установившиеся значение можно принять равным:
yss (t ) c eit c eit ,
где вычеты c и c могут быть вычислены следующим образом:
2 A
A W (i )
,
c lim ( s i ) W ( s) 2
2
s i
s
2i
2 A
A W (i )
c lim ( s i ) W ( s) 2
.
2
s i
s
2i
Так как W (i ) является комплексной величиной, можно, используя формулу Эйлера,
записать:
W (i ) | W (i ) | ei ,
где
arctan(
Im(W (i ))
).
Re(W (i ))
С другой стороны имеет место равенство:
W (i ) | W (i ) | ei . Тогда установившийся
выход равен:
ei ( t ) ei ( t )
yss (t ) A | W (i ) |
A | W (i ) | sin( t ) H sin( t ) ,
2i
где переменная H A | W (i ) | . Отсюда можно сделать вывод, что установившийся выход
также представляет собой синусоидальный сигнал той же частоты, что и входной синусоидальный
сигнал. Но амплитуда выходного сигнала и фаза (угол сдвига) будут, в общем случае, отличными от
входного сигнала.
Функция
W (i )
называется частотной передаточной функцией или амплитудно-фазовой
характеристикой (АФЧХ). Функция
| W (i ) |
амплитудно-частотной характеристикой системы
8
(АЧХ). Функция
() arg(W (i )) arctan(
Im(W (i ))
)
Re(W (i ))
называется фазо-частотной
характеристикой (ФЧХ) системы или частотным аргументом.
Пример.
Задана передаточная функция системы в форме Боде:
W ( s)
k
, где k , T
Ts 1
вещественные константы.
x(t ) A sin( t ) . Тогда амплитудно-фазовая характеристика имеет
k
k
вид: W (i )
. Амплитудно-частотная характеристика равна: | W (i ) |
и
2 2
1 i T
1 T
фазо – частотная характеристика: () arctan(T ) . Установившийся выходной сигнала будет
k
sin( t arctan(T )) .
описываться соотношением: yss (t ) A
1 T 2 2
Пусть задан синусоидальный вход:
Так как передаточная функция
W ( s)
определяется как:
W ( s) e st w(t )dt ,
где
w(t ) - импульсная функция, то формально можно записать, что: W (i ) eit w(t )dt .
Таким образом, амплитудно-фазовая характеристика является Фурье преобразованием от
импульсной функции.
Так как имеет место равенство:
eit cos( t ) i sin( t ) , то тогда можно записать
следующие соотношения:
Re{W (i )} w(t ) cos( t )dt , Im{W (i )} w(t )sin( t )dt .
Отсюда несложно получить:
Re{W (i )} Re{W (i )} , Im{W (i )} Im{W (i )}.
Таким образом, амплитудно-фазовая характеристика для частот из диапазона 0
является зеркальным отражением этой же характеристики для диапазона частот 0
относительно вещественной оси.
Если все полюса и нули передаточной функции W ( s ) имеют отрицательную вещественную
часть, то такая передаточная функция называется минимально – фазовой. Если передаточная
функция имеет хотя бы один нуль или один полюс с положительной вещественной частью, то такая
система называется неминимально – фазовой.
Если система является минимально – фазовой, тогда амплитудно – фазовая характеристика
W (i ) будет полностью определяться амплитудно – частотной характеристикой | W (i ) | .
Если система является неминимально – фазовой, то функция
| W (i ) |
, так и
()
.
W (i )
будет определяться, как
9
Правильная передаточная функция ( m n ) будет ограничена по амплитуде с увеличением
частоты, то есть | W (i ) | . Строго правильная передаточная функция ( m n ) будет
стремиться к нулю при увеличении частоты, то есть | W (i ) | 0 .
Пример.
1 Ts
, где 0 T T1 .
1 T1s
Передаточная функция Wa ( s ) является минимально – фазовой, а передаточная функция Wb ( s ) Пусть заданы две передаточные функции:
Wa ( s)
1 Ts
1 T1s
и
Wb ( s)
неминимально - фазовой.
1 T 2 2
Амплитудно-частотная характеристика обоих систем равна: A() | Wa (i ) || Wb (i ) |
1 T12 2
На
s
- плоскости представлены нули и полюса
Wa ( s) and Wb ( s) .
Фазо – частотные характеристики систем равны:
a () arctan(
(T1 T )
(T1 T )
,
)
(
)
arctan(
).
a
1 2T1T
1 2T1T
Построение частотных характеристик.
Амплитудно-частотная характеристика является комплексной функцией и зависит от параметра
. Имеется несколько вариантов графического представления такой характеристики, к которым
относятся следующие:
- Диаграммы Боде или логарифмические частотные характеристики;
- Диаграмма Найквиста;
- Диаграмма Николса.
.
10
Диаграммы Боде. Логарифмические частотные характеристики.
Диаграмма Боде состоит из двух графиков. Первый график представляет собой зависимость
логарифма амплитудно – частотной характеристики (ЛАЧХ), а второй график отражает зависимость
фазо-частотной характеристики (ЛФЧХ) от частоты входного синусоидального сигнала. Оба графика
строятся в логарифмической шкале частот.
Логарифмическая амплитудно-частотная характеристика строится в следующем масштабе:
Log ( A()) является
децибел, которая имеет аббревиатуру dB. Таким образом, изменение амплитуды в 10 раз при
LogA() 20log10 | W (i ) |
. То есть единицей измерения функции
изменении частоты в 10 раз будет приводить к изменению наклона логарифмической амплитудночастной характеристике на графике на 20dB / decade . Хотя график функции Log ( A())
невозможно построить во всем диапазоне частот, так как
log10 0
, это не создает серьезных
трудностей. Это обусловлено тем, что исследование реальных систем ведется на частотах
отличных от нуля, то есть 0 .
Передаточную функцию можно представить в виде произведения элементарных множителей
вида:
W ( s ) W j ( s )
,
j
Где
W j ( s)
(Ts 1) ,
- передаточные функции некоторых элементарных звеньев, таких как:
1
1
2
, (T s 2T 1) , 2
T s 2T 1
Ts 1
T1s 1
T s 2T2 1
и
2
2
K , sv ,
1
sv
,
. Если известны
характеристики таких элементарных передаточных функций, то по ним легко определить и
соответствующие характеристики передаточной функции W (s) . Действительно, логарифмические
амплитудно-частотную и фазо - частотную характеристики (ЛАЧХ и ЛФЧХ) можно вычислить с
помощью следующих соотношений:
Log ( A()) 20log10 | W j (i ) | 20 log10 | W j (i ) | Log ( Ai ()) ,
j
j
j
() arg(W (i )) arg(W j (i )) arg(W j (i )) .
j
Если ЛАЧХ
Log ( A())
и ЛФЧХ
()
j
известны, то их можно соответствующим образом
корректировать, добавляя в исходную передаточную функцию в качестве множителей передаточные
функции некоторых элементарных звеньев.
Построение диаграмм Боде в среде Matlab.
Функция bode(.) средыMatlab позволяет строить диаграммы Боде для LTI моделей, заданных в
виде передаточной функции. Данная функция имеет следующий синтаксис:
bode(W ) - автоматическое построение диаграммы Боде;
bode(W ,(m , M )) - построение диаграммы Боде в диапазоне частот (m , M ) ;
bode(W , ) - построение точек диаграммы Боде для частот заданных вектором ;
[ A, , ] bode(W ) - вывод данных диаграммы Боде в массив данных без построения графика;
bode(W1 ,' ',W2 ,'*',W3 .' : r ') - построение единой диаграммы Боде для различных систем;
11
bode( A, B, C, D)
- построение диаграммы Боде для системы, заданной моделью состояния:
y Cx D , где x n , y, u
Пример.
Задана передаточная функция:
W ( s)
s 8
s( s 2 0.2s 4)( s 1)( s 3)
. Построить диаграмму Боде.
Файл E10_BodeDg.
Построение асимптотических логарифмических характеристик.
Пусть передаточная функция системы представляет собой произведение передаточных функций
элементарных звеньев, то есть:
W ( s) W j ( s) . Тогда ЛАЧХ и ЛФЧХ такой передаточной
j
log
функции можно представить в виде: Log ( A()) 20
10
| W j (i ) | Log ( Aj ())
j
,
j
() arg(W j (i )) j () . Для оценки правильности расчетов, а также для
j
j
предварительного анализа такой системы, можно построить асимптотическую аппроксимацию ЛАЧХ
и ЛФЧХ системы с передаточной функцией W ( s ) . Для этого можно использовать следующие
правила:
1. Представить передаточную функцию
W ( s)
в форме Боде;
12
2. Определить значение коэффициента усиления в логарифмической шкале
также угловые (собственные) частоты элементарных звеньев
упорядочить эти частоты по возрастающему порядку, то есть
3. Построить точку ( 1,
20log10 K )
4. Провести через точку ( 1,
на диаграмме Боде
20log10 K )
j
1
Tj
20log10 K
,а
; при этом
1 2 .... ;
(log10 , dB) ;
первую асимптотическую линию с наклоном
20v
dB/decade, где v определяет порядок идеальных интегральных или дифференциальных
звеньев
sv
до угловой частоты
1 .
5. Построить вторую асимптотическую линию для диапазона частот
через точку ( 1,
1 2
, проходящую
20v плюс наклон асимптоты первого
элементарного звена, имеющего угловую (собственную) частоту, совпадающую с 1 .
6. Построить асимптотические линии для последующих участков j 1 j , пользуясь
правилом пункта 5; наклон j - ой асимптоты определяется наклоном итоговой ЛАЧХ в точке
j 1 плюс наклон асимптоты j - го элементарного звена;
20log10 K )
и имеющую наклон
7. Построить последнюю асимптоту, проходящую через последнюю угловую (собственную)
частоту, и продолжить ее до .
Пример.
100( s 1)
. Таким
s(10s 1)(0.01s 2 0.1s 1)
1
1
1 0.1, 2 1, 3
10 . ЛАЧХ системы
10
0.1
Пусть задана передаточная функция в форме Боде:
образом, получим следующие угловые частоты:
W ( s)
можно записать в виде:
LogA() 40 20log10 2 1 20log10 20log10 (10) 2 1
20log10 (1 0.012 )2 (0.1) 2
График асимптотической ЛАЧХ имеет вид:
Для вычисления ЛАЧХ и ЛФЧХ можно использовать следующий script – файл.
13
Файл E11_BodeDg.
Диаграмма Найквиста.
Диаграмма Найквиста для АФЧХ
W (i )
строится в координатах Re{W (i )}, Im{W (i)}
при изменении частоты от 0 до ( в среде Matlab обычно принимается, что изменяется от
до ). График W (i ) на диаграмме представляет кривую, которую описывает конечная
точка вектора ( Re{W (i )}, Im{W (i)} ) при изменении частоты
.
Таким образом, проекции точки на оси координат определяют вещественную и мнимую части
АФЧХ W (i ) .
14
Пример.
1
1
e arctan( T ) . Тогда можно
2
2
1 i T
1 T
1
записать: W (i ) U () i V () , где U () Re{W (i )}
,
1 2T 2
1 2
1 2
T
2
(
U
(
)
)
(
V
(
))
(
) . Таким
.
Можно
заметить,
что:
V () Im{W (i )}
2
2
1 2T 2
образом, в координатах (U ,V ) график W (i ) представляет собой окружность с центром в точке:
1
1
U , V 0 и радиусом .
2
2
Пусть АФЧХ системы имеет вид:
W (i )
Пример.
Задана передаточная функция:
W (i )
W ( s)
T
1
i
.
1 2T 2
(1 2T 2 )
1
. Тогда АФЧХ системы имеет вид:
s(Ts 1)
Отсюда получим, что асимптотой к графику АФЧХ при
0
является вертикальная линия, проходящая через точку (T ,0) оси абсцисс.
Построение диаграммы Найквиста в среде Matlab.
Для построения диаграммы Найквиста в среде Matlab используется функция nyquist (.) .
Команда ее выполнения может иметь следующий синтаксис:
nyquist (W ) - автоматическое (с параметрами по умолчанию) построение диаграммы Найквиста;
[Re, Im, ] bode(W )
- вычисляет координат точек диаграммы Найквиста без построения
графика;
nyquist (num, den)
W ( s)
- построение диаграммы Найквиста для передаточной функции
numerator polynomial in s
deno min ator polynomial in s
, заданной числителем и знаменателем;
15
nyquist (num, den, ) - построение диаграммы Найквиста для фиксированных частот, заданных
вектором ;
nyquist ( A, B, C, D, iu) - построение диаграммы Найквиста для модели состояния x Ax Bu ;
y Cx D , где x n , y, iu , параметр iu определяет вход системы.
Пример.
Задана передаточная функция:
W (s)
1
s 0.8s 1
2
. Построить диаграмму Найквиста.
Файл E12_Nyquist.
Анализ линейных систем при представлении в виде уравнений состояния.
Рассмотрим автономную линейную систему следующего вида:
где
x Ax Bu; y Cx Du; x n , u m , y p ,
A, B, C, D - действительные матрицы размерностей (n n),(n m),( p n),( p m)
соответственно.
Обычно предполагается, что начало координат x0 0 является положением равновесия
для свободной, или неуправляемой системы, для которой u (t ) 0 . При этом считается, что более
общая система с ненулевым положением равновесия
(x,u )
во многих случаях может быть
приведена к такому же виду с помощью параллельного переноса осей координат. Для такого
приведения необходимо уметь решать линейное уравнение, которое определяет точку x 0 при
значении u 0 , то есть уравнение вида Ax b .
Переход от уравнения состояния к передаточной функции.
Пусть модель состояния системы представляется соотношениями:
x Ax Bu; y Cx Du; x n , u m , y p .
Эквивалентная tf модель будет иметь матричную передаточную функцию вида:
16
где
W (s) C (s E A)1 B D ,
E - единичная диагональная матрица размера (n n) . Для вычисления такой передаточной
функции можно использовать метода Леверье – Фадеева. Обозначим характеристический полином
системы в виде:
P(s) det(s E A) n s n n 1s n 1 ... 0 ,
где
n 1 . Тогда обратная матрица (резольвента): ( s E A)1
может быть записана в форме:
Qn 1s n 1 Qn 2 s n 2 ... Q0
( s E A)
P( s )
1
где
Qn 1 E
,
. Метод Леверье – Фадеева позволяет одновременно вычислять коэффициенты
характеристического полинома
i , i 0,1, 2,..., n 1 и матрицы Q j , j 0,1, 2,..., n 2
используя следующие рекуррентные соотношения / 4 /:
Qn 1 E, n 1;
Fn 1 A Qn 1 , n 1 tr ( Fn 1 ) / 1; Qn 2 Fn 1 n 1 E ;
Fn 2 A Qn 2 , n 2 tr ( Fn 2 ) / 2; Qn 3 Fn 2 n 2 E
.........
F1 A Q1 , 1 tr ( F1 ) / (n 1); Q0 F1 1 E ;
F0 A Q0 , 0 tr ( F0 ) / n; F0 0 E 0 .
n
Здесь:
tr ( F ) Fii .
Последнее соотношение может быть использовано для вычисления
i 1
1
обратной матрицы A
Q0
0
.
Преобразования математических моделей в среде Matlab.
Matlab позволяет эффективно преобразовывать различные модели LTI систем . Команда
[ A, B, C, D] tf 2ss(num, den) позволяет преобразовать tf модель в модель состояния. При этом,
передаточная функция задается в виде векторов коэффициентов числителя и знаменателя:
W ( s)
numerator polynomial in s
.
deno min ator polynomial in s
Обратное преобразование от модели состояния (ss модели) к модели передаточной функции
(tf модели) осуществляется командой
[num, den] ss 2tf ( A, B, C, D, iu) , где параметр iu должен
указываться для систем, имеющих более, чем один вход. Если система имеет только один вход, то
команда может иметь синтаксис:
[num, den] ss 2tf ( A, B, C, D) .
17
Пример.
Файл E8_TFinSS. Преобразование tf – модели в ss – модель, и наооборот.
18
Внутренняя устойчивость. Устойчивость по Ляпунову.
Понятие внутренней устойчивости или устойчивости по Ляпунову основывается на оценке
количественного поведения систем при нулевом входном сигнале, то есть исследуется устойчивость
собственного движения системы исключительно от ненулевых начальных условий.
Рассмотрим нелинейную систему:
x F (t , x) , F (t ,0) 0 , x(t0 ) x0 , x n ,
которая имеет изолированную точку равновесия в начале координат. Будем считать, что функция
F (.)
является непрерывной и имеет непрерывные частные производные
t t0 . Тогда решение системы x x(t ) n
можно представить как траекторию в
пространстве в зависимости от независимого параметра
Определение.
Решение (t ) системы x
F (t , x)
F
, k 1, 2,..., n
xk
при
n - мерном
t.
называется устойчивым по Ляпунову при
t (или
просто устойчивым), если:
1. Для любого 0 существует такая величина
2.
(, t0 ) 0 , что все решения x(t )
системы x F (t , x) (включая (t ) ) удовлетворяют условию || x(t0 ) (t0 ) || ;
Все решения определены на всем интервале [t0 , ) и удовлетворяют неравенству
|| x(t ) (t ) || для t [t0 , ) .
Интуитивно понятно, что если решение
(t )
является устойчивым, то другие
решения x
x(t ) , которые были достаточно близки к решению
t [t0 , )
лежать в некоторой заданной узкой трубке, построенной вокруг решения
F (t ,0) 0
в момент
t0
, будут для
(t )
(t ) 0 (называемое также
невозмущенным движением системы) будет устойчивым, если для любого t0 существует такая
величина (, t0 ) 0 , из неравенства || x(t0 ) || будет вытекать следующее соотношение
|| x(t ) || для t [t0 , ) .
Следует понимать, однако, что если решение (t ) является устойчивым, то это не
значит, что все решения будут ограничены. Таким образом, если все решения системы x F (t , x)
В частности, когда
то тривиальное решение
являются ограниченными, то это не значит, что они устойчивы, и наоборот.
Если можно выбрать величину 0 такую, что она не будет зависеть от параметра t 0 , то есть
() , тогда такое решение (t )
называется равномерно устойчивым.
19
Определение.
Решение
(t ) называется неустойчивым по Ляпунову, если существует такое
решение x x (t ) и величины 0 , 0 такие, что || x (t0 ) (t0 ) || , и найдется такой
момент t1 t1 () t0 , при котором будет выполняться неравенство || x (t1 ) (t1 ) || .
В частности, когда F (t ,0) 0 , то тривиальное решение (t ) 0 будет неустойчивым, если
существует такое решение
x x (t )
и величины 0 , 0 такие, что || x(t0 ) || , и найдется
такой момент t1 t1 ( ) t0 , при котором будет выполняться неравенство
Определение.
Решение (t ) называется асимптотически устойчивым при
|| x(t1 ) || .
t , если выполняются
следующие условия:
- Решение (t ) является устойчивым;
-
(t0 ) 0 , t0 , что все решения x x(t ) , которые
удовлетворяют условию || x(t0 ) (t0 ) || , имеют свойство lim || x(t ) (t ) || 0
Существует такая величина
t
В частности, если решение
условия
|| x(t0 ) ||
(t ) 0
является асимптотически устойчивым, то при выполнении
будет выполняться свойство
lim || x(t ) || 0 .
t
Устойчивость LTI моделей систем.
Рассмотрим однородную LTI модель системы:
x Ax , x(t0 ) x 0 , x n .
Общее решение системы можно записать в виде x(t ) e
A(t t 0 ) 0
x . Построим преобразование
системы в эквивалентную жорданову каноническую форму с помощью преобразования z Tx ,
где T R n n несингулярная матрица. Тогда можно записать следующее соотношение:
A T 1 ( J1 (1 ) ... J s (s ))T T , где 1 ,..., s , s n - собственные значения матрицы A с учетом их
кратности. Обозначим через m1 ,..., ms порядки жордановых клеток (ячеек). Отсюда можно записать
следующее соотношение:
exp{(t t0 ) I p ( p )} e
p (t t 0 )
[ Em p
(t t0 )
(t t0 )
(t t0 )
I1, p
I 2, p ...
I m 1, p ] ,
1!
2!
(m p 1)! p
где I k , p ( I p ) k , I p - единичная косая матрица размера m p m p , Em p - единичная диагональная
матрица размера m p m p .
20
Теорема.
Если однородная LTI модель системы является x Ax устойчивой по Ляпунову, то тогда все
собственные значения 1 ,..., n матрицы A R n n имеют неотрицательные вещественные части, то
есть Re(k )) 0; k 1,2,..., n . Кроме того, алгебраическая кратность всех собственных значений,
которые имеют нулевую вещественную часть, должна совпадать с их геометрической кратностью.
Отсюда вытекает, что если LTI модель устойчива, то она будет равномерно устойчива.
Действительно, если система x Ax устойчива, то имеет место оценка || e At || M , где M является
константой t t0 . Тогда можно записать неравенство: || x |||| e
A(t t 0 )
|| || x 0 || M || x 0 || . Если
выбрать: || x 0 ||
, тогда получим: || x || . То есть 0 и не зависит от t 0 . Таким образом
M
тривиальное решение x(t ) 0 является равномерно устойчивым.
Теорема.
Однородная LTI модель системы x Ax будет асимптотически глобально устойчивой тогда,
когда все собственные значения матрицы A R n n имеют отрицательные вещественные части, то
есть Re(k ) 0; k 1,2,..., n .
Очевидно, что все решения LTI моделей систем x Ax являются или устойчивыми, или
неустойчивыми одновременно. К сожалению, такого нельзя сказать о нелинейных системах, где
некоторые решения могут быть устойчивыми, а другие неустойчивыми.
Устойчивость по Ляпунову неоднородных LTI моделей систем.
Рассмотрим неоднородную LTI модель системы:
x Ax Bu , x R n , u R m ,
где вход u u (t ) является ограниченной функцией, и который воздействует на систему на
u (t ), t [t0 , T ]
, || u(t ) || um . Общее
,
t
T
ограниченном интервале [t0 , T ] , T , то есть u (t )
решение системы имеет вид: x e
t
x e A(t )u ( )d . Предположим, что соответствующая
A(t t 0 ) 0
t0
однородная система x Ax является асимптотически устойчивой, то есть max(Re (k )) 0 , где
k
k , k 1,2,..., n собственные значения матрицы A . Тогда можно записать следующие соотношения:
|| x |||| e
t
A(t t 0 ) 0
x || || e
A(t )
(t t 0 )
u ( )d || e
|| x || u m || e A(t ) || d . Так как
t0
t
T
A(t )
|| d | e (t ) d | et | [
|| e
t0
t0
T
t0
e t 0 e T
] | , и 0 , то отсюда получим || x(t ) || 0 при t .
Можно также показать, что если собственные значения с максимальными вещественными
частями имеют одинаковую алгебраическую и геометрическую кратности, то если однородная
система является устойчивой, то и неоднородная соответствующая система с ограниченным входом
будет также устойчивой.
21
Таким образом, устойчивость неоднородной LTI модели системы с ограниченным входным
воздействием будет полностью определяться устойчивостью соответствующей однородной LTI
модели.
Оценка устойчивости автономной нелинейной системы по устойчивости линейной системы первого
приближения.
Пусть задана нелинейная автономная система x F (x) , F (0) 0 . Построим в окрестности
начала координат LTI модель первого приближения вида x Ax и найдем связь между ее
устойчивостью и устойчивостью исходной нелинейной системы.
Теорема (1-я теорема Ляпунова об устойчивости).
Если все собственные значения матрицы A линеаризованной стационарной системы первого
приближения лежат в левой части комплексной плоскости, то невозмущенное движение
соответствующей (то есть тривиальное решение) соответствующей (порождающей) нелинейной
автономной системы x F (x) будет устойчивым.
Если, по меньшей мере, имеется одно собственное значение матрицы A лежит в правой
части комплексной плоскости, то тогда невозмущенное движение соответствующей нелинейной
автономной системы x F (x) будет неустойчивым.
В случае, когда имеются собственные значения, лежащие на мнимой оси, но нет, ни одного
собственного значения, лежащего в правой полуплоскости, называется критическим. В критическом
случае нельзя оценить устойчивость исходной нелинейной системы на основе математической
модели линеаризованной стационарной системы первого приближения и требуется проведение
дополнительного исследования.
Пример.
Нелинейная система задана уравнениями:
x x(6 2 x y)
y y(4 x y)
y1 0 , x2 y2 2 . Построим линеаризованную
x 6 0 x
.
систему первого приближения для первой точки равновесия. Система имеет вид:
y 0 4 y
В данной системе имеется две точки равновесия x1
Собственные значения матрицы
6 0
равны 1 6 , 2 4 . Таким образом, линеаризованная система
0 4
является неустойчивой. Тогда по теореме Ляпунова исходная нелинейная система будет также неустойчива в
окрестности точки x1
y1 0 .
Построим теперь линеаризованную систему первого приближения в окрестности второй точки равновесия:
x 4 2 x
. Собственные значения матрицы
y 2 2 y
4 2
равны 1,2 3 5 . То есть все
2 2
собственные значения лежат в левой полуплоскости комплексной плоскости. Тогда в соответствии с теоремой
Ляпунова исходная нелинейная система будет устойчива в окрестности точки x1
Пример.
Заданы две нелинейные автономные системы:
x1 x2 x1 ( x12 x22 ); x2 x1 x2 ( x12 x22 ) ;
x1 x2 x1 ( x12 x22 ); x2 x1 x2 ( x12 x22 ) .
y1 2 .
22
Обе системы имеют точку равновесия
первого приближения: x1
( x1 , x2 ) (0,0)
и одни и те же уравнения линеаризованной системы
0 1
имеет собственные значения 1,2 0 .
1 0
x2 ; x 2 x1 . Матрица системы
То есть, линеаризованная система имеет точку равновесия ( x1 , x2 ) (0,0) , типа центр.
Введем полярные координаты x1
r sin и x2 r cos для обеих нелинейных систем. Тогда уравнения
первой системы в полярных координатах будут:
r r 3 ; 1 , а уравнения второй системы: r r 3 ; 1 .
Отсюда легко построить фазовые портреты обеих нелинейных систем и линеаризованной системы.
Если линеаризованная система первого приближения будет нестационарной, то необходимо
использовать другие методы оценки устойчивости.
Внешняя устойчивость. BIBO устойчивость.
Внешняя устойчивость еще называется устойчивостью по входу – выходу и определяет
зависимость выхода системы, при воздействии на нее ограниченного входа. То есть, это
принципиально отличается от внутренней устойчивости или устойчивости по Ляпунову, где в
качестве возмущения рассматривается только отклонение в начальных условиях. Так как при
внешней устойчивости исследуется ограниченность выхода при подаче на вход ограниченного
сигнала (Bounded -Input- Bounded-Output или BIBO), то часто такой тип устойчивости, особенно в
англоязычной литературе называют BIBO устойчивостью. Однако, хотя внутренняя и внешняя
устойчивости по основным признакам отличаются, они являются в определенном смысле
связанными, что будет показано позднее.
Говорят, что вход u (t ) является ограниченным, если u (t ) не стремится по амплитуде к
или, что эквивалентно, существует такая константа u m , что || u(t ) || um для t 0 .
Будем говорить, что система называется BIBO устойчивой, если для каждого ограниченного
входа существует свой ограниченный выход.
Данный тип устойчивости определяется при нулевом состоянии системы, то есть, когда все
начальные состояния ее координат равны нулю. Другими словами, система будет BIBO устойчива,
если существует такая положительная константа , что при подаче на вход системы ограниченного
сигнала u (t ) , при нулевом состоянии системы, выход системы будет удовлетворять неравенству
sup || y(t ) || sup || u (t ) || .
t 0
t 0
BIBO устойчивость SISO LTI систем.
Рассмотрим SISO LTI систему, заданную передаточной функцией W (s) . Тогда изображение
выхода системы будет равно Y (s) W (s)U (s) , где U (s) L(u(t )) , | u(t ) | um , или во временной
t
области y (t ) w(t )u ( )d .
23
Теорема.
t
SISO LTI система y (t ) w(t )u ( )d будет BIBO устойчивой, тогда и только тогда, когда ее
импульсная функция w(t ) будет абсолютно интегрируема на интервале [0, ) , или
| w(t ) | dt M , где M
- некоторая положительная константа.
Рассмотрим свойства BIBO устойчивости для SISO LTI систем.
Теорема.
Если система с импульсной функцией w(t ) является BIBO устойчивой, то тогда при t ,
получим:
1. Если вход равен u (t ) a для t 0 то выход определяется соотношением W (0) a ;
2. Если вход равен u(t ) sin(0t ) для t 0 , то выход определяется
соотношением | W (i 0 ) | sin(0t arg W (i 0 )) , где: W ( s )
w(t )e
st
dt .
Теорема.
SISO система, которая задается правильной дробно – рациональной передаточной функцией
W (s) будет BIBO устойчивой, тогда и только тогда, когда каждый полюс W (s) имеет отрицательную
вещественную часть или, что эквивалентно, лежит в левой части комплексной плоскости (s –
плоскости).
BIBO устойчивость MIMO LTI систем.
Теорема.
MIMO LTI система, с правильной дробно - рациональной передаточной матрицей G( s) ( g ij ( s)) ,
i 1,2,..., p; j 1,2,..., m , является BIBO устойчивой, тогда и только тогда, когда каждый полюс
каждого элемента g ij (s) имеет отрицательную вещественную часть.
Рассмотрим связь между полюсами передаточной матрицы G(s) и собственными значениями
C [ Adj ( sE A)]B
D , то очевидно, что каждый полюс G(s) является
det(sE A)
собственным значением матрицы A . То есть, если каждое собственное значение имеет
матрицы A . Так как G ( s)
отрицательную вещественную часть, то и каждый полюс будет иметь отрицательную вещественную
часть и, следовательно, система будет BIBO устойчивой. Однако, при вычислении передаточной
матрицы возможно сокращение полюсов и нулей, в том числе и положительных. Поэтому, если
исходная матрица A имеет собственные значения с нулевой или даже положительной частями, то
система x Ax Bu , y Cx Du может быть BIBO устойчивой . То есть, используя понятие BIBO
устойчивости, ничего нельзя сказать о внутренней устойчивости при не нулевых начальных
условиях и нулевом входе.
Управляемость линейных систем..
Интуитивно, понятие управляемости обозначает возможность перевода состояния системы
из одной точки пространства состояний в другую с помощью допустимого управления.
24
Событие (t0 , x0 ) линейной системы x A(t ) x B(t )u называется управляемым относительно
события (t1 , x1 ) , если существует момент времени t1 t0 и вход u (t ) , который определен на
конечном интервале [t0 , t1 ] и, обеспечивает перемещение системы из события (t0 , x0 ) в
событие (t1 , x1 ) .
Для линейных систем управляемость по состоянию, обычно, осуществляется по отношению к
точке x1 0 , так как если систему управляема по отношению к точке начала координат, то она будет
управляема и относительно любой точки x1 .
Теорема.
Линейная система x A(t ) x B(t )u является управляемой, если существует
управление u (t ) , которое позволяет переместить систему из точки x(t0 ) x0 в точку x(t1 ) x1 ,
когда вектор x0 (t0 , t1 ) x1 принадлежащую области линейного преобразования:
t1
W (t 0 , t1 ) (t 0 , t ) B(t ) BT (t )T (t0 , t )dt , где матричная функция (t , t0 ) 1 (t0 , t ) iявляется
t0
переходной функцией однородной линейной системы: x A(t ) . Кроме того, если x является
решением уравнения: W (t0 , t1 ) x x0 (t0 , t1 ) x1 , то тогда вход u(t ) BT (t )T (t0 , t1 ) x
является возможным управлением, которое позволяет переместить систему из точки
x(t0 ) x0 в точку x(t1 ) x1 за интервал времени [t0 , t1 ] .
Следствие.
Если в некоторые моменты времени t 0 и t1 , матрица W (t0 , t1 ) имеет ранг равный n , то тогда
система x A(t ) x B(t )u является управляемой.
Таким образом, управляемая система может быть преобразована в начало координат из
любой другой точки за любой малый интервал времени. Это будет возможно, если пара
матриц ( A, B) полностью известна и удовлетворяет условию W (t0 , t1 ) 0 . Кроме того, нет других
ограничений на сигнал управления. Однако, в этом случае, сигнал управления может быть очень
большим.
Рассмотрим теперь случай LTI модели системы x Ax Bu , x R n , u R m . Если LTI
модель системы является управляемой, то тогда говорят, что пара матриц ( A, B) является
управляемой. Для простоты, в дальнейшем примем, что t 0 0 .
Теорема.
Следующие утверждения являются эквивалентными:
1. Матричная пара ( A, B) является управляемой;
t
2. Матрица W (t ) e
A
T AT
BB e
t
d e A(t ) BB T e A
T
(t )
d является несингулярной для
любого t 0 ;
3. Матрица управляемости F ( B, AB,..., An1B) , размера (n nm) имеет ранг n (то есть
полный построчный ранг);
4. Если все собственные значения матрицы A имеют отрицательные вещественные части,
то тогда единственное решение WG уравнения: AWG WG AT BB T является
25
положительно определенным. Решение WG называется Грамианом управляемости и
T
задается соотношением: WG e A BB T e A d .
Следствие.
Пара матриц ( A, B), A R nn , B R nm будет управляемой тогда и только тогда, когда
матрица Fn r ( B, AB,..., An r B) имеет ранг равный n , где r rank (B) .
Пример.
Рассмотрим систему стабилизации платформы / 2 /.
0.5 0
0.5
, B . Начальное состояние
1
0
1
Система описывается уравнением is x Ax Bu , где A
платформы x1 (0) 10 , x2 (0) 1 . Найти управление, которое приведет платформу в состояние равновесия
за 2 секунды.
Ранг матрицы управляемости равен: rankF
0.5 0.25
2 , то есть система
rank ( B, AB ) rank
1
1
управляема. Вычислим Грамиан:
2
0.5
0.5
e 0.5
0
.
5
1
t 1
0
e
0.2162 0.3167
0
. Тогда управление имеет вид:
d
t
e
0.3167 0.4908
T
e 0.5( 2 t )
0 1 e 1 1 10
.
u (t ) BT e A (t1 t )W 1 (t1 )(e At1 x0 x1 ) 0.5 1
W (2)
2
e ( 2 t )
0 e 1
e
W (t 2)
0 0
Рассмотрим, изменяется ли управляемость LTI модели, если провести линейное
преобразование координат системы ~
x Tx , где T является несингулярной матрицей.
Теорема.
Управляемость LTI моделей является инвариантом по отношению ко всем невырожденным
эквивалентным преобразованиям.
Лемма.
SISO LTI модель, которая задается дифференциальным уравнением вида:
x (n) qn 1x ( n 1) ... q0 x p0u, p0 0 является управляемой.
Область пространства состояний, где все точки являются управляемыми, называется
областью управляемости LTI модели системы. Если линейная система является управляемой, то
это значит, что она полностью управляема во всем пространстве состояний. Если линейная
26
система является неуправляемой, то есть ранг матрицы управляемости rankF r , 0 r n , то
тогда говорят, что система является частично управляемой. Для частично управляемой системы
существует подпространство (область) управляемости, которое определяется базисом из
независимых столбцов матрицы F .
Функции оценки управляемости LTI моделей среды Matlab.
Функция ctrb(.) позволяет строить матрицу управляемости LTI модели. Синтаксис
команды ctrb(sys) позволяет формировать матрицу управляемости F на основе дескриптора
(указателя) системы sys . Синтаксис ctrb( A, B) формирует матрицу F на основе ссылки на пару
матриц ( A, B) .
Функция gram(.) вычисляет Грамиан управляемости путем решения матричного уравнения
Ляпунова.
Пример.
Наблюдаемость линейных систем.
Понятие наблюдаемости можно охарактеризовать, как возможность оценки состояния
системы по данным наблюдений за выходом. Проблема наблюдаемости системы отличается от
проблемы ее идентификации. Например, для LTI моделей, проблема идентификации сводится к
нахождению оценок матриц A, B, C, D на основе данных измерений входа и выхода системы. Задача
же наблюдения сводится к оценке текущего (или начального) состояния системы по измерениям
входа и выхода.
LTI модель состояния: x Ax Bu , y Cx Du называется наблюдаемой, если для любого
неизвестного начального состояния x(t0 ) , существует такой конечный момент времени t1 t 0 , что по
результатам измерений входа u (t ) и выхода y (t ) на интервале [t0 , t1 ] можно получить оценку
состояния x(t0 ) . В противном случае, модель называется ненаблюдаемой.
Выход модели: x Ax Bu , y Cx Du , который обусловлен влиянием начального
ненулевого состояния x(t0 ) и входным воздействием u (t ) , имеет вид:
27
t
y (t ) Ce A(t t 0 ) x(0) C e A(t ) Bu ( )d Du(t ) . Для простоты, будем в дальнейшем полагать, что:
t0
t0 0 . При этом предполагается, что вход u (t ) и выход y (t ) известны, а начальное состояние x(0) неизвестно. Введем обозначение функции ~
y (t ) Ce At x(0) , которую также можно определить с
t
помощью выражения: ~
y (t ) y (t ) C e A(t ) Bu ( )d Du(t ) .
Теорема.
Следующие утверждения эквивалентны:
1. Пара матриц ( A, C ) является наблюдаемой, где A R n n , C R p n ;
t
T
2. Матрица Wo (t ) e A t C T Ce At dt является несингулярной для любого t 0 ;
C
CA
3. Матрица наблюдаемости: H
, размера ( pn n) , имеет ранг равный n (полный
...
n 1
CA
ранг по столбцам).
4. Если все собственные значения матрицы A имеют отрицательные вещественные
части, то матричное уравнение Ляпунова AT Wo Wo A C T C имеет единственное
решение Wo , которое является положительно определенным. Решение Wo называется
T
Грамианом управляемости и имеет вид: Wo e A t C T Ce At dt .
Следствие.
Пара матриц ( A, C ) является наблюдаемой, тогда и только тогда, когда
матрица H n r (C T , AT C T ,..., ( AT ) n r C T )T имеет ранг равный n , где r rank (C ) .
Рассмотрим линейное эквивалентное преобразование ~
x Qx модели состояния, где Q несингулярная матрица.
Теорема.
Свойство наблюдаемости является инвариантным по отношению к невырожденному
линейному преобразованию.
Функции оценки наблюдаемость LTI систем в среде Matlab.
Функция obsv(.) позволяет вычислить матрицу наблюдаемости LTI модели системы. Синтаксис
команды obsv(sys) формирует матрицу наблюдаемости H на основе дескриптора (указателя) sys
системы. Команда obsv( A, C ) позволяет вычислить матрицу H по паре матриц ( A, C ) системы.
Функция gram(.) позволяет вычислить Грамиан наблюдаемости, путем решения соответствующего
матричного уравнения Ляпунова. Синтаксис команды имеет вид: gram(sys, ' o' ) .
28
Принцип дуальности (двойственности).
Построим следующую модель состояния вида: xˆ AT xˆ C T uˆ, yˆ BT xˆ DT uˆ , где:
xˆ R n , uˆ R p , yˆ R m . Такая модель системы называется двойственной по отношению к системе:
x Ax Bu , y Cx Du, x R n , u R m , y R p .
Теорема (принцип двойственности).
Пара матриц ( A, B) является управляемой, тогда и только тогда, когда пара матриц
( AT , BT ) является наблюдаемой и, соответственно, пара матриц ( A, C ) будет наблюдаемой,
тогда и только тогда, когда пара ( AT , C T ) является управляемой.
Декомпозиция линейных стационарных систем по Калману
Рассмотрим LTI модель: x Ax Bu , y Cx Du, x R n , u R m , y R p . Пусть задано
преобразование ~
x Px , где матрица P является несингулярной. Тогда, преобразованная система
~
~
~
~
~
~
~
~
имеет вид: ~
x A~
x Bu, y C~
x Du , где A PAP 1, B PB, C CP 1, D D .
Теорема (декомпозиция управляемости).
Пусть задана LTI модель: x Ax Bu , y Cx Du, x R n , u R m , y R p где матрица
управляемости имеет ранг: rank ( F ) rank ( B, AB,..., An 1B) r n . Тогда существует такая
матрица Q P 1 (q1, q2 ,..., qr , qr 1,..., qn ) , первые r столбцов которой состоят из r
независимых столбцов матрицы F , а оставшиеся (n r ) столбцов формируются произвольно,
таким образом, чтобы матрица P была несингулярной. При этом преобразование ~
x Px
позволяет записать преобразованную модель состояния системы в виде:
~
~
x ~ Ac
x ~ ~
c A ~ c B
u
~
xuc
xuc
0
~
~
~
xc Bc
A12 ~
~ ~ xc
u , y (Cc , Cuc ) ~ Du ,
~
xuc 0
Auc ~
xuc
~
~
~
где матрица управляемости F PF , Ac R r r , Auc R ( n r )( n r ) , а подмодель состояния:
~
~
~
~
x A ~
x B u, ~
y C ~
x Du, ~
x R r , u R m , y R p является управляемой и имеет
c
c c
c
c c
c
передаточную функцию, совпадающую с передаточной функцией исходной модели состояния
системы.
Пример.
Задана модель состояния x
0 1 0
0
1
1 2 0
1
0
Ax Bu , y Cx Du , где A 1 0 0 , B 0 , C T 0 , D 0 .
0 0 0
Вычислим матрицу управляемости F (b, Ab, A 2 b) 0 0 0 . Получим rank ( F ) 1 . Построим
1 0 0
0 1 1
матрицу P 0 1 0 , det(P) 1 , где первый столбец P является первым столбцом матрицы F , а другие
1 0 0
29
столбцы выбираются из условия, что det(P) 0 . Тогда преобразованная система может быть записана в
форме:
~
~
x A
~
~
~
x A~
x Bu или ~ c c
xuc 0
~
~
1
xc Bc
0 2 2 ~
A12 ~
~
1
,
где
,
u
A PAP 0 0 1 B PB 0 ,
~ ~
0 1 0
Auc xuc 0
0
~
~
~
~
C C P 1 0 0 1 , Ac 0 , Bc 1 , Auc 0 1 .
1 0
Теорема (декомпозиция наблюдаемости).
Пусть задана модель состояния: x Ax Bu , y Cx Du, x R n , u R m , y R p , матрица
C
CA
l n . Тогда существует такая
наблюдаемости которой имеет ранг: rank ( H ) rank
...
n 1
CA
p1
...
матрица P pl , где первые l строк состоят из l независимых транспонированных
...
pn
столбцов матрицы H , а оставшиеся строки матрицы P выбираются таким образом, чтобы
она была несингулярной. Тогда линейное преобразование ~
x Px позволяет записать
эквивалентную систему в форме:
~
~
Ao
~
xo
x ~
xo Bo
~
x ~ ~
0 ~
~
o A ~ c B
Du , где
y
(
C
,
)
,
u
u
~
~
~
o
~
~
A
~
xuo
xuc
xuo
21 Au 0 xuo Buo
~
~
Ao R l l , Auo R (n l )(n l ) , и подмодель состояния:
~
~
~
~
xo Ao ~
xo Bou, ~
y Co ~
xo Du, ~
xc R l , u R m , y R p является наблюдаемой и имеет
передаточную функцию, равную передаточной функции исходной системы.
Комбинируя декомпозицию управляемости и декомпозицию наблюдаемости можно перейти к
общей декомпозиции моделей состояния по Калману.
Теорема (Декомпозиция LTI моделей по Калману).
Каждая модель состояния может быть преобразована некоторым линейным
преобразованием к следующей форме:
~
~
xc 0 Aco
~
~
xcuo A21
~
xuco 0
~
xucuo 0
~
A13
~
A23
~
Auco
~
A43
~
~
xc 0 Bco
~ ~
xcuo Bcuo
~
~
xco
x Du , где вектор ~
u , y Cco 0 Cuco 0 ~
~
0 xuco 0
~
xucuo 0
Aucuo ~
x
является управляемым и наблюдаемым; вектор ~
является управляемым, но не наблюдаемым;
~
Acuo
~
A24
cuo
xuco является наблюдаемым, но не управляемым, и вектор ~
xucuo является не
вектор ~
управляемым и не наблюдаемым. Более того, исходная модель состояния эквивалентна, при
нулевых начальных условиях, управляемой и наблюдаемой ее редуцированной части:
~
~
~
~
xco Aco ~
xco Bcou, ~
y Cco ~
xco Du , и имеет передаточную функцию равную:
~
~
~
~
W (s) Cco (sE Aco ) 1 Bco D .
30
Результаты данной теоремы могут быть проиллюстрированы в графическом виде.
Из рисунка видно, что передаточная функция описывает только часть существующей модели
состояния системы.
Функции декомпозиции LTI моделей в Matlab.
Функция ctrbf (.) позволяет осуществлять преобразование моделей состояния, которые
задаются тройкой матриц ( A, B, C ) , в форму декомпозиции управляемости. Синтаксис команды
[ Ac, Bc, Cc,T , k ) ctrbf ( A, B, C) позволяет построить декомпозицию модели, заданную
~
A
0
0
~
~
матрицами Ac, Bc, Cc , представляемых в виде: Ac ~uc ~ , Bc ~ , Cc Cuc Cc .
A
Bc
21 Ac
Матрица преобразования T определяется как ортогональная несингулярная матрица.
Вектор k показывает, какое число итераций было проведено алгоритмов декомпозиции на каждом
шаге. Величина sum(k ) указывает на число независимых столбцов в матрице управляемости, а
также определяет размерность управляемой части модели состояния.
Пример.
Функция obsvf (.) позволяет преобразовать модель состояния, заданную тройкой
матриц ( A, B, C ) , в форму декомпозиции наблюдаемости. Синтаксис
команды [ Ao, Bo, Co,T , k ) ctrbf ( A, B, C ) позволяет вычислять матрицы Ao, Bo, Co , формы
31
~
~
~
Auo A12
Buo
~
, Bo ~ , Co 0 Co . Матрица
декомпозиции наблюдаемости в виде: Ao
~
0
Ao
Bo
преобразования T ищется как ортогональная и несингулярная. Вектор k определяет число итераций
алгоритма на каждом шаге. Величина sum(k ) определяет число независимых строк в матрице
наблюдаемости, а также размерность наблюдаемой части преобразованной модели.
Пример.
Минимальная реализация SISO систем. Понятие грубость линейных систем.
Нам уже известно, что система, заданная передаточной функцией, может иметь несколько
моделей состояния, заданных разными уравнениями одной размерности. Одновременно,
передаточная функция является в ряде случаев избыточным описанием, так имеет сокращаемые
нули и полюса. Это порождает проблему поиска такой реализации передаточной функции, чтобы
получить минимальную размерность уравнения состояния среди всех возможных.
В данном разделе мы рассмотрим задачу поиска минимальной реализации для SISO LTI систем,
которая тесно связана с проблемой грубости линейных систем.
Пусть система с одним входом и одним выходом задана правильной передаточной функцией
вида:
P( s) pm s m pm 1s m 1 ... p0
,m n.
Q( s )
s n qn 1s n 1 ... q0
Будем называть передаточную функцию W (s) несократимой, нет нулей и полюсов, которые можно
W ( s)
взаимно сократить, или, что эквивалентно, числитель и знаменатель не имеют общих корней или
множителей.
Рассмотрим теперь модель состояния SISO системы: x Ax Bu , y Cx , где x R n , u, y R . Эта
система имеет передаточную функцию W ( s) C ( sEn n A) 1 B . Предположим теперь, что один или
более нулей и полюсов можно взаимно сократить, то есть сократить общие множители в числителе
~
P ( s)
~
~
~ deg ree( P
(s)) m ,
P(s) и знаменателе Q(s) . Тогда можно записать: W ( s) ~ , где m
Q( s)
~
n~ deg ree(Q(s)) n .
~
Ясно, что модель состояния, построенная по передаточной функции W ( s) , имеет
~ , в то время, как исходная система имеет порядок n n~ , то есть не является
размерность n
~ ~ ~
минимальной. Обозначим через ( A, B , C ) матрицы модели состояния, построенной по
~
~ ~
~
~
передаточной функции W ( s) : W ( s) C ( sEn~ n~ A) 1 B . Возникает вопрос, является ли такая
реализация системы минимальной?
32
Теорема.
Пусть задана передаточная функция W ( s)
pm s m pm 1s m 1 ... p0
где n deg ree(Q(s)) s n qn 1s n 1 ... q0
степень знаменателя передаточной функции. Тогда n - мерная модель состояния, реализованная по
передаточной функции W (s) , будет минимальной тогда и только тогда, когда передаточная функция
W (s) является несократимой.
Пример.
Задана модель состояния системы x Ax Bu ,
0
0
1 ;
6 11 6
y Cx , где x R 3 , u, y R , A 0
1
1
0
6( s 1)
; T . Отсюда найдем передаточную функцию
. Очевидно, что
W ( s)
B 6 C 0
(
s
1
)(
s
2
)(
s
3
)
0
30
6
. Тогда минимальной реализацией
( s 2)(s 3)
1 ~ 0
~
~ 0
~
~
2
, B ,
системы будет модель состояния: x Ax Bu , y Cx , где x R , u, y R , A
6 5
6
несократимая передаточная функция равна
~
W ( s)
1
C T .
0
Пример.
4 1.25 1.5
2
0 , B 0 ,
Пусть задана система x Ax Bu , y Cx , где x R , u, y R , A 4
0
0
1
0
0
6( s 1)
T
C 0.75 . Отсюда найдем передаточную функцию W ( s)
, то есть
( s 1)(s 2)(s 3)
0.75
3
~
W ( s)
~
6
~
~
2
. Тогда минимальная реализация имеет вид: is: x Ax Bu , y Cx , where x R ,
( s 2)(s 3)
1 ~ 0 T 1
~ 0
, B , C .
u, y R , A
6 5
6
0
Таким образом, из анализа примеров, мы видим, что две различные системы имеют одну и ту же
минимальную реализацию. Рассмотрим теперь поведение исходных систем и их минимальной
реализации.
Все траектории системы первого примера и ее минимальной реализации стремятся к точке
равновесия, которая для обеих систем имеет один и тот же тип: устойчивый узел. Следовательно,
минимальная реализация имеет такое же поведение, что и исходная система.
Во втором примере часть траекторий исходной системы стремится к точке равновесия, а другая
часть траекторий удаляется от точки равновесия. То есть точка равновесия имеет тип: седло. У
системы минимальной реализации все траектории стремятся к точке равновесия типа устойчивый
узел. То есть поведение исходной системы и ее минимальной реализации различно.
Следовательно, если для дальнейшего проектирования выбрать минимальную реализацию, то при
33
проектировании можно совершить ошибку. Малейшее различие в параметрах исходной системы
приведет к невозможности сокращения нулей и полюсов, лежащих в правой части комплексной
плоскости. Это приведет к тому, что система станет неустойчивой.
Отсюда следует, что при построении несократимой передаточной функции можно взаимно
сокращать только нули и полюса, лежащие в левой части комплексной плоскости.
Понятие грубости систем.
Интуитивно понятно, что грубость систем, как математических моделей, заключается в том, что
небольшие изменения их параметров должны приводить к небольшим изменениям в их поведении.
Впервые понятие грубости систем было введено А. Андроновым и Л. Понтрягиным.
Физическое обоснование грубости систем заключается в следующем: если математическая
модель известна приблизительно, а такое имеет место в большинстве практических случаев, то
дискуссия о влиянии на поведение системы малых изменений ее параметров является
бессмысленной.
Математический аспект проблемы грубости заключается в том, что малые изменения параметров
системы не должны изменять топологический портрет ее траекторий, то измененная система и
исходная должны быть эквивалентными, за исключением особых критических случаев, требующих
дополнительного исследования (так называемые точки бифуркации параметров).
Основной вывод из приведенных рассуждений является следующим: проектируемая система
должна быть грубой, то есть малые изменения параметров не должны изменять ее качественное
поведение и переводить в другой топологический тип, за исключением критических случаев,
связанных с потерей устойчивости.
Управляемость, наблюдаемость и минимальная реализация MIMO систем.
Рассмотрим SISO LTI модель системы, которая задается строго правильной передаточной
функцией: W ( s)
D( s )
, deg ree( D( s)) deg ree( A( s)) . Предположим , что нули k , k 1,2,..., n
A( s)
полинома A(s) являются простыми и вещественными. То есть можно записать
соотношение: W ( s)
n
D (k )
D ( k )
. Очевидно, что величины bk
являются константами.
A(k )
k )( s k )
A(
k 1
Построим по передаточной функции модель состояния, используя преобразование Лурье:
n
xk k xk bk u, k 1,2,..., n , y xk . Тогда модель состояния примет вид: x Ax Bu , y Cx , где
k 1
b1
1 ... 0
B ... , A 0 ... 0 , C 1 ... 1 . Вычислим матрицу управляемости такой модели:
b
0 ...
n
n
b1
b
F ( B, AB..., An 1B) 2
.
b
n
1 1
n
1 2
det(F ) bk det
k 1
. .
1
n
... b11n 1
b2 2 ... b2 n2 1
. Отсюда найдем:
.
...
.
bn n ... bn nn 1
... 1n 1
... n2 1
.
...
.
... nn 1
b11
34
Второй множитель, в выражении определителя матрицы управляемости, является
определителем Вандермонда, который не равен нулю. Тогда определитель матрицы F будет равен
нулю, когда существует хотя бы один коэффициент bk 0 . НО это будет возможно, когда нуль
числителя D(s) равен s k . То есть, это возможно, когда передаточная функция имеет
сокращаемые нуль и полюс. В этом случае, координата xk модели состояния будет не зависеть от
входа u . Таким образом, передаточная функция, имеющая сокращаемые нули и полюса, не может
описывать полностью модель состояния.
Рассмотрим LTI модель x Ax Bu , y Cx Du, x R n , u R m , y R p . Примем, для
простоты, D 0 . Пусть передаточная функция системы W (s) C ( sE A) 1 B является
несокращаемой, то есть полиномы числителя и знаменателя не содержат общих множителей / 5, 6/.
Тогда можно сформулировать следующие условия полной наблюдаемости и управляемости.
Теорема.
Модель состояния x Ax Bu , y Cx, x R n , u R m , y R p является полностью
управляемой и наблюдаемой, тогда и только тогда, когда передаточная
функция W (s) C ( sE A) 1 B является несокращаемой.
Если начальное состояние модели состояния является нулевым, то выход y (t )
системы ( A, B, C, D) может быть получен, на основании декомпозиции по Калману, в виде:
t
~
~
~
y (t ) Cco e Aco (t ) Bcou ( )d Du(t ) . Данное выражение совпадает с выходом полностью
t0
управляемой и наблюдаемой части системы. Для простоты примем, что D 0 . Тогда
~
~
~
подсистема ( Aco , Bco , Cco ) будет совпадать с минимальной реализацией исходной системы.
Минимальная реализация модели всегда является полностью управляемой и наблюдаемой, и имеет
несокращаемую передаточную функцию.
Алгоритм получения минимальной реализации модели состояния может быть следующим (здесь
формы декомпозиции управляемости и наблюдаемости приняты совпадающими с
соответствующими формами среды Matlab):
1. Найти матрицу преобразования Tc1 , позволяющую выделить управляемую и
~
Auc
~
1
неуправляемую части (подсистемы): A Tc ATc ~
A
21
~
~
~
C CTc Cuc Cc .
0
0 ~
1
~ , B Tc B ~ ,
Ac
Bc
~ ~ ~
2. Найти матрицу преобразования To такую, что в управляемой подсистеме ( Ac , Bc , Cc )
~
~
Acuo Ac,12
~
1 ~
позволяет выделить наблюдаемую часть (подсистему): Ao To AcTo
~ ,
0
Aco
~
Bcuo ~
~
~
~
1 ~
Bo To Bc ~ , Co CcTo 0 Cco .
Bco
0
En rank( A~ )
c
. Тогда матрица T 1 Tˆo1Tc1 будет определять
3. Построить матрицу: Tˆo1
1
To
~ ~ ~
преобразование исходной системы в минимальную ее реализацию в виде ( Aco , Bco , Cco ) .
35
С помощью линейного преобразования с матрицей T исходная система преобразуется в форму
~
Auc
~
декомпозиции по Калману: z A21
0
~
Acuo
0 0
~
~ ~
Ac,12 z Bcuo u , y Cuc
~
~
Aco Bco
~
0 Cco z Du .
Пример.
Задана модель состояния x
5
4
A
0
0
0
,
0 0 4
0 2 6
8
7
Ax Bu , y Cx Du, x R 4 , u R, y R , где:
4
2
B , C 2 2 2 1 , D 0 .
2
1
Выделить управляемую и наблюдаемую подсистему (минимальную реализацию) .
Ранг матрицы управляемости F
( B, AB, A2 B, A3 B) равен 3. Ранг матрицы наблюдаемости
H (C T , AT C T , ( AT ) 2 C T , ( AT )3 C T )T равен 3. Передаточная функция исходной системы в zpk 10( s 2.6)(s 3)(s 4)
представлении имеет вид: W ( s)
. Таким образом, порядок минимальной
( s 1)(s 2)(s 3)(s 4)
10( s 2.6)
реализации системы равен 2. Передаточная функция минимальной реализации равна Wmr ( s)
.
( s 1)(s 2)
Файл E41_CalmDcm.
36
Каноническая форма наблюдаемости.
Каноническая форма наблюдаемости – минимальная реализация (метод 1)
37
Каноническая форма наблюдаемости – минимальная реализация (метод 2).
В среде Matlab имеется функция min real (.) , которая позволяет формировать минимальную
реализацию LTI моделей систем. С помощью этой функции из системы удаляются неуправляемая и
ненаблюдаемая части.
Пример.
x Ax Bu , y Cx Du, x R 4 , u R, y R , где:
4
2
B , C 2 2 2 1 , D 0 .
2
1
Задана LTI модель системы:
5
4
A
0
0
0
,
0 0 4
0 2 6
8
7
Построить наблюдаемую и управляемую подсистему (минимальную реализацию).
Файл E42_MinReal.
38
Анализ линейных нестационарных систем.
Асимптотическая устойчивость линейной стационарной системы полностью определяется
корнями характеристического многочлена ее коэффициентов. Поэтому, естественно, попытаться
связать асимптотическую устойчивость неавтономной системы x A(t ) x с корнями j (t )
уравнения:
det[A(t ) E ] 0 .
В первую очередь желательно ответить на вопрос: гарантирует ли неравенство:
max Re( j (t )) ,
j
где
0 - положительная постоянная для любого t [t0 , ) ,
асимптотическую устойчивость системы? Ответ на этот вопрос отрицательный.
Пример.
1 2 cos 4t 2 2 sin 4t
x . Собственными значениями ее
x
2 2 sin 4t 1 2 cos 4t
характеристического многочлена являются функции 1 (t ) 2 (t ) 1 . Однако эта система обладает
Пусть задана система:
неограниченным решением
1 (t ) e t sin 2t , 2 (t ) e t cos 2t .
39
Таким образом, в общем случае, неравенство
max Re( j (t ))
для корней
j
характеристического многочлена не гарантирует, в общем случае, асимптотическую
устойчивость неавтономной системы.
Переходная матрица автономной линейной системы вида x Ax , имеет следующий вид:
X (t ) exp( At ) . При этом, элементы фундаментальной матрицы, нормированной в точке t 0 ,
представляют собой линейные комбинации функций p j (t ) exp( j t ) , где p j (t ) - некоторые
многочлены,
j - собственные числа матрицы A . В терминах этих чисел (собственных значений),
как известно, можно полностью характеризовать свойство асимптотической устойчивости системы
x Ax . Оказывается, как показал А.М. Ляпунов, со всякой неавтономной линейной
дифференциальной системой можно связать конечный набор действительных чисел, которые, как и
в автономном случае, позволяют исследовать ее асимптотическую устойчивость.
Пусть f (t ) - числовая функция, заданная и непрерывная на полупрямой R . Число ( или
символы и ) определяемое формулой:
____
1
t t
[ f ] lim ln | f (t ) | , называется
характеристическим показателем или показателем Ляпунова функции f (t ) . Здесь под
____
обозначением lim (t ) понимается верхний предел функции
t
(t ) при t .
Рассмотрим теперь однородную неавтономную систему x A(t ) x с непрерывной на R
матрицей A(t ) R
n*n
.
Теорема Ляпунова (о характеристических показателях решений линейной неавтономной
системы) . Если матрица A(t ) ограничена, то есть || A(t ) || c const , то каждое
ненулевое решение системы x A(t ) x имеет конечный характеристический показатель.
Определение. Множество всех характеристических показателей ненулевых решений системы
x A(t ) x называется ее спектром.
Теорема. Спектр линейной системы x A(t ) x с ограниченной матрицей A(t ) состоит из не
более, чем n различных элементов.
Используя понятие спектра, можно достаточно детально описать структуру общего решения
уравнения x A(t ) x . Пусть совокупность функций x1 (t ), x2 (t ),..., xn (t ) образует фундаментальную
систему, а числа
1 , 2 ,..., n -полный спектр рассматриваемого уравнения, то есть
[ x j ] j , j 1,2,..., n . Положим j e
j t
x j (t ) . Отсюда получим [ j ] j [ x j (t )] 0 .
Отсюда, общее решение уравнения x A(t ) x может быть представлено в виде
n
c j j (t )e
jt
, где
j 1
c j - произвольные постоянные, а функции j имеют нулевые характеристические показатели.
Отсюда сразу же следует, что если полный спектр системы отрицательный, то есть
такая система асимптотически устойчива.
j 0 , то
40
Теорема. Неравенство Важевского. Для любого решения системы x A(t ) x , где A(t ) C ,
t
t
t0
t0
при t 0 t справедливо неравенство: || x(t 0 ) || exp{
где
( )d } || x(t ) |||| x(t 0 ) || exp{ ( )d } ,
(t ) и (t ) наименьший и наибольший характеристический корни эрмитово-
симметризованной матрицы: A (t )
H
1
[ A(t ) A* (t )] .
2
Следствие 1. Для асимптотической устойчивости системы x A(t ) x достаточно выполнения
условия (t ) h 0 , где h - положительная постоянная.
Нетрудно убедиться, что спектр системы x A(t ) x отрицательный, если квадратичная
T
H
форма x A (t ) x определенно отрицательна, то есть справедливо условие
x T A H (t ) x | x | 2 , где const 0 , а A H (t ) в случае действительной матрицы A(t ) имеет
вид A (t ) [ A(t ) A (t )] . Данное условие имеет место, когда главные миноры k (t ) матрицы
H
T
A H (t ) отделены от нуля, то есть det{ k (t )} const 0 и чередуют знаки. Причем первый из
определителей меньше нуля.
Критерий неустойчивости неавтономных систем. Пусть задана система x A(t ) x . Тогда
1t
lim tr ( A( )) d 0
t t
t0
____
из условия
следует неустойчивость ненулевых решений системы.
Методы и алгоритмы синтеза линейных стационарных систем
Частотные методы синтеза систем.
Можно сформулировать следующую схему решения задачи синтеза в частотном диапазоне.
1.
Построить ЛАЧХ Lnc ( ) и ЛФЧХ nc ( ) (диаграмму Боде) некомпенсированной разомкнутой
системы.
2.
На основе требований к точности и качеству замкнутой системы построить ЛАЧХ Ld ( ) и
ЛФЧХ d ( ) желаемой компенсированной разомкнутой системы.
Получит, методом графического вычитания, частотные характеристики компенсатора:
Lmc ( ) Ld ( ) Lnc ( ) , mc ( ) d ( ) nc ( ) .
3.
4.
Построить, на основе анализа частотных характеристик компенсатора: Lmc ( ), mc ( ) ,
передаточную функцию компенсатора.
5.
Вычислить передаточную функцию замкнутой системы при включении в ее состав
спроектированного компенсатора. Оценить устойчивость, точность и качество замкнутой системы.
6.
Сравнить полученные оценки с заданными требованиями к замкнутой системы. Если
требования выполняются, то приступить к проектированию физически реализуемого компенсатора.
Если требования не выполняются, то вернутся к шагу 2 алгоритма и откорректировать частотные
характеристик желаемой разомкнутой системы.
Как видно, из приведенных выше рассуждений и примеров, основной проблемой синтеза систем в
частотном диапазоне является построение желаемых частотных характеристик компенсированной
разомкнутой системы, которая будет удовлетворять заданным требованиям по точности и качеству.
41
Укажем несколько практических рекомендаций, которые облегчают процесс построения желаемых
частотных характеристик.
Прежде всего, необходимо определить частотный диапазон работы замкнутой системы путем
анализа частотных характеристик входа, выхода и шумов. Затем выделить низко – частотный,
средне – частотный и высоко - частотный диапазоны.
Низко – частотный диапазон определяется минимальной угловой частотой l наиболее
медленного звена, входящего в состав некомпенсированной разомкнутой системы. ЛАЧХ
разомкнутой системы имеет наклон 0 dB или кратный 20v dB/decade, где параметр v определяет
порядок астатизма системы. Амплитуда ЛАЧХ на частоте 1 определяет желаемый коэффициент
усиления системы.
Средне – частотный диапазон определяется частотами в интервале c1 c 2 , который
содержит частоту среза c ; c [c1, c 2 ] желаемой ЛАЧХ разомкнутой системы. В
соответствии с критерием Найквиста частота среза должна удовлетворять соотношению c b ,
где частота b определяется условием d (b ) 1800 . Обычно предполагается, что в окрестности
частоты: c желаемые частотные характеристики должны совпадать с соответствующими
частотными характеристиками устойчивых звеньев первого или второго порядка, то есть наклон
желаемой ЛАЧХ должен составлять -20 dB/decade или -40dB/decade. Величина c выбирается из
условия обеспечения требуемого времени t s затухания переходного процесса и максимального
P
перерегулирования . Для этого можно воспользоваться таблицей значений функции t sc ( max ) ,
P0
где Pmax , P0 - максимальное и начальное значения амплитудно – частотной характеристики Pc ( )
замкнутой системы.
P
Используя заданное значение перерегулирования r , можно по таблице найти величину max , а
P0
затем найти произведение t sc . Зная требуемое значение t s , отсюда легко найти частоту среза c
желаемой разомкнутой системы. Затем, по оценкам запасов устойчивости по амплитуде и фазе,
можно определить коэффициенты усиления: Ld (c1 ), Ld (c 2 ) . Частоты c1 , c 2 выбираются таким
образом, чтобы c1 10
c 2
10 , 2 c 4 .
c1
c1
Ld ( c1 )
20
c и c 2 10
Ld ( c 2 )
20
c , а также выполнялись неравенства:
42
Участок стыковки ЛАЧХ низко – частотного и средне – частотного диапазонов (то есть,
когда l c1 ) выполняется в виде прямых линий с наклоном -40 dB/decade или -60 dB/decade.
Высоко – частотный диапазон желаемой ЛАЧХ (то есть при c 2 ) выполняется с наклоном
превышающим -40 dB/decade.
Методы компенсации и желаемой передаточной функции.
Метод желаемого размещения полюсов (root – locus method) является одним из основных
способов улучшения динамики замкнутой системы во временной области, позволяющий достигать
требуемых прямых показателей качества. При этом используются различные схемы включения
компенсаторов (контроллеров) в систему.
На рисунках a), b) показаны последовательные схемы включения компенсаторов. Параллельная
схема включения компенсатора приведена на рисунке c).
Пример.
Рассмотрим структурную схему, приведенную на рисунке a). Пусть передаточная функция объекта
управления равна W p ( s)
2
s 3s 1
2
. Спроектируем контроллер, который будет гарантировать
время затухания переходного процесса не более t s 3 секунд, перерегулирование 0% и ошибку
в установившемся режиме e() 0 .
ki
s 2 as b
Для этого будем использовать PID контроллер (регулятор): Wcm ( s) k p k d s k d
,
s
s
kp
k
, b i . Выберем в качестве желаемой передаточной функции, удовлетворяющей
где a
kd
kd
1
. Замкнутая передаточная функция системы с PID
s 1
Wcm ( s)W p ( s)
указанным выше требованиям: Wd ( s)
контроллером равна: Wc ( s)
1 Wcm ( s)W p ( s)
. Выберем в качестве компенсатора динамики объекта
контроллер с передаточной функцией вида: s 2 as b s 2 3s 1 , i.e. a
kp
kd
3, b
ki
1 . Тогда
kd
43
передаточная функция замкнутой системы примет вид: Wc ( s)
2k d
1
s 2k d
1
2k d
s 1
Wd ( s) . Отсюда
найдем, что: k d 0.5 . Следовательно, получим: get k p 1.5, ki 1 и, окончательно, получим, что
передаточная функция контроллера равна: Wcm ( s) 1.5 0.5s
1
.
s
Рассмотрим случай, когда желаемая передаточная функция последовательного
контроллера Wd (s) замкнутой системы известна и удовлетворяет заданным требованиям. Тогда
передаточную функцию компенсатора (контроллера) Wcm (s) можно найти из
уравнения:
Wcm ( s)
Wcm ( s)W p ( s)
1 Wcm ( s)W p ( s)
Wd ( s) . Решение Wcm (s) данного уравнения имеет вид:
Wd ( s)
1
m( s )
. Перепишем передаточную функцию объекта в виде: W p ( s)
. Тогда
W p ( s) (1 Wd ( s))
p( s)
можно записать соотношение: Wcm ( s)
p( s) Wd ( s)
. Чтобы такая функция была физически
m( s) (1 Wd ( s))
w(t ), t 0
, то есть:
0, t 0
реализована, она должна удовлетворять соотношениям: L1{Wcm ( s)}
Wcm ( s)
rl s l ... r0
s k qk 1s k 1 ... q0
, l k . При этом, замкнутая система должна быть грубой, то есть при
компенсации не происходило сокращение положительных полюсов и нулей на s – плоскости.
Детализируем передаточную функцию объекта в виде: W p ( s)
m( s) m ( s)m ( s)
, где p ( s), m (s)
p( s) p ( s) p ( s)
полиномы, имеющие нули с отрицательной вещественной частью, а p ( s), m (s) - полиномы,
имеющие нули с положительной вещественной частью. Такое представление передаточной
функции называется ее факторизацией / 2 /. Тогда передаточную функцию компенсатора можно
записать в виде: Wcm ( s)
p ( s) p ( s)
Wd ( s)
. Отсюда, чтобы нули и полюса с положительной
m ( s)m ( s) (1 Wd ( s))
вещественной частью взаимно не сокращались, желаемая передаточная функция должна
m (s) K (s)
p ( s ) D( s )
удовлетворять условиям: Wd ( s)
, 1 Wd ( s)
, где полином N (s) - является
N ( s)
N (s)
знаменателем передаточной функции Wd (s) , а K (s), D(s) - некоторые полиномы, которые должны
удовлетворять уравнению: m (s) K (s) p (s) D(s) N (s) . Данное уравнение будет разрешимо
относительно K (s), D(s) , когда степень полинома N (s) будет меньше,
чем deg ree( K (s)) deg ree( D(s)) 2 , то есть:
deg ree( N (s)) deg ree( K (s)) deg ree( D(s)) 1 (a).
Если передаточная функция объекта имеет, в своем составе, v p интегральных звеньев, а
передаточная функция замкнутой системы должна иметь астатизм v - го порядка (тип системы
будет равен v ), то передаточная функция компенсатора будет содержать множитель s
v v p
, то есть:
44
1 Wd ( s)
p ( s ) D( s ) s
N ( s)
v v p
, или: Wcm ( s)
примет вид: m (s) K (s) p (s) D(s)s
компенсатора можно записать, как:
v v p
p (s)
K ( s)
m ( s ) D( s ) s v v p
. Тогда полиномиальное уравнение
N (s) . Тогда условие физической реализуемости
deg ree( p (s)) deg ree( K (s)) deg ree(m (s)) deg ree( D(s)) v v p (b).
Так как степень числителя передаточной функции Wd (s) не может превышать степень
знаменателя, то тогда функция (1 Wd ( s)) должна иметь числитель и знаменатель одного порядка.
То есть получим следующее условие:
deg ree( N (s)) deg ree( p (s)) deg ree( D(s)) v v p (c ).
Отсюда, из анализа условий (a ), (b ), (c ) получим неравенство / 2,3 /:
deg ree( N (s)) deg ree( p(s)) deg ree( p (s)) deg ree(m (s)) v v p 1 .
Таким образом, алгоритм определения передаточной функции компенсатора можно записать в
следующем виде:
1.
Определить желаемую передаточную функцию Wd ( s)
m (s) K (s)
замкнутой системы со
N (s)
знаменателем N (s) .
m ( s )m ( s )
2.
Провести факторизацию передаточной функции объекта: W p ( s)
3.
Записать условия (a ), (b ), (c ) и выбрать соответствующие минимальные степени nK , nD
p (s) p (s)
.
полиномов K (s), D(s) .
4.
Записать выражения полиномов K (s), D(s) с неопределенными коэффициентами и построить
полиномиальное уравнение: m (s) K (s) p (s) D(s) N (s) .
5.
Решить уравнение m (s) K (s) p (s) D(s) N (s) путем нахождения соответствующих
коэффициентов K (s), D(s) .
6.
Задать передаточную функцию компенсатора в виде выражения: Wcm ( s)
p (s)
K ( s)
m ( s ) D( s ) s v v p
.
Пример.
Пусть передаточная функция объекта управления равна: W p ( s)
1
. Провести синтез
( s 1)( s 1)
компенсатора, который будет обеспечивать замкнутой системе астатизм первого порядка и
характеристический полином (знаменатель) равный: N ( s) ( s 1) 3 .
Таким образом, получим, что: v 1 . Факторизованная передаточная функция объекта
имеет: m ( s) m ( s) 1 , p ( s) s 1, p (s) s 1 . Условия физической реализации
последовательного компенсатора имеют вид: deg ree( N (s)) 3 , deg ree(m (s)) deg ree(m (s)) 0 ,
deg ree( p (s)) deg ree( p (s)) 1 . Тогда условия (a ), (b), (c ) можно записать в виде: (a )
3 nK nD 1 , (b ) 1 nK nD 1 , (c ) 3 1 nD 1 . Отсюда найдем:
nD deg ree( D(s)) 1; nK deg ree( K (s)) 1 . Следовательно: D(s) d1s d 0 ; K (s) k1s k0 .
45
Полиномиальное уравнение тогда имеет вид: k1s k0 ( s 1)(d1s d 0 ) s s 3 3s 2 3s 1 . Отсюда
найдем: d1 1, d 0 4 , k1 7, k0 1 . Следовательно передаточная функция компенсатора может
быть записана в форме: Wcm ( s)
p (s)
K (s)
m ( s ) D( s ) s v v p
( s 1)(7 s 1)
, а желаемую передаточную
( s 4) s
m ( s) K ( s) 7 s 1
функцию можно записать, как: Wd ( s)
.
N ( s)
( s 1) 3
Как мы уже обсуждали, имеется тесная связь между прямыми показателями качества замкнутой
системы и расположением полюсов и нулей ее передаточной функции, которые определяются
видом желаемой передаточной функции. Ограничения на прямые показатели качества обычно
характеризуют не один, какой – то, конкретный переходный процесс, а целый класс временных
ответов системы на стандартный вход / 4 /. Возникает вопрос, как в этом случае определить такой
класс желаемых передаточных функций.
Теорема.
Пусть SISO LTI модель системы задана нормализованной передаточной функцией вида:
Wn ( s)
(bm s m ... b1s b0 ) ( p1 p2 ... pn )
, m n , где полюса pi являются вещественными и
n
b0
(s pi )
i 1
представлены в упорядоченном виде: 0 p1 p2 ... pn . Обозначим переходную функцию такой
системы Wn (s) , как xnm (t , p1 , p2 ,..., pn ) . Соответственно, обозначим переходную функцию
системы Wdn ( s)
(bm s m ... b1s b0 ) ( p n )
m
, как xdn
(t , p) . Тогда существует такое вещественной
n
b
( s p)
m
(t , p(t )) для любого t [0, ] ,
число p 0 , что будет выполняться равенство: xnm (t , p1, p2 ,..., pn ) xdn
при условии, что: p1 p pn .
Пример.
( s 2 3s 3) ( p 3 )
( s 2 0.5s 3) ( p 3 )
W
(
s
)
,
и
.
dn
3
3
( s p) 3
( s p) 3
Построить области переходных процессов для указанных передаточных функций при p [0.4,4] , и
p [1,6] , соответственно. Оценить область нахождения переходных функции систем относительно
Заданы две передаточные функции: Wdn ( s)
указанных областей с передаточными функциями:
46
Wdn ( s)
Bi ( s)
( p1 p 2 P3)
, p1, p 2, p3 0 , где B1 (s) s 2 3s 3 и
( s p1)(s p 2)( s p3)
3
B2 (s) s 2 0.5s 3 , соответственно.
Файл E53_DsMult.
p1 0.5; p2 1; p3 3;
p1 1.5; p2 3; p3 5;
Существует и более общий результат. Пусть задана SISO LTI модель системы, заданная
bm s m ... b1s b0 a0
, m n , которая имеет
передаточной функцией: Wn ( s) n
s an 1s n 1 ... a0 b0
47
полюса: pi , i 1,2,..., n с отрицательной вещественной частью, такими, что: Re( pi ) [ pmin , pmax ] , где:
pmin min Re( pi )) , pmax max Re( pi )) , и переходной функцией xnm (t , p1 , p2 ,..., pn ) . Обозначим
i
i
m
через: xdn
(t , p) переходную функцию системы с передаточной функцией с кратными полюсами:
bm s m ... b0 | p |2q
, n 2q
( s p) q ( s p ) q b0
, такими, что: pr [ pmin , pmax ] ,
Wdn ( s)
bm s m ... b0
pr | p |2 q
, n 2q 1
( s p )( s p) q ( s p ) q
b
r
Re( p) [ pmin , pmax ] , Im( p) [ f max , f max ] , где f max max Im( pi )) . Тогда существуют вещественное
i
число pr 0 , и комплексное число p, Re( p) 0 такие, что будет выполняться равенство:
m
xnm (t , p1, p2 ,..., pn ) xdn
(t , pr , p) для любого t [0, ] , и pr [ pmin , pmax ] , Re( p) [ pmin , pmax ] ,
Im( p) [ f max , f max ] .
Если определить связь между величинами: Re( p) , и Im( p) как:
определяет сектор с корневым углом: d , d
2
Re( p)
tan( d ) , где угол d
Im( p)
, то тогда можно задать только два параметра:
Re( p) [ pmin , pmax ] , pr [ pmin , pmax ] , для случая, когда: n 2q 1 , или один параметр вида:
Re( p) [ pmin , pmax ] , когда n 2q . Здесь угол d определяет колебательность системы.
Пример.
Задана система с передаточной функцией с кратными полюсами вида:
s 2 3s 3
r | p |2
Wdn ( s)
. Временная область ее переходных процессов на
( s r )(s 2 2 ps p 2 f 2 ) 2 3
входной единичный скачок (функцию Хэвисайда) определяется, при вариациях
параметров: r [4;0.2] ; p [4;0.2] и f [0;2] . Оценить, будет ли попадать в данную временную
область переходная функция системы с передаточной функцией:
Wdn ( s)
s 2 3s 3
7.5
.
( s 1.5)(s 2s 2)(s s 2.5) 3
2
Файл E54_DsCMult
2
48
Синтез систем методом обратной связи по состоянию.
Пусть задана SISO LTI система: x Ax Bu; y Cx; x R n ; u R, y R , и управление имеет
вид: u(t ) g Kx , где g R задающее воздействие, а вектор K имеет размер (1 n) . Тогда
замкнутая система будет описываться уравнением: x ( A BK ) x Bg ; y Cx / 3 /.
Теорема.
Пусть SISO LTI модель системы x Ax Bu; y Cx; x R n ; u R, y R является полностью
управляемой и, задан желаемый полином вида d ( ) n n 1n 1 ... 0 с вещественными
коэффициентами. Тогда существует вектор K R1 n обратной связи по состоянию такой, что
49
модель состояния x ( A BK ) x; y Cx замкнутой системы имеет характеристический
полином, совпадающий с желаемым полиномом, то есть det(E A BK ) d ( ) .
Лемма.
Пусть SISO LTI модель системы: x Ax Bu; y Cx; x R n ; u R, y R является полностью
управляемой и, задан желаемый полином вида: d ( ) n n 1n 1 ... 0 с вещественными
коэффициентами. Тогда существует вектор K R1 n обратной связи по состоянию такой, что
передаточная функция замкнутой системы Wc ( s) C ( sE A B K ) 1 B имеет знаменатель,
совпадающий с желаемым полиномом d (s) , а числитель совпадает с числителем передаточной
функции Wo ( s) C ( sE A ) 1 B разомкнутой системы.
Таким образом, использование обратной связи по состоянию позволяет управлять расположением
полюсов, но не нулей, передаточной функции замкнутой системы. Алгоритм построения обратной
связи по состоянию можно сформулировать в следующем виде.
1.
Проверить управляемость разомкнутой системы путем вычисления ранга матрицы
управляемости: Fo ( B, AB,..., An 1B) .
2.
Если разомкнутая система управляема, то определить характеристический
полином A ( ) n an 1n 1 ... a0 , матрицы A, B канонической формы Фробениуса - Калмана,
матрицу Fo ( B , A B ,..., A n 1B ) и передаточную
функцию Wo ( s) C ( sE A) 1 B
3.
4.
p s n 1 ... p0
P( s )
n n 1 n 1
.
A (s) s an 1s ... a0
1
Найти матрицу преобразования T Fo Fo .
Записать желаемый характеристический полином замкнутой системы:
d ( ) n n 1n 1 ... 0 .
5.
Определить элементы вектора обратной связи K для канонической формы, на основе
соотношений: k j j a j ; i 1,2,.., n .
6.
Вычислить вектор обратной связи для исходной системы: K KT . Записать передаточную
функцию замкнутой системы в форме: Wc ( s) C ( sE A B K ) 1 B
pn s n 1 pn 1s n 2 ... p1
.
d ( s)
Как мы видим, определение вектора K обратной связи сводится к поиску решения уравнения:
det(E A BK ) d ( ) .
Лемма (формула Аккермана).
Если пара матриц ( A, B) является полностью управляемой, то решение
уравнения det(sE A BK ) d (s) определяется следующим соотношением:
K [0,0,...,0,1][ B, AB,..., An 1B]1d ( A) .
MIMO LTI модель замкнутой системы с обратной связью по состоянию.
Рассмотрим MIMO LTI модель системы: x Ax Bu; y Cx; x R n ; u R m , y R p , m 1, p 1 .
50
Теорема.
Если пара матриц ( A, B) является полностью управляемой и, задан желаемый полином:
d (s) s n n 1s n 1 ... 0 с вещественными коэффициентами, то тогда существует
обратная связь по состоянию: u(t ) Kx(t ) , где K R m n - вещественная постоянная матрица,
которая, для замкнутой системы: x ( A BK ) x Bg ; y Cx , позволяет получить
характеристический полином равный желаемому, то есть: det(sE ( A BK )) d (s) .
Функции среды Matlab для синтеза обратной связи по состоянию.
В среде Matlab имеются две функции ac ker(.) и place (.) , которые позволяют рассчитывать
параметры обратной связи по состоянию.
Команда K ac ker( A, B, p) предназначена для вычисления вектора обратной связи K дляSISO LTI
полностью управляемых моделей, задаваемых парой матриц ( A, B) . Желаемый характеристический
полином задается множеством желаемых собственных значений матрицы ( A BK ) , которое
задается вектором - строкой p .
Пример.
s 0.5
Передаточная функция разомкнутой системы имеет вид: Wo ( s)
. Построить
s 3 2s 2 3s 0.5
обратную связь по состоянию, которая будет обеспечивать желаемые полюса (1;1.5 i;1.5 i)
для замкнутой системы.
Файл E43_Acker.
Функция place (...) предназначена для вычисления матрицы K обратной связи по состоянию для
MIMO LTI полностью управляемых моделей систем. Команда K place ( A, B, p) , где пара
матриц ( A, B) является полностью управляемой, а вектор – строка p , определяет желаемые
собственные значения матрицы ( A BK ) , позволяет вычислять матрицу K . Следует отметить, что
вектор p не может содержать кратные собственные значения (функция, в этом случае, формирует
сообщение об ошибке).
51
LTI модели систем, стабилизируемые обратной связью по состоянию.
Если модель состояния является полностью управляемой, то с помощью обратной связи по
состоянию можно расположить все собственные значения матрицы замкнутой системы в заданной
области. Теперь рассмотрим случай, когда модель состояния не является полностью управляемой.
Рассмотрим систему: x Ax Bu; y Cx; x R n ; u R m , y R p , которая имеет матрицу
управляемости F ( B, AB,..., An 1B), rank ( F ) r n . Тогда, используя декомпозицию по Калману,
zc Ac
zuc 0
можно уравнение состояния записать в эквивалентной форме:
A12 zc Bc
u ,
Auc zuc 0
где zc R r ; zuc R ( n r ) ; Ac R r r ; A12 R r ( n r ) ; Auc R ( n r )( n r ) ; Bc R ( r m) . Так как пара ( Ac , Bc )
является управляемой, то существует такой закон управления u (t ) K c zc , что матрица ( Ac Bc K c )
замкнутой управляемой части системы, будет иметь желаемое распределение собственных
значений. Тем не менее, необходимо, чтобы вся замкнутая система обладала внутренней
устойчивостью. То есть, матрица Auc должна быть гурвицевой, или характеристический
полином det(sE Auc ) должен иметь нули, имеющие отрицательные вещественные части.
Определение.
MIMO LTI модель системы: x Ax Bu; y Cx; x R n ; u R m , y R p , называется
стабилизируемой, если существует закон управления u(t ) Kx(t ) такой, что система
является асимптотические устойчивой.
Теорема.
Если MIMO LTI модель системы не является полностью управляемой, то она будет
стабилизируемой, тогда и только тогда, когда матрица Auc неуправляемой части системы
будет гурвицевой, то ее собственные значения будут иметь отрицательную вещественную
часть.
Напомним, что матрица обратной связи для исходной модели системы определяется с помощью
соотношения: K KP.
Пример.
Модель состояния вертолета описывается уравнением: x
Ax Bu; y Cx; x R 6 ; u R; y R . Построить
обратную связь по состоянию таким образом, чтобы замкнутая система была стабилизируемой.
52
Файл E45_StabNCtr.
Синтез LTI SISO систем с обратной связью.
Рассмотрим, следующую LTI SISO модель состояния разомкнутой системы вида:
x Ax Bu; y Cx; x R n , y, u R . Как известно, если пара ( A, B) является полностью
управляемой, то существует такой вектор K R1 n обратной связи по состоянию, что замкнутая
система: x ( A BK ) x Bg , где g R - задающее воздействие, будет иметь желаемый
характеристический полином d (s) . Кроме того, если пара ( A, C ) будет полностью наблюдаемой,
то передаточная функция замкнутой системы является несокращаемой и имеет вид:
Wc ( s) C ( sE ( A BK )) 1 B
Wo ( s) C ( sE A) 1 B
P( s )
, где полином P(s) равен числителю передаточной функции
d (s)
P( s )
разомкнутой системы.
det(sE A)
53
Таким образом, чтобы переходной процесс в замкнутой системе удовлетворял заданным
требованиям качества, необходимо исследовать передаточную функцию Wc ( s )
P( s )
, где
d ( s)
полином P(s) уже известен. Поэтому, для выбора желаемой области полюсов можно использовать
метод желаемой передаточной функции с кратными вещественными полюсами, которую можно
записать в форме: Wd ( s)
P( s )
n
( s ) n P(0)
. Отсюда, нетрудно определить возможный интервал
[ min, max] . Если выбрать некоторую внутреннюю точку интервала * [ min, max] , то
тогда желаемый полином будет определяться, как d ( s) ( s ) n .
Данное условие возможно расширить, если выбрать более общий вид желаемой передаточной с
кратными полюсами.
После нахождения желаемого полинома d (s) можно найти вектор коэффициентов обратной
связи путем использования функций ac ker(.) или place (.) среды Matlab.
Пример.
Рассмотрим LTI модель судна на подводных крыльях, заданную матрицами ( A, B, C ) , где:
Построим обратную связь по состоянию вида:
процесс по выходу
10%; t s 25s .
u3 Kx g такую, что, замкнутая система имела переходной
y2 (t ) при задающем воздействии g 1[t ] со следующими показателями качества:
54
Файл E57_FBhydro.
55
Синтез матрицы усилений замкнутой системы.
Рассмотрим обратную связь по состоянию для LTI MIMO системы в форме: u Kx Hg , где K матрица обратной связи соответствующей размерности, H - матрица усилений размера (m p) и
имеющая ранг rank ( H ) p , и g R p вектор входных задающих воздействий. LTI модель
разомкнутой системы имеет вид: x Ax Bu; y Cx; x R n , u R m , y R p . Тогда уравнение
замкнутой системы можно записать в виде: x ( A BK ) x BHg; y Cx . Передаточная функция
замкнутой системы равна: Wc ( s) C ( sE ( A BK )) 1 BH . Заметим, что обеспечить требования к
установившемуся значению каждого выхода можно только тогда, когда размерность задающего
вектора входных воздействий g имеет размерность равную p / 2 /. Поэтому, обычно, должно
выполняться ограничение m p , то есть в разомкнутой системе число управляющих входов должно
быть больше, чем число измеряемых выходов. Будем также предполагать, что LTI модель
разомкнутой системы будет полностью управляемой и наблюдаемой. Так как замкнутая система
должна быть асимптотически устойчивой, а преобразование Лапласа для вектора задающих
воздействий в виде единичной функции 1[t ] для j - ой компоненты, может быть записано в форме:
1
, где e j (0,0,.,0,1,0,..0)T R p , то тогда для соответствующего j - го выхода имеет
s
1
место соотношение: y j () lim sY j ( s) lim sWc ( s)e j Wc (0)e j , j 1,2,..., p , а матрица
s 0
s 0
s
L{g} e j
Wc (0) C ( A BK ) 1 BH будет определять коэффициенты усиления замкнутой системы в
установившемся режиме. Выберем матрицу H таким образом, чтобы выполнялось соотношение
Wc (0) E , то есть позиционная ошибка в установившемся режиме для каждой компоненты выхода
равнялась нулю (или замкнутая система имеет астатизм первого порядка по каждой
соответствующей компоненте выхода при воздействии соответствующей компоненты задающего
входного воздействия, когда другие компоненты входа равны нулю).
Предположим, что матрица коэффициентов усиления Wc (0) имеет размерность ( p p) . Тогда
найдем, что H [C ( A BK ) 1 B]1 . Если rank ( H ) p , то это означает, что коэффициент усиления
по какому – то выходу от соответствующего входного воздействия будет равен нулю.
Если же m p , то фактор [C ( A BK ) 1 B]1 будет иметь размерность ( p m) , то есть число
столбцов будет больше числа строк. В том случае, при выполнении условия rank ( H ) p можно
определить матрицу H с помощью процедуры псевдо - обращения Мора – Пенроуза, которая
задается соотношением: H [C ( A BK ) 1 B]T (C ( A BK ) 1 B [C ( A BK ) 1 B]T ) 1 . В этом случае
получим, что: Wc (0) E .
Пример.
Рассмотрим, повторно, LTI модель судна на подводных крыльях, которая задается матрицами
которая обсуждалась в предыдущем примере. Построим обратную связь по состоянию вида:
( A, B, C ) , и,
u3 Kx Hg ,
K . При этом потребуем, чтобы y2 () g , когда g L 1[t ], L const . Расчет
требуемого коэффициента (матрицы H ) приведен ниже.
с таким же вектором
56
Файл E58_Sthydro.
Восстановление фазовых координат объекта по измеряемому выходу.
Построение обратной связи по состоянию осуществлялось из условия предположения, что все
координаты модели состояния x(t ) R n доступны для измерения. К сожалению, для большинства
реальных технических систем, данное ограничение не является выполнимым, и доступным для
измерения является вектор выхода y(t ) R p . Поэтому возникает задача получения оценки x (t ) R n
вектора состояния по измерениям выхода y (t ) и входа u (t ) R m .
57
Динамическая подсистема, позволяющая, на основе измерения входа u (t ) и выхода y (t )
системы, получать оценку x (t ) R n вектора состояния, такую, что || x (t ) x(t ) || 0 при t ,
называется асимптотическим наблюдателем.
Если наблюдатель имеет такую же размерность, что и исходная модель состояния системы, то он
называется наблюдателем полной размерности.
Асимптотический наблюдатель полной размерности.
Рассмотрим SISO LTI модель исходной системы: x Ax Bu; y Cx; x R n ; u R; y R , у которой
невозможно измерение вектора состояния x(t ) R n . Если пара матриц ( A, B) известна, то можно
дублировать исходную модель, то есть построить систему вида x̂ Axˆ Bu . Если дублирующая
модель будет иметь такие же начальные условия, что и исходная система, то при любом входе u (t )
будет выполняться соотношение xˆ (t ) x(t ) для любого t 0 . Такое динамическое устройство для
оценки вектора состояния исходной системы называется наблюдателем с разомкнутым контуром
/ 3 /.
Однако, такая схема построения наблюдателя имеет существенные недостатки. Например, если
матрица A имеет собственные значения с положительной вещественной частью, то малейшее
различие в начальных условиях x(t0 ) и xˆ (t0 ) будет приводить к тому, что ошибка || x (t ) x(t ) || будет
расти экспоненциально во времени.
Поэтому рассмотрим задачу построения асимптотического наблюдателя полного порядка с
обратной связью по выходу, то основанного на ошибке между сигналом выхода исходной системы
(объекта) Cx(t ) и сигналом выхода наблюдателя: Cxˆ (t ) . Обозначим вектор коэффициентов такой
обратной связи, как L R n1 , тогда схема асимптотического наблюдателя с обратной связью будет
описываться уравнением: x̂ Axˆ Bu L( y Cxˆ ); xˆ R n ; u R; y R или x̂ [ A LC]xˆ Ly Bu .
58
Обозначим вектор ошибки оценки состояния, как e x xˆ . Тогда уравнения, описывающие
исходную систему (объект), и наблюдатель, можно записать в виде: x Ax Bu; e [ A LC]e . Таким
образом, если матрица [ A LC] имеет собственные значения с отрицательной вещественной
частью, то все компоненты вектора ошибки оценки состояния будут стремиться к нулю, то есть:
|| x (t ) x(t ) || 0 при t .
Теорема.
Если SISO LTI модель системы : x Ax Bu; y Cx; x R n ; u, y R является полностью
наблюдаемой, то тогда существует подсистема асимптотического наблюдателя полной
размерности с обратной связью, характеристический полином которой равен заданному
желаемому полиному, то есть det(sE A LC) o ( s) , где o ( s) s n n 1s n 1 ... 0 желаемый полином с вещественными коэффициентами.
Рассмотри теперь MIMO LTI модель: x Ax Bu; y Cx; x R n ; u R m , y R p ; m 1 . Рассмотрим
схему асимптотического наблюдателя полной размерности с обратной связью по выходу.
Уравнение наблюдателя имеет вид: x̂ Axˆ Bu L(Cxˆ Du y) ( A LC) xˆ ( B LD)u Ly ,
где xˆ (t ) - выход наблюдателя, L - матрица обратной связи наблюдателя, L R n p . Обозначим
ошибку оценки вектора состояния, как: e xˆ x . Тогда уравнение наблюдателя можно записать в
форме: e ( A LC)e . Отсюда найдем решение: e e
( A LC )(t t 0 )
e(t0 ) . Если матрица L выбрана
таким образом, что матрица ( A LC) является гурвицевой, то тогда: lim [ xˆ (t ) x(t )] 0 . Кроме того,
t
59
можно выбрать матрицу L таким образом, что характеристический полином: det(sE A LC) будет
совпадать с заданным желаемым полиномом: o ( s) s n n 1s n 1 ... 0 / 1 /.
Теорема.
Если LTI модель системы: x Ax Bu; y Cx Du; x R n ; u R m , y R p является полностью
наблюдаемой и, задан желаемый полином o ( s) s n n 1s n 1 ... 0 с вещественными
коэффициентами, то существует асимптотический наблюдатель полной размерности:
x̂ ( A LC) xˆ ( B LD)u Ly такой, что матрица [ A LC] имеет характеристический полином:
det[sE ( A LC)] , который совпадает с желаемым полиномом: o (s) s n n 1s n 1 ... 0
Асимптотический наблюдатель неполной размерности.
Запишем полностью наблюдаемую SISO LTI модель в виде канонической формы наблюдаемости:
b1
0 0 ... a0
1
...
a
b
1
, C (0, 0, ..., 1) , B 2 .
x A x B u; y C x; x R n ; u, y R , где A
. . ...
...
...
0 0 ... a
b
n 1
n
Отсюда видно, что выход y (t ) совпадает с компонентой xn (t ) вектора состояния x (t ) .
Следовательно, при измерении выхода, необходимо получать оценки только компонент
x1 (t ),...xn1 (t ) .
Рассмотрим уравнение асимптотического наблюдателя полной размерности:
x̂ A xˆ B u L ( y C xˆ ); xˆ R n ; u R; y R . Здесь выход наблюдателя равен: C xˆ xˆn , а выход
исходной системы (объекта) y C x xn xˆ n . Выделим составляющую y xˆ n и запишем уравнение
~
~
~
x A~
x Bu a~n 1 y , где ~
x ( xˆ1,..., xˆn 1 )T ; A наблюдателя относительно (n 1) - го вектора ~
x: ~
~
подматрица размера (n 1) (n 1) , которая расположена в левом верхнем углу матрицы A ; B -
~
вектор (b1 , b2 ,..., bn 1 )T ; a
n 1 - вектор размера (n 1) , который совпадает с последним столбцом
~
матрицы Â . Определим следующий вектор ошибки: e~ x ~
x , e~ R n 1 . Отсюда найдем: e~ A e~ .
Рассмотрим невырожденное линейное преобразование x Pxˆ такое, чтобы матрица A PA P 1
имела заданный характеристический полином: o ( s) s n 1 n 2 s n 2 ... 0 . Для этого построим
1
0
матрицу P в следующем виде: P .
0
0
0
1
.
. Тогда обратная матрица равна:
0 ... 1 n 2
0 ... 0
1
0 ... 0
1 ... 0
. ... .
60
1
0
1
P .
0
0
0 ... 0
1 ... 0
. ... .
0
1
. . Так как пара матриц ( A , C ) имеет каноническую форму представления,
0 ... 1 n 2
0 ... 0
1
0
1
1 .
то отсюда получим: PAP
0
0
0
0 ... 0
0
0 ... 0
1
. ... .
.
0 ... 0 n 3
0 ... 1 n 2
0 ... 0
1
(an 1 n 2 ) 0 a0
(an 1 n 2 ) 1 a1 0
...
.
(an 1 n 2 ) n 3 an 3 n 4
(an 1 n 2 ) n 2 an 2 n 3
an 1 n 2
bn 0b1
bn 1 1b1
Соответственно, первые (n 1) компонент вектора PB имеют вид:
. При этом,
...
b b
n2 1
2
характеристический полином матрицы PA P 1 совпадает с заданным полиномом:
o (s) s n 1 n 2 s n 2 ... 0 .
Теорема.
Пусть SISO LTI модель системы: x Ax Bu; y Cx; x R n ; u, y R является полностью
наблюдаемой, тогда существует асимптотический (n 1) - мерный наблюдатель, который
имеет заданный характеристический полином o ( s) s n 1 n 2 s n 2 ... 0 . Если пара
матриц ( A, C ) имеет представление в виде канонической формы наблюдаемости, то тогда
выход системы представляет собой n - ый компонент вектора состояния.
Такой асимптотический (n 1) - мерный наблюдатель называется наблюдателем Люенбергера.
Наблюдатель Люенбергера может быть также построен и для MIMO LTI модели.
Теорема.
Для полностью наблюдаемой MIMO LTI модели: x Ax Bu; y Cx; x R n ; u R m , y R p
существует асимптотический (n p) - мерный наблюдатель, который имеет заданный
характеристический полином (n p) - го порядка, где p rank (C )
Рассмотрим алгоритм построения наблюдателя Люенбергера (n p) - мерной размерности.
Пусть матрица C имеет ранг равный p . Выделим p линейных комбинаций y Cx , которые
определяются выходом системы. Оставшиеся (n p) независимых компонент вектора z определим
с помощью соотношения z Lx , где матрица L R (n p) n . То есть, определим линейное
z L
x Px . Очевидно, что в силу построения, матрица P
y C
преобразование соотношением:
61
z
y
является обратимой (несингулярной). Тогда можно записать следующие выражения: x P 1
z
z
PAP 1 PBu . Последнее уравнение можно записать, с помощью блочной структуры
y
y
и
z Azz
y
Ayz
матриц, в виде:
Azy z Bz
u . Первое матричное уравнение имеет вид:
Ayy y B y
z Azz z Azy y Bz u , где слагаемое ( Bzyu ) y Azy y Bz u является измеримым, так как выход y Cx
и вход u являются измеримыми. Тогда уравнение наблюдателя для получения оценки вектора z
можно записать в виде zˆ Azz zˆ Azy y Bz u; zˆ(0) 0 . В этом случае, вектор ошибки по
выходу ez (t ) z (t ) zˆ(t ) будет удовлетворять уравнению: ez Azzez ; ez (0) z (0) zˆ(0) . Оценка
zˆ
. Если матрица L выбрана таким
y
вектора состояния исходной системы определяется, как xˆ P 1
образом, чтобы матрица: Azz R ( n p)( n p) была гурвицевой, то тогда найдем:
zˆ
e x xˆ P 1 0 при t . Структура (n p) - мерного наблюдателя Люенбергера
0
приведена ниже на рисунке.
Обычно, при построении наблюдателя не полной размерности в качестве координат выхода y
используют соответствующие компоненты вектора состояния. Это приводит к упрощению расчетов.
Матричное уравнение TA DT S .
Рассмотрим LTI однородную модель системы: x Ax; y Cx; x R n ; y R p . Пусть выход системы
является входом наблюдателя, заданного в виде следующей динамической модели:
z Vz LCx; z R q , где L R q p , V R q q . Будем также считать, что пара матриц (V , L) является
управляемой, а матрица V является гурвицевой. Обозначим LC S и запишем уравнения
наблюдателя в виде: z Vz Sx . Предположим, что существует матрица линейного
преобразования T R q n такая, что имеет место равенство: TA VT S , то есть z TAx Tx .
Тогда, если матрица V не имеет общих собственных значений с матрицей A , то решение
уравнения z Vz Sx может быть записано в форме: z (t ) Tx(t ) eVt [ z (0) Tx(0)] . Так как
матрица V является гурвицевой, то тогда: z (t ) Tx(t ) при t для любых z (0), x(0) .
Если матрица C имеет ранг равный p , то тогда решение T уравнения TA VT S имеет ранг
равный (n p) .
62
Рассмотрим случай, когда матрица B 0 , то есть система имеет вид:
x Ax Bu; y Cx; x R n ; u R m , y R p и, пусть существует решение T уравнения TA VT S , где
z Tx , тогда уравнение асимптотического наблюдателя может быть записано в форме:
z Vz Sx TBu; z R q . При этом, решение z имеет тот же самый вид: z(t ) Tx(t ) eVt [ z(0) Tx(0)] .
Отсюда, можно сформулировать следующее утверждение / 1 /.
Теорема.
Пусть задана LTI модель системы x Ax Bu; y Cx; x R n ; u R m и асимптотический
наблюдатель, заданный уравнением: z Vz Sx TBu; z R q , u R m , x R n ,где V является
гурвицевой матрицей, не имеющей общих собственных значений с матрицей A ; матрица
S LC и пара матриц (V , L) является полностью управляемой. Тогда выход z наблюдателя
будет стремиться к вектору Tx при t для любых начальных условий z (0), x(0) , тогда и
только тогда, когда матрица T является решением матричного уравнения TA VT S .
Пусть пара матриц pair ( A, C ); p rank (C ) является полностью наблюдаемой. Тогда можно
детализировать алгоритм построения асимптотического наблюдателя редуцированной
размерности.
1.
Выберем, некоторым заданием, устойчивую матрицу V R ( n p)( n p) , которая не имеет
общих собственных значений с матрицей A .
2.
Построим такую матрицу L R ( n p) p , чтобы пара матриц (V , L) была управляемой.
3.
Найдем решение T R ( n p)n уравнения TA VT LC .
4.
Если квадратная матрица P порядка n является вырожденной (сингулярной), то
C
T
вернуться к шагу 2. Если матрица P является несингулярной, то тогда (n p) - мерный
1
C y
асимптотический наблюдатель z Vz TBu Ly , xˆ позволяет получать оценку вектора
T z
состояния системы x .
Матричное уравнение VX XA Q 0 , где матрица X является решением, часто называют
расширенным матричным уравнением Ляпунова.
Функции среды Matlab для синтеза наблюдателей.
Функция ac ker(.) позволяет рассчитывать вектор L асимптотического наблюдателя полной
размерности для SISO LTI полностью наблюдаемой модели системы. Для этого используется
следующий синтаксис команды: L accer ( AT , C T , p)T .
Функция place (.) позволяет рассчитывать матрицу L асимптотического наблюдателя полной
размерности для полностью наблюдаемой MIMO LTI модели системы. Синтаксис команды имеет
вид: L place ( AT , C T , p)T , где p - вектор – строка, определяющий желаемые нули
характеристического полинома матрицы ( A LC) наблюдателя.
Функция est estim(sys, L, sensors, known) позволяет построить модель наблюдателя на основе
заданной матрицы L , измеряемых компонент входа u , измеряемых компонент выхода y , не
измеряемых компонент выхода z и не измеряемых возмущений w . То есть модель системы
63
x Ax B1w B2u
z
D
D . Вектора sensors и known задают индексы
задается в виде: C1
x 11 w 12 u
y C2 D21
D22
измеряемых входов и выходов. Модель наблюдателя формируется в форме:
xˆ Axˆ B1w B2u L( y C2 xˆ D22u )
.
yˆ C2 D22
ˆ
x
u
xˆ E 0
Функция X lyap (V , A, Q) позволяет решать расширенное уравнение Ляпунова. Матрицы V , A, Q
должны иметь согласованные размеры, но не обязательно должны быть квадратными.
Пример.
LTI модель судна на подводных крыльях задана матрицами ( A, B, C ) :
Построить асимптотической наблюдатель полной размерности и наблюдатель редуцированной
размерности.
Матрица наблюдаемости H (C T , AT C T , ( AT ) 2 C T , ( AT )3 C T ) имеет ранг 4. Вектор желаемого
распределения нулей характеристического полинома: p [0.6 0.61 0.62 0.69] . Динамика
сигнала ошибки: Erf 1 xˆ1 x1 и Erf 2 xˆ3 x3 показана на рисунке 1.
Ранг матрицы C равен 2, поэтому наблюдатель редуцированной размерности имеет порядок
равный 2. Матрица Ld R 2 2 определяет обратную связь наблюдателя. Матрица преобразования
T ищется из решения матричного уравнения TA VT LC . Динамика сигнала
ошибки Erd z Tx R 2 показана на рисунке 2.
Файл E46_ObvRed.
64
65
Обратная связь на основе оценок состояния.
Пусть задана полностью управляемая и полностью наблюдаемая LTI система:
x Ax Bu; y Cx Du; x R n ; u R m , y R p и , построена матрица K обратной связи по
состоянию такая, что характеристический полином матрицы замкнутой системы ( A BK ) совпадает
с заданным полиномом d ( ) s n n 1s n 1 ... 0 . Кроме того, пусть построен асимптотический
наблюдатель полной размерности x̂ ( A LC) xˆ ( B LD)u Ly; u g Kx такой, что
характеристический полином матрицы ( A LC) совпадает с полиномом:
o (s) s n n 1s n 1 ... 0 . Подставим в обратную связь Kx вместо вектора x его оценку x̂ на
выходе наблюдателя. Тогда уравнения замкнутой системы с наблюдателем можно записать в
форме: x Ax BKxˆ Bg , x̂ ( A LC) xˆ B( g Kxˆ ) LCx LDg .
Теорема.
Пусть задана полностью управляемая и полностью наблюдаемая LTI модель системы ( A, B, C ) .
Предположим, что построена матрица K обратной связи такая, что матрица ( A BK ) имеет
характеристический полином равный: d ( ) s n n 1s n 1 ... 0 . Предположим также, что
построена матрица L асимптотического наблюдателя такая, что матрица ( A LC) имеет
характеристический полином равный: o ( s) s n n 1s n 1 ... 0 . Обозначим через 2n ( s)
характеристический полином замкнутой системы с наблюдателем, заданной уравнениями:
66
x Ax BKxˆ Bg , x̂ ( A LC) xˆ B( g Kxˆ ) LCx LDg . Тогда имеет место равенство:
2n (s) d (s)o (s) .
Доказательство.
Запишем уравнения замкнутой системы с наблюдателем в следующем виде:
x A
̂
x LC
BK
x B
g . Обозначим: ~
x x xˆ и запишем линейное
A LC BK xˆ B LD
x
x
E
преобразование: ~ M , где M
E
x
xˆ
0
, E R n n . В новых координатах уравнения
E
x A BK
0
замкнутой системы примут вид: ~
x
BK x B
g . Отсюда видно, что
A LC ~
x LD
характеристический полином такой системы имеет вид: 2n ( s) d ( s)o (s) .
Так как построение обратной связи по состоянию и асимптотического наблюдателя
осуществлялось независимо друг от друга, то не имеет значения и порядок выполнения расчетов
матриц K и L . Такой подход к проектированию систем управления называется принципом
разделения.
Рассчитаем передаточную функцию Wc (s) замкнутой системы с наблюдателем. Имеют место
1
BK
sE A BK
B
следующие соотношения: Wc ( s) C 0
sE A LC 0
( sE A BK ) 1 ( sE A BK ) 1 BK ( sE A LC ) 1 B
C ( sE A BK ) 1 B P( s)o ( s) ,
C 0
0
d ( s)o ( s)
( sE A LC ) 1
P( s )
где
C ( sE A BK ) 1 B передаточная функция замкнутой системы без наблюдателя. Таким
d ( s)
образом, матрица L оказывает влияние на переходной процесс в системе только при ненулевых
начальных условиях x(0), ~
x (0) 0 . Так как, в передаточной функции Wc (s) имеются сокращаемые
полюса и нули, то замкнутая система с наблюдателем будет неуправляемой и/или ненаблюдаемой.
Пример.
LTI модель судна на подводных крыльях задана матрицами ( A, B, C ) :
Построить регулятор, в виде обратной связи по состоянию, с заданным распределением собственных
значений p [0.6 0.61 0.62 0.69] , и асимптотического наблюдателя полной размерности, с заданным
распределением собственных значений q [0.7 0.71 0.72 0.79] . Вычислить передаточную функцию
замкнутой системы по входу
u3 .
67
Файл E47_ModCnt.
68
Синтез линейных оптимальных систем с интегральным квадратичным функционалом.
Рассмотрим следующую линейную систему:
x A(t ) x B(t )u; y C T (t ) x; x(t0 ) x0 .
Квадратичный критерий качества определим с помощью следующего соотношения:
tf
J {u T ( )u ( ) xT ( ) L( ) x( )}d xT (t f )Qx(t f ) .
Задача синтеза состоит в поиске
t0
управляющей функции
u (t ) , определенной на интервале t0 t t f , которая определена на
траекториях движениях системы и минимизирует функционал
J . Матрицы Q, L , обычно, полагают
симметричными. Данная задача является вариационной и, часто, рассматривается в разделе
теории оптимального управления. Однако, если в качестве управляющих функций рассмотреть m мерное пространство {C[t 0 , t f ]}
m
непрерывных функций, то можно решить указанную задачу
геометрическими способами, используемыми при исследовании линейных систем /Андреев с359/.
Можно выделить два вида решений приведенной задачи.
В первом случае оптимальное управление u (t ) определяется как функция от матриц A, B, Q, L ,
начального состояния x0 ,t 0 , временного параметра t и множества значений, которому должно
69
принадлежать состояние x1 , t f . То есть, управление рассматривается для разомкнутого контура
системы, полностью определяется априори заданными параметрами задачи и не зависит от
текущего состояния системы.
Во втором случае, управление ищется в виде линейной функции от текущего состояния ( x(t ), t )
при движении системы по оптимальной траектории. То есть ищется оптимальное управление в
виде u(t ) K (t ) x(t ) . В этом случае управление представлено в виде сигнала обратной связи
или, говорят, управление осуществляется по замкнутому контуру.
При решении поставленных задач будем предполагать, что состояние системы, либо ее выход
доступны для измерения. Рассмотрим предварительно задачу без ограничений на конечное
состояние системы. Определим матричное дифференциальное уравнение, называемое
уравнением Риккати, с помощью соотношения:
K (t ) AT (t ) K (t ) K (t ) A(t ) K (t ) B(t ) B T (t ) K (t ) L(t )
Теорема (о существовании оптимального управления для линейной системы без ограничений на
конечное состояние). Пусть заданы матрицы A(t ), B(t ), L(t ), Q(t ) , где L(t ), Q(t ) симметричные матрицы. Предположим, что на интервале t 0 t t f существует решение K (t )
дифференциального уравнения Риккати
K (t ) AT (t ) K (t ) K (t ) A(t ) K (t ) B(t ) BT (t ) K (t ) L(t ) с начальным условием K (t1 ) Q . Тогда
существует управление u (t ) , которое дает минимум критерию качества:
tf
J {u T ( )u ( ) x T ( ) L( ) x( )}d x T (t f )Qx (t f ) на траекториях системы
t0
x A(t ) x B(t )u; x(t0 ) x0 . Минимальное значение J равно xT (t 0 ) K (t ) x(t 0 ) . Оптимальное
управление в виде обратной связи определяется следующим соотношением:
u(t ) BT (t ) K (t , Q, t f ) x(t ) . Оптимальное управление по разомкнутому контуру имеет вид:
u(t ) BT (t ) K (t , Q, t f )(t , t 0 ) x0 , где (t , t 0 ) - переходная матрица для уравнения:
x [ A(t ) B(t ) BT (t ) K (t , Q, t f )] x .
Решение задачи в виде обратной связи называется построением регулятора нулевого состояния.
Структура замкнутой системы, при этом, будет иметь следующий вид:
Пример.
Пусть
x(0) 1 . Необходимо найти x(t )
на интервале
0 t T так, чтобы критерий, задаваемый
T
соотношением:
J {[ x ( )]2 x 2 ( )}d
t0
достигал минимума. В этом случае, уравнение движения можно
70
записать в следующем виде:
x (t ) u(t ) . Тогда уравнение Риккати для этой задачи имеет вид:
k(t ) k 2 (t ) 1 . Решение, которое проходит через нуль при t T , описывается функцией
k (t ) th(T t ) . Поэтому оптимальная траектория удовлетворяет уравнению: x (t ) th(T t ) x(t ) .
sh (T )
Интегрируя это уравнение, получим следующие соотношения: x(t ) ch(t ) (
) sh (t ) ,
ch(T )
u(t ) th(T t ) x(t ) . Минимальное значение квадратичной оценки качества равно:
x T (0) K (0,0, T ) x(0) x 2 (0)th(T ) .
Достаточно часто, требуемое ограничение на управление, которое в явной форме не учитывается
в постановке задачи оптимизации, может быть выбором соответствующей матричной функции
G(t ) .
Теорема. Пусть задана система
x A(t ) x B(t )u h(t ); y C T (t ) x; x(t0 ) x0 ; x R n ; u R m ; y R P и критерий оптимизации
tf
вида
J {u T ( )G ( )u ( ) x T ( ) L( ) x( )}d x T (t f )Qx(t f ) . Тогда решение задачи синтеза
t0
оптимального нестационарного линейного регулятора существует и единственно, а
оптимальное управление имеет вид u * (G 1 B T Kx
матрица
K
и
n - мерный вектор p
1 1 T
G B p) , где симметричная (n m)
2
определяются из уравнений
K KA AT K KBG 1BT K L ; p KBG 1 BT p AT p 2Kh при граничных условиях
K (t f ) Q ; p(t f ) 0 . При этом для любого t [t0 , t f ] справедливо соотношение
tf
x (t ) K (t ) x(t ) p (t ) x(t ) q(t ) x (t f )Qx(t f ) [( x * )T L( ) x * ( ) (u * )T G( )u * ( )]d , где
T
T
T
t
q(t )
- скалярная функция, которая определяется из уравнений
q
1 T
p BG 1 BT p p T h; q(t f ) 0 .
4
h(t ) 0 оптимальное управление имеет вид u * G 1BT Kx , где
(n n) матрица K определяется из уравнения
Следствие. При
симметричная
K KA AT K KBG 1BT K L ; K (t f ) Q .
Рассмотрим теперь случай линейной стационарной системы вида:
x Ax Bu; y C T x; x(t0 ) x0 ; x R n ; u R m ; y R p .
Теорема (о связи управляемости и идентифицируемости системы с решением алгебраического
уравнения Риккати ). Если система { A, B, C} управляема и идентифицируема, то существует
71
симметрическая положительно определенная матрица V , которая является решением
алгебраического уравнения Риккати:
(без доказательства).
AT V VA VBBT V CC T 0
Таким образом, управляемость и идентифицируемость тройки матриц
.
{ A, B, C} является
необходимым условием существования соответствующего алгебраического уравнения Риккати.
Теорема (об оптимальном стационарном регуляторе). Пусть стационарная система
управляема и идентифицируема. Пусть также
{ A, B, C}
Vr - положительно определенное решение
алгебраического матричного уравнения: A V VA VBB V CC 0 . Тогда существует
управление, которое минимизирует квадратичный критерий качества следующего вида:
T
J {u T ( )u ( ) x T ( )CC T x( )}d
для системы
T
T
x Ax Bu; y C T x; x(t0 ) x0 ; . Причем
t0
минимальное значение
имеет вид:
J
равно
x0T Vr x0 , а оптимальное
u(t ) BT Vr x(t ) , и в
Доказательство.
разомкнутой форме:
управление в виде обратной связи
T V ]t
r
u (t ) BT Vr e[ A BB
.
Сделаем замену переменной: v(t ) u (t ) B Vr x(t ) , называемую
T
преобразованием Риккати. Уравнения движения примут вид: x (t ) [ A BB Vr ]x(t ) Bv(t ) , и
T
функционал J будет определяться соотношением:
J [v( ) B T Vr x( )]T {[v( ) B T Vr , x( )] x T ( )CC T x( )}d . Преобразуем подынтегральное
t0
выражение, применив равенство:
AT Vr Vr A Vr BBT Vr CC T . Получим, после всех
преобразований, следующее соотношение:
J {[v T ( )v( ) [( A BBT Vr ) x( ) Bv( )]T Vr x( ) x T ( )Vr [( A BBT Vr ) x( ) Bv( )]}d
t0
или J v ( )v( )d x (t )Vr x(t ) |t0 . Очевидно, что минимум значения J достигается при
T
T
t0
v(t ) u(t ) BT Vr x(t ) 0 . Отсюда следует доказательство утверждения теоремы.
Структурная схема стационарного оптимального регулятора приведена на рисунке.
Имеет место также и более общий результат. Пусть линейная стационарная система описывается
уравнением: x Ax Bu , а критерий оптимальности имеет вид:
72
J {u T ( )Gu ( ) x T ( )Qx( )}d
(*)
Теорема (об оптимальном стационарном регуляторе). Задача (*) имеет решение тогда и
только тогда, когда стационарная система { A, B} стабилизируема, а оптимальное управление
имеет вид
u(t ) G 1BT Vr x(t )
, где
алгебраического матричного уравнения:
- положительно определенное решение
Vr
AT V VA VBG1BT V Q . При этом для любого
t
t 0 справедливо равенство x (t )Vr x(t ) {u T ( )Gu ( ) xT ( )Qx( )}d .
T
Линейные стохастические системы и их анализ
Линейные системы со случайной правой частью
Введенные понятия с.к. производной и с.к. интеграла позволяют дать корректное описание
линейного дифференциального уравнения со случайными возмущениями в правой части и
случайными начальными условиями: t a(t )t b(t )t , t t0 ,t 0 v , где t , t - с.к. непрерывные
случайные функции, t - с.к. производная функции t , a(t ), b(t ) - непрерывные неслучайные
функции, а v - некоторая случайная величина с известным распределением.
Случайная функция t является решением дифференциального уравнения при t t0 , если будет
t
t
t0
t0
удовлетворять уравнению: t v a( ) d b( ) d , где все интегралы понимаются в с.к.
смысле.
Для явного построения решения, введем вспомогательную неслучайную функцию (t , t0 ) ,
(t , t0 ) a(t ) (t , t0 ), t t0
. Тогда случайная
(t0 , t0 ) 1
удовлетворяющую следующим соотношениям:
t
функция t (t , t0 )v (t , )b( ) d удовлетворяет исходному дифференциальному уравнению
t0
со случайной правой частью.
Рассмотрим теперь векторный случай, когда система состоит из n линейных дифференциальных
уравнений, записанных в матричной форме: t A(t )t B(t )t , t t0 ,t 0 v , где t - с.к.
производная n - мерной случайной функции t R n , t R m - с.к. непрерывная при t t0
случайная функция, v R n - случайный вектор, а матричные функции A(t ) R n n , B(t ) R n m
имеют неслучайные кусочно - непрерывные компоненты .
t
Общее решение уравнения t A(t )t B(t )t ,t 0 v имеет вид: t (t , t0 )v (t , ) B( ) d ,
t0
где (t , t0 ) - матричная переходная функция, удовлетворяющая системе однородных
A(t )
дифференциальных уравнений
, t t0 . Интеграл в правой части решения понимается
(t 0 , t0 ) E
в с.к. смысле и применяется к случайной векторной подынтегральной функции по каждому ее
компоненту.
73
Преобразование стационарного процесса линейной стационарной системой.
Пусть стационарная случайная функция t имеет спектральное разложение вида:
t
e
it
dZ (d ) , где Z (d ) некоторая случайная мера с некоррелированными приращениями,
имеющая спектральную функцию F ( ) .
Определение.
Если случайная функция t , t R , допускает представление вида: t
e
it
( ) Z (d ) , где
комплексная функция ( ) удовлетворяет условию
| ( ) |
2
dF ( ) , то говорят, что
случайная функция t получена из случайной функции t с помощью стационарного линейного
преобразования.
Функция ( ) называется частотной характеристикой этого преобразования. Несложно
убедиться, что тогда будут справедливы следующие соотношения: D
| ( ) |
2
dF ( ) ,
F ( )
| (v ) |
2
dF (v); , v R . Если же, спектральная функция случайной функции t будет
иметь спектральную плотность S ( ) , то случайная функция t будет иметь спектральную
плотность S ( ) и, эти спектральные плотности связаны равенством: S ( ) | ( ) |2 S ( ) .
Процесс t будет стационарным, если выполняется неравенство
| ( ) |
2
dF ( ) . Тогда
корреляционная функция процесса t определяется с помощью следующего соотношения:
R (t )
e
it
| ( ) |2 dF ( ) .
Рассмотрим теперь стационарное линейное преобразование во временной области. Пусть задана
функция w(t ) такая, что
| w(t ) | dt . Рассмотрим следующее с.к. интегральное преобразование
стационарной центрированной случайной функции t
w(t s) s ds .
Очевидно, что будут
справедливы следующие соотношения:
M {t }
w(t s)M { s }ds const ,
M {| t |2 }
w(t s1 ) R ( s2 s1 ) w(t s2 )ds1ds2 D ( | w( s) | ds) 2 .
74
e
Если функция t имеет спектральное разложение t
it
dZ (d ) , то тогда будет также
выполняться следующее соотношение: t
( )
e
it
w(t s)( e
it
Z (d ))ds
e
it
( ) Z (d ) , где
w( s)ds . Причем будет справедливо неравенство
| ( ) |
2
dF ( ) M {| t |2 } .
Следовательно, преобразование t
w(t s) s ds является стационарным линейным
преобразованием с частотной характеристикой ( )
e
it
w( s)ds .
Пусть теперь стационарный процесс t является m раз с.к. дифференцируем. Пусть также
существует n раз с.к. дифференцируемый процесс t , связанный с t линейным
n
дифференциальным уравнением вида:
akt(k )
k 0
m
bk t(k ) . Найдем условия, при которых процесс
k 0
t будет стационарным случайным процессом, а также определим частотную характеристику
преобразования t в t .
Обозначим ( )
B(i )
. Если существует спектральная плотность S ( ) , то тогда, как было
A(i )
показано выше, будет справедливо соотношение S ( ) | ( ) |2 S ( ) . По построению
случайный процесс t существует и является стационарным, если его дисперсия будет конечной, то
есть M {| t |2 }
B(i )
| A(i ) |
2
dF ( ) . Данное соотношение будет выполняться, если полином
A(s) не будет иметь положительных и чисто мнимых корней.
Теперь рассмотрим процесс прохождения через физически реализуемый линейный фильтр с
частотной характеристикой ( ) стационарного процесса nt с постоянной спектральной
плотностью, равной S n ( ) 1 . Спектральная плотность S ( ) процесса t на выходе фильтра,
тогда определяется следующим образом: S ( ) | ( ) | , где ( )
2
e
it
w( s)ds , w(t ) -
импульсная функция фильтра, удовлетворяющая условию абсолютной интегрируемости
| w(t ) | dt .
Анализ стационарных линейных систем в пространстве состояний при воздействии белого шума.
Рассмотрим линейную стационарную систему, описываемую следующим уравнением:
x Ax Bu; y Cx; x R n , u R q , y R p , x(t0 ) xt 0 . Пусть входное воздействие
u (t ) (u1, u2 ,..., uq )T представляет собой стационарную случайную векторную функцию. В этом
случае запись уравнений состояния будем понимать в среднем квадратичном смысле. В
75
соответствие с условием задания случайных функций, такое входное воздействие будет
полностью характеризоваться многомерным совместным распределением вида
Fut1 ,ut 2 ,..., utk ( z1, z2 ,...., zkq ) . Для задания стационарных в широком смысле случайных функций
достаточно рассмотреть первые два момента такой векторной функции и ее компонентов: вектор
средних значений mu M {u (t )} (mu1 , mu 2 ,..., mu q )T и матрицу авто- и взаимно- корреляционных
функций Ru (t ) ( Ru i u j (t ) M {ui0u 0j }; i, j 1,2,..., q , где M {ui0u 0j } M {(u mu )(u mu )T } .
Общее решение уравнений состояния можно записать в следующем виде :
t
xt (t , t 0 )v (t , ) B( )u ( )d , где (t , t0 ) - матричная переходная функция. Интеграл в правой
t0
части решения понимается в с.к. смысле и применяется к случайной векторной подынтегральной
функции.
Для постоянной матрицы A переходная функция имеет вид (t , t0 ) e
A(t t 0 )
.
Если перейти к центрированным значениям вектора состояния x 0 x(t ) mx (t ) и вектора
входного воздействия u 0 n(t ) mu (t ) , то уравнение состояния системы можно записать в
следующем виде: x 0 Ax 0 Bu 0 . Тогда, предполагая, что вектор состояния x(t ) также является
стационарным процессом, можно записать следующее соотношение:
M {x 0 ( x 0 )T } D x (t ) M {Ax 0 ( x 0 )T Bu 0 ( x 0 )T } M {x 0 ( x 0 )T AT x 0 (u 0 )T BT }
ADx (t ) Dx (t ) AT BM {u 0 ( x 0 )T } M {x 0 (u 0 )T }BT .
Учитывая, что входной сигнал и начальные условия обычно не коррелированны, можно записать
следующее соотношение:
M {u ( x ) } M {u
0 T
( xt00 )T } T (t , t0 )
t
M {u
t
( )(u ( )) }B G (t )d Ru ( ) BT G T (t )d .
T
T
T
t0
t0
Здесь G(t ) (t , t0 ) 1 ( , t0 ) / 4 /.
Рассмотрим случай, когда входной сигнал представляет собой белый шум с независимыми
компонентами, то есть u 0 N 0 (t ) и RN (t ) V (t ) (t ) , где V (t ) диагональная матрица
интенсивностей векторного единичного белого шума N 0 (t ) . Тогда, с учетом данного допущения,
t
получим: M {u 0 ( x 0 )T } V ( ) ( ) BT GT (t )d
t0
1
1
1
[V ( ) BT ] t V (t ) BT . Множитель
2
2
2
учитывает тот факт, что момент приложения - функции совпадает с верхним пределом
интегрирования. Аналогично можно показать, что M {x 0 (u 0 )T }
1
BV (t ) .
2
Тогда уравнение для соотношения M {x 0 ( x 0 )T } может быть записано в виде:
D x (t ) ADx (t ) Dx (t ) AT BV (t ) BT . Данное матричное уравнение называется дисперсионным
уравнением и является уравнением Риккати. Так как матрица Dx (t ) является симметричной
матрицей, то оно эквивалентно системе
n(n 1)
скалярных уравнений.
2
76
Таким образом, учет того, что белый шум является случайным процессом с независимыми
приращениями, то есть не дифференцируемым случайным процессом, привел к тому, что в
уравнении для корреляционной функции Dx (t ) появилось дополнительное слагаемое BV (t ) BT .
Линейные стохастические системы.
Рассмотрим стохастическую систему, которая описывается с помощью уравнения
dxt a( xt , t )dt ( x, t )dwt на интервале t [0, T ] с начальным условием xt 0 v . Здесь
wt R m стандартный винеровский процесс. Данное уравнение может быть представлено в
скалярном виде: dxk , t (t ) ak ( xt , t )dt
обозначения | a( x, t ) |
1 m
ki ( x(t ), t )dwi,t , xk ,t 0 vk , k 1,2,..., n . Введем
2 i 1
n
n m
k 1
k 1i 1
ak2 ( x, t ) и || ( x, t ) || 2 ki2 ( x, t ) . Рассмотрим условия существования и
единственности решения приведенной стохастической системы.
Теорема / 3 /.
Пусть случайная величина v не зависит от {wt , t [0, T ]}, M {| v |2 } , а коэффициенты
уравнения dxt a( xt , t )dt ( x, t )dwt непрерывны по переменным t [0, T ], x R n . Пусть также:
- найдется такая величина K , что при всех t [0, T ], x R n будет выполняться соотношение
| a( x, t ) |2 || ( x, t ) || 2 K (1 | x |2 ) ;
- найдется такая величина C , что при всех t [0, T ], x, y R n будет выполняться соотношение
| (a( x, t ) a( y, t ) | 2 || ( x, t ) ( y, t ) || 2 C || x y || 2 .
Тогда на интервале t [0, T ] существует и единственно ( почти наверное) непрерывное решение
xt уравнения dxt a( xt , t )dt ( x, t )dwt , xt 0 v , причем M {| xt |2 } L(1 M {| v |2 }) , где константа
L зависит лишь от T , K .
Рассмотрим теперь следующее стохастическое дифференциальное уравнение:
dt A(t )t dt u(t )dt B(t )dwt ;t 0 v , где A(t ) R nn , B(t ) R n q , u(t ) R n непрерывные
неслучайные функции.
Можно показать, что процесс, определяемый равенством:
t
t
t0
t0
t (t , t0 )v (t , )u ( )d (t , ) B( )dw будет являться решением стохастического
дифференциального уравнения, если функция (t , t0 ) будет удовлетворять следующей системе:
(t , t0 ) A(t ) (t , t0 ), t t0
.
(t0 , t0 ) E
Второй интеграл в равенстве, определяющем решение, понимается в смысле Ито.
Уравнения моментов для линейных стохастических систем.
Рассмотрим линейную стохастическую систему следующего вида dxt A(t )dt u(t )dt B(t )dwt , где
xt 0 v , wt - случайный процесс с независимыми приращениям, mw 0; Rw (t , t ' ) Q (t t ' ) ; а
детерминированные матричные функции A(t ), u(t ), B(t ) являются непрерывными . Случайная
77
величина v является независимой от процесса wt . Решением этого уравнения, как уже
t
t
t0
t0
отмечалось, является процесс t (t , t 0 )v (t , )u ( )d (t , ) B( )dw , если матричная
(t , t0 ) A(t ) (t , t0 ), t t0
.
(t0 , t0 ) E
Пусть процесс xt в момент t [0, T ] имеет среднее значение mx (t ) M {x(t )} и дисперсионную
функция (t , t0 ) , удовлетворяет системе уравнений:
матрицу Dx (t ) M {( x(t ) mx (t ))( x(t ) mx (t ))T } , а начальное значение xt 0 v не зависит от wt .
Тогда функции mx (t ), Dx (t ) будут являться решениями следующих систем обыкновенных
x (t ) A(t )mx (t ) u(t ) , mx (0) mv ;
дифференциальных уравнений: m
D x (t ) A(t ) Dx (t ) Dx (t ) AT (t ) B(t )QBT (t ) , Dx (0) Dv , где Dv M {(v M {v})(v M {v})T } дисперсионная матрица вектора x(0) v . Данные уравнения называются уравнениями метода
моментов.
Пример.
Пусть скалярная случайная функция
xt удовлетворяет уравнению dxt adt bdwt , xt 0 v, a 0 .
t
Решением этого уравнения будет процесс
xt e at mv be at e a w( )d . Уравнения моментов будут иметь
вид:
x (t ) amx (t ) , mx (0) mv и D x (t ) 2aDx (t ) b , Dx (0) v2 . Отсюда получим следующие
m
2
соотношения
mx (t ) mv e at , Dx (t ) e 2at v2
mx (t ) 0, Dx (t )
b2
(1 e 2at ) . Следовательно, при t получим
2|a|
b2
2
для любых mV , Dv v .
2|a|
Различные формы представления стохастических систем.
Во многих, практически важных приложениях, стохастические системы описываются уравнением
вида xt a( xt , t ) ( xt , t )t , xt 0 v . Здесь xt - производная в среднеквадратическом (с.к.) смысле
случайного процесса xt R n ; t R m - с.к. – непрерывная при t 0 случайная функция;
a( x, t ), ( x, t ) - непрерывные неслучайные функции; v - случайная величина. Данная запись
является, по сути, записью следующего соотношения:
t
t
xt v a( x( ), )d ( x( ), ) d , где все интегралы понимаются в с.к. – смысле.
Если полоса пропускания системы значительно меньше полосы частот равномерности
спектральной плотности стационарного возмущающего процесса t , то этот процесс может
считаться белым шумом по отношению к данной системе, то есть имеющим постоянную
спектральную плотность в диапазоне частот существенном для динамики системы.
78
Таким образом, понятие белошумности случайного процесса достаточно условно. Один и тот же
случайный процесс может считаться белошумным по отношению к одной динамической системе и
«цветным» по отношению к другой, имеющей более широкую полосу пропускания. Для
формализации этого утверждения, была предложена параметризованная модель белошумного
случайного процесса. Если ввести малый параметр 0 и построить процесс t
1
t
, где
2
t R m - стационарный гауссов случайный процесс с корреляционной функцией R ( ) , такой, что
R (t , t )d E , и спектральной плотностью S ( ) , то процесс t
будет иметь корреляционную
функцию R ( ) и спектральную плотность S ( ) , такие, что при 0 будут выполняться
соотношения
R ( )
1
2
R (
) ( ) E , S () S ( 2) S (0) .
2
Таким образом, стохастическая система, для параметризованных случайных возмущающих
процессов, принимает вид xt a( x, t ) ( x, t )
процесс t
1
, x(t0 ) v . При достаточно малом значении
t
2
может считаться белошумным по отношению к данной динамической системе,
t
1
2
так как полоса равномерности спектра возмущающего процесса
1
t
неограниченно расширяется
2
при 0 , а полоса пропускания системы остается постоянной.
Если полоса пропускания системы не достаточно узкая по отношению к полосе равномерности
спектра возмущающего процесса (t ) , то белошумная аппроксимация этого процесса будет не
адекватна постановке задачи стохастической динамики. В этих случаях, необходим учет свойств
«цветности» процесса t . Для этого спектральную плотность процесса подвергают процедуре
факторизации, то есть представляют в следующем виде: S ( ) (i )(i )
B(i ) B(i )
, где
A(i ) A(i )
B(s), A(s) - многочлены конечных степеней порядка nB , n A ; nB n A . Причем многочлен A(s)
является гурвицевым. Отсюда можно построить уравнения соответствующего формирующего
фильтра: zt Dzt Gnt ; t C T zt , где nt - процесс типа белого шума; ( s) C T ( sE D) 1 G
B( s )
.
A( s)
79
Таким образом, достаточно универсальной математической моделью, описывающей стохастическую
динамику разнообразных систем, будет следующая: xt a( xt , t ) ( xt , t )t , x(t0 ) v .
Однако, необходимо отметить, что случайные процессы типа белого шума, имеют бесконечную
дисперсию. То есть, говорить о с.к. – смысле полученного соотношения надо с большой
осторожностью, так как все математические преобразования придется выполнять с обобщенными
функциями. Поэтому такую систему обычно приводят к форме стохастического уравнения Ито при
проведении аналитических вычислений, то есть к форме dxt a( xt , t ) ( xt , t )dwt , x(t0 ) v .
Решение такого уравнения представляется в виде соотношения
t
t
xt (t ) v a( x , )d ( x , )dw , где первый интеграл в правой части понимается в с.к. –
смысле, а второй – является интегралом Ито.
Восстановление координат состояния линейной стационарной модели объекта при
воздействии случайных шумов.
Пусть задана управляемая система
x Ax Bu n0 (t ); x(t0 ) x0 ; y Cx nн (t ) , где
матричные функции
A, B, C , в общем случае, являются детерминированными функциями времени,
x0 - случайная величина, а n0 (t ), nн (t ) - белые шумы со следующими характеристиками:
M {x0 } x0 ; M {( x0 x0 )( x0 x0 )T } P0 ; M {n0 (t )} 0; M {n0 (t )n0T (t )] Q0 (t ) (t t ) ;
M {nн (t )} 0; M {nн (t )nнT (t )] R0 (t ) (t t ) ; M {n0 (t )nнT (t )] W0 (t ) (t t ) . Здесь Q0 , P0
- положительно полуопределенные матрицы; R0 (t ) - положительно определенная матрица. При
этом случайная величина x0 не коррелированна с шумами n0 (t ), nн (t ) .
Требуется на основе наблюдения выхода y (t ) на интервале [t 0 , t ] определить несмещенную
оценку xˆ (t ) на выходе линейного фильтра, обеспечивающую минимум среднего квадрата ошибки:
J M [( x(t ) xˆ (t ))T ( x(t ) xˆ (t ))] (эффективная оценка).
Условие положительной
определенности матрицы R0 (t ) означает, что ни одна компонента выхода y (t ) не измеряется
точной. То есть задача оценивания является несингулярной (невырожденной).
Теорема . Если шумы объекта и наблюдения не коррелированны ( W0 (t ) 0 ), то оценка
xˆ (t )
является несмещенной и эффективной в том случае, когда она находится из уравнения
xˆ Axˆ Bu K 0 ( y Cxˆ ); xˆ (t0 ) x (t0 ) с матрицей коэффициентов усиления K 0 PCT R01 ,
AP PAT PCT R 1CP Q ; P(t ) P .
где P определяется из уравнения P
Матрица P является дисперсионной матрицей ошибки
e x(t ) xˆ (t ) , то есть
P M [e(t )eT (t )], а уравнение, которому удовлетворяет данная матрица, называется
дисперсионным уравнением.
Теорема . Если шумы объекта и наблюдения коррелированны ( W0 (t ) 0 ), то оценка
xˆ (t )
является несмещенной и эффективной в том случае, когда она находится из уравнения
xˆ Axˆ Bu K 0 ( y Cxˆ ); xˆ (t0 ) x (t0 )
с матрицей коэффициентов усиления
80
K 0 ( PCT W0 ) R01 , где P определяется из уравнения
P ( A W0 R01C ) P P( A W0 R01C )T PCT R01CP Q0 W0 R01W0T ; P(t0 ) P0 .
Полученные фильтры называются фильтрами (наблюдателями) Калмана-Бьюси. Фильтры
Калмана-Бьюси имеют такую же структуру, что и наблюдатели (идентификаторы) полного порядка в
детерминированном случае. Отличие состоит в том, что когда случайные воздействия не
учитываются (то есть в детерминированном случае), матрица коэффициентов усиления выбирается
достаточно произвольно. Когда случайные воздействия учитываются, эта матрица определяется
однозначно.
Пример.
Определить оптимальную оценку постоянной величины
x
M {x} x и дисперсией
- белый шум с интенсивностью r0 , математическим ожиданием
nн (t )
y x nн (t ) , где
по наблюдениям за выходом
M {( x x ) 2 } p0 . Очевидно, что уравнение объекта имеет вид x 0 , то есть
A 0; B 0; C 1; R0 r; Q0 0 . Тогда фильтр Калмана-Бьюси будет описываться уравнением
xˆ K ( y xˆ ); xˆ (t ) x , где K p / r . Параметр p будет определяться из уравнения
p p / r0 ; p(t 0 ) p0
2
. Отсюда коэффициент усиления получаем
k0 p0 /(r0 p0t ) .
Фильтр Калмана-Бьюси при цветном шуме объекта.
Рассмотрим теперь задачу линейного оптимального оценивания при условии, что шум
объекта является цветным, то есть не белым. В этом случае задача синтеза оптимального фильтра
сводится к следующей постановке. Управляемая система описывается уравнениями:
x (1) Ax(1) B1u x (2) (t ); x (1) (t0 ) x0(1) ; x (2) A2 x ( 2) n~0 (t ); x ( 2) (t0 ) x0( 2) ; y C1 x (1) nн (t ) ,
(1) ( 2)
~ (t ), n (t ) - белые шумы, с вероятностными
где x , x
- случайные величины, n
н
характеристиками:
M {x0(1) } x0(1) ; M {( x0(1) x0(1) )(x0(1) x0(1) )T } P01 ;
M {x0(1) } x0(1) ; M {( x0(1) x0(1) )(x0(1) x0(1) )T } P01 ;
~
M {n~0 (t )} 0; M {n~0 (t )n~0T (t )] Q0 (t ) (t t ) ; M {nн (t )} 0; M {nн (t )nнT (t )] R0 (t ) (t t ) ;
~
M {n~ (t )nT (t )] W (t ) (t t ) . Здесь Q , P - положительно полуопределенные матрицы;
н
R0 (t )
- положительно определенная матрица. При этом случайные величины
коррелированны с шумами
интервале
[t 0 , t ]
x0(1) , x0( 2)
n~0 (t ), nн (t ) . Требуется на основе наблюдения выхода y (t )
определить эффективную несмещенную оценку
xˆ (1) (t )
не
на
на выходе линейного
фильтра, обеспечивающую минимум среднего квадрата ошибки:
J M [(x (1) (t ) xˆ (1) (t ))T ( x (1) (t ) xˆ (1) (t ))] . Преобразуем данную задачу в задачу фильтрации с
белыми шумами. Введем следующие переменные: x x , x x 0 ; x x 0 ;
x ( 2) 0 x ( 2) 0 x ( 2)
0
0
A
A 1
0
0
P10
E ;
B
B 1 ; C [C1 ,0] ; G ; P0
A2
E
0
0
(1)
(1)
(1)
~ T
0
; Q0 GQG . Тогда задачу можно
P20
81
x Ax Bu n0 (t ); x(t0 ) x0 ; y Cx nн (t ) ,
J M [( x(t ) xˆ (t ))T ( x(t ) xˆ (t ))] min , где n0 (t ) Gn~0 (t ) . При этом фазовый вектор и шумы
переформулировать следующим образом
не коррелированны между собой и имеют вероятностные характеристики
M {x0 } x0 ; M {( x0 x0 )( x0 x0 )T } P0 ; M {n0 (t )} 0; M {n0 (t )n0T (t )] Q0 (t ) (t t ) ;
M {nн (t )} 0; M {nн (t )nнT (t )] R0 (t ) (t t ) ; M {n0 (t )nнT (t )] W0 (t ) (t t ) . Отсюда
T 1
найдем уравнения для оптимальной оценки: xˆ Axˆ Bu K ( y Cxˆ ); xˆ (t ) x (t ) , K PC R ,
P AP PAT PCT R01CP Q0 ; P(t 0 ) P0
Структурная схема фильтра Калмана-Бьюси, при этом, помимо модели объекта будет
включать еще и модель формирователя цветного шума.
Фильтр Калмана-Бьюси при цветном шуме наблюдения.
Рассмотрим теперь задачу линейной оптимальной фильтрации при условии, что шум объекта
является белым, а шум наблюдения – цветным. В этом случае управляемая система описывается
уравнениями
~
x Ax Bu n0 (t ); x(t0 ) x0 ; ~
y Cx z; z Dz n~н (t ) , где
M {x0 } x0 ; M {( x0 x0 )( x0 x0 )T } P0 ; M {n0 (t )} 0; M {n0 (t )n0T (t )] Q0 (t ) (t t ) ;
~
~
M {n~н (t )} 0; M {n~н (t )n~нT (t )] R0 (t ) (t t ) ; M {n0 (t )n~нT (t )] W (t ) (t t ) .
~
Здесь Q0 , P0 - положительно полуопределенные матрицы; R0 (t ) - положительно
определенная матрица. При этом случайная величина x0 не коррелированна с шумами
n0 (t ), n~н (t ) .
Требуется на основе наблюдения выхода
y (t )
на интервале
эффективную несмещенную оценку
xˆ (t )
минимум среднего квадрата ошибки:
J M [( x(t ) xˆ (t ))T ( x(t ) xˆ (t ))] .
[t 0 , t ]
определить
на выходе линейного фильтра, обеспечивающую
Преобразуем данную
задачу в задачу фильтрации с белыми шумами. Продифференцировав уравнение выхода по
~ ~
~
~
~
y (C CA) x CBu Dx Cn0 n~н .
~
y~
y CBu D~
y . Тогда после преобразований получим
времени, можно записать следующее соотношение
Введем новый вектор наблюдения
~
~ ~
~
y Cx nн , где C C CA DC ; nн Cn0 n~н . Шум nн имеет интенсивность
~ ~
~~ ~
~ ~
R0 CQ0C T W T C T CW R0 . Интенсивность взаимной корреляционной функции шума nн и
~T ~
шума объекта определяется соотношением W0 Q0C W . Очевидно, что шум объекта n0 (t ) и
шум nн (t ) преобразованного уравнения наблюдения будут коррелированны. Если матрица R0 (t )
положительно определена, то оптимальный фильтр будет описываться уравнениями
82
xˆ Axˆ Bu K 0 ( y Cxˆ ); xˆ (t0 ) x (t0 ) , K 0 ( PC T W0 ) R01 ,
P ( A W0 R01C ) P P( A W0 R01C )T PCT R01CP Q0 W0 R01W0T ; P(t 0 ) P0 .
Однако в это соотношение входит производная выхода ~
y , что нежелательно. Поэтому
введем следующий вектор ~
xK ~
y.
x , который определим следующим равенством xˆ ~
Продифференцировав это равенство по времени получим
~
~
x ( A K 0C ) xˆ ( B K 0CB)u ( K 0 K 0 D) ~
y или
~
~
x (t0 ) x0 K 0 (t0 ) ~
y (t0 ) . В
x ( A K 0C ) ~
x ( B K 0CB)u [( A K 0C ) K 0 K 0 K 0 D]~
y; ~
это уравнение производная ~
y не входит.
Это позволяет определить оценку ~
x без
x K0 ~
y найти искомую
дифференцирования выхода объекта, а затем с помощью соотношения xˆ ~
оценку
x̂ .
Пример.
x 0; ~y x z; M [ x(0)] 0; M [ x 2 (0)] p0 . Шум наблюдения является
1 | |
стационарным случайным процессом с характеристиками M [ z ] 0; M {z (t ) z (t )] e
. Тогда
2
~ , где M [n~ ] 0; M [n
~ (t )n
~ (t )] (t t ) . Таким
уравнение формирователя имеет вид z z n
н
н
н
н
Объект описывается уравнениями
образом, в данном случае, будут справедливы соотношения
~
~
~
A 0; B 0; C 1; Q0 0; D 1; R0 1;W 0
C 1; R0 1;W0 0 . Тогда
x k0 ~
x (k02 k0 k0 ) ~
y;
можно записать следующие уравнения оценивания ~
~
~
~
~
x (0) k0 y (0); xˆ x k0 y . Параметры k 0 , p определяются из результатов теоремы ( W0 0 ) с
помощью уравнений
. Из этого следует, что
k0 p; p p 2 ; p(0) p0 . Отсюда легко вычислить p k 0
p0
.
1 p0t
Линейные стохастические оптимальные системы при неполной информации о
состоянии.
Пусть уравнение вполне управляемого объекта описываются линейными уравнениями со
случайной правой частью вида
x Ax Bu 0 (t ); y Cx н (t ); x(t0 ) x 0 ; x R n ; u R m , а критерий оптимальности
имеет вид:
tf
J M {x (t f )x(t f ) [ x T (t ) Lx(t ) u T (t )Gu (t )]dt} ,
T
где
0 (t ) - гауссов белый шум
t0
M { 0 (t )} 0; M { 0 (t ) 0 (t )} Q0 (t t ) , Q0 - симметричная матрица; н (t ) гауссов белый шум на выходе (шум наблюдения), M { н (t )} 0; M { н (t ) н (t )} Rн (t t ) ,
на входе объекта,
(принимается, что процессы
0 (t ), н (t ) некоррелированные); x 0 - гауссова случайная величина
M {x 0 } x 0 ; M {( x 0 x 0 )(x 0 x 0 )T } H 0 , некоррелированная со случайными процессами
0 (t ), н (t ) . Матрицы A, B, L, G являются, в общем случае, функциями времени. Требуется
найти управление u u ( y( ), t 0 t ); t 0 t t f , при котором критерий оптимальности
принимает минимальное значение.
83
Теорема (принцип разделения). Стохастическое оптимальное управление с обратной связью
для объекта
x Ax Bu 0 (t ); y Cx н (t ); x(t0 ) x 0
при критерии оптимальности
tf
J M {x (t f )x(t f ) [ x T (t ) Lx(t ) u T (t )Gu (t )]dt} имеет вид u G 1BT Kxˆ , где K
T
-
t0
симметричная матрица, которая определяется из матричного уравнения Риккати
K KA AT K KBG 1BT K L
при граничном условии
K (t f ) ; а x̂
- оптимальная
оценка, которая определяется с помощью фильтра Калмана – Бьюси:
xˆ Axˆ Bu K 0 ( y Cxˆ ); xˆ (t0 ) x 0 ) с матрицей коэффициентов усиления K 0 PCT Rн1 ,
AP PAT PCT R 1CP Q ; P(t ) H .
где P определяется из уравнения P
н
o
Таким образом, в законе управления используется не сам фазовый вектор, а его оценка, которая
получается на выходе фильтра Калмана – Бьюси. То есть стохастический оптимальный регулятор
состоит из оптимального фильтра и детерминированного оптимального регулятора.
Этот результат носит название принципа разделения или принципа стохастической
эквивалентности. В соответствии с этим принципом задача синтеза стохастической оптимальной
системы управления при неполной информации сводится к двум задачам, решаемым раздельно.
Первая задача сводится к задаче синтеза фильтра Калмана – Бьюси. Вторая задача сводится к
синтезу детерминированной оптимальной системы. Если шумы и начальное состояние подчиняются
гауссову закону распределения, то в результате такого синтеза получим оптимальную
стохастическую систему.
Методы интегрирования линейных стохастических уравнений.
Большинство численных методов для обыкновенных дифференциальных уравнений строятся
с помощью разложения в ряд Тейлора. Поэтому попробуем применить такой же подход для
решения стохастических уравнений.
Рассмотрим скалярное стохастическое дифференциальное уравнение вида:
xt a( xt , t ) b( xt , t ) nt ,
где
xt
- решение уравнения;
белого шума,
a, b : [0, T ] ; nt - случайный процесс в виде гауссова
nt dt dwt .
Запишем это же уравнение в форме стохастического дифференциального уравнения Ито:
dxt a( xt , t ) b( xt , t ) dwt ,
где wt - стандартный винеровский процесс.
84
Применим обычную идеологию численных методов к данным стохастическим уравнениям.
Согласно формуле Тейлора при значении s t имеем:
dxt
d 2 xt ( s t )2
xs xt
(s t ) 2
O(( s t )3 ) ,
dt
dt
2
где запись
t O(s t ) )
понимается как:
lim
s t
M {t2 }
( s t )
C const . Используя первое
уравнение, найдем:
d
( s t )2
xs xt a( xt , t ) ( s t ) b( xt , t ) ( s t ) nt [a( xt , t ) b( xt , t ) nt ]
O((s t )3 ) .
dt
2
Правая часть данного соотношения неопределенна. Это связано с тем, что белый шум nt не
является дифференцируемым в каком – либо смысле случайным процессом.
Рассмотрим теперь применение этой же идеологии в отношении уравнения Ито. Используем
для случайных функций
такие же правила дифференцирования, как и для
a( xt , t ), b( xt , t )
неслучайных функций, то есть:
da( xt , t )
a( xt , t )
a( xt , t )
b( xt , t )
b( xt , t )
dxt
dt , db( xt , t )
dxt
dt ,
x
t
x
t
как это делается в обычных алгоритмах численного решения. Запишем решение уравнения в
интегральной форме:
s
s
t
t
xs xt a( x, ) d b( x , ) dw ,
s
a( x , ) a( x , )
a( x , )
a( xs , s) a( xt , t ) [a( x , )
] d b( x , )
dw
x
x
t
t
s
где:
s
b( x , ) b( x , )
b( x , )
b( xs , s ) b( xt , t ) [a ( x , )
] d b (x , )
dw
x
x
t
t
,
s
.
Подставляя данные выражения в разложение Тейлора, получим следующее разложение:
b( xt , t ) s
xs xt a ( xt , s ) ( s t ) b( xt , t ) dw b( xt , t )
dw1 dw
x
t
t t
a( xt , t ) s
b( xt , t )
b( xt , t ) s
b( xt , t )
dw1 d [a ( xt , t )
dxt
] d 1 dw
x
x
t
t t
t t
s 1
2
b( xt , t ) 2
b( xt , t )
b( xt , t ) [(
) b( xt , t )
] dw2 dw1 dw O(( s t ) 2 )
2
x
x
t t t
s
Однако, поскольку для случайных процессов
a( xt , t ), b( xt , t ) , правила дифференцирования
должны выполняться по формуле Ито, то с вероятностью 1 получим:
s
a( xs , s) a( xt , t ) [a( x , )
t
a( x , )
b( x , )
dw
x
t
s
a( x , ) a( x , ) 1 2
2 a( x , )
b ( x , )
] d
x
2
x 2
85
s
b( xs , s) b( xt , t ) [a( x , )
t
b( x , ) b( x , ) 1 2
2b( x , )
b ( x , )
] d
x
2
x 2
b( x , )
b( x , )
dw
x
t
s
Подставляя полученные выражения в соотношение:
s
s
t
t
xs xt a( x, ) d b( x , ) dw ,
получим с вероятностью 1 разложение:
b( xt t ) s
xs xt a( xt , s) ( s t ) b( xt , t ) dw b( xt , t )
dw1 dw
x
t
t t
s
a( xt , t ) s
b( xt , t )
dw1 d
x
t t
b( xt , t )
b( xt , t ) 1 2
2b( xt , t ) s s
.
[a( xt , t )
dxt
b ( xt , t )
] d 1 dw
x
t
2
x 2
t t
b( xt , t ) 2
2b( xt , t ) s 1
b( xt , t ) [(
) b( xt , t )
] dw2 dw1 dw O(( s t ) 2 )
2
x
x
t t t
Сопоставляя разложения при разных правилах дифференцирования, найдем что они
2b( xt , t )
1 2
b ( xt , t )
] d 1 dw . Это слагаемое имеет порядок
2
x 2
t t
s s
различаются на слагаемое
3
.
2
Таким образом, приходим к выводу, что при применении обычных правил дифференцирования
будут отсутствовать часть слагаемых. То есть формальное применение правил
дифференцирования для неслучайных функций к случайным процессам приводит, либо к
неопределенности, либо к неверному результату. Поэтому применение сильных и слабых
численных методов интегрирования на основе детерминированного разложения Тейлора является
малоэффективным.
При сильных методах нахождения решения предполагается, что реализации
аппроксимирующего процесса
xˆt
были близки к траекториям исходного процесса
xt . В этом случае
M {| xt xˆt |} .
Говорят, что аппроксимирующий процесс xˆt сходится в сильном смысле с порядком
критерием близости является абсолютная ошибка
[0, )
к решению
константа
0 0 ,
xt , если,
существуют конечная константа
такие, что
максимальным размером шага
M {| xt xˆt |} R
R
и положительная
для любой временной дискретизации с
(0, 0 ) .
Однако, во многих случаях, нужна не близость аппроксимирующей и истинной траекторий, а
только близость некоторой функции (и/или момента) от значений случайного сигнала. Тогда
методы, реализующие такое приближение, называются слабыми.
86
Говорят, что аппроксимирующий процесс
[0, ) , если для любого полинома g
константа
0 0 ,
такие, что
xˆt
сходится в слабом смысле с порядком
существуют конечная константа
| M {g ( xt )} M {g ( xˆt )}|} R
дискретизации с максимальным размером шага
R
и положительная
для любой временной
(0, 0 ) .
Рассмотрим метод интегрирования линейных стохастических уравнений на основе
численного вычисления интеграла Ито.
Интегрирование линейных стохастических дифференциальных уравнений.
Рассмотрим линейную систему, описываемую стохастическим дифференциальным
уравнением вида:
dxt A(t ) xt dt B(t ) f (t ) dt G(t ) dwt ; xt t0 v ,
где
xt n
условием
- случайный процесс, являющийся решением уравнения с начальным случайным
xt t0 v ; wt m
- винеровский случайный процесс с независимыми компонентами
wt (w1,t ,..., wm,t )T , wk ,t , M {wt } 0; M {wt wtT } 2w diag{2wk }km1 ;
f (t ) p - детерминированное возмущение;
A(t ), B(t ), G(t ) - соответствующие матричные функции, удовлетворяющие условиям
существования и единственности решения стохастического уравнения.
Приведенную систему можно рассматривать, как систему уравнений Ито или Стратоновича,
поскольку можно показать, что в силу своей структуры оно имеет одно и то же решение, в каком бы
смысле она не рассматривалась. Поэтому, часто применяемая в технической литературе форма
описания линейной стохастической системы в виде:
xt A(t ) xt B(t ) f (t ) G(t ) w nt ; xt t0 v ,
с белым шумом
nt с единичной дисперсионной матрицей в правой части не приводит к
неоднозначности математического определения решений системы.
Рассмотрим стационарный случай системы, то есть когда ее модель описывается
стохастическим уравнением:
dxt A xt dt B f (t ) dt G dwt ; xt t0 v .
В этом случае решение будет определяться с вероятностью 1 в виде:
xt e
A( t t0 )
t
v e
A( t s )
t
B f ( s) ds e A(t s ) G dws .
t0
t0
Для численного решения выберем регулярную сетку
tk 1 tk h; t0 0; k 0,1, 2,... , где h h
шаг сетки. Введем следующие обозначения
ˆ ) e A( h s ) B f (kh s) ds ,
Aˆ e Ah , ( Bf
k
h
nˆk 1 e A( h s ) Gdw( kh s ) . Тогда непрерывную систему можно аппроксимировать дискретной
моделью вида:
87
xk 1 Aˆ xk ( Bˆ f )k nˆk 1 ,
где случайные величины
ˆ ) , nˆ сохраняют свое значение на интервале h .
xk ,( Bf
k
k 1
Векторный гауссов дискретный случайный процесс n̂k имеет дисперсионную матрицу:
h
D{nˆk 1} Dnˆ (h) e A( h s )G 2w GT e A
T
(hs)
ds .
Несложно показать, что дисперсионная матрица
Dnˆ (t ), t [0, h]
является решением
уравнения:
Dnˆ (t ) A Dnˆ (t ) Dnˆ (t ) AT G 2w GT ; Dnˆ (0) 0 .
То есть, для моделирования случайной составляющей решения надо на каждом шаге
моделировать гауссов случайный процесс с нулевым средним и дисперсионной матрицей
Dnˆ (h) .
Если шаг является переменным, то надо вычислять, при каждом его изменении, новую
дисперсионную матрицу, что является вычислительно трудоемким.
Вычисление систематической составляющей решения линейной стохастической модели.
Оценка систематической составляющей сводится к вычислению интеграла:
h
Z h e A( h s ) B f (kh s ) ds .
Представим векторную функцию
f (kh s), s [0, h]
с точностью
следующей
полиномиальной функцией:
L
f (kh s) q j (k ) s j r (k , s) ,
j 0
где
| r (k , s) | s L 1 для величин k 0,1, 2,... ; а векторные коэффициенты q j (k ), j 1, 2,..., L
могут быть различными способами. Тогда можно записать:
L
Z h Z h, j B q j (k ) Rh, L 1 ,
j 0
h
где
|| Z h, L 1 |||| e
A( h s )
h
Br (k , s) ds || || Z h, L 1 || || B || , Z h, j s j e A( h s ) ds .
Для вычисления интегралов
Z h, j
можно использовать следующие очевидные рекуррентные
соотношения:
Z h,0 A1 (e Ah E ), Z h, j A1 (h j E j Z h, j 1 ); j 1, 2,..., L .
Отсюда найдем:
( j 1)
Z h, j A
или, отбрасывая остаточный член
j ! ( A Z h,0
hl Al
); j 1
l!
l 1
j
Rh, L 1 , получим следующее приближенное представление:
88
L
j
j 0
l 0
Z h, j A ( j 1) j ! (e Ah
При кусочно – постоянной аппроксимации
( L 0)
hl Al
).
l!
выполняется соотношение
q0 (k ) f (kh) f k . Тогда:
Z h A1 (e Ah E ) B f k Bˆ f k .
То есть:
( Bˆ f ) k Bˆ f k ,
где
Bˆ A1 (e Ah E ) B .
При линейной интерполяции ( L 1 ) получим:
q0 (k ) f k , q1 (k )
1
( f k 1 f k ) .
h
Отсюда найдем:
1
Z h A (e
Ah
A2
E) B fk
(e Ah E h A) B ( f k 1 f k )
h
Вычисление случайной составляющей решения линейной стохастической модели.
Рассмотрим дисперсионную матрицу:
h
Dnˆ (h) M {nˆk 1 nˆkT1} e A( h s ) G 2w GT e A
T
(hs)
ds .
Так как эта матрица является квадратной, симметричной и невырожденной, то ее SVD
разложение можно записать в виде:
Dnˆ (h) M {nˆk 1 nˆkT1} T 2D T 1 ,
где матрица
T
D diag (1 , 2 ,.., n ) , k , k 1,2,..., n - собственные значения матрицы Dnˆ (h) ;
- матрица ортонормированных собственных векторов матрицы
матрица
T
является ортогональной. Представим вектор
случайный вектор – столбец, такой, что его компоненты
nˆk 1
Dnˆ (h) . Очевидно, что, так как
в виде nˆk 1 T D nk 1 , где nk 1 -
n j ,k 1 , j 1, 2,.., n
являются случайными
величинами с нулевым средним с единичной дисперсией. Поэтому, чтобы смоделировать
случайный вектор
nˆk 1
надо найти матрицу
векторы, а также смоделировать
n
Dnˆ (h) , ее собственные значения и собственные
независимых гауссовых случайных величин с нулевым средним
и единичной дисперсией, для построения вектора
nˆk 1 .
Приближенный метод интегрирования линейных стохастических дифференциальных уравнений.
Вычисление дисперсионной матрицы
Dnˆ (h) является алгоритмически достаточно
трудоемким. Поэтому достаточно часто используют приближенный метод численного решения
стохастических уравнений. Для этого белый шум
[kh,(k 1) h]; k 0,1, 2,...
nt ,(dwt w nt dt )
кусочно постоянным случайным процессом
заменяется на промежутке
nt
с достаточно малым
89
шагом дискретизации . Запишем следующее приближенное представление стохастической
системы в следующей форме:
xt A xt B u (t ) G nt ; yt C xt ; xt t0 v ,
где
nt - кусочно постоянный на промежутках [kh,(k 1) h]; k 0,1, 2,...;
h
; N 1 , случайный
N
процесс с независимыми значениями на этих промежутках.
То есть имеют место соотношения:
nt zk ,r const
при значении
t [r,(r 1) ) ;
1
M {zk ,r } 0 ; M {zk ,r zkT,r } w .
Корреляционная функция случайного процесса
nt
имеет вид:
1
|t |
(1 ) 2w ,| t |
Rn (t )
.
0,| t |
0 корреляционная функция процесса nt
Таким образом, при
будет стремиться к
-
функции Дирака, а спектральная плотность к постоянной величине. Если выполняется условие
2
, где max - максимальная частота, наблюдаемая в системе, то будет справедливо
8
2
2
соотношение 0.95 w Sn () w . Таким образом, на интервале длиной h должно быть не
max
менее 8 подинтервалов
постоянства процесса
nt .
Построим дискретное представление линейной стохастической модели на регулярной сетке:
[kh,(k 1) h]; k 0,1, 2,...;
в виде:
h
xk 1 e xk e
Ah
A( h s )
h
B u (kh s) ds e A( h s ) G dn( kh s ) .
Оценим стохастический интеграл, используя кусочно постоянное приближение процесса nt
на более мелкой сетке с шагом
h
e
:
A( h s )
N
G dn( kh s )
1
A (e
A
e A( h s ) ds G n( kN j 1)
j 1 ( j 1)
N
A( h j )
1
j
E) e
j 1
Представляя дисперсионную матрицу
.
G zk , j
w T w T 1 в виде SVD представления,
получим:
h
A( h s )
e G dn(kh s )
где
1
N
A1 (e A E ) e A( h j) G zk , j
,
j 1
G G T w , zk , j - вектор-столбец независимых стандартных гауссовых случайных
величин.
90
Представляя систематическую составляющую с помощью кусочно – постоянной
аппроксимации, найдем, что:
h
e
A( h s )
B u (kh s) ds A1 (e Ah E ) Buk .
Тогда можно записать следующую рекуррентную формулу для приближенного
моделирования линейной стохастической системы:
xk 1 e Ah xk A1 (e Ah E ) B uk
1
N
A1 (e A E ) e A( h j) G zk , j ; yk C xk ,
j 1
где N 8 .
Алгоритм моделирования можно записать в следующей форме.
1.
Определение временного интервала моделирования
сеток
h
2.
Вычисление с заданной точностью матриц:
[0, t max] и задание шагов размеров
tmax
h
, k 1, 2,..., M ; , j 1, 2,..., N ; N 8 .
M
N
e A , e Ah , A1 (e Ah E ) B,
1
A1 (e A E ); G .
3.
k 1; x1 xt 0 v ; u1 u(t 0) .
4.
Численное моделирование N случайных независимых гауссовых векторов
zk , j ; j 1, 2,..., N размерности m .
5.
Вычисление с заданной точностью матриц
6.
Вычисление
xk 1 e Ah xk A1 (e Ah E ) B uk
1
e A( h j) G, j 1, 2,..., N .
N
A1 (e A E ) e A( h j) G zk , j ; yk C xk
j 1
Если k M , то k k 1 и переход к шагу 4.
Конец работы алгоритма
7.
8.
Расчет моментных характеристик решений линейных стохастических уравнений.
Пусть задана линейная стохастическая система:
dxt A xt dt B f (t ) dt G dwt ; yt C xt ; xt t0 v ,
где
xt n
условием
v
- случайный процесс, являющийся решением уравнения с начальным случайным
xt t0 v ; yt q
- выход системы;
- случайная гауссова величина с законом распределения
wt m
N (mv , Dv ) ;
- винеровский случайный процесс с независимыми компонентами
wt (w1,t ,..., wm,t )T , wk ,t , M {wt } 0; M {wt wtT } 2w diag{2wk }km1 ;
f (t ) p
A, B, G, C
- детерминированное возмущение;
- постоянные вещественные матрицы соответствующего размера.
91
Введем обозначения:
mx (t ) M {xt }; xt0 xt mx (t ) , Rx (t , s) M {xt0 xs0 }; Dx (t ) Rx (t , t ) .
Используя уравнения моментов можно записать следующие соотношения:
mx (t ) A mx (t ) B f (t ); mx (0) mv ;
Dx (t ) A Dx (t ) Dx (t ) AT G 2w GT ; Dx (0) Dv .
Если
f (t ) f const
и матрица
A
является гурвицевой, то приведенные соотношения
имеют стационарное решение, определяемое уравнениями:
A mx () B f 0 ;
A Dx () Dx () AT G 2w GT 0 ;
e A(t s ) Dx (), t s
Rx (t , s) Rx (t s)
AT ( s t )
,t s
Dx ()e
Таким образом, процесс
xt
является асимптотически стационарным.
Моментные характеристики выхода
yt
системы определяются соотношениями:
M { yt } my (t ) C mx (t ) ;
M { yt0 ( yt0 )T } Dy (t ) C Dy (t ) C T ;
Ry (t , s) C R x (t , s) C T .
Спектральная плотность стационарного выходного процесса
yt
определяется формулой:
1
S y ()
Ry (t ) eit dt C (iE A)1 G 2w GT (iE AT )1 C T ,
2
где
C (sE A)1 G w
nt ,(dwt w nt dt )
- матричная передаточная функция системы от входа
в виде белого шума с единичной дисперсионной матрицей.
