САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ
КАФЕДРА ТЕОРИИ УПРАВЛЕНИЯ
Богославец Александра Игоревна
Выпускная квалификационная работа бакалавра
Стабилизирующая роль запаздывания в
системах обыкновенных дифференциальных
уравнений
Направление 010400
Прикладная математика, фундаментальная информатика
и основы программирования
Научный руководитель,
кандидат физ.-мат. наук,
доцент
Егоров А. В.
Санкт-Петербург
2016
Содержание
1. Введение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
2. Постановка задачи . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
3. Основные используемые понятия и теоремы . . . . . . . . . . . . . . .
9
4. Необходимое и достаточное условие устойчивости . . . . . . . . . . . . 12
5. Метод оценки запаздывания . . . . . . . . . . . . . . . . . . . . . . . . 19
6. Заключение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
Список литературы . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2
1. Введение
В данной работе рассматривается система линейных дифференциальных
уравнений с запаздывающим аргументом на примере модели шара с желобом. Она стабилизируется PD контроллером. Сохраняя коэффициенты в нём,
записывается уравнение пропорционального контроллера с запаздыванием.
Целью работы является нахождение запаздывания, при котором эта система
остается экспоненциально устойчивой, с помощью матриц Ляпунова и функционалов Ляпунова-Красовского полного типа. Под системой дифференциальных уравнений с запаздыванием обычно понимают систему дифференциальных уравнений, в которых неизвестная функция и ее производные входят
при различных значениях аргумента [2]. И если для линейных стационарных
систем без запаздывания проблема устойчивости решена, то для систем с запаздываниями исследование устойчивости по-прежнему актуально.
При решении задач стабилизации линейных стационарных систем с запаздыванием находят широкое применение два метода [8]. Первый подход основывается на том, что для любой линейной системы с запаздыванием можно
построить характеристический квазиполином, по расположению корней которого можно сделать вывод об устойчивости системы. Второй подход называется методом Ляпунова-Красовского и состоит в обобщении классического второго метода Ляпунова на случай систем с запаздывающим аргументом [5]. При таком подходе для определения устойчивости системы используются функционалы, определённые на множестве вектор-функций (функционалы Ляпунова-Красовского).
Работа состоит из двух основных частей. Первая часть содержит в себе
необходимое условие устойчивости, которым является положительная опре3
делённость специальной матрицы, построенной исключительно по матрице
Ляпунова. Также строится область экспоненциальной устойчивости в пространстве параметров заданной системы. Для того чтобы были видны точные границы областей устойчивости и неустойчивости, используется метод
D-разбиений [2], основная идея которого состоит в разбиении пространства
параметров системы кривыми, соответствующими тем значениям параметров, при которых характеристическое уравнение имеет чисто мнимые корни. При этом области обладают следующим свойством: если одна из точек
этой области соответствует экспоненциально устойчивой системе, то и все
остальные точки области соответствуют экспоненциально устойчивым системам, и наоборот. Если область, ограниченная такими кривыми, содержит в
себе точки, в которых необходимое условие устойчивости не выполняется,
то эта область является областью неустойчивости. Из-за того, что численно
невозможно проверить условие для всех точек, проверка осуществляется для
конечного числа равномерно распределённых по плоскости точек. Также в
первой части рассматривается критерий экспоненциальной устойчивости для
систем дифференциально-разностных уравнений с запаздывающим аргументом. В критерии сказано, что для того чтобы система была экспоненциально
устойчивой, необходимо и достаточно показать положительную определенность матрицы специального вида некоторой размерности r. Проблема заключается в том, какой именно нужно её брать. Поэтому находятся условия,
при которых размерность r уменьшается и высчисление специальной матрицы с такой размерность уменьшается. И проводится анализ графика области
экспоненциальной устойчивости, в котором точки из областей, подозрительных на устойчивость, проверяются критерием устойчивости и находятся все
4
запаздывания, при которых система является экспоненциально устойчивой.
Вторая основная часть посвящена методу оценки запаздывания с помощью
функционалов полного типа [1]. В ней рассматривается система, где запаздывание входит только в одну матрицу
ẋ(t) = (hP + Q)x(t) + A1 x(t − 1),
t ≥ 0.
Предполагается, что эта система экспоненциально устойчива при достаточно малом h. Находится оценка на возмущенную матрицу такая, что система
остается устойчивой при некотором условии.
Стоит отметить, что задача поиска запаздывания может быть рассмотрена
и на других примерах кроме модели шара с желобом и к ним могут быть
применимы методы, описанные в данной работе.
5
2. Постановка задачи
Рассматриваем задачу на упрощенной модели движения шара по желобу:
ẍ = g sin α.
Здесь α — угол наклона желоба к горизонту; g — ускорение силы тяжести.
Записываем систему дифференциальных уравнений в матричном виде:
x˙
0 1
x
0
1 =
1 + u,
(1)
x˙2
0 0
x2
1
x1
y = 1 0 .
x2
(2)
Здесь x(t) = (s(t), v(t))Т — состояние системы; x1 = x, x2 = ẋ; u = g sin α.
В качестве контроллера, стабилизирующего систему (1)-(2), берем пропорционально дифференциальный контроллер, так как пропорциональный контроллер сам по себе не стабилизирует заданную систему.
Запишем PD контроллер:
u = Kp y + Kd ẏ.
(3)
Заменяем PD контроллер (3) управлением с запаздыванием:
u = Kp y + Kd
y(t) − y(t − h)
.
h
(4)
Коэффициенты Kp и Kd берутся из выражения (3), h - запаздывание. Из-за
невозможности вычисления производной на практике, будем находить коэффициенты Kp и Kd для PD контроллера из выражения (3).
6
Для этого рассмотрим систему (1), где вместо u подставляется управление
вида (3):
0
x˙
0 1
x
1 =
1 + (Kp y + Kd ẏ).
1
x˙2
0 0
x2
(5)
Уравнение измерительного устройства (2) подставляем в систему (5):
x˙
0 1
x
0
x
x˙
1 =
1 + Kp 1 0 1 + Kd 1 0 1 .
x˙2
0 0
x2
1
x2
x˙2
Раскрываем скобки и приводим подобные слагаемые, в результате получаем:
0 1
x
x˙
1 .
1 =
Kp Kd
x2
x˙2
(6)
По критерию Гурвица для уравнения второго порядка необходимым и достаточным условием устойчивости является положительность всех коэффициентов характеристического уравнения. Поэтому запишем характеристический полином для системы (6):
f (λ) = λ2 − Kd λ − Kp .
Отсюда видно, что для положительности коэффициентов нужно брать Kp
и Kd отрицательными:
Kp < 0, Kd < 0.
(7)
Теперь необходимо подставить их в систему (5), где вместо ẏ берется
y(t)−y(t−h)
h
:
x˙1
0 1
x1
0
y(t)
−
y(t
−
h)
=
+ Kp y + Kd
,
h
x˙2
0 0
x2
1
x1
y = 1 0 .
x2
7
(8)
(9)
Замыкаем систему (8) показаниями прибора (9):
0 0
0
1
x(t) +
x(t − h).
ẋ(t) =
Kd
Kd
−h 0
Kp + h 0
Делаем замену z(t) = x(ht), которая не оказывает влияния на экспоненциальную устойчивость системы и после которой будем работать только с
запаздыванием, равным единице:
0 0
0
h
z(t − 1).
z(t) +
ż(t) =
−Kd 0
Kp h + Kd 0
(10)
Теперь h — параметр, который пробегает различные значения и важно
отметить, что от h зависит только матрица A0 .
8
3. Основные используемые понятия и теоремы
Рассмотрим систему дифференциально-разностных уравнений вида:
ẋ(t) = A0 x(t) + A1 x(t − 1),
(11)
где x — n-мерный вектор, A0 и A1 — вещественные постоянные матрицы
размерности n × n.
Пусть ϕ : [−1, 0] → Rn — начальная функция, принадлежащая пространству кусочно-непрерывных вектор-функций P C([−1, 0], Rn ), заданных на отрезке [−1, 0].
Пусть x(t, ϕ) — решение задачи Коши при t ≥ 0, соответствующее начальной функции ϕ для системы (11) с начальными условиями:
x(θ, ϕ) = ϕ(θ),
θ ∈ [−1, 0].
Основным инструментом решения поставленной задачи является применение матрицы Ляпунова. Приведем ее определение.
Определение 1 [1]. Будем говорить, что U (τ ) — матрица Ляпунова для
системы (11), если она удовлетворяет следующим свойствам:
1) Динамическое свойство:
U 0 (τ ) = U (τ )A0 + U (τ − 1)A1 , τ > 0;
2) Свойство симметрии:
U (τ ) = U T (−τ );
3) Алгебраическое свойство:
T
U (0)A0 + AT
0 U (0) + U (−1)A1 + A1 U (1) = −W .
9
Для построения матрицы Ляпунова необходимо убедиться в ее существовании и единственности. Для этого воспользуемся условием Ляпунова.
Определение 2 [1]. Будем говорить, что система (11) удовлетворяет
условию Ляпунова, если спектр системы
Λ = s | det sI − A0 − e−s A1 = 0
не содержит точки s0 , такой что -s0 также принадлежит спектру.
Будем использовать равномерную норму:
kϕk∞ = sup kϕ(t)k.
t∈[−1,0]
Определение 3 [1]. Система (11) называется экспоненциально устойчивой, если существуют константы γ ≥ 1 и σ>0 такие, что
kx(t, ϕ)k ≤ γe−σt kϕk∞ ,
t ≥ 0,
для любой начальной функции ϕ.
Теорема 1 [4]. Система (11) экспоненциально устойчива тогда и только
тогда, когда нули характеристического квазиполинома
p(s) = det(sI − A0 − e−s A1 )
(12)
лежат в левой открытой комплексной полуплоскости.
Определение 4
[1]. Будем говорить, что матрица K(t) размерности
n × n — фундаментальная матрица, если она удовлетворяет матричному
уравнению:
d
K(t) = K(t)A0 + K(t − 1)A1 , t ≥ 0,
dt
K(t) = 0n×n для t < 0, K(0) = I. Здесь I — единичная матрица.
10
Определение 5 [3]. Будем говорить, что симметрическая матрица M
размерности n × n положительно определена, если для всех ненулевых векторов z ∈ Rn :
z T M z > 0.
Определение 6 [3]. Будем говорить, что симметрическая матрица M
размерности n × n положительно полуопределена, если для всех ненулевых
векторов z ∈ Rn :
z T M z ≥ 0.
Определение 7 [4]. Функционал
g : P C([−1, 0], Rn ) → Rp ,
где p — некоторое натуральное число, будем называть непрерывным в нуле,
если для любого ε > 0 существует δ = δ(ε) > 0 такое, что при kϕk∞ < δ
выполнено неравенство:
kg(ϕ) − g(0)k < 0.
Теорема 2 (Красовский, 1956.) [5]. Существование непрерывных в нуле скалярных функционалов υ(ϕ) и ω(ϕ), ϕ ∈ P C([−1, 0], Rn ), удовлетворяющих условию υ(0) = ω(0) = 0 и связанных на решениях системы (11)
соотношением:
dυ(xt )
= −ω(xt ),
dt
для которых найдутся числа α > 0, β > 0 такие, что
υ(ϕ) ≥ αkϕ(0)k2 ,
ω(ϕ) ≥ βkϕ(0)k2 ,
гарантирует экспоненциальную устойчивость системы (11).
11
4. Необходимое и достаточное условие устойчивости
Необходимым условием устойчивости является положительная определенность данной матрицы, построенной по матрице Ляпунова.
Теорема 3 [4]. Если система (11) экспоненциально устойчива, то для любого натурального r выполнено условие
r
j−i
Kr = U
> 0,
r−1
i,j=1
которое может быть переписано
1
)
U ( r−1
U (0)
U T ( 1 ) U (0)
r−1
U T ( 2 ) U T ( 1 )
r−1
r−1
...
...
U T (1) U T ( r−2
r−1 )
(13)
в виде:
2
U ( r−1
) ...
1
U ( r−1
) ...
U (0)
...
...
...
r−3
U T ( r−1
) ...
U (1)
U ( r−2
)
r−1
r−3 > 0.
U ( r−1 )
...
U (0)
Так как мы рассматриваем систему (10), то и условие (13) применяем к
ней. Для того чтобы построить матрицу Kr , сначала нужно построить матрицу Ляпунова U (τ ), ассоциированную с положительно определенной матрицей
W, для системы (10). Пусть матрица W — это единичная матрица. Берем
произвольные отрицательные значения для коэффициентов Kp и Kd . Для
определенности, Kp = −1 и Kd = −0.01. Алгоритм вычисления матрицы Ляпунова основан на полуаналитическом методе из работы [1] и реализован в
математическом пакете Matlab. В точке τ = 0 матрица Ляпунова выглядит
следующим образом:
U (0) =
119.7176 −0.5000
.
−0.5000 119.7176
12
После того как построили матрицу Ляпунова переходим к построению матрицы Kr . В качестве размерности для этой матрицы берем r = 4. Применим
метод D-разбиений [2] и найдем условия разрешимости уравнения f (iω) = 0
при вещественном ω.
Характеристическое уравнение для системы (10):
f (s) = s2 − hKd − h2 Kp + hKd e−s .
(14)
Кривые D-разбиения имеют вид:
Kd = 0,
Kd = −
π 2 k 2 Kp h
−
,
2h
2
(15)
(16)
где k нечетное положительное,
πk
h= p
,
| Kp |
(17)
где k четное положительное.
Построим сначала график для кривых D-разбиения (см. рис. 1). Вся область, где Kd принимает положительные значения, не является экспоненциально устойчивой.
13
Рис. 1: Кривые D-разбиения (16) (17) при фиксированном Kp = −1
Важно отметить, что коэффициенты Kp и Kd постоянные, то есть не меняют своих значений. Но для наглядности и дальнейших исследований мы
зафиксировали только Kp , а Kd мы варьируем. Пустая область дальше на
рисунке будет значить выполнимость необходимого условия.
14
Рис. 2: Области устойчивости системы (4) при фиксированном Kp = −1
Проверяем необходимое условие (13) для каждой точки области. Не будем учитывать те точки, в которых это условие не выполняется. Благодаря
методу D-разбиений, нам не надо проверять каждую точку на экспоненциальную устойчивость, достаточно взять по одной точке из каждой области
и проверить [7]. Для этой проверки мы будем использовать необходимое и
достаточное условие экспоненциальной устойчивости.
15
Теорема 4 Пусть α1 ∈ (0, α1∗ ), α = α2 /α1 . Система (11) экспоненциально
устойчива тогда и только тогда, когда выполняется условие Ляпунова, и
матрица
Kr − α1 Ar
(18)
p
r = 1 + de (M + L)(α + α(α + 1) − L)e,
(19)
положительно определена, где
L
α1 kϕ(0)k2 ≤ v1 (ϕ) ≤ α2 kϕk2∞ ,
r
Ar = K T (τi )K(τj ) i,j=1 ,
M = kA0 k + kA1 k,
β
,
m+1
α1∗ =
скобки d·e обозначают функцию округления сверху, m - количество запаздываний, β — положительное число, такое что матрица L(β) положительно
полуопределена
W 0 ... 0
A + A0 A1 . . . Am
0
0 W ... 0
A1
0 ... 0
+β
.
L(β) =
. . . . . . . . . . . .
...
. . . . . . . . .
0 0 ... W
Am
0 ... 0
0
Здесь Kr — матрица из формулы (13), L такое, что kK (t)k ≤ L на (0,h),
v1 (ϕ) — функционал, который имеет вид:
R0
R0 R0 T
v1 (ϕ) = ϕT (0)U (0)ϕ(0) + 2ϕT (0) U (−s − h)A1 ϕ(s)ds +
ϕ (s1 )AT1 ×
−1
× U (s1 − s2 )A1 ϕ(s2 )ds2 ds1 +
R0
ϕT (s)W ϕ(s)ds.
−1
16
−1 −1
То есть для того чтобы найти значение α1 , нужно сначала найти значение
β. Будем его искать из уравнения det(L(β)) = 0. Подставляем в L(β) наши
матрицы:
1
0
L(β) =
0
0
0
h + Kp h + Kd 0
h + K p h + K d
0
−Kd
1 0 0
+β
0
−Kd
0
0 1 0
0
0
0
0 0 1
0 0 0
0
0
.
0
0
Если за γ обозначить h + Kp h + Kd , тогда характеристическое уравнение
будет иметь вид:
1 − β 2 Kd2 − β 2 γ 2 = 0.
Наименьшее β, при котором определитель обращается в ноль, является
наибольшим возможным значением для β. Тогда можно брать любое β из
промежутка от нуля до получившегося числа
1
β=p 2
,
Kd + γ 2
Или при обратной подстановке γ:
1
β=p 2
.
Kd + (h + Kp h + Kd )2
(20)
Соответственно, для каждого h подбираем свое β.
Хотим найти условия, при которых размерность r матрицы Kr будет относительно мала, чтобы не строить матрицу большой размерности и не проверять её с матрицей Ar на положительную определенность, когда можно построить матрицу меньшей размерности. Этого можно добиться за счет уменьшения α = α2 /α1 . Будем брать α1 очень близким к α1∗ . При увеличении β увеличивается α1 , а увеличения β можно добиться, рассмотрев знаменатель (20).
17
Из условия Kd + h(Kp + 1) = 0 можно найти соотношение между Kp и Kd :
Kp = −
Kd
− 1.
h
(21)
Будем находить α2 по следующей формуле из книги [1]:
α2 = ν(1 + a)2 + kW1 + W2 k.
Пусть
ν = max kU (θ)k.
(22)
θ∈[0,1]
Здесь a = kA1 k,
W = W0 + W1 + W2 .
(23)
Чтобы α2 было наименьшим, надо чтобы норма сумм матриц W1 и W2 была наименьшей среди других норм. Так как матрица W0 нигде не участвует,
то можно сделать её любой положительно определенной, но чтобы выполнялось условие (23). Будем работать с единичными матрицами, умноженными
на некоторый коэффициент, чтобы легче определить собственные числа этой
матрицы. Поэтому выберем их следующим образом:
98
1
1
0
0
100
100
100
W0 =
, W1 =
, W2 =
98
1
0 100
0 100
0
0
1
100
.
Таким образом, при выполнении условий (20)-(21) размерность матрицы
можно получить сравнительно небольшую.
Например, при h = 1, α1 = 35.3553, α2 = 111.0862 получили r = 46.
18
5. Метод оценки запаздывания
Рассмотрим номинальную систему вида:
ẋ(t) = (εP + Q)x(t) + A1 x(t − 1), t ≥ 0,
(24)
где ε принимает достаточно малые значения, при которых система (24) экспоненциально устойчива.
P =
0
1
Kp 0
,Q =
0
Kd
0
,
0
Возмущенная система имеет вид:
ẏ(t) = (εP + Q + ∆)y(t) + A1 y(t − 1).
(25)
Здесь матрица ∆ такая, что выполняется условие:
k∆k = k(h − ε)P k ≤ ρ.
(26)
Пусть система (24) экспоненциально устойчива. Хотим найти условие на
ρ такое, что система (25) остается устойчивой для всех ∆, удовлетворяющих (26). С этой целью используем функционал полного типа v(ϕ), вычислить который можно по формуле [1]:
Z0
v(ϕ) = v0 (ϕ) +
[ϕ(θ)(W1 + (1 + θ)W2 )]ϕ(θ)dθ.
−1
Вычисляем производную по времени функционала в силу возмущенной
системы (25). Оцениваем производную сверху [1].
Теорема 5 [1]. Пусть система (11) экспоненциально устойчива. Заданы
положительно определенные матрицы W0 , W1 , W2 , система (25) остается
экспоненциально устойчивой для всех ∆, удовлетворяющих (26), если выполнены следующие неравенства:
19
1) λmin (W0 ) ≥ ν[2ρ + hρa];
2) λmin (W2 ) ≥ νρa;
Здесь a = kA1 k, ν строится по матрице Ляпунова, которая, в свою очередь,
строится по матрицам (εP + Q) и A1 . При ε = 10−3 система экспоненциально
устойчива.
Находим наибольшее возможное ρ, удовлетворяющее неравенствам из теоремы 5, считая матрицы W0 , W1 , W2 единичными, домноженными на число,
и удовлетворяющими уравнению:
ρ=
λmin (W2 ) λmin (W0 )
=
.
νa
ν(2 + a)
Задаем W1 с небольшими собственными числами, так как она не участвует
непосредственно в процессе нахождения ρ, но влияет на другие матрицы из
(23). Пусть λmin (W1 ) =
1
100 ,
тогда собственные числа W0 и W2 будут предста-
вимы как λmin (W0 ) + λmin (W2 ) =
99
100 .
И сами матрицы выглядят следующим
образом:
W1 =
99
1
I, W0 + W2 =
I.
100
100
Теперь подставляем вместо W0 её выражение через матрицу W2 :
99
I − W2 )
λmin (W2 ) λmin ( 100
ρ=
=
.
νa
ν(2 + a)
Преобразуем выражение:
2+a
99
λmin (W2 ) − λmin (
I − W2 ) = 0.
a
100
Следовательно,
λmin (W2 ) =
99a
.
200 + 200a
20
Тогда
λmin (W0 ) =
99
− λmin (W2 ).
100
Формируем матрицы как было указано в предыдущей главе:
0.9851
0
0.0049
0
, W2 =
.
W0 =
0
0.9851
0
0.0049
Отсюда
ρ = 4.8765−6 .
Выбираем h исходя из требований теоремы Красовского (Теорема 2) о том,
что производная должна быть отрицательно определена.
h≤
ρ + εkP k
ρ
=ε+
.
kP k
kP k
Тогда окончательно приходим к тому, что h = 0.001. Недостаток в том, что
h получается малым, а достоинство метода в том, что значение h вычисляется
быстро и просто.
21
6. Заключение
Нами были исследованы два метода для нахождения значения h критического, при котором система теряет экспоненциальную устойчивость. В
сравнении первый метод лучше, он дает хорошую оценку на h, но посчитать
матрицу и проверить её на положительную определенность с большой размерностью не всегда удается быстро. Со вторым методом таких трудностей
не возникает, не приходится долго ждать результата, но он получается консервативным.
В качестве направлений дальнейших исследований отметим, что во втором
методе можно попытаться найти матрицы W0 , W1 , W2 , которые не являлись
бы единичными, умноженными на некоторый коэффициент. Также можно
находить итерационно каждый раз новое h, подставляя его в систему для
определения h критического. Кроме того, возможно уменьшение размерности
r.
22
Список литературы
1. Kharitonov V. L. Time-delay systems: Lyapunov functionals and matrices //
Birkhauser, Basel. 2013.
2. Эльсгольц А. Э., Норкин С. Б. Введение в теорию дифференциальных
уравнений с отклоняющимся аргументом // М.: Наука. c.296. 1971.
3. Bhatia R. Positive definite matrices // Princeton Series in Applied
Mathematics, c.1-3. 2007.
4. Егоров А. В. Новые условия экспоненциальной устойчивости линейных
систем с запаздыванием // дисс. на соиск. уч. ст. канд. физ.-мат. наук.
2013.
5. Красовский Н. Н. О применении второго метода Ляпунова для уравнений
с запаздываниями времени // Прикладная математика и механика. 1956.
Т. 20. с. 315-327.
6. Bellman R., Cooke K. L. Differential-Difference Equations. // New
York/London 1963. Academic Press.
7. Воронов А. А. Теория автоматического управления // М.: Высшая школа,
Часть 1. c.153-165. 1977.
8. Чашников М. В. Анализ устойчивости линейных систем с запаздывающим
аргументом // дисс. на соиск. уч. ст. канд. физ.-мат. наук. СПб., 2010. с.94.
23
Отзывы:
Авторизуйтесь, чтобы оставить отзыв