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

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

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

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/мин, а 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 до 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'});

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

Повышение производительности вентилятора охлаждения

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

где AF - расход воздуха, а _Bi - коэффициент для термина 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')

Гистограмма показывает, что шаг и Pitch2 являются доминирующими факторами. Можно просмотреть взаимосвязь между несколькими входными переменными и одной выходной переменной, создав график поверхности отклика. Использовать plotSlice для создания графиков поверхности отклика для модели mdl в интерактивном режиме.

plotSlice(mdl)

График показывает нелинейное соотношение воздушного потока с шагом. Переместите синие пунктирные линии вокруг и просмотрите влияние различных факторов на поток воздуха. Хотя вы можете использовать plotSlice для определения оптимальных параметров коэффициента можно также использовать панель инструментов оптимизации для автоматизации задачи.

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

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

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

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

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

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

rng('default')

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

Рассчитайте воздушный поток для 10000 случайных комбинаций коэффициентов с помощью модели.

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 ^ft3/мин или ниже.

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

pfail =

     1.509289008603141e-07


pass =

  99.999984907109919

Конструкция, по-видимому, обеспечивает по меньшей мере 875 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-bar и 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 ^{фут3/мин}) более 99,999% времени.