Численные методы построения области достижимости динамической системы Текст научной статьи по специальности «Математика»

СИСТЕМЫ УПРАВЛЕНИЯ

J

УДК 519.6

Е. М. Воронов, А. П. Карпенко, О. Г. Козлова, В. А. Ф е д и н

ЧИСЛЕННЫЕ МЕТОДЫ ПОСТРОЕНИЯ ОБЛАСТИ ДОСТИЖИМОСТИ ДИНАМИЧЕСКОЙ СИСТЕМЫ

Рассмотрена задача приближенного построения области достижимости динамической системы. Исследованы подходы к решению этой задачи на основе последовательных и параллельных методов мультифиниша, аппроксимации векторного поля соответствующей системы обыкновенных дифференциальных уравнений, а также на основе нейронных сетей. В качестве примера рассмотрено применение указанных методов для приближенного построения границы области достижимости летательного аппарата.

E-mail: apkarpenko@mail.ru

Ключевые слова: динамическая система, область достижимости, аппроксимация векторного поля, нейронная сеть, летальный аппарат.

Во многих приложениях, использующих модели исследуемых объектов в виде обыкновенных дифференциальных уравнений (ОДУ), возникает задача построения области достижимости соответствующей динамической системы. Например, такая задача возникает при решении проблемы траекторной безопасности летального аппарата [1], при решении близкой задачи о посадке вертолета на подвижный носитель [2]. Важной особенностью задачи построения области достижимости является то, что ее часто приходится решать в режиме реального масштаба времени.

Аналитически построить области достижимости удается лишь в простейших случаях, не представляющих практического интереса, поэтому при рассмотрении реальных ситуаций приходится использовать численные методы, т.е. строить некоторую дискретную аппроксимацию области достижимости или ее границ.

Численное построение области достижимости можно выполнить путем многократного интегрирования модельной системы ОДУ при различных управлениях (метод мультифиниша). В практических задачах типичной является ситуация, когда основные вычислительные затраты при построении области достижимости методом мультифи-ниша обусловлены затратами на интегрирование указанной системы ОДУ. Поэтому ускорить вычисления можно за счет параллельного интегрирования модельной системы ОДУ.

К параллельным методам интегрирования ОДУ относятся блочные методы, а также методы на основе геометрической схемы декомпозиции системы ОДУ.

Методы первого класса основаны на использовании одношаговых и многошаговых блочных формул интегрирования ОДУ, например формулы Адамса-Башфорта [3]. В вычислительной практике обычно используют блочные формулы не выше четырехточечных. Поэтому потенциальный параллелизм этой группы методов невелик.

Идея методов второго класса состоит в разбиении модельной системы ОДУ на подсистемы, число которых равно числу используемых процессоров вычислительной системы, и интегрировании каждой подсистемы на своем процессоре. В силу, как правило, относительно невысокой размерности вектора фазовых координат исследуемой динамической системы данные методы также принципиально не могут обеспечить значительное ускорение.

В настоящей работе полагается, что используемая многопроцессорная вычислительная система (МВС) представляет собой однородную М1МО-систему с общей или распределенной памятью [4]. Поскольку метод мультифиниша использует интегрирование совокупности систем модельных ОДУ, естественным в данном случае является распараллеливание путем разделения этой совокупности на наборы, число которых равно числу процессоров в используемой МВС. В работе рассматривается проблема балансировки загрузки МВС при распараллеливании в методе мультифиниша по указанной схеме, получены оценки эффективности параллельных вычислений.

Даже при вычислениях на однопроцессорной ЭВМ значительное ускорение при построении области достижимости можно обеспечить, аппроксимируя векторное поле модельной системы ОДУ [2]. При этом область допустимых значений вектор-функции, соответствующей правым частям системы ОДУ, покрывается некоторой сеткой. Предварительно во всех узлах этой сетки вычисляются значения указанной вектор-функции и сохраняются в памяти ЭВМ. При интегрировании модельной системы ОДУ используются эти значения или их подходящая аппроксимация. В работе получены оценки объема требуемой памяти, а также оценки ускорения данного метода.

Очевидной является идея использовать метод мультифиниша с заменой модельной системы ОДУ ее некоторой нейросетевой аппроксимацией. Однако данный метод не дает кардинального сокращения времени построения области достижимости. Поэтому в настоящей работе используется другой подход к построению области достижимости с использованием нейронных сетей — подход, основанный на построении нейронной сети, выходами которой являются непосредственно

координаты точек, принадлежащих области достижимости или ее границе.

В качестве примера рассмотрена задача построения области достижимости для летательного аппарата, описываемого системой ОДУ шестого порядка. Рассмотрено решение задачи с помощью метода муль-тифиниша, метода аппроксимации векторного поля этой системы ОДУ, а также нейросетевого метода.

Постановка задачи. Рассмотрим динамическую систему

х = ^(г,х,и), X(0) = X0, (1)

