Sequential Decision Analytics and Modeling 第2版
Back to SDA site →

第8章:储能 I

本章概览

本章考察一个初看起来相当简单的库存问题,该问题出现在存储可从电网购买或向电网出售的能源时,而电网价格具有高度随机性。与前六章不同,我们使用一套更为丰富的模型来描述这些随机价格,这开始揭示出我们在对不确定性建模时可能遇到的复杂性。我们将概览用于价格过程的不同模型,包括经典时间序列模型、跳跃扩散模型(用于捕捉价格尖峰)、分位数分布,以及最后一种将分位数分布与标准正态分布相结合的混合模型,从而为使用依赖于正态性的方法打开了大门。

接着我们描述一系列策略,从一个基本的”低买高卖”策略(一种策略函数近似形式)开始,然后过渡到几种基于近似贝尔曼方程的方法,我们在第5章中首次见到了贝尔曼方程。我们从贝尔曼方程的基础形式开始(这对几乎所有问题在计算上都是不可行的),然后概览被称为反向近似动态规划(ADP)、正向ADP以及结合参数调优的正向ADP混合策略等各种变体。

叙述

新泽西州正寻求开发3,500兆瓦(MW)的离岸风力发电能力。一个挑战是风能(尤其是离岸风能)可能具有高度可变性。风能对电网的这种可变性影响,由于风能(在中间范围内)随风速的三次方增长这一特性而被放大。图8.1描绘了这种可变性。

来自五种不同风力发电容量水平的电力。
图8.1。 来自五种不同风力发电容量水平的电力。

风能在风力资源丰富的地区已变得越来越受欢迎,例如美国中西部、欧洲沿海地区、巴西东北部以及中国北方地区(仅举几例)。有时社区(和公司)投资可再生能源(风能或太阳能)是为了帮助减少其碳足迹,并尽量减少对电网的依赖。

然而,这些项目很少能让社区完全摆脱电网。通常的做法是让可再生能源(风能或太阳能)直接向电网出售电力,而公司则可以从电网购买电力。这可以作为一种对冲手段发挥作用,因为公司会在价格飙升期间(价格可能从每兆瓦时(mwh)20美元跃升至300美元甚至更高)赚取大量利润,从而抵消在这些时期购买电力的成本。

可再生能源面临的主要困难是应对其可变性。虽然一种解决方案是简单地将可再生能源产生的所有能量输入电网,并利用电网的容量来应对这种可变性,但人们一直对使用储能(特别是电池储能)来削峰填谷抱有浓厚兴趣。除了平滑可再生能源的可变性外,人们还对利用电池来利用价格尖峰产生兴趣,即在电价低廉时购买电力(价格甚至可能变为负值),并在电价高时卖出。利用电网电价的可变性在价格低时买入、价格高时卖出,被称为电池套利。

用于电力稳定和电池套利的电网到储能系统。
图8.2。 用于电力稳定和电池套利的电网到储能系统。

我们将使用图8.2所示的配置来说明储能中的若干建模和算法问题。该问题将为几乎任何库存/储存问题提供洞见,包括:

储能是一种特别丰富的库存问题形式。虽然我们不会考虑所有可能的变体(这些变体数不胜数),但我们的问题将呈现以下特征:

构建问题框架

我们对三个框架问题的回答是:

基本模型

我们将使用该问题的一系列变体来说明不同的建模问题,首先从使用电池向电网买卖电力以利用价格波动的基本系统开始。对于该应用,我们将以5分钟为时间增量向前推进,因为这是电网价格更新的频率(该时间增量因电网运营商而异)。

状态变量

对于我们的基本模型,我们只需追踪两个变量:$R_t$,即时刻$t$电池中储存的能量数量(以兆瓦时,即mwh计量);以及$p_t$,即电网上的电价。我们的状态变量则为

\[S_t = (R_t, p_t).\]

随着我们向模型中添加不同要素,状态变量会迅速变得更加复杂。例如,表示价格的状态变量取决于我们如何对价格过程建模,如下文转移函数中所述。

决策变量

我们唯一的决策是从电网购买还是向电网出售:$x_t$,即从电网购买($x_t > 0$)或向电网出售($x_t < 0$)的电力数量。

当我们将能量转入或转出电池时,我们将假设在转移过程中只能获得一部分$\eta$,这意味着有$1-\eta$的损失。为简化起见,我们将假设无论是给电池充电还是放电,这种损失都是相同的。

该决策受电池容量的限制,这意味着我们必须遵守以下约束

\[x_t \leq \frac{1}{\eta} (R^{max} - R_t), \qquad x_t \geq -\eta R_t,\]

其中第一个约束适用于我们从电网购买时($x_t > 0$),而第二个约束适用于我们向电网出售时($x_t < 0$)。

与往常一样,我们假设决策是通过策略$X^\pi(S_t)$做出的,具体将在下文确定。

外源信息

在我们的基本模型中,唯一的外源信息是价格的变化。我们可以假设每个时间段的价格都会被揭示出来,而不需要任何基于历史价格来预测价格的模型。在这种情况下,我们的外源信息$W_t$将是

\[W_{t+1} = p_{t+1}.\]

或者,我们可以假设观察到的是价格变化$\phat_t = p_t - p_{t-1}$,此时我们会写作

\[W_{t+1} = \phat_{t+1}.\]

转移函数

状态变量的演化由下式给出

