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

库存过剩不会得到任何补偿，缺货会受到严厉的惩罚（市价的200%）。因此，成熟度的值函数在I中具有非光滑的曲棍球棒形状，其斜率为零<1000，斜率为-2总的来说。提取和注射速率取决于库存，由CWDR（I）=-k√I和cinj（I）=krI+k-k、 这些功能规范源自储气设施的物理流体力学，见【34】。库存的结果动态为：Itk+1=Itk+cwdr（Itk）t、 如果mtk=-1（撤回）；如果mtk=0，则为0；（cinj（Itk）-k）t、 如果mtk=+1（注入）。（25）常数K测量注入成本，其表示为“气体损失”。这导致生产和注入之间存在一个适当的差距，因此不采取任何行动（m*= 0）如果价格“接近”平均水平，将是最佳措施，Pt≈ P在下面的数值实验中，我们将T离散为K=1000步，以便t=0.001T，其余参数列于表1中。切换成本取零K（i，j）≡0i、 j.无切换成本将（17）中连续函数q（tk，·）中的状态变量减少到（Ptk，Itk+1）。因此，在本例中，每个时间步都需要2D value函数的一个投影。α=2.38，σ=0.59，P=6k=2040.41，k=7.3·10，k=500，k=2500，k=620.5Imax=2000 MMcf，T=3，t=0.003，r=10%。表1：第6.6.2节“基准设置”中储气设施的参数，用于基准性能，我们比较了使用传统产品设计P×G（库存均匀离散化）的两种方案：3次全局多项式近似over（P，I）（PR-2D）和分段连续近似（在每个库存水平Ik，PR-1D上具有3次多项式回归）。

这些可以被视为具有联合回归的“经典”TvR方案【10，14】，以及【3，4，7，30】所使用的离散化-I版本。此外，在不同的设计上采用五种回归方法（PR-1D、GP-1D、PR-2D、黄土-2D和GP-2D），模拟预算为：N’10K、40K、100K：o通过LHS实施的随机空间填充设计（分段连续回归和Lotherwise为L×G）。我们使用一个大的保守输入域P∈ [2, 10].o 混合料-2D设计，40%的场地来自空间填充，其余60%来自接缝温度分布P，DM：=P（0.6N）∪L（0.4N）。为了实现P，我们需要估计pp（tk，·）和pI（tk，·）。这是通过首先以较小的预算和传统的产品设计运行算法来实现的。然后，我们生成前向路径（Pntk，^Intk）来估计tk处的关节（pP（tk，·），^pI（tk，·））。随着合同接近到期，边际分布^pI（tk，·）开始集中在I=1000附近，因此该领域其他部分的分辨率降低，并使用Lis来补偿这种影响自适应1D：对于1D回归，我们如上所述估计pP（tk，·），然后非均匀分解库存，以纳入最优控制库存过程集中在Imin、Imax和I=1000周围的事实。附录A中详细说明了每个模拟预算的离散化水平。动态时间相关设计改变了模拟预算Ntkacross tk。我们使用规范Ntk（N（1），N（2））：=N（1）{k<900}+N（2）{k≥900}这样（大约）0.9N（1）+0.1N（2）∈ {10K、40K、100K}。

附录中的表5给出了具体规格。我们使用了它的两种变体：–固定投影空间H（即GP-1D）和常规产品设计D=P×G。–时间相关投影和设计，即GP-2D和混合设计用于k<900，PR-1D和常规设计用于k≥ 上述动态方案的动机是通过投入更大的模拟预算以及使用-1D回归，更好地处理非光滑终端条件。对于GP，我们使用Matlab的内置实现fitrgp。对于黄土，我们再次使用Matlab中的curvefittingtoolbox，该工具箱基于三立方权重函数κ（x*, xn）=1.- （| x）*-xn |λ（d，x*)). λ（d，x）以上*) 是x的欧几里德距离*到跨度d中最远的xKw。我们使用默认跨度d=25%，将P保持在其原始比例上，并将比例I重新调整为范围[0，2]。为了减少HGPR和黄土的回归开销，我们在基础设计类型的基础上使用了丁腈橡胶（见附录A）的分批设计。为了比较不同设计/回归的性能，我们使用P=6，I=1000时的值函数^V（0，P，I）的估计值，使用N=10000的固定样本路径集（即固定Pn0:T）作为性能度量。6.3结果表2显示了不同设计和回归方法的性能。我们继续讨论三个不同方面的结果：（i）不同回归方案的影响，特别是参数PR与非参数GP和黄土方法；（ii）模拟设计的影响；（iii）接头-2D vs。

