Scientific journal
Modern problems of science and education
ISSN 2070-7428
"Перечень" ВАК
ИФ РИНЦ = 0,940

ALGORITHM FOR OPTIMAL DISTRIBUTION OF GAS FLOWS IN CLOSED THE GAS-TRANSPORT SYSTEM

Snezhin A.N. 1 Petrov S.V. 1
1 Company «Gazprom development»
In article the mathematical model and realizing her method for calculating the optimal allocation of gas flows in difficult close system of gas pipelines. The novelty of the proposed method consists in addition of a fictitious point in the center of a reversive arch and in a choice coefficient linear form of the target function that allows at introduction of additional characteristics (technically possible productivity, reversible arcs) automatically reformulate a new task for the description in the standard form, for example, the simplex method and to find the solution spectrum of flows tasks. In article the case when the target function is commodity-transport work is taken. The method was developed and successfully tested in Company «Gazprom development» in the process of creating a universal algorithm for solving flow control. This article also contains information received during the development of a method the theorem the decision-making tasks control gas flows on reversible arcs.
Gas-transport system (GTS); united system of gas supply; mathematical model; adjacency matrix; graph; bijection; the distribution of gas flows; commodity-transport work; technically possible productivity(TPP)
Введение

Одной из основных задач, стоящих перед большинством потоковых комплексов, является задача расчета оптимального распределения потоков, при котором некоторая целевая функция стремится к своему экстремуму. При вводе ограничений пропускной способности и управляющих воздействий эта задача по определению [3] становится задачей управления потоками, включающая в себя сразу несколько типов потоковых задач - транспортные задачи, задачи определения существования потока, задачи о поиске обобщенного потока и пр. Сегодня имеется целый ряд решений таких задач. Описание и алгоритмы решений можно найти во множестве учебников и научных трудов, например [1; 2; 5]. Каждое из них имеет свои достоинства, недостатки и условия применимости. Но в силу разных причин степень применимости существующих алгоритмов ограниченна. Поэтому поиск новых подходов, повышающих прикладную пригодность и улучшающих алгоритмы, ведется до сих пор.

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

  • система является сбалансированной (количество поставляемого газа равно количеству потребляемого газа), достигается это за счет изменения запаса газа в газопроводе;
  • газопровод может быть реверсивным и нереверсивным (нереверсивный - газ течет только в одном направлении, от начала трубопровода к его концу; реверсивный - поток газа может течь в любом направлении (рис. 1 случаи а и б), в том числе поток может быть встречным (случай в) или разнонаправленным (случай г));
  • в узлах системы возможно потребление и поступление газа (компрессорные станции, подземные хранилища газа, месторождения и т.д.);
  • в газопроводе возможно потребление и поступление газа;
  • газопровод имеет определенную пропускную способность;
  • возможно задание точек с заданными значениями потока фиксированно или в диапазоне (границы трансгазов, ГИС).

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

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

Рис. 1. Возможные направления потоков газа на реверсивной дуге (на рисунке i- начало газопровода, j - конец газопровода):

а) поток идет от начала к концу;

б) поток идет от конца к началу (реверс);

в) встречный поток (реверс на конце);

г) разнонаправленный поток (реверс на начале - газопровод становится донором за счет запаса газа).

Метод решения задачи

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

В систему входят уравнения баланса газа в узлах. Для каждого узла j составляется уравнение равенства входящих и выходящих потоков:

 (1)

Здесь i и j- вершины входящих в расчетный узел дуг, j и k- вершины дуг, выходящих из расчетного узла;x - величина потока на входах входящих в рассчитываемый узел дуг (xij - в левой части формулы) и выходящих из рассчитываемого узла дуг (xjk - в правой части формулы); g - потребление или распределение на входящих дугах; fj- известный приток или распределение газа в самом рассчитываемом узле (добыча для месторождений, закачка/отбор для ПХГ, экспорт или импорт и т.д.).

Замечание 1: потоки газа, входящие в данную замкнутую систему из других систем и выходящие из неё наружу, имеют вид fj при некоторых значениях индекса j.

Замечание 2: для тех индексов i и j, для которых дуги не существуют, суммирование не производится, переменные и величины не выписываются.

Формула (1) описывает систему ограничений для поставленной задачи. Необходимо учитывать, что задача решается для случая баланса поступления и потребления газа в замкнутой системе за единицу времени, например за один день, т.е. количество (общее) поступления газа в систему и отток из нее совпадают. Математически это означает, что:

 (2)

