VladimirDobrovolsky611/Алгоритм SDDP: различия между версиями

Материал из Алговики
Перейти к навигации Перейти к поиску
[непроверенная версия][непроверенная версия]
 
(не показана 131 промежуточная версия 2 участников)
Строка 1: Строка 1:
 +
{{Assignment|ASA}}
 +
 +
Автор описания [[Участник:VladimirDobrovolsky611|В.М.Добровольский]]
 +
 
=  Свойства и структура алгоритма =
 
=  Свойства и структура алгоритма =
 
== Общее описание алгоритма==
 
== Общее описание алгоритма==
Строка 22: Строка 26:
 
Считается, что на 1-м этапе задача детерминирована. Здесь присутствует всего одно состояние <math> \xi_1 =(A_1, \ b_1, \ c_1)  </math>.
 
Считается, что на 1-м этапе задача детерминирована. Здесь присутствует всего одно состояние <math> \xi_1 =(A_1, \ b_1, \ c_1)  </math>.
  
−
[[File:Graph.png |мини|400px|рис.1 Сценарная решетка и выделенный сценарий (синим)]]
+
[[File:Graph.png |thumb|center|600px|Рис.1 Сценарная решетка и выделенный сценарий (синим)]]
  
  
Строка 38: Строка 42:
 
\\x_1 \geqslant 0
 
\\x_1 \geqslant 0
  
−
\end{matrix}\right.</math>
+
\end{matrix}\right. \ \ \ (1)</math>
  
 
где
 
где
Строка 50: Строка 54:
 
\\x_t \geqslant 0
 
\\x_t \geqslant 0
 
\end{matrix}\right.
 
\end{matrix}\right.
−
t=2,...,T-1
+
t=2,...,T-1 \ \ \ (2)
 
</math>
 
</math>
  
Строка 64: Строка 68:
  
 
\end{matrix}\right.
 
\end{matrix}\right.
 +
(3)
 
</math>
 
</math>
 +
 +
При этом функцию <math>Q_t(x_{t-1},\xi_i^t)</math> называют функцией будущего спроса (cost-to-go function).
  
 
Данная постановка задачи позволяет решать задачи стохастического программирования, где в качестве метрики взято мат. ожидание по состояниям системы.
 
Данная постановка задачи позволяет решать задачи стохастического программирования, где в качестве метрики взято мат. ожидание по состояниям системы.
Строка 70: Строка 77:
  
 
===Прикладная область ===
 
===Прикладная область ===
−
Существует несколько прикладных задач, для которых используется данный алгоритм:
+
Существует множество прикладных задач, для которых используется данный алгоритм. далее приведены некоторые примеры:
  
 
1. Управление ГЭС
 
1. Управление ГЭС
Строка 85: Строка 92:
 
Инвестор обладает некоторым стартовым капиталом, который можно инвестировать в заданный набор ценных бумаг. Инвестор формирует финансовый портфель на несколько лет, но внутри данного периода возможны ребалансировки портфеля, проводимые, однако, без вывода и без привлечения дополнительных средств. Прирост цен активов портфеля, а также величина дивидендных выплат - случайны. При этом, транзакционные издержки, начальный капитал инвестора, а также набор возможных активов для вложения - фиксированные величины.  
 
Инвестор обладает некоторым стартовым капиталом, который можно инвестировать в заданный набор ценных бумаг. Инвестор формирует финансовый портфель на несколько лет, но внутри данного периода возможны ребалансировки портфеля, проводимые, однако, без вывода и без привлечения дополнительных средств. Прирост цен активов портфеля, а также величина дивидендных выплат - случайны. При этом, транзакционные издержки, начальный капитал инвестора, а также набор возможных активов для вложения - фиксированные величины.  
 
В данных условиях стоит задача минимизировать финансовые потери инвестора.
 
В данных условиях стоит задача минимизировать финансовые потери инвестора.
 +
 +
=== Алгоритм решения ===
 +
Алгоритм решения задачи (1) - (3) опирается на фундаментальное утверждение:
 +
 +
''cost-to-go function, опсанная уравнением'' (2) ''является кусочно-линейной относительно аргумента'' <math>x_{t-1}</math>
 +
 +
Данное утверждение позволяет сделать вывод, что задача (1) - (3) может быть сведена к задаче линейного программирования.
 +
 +
В ходе доказательства данного утверждения также показывается, что в задачу можно решить итерационно, причем каждую итерацию можно разбить на 3 этапа:
 +
 +
1. Инициализаци
 +
 +
2. Обратный ход
 +
 +
3. Прямой ход
 +
 +
На каждой итерации,в ходе прямого хода задачи достигается верхняя оценка решения, а в ходе обратного - нижняя. Причем при устремлении количества итераций к бесконечности, данные оценки сходятся к истинному решению.
 +
 +
==== Инициализация.====
 +
На данном этапе генерируются случайные сценарии на сценарной решетки задачи. Количество инициализируемых случайных сценариев может варьироваться в зависимости от структуры сетки. До сих пор не исследовано оптимальное с точки зрения скорости выполнения количество пробных сценариев и зависимость выбора этого количества от размерности решетки. Будем считать, что на каждой итерации формируется <math>K</math> случайных сценариев, которые представляют из себя наборы состояний <math>(\widehat{\xi}_{1}^{i_k},...,\widehat{\xi}_{T}^{i_k}) \ k = 1,...,K</math>
 +
 +
Для каждого пробного сценария на каждом шаге решается прямая задача линейного программирования:
 +
 +
<math>
 +
Q_t^{i_k}(\widehat{x}_{t-1}^{k},\widehat{\xi}_{t}^{i_k}) =
 +
