为什么蒙特卡罗?

蒙特卡罗方法在任何机械工程师的工具包中都是一种非常强大的计算技术。在机械工程中有许多常见的用例，蒙特卡罗模拟自然适用于这些用例，包括:

  • 统计公差分析
  • 可靠性设计
  • 测量不确定度分析
  • 机器人误差、位置和姿态分析
  • 优化设计

在本文中，我们将专注于构建对该方法的理解，建立一个可以应用于各种问题的分析框架，并将该框架应用于机械工程师可能经常面临的常见统计公差分析问题。

方法概述

蒙特卡罗方法的基本概念是使用随机性来解决确定性问题。例如，如果我掷三个六面骰子，骰子之和会大于七的频率是多少？问题本身是确定性的，这意味着它只有一个答案。无论骰子掷到哪里，掷到谁手中，骰子的行为都是一样的但是，也有随机因素。在掷骰子之前，你无法预测任何一个骰子会落在什么地方。

然而，很容易看到一个骰子的结果，并检查我们在问题陈述中列出的条件(它们加起来是否大于7 ?)。

这就是蒙特卡罗方法的魅力和力量。在许多问题中，直接计算答案是极其困难的，但检查“实验”(在本例中是掷骰子)的结果是否微不足道。这种方法依赖于进行足够多的实验或试验，这样计算每次试验的结果就能在一般意义上告诉我们最初问题的答案。

分析框架基础

概率密度函数（PDF）是任何蒙特卡罗模拟的核心组件。PDF完全描述了正在模拟的任何随机变量的行为。在掷骰子示例中，PDF很容易理解。当你掷一个骰子时，有一个相等的骰子(⅙) 任意数被滚动的概率。此均匀概率密度函数如下图所示。

为了构建有用的蒙特卡罗模拟，必须使用一个具有代表性的PDF，该PDF准确地捕获真实世界中每个试验的行为。与骰子不同，机械零件和其他工程过程的尺寸表现出连续的可变性。对于连续随机变量，即使在有限的范围内，也有无限个可能的值。因此，您需要一个诸如高斯(或正态)分布的函数来准确地表示连续的随机变量。

示例PDF：高斯函数

您的分析所需的确切PDF取决于您试图表示的现象，可能需要一些其他分配除了高斯(如幂律、威布尔、伯努利、二项式等)

分析框架

现在我们已经了解了试验的核心组件(概率密度函数)，我们可以围绕它构建一个分析框架。下面概述了框架中的五个步骤。这些步骤将使用伪代码以及谷歌Sheets、MATLAB和python函数的特定引用进行解释，您可以在自己的应用程序中使用这些函数。只要您理解框架并准确地应用它，您就可以在电子表格中或使用选择的编程语言构建有用的模拟。

  1. 构建pdf文件.将所有变量表示为概率密度函数，离散或连续。（例如：多个正态分布，每个正态分布具有各自独立的平均值和标准偏差值）
  1. 构造函数从这些分布中取样。在谷歌表格中，您可以结合使用NORMINV()函数和RAND()进行随机抽样。在python中，你可以使用np.random.normal从numpy包。
  1. 构造一个试用检查函数。此函数用于查看单个模拟试验，并根据某些确定性结果检查该试验的结果。例如，三个骰子的总和是否大于七（对还是错）？此函数（或多个函数）应聚合并存储您的单个试验结果，以供以后使用。在excel或g-sheets中，这通常意味着构建一个试验结果列。在python或MATLAB中，您应该构造一个列表、元组或数组（或其他数据结构），用于存储试验结果。
  1. 构造一个仿真结果检查函数。这个函数以这样一种方式提出了您的工程问题，它可以通过对个别试验结果进行一些计算来回答。这是你真正发挥创造力的地方，也是方法灵活和强大的地方。想出创造性的方法将试验数据与工程问题联系起来，正是这种方法能够在如此多的领域发挥作用。
  1. 使用下面的伪代码作为指导构造主模拟。
#从相关的概率分布函数中提取样本#使用您的试验检查函数检查样本#存储每个试验检查的结果#运行模拟结果检查函数来计算您的模拟结果

笔记:有许多其他的模拟框架，你可以使用蒙特卡洛模拟。一些简单的模拟受益于在主模拟循环中运行结果检查函数，每次试验后增加结果变量。一个框架相对于另一个框架的适用性取决于所解决问题的性质和复杂性以及您的个人偏好。

典型的收敛图显示随着试验计数的增加误差减少

公差分析的例子:O形环压缩问题

假设我们正在设计一个液压活塞和气缸执行器，使用丁腈橡胶o形环密封活塞(见下面的原理图)。

O形圈设计为在一定的干涉（压缩）范围内工作，我们称之为X。气缸、活塞和O形圈的尺寸和制造公差如下。

o形圈直径= 3mm±0.09 mm = DO

气缸内径= 25mm±0.1 mm = DC

活塞环槽直径=22.5mm±0.03 mm = DP

为了使活塞密封正常工作，o形圈压缩必须保持在0.3mm和0.6mm之间。有了上述制造公差，执行机构将不能正常工作的百分比是多少?

