УПРАВЛЕНИЕ ТЕХНОЛОГИЧЕСКИМ ПРОЦЕССОМ РЕКТИФИКАЦИИ НА ОСНОВЕ ВИРТУАЛЬНЫХ АНАЛИЗАТОРОВ

Олег Снегирев

УПРАВЛЕНИЕ ТЕХНОЛОГИЧЕСКИМ ПРОЦЕССОМ РЕКТИФИКАЦИИ НА ОСНОВЕ ВИРТУАЛЬНЫХ АНАЛИЗАТОРОВ

Во введении представлена актуальность рассматриваемого вопроса, обозначены цели и задачи, рассматриваемые в работе. В обзорной части рассмотрен процесс ректификации, существующие системы управления процессом ректификации. В исследовательской части проведено исследование наблюдаемости технологического объекта ректификации. Разработан адаптивный виртуальный анализатор по методу “движущегося окна”, предложены критерии выбора оптимальной ширины окна адаптивного виртуального анализатора. Рассмотрен процесс определения передаточной функции в PDS HONEYWELL. Построена модель ПИД-регулятора и получены его коэффициенты различными методами. В заключении приводятся результаты проделанных работ.

Автоматика. Вычислительная техника

Диссертации

Вуз: Дальневосточный федеральный университет (ДВФУ)

ID: 5ea7c9d6ca23590001c9b54b

UUID: 7b3f3750-6b45-0138-32f2-0242ac180006

Язык: Русский

Опубликовано: около 4 лет назад

Просмотры: 119

16.38

Олег Снегирев

Дальневосточный федеральный университет (ДВФУ)

Комментировать 0

Рецензировать 1

Скачать - 3,7 МБ

Поделиться работой

МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ Федеральное государственное автономное образовательное учреждение высшего образования «Дальневосточный федеральный университет» (ДВФУ) Инженерная школа Кафедра технологий промышленного производства СНЕГИРЕВ ОЛЕГ ЮРЬЕВИЧ УПРАВЛЕНИЕ ТЕХНОЛОГИЧЕСКИМ ПРОЦЕССОМ РЕКТИФИКАЦИИ НА ОСНОВЕ ВИРТУАЛЬНЫХ АНАЛИЗАТОРОВ МАГИСТЕРСКАЯ ДИССЕРТАЦИЯ по образовательной программе подготовки магистров по направлению 15.04.04 «Автоматизация технологических процессов и производств» магистерская программа «Автоматизация технологических процессов и производств (в промышленности)» г. Владивосток 2018

Автор работы ______________ Снегирев О.Ю. (подпись) Группа М3215 " _____ " ______________ 2018 г. Руководитель магистерской работы, профессор, док.тех.наук ___________________ А.Ю. Торгашов (подпись) «____» _________________ 2018 г. Рецензент к.т.н., доцент __________________________ Д.А. Назаров «____» _________________ 2018 г. Защищена в ГЭК с оценкой _____________________________ Секретарь ГЭК _________________ С.Е. Коровин «Допустить к защите» Заведующий кафедрой к.т.н., доцент _____________________________ К.В.Змеу " ___ " ____________ 2018г. (подпись) «____» _________________ 2018 г. Регистрационный № ________ ___________________ подпись И.О.Фамилия « _____» ___________ 20 г. Руководитель ОП, к.т.н., доцент ________________________К.В. Змеу (подпись) «____» _________________ 2018 г. 2

МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ Федеральное государственное автономное образовательное учреждение высшего образования «Дальневосточный федеральный университет» (ДВФУ) ИНЖЕНЕРНАЯ ШКОЛА Кафедра технологий промышленного производства УТВЕРЖДЕНО Руководитель ОП, к.т.н., доцент (должность, ученое звание) _____________ Змеу К.В. (подпись) (Ф.И.О.) « 21 » октября 2016 г. Заведующий кафедрой, к.т.н., доцент (должность, ученое звание) ____________________ (подпись) Змеу К.В. (Ф.И.О.) « 21 » октября 2016 г. ЗАДАНИЕ на выпускную квалификационную работу – магистерскую диссертацию Студенту Снегиреву Олегу Юрьевичу группы М3215 1) На тему: Управление технологическим процессом ректификации на основе виртуальных анализаторов 2) Основания для разработки: График учебного процесса и приказ директора. 3) Вопросы, подлежащие разработке (исследованию): Анализ существующих подходов к разработке систем управления процессом ректификации. Исследование наблюдаемости технологического объекта Разработка адаптивного виртуального анализатора Определение передаточной функции Моделирование и расчет коэффициентов ПИД-регулятора 4) Основные источники информации и прочее, используемые для разработки темы: Литературные источники Интернет источники 5) Технические требования (параметры): Среднеквадратическая ошибка – не более 0,5; время переходного процесса – не более 6 часов Перечень графического материала (с точным указанием чертежей): Презентация 3

КАЛЕНДАРНЫЙ ГРАФИК ВЫПОЛНЕНИЯ ВКР Наименование этапов магистерской работы Сроки выполнения этапов работы 1. Изучение процесса ректификации. 11.12.2016 2. Исследование наблюдаемости технологического объекта. 22.05.2017 3. Разработка адаптивного виртуального анализатора. 13.11.2017 4. Определение передаточной функции. 19.01.2018 5. Разработка модели ПИДрегулятора. 23.04.2018 6. Получение коэффициентов ПИД-регулятора различными методами. 4.06.2018 7. Оформление пояснительной записки. 15.06.2018 Срок представления работы 6 июля 2018 г. Дата выдачи задания 17 октября 2017 г. Руководитель ВКР _______________ Прим. Торгашов Андрей Юрьевич (подпись) Задание получил _______________ Снегирев Олег Юрьевич (подпись) 4

Оглавление АННОТАЦИЯ.......................................................................................................... 8 ОБОЗНАЧЕНИЯ И СОКРАЩЕНИЯ .................................................................... 9 ВВЕДЕНИЕ ............................................................................................................ 11 1. ОПИСАНИЕ ТЕХНОЛОГИЧЕСКОГО ПРОЦЕССА РЕКТИФИКАЦИИ 12 2. СУЩЕСТВУЮЩИЕ ПОДХОДЫ К РАЗРАБОТКЕ СИСТЕМ УПРАВЛЕНИЯ ПРОЦЕССОМ РЕКТИФИКАЦИИ ......................................... 14 3. МАТЕМАТИЧЕСКОЕ ОПИСАНИЕ ТЕХНОЛОГИЧЕСКОГО ОБЪЕКТА 16 4. ИССЛЕДОВАНИЕ НАБЛЮДАЕМОСТИ ТЕХНОЛОГИЧЕСКОГО ОБЪЕКТА............................................................................................................... 18 4.1. Формирование вектора состояний ....................................................... 19 4.2. Получение матрицы А пространства состояний ................................ 20 4.3. Получение матрицы С пространства состояний ................................ 22 4.4. Исследование влияния технологических параметров массообменного объекта управления на наблюдаемость объекта ............................................ 23 4.4.1. Влияние флегмового числа и расхода верхнего продукта.......... 23 4.4.2. Влияние номера ступени разделения, на которую подается сырье 29 4.5. 5. Наблюдаемость объекта........................................................................ 30 РАЗРАБОТКА ВИРТУАЛЬНОГО АНАЛИЗАТОРА ............................... 31 5.1. Выбор наиболее информативных входных переменных .................. 31 5.1.1. Матрица парных коэффициентов корреляции ............................. 31 5.1.2. Метод регрессии «лассо» ............................................................... 32 5

5.2. Построение модели виртуального анализатора: метод робастной регрессии (взвешенный МНК - WLS) ............................................................. 34 5.3. Адаптивный виртуальный анализатор ................................................ 37 5.3.1. Принцип работы разрабатываемого адаптивного ВА ................. 37 5.3.2. Выбор оптимальной ширины окна ВА ......................................... 39 5.3.2.1. Критерий ошибки прогнозирования ........................................ 39 5.3.2.2. Критерий изменчивости параметров модели.......................... 43 5.3.2.3. Критерий соответствия (адекватности) модели промышленным данным........................................................................... 46 5.3.2.4. Векторный критерий ................................................................. 50 5.3.2.5. Сравнение критериев ................................................................. 52 5.3.3. 6. Сравнение адаптивного и неадаптивного ВА .............................. 54 ИСПОЛЬЗОВАНИЕ ВИРТУАЛЬНОГО АНАЛИЗАТОРА В СИСТЕМЕ УПРАВЛЕНИЯ ТЕХНОЛОГИЧЕСКИМ ПРОЦЕССОМ РЕКТИФИКАЦИИ 56 6.1. Модель динамики технологического объекта (PDS HONEYWELL определение ПФ) ............................................................................................... 57 6.2. Синтез ПИД регуляторов (CAI_MSBE, CAI_DIME) на основе виртуальных анализаторов ............................................................................... 68 6.2.1. Базовое ПИД управление (PID) ..................................................... 70 6.2.2. Пропорциональный (PV) и дифференциальный алгоритм ПИД управления (I-PD).......................................................................................... 71 6.2.3. Дифференциальный алгоритм ПИД управления (PI-D) ............. 71 6.2.4. Автоматическое определение ........................................................ 72 6.2.5. Автоматическое определение 2 ..................................................... 72 6.3. Настройка ПИД-регуляторов ............................................................... 73 6.3.1. Метод Кохен-Куна .......................................................................... 73 6

6.3.2. Метод на основе целевой функции ............................................... 75 6.3.3. Метод на основе внутренней модели(IMC) ................................. 79 6.3.4. Метод оптимизации ........................................................................ 81 6.3.5. Сравнение методов ......................................................................... 85 ЗАКЛЮЧЕНИЕ ..................................................................................................... 87 Приложение A Лицензия продукта Honeywell Profit Suit ................................. 92 7

АННОТАЦИЯ Объем – 87 страниц, 52 рисунка, 16 таблиц, 1 приложение, список литературы: 29 источников, из них 11 англоязычных. Структура ВКР включает введение, обзорную часть, исследовательскую часть, заключение и список литературы. Во введении представлена актуальность рассматриваемого вопроса, обозначены цели и задачи, рассматриваемые в работе. В обзорной части рассмотрен процесс ректификации, существующие системы управления процессом ректификации. В исследовательской части проведено исследование наблюдаемости технологического объекта ректификации. Разработан адаптивный виртуальный анализатор по методу “движущегося окна”, предложены критерии выбора оптимальной ширины окна адаптивного виртуального анализатора. Рассмотрен процесс определения передаточной функции в PDS HONEYWELL. Построена модель ПИД-регулятора и получены его коэффициенты различными методами. В заключении приводятся результаты проделанных работ. Ключевые слова: наблюдаемость, виртуальный анализатор, ПИДрегулятор, передаточная функция, ректификация, система управления 8

ОБОЗНАЧЕНИЯ И СОКРАЩЕНИЯ В настоящей работе применены следующие обозначения и сокращения: АСУТП – процессом, автоматизированная система управления технологическим ББФ – бутан-бутиленовая фракция, ВА – виртуальный анализатор, МНК – метод наименьших квадратов, МТБЭ – метил-трет-бутиловый эфир, ПИД – пропорционально-интегрально-дифференцирующий, ТО – технологический объект, ТП – технологический процесс, FIR – finite impulse response, IAE – integral of the absolute value of error, IRLS – iteratively reweighted least-squares, ISE – integral of the squared error, LASSO – least absolute shrinkage and selection operator, MSE – mean squared error, WLS – weighted least squares estimation, c  количество компонентов на тарелке колонны, K  число ступеней разделения, ~ x j , i  концентрация i -ого компонента на j -ой ступени разделения, ~ xF , i  концентрация i -ого компонента в сырье, f  номер ступени подачи сырья, R  расход орошения, F  расход сырья, D  расход верхнего продукта,  i  относительная летучесть i -ого компонента, 9

