Компьютерный анализ многомерных гибридных систем явными методами в Исма Текст научной статьи по специальности «Математика»

Компьютерный анализ многомерных гибридных систем явными методами в ИСМА

Ю. В. Шорников, И. Н. Томилов, М. С. Денисов

Новосибирский государственный технический университет, 630092, Новосибирск, Россия

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

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

In this paper the high-dimensional problem of chemical kinetics has been considered. Performance analysis of explicit methods with variable order and stability control has been made. High effective algorithm for the problem solving has been detected.

Key words: hybrid system, visually-linguistic modeling, localization of switch points, array, the algorithmic determination of the system ODE.

Введение. Интерес к исследованию многомерных гибридных систем (ГС) явными методами в ИСМА обусловлен следующими факторами. В настоящие время сложность объектов моделирования постоянно увеличивается. Поведение сложного объекта или процесса часто описывается системой дифференциальных уравнений (ДУ) большой размерности. Существует класс задач, в которых задание правых частей системы ДУ носит алгоритмический характер, т. е. существует непосредственная зависимость между индексом дифференциального уравнения и его правой частью (подобная форма описания ДУ распространена в задачах химической кинетики). Как правило, системы ДУ большой размерности обладают свойством жесткости. Традиционно для решения жестких задач применяются неявные методы, так как они имеют лучшие показатели производительности по сравнению с классическими явными методами. Это обусловлено тем, что на участках установления [1] жесткой задачи шаг интегрирования явной схемы ограничен, прежде всего, условием устойчивости. В настоящей работе рассматривается применение явных методов переменного порядка и шага с контролем устойчивости к жестким задачам большой размерности. В качестве примера рассматривается гибридная задача химической кинетики, сформулированная в лаборатории Akzo Nobel Central Research [2].

Описание системы и математическая модель. Лаборатория Akzo Nobel Central Research сформулировала задачу о проникновении помеченных радиоактивной меткой антител в пораженную опухолью ткань живого организма. Исследование проводилось в диагностических и терапевтических целях. Рассматривается система одномерных уравнений реакции-диффузии [2]

ди д 2и

- kuv,

Ы дх2 ду

— = -киу.

дл

Дискретизации производных первого и второго порядков по пространственной переменной соответственно имеют вид

ди_ = _ - и}-1 ЫЧ = _ - 2и} + _ 1 < _ < н дС 2К ' д^2 (А^)2 ' '

