ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

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

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

問題の定義

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

冷却ファンの性能評価

標本データを読み込みます。

load(fullfile(matlabroot,'help/toolbox/stats/examples','OriginalFan.mat'))

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

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

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

データの中心は約 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')

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 =

    -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
     0     1    -1
     0     1     1
     0     0     0
     0     0     0
     0     0     0

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]';

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

disp({'Run Number','Distance','Pitch','Clearance','Airflow'})
disp(sortrows([runorder' RealValue TestResult]))
'Run Number'    'Distance'    'Pitch'    'Clearance'    'Airflow'

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

設計値と応答を table に格納します。

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

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

冷却ファンの性能の改善

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

AF = β0+β1*Distance+β2*Pitch+β3*Clearance+β4*Distance*Pitch+β5*Distance*Clearance+β6*Pitch*Clearance+β7*Distance2+β8*Pitch2+β9*Clearance2,

ここで、AF は気流速度、Bi は項 i の係数です。Statistics and Machine Learning Toolbox で関数 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')

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

plotSlice(mdl)

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

制約付き最適化関数 fmincon を使用して、最適な要因設定を求めます。

目的関数を記述します。

f = @(x) -x2fx(x,'quadratic')*mdl.Coefficients.Estimate;

目的関数は、データを近似する 2 次応答曲面です。fmincon を使用して負の気流を最小化することと、元の目的関数の最大化は同じになります。制約は、(コード化した値で) テストした上限と下限です。初期の開始点を実験計画テスト行列の中央に設定します。

lb = [-1 -1 -1]; % Lower bound					
ub = [1 1 1];    % Upper bound                       
x0 = [0 0 0];    % Starting point
[optfactors,fval] = fmincon(f,x0,[],[],[],[],lb,ub,[]); % Invoke the solver
Local minimum found that satisfies the constraints.

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

結果を最大化問題および実世界の単位に変換します。

maxval = -fval;
maxloc = (optfactors + 1)';
bounds = [1 1.5;15 35;1 2];
maxloc=bounds(:,1)+maxloc .* ((bounds(:,2) - bounds(:,1))/2);
disp('Optimal Values:')
disp({'Distance','Pitch','Clearance','Airflow'})
disp([maxloc' maxval])
Optimal Values:
    'Distance'    'Pitch'    'Clearance'    'Airflow'

            1       27.275            1       882.26

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

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

load(fullfile(matlabroot,'help/toolbox/stats/examples','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

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

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

    27

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

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

感度分析

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

因子実際の値コード化した値
ラジエーターからの距離1.00 +/- 0.05 インチ1.00 +/- 0.20 インチ
ブレード ピッチ角度27.3 +/- 0.25°0.227 +/- 0.028°
翼先端間隔1.00 +/- 0.125 インチ-1.00 +/- 0.25 インチ

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

モンテカルロ シミュレーションを使用して設計を検証できます。指定した許容誤差で 3 つの因子について 10,000 個の乱数を生成します。まず、異なる実行でも結果が一致するように乱数発生器の状態を設定します。

rng('default')

モンテカルロ シミュレーションを実行します。近似モデル mdl のノイズに比例するノイズ変数 (モデルの RMS 誤差) を含めます。モデルの係数はコード化された変数なので、distpitch および clearance はコード化された定義を使用して生成する必要があります。

dist = random('normal',optfactors(1),0.20,[10000 1]);
pitch = random('normal',optfactors(2),0.028,[10000 1]);
clearance = random('normal',optfactors(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');
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')

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

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

format long
pfail = cdf(pd,875)
pass = (1-pfail)*100
pfail =

     1.509289008603141e-07


pass =

  99.999984907109919

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

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

S = capability(simflow,[875.0 890])
pass = (1-S.Pl)*100
S = 

       mu: 8.822982645666709e+02
    sigma: 1.424806876923940
        P: 0.999999816749816
       Pl: 1.509289008603141e-07
       Pu: 3.232128339675335e-08
       Cp: 1.754623760237126
      Cpl: 1.707427788957002
      Cpu: 1.801819731517250
      Cpk: 1.707427788957002


pass =

  99.9999849071099

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

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

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

load(fullfile(matlabroot,'help/toolbox/stats/examples','spcdata.mat'))

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

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

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

[row,col] = size(spcflow);
S2 = capability(reshape(spcflow,row*col,1),[875.0 890])
pass = (1-S.Pl)*100
S2 = 

       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 =

  99.9999699106784

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