CC BY
46
11. Турчин, В.Ф. К математическому моделированию социальной интеграции. Моделирование социальных процессов I В.Ф. Турчин. - М.: Наука, 1970.
12. Моришма, М. Равновесие, устойчивость, рост I М. Моришма. - М.: Наука, 1972.
13. Магницкий, Н.А. Новые методы хаотической динамики I Н.А. Магницкий, С.В. Сидоров. - М.: Едиториал УРСС, 2004.
14. Волкова, В.Н. Устойчивость социально-экономических систем I В.Н. Волкова, А.А. Денисов II Системный анализ в экономике: сб. мат-лов Межвуз. конф. - Таганрог, 2000. - С. 4-12.
15. Huxley, J. Problems of Relative Growth, methurn I J. Huxley. - London,1932.
16. Александров, В.В. Структурный анализ диалога I В.В. Александров, А.В. Арсентьева. - Л.: ЛНИВЦ АН СССР, 1984.
УДК 004.10:531.1 Е.А. Новиков, Ю.В. Шорников
ОСОБЕННОСТИ КОМПЬЮТЕРНОГО МОДЕЛИРОВАНИЯ КИНЕМАТИКИ СЫПУЧИХ СРЕД В СИСТЕМЕ ИСМА*
Разработаны специализированные инструментальные средства машинного анализа (ИСМА), позволяющие сократить затраты на подготовку компьютерной модели, численный анализ и автоматическую обработки результатов решения. Особенностями ИСМА являются: 1) смешанное описание динамических систем, то есть одновременное применение дифференциальных уравнений и структурных схем; 2) учет разрывов первого рода в правой части системы; 3) применение явных методов переменного порядка, шага с переменным числом стадий для численного анализа жестких и нежестких задач. Инструментальная среда реализована с учетом простоты описания математических моделей. Класс исследуемых объектов ограничен системами с запаздыванием.
Большинство технологических процессов в сельском хозяйстве осуществляются с использованием бункерных устройств. Они позволяют выполнять важные технологические операции дозирования, хранения, транспортировки и переработки сыпучих материалов. Движение сыпучей среды является одной из интересных и сложных проблем в механике сплошных сред. Вопросам изучения динамики сыпучих материалов посвящено много работ (например, [2-3]). Нами использована имитационная модель описания движения массы сухих семян в вибрирующем лотке универсального высевающего устройства [1].
Истечение массы семян происходит через семяводы - отверстия в дне лотка. Наибольший интерес представляет показатель равномерности высева. В соответствии с законом Ньютона для /-го элемента сыпучего материала определена система обыкновенных дифференциальных уравнений первого порядка [1]
X, = и,
(1)
= Fx -Mu,,
У i = у,
У = Fy -му
(2)
Й): J =М,
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант №05-01-00579-а).
Здесь ^ и ^ - действующие на /-й элемент суммарные силы в направлении х и у; щ, V - компоненты скорости центра масс в направлении х и у соответственно; лщ и /лу{ - диссипативные (вязкие) члены; л > 0 - коэффициент вязкости; J - главный момент инерции элемента относительно оси, перпендикулярной рассматриваемой плоскости; М. - момент сил трения.
Программирование и реализация (1)-(2) имеет свои особенности и трудности. Здесь для понижения трудоемкости реализации предлагается использовать специализированные инструментальные средства ИСМА [4], которые, как показывает практика [5], позволяют сократить трудозатраты на подготовку компьютерной модели, численный анализ и автоматическую обработки результатов решения. Особенностями ИСМА являются следующие факторы:
Смешанное описание динамических систем, то есть одновременное применение дифференциальных уравнений и структурных схем, учет разрывов первого рода в правой части системы обыкновенных дифференциальных уравнений (ОДУ), применение явных методов переменного порядка, шага с переменным числом стадий для численного анализа жестких и нежестких систем ОДУ и т.д.
Инструментальная среда ИСМА реализована с учетом простоты описания математических моделей и ориентирована на предметного пользователя.
Класс исследуемых в ИСМА объектов ограничен системами ОДУ с запаздыванием
^ = Дг(0, г(г -©), *], * > 0, г(0) = гй, (3)
где г еЯ" - вектор состояния, 2(г) = у(1) при г е [-©,0); I- независимая переменная; \у(г) - т-мерная вектор-функция запаздывания, т<"; ©={г1,...,^}т- вектор чистых запаздываний;
/ = {/г,■■■,/п}т- нелинейная вектор-функция, удовлетворяющая условию Липшица; = {^,..,г0п}т
- вектор начальных условий.
Структурное описание. Структурное редактирование и дальнейший анализ в ИСМА осуществляется с помощью интерпретатора схем. Структуры компонуются из набора функциональных модулей или примитивов (интегратор, нелинейность, запаздывание и т.д.) на экране монитора средствами разработанного графического редактора. Структурная модель динамики сыпучей среды (1) представлена на рисунке 1.
Рис. 1. Структура модели (1) в ИСМА, к = л
Структурная схема для вертикальной составляющей силы (2), действующей на компонент сыпучей среды, имеет тот же вид, только в обратной связи вместо компоненты скорости центра масс щ появится
компонента V.. Макроблоки действующих на элемент сил Гх и Г учитывают возможность изменения правой части в (1)-(2) следующим образом. Если лоток приходит в движение, то элементы сыпучей среды за несколько шагов по времени принимают «плотную упаковку». В этом случае к силам, возникающим при сжатии элементов, необходимо добавить соответствующие компоненты сил трения, причем в зависимости от
СО (° < 0 или О — 0) меняются силы ^, Ру и момент М.. С учетом этого общая структура будет иметь вид, представленный на рисунке 2.
а б
Рис. 2. Гибридная структура модели сыпучей среды
На рисунке 2,а представлена развернутая модель, а в силу идентичности структур (Fx - x) и (Fy - y) на рисунке 2,б произведено агрегирование идентичных структур макроблоками.
Таким образом, особенностью предлагаемой модели динамики сыпучих сред является наличие дискретно-непрерывного поведения системы. Такие системы классифицируются как гибридные [5]. Их поведение характеризуется непрерывными локальными зависимостями от определенных параметров (в нашем случае о) и дискретных событий, когда соответствующий параметр меняет свое значение. В этом случае система приобретает новое локальное непрерывное поведение. Такие системы представляются разным формализмом. Например, в MVS [5] и HYTECH [6] гибридные системы формализованы конечными автоматами с графическим интерфейсом «карты состояний»; в Simulink [5] визуализация производится диаграммами в Stateflow; в Dymola [7] в качестве базового формализма используются временные сети Петри с визуализацией состояний и допустимых переходов. В ИСМА такие системы удалось описать введением бинарных компонент в правую часть (3) [8].
Временная зависимость бинарной составляющей предопределяет поведение системы в настоящем и будущем. Обратимся к нашему примеру. Зададим логический предикат р = true при со< 0 и р = false при со> 0. Пусть локальное непрерывное поведение некоторой гибридной системы описывается ОДУ с разными правыми частями
z = g(z, t,a), о > 0; z = q(z, t,о), о < 0.
(4)
Разное локальное описание таких систем обуславливает изменение поведения в каждом локальном состоянии, причем, следуя описанию, сами решения являются непрерывными. Переходы от одного состояния к другому осуществляются в дискретные моменты времени, обусловленные соответствующим значением предиката р(о). Предикат р(о) в данном случае рассматривается как бинарная компонента из множества {0,1} = {false,true}. Реализация р(о) в ИСМА имеет вид нелинейного блока с разрывной функцией. Способом задания начальных условий для анализа движения в таких системах является механизм припасо-вывания краевых значений.
Языковая модель. Переход от системы дифференциальных уравнений в форме Коши к компьютерному описанию происходит переписыванием правой части ОДУ либо во встроенном редакторе ИСМА, либо в любом другом текстовом процессоре, допускающим текстовый формат данных. Языковой процессор ИСМА интерпретирует введенные данные и придает им соответствующие типы, спецификации и т.д. В результате трудоемкость перехода к компьютерной модели существенно снижается. В известных системах моделирования при описании правой части системы ОДУ в грамматику языка вводят логические выражения. Так, для языка Modelika в системе моделирования Dymola [7] такое представление компьютерной модели для реализации (4) будет иметь вид der z = when о < 0 then q else g. Введение в язык явных логических конструкций типа when - then - else усложняют грамматику языка, и как следствие усложняется языковой процессор. Этого можно избежать, если воспользоваться следующими несложными преобразованиями [8]. Очевидно, объединением (4) через предикат р(о) будет уравнение
z = p(o)[q(z, t,о) -g(z, t, о)] + g(z, t, о) .
(5)
Можно также воспользоваться встроенными в ИСМА библиотечными функциями без универсальной нелинейной функции p(o). В этом случае достаточно использовать функцию sign, а уравнение (4) примет вид
z = sign (o)[q(z, t, о) - g(z, t, O)] + g(z, t, о), (6)
где sign (о) = 0 при o> 0 и sign(o) = 1 при o< 0. Встроенный в ИСМА текстовый редактор Redit с точностью до обозначений в соответствии с грамматикой G^CMA] допускает запись (6) практически без изменений. Программный фрагмент на специализированном языке системы ИСМА будет иметь вид: z = sign (о)* (q(z, t,o)-g(z, t,o)) + g(z, t,o) . Функции g (z, t,o) и g (z, t,o) задаются как макроопределения с арифметическими выражениями в правой части в соответствии с (1)-(2), а sign (о) является встроенной внутренней функцией библиотеки ИСМА.
Жесткость системы. При описании методов интегрирования задача (1)-(2) для сокращения записи ниже будет рассматриваться в виде
z = f z) z(to) = ^ to ^^ Ч, (7)
где z и f- гладкие вещественные л-мерные вектор-функции; t - независимая переменная. Особенностью рассматриваемой модели динамики сыпучих сред является жесткость системы ОДУ [1]. Поэтому для их численного анализа необходимо применять специальные численные методы решения жестких задач. Библиотека методов численного интегрирования ИСМА включает не только классические схемы типа Рунге-Кутта разного порядка точности
m i-1 i-1
zs+1 = zs +Z p>k, k = hf (ts + b'LPv , zs + Z jj ), (8)
i=1 j=1 j=1
но и оригинальные современные методы решения задач разной степени жесткости [9]. В частности, в состав библиотеки системы ИСМА включены известные методы Мерсона [11] и Фельберга [10], соответственно четвертого и пятого порядков точности. Для каждого метода на основе тех же самых стадий достроены численные формулы первого и второго порядков с расширенными областями устойчивости. С применением данных дополнительных численных схем построены алгоритмы интегрирования переменного порядка и шага с выбором наиболее эффективного метода на каждом шаге исходя из критерия устойчивости. В результате быстродействие методов Мерсона и Фельберга при решении жестких задач увеличилось соответственно в 14 и 20 раз. Следует отметить, что даже в случае нежестких задач за счет изменения порядка точности в процессе расчетов эффективность алгоритмов тоже возросла, хотя и не столь значительно.
В качестве примера конкретного алгоритма интегрирования изучим модифицированный метод Мерсона переменного порядка и шага, настроенный на эффективное решение жестких систем. Рассмотрим пятистадийный метод вида
6
zs+1 = zs +Z p<k , k1 = hf (tn , Уп X k 2 = hf (tn +^ К zn +T k1X
11
-h, zn +3 n 3
К = ¥& +1 h,zn +1 kl +1 k4 = ¥& +1 h,zn +1 kl + 3(9)
3 6 6 2 8 8
1 3
К5 = ¥(tn + h, Уп +-К -- К3 + 2К4)-
При значениях коэффициентов p1 = p5 = 1/6, p2 = Pз = 0, p4 = 2/3 численная формула (9) совпадает с методом Мерсона [11] и имеет четвертый порядок точности. Используя обозначения Є', =|| 2К - 9£3 + 8К4 - к II /150 и V = 6тах | (К3 - К ) /(к - К) I, Для контроля точности (9) мож-
1<г<п
но использовать неравенство ^ <є, а для контроля устойчивости - следующее: V <C [9]. Постоянную
C можно выбрать равной 3,5, то есть равной длине интервала устойчивости. Если использовать набор коэффициентов p1 = 5,248365568 -10-1 , p2 = 3,260928 -10-1, ^ = 1,395154944 -10-1 ,
p4 = 9,5158272 -10-3 и p5 = 3,93216 -10-5, то численная формула (9) имеет первый порядок точности,
i =1
а ее интервал устойчивости расширен до 50 по вещественной оси. Поэтому, если шаг выбирается из условия устойчивости, что имеет место на участке установления, то метод первого порядка будет примерно в 14 раз эффективнее метода Мерсона. Введем обозначения А" = 0,5(3 -6с2)|| £2 -||,
А" = 0,5(3 6°2) || Л+0 к1||. Тогда для контроля точности схемы первого порядка можно при-
менять неравенства А' <£ и А" <£ при с2 = 0,16 [9].
В силу того, что интервал устойчивости численной формулы первого порядка в 14 раз больше, чем у метода четвертого порядка, то переход с четвертого на первый порядок может сопровождаться увеличением шага в 14 раз. Заметим, что увеличение шага интегрирования в 14 раз осуществляется только за счет смены коэффициентов ,1 < / < 6, не приводит к увеличению вычислительных затрат. Однако, как показывают
расчеты, такой быстрый рост шага может приводить к невыполнению неравенства для контроля точности метода первого порядка. В результате осуществляется обратный переход. Чтобы избежать таких неоправданных переходов туда и обратно, в алгоритм интегрирования введена численная схема второго порядка точности. Переключение с первого порядка на четвертый и обратно осуществляется через метод второго порядка.
При коэффициентах р = 3,7789366573 2-10-1 , р2 =-9,30131004367-10-1,
р =-2,0390491435 8-10-2 , р = 1,5115746629 4, р = 6,1053167133 • 10-2 схема (9) имеет второй порядок точности и максимальный интервал устойчивости, равный приблизительно 8,6. Переход с одной схемы на другую осуществляется по следующему правилу. При расчетах по схеме четвертого порядка нарушение неравенства V < 3,5 вызывает переход на численную формулу второго порядка. При расчетах по
методу первого порядка выполнение неравенства V < 8,6 вызывает переключение на схему второго порядка. При расчетах по методу второго порядка выполнение неравенства V < 3,5 означает переход на четвертый порядок, а нарушение неравенства V < 8,6 вызывает переключение на первый порядок. Неравенства для контроля точности метода второго порядка имеют вид 0,17А" <е и 0,17А" <£, где величины А и А заданы при описании метода первого порядка точности.
Применение алгоритмов переменного порядка и шага приводит к существенному повышению эффективности алгоритмов интегрирования. Однако в случае применения таких алгоритмов к жестким системам точность вычислений получается всегда лучше требуемой. Это связано с тем, что на участке установления шаг интегрирования по точности может быть выбран достаточно большим, потому что на данном участке производные решения малы, но он ограничен неравенством для контроля устойчивости. На участках установления нет смысла использовать численные формулы высокого порядка точности. Быстродействие можно повысить за счет применения методов низкого порядка, но с большими областями устойчивости. Поэтому дальнейшее повышение эффективности достигается за счет построения алгоритмов интегрирования не только с переменным порядком и шагом, но и с переменным числом стадий. Известно, что максимальный порядок точности т -стадийного метода типа Рунге-Кутта первого порядка точности равен 2т2. В результате на каждое вычисление правой части приходится 2т единиц длины интервала устойчивости. Это означает, что если шаг ограничен по устойчивости, то с ростом т эффективность метода возрастает. В состав ИСМА включены методы с адаптивной областью устойчивости [9]. В результате эффективность повысилась более чем в 200 раз.
Выводы. Использование специализированных инструментальных средств позволяет предметному пользователю без дополнительных затрат достаточно просто перейти от математической модели сыпучей среды к компьютерной реализации через структурное и/или языковое описание математических абстракций. Адаптивный алгоритм выбора численной схемы интегрирования с переменным порядком, шагом и переменным числом стадий обеспечит необходимую точность и адекватность численных экспериментов с учетом жесткости. Средства автоматизации графической интерпретации результатов моделирования в ИСМА освобождают исследователя от рутинных задач документирования и хранения результатов решения.
Литература
1. Богульская, Н.А. Имитационная динамическая модель движения сыпучей среды / Н.А. Богульская, И.О. Богульский, А.А. Вишняков. - Красноярск: ИВМ СО РАН, 2005 (в печати).
2. Садовский, В.М. Задачи динамики сыпучих сред / В.М. Садовский // Мат. моделирование. - 2001. - Т. 13.
- № 5. - С. 62-74.
3. Бондаренко, Н.Ю. Параллельные алгоритмы решения задач динамики сыпучих сред / Н.Ю. Бондаренко, В.М. Садовский // Распределенные и кластерные вычисления: мат-лы II-й школы-семинара. - Красноярск: ИВМ СО РАН, 2002. - С. 166-176.
4. Инструментальные средства машинного анализа. Свидетельство официальной регистрации программы для ЭВМ № 2005610126 / Ю.В. Шорников [и др.]. - М.: Роспатент, 2005.
5. Бенькович, Е.А. Практическое моделирование динамических систем / Е.А. Бенькович, Ю.Б. Колесов, Ю.Б. Сениченков. - СПб.: БХВ, 2002. - 464 с.
6. http//www.eecs. Berkeley. edu/tah/HYTECH. - Henziger, T.A. A user Guide to HYTECH.
7. Elmgvist, H. Objеct-Oriented Modeling and Automatic Formula Manipulation in Dymola / Н. Elmgvist// SIMS'93.
- Scandinavian Simulation Society, 1993.
8. Шорников, Ю.В. Компьютерное моделирование билиарной системы специализированными средствами / Ю.В. Шорников // Науч. вестн. НГТУ. - 2004. - №3(18). - С. 31-42.
9. Новиков, Е.А. Явные методы для жестких систем / Е.А. Новиков. - Новосибирск: Наука, 1997.- 195 с.
10. Fhelberg, E. Klassische Runge-Kutta-Formeln funfter und siebenter Ordnung mit Schrittweitenkontrolle /
Е. Fhelberg // Computing. - 1969. - № 4. - P. 93-106.
11. Merson, R.H. An operational methods for integration processes / R.H. Merson // Proc. Symp. On Data Proc. Weapons Research Establishment. - Salisbury. - Australia. - 1957.
УДК 621.3: 519.876.5 (06) М.И. Соколов
МОДЕЛИРОВАНИЕ И ОПТИМИЗАЦИЯ ПРОЦЕССА ФУНКЦИОНИРОВАНИЯ ПРЕДПРИЯТИЙ ПРИ НЕПОЛНОЙ ИНФОРМАЦИИ
Предлагаются математические модели и методы оптимизации функционирования предприятий с позиции линейных динамических систем, включая условия неполной информации. Синтез моделей осуществляется последовательной конкретизацией за счет учета различных компонент и аспектов деятельности предприятия, что позволяет получить ее многопродуктовую модель, учитывающую длительность производственного цикла, временные запаздывания, взаимоотношения предприятия с финансовыми структурами, различные стратегии маркетинга. Решение задачи оптимизации процесса функционирования предприятия осуществляется с помощью модифицированного алгоритма стохастического квазиградиента.
Введение. Центральной проблемой современной теории моделирования и управления являются задачи исследования сложных систем в условиях различного уровня априорной информации о присущих им закономерностях. Подобная ситуация характерна также и для задач моделирования и оптимизации предприятий, что требует разработки адекватных математических и информационных средств.
Функционирование предприятия осуществляется в условиях неопределенности относительно внешних условий, что требует создания соответствующих математических моделей и алгоритмов управления. В подобных условиях применение традиционных аналитических и статистических методов приводит к упрощению закономерностей функционирования предприятия, а получаемые при этом результаты в основном имеют ограниченное практическое значение.
Перспективное направление выхода из создавшегося положения состоит в разработке математических средств исследования систем, сочетающих преимущества статистических и аналитических моделей, обеспечивающих наиболее полное использование априорной информации (экспериментальные данные и физические закономерности функционирования изучаемых объектов). В предлагаемой статье систематизи-
