考虑不确定性的数据驱动的动力学序列线性规划方法

1.本发明涉及数据驱动计算动力学技术领域，尤其是涉及一种考虑不确定性的数据驱动的动力学序列线性规划方法。


背景技术：

2.随着现代科学技术的发展，力学已经成为不同工程领域强有力的理论支撑，特别是经典的数值计算力学方法已经成为分析结构应力、位移响应强有力的工具，为实际工程提供理论指导。
3.经典基于本构关系的计算动力学方法(constitutive low-based computing in dynamics,简称clcd)实际上是建立在基于显式本构模型的框架上的。该框架根据从特定实验中获得的有限数据点(通常是在简单的拉/压/纯剪切加载状态下)建立相空间中状态变量之间的显式本构函数，然后利用该显式本构函数来驱动数值求解。但是显式本构关系的获取需要长期经验的积累和极高的人力物力成本，因此现有的本构关系常常以人名命名，例如hook定律、mooney rivlin模型等；而且目前不存在能够表达所有材料的本构关系的显式函数。随着科技的发展，新型材料层出不穷，现有的显式本构已难以满足当前的需求，然而为每种新材料都构建显式本构关系已难以满足当前产品研发周期的需求。同时随着实验方法的不断进步，获取的实验数据出现了爆炸式增长。而且继实验、理论、计算之后，数据已成为人类认知世界的第四科学范式。因此，直接利用材料的实验数据建立一套全新的框架去预测结构响应引起了很多力学工作者的关注。这种直接利用离散的实验数据，绕开建立解析形式的显式本构模型的全新范式被称为“数据驱动计算力学”，特别是将这一方法引入动力学领域，已经发展成为经典的范式“数据驱动计算动力学”(data-drivencomputingindynamics简称ddcd)。
4.该类方法不对本构数据做任何的拟合，且脱离经典计算力学框架，通过求解本构数据中的点与满足守恒定律的点的距离极值寻求问题的解。尽管如今的框架下能根据数据集求解出满足数据本构以及协调方程和平衡方程的响应曲线，但由于实验数据的不确定性，材料性质的波动，此单一曲线缺乏可信度，不能很好的满足工程需求。
5.实际上，当采用数据驱动框架时，由于实验中不可避免的不确定性因素，难以确保实验数据点真的表征了材料的本构行为，进而导致在此框架下求得的单一响应曲线缺乏实际的意义。


技术实现要素：

