Программные продукты и системы

Хорошо известный подход к математическому моделированию нелинейной динамической системы типа «вход-выход» основан на представлении отклика системы на внешнее воздействие в виде интегростепенного ряда (полинома) Вольтерра. Построить интегральную модель в виде полинома Вольтерра – значит решить задачу идентификации ядер Вольтерра, так как основная проблема при аппроксимации непрерывного отображения входного сигнала в выходной заключается в идентификации многомерных переходных характеристик системы.

В настоящее время разработано довольно много способов определения динамических характеристик, а универсальность применения этого математического аппарата позволяет создавать ПО для проведения экспериментов на компьютере. В частности, пакет VoltaireXL (компании Applied Wave Research) показал свою эффективность при описании электронных схем конечными суммами ряда Вольтерра. Один из наиболее известных способов решения задачи идентификации ядер Вольтерра основан на задании импульсных входных сигналов, однако многие физические процессы не допускают импульсных входов. Кроме того, как отмечено в [1], построение интегральных моделей в виде полиномов Вольтерра затрудняется неприемлемо большим объемом вычислений. В связи с этим разработка конструктивных методик моделирования нелинейной динамики с помощью полиномов Вольтерра по-прежнему является актуальной прикладной задачей.

В Институте систем энергетики им. Л.А. Мелентьева (ИСЭМ) Сибирского отделения РАН с начала 1990-х годов развивается методика идентификации во временной области, основанная на использовании физически реализуемых тестовых входных сигналов в виде функций Хевисайда с отклоняющимися аргументами. В работе [2] описана программная реализация алгоритмов идентификации, использующая эталонную модель теплообменного аппарата и созданная с целью построения квадратичных полиномов Вольтерра в случае скалярных и векторных входных возмущений. Проверка адекватности работы программы проводилась на реальных данных, полученных в ходе натурных экспериментов на высокотемпературном контуре ИСЭМ СО РАН.

Недостаток подхода [2] заключается в том, что многомерные интегральные уравнения Вольтерра I рода, к которым сводится задача восстановления ядер, имеют решения в нужных классах функций при весьма обременительных условиях разрешимости. Отметим, что для прогнозирования реакции системы на то или иное внешнее возмущение нет необходимости в знании ядер.

Рассмотрим новые алгоритмы идентификации, основанные на использовании метода интегрирования произведения [3]. Основное достоинство этого подхода в том, что исходная задача сводится к нахождению не самих ядер, а соответствующих интегралов от них.

Новые алгоритмы идентификации динамических характеристик

Пусть входной сигнал x(t) есть вектор-функция времени, состоящая из p³2 компонент xi(t),  отклик y(t) – скалярная (что не уменьшает общности) функция времени, непрерывно зависящая от x(t), кроме того, система стационарна в том смысле, что ее динамические характеристики не меняются за исследуемый промежуток времени tÎ[0, T]. Тогда математическая модель системы типа «вход-выход» может быть представлена в виде полинома Вольтерра N-й степени:

                              (1)

, (2)

где y(0)=0, , tÎ[0, T]. Ядра Вольтерра  в (2) симметричны лишь по тем переменным, которые соответствуют совпадающим индексам.

Ограничимся наиболее важным в приложениях случаем N=2. Вместо (1), (2) имеем

  (3)

                                              (4)

,    (5)

,(6)

причем ядра Вольтерра , , в (4) симметричны по переменным s1, s2, tÎ[0, T].

Введем равномерную сетку узлов ti=ih, nh=T. При достаточно малом шаге сетки h можно аппроксимировать многомерные свертки в (4)–(6) согласно методу интегрирования произведения [3] и перейти от (3)–(6) к сеточному аналогу квадратичного полинома Вольтерра:

,                                (7)

,                                              (8)

,                    (9)

,                   (10)

.

Задача идентификации для (7) состоит в определении массивов (8)–(10). Отметим, что векторность входного сигнала дает некоторую степень свободы, позволяющую варьировать формы тестовых сигналов для отдельных компонент вектора x(t). Выберем тесты, отличающиеся от введенных ранее в [2].

Для решения задачи идентификации gμ, gμμ (при фиксированном значении μ) зададим наборы сигналов:

                        (11)

                                (12)

где h – шаг дискретизации отрезка [0, T]; e(t) – функция Хевисайда; , . Вещественные числа , , характеризуют высоту возмущающих воздействий по входу xμ.

Отметим, что тестовые возмущения вида (11), (12) реализуют процедуру разделения отклика моделируемой системы на составляющие , r=1, 2, обусловленные индивидуальным влиянием m-й компоненты вектора x(t). Подстановка сеточного аналога (11), (12) в (7) приводит к системе линейных алгебраических уравнений:

   (13)

, 1≤j≤i≤n,                    (14)

относительно неизвестных gμ, gμμ. В (13), (14) , . Число уравнений в (13), (14) равно , что совпадает с количеством неизвестных (в силу симметричности Kμμ). Прямой подсчет дает, что

,

, ,

,

.

Перейдем к поиску оставшихся  неизвестных в (7), отражающих чувствительность системы к одновременному изменению входных возмущений xm(t) и xn(t). Восстановление  диагональных элементов  обеспечивают тестовые сигналы вида

                            (15)

