EringiTech’s blog

I am a mechanical engineer

【ワイブル分布】熱真空試験のサイクル数

1.概要

人工衛星は軌道上で動作するため，高温と低温の環境を交互にさらされる．

軌道上に打ち上げてからテストしても直すことが出来ないため，地上で熱真空環境を作り出しテストを行う． この試験を行うとき，高温低温を1サイクルとして，試験を何サイクル行うべきか迷う．

そんなときに参考文献[1]が役に立った．
文書中ではワイブル分布を用いており，最終的にスクリーニング率を99.9%にするためには8サイクル以上行う方が良いといったことが記載されている．

このブログでは，ワイブル分布の使い方に注目して導出を追ってみる． 今回は，なぜワイブル分布を使うのかというところまで記述している．

2.目的

  • ワイブル分布について復習する
  • ワイブル分布の応用例として参考文献[1]の導出を追う

3.ワイブル分布とは

3.1.式

尺度パラメータ \lambda と形状パラメータ c の2つのパラメータと確率変数 x を用いて確率密度変数(PDF)は

 
p(x)=
\begin{cases}
\frac{c}{\lambda} (\frac{x}{\lambda})^{c - 1} e^{-(\frac{x}{\lambda})^{c}} & x \ge 0 \\
0 & x \lt 0
\end{cases}

と表せる．

3.2.グラフ

とりあえず式だけみても意味がわからないので，図示をしてみる． 自分で書くのも面倒くさいのでChatGPTに頼み，各種パラメータにおけるプロットを1つの図に示した．

各種パラメータにおけるワイブル分布

グラフから

  • 形状パラメータ c が小さいときは単調減少

がわかる（まだ他はよく分かってないので割愛）

3.3.コード

念のため，コードを示しておく．

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import weibull_min

# ワイブル分布のパラメータを設定
c_values = [1, 2, 5]  # 形状パラメータ
lambda_values = [1, 2, 3]  # 尺度パラメータ

# グラフの設定
x = np.linspace(0, 5, 1000)  # x軸の範囲を設定

# パラメータを変えてワイブル分布をプロット
plt.figure(figsize=(10, 6))
for c in c_values:
    for lambd in lambda_values:
        pdf = weibull_min.pdf(x, c, scale=lambd)  # ワイブル分布の確率密度関数を計算
        label = f'Weibull (c={c}, λ={lambd})'
        plt.plot(x, pdf, label=label)

plt.title('Weibull Distribution with Different Parameters')
plt.xlabel('x')
plt.ylabel('Probability Density')
plt.legend()
plt.grid(True)
plt.show()

3.4.式の意味

式と図を書いても何となくしか分かりませんが，とりあえず以下のように考えればよさそうです．

ある時間 x で故障する確率 f(x)

4.信頼性モデル

