CC BY
20
»
Математическое моделирование: методы, алгоритмы, технологии
УДК 681.3.06
К.Ю. Алтунин, Ю.Б. Сениченков
О ВОЗМОЖНОСТИ СИМВОЛЬНЫХ ВЫЧИСЛЕНИЙ В ПАКЕТАХ ВИЗУАЛЬНОГО МОДЕЛИРОВАНИЯ СЛОЖНЫХ ДИНАМИЧЕСКИХ СИСТЕМ
В пакетах визуального моделирования сложных динамических систем решаются две основные задачи: автоматически строится математическая модель исследуемой или проектируемой системы в виде уравнений различного типа по заданному графическому описанию и численно находится решение полученной системы уравнений.
Для описания моделируемых систем в МуБшсНит [10] используется графический язык иерархических функциональных схем с блоками ориентированными (т. е. с входами-выходами) и неориентированными (с контактами и потоками), внутреннее поведение которых представляется картами поведения (гибридными автоматами). Математические модели пакета строятся в виде систем алгебро-дифференциальных, дифференциальных и алгебраических уравнений.
В современных пакетах визуального моделирования более свободными становятся формы описания проектируемых систем и все чаше автоматически строятся различные аппроксимации исходной системы, которые затем преобразуются к виду, облегчающему численное решение с заданной точностью. Иными словами, от принципа "строим только то, что удобно вычислительной машине", переходим к принципу "делаем то, что удобно разработчику, автоматизируя рутинные операции".
Пакет МуБшсИшп практически не ограничивает форму записи локальных систем уравнений в картах поведений и позволяет пользователю использовать широкий набор численных методов с различным внутрен-
ним, невидимым пользователю, интерфейсом, что часто требует автоматического преобразования исходной системы к заданному виду. При автоматическом построении итоговой системы для функциональной схемы, содержащей неориентированные блоки, часто возникает задача понижения индекса построенной системы алгебро-дифферен-циальных уравнений. Ее решение требует дифференцирования уравнений.
Преобразование систем от формы, удобной пользователю, к форме, необходимой пакету, невозможно выполнить без применения символьных вычислений. Наиболее яркая иллюстрация использования символьных вычислений в моделировании — пакет MapleSim. который можно рассматривать как пакет Maple, дополненный графическим языком "физического" моделирования. Символьные вычисления пакета Maple используются как на этапе построения системы, так и для ее преобразования и подготовки к численному решению. Можно сомневаться в необходимости столь радикального подхода к использованию символьных вычислений при моделировании сложных динамических систем, но целесообразность использования символьных вычислений в современных численных методах стала очевидной для разработчиков пакетов моделирования. В этой статье мы будем говорить только о символьном дифференцировании.
При проектировании пакета MvStudium предполагалось, что в нем будут использованы только хорошо известные численные методы, доступные пользователю на уровне
исходных текстов. К ним относятся библиотеки СЮЕРАСК, программы из книг [8—10]. Эти программы действительно высокоэффективны, но для них характерна "ручная" подготовка системы к решению, недопустимая в пакетах визуального моделирования. Например, в случае, если исходная система жесткая, приходится применять неявные методы, требующие построения матрицы Якоби.
Пакет МуБшсПит позволяет находить оптимальные параметры для моделей, что предполагает использование численных методов оптимизации, часть из которых требует вычисления матрицы Гессе.
Перечисленных примеров достаточно, чтобы оправдать попытку использовать символьное дифференцирование в пакетах моделирования.
Символьное вычисление производных
Вычисление значений производных функций может осуществляться численно на основе различных разностных схем (численное вычисление производных). Это наиболее часто используемый вариант в численных алгоритмах, и его недостатки хорошо известны. Производные функций могут быть построены явно с помощью символьных вычислений, а затем вычислены при заданных значениях аргументов (символьное вычисление производных). Наконец, можно проводить символьное построение выражения для производной, но не сохранять производную в виде выражения, а непосредственно вычислять значение производной при заданных значениях аргументов (автоматическое дифференцирование). Одним из преимуществ последнего подхода является возможность получения производной функции, для описания которой используется не математический, а алгоритмический язык. Эту разновидность автоматического дифференцирования называют "дифференцированием алгоритмов". При автоматическом дифференцировании могут использоваться обычные правила Лейбница, применяемые к исходной функции в любой форме записи, но возможно также и использование рядов Тейлора.
Пакет МуБШ<1шт разрабатывается уже 15 лет. постоянно модифицируется, и по-
этому речь может идти только о встраивании нового блока "символьное дифференцирование" в уже существующий программный продукт.
Использование таких пакетов, как Maple или Mathematica, в виде подпрограмм в пакетах моделирования уже обсуждалось на примере пакета MapleSim.
Можно было бы использовать в виде подпрограмм следующие уже созданные пакеты для автоматического дифференцирования.
Пакет для автоматического дифференцирования алгоритмов, написанных на языке С/ С++. [1]. Пакет использует метод, определяемый авторами как автоматическое дифференцирование с помощью переопределения операций в языке С++ (ADOL-C). Исходная функция представлена в виде кода на языке С++, для вычисления производной используется ее разложение в ряд Тейлора. Таким образом, в пакете ADOL-C производная в виде выражения недоступна.
Пакет Jacal [3]. Не обсуждая возможности символьных вычислений пакета, отметим. что он написан на языке программирования Scheme, что вызывает сложности при стыковке с другими программами. Символьные вычисления в современных пакетах моделирования используются как на этапе компиляции, так и на этапе исполнения. Использование кода для языка Scheme на этапе исполнения скорее всего приведет к медленной скорости работы.
Пакет DiffExpress [4]. Программа позволяет дифференцировать выражения в символьном виде и поддерживает все необходимые математические функции и операции. Пакет не позволяет вычислять значение функций и производных в конкретной точке и работать с системами уравнений. Нет возможности вычислять частные производные, вычисления возможны только для одной переменной.
Пакет Symbolic Differentiation [5]. Исходные выражения дифференцируемой функции хранятся в польской записи, где все операции, переменные и константы в свою очередь хранятся в стеке. На верху стека лежит операция, которую необходимо выполнить для получения значения выражения, следующие значения стека являются аргументами
этой операции. Вычисления производной (верхняя операция из стека) производится согласно заранее построенной таблице производных, а в качестве аргумента выступает следующий элемент стека. В предложенном здесь подходе не рассматривается возможность вычисления частных производных и отсутствует возможность вычисления значения производной в конкретной точке.
Даже этот неполный обзор показывает, что пока еще пакеты автоматического дифференцирования далеки от совершенства. Целесообразно разработать свой пакет "символьное дифференцирование", в котором будут учтены особенности пакета MvStudium — его методы хранения информации, алгоритмы построения итоговой системы уравнений и численные методы ее решения.
Блок "Символьные вычисления" пакета MvStudium
Блок символьного дифференцирования, для того чтобы его можно было применять в пакете MvStudium, должен обладать следующими возможностями:
уметь обрабатывать выражения, содержащие как скалярные переменные, так и векторно-матричные;
вычислять не только значения производной заданной функции, но и строить выражение для производной в символьном виде;
вычислять градиент, матрицу Якоби, матрицу Гессе;
обеспечивать высокую скорость вычислений. В существующих программах дифференцирования используется метод автоматического дифференцирования, который позволяет получить значение производной для конкретной точки, но не предусматривает получение производной в явном виде. Для применения же в пакете MVS в модуле дифференцирования необходимо вычислять не только значение производной для конкретных точек, но и получать производную в явном виде. Данная возможность необходима для понижения индекса системы. Построение в явном виде производной происходит по тому же принципу, по которому осуществляется автоматическое дифференцирование. Вычисление производной ведется по правилам Лейбница, за-
ранее определенным в системе. Выражение для дифференцирования в программе представлено как бинарное дерево, узлы которого определяют операции, константы и переменные. Каждый узел дерева является простейшим элементом выражения. Узел N содержит в себе два дочерних узла, N/ф и Nngh„ которые являются левым и правым операндом. Для бинарных операций узел N есть отображение N( N,efn Nright) -> {Nresull, ^'result)■> где /V Nright, Nresu„ - тройка непустых множеств и N с Nk,ft х Nrighl. Каждый узел воплощает в себе правило Лейбница для дифференцирования операции, представленной этим узлом, и правила для упрощения. В общем случае узел представляет собой элементы Nlefh Nrigh„ Rllif, и {Rsiinpii/y}, где Rdif — правило Лейбница для дифференцирования; Rcalc — операция; {^simplify} ~ множество правил упрощения. N = {N,ef,,
Nright, Rdif> Kale. Simplify}- Для унарных опе-раций Nieft = null, где null — пустое множество. Переменные и константы — частные случаи узла — представлены в следующей форме: NX = [null, null, NX' = I, x0, null}-, NC= {null, null, NC'= 0, value, null), где NX -узел переменной; NC — узел константы, Value — значение константы, х{) — значение переменной в заданной точке, которое используется при вычислении выражения.
В нашем модуле при вычислении производной в явном виде либо — при использовании подхода автоматического дифференцирования — в конкретной точке осуществляется обход дерева в глубину; при этом находится производная для каждого поддерева, и в итоге мы получаем производную для всего выражения. В существующих подходах применяется метод дифференцирования алгоритмов, где вводящий — новый тип данных, у которого переопределены все операции. Для каждой операции у нового типа данных прописан алгоритм дифференцирования, однозначно согласно правилам Лейбница определяющий производную от каждой операции. В нашем случае данный подход теряет всякий смысл, так как неизвестно выражение, которое необходимо дифференцировать до момента исполнения. И при переопределении операций для нового типа данных, как это делается сейчас
во многих пакетах, мы не можем формировать уравнения на этапе исполнения программы (это невозможно без изменения структуры внутренней работы MVS). Поэтому в нашем модуле применен тот же метод дифференцирования алгоритмов, что и во многих пакетах автоматического дифференцирования, но дифференцирование операций происходит не для конкретного типа, а каждая операция представляет собой новый тип данных. Такой метод позволяет строить уравнения в ходе выполнения программы и выполнять дифференцирование выражений, заранее не определенных в программе.
Существенный минус при вычислении производной по правилам Лейбница — это наличие множества подобных членов, а также констант, которые задают множество операций, поддающихся упрощению. Для того чтобы избежать громоздких выражений, в модуле определен блок, который упрощает полученное выражение по правилам, заранее заданным в системе. Так, например, происходит приведение подобных членов, осуществляются численные операции с константами там, где это возможно. В результате мы получаем выражение, при работе с которым в дальнейшем, вычислительные затраты будут минимальны. Работа модуля упрощения построена таким образом, что после вычисления производной в явном виде осуществляется обход дерева в ширину. И каждый узел дерева осуществляет проверку своих поддеревьев по определенным в узле правилам. Так, например, вычисляются все операции с константами, которые можно вычислить в выражении, происходит отсечение поддеревьев при умножении на ноль и так далее.
Таким образом, нахождение производной представляет собой обход дерева выражения. Если необходимо найти значение производной в заданной точке, то модулю указывается эта точка и осуществляется обход дерева, но в этом случае правила Лейбница применяются с учетом конкретных числовых значений. В результате мы получаем численное значение производной в конкретно указанной точке.
Модуль позволяет вычислять матрицу Якоби двумя способами.
1. Вычисление якобиана происходит каждый раз, когда это необходимо; операции производятся с числами по правилам Лейбница: упрощение выражений не выполняется.
2. Вычисление якобиана производится при инициализации системы, т. е. один единственный раз, в начале работы модуля. В таком случае в памяти приходится хранить не только исходную систему, но и матрицу Якоби, представленную символьно в виде деревьев выражений. При расчете матрицы Якоби после каждого вычисления производной производится упрощение дерева полученного выражения согласно правилам, описанным для каждого узла дерева. Когда в численном методе происходит вызов функции для подсчета якобиана, модуль не вычисляет производные, а использует уже найденную матрицу производных. В таком случае не приходится рассчитывать производные каждый раз, и мы экономим на скорости вычислений. Так как приходится хранить матрицу Якоби, то возрастают расходы оперативной памяти на хранение системы.
В модуле учитывается информация о структуре системы, которую MVS передает в блоки численных методов. Данная информация указывает, какие переменные входят в конкретные уравнения системы. По этой информации модуль строит матрицу А вхождения независимых переменных в уравнения системы:
А = [а0 [ / = 1... т. j = 1 ... и.
Матрица А в памяти компьютера хранится как пары значений [/', j\, для которых djj — 1. Таким образом, вычисления
о « ди,
значении для матрицы Якоби —- произ-
дх,
водятся лишь в том случае, если a,j = 1. В результате экономим время вычислений и расход памяти.
Для интеграции с пакетом MVS в модуле символьного дифференцирования предусмотрена загрузка уравнений несколькими способами:
загрузка из файла rnvl. который представляет собой описание моделируемой системы на языке Model Vision;
загрузка из строки, содержащей xml представление уравнения. В качестве узлов xml выступают функции, операции, переменные и константы. В результате возможен простой обмен уравнениями между MVS и модулем дифференцирования.
Экспериментальная проверка
Для проверки возможностей блока символьного дифференцирования была использована программа Radaub [9]. Предусмотренные в ней подпрограммы построения матрицы Якоби создавались автоматически с использованием блока "символьное дифференцирование".
Задача "Загрязнение окружающей среды" входит в набор тестов для решателей дифференциальных уравнений [6]. Модель использует химические реакции, приводящие к загрязнению воздуха. Она разработана Голландским национальным институтом здоровья и зашиты окружающей среды (RIVM) и описана Вервером [7]. Модель состоит из 25 реакций и 20 смесей, полученных в результате реакции. Уравнения модели являются жесткой системой из двадцати нелинейных обыкновенных дифференциальных уравнений.
Сравнивались следующие режимы работы программы:
Режим 1: матрица Якоби была вычислена вручную и оформлена в виде подпрограммы.
Режим 2: матрицы Якоби вычислялась численно программной Radau5
Режим 3: матрица Якоби вычислялась численно пакетом "символьное дифференцирование"
Режим 4: матрица Якоби вычислялась символьно блоком "символьное дифференцирование" в символьном виде
Режим 5: при вычислении матрицы Якоби блоком "символьное дифференцирование" учитывалась ее разреженность.
Результаты вычислений представлены в табл. 1.
Эти результаты показывают, что символьные вычисления могут приводить к большим потерям производительности (более пятнадцати раз для режима 4), но положение можно исправить, если учесть разреженность.
Таблица I
Символьные вычисления матрицы Якоби
Метод Время Относи- Абсолют- Относи-
вычис- вычис- тельная ная по- тельное
ления ления, погреш- грешность время
матрицы мс ность вы- вычисле- выполне-
Якоби числения, % ния, % ния, %
Режим 1 4,5515 0 0 100,00
Режим 2 4,8808 7,674-Ю-'5 9,39-10 17 107,24
Режим 3 4,9844 7,674-10 15 9,39-10 17 109,51
Режим 4 69,831» 0 0 1534,24
Режим 5 29,6289 0 0 650,97
В методе Яаиаи5 предусмотрены вычисления с разреженными матрицами. В этом случае в исходном коде либо матрица Якоби считается численно, либо требуется ее ручное вычисление в символьном виде. Мы повторили наши эксперименты для систем "разреженных" дифференциальных уравнений. Сравнение различных режимов вычислений проводилось на специально сконструированной тестовой системе линейных дифференциальных уравнений 400-го порядка, с трехдиагональной матрицей. Результаты вычислений представлены в табл. 2.
Таблица 2 Вычисления разреженной системы
Метод Время выполнения 100 шагов Radau5, с Время выполнения 100 шагов Radau5, % Расход памяти, Мб
Исходный вариант программы Ка11аи5 с численным нахождением матрицы Якоби 16,7 100 12,8
Использованы символьные вычисления при нахождении матрицы Якоби 12,4 75,3 13,7
Учёт разреженности системы. Численное вычисление матрицы Якоби 0,74 4,4 7,8
Символьное вычисление матрицы Якоби 1,2 7,1 8,5
Учёт разреженности системы. Символьное вычисление матрицы Якоби 1 6 8,7
Из табл. 2 видно, насколько можно увеличить быстродействие при учете разреженной структуры матрицы Якоби. На конкретном примере при решении системы обыкновенных дифференциальных уравнений при матрице Якоби, заполненной на 0,7 %, выигрыш в скорости относительно полного расчета матрицы составил 95,6 %. При использовании символьного вычисления матрицы Якоби прирост скорости составил 92,9 %.
Данную работу следует рассматривать как первый практический шаг на пути использования символьных вычислений производных в численных методах решения обыкновенных дифференциальных уравнений. Помимо совершенствования алгоритмов дифференцирования, будет продолжен поиск практически важных приложений, где оправдано использование данного метода.
СПИСОК ЛИТЕРАТУРЫ
1. Walther A„ Kowarz A., Griewank A. A Package for the Automatic Differentiation of Algorithms Written in C/C++ Version 1.10.0, July 2005.
2. Gesslein G. Mathomatic — automatic algebraic manipulator // http://mathomatic.org/math/.
3. JafTer A.. Thomas M.. Hedden J.D. Jacal // http://people.csail.mit.edu/jalTer/JACAL.
4. Aid Aim Software «DiffExpress» // http:// www.aidaim.com/products/free_components/ difr_spc.php# Description.
5. Hatem M. Symbolic Differentiation // http:// www.codeproject.com/KB/recipes/ Differentiation.aspx.
6. Houwen P.J. van der, Hoffmann W., Sommeijer B.P. Test Set for Initial Value Problems Solvers // http://ftp.cwi.nl/IVPtestset/index.litm.
7. Verwer J.G. Gauss-Seidel iteration for stiff ODEs from chemical kinetics. SIAM J. Sci. Comput. 1994. № 15(5): P. 1243-1259.
8. Хайрер Э., Нерсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи. М.: Мир. 1990.
9. Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи. М.: Мир. 1999. 685 с.
10. Колесов Ю.Б., Сениченков Ю.Б. Моделирование систем: Т. 1-3. СПб.: БХВ. 2007-2008.
УДК 004.94
Е.В. Бочкарева, Л.И. Сучкова, А.Г. Якунин
ПРИМЕНЕНИЕ ИМИТАЦИОННОГО МОДЕЛИРОВАНИЯ ДЛЯ ИССЛЕДОВАНИЯ ПРОЦЕССОВ СБОРА И ОБРАБОТКИ ДАННЫХ МИКРОКОНТРОЛЛЕРНЫМИ УСТРОЙСТВАМИ
В области разработки распределенных вычислительных систем с гетерогенной структурой перспективно предварительное моделирование их функционирования, позволяющее выявить "узкие места" и проверить применимость алгоритмов с различной временной сложностью для устройств с ограниченными возможностями. Событийно-ориентированное моделирование, основанное на логике параллельной работы процессов системы, служит инструментом для исследования работы систем диспетчер-
ского управления и сбора данных. Эти системы реализуют методы автоматизированного управления сложными динамическими объектами и включают подсистемы управления, сбора, обработки, передачи, хранения и отображения информации [1|.
Как правило, при проектировании гетерогенных систем сбора и обработки данных можно выделить три уровня их архитектуры [2|.
Нижний уровень включает различные датчики, собирающие и поставляющие информацию локальным интеллектуальным
