随机存储问题的模拟方法：一种统计方法

我们参考文献[37]中的一些自适应构造，包括等权重和等网格。4.1.1分段连续近似另一种方法是通过线性插值构造一维连续的分段近似。例如，将内生库存变量离散化为MI+1水平I，I，IMI，我们在P中为每个级别设置了一个独立的阶数MPmonomial，即优化^hj（P）：=Piβijφi（P），j=0，MI，给出总（MI+1）mP回归系数。（相反，我们也可以建立一个分段线性模型，在P维中包含多个子域。）任意I的最终插值预测∈ （Ij，Ij+1）则逐段定义为bhtk（P，I）：=δ（I）^hj（P）+（1- δ（I））^hj+1（P），（20），其中δ（I）=I-IjIj+1-Ij。4.2局部多项式回归（黄土）局部回归通过构造一个非参数函数，解决每个预测点的优化问题，最大限度地减少了对基函数选择的担忧。给定一个数据集{x，y}（作为提醒，x≡ （P，I）在本节中是二维的），使用x处的黄土进行预测*是h（x*) =Pri=1βi（x*)φi（x*) 其中局部系数βi（x*) 由加权最小二乘回归~β（x）确定*) = arg最小值~β∈RrNNXn=1κ（x*, xn）hyn-~βT~φ（xn）i.（21）权函数κ（x*, x）∈ [0，1]通过接近x的输入为yn提供更多权重*, 类似于tokernel回归。因为每个x都需要单独的优化*, 作出预测（21）的复杂性为O（NN）。[26]在采矿实物期权的背景下实施，Loesss通过一种“滑动技巧”得到增强，该技巧将复杂性降低到O（（N+N）log N）。4.3高斯过程回归（GPR）GPR与黄土非常相似，因为预测h（x*) 是采样输出y的加权平均值。

基本关系h:x 7→ y被视为阿高斯随机场的一种实现，即{h（xi）}Ni=1是来自多元高斯分布的样本，平均值{m（xi）}Ni=1，协方差函数{κ（xi，xj）}Ni，j=1。在文献中的许多选项中，可以说最流行的是平方指数核，κ（xi，xj）=σfexp-（xi）- xj）T∑-1（xi）- xj）,其中，∑是对角矩阵。传统上，矩阵的对角元素∑被称为长度尺度参数，σfas被称为信号方差。它们通常一起被称为超高参数。长度比例参数决定了相应尺寸表面的平滑度，而σf则决定了波动的幅度。在图2c中，我们观察到，对于固定价格，连续值函数在库存维度上是线性的，在价格维度上具有非线性行为。GPR通过其尺度参数捕捉到这种差异，对于库存而言，尺度参数“大”（相关性衰减缓慢，曲率很小），对于价格维度而言，尺度参数“小”（快速相关性衰减允许P的波动）。通过似然最大化估计超参数。对于先前的平均值，我们取m（x）=β，其中β与其他超参数一起学习。对于任何站点x*, h（x*) 是一个随机变量，其条件分布{x，y}为：h（x*)|y~ Nm（x*) + H*H-1（y- m（x）），H**（十）*) -H*（十）*)H-1小时*（十）*)T（22）其中N×N矩阵协方差矩阵H和N×1向量H*（十）*) areH公司：=κ（x，x）κ（x，x）。κ（x，xN）κ（x，x）κ（x，x）。κ（x，xN）。。。。。。。。。。。。κ（xN，x）κ（xN，x）。κ（xN，xN）, H*（十）*)电话：=κ（x*, x） κ（x*, x） 。。。κ（x*, xN）, H**（十）*) = κ（x*, x个*),（23）式中，κ（xi，xj）=κ（xi，xj）+σ。