Pj  давление на j -й тарелке,  j – вычисляется по формуле, RR  флегмовое число, QB  матрица наблюдаемости,  Q   число обусловленности матрицы Q, xy  среднее значение произведения переменных x и y, x  среднее значение переменной x,  x   среднеквадратическое отклонение переменной x, KR  матрица парных коэффициентов корреляции,   целевая функция,   кривая влияния, we  весовая функция, y k - значение выходной переменной в k-ый момент времени, Bh yˆ k k 1 -спрогнозированное значение выходной переменной в k-ый момент времени, bˆi , k , h - i-ый параметр регрессионной модели в k-ый момент времени при ширине окна виртуального анализатора h, bi , h - среднее арифметическое значений i-ого параметра регрессионной модели при передвижении виртуального анализатора с шириной окна h, J Bˆ , h  критерий изменчивости параметров модели, J RMSE  критерий ошибки прогнозирования, J R 2  критерий соответствия модели данным, ~ J  векторный критерий. 10

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

В соответствии с целью были поставлены следующие задачи:  исследование наблюдаемости процесса ректификации с использованием упрощенной математической модели процесса, с целью обоснования необходимости использования виртуальных анализаторов в системе управления процессом ректификации;  определение наиболее информативных входов для построения моделей виртуальных анализаторов;  построение моделей адаптивных виртуальных анализаторов, основанных на методе «движущегося окна», на основе промышленных данных;  выбор оптимальной длины обучающей выборки для адаптивных виртуальных анализаторов по концентрации заданных веществ в выходном продукте;  определение передаточных функций изменения концентраций заданных веществ в выходном продукте с использованием Honeywell Profit Design Studio;  построение и настройка ПИД регуляторов на основе виртуальных анализаторов. 1. ОПИСАНИЕ ТЕХНОЛОГИЧЕСКОГО ПРОЦЕССА РЕКТИФИКАЦИИ Процесс получения МТБЭ основан на реакции селективного взаимодействия изобутилена, входящего в состав ББФ, с метанолом: (1.1) Синтез и выделение МТБЭ осуществляется в реакционно- ректификационном аппарате, состоящем из двух ректификационных и одной реакционной зон. Исходная ББФ и метанол поступают в реактор Р-351, предназначенный для синтеза основного количества эфира. 12

В зависимости от концентрации изобутилена в исходной С4-фракции, тепло реакции в реакторе Р-351 расходуется на разогрев реакционной массы или дополнительно на испарение части реакционной массы. Процесс испарения контролируется давлением в реакторе. Реактор Р-351 представляет собой полый цилиндрический аппарат, заполненный катализатором. Реакционная масса выводится из реактора Р-351 с верха аппарата одним или двумя потоками: в паровой и жидкой фазе. Реакционно-ректификационный аппарат Кр-351/1 включает три зоны:  верхнюю ректификационную зону (для отделения непрореагировавших углеводородов С4, от метанола и эфиров);  среднюю реакционно-ректификационную зону, заполненную катализатором (для синтеза эфиров и их вывода из зоны реакции).  нижнюю ректификационную зону (для отделения МТБЭ от углеводородов С4 и метанола). Катализатор в аппарате Кр-351/1 расположен в виде слоев на опорнораспределительных тарелках специальной конструкции. Реакционная масса из реактора Р-351 поступает в аппарат Кр-351/1 под слой катализатора. Наверх катализатора в Кр-351/1 подается метанол. Сверху аппарата Кр-351/1 отбирается С4-фракция, которая подается в колонну Кр351/2 водной отмывки углеводородов от метанола. Кубовый продукт колонны Кр-351/1 - товарный МТБЭ выводится с установки. В верхнюю часть колонны Кр-351/2 подается вода. Сверху колонны Кр351/2 отбирается отработанная С4-фракция. Промывная вода с метанолом из куба колонны Кр-351/2 подается в качестве питания в колонну К-1, предназначенную для отгонки метанола от воды. Метанол, отбираемый с верха колонны К-1, возвращается в емкость со свежим метанолом. Вода из куба колонны К-1 подается в верхнюю часть колонны Кр-351/2. Принципиальная технологическая схема процесса получения МТБЭ представлена на рисунке 1.1. 13

Рисунок 1.1 – Принципиальная технологическая схема процесса получения МТБЭ Выходным продуктом определенные примеси. не является чистый МТБЭ, имеются Количество этих примесей должно быть ниже определенного уровня. Так для метил-втор-бутилового эфира и димеров стоит требование, чтобы концентрация данных веществ в выходном продукте была ниже 2%. 2. СУЩЕСТВУЮЩИЕ ПОДХОДЫ К РАЗРАБОТКЕ СИСТЕМ УПРАВЛЕНИЯ ПРОЦЕССОМ РЕКТИФИКАЦИИ Известен способ автоматического управления процессом ректификации путем изменения расхода орошения ректификационной колонны в зависимости от значения температуры верхней части колонны [2]. 14

Недостаток известного способа заключается в низкой точности регулирования качества дистиллята ректификационной колонны вследствие низкого быстродействия системы управления. Известен способ автоматического управления процессом ректификации путем изменения расхода орошения в зависимости от температуры и давления верхней части колонны, молекулярного веса нефтепродукта [3]. Недостатком регулирования известного качества способа нефтепродукта, является так как низкая точность нефтепродукты многокомпонентные смеси углеводородов с одинаковым молекулярным весом могут иметь различные показатели качества. Известен способ автоматического регулирования процесса ректификации [4] путем задания температурного профиля колонны и изменения подачи теплоносителя и орошения в зависимости от изменения температурного профиля ректификационной колонны, который реализует измерение давления в верхней и нижней частях колонны и, в зависимости от измеренных значений давлений, корректирование температурного профиля колонны, определение высоты участка колонны, на котором температура не меньше температуры кипения кубового продукта заданного состава, и высоты участка колонны, на котором температура не превышает температуры кипения дистиллята заданного состава, расчет скорости изменения температуры по высоте колонны и, в зависимости от значений определенных высот участков колонны и скорости изменения температуры по высоте колонны параметров, корректирование расходов теплоносителя и орошения, расхода и теплосодержания питающей смеси. Данный способ регулирования процесса ректификации является наиболее близким к предлагаемому решению. Недостатком способа является увеличение динамической погрешности процесса управления и время регулирования температуры питающей тарелки. Также в известном способе отсутствует компенсация инерционности процесса ректификации, что приводит к ухудшению качества регулирования. 15

Целью предлагаемого решения является создание способа управления технологическим процессом ректификации, достигающего высокой точности регулирования качества нефтепродукта, позволяющего оперативно реагировать на отклонения технологического процесса. Предлагаемый способ управления процессом ректификации заключается в использовании в структуре управления блока виртуального анализатора. Получая на вход показания датчиков температуры, давления и уровня сырья на ступенях разделения, а также результаты лабораторного анализа выходного продукта, в блоке виртуального анализатора строится математическая модель ректификационной колонны. В связи с тем, что лабораторный анализ выходного продукта производится с интервалом от 8 часов, управление процессом ректификации основывается на полученной математической модели. Данная модель считается актуальной до получения результатов лабораторного анализа выходного продукта, после которого проверяется актуальность модели, и, если математическая модель технологического процесса устарела, пересчитываем ее на основе последних данных. На основе математической модели задается требуемый температурный профиль колонны и изменение подачи теплоносителя и орошения. 3. МАТЕМАТИЧЕСКОЕ ОПИСАНИЕ ТЕХНОЛОГИЧЕСКОГО ОБЪЕКТА Рассматривается многоступенчатый процесс многокомпонентной ректификации, протекающий в технологическом аппарате колонного типа с jмя ступенями разделения и описываемый следующей системой: При j  1 : d~ x j,i dt  R  D  i ~x j 1,i c 1 c   ~ x j 1, i  i   c   R~ x j , i  D~ x j,i ; (3.1) i 1 16

При j  2... f  1: d~ x j,i dt  R~ x j 1, i  R  D  i ~x j 1, i c 1  c   ~x j 1, i  i   c   (3.2) i 1  R~ x j,i  R  D  i ~x j ,i c 1  c   ~x j , i  i   c  ; i 1 При j  f : d~ x j,i dt  R~ x j 1,i  R  D αi ~x j 1,i c 1 αc   ~ x j 1,i αi  αc   (3.3) i 1  R  F ~ x j,i  R  D αi ~x j,i c 1 αc   ~ x j,i αi  αc   F~ x F,i ; i 1 При j  f  1...K  1: d~ x j,i dt  R  F ~ x j 1,i  R  D αi ~x j 1,i c 1 αc   ~ x j 1,i αi  αc   (3.4) i 1  R  F ~ x j,i  R  D αi ~x j,i c 1 αc   ~ x j,i αi  αc  ; i 1 При j  K : d~ x j,i dt  R  F ~ x K 1,i  F  D ~ x K,i  R  D αi ~xK,i c 1 αc   ~ x K,i αi  αc  ; (3.5) i 1 где i  1...c; c  количество компонентов на тарелке колонны; K  число x  концентрация i -ого компонента на j -ой ступени ступеней разделения; ~ j,i разделения; ~ xF , i  концентрация i -ого компонента в сырье; f  номер ступени подачи сырья; R  расход орошения; F  расход сырья; D  расход верхнего продукта;  i  относительная летучесть i -ого компонента. 17

Изменение температуры на j -ой ступени разделения описывается выражением:  Pj  Bb  Cb lg    Cb Ab  j    Tj  ,  Pj  Ab  lg    j    (3.6) где Ab , Bb,Cb  константы Антуана для базового состава; Pj  давление на j -й тарелке;  j – вычисляется по формуле: c 1 (3.7)  j   c   ~x ji  i   c . i 1 Для упрощения последующего использования уравнений (3.1) - (3.7) в исследованиях введем следующее обозначение: d~ xl , p q  dt (3.8)  q  1  q  1    1 , p  q  c  1  где l    c  1 .  c  1  Правую часть функции изменения концентрации компонента можно представить в виде суммы из нескольких слагаемых:   d~ xl , p q   Sl 1, p  Sl 1, p  Sl , p  Sl , p dt где S  K~ x ,S  K~ x , l, p  Sl , p  l, p l 1, p l 1, p R  D αi ~xl , p c 1 αc   ~ xl , p αi  αc  i 1 (3.9)  , Sl 1, p  R  D αi ~xl 1, p c 1 αc   ~ xl 1, p αi  αc  . i 1 4. ИССЛЕДОВАНИЕ НАБЛЮДАЕМОСТИ ТЕХНОЛОГИЧЕСКОГО ОБЪЕКТА 18

