Software & Systems

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

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

Разработкой подходов к решению задачи структурной оптимизации ТС занимались многие оте- чественные ученые [1–3]. В исследованиях пос- ледних лет для решения такой задачи активное применение находят различные эвристические алгоритмы [4–8].

Для решения задачи схемно-структурной оптимизации ТС успешно применяется развиваемая в ИСЭМ СО РАН (г. Иркутск) методика избыточных проектных схем [9]. Избыточная схема представляет собой схему системы теплоснабжения, в которую включены всевозможные способы и направления прокладки сети трубопроводов от источников тепловой энергии до потребителей.

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

Постановка задачи

Структура избыточной схемы задается в виде графа G* = (V, E*), где V = {vj | j = 0, 1, …, m} –  множество узлов графа; E* = {ei = {vu, vk}| u ¹ k, i = 0, 1, …, n} – множество ветвей, соответствую- щих участкам схемы, m – число узлов, – число участков. Множество V определяется как объединение непересекающихся множеств V = V0 È V1 È V2, где V0 – множество потребителей; V1 множество источников; V2 – множество узлов разветвления. Избыточная схема дополняется следующими параметрами элементов схемы: Qj  – расход теплоносителя в узле j;  – верхнее и нижнее ограничения на давление в узле j; li  – длина участка i;  – верхнее и нижнее ограничения на скорость движения теплоносителя по участку i.

Структура выделяемой сети теплоснабжения определяется подграфом G = (V, E), где E Í E*. Оптимальная топологическая структура сети теплоснабжения с учетом экономического критерия находится среди множества древовидных схем избыточной схемы [3], поэтому известные алгоритмы ее поиска относительно минимальной функции затрат F(G) ориентированы на топологию сетей древовидной структуры. В общем виде алгоритм оптимальной сети включает несколько взаимосвязанных этапов, в числе которых нормализация структуры (определение и исключение ветвей с нулевым расходом), расчет потокораспределения сети, определение суммарных приведенных затрат. Содержание и последовательность реализации этого алгоритма могут быть представлены следующим образом.

1. Граф G = (V, E), где E Í E*, ориентируется относительно источников теплоносителя ei = (vu, vk), ei Î E, где vu – начальная вершина; vk – конечная вершина.

2. Для каждого узла vj определяется Mj  = = {m | em = (vj, vk)} – множество ветвей, направленных от узла vj.

3. По заданным нагрузкам, начиная с концевых узлов и исходя из материальных балансов в узлах сети, определяются расходы на участках:

.                        (1)

4. Производится нормализация сети, для каждой ветви с нулевым расходом " ei Î B = { ei = (vj, vk) Î E| xi = 0} положим hi = 0, Hi = 0, Mk = Æ, wi = = .

5. Из-за ограничения на скорость течения теплоносителя принято потери давления на ветвях рассчитывать относительно рекомендуемой скорости.

Скорость потока wi  жидкости на участке i рассчитывается по формуле

,               (2)

где wi – скорость течения жидкости, м/с; xi – расход жидкости (объемная скорость), т/ч; p – плотность жидкости, кг/м3; di – внутренний диаметр трубы, мм; p – число Пи.

Потери давления hi на участке i [4] рассчитыва- ются по формуле

,                  (3)

где ji – коэффициент, зависящий от эквивалентной шероховатости труб; ai – коэффициент местных потерь.

Выразив диаметр di  из (2) и заменив его в выражении (3), получим зависимость потерь давления от величины расхода теплоносителя:

.(4)

Потери давления hi на участке i при заданной (рекомендуемой) скорости течения теплоносителя wi рассчитываются по формуле (4) как функции от переменной расхода теплоносителя на участке сети.

6. Рассчитываются значения давления в узлах и действующего напора на участках.

Для концевых узлов давление принимается равным среднему между верхним и нижним ограничениями давления в этом узле:

.                                     (5)

Для остальных узлов vj давление Pj берется равным максимальному рассчитанному среди ветвей из Mj:

.               (6)

