基于数据的概率约束混合整数非线性双层优化问题及Copulas方法的应用

《Industrial & Engineering Chemistry Research》：Data-Driven Chance Constrained Mixed Integer Nonlinear Bilevel Optimization via Copulas

　　高分辨率图像
下载MS PowerPoint幻灯片

过程供应链的核心功能是计划和调度，以确保其盈利能力和响应性运营。同时，现实世界中的需求数据挑战，如多变量数据依赖性和相关性、噪声数据分布以及参数不确定性，显著增加了决策过程的复杂性，使得识别最优解决方案变得困难。这需要采用数据驱动的优化方法，有效考虑观察数据中固有的依赖结构，从而高效探索解决方案空间并识别高质量的结果。本研究提出了一种数据驱动的概率框架，该框架结合了机会约束编程（CCP）和copulas。CCP是一种基于优化的方法，在指定的风险阈值下强制执行概率约束。Copulas是一种统计技术，用于在捕捉多个不确定参数之间的依赖性的同时估计数据。所提出的框架能够准确建模不同场景下的变量依赖性，尤其是当数据表现出复杂分布或非平凡依赖性时。我们重点关注了三个来自过程供应链计划和调度以及原油调度的综合案例研究。我们将提出的估计框架与数据驱动优化双层混合整数非线性问题（DOMINO）框架相结合，后者是一种用于解决双层优化问题的数据驱动灰箱算法，以得出具有保证需求满足率的决策。计算实验证明，我们提出的基于copula的机会约束优化框架能够考虑需求相关性，并实现更高的联合需求满足率、更低的总成本和更高的效率。

1. 引言
随着大数据的兴起，数据量的快速增长为实际工业应用中的高级数据驱动决策提供了有利机会。传统的确定性方法依赖于过于简化的模型，这些模型往往无法充分捕捉不确定性、非线性动态和强系统依赖性，从而限制了它们利用数据资源的有效性。因此，对于许多工业相关的优化问题来说，找到最优或可行的解决方案仍然具有内在挑战性。几个关键障碍导致了这种困难，包括缺乏关于数据分布的先验知识、随机需求波动的存在、内部数据相关性及结构，以及参数不确定性，所有这些都会阻碍传统模型有效解释复杂问题结构并识别高质量解决方案以支持精细的决策过程。相比之下，数据驱动优化（DDO）利用了现实世界数据的粒度，提供了提取复杂数据模式的明显优势。因此，DDO越来越多地被用来识别数据价值，利用统计资源来支持更智能的运营。这里感兴趣的一类问题是企业范围内的优化范式，最初由Grossmann引入。不确定性是过程操作中的一个关键实际因素，这已成为许多领域研究的热点，包括微电网市场 bidding、在线调度、随机工业气体市场建模和氢基础设施规划。特别是计划和调度的整合，在过去二十年中引起了研究界的兴趣，因为决策者可以通过利用上述功能的相互依赖性获得好处。为此，已经提出了不同的框架来捕捉这种依赖性，同时平衡由集成问题的多尺度性质所带来的计算 effort。

机会约束编程（CCP）是一种在存在不确定性的情况下进行决策的优化方法，它保证了在一定的可靠性水平下满足随机约束，从而提供了控制风险和考虑决策过程中固有变异性的灵活性。许多研究已成功将CCP集成到各种优化问题中，应用于能源生产、机场运营管理、水资源管理、供应链管理、电力系统优化、随机项目网络调度和金融风险管理。为了确保CCP的可靠性，需要一种统计方法来建模和评估多个不确定参数之间的依赖结构。Copulas是一类专门设计用于捕捉具有内在相关性的多个随机变量之间依赖性的数据驱动耦合函数。将copula公式集成到CCP框架中，增强了其在不同场景下更好地建模变量依赖性的能力，并利用历史数据模式，特别是在底层数据表现出复杂分布、非平凡依赖性和相关风险的情况下。这种集成提高了在具有不确定参数的优化问题中的决策过程的估计准确性。

近年来，copulas和CCP的结合显示出巨大潜力，利用两者的优势开发了基于风险的框架，能够对不同来源的不确定性及其依赖性进行定量分析。Ding等人开发了一个基于copula的CCP框架，利用vine copulas进行风险管理，重点关注与管理效率相关的多维相关不确定性，旨在解决水质管理中的非点源污染问题。Khezri和Khodayifar应用联合机会约束来解决多商品网络流问题中的容量需求，使用Archimedean copula来估计不同容量系数变量之间的依赖性。Hosseini Nodeh等人利用copulas来建模随机资源受限最短路径问题中的资源消耗依赖性。Alizadeh等人结合copulas和CCP来解决来自水库灌溉区系统的多个不确定来源问题，旨在平衡农经济参数的最优值、目标灌溉供应可靠性水平和决策者的风险态度之间的相互关系。Chen等人提出了一种联合机会约束方法，该方法从一组广泛应用的bivariate copulas中选择表现最好的copula来表征不确定性概率分布和多周期、多区域住宅垃圾数据的相关性质。

在这项工作中，我们提出了一个数据驱动的机会约束框架，旨在将不确定性转化为定量的风险水平。该框架结合了集成copulas，以在过程调度优化问题的背景下估计多变量依赖性来源。基于机会约束的方法比传统的基于场景的方法具有明显优势，因为它消除了显式跟踪每个场景的需要，从而减少了资源强度和计算时间要求。此外，我们将提出的优化框架与DOMINO（数据驱动的双层混合整数非线性问题优化）框架相结合，后者是一种用于解决具有各种下层公式的双层优化问题的数据驱动灰箱算法。DOMINO框架已被证明可以解决确定性（6,39,40）和随机（7）集成计划和调度问题，后者通过基于场景的优化来处理。这项工作的目标是减少DOMINO中解决重复场景的计算负担，同时通过基于copula的机会约束优化方法保证需求满足。