Перепишем формулу в следующем виде:

 (3)

В уравнениях (1) для величины fj используем знак «-», если это приток газа, и знак «+», если это отток газа. Для величины gij- используем знак «+», если на дуге распределение - приток, и знак «-», если на дуге потребление - отток.

Первая задача, которую необходимо решить, - построить распределение газовых потоков в газотранспортной системе, при условии минимизации товарно-транспортной работы. Пусть cij- это длина газопровода, который соединяет вершину i с вершиной j. Тогда целевая функция будет иметь следующий вид:

 (4)

Уравнения и целевая функция , которую надо минимизировать, составляют постановку задачи линейного программирования. Будем решать поставленную задачу линейного программирования при помощи симплекс-метода [2]. Стандартный вид такой задачи следующий: пусть заданы n действительных чисел ;m действительных чисел ; и m x n действительных чисел , где  и . Требуется найти n действительных чисел x1, x2, ..., n, которые максимизируют линейную форму:

 (5)

при условиях:

 где i = 1, 2, ..., m, (6)

 , где j = 1, 2, ..., n. (7)

В матричной записи это можно записать в виде:

, (8)

, (9)

где: - матрица коэффициентов размера m x n, знаки неравенств должны выполняться для каждой координаты вектора, а скобки  обозначают скалярное произведение векторов.

Для решения задачи  требуется перейти от нумерации переменных и констант двумя индексами к нумерации одним индексом с сохранением взаимно-однозначного соответствия между ними [4].

Пусть m - число всех вершин графа, а n - число всех дуг. Пронумеруем вершины графа уникальными номерами от 1 до m . Тогда ориентированный граф - это множество вершин  и множество упорядоченных пар:

 (10)

Номер ij для переменной xij потока тоже является уникальным, и число таких номеров равно n. Построим взаимно-однозначное отображение пары (i,j) во множество чисел (1,2,...,n) в виде модифицированной матрицы смежности размерностью m x m по следующему правилу:

пусть i- номер вершины, откуда течёт поток, а j- номер вершины, куда втекает поток, тогда в ячейку на пересечении i-й строки и j- го столбца впишем номер дуги l. Для обратного преобразования понадобится матрица дешифратор размерностью 2n, где для каждого номера дуги l в строке l на первом месте будет стоять номер i начала этой дуги, а на втором - месте номер j конца дуги. Так построено взаимно-однозначное соответствие между двойной нумерацией дуг и одинарной.

Полученное соответствие позволяет переформулировать задачу для применения симплекс-метода. Обозначим описанное соответствие буквой T, тогда имеем:

, , , (11)

где , где z - либо переменная, либо заданная величина.

С учетом перенумерации задача (3 - 4) может быть записана следующим образом:

, (12)

, (13)

. (14)

Учитывая, что T - биекция, условие можно записать так:

(15)

Для решения задачи линейного программирования надо составить матрицу A, вектор ограничений b, и вектор с коэффициентами линейной формы c. После этого симплекс-методом вычисляются значения переменных, которые удовлетворяют условиям задачи. Переход от переменных величин с двумя индексами к переменным и величинам с одним индексом осуществляется по следующим правилам: если величина задана в двойной нумерации, её номер лежит в ячейке матрицы смежности, по нему находим значение величины с одним индексом. Если величина задана в одинарной нумерации, подставляем номер в матрицу дешифратор, находим два значения индекса по ним величину в двойной нумерации.

Матрица  заполняется значениями с помощью матрицы смежности. Для начала все элементы матрицы приравниваются к 0. Индекс  нумерует строки матрицы. Заполняется строка j для входящих в узел j дуг. Если значение в ячейке с номером ij равно l, то ajl = 1. Далее заполняется строка j для выходящих из узла j дуг. Если значение в ячейке с номером ji равно l, то ajl = -1. Построение матрицы А завершено.

Вектор  строим следующим образом. Индекс  нумерует координаты данного вектора b. Проведем инициализацию вектора . Далее каждого j перебираем в цикле значения i=1,2,...,m. Если значение в ячейке равно l, то вычитаем из предыдущего значения bj величину .