\[\begin{align} R_{t+1} &= \begin{cases} R_t + \eta x_t & x_t \geq 0, \\ R_t + \dfrac{x_t}{\eta} & x_t < 0. \end{cases} \label{eq:energytransition1}\\ p_{t+1} &= p_t + \phat_{t+1}. \label{eq:energytransition2} \end{align}\]

这种通过”观察”价格变化来对价格过程建模的方式有助于我们编写转移函数。在实践中,我们通常会直接观察到$p_{t+1}$(而不是其变化量),在这种情况下就不需要显式的转移方程。这两个方程构成了我们的转移函数$S_{t+1} = S^M(S_t,x_t,W_{t+1})$。

稍后,我们会发现建模决策后状态$S^x_t$很有用处,这是我们做出决策$x_t$之后、但新信息到达之前的状态。决策后资源状态为

\[R^x_t = \begin{cases} R_t + \eta x_t & x_t \geq 0, \\ R_t + \dfrac{x_t}{\eta} & x_t < 0. \end{cases}\]

由于储存变量的演化是确定性的,因此向下一个决策前状态的转移就是

\[R_{t+1} = R^x_t.\]

另一方面,价格$p_t$不受决策影响,因此决策后价格就是

\[p^x_t = p_t.\]

这意味着决策后状态为

\[S^x_t = (R^x_t, p_t).\]

目标函数

在任何一个时间段内,我们赚取或损失的金额由下式给出

\[C(S_t,x_t) = -p_t x_t.\]

我们的目标函数则是由下式给出的典范问题

\[\max_\pi \E \sum_{t=0}^T -p_t X^\pi(S_t),\]

其中$S_{t+1} = S^M(S_t,X^\pi(S_t),W_{t+1})$由方程$\eqref{eq:energytransition1}$和$\eqref{eq:energytransition2}$给出。我们随后还需要指定初始状态$S_0$(即$R_0$和$p_0$),并且需要一种生成$W_1, W_2, \ldots$的方法,我们将在下文对此加以描述。

对不确定性建模

第2章中的资产出售问题里,我们假设可以按照下式对价格建模

\[p_{t+1} = p_t + \varepsilon_{t+1},\]

然后我们假设$\varepsilon_{t+1}$服从均值为0、方差已知的正态分布。图8.3展示了一年内电网价格的走势,这类价格在电力行业术语中被称为”节点边际价格”(或LMP),它揭示了电网价格所表现出的巨大波动性。这种波动性源于负荷激增(或电力损失)可能造成短期短缺。由于需求缺乏弹性(电网被要求满足100%的负荷),价格可能在短时间内(价格以5分钟为增量更新)跃升20到50倍。

2010年PJM电网的节点边际价格(以5分钟为间隔)。
图8.3。 2010年PJM电网的节点边际价格(以5分钟为间隔)。

有多种方法可用于对电价建模。下面我们将描述已被用于该问题的四种方法。

时间序列模型

时间序列文献相当丰富,因此我们仅将说明一个基本模型,该模型将价格$p_{t+1}$表示为近期历史价格的函数。为便于说明,我们将使用最近三个时间段,这意味着我们的模型可以写作

\[\begin{align} p_{t+1} &= \thetabar_{t0} p_t + \thetabar_{t1} p_{t-1} + \thetabar_{t2} p_{t-2} + \varepsilon_{t+1}, \label{eq:energytimeseriesmodel}\\ &= \thetabar^T_t \phi_t + \varepsilon_{t+1}, \nonumber \end{align}\]

其中

\[\phi_t = \begin{pmatrix} p_t \\ p_{t-1} \\ p_{t-2} \end{pmatrix}\]

是我们的价格向量。我们假设对于给定的$\sigma^2_\epsilon$,噪声$\varepsilon \sim N(0,\sigma^2_\epsilon)$成立。

系数向量$\thetabar_t = (\thetabar_{t0},\thetabar_{t1},\thetabar_{t2})^T$可以递归地估计。假设我们从系数向量的初始估计$\thetabar_0$开始。我们还需要一个三乘三矩阵$M_0$,目前可以假设它是一个按比例缩放的单位矩阵(我们将在下文提供更好的想法)。

$\thetabar_t$的基本更新方程由下式给出

\[\thetabar_{t+1} = \thetabar_t - H_t\phi_t \varepsilon_{t+1},\]

误差$\hat{\varepsilon}_t$使用下式计算

\[\varepsilon_{t+1} = \thetabar^T_{t}\phi_t - p_{t+1}.\]

三乘三矩阵$H_t$使用下式计算

\[H_t=\frac{1}{\gamma_t}M_t,\]

其中矩阵 $M_t$ 使用以下递归方式计算

\[M_t = M_{t-1} - \frac{1}{\gamma_t} (M_{t-1} \phi_t (\phi_t)^T M_{t-1}).\]

变量 $\gamma_t$ 是使用以下公式计算的标量

\[\gamma_t = 1 + (\phi_t)^TM_{t-1}\phi_t.\]

这些方程需要 $\thetabar_0$ 和 $M_0$ 的初始估计值。一种方法是先收集一些初始数据,然后求解一个静态估计问题。假设我们观测到 $K$ 个价格。设 $Y_0$ 为观测价格 $p_3, p_4, \ldots, p_{K+3-1}$ 组成的 $K$ 元列向量(我们必须从第三个价格开始,因为模型需要用到之前三个价格)。

然后设 $X_0$ 为一个具有 $K$ 行的矩阵,其中每一行由 $p_k, p_{k-1}, p_{k-2}$ 组成。我们对 $\thetabar$ 的最佳估计由正规方程给出