Для осуществления управления необходимо иметь информацию о текущем состоянии системы, то есть о значениях вектора состояния 𝑋 в каждый момент времени. Состояние 𝑋 называется наблюдаемым, если в момент времени наблюдения можно однозначно определить 𝑋 по данным измерения входных u(t) и выходных y(t). переменных на конечном интервале времени. Для оценки наблюдаемости объекта необходимо знать матрицу системы A и матрицу выхода C пространства состояний [5]. Формирование вектора состояний 4.1. Матрица состояний включает в себя значения концентраций каждого из c компонентов на каждой из K тарелок. Тогда матрица состояний будет иметь размерность (K x c): ~ x1,1  ~ x1, c    ~ x      ~ ~   x K ,1  x K , c  K  c Для уменьшения количества (4.1) элементов вектора состояний Х воспользуемся тем, что на каждой ступени разделения сумма концентраций всех компонентов равна 1: c c c i 1 i 1 i 1  ~x1, i   ~x...,i   ~xK , i  1 (4.2) Тогда можно составить вектор состояний Х, содержащий Q  K  c  1 элементов:  X1        X  X q ,      X   Q (4.3) где q=1, …, Q. 19

Каждому элементу вектора состояний X будет соответствовать определенный элемент матрицы состояний ~ x: Xq  ~ x  q 1   q 1  1,  q  c 1        c 1     c 1    (4.4) , где q=1…Q – номер элемента вектора состояний Х, c – количество компонентов смеси, a  maxn  Z | n  a  функция пол, возвращающее наибольшее целое число, меньшее или равное a [6]. 4.2. Получение матрицы А пространства состояний Матрица системы А пространства состояний является матрицей частных производных функций системы от элементов вектора состояний. Функции системы представляют собой изменения концентрации каждого из (c-1) компонентов на каждой из K ступеней разделения. Всего мы имеем Q=K  (c1) функций. Тогда матрица системы А является квадратной матрицей размерности Q  Q:  1   X 1     k A  X 1     Q  X  1 1 X g    k  X g    Q  X g       1   X Q     k   X Q     Q  X Q  (4.5) Q Q Рассмотрим функции изменения концентраций компонентов i. Они различны в зависимости от номера ступени разделения j и от номера ступени f, на которую подается сырье. Исследуя функции изменения концентраций компонентов на ступени разделения, можно заметить, что концентрация компонента на j ступени зависит от концентрации компонентов на j ступени и концентраций компонентов на (j - 1) и (j + 1) ступенях разделения. 20

Для нахождения значений элементов матрицы системы необходимо найти частные производные функций изменения концентраций от каждой концентрации. Тогда значения коэффициентов матрицы системы Ak,g будут ненулевыми для элементов матрицы, индексы которых удовлетворяют системе неравенств:   j  1  c  1  k  j  c  1   j  2  c  1  g   j  1  c  1 (4.6) Значение ненулевых элементов матрицы А равны:    k Sl 1, p Sl 1, p Sl , p Sl , p Ak , g   ~  ~  ~  ~ X x x x x g m, n m, n m, n (4.7) m, n  k  1  k  1    1 , p  k  c  1  где l    c  1 ,  c  1   g  1  k  1   m  1 , n  g  c  1   c  1 .  c  1  Частная производная слагаемых Sl , p от ~ xm, n равна:  S l , p  K, m  l  n  p  ~  xm, n  S  l , p  0, m  l  n  p ~   xm, n (4.8)  Частная производная слагаемых S l , p от ~ xm, n равна: 21

c 1   R  D  p   c   ~xl , p  p   c S l , p p 1   2 ~ xm, n c 1    c   ~ xl , p  p   c   p  1   ~ R  D  p xl , p  p   c  2 c 1   ~   c   xl , p  p   c    p 1    R  D  p ~xl , p  p   c S l , p  2 ~ xm, n c 1   ~   c   xl , p  p   c    p 1    S l , p 0 ~ x         ml n p (4.9) ml n p (4.10) ml n p (4.11)      m, n 4.3. Получение матрицы С пространства состояний Матрица выхода С пространства состояний является матрицей частных производных функций изменения температуры от концентрации Q компонентов. Всего мы имеем K функций. Матрица выхода С является матрицей размерности K  Q:  T1   X 1    Th C   X 1    TK  X  1      T1 X g  Th X g  TK X g      T1   X Q    Th   X Q    TK  X Q  (4.12) K Q Рассмотрим функции изменения температуры на j ступени разделения. Исследуя функции изменения температуры на j ступени, можно заметить, что температура зависит только от концентрации компонентов на 𝑗 ступени. Тогда значения коэффициентов матрицы состояний 𝐶ℎ,𝑔 будут ненулевыми для элементов матрицы, индексы которых удовлетворяет системе: 22

h j    j  1  c  1  g  j  c  1 (4.13) Значение ненулевых элементов матрицы C равны: C h, g    Th  Bb ln 10  g  h 1 c 1   c  , 2 X g   P   Ab ln 10  ln  h    h      h   (4.14) c 1 где  h   c   X h 1c 1 i   i   c  i 1 4.4. Исследование влияния технологических параметров массообменного объекта управления на наблюдаемость объекта 4.4.1. Влияние флегмового числа и расхода верхнего продукта Проведем исследования наблюдаемости технологического объекта с числом ступеней разделения K=20 и сырьем, состоящим из c=4 компонентов. Изменение расхода верхнего продукта D в интервале {100 … 300} с шагом 10. Изменение флегмового числа RR в интервале {0.5 … 3} с шагом 0.1. Флегмовое число [7] представляет собой отношение количества жидкости R, стекающей с любой тарелки в концентрационной секции колонны, к количеству отбираемого верхнего продукта D. Флегмовое число записывается отношением: RR  (4.15) R D Оценить наблюдаемость объекта можно с помощью матрицы наблюдаемости:  C     CA  T QB   CA 2   C CA CA 2  CA n 1       CA n 1      (4.16) 23

Система наблюдаема, если и только если ранг матрицы наблюдаемости Rang(QB) [8] равен размерности вектора состояний Q: RangQB   Q На рисунке (4.17) 4.1 представлена поверхность рангов матрицы наблюдаемости в зависимости от расхода верхнего продукта и флегмового числа в заданных интервалах. Видно, что рассматриваемый объект ненаблюдаемый при заданных параметрах. Можно оценить степень наблюдаемости объекта. Используя критерий оценки наблюдаемости через ранг матрицы наблюдаемости, можно выделить область, на которой степень наблюдаемости будет выше. Для этого выделим область, где ранг матрицы наблюдаемости больше либо равен 9. На рисунке 4.2 представлена данная область. Рисунок 4.1 – Изменение ранга матрицы наблюдаемости при изменении значений RR и D 24

Рисунок 4.2 – Рекомендуемая область D и RR по критерию ранга матрицы наблюдаемости Существует способ оценки наблюдаемости объекта с использованием числа обусловленности матрицы наблюдаемости  QB  [9] [10] [11]. Чем меньше  QB , тем выше степень наблюдаемости рассматриваемого объекта. На рисунке 4.3 представлен график зависимости числа обусловленности матрицы наблюдаемости от параметров массообменного процесса. Для упрощения восприятия полученных значений числа обусловленности представлены в виде:  g QB   log 10  QB   100% maxlog 10  QB  (4.18) где maxlog10  QB  – максимальное значение десятичного логарифма числа обусловленности матрицы наблюдаемости на рассматриваемом интервале изменения флегмового числа RR и расхода верхнего продукта D. 25

Рисунок 4.3 – Изменение числа обусловленности матрицы наблюдаемости при изменении значений RR и D Выберем область изменения значений RR и D, которую можно рекомендовать к использованию в рассматриваемом технологическом объекте, где  g QB   90% . На рисунке 4.4 представлена выбранная область. Рисунок 4.4 – Рекомендуемая область D и RR по критерию числа обусловленности матрицы наблюдаемости 26

Еще один критерий, по которому можно оценить наблюдаемость объекта, использует понятие грамиана наблюдаемости [11] [9]. Согласно принципу дуальности Калмана в рассмотрение вводится дуальная система, переменные состояния которой образуют свою матрицу Грама в сходных рассмотренным выше условиях. Грамиан наблюдаемости можно найти из следующего уравнения Ляпунова: GA  A' G  c' c  0 (4.19) В оценке степени наблюдаемости объекта рассматривается число обусловленности грамиана наблюдаемости  GB  . Чем меньше  GB  , тем выше степень наблюдаемости рассматриваемого объекта. На рисунке 4.5 представлен график зависимости числа обусловленности грамиана наблюдаемости от параметров массообменного процесса. Для упрощения восприятия полученных значений числа обусловленности представлены в виде:  g GB   log 10  GB   100% maxlog 10  GB  (4.20) где maxlog10  GB  – максимальное значение десятичного логарифма числа обусловленности грамиана наблюдаемости на рассматриваемом интервале изменения флегмового числа RR и расхода верхнего продукта D. 27

Рисунок 4.5 – Изменение числа обусловленности грамиана наблюдаемости при изменении значений RR и D Выберем область изменения значений RR и D, которую можно рекомендовать к использованию в рассматриваемом технологическом объекте, где 𝜇𝑔 (𝐺𝑏 ) < 90%. На рисунке 4.6 представлена выбранная область. Рисунок 4.6 – Рекомендуемая область D и RR по критерию числа обусловленности грамиана наблюдаемости 28

Рекомендуемые области по трем представленным критериям оценки наблюдаемости объекта имеют схожие черты, но отличаются. Найдем такую область, которая будет удовлетворять всем вышеописанным критериям. На рисунке 4.7 представлена полученная область. Рисунок 4.7 – Рекомендуемая область D и RR 4.4.2. Влияние номера ступени разделения, на которую подается сырье Ранг матрицы наблюдаемости также зависит от номера f ступени разделения, на которую подается сырье. На рисунке 4.8 представлены графики изменения максимальных и минимальных значений ранга матрицы наблюдаемости для рассматриваемого технологического объекта в заданной области изменения расхода выходного продукта и флегмового числа. 29

Рисунок 4.8 – Изменение максимальных и минимальных значений ранга матрицы наблюдаемости при изменении номера ступени подачи сырья Из графиков рисунка 4.8 видно, что при увеличении номера ступени разделения, на которую подается сырье, уменьшаются значения ранга матрицы наблюдаемости. 4.5. Наблюдаемость объекта Исследование влияния технологических параметров непрерывного ТО многокомпонентной ректификации на наблюдаемость показало, что объект не полностью наблюдаемый. Количество наблюдаемых компонент вектора состояния различно в зависимости от значений флегмового числа и количества отбираемого верхнего продукта, а также номера ступени разделения, на которую подается сырье. Для массообменного объекта с 20 ступенями разделения и четырехкомпонентной исходной смеси веществ можно рекомендовать область значений флегмового числа и количества отбираемого верхнего продукта представленную на рисунке 4.7. Для достижения максимального количества наблюдаемых комбинаций элементов вектора состояния ТО рекомендуется 30

выбирать ступень разделения, на которую подается сырье, с наименьшим возможным номером. В связи с тем, что объект не полностью наблюдаемый использование наблюдателей невозможно, поэтому возникает необходимость в разработке виртуального анализатора. 5. РАЗРАБОТКА ВИРТУАЛЬНОГО АНАЛИЗАТОРА Выбор наиболее информативных входных переменных 5.1. Для создания модели процесса необходимо отобрать ряд независимых переменных, которые наиболее влияют на значение зависимой переменной. Для этого составим матрицу парных коэффициентов корреляции. 5.1.1. Матрица парных коэффициентов корреляции Корреляция или корреляционная зависимость – статическая взаимосвязь двух и более случайных величин. Математической мерой корреляции двух случайных величин служит коэффициент корреляции R. Для вычисления значения коэффициента корреляции переменных x и y можно воспользоваться формулой: Rxy  где xy  xy  x y ,  x   y  (5.1)  xi yi - среднее значение произведения переменных x и y при n числе n наблюдений, x   xi и y   yi - среднее значение переменных x и y при n n n 2 2 yi2 xi2   числе наблюдений соответственно,   x   y  x и  y  n n среднеквадратическое отклонение переменных x и y соответственно. 31