Значения ^ , им+1, полученные из граничных условий, имеют вид и0 = (р(1) и им+1 = им . В предположении у = (и1,У1,и2,у2,...,ин,Ун) и 1к = 20 , эта задача записывается в виде

йг

= / (г, у), у (0 ) = 2, у е Я 2" ,0 < г < 20

(1)

(N - задаваемый пользователем параметр). Функция / определяется формулами

г у 2 /+1 - у 2/-3 , о у 2/-3 - 2у 2 /-1 + у 2 /+1 ,

/2/-1 = \ ' + ^-7--^-~ - кУ2-У 2/,

2А^

/г / = -ку 2/у 2/-1,

(2)

где

а

= 2(А^ -1)3 с2, в/ = (/АС -1)4 с2,1 < / < N.

4 2

К = —, У-1 (г) = ф(г), У2N+1 = У2N

(3)

<(г )=

2 е Я2N,2 = (0,^,0,^,...,0,^) . Функция <р(г) задается в виде

[2, г е [0,5], [0, г е (5,20].

Для параметров к , У0, с подходящими значениями являются к = 100, У0 = 1, с = 4. Рассматривается случай для N = 200 , т. е. система (1) состоит из 400 уравнений.

Особенности спецификации задачи. Существование разрыва первого рода функции <(г) подразумевает наличие двух режимов функционирования системы. Переключение происходит при г > 5 : рг1 : 2(г) = г > 5 . Карта поведения системы показана на рис. 1.

В ИСМА возможны два варианта реализации смены режима: использование нелинейной функции и описание локального состояния на языке ЬШМЛ, представленные на рис. 2, 3 соответственно. Первый вариант является более очевидным, так как пользователь имеет возможность интерактивного задания требуемой нелинейности и оценки ее визуального отображения. Однако в случае, когда переход инициируется выходным сигналом нелинейной функции, применение алгоритмов точной локализации точек переключения ГС невозможно. Как показано ниже, этот момент является ключевым, поэтому при построении компьютерной модели используется второй вариант.

С

шк

(^""сЬавде

Рис. 1. Карта поведения

Рис. 2. Блок нелинейной функции, реализующий предикат рг1

Рис. 3. Фрагмент компьютерной модели, реализующий предикат рг1

Согласно математической модели начальные условия задаются математическим соотношением

У(0) = g, g = (0, v0,0,v0,...,0,v0)T .

Соответствующие начальные условия с применением новых языковых конструкций LISMA задаются в виде

array y [400]; i = [1 - 200]; y[2 * i - 1] = 0; y[2 * i - 1] = v0;

Как следует из (2), (3), при алгоритмическом задании правых частей системы ДУ необходимо учитывать, что при 1 < j < N неизбежно появятся переменные с индексами, выходящими за границу массива [1.. .400].

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

Язык спецификации. Версия языка LISMA, представленная в работе [3], не позволяет качественно специфицировать модели описанного класса задач. Данная возможность реализована путем добавления в язык конструкций, позволяющих описывать массивы фазовых и (или) алгебраических переменных и соответствующих им правых частей. Основная идея заключается в том, чтобы добавить в описание конструкций языка (начальные условия, уравнения, макросы) определенный параметр, изменяемый в заданных пределах. Такой подход позволяет для каждой параметризированной языковой конструкции генерировать соответствующие серии (массивы).

Для получения возможности спецификации моделей описанного класса в грамматику языка добавлены следующие конструкции:

1) оператор объявления массива <МАССИВ> ^ 'array' <ИД> '[' <ЦБЗ> ']' ; 'array' - ключевое слово,

<ИД> - имя массива (идентификатор),

<ЦБЗ> - размерность массива (целое положительное число);

2) оператор объявления счетчика

<СЧЕТЧИК> ^ <ИД>'=''[' <ЦБЗ>'-' <ЦБЗ> [ { ',' <ЦБЗ> '-' <ЦБЗ> } ]']' ; <ИД> - имя счетчика (идентификатор),

<ЦБЗ> - минимальное и максимальное значения счетчика (целые положительные числа). По умолчанию шаг изменения значения счетчика равен 1;

3) идентификатор с параметром <ПИД> ^ <ИД> '[' <АВ> ']' <ИД> - идентификатор,

<АВ> - индекс (арифметическое выражение).

Заполнение массива может производиться двумя способами: алгоритмически и явно. Пусть для алгоритмического задания элементов массива y[dim] в модели используются n независимых счетчиков или n -интервальный счетчик (выражения 1, 2 соответственно):

1. ci = [min. - max.], i e 1, n;

2. c = [min1 - max1, min2 - max2,..., minn - maxn].

В обоих случаях счетчики определяют интервалы алгоритмического задания элементов системы r = [min1, max1 ], ^ =[min2, max2],..., rn =[min„, max„]. При задании интервалов значений счетчика должны выполняться следующие ограничения:

mini < maxi, i e 1, n, min i < min i+1, i e 1, n -1, max i < maxi+1, i e 1, n -1, min1 > 1, max n < dim

Области алгоритмического A и явного B задания элементов массива описываются соответствующими выражениями

n n-1

A = U[minl, maxl],B = [l, minx) U U(maxl, minl+1) U (maxn, dim].

l=1 l=1

Обозначим через |r| количество целых чисел, принадлежащих интервалу r = [min; max ].

Тогда условие полноты системы при использовании массивов

принимает вид |A| + |B| = dim

(рис. 4).

Синтаксические формы записи для явного и алгоритмического задания элементов массива идентичны. Различие заключается в том, что при явном задании в качестве арифметического выражения указывается конкретное значение индекса.

С учетом изложенного построим компьютерную модель.

Компьютерная модель. Компьютерная модель системы (1), (2) приведена на рис. 5.

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

Область алгоритмического задания

Рис. 4. Условие полноты системы при использовании массива

Рис. 5. Компьютерная модель на языке LISMA

системы. Это, в свою очередь, приводит к неверному глобальному решению, и ограничивается использование неявных методов для решения рассматриваемого класса задач [4].

С другой стороны, решение жестких задач классическими явными схемами достаточно затруднено, так как обычно алгоритм управления шагом интегрирования строится на контроле точности численной схемы. Это естественно, потому что основным критерием является точность нахождения решения. Однако при применении алгоритмов интегрирования на основе явных формул для решения жестких задач данный подход приводит к потере эффективности и надежности, потому что на участке установления вследствие противоречивости требований точности и устойчивости шаг интегрирования раскачивается. В лучшем случае это приводит к большому количеству возвратов (повторных вычислений решения), а шаг выбирается значительно меньше максимально допустимого [1]. В связи с этим для решения данной задачи в системе ИСМА предлагается использовать алгоритмы на основе явных схем переменного порядка с алгоритмом управления шагом на основе не только критерия точности, но и устойчивости.

Рассмотрим пример контроля устойчивости на основе метода Фельберга 7-го порядка. Алгоритм управления шагом строится на основе соотношения [5]

h = max[h, min(hac, hst)], (4)

где h - текущий шаг интегрирования; hac - шаг по точности; hst - шаг по устойчивости, рассчитываемый по формуле

hst = (D / Vn)h, Vn = max |12K3 - 18K2 + 6K1)11 / | (K2 - Kl)l |,