Аналогично строится вектор . Производится двойной циклический перебор: . Если значение в ячейке равно l, то . Чтобы перейти от задачи в виде равенств (3) к задаче (6) в виде неравенств , следует взять вектор и матрицу . Новая матрица условий A" и новый вектор b" ограничений приобретают следующий вид:  и .

Чтобы найти минимум F(x), следует искать максимум линейной формы , т.е. необходимо произвести замену вектора c на вектор c", где .

Опишем теперь модель распределения потоков газа внутри конкретного газотранспортного общества с учётом ТВПС и требования удовлетворения контрактов, т.е. распределения или потребления на дугах.

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

, (16)

, (17)

, (18)

условия (16-18) только для тех индексов i и j, для которых существуют дуги.

Используя описанное выше соответствие индексов T, имеем:

. (19)

Поэтому можно записать в виде:

, где  (20)

В матричном виде это записывается следующим образом:

, (21)

где - единичная матрица размера , умноженная на -1, .

ТВПС выделенного участка или дуги будем обозначать символом . Эта величина обозначает, какой максимальный поток газа  может войти в дугу, которая начинается в вершине i и заканчивается в вершине j.

Замечание: нестрого - это поток в начале вершины i.

Ограничение ТВПС на дугах запишем так:

. (22)

Используя описанное выше соответствие индексов T, имеем:

, (23)

поэтому (23) можно записать в виде:

, где l=1,2,...,n. (24)

В матричном виде это можно записать:

, (25)

где - единичная матрица размера n x n, , .

Вектор  вычисляется по следующей схеме.

Перебираем в цикле . Используя матрицу дешифратор, получаем двойной индекс ij. По этому индексу находим qij и присваиваем его координате l. То же относится к вектору g.

Строим матрицу , вектор  и вектор линейной формы . Задача в виде (5) - (7) поставлена. Ее решение является решением задачи распределения потоков газа внутри газотранспортного общества с учётом ТВПС, потребления/распределения газа на дугах и в вершинах.

Сделаем очень важное предположение, что имеются реверсивные дуги. Это означает, что допустимы отрицательные потоки и .

Пусть имеется реверсивная дуга ij без потребления на ней самой. Тогда её можно заменить двумя встречными дугами ij и ji. Несложно понять, что два взаимных потока  и  на этих дугах моделируют поток на дуге ij, но с учётом знака, по формуле , где , но xij может быть и отрицательным. В случае потребления на реверсивной дуге ситуация более сложная. Предлагается следующее решение этой задачи. Пусть M- число реверсивных дуг M≥1. Для каждой реверсивной дуги ij добавим новую «фиктивную» вершину , где (рис. 1). После этого дугу ij заменим на две новые и разные дуги ik и kj. Заметим, что т.к. k > m, то новые дуги не будут совпадать со старыми дугами и i ≠ j, поэтому они сами разные. Будем считать, что потребление или распределение газа gij на дуге ij переходит в потребление или распределение газа fв узле , поэтому считаем, что . Дуга ik подлежит замене на две встречные дуги ik и ki, а дуга kj- на дуги kj и jk. Добавление новой вершины означает, что в матрице смежности появляется новая строка и новый столбец. В их объединении есть ячейка с новым номером новой дуги. Встречная дуга ji это транспонированная дуга ij. Длину каждой новой дуги ik, ki, kj и jk определим, как 1/2 длины дуги ij: По смыслу задачи ТВПС  дуги ij присваиваем дугам ik и kj, т.е. , для дуг ki и jk ограничение ставятся аналогично.

Здесь использовалось важное допущение, что потребление/распределение сосредоточено в середине реверсивной дуги.

Описание алгоритма вычислений

Произведённые выше операции с дугами графа приводят к изменению модифицированной матрицы смежности и матрицы дешифратора. Построим «новую» матрицу смежности по следующему правилу. Создадим новую матрицу смежности размера , где m- исходное число вершин графа, M- число реверсивных дуг. Присвоим . Перебираем в цикле все дуги начального графа . Если дуга не реверсивная, то, используя «старый» дешифратор, запишем значение дуги 1 в ячейку «новой» матрицы смежности с тем же номером, который используется в «старой» матрице смежности. Пусть дуга l реверсивная, тогда имеем некий двойной номер дуги: . Берем новый номер фиктивной вершины  и записываем 1 в ячейки с номерами . Одновременно вычисляется «новый» вектор распределения газа в узле . Записываем в специальную матрицу «биекции» размера m x m в ячейку с номером ij значение . Одновременно записываем в «специальную матрицу «дешифратор» размера M*2 в строку с номером k на первое место номер i, на второе место номер j. Нумерация должна начинаться с номера m+1. Увеличиваем k на 1. Действия повторяются со всеми реверсивными дугами до окончания цикла.