\[\thetabar_0 = [(X_0)^T X_0]^{-1} (X_0)^T Y_0.\]

最后设 $M_0 = [(X_0)^T X_0]^{-1}$,这表明矩阵 $M_t$ 是 $[(X_t)^T X_t]^{-1}$ 在时间 $t$ 的估计值。

存在整整一族用于捕捉变量随时间关系的时间序列模型。如果我们直接将这些方法应用于价格数据,效果会相当差。首先,价格并非正态分布。其次,虽然价格可能变为负值,但这种情况相当罕见。然而,如果将方差 $\sigma^2_\epsilon$ 校准为这类数据的高噪声水平,直接应用该模型极有可能产生负价格。最后,随时间推移的价格跳变行为也不会真实。

跳跃扩散

对上述线性模型的一个主要批评是,它无法很好地捕捉电力价格研究中常见的大幅尖峰。克服这一局限性的一个简单思路是使用所谓的跳跃扩散模型(jump diffusion model),即在方程 $\eqref{eq:energytimeseriesmodel}$ 中加入另一个噪声项,得到

\[p_{t+1} = \thetabar_{t0} p_t + \thetabar_{t1} p_{t-1} + \thetabar_{t2} p_{t-2} + \varepsilon_{t+1} + \mathbb{I}_t \varepsilon^J_{t+1}.\]

这里,指示变量 $\mathbb{I}_t = 1$ 以某个概率 $p^{jump}$ 取值,噪声 $\varepsilon^J_{t+1}$ 服从均值为 $\mu^{jump}$(通常远大于零)、方差为 $(\sigma^{jump})^2$(相当大)的正态分布。

我们需要估计跳跃概率 $p^{jump}$,以及均值和方差 $(\mu^{jump}, (\sigma^{jump})^2)$。这可以通过先建立一个基本模型(其中 $p^{jump} = 0$)来实现。我们利用这个基本模型来估计 $\sigma^2_\epsilon$。然后我们选定某个容差范围,例如三个标准差(即 $3 \sigma_\epsilon$),任何超出此范围的观测值都被认为是由另一种噪声来源引起的。设 $p^{jump}$ 为出现此类观测值的时间段所占的比例。然后,计算这些观测值的均值和标准差,得到 $(\mu^{jump}, (\sigma^{jump})^2)$。

我们不会止步于此。在从数据中剔除这些极端变化之后,我们应该在不包含这些观测值的情况下重新拟合我们的线性模型。标准做法是重复这一过程若干次,直到这些估计值不再变化。

跳跃扩散模型在复现尾部行为方面表现更好,但它仍然依赖于正态分布的尾部行为。通过认识到噪声的方差取决于温度,尤其是极端温度,可以获得更好的拟合效果。我们可以将温度划分为三个区间:低于冰点、高于90华氏度,以及介于这两者之间。引入温度依赖性,虽然会给状态变量集合增加一个变量,但也会引入额外的复杂度(这种影响取决于所采用的策略类别)。

分位数分布

尽管拟合其他参数分布也是可能的,但一种强有力的策略是根据数据数值化地计算累积分布,从而构建通常所称的分位数分布(quantile distribution)。为计算此分布,我们只需将价格从小到大排序。将该有序序列记为 $\ptilde_t$,其中 $\ptilde_{t-1} \leq \ptilde_t$。设 $T = 105,210$ 为一年中5分钟时间段的数量。价格小于 $\ptilde_t$ 的时间段所占的百分比即为 $t/T$。我们可以使用以下方法建立累积分布

\[F_P(\ptilde_t) = \frac{t}{T},\]

如图 8.4 所示。

价格的分位数分布。
图 8.4。 价格的分位数分布。

我们可以通过找到最大的 $\ptilde_t < p$ 并将 $F_P(p)$ 设为该值,从而为任意 $p$ 创建一个连续分布 $F_P(p)$,形成一个阶梯函数。函数 $F_P(p)$ 是一种非参数(nonparametric)分布形式,因为我们并未将该分布拟合到任何已知的参数形式。好消息是,它将与数据完美匹配,这意味着我们能够准确表示电力价格中出现的极端尾部。坏处是,我们需要一个良好的数据集来创建这些分布,并且必须保留该数据集以便计算分布,而不能像拟合参数模型那样仅存储少量参数。

我们可以通过生成一个在0到1之间均匀分布的随机变量 $U$ 从该分布中抽样。假设我们生成了 $U= 0.70$。然后我们想要找到对应于 $F_P(p^{.70}) = 0.70$ 的价格 $p^{.70}$,如图 8.4 所示。我们通过定义逆函数 $F^{-1}_P(u)$ 来在数学上表达这一过程,该函数返回产生 $F_P(p) = u$ 的价格 $p$。我们可以通过重复抽样均匀随机变量 $U$,然后观察对应的价格 $p=F^{-1}_P(U)$,从而反复从价格分布中抽样。

结合数据变换的混合时间序列方法

一种强有力的策略是将经验分布的使用与经典时间序列方法结合起来。我们首先对价格数据拟合一个经验分布,得到累积分布 $F_P(p)$。现在,设 $p_t$ 为一个价格,计算 $u_t = F_P(p_t)$,其中 $0\leq u_t \leq 1$。接下来,设 $\Phi(z)$ 为标准正态随机变量 $Z \sim N(0,1)$ 的累积分布,并设 $\Phi^{-1}(u)$ 为其逆函数。接着设 $z_t = \Phi^{-1}(u_t)$。映射 $p_t \rightarrow u_t \rightarrow z_t$ 的过程如图 8.5 所示。