1<i'< N

постоянная D = 5 ограничивает область устойчивости метода Фельберга 7-го порядка; Vn - оценка максимального собственного числа матрицы Якоби степенным методом; Kt (1 < l < 13) - стадии метода.

Применение соотношения (4) для управления шагом позволяет ограничить чрезмерный рост шага интегрирования по точности на участках установления жесткой задачи. Для метода Фельберга 7-го порядка с контролем устойчивости, интегрированного в систему ИСМА, также построена дополнительная численная схема 1-го порядка на основе первых семи стадий метода Фельберга 7-го порядка [5]. Интервал устойчивости схемы первого порядка D = 90 , что делает его в несколько раз эффективнее метода 7-го порядка на участках установления, где производные решений малы.

Локализация точек переключения гибридной системы. Проблема точки переключения событий является ключевой задачей в ГС [6-8]. Поиск точки переключения необходимо проводить при машинном анализе ГС вследствие того что выбор шага интегрирования по критериям точности и устойчивости может оказаться несогласованным с дискретными моментами переключения, т. е. обнаружение момента смены локальных поведений ГС не гарантировано. Наиболее опасными в этом смысле являются участки установления, когда производные решения невелики и как следствие шаг выбирается достаточно большим. Для решения данной проблемы в систему ИСМА встроены специальные алгоритмы локализации точки переключения.

Рассмотрим задачу Коши

У = f (t,y),y(tо) = Уо , (5)

где f : R x RN ^ RN - известная непрерывная вектор-функция; RN - N -мерное линейное вещественное пространство; Уо - решение в начальной точке to .

Решение задачи (5) будем искать явными методами, которые записываются в векторной форме

Уп+1 = Уп + hn+1Pn, n = 0,1,2,K

