シックスシグマ手法の計画によるエンジン冷却ファンの改善

この例では次を使用します。

この例では、定義、測定、分析、改善、管理 (DMAIC) を使用するシックスシグマ手法の計画によりエンジン冷却ファンの性能を改善する方法を示します。初期状態のファンでは、十分な空気がラジエーターを通して循環せず、困難な状況ではエンジンを冷却できません。この例では、ラジエーターからファンまでの距離、翼先端間隔およびブレードピッチ角度という 3 つの性能因子の効果を調査する実験を計画する方法を初めに示します。次に、テストデータを使用して、875 フィート $^{3}$ /分という目標を超える気流が発生する設計になるように、各因子に最適な値を推定する方法を示します。最後に、新しい設計では仕様に従って気流が生み出されるファンの比率が製造したファンの 99.999% を超えることを、シミュレーションを使用して検証する方法を示します。

問題の定義

この例では、ラジエーターから十分な量の空気を供給できないため困難な状況 (交通渋滞や高温) ではエンジンを冷却できないエンジン冷却ファンの設計に対処します。困難な状況でエンジンを冷却するには少なくとも 875 フィート $^{3}$ /分の気流が必要であると推定したとします。現在の設計を評価し、目標の気流を実現できる別の設計を開発する必要があります。

冷却ファンの性能評価

この例を実行すると提供される標本データを読み込みます。

load("OriginalFan.mat")

データは、既存の冷却ファンの性能に関する 10,000 個の測定値 (過去の生産データ) から構成されています。

冷却ファンの性能を分析するため、データをプロットします。

plot(originalfan)
xlabel("Observation")
ylabel("Max Airflow (ft^3/min)")
title("Historical Production Data")

Figure contains an axes object. The axes object with title Historical Production Data, xlabel Observation, ylabel Max Airflow ( f t Cubed baseline / m i n ) contains an object of type line.

データの中心は約 842 フィート $^{3}$ /分になっており、ほとんどの値は約 8 フィート $^{3}$ /分の範囲に収まっています。しかし、このプロットからは基礎となるデータの分布を理解できません。ヒストグラムをプロットし、正規分布をデータに当てはめます。

figure
histfit(originalfan) % Plot histogram with normal distribution fit
format shortg
xlabel("Airflow (ft^3/min)")
ylabel("Frequency (counts)")
title("Airflow Histogram")

Figure contains an axes object. The axes object with title Airflow Histogram, xlabel Airflow ( f t Cubed baseline / m i n ), ylabel Frequency (counts) contains 2 objects of type bar, line.

pd = fitdist(originalfan,"normal") % Fit normal distribution to data

pd = 
  NormalDistribution

  Normal distribution
       mu = 841.652   [841.616, 841.689]
    sigma =  1.8768   [1.85114, 1.90318]

fitdistは、正規分布をデータに当てはめ、データからパラメーターを推定します。平均気流速度の推定値は 841.652 フィート $^{3}$ /分、平均気流速度の 95% 信頼区間は (841.616, 841.689) です。この推定値から、現在のファンの気流速度は必要な 875 フィート $^{3}$ /分に近くないことが明らかになります。目標の気流を実現するようにファンの設計を改善する必要があります。

ファンの性能に影響を与える因子の判別

実験計画法 (DOE) を使用して、冷却ファンの性能に影響を与える因子を評価します。応答は冷却ファンの気流速度 (フィート $^{3}$ /分) です。次の因子を変更および制御できるとします。

ラジエーターからの距離
ピッチ角度
翼先端間隔

一般に、流体系の挙動は非線形です。したがって、応答曲面計画を使用して因子間の非線形交互作用を推定します。コード化 (正規化) した変数 [-1, 0, +1] でボックスベーンケン計画の実験実行を生成します。ボックスベーンケン計画を参照してください。

CodedValue = bbdesign(3)

CodedValue = 15×3

    -1    -1     0
    -1     1     0
     1    -1     0
     1     1     0
    -1     0    -1
    -1     0     1
     1     0    -1
     1     0     1
     0    -1    -1
     0    -1     1
      ⋮

1 列目はラジエーターからの距離、2 列目はピッチ角度、3 列目は翼先端間隔に対応しています。次の最小値および最大値で変数の効果をテストするとします。

ラジエーターからの距離: 1 ～ 1.5 インチ

ピッチ角度: 15 ～ 35°

翼先端間隔: 1 ～ 2 インチ

実行の順序を無作為化し、コード化した設計値を実世界の単位に変換し、指定された順序で実験を行います。

runorder = randperm(15);     % Random permutation of the runs
bounds = [1 1.5;15 35;1 2];  % Min and max values for each factor

