Документация

Определение проблемы

Этот пример обращается к проекту вентилятора механизма, который не может вытянуть достаточно воздуха через излучателя, чтобы сохранить механизм холодным во время трудных условий, таких как трафик останавливать-и-идти или жаркая погода). Предположим, что вы оцениваете, что вам нужен поток воздуха по крайней мере 875 ^ft3/min, чтобы сохранить механизм холодным во время трудных условий. Необходимо оценить текущий проект и разработать альтернативный дизайн, который может достигнуть целевого потока воздуха.

Оцените эффективность вентилятора

Загрузите выборочные данные.

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 ^ft3/min и большинство значений находятся в пределах области значений приблизительно 8 ^ft3/min. График не говорит очень о базовом распределении данных, как бы то ни было. Постройте гистограмму и соответствуйте нормальному распределению к данным.

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 ^ft3/min, и 95%-й доверительный интервал для средней скорости потока воздуха (841.616, 841.689). Эта оценка проясняет, что текущий вентилятор не близко к необходимым 875 ^ft3/min. Существует потребность улучшить проект вентилятора, чтобы достигнуть целевого потока воздуха.

Определите факторы, которые влияют на эффективность вентилятора

Оцените факторы, которые влияют на эффективность вентилятора с помощью проекта экспериментов (DOE). Ответ является уровнем потока воздуха вентилятора ^(ft3/min). Предположим, что факторы, которые можно изменить и управление:

В общем случае гидросистемы имеют нелинейное поведение. Поэтому используйте проект поверхности ответа, чтобы оценить любые нелинейные взаимодействия среди факторов. Сгенерируйте экспериментальные запуски для проекта Поля-Behnken в закодированных (нормированных) переменных [-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 к 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

Сохраните расчетные значения и ответ в a table.

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

D обозначает DistanceP обозначает Pitch, и C обозначает Clearance. На основе экспериментальных результатов испытаний уровень потока воздуха чувствителен к изменяющимся значениям факторов. Кроме того, четыре экспериментальных запуска выполняют или превышают целевой уровень потока воздуха 875 ^ft3/min (запускается 2, 4,12, и 14). Однако это не ясно, который, если таковые имеются, этих запусков оптимальный. Кроме того, не очевидно, насколько устойчивый проект к изменению факторов. Создайте основанное на модели на текущих экспериментальных данных и используйте модель, чтобы оценить оптимальные факторные настройки.

Улучшайте производительность вентилятора

Проект Поля-Behnken позволяет вам протестировать на нелинейные (квадратичные) эффекты. Форма квадратичной модели:

где AF является уровнем потока воздуха, и _Bi является коэффициентом для термина i. Оцените коэффициенты этой модели с помощью fitlm функция из Statistics and Machine Learning Toolbox.

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 и Pitch ² являются доминирующими факторами. Можно посмотреть на отношение между несколькими входными переменными и одной выходной переменной путем генерации объемной поверхностной диаграммы ответа. Используйте plotSlice сгенерировать объемные поверхностные диаграммы ответа для модели mdl в интерактивном режиме.

plotSlice(mdl)

График показывает нелинейное отношение потока воздуха с тангажом. Переместите синие пунктирные линии вокруг и смотрите влияние, которое различные факторы оказывают на поток воздуха. Несмотря на то, что можно использовать plotSlice чтобы определить оптимальные факторные настройки, можно также использовать Optimization Toolbox, чтобы автоматизировать задачу.

Найдите оптимальные факторные настройки с помощью ограниченной оптимизационной функции fmincon.

Запишите целевую функцию.

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

Целевая функция является квадратичной подгонкой поверхности ответа к данным. Минимизация отрицательного использования потока воздуха 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

Результат оптимизации предлагает разместить новый вентилятор в одном дюйме от излучателя с однодюймовым разрешением между советами лопастей вентилятора и оболочки.

Поскольку угол тангажа оказывает такое значительное влияние на поток воздуха, выполните дополнительный анализ, чтобы проверить, что 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

Результаты показывают, что квадратичная модель объясняет эффект тангажа на потоке воздуха хорошо.

Постройте угол тангажа против потока воздуха и наложите подобранную модель.

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 углов тангажа степени оптимальны.

Улучшенный проект вентилятора удовлетворяет требования потока воздуха. У вас также есть модель, которая аппроксимирует эффективность вентилятора хорошо на основе факторов, которые можно изменить в проекте. Убедитесь, что эффективность вентилятора устойчива к изменчивости в производстве и установке путем выполнения анализа чувствительности.

Анализ чувствительности

Предположим, что на основе исторического опыта производственная неопределенность следующие.

Проверьте, что эти изменения факторов включат, чтобы обеспечить устойчивый проект вокруг целевого потока воздуха. Философия Шести целей Сигмы скорость дефектообразования не больше, чем 3.4 на 1 000 000 вентиляторов. Таким образом, вентиляторы должны поразить 875 целей ^ft3/min 99,999% времени.

Можно проверить проект с помощью симуляции Монте-Карло. Сгенерируйте 10 000 случайных чисел для трех факторов с заданным допуском. Во-первых, установите состояние генераторов случайных чисел, таким образом, результаты сопоставимы через различные запуски.

rng('default')

Выполните симуляцию Монте-Карло. Включайте шумовую переменную, которая пропорциональна шуму в подобранной модели, mdl (то есть, ошибка RMS модели). Поскольку коэффициенты модели находятся в закодированных переменных, необходимо сгенерировать distТангаж, и 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 ^ft3/min и, кажется, лучше, чем 875 ^ft3/min для большинства данных.

Определите вероятность, что поток воздуха в 875 ^ft3/min или ниже.

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

pfail =

     1.509289008603141e-07


pass =

  99.999984907109919

Проект, кажется, достигает по крайней мере 875 ^ft3/min потока воздуха 99,999% времени.

Используйте результаты симуляции, чтобы оценить возможность процесса.

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. Cpk похоже на Cp значение, которое указывает, что процесс сосредоточен. Теперь реализуйте этот проект. Контролируйте его, чтобы проверить процесс проектирования и гарантировать, что вентилятор поставляет высококачественную эффективность.

Производство управления улучшенного вентилятора

Можно контролировать и оценить производственный процесс и процесс установки нового вентилятора, использующего диаграммы управления. Оцените первые 30 дней производства нового вентилятора. Первоначально, пять вентиляторов в день были произведены. Во-первых, загрузите выборочные данные от нового процесса.

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 значение. Однако только Cpk значение меньше чем 1,33, который указывает, что процесс, переключенный значительно к одному из пределов процесса, являются беспокойством. Процесс хорошо в пределах, и он достигает целевого потока воздуха (875 ^ft3/min) больше чем 99,999% времени.