CC BY
18
УДК 517.977
ПАРАЛЛЕЛЬНЫЙ ПРОГРАММНЫЙ КОМПЛЕКС РЕШЕНИЯ НЕГОЛОНОМНЫХ ЗАДАЧ УПРАВЛЕНИЯ
(Работа выполнена в рамках ФЦП «Исследования и разработки по приоритетным направлениям развития научно-технологического комплекса России на 2007-2013 годы», ГК№ 07.514.11.4033)
А.П. Маштаков
(Институт программных систем им. А.К. Айламазяна РАН, г. Переславль-Залесский,
а1ехеу. mashtako 1@дтаИ. сот)
Описан метод управления нелинейными системами с линейным управлением на основе нильпотентной аппроксимации. Представлен алгоритм приближенного решения конструктивной задачи управления пятимерными системами такого вида в классах кусочно-постоянных и оптимальных управлений для аппроксимирующей системы. Приведенный алгоритм был положен в основу параллельного программного комплекса МойопР1апт^235, предназначенного для решения поставленной задачи.
Ключевые слова нильпотентная аппроксимация, оптимальное управление, параллельные алгоритмы и программы.
Перевод механической системы в требуемую конфигурацию является одной из важнейших задач робототехники и инженерии: это вывод спутника на орбиту, управление мобильным роботом, вращение твердого тела в руке робота-манипулятора и др. В математике подобная задача формализуется как конструктивная задача управления.
Имеется динамика механической системы, описываемая системой дифференциальных уравнений с функциональными параметрами (управления). Требуется найти управляющие воздействия из заданного класса функций, при применении которых система переходит из начального в заданное целевое состояние.
В работе рассматривается двухточечная задача управления следующего вида:
д = и (() X (¿) + и (/) Х2 (¿), (1)
ч(0) = ч0, д(т) = ^, (2)
где пространство состояний qeQ - это связное пятимерное гладкое многообразие, dim(Q)=5; управление принимает значения на двухмерной плоскости, (иь и2)е^2; гладкие векторные поля Х^), Х2^) удовлетворяют условию полного ранга [1] на многообразии Q. Требуется найти управление и(0=(и1(0, и2(0), переводящее систему (1) за время Т>0 из начального состояния q0 в терминальное ql с любой ранее заданной точностью е>0, то есть в такое состояние ¿¡1, что (^, ¿¡1) < е , где dist -расстояние на многообразии Q (например, если
Q=ЭТ5, то dist(q1, $ ) = ^£ ($ - $ )2 ).
Системы вида (1) характеризуются тем, что размерность пространства управлений меньше размерности пространства состояний йт(^2)= =2<dim(Q)=5, однако любые точки пространства состояний могут соединяться траекторией системы. В теории управления такие системы называются вполне неголономными. Нелинейная система (1), линейно зависящая от управлений, количество которых меньше размерности пространства со-
стояний, характеризуется тем, что возможность передвижения по различным направлениям неодинакова. Величина смещения в направлении полей Х^), Х2^) за малое время t есть 0(0, в направлении коммутатора Х3^)=[ХЬ Х2]^) есть 0(/2), в направлении коммутаторов более высокого порядка Х4^)=[ХЬ Хз]^) и Х5^)=|Х2, Хз]^) есть 0(^). В силу этой анизотропии пространства состояний задача управления для таких систем весьма нетривиальна.
В последнее время конструктивная задача управления активно изучается в нелинейной теории управления. Удовлетворительное решение она имеет лишь для некоторых специальных классов систем [2]. Однако пятимерные системы общего положения с двумя управлениями имеют вектор роста (2, 3, 5), потому эти результаты неприменимы к таким системам.
В данной работе представлен способ отыскания приближенного решения задачи (1)-(2), основанный на построении нильпотентной аппроксимации. Идея метода в том, что исходная нелинейная система заменяется приближенной нильпотентной системой, для которой точно решается задача управления. Затем найденные управления подставляются в исходную систему. Если состояние, достигнутое после применения найденного управления, отличается от желаемого в пределах допустимой погрешности, задача считается решенной, иначе процедура повторяется. Управления, точно решающие нильпотентную систему, дают приближенное решение исходной задачи управления в малой окрестности целевой точки. Метод нильпотентной аппроксимации применим к задачам управления общего вида; существенно только умение решать задачу управления для нильпотентной аппроксимации. Этот метод предложен в [3].
В данной работе рассматривается задача управления нелинейными пятимерными системами общего вида, линейно зависящими от двухмерного управления. Конкретными примерами систем
такого вида, которые в робототехнике имеют самостоятельное значение из-за своего прикладного значения, являются система, моделирующая качение шара по плоскости без прокручивания и проскальзывания, и система, моделирующая движение машины с двумя прицепами по плоскости.
Алгоритм поиска приближенного управления на основе вычисления траекторий аппроксимирующей системы и его реализация программным комплексом (ПК) MotionPlanning235. Итерационный алгоритм приближенного решения задачи (1)-(2) основан на локальном приближении исходной нелинейной системы нильпотентной канонической системой (3), для которой задача управления решается точно на каждой итерации.
Итак, имеется исходная система в начальном состоянии q0. Требуется найти управление, переводящее систему в конечное состояние qi с предварительно заданной точностью е. Поиск приближенного решения осуществляется следующим образом.
1. В окрестности конечной точки qi строится аппроксимирующая система, для которой точно решается задача управления.
2. Найденные управления подставляются в исходную систему. Если достигнутое состояние не попадает в е-окрестность конечного состояния, значит, нужная точность не достигнута. Алгоритм повторяется снова для исходной системы, где в качестве начального состояния выбирается состояние, достигнутое на предыдущей итерации алгоритма, иначе вычисления прекращаются.
В среде Mathematica 8 (http://www.wolfram.com/ mathematica/) был реализован параллельный ПК MotionPlanning235 (см. рис. 1) решения задачи (1)-(2). Он является дополнительным пакетом для системы Wolfram Mathematica (MotionPlann-ing235.m) и состоит из следующих основных модулей:
- NilpotentApproximation235 строит нильпо-тентную аппроксимацию NA (см. (3)) системы (1) и возвращает замену переменных т=А°Ф, в которых NA имеет канонический вид;
- CPControl235 решает задачу (1)-(2) в классе кусочно-постоянных управлений;
- OptimalControl235 решает задачу (1)-(2) в классе оптимальных управлений.
Модуль NilpotentApproximation235 построения нильпотентной аппроксимации. Локальное приближение управляемой системы другой, более простой, широко используется в теории управления. Обычно в качестве локальной аппроксимации используется линеаризация управляемой системы. Однако для линейных по управлению систем вида (1) линеаризация дает слишком грубое приближение. Так как размерность управления меньше размерности состояния, линеаризация не может быть вполне управляемой. Естественной заменой линейной аппроксимации в этом случае является
нильпотентная аппроксимация - наиболее простая система, сохраняющая структуру управляемости исходной системы (в частности, сохраняется такой важный инвариант, как вектор роста).
В данной работе использован алгоритм вычисления нильпотентной аппроксимации для линейных по управлению систем, предложенный в [2]. Этот алгоритм был конкретизирован для систем вида (1), а именно, вычислена замена переменных для перехода в привилегированные координаты A:q^■z, а векторные поля нильпотентной аппроксимации Х1 (г) , Х2 (г) в привилегированных координатах zi, i=1, ..., 5, были выражены через векторные поля исходной системы Х^), X2(z) и их производные.
Кроме того, алгоритм Беллаиша дополнен следующим образом: выполняется переход в систему координат у (замена переменных Ф^^у), в которой нильпотентная аппроксимация имеет канонический вид
у = и1>
(3)
У 2 = и2 >
У Ъ = -2 (У 1й2 - У2и1 ) >
у4 = 1 (У + У2 ) и2 >
У5 =-1 (Л2 + У2 ) и1 >
а граничные условия представлены как
.у(0)=у0, У(Г)=(0, 0, 0, 0, 0). (4)
Модуль CPControl235 решения задачи в классе кусочно-постоянных управлений. В ПК Мой-опР1апп^235 реализован модуль СРСоп^о1235 решения задачи управления нильпотентной канонической системой (3) с граничными условиями (4) в классе кусочно-постоянных функций. При этом время можно перепараметризовать так, что 7=1. Требуется найти управления (щ((), и2(/)) в классе кусочно-постоянных функций на отрезке /е[0, 1], переводящие систему из произвольного начального состояния у0=Ш5 в начало координат. Доказано следующее утверждение.
Предложение. Для решения задачи управления (3)-(4) достаточно управлений с тремя точками переключения:
а, 2, t е [0, 0,25), t е [0,25, 0,5), t е [0,5, 0,75], t е (0,75, 1].
в 1,2 , Yl,2 > §, ,,
Коэффициенты а, в, у, и 5 определяются из системы пяти алгебраических уравнений с восемью неизвестными. Получается трехпараметриче-ское семейство решений, причем для любого начального состояния у0=Ш5 существует способ зафиксировать свободные параметры так, чтобы
Начало
xs (символьное обозначение переменных), Х1^), Х2^) (векторные поля при управлениях), q0, q1, (начальное и целевое состояния), eps (точность), Т(конечный момент времени), и (класс управлений)
NilpotentApproximation235 (Построение нильпотентной аппроксимации NA)
CPControl235 U=CP Выбран класс -кусочно-постоянных- управлений Switch U /-""Выбор класса функций/Х. в котором ищется U=OC ^_ Выбран класс оптимальных управлений OptimalControl235
управление
\ У У
Вывод результатов: ui(t), u2(t) (найденные управления), q(t) (соответствующая траектория) N '
q=qo
CPSolver
(Поиск кусочно-постоянных управлений ul(t), u2(t), переводящих NA из q в ql за время tl=1)
I ~~
Trajectory
(Вычисление траектории q(t) исходной системы, соответствующей найденным _ui(t), U2(t))_
q=q(ti)
Нужная точность достигнута
_к_
Конкатенация управлений и перепараметризация времени
Конец
л
OCSolver
(Поиск оптимальных управлений ul(t), u2(t), переводящих NAиз яв ql за время tl=topt)
Trajectory
(Вычисление траектории q(t) исходной системы, соответствующей найденным ui(t), U2(t))
Нужная точность достигнута
Конкатенация управлений и перепараметризация времени
q=qo
Рис. 1. Структура ПКMotionPlanning235
получалось решение без особенностей. В модуле СРСоп^о1235 фиксация свободных параметров осуществляется так, чтобы в соответствующей траектории не было больших отклонений, а именно, используется критерий тах | и, | ^-тт.
Модуль OptimalControl235 решения задачи в классе оптимальныгх управлений. Для управляемой системы (3) с граничными условиями (4) рассматривается задача оптимального управления с естественным интегральным критерием (функционалом субримановой длины)
U
U + u2 dt ^ min.
(5)
Эта задача известна в теории управления как обобщенная задача Дидоны, а в субримановой геометрии - как субриманова задача с вектором роста (2, 3, 5). Она имеет богатую историю и была детально теоретически изучена в [4]. Основным результатом теоретических исследований явилось описание структуры экспоненциального отображения и первых времен Максвелла вдоль экстремальных траекторий. На основе этого задача оптимального управления (3)-(5) была сведена к задаче решения пяти алгебраических уравнений в неэлементарных функциях с пятью неизвестными. Явно решить эту систему практически невозмож-
но из-за сложности получившихся формул, поэтому предложено использовать численные методы.
Приведем некоторые результаты из статьи [4], необходимые для дальнейшего изложения. Семейство экстремальных траекторий в задаче оптимального управления параметризуется экспоненциальным отображением ехр, переводящим пару (вектор сопряженных переменных, время) в соответствующую точку экстремальной траектории. Прообраз Ь и образ М экспоненциального отображения ехр: Ь^М разбиваются поверхностями раз-
4
реза на четыре области (Ь = и Ь1 в прообразе и
¿=1
4
М = иМв образе), таких, что Ь, переходит в М,
¿=1
и экспоненциальное отображение является диффеоморфизмом на этих областях. В свою очередь, каждая область Ь, разбивается на непересекающиеся множества С1 и С2, в которых экспоненциальное отображение задается разными формулами. Задача построения оптимального синтеза состоит в обращении экспоненциального отображения. При этом разбиение Ь на Ь, и М на М, известно, а разбиение на С, неизвестно. В задаче (3)-(5) была найдена двухмерная группа симметрий, состоящая из вращений и растяжений. Факторизация задачи по этой группе уменьшает размерность задачи с пяти до трех. Получается система из трех уравнений с тремя неизвестными:
' P(u, v, k ) = P1,
■Q(u, v, k ) = Q, (6)
R(u, v, k ) = R1.
Итак, для построения оптимального синтеза в задаче (3)-(5) требуется решить систему (6) относительно u, v, k при заданной правой части Pb Qb R\. Заметим, что функции P(u, v, k), Q(u, v, k) и R(u, v, k) выражены в эллиптических функциях Якоби и эллиптических интегралах первого и второго рода и имеют сложную аналитическую запись. Известно, что при любой правой части (за исключением особых множеств меньшей размерности) система (6) имеет единственный корень.
Для решения поставленной задачи был разработан модуль OptimalControl235. Схема его работы представлена на рисунке 2. Пользователь легко может организовать параллельное вычисление корня на нескольких ядрах процессора. Для этого требуется указать системе Mathematica 8 (функция недоступна в более ранних версиях) количество ядер и изменить при необходимости функцию SolveHen (выбрать, какие алгоритмы решения системы - AlgorithmA, AlgorithmB, Algorithme, AlgorithmD - будут выполняться одновременно и сколько именно). Кроме стандартного метода Ньютона (AlgorithmA) и метода хорд (AlgorithmB), пользователю предлагаются гибридные методы (Algorithme, AlgorithmD). Проведенное обширное тестирование показало, что в большинстве случа-
Рис. 2. Схема работы модуля Optimalcontrol235
ев первым результат выдает А^опШтА или А^о-пШтВ, хотя и бывают ситуации, когда гибридные методы справляются с задачей быстрее. Поэтому по умолчанию предлагается делить доступные ядра поровну между А^ойШтА и А^опШтВ. Отметим, что одновременный запуск одного и того же алгоритма на разных ядрах имеет смысл, так как во всех алгоритмах используется генератор случайных начальных приближений. Так как основное время вычисления занимает выбор удачного начального приближения, использование нескольких ядер дает существенное ускорение.
Дополнительные функции ПК Motionlan-ning235. Помимо основных модулей Nilpotent-АрргохтаШп235, СРСоп^о1235 и ОрйтаЮоп-^о1235, решающих задачу (1)-(2), пользователю предоставляются некоторые дополнительные инструменты для представления результатов и слежения за процессом вычислений. Речь идет о следующих функциях:
1. TrajectoryNatPar[X1, Х2, {и1,и2,Т}, q0, xs] возвращает траекторию д(0 системы (1), где ¿е[0, Т], соответствующую управлениям {и1, и2}, выходящую из точки q0.
2. PlotTrajectoryNatPar[trajectory, q0, q1, Т] строит разным цветом графики пяти компонент заданной траектории й^ейогу^) за время ¿е[0, Т], а также отмечает на нем заданные состояния q0 и q1. Эта функция может использоваться для визуальной оценки достигнутого состояния и найденных траекторий. Пользователь может включить режим работы ПК, в котором графики найденных траекторий будут выводиться на экран на каждой итерации.
3. PlotTrajectoryOtklNatPar[trajectory, q0, q1,
Т] строит разным цветом графики отклонений пяти компонент заданной траектории йа|ес1огу(0 за время ¿е[0, Т] от состояния q1.
4. AnimateCar[trajectory, Т, q0, q1, N {{xmin, хтах}, {ymin, утах}}] строит анимацию движения машины с двумя прицепами при движении по траектории й^ейогу^) за время ¿е[0, Т] из состояния q0. Область, по которой движется машина с прицепами, ограничивается прямоугольником {{хтш, хтах},{утт, >тах}}. В результате на жестком диске создается последовательность из N кадров. Последовательность кадров - это упорядоченный набор изображений .род, который впоследствии легко может быть собран в видеофайл стандартными утилитами. Размер изображения является параметром, доступным пользователю (по умолчанию 1024x768 пикселей). Помимо положения машины и прицепов, на каждый кадр пунктиром наносится состояние q1. Функция предназначена как для самостоятельного использования (для моделирования системы «машина с двумя прицепами» и ее исследования), так и для проверки решений, найденных с помощью ПК MotionPlanning235.
5. AnimateBall[trajectory, Т, q0, q1, N визуализирует качение без прокручиваний и проскальзываний сферы по траектории й^еСогу(0 (аналогично функции AnimateCar).
Пример использования ПК MotionPlan-т^235. Продемонстрируем работу ПК Мойоп-Р1ашод235 на задаче о перемещении по плоскости машины с двумя прицепами. Состояние системы описывается пятью координатами (х, у, 9, фь ф2)е0=^2х^3. Здесь (х, у)еЖ.2 - координаты центра машины на плоскости; 9 - угол, задающий ориентацию машины на плоскости; ф! - угол, задающий положение первого прицепа относительно машины; ф2 - угол, задающий положение второго прицепа относительно первого прицепа. Динамика системы имеет вид х = cos(9)и1, у = sm(9)u2,
ф 1 = - sin(ф1 )и1 + (-1 - cos(ф1 ))и2
(7)
ф 2 = ^п(Ф1 -ф2) + sin(фl ))и1 + + (cos(ф1 - ф2 ) + cos(ф1 ))и2.
Зададим начальное состояние qo=(0, 0, п/4, п/4, -п/4) и целевое q1=(-0,252, -0,339, 1,085, 0,514, -1,281), в которое система должна перейти с точностью е=10-3.
На рисунках 3 и 4 приведены результаты работы модулей CPControl235 и OptimalControl235. В случае применения кусочно-постоянных управлений алгоритму потребовалось шесть итераций для достижения требуемой точности, а для опти-
=и
2
мальных управлений - четыре. При использовании кусочно-постоянных управлений система совершает больший маневр (траектория имеет большую амплитуду отклонений от целевого состояния).
В заключение следует отметить, что рассматриваемый в работе метод приближенного решения двухточечной задачи управления (1)-(2) был реализован в параллельном ПК MotionPlaning235. Комплекс испытан на двух прикладных задачах (задача о качении шара по плоскости и задача управления машиной с двумя прицепами). В случаях, когда граничные условия не были слишком далеки, комплекс успешно решал поставленную задачу управления. При далеких граничных условиях алгоритм не сходился, что соответствует теоретическому обоснованию метода (нильпотентная аппроксимация является локальным приближением исходной системы). В будущем планируется расширить функционал ПК для решения глобаль-
ной задачи управления путем ее сведения к последовательности локальных задач. В настоящее время ПК MotionPlaning235 является удобным и надежным средством решения локальной задачи (1)-(2).
Литература
1. Аграчев А.А., Сачков Ю.Л. Геометрическая теория управления. М.: Физматлит, 2005. 391 с.
2. Bellaiche A. The tangent space in sub-Riemannian geometry // Sub-Riemannian Geometry, A. Bellaiche and J.J. Risler, Eds. Basel, Swizerland: Birkh auser, 1996, pp. 1-78.
3. Laferriere G., Sussmann H.J. A differential geometric approach to motion planning // Nonholonomic Motion Planning, Zexiang Li and J.F. Canny Eds. Basel, Swizerland: Kluwer, 1992.
4. Сачков Ю.Л. Полное описание стратов Максвелла в обобщенной задаче Дидоны // Мат. сб. 2006. Т. 197. № 6. С. 111-160.
5. Laumond J.P. Lecture Notes in Control and Information Science. Springer. 1998. № 229. 343 с.
6. Murray R.M., Sastry S.S. Steering controllable systems, Proc. 29th IEEE Conf. Dec. and Control. Honolulu, Hawaii, 1990, pp. 408-412.
УДК 004.942
МОДЕЛИРОВАНИЕ РАСТВОРЕНИЯ ТВЕРДЫХ ТЕЛ С ПОМОЩЬЮ КЛЕТОЧНЫ1Х АВТОМАТОВ
(Работа выполнена при поддержке Министерства образования и науки РФ, ГК № 16.552.11.7046)
Н.В. Меньшутина, д.т.н.; С.И. Иванов; Д.Д. Шипилова (Российский химико-технологический университет им. Д.И. Менделеева, г. Москва,
patephon2009@yandex.ru)
Статья посвящена моделированию растворения твердых тел с помощью вероятностных клеточных автоматов. Приведено описание математической модели растворения и рассмотрен алгоритм работы клеточного автомата, реализующего модель. Проведено сравнение экспериментальных (кинетика растворения таблетки аскорбиновой кислоты) и расчетных данных.
Ключевые слова: клеточные автоматы, растворение, моделирование, параллельные вычисления, фармацевтическая промышленность, многокомпонентные смеси, диффузия.
Высокие темпы развития пищевой, фармацевтической и косметической промышленности в мире обусловливают необходимость решения разнообразных задач моделирования процессов в этих предметных областях. Одним из важных является процесс растворения твердых тел в различных растворителях при разных условиях.
При растворении твердых тел стоит учитывать ряд факторов, которые могут повлиять на ход процесса растворения, например: твердое тело может состоять из нескольких компонентов, имеющих разную растворимость; оно может быть покрыто оболочкой, затрудняющей растворение, а также иметь несимметричную форму; в состав твердого тела могут входить вещества, способные увеличивать растворимость других входящих в него веществ.
На текущий момент уже существуют разнообразные математические модели растворения твер-
дых тел на основе дифференциальных уравнений ^еШи11, Hixon-Crowell, Korsemeyer-Peppas и т.д.).
Данная статья посвящена описанию модели растворения твердых тел на основе вероятностного клеточного автомата. При работе над моделью авторы попытались учесть все факторы, перечисленные выше.
Описание модели клеточного автомата
В основе модели лежит вероятностный клеточный автомат [1]. Клетки автомата образуют прямоугольное поле размером ЫхЫ клеток. Модель имеет следующие допущения:
- система представляется в виде поля, состоящего из квадратных клеток;
- каждая клетка имеет четыре соседние клетки;