6.本发明的目的是解决上述现有技术存在的技术问题。
7.为实现上述目的，本发明提供了一种考虑不确定性的数据驱动的动力学序列线性规划方法，具体步骤如下：
8.步骤s1：将动力学的二阶微分方程在时间域离散，形成多个时刻点，输入初始参数，并进行初始计算获得质量矩阵m和协调矩阵b，设置时间步δt、newmark算法参数以及初
始条件；
9.步骤s2：计算当前时刻等效右端载荷，为每个单元分配初始凸组合数据并得到线性规划求解；
10.步骤s3：判定线性规划求解是否可行，当线性规划求解可行时，自适应更新局部凸组合点，当线性规划求解不可行时，修正后自适应更新局部凸组合点，直到更新后局部凸组合点满足收敛条件；
11.步骤s4：通过newmark算法计算收敛后的速度和加速度，并将速度和加速度存储为当前时刻响应，返回步骤s2中进行下一时刻响应的计算；
12.步骤s5：直到计算完时间域所有时刻点的响应，输出相应曲线结果。
13.优选的，在步骤s1中，newmark算法参数包括常数参数α和常数参数β。
14.优选的，在步骤s2中，当前时刻tk等效右端载荷计算公式如下：
[0015][0016]
其中，为当前时刻等效右端载荷，pk为当前时刻外载荷向量，为当前时刻newmark预测值；
[0017]
为每个单元分配初始凸组合数据具体操作如下：
[0018]
从数据集中选择初始局部凸组合的应力应变数据：
[0019][0020]
其中，enc为构成第e个单元的第nc个数据点的标识，其中e＝1,
…
,m，m为划分单元的个数，表示从数据集d中选出构成局部凸组合的数据点；
[0021]
式中初始数据编号如下：
[0022][0023]
其中l
(1)
是(n
d-1)/nc的整数部分；
[0024]
根据每个单元分配初始凸组合数据计算得到线性规划求解和
[0025]
优选的，在步骤s3具体为：设l为更新次数，
[0026]
当线性规划求解可行时，采用如下公式获取单元应力应变状态
[0027][0028]
其中，表示为局部凸组合的系数，
[0029]
当线性规划求解不可行时，通过修正公式获取单元应力应变状态计算
公式如下：
[0030][0031]
其中，为修正系数和过渡应力应变状态参数；
[0032]
通过修正公式获取修正参数如下：
[0033][0034][0035][0036]
其中，we表示为单元体积，be表示为单元协调矩阵，为修正系数；
[0037]
通过上述得到的单元应力应变状态得到更新后的局部凸组合的数据点
[0038]
更新公式如下：
[0039][0040]
其中，t为nc/2的整数部分，l
(l+1)
＝l
(l)
+1，l
(l)
为大于等于1的整数，为距离当前状态最近的数据点的编号；
[0041]
收敛判断公式如下：
[0042][0043]
其中δ为收敛设定参数。
[0044]
因此，本发明采用上述一种考虑不确定性的数据驱动的动力学序列线性规划方法，具有以下有益效果：
[0045]
(1)局部凸组合描述本构关系的方法明确考虑了数据集中不可避免的不确定性(可能是由于测量误差、信息不足、建模不准确等原因造成的)；通过更改目标函数来求解结构响应的解集，可以衡量数据集不确定性对响应的影响，比经典ddcd方法的单一响应曲线更具有可信度，更有利于工程师的判断。
[0046]
(2)在数据点组成的凸包内求解结构响应，由于凸包的封闭性，使得求得的响应并不会远离数据点，从而使得该方法对于小样本数据集也具有有效性，降低了经典ddcd方法中对实验数据点数目的要求。
[0047]
下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。
附图说明
[0048]
图1为本发明一种考虑不确定性的数据驱动的动力学序列线性规划方法流程图；
[0049]
图2为本发明三杆桁架结构和局部凸包数据点更新数据示意图；
[0050]
图3为本发明实施例1不同目标函数下u3的位移响应曲线。
[0051]
图4为本发明三维桁架结构示意图；
[0052]
图5为本发明实施例2非线性本构关系数据集；
[0053]
图6为本发明实施例2应力应变状态与参考解误差随数据集总数变化趋势图；
[0054]
图7为本发明含噪声的非线性材料数据集(θ0＝0.02)图；
[0055]
图8为本发明2不同θ0的位移下界对应30组数据结果的直方分布图；
[0056]
图9为本发明不同θ0的位移上界对应30组数据结果的直方分布图；
[0057]
图10为本发明带有32个坏点的数据集图；
[0058]
图11为本发明0.8倍参考值坏点的位移下界对应30组数据结果的直方分布图；
[0059]
图12为本发明0.8倍参考值坏点的位移上界对应30组数据结果的直方分布图；
[0060]
图13为本发明1.2倍参考值的坏点的位移下界对应30组数据结果的直方分布图；
[0061]
图14为本发明1.2倍参考值的坏点的位移上界对应30组数据结果的直方分布图；
[0062]
图15为本发明不同噪声下的位移响应曲线；
[0063]
图16为本发明不同随机噪声数据集的上下界有效率统计图；
[0064]
图17为本发明64个不同倍数坏点的数据集得到的位移u
580
响应曲线；
[0065]
图18为本发明不同倍数异常值的上下界有效率统计图。
具体实施方式
[0066]
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。
[0067]
下面结合附图，对本发明的一些实施方式作详细说明。
[0068]
实施例1
[0069]
本实施例以二维三杆桁架，考虑无阻尼任意时刻t的动力学问题可以表示为如下优化列式，理论公式如下：
[0070]
findd＝(u(t)
t
,λ(t)
t
)
t
[0071]
mini＝s
t
u(t)(1.a)
[0072][0073][0074]
[0075][0076][0077][0078]
式中分别表示t时刻各节点的加速度和速度；表示t时刻的位移边界条件；表示t时刻的位移边界条件；表示各节点在初始t0时刻的位移、速度和加速度；d表示设计变量。s＝(0,0,
…
,
±1ith
,0,
…
,0
t
∈r
n+2m
为目标函数的系数向量，当s＝(0,0,
…
,-1
ith
,0,
…
,0)
t
∈r
n+2m
时，可求得di的上界；当s＝(0,0,
…
,+1
ith
,0,
…
,0)
t
∈r
n+2m
时，可求得di的下界。的下界。是第e个杆单元的nc个数据点的凸组合系数组成的列向量。
[0079]
求解(1)式中的微分方程本质上是积分过程，常用逐步积分法，把时间历程(0≤t≤t)划分为一系列时间段δt，根据一定的算法由前一时刻t
k-1
的位移u
k-1
、速度加速度计算tk＝t
k-1
+δt(k＝1,2,
…
,k)时刻的位移uk、速度加速度本发明采用newmark算法：
[0080][0081][0082]
其中α、β为newmark参数，通常取α＝0.25，β＝0.5。
[0083]
为了由前一时刻的位移u
k-1
、速度加速度表示下一时刻的速度加速度由(2)式得：
[0084][0085][0086]
其中其中
[0087]
对(1)式在时间域离散后可将式(3.a)中k时刻的加速度代入动力平衡方程，得到tk＝t
k-1
+δt时刻每次迭代的动力学线性规划数学列式+δt时刻每次迭代的动力学线性规划数学列式
[0088]
[0089][0090][0091][0092][0093][0094][0095]
在每个时刻迭代求解的线性规划问题唯一的不同就是这nc个数据点因此这些局部数据点必须在求解之前确定，并且会随着迭代过程中的dk自动更新。
[0096]
如图2所示，本实施例杆件长度l1＝l3＝10,横截面积a＝1；材料密度ρbar＝2.5
×
10-6
；外荷载激励其中为外荷载激励频率，t＝0.1为外荷载激励周期。由如下公式生成带有噪音的数据集，共计501个数据点，
[0097][0098]
其中，e＝1，当时，当时，θ＝0.05。为在0,1区间上的均匀分布的随机数。
[0099]
步骤s1：将动力学的二阶微分方程在时间域离散，形成多个时刻点，输入初始参数，并进行初始计算获得质量矩阵m和协调矩阵b，设置时间步δt、newmark算法参数以及初始条件。newmark算法参数包括常数参数α和常数参数β，设置初始条件
[0100]
步骤s2：计算当前时刻等效右端载荷，为每个单元分配初始凸组合数据并得到线性规划求解。
[0101]
当前时刻tk等效右端载荷计算公式如下：
[0102][0103]
其中，为当前时刻等效右端载荷，pk为当前时刻外载荷向量，为当前时刻newmark预测值。
[0104]
为每个单元分配初始凸组合数据具体操作如下：
[0105]
从数据集中选择初始局部凸组合的应力应变数据：
[0106][0107]
其中，其中，enc为构成第e个单元的第nc个数据点的标识，其中e＝1,
…
,m，m为划分单元的个数，表示从数据集d中选出构成局部凸组合的数据点；
[0108]
式中初始数据编号如下：
[0109][0110]
其中l1是n
d-1/nc的整数部分。
[0111]
步骤s3：判定线性规划求解是否可行，当线性规划求解可行时，自适应更新局部凸组合点，当线性规划求解不可行时，修正后自适应更新局部凸组合点，直到更新后局部凸组合点满足收敛条件。
[0112]
设l为更新次数，
[0113]
当线性规划求解可行时，采用如下公式获取单元应力应变状态
[0114][0115]
其中，表示为局部凸组合的系数。
[0116]
当线性规划求解不可行时，通过修正公式获取单元应力应变状态计算公式如下：
[0117][0118]
其中，为修正系数和过渡应力应变状态参数；
[0119]
通过修正公式获取修正参数如下：
[0120][0121][0122][0123]
其中，we表示为单元体积，be表示为单元协调矩阵，为修正系数；
[0124]
通过上述得到的单元应力应变状态得到更新后的局部凸组合的数据
点
[0125]
更新公式如下：
[0126][0127]
其中，t为nc/2的整数部分，l
(l+1)
＝l
(l)
+1，l
(l)
为大于等于1的整数，为距离当前状态最近的数据点的编号；
[0128]
收敛判断公式如下：
[0129][0130]
其中δ为收敛设定参数。如图2所示，在图2中：(a)三杆桁架示意图；(b)初始数据点构成的局部凸组合；(c)-(e)迭代过程中1号杆局部凸组合自动识别过程；(f)目标函数为i(u＝p
t
u收敛时，局部凸包数据点。
[0131]
当可行时，但是此时局部凸包占据的空间过大，即凸包的形状与数据集中包含的本构关系相差甚远，如图2(b)所示。因此，为了得到更紧的上下界，或者说让得到的上下界更加贴近数据集中包含的本构关系，需要在下一步中缩小局部凸包的大小，即其占据的空间。在当前的算法中，首先计算了数据集中所有数据点与当前状态的欧拉距离，即随后便可根据得到离当前状态最近的数据点在数据集中的编号因为在第一步中数据集中的数据点便以按sign(εd)||(εd,σd)||2的大小进行了排序(即根据εd正负定义的向量范数的大小)。因此算法只需要引进一个大于等于1的整数l
(l)
便能很好的衡量局部凸包的大小，缩小局部凸包时也只需要减小l
(l)
的值即可，可让l
(l+1)
取l
(l)
/ρ的整数部分，进行更新，其中ρ为初始设定的大于1的参数。一旦与l
(l+1)
确定，第l+1步的局部凸包数据点便可确定，如果nc为奇数，相应的表达式可如下所示：
[0132][0133]
其中t为nc/2的整数部分。
[0134]
相关参数分别设定为：l
(1)
＝72，ρ＝3，nc＝7，δ＝0.5。
[0135]
当不可行时，此时意味着局部凸包的位置不合适或者过小，因此相比缩小局部凸包的范围，此时更应该扩大局部凸包的范围，可以让l
(l+1)
＝l
(l)
+1或者将当前局部凸包的数据点编号向两端扩展。一旦l
(l+1)
确认，便可按上述可行解的方式更新局部凸包的数据点。
[0136]
任取一个时间步t＝0.06为例说明本文算法的迭代收敛过程。每个时间步迭代开
始前从全局数据集中按照如下规律所有单元分配局部数据集：
[0137]
构成初始凸组合，并求解线性规划得到初始的可行解，如图2b所示；随后逐步缩小构成凸组合数据点之间的间距l
l+1
＝l
l
/ρ，在第l步之后更新第l+1步构成局部凸组合的数据点，如图2b-c所示，数值算例显示，即使在每个时刻构成每根杆单元初始凸组合的数据点一致，但是随着局部凸组合的自适应更新，自动识别每个杆单元对应的应力应变所在的大致位置，并在其附近选择局部数据点构成凸组合，图2b-e便很好地显示了这一点。当l
(l)
＝l
min
时，此时已成功识别每个杆单元的应力应变，并收敛，如图2f所示，此时基本得到t＝0.06时刻解答。
[0138]
为了获取2号节点水平位移u3的在某时刻t＝0.06的上下界，分别取目标函数为i(u)＝-u3和i(u)＝u3，并取目标函数i(u)＝p
t
u的特殊解和精确本构σ＝eε的有限元解答作为参考。下表列出了2号节点水平位移u3上下界以及特殊解在迭代中变化。
[0139]
三杆桁架算例在三个目标函数下节点位移的迭代过程
[0140][0141]
可以看出，随着局部凸组合的逐渐收敛，u3的上界(i(u)＝-u3)从第一次迭代的3.0480下降到第5次迭代的3.0437。同样u3的下界(i(u)＝u3)从第一次迭代的2.1046上升到第5次迭代的2.1083。说明随着迭代过程的进行，数据集中不确定量化的置信度有了一定的提高。另外，特殊解(i(u)＝p
t
u)位移u3＝2.2176和有限元解u3＝2.5781，均落在了u3上下界所在的区间[2.1083,3.0437]内，上述过程和结果清晰地证明了本方法实现过程以及数值实现的有效性。
[0142]
步骤s4：通过newmark算法计算收敛后的速度和加速度，并将速度和加速度存储为当前时刻响应，返回步骤s2中进行下一时刻响应的计算。
[0143]
步骤s5：直到计算完时间域所有时刻点的响应，输出相应曲线结果。
[0144]
图3给出了外荷载激励作用总时间为t＝0.1，时间步长δt＝0.004时2号节点水平位移u3的上下界和相关参考解响应随时间变化曲线。从曲线可以看出，u3上界(upper)和下界(lower)构成的曲线恰好可以将以精确本构关系(σ＝eε)计算得到响应(fem)包罗在内，而以目标函数为i(u)＝p
t
u的特殊解(pu)只有一个时间节点不在上下界之间。说明本发明的有效性。
[0145]
实施例2
[0146]
本实施例以三维桁架为例，如图4所示，其中包含1194个杆单元，1002个自由度，杆的横截面积a＝1，密度ρbar＝2.5
×
10-5；在桁架底部节点施加10个周期的正弦位移激励其中um＝1为位移激励幅值，为位移激励频率，t＝0.1为位移激励周期。假设各杆件处于小变形弹性状态，精确非线性本构关系数据集如图5所示，表达
式如下：
[0147][0148]
其中，e表示为弹性模量。
[0149]
取目标函数分别取i(u)＝u
580
、i(u)＝-u
580
获得194号节点水平位移下界、上界以及i(u)＝p
t
u的特殊解，使用精确数据集进行求解，为了度量不同目标函数下数据驱动解和参考解之间的误差，提出度量整个时间域上应力应变状态之间误差error2和u
580
位移误差的公式，如下所示：
[0150][0151][0152][0153]
式中，t1为第一个时间步之后的时刻，tf为最终时刻；(εe,σe)表示本方法求解的应力应变状态，表示参考解的应力应变状态。
[0154]
分别取数据集总数为nd＝101,317,1001,3429,10001，取ρ＝2,取ρ＝2,δ＝0.1，下表给出本发明方法在不同数量数据集下个输出节点位移误差每个时间步线性规划的平均次数、不可行线性规划总次数、不可行线性规划总时间以及求解总时间。
[0155][0156][0157]
由上述可知，在精确数据集下，当数据点逐渐增大时，位移误差和应力应变误显著降低，说明本发明计算得到的上下界都趋近于参考解。这是因为数据点越密集，最终收敛时算法自适应得到的局部凸包就越小，越接近参考本构；而对于数据集总数从nd＝3429增大到nd＝10001的过程中，误差并没有显著降低，说明增大数据集并没有效降低误差，说明本发明算法更适用于稀疏数据样本。
[0158]
另外值得注意的是，局部凸组合数不变时，随着数据集增大，在从101增大到10001过程中，本发明算法总时间有了明显增加，主要有两方面的原因：一方面，在数据集nd增大过程中，ρ＝2为固定值，每个时间步的迭代次数也在增加，导致计算总时间增加；另一方面，由于局部凸组合数nc＝5为固定值，这对于较大数据集个数而言，在局部凸组合缩小过程中出现局部图组合范围过小，无法包罗上一步求解的可行解，导致该迭代步线性规划出现不可行解的，增加求解时间。为了验证这一原因，分别取ρ＝3、局部凸组合数nc＝7以及其组合的求解时间进行对比，根据表3求解时间可知，适当增大ρ和局部凸组合数nc显著降低了每个迭代步的平均线性规划次数、不可行解时间和总求解时间，特别是对于上下界的求解。
[0159]
因此，本文算法关于数据集具有收敛性，即上下界随着数据点的密集而接近于参
考解；从求解误差和求解效率来看，本发明算法更适用于稀疏数据样本，同时对于密集数据样本，选择合适的ρ和局部凸组合数可以显著降低求解时间。
[0160]
本发明方法的鲁棒性分析。
[0161]
对本构关系添加不同的随机噪声进行计算。其中，nd＝1001，当时，θ＝θ0，否则，图5所示为为θ0＝0.02，nd＝1001的数据集。计算三组不同θ0(0.01,0.02,0.04)的数据集，每组包含30个数据集，每个数据集计算中其他参数设置为ρ＝2，nc＝9，δ＝0.1。对应的位移和应力应变状态误差的平均值以及均方差如表所示：
[0162][0163]
在同一θ0下，表中上下界对应的位移和相状态的均方差很好地说明了算法的鲁棒性，特别是相同随机噪声θ0的上下界对应的位移、相状态的均值和方差接近，说明不同随机噪声对上下界不确定性的影响接近。
[0164]
为了进一步验证本发明的鲁棒性，在θ0＝0.02的数据集中随机选取16,32,64,128个数据点将其数值变为参考值的0.8(1.2)倍，如图10所示，同时统计了30组数据误差的均值以及均方差，如下表所示，
[0165]
不同坏点数目(0.8倍参考值的坏点)下u
re
与error2的取值
[0166][0167]
不同坏点数目(1.2倍参考值的坏点)下u
re
与error2的取值
[0168]
[0169][0170]
从表格数据中误差均值可知，本发明对于坏点的位置具有很好的鲁棒性，同时也发现坏点的数目对于位移的影响很小，但其对应力的影响却随坏点数目的增大而增大。这是由于应力只在数据点组成的凸包中取值，每次线性规划仅满足平衡方程，而每个时间步的位移不仅需要满足位移协调方程和小量位移约束，也需要满足动力学平衡方程，因此能很好地保证其稳定性，这和静力学平衡方程中不包含位移项的特点不同，因此动力学算法更有利于在实际问题中对位移的控制。
[0171]
从图11-14的直方图来看，数据点的噪音和坏点对于应力的影响远大于位移，更进一步的说明数据驱动求解存在很大的波动性，这也就意味着在数据驱动的框架下，由于数据点的不确定性，追求单一的解是缺乏实际意义的。
[0172]
本发明的上下界分析。
[0173]
在实际实验中由于测量误差，实验材料的波动等等都会导致数据不确定性增加，难以确保实验数据点真的表征了材料的本构行为，更无法预测材料本构不确定性对结构动力响应带来的影响，因此在这种情况下，相比求出单一解，给出位移的上下界更有实际意义，也更利于工程师作出判断。因此，本节主要讨论了在噪声和异常值下，上下界表现的有效性。
[0174]
首先，计算了三种不同噪声θ0(0.01,0.02,0.04)时u
580
的上下界以及i(u)＝p
t
u的特殊解，相应位移响应曲线如图15所示，由图可知，在求解时间域上，上界均严格大于下界，且均能包络参考解和特殊解。且注意到随着数据噪声的增大(θ0＝0.04)，响应上下界曲线的变化趋势、峰值出现的先后和参考解有了明显不同，如图15中的(c)所示，这因为随着数据集的不确定性增加，结构响应变化趋势也会存在不确定性，并和参考解以及特殊解(i＝pu)产生差异，且随着响应时间延长，这种不确定性体现的也更加明显。这更加体现了上下界的求解不是对参考解简单的上移和下移，需要考虑数据集不确定性对动力学响应的影响，对应的上下界对工程更具有参考价值，而不是针对不确定性的数据集只求解单一的响应曲线。
[0175]
统计给定时间域内位移响应曲线上满足上界大于下界、上下界包络参考解或者特殊解的时间步个数占总时间步的比例，即有效率公式：
[0176][0177]
其中，上界大于下界时，r＝r
u》l
,n＝n
u》l
；上下界包络参考解时，r＝r
u》r》l
,n＝n
u》r》l
；上下界包络特殊解时，r＝r
u》p》l
,n＝n
u》p》l
；本算例中响应时间t＝1，响应步长δt＝0.004，总时间步n
total
＝t/δt＝250。
[0178]
为了进一步说明本发明上下界的有效性，在不同噪声θ0(0.01,0.02,0.04)时，取nc＝9,nd＝3163，分别统计30组随机噪声数据集上下界有效率情况，统计图如图16所示。在θ(0.01,0.02)，每组位移响应曲线的上界均大于下界，在位移响应曲线图中表现为虚线表示的上界均在点线表示的下界上方。在30组数据集中，每组数据集求出的上下界位移曲线包
络参考解和特殊解曲线的比例均接近100％，在位移响应曲线中表现为实线表示的参考解和点画线表示的特殊解位移上下界之间。综合来看，在这30组数据集中，上界均严格大于下界，充分说明了本文算法的有效性，上下界之间的差距也体现了考虑不确定性的必要性，同时上下界也能将所求的参考解和特殊解接近100％包含在内，说明了本文算法所求上下界的可靠性。对于θ0(0.04)，在30组响应中，位移响应上界严格大于下界的有27组(即每组的r
u》l
＝100％)，而对于其余3组，即使没有在所有时间步(250steps)中实现上界严格大于下界，但其上界大于下界的有效率分别为99.2％(248steps)，99.6％(249steps)，99.6％(249steps)；另外，提高局部凸组合数(nc＝11)后，30组数据集的上界大于下界的有效率均达到100％，如图16d所示，并且上下界包络参考解的有效率有明显提升，充分说明适当增加局部凸组合数提高了上下界的有效性。
[0179]
随后在噪声的基础上随机添加了一定数目的坏点，参考图10，计算了具有64个固定位置坏点，大小分别为0.8倍异常值和1.2倍异常值的含噪声数据集，取nd＝1001,nc＝9，θ0＝0.02，图17为不含异常值和含不同倍数异常值的数据集计算得到的位移u
580
响应曲线。
[0180]
对比图17a中无坏点的位移响应，对于含有0.8倍异常值得到的位移响应图17b，上下界和特殊解位移响应峰值有明显增大趋势，且峰值到达的时间有明显延迟现象；而对于含有1.2倍异常值得到的位移响应图17c，上下界和特殊解位移响应峰值有明显降低趋势，且峰值到达的时间有明显提前现象。这是由于具有0.8倍异常值的坏点具有更弱的材料强度，算法在选择这些坏点后导致单元在同等受力条件下变形更大，计算的位移也更大，而1.2倍异常值的坏点的结果则相反。且对于含噪声数据集中具有的坏点，更加增大了数据集的不确定性，因此，相比无坏点的位移响应曲线图17a，图17b、图17c位移响应上下界曲线的变化趋势、峰值出现的先后和参考解有了明显不同，这和前面所述的提高数据集不确定性后，结构的响应变化趋势存在不确定性一致，更加说明上下界不是对参考解简单的上移和下移，考虑数据集不确定性后求解结构响应上下界更加合理。
[0181]
另外，对于含坏点的数据集，同样测试不同坏点位置的30组数据集，上下界包络参考解以及特殊解的情况，有效率统计图如图18所示。由图可知，无论对于多少倍异常值的坏点，位移上界均严格大于位移下界，即r
u》l
＝100％，且基本能包络参考解和特殊解，充分说明本算法对于含异常值的上下界的有效性。
[0182]
基于上述结果，可以得出，考虑到数据集中不可避免的不确定性，本发明的方法针对动力学响应可以给出严格的上下界曲线，并且具有覆盖参考解的能力，而且算法得到的上下界响应曲线对数据集中的噪声和异常值也具有鲁棒性，这不仅展现了算法的鲁棒性，也具有上下界的有效性，提高了实际意义，同时与传统的ddcd框架下得到的单一解曲线相比，本方法还绕过了一些数值上的困难，求解更简单，对于稀疏的小样本数据集也具有有效性。

技术特征：
1.一种考虑不确定性的数据驱动的动力学序列线性规划方法，其特征在于，具体步骤如下：步骤s1：将动力学的二阶微分方程在时间域离散，形成多个时刻点，输入初始参数，并进行初始计算获得质量矩阵m和协调矩阵b，设置时间步δt、newmark算法参数以及初始条件；步骤s2：计算当前时刻等效右端载荷，为每个单元分配初始凸组合数据并得到线性规划求解；步骤s3：判定线性规划求解是否可行，当线性规划求解可行时，自适应更新局部凸组合点，当线性规划求解不可行时，修正后自适应更新局部凸组合点，直到更新后局部凸组合点满足收敛条件；步骤s4：通过newmark算法计算收敛后的速度和加速度，并将速度和加速度存储为当前时刻响应，返回步骤s2中进行下一时刻响应的计算；步骤s5：直到计算完时间域所有时刻点的响应，输出相应曲线结果。2.根据权利要求1所述的一种考虑不确定性的数据驱动的动力学序列线性规划方法，其特征在于：在步骤s1中，newmark算法参数包括常数参数α和常数参数β。3.根据权利要求2所述的一种考虑不确定性的数据驱动的动力学序列线性规划方法，其特征在于：在步骤s2中，当前时刻t
k
等效右端载荷计算公式如下：其中，为当前时刻等效右端载荷，p
k
为当前时刻外载荷向量，为当前时刻newmark预测值；为每个单元分配初始凸组合数据具体操作如下：从数据集中选择初始局部凸组合的应力应变数据：其中，en
c
为构成第e个单元的第n
c
个数据点的标识，其中e＝1,
…
,m，m为划分单元的个数，表示从数据集d中选出构成局部凸组合的数据点；式中初始数据编号如下：其中l
(1)
是(n
d-1)/n
c
的整数部分；根据每个单元分配初始凸组合数据计算得到线性规划求解和4.根据权利要求3所述的一种考虑不确定性的数据驱动的动力学序列线性规划方法，其特征在于，在步骤s3具体为：设l为更新次数，当线性规划求解可行时，采用如下公式获取单元应力应变状态
其中，表示为局部凸组合的系数，当线性规划求解不可行时，通过修正公式获取单元应力应变状态计算公式如下：其中，为修正系数和过渡应力应变状态参数；通过修正公式获取修正参数如下：通过修正公式获取修正参数如下：通过修正公式获取修正参数如下：其中，w
e
表示为单元体积，b
e
表示为单元协调矩阵，为修正系数；通过上述得到的单元应力应变状态得到更新后的局部凸组合的数据点更新公式如下：其中，t为n
c
/2的整数部分，l
(l+1)
＝l
(l)
+1，l
(l)
为大于等于1的整数，为距离当前状态最近的数据点的编号；收敛判断公式如下：其中δ为收敛设定参数。

技术总结
本发明公开了一种考虑不确定性的数据驱动的动力学序列线性规划方法，具体步骤如下：输入初始参数并进行初始计算，计算当前时刻等效右端载荷，分配初始凸组合数据并得到线性规划求解，当线性规划求解可行时，自适应更新局部凸组合点，当线性规划求解不可行时，修正后自适应更新局部凸组合点，直到更新后局部凸组合点收敛；通过Newmark算法计算收敛后的速度和加速度，并将速度和加速度存储为当前时刻响应，进行下一时刻响应的计算；直到计算完时间域所有时刻点的响应，输出相应曲线结果。采用上述一种考虑不确定性的数据驱动的动力学序列线性规划方法，明确考虑了数据集中不可避免的不确定性，得到响应曲线更具有可信度，更有利于工程师的判断。利于工程师的判断。利于工程师的判断。


技术研发人员：郭旭 杜宗亮 孙培勇 黄孟成 刘畅 张维声
受保护的技术使用者：大连理工大学
技术研发日：2023.03.10
技术公布日：2023/10/8