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

Термическая обработка является одним из главных этапов в технологии приготовления продукта, его сохранении и консервации. При этом разработка алгоритмов функционирования управляющей части устройств, связывающих электро- и теплофизические процессы, имеет давнюю историю. Если объект включается в контур программного управления и его динамические свойства учитываются при синтезе системы автоматического регулирования (САР), то необходимо подобрать модель процесса по степени сложности, приемлемой для такой цели. Наиболее детальные и полные модели объекта с описанием температурных полей нагрева его элементов имеют в своей основе расчетные схемы с распределенными параметрами. Для решения конечно-разностных уравнений, описывающих тепловой процесс, привлекаются специализированные пакеты, идет анализ температурных полей [1–3]. Такие модели не могут быть взяты за основу при синтезе САР для реализации на микро- контроллере (МК) программных алгоритмов прямого цифрового управления. Их цель – визуализация тепловых процессов, разработка конструкторских решений и т.п.

Для задач адаптивного управления [4] более продуктивным является использование моделей с сосредоточенными параметрами, упрощенных приемами декомпозиции [5, 6]. Но и при декомпозиции надо контролировать сложность модели, программно-аппаратной части МК, обеспеченность контроллера математическими пакетами при использовании, например, принципов оптимальной фильтрации Кальмана [7]. В худшем случае, при отсутствии такого контроля, стоимость системы на базе МК может оказаться неоправданно высокой и соизмеримой со стоимостью самого устройства.

Программная реализация алгоритма управления для рассматриваемого устройства на­грева жидкой среды средней мощности предполагает ряд допущений при переходе к модели с сосредоточенными параметрами, декомпозиции процесса. Включение в программу работы МК алгоритмов идентификации ненаблюдаемых параметров позволяет повысить потребительские свойства устройства и остаться в приемлемой ценовой нише: используются обычные контроллеры общепромышленного назначения [8].

Процесс разработки устройства нагрева жидкости (аппаратной части и ПО управляющего МК) представлен на рисунке 1.

Статья детально касается блока 2 (рис. 1) как части общего алгоритма разработки устройства нагрева жидкости. А именно, рассматривается задача идентификации параметров установки и жидкости с целью контроля качества работы САР, соблюдения технологии приготовления продукта. Материал строится так: предлагается математическая модель САР устройства, и на ее базе формулируется алгоритм идентификации параметров процесса, оцениваются эффективность и управляемость алгоритма, даются рекомендации по отладке ПО блока идентификации.

Математическое описание процесса нагрева жидкости и модель устройства

Рассмотрим алгоритм управления САР и идентификации параметров процесса для устройства нагрева жидкости, которая может быть как конечным продуктом, так и средой для опосредованного нагрева (процесс пастеризации) основного продукта, его консерва- ции. В основе САР лежит модель устройства нагрева жидкой среды, состоящая из индуктора, первичного источника переменного магнитного поля, вторичного тела, нагревающегося токами Фуко и являющегося сосудом для жидкости, и самой жидкой среды. При составлении модели процесса учитывается теплообмен с внешней средой. Взаимодействующие при нагреве элементы устройства изображены на рисунке 2. На блок-схеме указаны номера уравнений в следующей системе, описывающей тепловой и электромагнитный баланс:

(1)

Введены обозначения: изображения напряжения источника u1(p), тока индуктора i1(p), тока вторичного тела (тока Фуко) i2(p). Используются превышения температур τi(p) над температурой окружающей среды t4 (ti(p) – абсо- лютные значения температур), причем τi(p) = = ti(p) – t4, где индексы соответствуют индуктору (i = 1), вторичному телу (i = 2), жидкой среде (i = 3), окружающей среде (i = 4) [3, 4]. Параметры электромагнитных цепей индуктора и вторичного тела: активные (R1, R2) и индуктивные (L1, L2) сопротивления, коэффициент взаимной индукции (M12). Теплоемкости отдельных элементов схемы обозначены через Сj при j = 1, 2, 3, а коэффициенты теплопередачи – через Аjk при j=1, 2, 3, k = 2, 3, 4 и j ≠ k. Модель контура токов Фуко принята однослойная, с сосредоточенными параметрами. Учтено изменение активных сопротивлений индуктора (R1) и вторичного тела (R2) при нагреве (температурный коэффициент αi при i = 1, 2, возможно, для разных материалов изготовления индуктора и вторичного тела, базовая температура 20 °С).

