AOKI's copy&paste archive

高専から駅弁大学から東工大を経て大企業へ 浅く広い趣味とかキャリアの日記を

Monte Carlo Simulation to Calculate Reliability Index

学部の卒論に用いたメインモジュール

# -*- coding: utf-8 -*-
"""
Created on Mon Nov 19 15:15:20 2018

@author: Aoki
"""
#span center bending tention
#import
import numpy as np
import matplotlib.pyplot as plt
import runpy
import ini
#%%
#initialize
i=0
rlist=[]
slist=[]
zlist=[]
num=100000#simulation time
#%%
while i<num:
    if i%5000==0:#advance check
        print('i',i)
    #run all times
    #cf:https://teratail.com/questions/165243
    stress_mod = runpy.run_module('stress')
    #calcurate z(=R-S)
    #cf:manual(shihosho)
    #dim:N/mm2
    r=0.9*0.85*stress_mod["yk"]
    s=1.05*(stress_mod["sd"]+stress_mod["vd"]+ini.others)+1.25*stress_mod["l"]
    z=r-s
    rlist.append(r)
    slist.append(s)
    zlist.append(z/ini.strength_s)#to no dim
    i+=1
    #end while loop
#%%
if __name__ == '__main__':
    #graph drow
    #z hist
    plt.hist(zlist,bins=100,density=True)
    plt.title("",fontsize=20)
    plt.xlabel('${Z/{\sigma}_{yk}}$',fontsize=16)
    plt.ylabel('Probability density',fontsize=16)
    plt.xlim(-0.2,0.6)
    plt.show()
    #R&S hist
    plt.hist(rlist,bins=100,density=True,alpha=0.5,label='Resistance')
    plt.hist(slist,bins=100,density=True,alpha=0.5,label='Load effect')
    plt.xlabel('stress($N/{mm}^2$)',fontsize=16)
    plt.ylabel('Probability Density',fontsize=16)
    plt.legend()
    plt.xlim(300,550)
    plt.show()
    #
    print('r std',np.std(rlist))
    print('s std',np.std(slist))
#%%
#calculate & output reliability index
print('mu ',"{:.4f}".format(np.mean(zlist)))
print('std',"{:.4f}".format(np.std(zlist)))
beta=np.mean(zlist)/np.std(zlist)
print('Reliability Index beta',"{:.2f}".format(beta))

今思うとWhile文じゃなくてFor文でよかったじゃんとか色々思うところはあるが
これは複数モジュールで構成されていて初期値というかパラメータをini.pyから読み込んでいる
また先日のモジュールも直接は出てこないが孫請け的にShape.pyとして参照されている
pytho.hatenablog.com
これを受けて応力をstress.pyで算出しrunpyでシミュレーションごとに参照するコードの構成である
このrunpyが味噌で当初はバラツキが与えられず困ったのだが質問サイトの答えから参照の仕方を工夫することで乗り越えた
1つのモジュールに詰め込むのもありだったが役割ごとに分割した方がやはり可読性が高いだろうということで敬遠した
特に研究の引継ぎもあってエゴでは書けなかったし
ただ所詮工学系が1年間の付け焼刃で書いたコードなので高は知れているかと
あとはギリシャ文字や添え字がやや面倒だったがソースを見ればその通りという感じ