插值-1D方法。首先，我们的结果证实，插值1D方法在这个经典示例中表现非常好，完美地利用了具有一维库存变量的二维设置。从这个意义上讲，现有的最先进技术已经非常优秀了。这有两个重要原因。第一个PR-1D高度灵活，有很多自由度，允许良好的安装（由于1D设置，过度安装的危险降至最低）。其次，PR-1D完美地利用了价值函数在库存中几乎（分段）线性的事实。我们通过用GP-1D替换PR-1D获得了轻微的改进；另一个微小的收获是用自适应方法代替等间距库存离散化，使更多的级别接近库存边界。作为进一步的增强，动态设计利用了一种依赖于步长的模拟预算（由于非光滑的“曲棍球棒”终端条件，需要更多的努力才能从统计上学习，从而获得边界层效应），从而带来了显著的改善，突出了跨时间步长混合和匹配近似策略的潜在好处。通过采用依赖时间的方法，可以有效地节省仿真预算（例如，动态GP-1D的估值N=10与自适应GP-1D的估值N=2·10相当）。然而，正如我们反复强调的那样，-1D方法的范围必然有限，尤其是当添加更多的因素/库存变量时，不能进行缩放。1D方法的另一个限制是它们对产品设计的要求。更通用的二元2D回归对所有状态变量进行了统计上的公平处理，因此允许任意模拟设计。由此产生的潜在实现的巨大范围既是福也是祸。

因此，我们记录了在选择回归方案和选择模拟设计方面的一些好的和一些坏的选择。我们发现，PR-2D往往表现不佳，这并不奇怪，因为它对连续值实施了严格的参数化形状，灵活性空间不足。同样，我们观察到黄土的中等性能；另一方面，探地雷达通常工作得很好。将我们的注意力转向不同的2D模拟设计，我们将传统的P×G选择与3种备选方案进行比较：在大输入域上填充保守的空间；联合概率设计P；混合前两者的混合物设计。我们发现，普通空间和联合概率都不能很好地工作；第一个目标不够明确，在对^V贡献不大的地区花费了太多预算；第二种方法过于激进，通常在计算^m时，requiresextrapolation会产生不准确的预测。相比之下，混合料设计（我们发现60/40的混合料效果良好）是一个赢家，显著优于传统的混合料设计。特别是，混合设计的GP-2D是唯一一个与GP-1D并驾齐驱的二元回归方案。这一点很重要，因为与GP-1D不同，GP-2D可以直接扩展到更高的维度，并且不需要产品设计（或通常速度较慢的任何插值）。这些发现突出了设计选择中提案密度的重要性。为了突出我们算法的灵活性，我们还提出了另一种结合自适应PR-1D和混合GP-2D的动态设计。

这种组合保持了相同的精确度，但由于PR-1D的回归开销较低，运行速度快了约10-15%。回归模拟预算设计方案中低规模常规PR-1D 4965 5097 5231GP-1D 4968 5107 5247PR-2D 4869 4888 4891黄土-2D 4910 4969 5011GP-2D 4652 5161 5243空间-fillingPR-1D 4768 4889 5028GP-1D 4854 5064 5224PR-2D 4762 4789 4792黄土-2D 4747 4912 4934GP-2D 4976 5080 5133自适应1DPR-1D 5061 5187 5246GP-1D 5079 5195 5245动态CGP-1D 51325225 5266混合5137 5205 5228混合物2DPR-2D 4820 4835 4834黄土-2D 4960 4987 5003GP-2D 5137 5210 5233表2：使用不同的设计回归对和三个模拟预算进行估值^V（0，61000）（千）：低N’10K，中N’40K，大N’100K，参见附录A。估价是每个方案10次运行的平均值，但预算较大的黄土-2D除外：由于黄土的管理费用过高，仅进行了一次运行。P设计对-1D方法的影响。在注重P回归和I离散+插值的-1D方法中，将设计从P×G改为空间填充L×G对所有方法都有不利影响，例如，PR-1D观察到估值下降了200K以上（参见右图4）。此外，使用空间填充但“粗略”模仿Pt分布的中间设计（例如，范围内60%的站点∈ [4.5，7.5]，P=4.5以下20%，P=7.5以上20%）在L×G和P×G之间插入估值，确认P维度设计的重大影响。I维模拟设计的形状在评估中具有同样重要的作用。

通过比较传统PR-1D和自适应PR-1D（表2）的性能可以观察到这一点，它们基于PTP的边际分布共享相同的P-P设计，但使用不同的方法离散库存。自适应设计中的非均匀离散化提高了接近Imin、Imax和I=1000的精度，并导致相对于传统设计的更高估值，还有一个额外的优势，即我们不需要在内存中存储数十万条路径。使用GPR可视化设计影响。除了更高的估值外，GPR还提供了一种通过后验标准差量化设计效果的自然方法。在图5中，我们绘制了基于混合和空间填充2D设计的后验GP标准偏差s（x）（局部估计标准误差的代理）。这显示了设计形状和探地雷达学习精度之间的直接联系。对于空间填充设计，我们观察到一个恒定的后验方差，即GPR在回归域内部学习^q（tk，·）同样好。相比之下，[Pmin，Pmax]×[Imin，Imax]边缘的后验标准差非常高，这对于在库存几乎为满或几乎为空时正确识别策略是个问题。对于混合料设计，与图3c中的联合概率设计销相比，后验方差反映了对角线上输入点的浓度。反过来，这有利于学习控制图，因为GPR预测沿切换边界最准确，见图2a。切换边界周围的更高精度允许混合料设计将模拟预算分配到最需要精度的区域，并且部分解释了与LHS设计相比更高的估值。