这项工作的新颖之处在于将机会约束和DOMINO结合起来，形成了一种高效的数据驱动优化方法，该方法通过copula技术保留了依赖结构。具体来说，将机会约束copula集成到DOMINO框架中，产生了一个更高效的风险基础决策过程，在该过程中，不确定性参数是根据预定义的风险水平和内部数据依赖性来估计的，而不是通过显式计算所有场景来估计的，从而大大减少了相对于现有基于随机编程框架的计算时间。我们在一个激励性的单阶段连续搅拌罐反应器案例研究和一个更高维度的多产品原油规划与调度案例研究中测试了我们的优化框架。我们的结果表明，所提出的基于copula的数据驱动机会约束模型在产生具有更高利润和更低 penalty 的稳健决策方面表现出色，与具有相似风险水平的传统方法相比。计算实验进一步证明，所提出的机会约束copula方法能够有效识别稳健的解决方案，实现更高的联合满足率和接近最优的性能，同时显著减少了相对于现有随机优化方法的计算时间。

本文的其余部分组织如下。第2节介绍了copula理论、机会约束公式以及两者结合形成数据驱动机会约束copula框架的过程。第3节通过一个激励性研究和两个案例研究展示了我们提出的数据驱动机会约束copula框架。第4节包括结论和未来工作。