Численные значения параметров и коэффициентов рассчитаны для нагревателя средней мощности (10–15 кВт): А12 = А23 = 200 Вт/°С; А14 = А24 = А34 = 20 Вт/°С; R1 = 0,3 Ом; L1 = = M12 = 40∙10-5 Гн; R2 = 1 Ом; L2 = 40∙10-4 Гн; С1 = 500 Вт∙с/°С; С2 = 1 000 Вт∙с/°С; С1 = = 21 000 Вт∙с/°С; α = 0,005 Ом/°С. Структура модели статической САР нагревателя приведена на рисунке 3. Среда моделирования – MATLAB [9, 10]. Показана подсистема модели нагревателя (блок 3), включенная в САР с коэффициентом обратной связи (блок 4), блоком 1, обеспечивающим нужный по технологии временной закон изменения температуры жидкости τ3ref(t). Подключаемая подсистема цепей коррекции (блок 2) обеспечивает для САР нужное качество в динамике [4]. Из критически значимых элементов для исследования отметим подсистему оценки качества идентификации (блок 9) и подсистему задания периода квантования (блок 8), вызванного конечным временем работы ПО МК. Дополнительные элементы модели выполняют сервисные функции задания коэффициентов контуров невязки (блок 7), формирования нестабильности значения теплоемкости C3ref.

На цифровой модели исследуются процессы идентификации двух ненаблюдаемых параметров: теплоемкости жидкости (С3) и активной мощности, выделяющейся во вторичном теле от протекания токов Фуко (i22R2(τ2)), где сопротивление вторичного тела есть функция превышения его температуры: R2(τ2) по (1). Ненаблюдаемыми эти параметры называем потому, что трудно обеспечить систему датчиками, напрямую измеряющими указанные величины. Тем не менее, они важны для обеспечения контроля над технологией изготовления продукта и техническим состоянием нагревателя. Предлагаемые в статье алгоритмы могут быть распространены на идентификацию и другие интересующие разработчика параметры. Качество идентификации будет оцениваться при скачкообразном изменении параметра на величину, равную 30–50 % от его первоначального значения.

Алгоритм идентификации параметров процесса нагрева жидкости

Построим идентификатор, используя идеи Эйкхоффа [4]. Метод использует математическое описание САР, прогноз значения идентифицируемого параметра (или параметров) с привлечением сигнала невязки (погрешности идентификации ε) и интегральных звеньев в каналах идентификации. Применительно к нагревателю жидкости, как отмечалось, идентифицируем два параметра: активные потери в контуре вихревых токов Фуко (i22R2(τ2)) и теплоемкость жидкости (С3). Остальные параметры устройства (теплоемкости, коэффициенты теплопередачи) определяются конструкцией нагревателя и считаются известными. Также даны параметры и переменные первичной электрической цепи, коэффициент взаимной индукции и комплекс температур. Необходимые для расчетов производные вычисляются программно-аппаратно. Решая в системе урав- нения 5 и 4 из (1), получим аналитическую связь параметров, предложенных к идентификации в виде:

,                 (2)

где F1, F2, F3 – полиномы, не содержащие идентифицируемые параметры (F1 = 1, F2 = dτ3/dt и т.д.).

Процесс идентификации по Эйкхоффу состоит из следующих этапов.

1.    Составляем модель исходного уравнения:

                (3)

2.    Умножаем левую и правую части (3) на ε:

      (4)

3.    Принимаем составляющую F1·ε (4) как производную от оценки первого идентифицируемого параметра, то есть

F1·ε = d(i22·R2(τ2))id/dt.

4.    Принимаем составляющую F2·ε (4) как производную от оценки второго идентифицируемого параметра, то есть

F2·ε = dC3id/dt.

5.    Получаем оценку идентифицируемых параметров на выходе соответствующих интегрирующих звеньев:

6.    Уравнение идентификатора аналогично исходному (3) и имеет вид:

           (5)

Блок-схема идентификатора, составленная по уравнению (5), приведена на рисунке 4а.

