CC BY
66
Вычислительные технологии
Том 18, Специальный выпуск, 2013
Моделирование жёстких гибридных систем с односторонними событиями в инструментальной
Е.А. Новиков1, Ю.В. Шорников2 1 Институт вычислительного моделирования СО РАН, Красноярск, Россия
2Новосибирский государственный технический университет, Россия e-mail: novikov@icm.krasn.ru, shornikov@inbox.ru
Рассмотрены особенности компьютерного анализа гибридных систем в инструментальной среде ИСМА. Приведены классы математических моделей непрерывного поведения гибридных систем, символьная и графическая спецификации обозначенного класса и особенности численной реализации с учётом жёсткости и односторонности событий.
Ключевые слова: гибридная система, структурная схема, импорт данных, жёсткие режимы, событийная функция, дифференциально-алгебраические уравнения, неявные задачи.
Введение
Для качественного описания довольно большого класса практических задач требуется учитывать как непрерывное, так и дискретное поведение динамических систем. Современная теория гибридных систем (ГС) является универсальным аппаратом математического описания сложных дискретно-непрерывных процессов различной физической природы. Рассматриваемые задачи спецификации и эффективного численного анализа динамических и гибридных систем относятся к категории фундаментальных [1, 2]. Программная система ИСМА [3] для эффективной реализации моделей из обозначенного класса наполнена графическими и текстовыми входными языками предметного описания таких приложений, как электромеханика, электроэнергетика, автоматика, химическая кинетика, биология и др. Эффективные оригинальные алгоритмы машинного анализа в идеологии гибридного моделирования рассматриваются для различных приложений в соответствии с выбранными классами жёстких и нежёстких задач в условиях односторонних событий.
1. Задача Коши
Класс динамических систем и непрерывные поведения ГС представлены в ИСМА задачей Коши с запаздывающим аргументом
*Работа выполенена при частичной финансовой поддержке РФФИ (грант № 11-01-00106-а и № 14-01-00047) и Президиума РАН (проект № 15.4).
среде ИСМА
*
X = f [x (t) , x (t — t), t] , x (t0) = x0,
(1)
где x(t) = p(t) при t G [t0 — т, to), x G Rn — вектор состояния, ^ G Rr — вектор-функция запаздывания; r < n; t — независимая переменная; t G [t0, ]; т = {ti, ..., тг}T — вектор чистых запаздываний; f : R x Rn ^ Rn — нелинейная вектор-функция, удовлетворяющая условию Липшица. В отличие от передовых мировых аналогов, например системы Simulink/MATLAB, особенностью ИСМА является возможность импорта данных непосредственно в программную модель на оригинальном графическом языке структурных схем [4]. Многие пользователи подготавливают и хранят данные в файлах внешних приложений формата MS Excel. Это связано в первую очередь с удобством представления и обработки данных. В ИСМА реализованы средства импорта массива точек из Excel, что позволяет свободно манипулировать данными. Импорт осуществляется через интерфейс нелинейного блока. В качестве иллюстрации импорта данных рассмотрим задачу моделирования электропогрузчика [5]. Программная модель тягового электропривода электропогрузчика представлена на рис. 1. Здесь выделенный нелинейный блок f (x) реализует изменение момента двигателя. Изменения регистрировались с помощью специальной аппаратуры.
После обработки данных в Excel потребовалось в нелинейную функцию модели ввести массив из 5000 точек в формате внешнего приложения. Массив экспериментальных точек из внешнего приложения MS Excel введён программно [4] как последовательность
Рис. 1. Программная модель электропривода
б
а
Рис. 2. Изменение момента двигателя (а) и переходные процессы двигателя (б)
значений динамики момента двигателя (рис. 2, а). Результаты машинных экспериментов с импортированными в компьютерную модель данными приведены на рис. 2, б.
Под этот же класс систем (1) подпадают задачи химической кинетики. Для спецификации задач химической кинетики разработан язык Ы8МА+ [6], являющийся расширением базового языка спецификации динамических и гибридных систем ЫБМА [7]. Представим порождающую грамматику химических реакций О [С] в виде С ^ Б С |Б и Б ^ EsaEs, где Б — стадии химической реакции; Es — подмножество арифметических выражений; а Е X* — символ итерации терминального алфавита основной грамматики О [Е]. Символ итерации однозначно определяется терминальной цепочкой а ^ ¿^ > |е с идентификатором ¿^ соответствующей скорости стадий. Значение Ев выбирается так, чтобы выполнялось условие О [Ев] С О[Е]. Тогда необходимо, чтобы продукции новой грамматики имели вид Ез ^ Т |Т + Ез, Т ^ О |О * Т и О ^ ¿^ |с. Здесь г^ — идентификатор реагента химической реакции, представляет собой запись переменной без индекса и поэтому не противоречит общепринятой в О[Е] записи простых переменных; с = сош1 или целое без знака, которое означает число реагентов в реакции. С учётом вложенности для расширенной грамматики О [Ев] используются наследуемые однозначные методы анализа, что и для базовой грамматики О[Е]. Поэтому не требуется доработки языкового процессора системы ИСМА в целом. Этим примером показано важное решение задачи унификации математического и программного обеспечения в рамках разработанной инструментальной среды.
2. Дифференциально-алгебраические задачи
При описании ГС с непрерывным поведением в классе дифференциально-алгебраических уравнений, разрешённых относительно производной, ограничимся уравнениями вида
у' = / (¿,Х,У) , ж = V (¿,Х,У) , Рг : 9 (¿,Х,У) < О,
Ь Е [¿с,4], х Е ЯМх, у Е ЯМу, х (¿о) = же, у (¿о) = уо,
/ : Я х ЯМх х ЯМг ^ ЯЫу, V : Я х Ямх х ЯЫу ^ ЯМх,
9 : Я х ЯИх х ЯИу ^ Я3, Б < . (2)
В (2) использованы принятые при описании ГС обозначения предиката рг, событийной функции 9(Ь,ж,у) и т.д. Существует множество приложений в этом классе, в том числе задачи диффузии, описываемые уравнениями в частных производных. Например, модель конкуренции Лотки — Вольтерра на основе системы реакции — диффузии в двухмерном пространстве представлена уравнениями [8]
Ж = *®т + Ш + f (c1-2), (3)
где di = 0.05, d2 = 1.0, an = 106, ai2 = 1, «21 = 106 — 1, «22 = 106, bi = 62 = 106 — 1 + 10-6, f1 (c1, c2) = c1 (61 — a11c1 — a12c2), f2 (c1, c2) = c2 (62 — a21c1 — a22c2). Граничные условия SC/Sx = 0 при x = 0, x =1 и dc/dz = 0 при z = 0, z =1.8. Начальные условия имеют вид
c1 (x, z, 0) = 500 + 250 cos (nx) cos (10nz/1.8), c2 (x, z, 0) = 200 + 150 cos (10nx) cos (nz/1.8).
// Начальные условия у1 = 7.5000Е+2; у2 = 3.5000Е+2; // ...
// Дифф. уравнения
у11 = 8.0000Е-1*(уЗ -2.О *у1+уЗ)+2.4691Е-1*(у11-2.О *у1+у11) +
у1*( (1.0Е6-1.0+1.0Е-6)-1.0Еб*у1-1.0*у2) ; у21 = 1.бОООЕ+1*(у4-2.0*у2+у4)+4.9383Е+0*(у12-2.0*у2+у12)+ у2*( (1.0Е6-1.0+1.0Е-6)-(1.0Е6-1.0)*у1-1.ОЕ6*у2) ;
// ...
Рис. 3. Фрагмент текстовой модели
С учётом приведённых значений параметров и начальных условий решение этой реакционно-диффузионной системы сходится при £ ^ то к решению с1 = е\ = 1 — 10-6, с2 = с2 = 10-6.
Перейдём к сетке размером 3 х К, соответственно по х и г получаем шаги Дх = 1/(3 — 1) и Дг = 1.8/(К — 1), где с.к — аппроксимация сг (х., гк,£), х. = (3 — 1) Дх, гк = (к — 1)Дг, 1 < 3 < 3, 1 < к < К. В результате получим систему уравнений размерности N = 23К, причем
в в
в (сг _ 2С + С А I Дх2 \С^+1,к 2С.к + 4-1,^ + Дг2
где 1 < г < 2, 1 < 3 < 3, 1 < к < К, /]к = /г (с]к, с2к). Граничные условия на сетке
имеют вид с0,к = с2,к, с!?+1,к = сЬ-1,к для 1 < к < К и . = с5,2 , с5,к+1 = с5,к-1 для 1 < 3 < 3.
При построении компьютерной модели от тройного индекса г, 3 и к перейдём к индексу т = г + 2 (3 — 1) + 23 (к — 1). Получим текстовую модель (фрагмент на рис. 3). Система (3) рассчитывалась при размерностях от 50 до 1800 уравнений. За разумное время решение удалось получить только явными и полуявными схемами из библиотеки оригинальных методов (см. таблицу).
С]к д г2 (С.?+1,к 2с.к + с.-1,к) + Л г2 (С.?,к+1 2с?к + Сз,к-\) + .,
Таблица 1. Численные методы в системе ИСМА
Алгоритм Характеристика
БШРР Алгоритм переменного порядка с контролем устойчивости (максимальный порядок точности 5)
ИЛБЛЩ Неявный метод анализа жёстких режимов
БШРР1_, ИЛБЛЩ Метод БШРР в комбинации с ИЛБЛЩ
БТЕКБ Явный метод четвёртого порядка с контролем устойчивости на основе метода Мерсона
БР78БТ Явный метод восьмого порядка с контролем устойчивости на основе метода Дорманда—Принса
БШРР Алгоритм переменного порядка с контролем устойчивости (максимальный порядок точности 3)
ИК28Т, ИКЗБТ Явные методы 2-го и 3-го порядка с контролем точности и устойчивости
МК22, МК21 Второй порядок точности, "замораживание" матрицы Якоби, жёсткие режимы
Особенность решения состоит в том, что в начале интервала интегрирования оно меняется быстро, а далее — очень медленно. Из явных схем переменного порядка и шага с контролем точности и устойчивости вычислений наилучшие результаты показал алгоритм Б18РЕ, который использует метод пятого порядка в начале интегрирования и ше-стистадийный метод первого порядка с расширенной областью устойчивости на участке с медленно меняющимся решением. Полуявные ¿-устойчивые методы МК22 и МК21 более эффективны, чем явные схемы. Их эффективность обусловлена ¿-устойчивостью и замораживанием матрицы Якоби.
Таким образом, традиционно используемые неявные схемы для исследования жёстких режимов могут оказаться бесполезными по двум причинам. Во-первых, в задачах высокой размерности для декомпозиции матрицы Якоби может потребоваться недоступно большое время для современных процессоров. Во-вторых, как показали проведённые исследования [9], если жёсткая модель высокой размерности является гибридной, применение неявных схем приводит к неверному глобальному решению. В связи с этим для эффективного анализа таких систем в библиотеку методов ИСМА включены разработанные оригинальные явные схемы с контролем устойчивости и ¿-устойчивые полуявные методы, которые дополнены алгоритмом корректного обнаружения событий.
3. Неявные задачи
При моделировании электрических цепей, электроэнергетических процессов и во многих других приложениях возникает необходимость численного решения жёстких систем дифференциально-алгебраических уравнений, неразрешённых относительно производной:
^ (ж, ж', ¿) = О, рг : 9 (ж, ж', ¿) < О, Ь Е [¿о,4] , ж (¿о) = жо, ^ : Ям х Ям х Я ^ Ям, 9 : Ям х Ям х Я ^ Я. (4)
Современные численные методы обычно предполагают задание явной зависимости производной от решения. Приведение неявной задачи к разрешённому виду порождает дополнительные вычислительные затраты. В библиотеку ИСМА включён оригинальный ¿-устойчивый алгоритм решения неявных задач [10] на основе схемы типа Розен-брока. В данном алгоритме разрешение задачи относительно производной и обеспечение ¿-устойчивости осуществляется одновременно. Задачу (4) можно записать в виде
•' = у, ^(ж,у,Ь) = 0, ж(Ьо) = жо, у(Ьо) = Уо.
ж
Дополнительное условие у(Ьо) = уо можно вычислить, например, решая задачу ^(жо, у) = 0 на установление. Тогда метод типа Розенброка для решения этой задачи записывается в виде
жп+1 жп + кХ, уп+1 = уп + ку ,
ад = ^¿пу - ¿п), ку = О- (кх - ^уп),
^п ¿пу + аh¿nж, ¿пу д¿(yn, жn, Ьп)/ду
¿пх = ^(уп, жп, ¿п)/дж. (5)
Для контроля точности схемы (5) на каждом шаге проверяется неравенство ||кХ11 — £, где £ — требуемая точность расчётов; || • || — некоторая норма Ям. В отличие от методов типа Розенброка применительно к решению разрешённой задачи производная решения в (5) вычисляется приближённо. Поэтому при выборе величины шага интегрирования дополнительно проверяется неравенство 1^га|| — £.
Важной проблемой в моделировании ГС является обнаружение смены режимов. В библиотеку ИСМА включён алгоритм корректного обнаружения событий, в котором наряду с точностью и устойчивостью вычислений учитываются динамика событийной функции и односторонность событий, которая имеет место в большинстве практических задач. Оригинальный метод локализации точек переключения основывается на доказанной теореме о выборе шага интегрирования по формуле [10]
^+1 = (7 - 1) 9п ! (^ПVп + ^^ , 7 е (0,1).
Такой способ выбора шага с одновременным учётом величины шага по точности и устойчивости обеспечивает асимптотическое приближение к границе режима.
Для иллюстрации работы алгоритма корректного обнаружения событий приведём результаты моделирования простой гибридной системы — прыгающего с неупругим отскоком мячика. Режимное поведение зададим неявной системой дифференциальных уравнений из класса (4) в виде
у' — V = 0, V + а = 0, рг : —у < 0,
где y — расстояние мячика от поверхности отскока; v — скорость падения мячика; a — ускорение свободного падения. В момент отскока y = 0 и v = —av, где а — коэффициент сохранения скорости. На рис. 4 представлены зависимости координаты мячика от времени, обработанные графическим интерпретатором GRIN среды ИСМА.
При расчётах без контроля динамики событийной функции (рис. 4, а) возникает существенная ошибка в обнаружении события. Это приводит к нарушению условия односторонности и, как следствие, ошибочному глобальному решению. Использование алгоритма для асимптотического приближения к границе режима (рис. 4, б) обеспечивает точное обнаружение момента смены режима. При приближении к поверхности y = 0 происходит уменьшение шага интегрирования, а при удалении от границы режима шаг
а б
4 5 6 t 4 5 б t
Рис. 4. Результаты моделирования: а — без учёта динамики событийной функции, б — с использованием алгоритма обнаружения событий
определяется по критериям точности расчётов и устойчивости численной схемы. Этим свойством представленный алгоритм детекции выгодно отличается от разработанных ранее алгоритмов без учёта динамики событийной функции [6].
Заключение
Инструментальная среда ИСМА с оригинальной библиотекой численных методов и предметно-ориентированными средствами символьного и графического описания моделей из класса (1), (2) и (4) адекватно настроена для исследования гетерогенных гибридных систем в довольно широком диапазоне математических моделей с учётом жёсткости и в условиях односторонних событий.
Список литературы
[1] Моисеев Н.Н. Математические задачи системного анализа. М.: Наука, 1981. 488 с.
[2] Яненко Н.Н., Карначук В.И., Коновалов А.Н. Проблемы математической технологии // Числ. методы механики сплошной среды. Новосибирск: ВЦ СО АН СССР. 1977. Т. 8, № 3. С. 129-157.
[3] Шорников Ю.В., Дружинин В.С., Макаров Н.А. и др. Инструментальные средства машинного анализа: Свидетельство об официальной регистрации программы для ЭВМ № 2005610126. М.: Роспатент, 2005.
[4] Шорников Ю.В., Дружинин В.С. Импорт данных в программной среде: Свидетельство об официальной регистрации программы для ЭВМ № 50200600117. М.: ВНТИЦ, 2006.
[5] Аносов В.Н., Кавешников В.М., Шорников Ю.В. Характеристики управляющих воздействий тягового электропривода автономного напольного транспортного средства // Науч. вест. НГТУ. 2005. № 3(21). С. 37-44.
[6] Новиков Е.А., Шорников Ю.В. Компьютерное моделирование жёстких гибридных систем. Новосибирск: Изд-во НГТУ, 2012. 451 с.
[7] Шорников Ю.В., Томилов И.Н. Программа языкового процессора с языка LISMA (Language of ISMA): Свидетельство об официальной регистрации программы для ЭВМ № 2007611024. М.: Роспатент, 2007.
[8] Brown P., Hindmarsh A. Matrix Free Methods in the Solution of Stiff Systems of ODEs. Lawrence Livermore National Laboratory, 1983. 38 p.
[9] Novikov E.A., Shornikov Yu.V., Tomilov I.N. et al. Modeling stiff hybrid systems of high dimension in ISMA // Proc. IASTED Intern. Conf. on Automation, Control and Information Technology (ACIT 2010). Novosibirsk, Russia: ACTA Press, 2010. P. 256-260.
[10] Shornikov Yu.V., Dostovalov D.N., Myssak M.S. Simulation of hybrid systems with implicitly specified modal behavior in the ISMA Environment // Универ. науч. журн. 2013. № 5. С. 171-178.
Поступила в 'редакцию 29 ноября 2013 г.