将经验分布转换为正态分布(并再转换回来)。
图 8.5。 将经验分布转换为正态分布(并再转换回来)。

我们可以利用此方法将高度非正态的价格 $p_t$ 转换为均值为0、方差为1的正态分布数值序列 $z_t$,同时我们还可以捕捉相关性。然后,我们可以对序列 $z_t$ 执行任意时间序列建模。之后,从该标准化模型得出的任何估计值都可以通过沿图 8.5 的路径反向追踪转换回价格:$u_t = \Phi(z_t)$,然后 $p_t = F^{-1}_P(u_t)$。

在处理非正态分布数据时,这种策略非常有效,其表现远优于金融领域中流行的跳跃扩散模型。

设计策略

我们将通过两类策略再加一种混合策略来说明如何求解该问题:

基于贝尔曼方程的策略需要计算(或近似)值 $V_{t+1}(S_{t+1})$,该值是处于状态 $S_t$、采取决策 $x_t$、然后观测到随机外源信息 $W_{t+1}$ 所产生的结果。我们最早是在最短路径问题的背景下见到这些方法的。现在一个重要的区别是,给定 $S_t$ 和 $x_t$ 时状态 $S_{t+1}$ 是随机的(在最短路径问题中,只有成本 $\chat_t$ 是随机的)。此外,我们的状态变量现在具有两个连续维度,而不仅仅是离散节点。

我们先描述低买高卖策略,然后介绍三种基于近似贝尔曼方程的方法:

最后,我们将描述一种混合策略,将来自贝尔曼方程的值函数逼近与某种形式的策略搜索结合起来。

低买高卖

低买高卖策略遵循一个简单原则:当价格低于某个下限时为电池充电,当价格高于某个上限时卖出。该策略可以写成

\[X^{low-high}(S_t\vert \theta) = \begin{cases} -1 & \text{if } p_t \leq \theta^{buy}, \\ 0 & \text{if } \theta^{buy} < p_t < \theta^{sell}, \\ +1 & \text{if } p_t \geq \theta^{sell}. \end{cases}\]

现在我们需要调整 $\theta = (\theta^{buy}, \theta^{sell})$。我们通过沿着一条价格样本路径 $p_t(\omega)$ 来评估我们的策略(或者我们观测到的可能是价格的变化 $\phat(\omega)$)。假设我们是从某个数学模型生成样本路径,我们可以生成样本路径 $\omega^1, \ldots, \omega^N$。然后我们可以在每条样本路径上模拟该策略的表现,并使用以下方式取平均值

\[\Fbar^{low-high} = \frac{1}{N} \sum_{n=1}^N C\big(S_t(\omega^n),X^{low-high}(S_t(\omega^n)\vert \theta)\big).\]

调整 $\theta$ 需要求解以下问题

\[\begin{align} \max_\theta \Fbar^{low-high}(\theta\vert S_0). \label{eq:buylowpolicysearch} \end{align}\]

由于 $\theta$ 仅具有两个维度,一种策略是通过对每个维度进行离散化,然后搜索这两个维度所有可能的取值组合,从而进行全面的网格搜索。一种常见的离散化方式是将区域划分为5%的增量。包含边界值的情况下,这意味着我们必须表示每个参数的21个取值,形成一个大小为441个点的网格,对大多数问题来说这是可以处理的(尽管并非易事)。

我们注意到,暴力网格搜索只有在我们运行足够多的模拟次数 $N$ 使得估计值 $\Fbar^\pi(\theta)$ 的方差相对较小时才可行。然而,即使在策略性能估计存在噪声的情况下,我们也有方法来搜索 $\theta$,这将在第7章中介绍。

我们将反复看到由 $\eqref{eq:buylowpolicysearch}$ 给出的优化问题,因为最简单的策略总是具有可调参数。如果我们能够计算(或近似)$\Fbar^{low-high}(\theta\vert S_0)$ 关于 $\theta$ 的导数,那么问题 $\eqref{eq:buylowpolicysearch}$ 可以使用基于导数的方法来求解。当这不可行时,我们必须使用无导数方法,这正是我们在第2章中面对的问题。

逆向动态规划

逆向动态规划涉及直接求解贝尔曼方程

\[\begin{align} V_t(s_t) = \max_{x_t} \left(C_t(s,x_t)+ \E\{V_{t+1}(S_{t+1})\vert S_t,x_t\} \right), \label{eq:energystoragebellman} \end{align}\]

其中 $S_{t+1} = S^M(s_t,x_t,W_{t+1})$,期望是对随机变量 $W_{t+1}$ 取的。假设 $W_{t+1}$ 是离散的,取值于 $\Wcal = \lbrace w_1, w_2, \ldots, W_M\rbrace $,并使用以下方式表示概率分布

\[f^W(w\vert s_t,x_t) = Prob[W_{t+1} = w\vert s_t,x_t].\]

我们将该分布写成依赖于状态 $s_t$ 和决策 $x_t$ 的形式,但这取决于具体问题。例如,我们可能合理地假设价格变化 $p_{t+1}-p_t$ 取决于当前价格 $p_t$(如果价格非常高,它们更有可能下跌),这就是以 $s_t$ 为条件的理由。如果从电网购买大量电力会推高价格,我们甚至可能需要考虑对 $x_t$ 的依赖性。

然后我们可以将方程 $\eqref{eq:energystoragebellman}$ 重写为

