Правительство Российской Федерации
Федеральное государственное бюджетное образовательное учреждение
высшего образования «Санкт-Петербургский государственный
университет»
Математическое обеспечение и администрирование информационных
систем
Кафедра информационно-аналитических систем
Новиков Георгий Сергеевич
Разработка алгоритма выпуклой
оптимизации для минимизации всплесков в
динамической системе управления
Бакалаврская работа
Научный руководитель:
д. ф.-м. н., профессор Граничин О. Н.
Рецензент:
к.ф.-м.н., доцент Михайлова Е.В.
Санкт-Петербург
2016
SAINT-PETERSBURG STATE UNIVERSITY
Sub-Department of Analytical Information Systems
Novikov Georgii
Development of convex optimisation
algorithm for deviation minimisation in
dynamical control system
Graduation Thesis
Scientific supervisor:
Prof. Oleg Granichin
Reviewer:
associate professor, PhD Elena Mikhailova
Saint-Petersburg
2016
Оглавление
1. Введение
1.1. Актуальность . . . . . . . . . . . . . . . . . . . . . . . . .
1.2. Предварительный пример . . . . . . . . . . . . . . . . . .
1.3. Структура работы . . . . . . . . . . . . . . . . . . . . . . .
4
4
5
5
2. Задача о минимизации всплесков в системе управления
7
3. Случай нулевого собственного числа
3.1. Поведение решения . . . . . . . . . . . . . . . . . . . . . .
3.2. Постановка задачи . . . . . . . . . . . . . . . . . . . . . .
3.3. LMI для инвариантного подпространства . . . . . . . . .
9
9
10
10
4. Движение группы мобильных агентов
4.1. Задача о переводе в заданное расположение . . . . . . . .
4.2. Теоретические основы . . . . . . . . . . . . . . . . . . . . .
4.3. Ход Алгоритма . . . . . . . . . . . . . . . . . . . . . . . .
12
12
13
14
5. Программное обеспечение
15
6. Пример использования
6.1. Случай одного мобильного агента . . . . . . . . . . . . . .
6.2. Случай трех мобильных агентов . . . . . . . . . . . . . .
17
17
17
7. Заключение
23
Список литературы
25
3
1. Введение
1.1. Актуальность
Теория управления - молодая наука, находящаяся в процессе интенсивного развития. При этом существенно меняются взгляды на предмет
и основные проблемы данной дисциплины, равно как и используемый
математический аппарат. В XIX веке было введено понятие устойчивости регулируемого процесса и получены первые критерии устойчивости.
В 30-е годы в основном рассматривались частотные методы и, соответственно, частотные критерии устойчивости. В 50-е годы произошло
обновление в теории управления - появился новый аппарат описания
систем - в пространстве состояний. Иначе говоря, движение системы
подчиняется обыкновенному дифференциальному уравнению, в правой
части которого стоит функция, которая может выбираться проектировщиком (управление). Появились различные требования к функции
управления - возможность стабилизации системы при внешних возмущениях [1], робастная стабилизация и т.д., а также различные критерии
оптимальности функции управления. Один из таких критериев - скорость сходимости системы к своей точке равновесия. Этот критерий не
учитывает поведение решения на ранних этапах сходимости, а смотрит
лишь на асимптотическое поведение уравнения. Оказывается, это не
всегда допустимо - при определенных условиях, на начальных этапах
решение может сильно отклоняться от точки сходимости, причем, как
правило, чем больше скорость сходимости - тем больше отклонения.
Минимизации таких отклонений и посвящена эта статья.
Основная техника исследования задач в этой статье связана с линейными матричными неравенствами (С. Бойд с соавторами [5], Д. В.
Баландин и М. М. Коган [2]) и задачей полуопределенного программирования [11]. Для решения полученных задач разработаны мощные
оптимизационные процедуры, основанные на методе внутренней точки
[3], [6], [7]. В качестве вычислительного средства были использованы две
платформы - система Matlab вместе с программным пакетом SeDuMi
4
[8] и пакетом cvx [9], [10]. А также система ipython notebook для языка
python [13] с использованием программных пакетов numpy [14] и picos
[15].
1.2. Предварительный пример
Рассмотрим систему управления
0
1
0
0
0
0
1
0
ẋ =
0
0
0
1
−100 −20000.02 −1000040.0001 −2000.02
x
(1)
Собственные числа матрицы этой системы равны λ1 = λ2 = −1000
и λ3 = λ4 = −0.01. Норма Евклида решения с начальным положением
(0, 0, 1, 0) показана на рисунке 1, и с начальным положением (0, 1, 0, 0)
на рисунке 2. Обратите внимание, что сильное отклонение происходит в
разные моменты времени t. И величина отклонения, и момент времени,
в который оно происходит, зависят от начального положения системы.
1.3. Структура работы
В этой работе в разделе (2) сформулированы ключевые понятия
и теоремы. Затем в разделе (3) рассмотрен случай нулевого собственного числа матрицы и способы работы с такими матрицами. Далее в
разделе (4) построен алгоритм сходимости в требуемую конфигурацию
для агентов, закон движения которых описывается формулой ẋ = v.
В разделе (5) описано программное обеспечение, разработанное в ходе
выполнения работы. В разделе (6) находятся результаты моделирования для задачи о перемещении группы мобильных агентов в заданное
расположение.
5
Рис. 1: Система (1), x(0) = (0, 0, 1, 0)
Рис. 2: Система (1), x(0) = (0, 1, 0, 0)
6
2. Задача о минимизации всплесков в системе управления
Рассмотрим линейную непрерывную стационарную систему управления с одномерным управлением
ẋ = Ax + bu, A ∈ Rn×n , b ∈ Rn , u ∈ R
(2)
В случае, если вектора A, Ab, . . . , An−1 b линейно независимы, то система (2) является управляемой и по теореме о размещении полюсов существует линейная обратная связь, задаваемая вектором K такая, что
точка x = 0 является асимптотически устойчивой в уравнении
ẋ = Ax + bK T x.
(3)
Более того, выбором вектора K можно получить систему с любым набором собственных чисел L ∈ C, тем самым получив любую наперед
заданную скорость сходимости. Назовем величиной отклонения решения системы (3) от точки равновесия величину
sup sup |x(t, x0 )|.
(4)
0≤t≤∞ |x0 |=1
Cистемы с большой скоростью сходимости, как правило, сильно отклоняются от точки сходимости на начальном этапе, что было показано Р.
Измаиловым в [12]. А именно, им была доказана следующая теорема:
Theorem 1 ∃γ = γ(A, b) > 0 такая, что если {λ1 , . . . , λn } - собственные значения системы (3) такие, что Reλi ≤ −σ < 0, i = 1, n, то верно
неравенство
sup
|x(t, x0 )| ≥ λγ n−1 .
(5 )
0≤t≤1/σ,|x0 |=1
Встает вопрос нахождения вектора K не просто с условием получения требуемой скорости сходимости (требуемых действительных частей
собственных чисел матрицы), но также и минимизации величины возможных отклонений от точки сходимости. Данная задача исследована
7
в статье Б. Т. Поляка и Г. В. Смирнова [4]. В ней дан метод получения
хорошего кандидата на вектор обратной связи. Его можно получить
при помощи LMI подхода, решив задачу полуопределенного программирования.
Рассмотрим задачу минимизации
min α,
(6)
P (A + bK T )T + (A + bK T )P −2σP.
(7)
I P αI.
(8)
Здесь α ≥ 1, K ∈ Rn , P = P T ∈ Rn×n - переменные оптимизации. Неравенство Ляпунова (7) гарантирует, что действительные части всех собственных чисел системы (скорость сходимости всей системы) не превосходит σ и одновременно то, что эллипсоид E = {x ∈ Rn : hx, P −1 x ≤ 1i}
является инвариантным множеством данной системы. Условие (8) гарантирует нам минимальность эксцентриситета эллипсоида E (а так
как E инвариантное множество, то и минимальность величины отклонений, оцененных в таком виде). Заметим, что в таком виде задача
минимизации не линейна. Это можно преодолеть произведя замену переменных y = P K. Тогда неравенство (7) приобретет вид
AP + P AT + by T + ybT ≤ −2σP
(9)
В таком виде это линейная задача полуопределенного программирования, а ее решение является хорошим кандидатом в вектор обратной
связи. Более того, в [4] показано, что это решение асимптотически (при
больших σ) дает величину отклонений того же порядка, что и нижняя
оценка Р. Измаилова (5).
8
3. Случай нулевого собственного числа
3.1. Поведение решения
Рассмотрим матрицу A с простым собственным числом, равным нулю, и соответствующим собственным вектором x∗ . Для начала посмотрим, как будет вести себя решение с начальными данными x0 . Далее
будем рассматривать разложение пространства в прямую сумму двух
инвариантных подпространств матрицы А - собственного подпространства, соответствующего с.ч. 0 матрицы А и инвариантного подпространства, включающего в себя все остальные присоединенные вектора матрицы А. Назовем их нулевое и ненулевое инвариантные подпространства соответственно. Пусть x0 раскладывается в сумму x0 = αx∗ + βx⊥ ,
где x⊥ принадлежит ненулевому инвариантному подпространству. Тогда
(10)
x(t, x0 ) = eAt x0 = αeAt x∗ + βeAt x⊥
At ∗
e x =
∞ i i
X
tA
i!
i=0
x∗ = x∗
то
(11)
x(t, x0 ) = αx∗ + βeAt x⊥
Рассмотрим, как действует матричная экспонента на вектор из ненулевого инвариантного подпространства. Для этого рассмотрим
разло!
0 0
жение матрицы A в жорданов базис A = S −1
S. Матрица S
0 AJ
!
0
действует на x⊥ следующим образом: Sx⊥ =
из чего следует:
x
S −1
eAt x⊥ = e
0 0
0 AJ
S
x⊥ = S −1
1 0
0 e AJ t
9
!
Sx⊥ = S −1
0
eAj t x
!
(12)
Из этого следует, что решение из начальной точки x0 не меняет своей
составляющей x∗ , а изменяет только свою составляющую в ненулевом
инвариантном подпространстве, причем по закону ẋ = Aj x. Следующей задачей будет научиться работать со сходимостью в инвариантном
подпространстве.
3.2. Постановка задачи
Дана система ẋ = Ax + bK T x. Матрица А имеет простое собственное
число, равное нулю, соответствующее собственному вектору x∗ . Необходимо найти такой вектор K, что матрица (A+bK T ) не поменяет нулевое
и ненулевое инвариантные подпространства, собственные числа ненулевого инвариантного подпространства будут иметь вещественные части,
меньшие −σ для данного σ и величина отклонений системы от точки
сходимости будет минимальной.
3.3. LMI для инвариантного подпространства
Мы хотим решить LMI задачу для матрицы Aj , получив после этого
каким то образом вектор K. Наложим условия на вектора b, K, чтобы
они не испортили нулевое собственное число и соответствующие два
инвариантных подпространства:
A + bK T = S −1 JS + bK T = S −1 (
0 0
0 AJ
!
+ SbK T S −1 )S
(13)
Из этого следует, что:
n
X
(S −1 b)j Ki Si1 = 0, j = 1, n
(14)
(S −1 b)1 Ki Sij = 0, j = 1, n
(15)
i=1
n
X
i=1
Из (14) следует, что K T Si1 = 0 (иначе S −1 b ≡ 0 ⇔ b ≡ 0), из (15)
следует, что (S −1 b)1 = 0 (иначе K T S ≡ 0 ⇔ K ≡ 0). Вектор b является
10
входным параметром, поэтому условие на него должно быть соблюдеb
но изначально. Пусть LMI метод как результат выдал нам вектор K
(вспомним, что он размера n − 1, так как матрица Aj ∈ Rn−1×n−1 , а
также то, что он записан для нестандартного базиса матрицы А). Посмотрим на то, как должен выглядеть этот вектор в стандартном базисе
с учетом условия (14):
KT S1 = 0
K T S 2 = K
b1
...
K T S n = K
b n−1
⇔ ST K =
0
b
K
!
⇔ K = (S T )−1
0
b
K
!
(16)
В итоге получается, что нужно разложить матрицу A в базис, выделив нулевое собственное подпространство (например, подойдет жорданов базис). Подать на вход задаче LMI матрицу Aj , вектор S −1 b без
первого элемента (который по условию (15) должен быть равен нулю)
и необходимую скорость сходимости σ, и по результирующему
вектору
!
0
b получить вектор K по формуле K = (S T )−1
K
.
b
K
11
4. Движение группы мобильных агентов
4.1. Задача о переводе в заданное расположение
Рассмотрим систему мобильных агентов, движение которых задается уравнениями
ẋ = v
v̇ = q
(17)
x ∈ Rn - координаты системы, v ∈ Rn - скорости, q ∈ Rn - ускорения, которые являются нашим воздействием на эту систему. Я буду
рассматривать случай, когда q линейно зависит от x и v. Процесс будет
описываться следующей системой:
˙ !
x
=
v
0 I
Av Bv
!
x
v
!
(18)
Алгоритм состоит в следующем: наложим условия на Av и!Bv таким
x∗
образом, что для наперед заданного вектора x∗ вектор
являет0
ся собственным вектором системы (18), соответствующим собственному
числу ноль. x∗ - это некоторая конфигурация, в которую система долж!
x0
на сойтись. Но в соответствии с (11) для начального положения
v0
точкой сходимости является ее составляющая
в нулевом инвариантном
!
x∗
подпространстве, а не само
. Поборем это следующим способом:
0
добавим к системе новую размерность - управляющий параметр u ∈ R,
изменение которого также линейно зависит от полного состояния. Тогда
система примет следующий вид:
˙
x
0 I 0
x
v = Av Bv cv v ,
u
au bu cu
u
(19)
где au ∈ Rn , bu ∈ Rn , cu ∈ R. Использовать управляющий параметр мы
12
x∗
0
!
будем следующим образом: во-первых, вместо вектора
собствен
x∗
ным вектором системы будет вектор 0 . Во-вторых, по начальному
1
x0
v0
положению
!
можно вычислить (см. (4.2)) такой управляющий
∗
x0
x
параметр u0 , что v0 = 0 + yb, где yb принадлежит ненулевоu0
1
му инвариантному подпространству,
а то-есть система асимптотически
∗
x
сойдется в состояние 0
1
4.2. Теоретические основы
Сначала сформулируем условия на Av , Bv , cv , au , bu , cu , чтобы получить необходимый собственный вектор:
0 I 0
x∗
A x∗ + c = 0
v
v
Av Bv cv 0 = 0 ⇔
au x∗ + cv = 0
au b u c u
1
(20)
Теперь по x0 и v0 найдем необходимое начальное значение управляющего параметра. Пусть матрица S - матрица перехода в базис, который разделяет два инвариантных подпространства (подпространство
собственного числа 0 и остальных собственных чисел). Тогда
α
x
1 0
x0
S
. . . = v0
−1
⇔ S1 v0
(21)
=1
α
u
n
0
u0
α = 1
1
13
4.3. Ход Алгоритма
В алгоритме можно выделить три части:
1. Построение матрицы
2. Подготовка начального положения
3. Процесс сходимости
Первый этап - построение матрицы с нужным собственным вектором
и собственным числом и последующая оптимизация части матрицы,
действующей в нужном инвариантном подпространстве, при помощи
метода LMI. Это нужно сделать только один раз для конкретной конфигурации x∗ . Второй этап - вычисление u0 по заданным x0 и v0 . Третий
этап - процесс сходимости по закону (19). При построенной матрице на
первом этапе, второй и третий этапы можно повторять сколько угодно
раз, перестраивать матрицу каждый раз не нужно.
14
5. Программное обеспечение
В ходе работы было написано две версии вычислительной программы с идентичными интерфейсами - на языке MATLAB и на языке
Python в среде разработки jupyter notebook. Всего можно выделить
несколько базовых функций
1. Функция решения задачи полуопределенного программирования в
форме (9). На вход она получает матрицу A ∈ Rn×n , вектор b ∈ Rn
и число σ ∈ R+ . Она приводит задачу к необходимому формату,
после чего, используя внешний алгоритм минимизации, решает ее
и возвращает искомый вектор K ∈ Rn . Как расширение данная
функция может принимать на вход вместо вектора b матрицу B ∈
Rn×m и возвращать матрицу K ∈ Rn×k решая тем самым более
общую задачу поиска управления для системы ẋ = Ax + BK T x.
2. Функция, находящая базис, разделяющий нулевое и ненулевое инвариантные подпространства для данной матрицы А. В общем
случае результатом является овеществленный жорданов базис матрицы А. В нашем же случае известно, что у матрицы А имеется
простое собственное число 0, соответствующее собственному вектору x∗ . Так что процедуру поиска можно сильно упростить. Заметим, что для любого x = αx∗ +βb
x, где x
b принадлежит ненулевое
инвариантное подпространство, Ax = βAb
x. Поэтому базис можно
набрать из случайно сгенерированного набора векторов Ax.
3. Функция, моделирующая поведение системы мобильных агентов.
На вход принимает матрицу A ∈ R2n∗m+1×2n∗m+1 , собственный вектор y ∈ R2n∗m+1 , вектора x0 ∈ Rn∗m и v0 ∈ Rn∗m - начальное состояние системы и массив моментов времени t. n и m здесь - количество
агентов и размерность пространства соответственно. Результатом
является массив векторов состояний системы в данные моменты
времени, полученный в соответствии с законом (19).
4. Функция, находящая по заданному набору собственных чисел для
15
системы (2) необходимый вектор K. Метод описан в статье Б. Поляка [4]
16
6. Пример использования
Рассмотрим задачу перевода системы мобильных агентов в заданное расположение и для заданной скорости сходимости σ сравним поведение двух систем - системы с одинаковыми собственными числами
равными −σ и системы, полученной из нее применением подхода LMI.
6.1. Случай одного мобильного агента
Рассмотрим случай одного мобильного агента на!плоскости. Необ0
ходимое расположение задается вектором x∗ =
. Матрица А сге10
нерирована случайно с точностью до удовлетворения
условиям
(20).
!
!
−5
−7
Начальное состояние системы x0 =
, v0 =
, σ = 0.1. На
−5
−3
рисунках 3 и 4 изображены траектории агента для системы без оптимизации и с оптимизацией соответственно. для. На рисунках 5 и 6
изображены графики зависимости расстояния состояния системы до
точки сходимости от времени. Четко видно, что из-за недостаточной
силы воздействия на агента в неоптимизированной системе, в процессе
положение сильно отклоняется от точки сходимости. Тогда как в случае оптимизированной системы агент очень быстро с небольшими отклонениями сходится в нужную точку. Ненулевые собственные числа
матрицы A + bK T примерно равны −6.57, −0.36, −0.43 ± 0.85i.
6.2. Случай трех мобильных агентов
Более интересен случай трех и более мобильных агентов. Матрица
системы (19) будет иметь размер хотя бы 12. В этом случае алгоритм построения матрицы с данным набором собственных чисел перестает работать - из-за возведения матрицы в большую степень погрешность вычислений превосходит необходимую точность. LMI подход имеет более
хорошие показатели устойчивости, из-за чего может быть использован
для построения матриц для более высокоразмерных задач. Рассмотрим
17
500
0
-500
-1000
-1500
-2000
-2500
-2000
-1800
-1600
-1400
-1200
-1000
-800
-600
-400
-200
Рис. 3: Траектория движения агента в неоптимизированной системе
18
0
12
10
8
6
4
2
0
-2
-4
-6
-6
-5
-4
-3
-2
-1
Рис. 4: Траектория движения агента в оптимизированной системе
19
0
3000
2500
2000
1500
1000
500
0
0
10
20
30
40
50
60
70
80
90
100
Рис. 5: Зависимость до точки сходимости от времени в неоптимизированной системе
20
18
16
14
12
10
8
6
4
2
0
0
2
4
6
8
10
Рис. 6: Зависимость расстояния до точки сходимости от времени в оптимизированной системе
100
−5
−7
0
−5
−3
0
5
2
следующую систему: x∗ =
, x0 =
, v0 =
,σ =
0
−5
3
−100
0
0
0
0
0
0.01. Это задача выстроить три мобильных агента в линию в точки с
координатами (-100, 0), (0, 0) и (100, 0). На риcунке 7 изображены траектории агентов. По сравнению со случаем одного агента траектории
сильно усложнились. На рисунке 8 изображена зависимость расстояния
от состояния системы до точки сходимости от времени - отклонение не
превышает 4 норм равновесного состояния.
21
300
200
100
0
-100
-200
-300
-400
-400
-300
-200
-100
0
100
200
Рис. 7: Траектории агентов в системе из трех агентов
22
300
500
450
400
350
300
250
200
150
100
50
0
0
10
20
30
40
50
60
70
80
90
100
Рис. 8: Зависимость расстояния до точки сходимости от времени в системе из трех агентов
23
7. Заключение
В работе было проведено несколько исследований. Разработан способ применения метода поиска оптимального управления в смысле минимизации величины отклонения системы от точки равновесия к случаю матриц с простым нулевым собственным числом. Также разработана система, решающая задачу перемещения группы мобильных агентов в заданное расположение и применен к ней разработанный метод
минимизации отклонений. Написал программа, производящая необходимые вычисления на языке MATLAB и языке python и произведен ряд
экспериментов по сравнению движений агентов в соответствии с придуманной системой с и без применения метода минимизации отклонений.
Получены улучшения для оптимизированных систем как по скорости
сходимости, так и по величине всплесков.
24
Список литературы
[1] Б. Т. Поляк, М. В. Хлебников, and П. С. Щербаков. Управление
линейными системами при внешних возмущениях: Техника линейных матричных неравенств. М.: ЛЕНАНД, 2014.
[2] Д. Баландин and М. Коган. Синтез законов управления на основе
линейных матричных неравенств. ЛитРес, 2016.
[3] Нестеров Ю. Е. Методы выпуклой оптимизации. М.: МЦНМО,
2009.
[4] Boris T. Polyak and Georgi V. Smirnov.
Large deviations
in continuous-time linear single-input control systems.
{IFAC}
Proceedings Volumes, 47(3):5586 – 5591, 2014. 19th {IFAC} World
Congress.
[5] S. Boyd, L. El Ghaoui, E. Feron, and V. Balakrishnan. Linear Matrix
Inequalities in System and Control Theory, volume 15 of Studies in
Applied Mathematics. SIAM, Philadelphia, PA, June 1994.
[6] Stephen Boyd and Lieven Vandenberghe. Convex Optimization.
Cambridge University Press, New York, NY, USA, 2004.
[7] Yurii Nesterov, Arkadii Nemirovskii, and Yinyu Ye. Interior-point
polynomial algorithms in convex programming, volume 13. SIAM, 1994.
[8] Jos F. Sturm. Using sedumi 1.02, a matlab toolbox for optimization
over symmetric cones, 1998.
[9] Michael Grant and Stephen Boyd.
CVX: Matlab software for
disciplined convex programming, version 2.1. http://cvxr.com/cvx,
March 2014.
[10] Michael Grant and Stephen Boyd.
Graph implementations for
nonsmooth convex programs. In V. Blondel, S. Boyd, and H. Kimura,
editors, Recent Advances in Learning and Control, Lecture Notes
25
in Control and Information Sciences, pages 95–110. Springer-Verlag
Limited, 2008. http://stanford.edu/~boyd/graph_dcp.html.
[11] A. Ben-Tal and A. Nemirovski.
Lectures on Modern Convex
Optimization: Analysis, Algorithms, and Engineering Applications.
MPS-SIAM Series on Optimization. Society for Industrial and Applied
Mathematics, 2001.
[12] R. N. Izmailov. Peak effect in stationary linear-systems in scalar inputs
and outputs. Automation and Remote Control, 48(8):1018–1024, 1987.
[13] Fernando Pérez and Brian E Granger. Ipython: a system for interactive
scientific computing. Computing in Science & Engineering, 9(3):21–29,
2007.
[14] Travis E Oliphant. A guide to NumPy, volume 1. Trelgol Publishing
USA, 2006.
[15] Guillaume Sagnol. Picos, a python interface to conic optimization
solvers.
Technical report, Technical Report 12-48, ZIB, 2012.
http://picos. zib. de, 2012.
26
Отзывы:
Авторизуйтесь, чтобы оставить отзыв