Матрица парных коэффициентов корреляции [12] для зависимой переменной Y и ряда независимых переменных X1 , X 2 ,, X m будет иметь вид: R yy R yx1 KR   R yx m  R yx m  Rx1 x m    Rx m x m R yx1 Rx1 x1  Rx1 x m (5.2) Для нахождения матрицы парных коэффициентов корреляции KR составим матрицу A: 1 Y1 X 11  1 Y2 X 12  X m 2 A  1 Y3   X 13  X m 3 1 Yn X 1n  X m n   X m1  (5.3) n  m  2  где Yi , X 1i ,...,X m i - значение переменных при i  1,...,n наблюдении. В результате умножения транспонированной матрицы AT на матрицу A получим матрицу, значения элементов которой будут соответствовать: С использованием формулы (5.1) и значений матрицы AT A можно вычислить значения элементов матрицы парных коэффициентов корреляции KR. По первому столбцу матрицы KR можно выявить такие независимые переменные, которые наиболее сильно связаны с зависимой переменной Y. Чем ближе значение коэффициента корреляции к 1 (или -1), тем сильнее взаимосвязаны переменные [13]. 5.1.2. Метод регрессии «лассо» Также для отбора информативных переменных из ряда данных можно воспользоваться регуляризацией. Регуляризация – общий метод, 32

заключающийся в наложении дополнительных ограничений на искомые параметры, которые могут предотвратить излишнюю сложность модели. Метод регрессии «лассо» (LASSO, Least Absolute Shrinkage and Selection Operator) [14] заключается во введении дополнительного слагаемого регуляризации в функционал оптимизации модели, что позволяет получать более устойчивое (сглаженное) решение, то есть менее зависящее от обучающих данных. Условие минимизации квадратов ошибки параметров ˆ выражается следующей формулой: 2  n m   ˆ  arg min    yi    j xij     j 1   i 1  ,   (5.4) где  - параметр регуляризации, имеющий смысл штрафа за сложность. При значении параметра   0 , лассо-регрессия сводится к методу наименьших квадратов, а при увеличении  формируемая модель становится все более «лаконичной», пока не превратится в нуль-модель. Оптимальной величине  соответствует минимальная ошибка прогноза ŷi на наблюдениях, не участвовавших в построении модели. Таким образом достигается компромисс между размерностью используемого пространства переменных и ошибкой регрессии, выраженного суммой абсолютных значений коэффициентов  . В ходе минимизации некоторые коэффициенты становятся нулевыми, что и определяет отбор информативных переменных. В среде Matlab реализован метод «лассо». Команда [ B , FitInfo ] = lasso ( X , Y ) возвращает коэффициенты регрессионной модели, полученных методом наименьших квадратов, для набора значений параметра регуляризации  и информацию о состоянии модели (значение параметра регуляризации  , значение среднеквадратической ошибки MSE, количество переменных участвующих в модели и др.) для зависимой переменной Y и ряда независимых переменных X. 33

Выбрав из полученных данных минимальное значение среднеквадратической ошибки MSE и соответствующие ему коэффициенты регрессионной модели для всех независимых переменных, можно получить набор наиболее информативных независимых переменных. Для этого необходимо отобрать переменные, коэффициенты регрессионной модели которых ненулевые. 5.2. Построение модели виртуального анализатора: метод робастной регрессии (взвешенный МНК - WLS) После того, как были отобраны наиболее информативные входные переменные, необходимо построить регрессионную модель с помощью методов линейной оценки зависимости выходной переменной от входных. Можно воспользоваться наиболее простым и часто используемым методом – методом наименьших квадратов [15]. Но линейная оценка методом наименьших квадратов может вести себя неудовлетворительно, когда ошибка не нормально распределена, в частности когда ошибки сильно разбросаны. Одним из средств является удаления влиятельных наблюдений из подмножества. Другой подход к решению данной проблемы, называемый робастной регрессией [16], заключается в использовании критерия подгонки, который менее чувствительный к необычным данным, чем метод наименьших квадратов. Наиболее часто используемым методом робастной регрессии является метод M-оценок. Этот класс оценок можно рассматривать как обобщение оценки максимального правдоподобия. Рассматривается линейная модель yi    1 xi1   2 xi 2  ...   k xik   i  X i '    i (5.5) для i-ой из n наблюдений. Давая оценку b для  , получим модель yˆ i    b1 xi1  b2 xi 2  ...  bk xik  ei  X i ' b (5.6) и разность 34

ei  yi  yˆ i (5.7) При M-оценке оценка b определяется минимизацией конкретной целевой функции над всеми b n n i 1 i 1 (5.8)   ei      yi  X i ' b  где функция  приносит вклад каждого остатка в целевую функцию. Подходящая функция  должна обладать следующими свойствами:  Неотрицательность,  e  0  Равенство 0 при равенстве аргумента 0,  0  0  Симметричность,  e    e  Монотонность в ei ,  ei    ei ' для ei  ei ' Пусть    ' является производной  .  называется кривой влияния. Дифференцируя целевую функцию относительно коэффициентов b и приняв частные производные равные 0, получим систему k+1 оценочных уравнений для коэффициентов: (5.9) n   yi  X i ' b X i '  0 i 1 Определим весовую функцию we   e / e , и wi  wei . Оценочные уравнения можно записать (5.10) n  wi  yi  X i ' b X i '  0 i 1 Решение этих оценочных уравнений эквивалентно взвешенному МНК, минимизирование  wi2 ei2 . Веса зависят от остатков, остатки зависят от оцененных коэффициентов, оцененные коэффициенты зависят от весов. Поэтому требуется итерационное решение называемое итеративно редуцированные наименьшие квадраты (iteratively reweighted least-squares, IRLS) 35

Выбрать начальные оцененные коэффициенты b (0) , по методу a) наименьших квадратов. В каждой итерации t вычислять остатки ei(t 1) и связанные веса b)   wi(t 1)  w ei(t 1) с предыдущей итерации. Вычислить новые оцененные коэффициенты взвешенной МНК по c) формуле:  b (t )  X 'W t 1 X  1 (5.11) X 'W (t 1) y   где X – матрица состояний модели, W (t 1)  diag wi(t 1) - текущая матрица весов. Этапы b и c повторять до тех пор, пока оцененные коэффициенты не прекратят изменяться. В качестве целевой функции можно использовать функцию Huber и функцию Bisquare, целевая функция и соответствующая ей весовая функция которых представлены в таблице 5.1. Таблица 5.1 – Целевые функции и соответствующие им весовые функции для метода робастной регрессии Метод Целевая функция Huber  0.5e 2  H e    k e  0.5k 2 Bisquar e Весовая функция  1 wH e    k / e e k e k  2  2 3  k e      1  1       B e    6    k       2 k 6  e k e k e k e k  2 2 e    1 wH e     k      0  e k e k Значение k называется настроечной константой. Меньшие значения k дают большую устойчивость к выбросам, но за счет более низкой эффективности, когда ошибки обычно распределяются. Константа настройки обычно выбирается для обеспечения достаточно высокой эффективности в 36

нормальном случае, в частности k  1.345 для Huber и k  4.685 для Bisquare (где  - стандартное отклонение ошибок) производят 95-процентную эффективность, когда ошибки нормально распределены и по-прежнему обеспечивают защиту от выбросов. 5.3. Адаптивный виртуальный анализатор Разрабатываемый виртуальный анализатор является адаптивным работает по методу «движущегося окна» [17]. 5.3.1. Принцип работы разрабатываемого адаптивного ВА Разбиваем интервал наблюдения во времени от 1 до N на интервал инициализации и интервал работы виртуального анализатора. Длина интервала инициализации равна значению ширины окна виртуального анализатора. Рисунок 5.1 – Принцип работы адаптивного виртуального анализатора На основе k-ой обучающей выборки: Ak  Yk , X k , (5.12) 37

где  y k  h 1    Yk     - матрица h наблюдений выходной переменной,  y  k   h1  X k  h 1   ~  Xk   - матрица h измерений входных переменных, h – длина   X  k   h1 обучающей выборки, строится регрессионная модель методом робастной регрессии, описанным ранее, с параметрами:  Bˆ kh  aˆ k , bˆ1, k , ,...,bˆm, k  (5.13) Затем полученная модель используется для прогнозирования значения выходной переменной yˆ k 1 Bˆ h по полученным значениям входных переменных k xk 1,1,...,xk 1, m : yˆ (5.14) m k 1 Bˆ kh  bˆ0, k   bˆi , k xk 1,i i 1 При получении истинного значения выходной переменной yk 1 , расчитывается ошибка прогнозирования значения моделью: Ek 1  yk 1  y (5.15) k 1 Bˆ kh После того как ошибка прогнозирования была рассчитана, в обучающую выборку A виртуального анализатора включают (k+1) наблюдение и исключают (k+1-h) наблюдение выходной и входных переменных: Ak 1  Yk 1 , X k 1,  xk  2  h,1  xk  2  h, m   yk  2  h         где Yk 1     , X k 1     x   y   k 1  h1  k 1,1  xk 1, m  (5.16) h m Расчитываются параметры регрессионной модели для новой обучающей выборки Bˆ kh1 и выполняется последовательность действий описанная выше. В ходе процесса, получая (k,…,N) наблюдения выходной и входных переменных, получим вектор ошибок прогнозирования: 38

 E h  Ekh , Ekh1 ,...,E Nh  (5.17) Показателем качества адаптивного виртуального анализатора примем среднеквадратическую ошибку [18], вычисляемую по формуле:  Eih  N Eh  2 ik (5.18) N  k 1 5.3.2. Выбор оптимальной ширины окна ВА Различным значениям ширины окна виртуального анализатора соответствует различные значения среднеквадратической ошибки. Поэтому возникает необходимость выбора значения оптимальной ширины окна. Для этого рассмотрим различные критерии. На практике принято следующее отношение длины обучающей выборки h0 к количеству входных переменных модели m: h0  m  10 (5.19) 5.3.2.1. Критерий ошибки прогнозирования Для выбора оптимального значения ширины окна виртуального анализатора можно использовать критерий отражающий качество адаптивного виртуального анализатора через ошибку прогнозирования значения выходной переменной. Показателем качества адаптивного виртуального анализатора может являться среднеквадратическая ошибка MSE [18]: 1 MSE  N M 2     yk  yˆ k Bˆ h  , k 1  k  M 1 N (5.21) Bh где y k - значение выходной переменной в k-ый момент времени, yˆ k k 1 спрогнозированное значение выходной переменной в k-ый момент времени. 39

Качество адаптивного виртуального анализатора будет тем лучше, чем меньше значение среднеквадратической ошибки. Более легко интерпретируемым показателем качества адаптивного виртуального анализатора является среднеквадратичная ошибка RMSE: 1 RMSE  N M     yk  yˆ k Bˆ h  k 1  k  M 1 N 2 (5.22) На практике используется упрощенный показатель оценки качества: AMSE  1 N M N  yk  yˆ k Bˆ h k  M 1 (5.23) k 1 На графиках рисунка 5.2 и рисунка 5.2 представлены графики зависимостей показателей lg(MSE), lg(RMSE) и lg(AMSE) от значения ширины окна адаптивного виртуального анализатора для метил-втор-бутилового эфира и димеров соответственно. Из графиков видно, что поведения показателей MSE, RMSE и AMSE идентичные, отличающиеся только порядком значений. Рисунок 5.2 – Сравнение MSE, AMSE и RMSE от ширины обучающей выборки для метил-втор-бутилового эфира 40