\[V_t(s_t) = \max_{x_t} \left(C_t(s_t,x_t)+ \sum_{w\in\Wcal} f^W(w\vert s_t,x_t)V_{t+1}(S^M(s_t,x_t,w)) \right).\]

逆向动态规划的一个朴素实现表现出四层循环:

  1. 从时间 $T$ 到时间 $0$ 向后逐步递推的循环。
  2. 遍历所有可能状态 $s_t\in\Scal$ 的循环(更准确地说,这是时间 $t$ 时状态 $S_t$ 的所有可能取值的集合)。
  3. 为求解最大化问题而需要遍历所有可能决策 $x_t$ 的循环。
  4. 遍历计算 $V_t(s)$ 所需求和中包含的随机变量 $W$ 所有可能取值的循环。

后向动态规划

第0步. 初始化:对所有状态$S_{t+1}$初始化终端贡献$V_{T+1}(S_{T+1})=0$。

第1步. 对$t=T, T-1, \ldots, 1, 0$执行:

第2步. 对所有$s\in\Scal$,计算

$$V_t(s_t) = \max_{x_t} \left(C_t(s_t,x_t)+ \sum_{w\in\Wcal} f^W(w\vert s_t,x_t)V_{t+1}(S^M(s_t,x_t,w)) \right)$$

考虑每个循环可能取值的范围是有益的。对于能源问题,我们可能以小时为增量在一天内优化一个储能设备,这就给出了24个时间步。如果我们使用5分钟的时间步(一些电网运营商每5分钟更新一次价格),那么24小时的时域将意味着288个时间段(如果我们想在一周内进行规划,则要乘以七)。如果我们在做调频,那么我们必须每2秒做一次决策,这相当于一天内有43,200个时间段。

我们的状态变量由$S_t = (R_t,p_t)$组成,这意味着我们必须用对所有$R_t$取值以及所有$p_t$取值的嵌套循环,来替代对所有状态的循环。由于两者都是连续的,每一个都必须被离散化。资源变量$R_t$必须根据我们在单个时间增量内可能充电或放电的量来划分为若干增量。然后我们必须离散化电网价格$p_t$。电网价格可以低至-$100,也可以高达$10,000(在极端情况下)。一个合理的策略可能是构建一个经验分布,然后表示与累积分布每2%的增量对应的价格,这样就给出了50个可能的价格。

充放电决策的数量可能小到只有三个(充电、放电或不操作),如果我们可以以不同的速率充电或放电,则数量可能会大得多。

最后,概率分布$f^W(w)$将是价格随机变化的分布$\phat_{t+1}$。同样,我们建议构建$\phat_{t+1}$变化的经验分布,然后将累积分布离散化为,比如说,2%的增量。

如果我们有一个二维状态变量(正如我们基本模型的情形),那么我们已经有了五个循环(时间、两个状态变量、对$x$的max算子,以及对$W$结果的求和)。这可能会变得代价高昂,而我们才刚刚开始。现在设想我们使用方程$\eqref{eq:energytimeseriesmodel}$中的时间序列模型,此时我们还必须追踪价格$(p_t, p_{t-1}, p_{t-2})$。在这种情况下,我们的状态变量将是

\[S_t = (R_t, p_t, p_{t-1}, p_{t-2}).\]

在这种情况下,我们现在将有七个嵌套循环。虽然复杂度取决于连续变量的离散化程度,运行这个后向动态规划算法很容易需要一年(甚至更长)的时间。

鉴于即使对于这个相对简单的问题,使用贝尔曼方程也如此困难,令人惊讶的是这种具体方法仍在课堂上被教授。鉴于这种复杂性,已经有大量研究致力于近似贝尔曼方程的方法,这些方法被归类为近似动态规划强化学习等名称。我们将描述两种近似贝尔曼方程的策略,称为后向近似动态规划(backward ADP)和前向近似动态规划(forward ADP)。

后向近似动态规划

一个强大的算法策略被称为”后向近似动态规划”。这种方法的进展与我们上面刚做的完全一样,只有一处不同。我们不是遍历所有状态,而是选择一个随机样本$\Shat$。然后我们像最初那样计算处于状态$s\in\Shat$的价值,并计算相应的值$\vhat$。假设我们重复此过程$N$次,并获得一个数据集$(\shat^n, \vhat^n), n=1, \ldots, N$。然后我们用它来拟合一个统计模型,比如由以下给出的线性模型:

\[\begin{align} \Vbar(s) = \theta_0 + \theta_1 \phi_1(s) + \theta_2 \phi_2(s) + \ldots + \theta_F \phi_F(s), \label{eq:energylinearvfa} \end{align}\]

其中$\phi_f(s), f=1, \ldots, F$是一组适当选择的特征。特征的例子可能是

\[\begin{align*} \phi_1(s) &= R_t, \\ \phi_2(s) &= R^2_t, \\ \phi_3(s) &= p_t, \\ \phi_4(s) &= p^2_t, \\ \phi_5(s) &= p_{t-1}, \\ \phi_6(s) &= p_{t-2}, \\ \phi_7(s) &= R_t p_t. \end{align*}\]

请注意,包含常数项$\theta_0$的情况下,该模型只有八个需要估计的系数。抽样几百个状态应该足以获得一个良好的统计近似。这种方法对状态变量的数量相对不敏感,当然如果任何变量是连续的也没有问题。