Производится проверка соблюдения ограничений на давление [] в узлах, и, если требуется, осуществляется установка насосных (дросселирующих) станций (НС, ДС).

Правила выбора действующего напора НС на участке ei = (vj, vk) Î E/B включают

                       (7)

7. Осуществляется увязка давления в узлах, общих для нескольких участков. Производится корректировка давлений таким образом, чтобы для всех ветвей, выходящих из него, значение давления было равно максимальному из рассчитанных:

Hi = Pk  + hi – Pj, ei = (vj, vk) Î E/B.                 (8)

8. Производится расчет общих суммарных затрат на основе массивов переменных x, h, H для сети G:

,

где – расчетные затраты на сооружение и эксплуатацию трубопровода на участке с но- мером i;  – расчетные затраты на со- оружение и эксплуатацию НС на i-м участке; – расчетные затраты на электроэнергию, расходуемую на перекачку теплоносителя по i-му участку и подачу его потребителю.

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

 (9)

при соблюдении следующих ограничений:

                                          (10)

                        (11)

Эта задача относится к задачам дискретного нелинейного программирования. Среди алгоритмов ее решения хорошо зарекомендовал себя алгоритм целенаправленного перебора деревьев (ПД) [6].

Целенаправленный ПД

Алгоритм ПД заключается в поэтапном улучшении дерева начального приближения путем замены входящих в него участков хордами (участками, не входящими в текущее решение). При включении хорды в дерево образуется цикл. При исключении любого участка, входящего в цикл, будет определено новое дерево [10].

В качестве дерева начального приближения (G0) берутся деревья минимальной суммарной длины, максимальной суммарной длины, минимальных расстояний от потребителей до источника, максимальных расстояний от потребителей до источника, случайно сгенерированные остовные деревья. На каждой t-й итерации алгоритма целенаправленного ПД на базе принятого дерева начального приближения генерируется новое дерево (Gi).

Алгоритм целенаправленного ПД представляет последовательность следующих операций.

1. Определяется начальное приближение G0 = = (V, E0).

2. Формируется массив хорд U  = {u Î E*/E0}.

3. Берется первая хорда u = (vj, vk) из массива хорд u Î U. Хорда u исключается из U. Определяется образовавшийся в результате включения хорды u контур C = Phj D Phk, где Phj Ì E – путь от vj до источника, Phk Ì E  – путь от vk до источника.

4. Генерируются новые решения G Î A, A = = {G = (V, E) | E = {u} + Et / {c}} путем поочередного исключения одного из участков выделенного контура c Î C и включения хорды u.

5. Дерево Gi+1 определяется как дерево с наименьшими расчетными затратами:

                        (12)

" G Î A: (F(G)£ F(Gmin) Þ G = Gmin).

6. Если массив хорд пуст (U = Æ), работа алгоритма завершается, иначе осуществляется переход к шагу 3.

Повышение эффективности алгоритма ПД возможно путем его модификации, заключающейся в динамическом изменении массива хорд, или на основе перехода к новому решению аналогично алгоритму имитации отжига [11].

Модификация динамического массива хорд

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

Модифицированный алгоритм перебора деревьев.

1. Принимается, что n0 соответствует количеству хорд, а n1 представляет счетчик итераций без изменения решения.

2. Определяется начальное приближение G0 = = (V, E0).

3. Формируется массив хорд U  = {u Î E*/E0}.

4. Определяется n0 = |U|, n1 = 0.

5. Берется первая хорда u = (vj, vk)  из массива хорд u Î U. Хорда исключается из U. Определяется образовавшийся в результате включения хорды u контур C = Phj D Phk, где Phj Ì E – путь от vj до источника, Phk Ì E – путь от vk до источника.

6. Генерируются новые решения G Î A, A = = {G = (V, E) | E = {u} + Et / {c}} путем поочередного исключения одного из участков выделенного контура c Î C  и включения хорды u.

7. Дерево Gt+1 определяется как дерево с наименьшими расчетными затратами (12).

8. Если Gt+1 ¹ Gt, то исключенная ветвь добавляется в конец массива хорд U, обнуляется счетчик n1 = 0, иначе рассмотренная хорда добавляется в конец массива U и увеличивается значение счетчика n1 + 1.

