Медвежонков Д. С. УДК [51.001.57+519.853]:656
ТРАНСПОРТНАЯ МОДЕЛЬ С КУСОЧНО-ЗАДАННЫМИ НЕЛИНЕИНЫМИ ИЗДЕРЖКАМИ1
Рассматривается математическая модель транспортной системы, в которой однородный продукт передается по транспортным ветвям из пунктов производства в пункты потребления. Издержки на передачу линейно зависят от объема передачи по ветви, пока объем ниже некоторого критического для каждой ветви значения. Когда объем выше этого значения, издержки возрастают нелинейно.
Модель в виде задачи оптимизации, описывает оптимальное по суммарным издержкам на транспорт распределение потоков по ветвям при заданных предельно возможных объемах снабжения узлов-источников и заданных предельно необходимых объемах снабжения узлов -потребителей.
Предельно возможный объем снабжения узлов-источников не совпадает с предельно необходимым объемом снабжения узлов-потребителей, в связи с этим в целевую функцию входят слагаемые, отражающие штрафы за недопоставку требуемого продукта в узлы-потребители.
Эта модель возникла при проведении исследований в области проблем живучести больших трубопроводных систем - транспортных составляющих систем газо-, нефте- и нефтепродукто-снабжения [3]. В ситуации, когда объем передачи на некоторых ветвях сети выходит на границы пропускной способности и может на некоторую величину превышать их, требуется определить наиболее безопасные маршруты - маршруты, где превышение минимально.
Описание транспортной модели
Транспортная сеть задана направленным графом, имеющим т узлов (источники, потребители или узлы разветвления), и п ветвей (транс-
портные связи). Матрица инциденций А графа содержит элементы а7/, которые равны:
0, если узел 7 не связан с ветвью /, -1, если узел 7 является началом ветви /, +1, если узел 7 является концом ветви /. Для 7-го узла заданы величины й, и Ъг-, г = 1, ..., т . Если Ъг = 0 , а Ъг > 0, то в 7-м узле расположен потребитель с предельно необходимым объемом снабжения Ъг . Если Ъг < 0, а Ъ = 0, то в узле расположен источник, максимально допустимый объем снабжения которого равен |Ъг |. Если й, = 0 и Ъ{ = 0, то в 7-м узле находится узел разветвления без источника и без потребителя. Обозначим 1сош множество номеров узлов, в которых расположены потребители, а - множество номеров узлов, в которых расположены источники. Обозначим переменной Ь вектор
(Ъ,. ., Ът )Т , состоящий из объемов снабжения в узлах сети.
Обозначим переменной х вектор (хг,..., хи )Т, состоящий из объемов перевозок по ветвям сети. Каждая компонента х, вектора х в связи с наличием ограничений снизу на объем потока должна удовлетворять условию хг > хг, где
X! > 0 - заданная величина. Для /-ой ветви задана величина х}-, где х} > , } = 1,..., п . Если объем перевозки X} удовлетворяет условию х}- е [х}-, х}- ] , то издержки на этом интервале заданы линейной функцией от объема. Если х}- > х ■, то издержки
заданы нелинейной функцией от объема.
Обозначим Z - множество непрерывно дифференцируемых функций вещественного аргумен-
1 Работа поддержана грантом РФФИ 09-01-00306-а
та, равных нулю в нуле, производные которых монотонно возрастают и равны нулю в нуле.
Переменные издержки на передачу продукта по ветви ] представимы в виде суммы линейной и кусочно-заданной функций от объема передачи:
+ (Х} ) , (1)
где Я] - коэффициент линейной части, а функция Г; (X;) представима в виде:
F (X) =
(О, Xj < Xj < Xj,
[Fj (xj - Xj), xj > Xj,
где Fj ( • ) - функция из множества Z.
Исходная задача оптимизации
Рассмотрим задачу поиска объемов перевозок по ветвям и снабжения узлов:
P(x,b) = F(Xj) + SjXj)+ Xg(b -b) ^ min (2)
j=i
при ограничениях
Ax = b, x > x,
b < ъ, < о, i g is
0 < b < ъ, i G T
(3)
(4)
(5)
(6)
fj (X) =
I0, Xj <X <xj,
If (xj- xj), xj > xj
(7)
Пусть и = (щ, ..., ит)Т - вектор множителей Лагранжа ограничений (3) исходной задачи. Обозначим п = (Л, •■■, Лп)Т - вектор множителей Ла-
гранжа
для
ограничении-неравенств
(4),
пь = (лЬ, •, ЛЬт)Т и уь = V,..., уьт)Т - векторы множителей Лагранжа для нижних и верхних ограничений-неравенств из (5) и из (6). Построим функцию Лагранжа для задачи (2)-(6):
L(x, b, u, п, Пъ, vъ) = Z Fj (Xj)
j=i
+sjXj + Z (ъ- ъ) -
iGI™n?
-z u ({ax\ - b)-zij (xj - xj) - zvb (b - b>)+
,=1
+
j=1 b
iGl„,
ZV (ъ) - Zvb (ъ) - Zvb (ъ - ъ).
iGl„
iGl„
iG Ir
Целевая функция представляет собой сумму издержек на передачу продукта по ветвям и суммарный штраф за недопоставку продукта в узлы-потребители. Здесь gi - коэффициент штрафа за недопоставку в узел 7. Условие (3) описывает материальный баланс входящих и выходящих потоков во всех узлах. Условия (3), (4) накладывают ограничения на объем перевозок по каждой дуге. Условия (5), (6) накладывают ограничения на объемы снабжения узлов-источников и узлов-потребителей. Векторы х и Ь, для которых выполняются условия (3)-(6), будем называть допустимыми. Будем называть задачу (2)-(6) исходной задачей оптимизации.
Эквивалентная система уравнений к исходной задаче оптимизации
Обозначим fj (•) производную по своему
аргументу функции Г; (•), а /; (X;) - производную по х ■ функции Г; (X;). Заметим, что
Построим систему уравнений и неравенств от переменных х, Ь,и, п, П, Vй . Вначале найдем частные производные от функции Лагранжа по переменным х, Ь, и и приравняем эти производные к нулю, затем составим условия дополняющей нежесткости для ограничений (3)-(6) и перемен-
ь ь ных п, П , V .
]}(х;) + sJ-[лТи[.-Л; = 0, ; = 1,..., п, (8)
-v, = о
U, -ПЪг + Vh, = 0 , i G Ts
(9)
u -g,-vV + V = 0, i gicons, (10)
{Ai], -ъ = 0, i = 1,...,m, (11)
Xj - Xj > 0, Vj (Xj - Xj) = 0, Vj > 0, j = 1,..., n, (12)
ъ -bt > 0, V(ъ -b,) = 0, VV >0, i GIsrc, (13)
ъ < 0, ъ (vb ) = 0, vb > 0 , i G Is
(14)
b - b, > 0, vb (b, - b,) = 0, vb > 0, i e Icons, (15)
b, > 0, b (rtb ) = 0, v, > 0 , i e IconS • (16) Теорема 1. Для того чтобы векторы x, b составляли решение задачи (2)-(6), а векторы u, п, П b, vb являлись множителями Лагранжа ограничений (3)-(6) исходной задачи необходимо и достаточно, чтобы векторы x, b, u, п , П b, v b составляли решение системы (8)—(16).
Доказательство. Система (8)-(16) является записью условий оптимальности Куна-Таккера для задачи (2)-(6). Из выпуклости целевой функции задачи (2)-(6) и линейности её ограничений по теореме Куна-Таккера [1] следует необходимость и достаточность существования решения системы
n
иркутским государственный университет путей сообщения
(8)-(16) для существования решения и оптимальных множителей Лагранжа задачи (2)-(6). Теорема доказана.
Замечание. Переменные rf, vb можно исключить из системы (8)-(14) заменив их следующими выражениями:
rf = (u,)+, vb = (-ut)+, i 6Isrc, (17)
rf =(-ui + gi)+, vb =(ut -gi)+, i 6Ic0ns , (18)
где символом ( )+ обозначена неотрицательная срезка величины в скобках:
(а)+ = max (0, а). Экономическая интерпретация исходной задачи и эквивалентной системы
Вектор переменных u в системе (8)—(16) имеет смысл вектора цен на продукцию в узлах.
Величина (ATu). для j = 1, ..., n будет представлять собой увеличение цены в конечном узле ветви j по сравнению с начальным.
Величина fj (Xj) + sj - rf. совпадает с увеличением цены на ветви j, её будем считать тарифом на перевозку по ветви j=1, ..., n. Согласно (8) увеличение цены по ветви j происходит за счет
предельных издержек f (x.) + s. и дополнительного слагаемого - rf, которое представляет собой индикатор достижения ограничений снизу на величину потока. Если rf > 0 на ветви j, то Xj = Xj
(т.е. достигнуто ограничение снизу), а если rf = 0, то X — Xj. Абсолютное значение rf показывает
насколько тариф ниже предельных издержек на ветви j.
Аналогично, компоненты векторов П, vb - есть индикаторы достижения предельного снабжения в узлах-источниках и узлах-потребителях.
Если в 7-м узле-источнике rf > 0, то b = b (т.е. достигнут уровень предельно допустимого снабжения). Иначе, если rf = 0, то b — bi ■ Если в 7-м узле-потребителе vb > 0, то bt = bt (достигнут уровень предельного необходимого снабжения). Иначе, если vb = 0, то b1 < b1.
Двойственная задача оптимизации
Две дифференцируемые функции Fj (%) и
Ф ■ (у) из Z будем называть сопряженными, если
их производные (%) = (%)/ dx и
р.(у) = dФ■ (у)!dy являются взаимно обратными:
(Р} (У)) = У, (} (1} (1)) = 1.
Это понятие сопряженных является частным случаем сопряженных функций Лежандра и Фенхеля [4, 5].
Введем вектор у = (у, ..., уп )Т, каждая компонента у которого связана с компонентой х. соотношением:
у} = х ) , для } = 1,..., п . (19)
Обратная функция к / . обозначена р., таким образом из (7) и (19) следует, что каждая компонента вектора х связана с компонентой вектора у соотношением:
[х,х], У} = 0 ] [(}(У}) + х, У} > ^
для } = 1,..., п .
Введем функцию Ф ■ (у) из множества 2 сопряженную к ¥■ (%) из 2:
x, =Ф, (У,) = ■
(20)
Ф
у
(у) = jV/ (*)dT , для j
= 1,..., n.
Обозначим символом ф следующую функцию:
ф(У})-ф}(У}) + х}У}, } = 1... п . Функция Ф} (У}) является сопряженной к
функции (х} ) .
Рассмотрим задачу оптимизации с перемен-
Ъ Ъ
ными у, и, п, П , V :
DP(y,n, ць , v ь ) = £ Ф; (y}.) - £ xrfj
j=i
j=i
brfb + £ bvb ^ min (21)
i6^src i6^cons
y.j + Sj-\atu\ - n. = 0, j = 1,..., n, (22)
u - nb + vb = 0, i 6 Isrc, (23)
u, - g, -rf + vb = 0, i 6 Icons, (24)
rfj > 0, j = 1,..., n,
rf — 0, i = 1, ..., m ,
(25)
0
n
n
m
m
V6 > 0, i = 1, ..., m . (27)
Как будет показано далее, задача (21)-(27) имеет с задачей (2)-(6) одинаковые условия оптимальности и переменные в каждой из этих задачах являются множителями Лагранжа для ограничений другой из них. Поэтому задачу (21)-(27) будем называть двойственной задачей оптимизации к исходной задаче (2)-(6).
Теорема 2. Для того чтобы векторы у, и, п , П Ь, V Ь составляли решение задачи (21)-(27), а векторы х, Ь являлись оптимальными множителями Лагранжа для ограничений (22)-(24) двойственной
задачи необходимо и достаточно, чтобы векторы
. ь ь
х, Ь, у, и, п, П , V составляли решение системы
(8)—(16), (19).
Доказательство. Запишем функцию Лагранжа для задачи (21)-(27)
п П
L(у,и,л,Л,у6,х,Ъ) = 1 (Уз)"XХЛз -
з=1 1=1
т т п
- XЪЛ + XЫ "XХз(Уз + -[АТи]з-Л])
-г п
1=1
+ Хъ(иг-лЪ + уЬ) + Хъ(и - gI-лЬ + уЬ) .
Условия Куна-Таккера для задачи (21)-(27) построенные с использованием функции Лагранжа отличаются от системы (8)-(16) следующими условиями:
~(У]) - Х = 0, ] = 1,..., п , (28)
У] + -
к 4-лз=о, 1=1,..., п.
(29)
где условия (28) - дополнительные, а условия (29) используются вместо условий (8).
Замечаем, что функции (х) и ^ (/) взаимообратные. Поэтому равенства (28) перепишем в виде (19). Следовательно, равенства (29) можно переписать в виде (8). По теореме Куна-Таккера с учетом выпуклости целевой функции (21) и линейности ограничений (22)-(27) получаем необходимость и достаточность существования решения системы (8)-(16), (19) для существования решения задачи (21)-(27). Теорема доказана.
Экономическая интерпретация двойственной задачи
Величина Ф ■ (у ■) имеет следующую интерпретацию:
1) если Л]= 0 , то Ф. (у.) - прибыль перевозчика от транспортировки продукта по ветви у, связанная с формированием тарифа по предельным издержкам,
2) если Л]> 0 , то Ф . (у ) - возможная прибыль перевозчика от транспортировки продукта по ветви у, которую он получил бы, если бы на ветви тариф был равен предельным издержкам.
Если Л]> 0 , то х. = Х] . Тариф меньше предельных издержек на величину л. Величина ХЛ представляет собой нехватающий объем
средств, чтобы плата за передачу по ветви у была равна плате по тарифу по предельным издержкам.
Величина - ЬлЪ (которая равна - Ьг (щ X) представляет собой стоимость снабжения /-го узла-источника максимально допустимым объемом продукта при неотрицательной цене щ . Величина
Ъуъ равна Ъi (gi - иг)+ и представляет собой возможные затраты на выплату штрафа за недопоставку продукта в /-й узел в объеме, равном максимально необходимому, за вычетом платы за снабжение узла-потребителя тем же объемом продукта, при цене продажи не превышающей цену штрафа.
Решение исходной задачи алгоритмом внутренних точек
Для решения используем алгоритм внутренних точек [2]. При реализации алгоритма будем считать:
Р} X) = й]Хг.
Выберем начальное приближение х°, Ь° к решению задачи (2)-(6), удовлетворяющее ограничениям-неравенствам (4)-(6). Реализуем итерационный процесс, в котором на каждой итерации будет выбираться направление корректировки из решения вспомогательной задачи и с некоторым шагом вдоль выбранного направления осуществляться корректировка текущего приближения. Если текущее приближение не удовлетворяет ограничениям-равенствам (3), то процесс переходит в первую фазу - фазу ввода в допустимую область. Иначе процесс переходит во вторую фазу - фазу оптимизации в допустимой области.
иркутским государственный университет путей сообщения
1 п
0(Дх, ДЪ) = - ^
(Ах, )2
+
1 т
1 I-
2 }1(х) - х, )2 (дъ, )2
+
2 ^ттЪ - ъ},ъ} - ъ))
->ш1п, (30)
при ограничениях
где
[ЛДх^-ДЪ, -г) = 0, г = 1,..., т, (31)
гЧ . ЪЧ -[лхк
г = 1,..., т.
А для второй фазы - в следующем виде:
п — 1 п _ / ч
I ~ (х) ^ +11 } хК- )2 +
}=1
}=1
+
Шаг вдоль направления корректировки будет осуществляться так, чтобы не выйти за границы допустимой области и при этом:
1) на первой фазе - уменьшить невязки ограничений-равенств (3),
2) на второй фазе - оптимизировать значение целевой функции (2).
Вспомогательную задачу поиска направления корректировки текущего приближения для первой фазы алгоритма представим в виде:
Дх} =
[лТи ] (хЧ - х} )2, }
= 1,
ДЪ, = (-и,)Ш1П(Ъ} -,Ъ} -Ъ Ч)2,г = 1, ..., т
I^Ах] - Ig1 ДЪ +©(Дх, ДЪ) ^ ш1п, (32)
}=! ,е1соп,
при ограничениях
[ЛДх], - ДЪ = 0, , = 1, ..., т . (33) На первой фазе, в отличие от второй, в целевой функции вспомогательной задачи не используются слагаемые аппроксимирующие целевую функцию исходной задачи (2)-(6), а только слагаемые представляющие штрафы за выход переменных на границы ограничений (4)-(6). А в огра-
к ■ 1
ничениях присутствует слагаемые г, , ' = 1, ..., т , отражающие невязки текущего приближения.
Итерационный процесс состоит из следующих шагов.
1 шаг. Поиск направлений корректировки Дх, ДЬ и вектора множителей Лагранжа и. В зависимости от фазы решаем вспомогательную задачу (30)-(31), либо (32)-(33). Для этого с использованием правила множителей Лагранжа записываем условия оптимальности для вспомогательной задачи и выражаем аналитически переменные Дх, ДЬ через неопределенные множители Лагранжа - вектор и . Для задачи (30)-(31) эти выражения следующие:
Для задачи (32)-(33) получаем выражения:
1+г(х}м-х})2 ! }"••..• !
ДЪ, = (-и,)ш1п(ЪЧ -Ъ,,Ъ} -Ъ))2, , е 1„с, ДЪ, = (-и, + g, )ш1п(ЪЧ - Ъ,, Ъ - ЪЧ )2, , е ¡со,,,.
Подставляем эти выражения на первой фазе в (31), а на второй - в (33). Получаем систему линейных уравнений относительно переменных и7, , = 1, ..., т . Решим систему относительно вектора переменных и методом Холецкого. Зная и, вычислим значения переменных Дх, ДЬ .
2 шаг. Проверка критерия останова. Для допустимого приближения хк, Ьк проверим выполнение с заданной точностью условий оптимальности (8)-(16). Если каждое из условий выполнено с точностью е , где е - параметр точности оптимизации, то завершаем работу алгоритма, иначе переходим к следующему шагу.
3 шаг. Расчет величины Я . Найдем величины Я, Я, А3:
Я\ = ш1п {х, -хЧ Дх, : Дх, <0}
11 }'е{1,...п} } ]< ] } ;
А12 = Шт г{Ъ -Ъ-1 ДЪ, : ДЪ, < 0}, Аз = ш1п { Ъ - Ъ,Ч/ДЪ, : ДЪ > 0}.
Найдем величину Я = ш1п{Я11, А2, Я3},
со-
ответствующую кратчайшему шагу до границы области, формируемой ограничениями-
неравенствами (8)-(16).
4 шаг. Выбор величины шага Я .
Если алгоритм находится в первой фазе (режиме ввода в допустимую область), то:
1) если кратчайший шаг Я до границы области, формируемой ограничениями-неравенствами (8)-(16) меньше единицы, то назначаем Я = у\, где у - параметр метода, у е (0;1) (можно, например, взять 0.9).
2) если шаг А > 1, то назначаем Я = 1, переходим к шагу 7 и переводим алгоритм во вторую фазу.
п
п
Если алгоритм находится во второй фазе (режиме оптимизации в допустимой области), то рассчитываем величину Л как решение задачи оптимизации целевой функции зависящей от одной переменной Лк - величины шага вдоль направления корректировки текущего приближения:
Л = аг§тт
Лк е[0Д]
з=1
хк + Лк Ах,) +
+ 5,.(хк + Лк Ах,)
X \ Ъ + Лк АЪ )
. (34)
При решении задачи (39) используем метод дихотомии с критерием останова по длине текущего вложенного отрезка (длина отрезка меньше параметра оптимизации 3) и по удаленности от нуля производной целевой функции (производная по модулю меньше 3).
Далее назначаем Л = тт{Л, Л }. Если после этой операции шаг Л равен кратчайшему шагу Л до границы области, формируемой ограничениями-неравенствами (8)-(16), то назначаем Л = у\ .
6 шаг. Проверка условия выхода из допустимой области. Если во второй фазе алгоритма выполняется условие:
8
тах {г,
Г }>8,
\ г ) 1 АП
100
(35)
Таблица 1
Результаты расчетов
Размер Кол-во Невязки Точность
задачи итера- ограниче- выполне-
(число ций ний ра- ния усло-
узлов, метода венств в вий опти-
число вет- внут- точке ре- мальности
вей) ренних точек шения
(7, 10) 13 1.131*10-10 0.00437
(7, 22) 16 5.222*10-12 0.00571
(25, 30) 23 2.469*10-9 0.00144
(50, 67) 41 2.381*10-/ 0.00332
(100, 116) 68 1.865*10-/ 0.00931
(200, 218) 81 1.484*10-7 0.00641
(200, 240) 93 8.749*10-7 0.00689
то переводим алгоритм в первую фазу и переходим к следующему шагу.
7 шаг. Корректировка текущего приближения с шагом Л.
Вычисляем:
х^1 = хк + ЛАх, (36)
Ь^1 = Ьк + ЛАЬ (37)
и переходим к первому шагу.
Практические расчеты
С использованием алгоритма внутренних точек, реализованного на языке программирования С++, были проведены расчеты на тестовых задачах. Для формирования тестовых задач была разработана специальная программная среда, позволяющая задавать исходные данные модели. Результаты вычислений для рассчитанных примеров приведены в таблице 1.
Проведенные расчеты показали равномерное увеличение числа итераций метода внутренних точек с ростом размера задачи. Реализованный алгоритм позволил успешно решить задачи среднего размера для описанной модели.
БИБЛИОГРАФИЯ
1. Базара М., Шетти К. Нелинейное программирование. Теория и алгоритмы : пер. с англ. М. : Мир, 1982. 583 с.
2. Дикин И. И., Зоркальцев В. И. Итеративное решение задач математического программирования: алгоритмы метода внутренних точек. Новосибирск : Наука. Сиб. отд-ние, 1980. 144 с.
3. Еделев А. В. Разработка специализированной инструментальной среды для исследования проблем живучести больших трубопроводных систем : дис. ... канд. техн. наук : 05.13.18. Иркутск, 2001. 107 с.
4. Зоркальцев В. И., Хамисов О. В. Равновесные модели в экономике и энергетике. Новосибирск : Наука. Сиб. отд-ние, 2006. 221 с.
5. Рокафеллар Р. Выпуклый анализ. М. : Мир, 1973. 470 с.