О внутренних оценках множеств достижимости многошаговых систем при групповом движении Е.К. Костоусова kek@imm.uran.ru Институт математики и механики им. Н.Н. Красовского УрО РАН (Екатеринбург) Аннотация Рассмотрена задача достижимости для линейных многошаговых управляемых систем, описывающих так называемое групповое дви- жение, при котором объекты попарно не сближаются, но и не слишком отдаляются друг от друга. Условия группового дви- жения приводят к невыпуклым фазовым ограничениям специ- ального типа, в результате чего множества достижимости могут быть невыпуклыми и даже несвязными. Получающиеся трубки достижимости могут быть представлены в виде ветвящихся тру- бок, сечения которых являются выпуклыми политопами, а мно- жества достижимости представимы в виде объединения выше- упомянутых сечений. В работе предлагаются некоторые алгорит- мы построения внутренних полиэдральных оценок множеств до- стижимости в виде объединений параллелепипедов, построенных на основе “элементарных” параллелотопозначных оценок для ре- зультатов теоретико-множественных операций с параллелелотопа- ми/параллелепипедами и полосами. Ключевые слова: множества достижимости; многошаговые си- стемы; групповое движение; внутренние полиэдральные оценки. 1 Введение Решение многих задач гарантированного управления и оценивания в условиях неопределенности связано с построением трубок траекторий (многозначных функций, описывающих, например, динамику множеств достижимости, разрешимости, информационных областей) [1–4]. Практическое построение вышеупомя- нутых трубок может быть достаточно затруднительно даже для линейных систем, особенно для систем с фазовыми ограничениями, систем большой размерности. Среди многих численных методов решения таких задач активно развиваются методы, основанные на аппроксимации множеств более простыми областями фиксированной формы, такими как эллипсоиды, параллелепипеды, параллелотопы, зонотопы (см., напри- мер, [1–3,5–11] и библиографию, приведенную там; здесь и ниже для примера дано только несколько ссылок из большого числа работ). Такие методы позволяют находить и эффективно использовать аппроксимации исследуемых трубок и множеств с помощью не слишком затратных вычислительных средств и развиваются как некоторая альтернатива к разработке других численных методов, предназначенных для получения как c by the paper’s authors. Copying permitted for private and academic purposes. Copyright ° In: G.A. Timofeeva, A.V. Martynenko (eds.): Proceedings of 3rd Russian Conference "Mathematical Modeling and Information Technologies" (MMIT 2016), Yekaterinburg, Russia, 16-Nov-2016, published at http://ceur-ws.org 50 можно более точных аппроксимаций (например, путем аппроксимации множеств политопами с большим количеством вершин и граней [12, 13]), но могущих потребовать большого объема вычислений, особенно для систем большой размерности. Упомянем, что кроме указанного выше гарантированного подхода к ре- шению задач управления и оценивания в условиях неопределенности существует стохастический подход, при котором используются вероятностные характеристики неизвестных входных данных. Решения задач оценивания состояния линейных многошаговых систем при одновременном наличии детерминированных и случайных помех описано, например, в [14, 15]. Множества достижимости по начальным данным и их аппроксимации для стохастических систем рассматривались, например, в [16]. В настоящей работе рассматривается задача достижимости для линейных многошаговых управляемых систем, описывающих так называемое групповое движение. А именно, рассматривается конечное множе- ство управляемых объектов при условии, что их траектории должны не слишком сближаться, но и не слишком отдаляться друг от друга. Мотивация задач такого рода вытекает из активно развиваемой в последнее время теории управления групповым движением (см., например, [5, 17–19]). Предполагается, что начальные состояния и управления объектов стеснены данными параллелепипедозначными ограниче- ниями. Вышеупомянутые условия группового движения приводят к невыпуклым фазовым ограничениям специального типа, в результате чего множества достижимости могут быть невыпуклыми и даже несвяз- ными. Получающиеся трубки достижимости могут быть представлены в виде ветвящихся трубок, сечения которых являются выпуклыми политопами, а множества достижимости представимы в виде объединения вышеупомянутых сечений [20]. В работе предлагаются некоторые алгоритмы построения внутренних поли- эдральных оценок множеств достижимости в виде объединений параллелепипедов, построенных на основе “элементарных” параллелотопозначных оценок [9–11] для результатов теоретико-множественных операций с параллелелотопами/параллелепипедами и полосами. Ниже используются следующие обозначения: Rn и Rn×m — линейные пространства вещественных n-векторов и n×m-матриц соответственно; > — знак транспонирования; kxk2 = (x> x)1/2 и kxk∞ = max1≤i≤n |xi | — разные нормы вектора x = (x1 , . . . , xn )> ; e = (1, 1, . . . , 1)> ; A = {aji } = {a1 ; a2 ; . . . an } — матрица с элементами aji и со столбцами aj (верхним индексом нумеруются столбцы, нижним — ком- поненты векторов); 0 — нулевая матрица (вектор) произвольной размерности; I — единичная матрица; A = diag {Aj } — диагональная Pm матрица A с блоками Aj на диагонали; det A — определитель матрицы A; j kAk = max1≤i≤n j=1 |ai | — норма матрицы A∈Rn×m , индуцированная нормой kxk∞ ; int X — совокуп- ность внутренних точек множества X ⊂ Rn ; M (H) — мощность (число элементов) конечного множества H; vol P — объем параллелепипеда P . 2 Постановка задачи > > Предполагается, что система состоит из m подсистем с фазовыми координатами xj = (xj,1 , xj,2 )> ∈ Rn , где n = n1 + n2 , xj,1 ∈ Rn1 , xj,2 ∈ Rn2 , j = 1, 2, . . . , m, и описывается соотношениями xj [k] = Aj [k] xj [k−1] + B j [k] uj [k], k = 1, 2, . . . , N, j = 1, 2, . . . , m; (1) xj [0] ∈ P0j ; uj [k] ∈ Rj [k], k = 1, 2, . . . , N ; j = 1, 2, . . . , m; (2) kxj,1 [k] − xl,1 [k]k∞ ≥ µ, k = 0, 1, . . . , N, j, l = 1, 2, . . . , m, l > j; (3) kxj,1 [k] − xl,1 [k]k∞ ≤ µ, k = 0, 1, . . . , N, j, l = 1, 2, . . . , m, l > j. (4) Здесь n1 > 0, n2 ≥ 0; Aj [k] и B j [k] — известные матрицы, det Aj [k] 6= 0; uj [k] ∈ Rnu — управления, стеснен- ные геометрическими ограничениями из (2); множества P0j ⊂ Rn и Rj [k] ⊂ Rnu считаем параллелепипеда- ми. Формализация условий группового движения производится с помощью (3), (4), где 0 < µ < µ ≤ ∞ — заданные значения, и кубической нормы kxk∞ (отличающейся от используемой в [17] евклидовой нормы). Условия (3) (нижние ограничения) означают, что траектории подсистем попарно не должны сближаться ближе, чем на положительное значение µ в смысле первой группы координат xj,1 , а условия (4) (верхние ограничения) — что, кроме того, не должны расходиться более, чем на µ, если это число принято конечным. Если, например, рассматриваемая система получена дискретизацией по схеме Эйлера из дифференциаль- ной системы, описывающей механическое движение системы материальных точек, можно считать, что n1 = n2 , и интерпретировать xj,1 как координаты, а xj,2 — как скорости этих точек. Этим объясняется название доклада, хотя рассматриваемая система не обязательно должна иметь механический смысл. 51 > > Введем расширенные векторы фазового состояния z = (x1 , x2 , . . . , xm > )> ∈ Rmn и управления u = 1> > (u , u2 , . . . , um > )> ∈ Rmnu полной системы. Множеством достижимости (МД) Z[i] системы (1)–(4) в момент i ∈ {0, 1, . . . , N } называем множество всех тех точек z ∈ Rmn , для каждой из которых существуют такие z[0] и u[·], удовлетворяющие (2), что порождаемое ими в силу (1) решение z[·] будет удовлетворять условиям z[i] = z и (3)–(4) для k = 0, 1, . . . , i. Многозначную функцию Z[k], k = 0, 1, . . . , N , называем по аналогии с [2] трубкой достижимости или трубкой траекторий Z[·]. В [20] были исследованы некоторые свойства вышеупомянутых множеств достижимости и дано описание трубок достижимости в виде ветвящихся трубок. В настоящей работе рассматривается задача нахождения внутренних полиэдральных оценок для них с использованием “элементарных” параллелепипедозначных оценок для множеств. Напомним некоторые упомянутые выше и используемые далее понятия. Называем множество P ⊂ Rn внутренней оценкой множества Q ⊂ Rn , если Q ⊆ P . Pn Параллелепипедом P(p, P , π) ⊂ Rn называем множество P = P(p, P , π) = {x | x = p + i=1 pi πi ξi , |ξi | ≤ 1}, где p ∈ Rn ; P = {pi } ∈ Rn×n — неособая матрица (det P 6= 0) со столбцами pi , kpi k2 = 1 (приведенное условие нормировки в терминах евклидовой нормы может быть опущено с целью упрощения формул); π ∈ Rn , π ≥ 0 (векторные неравенства понимаем покомпонентно). Можно сказать, что p — центр параллелепипеда, P — матрица ориентации, pi — направления, πi — величины его “полуосей”. Параллелотопом P[p, P̄ ] ⊂ Rn называем множество P = P[p, P̄ ] = {x| x = p + P̄ ζ, kζk∞ ≤ 1}, где p ∈ Rn , а матрица P̄ = {p̄i } ∈ Rn×r может быть особой и необязательно квадратной (r ≤ n). Называем параллелотоп P невырожденным, если r = n и det P̄ 6= 0. Полосой или m-полосой S = S(c, S, σ, m) ⊂ Rn называем пересечение m гиперполос Σi , где 1 ≤ m ≤ n T m (здесь m не связано с числом подсистем в (1)): S = S(c, S, σ, m) = Σi , где Σi = Σ(ci , si , σi ) = {x | |x> si − i=1 ci | ≤ σi }. Здесь c ∈ Rm ; S = {si } ∈ Rn×m — матрица ранга m со столбцами si ∈ Rn ; σ ∈ Rm , σ ≥ 0. Векторы ±si определяют нормали к гиперплоскостям, ограничивающим гиперполосу Σi . Каждый параллелепипед — это параллелотоп: P(p, P , π) = P[p, P̄ ], P̄ = P diag {πi }; каждый невырож- денный параллелотоп — это параллелепипед с P = P̄ diag {kp̄i k−1 i 2 }, πi = kp̄ k2 или, иначе, с P = P̄ , π = e, > где e = (1, 1, . . . , 1) . Известно также, что при m = n полоса превращается в параллелепипед (см. формулы в [8]). Везде ниже считаем, что выполнено следующее предположение. Предположение 1. Множества P0j = P(pj0 , P0j , π0j ) = P[pj0 , P̄0j ] ⊂ Rn и множества Rj [k] = P(rj [k], Rj [k], ρj [k]) = P[rj [k], R̄j [k]] ⊂ Rnu в (2) — это параллелепипеды; все матрицы Aj [k] неособые. В следующем разделе приводится описание множеств достижимости и трубок достижимости рассматри- ваемых систем. Четвертый раздел посвящен алгоритмам построения искомых внутренних полиэдральных оценок. 3 Описание множеств достижимости и трубок достижимости Несложно видеть, что полная система может быть записана в виде: z[k] = A[k] z[k−1] + B[k] u[k], k = 1, 2, . . . , N ; z[0] ∈ P0 ; u[k] ∈ R[k]; k = 1, 2, . . . , N ; (5) z[k] ∈ Y, k = 0, 1, . . . , N. Здесь A[k] = diag {Aj [k]} и B[k] = diag {B j [k]} — блочные матрицы, составленные из матриц, фигуриру- ющих в (1); параллелепипеды P0 = P[p0 , P̄0 ] и R[k] = P[r[k], R̄[k]] (записанные для простоты обозначе- ний в виде параллелотопов) определяются множествами P0j = P[pj0 , P̄0j ] и Rj [k] = P[rj [k], R̄j [k]] из (2): > > > > j p0 = (p10 , p20 , . . . , pm 1 > 2 > m > > j 0 ) , P̄0 = diag {P̄0 }, r[k] = (r [k] , r [k] , . . . , r [k] ) , R̄[k] = diag {R̄ [k]}; множе- ство Y определяется нижними и верхними ограничениями из условий группового движения: T Y = Y̌ Ŷ; > > Y̌ = {z = (x1 , x2 , . . . , xm > )> | kxj,1 − xl,1 k∞ ≥ µ, j, l = 1, 2, . . . , m, l > j}; (6) 1> 2> Ŷ = {z = (x ,x , . . . , xm > )> | kxj,1 − xl,1 k∞ ≤ µ, j, l = 1, 2, . . . , m, l > j} 52 и задает фазовые ограничения в системе (5). Известно [3, § 18], что МД системы типа (5) удовлетворяют следующим рекуррентным соотношениям: Z (0) [0] = P0 ; Z (0) [k] = A[k] Z[k−1] + B[k] R[k], k = 1, 2, . . . , N ; (7) T Z[k] = Z (0) [k] Y, k = 0, 1, . . . , N, которые включают операции с множествами: линейное преобразование, сумму Минковского и пересечение. Точное построение МД для систем большой размерности даже в случае выпуклых фазовых ограничений может быть достаточно затруднительно. В рассматриваемом же случае дополнительная трудность связана с невыпуклостью множеств Y̌ и Y. Это приводит к тому, что МД могут быть, вообще говоря, не только невыпуклыми, но и несвязными. Рассмотрим, что представляет собой пересечение выпуклого множества Z (0) с Y. Несложно заметить, что невыпуклое множество Y̌ представимо в следующем виде: T S T S Y̌ = ( Y̌ijl ) = ( Y̌iγ ), Y̌ijl = {z | |xj,1 l,1 i − xi | ≥ µ}, j,l=1,2,...,m, l>j 1≤i≤n1 γ=1,2,...,J 1≤i≤n1 S jl S γ Y̌ijl = Y̌(i,0) jl γ Y̌(i,1) = Y̌(i,0) Y̌(i,1) = Y̌iγ , Y̌(i,0) jl ={z | xj,1 l,1 γ jl j,1 l,1 i −xi ≤−µ}=Y̌(i,0) , Y̌(i,1) ={z | xi −xi ≥µ}=Y̌(i,1) , γ (8) пересечение взято по всем различным парам подсистем, определяемых индексами j и l. Каждая из ком- понент пересечения есть объединение n1 множеств Y̌ijl = Y̌iγ , каждое из которых есть в свою очередь jl γ jl γ объединение двух полупространств Y̌(i,0) = Y̌(i,0) и Y̌(i,1) = Y̌(i,1) в Rmn с обратно коллинеарными норма- лями. Здесь и далее для упрощения обозначений вместо двух индексов j и l используем один индекс γ, которым нумеруем различные пары подсистем; число таких пар равно J = m(m − 1)/2 штук, где m — общее число подсистем в (1). Аналогичную структуру имеет множество Y: T S T S Y= ( Yijl ) = ( Yiγ ), Yijl = {z | |xj,1 l,1 i − xi | ≥ µ, kx j,1 − xl,1 k∞ ≤ µ}, j,l=1,2,...,m, l>j 1≤i≤n1 γ=1,2,...,J 1≤i≤n1 S jl S γ Yijl = Y(i,0) jl γ Y(i,1) = Y(i,0) Y(i,1) = Yiγ . (9) γ jl γ jl При этом можно считать, что пары множеств Y(i,0) = Y(i,0) и Y(i,1) = Y(i,1) в (9) с одинаковыми индексами γ и i — это пары непересекающихся n3 -полос, имеющих одинаковые параметры S и σ (будем называть такие полосы парными), где n3 = n1 при µ < ∞ и n3 = 1 при µ = ∞ (т.е. когда условия (4) отсутству- ют). Несложно заметить, что матрицы S имеют столбцы типа (0, . . . , 0, 1, 0, . . . , 0, −1, 0, . . . , 0)> с двумя γ ненулевыми элементами. Формулы для Y(i,ω) , ω ∈ {0, 1}, приведены в [20, Замечание 1]. Заметим также, что нетрудно распространить результаты статьи на случай, где неравенства в условиях группового движения (3), (4) заменены более общими типа |xj,1 l,1 j l i − xi | ≥ µi , i = 1, 2, . . . , n1 ; |xi − xi | ≤ µi , i = 1, 2, . . . , n, где 0 < µi < µi ≤ ∞. Отличие будет в конкретных формулах для полос, определяющих Y. 8 2 0 x −8 −8 0 8 x1 Рис. 1: Пример множеств Z (0) = P0 , Y̌, Ŷ и Y для случая m = 2, n1 = 1, n2 = 0 На рис. 1 дан пример для случая m = 2, n1 = 1, n2 = 0, показывающий, что может получиться в Rm(n1 +n2 ) при пересечении параллелепипеда Z (0) = P0 с множествами Y̌ и Y. Здесь Y̌ состоит из двух 53 непересекающихся полупространств, ограниченных сплошными линиями, Ŷ — это гиперполоса между штрихпунктирными T линиями, Y состоит из двух непересекающихся параллельных гиперполос. Множе- ство P0 Y̌ состоит из двух политопов, а P0 ∩ Y тоже несвязно и состоит из двух выпуклых компонент — частей этих политопов. S T TВ общем S T случае в силу взаимной дистрибутивности операций объединения и пересечения ((A B) C = (A C) (B C)) справедлива следующая лемма [20]. T Лемма 1. Множество Z = Z (0) Y представимо в виде объединения (2n1 )J компонент: S S ³ T¡ 1 T 2 T T J ¢´ Z= Z (0) Y(i1 ,ω1 ) Y(i 2 ,ω2 ) . . . Y (iJ ,ωJ ) 1≤i1 ,...,iJ ≤n1 0≤ω1 ,...,ωJ ≤1 ³ ´ (10) S S T¡ T J γ ¢ = Z (0) Y(i γ ,ωγ ) , 1≤i1 ,...,iJ ≤n1 0≤ω1 ,...,ωJ ≤1 γ=1 где J = m(m − 1)/2, а объединения берутся по всевозможным индексам i1 , i2 , . . . iJ , изменяющимся от 1 до n1 , и всевозможным индексам ω1 , ω2 , . . . ωJ , принимающим значения 0 и 1. TJ γ Замечание 1. Обозначим через Wh множества Wh = γ=1 Y(i γ ωγ ) , которые мы “занумеровали” с помощью векторов типа h = ((−1)ω1 +1 i1 , (−1)ω2 +1 i2 , . . . , (−1)ωJ +1 iJ )> ∈ RJ , (11) где iγ ∈ {1, 2, . . . , n1 }, ωγ ∈ {0, 1}, γ = 1, 2, . . . , J. Обозначим символом H множество всевозможных век- торов h такого типа; его мощность равна M (H) = K = (2n1 )J . Тогда соотношения (10) можно записать в виде S ¡ (0) T ¢ Z= Z Wh . (12) h∈H С учетом рекуррентных соотношений (7), леммы 1 и замечания 1 несложно видеть, что трубка дости- жимости Z[·] может быть представлена в виде ветвящейся трубки с ветвями, которые обозначаем через Z[·; H[·]] (здесь можно заметить некоторое сходство с трубками траекторий гибридных систем [21], [2, разд. 11.2]). Здесь матрица H[k] = [h[0]; h[1]; . . . ; h[k]] = {h[κ]}kκ=0 ∈ RJ×(k+1) увеличивающейся от шага к шагу размерности определяет “историю” построения ветви до шага k. Ее столбцы h[κ] (0 ≤ κ ≤ k) типа γ (11) хранят информацию о значениях индексов i[κ] и ω[κ] тех полос Y(i,ω) , которые определяют множество Y в (9) и используются при построении ветви Z[·; H[·]] на шаге κ в соответствии с приведенными ниже формулами (13). Утверждение 1. (См. [20]). Трубка достижимости Z[·] системы (1)–(4) (иными словами, системы (5)–(6)) может быть представлена в виде ветвящейся трубки с ветвями Z[·; H[·]], сечения которых Z[k; H[k]] в моменты k ∈ {0, 1, . . . , N } являются выпуклыми политопами. Каждая ветвь Z[·; H[·]] удо- влетворяет следующим соотношениям: T TJ γ T Z[0; H[0]] = P0 ( γ=1 Y(i γ [0],ωγ [0]) ) = P0 Wh[0] ; H[0] = [h[0]]; Z (0) [k; H[k−1]] = A[k] Z[k−1; H[k−1]] + B[k] R[k], k = 1, 2, . . . , N ; H[k] = [H[k−1]; h[k]] = [h[0]; h[1]; . . . ; h[k]], k = 1, 2, . . . , N ; ω1 [k]+1 ω2 [k]+1 (13) h[k] = ((−1) i1 [k], (−1) i2 [k], . . . , (−1)ωJ [k]+1 iJ [k])> , k = 0, 1, . . . , N ; T Z[k; H[k]] = Z (0) [k; H[k−1]] Wh[k] , k = 1, 2, . . . , N ; TJ γ Wh[k] = γ=1 Y(i γ [k],ωγ [k]) , k = 0, 1, . . . , N. При этом множества S достижимости Z[k] представимы в виде объединения вышеупомянутых сече- ний: Z[k] = H[k] Z[k; H[k]] (назовем эти сечения Z[k; H[k]] компонентами множеств достижимости), верхняя оценка числа которых дается большим числом L[k] = K k+1 , где K = (2n1 )J , J = m(m − 1)/2. Таким образом, для каждой ветви Z[·; H[·]] на каждом k-м шаге сначала производятся операции ли- нейного преобразования сечения и суммирования множеств, затем выбирается J штук n3 -полос из тех, которые определяют множество Y, конструируется соответствующее множество Wh[k] и производится опе- рация пересечения с ним. Индексы, определяющие эти полосы, запоминаются в последнем столбце h[k] матрицы H[k]. 54 Описание ветвящейся трубки достижимости Z[·] в форме с ветвями Z[·; H[·]] с использованием матриц H[k] отражает суть конструкций, но может быть не очень эффективно при компьютерной реализации. Заметим, что трубка Z[·] может быть представлена как позиционное дерево, а именно, как K-арное дерево, так что степень каждого узла не превышает K (см., например, [22, часть VIII, разд. Б.5.3]). Верхняя оценка числа получающихся таким образом ветвей и, соответственно, компонент МД дается очень большим числом L[k] = K k+1 = (2n1 )J(k+1) . При некоторых предположениях (типа приведенного ниже предположения 2) эту оценку можно немного уменьшить. Предположение 2. Известна следующая оценка: kxj [k]k∞ ≤ Cx , j = 1, 2, . . . , m, k = 0, 1, . . . , N , и выполняется следующее неравенство: CA · Cx + Cu < µ, где CA = max max kAj [k] − Ik, Cu = 1≤j≤m 1≤k≤N max max maxj kB j [k]uj [k]k∞ . 1≤j≤m 1≤k≤N u [k]∈R [k] j Здесь величина Cx может быть вычислена, например, как некоторая оценка сверху для точек множеств достижимости подсистем (1)–(2) без учета условий (3)–(4). Заметим, что предположение 2 выполняется, в частности, для многошаговых систем, полученных из дифференциальных систем с помощью аппроксимаций Эйлера при достаточно малом шаге дискретизации. Элементарными оценками несложно получить, что при указанных условиях траектория не может за один шаг перескочить с одной полосы на другую, парную ей, полосу, и имеет место следующее утверждение. Следствие 1. (См. [20]). В условиях предположения 2 множество Z[k; H[k]] получается пустым для такой “истории” H[k], у которой хотя бы для одного γ ∈ {1, 2, . . . , J} выполняются условия: iγ [k] = iγ [k−1], но ωγ [k] 6= ωγ [k−1]. Следствие 2. В условиях предположения 2 число непустых компонент множества достижимости Z[k; H[k]] на шаге k ≥ 1 не превосходит числа L0 [k] = (2n1 )J (K 0 )k , где K 0 = (2n1 −1)J . В частности, если m = 2 и n1 = 1, то L0 [k] = 2, т.е. получается не более двух ветвей Z[·; H[·]]. В таблицах 1 и 2 приведены сводки значений некоторых величин, упомянутых выше, для ряда значений m, n1 и n2 . В первой таблице для каждого их приведенных величин m указаны соответствующее значение J (число различных пар объектов/подсистем), размерность nz = m · (n1 +n2 ) расширенного пространства Rm(n1 +n2 ) и число J · n3 , равное числу гиперполос, определяющих множества типа Wh , пересечение с которыми требуется проводить на каждом k-м шаге. Заметим, что при увеличении m, начиная с некоторого значения это число гиперполос начинает превосходить размерность расширенного пространства, что еще более осложняет задачу построения множеств достижимости. Во второй таблице указаны соответствующее значения K и K 0 из следствия 2. Таблица 1: Сводка некоторых значений величин m, J, nz и числа гиперполос в Wh Число пар Размерность расширенного пространства Число гиперполос в Wh (= J · n3 ) объектов nz = m · (n1 +n2 ) m nz = 2mn1 nz = mn1 Нет внеш. огранич. Есть внеш. огранич.: m(m−1) J= при n2 = n1 при n2 = 0 (т.е. µ = ∞): n3 = n1 2 n1 =1 n1 =2 n1 =3 n1 =1 n1 =2 n1 =3 n3 = 1 n1 =1 n1 =2 n1 =3 2 1 4 8 12 2 4 6 1 1 2 3 3 3 6 12 18 3 6 9 3 3 6 9 4 6 8 16 24 4 8 12 6 6 12 18 5 10 10 20 30 5 10 15 10 10 20 30 6 15 12 24 36 6 12 18 15 15 30 45 4 Внутренние полиэдральные оценки множеств достижимости Рассмотрим вопрос о построении внутренних полиэдральных оценок множеств достижимости для случая группового движения. В силу (7) внутренние оценки для МД Z[k] будут найдены, если построим элементарные полиэдраль- ные оценки для результатов фигурирующих в (7) операций с множествами. В [9–11] указаны способы постороения ряда элементарных оценок. Кратко напомним их. 55 Таблица 2: Сводка некоторых значений величин m, J, K и K 0 K = (2n1 )J K 0 = (2n1 − 1)J m(m−1) m J= n1 = 1 n1 = 2 n1 = 3 n1 = 1 n1 = 2 n1 = 3 2 (K = 2J ) (K = 4J ) (K = 6J ) (K 0 = 1J ) (K 0 = 3J ) (K 0 = 5J ) 2 1 2 4 6 1 3 5 3 3 8 64 216 1 27 125 4 6 64 4096 46656 1 729 15625 5 10 1024 > 106 1 ≈ 104 Пусть P k = P[pk , P̄ k ], k = 1, 2, P̄ 1 ∈ Rn×n , P̄ 2 ∈ Rn×r . Внутренние параллелотопозначные оценки для Q = P 1 + P 2 могут быть найдены [9] в виде P − 1 2 1 2 Γ (Q) = P[p +p , P̄ +P̄ Γ], где Γ ∈ G r×n = {Γ = {γij } ∈ r×n R | kΓk ≤ 1}. Матрицы Γ можно рассматривать как параметры, определяющие целые семейства оценок. TΥ Пусть Q — это ограниченный политоп, определяемый пересечением Υ ≥ n + 1 гиперполос: Q = j=1 Σj , Σj = Σ(cj , sj , σj ) = {x| | (x> sj − cj | ≤ σj }, и пусть Q имеет непустую внутренность int Q. Внутренние параллелепипедозначные оценки P − p− ,P − (Q) для Q с произвольным фиксированным центром p − ∈ int Q − n×n и неособой матрицей ориентации P ∈ R могут быть вычислены по явным формулам [10, 11]; p− и P − выступают в роли параметров оценок. Один из возможных способов нахождения центра p− при вы- бранной матрице ориентации P − — это решение задачи максимизации объема оценки, т.е. нахождение: p− ∈ Argmax {vol P − − p− ,P − (Q) | p ∈ Q}. Для численного решения может быть использован, например, сим- плексный метод Нелдера и Мида [23]. Некоторые примеры построения полиэдральных оценок трубок траекторий (для систем без условий группового движения) с использованием элементарных оценок приведены в [8–11]. Заменяя в рекуррентных соотношениях для нахождения МД результаты операций с множествами со- ответствующими элементарными оценками, рассмотрим ветвящуюся полиэдральную трубку с ветвями P − [·; H[·]], где сечения P − [k; H[k]] каждой ветви P − [·; H[·]] удовлетворяют следующим соотношениям: ( P0 , если P0 ⊆ Wh[0] , − P [0; H[0]] = T H[0] = h[0]; P−p− [H[0]],P − [H[0]] (P0 Wh[0] ) в противном случае; P (0)− [k; H[k−1]] = P − − Γ[H[k−1]] (A[k] P [k−1; H[k−1]] + B[k] R[k]), k = 1, 2, . . . , N ; H[k] = [H[k−1]; h[k]] = [h[0]; h[1]; . . . ; h[k]], k = 1, 2, . . . , N ; (14) ( − P (0)− [k; H[k−1]], если P (0)− [k; H[k−1]] ⊆ Wh[k] , P [k; H[k]] = T k = 1, 2, . . . , N ; P−p− [H[k]],P − [H[k]] (P (0)− [k; H[k−1]] Wh[k] ) в противном случае, TJ γ Wh[k] = γ=1 Y(i γ [k],ωγ [k]) , k = 0, 1, . . . , N. Здесь матрица H[k] = [h[0]; h[1]; . . . ; h[k]] = {h[κ]}kκ=0 ∈ RJ×(k+1) определяет “историю” построения вет- ви аналогично утверждению 1. Матричные функции Γ[·] и P − [·] и векторные функции p− [·], удовле- творяющие следующим условиям: Γ[H[k]] ∈ Rmnu ×mn , kΓ[H[k]]k ≤ 1, k = 0, 1, . T . . , N −1; P − [H[k]] ∈ mn×mn R , det P [H[k]] 6= 0, k = 0, 1, . . . , N ; p [H[k]] ∈ R ; p [H[0]] ∈ int (P0 Wh[0] ); p− [H[k]] ∈ − − mn − (0)− T int (P [k; H[k−1]] Wh[k] ), k = 1, 2, . . . , N , представляют собой допустимые параметры оценок. Утверждение 2. Если ветвящаяся полиэдральная трубка P − [·] такова, что ее ветви P − [·; H[·]] удо- влетворяют рекуррентным соотношениям (14) при каких-либо допустимых параметрах Γ[H[·]], P − [H[·]] и p− [H[·]], то она является внутренней оценкой для трубки достижимости Z[·] системы (5)–(6). При этом внутренними оценками для множеств достижимости Z[k] служат следующие объединения па- раллелепипедов P − [k; H[k]] — сечений ветвей: S − H[k] P [k; H[k]] ⊆ Z[k], k = 0, 1, . . . , N. (15) Верхняя оценка числа непустых множеств в объединении (15) в момент k дается тем же самым чис- лом L[k] = K k+1 , что и в утверждении 1. (Для упрощения обозначений не исключаем ситуацию, при которой какие-либо множества P − [k; H[k]] в (14), (15) могут быть пусты). 56 При построении полиэдральных трубок по формулам (14) в качестве эвристического способа выбора допустимых параметров можно предложить брать матрицы ориентации P − [H[k]] такими же, как матрицы ориентации параллелелотопов P (0)− [k; H[k−1]] (если эти параллелотопы получились невырожденными), и находить векторы p− [H[k]] путем максимизации объема оценки P − [k; H[k]] при выбранной матрице ори- ентации P − [H[k]]. Если рассматриваемые многошаговые системы получены из дифференциальных систем с помощью аппроксимаций Эйлера при достаточно малом шаге дискретизации, то матричные параметры Γ[H[k−1]] можно находить, используя формулы из [9, разд. 6]. Заметим, что описанные рекомендации, вообще говоря, не гарантируют непустоту сечений всех ветвей вычисляемой полиэдральной трубки. Опять получается, фактически, K-арное дерево, где K = (2n1 )J , J = m(m − 1)/2. Вычисление всех ветвей полного K-арного дерева [22, часть VIII, разд. Б.5.3] может быть практически невозможно для m > 2, n1 > 1 и не очень малых значений высоты дерева (равной k + 1 на шаге k). Предлагается строить не всю систему ветвей P − [·; H[·]], а зафиксировать некоторое число L и строить не более L внутренних трубок (ветвей P − [·; H[·]]), оставив при k = 0 не более L параллелотопов P − [0; H[0]], а затем последовательно оставляя на каждом k-м шаге из потенциально возможных L · K трубок не более L штук, а именно, те, которые имеют большие (по сравнению с другими) объемы сечений — параллелепипедов P − [k; H[k]]. При этом строится система L параллелепипедозначных трубок P (α)− [·], α = 1, 2, . . . , L, вычисление которых можно описать соотношениями типа использовавшихся в [10]: ( (α)− P (α),(0)− [k], если P (α),(0)− [k] ⊆ Wh(α) [k] , P [k] = T k = 0, 1, . . . , N, α = 1, 2, . . . , L; P−p(α)− [k]],P (α)− [k] (P (α),(0)− [k] Wh(α) [k] ), в противном случае, P (α),(0)− [k] = P −Γ(α) [k−1] (A[k] P (α)− [k] + B[k] R[k]), k = 1, 2, . . . , N, P (α),(0)− [0] = P0 , α = 1, 2, . . . , L; TJ Wh(α) [k] = γ=1 Y γ (α) (α) , k = 0, 1, . . . , N, α = 1, 2, . . . , L. (iγ [k],ωγ [k]) (16) Но при этом вычисление всех трубок P (α)− [·], α = 1, 2, . . . , L, вообще говоря, взаимосвязано, поскольку выбор значений h(α) [k] для разных трубок на k-м шаге согласовывается. Опишем его подробнее, считая для простоты обозначений, что L ≤ K. Введем следующие обозначения. Пусть f (h) — функция аргумента h, принимающего значения из неко- торого непустого конечного множества H, которое имеет мощность (число элементов) M (H) = K ≥ L ≥ 1. Записью H̃ = argmax(L) {f (h) | h ∈ H} будем обозначать процедуру нахождения множества H̃ = {h̃1 , h̃2 , . . . h̃L } таких L элементов множества H, которые дают наибольшие, по сравнению с другими, значения функции f . А именно, при K > L упорядочим все значения f (h) при h ∈ H по величине: f (hβ1 ) ≥ f (hβ2 ) ≥ . . . ≥ f (hβL ) ≥ f (hβL+1 ) ≥ . . . ≥ f (hβK ). Если выполнено строгое неравенство f (hβL ) > f (hβL+1 ), то полагаем H̃ = argmax(L) {f (h) | h ∈ H} = {h̃1 , h̃2 , . . . h̃L } = {hβ1 , hβ2 , . . . hβL }; множество H̃ определяется однозначно. Если же f (hβL ) = f (hβL+1 ), то либо просто полагаем H̃ = {hβ1 , hβ2 , . . . hβL }, либо находим индексы K1 и K2 , такие что 1 ≤ K1 ≤ L < K2 ≤ K и f (hβK1 ) = f (hβK1 +1 ) = . . . = f (hβK2 ), причем K1 таков, что либо K1 = 1, либо f (hβK1 −1 ) > f (hβK1 ), а K2 таков, что либо K2 = K, либо f (hβK2 ) > f (hβK2 +1 ); из векторов hβK1 , hβK1 +1 , . . . hβK2 случайным образом выбираем L − (K1 − 1) штук и вносим их в H̃ после внесения туда hβ1 , hβ2 , . . . hβK1 −1 (которое выполняется при K1 > 1). При K = L просто полагаем H̃ = H. Далее в качестве H будем рассматривать множество всевозможных векторов h типа (11). Векторы h(α) [k], определяющие в (16) множества Wh(α) [k] , строятся следующим образом. При k = 0 задаем некотороеT множество H̃ = {h̃1 , h̃2 , . . . h̃L } ⊆ H; например, находим H̃ = argmax(L) {vol (P − p(α)− [0],P (α)− [0] (P0 Wh )) | h ∈ H}. Полагаем h(α) [0] = h̃α , α = 1, 2, . . . , L. (17) Для вычислений при k ≥ 1 задаем некоторые множества H̃(α) [k] ⊆ H с числом элементов M (H̃(α) [k]) = K̃[k], α = 1, 2, . . . , L, где 1 ≤ K̃[k] ≤ K. Можно рассмотреть, в частности, два следующих “крайних” случая. (1) Все H̃(α) [k] = H (α = 1, 2, . . . , L, k = 1, 2, . . . , N ), т.е. берутся наибольшими возможными. При этом K̃[k] = K; для больших значений m объем вычислений может оказаться слишком большим. (2) Все множества H̃(α) [k] (α = 1, 2, . . . , L, k = 1, 2, . . . , N ) — одноточечные, определяемые равенствами H̃ [k] = h̃α . При этом K̃[k] = 1 и в результате вычислений будут найдены “стационарно расположенные” (α) внутренние оценки компонент множеств достижимости, т.е. те оценки, которые лежат в зафиксированных вначале и не меняющихся от шага к шагу множествах. 57 На k-м шаге для перехода от P (α),(0)− [k] к P (α)− [k] T для каждого α ∈ {1, 2, . . . , L} строится K̃[k] штук множеств P (α,β)− [k] = P − p(α,β)− [k]],P (α,β)− [k] (P (α),(0)− [k] Whβ ), где индекс β ∈ {1, 2, . . . , K̃[k]} пробегает множество значений, позволяющих перебрать все элементы множества H̃(α) [k] возможных (допускаемых алгоритмом) на k-м шаге значений векторов hβ типа (11). Таким образом, на k-м шаге получаем L · K̃[k] параллелепипедозначных множеств P (α,β)− [k], α = 1, 2, . . . , L, β = 1, 2, . . . , K̃[k] (некоторые из них могут оказаться пустыми, но для простоты обозначений не будем это как-то выделять специальным образом). При K̃[k] = 1 полагаем P (α)− [k] = P (α,1)− [k], α = 1, 2, . . . , L, и переходим к следующему шагу. При K̃[k] > 1 находим argmax(L) {vol P (α,β)− [k] | α ∈ {1, 2, . . . , L}, β ∈ {1, 2, . . . , K̃[k]}} = {(α1∗ [k], β1∗ [k]), (α2 [k], β2∗ [k]), . . . (αL ∗ ∗ [k], βL∗ [k])}. И для простоты обозначений полагаем сечения искомых трубок P (α)− [·], α = 1, 2, . . . , L, и соответствующие им векторы h(α) [k] равными ∗ ∗ ∗ P (δ)− [k] = P (αδ [k],βδ [k]))− [k], δ = 1, 2, . . . , L; h(δ) [k] = hβδ [k] , δ = 1, 2, . . . , L. При этом, очевидно, может получиться, что для построенного таким образом сечения P (δ)− [k] имеем αδ∗ 6= δ (т.е. требуется смена нумерации трубок), и может быть, что не все значения α1∗ , α2∗ , . . . , αL ∗ различны (какие- то трубки “оборвались”, а какие-то “разветвились”). Поэтому, чтобы остаться в рамках описания трубок с помощью соотношений (16), надо, построив L сечений P (δ)− [k] на шаге k, скорректировать нумерацию начальных участков трубок P (δ)− [·]. А именно, можно считать, что ∗ ∗ P (δ)− [κ] = P (αδ [k])− [κ], κ = 0, 1, . . . , k − 1, h(δ) [κ] = h(αδ [k]) [κ], κ = 0, 1, . . . , k − 1, δ = 1, 2, . . . , L. Следствие 3. Пусть трубки P (α)− [·], α = 1, 2, . . . , L, построены описанным выше способом. Тогда внутренними оценками S для множеств достижимости Z[k] будут служить объединения параллелепи- педов P (α)− [k]: 0≤α≤L P (α)− [k] ⊆ Z[k], k = 0, 1, . . . , N . Заметим, что при численной реализации вычисления трубок P (α)− [·] на компьютере вместо работы с L массивами, содержащими описания вышеупомянутых трубок, эффективнее использовать указатели и структуры типа дерево (чтобы для двух трубок с одинаковыми начальными участками “историй” соответ- ствующие участки хранились в одном экземрляре). Ввиду большого объема требуемых вычислений может быть полезно проводить их на многопроцессор- ном вычислительном комплексе. В качестве простейшего способа распараллеливания (вообще говоря, не самого эффективного с точки зрения получения трубок с наибольшими объемами сечений) можно предло- жить простейшую “процессорную ферму”, когда “рабочие” процессоры рассчитывают отдельные системы внутренних параллелепипедозначных трубок (каждый — не более L штук аналогично тому, как описано выше), и операция argmax(L) применяется к трубкам, рассчитываемым только на этом процессоре. Численная верификация предложенных алгоритмов является предметом дальнейшей работы. Благодарности Работа выполнена при частичной финансовой поддержке РФФИ (Проект 15-01-02368a). Список литературы [1] A.B. Kurzhanski, I. Vályi. Ellipsoidal Calculus for Estimation and Control. Boston: Birkhäuser, 1997. [2] A.B. Kurzhanski, P. Varaiya. Dynamics and Control of Trajectory Tubes: Theory and Computation. (Systems & Control: Foundations & Applications, Book 85). Birkhäuser Basel, 2014. [3] F.L Chernousko. State Estimation for Dynamic Systems. Method of Ellipsoids. Moscow: Nauka, 1988. (in Russian) = Ф.Л. Черноусько. Оценивание фазового состояния динамических систем. Метод эллипсоидов. М.: Наука, 1988; or: F.L Chernousko. State Estimation for Dynamic Systems. Boca Raton: CRC Press, 1994. [4] M.I. Gusev. Application of penalty function method to computation of reachable sets for control systems with state constraints. AIP Conf. Proc., 1773:050003-1–050003-9, 2016; doi: 10.1063/1.4964973. [5] A.B. Kurzhanski, P. Varaiya. On synthesizing team target controls under obstacles and collision avoidance. J. Franklin Inst., 347(1):130–145, 2010. 58 [6] T.F. Filippova. Estimates of reachable sets of impulsive control problems with special nonlinearity. AIP Conf. Proc., 1773:100004-1—100004-10, 2016; doi: 10.1063/1.4964998. [7] T.F. Filippova, O.G. Matviychuk, E.K. Kostousova. Estimation techniques for uncertain dynamical systems with bilinear and quadratic nonlinearities. Dynamical Systems: Control and Stability. Proc. of the 13th International Conference on Dynamical Systems: Theory and Applications (DSTA 2015), Lodz: ARSA Druk i Reklama: 185–196, 2015. [8] E.K. Kostousova. State estimation for dynamic systems via parallelotopes: optimization and parallel computations. Optimiz. Methods & Software, 9(4):269–306, 1998. [9] E.K. Kostousova. Control synthesis via parallelotopes: optimization and parallel computations. Optimiz. Methods & Software, 14(4):267–310, 2001. [10] E.K. Kostousova. On internal polyhedral estimates of reachable sets of linear systems with state constraints. Algoritmy i programmnye sredstva parallelp1nykh vychislenii. Vol. 5. Ekaterinburg: Rossiiskaya Akademiya Nauk, Ural’skoe Otdelenie, Institut Matematiki i Mekhaniki: 167–187, 2001. (in Russian) = Е.К. Костоусова. О внутренних полиэдральных оценках множеств достижимости линейных систем с фазовыми ограничениями. Алгоритмы и программные средства параллельных вычислений: [сб.науч.тр.]. Екатеринбург, УрО РАН, 2001. Вып.5. С.167–187. [11] E.K. Kostousova. On control synthesis for uncertain dynamical discrete-time systems through polyhedral techniques. Discrete Contin. Dyn. Syst. 2015, Dynamical systems, differential equations and applications. 10th AIMS Conference. Suppl.: 723–732, 2015. [12] A.V. Lotov. A numerical method for constructing sets of attainability for linear controlled systems with phase constraints. USSR Comput. Math. Math. Phys., 15(1):63–74, 1975. = А.В. Лотов. Числен- ный метод построения множеств достижимости для линейных управляемых систем с фазовыми ограничениями. Журн. вычисл. математики и мат. физики, 15(1):67–78, 1975. [13] A.F. Shorikov. Minimax Estimation and Control in Discrete Dynamical Systems. Ekaterinburg: Izdat. Ural. Univ., 1997. (in Russian) = А.Ф. Шориков. Минимаксное оценивание и управление в дис- кретных динамических системах. Екатеринбург: Изд-во Урал. гос. ун-та, 1997. [14] I.Ya. Kats, A.B. Kurzhanski. Minimax multistep filtering in statically uncertain situations. Autom. Remote Control, 39:1643–1650, 1979. = И.Я. Кац, А.Б. Куржанский. Минимаксная многошаговая фильтрация в статистически неопределенных ситуациях. Автомат. и телемех., (11):79–87, 1978. [15] I.Ya. Kats, G.A. Timofeeva. A modified residual method in a statistically indefinite filtration problem. Autom. Remote Control, 55(2, Pt.1): 229–236, 1994. = И.Я. Кац, Г.А. Тимофеева. Модифицирован- ный метод невязки в статистически неопределенной задаче оценивания. Автомат. и телемех., (2):100–109, 1994. [16] B.I. Ananyev. Finitely approximable random sets and their evolution via differential equations. AIP Conf. Proc., 1789:040012-1—040012-9, 2016; doi: 10.1063/1.4968465. [17] A.B. Kurzhanski. The problem of control for multi-agent motion. Doclady Mathematics, 79(3):314–318, 2009. = А.Б. Куржанский. Задача управления групповым движением. Общие соотношения. Докл. РАН, 426(1):20–25, 2009. [18] V. Gazi, B. Fidan. Coordination and control of multi-agent dynamic systems: models and approaches. Swarm Robotics Ws, LNCS, 4433:71–102, 2007. [19] R. Olfati-Saber. Flocking for multi-agent dynamic systems: algorithms and theory. IEEE Trans. on Automatic Control, 51(3):401–420, 2006. [20] E.K. Kostousova. On state estimation for multi-agent motion: discrete-time systems. Bulletin of the South Ural State University. Series: Computational Mathematics and Software Engineering, 5(1), 13–23, 2016; doi: 10.14529/cmse160102. Also available at: http://vestnik.susu.ru/cmi/article/view/3888/4226 (accessed 25.12.2016). (in Russian) = Е.К. Костоусова. Об оценивании состояний многошаговых 59 систем при групповом движении. Вестник ЮУрГУ. Серия: Вычислительная математика и ин-форматика, 5(1):13–23, 2016. URL: http://vestnik.susu.ru/cmi/article/view/3888/4226 (Дата об- ращения 25.12.2016). [21] A.B. Kurzhanski, P.A. Tochilin. Weakly invariant sets of hybrid systems. Differential Equations, 44(11), 1585–1594, 2008. = А.Б. Куржанский, П.А. Точилин. Слабоинвариантные множества гибридных систем. Дифференциальные уравнения, 44(11):1523–1533, 2008. [22] T.H. Cormen, C.E. Leiserson, R.L. Rivest, C. Stein. Introduction to Algorithms (2nd ed.) Cambridge, MA: MIT Press, 2001. = Т.Х. Кормен, Ч.И. Лейзерсон, Р.Л. Ривест, К. Штайн. Алгоритмы: построение и анализ, 2-е изд. М.: Издательский дом “Вильямс”, 2005. [23] J.C. Lagarias, J.A. Reeds, M.H. Wright, P.E. Wright. Convergence properties of the Nelder-Mead simplex method in low dimensions, SIAM Journal of Optimization, 9(1):112–147, 1998. 60 On internal estimates for reachable sets of discrete-time systems describing multi-agent motion Elena K. Kostousova Krasovskii Institute of Mathematics and Mechanics (Ekaterinburg, Russia) Abstract. The reachability problem for linear discrete-time control systems which describe so called multi- agent motion is considered. Namely, we consider a finite set of control subsystems under a condition that the trajectories of the subsystems should be pairwise not very close to and not very far away from each other. The conditions of multi-agent motion lead to nonconvex state constraints of special type. As a result the reachable sets can be nonconvex and even disconnected. Corresponding reachable tubes may be presented in the form of branching trajectory tubes whose cross-sections are convex polytopes; the reachable sets may be presented in the form of unions of the mentioned cross-sections. In the paper, there are proposed some algorithms for constructing internal polyhedral estimates for the reachable sets in the form of unions of parallelepipeds, which are constructed on the base of primary parallelotope-valued estimates for results of set operations with parallelotopes/parallelepipeds and zones. Keywords: reachable sets, discrete-time systems, multi-agent motion, internal polyhedral estimates. 61