每当我们使用如上述线性模型这样的参数模型时,一个挑战是我们必须指定特征$\phi_f(S_t)$。随着神经网络的流行,研究人员开始使用这种方法,包括可能需要估计数百万个参数的深度神经网络。这种方法的优势在于它消除了指定模型结构的需要,但代价是你需要更多的观测数据。深度神经网络提供了能够逼近任何函数的诱人特性,但这也意味着它们可能会拟合噪声。神经网络在复制已知的问题结构(如单调性——库存越大,价值越大——或凸性)方面也存在困难。

我们发现后向近似动态规划在一小部分问题上表现异常出色(有关后向近似动态规划与基准方法比较的总结,请参见Reinforcement Learning and Stochastic Optimization第15.4节),但并没有保证,其表现明显取决于选择一组有效的特征。在一个应用中,我们将标准后向马尔可夫决策过程算法的30天运行时间,缩短到了20分钟,并且所得到的解与最优解(由一个月的运行时间产生)相差在5%以内。但再次强调,这种表现没有任何保证。

前向近似动态规划

前向近似动态规划以一种直观的方式运作。设想我们从围绕决策后状态$S^x_t$的值函数近似$\Vbar^{x,n-1}_t(S^x_t)$开始,这是我们从算法的前$n-1$次迭代中计算出来的。我们首次在第1章中引入了决策后状态的概念,这是我们做出决策之后、但在任何新信息到达之前的状态。

现在,设想我们在算法的第$n$次迭代中处于某个特定状态$S^n_t$,沿着一条指导我们随时间向前推进采样的样本路径$\omega^n$。假设我们有一个函数$S^{x,n}_t = S^{M,x}(S^n_t,x)$,能将我们带到决策后状态。对于我们的$S^n_t = (R^n_t,p^n_t)$的能源问题,决策后状态将是

\[S^{x,n} = (R^n_t+x^n_t, p^n_t).\]

然后我们使用以下方式做出决策

\[x^n_t = \argmax_x \big(C(S^n_t,x) + \Vbar^{x,n-1}_t(S^{x,n}) \big).\]

给定$S^n_t$和我们的决策$x^n_t$,我们接着抽样$W_{t+1}(\omega^n)$,这转化为价格的变化$\phat^n_{t+1}$。然后我们模拟前进到下一个状态

\[S^n_{t+1} = (R^n_t+x^n_t, p^n_t + \phat^n_{t+1}(\omega)).\]

因此,我们只是在时间上向前模拟,这意味着我们不必关心状态变量有多复杂。有不同的策略可用于更新值函数近似$\Vbar^{n-1}_t$:

前向近似动态规划

第0步. 初始化:初始化$V^{\pi,0}_t,~t\in\Tcal$。设置$n = 1$。初始化$S^1_0$。

第1步. 对$n = 1, 2, \ldots, N$执行:

第2步. 对$m = 1, 2, \ldots, M$执行:

第3步. 选择一条样本路径$\omega^m$。

第4步. 初始化$\vhat^m = 0$。

第5步. 对$t = 0, 1, \ldots, T$执行:

第5a步. 求解:

$$x^{n,m}_t = \argmax_{x_t\in\Xcal^{n,m}_t} \big(C_t(S^{n,m}_t,x_t) + V^{\pi,n-1}_t(S^{M,x}(S^{n,m}_t,x_t))\big)$$

第5b步. 计算:

$$S^{x,n,m}_t = S^{M,x}(S^{n,m}_t,x^{n,m}_t), \qquad S^{n,m}_{t+1} = S^M(S^{x,n,m}_t,x^{n,m},W_{t+1}(\omega^m)).$$

第6步. 对$t = T-1,\ldots, 0$执行:

第6a步. 累积路径成本(以$\vhat^m_{T} = 0$):

$$\vhat^m_t = C_t(S^{n,m}_t,x^m_t) + \vhat^m_{t+1}$$

第6b步. 更新从时间$t$开始的策略近似值:

$$\Vbar^{n,m}_{t-1} \leftarrow U^V(\Vbar^{n,m-1}_{t-1}, S^{x,n,m}_{t-1}, \vhat^m_t)$$

其中我们通常使用$\step_{m-1} = 1/m$。

第7步. 对所有$t = 0, 1, \ldots, T$更新策略值函数$V^{\pi,n}_t(S^x_t) = \Vbar^{n,M}_t(S^x_t)$。

第8步. 返回值函数$(V^{\pi,N}_t)_{t=1}^T$。

这将实际的更新留在了一个更新函数$U^V(\cdot)$中,因为这取决于我们如何近似值函数。

前向近似动态规划之所以吸引人,是因为它能够扩展到高维问题。我们在任何时候都不会遍历所有可能的状态或结果。事实上,如果我们适当地近似值函数,从而能够利用强大的算法,我们甚至可以处理高维决策$x$。然而,前向近似动态规划(与后向近似动态规划一样)在性能保证方面几乎没有什么保障。

一种混合策略搜索-值函数近似策略

无论我们选择哪种方式来近似值函数,我们的策略都由以下给出

\[X^{VFA}(S_t) = \argmax_x \big(C(S_t,x) + \Vbar^x_t(S^x_t)\big).\]

我们向前模拟策略,其中我们让$\omega^n$表示外源信息的一条样本路径(即价格变化的集合)。经常出现的情况是,我们将在历史数据上测试我们的策略,在这种情况下只有一条样本路径。然而,如果我们已经开发了一个不确定价格的数学模型,我们可以创建一条样本路径$\omega$,用来近似一个策略的价值(我们也可以创建多条样本路径并取平均值):

\[\Fbar^{VFA}(\omega\vert S_0) = \sum_{t=0}^T C\big(S_t(\omega),X^{VFA}(S_t(\omega))\big).\]