\left\{\begin{matrix}
 +
min \ c_t^{i_k}x_t
 +
\\ subject  \ to
 +
\\ A_t^{i_k}x_t + B_t^{i_k}\widehat{x}_{t-1}^k\geq b_t^{i_k}
 +
 +
\end{matrix}\right.
 +
\ \ \ (2)
 +
 +
</math>
 +
 +
Получается <math>K</math> пробных решений:
 +
 +
<math>
 +
\widehat{x}_t^k = \left\{\begin{matrix}
 +
argmin \ c_t^{i_k}x_t
 +
\\ subject  \ to
 +
\\ A_t^{i_k}x_t + B_t^{i_k}\widehat{x}_{t-1}^k\geq b_t^{i_k}
 +
 +
\end{matrix}\right. \ \  \ (3)
 +
 +
</math>
 +
 +
<math>k = 1,...,K</math>
 +
 +
====Обратный ход ====
 +
На обратном ходе с помощью каждого пробного решения формируется отсекающая гиперплоскость, ограничивающая решение снизу.
 +
[[File:Отсекающая гиперплоскость.bmp.jpg|thumb|center|600px|Рис.2 Отсекающие гиперплоскости ограничивают истинную cost-to-go function снизу]]
 +
 +
Уравнение гиперплоскости получается в ходе решения прямой и обратной задач линейного программирования:
 +
 +
На шаге <math>T</math>:
 +
 +
Прямые задачи:
 +
 +
<math>
 +
Q_T(\widehat{x}_{T-1}^k,\xi_T^i) =
 +
\left\{\begin{matrix}
 +
min \ c_T^{i_k}x_T
 +
\\ subject  \ to
 +
\\ A_T^{i_k}x_T + B_T^{i_k}\widehat{x}_{T-1}\geq b_T^{i_k}
 +
\end{matrix}\right.
 +
(4)
 +
</math>
 +
 +
 +
Двойственные задачи:
 +
 +
<math>
 +
Q_T(\widehat{x}_{T-1}^k,\xi_T^i) =
 +
\left\{\begin{matrix}
 +
max \ \pi_T^{i_k}(b_t^i - B_T^i\widehat{x}_{T-1}^{i_k})
 +
\\ subject  \ to
 +
\\ \pi_T^{i_k}A_T^{i_k}\leq c_T^{i_k}
 +
\end{matrix}\right.
 +
(5)
 +
</math>
 +
 +
Уравнение отсекающей гиперплоскости:
 +
 +
<math>
 +
l_T^k(x_{T-1}^k) = \sum_{j=1}^{m_T} \ p_{Tj}c_T^jx_t^{i_kj}+ (-\sum_{j=1}^{m_T}p_{tj}\pi_t^{i_kj}B_t^j)(x_{T-1}^{i_k} - \widehat{x}_{T-1}^{i_k})
 +
\ \ (6)
 +
</math>
 +
 +
На этапах <math>t = 1,...,T-1</math>:
 +
 +
Прямые задачи:
 +
 +
<math>
 +
Q_t(x_{t-1}^k,\xi_t^i) =
 +
\left\{\begin{matrix}
 +
min \ c_t^{i_k}x_t + \gamma_t
 +
\\ subject  \ to
 +
\\ A_t^{i_k}x_t + B_t^{i_k}\widehat{x}_{t-1}\geq b_t^{i_k}
 +
\\ \gamma_t\geq l_{t+1}(x_t^{i_k})
 +
\end{matrix}\right.
 +
(7)
 +
</math>
 +
 +
 +
Двойственные задачи:
 +
 +
<math>
 +
Q_t(\widehat{x}_{t-1}^k,\xi_t^i) =
 +
\left\{\begin{matrix}
 +
max \ \pi_T^{i_k}(\widetilde{b}_t^i - \widetilde{B}_t^i\widehat{x}_{t-1}^{i_k})
 +
\\ subject  \ to
 +
\\ \pi_t^{i_k}\widetilde{A}_t^{i_k}\leq \widetilde{c}_t^{i_k}
 +
\end{matrix}\right.
 +
(8)
 +
</math>
 +
 +
Уравнение отсекающей гиперплоскости:
 +
 +
<math>
 +
l_t^k(x_{t-1}^k) = \sum_{j=1}^{m_t} \ p_{tj}\widetilde{c}_t^jx_t^{i_kj}+ (-\sum_{j=1}^{m_t}p_{tj}\pi_t^{i_kj}\widetilde{B}_t^j)(x_{t-1}^{i_k} - \widehat{x}_{t-1}^{i_k})
 +
\ \ (9)
 +
</math>
 +
 +
В результате, значение функционала при оптимальном решении в уравнении (7) является нижней оценкой решения <math>\underline{s}</math>
 +
 +
==== Прямой ход ====
 +
На прямом ходе алгоритма производятся те же действия, что и при инициализации, однако при решении прямых задач линейного программирования дополнительными ограничениями выступают отсекающие гиперплоскости, полученные на обратном ходе алгоритма.
 +
 +
На прямом ходе алгоритма генерируются пробные сценарии <math>\widehat{\xi}_t^k; \ k = 1,...,K</math>
 +
 +
Для каждого из K сценариев решаются задачи (7), в результате получается верхняя оценка <math>\overline{s}</math>
 +
 +
Если <math>|\overline{s} - \underline{s}| \leq \varepsilon</math> то решение найдено. Иначе, переход на следующую итерацию.
  
 
== Вычислительное ядро алгоритма==
 
== Вычислительное ядро алгоритма==
 +
С вычислительной точки зрения, алгоритм SDDP является композицией большого количества вычислений задач линейного программирования сравнительно небольших размерностей.
 +
С точки зрения компьютерной реализации, алгоритм состоит из двух фаз:
 +
1. Вычисление прямой или двойственной задачи линейного программирования. Получения значения функционала и оптимального решения
 +
2. Запись в пространство гиперплоскостей коэффициентов задачи.
 +
Следует отметить, что при каждой новой итерации, пространство гиперплоскостей пополняется новыми отсекающими гиперплоскостями, полученными посредством решения более приближенных задач, чем на предыдущих итерациях. Поэтому, размерность задач линейного программирования растет(т.к. растет количество ограничений) в то время, как эффективность данных задач не увеличивается. Для этого был введен также алгоритм по удалению ненужных отсекающих гиперплоскостей.
 +
 +
[[File:Up down.jpg|thumb|center|600px|Рис.3 Типичное поведение верхней и нижней границы решения с увеличением итераций]]
  
 
== Макроструктура алгоритма==
 
== Макроструктура алгоритма==
 +
Алгоритм SDDP является композицией задач линейного программирования.
 +
Пусть выполнены следующие условия:
 +
 +
<math>x = (x_1,...,x_N);</math>
 +
 +
<math>\pi = (\pi_1,...,\pi_M);</math>
 +
 +
<math>A \in R^{N \ X \ M};</math>
 +
 +
<math>c \in R^N
 +
 +
</math>
 +
 +
Тогда, задачей линейного программирования называется задача следующего вида:
 +
 +
<math>
 +
\left\{\begin{matrix}
 +
min \ c^Tx
 +
\\ subject to
 +
\\ Ax \geq b
 +
\\ x_i \geq 0; \ i = 1,...,N
 +
\end{matrix}\right. (10)
 +
</math>
 +
 +
Двойственной задачей к задаче (10) является задача следующего вида:
 +
 +
<math>
 +
\left\{\begin{matrix}
 +
max \ \pi^Tb
 +
\\ subject to
 +
\\ \pi^TA \leq c
 +
\\ \pi_i \geq 0; \ i = 1,...,M
 +
\end{matrix}\right. (11)
 +
</math>
 +
 +
Задача (10) и задача (11) могут быть решены Симплекс методом
 +
 +
Обозначим основные свойства данных задач и Симплекс метода:
 +
 +
Пусть <math> M</math> - количество ограничений, <math>N</math> - количество переменных задачи линейного программирования. Тогда:
 +
 +
1. Практика показывает, что в среднем вычислительная сложность алгоритма составляет <math> M^3 </math> операций
 +
 +
2. Количество итераций знаключение между <math> M </math>и <math>2M </math>
 +
 +
3. В худшем случае ввиду комбинаторных особенностей алгоритма, симплекс метод имеет экспоненциальную сложность, однако в среднем, задача решается с полиноминальной сложностью.
 +
 +
4. Количество ограничений влияет на время выполнения метода гораздо сильнее, чем количество переменных
 +
 +
5. Оптимальное значение функционала задачи (10) совпадает с оптимальным значением функционала (11)
 +
 +
6. Для решения алгоритма SDDP использовался готовый солвер для решения задач линейного программирования - lpSolve
  
 
== Схема реализации последовательного алгоритма==
 
== Схема реализации последовательного алгоритма==
 +
Алгоритм SDDP начинается с инициализации пробных сценариев, данный этап описан в разделе 1.2.3.1
 +
 +
После генерации пробных сценариев и вычисления пробных решений, производится обратных ход алгоритма, описанный в 1.2.3.2
 +
 +
На данном этапе на каждом этапе формируются поддерживающие гиперплоскости, а также составляется нижняя оценка решения
 +
 +
После проведения обратного хода в алгоритме запускается прямой ход, описанный в 1.2.3.3, вычисляется верхняя оценка решения.
 +
 +
В случае, если норма разности верхней и нижней оценки оказывается меньше заранее заданного числа, алгоритм заканчивает работу, любая из двух оценок является оптимальным решением.
 +
 +
В случае, если норма разности верхней и нижней оценки оказывается выше заранее заданного числа, алгоритм повторяет итерацию.
 +
[[File:Scheme SDDP.jpg|thumb|center|600px|Рис.4 Схема SDDP]]
  
 
== Последовательная сложность алгоритма==
 
== Последовательная сложность алгоритма==
 +
Обозначим вычислительную решения задачи линейного программирования симплекс-методом с размерностью переменных оптимизации N и количеством ограничений M за <math>lp(N,M)</math>
 +
Пусть задача SDDP имеет T этапов, размерность вектора переменных - N, размерность ограничений - M
 +
Количество состояний на каждом этапе - <math>S_t; \ t = 1,...,T</math>
 +
Пусть также на каждой итерации генерируется k пробных сценариев
 +
Тогда для первой итерации алгоритма SDDP вычислительная сложность соответствует значению
 +
 +
<math>
 +
kTlp(N,M) + k(lp(N,M)+lp(M,N)) + k(T-2)(lp(N,M+k)+lp(M+k,N)) + lp(N,M+k)
 +
</math>
 +
 +
С увеличением номера итерации, вычислительная сложность увеличивается в худшем случае линейно, в среднем, пропорционально корню номера итерации (это происходит потому, что некоторые отсекающие гиперплоскости удаляются)
 +
 +
На j-той итерации:
 +
<math>
 +
kTlp(N,M) + k(lp(N,M)+lp(M,N)) + k(T-2)(lp(N,M+k\sqrt{j})+lp(M+k\sqrt{j},N)) + lp(N,M+k\sqrt{j})
 +
</math>
 +
 +
Таким образом, сложность алгоритма при K итерациях составляет:
 +
 +
<math>
 +
P = \sum_{j=1}^{K} \ kTlp(N,M) + k(lp(N,M)+lp(M,N)) + k(T-2)(lp(N,M+k\sqrt{j})+lp(M+k\sqrt{j},N)) + lp(N,M+k\sqrt{j})
 +
</math>
 +
 +
 +
Вообще говоря на решетке с такими данными возможно <math>\prod_{t=1}^TS_t</math> сценариев. Также, в худшем случае, в ходе решения не была удалена ни одна отсекающая гиперплоскость
 +
Тогда в худшем случае, решение будет достигнуто с вычислительной сложностью:
 +
 +
<math>
 +
P = \sum_{j=1}^{\frac{\prod_{t=1}^TS_t}{k}} \  kTlp(N,M) + k(lp(N,M)+lp(M,N)) + k(T-2)(lp(N,M+kj)+lp(M+kj,N)) + lp(N,M+kj)
 +
</math>
 +
 +
Так как средняя вычислительная сложность симплекс метода <math>lp(N,M) = M^3</math>, вычислительная сложность алгоритма SDDP в худшем случае составляет
 +
 +
<math>
 +
P = \sum_{j=1}^{\frac{\prod_{t=1}^TS_t}{k}} \  kTM^3 + k(M^3+N^3) + k(T-2)((M+kj)^3+N^3) + (M+kj)^3
 +
</math>
 +
 +
В терминах эквивалентности:
 +
 +
<math>
 +
P \sim O(\overline{S}^{2T}*(M^3+N^3))
 +
</math>
 +
 +
Как было показано, в худшем случае сложность алгоритма экспоненциально зависит от количества этапов. Однако весь замысел алгоритма заключается в том, чтобы не перебирать все возможные сценарии.
  
 
== Информационный граф==
 
== Информационный граф==
 +
Информационный граф 1 итерации алгоритма.
 +
 +
На Рис.5 изображен обратный ход алгоритма, проход от последних состояний к корневым
 +
 +
На Рис.6 изображен прямой ход алгоритма по сценарию, выделенному зеленым
 +
 +
Под ОП1 следует понимать решение прямой и двойственной задачи ЛП, составление коэффициентов отсекающей гиперплоскости.
 +
 +
Под ОП2 следует понимать решение прямой задачи ЛП с использованием составленных на обратном ходе отсекающих гиперплоскостей.
 +
 +
[[File:Reverse step.jpg|thumb|center|600px|рис.5 Информационный граф 1 итерации алгоритма SDDP (обратный ход)]]
 +
 +
 +
[[File:Forward step.jpg|thumb|center|600px|рис.6 Информационный граф 1 итерации алгоритма SDDP (прямой ход)]]
  
 
== Ресурс параллелизма алгоритма ==
 
== Ресурс параллелизма алгоритма ==
 +
При параллельной реализации алгоритма выгодно с точки зрения быстродействия увеличивать количество пробных сценариев на каждой итерации, так как это уменьшает количество необходимых итераций алгоритма. При условии, что в расположении имеется вычислительный кластер бесконечного размера, решение задачи можно гарантированно получить за 1 итерацию. Для этого потребуется <math> \prod_{t=1}^T \ S_t </math> вычислительных узлов 1-го уровня, каждый из которых будет рассчитывать коэффициенты отсекающей гиперплоскости для каждого сценария. Внутри каждого узла производится последовательное, от этапа к этапу, вычисление прямых и двойственных задач линейного программирования. В рамках данной статьи не рассматривается рессурс параллельного алгоритма симплекс метода, к тому же lpSolve доступен только в скалярной реализации. Поэтому, в рамках данной статьи будем считать, что решение задачи линейного программирования - задача для одного процессора. С учетом таких допущений, параллельная вычислительная мощность может составить <math>P \sim O(max(M^3,N^3)T)</math>
  
 
== Входные и выходные данные алгоритма ==
 
== Входные и выходные данные алгоритма ==
 +
Параметры задачи:
 +
 +
<math>
 +
T \in \mathbb{N}
 +
</math> - количество этапов задачи
 +
 +
<math>
 +
N \in \mathbb{N}
 +
</math> - размерность вектора переменных задачи
 +
 +
<math>
 +
S \in \mathbb{N}^T
 +
</math> - вектор количества состояний на каждом этапе
 +
 +
Входные данные:
 +
 +
<math>p</math> - вероятности перехода от состояния на предыдущем этапе в состояние на следующем
 +
 +
<math>
 +
p_t \in [0,1]^{S_t \times S_{t+1}}; \ t = 1,...,T
 +
</math>
 +
 +
<math>
 +
\sum_{i=1}^{S_{t+1}} \ p_t^{j \rightarrow i} = 1; \ \forall j = 1,...,S_t
 +
</math>
 +
 +
<math>\xi</math> - состояние
 +
 +
Каждое состояние описывается набором данных
 +
 +
<math> \xi_t^i = \{A_t^i,B_t^i,c_t^i,b_t^i\}; \ t=1,...,T; \ i = 1,...,S_t</math>
 +
 +
где:
 +
 +
<math> A_t^i \in \mathbb{R}^{N\times M}</math>
 +
 +
<math> B_t^i \in \mathbb{R}^{N\times M}</math>
 +
 +
<math> c_t^i \in \mathbb{R}^{N}</math>
 +
 +
<math> b_t^i \in \mathbb{R}^{M}</math>
 +
 +
SDDP применяется для задач с количеством этапов порядка 100, количеством состояний порядка 20 и размерностью решения и ограничений порядка 50 на 50. Таким образом, в памяти на протяжении всей программы должно храниться 100*20*(2*50+2*50*50) = 1 000 000 000 чисел типа double, т.е. около 8 гб.
 +
 +
Выходные данные алгоритма:
 +
 +
<math>X \in \mathbb{R}^{N\times T}</math> - матрица решения.
  
 
== Свойства алгоритма ==
 
== Свойства алгоритма ==
 +
1. К отрицательным свойствам алгоритма необходимо отнести его недетерминированность. В зависимости от входных данных а также от выбранных пробных сценариев зависит количество итераций, которое необходимо проделать для того, чтобы добиться удовлетворительной погрешности.
 +
 +
2. Вычислительная мощность алгоритма в худшем случае:
 +
<math>
 +
P \sim O(\overline{S}^{2T}*(M^3+N^3))
 +
</math>
 +
 +
3. Вычислительная мощность алгоритма в среднем:
 +
<math>
 +
P \sim O(\overline{S}*2T^3*(M^3+N^3))
 +
</math>
 +
 +
= Программная реализация алгоритма =
 +
 +
== Особенности реализации последовательного алгоритма ==
 +
Важной особенностью реализации последовательного алгоритма SDDP является выбор количества пробных сценариев на каждой итерации. С одной стороны, увеличив количество пробных сценариев на каждой итерации, задача быстрее усложняется с вычислительной точки зрения от итерации к итерации. С другой стороны, увеличивается и вероятность нахождения оптимального решения на более ранних итерациях. До сих пор нет единого мнения о том, какое количество пробных сценариев - оптимальное. Авторы первых статей об SDDP применяли всего 1 пробный сценарий на итерацию, аргументировав это тем, что каждое следующее вычисление новой гиперплоскости использует уже все существующие, тем самым повышая скорость сходимости.
 +
 +
==Локальность данных и вычислений ==
 +
На практике алгоритм SDDP необходим для решения задачи с количеством этапов порядка 100, количеством состояний порядка 20 и размерностью решения и ограничений порядка 50 на 50. Таким образом, в памяти на протяжении всей программы должно храниться 100*20*(2*50+2*50*50) = 1 000 000 000 чисел типа double, т.е. около 8 гб. Также, каждому состоянию в ходе программы привязывается набор коэффициентов гиперплоскости, с каждой итерацией добавляются новые наборы. Так, на 100 - итерации, в памяти должно храниться около 100*20*51*100 = 10 200 000 т.е. около 80мб. Обращение к входным данным имеет временную неравномерность. Если количество пробных сценариев = k, то один раз за итерацию к каждой матрице каждого этапа алгоритм обратится k раз, причем почти одновременно. Таким образом, можно ускорить выполнение алгоритма в условиях недостаточной быстрой памяти с помощью загрузки данных нескольких этапов в кэш. При параллельной реализации можно предоставить памяти, выделенной под исходные данные, общий доступ, так как запись в эту область не предусматривается.
 +
 +
Другая область памяти, необходимая алгоритму - пространство коэффициентов гиперплоскости. В отличие от области входных данных, данная область растет с увеличением количества итераций алгоритма пропорционально корню номера итерации. Обращение к данному участку памяти происходит также неравномерно, т.е. k почти одновременных обращений раз в итерацию. Также нельзя без дополнительных условий организовать общий доступ к данной области памяти, т.к. в ней предусмотрены как запись новых значений, так и удаление старых.
 +
 +
==Возможные способы и особенности параллельной реализации алгоритма ==
 +
Параллельная реализация заключается в том, чтобы давать разным вычислительным узлам высчитывать обратный ход для разных пробных сценариев. Внутри же каждого сценария на каждом этапе решается прямая и двойственная задачи, которые могут решаться параллельно. Таким образом, предлагается следующая структура параллельного исполнения алгоритма:
 +
 +
Задается некоторое число вычислительных узлов. В параметрах алгоритма указывается количество пробных сценариев, кратное количеству вычислительных узлов. Каждый вычислительный узел берет на себя вычисления обратного хода по нескольким сценариям. Синхронизация происходит на каждом этапе, т.к. в этот момент загружаются новые гиперплоскости, полученные от решения задач ЛП в каждом пробном сценарии. Данная стратегия повторяется на каждом этапе вплоть до 1-го. Далее начинается прямой ход, в котором каждый вычислительный узел решает задачи ЛП по набору состояний на следующем этапе. Далее проводится процедура удаления лишних гиперплоскостей и инициализация следующей итерации.
 +
Барьерные функции должны отрабатывать на каждом этапе обратного хода.
 +
 +
Возможно ускорение также при параллельной реализации задач ЛП, однако как показывают реальные опыты, задачи ЛП решаются за незначительное время и не время их выполнения на одном процессоре сопоставимо с обменом данными между несколькими при параллельной реализации.
 +
 +
==Масштабируемость алгоритма и его реализации ==
 +
Ввиду того, что существуют реализации алгоритма только на языке R и Matlab, которые не установлены ни на Blue gene (установить невозможно), ни на Lomonosov(в процессе установки) Результаты расчетов были проведены лишь для количества задействованных ядер многоядерного домашнего компьютера.
 +
Параметры запуска:
 +
 +
Процессор Intel Core i5 (4 ядра)
 +
 +
Программа написана на языке R, программа запуска - R.exe (R - не компилируемый, а интерпритируемый язык, изначально приложение R.exe написано на С и было откомпилировано для операционной системы Windows до установки на ПК)
 +
 +
Версия исполняемого скрипта - мой собственный код, который я пишу в рамках магистерской диссертации.
 +
 +
Используемые библиотеки: lpSolve, Rmpi, snow, Rcluster
 +
 +
Результаты тестов:
 +
 +
{|border = "1"
 +
!Размерность данных (T;S;N;M)
 +
!1
 +
!2
 +
!4
 +
|-
 +
!(5;5;10;10)
 +
|279
 +
|159,064
 +
|86,24
 +
|-
 +
!(10;10;10;10)
 +
|1120,8
 +
|608,402
 +
|323,239
 +
|-
 +
!(28;59;6;6)
 +
|62344363,245
 +
|Na
 +
|Na
 +
|}
 +
 +
Выводы о масштабируемости:
 +
 +
Данный алгоритм хорошо масштабируется по измерению T, т.е. в случае, когда количество этапов T - достаточно велико, использование дополнительных процессоров заметно ускоряет решение задачи.
 +
 +
==Динамические характеристики и эффективность реализации алгоритма==
 +
Те немногие данные, которые удалось получить посредством экспериментов запуска алгоритмов на  многоядерном процессоре домашнего компьютера не позволяют сделать глубоких и исчерпывающих выводов. Однако, судя по таблице, указанной в параграфе 2.6 данной статьи можно сделать вывод, что при бОльшем количестве вычислительных узлов, время выполнения алгоритма алгоритма растет медленнее при росте числа этапов.
 +
 +
==Выводы для классов архитектур==
 +
К сожалению, ввиду того, что не существует универсальных реализаций на низкоуровневых языках программирования, трудно дать какие либо рекомендации касательно того, какие архитектуры лучше подходят для данного алгоритма.
 +
 +
==Существующие реализации алгоритма==
 +
На данный момент существует пакет Fast в матлабе, который позволяет запускать алгоритм SDDP с мат. ожиданием в качестве метрики на скалярном компьютере.
 +
Также в разработке находится проект, написанный на R, позволяющий использовать в качестве метрики как мат.ожидание, так и VaR, так и CVaR, разрабатывается параллельный алгоритм вычисления с помощью пакета Rmpi.
 +
 +
=Литература=
 +
1. M.V.F. Pereira and L.M.V.G. Pinto "Multi-stage stochastic optimization applied to energy planning" Mathematical Programming 52 (1991) 359-375
 +
 +
2. Alexander Shapiro, Darinka Dentcheva, Andrzej Ruszczy´nski "Lectures on Stochastic Programming: Modeling and Theory"

Текущая версия на 15:32, 8 февраля 2017

Symbol wait.svgЭта работа прошла предварительную проверку
Дата последней правки страницы: 08.02.2017
Данная работа соответствует формальным критериям.
Проверено ASA.


Автор описания В.М.Добровольский

Содержание

1 Свойства и структура алгоритма

1.1 Общее описание алгоритма

Стохастическое двойственное динамическое программирование (SDDP) – это метод оптимизации, предназначенный для решения динамических задач в условиях неопределенности, то есть в случае, когда некоторые параметры задачи не являются детерминированными. В силу фундаментальности постановки задачи, данный алгоритм может быть применен в самых различных прикладных областях. Например, на сегодняшний день, стохастическое двойственное динамическое программирование активно используется для управления ГЭС в Норвегии, а также вводится в банках для управления рыночными рисками. На сегодняшний день также широко распространены альтернативные динамические методы поиска решений в условиях неопределенности, например, методы, работающие на принципах построения дерева возможных исходов, или методы, работающие на принципах управляющих правил. Однако, методы, работающие по принципу построения дерева, неизбежно сталкиваются с широко известным «проклятием размерности» (curse of dimensionality), а методы, построенные на принципах управляющих правил, как правило, требуют серьезные ограничения на тип управляющих правил, а также на свойства стохастических параметров задачи. Также, в задачах динамического управления присутствует проблема тайм-консистентности решения (time-consistence solution).

Алгоритм SDDP (Stochastic dual dynamic programming) впервые был предложен в статье M.V.F. Pereira и L.M.V.G. Pinto "Multi-stage stochastic optimization applied to energy planning" в 1991 году. Далее алгоритм претерпел множество модернизаций и спецификаций, описанных в труде Alexander Shapiro, Darinka Dentcheva, Andrzej Ruszczynski "Lectures on Stochastic Programming: Modeling and Theory", Теперь под аббревиатурой SDDP подразумевается целое семейство алгоритмов.

1.2 Математическое описание алгоритма

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

Исходные данные:

1. Количество этапов T, количество состояний на каждом этапе [math]m_t[/math], [math]t = 1,...,T[/math]

2. Размерность задачи N (размерность управляющего правила)

3. Вероятности переходов [math]p_{nt}; t = 1,...,T; n = 1,...,m_t[/math]

4. матрицы и векторы, характеризующие каждое состояние системы [math] \xi_i^t=(A_i^t, \ B_i^t, \ b_i^t, \ c_i^t) [/math]

Совокупность входных параметров в пунктах 1 - 4 формируют сценарную решетку задачи (см. рис. 1)

Считается, что на 1-м этапе задача детерминирована. Здесь присутствует всего одно состояние [math] \xi_1 =(A_1, \ b_1, \ c_1) [/math].

Рис.1 Сценарная решетка и выделенный сценарий (синим)


Вычисляемые данные: Матрица управляющих действий [math]X[/math] (элементы [math]x_{it}; \ t \in 1,...,T; \ i \in 1,...,N[/math] - управления для i-го элемента на шаге t)


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

[math] \left\{\begin{matrix} min \ c_1x_1 + \sum_{i=1}^{m_2}p_{1i}^1Q_2(x_1,\xi_i^2) \\subject \ to \\A_1x_1\geqslant b_1 \\x_1 \geqslant 0 \end{matrix}\right. \ \ \ (1)[/math]

где

[math] Q_t(x_{t-1},\xi_i^t) = \left\{\begin{matrix} min \ c_t^ix_t + \sum_{j=1}^{m_{t+1}}p_{ji}^tQ_{t+1}(x_t,\xi_j^{t+1}) \\subject \ to \\A_t^ix_t + B_t^ix_{t-1}\geqslant b_t^i \\x_t \geqslant 0 \end{matrix}\right. t=2,...,T-1 \ \ \ (2) [/math]

...

[math] Q_T(x_{T-1},\xi_i^T) = \left\{\begin{matrix} min \ c_T^ix_T \\subject \ to \\A_T^ix_T + B_T^ix_{T-1}\geqslant b_T^i \\x_t \geqslant 0 \end{matrix}\right. (3) [/math]

При этом функцию [math]Q_t(x_{t-1},\xi_i^t)[/math] называют функцией будущего спроса (cost-to-go function).

Данная постановка задачи позволяет решать задачи стохастического программирования, где в качестве метрики взято мат. ожидание по состояниям системы. Существуют также и другие расширения данного алгоритма, где в качестве метрики взято значение Value at Risk или Expected Shortfall, широко используемые в оценки рисков финансовых потерь.

1.2.2 Прикладная область

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

1. Управление ГЭС

Моделируется следующая ситуация:

Существует водохранилище с ограниченным объемом воды. водохранилище наполняется за счет осадков, а опустошается с помощью шлюзов ГЭС. Также существует населенный пункт, потребляющий электроэнергию. ГЭС может генерировать электричество 2-мя способами: спускать воду из водохранилища или использовать топливные генераторы. Себестоимость электроэнергия, полученная ГЭС за счет спуска воды можно считать нулевой, в то время, как электроэнергия, выработанная за счет топливных генераторов имеет значительную себестоимость.

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

2. Управление портфелем ценных бумаг.

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

1.2.3 Алгоритм решения

Алгоритм решения задачи (1) - (3) опирается на фундаментальное утверждение:

cost-to-go function, опсанная уравнением (2) является кусочно-линейной относительно аргумента [math]x_{t-1}[/math]

Данное утверждение позволяет сделать вывод, что задача (1) - (3) может быть сведена к задаче линейного программирования.

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

1. Инициализаци

2. Обратный ход

3. Прямой ход

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

1.2.3.1 Инициализация.

На данном этапе генерируются случайные сценарии на сценарной решетки задачи. Количество инициализируемых случайных сценариев может варьироваться в зависимости от структуры сетки. До сих пор не исследовано оптимальное с точки зрения скорости выполнения количество пробных сценариев и зависимость выбора этого количества от размерности решетки. Будем считать, что на каждой итерации формируется [math]K[/math] случайных сценариев, которые представляют из себя наборы состояний [math](\widehat{\xi}_{1}^{i_k},...,\widehat{\xi}_{T}^{i_k}) \ k = 1,...,K[/math]

Для каждого пробного сценария на каждом шаге решается прямая задача линейного программирования:

[math] Q_t^{i_k}(\widehat{x}_{t-1}^{k},\widehat{\xi}_{t}^{i_k}) = \left\{\begin{matrix} min \ c_t^{i_k}x_t \\ subject \ to \\ A_t^{i_k}x_t + B_t^{i_k}\widehat{x}_{t-1}^k\geq b_t^{i_k} \end{matrix}\right. \ \ \ (2) [/math]

Получается [math]K[/math] пробных решений:

[math] \widehat{x}_t^k = \left\{\begin{matrix} argmin \ c_t^{i_k}x_t \\ subject \ to \\ A_t^{i_k}x_t + B_t^{i_k}\widehat{x}_{t-1}^k\geq b_t^{i_k} \end{matrix}\right. \ \ \ (3) [/math]

[math]k = 1,...,K[/math]

1.2.3.2 Обратный ход

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

Рис.2 Отсекающие гиперплоскости ограничивают истинную cost-to-go function снизу

Уравнение гиперплоскости получается в ходе решения прямой и обратной задач линейного программирования:

На шаге [math]T[/math]:

Прямые задачи:

[math] Q_T(\widehat{x}_{T-1}^k,\xi_T^i) = \left\{\begin{matrix} min \ c_T^{i_k}x_T \\ subject \ to \\ A_T^{i_k}x_T + B_T^{i_k}\widehat{x}_{T-1}\geq b_T^{i_k} \end{matrix}\right. (4) [/math]


Двойственные задачи:

[math] Q_T(\widehat{x}_{T-1}^k,\xi_T^i) = \left\{\begin{matrix} max \ \pi_T^{i_k}(b_t^i - B_T^i\widehat{x}_{T-1}^{i_k}) \\ subject \ to \\ \pi_T^{i_k}A_T^{i_k}\leq c_T^{i_k} \end{matrix}\right. (5) [/math]

Уравнение отсекающей гиперплоскости:

[math] l_T^k(x_{T-1}^k) = \sum_{j=1}^{m_T} \ p_{Tj}c_T^jx_t^{i_kj}+ (-\sum_{j=1}^{m_T}p_{tj}\pi_t^{i_kj}B_t^j)(x_{T-1}^{i_k} - \widehat{x}_{T-1}^{i_k}) \ \ (6) [/math]

На этапах [math]t = 1,...,T-1[/math]:

Прямые задачи:

[math] Q_t(x_{t-1}^k,\xi_t^i) = \left\{\begin{matrix} min \ c_t^{i_k}x_t + \gamma_t \\ subject \ to \\ A_t^{i_k}x_t + B_t^{i_k}\widehat{x}_{t-1}\geq b_t^{i_k} \\ \gamma_t\geq l_{t+1}(x_t^{i_k}) \end{matrix}\right. (7) [/math]


Двойственные задачи:

[math] Q_t(\widehat{x}_{t-1}^k,\xi_t^i) = \left\{\begin{matrix} max \ \pi_T^{i_k}(\widetilde{b}_t^i - \widetilde{B}_t^i\widehat{x}_{t-1}^{i_k}) \\ subject \ to \\ \pi_t^{i_k}\widetilde{A}_t^{i_k}\leq \widetilde{c}_t^{i_k} \end{matrix}\right. (8) [/math]

Уравнение отсекающей гиперплоскости:

[math] l_t^k(x_{t-1}^k) = \sum_{j=1}^{m_t} \ p_{tj}\widetilde{c}_t^jx_t^{i_kj}+ (-\sum_{j=1}^{m_t}p_{tj}\pi_t^{i_kj}\widetilde{B}_t^j)(x_{t-1}^{i_k} - \widehat{x}_{t-1}^{i_k}) \ \ (9) [/math]

В результате, значение функционала при оптимальном решении в уравнении (7) является нижней оценкой решения [math]\underline{s}[/math]

1.2.3.3 Прямой ход

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

На прямом ходе алгоритма генерируются пробные сценарии [math]\widehat{\xi}_t^k; \ k = 1,...,K[/math]

Для каждого из K сценариев решаются задачи (7), в результате получается верхняя оценка [math]\overline{s}[/math]

Если [math]|\overline{s} - \underline{s}| \leq \varepsilon[/math] то решение найдено. Иначе, переход на следующую итерацию.

1.3 Вычислительное ядро алгоритма

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

Рис.3 Типичное поведение верхней и нижней границы решения с увеличением итераций

1.4 Макроструктура алгоритма

Алгоритм SDDP является композицией задач линейного программирования. Пусть выполнены следующие условия:

[math]x = (x_1,...,x_N);[/math]

[math]\pi = (\pi_1,...,\pi_M);[/math]

[math]A \in R^{N \ X \ M};[/math]

[math]c \in R^N [/math]

Тогда, задачей линейного программирования называется задача следующего вида:

[math] \left\{\begin{matrix} min \ c^Tx \\ subject to \\ Ax \geq b \\ x_i \geq 0; \ i = 1,...,N \end{matrix}\right. (10) [/math]

Двойственной задачей к задаче (10) является задача следующего вида:

[math] \left\{\begin{matrix} max \ \pi^Tb \\ subject to \\ \pi^TA \leq c \\ \pi_i \geq 0; \ i = 1,...,M \end{matrix}\right. (11) [/math]

Задача (10) и задача (11) могут быть решены Симплекс методом

Обозначим основные свойства данных задач и Симплекс метода:

Пусть [math] M[/math] - количество ограничений, [math]N[/math] - количество переменных задачи линейного программирования. Тогда:

1. Практика показывает, что в среднем вычислительная сложность алгоритма составляет [math] M^3 [/math] операций

2. Количество итераций знаключение между [math] M [/math]и [math]2M [/math]

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

4. Количество ограничений влияет на время выполнения метода гораздо сильнее, чем количество переменных

5. Оптимальное значение функционала задачи (10) совпадает с оптимальным значением функционала (11)

6. Для решения алгоритма SDDP использовался готовый солвер для решения задач линейного программирования - lpSolve

1.5 Схема реализации последовательного алгоритма

Алгоритм SDDP начинается с инициализации пробных сценариев, данный этап описан в разделе 1.2.3.1

После генерации пробных сценариев и вычисления пробных решений, производится обратных ход алгоритма, описанный в 1.2.3.2

На данном этапе на каждом этапе формируются поддерживающие гиперплоскости, а также составляется нижняя оценка решения

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

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

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

Рис.4 Схема SDDP

1.6 Последовательная сложность алгоритма

Обозначим вычислительную решения задачи линейного программирования симплекс-методом с размерностью переменных оптимизации N и количеством ограничений M за [math]lp(N,M)[/math] Пусть задача SDDP имеет T этапов, размерность вектора переменных - N, размерность ограничений - M Количество состояний на каждом этапе - [math]S_t; \ t = 1,...,T[/math] Пусть также на каждой итерации генерируется k пробных сценариев Тогда для первой итерации алгоритма SDDP вычислительная сложность соответствует значению

[math] kTlp(N,M) + k(lp(N,M)+lp(M,N)) + k(T-2)(lp(N,M+k)+lp(M+k,N)) + lp(N,M+k) [/math]

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

На j-той итерации: [math] kTlp(N,M) + k(lp(N,M)+lp(M,N)) + k(T-2)(lp(N,M+k\sqrt{j})+lp(M+k\sqrt{j},N)) + lp(N,M+k\sqrt{j}) [/math]

Таким образом, сложность алгоритма при K итерациях составляет:

[math] P = \sum_{j=1}^{K} \ kTlp(N,M) + k(lp(N,M)+lp(M,N)) + k(T-2)(lp(N,M+k\sqrt{j})+lp(M+k\sqrt{j},N)) + lp(N,M+k\sqrt{j}) [/math]


Вообще говоря на решетке с такими данными возможно [math]\prod_{t=1}^TS_t[/math] сценариев. Также, в худшем случае, в ходе решения не была удалена ни одна отсекающая гиперплоскость Тогда в худшем случае, решение будет достигнуто с вычислительной сложностью:

[math] P = \sum_{j=1}^{\frac{\prod_{t=1}^TS_t}{k}} \ kTlp(N,M) + k(lp(N,M)+lp(M,N)) + k(T-2)(lp(N,M+kj)+lp(M+kj,N)) + lp(N,M+kj) [/math]

Так как средняя вычислительная сложность симплекс метода [math]lp(N,M) = M^3[/math], вычислительная сложность алгоритма SDDP в худшем случае составляет

[math] P = \sum_{j=1}^{\frac{\prod_{t=1}^TS_t}{k}} \ kTM^3 + k(M^3+N^3) + k(T-2)((M+kj)^3+N^3) + (M+kj)^3 [/math]

В терминах эквивалентности:

[math] P \sim O(\overline{S}^{2T}*(M^3+N^3)) [/math]

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

1.7 Информационный граф

Информационный граф 1 итерации алгоритма.

На Рис.5 изображен обратный ход алгоритма, проход от последних состояний к корневым

На Рис.6 изображен прямой ход алгоритма по сценарию, выделенному зеленым

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

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

рис.5 Информационный граф 1 итерации алгоритма SDDP (обратный ход)


рис.6 Информационный граф 1 итерации алгоритма SDDP (прямой ход)

1.8 Ресурс параллелизма алгоритма

При параллельной реализации алгоритма выгодно с точки зрения быстродействия увеличивать количество пробных сценариев на каждой итерации, так как это уменьшает количество необходимых итераций алгоритма. При условии, что в расположении имеется вычислительный кластер бесконечного размера, решение задачи можно гарантированно получить за 1 итерацию. Для этого потребуется [math] \prod_{t=1}^T \ S_t [/math] вычислительных узлов 1-го уровня, каждый из которых будет рассчитывать коэффициенты отсекающей гиперплоскости для каждого сценария. Внутри каждого узла производится последовательное, от этапа к этапу, вычисление прямых и двойственных задач линейного программирования. В рамках данной статьи не рассматривается рессурс параллельного алгоритма симплекс метода, к тому же lpSolve доступен только в скалярной реализации. Поэтому, в рамках данной статьи будем считать, что решение задачи линейного программирования - задача для одного процессора. С учетом таких допущений, параллельная вычислительная мощность может составить [math]P \sim O(max(M^3,N^3)T)[/math]

1.9 Входные и выходные данные алгоритма

Параметры задачи:

[math] T \in \mathbb{N} [/math] - количество этапов задачи

[math] N \in \mathbb{N} [/math] - размерность вектора переменных задачи

[math] S \in \mathbb{N}^T [/math] - вектор количества состояний на каждом этапе

Входные данные:

[math]p[/math] - вероятности перехода от состояния на предыдущем этапе в состояние на следующем

[math] p_t \in [0,1]^{S_t \times S_{t+1}}; \ t = 1,...,T [/math]

[math] \sum_{i=1}^{S_{t+1}} \ p_t^{j \rightarrow i} = 1; \ \forall j = 1,...,S_t [/math]

[math]\xi[/math] - состояние

Каждое состояние описывается набором данных

[math] \xi_t^i = \{A_t^i,B_t^i,c_t^i,b_t^i\}; \ t=1,...,T; \ i = 1,...,S_t[/math]

где:

[math] A_t^i \in \mathbb{R}^{N\times M}[/math]

[math] B_t^i \in \mathbb{R}^{N\times M}[/math]

[math] c_t^i \in \mathbb{R}^{N}[/math]

[math] b_t^i \in \mathbb{R}^{M}[/math]

SDDP применяется для задач с количеством этапов порядка 100, количеством состояний порядка 20 и размерностью решения и ограничений порядка 50 на 50. Таким образом, в памяти на протяжении всей программы должно храниться 100*20*(2*50+2*50*50) = 1 000 000 000 чисел типа double, т.е. около 8 гб.

Выходные данные алгоритма:

[math]X \in \mathbb{R}^{N\times T}[/math] - матрица решения.

1.10 Свойства алгоритма

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

2. Вычислительная мощность алгоритма в худшем случае: [math] P \sim O(\overline{S}^{2T}*(M^3+N^3)) [/math]

3. Вычислительная мощность алгоритма в среднем: [math] P \sim O(\overline{S}*2T^3*(M^3+N^3)) [/math]

2 Программная реализация алгоритма

2.1 Особенности реализации последовательного алгоритма

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

2.2 Локальность данных и вычислений

На практике алгоритм SDDP необходим для решения задачи с количеством этапов порядка 100, количеством состояний порядка 20 и размерностью решения и ограничений порядка 50 на 50. Таким образом, в памяти на протяжении всей программы должно храниться 100*20*(2*50+2*50*50) = 1 000 000 000 чисел типа double, т.е. около 8 гб. Также, каждому состоянию в ходе программы привязывается набор коэффициентов гиперплоскости, с каждой итерацией добавляются новые наборы. Так, на 100 - итерации, в памяти должно храниться около 100*20*51*100 = 10 200 000 т.е. около 80мб. Обращение к входным данным имеет временную неравномерность. Если количество пробных сценариев = k, то один раз за итерацию к каждой матрице каждого этапа алгоритм обратится k раз, причем почти одновременно. Таким образом, можно ускорить выполнение алгоритма в условиях недостаточной быстрой памяти с помощью загрузки данных нескольких этапов в кэш. При параллельной реализации можно предоставить памяти, выделенной под исходные данные, общий доступ, так как запись в эту область не предусматривается.

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

2.3 Возможные способы и особенности параллельной реализации алгоритма

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

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

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

2.4 Масштабируемость алгоритма и его реализации

Ввиду того, что существуют реализации алгоритма только на языке R и Matlab, которые не установлены ни на Blue gene (установить невозможно), ни на Lomonosov(в процессе установки) Результаты расчетов были проведены лишь для количества задействованных ядер многоядерного домашнего компьютера. Параметры запуска:

Процессор Intel Core i5 (4 ядра)

Программа написана на языке R, программа запуска - R.exe (R - не компилируемый, а интерпритируемый язык, изначально приложение R.exe написано на С и было откомпилировано для операционной системы Windows до установки на ПК)

Версия исполняемого скрипта - мой собственный код, который я пишу в рамках магистерской диссертации.

Используемые библиотеки: lpSolve, Rmpi, snow, Rcluster

Результаты тестов:

Размерность данных (T;S;N;M) 1 2 4
(5;5;10;10) 279 159,064 86,24
(10;10;10;10) 1120,8 608,402 323,239
(28;59;6;6) 62344363,245 Na Na

Выводы о масштабируемости:

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

2.5 Динамические характеристики и эффективность реализации алгоритма

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

2.6 Выводы для классов архитектур

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

2.7 Существующие реализации алгоритма

На данный момент существует пакет Fast в матлабе, который позволяет запускать алгоритм SDDP с мат. ожиданием в качестве метрики на скалярном компьютере. Также в разработке находится проект, написанный на R, позволяющий использовать в качестве метрики как мат.ожидание, так и VaR, так и CVaR, разрабатывается параллельный алгоритм вычисления с помощью пакета Rmpi.

3 Литература

1. M.V.F. Pereira and L.M.V.G. Pinto "Multi-stage stochastic optimization applied to energy planning" Mathematical Programming 52 (1991) 359-375

2. Alexander Shapiro, Darinka Dentcheva, Andrzej Ruszczy´nski "Lectures on Stochastic Programming: Modeling and Theory"

Источник — https://algowiki-project.org/w/ru/index.php?title=VladimirDobrovolsky611/Алгоритм_SDDP&oldid=22470