9. Если n1 = n0, то нужно завершить алгоритм, иначе перейти к шагу 5.

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

Метод отжига

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

(13)

Начальное приближение определяется как Gb, и производится выбор плана Gn из окрестности текущего плана Gb. Вероятность перехода (13) от плана Gb к Gn регулируется коэффициентом t°, имитирующим падение температуры. Чем выше показатель, тем больше вероятность перехода к плану, хуже текущего F(Gn) > F(Gb). Такой подход позволяет выйти из локального минимума функции F. Если p(Gn, Gb) больше случайного числа из интервала (0, 1), то на следующей итерации рассматривается Gb = Gn. С течением времени значение t° убывает, что обеспечивает сходимость алгоритма к локальному оптимуму.

Необходимо экспериментально выбрать начальное tstart и конечное tstop значения температуры и закон T его изменения. В качестве T предлагается использовать один из следующих подходов:

-     отжиг Коши , где it – номер итерации алгоритма;

-     метод тушения T(t°) = c × t°, где c – некоторая константа, 0,7 < c < 1.

Алгоритм ПД и колец с имитацией отжига

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

1. Принимается, что n0 соответствует количеству хорд, а n1 представляет счетчик итераций без изменения решения.

2. Определяется начальное приближение G0 = = (V, E0).

3. Формируется массив хорд U = {u Î E*/E0}.

4. Определяется n0 = |U|, n1 = 0, t° = tstart.

5. Берется первая хорда u = (vj, vk) из массива хорд u Î U. Хорда u исключается из U. Определяется образовавшийся в результате включения хорды u контур C = Phj D Phk, где Phj Ì E – путь от vj до источника, Phk Ì E – путь от vk до источника.

6. Генерируются новые решения G Î A, A = = {G = (V, E) | E = {u} + Et/{c}} путем поочередного исключения одного из участков выделенного контура c Î C и включения хорды u, Gt+1 = Gt.

7. " G Î A, если p(G, Gt+1) > r при (0 < r < 1), где r – случайное число, то Gt+1 = G.

8. Если Gt+1 ¹ Gt, то исключенная ветвь c добавляется в конец массива хорд U, обнуляется счетчик n1 = 0, иначе рассмотренная хорда u добавляется в конец массива U и увеличивается значение счетчика n1 + 1.

9. Если n1 ³ n0 и t° £ tstop, то завершить алгоритм, иначе t° = T(t°, t), перейти к шагу 5.

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

Сравнительный анализ методов схемно-структурной оптимизации ТС

С целью более компактного изложения результатов исследования введем сокращенные названия алгоритмов:

ПД – исходный алгоритм целенаправленного ПД;

ПДХ – модификация алгоритма ПД с динамическим добавлением хорд в массив хорд;

МОК – модификация ПДХ с законом перехода метода имитации отжига, в качестве закона изменения параметра температуры принят отжиг Коши;

МОТ – модификация ПДХ с законом перехода метода имитации отжига, в качестве закона изменения параметра температуры принят метод тушения;

СП – случайный поиск (генерация и расчет слу- чайных деревьев).

В рамках вычислительного эксперимента сравнивались результаты работы алгоритмов для четырех случайно сгенерированных избыточных и различных по масштабам и параметрам ТС (табл. 1).

Число возможных вариантов решений определяется числом остовных деревьев, связывающих все вершины графа. Структура сети возможного решения является древовидной, но необязательно может быть представлена в виде единого дерева (это может быть лес деревьев). Вместе с тем, если узлы-источники с помощью фиктивных ветвей свести в один узел, то такой граф будет деревом. Для определения количества возможных решений избыточной схемы воспользуемся матричной теоремой о деревьях. Для этого найдем матрицу Кирхгофа K для графа с совмещенными узлами-источниками G*. Число вариантов остовных деревьев равно определителю матрицы, получающейся из матрицы K путем удаления 1-й строки и 1-го столбца.

Характеристики рассчитываемых сетей, включая количество возможных различных деревьев, сведены в таблице 1.