解决方案

使用我们的分析框架，我们必须首先选择一个PDF来表示每个加工过程和装配的关键尺寸。我们将使用一个正态分布，其中平均值是目标尺寸，而标准偏差是由指定的加工公差导出的。对于这些类型的公差问题，一个常见的假设是，在99.7%的零件上达到设计中规定的公差，或三个标准差(西格玛)的覆盖范围。

在为公差堆栈中的每个维度选择pdf文件之后，我们将对每个分布进行抽样，以便为每个模拟试验“构建一个组件”。最后，我们将根据我们的功能o形环规格检查每个组件，并汇总结果。要查看Python代码，请展开下面的部分。

Python示例代码 
导入numpy作为np导入matplotlib。pyplot plt n = 1000000 # #组件规范σ= 3 #的试验假设99.7%的公差组件规范dp = 22.5 #活塞直径tol_dp = 0.03 = 3 # #活塞直径公差做o形圈直径tol_do = 0.09 # o形圈直径公差dc = 25 #缸直径tol_dc = 0.1 #缸直径公差defcalculateInterference (dp,特区):“”:直径:o型环干涉”“干扰= dp +做-直流干扰def checkInterference回报(干扰,低,上):“”:o形圈干扰:Bool,通过T, F为失败”“如果干扰<低或干扰>上:返回错误信息:还真#设置干扰限制低= 0.3 #最小容许干扰上= 0.6 #最大容许干扰interference_list =[] #初始化列表存储每个试验结果的干扰值=[]#初始化列表来存储我们的审判结果#主仿真循环试验的范围(n):# Sample PDF for each component piston_sample = np.random。正常(dp, tol_dp/sigma) oring_sample = n .随机。Normal (do, tol_do/sigma)正常(dc, tol_dc/sigma) #计算和存储干涉干涉= calculateInterference(piston_sample, oring_sample, cylinder_sample) interference_list.append(interference) #试验检查:记录干涉检查结果的结果。append(checkInterference(interference, lower, upper)) # results check goodAssemblyCount = sum(results) failurePercentage = 100*((n - goodAssemblyCount) / n) print(' failed piston assembly: '， failurePercentage， "%")

我们的100万次模拟试验显示，活塞装配失败率为1.4738%。

然而，我们不知道故障的分布情况，也不知道如何修改设计以防止这些故障。为了可视化o形环干涉的分布，我们可以用下面的python代码绘制模拟数据的直方图。

Python示例代码 
#设置直方图参数，绘制n_bins = 500 heights, bins， _ = plt。Hist (interference_list, bins=n_bins) # get position and heights of bars bin_width = np.diff(bins)[0] bin_pos = bins[:-1] + bin_width / 2 # mask for coloring failures differently maskleft = (bin_pos <= lower) maskright = (bin_pos >= upper) # plot data in three steps plt. txt . txt4)图(figsize =(14日)plt。text ='#281E78'， label='Good Assemblies'， text ='好程序集'，text ='好程序集'plt. Bar (bin_pos[maskleft]， heights[maskleft]， width=bin_width, color='#ea4228'， label='left tail failures')bar(bin_pos[maskright]， height [maskright]， width=bin_width, color='#49C6E5'， label='right tail failures') #轴标签和右尾失败的垂直标记xlabel(“O形环干涉(mm)”)ylabel plt(“#样品”)。Axvline (low, linestyle='——'，linewidth=1, color='grey')axvline(upper, linestyle='——'，linewidth=1, color='grey') # remove上方和右侧的spine, set ticks ax = plt.gca() ax.spine ['right'].set_visible(False) ax.spine ['top'].set_visible(False) ax.xaxis.set_ticks_position('bottom') ax.yaxis.set_ticks_position('left') plt.legend()

从o形环干涉直方图上可以看出，所有的故障都发生在分布的右尾，意味着干涉过多。

假设在分布的左尾部有空间，一个简单的解决方法是调整堆栈中的一个维度，将组装分布向左移动。我们可以通过将活塞环槽直径从22.5mm减小到22.45mm来进行模拟，同时保持所有公差相同。通过对o形环干涉分布在左右尾端截止阀之间的定心，模拟该工况可将失效率降低至0.1118%。如下图所示

中心o形圈干涉分布
正态分布中的左尾失效细节

外卖

我们希望你能从这篇文章中学到的主要内容是，蒙特卡罗方法并不难实现，而且对于各种常见的机械工程问题可以非常灵活和有用。了解基本框架允许您在几分钟内快速创建有用的电子表格模拟或python代码。在评估模拟的准确性时，记住以下关键事实是很重要的:

  • 概率密度函数必须能代表现实世界中实际的系统行为。
  • 你从概率密度函数中抽取的样本必须是随机样本，否则你的模型就不能准确地代表总体。
  • 大数定律表明，随着样本(试验)计数的增加，模型的准确性也会增加。
  • 对于需要高度精确答案的模型，应增加试验计数，直到模拟的准确性不再改变。检查收敛性并不能保证模型是准确的，但是通过允许您平衡计算成本和模型精度的边际增长，它确实产生了有效的结果。