2. 方法论
2.1. Copula理论
Copulas是全面分析多变量数据分布中所示依赖结构的重要统计工具。Copulas可以分为参数化、半参数化或非参数化类型，允许灵活指定边际分布。一个n维copula函数C将来自未知底层分布的多个边际累积分布连接或“耦合”起来，形成[0, 1]n上的多变量联合分布，同时保持分布间的依赖结构。Sklar定理（42）形式化了任何多变量分布函数F(x1, x2,···, xn)可以用两个不同组件来表示的数学基础：一组单变量边际分布函数Fi(xi)（i ∈ {1,···, n}），每个函数代表单个变量的行为，以及预定的n维copula函数C，它封装了边际分布之间的依赖结构，如方程1所示。
??(??1,??2,...,????)=??(??1(??1),??2(??2),...,????(????))
F(x1,x2,...,xn)=C(F1(x1),F2(x2),...,Fn(xn)
(1)
其中x1, x2, ..., xn是从边际分布F1, F2,···,Fn中获得的随机变量(x)的实现值。这使得可以使用选定的copula和指定的边际分布根据观察数据构建任何多变量分布，无论底层真实分布是否完全定义。当底层边际分布是连续的时候，copula函数是有唯一定义的。一些常用的copula模型包括高斯copulas（33）、t-copulas（11）和Archimedean copulas（21），每种都针对特定类型的依赖结构进行了定制。

当涉及两个随机变量时（n = 2），copula函数被称为双变量copula。给定概率密度函数??(??)=???(??)???的计算，双变量copula捕捉了两个变量之间的成对依赖结构，每个变量都与一个边际分布相关联。因此，我们可以使用copula密度函数c12和两个边际密度函数按照下面的方程2推导出联合密度函数。
??(??1,??2)=?2??(??1,??2)???1???2=?2??12(??1(??1),??2(??2))???1???2=?2C12(??1(??1),??2(??2))·??1(??1)·??2(??2)
f(x1,x2)=?2F(x1,x2)?x1?x2=?2C12(F1(x1),F2(x2))?x1?x2=c12(F1(x1),F2(x2))·f1(x1)·f2(x2)
(2)
换句话说，任何双变量copula密度函数都可以使用两个底层随机变量的边际和联合密度函数来推导。因此，我们可以将方程2推广到具有多个变量和相关边际分布的更高维度，如方程3所示。
??(??1,??2,...,????)=??(??1(??1),??2(??2),...,????(????))·??1(??1)...????(????)
f(x1,x2,...,xn)=c(F1(x1),F2(x2),...,Fn(xn))·f1(x1)...fn(xn)
(3)

Vine copula是一种构造方法，它将高维多变量依赖结构（n ≥ 3）分解为一系列成对的双变量copulas和边际密度函数。给定条件密度函数f(x1|x2) = ??(??1,??2)??(??2)，我们可以基于单变量条件密度函数和双变量copulas构建包含多个变量的更高维vine copula函数。假设示例中包括三个变量，应用条件密度函数方程4，
??(??1|??2,??3)=??(??1|??2,??3)·??(??2|??3)·??(??3)
f(x1,x2,x3)=f(x1|x2,x3)·f(x2|x3)·f(x3)
(4)
其中条件项可以通过方程5和6进一步推导。
??(??1|??2,??3)=??1,2|3(??1|3(??1|??3),??2|3(??2|??3))·??(??1|??3)=??1,2|3(??1|3(??1|??3),??2|3(??2|??3))·??(??1|??3)=??1,2|3(??1|3(??1|??3),??2|3(??2|??3))·??1,3(??1(??1),??3(??3))·??(??1)
f(x1|x2,x3)=c1,2|3(F1|3(x1|x3),F2|3(x2|x3))·f(x1|x3)=c1,2|3(F1|3(x1|x3),F2|3(x2|x3))·c1,3(F1(x1),F3(x3))·f(x1)
(5)
??(??2|??3)=??2,3(??2(??2),??3(??3))·??(??2)
f(x2|x3)=c2,3(F2(x2),F3(x3))·f(x2)
(6)
鉴于条件于X3 = x3的(X1, X2)分布对应的copula密度表示为c1,2|3，我们可以通过双变量copulas和边际分布推导出联合密度函数，如方程7所示。
??(??1,??2,??3)=??(??1)·??(??2)·??(??3)·??1,3(??1(??1),??3(??3))·??2,3(??2(??2),??3(??3))·??1,2|3(??1|3(??1|??3),??2|3(??2|??3))
f(x1,x2,x3)=f(x1)·f(x2)·f(x3)·c1,3(F1(x1),F3(x3))·c2,3(F2(x2),F3(x3))·c1,2|3(F1|3(x1|x3),F2|3(x2|x3))
(7)
通过关注边际分量的联合分布，vine copulas保留了变量之间的依赖性，同时允许我们抽象出各个边际特征。由于选择哪些变量进行条件化直接决定了在方程7中应用哪一组双变量copulas，从而决定了变量之间的潜在依赖结构。因此，如果多变量函数中的变量顺序发生变化（例如，假设分解函数为f(x1, x3, x2)或f(x3, x2, x1)，则copula结构的灵活性非常大。可能的组合枚举影响了哪些变量对应该直接耦合，条件函数的选择以及特定的双变量copula函数的选择，特别是在推广到更高维度时影响更大。因此，变量顺序的选择在确定数据中的强相关性以最大化依赖性方面起着关键作用，是该方法的关键超参数之一。

与将单个多变量copula拟合到高维数据集相比，vine copula由双变量copulas构建，随着变量数量的增加，它以更好的灵活性和可处理性分解联合分布来估计复杂的相关性。此外，与其他假设对称或椭圆依赖性的多变量Copula不同，vine Copula能够更好地捕捉复杂的、非线性的以及尾部依赖性（极值事件），而这些往往是传统相关性指标所忽视的。（9）在这项工作中，我们使用vine Copula来确定Copula家族和结构，从而最大化依赖性并保持变量间的强相关性。

2.2. 机会约束编程（Chance Constrained Programming, CCP）
机会约束编程（CCP）是一种基于优化的技术，它通过制定概率约束来模拟决策问题中的不确定性。（30）这种方法通过允许概率约束在预定可靠性水平下得到满足来解决不确定性。（10）单个概率约束，称为单个机会约束（Individual Chance Constraint, ICC），规定了该约束必须在一定的置信水平下得到满足，该置信水平用概率项（1?α）%表示，其中α代表给定约束的风险水平。（10）这种概率解释使得在严格可行性和约束遵循的灵活性之间可以进行 controlled trade-off，这在动态或不确定的环境中特别有用。例如，将风险水平设定为α = 0.05意味着该约束预计至少有95%的概率成立，同时对潜在违反该约束的情况容忍5%的风险。ICC的最一般形式如下所示为方程8a。每个ICC也可以用不确定性的逆累积分布函数（CDF）来等价表达，如下所示为方程8b。
??(????(??,??)≤0)≥1?????,???∈?
P(gi(x,ξ)≤0)≥1?αi,?i∈I
(8a)
inf{??:?(????(??,??)≤??)≥1?????}≤0,???∈?
inf{z:P(gi(x,ξ)≤z)≥1?αi}≤0,?i∈I
(8b)
在一般形式方程8a中，决策向量x和随机向量ξ位于约束函数g(x, ξ)内，同时还有违规概率α，这表明每个约束都必须以其自身的概率水平得到满足，以确保每个单独情景或实现中的指定需求要求约束的概率不小于（1?αi）%。每个重新表述的逆CDF约束，即分位数约束，强制随机变量的1?αi分位数不超过零，从而确保原始的概率要求得到满足。

联合机会约束（Joint Chance Constraint, JCC）是ICC的一种变体，其特点是在单个给定置信水平下同时考虑多个约束。类似地，JCC的数学表述可以表示为一般形式方程9a，并可以重新表述为逆CDF（分位数）约束方程9b，如下所示：
??(????(??,??)≤0,???∈?)≥1???
P(gi(x,ξ)≤0,?i∈I)≥1?β
(9a)
inf{??∈?:?(max??∈?????(??,??)≤??)≥1???}≤0
inf{z∈R:P(maxi∈I?gi(x,ξ)≤z)≥1?β}
(9b)
其中逆CDF重新表述方程9b规定最大约束违规随机变量的1?β分位数必须小于或等于0。因此，JCC方法对解空间的可行区域施加了更严格的限制，使得实际解决问题时通常更具挑战性。

一般来说，CCP提供了一种基于风险水平的简化不确定性估计方法，与随机规划等替代方法相比更为高效，因为它需要显式表示所有情景，并且在处理高维问题时会产生计算负担。在这项工作中，我们使用Copula来捕捉ICC实体之间的依赖性，并在各自的风险水平上执行约束。选择ICC而不是JCC的原因是，ICC的选择通过允许不同的需求遵循不同的不确定性水平来提高灵活性。通过分解为边际或低维概率，ICC避免了计算高维非线性联合概率的需要。因此，优化过程在计算上更为可行，同时通过Copula仍然保留了它们相关的依赖行为。

2.3. 基于Copula的数据驱动机会约束编程（Copula-based Data-driven Chance Constrained Programming）
由于历史数据的概率分布存在不确定性，构建和解决机会约束问题本质上具有挑战性。当只有有限的观测数据可用时，这些分布往往是未知的或难以准确建模的。基于现有参数分布（例如高斯分布）估计未知分布在高维或非线性情境中变得具有挑战性，因为这些分布没有封闭形式的表达式，通常需要近似或数据驱动的实证方法来生成能够更准确地捕捉实际分布的估计器。当考虑到相关性时，这一挑战进一步加剧，因为需要同时考虑边际概率分布之间的依赖结构以确保对未知分布的准确建模。我们提出了基于Copula的个体机会约束编程（Copulas-based Individual Chance Constrained Programming, C3IP），这是一个将Copula集成到数据驱动机会约束优化中的框架，以建模变量之间的依赖性，如图1所示。
图1
图1. 所提出的基于Copula的数据驱动C3IP框架的逻辑流程

将Copula技术集成到C3IP的建模中提供了一种基于实证数据的方法。通过这种方式，我们利用不确定变量的边际分布及其依赖结构来构建联合分布，从而考虑了各个分布之间的相关动态。该框架首先使用pyvinecopulib（38）这个Python接口来拟合数据中的Copula，它涵盖了一类依赖模型和双变量Copula构建块，用于建模和模拟复杂的多变量依赖性。在第一步中，应用了最大似然估计来识别最兼容的依赖结构并估计最合适的Copula家族及其相关参数。目的是捕捉强烈的成对双变量Copula依赖性，以最大化所使用的vine Copula中的整体依赖性。随后，我们基于拟合的Copula算法构建一个替代模型，模拟相关的伪观测值。在完成Copula拟合后，从拟合的Copula模型生成合成样本。拟合的Copula生成的伪观测值是从样本数据中得出的，并转换为均匀分布，同时保留了原始数据的排序，从而保持了变量之间的原始依赖结构。转换过程通过重构具有推断出的依赖结构的完整数据集来提高可扩展性和鲁棒性，同时去除了边际分布。之后，可以使用非参数方法根据实证数据估计未知的真实连续概率分布函数。在我们的方法中，我们调整了pyvinecopulib中的内置函数来根据底层边际分布推断联合CDF。最后，一旦拟合了vine Copula和边际分布，就可以通过方程1结合拟合的Copula和边际分布来构建联合分布。随后可以利用它来近似分位数函数（逆CDF）及其相关的风险水平。算法1总结了所采用的步骤。

在上面概述的数据驱动Copula拟合和模拟过程的基础上，我们通过蒙特卡洛（Monte Carlo, MC）模拟来评估其性能。在这个过程中，从基于Copula的分位数函数生成的每个情景都作为问题实例的输入，允许事后计算所有实现中约束满足的概率。

2.4. 通过DOMINO的数据驱动混合整数非线性双层优化（Data-Driven Mixed Integer Nonlinear Bilevel Optimization via DOMINO）
所提出的C3IP-DOMINO框架是一种用于双层优化问题的通用解决方案方法，其特点是具有层次化的决策结构，其中上层问题包含作为其约束部分的下层优化问题。在双层集成规划和调度（iPS）问题的背景下，利用历史需求数据来计算与特定风险水平相关联的Copula估计需求集，遵循图1中所示的程序。然后可以使用DOMINO框架来解决iPS优化问题（12,14,17）中的需求集。DOMINO（5）是一个用于解决通用约束双层问题的数据驱动框架（36,37,45），这一主题已经引起了重要的研究兴趣（3,41,48）。在这里，我们简要概述了DOMINO如何处理此类问题以及参与优化过程的求解器的具体角色，使用了在（12）中介绍的规划和调度制造过程以及（40）中介绍的原油精炼过程中的两个具体iPS问题。双层iPS问题的数学形式可以在附录A中找到。该框架通过采样上层决策变量（例如规划目标），使用确定性求解器（如BARON，44）解决相应的底层问题（例如调度）以达到全局最优，然后使用这些评估来构建一个灰盒优化模型。然后使用数据驱动的优化器，如局部搜索算法NOMAD（Nonlinear Optimization by Mesh Adaptive Direct Search, 29）或全局算法来搜索最优的上层决策。这个采样、底层优化和数据驱动的上层优化的循环不断迭代进行。该过程持续到满足收敛标准为止，这些标准包括达到最大样本数量、多次迭代后没有改进或CPU时间限制。有关DOMINO框架的全面描述，请参阅（5）。请注意，我们使用数据驱动优化作为解决方案策略，与其他此类数据驱动优化方法一样，它不能保证总是在具有有限计算预算的实证数据集中找到可行的解决方案。尽管如此，一个重要的区别是，每当找到可行的解决方案时，其可行性都是由DOMINO求解器保证的。从这个意义上说，虽然不能事先保证全局发现，但任何报告的解决方案都是可证明的可行。

图2展示了构成C3IP-DOMINO框架的Copula和DOMINO之间的工作流程。上方的子图描绘了传统的ICC模型，其中不确定性通过假设所有需求都是独立且不相关的个别风险阈值来表征。相比之下，下方的子图展示了C3IP-DOMINO模型，在该模型中通过Copula估计参数相关性。所有情景首先在多元正态分布下生成，然后分成两组：训练组和评估组。训练组用于拟合Copula，以捕捉相关变量的联合分布。一旦Copula模型训练完成并且分布估计出来，就应用分位数估计来重新表述机会约束。评估组与训练集无关，随后在案例研究中使用，特别是蒙特卡洛模拟中，来评估和展示所提出框架的性能。

3. 案例研究
3.1. 实验设置
我们在过程系统工程文献中的问题上测试了我们提出的C3IP-DOMINO框架，特别关注来自多产品制造过程（13）和原油精炼调度问题（26,39）的三个案例研究。我们使用DOMINO框架来解决双层优化问题。具体来说，集成问题的非线性方面源于序列依赖的切换以及混合操作产生的双线性项。激励示例使用集成Copula模型制定了CCP，对于其余两个案例研究，我们进一步将模型集成到DOMINO框架中，通过数据驱动的方法来解决由此产生的双层MINLP。在多产品制造和原油精炼案例研究的iPS问题背景下，不确定性出现在规划层面，并与不同情景下的需求变化相关。此外，在原油精炼过程中，不确定性也出现在调度层面，表现为容器到达时间以及每个储罐中每种组分的最低、最高和初始浓度水平。
下面详细介绍了实验的具体问题和数据设计。对于激励示例和案例研究1，我们考察了一个包含连续多产品过程中需求不确定性的iPS问题。具体来说，我们考虑了一个在连续搅拌罐反应器（CSTR）中进行的连续时间、单阶段生产过程，涉及六个不同的产品，时间跨度为两周（12,39）。对于案例研究2，我们考察了原油调度问题（26）。在整个实验过程中，系统地生成了两组实现。每个实现对应于一个独立的数据运行，由所有产品和时期的固定需求场景定义。第一组实现用于训练Copula模型，使其能够估计经验CDF并构建联合CCP所需的分位数函数，从而将每个需求量映射到相关的风险水平。第二组实现结果作为基准，用于评估蒙特卡洛模拟过程中需求满足率的性能。作为示例，生成了5000个实现结果来训练在常规藤蔓copula配置（R-vine）和0.05风险水平下的copula模型。对于这两个案例研究，我们应用了包含10,000个实现结果的更大集合，以及相同的R-vine copula类型和0.01的风险水平。对于涉及DOMINO的两个案例研究，我们首先定义所需的风险水平，并使用预训练的copula模型生成所有产品的相应估计不确定需求，这些需求被输入到DOMINO框架中，以优化基于输入需求的数据驱动生产水平的确定过程，并通过基于模拟的评估过程评估需求满足率。在DOMINO框架内，案例研究在配备Red Hat Enterprise Linux 8.9（Ootpa）的高性能计算（HPC）机器上执行了10次。DOMINO框架中的NOMAD算法是串行执行的，每个节点使用一个核心和48 GB的RAM，从随机初始点开始NOMAD算法搜索。这个点使用决策变量允许范围内的随机生成的值初始化搜索过程的种子，这些值不一定需要位于可行区域内。

对于C3IP-DOMINO框架的计算时间，我们将离线计算时间定义为包括机会约束的构建和copula拟合的时间，在示例中评估了多种样本大小的运行时间。离线计算还包括蒙特卡洛模拟。相比之下，运行DOMINO的在线计算时间在案例研究1和2中分别进行评估。对于这些案例研究，我们没有明确报告单个较低层次（调度）解决的计算时间，因为它在很大程度上取决于DOMINO选择的输入样本。调度问题表现出强烈的样本依赖性，有些实例的收敛速度显著快于其他实例。因此，报告单一代表性的调度解决时间并不能提供有意义的生产性能度量。

3.2. 动机示例
在这个动机示例中，我们研究了一个iPS问题，以确定在不确定需求下每个产品的最佳生产策略。问题配置和混合整数非线性规划（MINLP）公式在(12)中给出。iPS问题基于涉及6个产品和2个计划周期（周）的生产过程，目标是在时间范围内最大化总利润，同时最小化库存和积压的惩罚。客户需求的不确定性引入了生产目标和决策者的计划。iPS问题被构建为一个具有108个二进制变量、12个整数变量、132个连续变量和总共280个约束的MINLP。

这个动机示例旨在展示使用copulas涉及内部数据优势的重要性，强调作为DOMINO输入的copula生成的生产计划在0.01风险水平下或等效于99%的时间内对于每个可能的情景实现都是可行的。在不同的方法中，(4,7,34)提出了通过双层视角研究iPS的方法。虽然将iPS问题以双层方式构建自然反映了决策层次结构，但所涉及的计算工作可能是额外的。(7)中关于不确定性的双层iPS问题的研究仍然有限，特别是涉及非凸较低层次的情况研究更少。因此，我们提出了一种新的基于机会约束的方法。因此，我们不采用(7)中提出的随机情景分析方法，该方法以迭代方式处理每个情景并且可能非常耗资源和时间，而是利用copulas来识别产品需求之间的多变量依赖性。然后我们在决策阶段之前以离线方式计算copula估计的CDF。这种方法使得在优化过程中可以即时、动态地进行需求估计，从而大幅减少总体计算时间。

在我们的实验设计中，主要目标是在生成的生产计划中实现99%的需求满足率（或等效的风险水平0.01）。为了评估我们方法的有效性，我们评估实际性能（事后）在随机环境中的表现与最初基于copula生成的生产目标假设和建立的目标可靠性水平（事前）进行比较。具体来说，我们首先根据copula估计的需求水平计算一个固定计划，然后将其设置为所有需求情景的基准生产策略。在此基础上，我们结合copula来模拟不同需求条件下的需求，生成一系列需求情景，并系统地评估模拟生产目标满足所有预定义运营情景和约束的频率。这种方法衡量了copula估计需求在处理生产计划不确定性方面的鲁棒性。

事前和事后性能之间的差异可能是由于copula方法是数据驱动和非参数的。由于我们的方法不对需求数据的数据结构或参数设置做任何先验假设，我们依赖于观察到的数据来训练最佳拟合的copula模型来估计联合和个别需求的CDF。为了检查这种差异的影响，我们进行了一项实验，旨在隔离和评估在copula模型训练过程中引入的情景大小对预测需求满足率准确性的影响程度。具体来说，我们使用联合和个别需求满足率作为关键定量指标，来评估在预定义风险水平下copula生成的需求的质量和可靠性，以及copula方法在不同情景大小下与预期结果的一致性。联合满足率表示所有需求目标同时得到满足的模拟情景的比例，或者等效地，没有积压发生的情景的比例。相比之下，个别满足率独立评估每个产品的需求，不考虑它们之间的相关性或相互依赖性。将这些单独的比率相乘可以得到独立联合满足率，这是一个假设所有需求都单独满足的情况下的度量，不考虑任何相关性。联合和个别满足率可以分别使用公式10和11计算，其中δij表示实例i中产品索引j的积压水平。指示函数在条件满足时返回值1，不满足时返回0。

图3展示了在不同历史数据集大小下训练copula所需的满足率和计算时间。结果表明，假设需求之间相关性的copula模型始终产生比独立假设下个别比率乘积更高的联合满足率，突出了copula模型捕捉需求数据集体行为的能力。此外，copula产生的更高联合满足率表明需求水平之间的相互作用可能带来额外的价值。结合这些相互作用对于提高所有产品的整体满足率可能至关重要。此外，图3还展示了从10轮不同数量生成的情景中计算的平均满足率。我们观察到，随着训练情景总数从1000增加到10,000，联合满足率始终较高。这表明机会约束copula模型在估计CDF方面具有很高的准确性，从而生成了一个固定的生产水平，最终增加了需求满足的可能性，与个别边际满足率的乘积相比。

图3
图3. 对于6个产品和2周的实例大小，联合和个别满足率（左侧y轴实线）及copula训练时间（右侧y轴虚线）与不同数量的生成情景的关系。所有结果均为10次计算轮次的平均值。

高分辨率图片
下载MS PowerPoint幻灯片

DOMINO解决集成规划和调度问题所需的平均计算时间大约为10小时，对于6个产品和2个计划周期的实例大小。这强调了以离线方式管理不确定性的重要性，正如我们提出的数据驱动copulas框架所倡导的。此外，我们提出的copula方法只需要估计一次CDF并解决模型，从而进一步提高其效率。

我们还评估了copula估计的生产序列在所有需求情景下的质量和鲁棒性，与单独处理这些情景相比。在图4a中，我们的平均积压成本为5.08，联合满足率为99.9%，而通过解决MINLP公式创建的最优解的积压成本为0.01，联合满足率为100%（见图4b）。一旦copula模型训练完成，它就被用来生成用于所有评估实现的生产序列。

图4
图4. 通过机会约束copula框架解决的100个解决方案，包括固定估计的生产序列（左）和作为单独问题的不同生产序列（右）。

高分辨率图片
下载MS PowerPoint幻灯片

3.3. 多产品规划和调度案例研究
在这个案例研究中，我们继续研究不确定需求下的集成规划和调度问题，建立在动机研究的发现之上。该案例研究的配置涉及6个产品、2个计划周期（周）和10,000个实现结果，风险水平为1%，用于机会约束copula模型C3IP及其确定性等效模型，后者是一种传统ICC公式，通过简单地将生产目标设置为来自6个产品和2周的99%最高总需求对应的需求情景来确定生产目标。然后分别将C3IP和确定性模型的估计需求提供给DOMINO框架，在那里解决集成规划和调度问题。

将copula估计的需求集输入DOMINO框架后，应用NOMAD算法解决上层（数据驱动）问题，而BARON用于解决下层问题以达到全局最优。一旦需求生成并输入DOMINO，就在DOMINO内部进行了10次运行，所有运行的结果列在表1中。单个较低层次调度问题的解决时间可以通过将总计算时间除以评估的样本数量来获得。根据这个计算，传统ICC模型的平均解决时间约为每样本轮次4秒，提出的C3IP模型为4.3秒。图5展示了10次运行中最佳情景的生产目标和需求水平结果的具体结果。根据结果，2周内所有产品的库存水平均为零。

图5
图5. DOMINO解决的6个产品2周实例的最佳解决方案，包含copula生成和确定性需求水平，在10k个情景下。所有库存水平均为零，因此没有在图中显示。

高分辨率图片
下载MS PowerPoint幻灯片

表1. 机会约束copula模型（C3IP）与其确定性等效物（传统ICC模型）的性能比较
传统ICC模型（99%）
C3IP模型
运行时间（小时）
利润（千元）
样本时间（小时）
利润（千元）
样本
1
12.5
45.2
1996
39.6
14.9
385
88
28.4
24.0
87
560
10.3
34.6
483
94
36.5
75
60
10.1
64.9
182
124
7.9
44.8
270
98
64.6
169
73
58.8
33.9
79
56
15.0
74.9
41
3,217
68.5
5.0
577
95
7.7
4.9
46
500
76.8
65.2
162
60
6.8
14.9
356
678
7.7
5.2
168
89
9.0
4.9
378
99
7.5
5.2
167
526.9
34.9
359
69
10
11.2
75.2
105
65
18.2
64.9
3
13.5
43
平均
8.6
24.9
276
68.8
10.2
54.8
78
496.20
标准差
1.7
90.4
61
44
3.2
23.4
80.1
226
25.04

实际上，风险规避的需求估计方法通常与更高的生产目标相关联。这是因为风险规避倾向于高估需求概况，这可能导致在潜在积压的情况下增加生产水平。其背后的原理是，通过预期更极端的情景（其中需求水平显著更高），决策者可以减少生产不足的风险，避免需求违规，从而提高整体需求满足率。然而，过度高估生产目标以防止积压可能会产生额外的库存成本。此外，虽然这种策略可能提高某些产品组的短期韧性，但它可能会损害生产计划的鲁棒性和适应性。当总生产时间或资源预算有限时，这一点尤为明显。当时间预算非常紧张时，这些限制就尤为明显，因为某些产品类型可能由于其他产品的过度生产而无法按时生产。从表1可以看出，在99%的置信水平下，确定性需求情况下当前的平均生产目标更高，因此平均销售利润也更高。所提出的C3IP模型作为一种风险厌恶型方法，虽然名义上的利润可能会略有降低，但其目的是为了在生产目标和销售方面规避风险。尽管如此，该模型专门设计用于捕捉产品需求之间的相关性，并生成更现实和平衡的生产目标，从而减少过度生产的风险以及随之而来的库存成本增加。因此，由copula生成的需求通常低于其确定性对应值，更准确地反映了需求状况。表1中的结果表明，尽管copula估计的需求较低，导致目标值低于99%置信水平下的确定性模型估计值，但利润差异很小。此外，我们观察到copula生成的需求导致目标值的波动性更小，表明在不同运行中需求状况的估计更为稳定。因此，copula技术的整合增强了模型的鲁棒性及其在现实世界应用中的实用性。

接下来，我们通过蒙特卡洛模拟来评估这两种模型的性能，以研究我们提出的C3IP模型如何缓解需求的高估，从而在保持有希望的利润水平的同时减少 penalty。蒙特卡洛模拟从生成一组基于多元正态分布假设的样本情景开始。具体来说，使用了5000个情景来训练C3IP模型。这为copula依赖结构和机会约束中使用的尾部分位数提供了稳定的估计。训练过程遵循图1中详细说明的步骤，这些步骤指导了copula家族和参数的选择和调整，以适应各种情景并确保模型的鲁棒性。基于生成的训练情景，所提出的C3IP模型构建了一个95%分位数的需求集，该需求集在2周的时间内满足所有产品的预期需求的95%置信水平。另外，我们还生成了与传统的ICC公式（称为“确定性模型”）进行比较的需求，其中生产目标是通过选择对应于所有产品和时期总需求α百分位数的需求情景来确定的。这些需求是为每个目标置信水平生成的，包括90%、95%和99%的阈值。这样做是为了使用类似的生产满意度标准来比较不同确定性模型的性能。我们将机会约束的95%分位数C3IP模型称为C3IP，并将确定性模型分别表示为ICC90、ICC95和ICC99，对应于它们预定义的置信水平。然后，将每个模型估计的需求输入到（12）中提出的计划和调度优化问题中，并作为MINLP公式求解以获得解决方案。对于每个估计的需求，相应的MINLP公式被求解到最优解，以评估目标函数。我们从解决方案中获得关键性能指标，包括预期利润、预期库存和预期积压量。这些性能指标用于量化与特定模型估计的需求集相关的总体利润和 penalty，从而评估某个模型的性能。理想情况下，每个模型的需求估计应能够代表广泛潜在情景的特征，在面对训练数据之外的实际需求实现时，能够一致地保持高总体利润。具体来说，我们生成了另一个与训练集无关的评估情景样本，其中包含多元分布的产品需求。确定2000个情景的样本量，以确保样本平均值的可靠性收敛。对于copula和确定性模型生成的每个估计需求实例，MINLP公式解决一次。MINLP公式计算出的预期总利润、库存和积压量分别作为图6中每个模型需求的虚线存储。此外，还存储了与生产目标和时间调度相对应的决策变量。然后将这两组变量固定并在评估情景中重复使用，以比较copula模型和确定性模型的性能。计算每个实现的总利润、库存和积压量水平。图6中的实线反映了每个实现的实际需求情况，并说明了随着考虑更多情景结果的发展情况。结果使用以下公式生成：CumProf(??)=CumProf(???1)+Prof(??) EV(??)=CumProf(??)???? CumProf(i)=CumProf(i?1)+Prof(i)EV(i)=CumProf(i)Ni（12），其中Prof(i)表示来自实现i的利润，而CumProf(i)表示所有先前实现的总累积利润。期望值EV是总累积利润除以考虑累积的实现次数Ni，其中2 ≤ i ≤ N且1 ≤ Ni ≤ N。积压量和库存 penalty的期望值使用基于方程12的相同方法计算，其中累积利润和迭代利润被替换为相应的积压量和库存。对于每个实现，目标值计算为利润减去积压量和库存的 penalty之和。为了评估所有模型的预测准确性，我们通过比较基于估计需求的MILP模型得出的目标值与每个实现的实际需求之间的差异来跟踪所有模型的差异。

我们在图6中观察到，当考虑足够的评估情景时，copula和确定性模型的性能趋于收敛，实线在所有图表中都趋于平稳。通过从期望值中减去蒙特卡洛模拟生成的实现值，我们得到了预期结果与实际结果在目标函数之间的差异。这种差异可以在图6a中看到，以突出显示性能差距。图表显示，C3IP模型在评估情景中产生的干运行期望值（虚线）与实际利润实现（实线）之间的差距要小得多，表明与确定性模型相比，C3IP模型的需求预测更准确，即使是在99%的更高置信水平下也是如此。此外，copula估计的预期生产目标通常产生的总需求值低于其确定性等效目标，主要是由于利用了数据间的相关性。由于这一优势，C3IP模型产生的库存水平更低，因此与确定性模型相比，所承受的库存成本也更低，如图6c所示。在高置信水平下，确定性模型生成的预期生产目标更可能导致过度生产，因为大多数情况下的实际需求可以得到充分满足，从而导致更高的库存成本。从图6d可以看出，copula估计的生产目标在所有评估情景中产生的预期积压量最小。这一结果可以归因于确定性模型的内在局限性。对于置信水平较低的确定性模型（ICC90），估计的生产目标可能不足以满足大多数情景的实际需求，从而导致显著的积压量。相反，置信水平较高的确定性模型（ICC99）虽然产生了更高的预期生产目标，但由于每个时期的固定生产时间容量可能无法及时支持全部生产范围，因此也会产生更高的预期积压量。

总体而言，C3IP模型表现出更大的鲁棒性，并提供可靠的需求估计，生成更高的平均利润以及更低的积压量和库存 penalty。此外，C3IP模型在预期需求与实际需求实现之间的目标值差异最小。

3.4. 原油规划和调度案例研究
在这个案例研究中，我们考察了一个包括卸货船、储存罐、加注罐和原油蒸馏装置（CDU）的炼油厂调度问题，如图7所示。系统配置的完整细节在我们之前的工作中已经提供。（40）调度问题的目标是最小化总生产成本，而规划层次的目标函数是整个时间范围内的成本最小化。

这个例子涉及一个包含两种原油产品的15天系统的运行。在几个关键参数中引入了不确定性，包括：（i）每种产品的需求；（ii）船舶到达时间；（iii）每个储存罐中组分k的最低、最高和初始浓度；（iv）每个加注罐中原油混合物中组分k的初始浓度；（v）加注罐中组分k的上下浓度限制；（vi）每个加注罐中原油混合物的初始体积。我们从参考文献（40）中的原油炼油厂操作数据中选择了一组参数，并使用标准差为相应平均值10%的多元正态分布生成了框架的不确定性版本。这些不确定参数及其对应平均值的详细列表见表2。这些不确定参数是使用copula框架生成的，并被采样并输入到DOMINO框架中以评估韧性和性能。上层（规划问题）有28个决策变量和61个非线性灰箱约束。调度问题包括24个二进制变量和134个连续变量，总共有386个约束。我们使用BARON解决了下层MINLP，而利用DOMINO框架内的NOMAD求解器优化了双层结构。这种集成通过将基于copula的机会约束估计与DOMINO框架在规划和调度层面结合起来，从而有效地管理不确定性。

在这个案例研究中，我们选择了来自参考文献（40）中的原油炼油厂操作数据的一组参数，并使用标准差为相应平均值10%的多元正态分布生成了框架的不确定性版本。表2提供了这些不确定参数及其相应平均值的详细列表。这些不确定参数利用copula框架生成，并被采样并输入到DOMINO框架中以进行评估。上层（规划问题）有28个决策变量和61个非线性灰箱约束。调度问题包括24个二进制变量和134个连续变量，总共有386个约束。我们使用BARON解决了下层MINLP，而使用NOMAD求解器在DOMINO框架内优化了综合双层结构。

这个例子涉及一个15天周期内的系统运行，涉及两种原油产品。在不確定因素包括：（i）每种产品的需求；（ii）船舶的到达时间；（iii）每个储存罐中组分k的最低、最高和初始浓度；（iv）每个加注罐中原油混合物中组分k的初始浓度；（v）加注罐中组分k的上下浓度限制；（vi）每个加注罐中原油混合物的初始体积。我们从参考文献（40）中选择了原油炼油厂操作数据的一组参数，并使用多元正态分布生成了框架的不确定性版本，标准差为相应平均值的10%。表2提供了这些不确定参数的详细列表及其相应的平均值。这些不确定参数是利用copula框架生成的，并被采样并输入到DOMINO框架中以评估韧性和性能。上层（规划问题）有28个决策变量和61个非线性灰箱约束。调度问题包括24个二进制变量和134个连续变量，总共有386个约束。我们使用BARON解决了下层MINLP，而使用NOMAD求解器在DOMINO框架内优化了综合双层结构。

图7展示了原油装卸和精炼网络的配置。该图改编自参考文献（40）。本文根据知识共享署名（CC BY）许可协议发布。

我们使用CCP模型根据下表中给出的上下限生成不确定参数。我们考虑了10,000个情景和0.01的风险水平。C3IP-DOMINO框架被运行10次，以生成目标值和样本大小。我们在高性能计算（HPC）上设置了144小时的上限计算时间。结果分别展示在表3和图8中。在表3中，双层问题的目标函数是在较低调度层最小化生产成本，在较高规划层最小化整个时间范围内的规划成本。报告的目标值是上层规划成本，其中包括调度层的目标（生产成本）。样本指的是在优化过程中生成和评估的试验点数量。相应的一轮解决时间分别为：传统ICC模型为157.2秒，C3IP模型为168.5秒。

图8：通过DOMINO解决的双产品15天实例的最佳解决方案，该实例在10k种情景下具有由copula生成和确定性的需求水平。高分辨率图像|下载MS PowerPoint幻灯片

表3：Chance Constrained C3IP模型与其99%分位数确定性等价模型（使用ICC）的性能比较

| 传统ICC模型（99%） | C3IP模型 |
| --- | --- |
| 运行时间（小时） | 目标值（k） | 样本时间（小时） | 目标值（k） | 样本 |
| 1 | 144 | 32.69 | 144 | 57.11 |
| 2 | 144 | 57.11 | 144 | 58.58 |
| 3 | 144 | 164 | 144 | 50.77 |
| 4 | 144 | 318 | 144 | 44.20 |
| 5 | 144 | 653 | 144 | 72.60 |
| 6 | 144 | 185 | 144 | 81.98 |
| 7 | 144 | 93.8 | 144 | 345 |
| 8 | 144 | 94.8 | 144 | 33.97 |
| 9 | 144 | 593 | 144 | 23.85 |
| 10 | 144 | 671 | 144 | 50.77 |
| 11 | 144 | 128 | 144 | 58.63 |
| 12 | 144 | 189 | 144 | 138.09 |
| 13 | 144 | 968 | 144 | 50.72 |
| 14 | 144 | 69.96 | 144 | 187 |
| 15 | 144 | 469 | 144 | 46.43 |
| 16 | 144 | 377 | 144 | 144 |
| 17 | 144 | 186 | 144 | 52.75 |
| 18 | 144 | 345 | 144 | 611 |
| 19 | 144 | 346 | 144 | 286 |
| 20 | 144 | 507 | 144 | 189 |
| 平均 | 144 | 3.71 | 30 | 1.50 |
| 标准差 | 0.20 | 20.23 | 0.39 | 0.24 |

注：目标值是规划成本，数值越低越好。

总体而言，在CCP计划下预测的需求水平低于确定性计划，这对应于两种产品的生产水平都较低。特别是在第7次运行中，我们观察到一个异常值，这是由于只抽取了96个样本，且由于DOMINO初始化问题，GAMS提前终止了搜索过程，因此达到了144小时的时间限制。第5次运行也观察到了类似的情况，其中99%分位数确定性等价模型的DOMINO返回的规划成本远高于Chance Constrained C3IP模型。此外，基于数据驱动的公式化方法的表现优于基于情景的方法，因为后者扩大了数据驱动优化器的搜索空间。我们注意到，与案例研究1中的规划和调度问题相比，原油问题通常更具挑战性，并且对DOMINO中的搜索初始化更加敏感。

在图8中，我们观察到在确定性计划下，产品1在第3天和第11天的生产水平较高，导致第4天和第11天的库存水平上升。相反，在CCP计划下，产品1仅在第10天的生产水平较高，从而导致第11天的库存水平达到峰值。此外，预测第6天的需求水平在CCP计划下低于确定性计划，而第7天的需求水平则较高。这导致CCP计划下第7天的库存水平更高。

**结论**：
我们提出了一种带有copula函数的机会约束优化框架，有效捕捉和模拟了多变量需求数据的依赖结构，并将其应用于集成规划和调度问题。我们的方法确保了在预定义的风险阈值内进行决策是可行的。我们使用DOMINO框架验证了由copula生成的生产目标，以解决由此产生的数据驱动的双层最小化问题（MINLP）。我们的计算实验表明，所提出的方法能够识别出稳健的解决方案，这些解决方案在提高产品联合满意度方面优于传统的单个机会约束公式化方法，并且性能接近最优，同时保持了合理的计算时间。