Рисунок 5.3 – Сравнение MSE, AMSE и RMSE от ширины обучающей выборки для димеров Таким образом, для выбора оптимального значения ширины окна адаптивного виртуального анализатора по критерию ошибки прогнозирования можно использовать либо MSE, либо AMSE, либо RMSE. Воспользуемся среднеквадратичной ошибкой RMSE. Критерий ошибки прогнозирования на основе показателя среднеквадратичной ошибки RMSE: J RMSE 1  N M     yk  yˆ k Bˆ h  k 1  k  M 1 N 2 (5.24) На рисунке 5.4 и рисунке 5.5 представлены графики изменения значения критерия ошибки прогнозирования в относительных единицах от ширины окна виртуального анализатора по концентрации МВБЭ и димеров соответственно. 41

Рисунок 5.4. Значение критерия ошибки прогнозирования от ширины окна виртуального анализатора по концентрации МВБЭ Рисунок 5.5. Значение критерия ошибки прогнозирования от ширины окна виртуального анализатора по концентрации димеров 42

Оптимальной ширине окна виртуального анализатора будет соответствовать ширина окна, при которой значение критерия ошибки прогнозирования будет минимальным. В таблице 5.2 представлены значения ширины окна виртуального анализатора по критерию ошибки прогнозирования для метил-вторбутилового эфира и димеров и значения среднеквадратичной ошибки и коэффициент детерминации адаптивного виртуального анализатора, соответствующие оптимальной ширине окна. Таблица 5.2 – Оптимальная ширина окна, среднеквадратичная ошибка и коэффициент детерминации по критерию ошибки прогнозирования Метил-втор-бутиловый Димеры эфир Оптимальная ширина окна 890 220 2.5187е-05 0.0011715 0.42087 0.0035104 hJ RMSE Среднеквадратичная ошибка RMSE Коэффициент детерминации R2 5.3.2.2. Критерий изменчивости параметров модели Для выбора оптимальной ширины окна адаптивного виртуального анализатора можно использовать критерий отражающий изменчивость параметров регрессионной модели в течение передвижения обучающей выборки по ряду данных. Для определения изменчивости параметров регрессионной модели используется среднеквадратическое отклонение оценок параметров модели [19]: Sbˆ i ,h  1 N M  bˆi, k , h  bi, h  , N 2 (5.25) k  M 1 43

где bˆi , k , h - i-ый параметр регрессионной модели в k-ый момент времени при ширине окна виртуального анализатора h, bi , h - среднее арифметическое значений i-ого параметра регрессионной модели при передвижении виртуального анализатора с шириной окна h. Критерий выбора оптимального значения ширины окна виртуального анализатора использует комплексную оценку среднеквадратических отклонений оценок параметров модели: (5.26) m J Bˆ , h   Sbˆ , h i i 0 На рисунке 5.6 и рисунке 5.7 представлены графики изменения значения критерия изменчивости параметров модели в относительных единицах от ширины окна виртуального анализатора для метил-втор-бутилового эфира и димеров соответственно. Оптимальная ширина окна виртуального анализатора будет соответствовать минимальному значению критерия. По графикам рисунков 5.6 и 5.7 видно, что функция изменения значения критерия изменчивости параметров модели от ширины окна является убывающей, содержащей локальные минимумы. Так как минимальному значению критерия будет соответствовать максимальное значение ширины окна на рассматриваемом интервале, что нецелесообразно для использования на практике, то будем считать, что оптимальной ширине окна виртуального анализатора будет соответствовать ширина окна локального минимума функции изменения значения критерия изменчивости параметров модели. 44

Рисунок 5.6 – Значение критерия изменчивости параметров от ширины окна виртуального анализатора для метил-втор-бутилового эфира Рисунок 5.7 – Значение критерия изменчивости параметров от ширины окна виртуального анализатора для димеров В таблице 5.3 представлены значения ширины окна виртуального анализатора по критерию изменчивости параметров для метил-втор45

бутилового эфира и димеров и значения среднеквадратичной ошибки и коэффициент детерминации адаптивного виртуального анализатора, соответствующие оптимальной ширине окна. Таблица 5.3 – Оптимальная ширина окна, среднеквадратичная ошибка и коэффициент детерминации по критерию изменчивости параметров модели Метил-втор-бутиловый Димеры эфир Оптимальная ширина окна 900 720 5.1915е-05 0.0018287 0.41925 0.11112 hJ ˆ B Среднеквадратичная ошибка RMSE Коэффициент детерминации R2 5.3.2.3. Критерий соответствия (адекватности) модели промышленным данным При выборе оптимального значения ширины окна виртуального анализатора можно использовать критерий соответствия модели данным. Для этого при оценке регрессионных моделей используют коэффициент детерминации [20]. Коэффициент детерминации показывает, какая доля вариации выходной переменной учтена в модели и обусловлена влиянием на нее входных переменных: 2     yk  yˆ k Bˆ h  k 1  R 2  1  k  M 1N , 2   yk  y  N (5.27) k  M 1 где yk y 1 N M - значение выходной переменной в k момент времени, N  yk - среднее арифметическое значений выходной переменной, k  M 1 46

yˆ k Bˆ kh1 - прогноз значения выходной переменной в k-ый момент времени на основе модели с вектором параметров h Bˆ k 1 . На рисунках 5.8 и 5.9 представлены графики зависимости коэффициента детерминации адаптивного виртуального анализатора от ширины окна. Оптимальному значению ширины окна будет соответствовать значение коэффициента детерминации наиболее близкое к 1 на рассматриваемом интервале. Рисунок 5.8 – Значение коэффициента детерминации от ширины окна виртуального анализатора по концентрации МВБЭ 47

Рисунок 5.9 – Значение коэффициента детерминации от ширины Окна виртуального анализатора по концентрации димеров При выборе оптимальной ширины окна виртуального анализатора по критерию коэффициента детерминации, решается задача оптимизации (минимизации) по критерию коэффициента детерминации, который будет иметь вид :     yk  yˆ k Bˆ h  k 1   1  R 2  k  M 1N 2   yk  y  N J R2 2 (5.28) k  M 1 На рисунках 5.10 и 5.11 представлены графики изменения значения критерия соответствия модели данным от ширины виртуального анализатора по концентрации МВБЭ и димеров соответственно. 48

Рисунок 5.10 – Значение критерия соответствия модели данным от ширины окна виртуального анализатора по концентрации МВБЭ Рисунок 5.11 – Значение критерия соответствия модели данным от ширины окна виртуального анализатора по концентрации димеров 49

Оптимальной ширине окна виртуального анализатора будет соответствовать ширина окна, при которой значение критерия соответствия модели данным будет минимальным. В таблице 5.4 представлены значения ширины окна виртуального анализатора по концентрации МВБЭ и димеров из критерия соответствия модели данным и значения среднеквадратичной ошибки и коэффициента детерминации адаптивного виртуального анализатора, соответствующие оптимальной ширине окна. Таблица 5.4 – Оптимальная ширина окна, среднеквадратичная ошибка и коэффициент детерминации по критерию соответствия модели данным Метил-втор-бутиловый Димеры эфир Оптимальная ширина окна 640 690 0.0012015 0.0018585 0.6388 0.12141 hR 2 Среднеквадратичная ошибка RMSE Коэффициент детерминации R2 5.3.2.4. Векторный критерий Изучив таблицы 5.2, 5.3 и 5.4, можно заметить, что, выбирая оптимальную ширину окна только по какому-то одному из вышеописанных критериев, будет получен виртуальный анализатор, у которого ошибка прогнозирования будет удовлетворять требованиям, а коэффициент детерминации не будет удовлетворять, и наоборот. Поэтому для адаптивного виртуального анализатора предложен векторный критерий (свертка критериев) [21]: ~ P J   Kp  Jp , (5.29) p 1 50

 где J  J RMSE J R2  J Bˆ , h - матрица значений вышеописанных критериев от  ширины окна виртуального анализатора, K  K RMSE KR2  K Bˆ , h - вектор весовых коэффициентов критериев, показывающих, какой из критериев наиболее важен при выборе оптимальной ширины окна виртуального анализатора. На рисунках 5.12 и 5.13 представлены графики изменения значения векторного критерия от ширины виртуального анализатора по концентрации МВБЭ и димеров соответственно, при условии, когда критерий ошибки прогноза, критерий коэффициента детерминации и критерий изменчивости параметров в одинаковой мере учитываются при выборе оптимальной ширины окна виртуального анализатора. В таблице 5.5 представлены значения ширины окна виртуального анализатора по концентрации МВБЭ и димеров из векторного критерия и значения среднеквадратичной ошибки и коэффициента детерминации адаптивного виртуального анализатора, соответствующие оптимальной ширине окна. Рисунок 5.12 – Значение векторного критерия от ширины окна виртуального анализатора по концентрации МВБЭ 51

Рисунок 5.13 – Значение векторного критерия от ширины окна виртуального анализатора по концентрации димеров Таблица 5.5 – Оптимальная ширина окна, среднеквадратичная ошибка и коэффициент детерминации по векторному критерию Метил-втор-бутиловый Димеры эфир Оптимальная ширина окна 750 510 0,00023305 0.0014992 0.53215 0.090857 hJ~ Среднеквадратичная ошибка RMSE Коэффициент детерминации R2 5.3.2.5. Сравнение критериев В таблицах 5.6 и 5.7 представлено сравнение среднеквадратичной ошибки и коэффициента детерминации для ширины окна виртуальных 52

анализаторов по концентрации МВБЭ и димеров, используемой на практике и выбранной по критериям оптимальности и векторному критерию. Таблица 5.6 – Сравнение оптимальной ширины окна, среднеквадратичной ошибки и коэффициента детерминации виртуального анализатора по концентрации МВБЭ Ширина Применяе окна тся на практике, По По По ошибки критерию критерию векторно прогнозирован коэффициен изменчивос h0 Значени По критерию ия hRMSE му та ти критерию детерминац параметров , hJ~ ии, hR 2 , hBˆ , h 900 430 890 640 0.0001929 0.000025187 0.0012015 750 е RMSE 4 R2 0.23001 0.42087 0.6388 0.00005191 0,0002330 5 5 0.41925 0.53215 Таблица 5.7 – Сравнение оптимальной ширины окна, среднеквадратичной ошибки и коэффициента детерминации виртуального анализатора по концентрации димеров Ширина Применяетс По критерию окна Значени я на ошибки практике, прогнозиров h0 ания, hRMSE По По По критерию критерию векторно коэффициен изменчивос му та ти критерию детерминац параметров , hJ~ ии, hR 2 , hBˆ , h 120 220 690 720 510 0.0014359 0.0011715 0.0018585 0.0018287 0.0014992 е RMSE 53

R2 0.0077663 0.0035104 0.12141 0.11112 0.090857 Из таблиц 5.6 и 5.7 видно, что векторный критерий выбора оптимального значения ширины окна выбирает значение ширины окна виртуального анализатора, при котором достигается компромисс между ошибкой прогнозирования и мерой соответствия регрессионной модели данным. 5.3.3. Сравнение адаптивного и неадаптивного ВА Разрабатываемый виртуальный анализатор является адаптивным из-за того, что во время технологического процесса ректификации могут измениться параметры колонны. То есть модель неадаптивного ВА, разработанного в начале технологического процесса может не удовлетворять данным через какой-то промежуток времени. На Рисунках 5.14 и 5.15 представлено сравнение значений десятичного логарифма от MSE для метил-втор-бутилового эфира и димеров соответственно. По графикам Рисунков 5.14 и 5.15 видно, порядок MSE для неадаптивного ВА больше по сравнению с порядком MSE для адаптивного ВА, то есть адаптивный виртуальный анализатор более точно предскажет значение выходной переменной при определенных значениях входных переменных. 54