где , , , v¹m¹l.

Для идентификации  (при i≠k) из (10) введем две группы тестовых сигналов с носителем ширины h, имеющих следующее аналитическое представление:

            (16)

           (17)

где c1,v≠0, c2,m≠0 – высота входных воздействий по входам xn и xm с разной физической природой, , , , v≠m≠λ. Согласуя точки разрывов в тестовых сигналах с узлами равномерной сетки  переходим к сеточным аналогам (15)–(17) и подставляем их в (7)–(10). В результате приходим к следующим формулам:

,

,                                                                    (18)

, ,                                 (19)

, ,                    (20)

определяющим элементы массива gvμ. В (18)  есть отклик системы на входные возмущения вида (15). В (19), (20) отклики  и  получены при входных сигналах (16) и (17) соответственно.

Замечание. Модификация изложенного в данной статье подхода для скалярной функции x(t) была с успехом применена в работе [4] при моделировании нелинейной динамики теплообмена в производственном процессе.

Программно-вычислительный комплекс для моделирования нелинейной динамики

С целью построения и тестирования интегральных моделей нелинейной динамики теплообмена был создан программно-вычислительный комплекс, использующий эталонную модель. Эталоном послужило описание процесса теплообмена в элементе теплообменного аппарата (теплообменнике) с независимым подводом тепла, представленное в [5]. Рассматривалось изменение энтальпии ∆i(t) (кДж/кг) на выходе при произвольных законах возмущений расхода вещества ∆D(t) (кг/с), при полной тепловой нагрузке ∆Q(t) (кВт) и отклонении энтальпии ∆iвx(t) (кДж/кг) на входе теплообменника (t – время, D – приращение к соответствующему параметру начального стационарного режима).

В данном комплексе реализованы алгоритмы, разработанные ранее в [2], а также алгоритмы, базирующиеся на использовании метода интегрирования произведения [3] и задании новых типов тестовых возмущений (11)–(12), (15)–(17). Вычислительный комплекс создан в объектно-ориенти­рованной среде программирования Borland C++ Builder и основан на функционально-модульном принципе. Комплекс включает модули построения интегральных моделей в виде линейных и квадратичных полиномов Вольтерра для скалярных и векторных входных сигналов, вычисления откликов эталонной и интегральных моделей при произвольном возмущающем воздействии, а также модуль подсчета управляемого входного воздействия, обеспечивающего желаемый (заданный) отклик интегральной модели (такая задача возникает в связи с задачами автоматического регулирования техническими объектами).

Архитектура системы показана на рисунке 1.

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

В блоке настройки входных данных пользователь имеет возможность изменять в диалоговом режиме параметры эталонной модели, задавать шаг дискретизации, длину временного отрезка [0, T], высоту тестовых входных сигналов.

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

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

Взаимодействие пользователя и системы происходит через интуитивно понятный интерфейс, который отправляет запросы блокам системы. Основную часть пользовательского интерфейса составляют окна ввода исходных данных для эталонной и интегральных моделей (рис. 2) и основное окно для проведения вычислительного эксперимента (рис. 3).

На рисунке 3 приведены графики входных возмущений расхода воды DD и теплоподвода DQ, а также результаты моделирования энтальпии Di по эталонной модели (I) и интегральным моделям (II, III) на заданные входные сигналы. Процедуры вычисления откликов интегральных моделей основаны на использовании метода средних прямоугольников и метода интегрирования произведения. Расчет выходных значений эталонной модели проводится с помощью метода трапеций. Достоверность результатов моделирования проверялась на наборе сигналов, использованных для построения модели (7)–(10).

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

В заключение отметим, что автором разработан программно-вычислительный комплекс, использующий эталонную модель теплообмена. Внедрены новые алгоритмы для моделирования, основанные на использовании метода интегрирования произведения. Описанные в работе алгоритмы в силу их универсальности могут использоваться при математическом моделировании самых разнообразных нелинейных систем, до- пускающих задание тестовых сигналов вида (11)–(12), (15)–(17). Для адаптации ПО к задаче моделирования других объектов потребуется введение дополнительных процедур, вычисляющих отклики исследуемых динамических систем на наборы тестовых входных сигналов.

Литература

1.     Пупков К.А., Цибизова Т.Ю. Реализация фильтра Вольтерра второго порядка для идентификации нелинейных систем управления // Наука и образование. 2006. № 6. URL: http://technomag.edu.ru/doc/58741.html (дата обращения: 19.01.2012).

2.     Солодуша С.В. Построение интегральных моделей нелинейных динамических систем с помощью рядов Вольтерра: дисс… канд. физ.-мат. наук. Иркутск, 1996. 153 с.

3.     Linz P. Product Integration Method for Volterra Integral Equations of the First Kind // BIT. 1971. Vol. 11, pp. 314–421.

4.     Щербинин М.С. Оптимизация потребления энергоресурсов турбокомпрессором М-1 ЭП-300 с использованием программно-вычислительного комплекса // Науч.-технич. вестн. ОАО «НК «Роснефть». 2010. № 3. С. 36–39.

5.     Таиров Э.А. Нелинейное моделирование динамики теплообмена в канале с однофазным теплоносителем // Изв. АН СССР: Энергетика и транспорт. 1989. № 1. С. 150–156.