с заданными начальными условиями y0 . Здесь hn+1 - шаг интегрирования, (рп - заданная гладкая N -мерная вектор-функция, зависящая от правой части задачи (5).

Рассмотрим режим односторонней ГС [4] в виде задачи Коши с ограничениями

У = f(y), У(to) = Уо, g(У,t)<0,

где £(у, ^) - событийная функция или нелинейный предохранитель [6, 9]. Поскольку многие модели £(у, ^) , представляющие интерес, являются линейными, будем рассматривать их как наиболее важный класс событийных функций.

В работе [6] показано, что с учетом точности, устойчивости и динамики событийной функции £ (у, ^) формула управления шагом принимает вид

где

hp+l = min (hg ,max

hg =(y-1)

hn,min (hac, hst)

g n

dgn m , dg

(6)

(7)

Pn +■

ду д

у е [0, 1) обеспечивает поведение событийной динамики как устойчивой линейной системы, которая приближается к поверхности переключения £ (у, ^) = 0.

Сформулируем алгоритм интегрирования с учетом точности, устойчивости и динамики событийной функции £ (у, ^) . Пусть решение уп в точке ¡п вычислено с шагом Ип . Тогда управление шагом можно выполнить по следующему алгоритму:

Шаг 1. Вычисляется /п = / (уп, ¡п ) .

Шаг 2. Вычисляются ^ = £ (уп, ¡п ) , д8п/ду = д£ (уп, ¡п , д8п/д = д£ (уп, ¡п . Шаг 3. По формуле (7) вычисляется шаг , причем (рп = /п . Шаг 4. По формуле (6) вычисляется новый шаг Ъп+Х. Шаг 5. Выполняется следующий шаг интегрирования.

Численное решение задачи. Для численного эксперимента будем использовать два метода: неявный метод RADAU5 и явный метод Фельберга переменного порядка с контролем устойчивости. Расчеты проводились на ЭВМ Intel Pentium D 3.2 ГГц, 2 Гб DDR в среде ИСМА, требуемая точность решения s = 1e-6, начальный шаг интегрирования h0 = 1e 9. Прежде чем проводить сравнительный анализ численных методов, приведем результаты решения некоторых компонент системы ДУ (рис. 6).

В табл. 1 приведены результаты, полученные двумя методами. Анализ результатов показывает, что, во-первых, явный метод Фельберга выигрывает по показателю времени, что критически важно при решении задач большой размерности. Во-вторых, количество точек, рассчитанных явным методом, значительно больше, чем у неявного метода, что не всегда удобно при анализе результатов решения. Поэтому в ИСМА проводится предварительная обработка полученного решения. При этом время, затраченное на оптимизацию, существенно меньше времени, затраченного на получение результатов. Поэтому обработка в GRIN выполняется без потери эффективности.

Рассмотрим решения, полученные неявным методом RADAU5.

Таблица 1

Результаты расчета модели с использованием методов RADAU5 RKF78STP

RADAU5 RKF78STP

Количество точек 426 15 750

Расчет правой части 86 446 194 798

Время (CPU), с 50,9 31,2

Рис. 6. Результаты решения некоторых компонент системы ДУ: а - компонента >79; б - компонента >199; в - компонента >172; г - гибридная компонента phi

На рис. 7, а видно, что гибридная компонента меняет значение с запаздыванием, переключение происходит в момент времени t = 5,0157. Погрешность при переключении состояния, равная 0,015, приводит к существенной глобальной погрешности (рис. 7, б). Данный пример наглядно демонстрирует неприспособленность неявных методов к решению гибридных систем. При использовании выражения (7) для интегрирования ГС неявным методом получим новое значение времени переключения t = 5,00000013; теперь решение будет соответствовать результатам, приведенным на рис. 6. Однако включение алгоритма точной лока-

б

Рис. 7. Результаты решений, полученные неявным методом RADAU5: а - гибридная компонента phi; б - компонента решения y199 (пунктирная линия - верное решение)

лизации негативно влияет на производительность неявного метода, с другой стороны, применение подобного алгоритма к явному методу не приводит к серьезным временным потерям. Приведем окончательные результаты с учетом алгоритма локализации и оптимизации, а также результаты расчета данной задачи классическим методом Фельберга 7-го порядка (табл. 2).

б

а

а

Из результатов, приведенных в Таблица 2

табл. 2, следует, что при решении гибридных задач большой размерности явные методы с контролем устойчивости и переменным порядком более эффективны, чем неявные методы. При использовании явного метода Фельберга переменного порядка с контролем устойчивости для решения данной задачи требуется в два раза меньше времени, чем при использовании неявного метода RADAU5, и почти в три раза меньше, чем при использовании классического метода Фельберга 7-го порядка.

В системе ИСМА наряду с методом Фельберга 7-го порядка реализованы другие явные методы с контролем устойчивости и переменным порядком - Фельберга 5-го порядка, Дормана - Принса 8-го порядка, Мерсона и др. Однако для решения задач повышенной жесткости применение явных методов, в том числе методов с контролем устойчивости, может быть неприемлемым, вследствие того что для получения решения на переходном участке шаг интегрирования очень мал. Это может привести к численной неустойчивости, поскольку накопленная ошибка округления может превысить допустимую погрешность. Кроме того, такой процесс может потребовать гораздо больше временных затрат, чем применение неявного метода к данному участку задачи. С другой стороны, производительность (по количеству расчетов правой части) явного метода с контролем устойчивости на переходном участке в случае средней и малой жесткости может быть сравнима с производительностью неявного метода, так как шаг интегрирования выбирается достаточно большим.

Для анализа режимов ГС повышенной жесткости предлагается применять адаптивный алгоритм с автоматическим контролем жесткости [10].

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

Результаты компьютерного анализа показали, что комбинирование методов низкого и высокого порядков с помощью неравенства для контроля устойчивости позволяет существенно повысить эффективность расчетов. Это делает возможным применение явных методов переменного порядка и шага с контролем устойчивости для решения жестких задач большой размерности. Также показано, что применение неявных методов для интегрирования гибридных систем без учета динамики событийной функции g(y,t) приводит к неверному глобальному решению.

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

Список литературы

1. Новиков Е. А. Явные методы для жестких систем. Новосибирск: Наука. Сиб. предпр. РАН, 1997.

2. Новиков А. Е., Новиков Е. А., Шорников Ю. В. Аппроксимация матрицы Якоби в (2,2)-методе решения жестких сис-

тем // Докл. АН ВШ РФ. 2008. № 1 (10). С.30-44.

3. Свидетельство об офиц. регистрации программы для ЭВМ № 2007611024. Программа языкового процессора с языка

LISMA (Language of ISMA) / Ю. В. Шорников, И. Н. Томилов // Программы для ЭВМ. Базы данных. Топологии интегральных микросхем: Офиц. бюл. Рос. агентства по патентам и товар. знакам. 2007.

4. Lee E. A., Zhenq H. Operational semantics of hybrid systems // Proc. of Hybrid systems: computational and control (HSCC),

Zürich (Swenzeland), March 9-11, 2005. LNCS 3414.

5. Новиков Е. А., Шорников Ю. В. Контроль устойчивости метода Фельберга седьмого порядка точности // Вычисл.

технологии. 2006. Т. 11, № 4. С. 65-72.

6. Esposito J., Kumar V., Pappas G. J. Accurate event detection for simulating hybrid systems / Hybrid systems: computation

and control (HSCC). V. LNCS 2034. Springer-Verlag, 1998.

7. Бахвалов Н. С. Численные методы / Н. С. Бахвалов, Н. П. Жидков, Г. М. Кобельков. М.: Наука, 2000.

Результаты расчета модели с использованием различных методов

RADAU5 RKF78STP RKF78

Количество точек 561 15 779 35 566

Количество точек 561 315 351

с оптимизацией

Расчет правой части 90 423 195 175 501 468

Время (CPU), с 66,8 31,6 84,15

8. Сениченков Ю. Б. Численное моделирование гибридных систем. СПб.: Изд-во С.-Петерб. политехн. ун-та, 2004.

9. Esposito J. M., Kumar V. A state event detection algorithm for numerically simulating hybrid systems with model singulari-

ties // ACM Transactions on Modeling and Computer Simulation. 2007. V. 17, iss. 1. Р. 1-22.

10. Свидетельство об офиц. регистрации программы для ЭВМ №2007611459. Адаптивный алгоритм численного анализа жестких систем / Ю. В. Шорников, Е. А. Новиков, М. С. Денисов // Программы для ЭВМ. Базы данных. Топологии интегральных микросхем: Офиц. бюл. Рос. агентства по патентам и товар. знакам. 2007.

11. Шорников Ю. В. Моделирование гибридных систем явными методами // Системы управления и информационные технологии. 2007. № 4.2(30). С. 307-311.

Юрий Владимирович Шорников - канд. техн. наук, доц. кафедры автоматизированных систем

управления Новосибирского гос. техн. университета, (383) 346-11-00, e-mail: shornikov@asu.cs.nstu.ru, shornikov@inbox.ru; Иван Николаевич Томилов - магистр техники и технологий, аспирант кафедры автоматизированных систем управления Новосибирского гос. техн. университета,

e-mail: studnstu@ngs.ru;

Михаил Сергеевич Денисов - магистр техники и технологий, аспирант кафедры автоматизированных систем управления Новосибирского гос. техн. университет, e-mail: miha.x56@gmail.com