《1 引言》
1 引言
在工程结构的可靠性分析中, 一次二阶矩法 (包括中心点法、JC法、映射变换法和实用分析法) 以及二次二阶矩法得到了广泛应用
[1 ,2 ] 。一次 (二次) 二阶矩法主要利用随机变量的均值和方差信息以及分布概型来计算结构的失效概率, 属于低阶矩法。使用低阶矩法的前提是:采用的随机变量的分布概型是正确的, 随机变量的有关统计参数是准确的。然而, 在实际工程中, 获得的样本数量往往很小, 在一定的置信度下, 2个或者几个不同的分布概型都有可能被接受, 随机变量的统计参数可能与真实值差别较大。此外, 目前采用的随机变量分布概型几乎都是理想的数学模型, 现实中的问题由于受多种因素影响可能并不严格符合理想的数学模型, 而是服从某种未知的分布或几种分布组合而成的分布
[3 ,4 ] 。为了全面反映和利用随机变量的统计信息, 利用随机变量的高阶矩对工程结构进行可靠性分析是一条新思路
[5 ] 。四阶矩分析法用到了随机变量的前四阶矩信息, 对结构功能函数非线性程度不高的结构可靠性分析问题, 比较适用
[6 ] 。以{1, x , \:, x n }为基的多项式逼近法在多项式的阶数不是很高的情况下, 可以较好地逼近功能函数的概率密度函数
[4 ] 。当阶数较高时容易产生龙格现象, 计算结构失效概率时误差较大。因此, 在实际工程中, 如何利用已有的统计信息, 对结构的失效做出符合实际情况的概率测度, 在结构可靠性研究中很有必要。作者基于数值逼近原理, 充分考虑随机变量的高阶矩信息, 利用切比雪夫正交多项式来逼近计算功能函数的概率密度函数, 进而可以利用数值积分计算结构的可靠度, 并给出了计算实例。
《2 随机变量概率密度函数的正交多项式逼近》
2 随机变量概率密度函数的正交多项式逼近
根据工程结构可靠度的定义
[1 ,2 ] , 结构功能函数可以用Z = (X 1 , X 2 , \:, X n ) 表示, X 1 , X 2 , \:, X n 表示结构的基本随机变量。假定结构的抗力随机变量为R 和荷载效应随机变量为S 是两个独立的连续随机变量, 其相应的概率密度函数为f R (r ) , f S (s ) , 则结构的失效概率p f 为:
p f = P ( Z < 0 ) = ∫ ∫ r < s f R ( r ) f S ( s ) d r d s ( 1 ) p f = Ρ ( Ζ < 0 ) = ∫ ∫ r < s f R ( r ) f S ( s ) d r d s ( 1 )
在实际问题中, 有时功能函数是非线性的, 比较复杂, 对f R (r ) 或f S (s ) 积分难以得到显式表达式, 因此直接通过数值积分计算结构的失效概率, 在实际工程中难以实现。正交多项式逼近法是采用正交函数族做基, 利用功能函数的高阶矩信息, 通过计算得到功能函数概率密度函数的近似表达式, 进而可以利用数值积分计算结构的失效概率。其具体推导过程如下:
设f (x ) , g (x ) 是闭区间[a , b ]上关于x 的连续函数, η (x ) 是[a , b ]上的权函数, 积分
( f ( x ) , g ( x ) ) = ∫ b a η ( x ) f ( x ) g ( x ) d x ( 2 ) ( f ( x ) , g ( x ) ) = ∫ a b η ( x ) f ( x ) g ( x ) d x ( 2 )
称为函数f (x ) 与g (x ) 在[a , b ]上的内积。满足内积定义的函数空间C [a , b ]称为内积空间。并且在C [a , b ]上首项系数a n ≠0的n 次多项式满足
∫ b a ω ( x ) φ i ( x ) φ j ( x ) d x = { 0 , i ≠ j A k > 0 , i = j ( i , j = 0 , 1 , \ : ) ( 3 ) ∫ a b ω ( x ) φ i ( x ) φ j ( x ) d x = { 0 , i ≠ j A k > 0 , i = j ( i , j = 0 , 1 , \ : ) ( 3 )
称{φ n (x ) }是C [a , b ]上带权ω (x ) 的正交多项式
[7 ] 。设f (x ) ∈C [a , b ], φ 0 (x ) , φ 1 (x ) , …, φ n (x ) 为空间C [a , b ]中的一个线性无关系, a 0 , a 1 , \:, a n 为任意实数, 则S ( x ) = ∑ i = 0 n ω ( x ) a i φ i ( x ) S ( x ) = ∑ i = 0 n ω ( x ) a i φ i ( x ) 是C [a , b ]中是的一个子集。根据函数最佳平方逼近原理
[7 ] , 存在最佳平方逼近函数S * (x ) , 使得
∥ f ( x ) − S ∗ ∥ 2 2 = inf S ∈ ∥ f ( x ) − S ∗ ∥ 2 2 = inf S ∈ C ∫ b a η ( x ) [ S ( x ) − f ( x ) ] 2 d x ( 4 ) ∥ f ( x ) - S * ∥ 2 2 = inf S ∈ ∥ f ( x ) - S * ∥ 2 2 = inf S ∈ C ∫ a b η ( x ) [ S ( x ) - f ( x ) ] 2 d x ( 4 )
即求I = ∫ b a η ( x ) [ ∑ i = 0 n ω ( x ) a i φ i ( x ) − f ( x ) ] 2 d x Ι = ∫ a b η ( x ) [ ∑ i = 0 n ω ( x ) a i φ i ( x ) - f ( x ) ] 2 d x 的最小值。
由于{φ i (x ) }为C [a , b ]上的一个无关系, 利用多元函数求极值的方法, 可得方程:
∂ I ∂ α k = 2 ∫ b a η ( x ) [ ∑ i = 0 n ω ( x ) a i φ i ( x ) − f ( x ) ] ⋅ φ k ( x ) d x = 0 ( k = 0 , 1 , 2 , \ : , ) ∂ Ι ∂ α k = 2 ∫ a b η ( x ) [ ∑ i = 0 n ω ( x ) a i φ i ( x ) - f ( x ) ] ⋅ φ k ( x ) d x = 0 ( k = 0 , 1 , 2 , \ : , )
∫ b a [ ∑ i = 0 n ω ( x ) a i φ i ( x ) ] φ k ( x ) d x = ∫ b a f ( x ) φ k ( x ) d x ( k = 0 , 1 , 2 , \ : , ) ( 5 ) ∫ a b [ ∑ i = 0 n ω ( x ) a i φ i ( x ) ] φ k ( x ) d x = ∫ a b f ( x ) φ k ( x ) d x ( k = 0 , 1 , 2 , \ : , ) ( 5 )
ω ( φ 0 , φ 0 ) a 0 + ω ( φ 0 , φ 1 ) a 1 + \ : + ω ( φ 0 , φ n ) a n = ( φ 0 , f ) ω ( φ 1 , φ 0 ) a 0 + ω ( φ 1 , φ 1 ) a 1 + \ : + ω ( φ 1 , φ n ) a n = ( φ 1 , f ) ⋯ ⋯ ω ( φ n , φ 0 ) a 0 + ω ( φ n , φ 1 ) a 1 + \ : + ω ( φ n , φ n ) a n = ( φ 0 , f ) ⎫ ⎭ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ( 6 ) ω ( φ 0 , φ 0 ) a 0 + ω ( φ 0 , φ 1 ) a 1 + \ : + ω ( φ 0 , φ n ) a n = ( φ 0 , f ) ω ( φ 1 , φ 0 ) a 0 + ω ( φ 1 , φ 1 ) a 1 + \ : + ω ( φ 1 , φ n ) a n = ( φ 1 , f ) ⋯ ⋯ ω ( φ n , φ 0 ) a 0 + ω ( φ n , φ 1 ) a 1 + \ : + ω ( φ n , φ n ) a n = ( φ 0 , f ) } ( 6 )
这是关于a 0 , a 1 , \:, a n 的线性方程组。因为φ 0 (x ) , φ 1 (x ) , \:, φ n (x ) 线性无关, 故系数行列式G (φ 0 , φ 1 , \:, φ n ) ≠0, 方程组 (6) 有唯一解。
根据函数族{φ n (x ) }的正交性, 方程组 (6) 左边系数矩阵除对角线上元素外, 其余都为0, 所以得到方程组 (6) 的解, 即:
a k = ( φ k , f ) ω ( φ k , φ k ) ( 7 ) a k = ( φ k , f ) ω ( φ k , φ k ) ( 7 )
S k = ∑ k = 0 n ( φ k , f ) ω ( φ k , φ k ) φ k ( x ) = ∑ k = 0 n a k φ k ( x ) ( 8 ) S k = ∑ k = 0 n ( φ k , f ) ω ( φ k , φ k ) φ k ( x ) = ∑ k = 0 n a k φ k ( x ) ( 8 )
考虑f (x ) ∈C [-1, 1], 利用切比雪夫正交多项式{T k (x ) }逼近f (x ) , 其中权函数ω ( x ) = 1 1 − x 2 √ ω ( x ) = 1 1 - x 2 , 并且满足如下关系:
T 0 ( x ) = 1 , T 1 ( x ) = x , T n + 1 ( x ) = 2 x T n ( x ) − T n − 1 ( x ) . ∫ − 1 1 T m ( x ) T n ( x ) d x 1 − x 2 − − − − − √ = ⎧ ⎩ ⎨ ⎪ ⎪ 0 , m ≠ n π 2 , m = n ≠ 0 π , m = n = 0 ( 9 ) Τ 0 ( x ) = 1 , Τ 1 ( x ) = x , Τ n + 1 ( x ) = 2 x Τ n ( x ) - Τ n - 1 ( x ) . ∫ - 1 1 Τ m ( x ) Τ n ( x ) d x 1 - x 2 = { 0 , m ≠ n π 2 , m = n ≠ 0 π , m = n = 0 ( 9 )
( φ k , f ) = ( T k , f ) = ∫ 1 − 1 T k ( x ) f ( x ) d x ( 1 0 ) ( φ k , f ) = ( Τ k , f ) = ∫ - 1 1 Τ k ( x ) f ( x ) d x ( 1 0 )
其中T k (x ) 为一般多项式, 设其形式为B k 0 +B k 1x +…+B kk x k , 代入式 (10) , 有
∫ 1 − 1 T k ( x ) f ( x ) d x = ∫ 1 − 1 ( B k 0 + B k 1 x + \ : + B k k x k ) f ( x ) d x = B k 0 ∫ 1 − 1 f ( x ) d x + B k 1 ∫ 1 − 1 x f ( x ) d x + ⋯ + B k k ∫ 1 − 1 x k f ( x ) d x ∫ - 1 1 Τ k ( x ) f ( x ) d x = ∫ - 1 1 ( B k 0 + B k 1 x + \ : + B k k x k ) f ( x ) d x = B k 0 ∫ - 1 1 f ( x ) d x + B k 1 ∫ - 1 1 x f ( x ) d x + ⋯ + B k k ∫ - 1 1 x k f ( x ) d x
μ = ∫ 1 − 1 x f ( x ) d x μ k = ∫ 1 − 1 x k f ( x ) d x ⎫ ⎭ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ( 1 1 ) μ = ∫ - 1 1 x f ( x ) d x μ k = ∫ - 1 1 x k f ( x ) d x } ( 1 1 )
( φ k , f ) = ( T k , f ) = ∫ 1 − 1 T k ( x ) f ( x ) d x = B k 0 μ 0 + B k 1 μ 1 + \ : + B k k μ k = ∑ i = 0 k B k i μ i ( 1 2 ) ( φ k , f ) = ( Τ k , f ) = ∫ - 1 1 Τ k ( x ) f ( x ) d x = B k 0 μ 0 + B k 1 μ 1 + \ : + B k k μ k = ∑ i = 0 k B k i μ i ( 1 2 )
ω ( φ k , φ k ) = ω ( T k , T k ) = ∫ 1 − 1 T k ( x ) T k ( x ) 1 − x 2 − − − − − √ d x = { π , k = 0 π 2 , k ≠ 0 ( 1 3 ) ω ( φ k , φ k ) = ω ( Τ k , Τ k ) = ∫ - 1 1 Τ k ( x ) Τ k ( x ) 1 - x 2 d x = { π , k = 0 π 2 , k ≠ 0 ( 1 3 )
α 0 = ( T 0 , f ) ω ( T 0 , T 0 ) = B 0 0 μ 0 π , a 1 = ( T 1 , f ) ω ( T 1 , T 1 ) = B 1 0 μ 0 + B 1 1 μ 1 π / 2 , \ : , a k = ( T k , f ) ω ( T k , T k ) = ∑ i = 0 k B k i μ i π / 2 , \ : ( 1 4 ) α 0 = ( Τ 0 , f ) ω ( Τ 0 , Τ 0 ) = B 0 0 μ 0 π , a 1 = ( Τ 1 , f ) ω ( Τ 1 , Τ 1 ) = B 1 0 μ 0 + B 1 1 μ 1 π / 2 , \ : , a k = ( Τ k , f ) ω ( Τ k , Τ k ) = ∑ i = 0 k B k i μ i π / 2 , \ : ( 1 4 )
得到概率密度函数f (x ) 的多项式逼近表达式f n (x ) , 即f ( x ) ≈ f n ( x ) = ∑ k = 0 n ω ( x ) a k T k ( x ) f ( x ) ≈ f n ( x ) = ∑ k = 0 n ω ( x ) a k Τ k ( x ) 。
一般情况下, f (x ) ∈C [a , b ]。因为切比雪夫多项式定义在区间[-1, 1]上, 所以求解在[a , b ]上的逼近多项式, 需做区间变换
x = b − a 2 t + b + a 2 ( − 1 ≤ t ≤ 1 ) x = b - a 2 t + b + a 2 ( - 1 ≤ t ≤ 1 )
x —t 之间为线性关系, 当x 在区间[a , b ]上变化时, 对应的t 在区间[-1, 1]上变化, 则
T 0 ( t ) = T 0 ( 2 x − ( b + a ) b − a ) = 1 ‚ T 1 ( t ) = t = 2 x − ( b + a ) b − a , T 2 ( t ) = 3 t 2 − 1 2 = 3 ( 2 x − ( b + a ) b − a ) 2 − 1 2 ‚ ⋯ Τ 0 ( t ) = Τ 0 ( 2 x - ( b + a ) b - a ) = 1 ‚ Τ 1 ( t ) = t = 2 x - ( b + a ) b - a , Τ 2 ( t ) = 3 t 2 - 1 2 = 3 ( 2 x - ( b + a ) b - a ) 2 - 1 2 ‚ ⋯
最后, 得到f (x ) 的多项式逼近表达式f n (x ) , 即
f ( x ) ≈ f n ( x ) = ∑ k = 0 n ω ( t ) a k T k ( t ) = ∑ k = 0 n a k ω ( 2 x − ( b + a ) b − a ) T k ( 2 x − ( b + a ) b − a ) ( 1 5 ) f ( x ) ≈ f n ( x ) = ∑ k = 0 n ω ( t ) a k Τ k ( t ) = ∑ k = 0 n a k ω ( 2 x - ( b + a ) b - a ) Τ k ( 2 x - ( b + a ) b - a ) ( 1 5 )
在结构可靠性分析中, 以切比雪夫正交多项式做基, 利用功能函数的高阶矩 (即f (x ) 的前k 阶原点矩) , 求解方程组 (6) , 得到各项系数a k , 即可求得功能函数概率密度函数的逼近表达式, 对结构的可靠性进行分析。
《3 结构可靠性分析的高阶矩法》
3 结构可靠性分析的高阶矩法
基于数值逼近的结构可靠性高阶矩法是以切比雪夫正交多项式{T k (x ) }做基, 利用功能函数f (x ) 的高阶原点矩, 求解线性方程组, 得到功能函数概率密度函数的逼近表达式f n (x ) , 然后估算结构的失效概率。具体计算步骤如下:
1) 分析工程结构的实际情况, 建立功能函数Z =g (X ) = (X 1 , X 2 , …, X n ) 。
2) 根据结构随机变量的统计特征或分布概型, 对g (X ) 进行抽样。
3) 考虑正交多项式性质, 区间做适当变换, 计算功能函数g (X ) 的各阶样本矩。
4) 以切比雪夫正交多项式{T k (x ) }做基, 利用得到的各阶样本矩, 求解功能函数概率密度函数f (x ) 的近似多项式f n (x ) 。
5) 计算结构的失效概率, 对工程结构进行可靠性分析。
在计算失效概率时
[1 ,2 ] , 引入标准化随机变量z , 即 (μ z =0, σ z =1) :
把f (x ) 的近似多项式f n (x ) 代入式 (1) , 得
p f = P ( Z < 0 ) = P ( z < − μ z / σ z ) = ∫ − μ z σ z − ∞ f n ( x ) ( 1 5 ) p f = Ρ ( Ζ < 0 ) = Ρ ( z < - μ z / σ z ) = ∫ - ∞ - μ z σ z f n ( x ) ( 1 5 )
在现实情况中, 随机变量不可能无限制地趋近于-∞, 因此, 对上式积分时, 积分下界可根据工程结构的实际情况估计一个界值, 使误差在可接受的范围之内。
《4 结构可靠性高阶矩法的检验》
4 结构可靠性高阶矩法的检验
《4.1 经典分布函数的数值检验》
4.1 经典分布函数的数值检验
以工程结构随机变量最常服从的正态分布和对数正态分布为例, 检验高阶矩法拟合各种经典分布曲线的优良性。首先根据经典概率密度函数的解析表达式, 画出其理论分布曲线, 然后用积分的方法计算出理论分布的各阶原点矩, 采用切比雪夫正交多项式为基, 做拟合曲线并进行比较, 观察后者拟合性能的优劣。
图1, 图2分别是用切比雪夫正交多项式为基对正态分布、对数正态分布 (具体表达式见表1) 的拟合曲线。从图中可以看出, 拟合曲线与两种经典分布已经基本重合。通过数值检验说明, 只要已知随机变量的前n 阶原点矩, 采用切比雪夫正交多项式为基生成的概率密度函数具有良好的拟合性能。
《图1》
图1 正态密度函数的数值检验
Fig.1 Numerical test of normal density function
《图2》
图2 对数正态密度函数的数值检验
Fig.2 Numerical test of log-normal density function
表1 两种经典分布曲线及其表达式
Table 1 Two classical distributions and their expressions
《表1》
名 称
正态分布
对数正态分布
表达式
y = 1 2 π − − √ ⋅ 0 . 3 e − x 2 2 ⋅ ( 0 . 3 ) 2 y = 1 2 π ⋅ 0 . 3 e - x 2 2 ⋅ ( 0 . 3 ) 2
y = 1 0 . 1 2 π − − √ x e { − 1 2 [ l n x − 1 0 . 1 ] 2 } y = 1 0 . 1 2 π x e { - 1 2 [ l n x - 1 0 . 1 ] 2 }
置信区间
[-1, 1]
[2, 4]
置信度β /%
99.91
99.89
《4.2 工程结构失效概率的计算》
4.2 工程结构失效概率的计算
某结构构件正截面强度的功能函数为Z = (R , S ) =R -S , 其中抗力R 和荷载效应S 均服从正态分布, 平均值μ R =100 kN·m, μ S =50 kN·m, 变异系数δ R =0.12, δ S =0.15。试求结构构件的失效概率。
因为抗力R 和荷载效应S 均服从正态分布, 所以σ R =μ R ·δ R =12, σ S =μ S ·δ S =7.5。功能函数Z 也服从正态分布, 并且μ z =50 kN·m, σ 2 z z 2 =200.25。-μ z /σ z =-3.533 3, 因此在区间[-5, -3]上进行概率密度函数的逼近, 同时将功能函数随机变量Z 标准化为z 。根据正交多项式逼近法, 利用样本数据的前n 阶矩, 计算得到在区间[-5, -3]上标准化功能函数z 的概率密度函数各阶近似表达式f n (t ) 。近似表达式的各项系数见表2, 把t =x +4代入, 即可求得f n (x ) 的具体表达式。
表2 切比雪夫多项式生成函数各项系数及计算结果
Table 2 The coefficients of generated probability density functions and result of p f
《表2》
各项系数
f 7 (x )
f 8 (x )
f 9 (x )
f 10 (x )
a 0
0.000 429 6
0.000 429 6
0.000 429 6
0.000 429 6
a 1
0.000 616 3
0.000 616 3
0.000 616 3
0.000 616 3
a 2
0.000 144 9
0.000 144 9
0.000 144 9
0.000 144 9
a 3
-0.000 167 6
-0.000 167 6
-0.000 167 6
-0.000 167 6
a 4
-0.000 236 3
-0.000 236 3
-0.000 236 3
-0.000 236 3
a 5
-0.000 184 4
-0.000 184 4
-0.000 184 4
-0.000 184 4
a 6
-0.000 120 5
-0.000 120 5
-0.000 120 5
-0.000 120 5
a 7
-7.812×10-5
-7.812×10-5
-7.812×10-5
-7.812×10-5
a 8
-5.428×10-5
-5.428×10-5
-5.428×10-5
a 9
-4.044×10-5
-4.044×10-5
a 10
-3.166×10-5
p f
2.032×10-5
2.053×10-5
2.058×10-5
2.058×10-5
根据结构可靠性分析失效概率的定义, 利用样本数据的前n 阶矩计算得到的f n (x ) , 分别进行数值积分求解失效概率p f , 计算结果见表2第13行。可以看出, 当n ≥7时, 计算结果已经基本趋于稳定。利用中心点法、验算点法以及蒙特卡罗模拟
[1 ] 的计算结果分别为2.078×10-4 , 2.05×10-4 和2.05×10-4 , 从中可以看出该方法的正确性和实用性。
《5 结论》
5 结论
1) 提出了工程结构可靠性分析的高阶矩法。该方法基于数值逼近原理, 以切比雪夫正交多项式做基, 利用功能函数的高阶矩信息, 通过求解线性方程组, 得到功能函数的逼近概率密度函数表达式f n (x ) , 然后利用工程结构可靠性分析的一般式计算结构的失效概率, 对工程结构的可靠性进行分析。
2) 高阶矩分析法能充分利用结构随机变量的统计信息, 避开了复杂的数值寻优计算, 很大程度地提高了计算效率, 并且可以按任意指定的高阶矩信息, 达到满意的计算精度。通过计算, 可以依次得到多项式各项系数, 计算难度和复杂性不随矩阶数的增加而增加, 数值计算结果稳定。根据结构构件失效概率的计算, 表明可靠性分析的高阶矩数值逼近法是正确而且实用的。