Рисунок 5.14 – Сравнение MSE адаптивного и неадаптивного виртуального анализатора по концентрации метил-втор-бутилового эфира. Рисунок 5.15 – Сравнение MSE адаптивного и неадаптивного виртуального анализатора по концентрации димеров. 55

6. ИСПОЛЬЗОВАНИЕ ВИРТУАЛЬНОГО АНАЛИЗАТОРА В СИСТЕМЕ УПРАВЛЕНИЯ ТЕХНОЛОГИЧЕСКИМ ПРОЦЕССОМ РЕКТИФИКАЦИИ Виртуальный анализатор по концентрации компонента в выходном продукте является одной из основных частей системы управления. Для управления технологическим объектом в режиме близком к режиму реального времени необходимо иметь данные о состоянии объекта и качестве выходного продукта, интервал получения которых минимален, например, 1 минута. В связи с тем, что невозможно получить результаты лабораторных исследования так часто, используется виртуальный анализатор, разработанный в предыдущей главе, для прогноза значений концентраций метил-вторбутилового эфира и димеров в выходном продукте по данным полученным при опросе датчиков объекта с интервалом в 1 минуту. Разрабатываемая система управления процессом ректификации предназначена для поддержания концентрации метил-втор-бутилового эфира или димеров в выходном продукте в заданном диапазоне. На рисунке 6.1 представлена схема колонны и системы ПИД-регулятора. Получая данные с датчиков, виртуальный анализатор прогнозирует значение концентрации МВБЭ (димеров). На основе данного значения, входящего в блок ПИД-регулятора “CAI_MSBE” (“CAI_DIME”) в качестве переменной процесса, и задания “SV_MSBE” (“SV_DIME”) ПИД-регулятор выдает задание для ПИД-регулятора расхода метанола на синтез “PID FIRC3654_2”. 56

Рисунок 6.1 – Схема колонны для ПИД регулирования 6.1. Модель динамики технологического объекта (PDS HONEYWELL определение ПФ) Для управления процессом необходимо знать динамику изменения концентрации МВБЭ (димеров) в зависимости от расхода метанола на синтез. 57

Динамику исследуемого процесса можно описать передаточной функцией первого порядка с запаздыванием: Gs   K p e s 1  s Для определения (6.1) параметров передаточной функции можно использовать программу Honeywell Profit Design Studio (лицензия на использование представлена в Приложении). Для вычисления параметров передаточной функции необходимо загрузить данные о состоянии объекта и рычагов управления, а также данные, показывающие качество выходного продукта. Загрузив данные о состоянии системы и о качестве выходного продукта в модель, используя функцию “Import From Excel” (рисунок 6.2) получим графики изменения параметров от времени (рисунок 6.3). Рисунок 6.2 – Чтение минутных данных из файла Excel 58

Рисунок 6.3 – Графики изменения значений показателей технологического объекта во времени Из всего набора полученных графиков изменения значений параметров во времени выберем те, которые входят в разрабатываемый контур управления: FIRC3654_2 (задание по расходу метанола на синтез), VA_MSBE (спрогнозированные значения концентрации метил-втор-бутилового эфира в выходном продукте), LAB_MTBE_MSBE измерений концентрации (результаты метил-втор-бутилового эфира лабораторных в выходном продукте), VA_DIME (спрогнозированные значения концентрации димеров в выходном продукте), LAB_MTBE_DIME (результаты лабораторных измерений концентрации димеров в выходном продукте). Также на изменение концентраций метил-втор-бутилового эфира и димеров в выходном продукте влияют: FIRC3656 (расход метанола на форконтакт) и FIRC3651 (загрузка блока МТБЭ). 59

Выведем графики изменения вышеописанных величин отдельно (рисунок 6.4 и 6.5). Рисунок 6.4 – Изменение LAB_MTBE_MSBE, VA_MSBE, FIRC3656, FIRC3654_2 и FIRC3651 во времени 60

Рисунок 6.5 – Изменение LAB_MTBE_DIME, VA_DIME, FIRC3656, FIRC3654_2 и FIRC3651 во времени С помощью команды “Identify -> Set Overall Options …” задаем количество испытаний “Number of Trials” и параметры конечной импульсной характеристики “FIR” (рисунок 6.6). 61

Рисунок 6.6 – Задание параметров для построения конечной импульсной характеристики Затем, с помощью команды “Identify -> Fit FIR/PEM/CLid Models...” выбираем входную и выходную величину передаточной функции (рисунок 6.7 и 6.8) и интервал времени, на основании которого будет вычисляться параметры передаточной функции (рисунок 6.9 и 6.10). При выборе интервала времени следует учесть то, что нам необходимо выявить передаточную функцию изменения концентраций метил-втор-бутилового эфира и димеров в выходном продукте от FIRC3654_2 (задание по расходу метанола на синтез). То есть следует выбрать интервал, на котором FIRC3654_2 изменяется, а FIRC3656 и FIRC3651 постоянны. 62

Рисунок 6.7 – Выбор входной и выходной переменной для вычисления параметров передаточной функции для МВБЭ Рисунок 6.8 – Выбор входной и выходной переменной для вычисления параметров передаточной функции для димеров 63

Рисунок 6.9 – Выбор интервала времени для вычисления передаточных функций для МВБЭ Рисунок 6.10 – Выбор интервала времени для вычисления передаточных функций для димеров После выбора интервала времени и входной и выходной величин, нажатием на “Fit FIR”, строится конечная импульсная характеристика. Конечные импульсные характеристики для МВБЭ и димеров от FIRC3654_2 представлены на рисунке 6.11. 64

Затем, с помощью команды “Identify -> Fit Parametric Models...” вычисляем параметры передаточной функции. Выбрав дискретный метод, устанавливаем индивидуальные параметры (рисунок 6.12). Затем, с помощью кнопок “Previous CV”, “Next CV”, “Prev. MV/DV” и “Next MV/DV” выбираем ячейку, в которой находится построенная ранее модель FIR (рисунок 6.13). Выбрав номер испытания (Trial) и нажав на “Options”, переходим к выбору порядка передаточной функции. В появившемся окне выставляем значение “Order” равным 1 и нажимаем “Do It”. Рисунок 6.11 – Конечные импульсные характеристики для МВБЭ и димеров 65

Рисунок 6.12 – Выбор метода построения параметрической модели Рисунок 6.13 – Выбор ячеек с полученными FIR моделями 66

После вышеописанных действий имеем модель FIR и параметрическую модель зависимости МВБЭ и димеров от FIRC3654_2 представленную на рисунке 6.14. Рисунок 6.14 – Построение параметрических моделей Затем, с помощью команды “Identify -> Select Final Models...” строятся модели и выписывается передаточные функции. Рисунок 6.15 – Конечные модели и передаточные функции 67

Таким образом с использование программы Honeywell Profit Design Studio были получены передаточные функции для МВБЭ (6.2) и димеров (6.3) от FIRC3654_2:  0.0951e 155s GMSBE s   1  20.9s (6.2)  0.0534e 122 s GDIME s   1  50.1s (6.3) При использовании данных передаточных функций стоит учесть, что переменные  и  имеют размерность минут. Это необходимо для использования в других программных продуктах, где данные параметры передаточных функций имеют размерность секунд. 6.2. Синтез ПИД регуляторов (CAI_MSBE, CAI_DIME) на основе виртуальных анализаторов Зная передаточную функцию процесса можно приступать к разработке системы управления. На практике используют ПИД-регулятор, так как он обеспечивает выполнение наиболее общих управляющих функций с использованием пропорционального, интегрального и дифференциального регулирования по отклонению переменной процесса (PV) от значения задания (SV). На практике при управлении процессом ректификации используют ПИД-регуляторы компании Yokogawa. Расчет ПИД управления представляет собой основу обработки ПИД управляющих воздействий, выполняющую вычисление изменения управляющего выхода ( MV ), используя алгоритмы ПИД управления [22]. ПИД регулятор использует следующие пять алгоритмов ПИД управления для выполнения расчета ПИД управления. Используемый вариант зависит от характеристик управляемой системы и цепи управления.  Базовое ПИД управление (PID) 68

 ПИД управление, пропорциональное PV и по производной PV (IPD)  ПИД управление по производной PV (PI-D)  Автоматическое определение  Автоматическое определение 2 На рисунке 6.15 приведена блок схема расчета ПИД управления. Рисунок 6.15. Блок-схема расчета ПИД управления При расчете ПИД управления входные параметры зависят от алгоритма ПИД управления. В таблице 6.1 перечислены алгоритмы расчета ПИД управления и входные параметры: Таблица 6.1. Алгоритмы расчета ПИД управления и входные параметры Алгоритм ПИД Входные переменные полинома управления Пропорциональны Дифференциальны Интегральны й член й член й член PID En En En I-PD PV PV En PI-D En PV En Автоматическо Аналогично I-PD в режиме AUT е определение Аналогично PI-D в режиме CAS или RCAS Автоматическо Аналогично I-PD в режиме AUT или RCAS е определение 2 Аналогично PI-D в режиме CAS 69

6.2.1. Базовое ПИД управление (PID) Базовое ПИД управление предполагает выполнение управляющих действий – пропорциональных, интегральных и дифференциальных, отслеживающих изменения значения задания. Данный алгоритм используется, когда константа времени процесса обработки велика, а управление предполагает немедленную реакциюна изменение значения задания. Расчетное выражение базового ПИД управления (PID):   T T MVn  K P  K S En  En  D En , TI T   (6.4) где MVn - приращение управляющего воздействия; En - отклонение En  PVn  SVn ; PVn - переменная процесса; SVn - значение задания; En - приращение отклонения En  En  En 1 ; T - период управления; K P  KS  100 , PB MSH  MSL , MSH – верхний предел шкалы MV, MSL – нижний предел SH  SL шкалы MV, SH – верхний предел шкалы PV, SL – нижний предел шкалы PV PB - диапазон пропорциональности (%); TI - время интегрирования; TD время дифференцирования. Индексы “n” и “n-1” означают номер выборки в ходе периода управления. Дифференциальное уравнение для расчета ПИД управления вычисляет изменение управляющего воздействия (приращение). Новое значение выхода получается путем прибавления текущего приращения управляющего воздействия ( MVn ) к предыдущему значению управляющего выхода ( MVn 1 ). Переменная процесса (PV) и значение задания (SV), используемые в расчете, задаются в физических единицах. Приращение управляющего воздействия 70