Пусть искомая матрица построена, найдём новую сквозную нумерацию дуг, общее число которых . Присвоим r = 1. Будем перебирать ячейки матрицы двойным циклом . Если в ячейке с индексом ij ненулевое значение, тогда ячейке присваиваем значение r и увеличиваем r на 1. В результате получена новая сквозная нумерация ребер  изменённого графа и построена «новая» матрица смежности. Построим «новую» матрицу дешифратора размера . В том же двойном цикле для каждого присвоения индекса r одновременно на первое место в строке r записываем номер i-й строки, на второе - номер j-го столбца. Таким образом, задана новая биекция T':

, (26)

где .

Построим «новый» вектор приток/распределение . Для номеров  имеем . «Новый» вектор (f'j) для номеров  вычислен ранее.

Построим «новый» вектор . Для этого производится перебор значений по . Если дуга с индексом l в старой нумерации не реверсивная, тогда используя «старую» матрицу шифратор находится двойной номер , и в «новой» матрице смежности в ячейке с номером ij берётся значение r. После этого присваиваем . В итоге переберём n-M «старых» дуг, и n-M новых дуг. Если дуга в старой нумерации реверсивная, тогда, используя двойной номер , по матрице «биекции» находим номер k=k(ij) фиктивной вершины. В «новой» матрице смежности из ячеек с номерами  берем последовательно значения r и присваиваем . В результате заданы все  значений для вектора c'. 

Построим «новый» вектор . Производится перебор значений по . Если дуга в старой нумерации не реверсивная, то, используя «старую» матрицу шифратор, находим двойной номер ij=ij(l) и в «новой» матрице смежности в ячейке с номером ij берём значение r. После этого присваиваем . В итоге перебираем n - M «старых» дуг, и n - M новых дуг. Потребление на остальных 4 * M дугах равны 0, так как с реверсивных дуг их значения перенесены на фиктивные вершины. В результате заданы n + 3 * M все значений для вектора g'.

Построим «новый» вектор . Производится перебор значений по l =1,2,...,n. Если дуга в старой нумерации не реверсивная, тогда, используя «старую» матрицу шифратор, находим двойной номер ij=ij(l) и в «новой» матрице смежности в ячейке с номером ij берём значение r. После этого присваиваем . В итоге перебираем n - M «старых» дуг, и n - M новых дуг. Если дуга в старой нумерации реверсивная, тогда, используя двойной номер ij=ij(l), по матрице «биекции» находим номер k = k(ij) фиктивной вершины. В «новой» матрице смежности из ячеек с номерами ik,kj берём последовательно значения  r и присваиваем . Для ячеек с номерами jk,ki поступаем аналогично. В результате заданы все n + 3 * M значений для вектора q'.

Результатом является новая постановка задачи, сведенная к условиям (5)-(7), которую решаем описанным выше методом (без реверсивных дуг). Решение принимает вид , . Далее необходимо вернуться к старым переменным.

Пусть найдено некоторое решение выше поставленной задачи. У нас было n «старых» переменных, стало n + 3 * M «новых» переменных, из них n - M не реверсивных и 4*M относятся к изначально реверсивным дугам. Для нереверсивных дуг (переменных) «новая» переменная означает поток в начале дуги. Поток в конце дуги вычисляется вычитанием из начального потока величины потребления/распределения на дуге. В случае реверсивных дуг поступим следующим образом. Пусть k = k(ij) номер фиктивной вершины. Сложим значения xik и , в результате получим поток в начале дуги ij, причем, если знак положительный - поток течет из вершины i, если знак отрицательный - поток течёт в вершину i. Так вычисляем поток в начале дуги. Сложим значения xkj и . В результате получим поток в конце дуги ij, причём, если знак положительный - поток течет в вершину j, если отрицательный - поток течёт из вершины j в обратную сторону. Так вычисляется поток в конце дуги. Посчитаем количество переменных после такой операции. Имеем n - M переменных для нереверсивных дуг и  переменных для реверсивных дуг. Далее удваиваем число n - M переменных для нереверсивных дуг  и складываем его с числом 2 * M переменных для реверсивных дуг. Получаем 2 * n переменных (поток в начале, поток в конце). Составим «Таблицу решений» задачи. Инициализируем значения таблицы нулевыми значениями. Номер столбца совпадает с номером дуги - переменной, первая строка - поток в начале дуги, вторая строка - поток в конце дуги.

Производится циклический перебор . Используя «новую» матрицу дешифратора определяется двойной номер . Если оба номера , тогда дуга не реверсивная (это признак). Тогда . По «старой» матрице смежности находится «старый» номер дуги  и . Присваиваются значения   и вычисляются значения , так как отток со знаком « — ». В «Таблицу решений» в столбец с номером l записываем в первую строку xl, во вторую yl.

Пусть один из индексов i' или j' больше m. Это значит, что i'j' есть «часть» первоначально реверсивной дуги. Разберем два случая (третьего по построению быть не может).

Случай 1: j'>m. Обозначим через k = j', тогда имеем дугу i'k. Находим по специальной матрице «дешифратору» двойной номер . По «старой» матрице смежности находим «старый» номер дуги . Случай 1 разбивается на 2 случая.

Случай 1.1:i' = i совпало с началом. Тогда в «Таблицу решений» к значению в ячейке пересечения столбца l и первой строки прибавляем x'r (рис. 1, дуга ik).

Случай 1.2: i' = j совпало с концом. Тогда в «Таблицу решений» к значению в ячейке пересечения столбца l и второй строки прибавим (рис. 1, дуга jk).

Случай 2: i' > m. Обозначим через k = j', тогда имеем дугу kj'. Находим по специальной матрице «дешифратору» двойной номер ij = ij (k). По «старой» матрице смежности находим «старый» номер дуги l - l(ij). Случай 2 разбивается на 2 случая.

Случай 2.1:j' = j совпало с концом. Тогда в «Таблицу решений» к значению в ячейке пересечения столбца l и второй строки прибавим x'r. (Поток в начале kj' суть поток в конце ij.) (Рис. 1, дуга kj).

Случай 2.2: j' =i совпало с началом. Тогда в «Таблицу решений» к значению в ячейке пересечения столбца l и первой строки прибавим -x'r. (Поток в начале ki' суть поток в начале ij.) (Рис. 1 дуга ki).

Задача распределения потоков газа внутри газотранспортного общества с учётом ТВПС, потребления/распределения газа на дугах и в вершинах, и с учётом наличия реверсивных дуг решена.

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

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

Существенным предположением для доказательства этой теоремы является условие, что четыре новые дуги, которые заменяют одну реверсивную дугу, получают каждая длину  старой дуги (идея симметрии) и условие переноса потребления/распределения газа со старой дуги на новую «фиктивную» вершину. Актуальность этой теоремы связана с тем, что при практическом использовании описанного выше метода для получения оптимального распределения газа в системе, реверс на конкретной дуге может задаваться пользователем дополнительно, так как является признаком дуги. Теорема утверждает, что такие изменения условий задачи заведомо не ухудшат результат расчета системы, а в практическом смысле часто видно улучшение в смысле значения целевой функции. Добавляя последовательно новые условия, можно моделировать управление потоков газа в системе. Общая формулировка теоремы была получена в результате решения конкретной задачи управления потоками газа, при условии, что задача ставится линейная.

Заключение

Вышеописанный метод решения задачи распределения потоков газа в сложной системе газопроводов с учетом реверсивных дуг, попутного поступления/ потребления в узлах и дугах с учетом ТВПС был реализован в ООО «Газпром развитие» с помощью RAD Delphi 2010. Успешная апробация метода показала его эффективность при расчёте ГТС ЕСГ, что может служить рекомендацией к использованию метода в потоковых комплексах.

Решение с управляющими значениями (границы ГТО, ГИС) и доказательство приведенной теоремы подробно будут рассмотрены в следующей статье.

Выражаем благодарность кандидату физико-математических наук В.М. Простокишину за внимание к работе и конструктивное обсуждение.

Рецензенты

  • Мирошин Николай Васильевич, доктор физико-математических наук, профессор, НИЯУ «МИФИ» Министерства образования и науки РФ, г. Москва.
  • Григорьев Леонид Иванович, доктор технических наук, профессор, РГУ нефти и газа им. Губкина Министерства образования и науки РФ, г. Москва.