蒙特卡洛模拟极限状态函数代码报错求助

问题详情 #

• 耗时两天编写蒙特卡洛模拟极限状态函数的VBA代码，始终出现报错
• 手动计算得出失效概率为2.3e-4、安全系数Beta值为3.5
• 代码中计算Beta的语句Beta = -Application.WorksheetFunction.Norm_S_Inv(Pf)无法处理小于0的数值，特此求助

代码示例 #

Private Sub CommandButton1_Click()
    ' worked example 6.2b
    Dim NumFails As Long
    Dim Ntrials As Long
    Dim Pi_2 As Double
    Dim Pf As Double
    Dim Beta As Double
    Dim u1 As Double, u2 As Double, u3 As Double
    Dim fb As Double, h As Double, b As Double, w As Double, q As Double
    Dim S As Double, Gx As Double
    Dim i As Long
    
    ' Initialize variables
    NumFails = 0
    Ntrials = 10000 ' Number of trials
    Pi_2 = 2 * 3.14159265358979
    Randomize
    
    ' Statistical properties
    Dim mean_fb As Double, cov_fb As Double, std_fb As Double
    Dim mean_h As Double, cov_h As Double, std_h As Double
    Dim mean_b As Double, cov_b As Double, std_b As Double
    Dim mean_w As Double, cov_w As Double, std_w As Double
    Dim mean_q As Double, cov_q As Double, std_q As Double
    Dim L As Double
    
    mean_fb = 30 ' N/mm²
    cov_fb = 0.15
    std_fb = mean_fb * cov_fb
    
    mean_h = 757 ' mm
    cov_h = 0.04
    std_h = mean_h * cov_h
    
    mean_b = 250 ' mm
    cov_b = 0.03
    std_b = mean_b * cov_b
    
    mean_w = 20 ' kN/m
    cov_w = 0.1
    std_w = mean_w * cov_w
    
    mean_q = 15 ' kN/m
    cov_q = 0.18
    std_q = mean_q * cov_q
    
    L = 8000 ' Beam span in mm (8m)
    
    ' Monte Carlo Simulation
    For i = 1 To Ntrials
        ' Generate random variables
        u1 = Rnd
        u2 = Rnd
        u3 = Rnd

        ' Generate lognormal random variable
        Dim mu_fb As Double, sigma_fb As Double
        mu_fb = Log(mean_fb / Sqr(1 + (std_fb / mean_fb) ^ 2))
        sigma_fb = Sqr(Log(1 + (std_fb / mean_fb) ^ 2))
        h = mean_h + std_h * Sqr(-2 * Log(Rnd)) * Cos(Pi_2 * Rnd)
        b = mean_b + std_b * Sqr(-2 * Log(Rnd)) * Cos(Pi_2 * Rnd)
        w = mean_w + std_w * Sqr(-2 * Log(Rnd)) * Cos(Pi_2 * Rnd)
        q = mean_q - (std_q / 0.5772) * Log(-Log(u3))
        
        ' Section modulus
        S = (b * (h ^ 2)) / 6
        
        ' Evaluate limit state function
        Gx = fb * S - ((w + q) * (L ^ 2)) / 8
        
        ' Count failures
        If Gx < 0 Then
            NumFails = NumFails + 1
        End If
    Next i
    
    ' Probability of failure and safety index
    Pf = NumFails / Ntrials
    Beta = -Application.WorksheetFunction.Norm_S_Inv(Pf)
    
    ' Output results
    Range("C8").Value = Ntrials
    Range("C10").Value = NumFails
    Range("C12").Value = Pf
    Range("C14").Value = Sqr((1 - Pf) / (Ntrials - 1) / Pf) ' COV of Pf
    Range("C16").Value = Beta
End Sub

核心问题 #

代码中Beta = -Application.WorksheetFunction.Norm_S_Inv(Pf)语句无法处理小于0的数值，导致报错，而手动计算的Beta值应为3.5。

内容的提问来源于stack exchange，提问作者martin conroy