这意味着该方法非常适合近似物流中可能出现的甚至是高维的问题。

当我们在构建一个值函数近似(属于前瞻类策略)时,我们通常不再有调整策略的步骤。然而,这并不意味着我们不能尝试。假设我们的值函数由方程$\eqref{eq:energylinearvfa}$中的线性模型给出。我们现在可以使用以下方式写出我们的策略

\[X^{VFA}(S_t\vert \theta) = \argmax_x \left(C(S_t,x) + \sum_{f=1}^F \theta_f \phi_f(S_t)\right).\]

使用我们的后向或前向近似动态规划算法之一来获得$\theta$的初始估计是有意义的,但正如我们上面所指出的,无法保证所得到的策略是高质量的。然而,我们总是可以将其作为起点,使其变得更好,

\[\Fhat^{VFA}(\theta,\omega\vert S_0) = \sum_{t=0}^T C\big(S_t(\omega),X^{VFA}(S_t(\omega)\vert \theta)\big).\]

现在,我们只需要求解我们可能提出的策略搜索问题,如下所示

\[\begin{align} \max_\theta \Fbar^{VFA}(\theta,\omega\vert S_0). \label{eq:optthetavfa} \end{align}\]

在我们之前的策略搜索示例中,$\theta$是一个标量,这使得这个问题相对容易。现在,$\theta$是一个可能有几十个维度的向量。我们将在稍后回到这个问题。

关于近似动态规划的一些警示说明

我们利用这个问题背景,对基于近似处于某状态之价值这一思想的方法进行了相对深入的介绍。这已在诸如”近似动态规划”或”强化学习”等术语下被研究。这些方法从学术研究界吸引了相当多的关注,但在实践中,这些方法并不容易。序贯决策问题无处不在,但实践中成功的应用相对罕见。

查找表近似——即我们为每个离散(或离散化)状态估计价值——在状态变量维度超过三维时无法扩展。使用诸如我们的线性近似这样的近似策略通常不起作用,因为这些近似必须在全局范围内准确,因为我们可能访问任何状态。与此同时,局部近似(这是非参数模型的一种形式)可能会遇到困难,因为局部近似的灵活性会带来不稳定性。

使这一过程更加复杂的是,我们依赖我们的近似值函数来做决策,这就产生了一个恶性循环。我们最初的近似不是很好,因此导致糟糕的决策。然后这些糟糕的决策又被用来更新值函数近似,由此你可以看出这种每况愈下的螺旋式下降。

将值函数近似进行调优的想法,正如我们在方程$\eqref{eq:optthetavfa}$中所做的那样,很有前景,因为它直接优化了策略的性能。奇怪的是,这个想法并没有被广泛使用。我们只需指出,使用这些模拟来调优$\theta$并不容易。因此,我们提醒任何打算尝试这些方法的读者要谨慎。

我们学到了什么?

习题

复习题

  1. 鉴于电力价格的重尾特性,使用时间序列模型有什么问题?
  2. 简要描述我们如何将落在正常波动范围内(三个标准差以内)的价格与更极端的观测值区分开来。
  3. 使用历史数据拟合经验分布应该能给我们一个与历史相符的概率分布。价格的随机模型中还可能存在哪些其他误差?
  4. 用文字描述混合时间序列部分中数据变换所完成的工作。
  5. 经典的反向动态规划很快会因维数灾难而崩溃。用文字描述反向近似动态规划如何克服维数灾难。例如,如果我们将状态变量的维度数量加倍,描述这会如何使反向ADP变得更加复杂。

问题求解题

  1. 写出文中给出的储能问题基本模型的五个要素。假设策略为低买高卖策略,写出目标函数 $$ X^{low-high}(S_t\vert \theta) = \begin{cases} -1 & \text{if } p_t < \theta^{buy}, \\ 0 & \text{if } \theta^{buy} \leq p_t \leq \theta^{sell}, \\ +1 & \text{if } p_t > \theta^{sell}. \end{cases} $$ 用对策略参数进行搜索的形式写出目标函数。同时,用反映每个随机变量的嵌套形式写出目标函数中的期望。
  2. 在建模不确定性的一节中,我们引入了如下给出的价格时间序列模型 $$ p_{t+1} = \thetabar_{t0} p_t + \thetabar_{t1} p_{t-1} + \thetabar_{t2} p_{t-2} + \varepsilon_{t+1}. $$ 书中描述了系数向量$\thetabar_t = (\thetabar_{t0},\thetabar_{t1}, \thetabar_{t2})$的更新方程。请记住,状态$S_t$包含了从时间$t$起对系统进行建模所需的*全部*信息,给出用于处理该价格过程的更新状态变量和转移函数。

编程题

