Стохастическое двойственное динамическое программирование (SDDP)

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

SDDP
Последовательный алгоритм
Последовательная сложность	[math]O(\overline{S}^{2T}\cdot(M^3+N^3))[/math]
Объём входных данных	размерность задачи [math]N[/math], количество ограничений [math]M[/math], количество этапов [math]T[/math], среднее количество состояний [math]\overline{S}[/math]

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

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

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

Алгоритм SDDP (Stochastic Dual Dynamic Programming) впервые был предложен в 1991 году в статье . Далее алгоритм претерпел множество модернизаций и спецификаций, описанных в труде . В настоящее время под аббревиатурой SDDP подразумевается целое семейство алгоритмов.

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

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

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

1. Количество этапов T, количество состояний на каждом этапе [math]m_t, \ 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-м этапе задача детерминирована. Здесь присутствует всего одно состояние [math] \xi_1 =(A_1, \ b_1, \ c_1) [/math].

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

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

[math] \left\{\begin{matrix} \min \ c_1x_1 + \sum_{i=1}^{m_2}p_{1i}^1Q_2(x_1,\xi_i^2) \\{\rm subject \ to} \\A_1x_1\geqslant b_1 \\x_1 \geqslant 0 \end{matrix}\right. \qquad (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}) \\{\rm 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 \qquad (2) [/math]

...

[math] Q_T(x_{T-1},\xi_i^T) = \left\{\begin{matrix} \min \ c_T^ix_T \\{\rm subject \ to} \\A_T^ix_T + B_T^ix_{T-1}\geqslant b_T^i \\x_t \geqslant 0 \end{matrix}\right. \qquad (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 \\ {\rm 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. \quad (2) [/math]

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

[math] \widehat{x}_t^k = \left\{\begin{matrix} {\rm argmin} \ c_t^{i_k}x_t \\ {\rm 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. \quad (3) [/math]

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

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

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

Рис.1. Отсекающие гиперплоскости ограничивают истинную 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 \\ {\rm subject \ to} \\ A_T^{i_k}x_T + B_T^{i_k}\widehat{x}_{T-1}\geq b_T^{i_k} \end{matrix}\right. \qquad (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}) \\ {\rm subject \ to} \\ \pi_T^{i_k}A_T^{i_k}\leq c_T^{i_k} \end{matrix}\right. \qquad (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}) \qquad (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 \\ {\rm 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. \qquad (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}) \\ {\rm subject \ to} \\ \pi_t^{i_k}\widetilde{A}_t^{i_k}\leq \widetilde{c}_t^{i_k} \end{matrix}\right. \qquad (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}) \qquad (9) [/math]

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

1.2.3.3 Прямой ход

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

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

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

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

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

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

1. Вычисление прямой или двойственной задачи линейного программирования. Получения значения функционала и оптимального решения.

2. Запись в пространство гиперплоскостей коэффициентов задачи.

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

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

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

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

[math]x = (x_1,...,x_N), \quad \pi = (\pi_1,...,\pi_M), \quad A \in R^{N \ X \ M}, \quad c \in R^N[/math].

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

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

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

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

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

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

Пусть [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, вычисляется верхняя оценка решения.

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

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

Рис.3. Схема алгоритма 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}^{{\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}^{{\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}\cdot(M^3+N^3)). [/math]

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

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

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

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

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

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

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

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

При параллельной реализации алгоритма выгодно с точки зрения быстродействия увеличивать количество пробных сценариев на каждой итерации, так как это уменьшает количество необходимых итераций алгоритма. При условии, что в расположении имеется вычислительный кластер бесконечного размера, решение задачи можно гарантированно получить за 1 итерацию. Для этого потребуется [math]\textstyle\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}}, \quad t = 1,...,T [/math]

[math] \sum_{i=1}^{S_{t+1}} \ p_t^{j \rightarrow i} = 1, \quad \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\}, \quad t=1,...,T, \quad i = 1,...,S_t[/math]

где

[math] A_t^i \in \mathbb{R}^{N\times M}, \quad B_t^i \in \mathbb{R}^{N\times M}, \quad c_t^i \in \mathbb{R}^{N}, \quad b_t^i \in \mathbb{R}^{M}[/math].

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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