因此，在x处的预测*is^h（x*) = m（x*)+H*（十）*)H-1（y-m（x））和后验GP方差s（x*) := H**（十）*)-H*（十）*)H-1小时*（十）*)t提供此预测的不确定性度量（类似于标准误差）。GPR一般具有O（N+NN）复杂性，类似于内核回归。GPR通常表现得非常好，但对于N而言，其成本过高 1000、备注5。有一个广泛的GP生态系统，允许对所描述的GPR进行多次扩展；特别相关的方面包括更高级的m（·）先验平均值规范；其他核家族κ（·，·）；能够处理状态相关仿真噪声σ（·）的异方差模型；选择GP超参数的进一步技术；的分段模型允许空间非平稳协方差。由于我们的重点是说明性的演示，而不是寻找性能最佳的方法（这必然取决于问题），因此我们将注意力限制在基本的GP方法上，并提供现成的实现。5模拟设计算法1的第二个核心部分涉及设计Dk。与回归子问题一样，该模板为用户在给定模拟预算的情况下选择DK提供了很大的自由度，并提供了许多可行的方法。我们确定了潜在设计的四个相关方面：jointvs。产品自适应与空间填充；确定性与随机性；独特而不重复。最后但并非最不重要的是，我们讨论了设计尺寸N。为了设置舞台，让我们总结一下“常规”设计【4、7、30】。传统上，解决存储问题的设计依赖于全局路径与库存离散化的混合。在价格维度中，它包括为价格过程选择NP初始条件Ptattime和抽样n=1，NPpaths Pntk+1遵循条件密度p（tk+1，.| tk，Ptk），直到终止日期T。

结果集合（Pntk）用于每个tk处的设计“网格”。对于内生库存维度，I被离散为NIlevels：Il=l我I=Imax-伊曼尼-1，l=0，1，镍-1、总体设计DKN=NPNIsites，构造为Cartesianproduct{Ptk，Ptk，…，PNPtk}×{I，I，…，INI-1}. 在我们的术语中，这是一种产品设计，在P中是自适应和随机的，在I中是空间填充和确定性的，并且没有复制。其形状便于第4.1.1节的分段连续回归方案，该方案分别处理P和I坐标。同时，这种D可以被视为联合（P，I）域中的一个站点云，并与二元回归方案一起使用。为了总结模拟设计的范围，我们使用了一个简单的术语。产品设计标识为DP×DI，而接头设计用单个符号表示。必要时可使用下标来确定各自设计的尺寸。不同的设计类型由不同的字母标识。例如，我们将G用于网格设计，将P用于基于密度的“概率”设计。因此，上述常规设计得到了labelledas P×G.5.1空间填充设计，以实现学习（P，I）7的目标→ q（tk，P，I，m）我们需要探索整个输入域的连续值。实现这一点的一个简单机制是展开设计场地以填充空间。网格设计，使用网格大小统一选择设计场地, 上文所述的传统方法就是此类填空序列的一个例子。空间填充可以是确定性的，也可以是随机的。对于确定性情况，除了上述网格G之外，还可以使用各种准蒙特卡罗（QMC）序列，例如，Sobol序列S。

Sobol序列在维度d>1时很有用，在维度d>1时，它们可以产生任意大小N的d变量空间填充设计，而网格仅限于大小NP×NI的矩形构造。从理论上讲，QMC序列也可以保证对特定矩形域提供良好的“统一”覆盖。我们试验了以下两种设置：o由于我们只有一个外生变量P，我们在每个时间步长tk，k=0，…，采用了一个一维Sobol序列，大小N限制为[Pmin，Pmax]，K、 然后，我们将inventorydimension离散为{I，I，…，INI}，类似于传统设计，最终D是产品×G。或者，我们从受限域[Pmin，Pmax]×[Imin，Imax]上的二维Sobol序列S={Pn，In}Nn=1生成N个设计位点。随机空间填充设计的一个例子是Pntk~ Unif（Pmin，Pmax）i.i.d.由于均匀样本倾向于聚类，因此有方差减少的版本，如拉丁超立方体样本（LHS）。在二维中，LHS设计将输入空间划分为矩形阵列，并确保每一行和每一列正好有一个设计站点。请注意，如果跨时间步使用相同的确定性空间填充方法，则设计Dtk≡ D在tk中变得相同。这可能会从回归模式中产生“混叠”效应，即由于重复回归和各自的误差反向传播，在（Pn，In）周围产生近似伪影。更改或随机化Dtkacross-tk是一种补救方法，通常作为实现默认值首选。