这些习题使用位于tinyurl.com/sdagithub的Python模块EnergyStorage_I

  1. 使用该Python模块,对参数向量$\theta = (\theta^{buy}, \theta^{sell})$进行网格搜索,将价格中的$\theta^{sell}$从1.0变化到100.0,增量为$1,同时将$\theta^{buy}$从1.0变化到$\theta^{sell}$,增量同样为$1。价格将采用8天期间的实际历史每小时价格。
  2. 使用上文描述的反向动态规划策略求解最优策略(该算法已在Python模块中实现)。假设价格过程按如下方式演化 $$ p_{t+1} = p_t + \varepsilon_{t+1}, $$ 其中$\varepsilon_{t+1}$遵循基于实际历史价格差异的经验分布。
    1. 分别以$1、$0.50和$0.25的增量对价格进行离散化并运行该算法。计算三种离散化程度下状态空间的大小,并绘制运行时间与状态空间大小之间的关系图。
    2. 使用$1离散化下的最优值函数,将其性能与你在(a)部分中找到的最佳买卖策略进行比较。
  3. 从[tinyurl.com/sdamodelingsupplements](https://tinyurl.com/sdamodelingsupplements/)下载电子表格"Chapter8_electricity_prices"。使用"electricity prices"选项卡中的数据完成以下问题:
    1. 构建一个经验累积分布$F_P(p) = Prob[P \leq p]$,其中$P$是数据集中一周期间内某特定小时随机选取的价格。
    2. 令$F^{-1}_P(u)$为逆累积分布,其中$u$介于0和1之间。求出对应于$u = 0, 0.1, 0.2, \ldots, 0.9, 1.0$的价格$p(u) = F^{-1}_P(u)$。给这些价格中的每一个赋予1/11的概率,求出累积分布,并将其与你在(a)部分创建的累积分布进行比较。它们看起来匹配吗?
  4. 使用电子表格"Chapter8_electricity_prices",拟合如下形式的均值回归模型 $$ \begin{align} p_{t+1} = p_t + \beta (\mubar_t - p_t) + \varepsilon_{t+1} \label{eq:priceexercise} \end{align} $$ 其中 $$ \mubar_t = (1-\alpha)\mubar_{t-1} + \alpha p_t. $$ 假设$\alpha = 0.15$。求出使下式最小化的$\beta$ $$ G(\beta) = \sum_{t=0}^T \big(p_{t+1} - (p_t + \beta (\mubar^t-p_t))\big)^2. $$
    1. 通过进行简单的一维搜索来拟合均值回归模型(例如,以0.1为增量尝试0到1之间的值)。
    2. 根据你的样本计算$\varepsilon$的标准差$\sigma$(我们假设均值为0)。请注意,我们假设存在单一的常数标准差,尽管我们允许均值$\mubar_t$随时间变化。
    3. 使用你在(a)中找到的$\beta$的值,通过从均值为0、标准差为$\sigma$的正态分布中抽样$\varepsilon_{t+1}$,利用方程$\eqref{eq:priceexercise}$生成10条样本路径。在图中绘制这些样本路径,并将样本路径的行为与历史价格进行比较。它们看起来相似吗?
  5. 现在你将拟合如下给出的跳跃扩散模型 $$ p_{t+1} = p_t + \beta(\mubar_t - p_t) + \varepsilon_{t+1} + J_{t+1} \varepsilon_{t+1}, $$ 其中$J_{t+1} = 1$以某个跳跃概率(我们在下文计算)取值,否则为0,$\varepsilon^J_{t+1}$是跳跃发生时的随机跳跃幅度。 按照以下步骤拟合跳跃扩散模型,并将结果与历史进行比较。
    1. 使用习题11中的$\beta$值,遍历数据并识别所有落在$[\mubar_t \pm 3 \sigma]$范围之外的数据点。
    2. 使用你在习题11中找到的相同$\beta$值,对数据再次运行均值回归,但这次只包括未在(a)部分中被排除的数据点。求出未被排除数据的$\sigma$的新均值和标准差。
    3. 对保留下来的数据点再重复一次(b),同样排除位于$\pm 3 \sigma$范围之外的数据点。
    4. 将完成(c)部分时被排除的点所占的比例计算为跳跃概率(到目前为止,你已经运行了两次排除数据点的过程)。同时,计算被排除点的均值和标准差。
    5. 现在,使用保留点和排除点的最终均值和方差,运行你的跳跃扩散模型的10次模拟,并使用跳跃扩散概率对跳跃进行抽样。将这些模拟结果与历史进行比较,并讨论所得到的价格路径是否比你在习题11中找到的结果更真实,并将价格路径与实际历史进行比较。
  6. 我们将通过使用变换后的价格重复习题11和12的部分内容,再次尝试很好地拟合价格。
    1. 使用习题10中的累积分布,利用恒等式$U_t = F_P(p_t)$将每个价格转换为均匀分布的随机变量。
    2. 接下来,使用$Z_t = \Phi^{-1}(U_t)$将你的均匀分布随机变量$U_t$转换为均值为0、方差为1的正态分布随机变量,其中$\Phi(z)$是标准正态分布的累积分布,$\Phi^{-1}(U_t)$是该分布的逆函数。这在Excel中由norm.s.inv(p)函数体现,该函数返回对应于概率$p$(由变量$U_t$给出)的$Z$值。
    3. 我们现在必须在习题11的均值回归方程中重新拟合$\beta$。这一次,我们不使用价格$p_t$,而是使用归一化后的量$Z_t$;其余部分都相同(因此你可以直接遵循为原始价格拟合$\beta$时的过程)。
    4. 现在使用(c)中的模型(带有$\beta$的新值)创建一条$Z_t$值的样本路径。然后,使用$U_t = \text{norm.s.dist}(Z_t,1)$得到$Z_t \leq z$(在0和1之间均匀分布)的概率。最后,使用你在习题11中找到的累积分布,将$U_t$值映射回价格。绘制一条样本路径(如果你精通Excel,这并不太难——否则最痛苦的部分就是最后这一步)。
    5. 将所得样本路径的行为与历史分布进行比较。请注意,价格的分布应该是完美的,但价格序列可能仍然看起来拟合得不好。我们可能仍在犯哪些错误?