我们强调，这种局部推理质量仅适用于非参数工具；像PR这样的全球计划不能直接从专注的设计中受益。2价格02004000600800100012001400160018002000库存6.56.56.588889.59.59.511111111313131515151515171717172 4 6 8价格02004000600800100012001600108002000库存9.59.51111111131313151515151717191919图5：设计对t=1.5年时探地雷达后验标准误差s（x）的影响。我们展示了x 7的轮廓→ s（x），带点，表示N=500个站点的基础各自设计。左面板：图3a的空间填充设计（Sobol QMC S）导致S（·）的矩形水平集。右面板：在混合设计下，s（·）类似于图3c中的销的形状。使用GPR复制设计。GPR的实施还需要管理设计站点数量NSA和复制量Nb之间的权衡。使用复制是必要的，因为拥有超过2000个不同的站点非常耗时，至少对于我们使用的GPR的有效实施来说是如此。在图6a的左窗格中，我们考虑了不同的Ns和Nb（因此N）。虽然较大的模拟预算明显改善了结果，但我们注意到，最终通过Ns固定增加NB并不能改善回归质量（尽管它仍能减少标准误差）。在图6b的右侧面板中，我们展示了Nsfor已执行总预算的影响N=Nb·Ns。我们发现，当NSI越大（即NBI越小）时，复制总体上是次优的，结果越好。

为了管理最终的速度/精度权衡，我们建议采用Nb∈ [20, 50]; 例如，当N=10时，我们使用Ns=2000，Nb=50，当N=10时，我们使用Ns=500，Nb=20.10 20 40 100 160复制数量490049505005105150520020550评估设计站点=1000（a）复制的影响Nb250 500 1000 2000个设计站点485049505000505051051505205520550评估预算=100000（b）不同设计站点的影响NSFix图6：左面板：混合GP-2D的性能作为N的函数。我们确定唯一设计站点的数量Ns=1000，并逐步增加增加重复次数Nb，在P=6，I=1000时报告结果^V（0，P，I）。右图：增加固定N=10的NSF的影响。结果是每个算法运行20次。在每个方框图中，中心标记表示中值，方框表示第25和75个百分位，胡须表示最极端的跑步。值得注意的是：我们的实验表明了以下关键观察结果：（i）在inventorydiscretized-1D方法中，高斯过程回归在所有情况下都优于标准多项式回归，并且对“差”设计或低模拟预算更为稳健。（ii）在联合2D方案中，我们继续观察到GPR的性能优于交替双变量回归方案（黄土和PR）。此外，GP-2D与性能最好的GP-1D不相上下。我们强调，GPR和黄土的高效实施依赖于批量设计，这是我们DEA实施中的另一项创新。（iii）我们还发现，设计选择与性能之间存在着强烈的依赖性。我们确认，最佳结果来自于设计，例如我们实施的混合和动态版本，这些设计平衡了输入空间的填充和针对大多数（Pt，^It）轨迹所在的领域。否则，纯空间填充或逆向渐进边界会降低性能。

（iv）在两种参数化方法PR1D和PR-2D之间，我们发现无论设计和模拟预算如何，PR-1D都显著优于PR-2D，表明分段连续回归的优势。6.4具有转换成本的储气库模型化我们通过合并转换成本来概括前面的示例。切换成本使控制图取决于当前的管理模式，并产生惯性，即倾向于继续使用相同的管理模式，以降低总体成本。为了处理离散的m维，我们将其视为| J |不同的连续函数，通过| J |不同的回归进行估计。否则，算法将以与之前完全相同的方式进行。这说明了DEA处理一系列问题公式的灵活性。除了通过k的注入损耗外，我们还添加了转换成本k（i，j），具体如下：k(-1，1）=K（0，1）=15000；K（1，-1） =K（0，-1) = 5000; K（1，0）=K(-1，0）=0，（26），即切换成本仅取决于控制器决定切换到的状态，切换到注入成本最高，切换到无动作。在图7中，我们展示了控制器针对不同制度的政策。由于相应的控制图具有最宽的存储区域，因此处于mtk=0状态的惯性很明显。当比较左侧和中间面板的存储区域时，K（i，j）的影响也很明显。如果控制器从注入模式切换到存储模式，由于切换成本，她会在尝试返回注入模式时发现更多阻力。（a） m=+1（注入）（b）m=0（存储）（c）m=-1（撤销）图7:t=2.7年时模型的控制图^m（t，P，I，m），转换成本为（26），对于m∈ {+1, 0, -1}. 颜色为^mt+t=+1（注入，浅黄色），^mt+t=0（储存，中青色），^mt+t=-1（撤回，深蓝色）。