空间填充设计的另一个问题是需要指定边界框【Pmin，Pmax】×【Imin，Imax】——这在I中很容易，但在P维度中并不明显，因为PTI是无界的。5.2自适应设计与旨在覆盖输入空间（并以全局L（Leb）近似误差为统计目标）的空间填充设计相比，自适应设计试图利用手头问题的特定结构。他们的出发点是观察到^vd的质量取决于沿受控路径（Pntk，^Intk）正确预测存储行为。关于系统状态分布的信息已经在传统方法中得到利用，该方法使用分布式概率设计进行DP。类似地，[30]在I中使用了非均匀离散化，以获得更多更接近Iminand Imax的网格点。上述讨论激发了更一般的适应性设计，旨在将DTA的形状匹配到建议密度p（·）（p和I为双变量或单变量）。例如，随机联合概率设计样本来自（P，^I）的二元分布。另一种方法是使用各自的边际密度构建概率产品设计P×P，采用PNTK~ pP（tk，·）和库存Intk~ pI（tk，·）。还可以通过量化建议密度p（类似于数值求积方法）来构建确定性自适应设计，以返回带有NPsites的离散表示。概率设计确保更好地估计（Pn0:T，^In0:T）轨迹最可能出现的区域中的连续函数。备注6。构建内生变量^itk的分布pI（tk，·）需要一些远见。

一种策略是使用模拟预算的一小部分，使用传统/空间填充设计创建策略，并在此策略上运行正向模拟，将正向模拟的结果保存为创建最优控制图设计的联合分布代理。5.3分批设计GPR和黄土等非参数技术通常会提高精度，但会增加回归开销。事实上，对于这两种方法，复杂度至少是站点数量N的二次方，这可能会对N造成阻碍 克服这一障碍的一个解决方案是使用复制，即重复使用同一设计站点进行多个模拟。因此，虽然有数千个不同的设计站点相当于模拟预算N，但我们只选择了几百个不同的站点NSA，并从每个设计站点生成Nb：=N/NSPATH。形式上，这意味着我们区分了N个初始条件（Pntk，Intk）Nn=1，用于模拟路径延拓值和唯一的设计站点Ns N构成设计D。后者可以是任何类型、空间填充、自适应等。给定设计场地（Pntk，Intk+1）Nsn=1，我们绘制Nbdraws Pn，（m）tk+1，并评估相应的路径延拓值v（m）k+1（Pn，（m）tk+1，Intk+1），m=1，注意。对于基于核的技术，如GPA和黄土，则可以使用预平均值，即首先评估经验平均值：(R)vnk+1（Pntk，Intk+1）=NbNbXm=1v（m）k+1，跨越Nbreplicates。然后将得到的回归数据集（Pntk、Intk+1、？vnk+1）Nsn=1输入回归方程，以估计连续函数。除了减少回归所花费的总时间（可以很容易地达到几个数量级），分批设计还具有在每个设计地点减少模拟方差'v的优势，从而提高信噪比。