4.1.ワイブル分布(再掲）

参考文献1の内容に入っていく．
しかし，ワイブル分布のパラメータ変数が違ったのでもう一度ワイブル分布を示す．

尺度パラメータ \eta と形状パラメータ mを用いて確率密度関数 f(x) は

 
f(x)=\frac{m}{\eta} (\frac{x}{\eta})^{m - 1} e^{-(\frac{x}{\eta})^{m}}

と表せる．また，累積分布関数(つまり確率密度関数の積分） F(x)は

 
F(x)=\int^{x}_{-\infty} {f(x)dx}=\int^{x}_{0} {f(x)dx} = (-{e^{-(\frac{x}{\eta})^{m}}} )_0^x = -{e^{-(\frac{x}{\eta})^{m}}} + 1

と求めることが出来ます． なお， x\lt 0 で f(x)=0 であることと，以下の微分を使って累積分布関数 F(x) を求めました．

 
\frac{d}{dx} {e^{-(\frac{x}{\eta})^{m}}} = \frac{d}{dx} ({-(\frac{x}{\eta})^{m}}) \times  {e^{-(\frac{x}{\eta})^{m}}} =  \frac{m}{\eta} \times ({-(\frac{x}{\eta})^{m-1}}) \times {e^{-(\frac{x}{\eta})^{m}}}

4.2.故障しない確率

ある時間 x まで故障しない確率は， x以上で故障する確率を足し合わせたものと考えられそうです．
つまり， x\to\infty で確率密度関数を足し合わせたものと考えられそうです．

参考文献[1]の式(4.2.2-2)がこの考えから求めることが出来ます．
信頼度 R(x) は，

 
R(x)=\int^{\infty}_{x} {f(x)dx}= (-{e^{-(\frac{x}{\eta})^{m}}} )^{\infty}_x =0 - (-{e^{-(\frac{x}{\eta})^{m}}} ) = e^{-(\frac{x}{\eta})^{m}}

と求められます．

この考え方だと， -\infty\to x を計算した F(x) というのは，全体確率から信頼度を引いたものと言えそうです．
なので， F(x) を不信頼度 と呼ぶそうです．

4.3.故障率

信頼度 R(x) は， x=0 で R(0)=1 であり，時間とともに徐々に減少していく．
逆に不信頼度 F(x) は， x=0 で F(0)=0 である．

ここで，単位時間の故障率(瞬間故障率） \lambda というものを

 
\lambda (x) =\frac{f(x)}{R(x)} = \frac{\frac{m}{\eta} (\frac{x}{\eta})^{m - 1} e^{-(\frac{x}{\eta})^{m}}}{e^{-(\frac{x}{\eta})^{m}}}=\frac{m}{\eta} (\frac{x}{\eta})^{m - 1}

とする．式からわかるようように，指数関数についていた係数部分が瞬間故障率となっている．

瞬間故障率から以下のことが分かる

  1. x \lt 1 のとき， xのべき乗は負の値となるため，単調減少する．これは，瞬間故障率が時間とともに減少していくことを示す．『初期故障』に相当する．
  2. x = 1 のとき，[tex: x0=1] となり，時間に関わらず一定の値となる．『偶発故障』に相当する．
  3. x \gt 1 のとき， x のべき乗は正の値となるため，単調増加する．つまり，故障率が時間とともに増加していく．『摩耗故障』に相当するらしい．

4.4.なぜワイブル分布を信頼性モデルとして使うのか

人工衛星の試験を行うとき，体感的によくあるミスは単純な人の手によるものである．
衛星を打ち上げる前に，受入試験(AT)などを行う．
これはロケット側に影響がないかや本当に起動で動作できるのかを確認する．（例えば，ロケットで急に爆発でもしたら大変である）

このロケットに搭載してよいかなどを試験するときに，どんなエラーを確認したらよいであろうか．
* 上述したようなワークマンシップエラーや材料や部品の不良（欠陥）→いわゆる初期不良 を確認したくなる．

また，参考文献[1]の図4.2.2-4のように軌道上での運用日数と不具合の累積検出率がワイブル分布の累積分布関数に類似している．

これらのことからワイブル分布を使うようである．（自動車などでも同じだと思うが，どうなのかまた調べてみる）

5.まとめ

  • ワイブル分布は，パラメータによって様々な形状を取ることが分かった．
  • 時間とともに故障率がどのように変化するかを考えた
  • 過去の実績とヒューマンエラーなどの初期不良を検出するという目的がワイブル分布の特性とマッチしている

6.次回

  • Test Effectivenessの概念を学んでいく

参考文献

[1] JAXA, "環境試験信頼性要求ハンドブック”，https://sma.jaxa.jp/TechDoc/Docs/JAXA-JERG-2-130-HB006_N1.pdf [2] https://engineeeer.com/hatena-markdown-python/ [3]LATEXチートシート - 数式記号の読み方・表し方 - SHOYAN BLOG [4] はてなブログでTex記法を使って行列を書く時の注意点 - 医療職からデータサイエンティストへ

おわりに

ご覧いただきありがとうございます．

もしよければ，X (Twitter)とnoteのフォローをお願いします．

twitter.com

  • note

はてなブログの内容をまとめております．

note.com