где X = X(£) = (ж!(£), х2(£),..., хп(£))т — п-мерный вектор фазовых переменных системы; и = и(£) = (м^), и2 (£),..., ит(£))т — т-мерный вектор управлений; ^(¿,Х, и) = (/!(^,Х, и),/2(^,Х, и), ..., /п(£, X, и))т — п-мерная вектор-функция; X0 = (х!, х2,..., хП)т — п-мерный вектор начальных условий, £ € [0,Т]. На вектор фазовых переменных X и вектор управления и наложены ограничения

X € Бх, и € Бц С Ьц[0,Т], (2)

где Ьц [0, Т] — некоторое пространство т-мерных функций, определенных на интервале [0,Т], например, пространство функций, интегрируемых с квадратом на этом интервале.

Среди фазовых переменных ж!, х2,..., хп выделим V ^ п переменных; не ограничивая общности, положим, что эти переменные образуют v-мерный вектор У(£) = (х!(£), х2(£),..., ху(£))т.

Областью достижимости Бу = Бу(Т, X0) системы (1) назовем множество всех возможных значений вектора У (Т), которые достигаются на решениях системы (1) при начальных условиях X0 и выполнении ограничений (2) на фазовые переменные и управления.

Ставим следующую прямую задачу: при заданных векторе начальных условий X0 и конечном времени Т построить область достижимости Бу системы (1). Также рассмотрим обратную задачу: при заданных векторе начальных условий X0, конечном времени Т и точке У(Т), принадлежащей области достижимости Бу, найти управление и € Бц, переводящее систему (1) в эту точку.

Для решения прямой задачи, очевидно, достаточно построить границу Гу области достижимости Бу. В некоторых случаях удается найти множество допустимых управлений БЦ С Бц, принадлежащих классу управлений ЬЦ [0, Т] С Ьц [0, Т], которые приводят систему (1) на эту границу [1]. В этом случае прямая задача сводится к построению границы Гу.

Сделаем следующие допущения. Во всех случаях при интегрировании системы ОДУ (1) используется алгоритм с постоянным шагом интегрирования Д£ = Т/К, где К — число шагов интегрирования.

Каждый шаг интегрирования требует l вычислений значений функции F(t,X, U), вычислительная сложность которой (число арифметических операций, необходимых для однократного вычисления значения этой функции) не зависит от аргументов t, X, U и равна CF. Вычислительная сложность затрат на реализацию алгоритма интегрирования на одном шаге интегрирования равна Cj = Cj (l).

Метод мультифиниша. Покроем множество Du некоторой сеткой с узлами Ui = U(t) = {uij(t), i G [1 : M], j G [1 : m]}, где M — общее число узлов сетки. Поставим в соответствие системе (1) совокупность M систем ОДУ с указанными управлениями:

X1 = F(t,Xi,Ui), Xi(0) = X0; ................o . (3)

XM = F(t XM, UM),XM(0) = X .

Тогда схему приближенного решения прямой задачи методом мультифиниша можно представить в следующем виде.

1. Путем интегрирования совокупности систем ОДУ (3) находим множество точек {Yi(T), i G [1 : M]}, представляющее собой дискретную аппроксимацию области Dy.

2. Сохраняем полученные наборы значений (Ui,Yi(T)) в памяти используемой ЭВМ.

3. Во множестве {Y (T)} находим граничные точки {Zij (T), ij G [1:ç]}, представляющие собой дискретную аппроксимацию границы Гу области достижимости; ç — число граничных точек, ç ^ M.

4. Сохраняем в памяти ЭВМ соответствующие наборы значений

(Uij, Zj (T)).

5. На основе множества граничных точек |Zij (T)} строим подходящую непрерывную аппроксимацию Г у границы Гу.

Задача построения дискретной или непрерывной аппроксимации границы Гу представляет собой классическую задачу вычислительной геометрии — задачу построения оболочки облака точек. Методы решения этой задачи хорошо разработаны [5]. Известно много программных систем, реализующих такие методы [6]. Поэтому в настоящей работе рассматривается только шаг 1 приведенной схемы приближенного решения прямой задачи.

Обратная задача в методе мультифиниша может быть решена по следующей схеме.

1. Во множестве точек {Yi(T)} находим точку Y*(T), которая в используемой метрике является ближайшей к заданной точке Y (T).

2. В качестве искомого управления принимаем управление из набора (Uij, Zij (T)), соответствующее точке Y*(T).

Схема распараллеливания вычислений при приближенном решении прямой задачи методом мультифиниша с использованием NN-процессорной МВС имеет вид:

1) процессор Pj, i G [1 : NN], получает от host-процессора векторы управлений U^-^+i, U(»-i)r+2,..., Ur;

2) интегрируя при каждом из этих управлений систему ОДУ (3), процессор находит координаты множества точек {Yj(T), j G [(i — 1)r + + 1 : ir]};

3) процессор Pi передает координаты полученных точек host-процессору;

4) host-процессор строит сначала дискретную, а затем непрерывную аппроксимацию границы rY области достижимости.

Здесь r = [M/NN] — число векторов управления, обрабатываемых процессором P^; [*] — символ ближайшего большего целого.