虽然复制设计是次优的（在最大化统计近似质量方面），但在实践中，对于大N（通常为几十万），精确性损失很小，并且超过了保证的计算收益。除了在每个站点使用固定数量的复制品的统一批处理方法外，还可以在切换边界周围使用更多复制品的自适应批处理（即Nbas站点特定），以获得更好的精度[6]。5.4动态设计DEA让我们能够在不同的时间步轻松地组合不同的设计。例如，可以改变分步模拟预算，使用更大的接近成熟度的N来有效地捕捉连续函数中终端条件的影响，然后使用更小的预算。一些回归方法可能在计算上很昂贵，但更“精确”（GPR），而其他方法可能不太“精确”，但计算效率很高（多项式回归）。明智地将两者结合起来，可以节省时间，而不会造成太多的损失。我们用两个插图来结束本节。在图3中，我们展示了四种具有代表性的设计：2D S中的Sobol空间填充序列、2D L中的LHS、联合概率P和常规设计P×G。这些也将在下面的数值示例中使用。虽然Sobol序列以对称模式填充输入空间，但LHS是随机的。联合概率设计模拟了状态变量（Pt，^It）的分布，将大多数站点放在库存Imin、Imax和价格平均值附近的边界。请注意，这种设计非常“激进”，并且没有探索足够的输入空间，这会使回归估计变得不稳健。因此，在以下数值示例中，我们将概率设计（通过统计混合）与其他类型混合。

最后，传统设计将库存离散化，同时保持价格维度的自适应分布。2 4 6 8 10 12价格02004000600800100012001400108002000库存（a）2 d Sobol QMC序列2 4 6 8 10价格02004000600800100012001600108002000库存（b）2 d LHS设计L2 6 8 10价格02004000600800100012001600108002000库存（c）联合概率设计P2 4 6 10 12价格02004000600800100012001600108002000库存（d）常规产品设计P×G图3：插图不同的模拟设计D。在所有情况下，我们取N=500个设计点。接下来，在图4中，我们显示了回归和设计对连续值函数和相应控制图的影响。左侧面板比较了成熟度t=t之前最后一步的GP-2D和PR-2D的连续值函数-t、 我们观察到，PR-2D与GP-2D相比表现不佳，后者几乎完全符合“曲棍球棒”的处罚。此外，由于I中的分段回归，1D回归具有“非光滑”切换边界。这在图的中间面板中很明显，其中控制图显示I处的跳跃∈{500，1000，1500}，库存的中间离散化水平（我们使用NI=5I=500）。当NI增加时，这种行为变得不那么突出。图的右面板显示了PR-1D的两个控制图，分别为LHS和常规设计。

保守域（我们使用[Pmin，Pmax]=[2，20]）对LHS的影响是显而易见的，因为我们注意到LHS变体中的存储区域太宽，当价格为∈ [7，7.5]导致估值较低。500 1000 1500 2000库存-6000-5000-4000-3000-2000-100001000连续值数据GP-2DPR-2D0200400600800100012001400108002000库存2价格BWDR2DBINK1D0200600800100012001400108002000库存2价格BWDR1DBINK1D02004006008001000120016018002000库存2价格BWDR常规BINJConventionalBWDRLHSBINJLHSF图4：左面板：不同回归的连续函数q（t，P，·）比较步骤到期日t=t- t和P=6。中央面板：t=2.7年时的控制图^m（t，P，I），对应于GP-1D和GP-2D回归。对于GP-1D，我们使用Ns=100、Nb=20和NI=5。对于GP-2D，我们使用Ns=500和Nb=20。水平线表示GP-1D设计的投资离散化水平。右面板：与常规和LHS设计相对应的控制图，带有PR-1D回归。我们使用NP=2000，NI=21和pricedomain来表示LHS[Pmin，Pmax][2，20]。6天然气储存设施6.1市场和储存说明在本节中，我们重新讨论了Chen和Forsyth[11]提出的天然气储存问题。异源状态过程P遵循对数平均值还原动力学Tk+1- Ptk=α（P- Ptk）t+σPtkWtk，（24），其中PTI为每单位天然气的现货价格，报价单位为“百万美元/英国热量单位（$/百万英热单位）”。存货报价单位为百万立方英尺（MMCf）。由于粗略的1000mmbtu=1MMcf，在计算收入或利润时，我们将PT乘以10。到期时的惩罚函数W（PT，IT）为W（PT，IT）=-2PTmax（1000- IT，0）。因此，目标库存是IT=1000或50%容量。