Таблица 1

Характеристики расчетных сетей

Table 1

Characteristics of calculated networks

Каждому расчету соответствует определенное дерево начального приближения. Последовательность расчетов на заранее определенных начальных приближениях представлена в таблице 2.

Таблица 2

Деревья начального приближения, принятые в расчетах, и последовательность их использования в предложенных алгоритмах

Table 2

Initial approximation trees in the calculations and a sequence of their use in the proposed algorithms

Для каждой сети выполнялась серия расчетов по каждому алгоритму. Одному проходу алгоритма соответствует набор входных параметров, включая заранее определенное дерево начального приближения. В результате одного расчета алгоритма определяются структура сети, суммарные затраты и время, потраченное на поиск решения. Серия расчетов алгоритма для сети 4 ограничивалась заданным временем вычислений  для остальных се- тей – количеством расчетов алгоритма  Параметры алгоритмов для тестовых сетей приведены в таблице 3.

Таблица 3

Параметры и ограничения алгоритмов для расчетных сетей

Table 3

Parameters and constraints of the algorithms for calculated networks

Для сопоставления результатов работы алгоритмов аналогично работе [9] определялись следующие показатели для каждой серии расчетов:

minF(G), midF(G), maxF(G) – соответственно минимальное, среднее арифметическое и максимальное значения F(G) среди всех результатов расчетов по всему набору алгоритмов для данной сети (меньше – лучше), тыс. руб.;

midt – среднее арифметическое значение затраченного времени на один расчет рассматриваемого алгоритма (меньше – лучше), cек.;

midG – среднее арифметическое количество рассмотренных сетей за один расчет рассматриваемого алгоритма;

midr – среднее арифметическое отклонение найденных решений от лучшего найденного решения среди всех алгоритмов (меньше – лучше), %;

l – отношение количества расчетов (best) по соответствующему алгоритму, в которых обнаружено лучшее решение F(G) = minF(G), к общему количеству расчетов (all) алгоритма (больше – лучше), %: ;

tmin – расчетное время, необходимое для поиска наилучшей сети с минимальным значением F(G) = = minF(G) (меньше – лучше), cек.: ;

rmin – расчетное количество расчетов алгоритма, необходимое для нахождения сети с F(G) = = minF(G) (меньше – лучше): .

В таблице 4 отражены результаты работы рассмотренных выше алгоритмов. Экспериментальные расчеты выполнялась на компьютере на базе процессора Intel® Core™ i5-4670K.

Визуально результаты всех проведенных исследований согласно рассматриваемым схемам систем теплоснабжения представлены на рисунках 1–4 в виде точечных графиков. Ось абсцисс соответствует значениям функции затрат F(G) для полученных решений, ось ординат отражает требуемое на их поиск время. Каждая точка представляет полученное решение работы алгоритма, а именно значение функции затрат найденной сети и время, затраченное на поиск. Чем точка ближе к началу координат, тем меньше времени потрачено на вычисление и выгоднее решение. Выделенные области решений соответствуют распределению результатов, полученных определенным алгоритмом.

Представленная на рисунках 1–4 область решений, полученная с помощью ПД, имеет наибольший разброс результатов относительно других рассматриваемых алгоритмов. Это объясняется сильной зависимостью получаемого результата от начального приближения. Множества локальных минимумов для первой и второй сетей, найденных алгоритмами МОК, МОД и ПДХ, как видно из рисунков 1 и 2, принимают форму столбцов из точек. Их графическое сопоставление отражает наличие небольшого количества локальных минимумов, при этом, чем столбец плотнее, тем больше начальных приближений приводят к решениям с равными значениями функции затрат в процессе оптимизации. Крайний левый столбец определяет решение, соответствующее условно-глобальному минимуму функции затрат. Высота столбца отражает распределение затраченного времени на поиск решения. Различие затрат времени на поиск решения объясняется отличием начальных приближений. Для каждого начального приближения алгоритмом для поиска оптимального решения выполняется различное количество итераций, на выполнение которых затрачивается различное время расчета. Для сетей 3 и 4 локальные минимумы представлены размытыми областями решений. Это свидетельствует о наличии большого числа минимумов, из которых исследуемые алгоритмы не способны выйти, используя предложенную для них логику вычислений. Вместе с тем распределение точек области решений для алгоритмов МОК, МОД и ПДХ смещено в сторону оптимума функции затрат. Наличие значительного количества обнаруженных локальных минимумов функции затрат свидетельствует о высокой сложности задачи оптимизации для данной сети в отличие от случаев, когда число обнаруживаемых минимумов незначительное, например, для сетей 1 и 2.

