Проведен анализ методов построения детерминированных и вероятностных моделей сложных технологических объектов для их дальнейшего применения при разработке и эксплуатации систем автоматизированного управления. Реализовано процедуры сбора и обработки экспериментальных данных для создания базы данных среды моделирования MATLAB Simulink. Определены основные технологические параметры необходимые для реализации модели (входные и выходные температура и давление газа, частота вращения турбины нагнетателя, перепад давления на конфузоре). Для расширения функциональных возможностей имитационной модели нагнетателя ГПА (газоперекачивающий агрегат) реализовано как стационарный режим работы, так и режимы переходных процессов (пуск и остановка ГПА). В качестве алгоритмического обеспечения применен метод группового учета аргументов для расчета коэффициентов аппроксимирующих полиномов. Проведен сравнительный анализ зависимостей основных технологических параметров нагнетателя ГПА на основе экспериментальных данных и расчетных значений при различных режимах работы. Предложены и апробированы процедуры преобразования компонентов имитационных моделей в программные функциональные блоки PLC на основании стандарта IEC 61131 для их интеграции на уровне аппаратных средств для систем управления технологическим оборудованием компрессорных цехов газотранспортной системы.
Ключевые слова: газоперекачивающий агрегат, система автоматизированного управления, аппроксимирующий полином, имитационная модель, режимы работы ГПА, стандарт IEC 61131, функциональные блоки PLC.
Введение
В настоящее время проведение экспериментальных исследований на объектах нефтегазовой отрасли связано со значительными сложностями обусловленными необходимостью оформления соответствующих согласований и режима экономии энергоресурсов (изменение режимов работы ГПА требует дополнительных затрат энергоресурсов), поэтому разработка алгоритмического обеспечения и создания инструментальных средств для исследований на основе имитационных моделей объектов и средств управления является актуальной научно-технической задачей.
Анализ исследований и публикаций
Сравнительный анализ методов построения детерминированных и вероятностных моделей сложных технологических объектов указывает на возможность и эффективность применения метода группового учета аргументов (МГУА) для решения задач моделирования сложных технологических объектов [1, с.3–4, 90–93; 2, с.64].
В математических моделях, построенных на основе МГУА, учитываются вероятностные взаимосвязи состояний и признаков состояний технологических объектов. Это позволяет заменить процесс изучения особенностей, например, с точки зрения заданного объекта моделирования, газодинамических процессов протекающих в ГПА, статистической обработкой технологической информации, собираемой системой автоматизированного управления (САУ) ГПА. Поэтому модели, построенные на основе МГУА, позволяют оценивать изменение технического состояния объекта управления используя, например, значение абсолютного отклонения между измеренным технологическим параметром и значением, рассчитанным на основе модели.
МГУА можно использовать для решения следующих задач [2, с.64; 3 с.1]:
- идентификации физических закономерностей;
- аппроксимации многомерных процессов;
- краткосрочного пошагового прогнозирования процессов и событий;
- долгосрочного пошагового прогнозирования процессов и событий;
- экстраполяции физических полей;
- нормативного векторного прогнозирования процессов.
Основной результат теории МГУА состоит в том, что при неточных зашумленных данных с объекта управления и коротких выборках, минимум критерия указывает на нефизическую модель (решающее правило), точность которой выше, а структура более простая в сравнении с структурой полной физической модели.
Выделениенерешенныхпроблем
Математическое моделирование сложных объектов имеет следующие отличительные особенности:
- высокая размерность вектора входных координат 
;
- наличие значительного количества внутренних источников случайных помех, статистические характеристики которых, как правило, неизвестны;
- неизученность механизмов, определяющих направление и особенности протекания процессов в объектах;
- сложность постановки экспериментов для получения значительного количества сигналов 
(
– одномерный выход объекта).
Указанные особенности усложняют построение неформальных математических моделей сложных объектов. Как следствие, статистические режимы их функционирования приходится описывать уравнениями вида:
(1)
где 
.
Формированиецели
Целью данной работы является разработка модели нагнетателя ГПА на основе метода группового учета аргументов, исследование режимов работы технологических объектов и интеграция результатов моделирования в программные функциональные блоки PLC системы управления технологическим оборудованием компрессорных цехов газотранспортной системы.
Результаты
Ниже приведены основные процедуры и алгоритм применения метода группового учета аргументов для построения математических моделей технологических параметров.
Алгоритм МГУА включает следующие процедуры (рис. 1):
- — создание вектора 
, в который заносятся n первых значений одного из выходных технологических параметров (рис. 1, блок 1);
- — создание матрицы 
размерностью n × m, в столбцах которой размещаются первые n значений m входных технологических параметров (рис. 1, блок 3);
- — последовательный выбор из матрицы всех возможных комбинаций номеров столбцов — переменные 
(рис. 1, блок 4);
- — нахождение значений коэффициентов аппроксимирующего полинома 
:
, (2)
где 
для каждой из возможных комбинаций номеров столбцов матрицы и вектора 
с помощью метода наименьших квадратов (МНК) (рис. 1, блок 5);
- вычисление, используя найденные наборы коэффициентов аппроксимирующих полиномов , их значений в каждой из точек k=1,2,...N:
(3)
где (рис. 1, блок 6);
- вычисления среднеквадратического отклонения (СКО) остатков-разностей между вектором , k=1,2,...N и соответствующим значением каждого аппроксимирующего полинома:
, (4)
где (рис. 1, блок 7);
- отбор полиномов с минимальной дисперсией, которые будут использованы для построения аппроксимирующих полиномов второго и выше уровней (каждый уровень соответствует номеру итерации) (рис. 1, блок 10) и переход на следующий шаг (рис. 1, блок 12);
- сравнение СКО на предыдущем и текущем шагах (рис. 1, блок 9).
Если СКО на текущем шаге превышает СКО на предыдущем, то из всех полученных полиномов выбирается полином с минимальной дисперсией (рис. 1, блок 11).
Учитывая значимость влияния входных и выходных параметров на адекватность модели нагнетателя ГПА, предлагается применить следующие технологические параметры.
Входные параметры:
- Рвх — давление на входе нагнетателя;
- Твх — температура на входе нагнетателя;
- N — обороты турбины нагнетателя.
Выходные параметры:
- Рвых — давление на выходе нагнетателя;
- Твых — температура на выходе нагнетателя;
- dP — перепад давления на конфузоре нагнетателя.
Для построения математической модели применены экспериментальные выборки значений определенных технологических параметров САУ ГПА, которые занесены в электронную таблицу Excel.
Данные выборки получены за период с 10:00 04.02.2014 до 15:00 08.02.2014, в различных режимах работы ГПА компрессорной станции (КС) «Бердичев» управления магистральных газопроводов «Киевтрансгаз» (рис. 2).
Рис. 1. Блок-схема алгоритма МГУА
Рис. 2. Выборки технологических параметров нагнетателя ГПА КС «Бердичев»
Для трансляции данных из MS Excel в среду MATLAB апробирована процедура на основании функций импорта данных:
(5)
где 
— переменная, в которую записывается результат;
— переменная, содержащая путь к файлу;
— переменная, содержащая название листа в электронной таблице Excel;
— переменная, содержащая диапазон полей таблицы, который необходимо считать.
В результате после выполнения заданных функций данные из таблицы Excel будут записаны в указанную переменную в соответствии с заданным листом и диапазоном [4, с.1].
Для решения задачи поиска коэффициентов аппроксимирующего полинома необходимо применить функцию:
(6)
где 
— переменная, в которую записываются рассчитанные коэффициенты;
— аппроксимирующая функция;
— начальное приближение для нахождения коэффициентов аппроксимирующей функции;
— значения аргументов аппроксимирующей функции;
— значения аппроксимирующей функции.
Данная функция предназначена для расчета коэффициентов аппроксимирующих полиномов, решая при этом следующую оптимизационную задачу (минимизация квадрата разницы между экспериментальными и расчетными значениями аппроксимирующего полинома):
(7)
где — вектор коэффициентов аппроксимирующего полинома (целевая величина);
— входные данные;
— выходные данные;
— аппроксимирующая функция.
Таким образом, функция (6) рассчитывает коэффициенты для аппроксимирующей функции с помощью МНК [5, с.1].
Ниже приведен фрагмент программы (m-файл MATLAB), реализующей расчет коэффициентов аппроксимирующих полиномов согласно МГУА.
N = 363480 % количество значений параметров
x1 = transpose(P_in_1s); % давление до нагнетателя
x2 = transpose(T_in_1s); % температура до нагнетателя
x3 = transpose(n_1s); % обороты турбины нагнетателя
y = transpose(P_out_1s); % давление после нагнетателя
z = transpose(T_out_1s); % температура после нагнетателя
d = transpose(delta_P_1s); % перепад давления на конфузоре
% Комбинации входных параметров
x12 = [x1(1:N) x2(1:N)];
x13 = [x1(1:N) x3(1:N)];
x23 = [x2(1:N) x3(1:N)];
% Определение коэффициентов аппроксимирующих полиномов для давления на выходе
a12 = lsqcurvefit(@myfun,zeros(1,6),x12,y(1:N))
a13 = lsqcurvefit(@myfun,zeros(1,6),x13,y(1:N))
a23 = lsqcurvefit(@myfun,zeros(1,6),x23,y(1:N))
% Вычисление значений аппроксимирующих полиномов для давления на выходе
y12 = myfun(a12, [x1 x2])
y13 = myfun(a13, [x1 x3])
y23 = myfun(a23, [x2 x3])
% Среднеквадратичные отклонения полиномов от эталонных значений
sigma_y12 = sqrt(sum((y — y12).^2)/length(y));
sigma_y13 = sqrt(sum((y — y13).^2)/length(y));
sigma_y23 = sqrt(sum((y — y23).^2)/length(y));
%------------2-я итерация---------------
% Определение коэффициентов аппроксимирующих полиномов для давления
a1 = lsqcurvefit(@myfun,zeros(1,6), [y12(1:N) y23(1:N)],y(1:N))
a2 = lsqcurvefit(@myfun,zeros(1,6), [y12(1:N) y13(1:N)],y(1:N))
a3 = lsqcurvefit(@myfun,zeros(1,6), [y13(1:N) y23(1:N)],y(1:N))
% Вычисление значений аппроксимирующих полиномов
y1 = myfun(a1, [y12 y23])
y2 = myfun(a2, [y12 y13])
y3 = myfun(a3, [y13 y23])
% Среднеквадратичные отклонения полиномов от эталонных значений
sigma_y1 = sqrt(sum((y — y1).^2)/length(y))
sigma_y2 = sqrt(sum((y — y2).^2)/length(y))
sigma_y3 = sqrt(sum((y — y3).^2)/length(y))
% Задание аппроксимирующего полинома
function F = myfun(a, x)
F = a(1) + a(2) * x(:,1) + a(3) * x(:,2) + a(4) * x(:,1).* x(:,2) + a(5) * x(:,1).^2 + a(6) * x(:,2).^2;
Полученные среднеквадратические отклонения полиномов на первой и второй итерациях по каждому выходному параметру приведены в табл. 1, 2, при этом комбинации входных параметров обозначены (PT — давление-температура, PN — давление-обороты, TN — температура-обороты).
Таблица 1
Значение СКО на первой итерации
|
Зависимость PT |
Зависимость PN |
Зависимость TN |
|
|
Pвых |
0.6785 |
0.2585 |
2.3947 |
|
Tвых |
1.1890 |
2.7310 |
0.5957 |
|
dP |
0.0384 |
0.0260 |
0.0295 |
Таблица 2
Значение СКО на второй итерации
|
Зависимость PT, TN |
Зависимость PT, PN |
Зависимость PN, TN |
|
|
Pвых |
0.5438 |
0.2565 |
0.2561 |
|
Tвых |
0.5478 |
1.0326 |
3.5174 |
|
dP |
0.0261 |
0.0240 |
0.0301 |
На рис. 3–5 изображены графические зависимости выходных параметров модели нагнетателя ГПА (экспериментальные и расчетные).
Рис. 3. Зависимости давления на выходе нагнетателя ГПА от времени (на основании экспериментальных данных и расчетная)
Рис. 4. Зависимости температуры на выходе нагнетателя ГПА от времени (на основании экспериментальных данных и расчетная)
Рис. 5. Зависимости перепада давления на конфузоре нагнетателя ГПА от времени (на основании экспериментальных данных и расчетная)
В табл. 3, 4 приведены значения рассчитанных коэффициентов аппроксимирующих полиномов первого и второго уровней.
Таблица 3
Значения коэффициентов аппроксимирующих полиномов первого уровня
|
Комбинации входных параметров |
a0 |
a1 |
a2 |
a3 |
a4 |
a5 |
|
|
Pвых |
PT PN TN |
3.5658 0.0551 3.5013 |
0.5954 0.9855 —0.0523 |
0.3079 —0.0080 0.0156 |
0.0218 0.0001 —0.0002 |
0.0042 0.0007 0.0149 |
‑0.0208 0.0000 0.0000 |
|
Tвых |
PT PN TN |
8.1046 3.7478 0.0649 |
‑0.6373 2.6938 1.0407 |
1.7469 —0.0254 —0.0045 |
0.0529 0.0005 0.0000 |
‑0.0008 —0.0560 —0.0005 |
‑0.0508 0.0000 0.0000 |
|
dP |
PT PN TN |
0.0465 0.0041 0.0169 |
0.0151 0.0018 0.0016 |
0.0056 0.0005 0.0004 |
0.0005 0.0000 0.0000 |
‑0.0004 —0.0001 —0.0001 |
‑0.0006 0.0000 0.0000 |
Таблица 4
Значения коэффициентов аппроксимирующих полиномов второго уровня
|
Комбинации входных параметров |
a0 |
a1 |
a2 |
a3 |
a4 |
a5 |
|
|
Pвых |
PT, TN PT, PN PN, TN |
2.3621 —0.0386 0.3073 |
1.1357 —0.1121 1.0436 |
‑0.3734 1.1195 —0.0651 |
0.0139 0.0006 0.0024 |
‑0.0087 0.0005 —0.0018 |
‑0.0010 —0.0013 —0.0003 |
|
Tвых |
PT, TN PT, PN PN, TN |
0.1057 4.9314 2.6068 |
0.6366 1.1236 —0.3798 |
0.3434 —0.7817 1.0993 |
0.0732 —0.0149 —0.0077 |
‑0.0414 0.0010 0.0120 |
‑0.0313 0.0288 0.0016 |
|
dP |
PT, TN PT, PN PN, TN |
0.0273 0.0009 0.0064 |
‑0.7374 0.4560 —0.9322 |
‑0.7181 0.3882 1.6602 |
3.2343 0.8081 13.7545 |
2.9355 —1.2146 —3.7856 |
1.8105 0.9161 —9.1155 |
На рис. 6 приведена Simulink-модель нагнетателя ГПА (блок «Subsystem»). В результате апробации модели получены результаты, позволяющие оценить взаимосвязь между входными и выходными параметрами и исследовать режимы работы ГПА путем применения различных тестовых сигналов, включая моделирование аварийных режимов работы.
Рис. 6. Simulink-модель нагнетателя ГПА (блок «Subsystem») и результаты симуляции
На рис. 7 приведена Simulink-модель нагнетателя ГПА (внутренняя структура блока «Subsystem»), которая включает входные и выходные порты сигналов (технологических параметров), блоки индикации и расчетный блок с реализацией алгоритма расчета выходных технологических параметров модели в соответствии с входными параметрами.
Рис. 7. Simulink-модель нагнетателя ГПА (внутренняя структура блока «Subsystem»)
Ниже приведен фрагмент программы, которая реализована в блоке «Расчет» (рис. 7).
function [Pvyh, Tvyh, dP] = F(Pvh, Tvh, n)
% Используется метод группового учета аргументов (МГУА)
% Pvh — Давление газа на входе в нагнетатель
% Tvh — Температура газа на входе в нагнетатель
% n — Обороты нагнетателя
% Давление
%----------------------------------------------
aPN = [0.055057216854257 0.985451836548602 -0.007992668494920...
0.000109646472252 0.000685139307610 0.000001410140111];
aTN = [3.501338767183844 -0.052267370075987 0.015585031034252...
—0.000195412507309 0.014873941618501 -0.000000678595129];
aPNTN = [0.307280744812314 1.043581754957738 -0.065098002123334...
0.002440332212594 -0.001792166499577 -0.000320084897135];
PvyhPN = aPN(1) + aPN(2)*Pvh + aPN(3)*n + aPN(4)*Pvh*n +...
+ aPN(5)*Pvh^2 + aPN(6)*n^2;
PvyhTN = aTN(1) + aTN(2)*Tvh + aTN(3)*n + aTN(4)*Tvh*n +...
+ aTN(5)*Tvh^2 + aTN(6)*n^2;
Pvyh = aPNTN(1) + aPNTN(2)*PvyhPN + aPNTN(3)*PvyhTN +...
+ aPNTN(4)*PvyhPN*PvyhTN + aPNTN(5)*PvyhPN^2 + aPNTN(6)*PvyhTN^2;
Для интеграции компонентов имитационных моделей (например, блок «Расчет») в аппаратные средства систем управления на базе PLC Simatic S7 [6, с.180–186] или другие, апробирован программный модуль MATLAB PLC Coder, что обеспечивает автоматическую генерацию программных функциональных блоков PLC на основе стандарта IEC 61131–3 [7, c.7–9].
Выводы
В результате выполнения проектно-исследовательских работ предложена и реализована модель для исследования и имитации режимов работы нагнетателя ГПА с целью использования результатов исследований при построении и внедрении современных систем автоматизированного управления технологическим оборудованием компрессорных цехов газотранспортной системы. Предлагается использование компонентов имитационных моделей путем их преобразования в программные функциональные блоки PLC на основании стандарта IEC 61131–3.
При построении модели, кроме реализации «стационарного» режима работы нагнетателя предусмотрены также переходные процессы («пуск» и «остановка» ГПА), что обеспечивает расширенные функциональные возможности и большее соответствие модели технологическому объекту.
В результате апробации и тестирования разработанной модели на основании метода группового учета аргументов выполнен сравнительный анализ зависимостей основных технологических параметров нагнетателя ГПА на основании экспериментальных данных и расчетных значений, свидетельствующий об адекватности разработанной модели.
Литература:
1. Ивахненко, А. Г., Юрачковский, Ю. П. Моделирование сложных систем по экспериментальным данным. М.: Радио и связь, 1987, 120 с.
2. Владимиров, В. А. Алгоритм выявления предвестников аварийных остановов газоперекачивающих агрегатов: дис. … кандидата тех. наук: 05.13.01. Тюмень, 2011, 202 с.
3. Group Method of Data Handling (GMDH) for data mining, forecasting algorithms optimization, fuzzy models analysis, statistical learning networks and modeling software systems [Электронный ресурс]: Available at: http://gmdh.net/.
4. Read Microsoft Excel spreadsheet file — MATLAB xlsread [Электронный ресурс]: xlsread. Available at: http://www.mathworks.com/help/matlab/ref/xlsread.html.
5. Solve nonlinear curve-fitting (data-fitting) problems in least-squares sense — MATLAB lsqcurvefit. Available at: http://www.mathworks.com/help/optim/ug/lsqcurvefit.html.
6. Назаренко, І. В., Николайчук, М. Я. Побудова систем управління об’єктами газотранспортної системи на базі уніфікованої технології генерування функціональних блоків з їх математичних моделей, Розвідка та розробка нафтових і газових родовищ ІФНТУНГ, 2012, № 3 (44), С. 177–186.
7. International Standard IEC (International Electrotechnical Commission) 61131–3, 3.0 Ed., Part 3, Programmable controllers — Programming languages, Geneva, 2013, 218 p.