RealValue = zeros(size(CodedValue));
for i = 1:size(CodedValue,2) % Convert coded values to real-world units
    zmax = max(CodedValue(:,i));
    zmin = min(CodedValue(:,i));
    RealValue(:,i) = interp1([zmin zmax],bounds(i,:),CodedValue(:,i));
end

実験が終わったときに次の応答値を変数 TestResult に格納するとします。

TestResult = [837 864 829 856 880 879 872 874 834 833 860 859 874 876 875]';

設計値と応答を table に保存します。

Expmt = table(runorder', CodedValue(:,1), CodedValue(:,2), CodedValue(:,3), ...
    TestResult,'VariableNames',{'RunNumber','D','P','C','Airflow'});

設計値と応答を表示します。

disp(Expmt)

    RunNumber    D     P     C     Airflow
    _________    __    __    __    _______

        6        -1    -1     0      837  
        3        -1     1     0      864  
       11         1    -1     0      829  
        7         1     1     0      856  
       14        -1     0    -1      880  
        8        -1     0     1      879  
        5         1     0    -1      872  
       15         1     0     1      874  
        1         0    -1    -1      834  
        2         0    -1     1      833  
        4         0     1    -1      860  
       13         0     1     1      859  
        9         0     0     0      874  
       10         0     0     0      876  
       12         0     0     0      875

D は Distance、P は Pitch、C は Clearance を表します。実験テストの結果から、気流速度は因子の値の変化に大きく影響されることがわかります。また、4 つの実験実行では、875 フィート $^{3}$ /分という目標の気流速度を満たしているか超えています (実行 2、4、12、14)。しかし、どの実行が最適であるかは明確ではありません。さらに、因子の変動に対する設計の耐性がどの程度であるか不明確です。現在の実験データに基づいてモデルを作成し、このモデルを使用して最適な要因設定を推定します。

冷却ファンの性能の改善

ボックスベーンケン計画では、非線形 (2 次) 効果をテストできます。2 次モデルの形式は次のようになります。

$\begin{array}{l} A F & = β_{0} + β_{1} * D i s t a n c e + β_{2} * P i t c h + β_{3} * C l e a r a n c e + β_{4} * D i s t a n c e * P i t c h \\ + β_{5} * D i s t a n c e * C l e a r a n c e + β_{6} * P i t c h * C l e a r a n c e + β_{7} * D i s t a n c e^{2} \\ + β_{8} * P i t c h^{2} + β_{9} * C l e a r a n c e^{2}, \end{array}$

ここで、AF は気流速度、 $β_{i}$ は項 $i$ の係数です。関数fitlmを使用して、このモデルの係数を推定します。

mdl = fitlm(Expmt,"Airflow~D*P*C-D:P:C+D^2+P^2+C^2");

(値を正規化した) 係数の大きさを棒グラフに表示します。

figure
h = bar(mdl.Coefficients.Estimate(2:10));
set(h,"facecolor",[0.8 0.8 0.9])
legend("Coefficient")
set(gcf,"units","normalized","position",[0.05 0.4 0.35 0.4])
set(gca,"xticklabel",mdl.CoefficientNames(2:10))
ylabel("Airflow (ft^3/min)")
xlabel("Normalized Coefficient")
title("Quadratic Model Coefficients")

Figure contains an axes object. The axes object with title Quadratic Model Coefficients, xlabel Normalized Coefficient, ylabel Airflow ( f t Cubed baseline / m i n ) contains an object of type bar. This object represents Coefficient.

この棒グラフは、"Pitch" と "Pitch $^{2}$ " が支配的な因子であることを示しています。応答曲面プロットを生成すると、複数の入力変数と 1 つの出力変数の関係を調べることができます。plotSliceを使用して、モデル mdl の応答曲面プロットを対話的に生成します。

plotSlice(mdl)

Figure Prediction Slice Plots contains 3 axes objects and other objects of type uimenu, uicontrol. Axes object 1 contains 5 objects of type line. Axes object 2 contains 5 objects of type line. Axes object 3 contains 5 objects of type line.

このプロットは、気流とピッチの非線形関係を示しています。青い破線を動かして、さまざまな因子が気流に与える影響を確認します。plotSlice を使用して最適な要因設定を判別できますが、Optimization Toolbox™ を使用してタスクを自動化することもできます。

要因設定の最適化

問題ベースの最適化ワークフロー (Optimization Toolbox)を使用して、最適な要因設定を求めます。最初に、最大化の最適化問題を定義します。

prob = optimproblem("ObjectiveSense","max");

x2fxを使用して予測子行列を計画行列に変換することで目的関数を記述します。結果をモデル係数の推定値で乗算します。

fun = @(x) x2fx(x,"quadratic")*mdl.Coefficients.Estimate;

factors という名前の最適化変数を作成します。下限は –1、上限は 1 で、3 つの因子を表す 3 つの成分で構成されます。