Отметим, что при большом числе узлов M задача построения дискретной и непрерывной аппроксимации границы rY области достижимости может иметь высокую вычислительную сложность и требовать распараллеливания.

При параллельном построении множества достижимости по рассмотренной схеме проблема балансировки загрузки МВС не возникает. При построении границы rY области достижимости DY ситуация может быть иной.

Равномерная декомпозиция точек переключения. Положим, что одну из границ множества достижимости формируют управления класса LU [0, T], в которых все компоненты вектора управления, кроме одной, постоянны, а компонента u1(t), например, имеет одну точку переключения:

( 1, t G [0,tS], u[(t, tS) = < ur (t) = COnst2, . . . , um(t) = COnstm.

( — 1, t G (tS,T],

Здесь tS G [0, T] — момент времени, в который происходит переключение управления u^(t).

Покроем интервал [0, T] равномерной сеткой с шагом AtS = T/M и узлами tS, i G [1 : M]. Положим, что шаг AtS кратен шагу At, так что AtS = qAt, где q = K/M ^ 1 — целое число. При этом естественно принять Ur = (ur(t,tS),ur,...,um) . Таким образом, схема распараллеливания в данном случае имеет вид, показанный на рис. 1.

В соответствии с этой схемой на процессоре P1 выполняется интегрирование системы ОДУ (3) при управлениях Uf,...,U,F, на процессоре P2 — при управлениях Ur+1,..., U2r и далее до процессора PNN, который выполняет интегрирование при управлениях

U(NN-1)r+1, . . . , UM.

Рис. 1. Схема распараллеливания при равномерной декомпозиции точек переключения

Вычисления на процессоре р организуем следующим образом.

Шаг 0 (этап разгона). Исходя из начальных условий Х0, выполняем интегрирование системы (3) при управлении (1, , ,..., и^) от момента времени 0 до момента времени 1 и сохраняем в памяти ЭВМ значения компонент векторов X) = Х((^-1)г+1,

Х (^г-1)г+2) = Х(0г-1)г+2, . . . , Х (^-1) = Хг>-1.

Шаг 1. Исходя из начальных условий! Х((^-1)г+1, выполняем интегрирование системы (3) при управлении (—1, , ,..., и^) от момента времени ¿^¿-1)г+1 до момента времени Т.

Шаг 2. Исходя из начальных условий Х((^-1)г+2, выполняем интегрирование системы (3) при управлении (-1,^2, от момента времени 1)г+2 до момента времени Т.

Шаг r. Исходя из начальных условий 1, выполняем интегрирование системы (3) при управлении (—1, , ,..., u^) от момента времени tfr_1 до момента времени T.

Замечание. В рассмотренной схеме каждый из последующих процессоров повторяет этап разгона своего предыдущего процессора и только затем выполняет интегрирование системы (3) на своем интервале [tg-1)r ,tf-1 ]. Если вычислительные затраты на интегрирование системы (3) велики, то может оказаться целесообразной более тонкая организация вычислений на этапе разгона — однократное интегрирование системы (3) на интервале [0, T] при управлении (1, u^, ,..., u^) на одном процессоре и передача необходимых начальных условий X°, X°,... остальным процессорам вычислительной системы. Однако, если r ^ 1, что можно считать типичной ситуацией, выигрыш от такой организации вычислений не может быть существенным.

В нашей работе [7] показано, что при равномерной декомпозиции точек переключения время последовательного решения задачи определяется выражением

M 1

Ti = tcalq (lCF + CI) ^ (M - (j - 1)) ~ ^tcalq (ICf + C7) M2,

j=i

а время параллельного решения задачи — выражением

TNN = Vtca;q(M — j)(1Cf + Cj) « tcaiq(1CF + Cj)

M2

Здесь и далее ¿са1 — время выполнения арифметической операции на одном процессоре используемой МВС.

Утверждение 1. Ускорение Б параллельного метода на основе равномерной декомпозиции точек переключения оценивается величиной

Б Ш

2

Из утверждения 1 следует, что равномерная декомпозиция точек переключения не может обеспечить ускорение, превышающее 50% максимального потенциально возможного ускорения, равного NN. Лучшей балансировки загрузки МВС можно добиться за счет неравномерной декомпозиции точек переключения.

Неравномерная декомпозиция точек переключения. Положим, что процессор р выполняет интегрирование системы ОДУ (3) при управлениях ..., а процессор Д+1 — при управлениях иь1+1,..., ЦТ, 1 ^ а < 6 < с ^ М. Потребуем, чтобы времена решения задачи процессорами р, Р+ были равны, т.е. чтобы выполнялось равенство Т = Т+1, г € [1 : (NN — 1)]. В работе [7] показано, что для этого величины а, 6, с должны удовлетворять рекуррентному уравнению

с2 — 2(М + 1)с + (а2 — 262 — 2Ма + 2(2М + 1)6) = 0. (4)

Если для процессора Р1 задать величины а = а1 = 1, 6 = 61, то из выражения (4) последовательно легко найти управления, которые должны обрабатывать процессоры Р2, Р3,..., .

Рассмотрим для примера ситуацию, когда М = 1000, 61 = 100. Из уравнения (4) вытекает, что в данном случае можно обеспечить приблизительно равномерную загрузку шести процессоров (таблица). Отметим, что для процессора число точек переключения вычислено по остаточному принципу и не удовлетворяет уравнению (4).

Как показано в работе [7], при неравномерной декомпозиции точек переключения время последовательного решения задачи можно оценить выражением

м

Т1 = tcaiq (ICf + Cj) M — 1 + ^ (M — i)

i=1

= ¿саг Я (1СР + С) (М — 1) (1 + 0,5М) а время параллельного решения — выражением

Т^м = ¿сагЯ (Юр + С) ((М + 1)61 — 0,561)

Таблица

Пример неравномерной декомпозиции точек переключения

Процессор Pi ts ts ts - ts tb

Pi 1 100 99

P2 101 212 111

P3 213 342 129

P4 343 504 151

P5 505 756 251

Pe 757 1000 243

Таким образом, имеем утверждение 2.

Утверждение 2. Ускорение Б параллельного метода на основе неравномерной декомпозиции точек переключения оценивается величиной

_ (М - 1)(1 + 0,5М) _ М2 = (М + 1)61 - 0,56? ~ 2МбГ-б|'

Из утверждения 2 следует, что в условиях рассмотренного примера ускорение можно записать как

10002

Б ~ -- ~ 5,26'

2 • 1000 • 100 - 1002

Заметное отличие ускорения от максимально возможного (равного 6) объясняется тем, что по причинам, указанным ранее, загрузка последнего процессора Р6 далека от оптимальной.

Явно решить рекуррентное уравнение (4) не удается. Поэтому число точек переключения 6?, обрабатываемых процессором Р?, приходится подбирать таким образом, чтобы вычислительная загрузка всех процессоров, включая последний, была примерно одинакова. При этом критерий окончания перебора удобно строить, исходя из того, что при несбалансированности загрузки процессора корни уравнения (4) оказываются комплексными.

Метод на основе аппроксимации векторного поля системы ОДУ. Идея данного метода состоит в следующем.

1. Предварительно множество [0, Т] х х покрываем некоторой сеткой П с узлами , мк,7), где ] € [1 : п], к € [1 : т],