Как показывает эксперимент, идентификатор по Эйкхоффу работает неудовлетворительно по нескольким причинам. Во-первых, два контура взаимосвязаны и возмущения в одном тут же вызывает возрастание погрешности работы другого контура. Величиной коэффициента обратной связи по невязке (–Kε) этот недостаток устранить нельзя. Во-вторых, первая особенность усугубляется тем, что идентифицируемые параметры различны по природе и абсолютным величинам значений: теплоемкость составляет величину, равную десяткам тысяч, температура жидкости увеличивается с очень малым градиентом, а джоулевые потери во вторичном теле измеряются сотнями ватт.

Была предложена модификация (упрощение) методики Эйкхоффа. Блок-схема, поясняющая алгоритм, представлена на рисунке 4б. Идентификация каждого параметра проводится по своему уравнению, появляются две невязки для каждого параметра (εW и εC) и два независимых коэффициента обратной связи по невязке (KεW и KεC). Последнее дает больше свободы при настройке процесса. Ценой нововведения является приемлемое увеличение числа математических операций ПО (примерно на 30 %). Этому соответствует программная реализация вычисления второго слагаемого правой части четвертого уравнения следующей системы:

        (6)

Входными величинами к модифицированному алгоритму Эйкхоффа по (6) являются F11, F31, F12, F32.

Таким образом, для идентификации двух параметров используются два уравнения системы (1). Отметим, что прямая цифровая идентификация (прямое решение уравнений 4 и 5 системы (1) относительно неизвестных параметров) по результатам двух шагов цифрового управления, необходимых для вычисления производных, проигрывает предложенному алгоритму. Наличие в алгоритме интегрирующих звеньев (рис. 4б) вносит стабильность и устойчивость в процесс. При идентификации активной мощности интегратор проявляет фильтрующие свойства для высокочастотной помехи, источником которой является пульсация активной мощности с двойной частотой питающей сети. Наличие интегрирующего звена в блоке алгоритма идентификации теплоемкости позволяет исключить из ПО операцию деления (первое слагаемое правой части уравнения 5 системы (1) предполагает выделение теплоемкости C3 из произведения), а значит, не надо программно обслуживать возможную ситуацию «деления на ноль» при прямой цифровой идентификации.

Благодаря разделению алгоритма идентификации на два независимых контура появляется возможность настройки контуров оптимальным образом путем подбора коэффициентов обратных связей KεW и KεC.

Так как устройством управления в САР является контроллер, третьим параметром настройки служит период квантования (обеспечен блоком 8, рис. 2) – время выполнения программы (интервал времени между обновлениями входных величин к модифицированному алгоритму Эйкхоффа (F11, F31, F12, F32)). Причем этот период может быть различным для идентификации величин активной мощности, выделяемой во вторичном теле, и теплоемкости жидкости в нагревателе.

Оптимальная настройка алгоритма идентификации параметров нагревателя

Демонстрацию работоспособности алгоритмов идентификации проведем для двух параметров нагревателя: активная мощность, выделяющаяся в виде тепла во вторичном теле нагревателя (i22R2(τ2))id = W2id, и теплоемкость жидкости C3id = c3∙m3id, где c3 – удельная теплоемкость жидкой среды, а m3id – масса жидкости. Эти параметры могут показать эффективность работы нагревательного элемента (индуктора), необходимость восстановления его работоспособности и дать возможность вычисления массы жидкой среды, что определит соответствие процесса технологическому заданию на приготовление конечного продукта.

Рассматривается квазиустановившийся процесс подогрева жидкости со скачкообразным изменением идентифицируемых параметров на 30–50 % (W2ref ≈ (100 → 150) [Вт], C3ref = (14 000 → 21 000) [Вт∙сек./°С]). Возмущение в канале идентификации мощности (W2ref) дано приближенной величиной: в квазиустановившемся процессе работы нагревателя (рассматривается процесс с начальным временем в 75 сек. после начала нагрева) текущее значение активной мощности, выделившейся в электрической цепи вторичного тела, все время вычисляется блоком 3 (рис. 3) и может меняться. Время рассмотрения процесса идентификации Tid = 5 сек. Целевые функции, по значениям которых будет определяться качество идентификации, есть относительные значения среднеквадратических отклонений идентифицируемых параметров (σ*W2, σ*C3) в псевдослучайном процессе отработки возмущения, базовые величины – максимальные среднеквадратические отклонения (σW2max, σC3max) в эксперименте:

                (7)