factors = optimvar("factors",1,3,LowerBound=-1,UpperBound=1);

関数fcn2optimexpr (Optimization Toolbox)を使用して、目的関数を factors の最適化式に変換します。

objective = fcn2optimexpr(fun,factors);

目的関数の式を問題 prob に配置します。

prob.Objective = objective;

初期点を実験計画テスト行列の中央、つまりベクトル [0 0 0] に設定します。問題ベースのアプローチでは、初期点は変数名を名前のフィールドとしてもつ構造体でなければなりません。

x0 = struct("factors",[0 0 0]);

最適な設計を求めます。

[sol,fval] = solve(prob,x0);

Solving problem using fmincon.

Feasible point with lower objective function value found.


Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

<stopping criteria details>

結果を実世界の単位に変換します。

maxloc = (sol.factors + 1)';
maxloc = bounds(:,1)+maxloc .* ((bounds(:,2) - bounds(:,1))/2);
fprintf("Optimal Values:\n" + ...
    "Distance  Pitch    Clearance  Airflow\n" + ...
    " %g       %g    %g         %g\n",maxloc',fval);

Optimal Values:
Distance  Pitch    Clearance  Airflow
 1       27.2747    1         882.257

最適化の結果から、新しいファンをラジエーターから 1 インチの位置に配置し、ピッチ角度を 27.3 にして、ファンの翼先端とシュラウドの間隔を 1 インチにすればよいことがわかります。

ピッチ角度は気流に大きい影響を与えるので、さらに分析を行って、27.3°というピッチ角度が最適であるか検証します。

load("AirflowData.mat")
tbl = table(pitch,airflow);
mdl2 = fitlm(tbl,"airflow~pitch^2");
mdl2.Rsquared.Ordinary

ans = 
      0.99632

この結果は、気流に対するピッチの影響が 2 次モデルによって十分に説明されることを示しています。

ピッチ角度と気流の関係をプロットし、近似モデルを重ねます。

figure
plot(pitch,airflow,".r") 
hold on
ylim([840 885])
line(pitch,mdl2.Fitted,"color","b") 
title("Fitted Model and Data")
xlabel("Pitch angle (degrees)") 
ylabel("Airflow (ft^3/min)")
legend("Test data","Quadratic model","Location","se")
hold off

Figure contains an axes object. The axes object with title Fitted Model and Data, xlabel Pitch angle (degrees), ylabel Airflow ( f t Cubed baseline / m i n ) contains 2 objects of type line. One or more of the lines displays its values using only markers These objects represent Test data, Quadratic model.

最大の気流に対応するピッチの値を求めます。

pitch(find(airflow==max(airflow)))

ans = 
    27

さらに分析した結果、27.3°というピッチ角度が最適であることが確認されました。

改善した冷却ファンの設計は、気流の要件を満たしています。また、設計で変更できる因子に基づいてファンの性能を十分に近似するモデルが得られました。感度分析を実行して、ファンの性能が製造と設置の変動性の影響を受けないことを確認します。

感度分析

過去の経験に基づくと、製造の不確実性が次のようになっているとします。

table(["Distance from radiator";"Blade pitch angle";"Blade tip clearance"],...
    ["1.00 +/- 0.05 inch";"27.3 +/- 0.25 degrees";"1.00 +/- 0.125 inch"],...
    ["1.00 +/- 0.20 inch";"0.227 +/- 0.028 degrees";"-1.00 +/- 0.25 inch"],...
    'VariableNames',{'Factor' 'Real Values' 'Coded Values'})

ans=3×3 table
             Factor                   Real Values                Coded Values       
    ________________________    _______________________    _________________________

    "Distance from radiator"    "1.00 +/- 0.05 inch"       "1.00 +/- 0.20 inch"     
    "Blade pitch angle"         "27.3 +/- 0.25 degrees"    "0.227 +/- 0.028 degrees"
    "Blade tip clearance"       "1.00 +/- 0.125 inch"      "-1.00 +/- 0.25 inch"

このように因子が変動しても目標の気流の付近で設計が堅牢であるかどうか検証します。シックスシグマの指針では、1,000,000 個のファンごとに欠陥品が 3.4 個を超えないという欠陥率を目標にします。つまり、ファンは 99.999% の確率で 875 フィート $^{3}$ /分という目標を達成しなければなりません。

モンテカルロシミュレーションを使用して設計を検証できます。指定した許容誤差で 3 つの因子について 10,000 個の乱数を生成します。近似モデル mdl のノイズに比例するノイズ変数 (モデルの RMS 誤差) を含めます。モデルの係数はコード化された変数なので、dist、pitch および clearance はコード化された定義を使用して生成する必要があります。

dist = random("normal",sol.factors(1),0.20,[10000 1]);
pitch = random("normal",sol.factors(2),0.028,[10000 1]);
clearance = random("normal",sol.factors(3),0.25,[10000 1]);
noise = random("normal",0,mdl2.RMSE,[10000 1]);

モデルを使用して、10,000 個の無作為な因子の組み合わせについて気流を計算します。

simfactor = [dist pitch clearance];
X = x2fx(simfactor,"quadratic");

ノイズ (モデルで考慮しないデータの変動) をモデルに追加します。

simflow = X*mdl.Coefficients.Estimate+noise;

ヒストグラムを使用して、モデルから予測した気流における変動を評価します。平均と標準偏差を推定するため、正規分布をデータに当てはめます。

pd = fitdist(simflow,"normal");
figure
histfit(simflow) 
hold on
text(pd.mu+2,300,["Mean: " num2str(round(pd.mu))])
text(pd.mu+2,280,["Standard deviation: " num2str(round(pd.sigma))])
hold off
xlabel("Airflow (ft^3/min)")
ylabel("Frequency")
title("Monte Carlo Simulation Results")
hold off

Figure contains an axes object. The axes object with title Monte Carlo Simulation Results, xlabel Airflow ( f t Cubed baseline / m i n ), ylabel Frequency contains 4 objects of type bar, line, text.

結果は有望なようです。平均気流は 882 フィート $^{3}$ /分なので、ほとんどのデータについて 875 フィート $^{3}$ /分より優れているようです。

気流が 875 フィート $^{3}$ /分以下になる確率を求めます。

format long
pfail = cdf(pd,875)

pfail = 
     1.454076873131737e-07

pass = (1-pfail)*100

pass = 
  99.999985459231269

この設計は、99.999% の確率で少なくとも 875 フィート $^{3}$ /分の気流を実現しているようです。

シミュレーション結果を使用して、工程能力を推定します。

S = capability(simflow,[875.0 890])

S = struct with fields:
       mu: 8.822983078353724e+02
    sigma: 1.422865134786495
        P: 0.999999823569868
       Pl: 1.454076873131737e-07
       Pu: 3.102244433021162e-08
       Cp: 1.757018243598422
      Cpl: 1.709768001886212
      Cpu: 1.804268485310632
      Cpk: 1.709768001886212

pass = (1-S.Pl)*100

pass = 
  99.999985459231269

Cp の値は 1.75 です。Cp が 1.6 以上の場合、工程は高品質であると考えられます。Cpk は Cp の値と似ており、工程の特性値が目標範囲の中心付近にあることを示しています。次に、この設計を実装します。監視して設計プロセスを検証し、冷却ファンが高い性能を発揮することを確認します。

改善した冷却ファンの製造の管理

管理図を使用すると、新しいファンの製造工程と設置工程を監視できます。初めの 30 日について新しい冷却ファンの生産を評価します。当初は、1 日に 5 枚の冷却ファンが生産されていました。まず、新しい工程の標本データを読み込みます。

load("spcdata.mat")

"X" バーと "S" の管理図をプロットします。

figure
controlchart(spcflow,"chart",{'xbar','s'}) % Reshape the data into daily sets
xlabel("Day")

Figure contains 2 axes objects. Axes object 1 with title Control charts, xlabel Day, ylabel XBAR contains 4 objects of type line. One or more of the lines displays its values using only markers These objects represent Data, Violation, Center, LCL/UCL. Axes object 2 with ylabel S contains 4 objects of type line. One or more of the lines displays its values using only markers

この結果によると、管理限界の違反はなく、データの時間変化に含まれているパターンは非無作為なので、製造工程は統計的に管理されています。また、データに対して能力分析を実行して工程を評価することもできます。

[row,col] = size(spcflow);
S2 = capability(reshape(spcflow,row*col,1),[875.0 890])

S2 = struct with fields:
       mu: 8.821061141685465e+02
    sigma: 1.423887508874697
        P: 0.999999684316149
       Pl: 3.008932155898586e-07
       Pu: 1.479063578225176e-08
       Cp: 1.755756676295137
      Cpl: 1.663547652525458
      Cpu: 1.847965700064817
      Cpk: 1.663547652525458

pass = (1-S.Pl)*100

pass = 
  99.999985459231269

Cp 値の 1.755 は、推定値の 1.73 と非常に近い値です。Cpk 値の 1.66 は、Cp 値より小さい値です。しかし、1.33 未満の Cpk は工程が一方の工程限界に向かって大幅にシフトしていることを示しており、唯一の懸念事項になっています。この工程は十分に範囲内に収まっており、99.999% を超える確率で目標の気流 (875 フィート $^{3}$ /分) を実現しています。

参考

シックス シグマ手法の計画によるエンジン冷却ファンの改善