а € [1 : N4, в € [1 : ], 7 € [1': ОД.

2. Предварительно во всех узлах сетки П вычисляем значения функции Р (£, X, и) и сохраняем в памяти ЭВМ.

3. Тем или иным методом численно интегрируем систему ОДУ (3) при управлении и», г € [1 : М], используя в качестве требуемых значе-

ний функции Р(£, Х, и) соответствующие интерполяции или аппроксимации сохраненных в памяти ЭВМ значений этой функции.

Здесь ИХ — число узлов сетки П по фазовой переменной х^; Ик — число узлов той же сетки по компоненте ик вектора управления.

Ускорение вычислений в данном методе достигается за счет того, что этапы 1 и 2 выполняются предварительно, так что при выполнении этапа 3 не требуется вычислять значения правых частей системы ОДУ (3).

Для того чтобы использовать весь потенциал рассматриваемого метода, значения функции Р(£, X, и) в узлах сетки П следует хранить в оперативной памяти компьютера. Примем, что длина вещественного

п т

числа равна 4 байтам. Пусть также NХ = П ИХ, Nи = П ИЦ. Тогда

3 = 1 к = 1 с учетом того факта, что общее число узлов сетки П равно Ии,

имеем утверждение 3.

Утверждение 3. Требуемый объем памяти ЭВМ для хранения значений функции Р(£, Х, и) в узлах сетки П равен

W = 4пИ4ИХ Ии байт.

Положим, что п = 5, т = 3; N4 = 10; ИХ = 10, € [1 : п]; Ик7 = 10, к € [1 : т]. Из утверждения 3 следует, что при этом требуемый объем памяти равен

W = 4 • 5 • 109 = 2 • 1010 байт и 20 Гбайт.

В простейшем случае в качестве значения Р(£, Х, и) используется значение этой функции в ближайшем узле сетки П. Такая схема вычислений применима при построении области в грубом приближении и означает использование кусочно-постоянной интерполяции функции

р (*,Х, и).

В работе [7] показана справедливость следующего утверждения.

Утверждение 4. Ускорение Б последовательного метода на основе кусочно-постоянной интерполяции функции Р(£, X, и) равно

+ С7 Ср ^ Б - —~-- » 1 (5)

С1 С1

Дальнейшего уменьшения времени решения задачи можно добиться путем организации параллельных вычислений. Схема распараллеливания при этом принципиально не отличается от схемы, рассмотренной ранее.

Кусочно-линейная интерполяция функции Р(£, X, и). При том же числе узлов сетки П большую точность вычислений можно обеспечить путем использования линейных интерполяций компонент функции Р(£, X, и). При этом рассмотренная ранее идея последовательного

метода существенно меняется. Дело в том, что построение интерполирующей функции, даже в том случае, когда она является линейной, требует значительных вычислительных затрат. Поэтому в данном случае целесообразно хранить не значения функции Р(£, X, и) в узлах сетки П, а заранее строить интерполирующие функции для каждой ячейки этой сетки и хранить в памяти ЭВМ коэффициенты полученных функций.

Рассмотрим методику линейной интерполяции функции Р(£, X, и) на примере г-й компоненты этой функции /¿(£, X, и) = / (£, X, и).

В ячейке П7, 7 € [1 : Nп], сетки П построим для функции /(£,Х, и) линейную интерполирующую функцию

/ (£, X, и) = аЛ0 + а_7д£ + «7,2X1 + ... + а7,п+1 Хп+

+ «7,п+2и1 + ... + Й7,п+т+1ит = «7,0 + (А, V), (6)

где А 7 = (07,1, а7,2,..., а7,п+т+1)т - (п + т + 1)-мерный вектор неизвестных коэффициентов, V = (£, х1,..., хп, и1,..., ит)т — вектор такой же размерности аргументов функции Р(£, X, и); (*, *) — символ скалярного произведения векторов;

пт

N° = (И - 1) П(ИХ - 1)П(Ики - 1) —

¿=1 к=1

общее число ячеек в сетке П.

Коэффициенты а7,0, а71, а7,2,..., а7,п+т+1 функции (6) находят из системы (п + т + 2) линейных алгебраических уравнений (СЛАУ):

07,0 + (А ,К7,1) = / (К7,1);

Д7,.0 + (¿7 ,К7,2) = / 7 ), (7)

а7,0 + (А7, ^7,п+т+2) = /(^7,п+т+2),

где (К^, г € [1 : (п + т +2)]} — совокупность значений аргументов функции Р(£,Х, и) в рассматриваемой части узлов ячейки П7. Отметим, что общее число узлов этой ячейки равно, очевидно, 2п+т+1.

Утверждение 5. Объем памяти ЭВМ, необходимый для хранения коэффициентов а70, а71, а72,..., а7,п+т+1 линейной интерполяции компонент функции Р(£, X, и) на сетке П, равен

W = 4п(п + т + 2)И° байт.

Из данного утверждения следует, что в условиях рассмотренного примера

W = 4 • 5 • 10 х 99 « 80 • 109 байт и 80 Гбайт.

Ускорение последовательного метода на основе кусочно-линейной интерполяции функции Р(¿, X, и) определяется так же, как в методе на основе кусочно-постоянной интерполяции.

Кусочно-линейная аппроксимация функции Р(¿, X, и). Известной проблемой интерполяции многомерных функций является невозможность использования произвольных сеток. Поэтому в вычислительной практике обычно используют не интерполяцию функции Р (¿, X, и), а ее аппроксимацию на основе метода наименьших квадратов (МНК).

Аналогично предыдущему пункту методику линейной аппроксимации функции Р(¿, X, и) рассмотрим на примере г-й компоненты этой функции /¿(¿, X, и) = /(¿, X, и). Для кусочно-линейной аппроксимирующей функции в этом случае имеем следующую СЛАУ относительно неизвестных коэффициентов а^о, а^, а^2,..., а7,п+т+1:

а/,о + (А,- ,К/д) = / (V},?);

