CC BY
65
2012
ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА Управление, вычислительная техника и информатика
№ 3(20)
ОБРАБОТКА ИНФОРМАЦИИ
УДК 519.676
Т.А. Аверина
АЛГОРИТМ АНАЛИЗА СИСТЕМ УПРАВЛЕНИЯ АНСАМБЛЕМ ТРАЕКТОРИЙ С УЧЕТОМ СЛУЧАЙНОГО ИЗМЕНЕНИЯ СТРУКТУРЫ И СКАЧКОВ
Рассматривается задача анализа нелинейных систем управления ансамблем траекторий с учетом случайного изменения структуры системы и скачков (случайных импульсных воздействий). Для решения рассматриваемой задачи предлагается алгоритм статистического моделирования.
Ключевые слова: система со случайной структурой, метод Монте-Карло.
Под системами со случайной (переменной) структурой понимаются динамические системы, поведение которых на случайных интервалах времени характеризуется различными структурами и описывается различными уравнениями [1]. К системам с переменной структурой относятся:
- системы поиска, захвата и сопровождения сигнала;
- системы, имеющие срыв управления или слежения;
- системы с периодическим случайным повторением процессов поиска, захвата, срыва управления;
- системы управления, в которых связи между функциональными элементами меняются в зависимости от состояния (системы с переменной структурой управляющего устройства);
- системы с возможными нарушениями, в которых возможен случайный отказ функционирования некоторых цепей в случайные моменты времени.
В настоящей работе рассматривается задача анализа многомерных нелинейных систем управления ансамблем траекторий [2, 3] при импульсных воздействиях и случайном изменении структуры. Для решения этой задачи применяется метод статистического моделирования, который позволяет оценивать вероятностные характеристики выходных процессов, в том числе и плотность вероятности.
1. Постановка задачи анализа нелинейных систем управления ансамблем траекторий
Рассмотрим процесс (у(0, s(t))T, где s(f) - дискретный случайный процесс с конечным множеством состояний {1, 2,...,£}; £ - число структур, а у(/) - ^-мерный случайный процесс, описываемый при условии 5(/) = I системой обыкновенных дифференциальных уравнений (ОДУ):
= /(1 )(1,У(0), Ж) = У0 к",
ш
где t e [t0, T]; f(l -1 (t, y) - вектор-функции размера d; l - номер структуры системы, множество Q ограничено, оно характеризует неопределенность задания начальных данных [3].
Вероятность перехода дискретного случайного процесса s(t) удовлетворяет следующим условиям [1]:
P{s(t + At) = r | s(t) = l, y(t) = y} = vlr (t, y)At + o(At),
P{s(t + At) = 11 s(t) = l, y(t) = y} = 1 — vii (t, y)At + o(At),
s(to) = So, l,r = 1,2,...,S, l ^ r ,
где функция Vlr (t, y):[to, T] x Rd ^ [0, го] называется интенсивностью перехода, при этом
vll y) = X vlr(t, y).
r=1*l
Это условие обеспечивает при любом фиксированном y e Rd отсутствие нескольких переключений процесса s(t) за малый интервал времени At. Предполагается, что в моменты переключения траектории процесса y(t) могут быть как непрерывными (случай точного восстановления реализаций), так и иметь разрывы [1]. В последнем случае дополнительно задаются законы, определяющие величину скачка в момент смены структуры. Например, в качестве такой характеристики удобно взять ц!1г (t, y | y) (условную плотность вероятности смены вектора состояния из положения y в положение y, при переходе из l-й структуры в r-ю в момент переключения t) или плотность вероятности ц!1г (t, 9) величины скачка 9 = y — y . При точном восстановлении реализаций
Vlr y1 у) = xVlr ^ y — у) = 5(y — у).
Систему, описываемую соотношениями (1) и (2), можно рассматривать как стохастическую мультиструктурную систему с распределенными переходами (при отсутствии случайных возмущений, действующих на объект управления).
Наиболее полной вероятностной характеристикой расширенного вектора состояния является упорядоченная совокупность взвешенных плотностей распределения p*(I )(t, y) вектора состояния (здесь и в дальнейшем * означает, что отсутствует нормировка), удовлетворяющих следующему условию:
X i p*(1)(t, y)dy = 1, t e [to, T]. (3)
l=l Rd
Взвешенные плотности распределения удовлетворяют системе обобщенных уравнений Фоккера - Планка - Колмогорова [1]. Вероятность того, что в момент времени t структура системы имеет номер l, т.е. s(t) = l, задается выражением
P(l)(t) = f P*Q)(t,y)dy, l = 1,2...,5 ,
Rd
а плотность вероятности p(l) (t, y) вектора состояния при условии s(t) = l связана с p*(I )(t, y) следующим соотношением:
р()(Г,у) = Р()(Г)р()(Г, у), I = 1,2,..., 5 .
Для безусловной плотности вероятности р($, у) вектора состояния выполняется
р&у) = х р*1^, у)
l=1
5
p(t, У) = X p(l) (t)P*(I) (t, y) •
или
l=1
Задача анализа систем управления ансамблем траекторий состоит в нахож-
*(i)
дении ненормированных плотностей распределения p к J(t, у) вектора состояния
по заданным функциям f(l)(t, у), интенсивностям vlr (t, у), законам поведения траектории в момент смены структуры (плотностям ц!1г (t, 9)) и ненормированным
плотностям распределения p*(l) (t, у), l, r = 1,2,..., S .
Наряду с нахождением функций p*(1 )(t, у) можно рассматривать задачу нахождения маргинальных плотностей вероятности и моментных характеристик вектора состояния, а также задачу определения вероятностных характеристик времени перехода из одной структуры в другую [1-5]. Некоторые вероятностные характеристики и их оценки приведены в 3-м разделе статьи.
2. Анализ систем управления ансамблем траекторий методом статистического моделирования
Опишем статистический алгоритм решения задачи анализа систем управления ансамблем траекторий. Статистический алгоритм должен в себя включать: численное решение ОДУ (1) некоторым численным методом [6], а также моделирование моментов смены структуры, номера новой структуры и величины скачка (при разрыве траекторий). В рассматриваемом случае распределение моментов смены структуры определяется интенсивностями переходов (2). Так как интенсивности переходов зависят от вектора состояния, то моделирование моментов смены структуры будет осуществляться по методу максимального сечения [7]. Применение этого метода требует выполнения условий
vlr (t, у) <Vlr = const, i = 1,2,..., S, i Ф l, на всем интервале интегрирования [t0, T].
Алгоритм моделирования траекторий процесса (у(ї), s(t))T :
0) k: = 0; моделируем (у0,s0)T согласно заданным ненормированным плотностям распределения p*(l) (t, у), l = 1,2,..., S;
1) l := sk; моделируем возможный момент выхода из l-й структуры
tk+1 = tk +т , где т - случайная величина с плотностью p(x) = vl exp(-vl x),
vl = X Vli , (по формуле т = -^na, где a - равномерно распределенная на ин-
i=М vl
тервале (0,1) случайная величина, которая моделируется с помощью датчика псевдослучайных чисел, например RAND [8]); если tk+1 > T , то tk+1 := T ;
2) моделируем номер r (возможный номер новой структуры) с вероятностью pr = Vlr /Vl, гф-l, r = 1,2,...,S;
3) решаем уравнение (1) для l-й структуры на интервале [tk, tk+1] численным методом с шагом h и находим yk+1 - вектор состояния системы в момент tk+1, при этом шаг должен быть согласован с интенсивностью перехода, например h < 0,1/ vl;
4) k : = k + 1;
5) проверяем условие смены структуры: моделируем равномерно распределенную на интервале (0,1) случайную величину a1: если a1 <vlr (tk, yk)/ Vlr, то переходим к п. 6, иначе переходим к п. 7;
6) меняем номер структуры на r-й: sk := r ; если задано условие разрыва траектории при переходе из l-й структуры в r-ю, то моделируется новая случайная величина yk согласно плотности уlr (t, y | yk) или величина скачка 9 согласно плотности ylr (t, 9), тогда yk := yk +9 ;
7) если tk < T , то переходим к п. 1, иначе процесс моделирования завершается.
3. Вероятностные характеристики решения и их статистические оценки
При решении задач со случайной структурой возникает потребность в вычислении вероятностных характеристик решения [1]. Некоторые из них мы приведем ниже.
1. Вероятность l-й структуры (или вероятности состояния, l = 1,...,S)
S
P(l) (t) = P(s(t) = l), X P(l) (t) = 1.
l=1
2. Плотность условного распределения y(t) при s(t) = l
p(l)(t, y) = p(y(t)| s(t) = l), J p(l)(t, y)dy = 1.
Rn
3. Взвешенная плотность условного распределения y(t) при s(t) = l
p*(l)(t, y) = p(y(t), s(t) = l) = p(l)(t, y)P(l)(t) .
4. Плотность вероятности y(t)
S SS
p(t, y) = X p(t, y(t) I s(t) = l)P(s(t) = l) = X p(l)(t, y)P(l)(t) = X p*(l)(t, y).
l=1 l=1 l
5. Условные начальные моменты k-го порядка при s(t) = l (k = ^ +... + kd)
rokl)(t) = ®kll.,kd(t) = ••• y^d)(l)(t) = J yf1 ••• yddp(ly)dy.
Rd
6. Взвешенные условные начальные моменты k-го порядка при s(t) = l
®k(l) (t) = roit(i).,kd (t) = E(yf1 ••• yd )*(l) (t) = J yf1 ••• yd p*(l) (t, y)dy .
1 d Rd
7. Безусловные начальные моменты к-го порядка
»к (1) = »к1,...,кё (1) = Е(Уікі ••• Уй ХО = Р Уікі ••• Ук/ р(1, у)іїу =
кё
=Х »к' >« )р(' )(о=х»?')(»).
1=1 '=1
8. Условный центральный второй момент при s(t) = І
©(')(1) = [ (У - т(І) (1))(у - т(І) (1 ))Т р(І) (1, у)ёу, т(І )(1) = [ ур(І) (/, у)ёу .
9. Взвешенный условный центральный второй момент при 5(1) =І
©*(І) (1) = [ (у - т*(1) (1))(у - т*(І) (1 ))Т р*(І) (1, у )ёу, ті*1' (1) = [ ур*(І) (/, у )р .
10. Безусловный центральный второй момент
©(1) = Р (у - т(1))(у - т(1))Тр(1, у)ёу, т(1) = р ур(1,у)ёу.
Пусть Ж - объем выборки (общее число моделируемых траекторий); Ык) -число траекторий, находящихся в 1-м состоянии в момент времени ґк
(X Ык? = N для всех 1к); уіХ обозначает і-ю компоненту численного решения
І=1
ук в момент времени 1к, на п-й траектории в 1-й структуре; уіпк обозначает і-ю компоненту численного решения ук в момент времени ґк, на п-й траектории.
Нас будут интересовать следующие вероятностные характеристики решения в момент времени ґк:
1. Безусловные:
а) вектор математического ожидания т(ґк), компоненты которого обозначим
тгк = тг (1к) = Р у,р(ґк, у)ёу; і = 1,..., ё;
с компонентами
b) матрица вторых моментов Ф(/к) = [ ууТ р(/к, у)ёу
к
Ф (Ч) = %к; ^ 1 = ё;
c) матрица центральных вторых моментов ©(/к) = Ф(/к) - ш(/к)тТ (/к) с компонентами ©у (/к) = ®ук ; 1,1 = 1,..., ё.
Их статистические оценки имеют вид
N N N _ _
X X у1пк X- утку1ик X (у,пк - тк )(у}пк - т1к )
Шк = П=-, Фук = П=----------, ©ук = п=1--------------------------•
N N N-1
2. Условные:
ш(1^) , Ф(г)(^к), ©(')(^к) = Ф('Ч) - т(1 \/к )(т(1 ))Т (/к).
Их статистические оценки имеют вид
< < N — —
_ ХУПк _ ХУПк_ Xой -тк})(УМ -т£>)
™(0 = и=1 ф(1) = _И=1___________ 0(1) = П=1______________________
ік № ’ ук ыЦ) ’ ук ы(-) -1
3. Взвешенные:
да
( )
(/к), Ф*(1Ч), ©*(1)(/к) = Ф*(1Ч)-т*(1)(^кХт™/(/к).
Их статистические оценки имеют вид
ык1) ^ N
*(0 = _и=______ ф*(1) = _П=1__________ ©*(') = П=______________________________
,к N ’ ']к N ’ ']к N -1 .
Замечание [9]. Пусть х1,...,xN - выборка независимых одинаково распределенных случайных величин, с функцией распределения Е(х), математическим ожиданием т, центральным вторым моментом ц2 (дисперсией), центральным четвертым моментом ц4. Тогда статистические оценки математического ожидания х и центрального второго момента т2
_ 1 N
і N і N
= N 'X х , т2 = NXх'- Х)2
имеют следующие математические ожидания и дисперсии:
Е (х) = да, Е(х - да)2 =“^> Е (да2) = ^ ц2,
Е(да22) - (Е(да2))2 =^^1 - 2(^4 ~2^2) + ^М..
22 N Ж2 Ж3
Согласно этому замечанию, можно найти первые два момента рассматриваемых выше оценок вероятностных характеристик решения систем со случайной структурой.
4. Вероятности нахождения системы в 1-й структуре в момент времени оцениваются по формуле
N .
5. Гистограмма маргинальной плотности условного распределения р(1) ^к, yi)
строится следующим образом:
a) [а/, Ь/ ] - область изменения у/ ^к) равномерно разбивается на ч частей с
Ь — а
шагом к = —------; обозначим 5у -у-й интервал области,у = 1,...,ч;
q ■'
b) моделируем п = 1,...,Ж траекторий и в момент tk на каждой траектории для /-й компоненты решения, находящейся в 1-й структуре, считаем сколько у/ППк попало в у-й интервал; полученное значение обозначим N (к (X Ч =1N $ = жк));
с) в качестве приближенного значения для p(1 \tk, y ) на интервале 5берется
Nk
величина Pj =—.
NhP( )
6. Гистограмма маргинальной плотности распределения p(tk,y ) строится следующим образом:
a) [a, bi ] - область изменения y (tk ) равномерно разбивается на q частей с
b — a-
шагом h = —-----L ; обозначим 5 - j-й интервал области, j = 1,.. .,q;
q
b) моделируем n = 1,...,N траекторий и в момент tk на каждой траектории для i-й компоненты решения считаем сколько yink попало в j-й интервал; полученное
q
значение обозначим N jk ( Е N]k = N );
j=i
c) в качестве приближенного значения для p (tk, y ) - плотности распределе-
Njk
ния yi (tk ) на интервале 5j - берется величина Pj = •
При вычислении математического ожидания некоторого функционала от решения рассматриваемых задач построенным статистическим алгоритмом важной является проблема оптимального (согласованного) выбора параметров статистического алгоритма: шага численного метода h и размера выборки N. При вычислении гистограммы решения также возникает проблема оптимального выбора параметров статистического алгоритма: шага численного метода h, размера выборки N и шага гистограммы hg . Используя методику доказательства, рассмотренную в
[10], для решения задачи условной оптимизации статистического алгоритма, можно доказать, что минимум трудоемкости вычисления гистограммы достигается при ng » у-'', N » у-3, h » у17 p , где у - требуемая точность вычислений гистограммы в норме пространства L2([a,b]), ng - число узлов гистограммы, p - порядок сходимости выбранного численного метода решения ОДУ.
4. Расчет модельной задачи анализа системы стабилизации малого искусственного спутника
Актуальными являются разработки и применения малых искусственных спутников стандарта CubeSat (нано-спутников и пико-спутников). В качестве примеров можно привести международную систему Disaster Monitoring Constellation, позволяющую производить мониторинг катастроф по всему миру; норвежский спутник nCube, отслеживающий перемещение кораблей по территориальным водам Норвегии; проект AAUSAT, предназначенный для получения детальных изображений Земли; спутник UWE-1 (Вюрцбургский университет), позволяющий анализировать использование технологий TCP/IP для телеметрических и телеко-мандных данных с учетом проблем задержек и помех; японский проект XI, созданный для демонстрации и тестирования систем спутниковой платформы с ис-
пользованием готовых элементов, включая проверку аппаратуры спутника в условиях реального орбитального полета.
Возмущенное движение спутника, находящегося под действием гравитационного и управляющего моментов, в плоскости орбиты описывается следующими уравнениями:
d 9(т) = д(х), dq(T) = -3Q2p sin 9(т) cos(9(x)) + Lv(т), d т d т
где 9 - угол отклонения оси спутника по отношению к радиус-вектору центра масс, q - угловая скорость вращения вокруг центра масс, Q - угловая скорость обращения спутника по круговой орбите, L и р - константы, зависящие от конструкции спутника, v - управление [11]. В [12] показано, что при малых колебаниях спутника с учетом возможного отказа управляющего устройства его движение приближенно описывается уравнениями
= У2(0, = - y\(t) + (2 - k )u(t), k = 1,2,
где t = т / а , y1 = 9/ у, y2 = q / 5 . Числа a, y, 5 выбраны таким образом, что
—=1, 3«2Éai=1, a£=1, t € [0,1].
Y 5 5
Случай k = 1 соответствует режиму нормального функционирования, а случай k = 2 - режиму отказа, т.е. срыву стабилизации. Соотношения для определения оптимального управления получены в [12]. Управляющее воздействие
u(t) = -0,2(1 - tg1 - 2tg21)cos t - 0,2(3 + 3tg1 + 2tg21)sin t
в среднем обеспечивает минимальный расход энергии, а также стабилизацию спутника в момент времени t = 1, т.е. математические ожидания величин (1) и y2(1) равны нулю. Пусть начальные значения y10 и y20 являются независимыми случайными величинами, имеющими усеченное нормальное распределение с параметрами
m10 = -0, 3, D10 = 1 m20 = 0Д, D20 = 1.
Интенсивности отказа и восстановления заданы функциями v12(t) = 0,2 и v 21 (t) = 0,1, соответственно. В начальный момент времени система функционирует нормально с вероятностью P0(1) = 0,95 .
В качестве конкретного численного примера была рассмотрена задача нахождения вероятностей работы системы в режимах нормального функционирования и срыва стабилизации, маргинальных плотностей вероятности, математических ожиданий и вторых начальных моментов координат вектора состояния.
Задача была посчитана приведенным выше алгоритмом статистического моделирования. При решении задачи построенным алгоритмом использовался обобщенный одностадийный метод типа Розенброка [6] с шагом h = 0,01; число моделируемых траекторий N = 106; для построения гистограммы отрезок [-4, 4] изменения координат вектора состояния равномерно разбивался на 100 частей (ng = 100 ). Оценки функционалов от решения и гистограммы вычислялись одно-
временно. Время счета составило 3 мин (для расчетов использовался РС с процессором Intel Celeron 2ГГц и 256 МБ оперативной памяти). Погрешность оценок функционалов от решения составляет O(h). Выбранные параметры задачи (N = 106, h = 10-2, ng = 102) являются оптимальными параметрами для получения
наилучшей оценки гистограммы при выбранном шаге h = 10~2 и гарантируют погрешность вычисления гистограммы в норме пространства L2([-ro, да]) порядка O(h). Дальнейшее уменьшение шага гистограммы или шага численного метода точность вычисления гистограммы не увеличивает. Дополнительные расчеты проводились при h = 10-2 и ng = 5 • 102, а также при h = 10-3 и ng = 5 • 102. Погрешность гистограммы не уменьшается, хотя время вычислений увеличивается (при h = 10-3 в 10 раз).
Заключение
Рассмотренная математическая модель позволяет описывать сложные системы управления, подверженные импульсным воздействиям, имеющие различные режимы функционирования и неопределенность в задании начальных данных. Построенным алгоритмом статистического моделирования была посчитана модельная задача анализа системы стабилизации малого искусственного спутника, находящегося под действием гравитационного и управляющего моментов, с учетом возможного отказа управляющего устройства. Точность оценок вероятностных характеристик решения, полученных методом статистического моделирования, зависит от шага интегрирования h, числа моделируемых траекторий N и шага гистограммы hg . Проведенные численные расчеты показали, что построенный алгоритм позволяет с высокой точностью вычислять вероятностные характеристики решения.
ЛИТЕРАТУРА
1. Казаков И.Е., Артемьев В.М., Бухалев В.А. Анализ систем случайной структуры. М.: Физматлит, 1993. 272 с.
2. Овсянников Д.А. Математические методы управления пучками. Л.: Изд-во ЛГУ, 1980. 228 с.
3. Пантнлеев А.В., Руденко Е.А., Бортаковский А. С. Нелинейные системы управления: описание, анализ и синтез. М.: Вузовская книга, 2008. 312 с.
4. Аверина Т.А. Статистический алгоритм моделирования динамических систем с переменной структурой // Сиб. журн. вычисл. матем. 2002. Т. 5. № 1. С. 1-10.
5. Averina Т.А. Algorithm for statistical simulation of two types of random-structure systems // Russ. J. Numer. Anal. Modelling. 2001. V. 16. No. 6. P. 467-482.
6. Artemiev S.S., Averina Т.А. Numerical analysis systems of ordinary and stochastic differential equations. Utrecht: VSP, 1997. 176 p.
7. Михайлов Г.А., Аверина Т.А. Алгоритм «максимального сечения» в методе Монте-Карло // ДАН. 2009. Т. 428. № 2. С. 163-165.
8. Михайлов Г.А., Ермаков С.М. Курс статистического моделирования. М.: Наука, 1976.
320 с.
9. Крамер Г. Математические методы статистики. М.: Мир, 1975. 648 с.
10. Михайлов Г.А. Весовые методы Монте-Карло. Новосибирск: Изд-во СО РАН, 2000. 248 с.
11. Гурман В.И. Вырожденные задачи оптимального управления. М.: Наука, 1977. 250 с.
12. Rybakov K.A., Sotskova I.L., Yudin M.A. Synthesis of optimal control algorithms for small satellite subject to the possible failure of control unit // Theoretical questions of computers and software. M.: MIREA, 2006. P. 98-103.
Аверина Татьяна Александровна
Институт вычислительной математики и математической геофизики СО РАН,
Новосибирский государственный университет,
E-mail: [email protected] Поступила в редакцию 22 мая 2012 г.
Averina Tatyana A. (Institute of Computation Mathematics and Mathematical Geophysics SB RAS, Novosibirsk State University). Algorithm for analysis of ensemble paths control systems subjected to the random change of structure and jumps.
Keywords: systems with a random structure, Monte Carlo method.
We consider the problem of analysis of ensemble path control systems subjected to the random change of structure and jumps. If S is the number of structures in the system, s(t) is a discrete random process with a finite set of values {1, 2,.. .,S}, then state vector of the system y(t) is a n-dimensional random process described under the condition s(t) = l by the following ordinary differential equations:
^ = f (l\t, y(t)), y(t0) = У0 eQc Rn, t e [t0,T ] dt
where; f(l)(t,y) is the vector function of the dimension n; l is a number of structures of the system. The transition probability of the discrete random process s(t) satisfies the following conditions:
P{s(t + At) = r | s(t) = l,y(t) = y} = vlr (t,y)At + o(At),
P{s(t + At) = 11 s(t) = l,y(t) = y} =1 - vll(t,y)At + o(At),
S
s(t0) = so, l,r = 1,2,...,S, l *r; vii<xy) =Z vr(t,y).
r=1*l
Statistical algorithm was constructed for solving the problem of analysis of ensemble path control systems subjected to the random change of structure and jumps. Developed algorithm is applied to analysis of the satellite stabilization system.
