Автоматика и телемеханика, № 1, 2014
© 2014 г. C.B. СОЛОДУША, канд. физ.-мат. наук (Институт систем энергетики им. Л.А. Мелентьева СО РАН, Иркутск)
К ЗАДАЧЕ МОДЕЛИРОВАНИЯ ДИНАМИКИ ТЕПЛООБМЕННИКОВ КВАДРАТИЧНЫМИ ПОЛИНОМАМИ ВОЛЬТЕРРА1
Применительно к проблеме автоматического управления нелинейной динамической системой проведено сопоставление различных методик идентификации квадратичных полиномов Вольтерра, основанных на введении специальных семейств кусочно-постоянных тестовых входных сигналов. Приводятся результаты вычислительных экспериментов для эталонной модели теплообмена.
1. Введение
Задача моделирования нелинейных динамических систем с помощью аппарата функциональных рядов Вольтерра имеет богатую историю. Применение конечного отрезка [1]
N * * т
(1) у(*) = Р^х(г)) / ... Кт(М1,...,5т)П * € [0,Т],
т=1 0 0 г=1
для моделирования нелинейной динамики различной природы базируется на теореме Фреше [2] и различных ее обобщениях. Проблема использования (1) для моделирования технических (в частности, электро- и теплотехнических) систем рассматривалась в монографиях [3-6], в которых приведена обширная библиография.
Работы Института систем энергетики им. Л.А. Мелентьева (ИСЭМ) СО РАН в области построения интегральных моделей на базе полиномов Вольтерра были начаты в 90-х гг. прошлого века. В публикациях [7, 8] предложена оригинальная методика идентификации ядер Вольтерра, основанная на задании многопараметрических семейств кусочно-постоянных тестовых входных сигналов. При этом задача идентификации сводится к решению линейных многомерных уравнений Вольтерра I рода с переменными верхними и нижними пределами. Данный подход был развит (см. библиографию в [9]) и применен для описания реальных процессов теплообмена на установке ВТК (высокотемпературный контур) ИСЭМ СО РАН [10], а также эталонной модели теплообмена [11-13].
Дальнейшие исследования коллектива связаны со следующим этапом математического моделирования. Если ядра Кт идентифицированы, то при
1 Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (проект № 12-01-00722-я).
заданном у(Ь) (1) является полиномиальным относительно ж(£) интегральным уравнением, которое возникает в задаче автоматического управления динамическим объектом. Исследованию специфики полиномиальных уравнений Вольтерра посвящен цикл публикаций А.С. Апарцина (см. библиографию в [14]). Численному решению подобных уравнений, связанных с проблемой автоматического управления динамикой теплообмена, посвящены работы [15, 16]. Разработанные алгоритмы были реализованы в программном продукте [17], на который получено авторское свидетельство [18].
Данная статья посвящена сравнению имеющихся подходов к моделированию нелинейной динамики в наиболее интересном для приложений квадратичном случае, когда N = 2 в (1). В качестве эталонной динамической системы рассмотрена математическая модель переходного процесса в элементе теплообменного аппарата (теплообменнике) с независимым подводом тепла, представленная в [19]:
А-До )( О0 \ / -Л11^К* -л2}
(2) АШ = ^^ J уА(2(г]) - ^ЧЛ)) Iе 17 -е "
ь е [0,Т].
В (2) Ь — время, А- и А2 — некоторые константы, индексами "0" обозначены параметры начального стационарного режима, Р = 0,16 (кг/с), Qo = = 100 (кВт), г0 = 434 (кДж/кг), А — приращение, например Р(Ь) = Р0 + + АР(Ь). Числовые характеристики, входящие в (2), принимались соответствующими реальной установке ВТК.
2. Построение интегральных моделей в виде квадратичных полиномов Вольтерра
Исследуется стационарная динамическая система, такая что ее динамические характеристики остаются неизменными за время Т переходного процесса. Рассмотрим случай векторного входа ж(Ь) = (АР(Ь), А2(Ь)) и скалярного отклика у(Ь) = Агег(Ь). Учитывая линейность (2) по А2(Ь), вместо (1) при N = 2 имеем
(3) Аг^(Ь) = у КТ1(в)АР(* - з)^ +
0
г г г
+ / / К^11(в1,в2)АР(Ь - в-)АР(Ь - «2)^1^2 + I ККо(в)Ад(Ь - в)^ +
0 0 0 г г
+ У J К^12(в1,в2)АР(Ь - в1)А2(Ь - в2)^Мв2, Ь е [0,Т].
00
В (3) (0) = 0, ¿1то^(Ь) е С(01)Т], ядро К11 симметрично по перемен-
ным 81,82. Помимо (3) рассмотрим квадратичную модель с комбинацией ли-
г
нейной нестационарной и билинейной стационарной составляющих [20]:
t
(4) Ai2mod (t) = j Ki(t,s)AD(s)ds +
t t
+ / / Kin(si,S2)AD(t - si)AD(t - s2)dsids2 + K2(t,s)AQ(s)ds +
0 0
t t
+ У ] К Фи $2^(1 - $г)Аа(Ь - $2)с$ф2, Ь € [0, Т].
0 0
Введем сетку узлов г = 1,?г, пЛ = Т. Всюду далее /г = 1 (с) (такой выбор шага связан с данными, полученными в ходе реального эксперимента на ВТК ИСЭМ СО РАН). Выбор устойчивого метода построения разностных аналогов (3), (4) основан на саморегуляризующем свойстве процедуры дискретизации [21]. Простейшими из кубатур, которые обладают данным эффектом, являются формулы прямоугольников — правых, левых и средних (формула средних прямоугольников предпочтительнее, так как имеет алгебраический порядок точности, равный двум). Методы типа Рунге-Кутта позволяют повысить порядок сходимости, однако при этом существенно увеличивается трудоемкость вычислений [22]. Следуя [22], в качестве "базовой" используем формулу средних прямоугольников.
Для построения квадратичных интегральных моделей использовались алгоритмы из [11, 20]. Согласно [11] идентификацию К\, К\\ обеспечивают тестовые сигналы вида
0
t
0
(5)
(6)
х
1,ш
хаа (t) = aDoI (t), X2 (t) = 0,
(t) = aDo (I (t) - I (t - u))
x2 (t) = 0
соответственно, где I(t) — функция Хевисайда, i = 1, 2, a1 + a2 = 0, 0 ^ u ^ ^ t ^ T. Согласно [20] восстановление Ki, K11 из (4) выполняют сигналы вида (6).
Во многих физических приложениях важно значение выходного сигнала динамического объекта в конце переходного процесса. В качестве критерия точности моделирования выбираем значение коэффициента "Mean Absolute Percentage Error" (MAPE) при ti = T:
(7)
MAPEj (t)
t=T
l_^\Aihet(T)-Ail В
jimo(i
(T )|
в=1
A ihet(T)
100%, j = 1, 2.
1 5 10 15 20 1 5 10 15 20 1 5 10 15 20
ю, с ю, с ю, с
□ если МЛРЕ^Т) > МАРЕ2(Т); если МЛРЕ](Т) < МЛРЕ2(Т); если МЛРЕ](Т) = МЛРЕ2(Т)
Рис. 1. Результаты расчетов для В = В1, В2, В3 с точностью 6 =10 4.
|а|, %
25
20
5 = 10-
¡ИР
15
10 15 20
ю, с
|а|, %
25
20
5 = 10-2
15
10 15 20
ю, с
|а|, %
25 20
5 = 10-1
И 'АО. 1 1 5 2
е>я>. ум 9. 2
УМУМХМ 2 2 5 2
УУМОМ9. 5
ъжтшя
умуохгл
15
10 15 20
ю, с
Рис. 2. Результаты расчетов для В = 25 с точностью 6 = 10 1 — 10
1-3
Проиллюстрируем результаты тестирования квадратичных моделей (3), (4) на допустимом множестве X(В,Т) входных сигналов
(8)
Арв(Ь) = вРо(1 (Ь) - 2/(Ь - ш) - I(Ь - Т)), А2(Ь) = 0,
где 0,0016 < вР < 0,0016В, 0 < ш < Ь < Т. Расчеты МАРЕ1(Т), МАРЕ2(Т) по (7) проводились для 1 ^ ш ^ 20, Т = 30, В = 35, 30, 25. Отклики Аг^ (¿¿), Аг^ ^ (¿¿) на (8) вычислялись по разностным аналогам (3), (4),
в которых ядра Вольтерра К1, К1, К11 настроены на тестовые сигналы с амплитудами 0,032 ^ |а|р ^ 0,0016В. На рис. 1 сравниваются МАРЕ1(Т), МАРЕ2(Т), полученные с точностью 6 = 10-4 при В = 35, 30, 25. Результа-
ты расчетов для В = 25 с точностью 6 = 10 — 10 иллюстрирует рис. 2 (обозначения на рис. 2 соответствуют принятым ранее на рис. 1). Области, представленные на рис. 1, 2, не являются универсальными, так как выявлены на основе анализа результатов вычислительных экспериментов, проведенных с использованием одной эталонной модели (2).
Отметим, что на точность моделирования существенно влияет выбор амплитуд тестовых сигналов, используемых для построения (3), (4). В частности, при согласовании амплитуд тестовых сигналов с величиной действующих возмущений отклики интегральных моделей совпадают с откликами эталонной модели, так что МАРЕ1(Т) = МАРЕ2(Т) = 0%.
Приведем значения невязок для наиболее "неприятных" входных сигналов из множества допустимых. Было получено, что для В = 25:
МЛРЕ1(Т) < 5,62%, МЛРЕ2(Т) < 5,18%, для В = 30: МЛРЕ1(Т) < 6,89%, МЛРЕ2(Т) < 6,38%, для В = 35: МЛРЕ1(Т) < 8,41%, МЛРЕ2(Т) < 7,93%.
Вычислительный эксперимент показал, что области предпочтительности использования той или иной интегральной модели зависят от длины отрезка Т, амплитуды а, используемой при построении интегральных моделей, а также точности расчетов 5.
До сих пор центральной в работе была проблема восстановления ядер Вольтерра в (3), (4). Допустим, что с помощью того или иного подхода эта проблема решена. В качестве следующего этапа математического моделирования можно рассмотреть задачу стабилизации (регулирования), связанную с поиском управляющего воздействия х(£), поддерживающего выходной сигнал у(£) на заданном уровне у*. Такая постановка возникает в связи с задачами автоматического управления техническими объектами. Выберем в качестве управляющего воздействия в (3), (4) сигнал Х\(Ь) = ДР(£). Для определенности входное возмущение х2(£) = ДQ(t), £ € [0, Т], считаем заданным. В этом случае уравнения (3), (4) являются полиномиальными интегральными уравнениями Вольтерра I рода, непрерывное решение которых, вообще говоря, носит локальный характер. В [16] рассмотрен вопрос получения управляющего воздействия Х1 (£) с учетом апостериорных данных об отклонении выходной переменной у (£) от желаемого значения у* = 0, так что х1 (£) = и (£ — К), и (£) = 0, £ € [-Ь, 0], К — известное постоянное запаздывание. В этом случае задача регулирования нелинейного динамического объекта сводится к поиску непрерывных решений и1(£), и*(£) полиномиальных уравнений
3. Полиномиальные уравнения Вольтерра в задаче автоматического управления
ь
(9)
о
ь ь
оо
ь ь
оо
ь
о
(10) У"к1(*,в)и2(в - Л)^ +
0
г г
+ J ! ^11(51, «2)«2(£ - «1 - Л)«2(^ - «2 - +
00 г г
+ / / К^12(«1,52)и2(^ - «1 - Л)Д2(£ - 82)^1^2 =
00
г
= /2(«) -/ ад^Дд^
соответственно. В (9), (10) /¿(¿) = - - Л), г = 1, 2. Сигнал ^(£) = у* -
- Уг(^)> где У1СО = Дг1то^ У2(£) = Дг2то^ £г(0 = 0, £ € [-Л, 0], считаем рассогласованием или ошибкой управления. При малом £ численное решение (9), (10) заведомо существует [14].
Пусть и = Иг, = (г — г = 1 ,п, п!г = Т. Обоснование сходимости численного решения полиномиального
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.