Анализ результатов экспериментальных исследований (табл. 4, рис. 1–4) показывает, что предложенная модификация динамического массива хорд обеспечивает прирост эффективности работы алгоритма ПД в 2–5 раз в процентном отношении (midr), но при этом затрачивается больше времени (в 2–3 раза) на один расчет (midt). Однако ПДХ при достаточном количестве расчетов (rmin) находит условно-глобальный минимум (сети 1 и 2) или точки, близкие к нему (сеть 4), в то время как ис- ходный алгоритм ПД не находит глобальные опти- мумы для рассчитываемых сетей.

Предложенные модификации имитации отжига имеют большую вероятность попадания в глобаль- ный оптимум (λ) по сравнению с ПДХ при достаточном количестве расчетов согласно предложенному алгоритму и подходящих параметрах.

Для больших сетей (сеть 4) при ограниченном времени расчета  модификации имитации отжига находят более выгодные решения, чем ПД и ПДХ, несмотря на меньшее количество расчетов. В общем случае для МОК и МОТ время tmin, обеспечивающее нахождение глобального минимума функции (9) или близкого к нему, меньше чем у ПДХ.

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

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

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

Заключение

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

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

Литература

1.     Меренков А.П., Хасилев В.Я. Теория гидравлических цепей. М.: Наука, 1985. 175 с.

2.     Некрасова О.А., Хасилев В.Я. Оптимальное дерево трубопроводной системы // Экономика и мат. методы. 1970. Т. 4. № 3. С. 427–432.

3.     Чупин В.Р., Майзель И.В., Чупин Р.В., Житов А.В. Схемно-структурная оптимизация систем водоотведения поверхностного водостока и ливневой канализации // Изв. вузов. Инвестиции. Строительство. Недвижимость. 2014. № 4. С. 97–105.

4.     Keedwell E., & Khu S.T. A hybrid genetic algorithm for the design of water distribution networks. Eng. Appl. Artif. Intell., 2005, vol. 18, no. 4, pp. 461–472.

5.     Kadu M.S., Gupta R., & Bhave P.R. Optimal design of water networks using a modified genetic algorithm with reduction in search space. Jour. of Water Resources Planning and Management, 2008, vol. 134, no. 2, pp. 147–160.

6.     Zecchin A.C., Simpson A.R., Maier H.R., Marchi A., & Nixon J.B. Improved understanding of the searching behavior of ant colony optimization algorithms applied to the water distribution design problem. Water Resources Research, 2012, vol. 48, no. 9.

7.     Johns M.B., Keedwell E., & Savic D.A. Pipe smoothing genetic algorithm for least cost water distribution network design. Proc. GECCO, 2013, pp. 1309–1316.

8.     McClymont K., Keedwell E., & Savic D., Randall-Smith M. A general multiobjective hyper-heuristic for water distribution network design with discolouration risk. Jour. Hydroinform, 2013, vol. 15, no. 3, pp. 700–716.

9.     Системные исследования в энергетике: ретроспектива научных направлений СЭИ-ИСЭМ; [отв. ред. Н.И. Воропай]. Н.: Наука, 2010. 686 с.

10.   Наумов И.В., Ямщикова И.В. Математическое обоснование выбора оптимизационной модели трассировки электрической сети // Евразийский союз ученых. 2015. № 7-2. С. 123–127.

11.   Лопатин А.С. Метод отжига // Стохастическая оптимизация в информатике. 2005. Т. 1. № 1-1. С. 133–149.