Для формирования рекомендаций по настройке САР поставим два эксперимента (по количеству независимых контуров идентификации (рис. 4б)). В качестве первого фактора (x1) принимаем периоды квантования для каждого контура (T0W2, T0C3), а в качестве второго (x2) – коэффициенты усиления сигналов невязки (KεW, KεС, рис. 4б).

Получим рекомендации по настройке контуров идентификации в виде регрессионных зависимостей из теории планирования эксперимента. Значения кодированных факторов (Xi при i = 1–2) для двух экспериментов связаны с факторами (xi при i = 1–2) формулой нормирования Xi = (xi – xi mid) ⁄ (xi mid – xi min), где xi mid – среднее значение фактора; xi min – минимальное значение фактора. Выбираем ортогональный центральный композиционный план (ОЦКП) проведения эксперимента и получения регрессионных зависимостей. Общее число опытов одного эксперимента для двух факторов (n = 2) будет N = 2n + 2∙n + 1 = 9. Добавляются эксперименты в центре плана. Для двух экспериментов число опытов равно 18. В таблице 1 представлены уровни факторов и их кодированные значения.

Таблица 1

Кодировка факторов для двух ОЦКП

Table 1

Factor encoding table for two OCCP (Orthogonal Central Composite Planning)

Интервалы изменения факторов показательны и получены по результатам предварительного моделирования. Результаты экспери- ментов обрабатывались авторской програм- мой, созданной в среде Excel [11].

В базе есть возможность проводить регрессионный анализ результатов: контролировать адекватность модели по F-критерию (распределение Фишера), исключать незначимые коэффициенты с использованием критерия Стьюдента, оценивать точность регрессии.

Для каждого эксперимента в таблице 2 приводятся по три показательных опыта с указанием принятых уровней факторов [X1, X2] и полученных значений целевых функций (σ*W2, σ*C3).

Сопоставляя результаты отдельных экспериментов (табл. 2), отметим, что при идентификации W2 рост периода квантования X1 ведет к зашумлению результата из-за неточного баланса в уравнении 4 системы (1): нарушается баланс высокочастотной составляющей (двойная частота сети) левой части и суммы изменяющихся величин правой части уравнения. В результате на W2id накладывается высокочастотная помеха, усиливающаяся еще и ростом второго фактора – коэффициента усиления канала KεW невязки εW. С другой стороны, благоприятное и малое значение периода квантования T0W2 (эксперимент в точке [–1, –1]) дает возможность проявить фильтрующие способности контуру идентификации W2 за счет наличия в прямом канале алгоритма интегрирующего звена (верхняя часть схемы рис. 4б).

Эксперименты по идентификации C3 показывают более терпимое отношение контура идентификации к росту периода квантования T0C3 и коэффициента усиления KεC. В контуре идентификации нет явных источников шума – сказывается инерционность тепловых процессов. Однако и здесь излишнее увеличение величин T0C3 и KεC начинает давать негативный результат (эксперимент в точке [+1, +1] против эксперимента [0, 0]).

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

Поверхность отклика σ*W2 (рис. 5) дает неочевидный полезный результат: лучшая настройка достигается при минимальном периоде квантования T0W2 (X1) и коэффициенте KεW (X2), близком к среднему значению. При минимальном значении коэффициента KεW (X2 = –1) есть два возможных альтернативных приемлемых значения периода квантования T0W2 (X1 = ±1).

Для второго эксперимента и поверхности отклика σ*С3 оптимальная настройка очевидна. Она достигается при максимальном значении KεC (X2) и минимальном T0C3 (X1). Заметно то, что уменьшение периода квантования T0C3 улучшает процесс идентификации более энергично, чем увеличение коэффициента KεC.

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

Заключение

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

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

Поставлен эксперимент, и получены регрес- сионные зависимости, связывающие периоды квантования (интервалы работы ПО), коэффициенты усиления контуров алгоритма идентификации с качеством идентификации ненаблюдаемых показателей процесса нагрева: активными потерями во вторичном теле устройства и теплоемкостью жидкой среды. Регрессионные зависимости оказываются востребованными при отладке ПО МК.

Литература

1.    Рапопорт Э.Я., Плешивцева Ю.Э. Оптимальное управление температурными режимами индукционного нагрева. М.: Наука, 2012. 309 с.

2.    Шарапова О.Ю. Создание численной модели индукционно нагревательной установки периодического действия в среде наукоемкого расчетного программного комплекса FLUX // Технич. науки: проблемы и перспективы: матер. Междунар. науч. конф. СПб, 2011. С. 191–195.