MV  , получаемое в физических единицах, преобразуется с использованием коэффициента преобразования шкалы K S  . 6.2.2. Пропорциональный (PV) и дифференциальный алгоритм ПИД управления (I-PD) Данный алгоритм расчета ПИД управления отличается от базового тем, что при изменении значения задания предполагается только интегральные действия. Данный алгоритм обеспечивает устойчивое управление даже при резком изменении значения задания (SV) при его установке путем ввода численного значения. Одновременно данный алгоритм обеспечивает адекватную реакцию на характерные изменения управляемого процесса, колебания нагрузки и возмущения с использованием пропорционального, дифференциального и интегрального управляющих действий соответственно. Расчетное выражение для пропорционального и дифференциального алгоритма ПИД управления (I-PD):   T T MVn  K P  K S PVn  En  D PVn , TI T   (6.5) где PVn  PVn  PVn1. 6.2.3. Дифференциальный алгоритм ПИД управления (PI-D) В отличии от базового данный алгоритм расчета ПИД управления предполагает только пропорциональное и интегральное действие при изменении значения задания (без действия по производной). Данный алгоритм используется, когда необходимо лучшее отслеживание изменений значения задания, например, при наличии вторичного блока управления в каскадном контуре управления. Расчетное выражение дифференциального алгоритма ПИД управления: 71

  T T MVn  K P  K S En  En  D PVn , TI T   (6.6) 6.2.4. Автоматическое определение При работе ПИД регулятора в каскадном (CAS) или внешнем каскадном (RCAS) режиме для выполнения расчетов используется алгоритм ПИД управления по производной PV (PI-D), что обеспечивает лучшее отслеживание изменения значения задания. При работе блока в автоматическом режиме (AUT) для выполнения расчетов используется пропорциональный и дифференциальный алгоритм ПИД управления (I-PD), что обеспечивает устойчивое управление даже в случае резкого изменения значения задания, обусловленного вводом численного значения. 6.2.5. Автоматическое определение 2 При работе ПИД регулятора в каскадном (CAS) для выполнения расчетов используется алгоритм ПИД управления по производной PV (PI-D). При работе блока во внешнем каскадном (RCAS) или автоматическом (AUT) режиме для выполнения расчетов используется пропорциональный и дифференциальный алгоритм ПИД управления (I-PD). В каскадном режиме (CAS) автоматическое определение 2 предполагает возможность отслеживания изменения значения задания (CSV). Во внешнем каскадном режиме (RCAS) данный алгоритм предотвращает резкое изменение выходного значения, обусловленное резким изменением внешнего значения задания (RSV). На практике в основном используется пропорциональный (PV) и дифференциальный алгоритм ПИД управления (I-PD). На рисунке 6.16 представлена модель ПИД-регулятора реализованная в среде MATLAB Simulink на основе формулы (6.5). 72

Рисунок 6.16 – Модель ПИД регулятора 6.3. Настройка ПИД-регуляторов Известно множество методов настройки ПИД-регуляторов. В данной работе будут рассмотрены метод Кохен-Куна, метод на основе целевой функции, метод IMC и метод оптимизации. 6.3.1. Метод Кохен-Куна Известен метод настройки ПИД регулятора, который был разработан Cohen и Coon [23]. Метод использует четвертичный коэффициент упадка, согласно которому отклик PV на изменение задания SP должен иметь небольшие колебания с высотой второго пика равной четверти высоты первого, и представляет собой набор формул, основанный на параметрах передаточной функции процесса. Для П-регулятора: Kc  1   1 K P   3  (6.7) Для ПИ-регулятора: 3   30   1      Kc  0 . 9  T     i K P   12   9  20     (6.8) Для ПИД-регулятора: 73

6     32      1  4   4  Kc   T T         i d 8 2   K P   3 4   13  11          (6.9) При настройке ПИД-регуляторов управления концентрациями МВБЭ и димеров в выходном продукте методом Cohen-Coon были получены коэффициенты, представленные в таблице 6.2. На рисунках 6.17 и 6.18 представлены графики переходных процессов изменения концентраций метил-втор-бутилового эфира и димеров соответственно. Таблица 6.2. Коэффициенты ПИД-регулятора полученные методом CohenCoon Коэффициент МВБЭ Димеры Kc -4,519 -14,94 Ti 9836 10500 Td 1440 1845 Рисунок 6.17 – Переходный процесс изменения концентрации МВБЭ в выходном продукте. Коэффициенты ПИД-регулятора получены методом Cohen-Coon 74

Рисунок 6.18 – Переходный процесс изменения концентрации димеров в выходном продукте. Коэффициенты ПИД-регулятора получены методом Cohen-Coon 6.3.2. Метод на основе целевой функции Существует метод настройки ПИД-регулятора, использующий штрафные функции [24]. Суть метода настройки ПИД регулятора заключается в минимизации значения целевой функции. Есть три основных настроичных критерия используемые Lopez, Miller, Smith и Murrill: интеграл абсолютного отклонения (IAE) (6.10), интеграл квадрата отклонения (ISE) (6.11), интеграл от времени абсолютного отклонения (ITAE) (6.12). Также существует еще один критерий, неиспользуемый Smith и Murrill, интеграл во времени квадрата отклонения (ITSE) (6.13). Каждый из этих критериев является целевой функцией, представляющей размер и продолжительность отклонения.  IAE   E .dt (6.10) 0  ISE   E 2 .dt (6.11) 0 75

 (6.12) ITAE   E t.dt 0  (6.13) ITSE   E 2 t.dt 0 При вычислении значения целевой функции мы не учитываем знак отклонения из-за того, что при устойчивом колебании с учетом знака отклонения значение целевой функции будет стремиться к нулю. Минимизация ISE эквивалентна минимизации дисперсии, а также стандартному отклонению ошибки регулятора. Но в результате настройки ПИД регулятора минимизацией ISE или IAE можно получить регулятор, который устраняет как можно быстрее большую ошибку, возникающую сразу после истечения времени ожидания, за счет чего возникает некоторое колебательное поведение в течение некоторого времени после нарушения. Добавление времени в ITAE и ITSE предоставляет такой весовой фактор, что малое отклонение в течении длительного времени после нарушения вносит такой же вклад в целевую функцию, что большое отклонение в начале нарушения. Поскольку абсолютное значение ошибки не превышает 1, то возведение в квадрат малой ошибки дает значение штрафа близкое к нулю, тем самым подрывает преимущества времени, как весового фактора. Исходя из данных замечаний на практике принято использовать ITAE в качестве целевой функции. В таблице 6.3 представлены коэффициенты, используемые при вычислении коэффициентов регулятора методом на основе целевой функции, разработанные Lopez, Miller, Smith и Murrill, для регуляторов различной структуры и алгоритма управления. Таблица 6.3 – Коэффициенты, используемые при вычислении коэффициентов регулятора   A   Kc    KP B Ti     A      Td  A     B B 76

A 0.902 0.984 1.435 1.411 1.305 1.495 0.490 0.859 1.357 B -0.985 -0.986 -0.921 -0.917 -0.959 -0.945 -1.084 -0.977 -0.947 A B A B IAE, P load PI 0.608 -0.707 change PID 0.878 -0.749 0.482 1.137 ISE, P load PI 0.492 -0.739 change PID 1.101 -0.771 0.560 1.006 ITAE, P load PI 0.674 -0.680 change PID 0.842 -0.738 0.381 0.995 IAE, P SP PI 0.758 -0.861 1.020 -0.323 change PID 1.086 -0.869 0.740 -0.130 0.348 0.914 ITAE, P SP PI 0.586 -0.916 1.030 -0.165 change PID 0.965 -0.855 0.796 -0.147 0.308 0.929 Вычислим коэффициенты ПИД регуляторов концентраций метил-вторбутилового эфира и димеров в выходном продукте. В качестве целевой функции будем использовать ITAE. Уравнения для расчета коэффициентов ПИД регулятора будут иметь вид:      Ti  Td  0.308   0.965  0.147 (6.14)       0.796  Kc  KP   При настройке ПИД-регуляторов управления концентрациями МВБЭ и 0.855 0.929 димеров в выходном продукте методом на основе целевой функции были получены коэффициенты, представленные в таблице 6.4. На рисунках 6.19 и 6.20 представлены графики переходных процессов изменения концентраций метил-втор-бутилового эфира и димеров соответственно. Таблица 6.4 – Коэффициенты ПИД-регулятора полученные методом на основе целевой функции Коэффициент Kc Ti Td МВБЭ -1,83 2115 2485 Димеры -8,443 4304 2117 77

Рисунок 6.19. Переходный процесс изменения концентрации МВБЭ в выходном продукте. Коэффициенты ПИД-регулятора получены методом настройке на основе целевой функции Рисунок 6.20. Переходный процесс изменения концентрации димеров в выходном продукте. Коэффициенты ПИД-регулятора получены методом настройке на основе целевой функции 78

6.3.3. Метод на основе внутренней модели(IMC) IMC (internal model control) tuning [24] [25] – настройка внутреннего регулятора модели. Метод разработан с использованием техники, известной как прямой синтез. Прямой синтез - это метод, принцип которого заключается в синтезировании регулятора, который будет отвечать на изменение SP по определенной траектории. Траектория определяется как приближение первого порядка к окончательному значению с пользовательской постоянной времени  . Низкое значение  даст более агрессивное управление, чем робастное. При высоком – более робастное. Прямой синтез может быть применен для процессов высокого порядка и для всех типов регуляторов. Однако, результат может не иметь форму ПИД алгоритма, поэтому при настройке ПИД регулятора требуется сделать некоторые приближения. Например, пренебречь составляющими высокого порядка. Полученный результат приводят к небольшим различиям между полученными формулами настройки. Значение  обычно подбирается методом “проб и ошибок” до получения разрешенного значения MV перерегулирования. Начальной точкой при подборе  можно принять:    (6.15) В таблице 6.5 представлены наиболее используемые формулы, но в литературе можно встретить и другие версии. Таблица 6.5. Уравнения для определения коэффициентов ПИД-регулятора методом IMC с помощью параметров передаточной функции процесса Self-regulating Ti Td Kc PID (noninteractive)  1 2 KP      2  2   Integrating Ti Kc 2   2   1 2 KP      2  Td   2 4 2   79

PID (interactive) 1  KP     2   2    1 2 2 2 2 KP      2  PI 2    1  1 2   KP   K P    2 При настройке ПИД-регуляторов управления концентрациями МВБЭ и 2 димеров в выходном продукте методом IMC были получены коэффициенты, представленные в таблице 6.6. На рисунках 6.21 и 6.22 представлены графики переходных процессов изменения концентраций метил-втор-бутилового эфира и димеров соответственно. Таблица 6.6 – Коэффициенты ПИД-регулятора полученные методом IMC Коэффициент Kc Ti Td МВБЭ -3,127 5904 987,7 Димеры -7,074 6666 1650 Рисунок 6.21. Переходный процесс изменения концентрации МВБЭ в выходном продукте. Коэффициенты ПИД-регулятора получены методом IMC 80

Рисунок 6.22. Переходный процесс изменения концентрации димеров в выходном продукте. Коэффициенты ПИД-регулятора получены методом IMC 6.3.4. Метод оптимизации Для настройки ПИД регулятора можно использовать методы оптимизации. Используя блок “Check Step Response Characteristics” раздела “Simulink Design Optimization” программы “MATLAB”, можно подобрать коэффициенты ПИД-регулятора. Для этого на вход блока подается значение настраиваемой функции. Далее необходимо выбрать ограничения переходного процесса (рисунок 6.23). Устанавливают интервал времени, в который процесс должен достичь определенного уровня от значения задания, момент времени, в котором процесс достиг значения задания, допустимое значение перерегулирования и недорегулирования(недолета) выходной величины. Для химической промышленности нежелательно наличие перерегулирования. 81

Рисунок 6.23 – Задание ограничений процесса Затем необходимо выбрать из списка всех переменных участвующих в модели те, значение которых будет оптимизироваться (рисунок 6.24). Рисунок 6.24 – Выбор оптимизируемых параметров 82