а/,.о + (А} ,К/,2) = / (К/,2);

а},о + (А} ) = / (/*),

где й — число используемых узлов ячейки П/; (п+т+2) < й ^ 2п+т+1.

В соответствии с методикой МНК решение этой системы получаем в виде [7]

А/ = (УV/)-1 V} Г/, (8)

где А/ = (а/,о,а/,1,''' ,а/,п+т+1 )т — (п + т + 2)-мерный вектор; Г/ = (/(К/д),/(К/,2),''',/(К/,й))т — й-мерный вектор; V/ =

= (К/,1, К/,2,......, К/,^ — й х (п+т+2)- матрица, которая в качестве

столбцов имеет (п + т + 2)-мерные векторы К/,» = (1, , , - - -, /, ,''',мт/)т. Важно, что, кроме вырожденного случая, когда все точки /(К/д),/(К/,2),''', /(У/,й) лежат в одной гиперплоскости, решение СЛАУ (8) существует и единственно.

Объем памяти ЭВМ, необходимый для хранения коэффициентов а^о, а/,1, а^2,''', а7,п+т+1 кусочно-линейной аппроксимации функции Ри) на сетке П, определяется утверждением 5. Ускорение последовательного метода на основе кусочно-линейной аппроксимации функции Р(¿, X, и) определяется утверждением 4.

Большое число вариантов рассмотренного подхода к кусочно-линейной аппроксимации функции Р(¿, X, и) можно получить, привлекая результаты теории планирования эксперимента. Так, каждый узел сетки П можно использовать в качестве центральной точки для построения регрессионного плана эксперимента первого порядка, удовлетворяющего критерию Б-оптимальности, А-оптимальности, Е-оптимальности и т.д. [8]. СЛАУ (8) при этом сохраняет свою структуру.

а б

Рис. 2. Нейросетевая аппроксимация области достижимости

Замечание. Повысить точность аппроксимации функции Р(£, X, и) можно, используя квадратичную аппроксимацию. При этом естественно ориентироваться на планы эксперимента второго порядка, такие как центрально-композиционные планы, композиционные планы Хартли и Вестлейка [9]. В алгоритмическом отношении для аппроксимации удобно использовать композицию одномерных полиномов (метод Брандона) [10].

Нейросетевые методы. Для аппроксимации многомерных функций чаще всего используют многослойные нейронные сети прямого распространения (многослойные персептроны) [11]. В качестве аппроксиматора функций могут использоваться также радиально-базисные нейронные сети. При прочих равных условиях аппроксимация этими сетями требует значительно большего числа нейронов, чем аппроксимация многослойными персептронами. Поэтому радиально-базисные нейронные сети в работе не рассматриваются.

Для решения прямой и обратной задач построения области достижимости используем нейронные сети НС1, НС2 соответственно (рис.2). Входами нейронной сети НС1 являются вектор начальных условий X0, конечный момент времени Т и совокупность управлений (и^,г € [1 : М]}, а выходами — вектор У(Т). Нейронная сеть НС2 имеет в качестве входов величины X0, Т, У (Т), а в качестве выходов — вектор и(£) € ^.

В случае, когда известно множество допустимых управлений С , приводящих систему (3) на границу области достижимости Гу, могут быть использованы сети, аналогичные сетям НС 1, НС2. Отличие состоит в том, что здесь управления (и} должны принадлежать множеству , а вектор У (Т) — границе Гу множества достижимости. В случае сложной топологии границы Гу может оказаться необходимым использовать различные нейронные сети для аппроксимации разных фрагментов этой границы.

Пример. Рассмотрим летательный аппарат, уравнения движения центра масс которого в нормальной земной системе координат Оху^ описываются системой нелинейных дифференциальных уравнений:

v = g (uT — sin O);

O = g/v (u cos Yc — cos O);

Ф = —gu sin yc /v cosO;

x = v cos O cos Ф;

y = v sin O;

z = —v cos O sin Ф,

где v — скорость летательного аппарата; O — угол наклона траектории; Ф — угол поворота траектории; y — высота летательного аппарата; uT — тангенциальная перегрузка; u — нормальная перегрузка; yC — скоростной угол крена; g — ускорение свободного падения [1].

Управлениями летательного аппарата являются тангенциальная перегрузка, нормальная перегрузка и скоростной угол крена, так что U = (uT, u, yC)т. На управления наложены ограничения

и™ ^ UT ^

max.

U

|u| ^ u

= -1,6;

U =

U

= 0,6;

|Yc | ^ n.

Дальняя, ближняя и боковая границы области достижимости системы (9) формируются управлениями, принадлежащими классу кусочно-постоянных управлений с одной точкой переключения [1]. Ограничимся рассмотрением дальней границы области достижимости. Для каждого допустимого угла крена yc (t) = const структура управлений, формирующих эту границу, приведена на рис. 3.

Метод мультифиниша. Для экспериментального исследования эффективности метода мультифиниша разработана MPI-программа [4], реализующая этот метод. Эксперименты выполнены на виртуальном кластере, созданном с помощью программной системы VMware [12] и функционирующим под управлением свободно распространяемой операционной системы Ubuntu [13]. Результаты экспериментов приведены на рис. 4, из которого следует, что при числе процессоров от 2 до 10 достигается ускорение, близкое к расчетному.

Рис. 3. Структура управлений, формирующих дальнюю границу области достижимости

Рис. 4. Ускорение метода мультифиниша

Приведенные на рис.4 экспериментальные результаты получены при следующих значениях параметров задачи: V(0) = 330 м/с; Э(0) = 0; Ф(0) = 0; х(0) = 0; у(0) = 10000м; г(0) = 0; Т =15 с.

Аппроксимация векторного поля системы ОДУ. Рассмотрим плоскую задачу (9), когда Ф(£) = 0, £ € [0,Т], Т = 30 с. Предположим, что область фазового пространства системы (9) определяется неравенствами

1Мв ^ V(£) ^ 1,9vзв; -90° ^ в(£) ^ 90°; Ф(£) = 0;

0 ^ х(£) ^ 5500 м; 0 ^ у(£) ^ 2700 м; г(£) = 0.

Здесь vзв = 331,46 м/с — скорость звука в воздухе. Принимаются следующие начальные условия летательного аппарата: V(0) = 132,584м/с; в(0) = 0; Ф(0) = 0; х(0) = 0; у(0) = 1000м; г(0) = 0; 7С(0) = 0.

В соответствии с методикой покроем интервал времени [0,Т] равномерной сеткой с 20 узлами. Указанную область фазового пространства системы (9) покроем равномерной сеткой так, что число шагов этой сетки равно 18, 20, 20, 18 по измерениям V, в, х, у соответственно. Наконец, в качестве моментов времени переключения управления и рассмотрим М = 30 узлов равномерной сетки, покрывающей интервал времени [0, Т] (рис. 3, а).

Численные эксперименты выполнены для кусочно-постоянной аппроксимации компонент функции Р(£,Х, и) на указанной сетке (рис. 5).

Эксперименты показали, что метод обеспечивает погрешность в определении расстояния от начала координат до дальней границы области достижимости, не превышающую 4%. В качестве точного положения этой границы принималось положение, полученное путем численного интегрирования системы (9) при соответствующем управлении.

Нейросетевая аппроксимация. В данном случае начальный угол поворота траектории Ф(0) и начальные координаты х(0), у (0), г (0) ле-

Рис. 5. Дальняя граница области достижимости ("плоская" задача):

* — точная граница; + — приближенная граница

тательного аппарата приняты равными нулю. Так что исходное положение летательного аппарата определяется начальной скоростью V (0)

и начальным углом наклона траектории ©(0): 1,6^в ^ V (0) ^ 1,9^в; п п

--^ ©(0) ^ —. Варьируемыми параметрами являются точки пере-

8 8

ключения нормальной перегрузки и и значения скоростного угла крена 7с. Полагается, что е [0,Т], где Т е [5,15] с. Поскольку область достижимости симметрична относительности плоскости Оху, принято, что 7с е [0,п].

Для указанной задачи были построены нейронные сети НС 1, НС2. Во всех случаях обучающая выборка состояла из 4500 элементов, из которых 20% случайно выбранных точек использовались в качестве проверочных для определения момента окончания обучения. Начальные веса и пороговые значения нейронов устанавливались в соответствии с алгоритмом инициализации Нгуен-Видроу (Nguyen-Widrow) [11]. Обучение каждой сети выполнялось в среднем 10 раз с различными начальными весами и пороговыми значениями, на основе чего выбиралась сеть с минимальным значением ошибки. Для исследования точности обученной сети использовалась тестовая выборка, содержащая более 20000 элементов. Все эксперименты выполнены в программной среде Ма1;ЬаЪ 6 [11].

Сеть НС1. Входной слой сети НС1 состоит из пяти элементов, соответствующих величинам V(0), ©(0), Т, 7С, . Сеть содержит два скрытых слоя, содержащих одинаковое число нейронов с функцией активации в виде гиперболического тангенса. Выходной слой сети состоит из трех нейронов с линейной функцией активации, соответствующих координатам х(Т), у (Т), г (Т) конечного положения летательного аппарата.

Исследование показало, что при числе нейронов в сети, равном 22, максимальная относительная ошибка определения нейронной сетью

Рис. 6. Нейросетевая аппроксимация дальней границы области достижимости

НС1 координат точек границы области достижимости не превышает 12 %. С ростом числа нейронов эта ошибка быстро уменьшается и при 40 нейронах не превышает 4 %. Результаты исследования приведены на рис. 6, на котором показана аппроксимация дальней границы области достижимости нейронной сетью НС1 с 32 нейронами для случая Т =15 с, V(0) = 1^в, в(0) = п/10.

Сеть НС2 аналогична сети НС1. Входной слой сети НС2 состоит из пяти элементов, соответствующих величинам V(0), в(0), Т, а, в, где а, в — углы, задающие требуемое направление в пространстве. Выходной слой состоит из двух нейронов, соответствующих положению точки переключения и значению скоростного угла крена 7С.

Для оценки точности аппроксимации сети определяются координаты точек границы области достижимости, соответствующие управлению, получаемому на выходе нейронной сети, и сравниваются с эталонными координатами этих точек, вычисленными путем интегрирования системы ОДУ (9).

Исследование показало, что максимальная относительная ошибка определения нейронной сетью НС2 координат точек границы области достижимости при 22 нейронах в сети не превышает 9 %, а при 40 нейронах составляет ^3 %.

Детально результаты исследования метода нейросетевой аппроксимации границы области достижимости летательного аппарата приведены в работах [14, 15].

Заключение. В работе рассмотрены три метода построения области достижимости динамической системы — метод мультифиниша, метод на основе аппроксимации векторного поля соответствующей системы ОДУ, а также метод на основе нейронных сетей. Для первых

двух методов предложены параллельные алгоритмы, ориентированные на вычислительные системы класса MIMD. Эффективность всех рассмотренных методов исследована на примере построения границы области достижимости летательного аппарата.

Проведенное исследование показало перспективность практического использования всех рассмотренных методов. Тот или иной метод может быть рекомендован в зависимости от вычислительной сложности правых частей системы ОДУ, описывающей исследуемую динамическую систему, а также от располагаемых вычислительных ресурсов.

Кроме рассмотренных методов интерес представляет использование для приближенного построения области достижимости нечетких множеств. Отметим также целесообразность исследования метода, основная идея которого состоит в следующем: для представительного множества начальных условий динамической системы тем или иным образом строим совокупность областей достижимости; границу каждой области аппроксимируем подходящей функцией; функциональную зависимость коэффициентов этой функции от начальных условий аппроксимируем, например, с помощью нейронной сети. Можно ожидать, что такой метод позволит получить приемлемую точность аппроксимации множества достижимости для широкого набора начальных условий.

Авторы благодарят А.Г. Трофимова за идею использования нейронных сетей для аппроксимации границ области достижимости, а также К.О. Вишневецкого и Д.С. Кулеша за проведение вычислений с использованием методов мультифиниша и аппроксимации векторного поля.

СПИСОК ЛИТЕРАТУРЫ

1. Воронов Е. М., Карпунин А. А. Алгоритм оценки границ области достижимости летательного аппарата с учетом тяги // Вестник МГТУ. Сер. "Приборостроение". - 2007. - № 4(69). - С. 81-99.

2. Г у р м а н В. И., К в о к о в В. И., У х и н М. Ю. Приближенные методы оптимизации управления летальным аппаратом // Автоматика и телемеханика. -2008. -№ 4. - С. 191-201.

3. ОртегаДж., Пул У Введение в численные методы решения дифференциальных уравнений / Пер. с англ., под ред. А.А.Абрамова. - М.: Наука, 1986. -288 с.

4. В о е в о д и н В. В. Параллельные вычисления. - СПб.: БХВ-Петербург, 2004. - 608 с.

5. С к в о р ц о в А. В. Обзор алгоритмов построения триангуляции Делоне // Вычислительные методы и программирование. - 2002. - Т. 3. - С. 14-39.

6. Directory of Computational Geometry Software [Электронный ресурс]. /http://www.geom.uiuc.edu/software/cglist/

7. Воронов Е. М. Численные методы построения области достижимости динамической системы [Электронный ресурс] // Электронное научно-техническое издание: наука и образование. - 2010. - № 1. http://www.technomag.edu.ru

8.Карпенко А. П. Аппроксимация функции предпочтений лица, принимающего решения, в задаче многокритериальной оптимизации. 1. Методы на основе планов первого порядка [Электронный ресурс] // Электронное научно-техническое издание: наука и образование. - 2008. - № 3. http://www.technomag.edu.ru

9.Карпенко А. П. Аппроксимация функции предпочтений лица, принимающего решения, в задаче многокритериальной оптимизации. 2. Методы на основе планов второго порядка [Электронный ресурс] // Электронное научно-техническое издание: наука и образование. - 2008. - № 3. http://www.technomag.edu.ru

10. Х о х л о в С. Ф., Ш к о л а О. И. Применение метода Брандона для обработки экспериментальных данных // Вопросы химии и химической технологии. - 1973. -Вып.28.-C.204-207.

11. Медведев В. С., Потемкин В. Г. Нейронные сети. Matlab 6. - М.: Диалог-МИФИ, 2001. - 496 с.

12. VM ware [Электронный ресурс]. http://www.vmware.com/

13. Ubuntu [Электронный ресурс]. http://www.ubuntu.com/

14. К о з л о в а О. Г. Нейросетевая аппроксимация границы области достижимости летательного аппарата [Электронный ресурс] // Электронное научно-техническое издание: наука и образование. - 2009. - № 7. http://www.technomag.edu.ru/doc/129990.html

15.Козлова О. Г. Нейросетевая аппроксимация границы области достижимости летательного аппарата. Трехмерный случай [Электронный ресурс] // Электронное научно-техническое издание: наука и образование. // - 2009. - № 8. http://www.technomag.edu.ru/doc/130282.html

Статья поступила в редакцию 28.01.2010

Евгений Михайлович Воронов родился в 1940 г., окончил в 1963 г. МВТУ им. Н.Э. Баумана и в 1969 МГУ им. М.В. Ломоносова. Д-р техн. наук, профессор кафедры "Системы автоматического управления" МГТУ им. Н.Э. Баумана. Автор 150 научных работ в области теории управления и ее приложений.

Ye.M. Voronov (b. 1940) graduated from the Bauman Moscow Higher Technical School in 1963 and the Lomonosov Moscow State University in 1969. D. Sc. (Eng.), professor of "Automatic Control Systems" department of the Bauman Moscow State Technical University. Author of 150 publications in the field of control theory and its applications.

Анатолий Павлович Карпенко родился в 1947 г., окончил МВТУ им. Н.Э. Баумана в 1975 г. Автор 141 научной работы в области информационных технологий.

A.P. Karpenko (b. 1947) graduated from the Bauman Moscow Higher Technical School in 1975. Author of 141 publications in the field of information technologies.

Оксана Геннадьевна Козлова родилась в 1984 г., окончила МВТУ им. Н.Э. Баумана в 2009 г. Автор пяти научных работ в области информационных технологий. O.G. Kozlova (b. 1984) graduated from the Bauman Moscow Higher Technical School in 2009. Author of 5 publications in the field of information technologies.

Владимир Андреевич Федин родился в 1986 г., окончил МВТУ им. Н.Э. Баумана в 2008 г. Автор пяти научных работ в области информационных технологий. V.A. Fedin (b. 1986) graduated from the Bauman Moscow Higher Technical School in 2008. Author of 5 publications in the field of information technologies.