3.    Данилушкин A.И., Батищев A.M., Зиннатуллин Д.А. Синтез квазиоптимальных алгоритмов пространственного распределения мощности индукционного нагревателя жидкости // Матем. моделирование и краевые задачи: сб. науч. тр. III Всерос. науч. конф. Самара, 2006. С. 46–49.

4.    Ким Д.П. Теория автоматического управления. Т. 2. Многомерные, нелинейные оптимальные и адаптивные системы М.: Физматлит, 2016. 440 с.

5.    Antipin A.S., Frizen V.E., Sarapulov F.N., Kondratev M.A. Power and energy parameters of the inductor with an electromagnetic shield. Proc. 15th Int. Conf. ELMA, 2017, pp. 406–409.

6.    Dolgikh I., Korolev A., Zakharov V. Temperature processes of induction heating simulation. Scientific Enquiry in the Contemporary World: Theoretical Basic and Innovative Approach, 2014, vol. 5, pp. 68–74. DOI: 10.15350/L_26/6.

7.    Задорожная Н.М., Лунев А.А. Адаптивные системы управления теплопрочностными испытаниями // Современные наукоемкие технологии. 2016. № 2-2. С. 216–219.

8.    Мюллер Д. Микроконтроллеры от Microchip выходят на новый уровень производительности и гибкости // Электроника: НТБ. 2018. № 5. С. 128–132. DOI: 10.22184/1992-4178.2018.176.5.128.132.

9.    Новожилов М.А., Пионкевич В.А. MATLAB в электротехнике. Иркутск: Изд-во ИНИТУ, 2016. 246 с.

10. Дьяконов В.П. MATLAB. Полный самоучитель. М.: ДМК Пресс, 2012. 768 с.

11. Льготчиков В.В. Регрессионный анализ результатов экспериментов. Свид. о гос. регистр. БД № 2005620067. М.: Роспатент, 2005.

References

1. Rapoport E.Ya., Pleshivtseva Yu.E. Optimal Temperature Control of Induction Heating. Moscow, Nauka Publ., 2012, 309 p. (in Russ.).
2. Sharapova O.Yu. Creation of a numerical model of a batch induction heating installation in the environment of the high-tech computational software complex FLUX. Proc. Intern. Conf. Tech. Sci.: Problems and Prospects. St. Petersburg, 2011, pp. 191–195 (in Russ.).
3. Danilushkin A.I., Batishchev A.M., Zinnatullin D.A. Synthesis of quasi-optimal algorithms of spatial power distribution of induction liquid heater. Proc. Conf. Math. Modell. and Boundary Problems. Samara, 2006, pp. 46–49 (in Russ.).
4. Kim D.P. Theory of automatic control. Vol. 2: Multidimensional, Nonlinear Optimal and Adaptive Systems. Мoscow, Fizmatlit Publ., 2016, 440 p. (in Russ.).
5. Antipin A.S., Frizen V.E., Sarapulov F.N., Kondratev M.A. Power and energy parameters of the inductor with an electromagnetic shield. Proc. 15th Intern. Conf. ELMA. 2017, pp. 406–409.
6. Dolgikh I., Korolev A., Zakharov V. Temperature processes of induction heating simulation. Scientific Enquiry in the Contemporary World: Theoretical Basiсs and Innovative Approach. 2014, vol. 5, pp. 68–74. DOI: 10.15350/L_26/6.
7. Zadorozhnaya N.М., Lunev А.А. Adaptive control system of heat-strength tests. Modern High Technologies. 2016, no. 2, part 2, pp. 216–219 (in Russ.).
8. Mueller D. Microchip microcontrollers reach a new level of performance and flexibility Electronics: STB, 2018, no. 5, pp. 128–132 (in Russ.). DOI: 10.22184/1992-4178.2018.176.5.128.132.
9. Novozhilov М.А., Pionkevich V.А. MATLAB in Electrical Engineering. Irkutsk, 2016, 246 p. (in Russ.).
10. Dyakonov V.P. MATLAB. Full tutorial. Moscow, DMK Press, 2012, 768 p. (in Russ.).
11. Lgotchikov V.V. Regression Experimental Results Mining. Certif. of database registration no. 2005620067, Moscow, Rospatent, 11.03.2005 (in Russ.).