Затем при нажатии кнопки “Optimize”, запускается процесс вычисления значений коэффициентов ПИД-регулятора методом оптимизации (рисунок 6.25). В окне “Response Optimization” показаны графики переходных процессов и установленные ограничения. В окне “Optimization Progress Report” представлены такие сведения процесса оптимизации, как значение, показывающего, на сколько переходный процесс вышел за рамки ограничений. Рисунок 6.25 – Процесс оптимизации После завершения процесса мы получаем переменные, значение которых получено с помощь метода оптимизации. При настройке ПИДрегуляторов управления концентрациями МВБЭ и димеров в выходном продукте методом оптимизации были получены коэффициенты, представленные в таблице 6.7. На рисунках 6.26 и 6.27 представлены графики переходных процессов изменения концентраций метил-втор-бутилового эфира и димеров соответственно. Таблица 6.7 – Коэффициенты ПИД-регулятора полученные методом оптимизации Коэффициент Kc Ti Td МВБЭ -0.4785 1128 1000 Димеры -3.541 3489 1000 83

Рисунок 6.26 – Переходный процесс изменения концентрации МВБЭ в выходном продукте. Коэффициенты ПИД-регулятора получены методом оптимизации Рисунок 6.27 – Переходный процесс изменения концентрации димеров в выходном продукте. Коэффициенты ПИД-регулятора получены методом оптимизации 84

6.3.5. Сравнение методов Что бы из пропорционального коэффициента K c перейти к значению PB, используемого при настройке ПИД-регуляторов компании Yokogawa, воспользуемся следующей формулой [22]: PB  (6.16) 100 Kc В таблице 6.8 и 6.9 представлены параметры, используемые ПИД регуляторами компании Yokogawa, полученные вышеописанными методами настройки, для ПИД регуляторов управляющих изменением концентраций метил-втор-бутилового эфира и димеров в выходном продукте соответственно. Также в таблице 6.8 и 6.9 приведена оценка качества настройки ПИДрегулятора с использованием критерий интеграла квадрата отклонения (6.11). Можно заметить, что оценка качества настройки ПИД-регулятора минимальна для метода на основе целевой функции, но переходный процесс при данном методе настройки имеет перерегулирование, что нежелательно для химических процессов. Остальные методы зависят от параметров настройки (  в методе IMC) и от ограничений наложенных на переходный процесс (метод оптимизации), которые подбираются индивидуально для каждого процесса. Подобрав определенные значения параметров настройки и ограничений можно добиться такого качества настройки ПИД-регулятора, что критерий ISE будет меньше чем при настройке ПИД-регулятора на основе целевой функции. Таблица 6.8 – Параметры настройки ПИД-регулятора компании Yokogawa, управляющего изменением концентрации МВБЭ Метод PB TI TD ISE Cohen-Coon -22.13 9836 1440 21130 Метод на основе -54.64 2115 2485 14950 целевой функции 85

IMC -31.98 5904 987.7 18270 Метод -208.98 1128 1000 19410 оптимизации Таблица 6.9 – Параметры настройки ПИД-регулятора компании Yokogawa, управляющего изменением концентрации димеров Метод PB TI TD ISE Cohen-Coon -6.69 10500 1845 16110 Метод на основе -11.84 4304 2117 13040 IMC -14.136 6666 1650 17210 Метод -28.24 3489 1000 16730 целевой функции оптимизации На рисунке 6.28 и 6.29 представлено сравнение переходных процессов изменения концентраций метил-втор-бутилового эфира и димеров при настройке ПИД-регуляторов различными методами. Рисунок 6.28 – Переходные процессы изменения концентрации МВБЭ при настройке ПИД-регулятора различными методами 86

Рисунок 6.29 – Переходные процессы изменения концентрации димеров при настройке ПИД-регулятора различными методами Из таблиц 8.8, 8.9 и рисунков 6.28, 6.29 видно, что наименьшее значение критерия интеграла квадрата отклонения имеет ПИД-регулятор настроенный методом на основе целевой функции, но переходный процесс имеет перерегулирование, наличие которого неприемлемо в химической промышленности. Из методов настройки, представленных выше, при которых переходный процесс не имеет перерегулирования, наиболее быстро достигается заданное значение при настройке ПИД-регулятора методом оптимизации. ЗАКЛЮЧЕНИЕ В ходе данной работы рассмотрена ректификационная колонна производства МТБЭ в качестве объекта управления. Приведена технологическая схема, разработана упрощенная математическая модель объекта управления. Проведено исследование наблюдаемости технологического объекта в зависимости от его параметров для подтверждения необходимости использования виртуальных анализаторов в системе управления. Выяснено, 87

что объект не полностью наблюдаемый. Количество наблюдаемых компонент вектора состояния различно в зависимости от значений флегмового числа и количества отбираемого верхнего продукта, а также номера ступени разделения, на которую подается сырье. Для определенных значений параметров ТО выделена область значений флегмового числа и количества отбираемого верхнего продукта, рекомендуемых к использованию на производстве. Были отобраны наиболее информативные входные переменные для виртуальных анализаторов по концентрации метил-втор-бутилового эфира и димеров, построена модель виртуального анализатора методом робастной регрессии (взвешенного МНК). Изучен принцип работы адаптивного виртуального анализатора по методу «движущегося окна» и применен к технологическому объекту ректификации. Показано, что значение MSE адаптивного виртуального анализатора является 6  10 9 % значения MSE неадаптивного виртуального анализатора. В рамках данной работы были исследованы различные подходы к выбору оптимальной ширины окна виртуального анализатора. Разработаны критерий ошибки прогнозирования, критерий изменчивости параметров модели, критерий соответствия модели данным, векторный критерий. С использованием программного продукта PDS HONEYWELL определены передаточные функции изменения концентраций МТБЭ и димеров в выходном продукте от изменения расхода метанола на синтез. Изучены алгоритмы ПИД управления, используемые ПИД- регуляторами компании Yokogawa. В среде MATLAB Simulink построена модель ПИД-регулятора с пропорциональным PV и дифференциальным алгоритмом ПИД управления. На базе полученной модели были изучены такие методы настройки ПИД-регуляторов как метод Кохен-Куна, метод на основе целевой функции, метод IMC и метод оптимизации. Было проведено их сравнение. 88

СПИСОК ИСПОЛЬЗОВАННЫХ ИСТОЧНИКОВ 1. А.А. М. Виртуальные анализаторы: концепция построения и применения в задачах управления непрерывными ТО // Автоматизация в промышленности, август 2003. pp. 28-44. 2. Шувалов В.В., Огаджанов Г.А., Голубятников В.А. Автоматизация производственных процессов в химической промышленности. Москва: Химия, 1991. 439-442 pp. 3. Мустафин А.И. К.С.Н..Е.А.Н. способ автоматического управления процессом ректификации, 2092222, Октябрь 10, 1997. 4. Шевчук В.П. А.С.А..Ф.С.О. Устройство автоматического регулирования процессом ректификации, 2449827, Май 10, 2012. 5. Поляк Б.Т., Щербаков П.С. Робастная устойчивость и управление. Москва: Наука, 2002. 303 pp. 6. Целая часть числа // Академик. URL: https://dic.academic.ru/dic.nsf/ruwiki/ 1187086 7. Филиппов В.В. Ректификация бинарной смеси. Метод. указ. к лабораторной работе. Самара: Самар. гос. тех. ун-т, 2016. 11 pp. 8. Определение ранга матрицы. Вычисление ранга матрицы по определению. // Math1.ru. URL: http://math1.ru/education/matrix/rang.html 9. Qi J., Sun K., Kang W. Optimal PMU placement for power system dynamic state estimation by using empirical observability Gramian // IEEE Transactions on Power Systems, Vol. 4, No. 30, 2015. pp. 2041-2054. 10. Цей Р., Шумафов М.М. Число обусловленности матрицы как показатель устойчивости при решении прикладных задач // Труды ФОРА, No. 16, 2011. 89

11. Singh A., Hahn J. On the use of epirical Gramians for controllability and observability analysis // Proc. of American Control Conference. Portland. 2005. pp. 140-146. 12. Индекс множественной корреляции // Semestr.ru. URL: https:// math.semestr.ru/regress/regres3.php 13. Корреляция, корреляционная зависимость // Математическая статистика для психологов. URL: http://statpsy.ru/correlation/correlation/ 14. Hans C. Article Navigation // Biometrika, Vol. 4, No. 96, Dec 2009. pp. 835845. 15. Luiz Mountinho, Graeme Hutcheson. The Sage Dictionary of Quantitative Management Research. London: Sage, 2011. 224-227 pp. 16. John Fox S.W. Robust Regression, October 2013. pp. 1-16. 17. Petr Kadlec, Ratko Grbic, Bogdan Gabrys. Review of adaptation mechanisms for data-driven soft sensors // Computers and Chemical Engineering, No. 35, 2011. pp. 5-7. 18. Среднеквадратическая ошибка модели // StatSoft. URL: http://statistica.ru/ glossary/general/srednekvadraticheskaya-oshibka/ 19. Дисперсия, среднеквадратичное (стандартное) отклонение, коэффициент вариации // statanaliz.info. URL: https://statanaliz.info/metody/opisaniedannyx/11-dispersiya-standartnoe-otklonenie-koeffitsient-variatsii 20. Коэффициент детерминации // Machinelearning.ru. URL: http:// www.machinelearning.ru/wiki/index.php?title=Коэффициент_детерминации 90

21. Ногин В.Д. ЛИНЕЙНАЯ СВЕРТКА КРИТЕРИЕВ В МНОГОКРИТЕРИАЛЬНОЙ ОПТИМИЗАЦИИ // ИСКУССТВЕННЫЙ ИНТЕЛЛЕКТ И ПРИНЯТИЕ РЕШЕНИЙ, No. 4, 2014. pp. 73-82. 22. Справочное руководство CENTUM VP. Часть D. 1st ed. 2008. 57-67 pp. 23. W.K. Ho, O.P. Gan, E.B. Tay. Performance and gain and phase margins of wellknown PID tuning formulas // IEEE Transactions on Control Systems Technology, Vol. 4, No. 4, Jul 1996. pp. 473-477. 24. King M. Process Control A Practical Approach. John Wiley & Sons Ltd, 2011. 51-65 pp. 25. Daniel E. Rivera M.M.S.S. Internal model control: PID controller design // Industrial &Engineering Chemistry Process Design and Development, No. 25, 1986. pp. 252-265. 26. Шляпников А.М. Стряхилева М.Н. Смирнов В.А. СИНТЕЗ МЕТИЛТРЕТ-БУТилового компонентов эфира товарных и других бензинов // высокооктановых МИР эфирных НЕФТЕПРОДУКТОВ. ВЕСТНИК НЕФТЯНЫХ КОМПАНИЙ, No. 3, 2008. pp. 9-12. 27. Mohd Ali J., Hoang N. Ha, Hussain M.A., Dochain D. Review and classification of recent observers applied in chemical process systems // Comp Chem Eng., No. 76, 2015. pp. 27-41. 28. Лебедев С.К. Лекция 14. Математические модели в пространстве состояний // Кафедра ЭП и АПУ. URL: http://drive.ispu.ru/elib/lebedev/ 14.html 29. Lee J.H. Model predictive control: Review of the tree decades of development // International Journal of Control, Automation and Systems, Vol. 3, No. 9, 2011. pp. 415-424. 91

Приложение A Лицензия продукта Honeywell Profit Suit 92

Рецензии:

Авторизуйтесь, чтобы добавить рецензию

Рецензия от Олег Снегирев

Представлена рецензия на выпускную квалификационную работу – магистерскую диссертацию научным сотрудником Института автоматики и процессов управления ДВО РАН

около 4 лет назад

Отзывы:

Авторизуйтесь, чтобы оставить отзыв