Планирование экспериментов и обработка данных активного эксперимента в Excel

ФЕДЕРАЛЬНОЕ АГЕНТСТВО ПО ОБРАЗОВАНИЮ РОССИЙСКОЙ ФЕДЕРАЦИИ РОССИЙСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ НЕФТИ И ГАЗА имени И.М. Губкина Кафедра прикладной математики и компьютерного моделирования

Калинина Э.В., Ретинская И.В. Планирование экспериментов и обработка данных активного эксперимента в Excel Методическое пособие

Москва 2007

1

Содержание 1. Регрессионные модели и организация экспериментов для получения данных ............................................................................................................................. 3 1.1 Независимые и зависимые переменные ............................................................. 3 1.2 Регрессионный анализ. Основная постановка задачи .................................... 5 1.3 Пассивный и активный эксперименты. Планирование эксперимента .. 9 1.4 Полный факторный эксперимент и дробные реплики ................................. 14 1.4.1. Критерий Стьюдента для проверки значимости параметров модели... 17 1.4.2. Проверка адекватности модели по критерию Фишера ............................ 19 1.5 Дробные реплики полного факторного эксперимента................................ 22 1.6 Примеры применения факторных планов и дробных реплик .................. 27 1.6.1 ......... Планирование эксперимента и математическое моделирование при поиске оптимального состава консервационного масла. ................................... 27 Полный факторный план .......................................................................................... 27 1.6.2 Разработка композиции присадок к консервационному маслу. Дробный факторный эксперимент .......................................................................... 34 1.7 Планирование эксперимента на диаграммах................................................. 41 состав - свойство.......................................................................................................... 41 1.7.1 Симплексная система координат.................................................................... 41 1.7.2 Симплекс - решетчатые планы ...................................................................... 44 1.7.3 Расчет параметров моделей при использовании симплекс - решетчатых планов ............................................................................................................................ 49 1.8. Информатизация обработки данных активного эксперимента в Excel... 51 Лабораторная работа по планированию эксперимента.................................... 60 Задача 1......................................................................................................................... 60 Задача 2......................................................................................................................... 61

2

1. Регрессионные модели и организация экспериментов для получения данных 1.1

Независимые и зависимые переменные

Экспериментальные

данные

широко

используются

в

науке

и

в

промышленности. При этом цели исследования в этих двух областях могут отличаться. Основная цель обработки экспериментальных данных в рамках научного исследования - описание статистической зависимости воздействия определенного фактора на изучаемую зависимую переменную. Целью и результатом является построение адекватной математической модели. Основная цель дорогостоящих экспериментов в промышленности обычно заключается в извлечении максимального количества

объективной

информации

о

влиянии

изучаемых

факторов

на

производственный процесс. При этом, очевидно, что число опытов должно быть как можно меньшим. В промышленности учет эффектов взаимодействия факторов часто считается излишним, при этом пытаются выявить лишь существенно влияющие факторы. Выделяют два основных типа переменных: независимые переменные (или факторы, или управляющие параметры, или переменные

(или

переменные-отклики,

входные переменные) и зависимые

или

выходные

переменные).

Под

независимыми переменными понимают такие переменные, для которых можно устанавливать желаемые значения, либо те, которые можно только наблюдать, но не управлять ими. В результате изменений, проведённых намеренно или происшедших с независимыми переменными случайно, появляется эффект, который передается на другие переменные-отклики. Разделение на факторы и отклики не всегда вполне чётко и зависит от целей исследования. Практически же роли переменных легко различимы. Примеры независимых и независимых переменных: • количество удобрений и количество выпавших осадков (факторы) при изучении урожайности (отклик); 3

• давление, температура кипящего слоя (факторы), при изучении толщины поверхностного слоя кремния на подложке при производстве микросхем (отклик); • давление, под которым производится розлив (фактор), на потери сиропа при изготовлении безалкогольных напитков (отклик). При изучении взаимосвязей между переменными величинами независимые входные переменные или факторы принято обозначать

как

x1 ,.., x n , (n - их

количество). Выходные зависимые переменные или отклики обозначаются обычно как

y1 ,... y p , (p - их количество).

Каждый набор чисел { x1 ,.., x n } можно

рассматривать как координаты точки в n- мерном факторном пространстве Rn или как вектор. Ограничения на величины входов xi , i=1÷n выделяют некоторую область в этом же пространстве. В частности, если ограничения имеют вид ci ≤ xi ≤ bi , то область представляет собой брус (прямоугольник, параллелепипед) в пространстве Rn. Рассмотрим пример, когда объектом исследования является смесь q различных компонентов, и изучается взаимосвязь некоторого свойства смеси с распределением в ней пропорций между компонентами, тогда область задания независимых переменных представляет собой геометрическую фигуру – симплекс. Переменные xi (i=1,2, . . ., q) таких систем являются пропорциями или концентрациями (доли единицы) i -тыx компонентов смеси и удовлетворяют условию нормирования:

∑ xi = 1, 1≤i ≤q

(xi ≥0).

Геометрическое место точек, удовлетворяющее условию нормирования, представляет собой

(q—1) - мерный правильный симплекс. (Правильный

симплекс с q вершинами в (q — 1)-мерном пространстве, треугольник для q = 3, тетраэдр для q = 4 и т. д.). Каждой точке такого симплекса соответствует смесь определенного состава, и, наоборот, любой комбинации относительных содержаний q компонентов соответствует определенная точка симплекса. 4

1.2

Регрессионный анализ. Основная постановка задачи

Регрессионный анализ позволяет установить связи между факторами { x1 ,.., x n } и

выбранным

откликом yi. Факторы в регрессионном анализе считаются

неслучайными величинами, а зависимая переменная

yi считается величиной

случайной. Природа случайности связана либо со случайным характером ошибок измерения отклика, либо значения отклика могут зависеть ещё и от ряда неконтролируемых факторов. Регрессионный анализ [9,11,12] основывается на предположении о том, что величина y является случайной величиной, распределение которой зависит от величин { x1 ,.., x n }. Измеряемая величина y

может быть

представлена в виде двух слагаемых, первое из которых – неслучайная часть – определяется некоторой функцией от

x1 ,.., x n , а второе – случайная часть –

учитывает случайный характер y:

y = f ( x1 ,..., x n , a 0 , a1 ,...) + ε ,

(1.1.1)

где a0 , a1,.. - некоторые числовые коэффициенты, а ε- случайная величина, у которой математическое ожидание равно нулю, а дисперсия не зависит от величин факторов, постоянна и равна

σ2

(предположения

независимости измерений выходной

характеристики и их равной точности в факторном пространстве). Вид функции f считается известным, (говорят, модель идентифицирована)1, неизвестны лишь параметры модели a0 , a1,.. . Цель регрессионного анализа – найти по результатам экспериментов оценки параметров и оценить степень точности построенной зависимости. Оценки параметров находятся с помощью известного

метода наименьших

квадратов (МНК) [8,9,10,11]. Будем считать, что аппарат регрессионного метода и его составные части, такие как проверка значимости параметров модели по критерию Стьюдента

1

и проверка адекватности модели по критерию Фишера,

чаще всего на начальном этапе исследования модель идентифицируют как линейную

5

изучены

читателем. Тем не менее в последующих параграфах приводятся основные понятия и формулы. На примере линейной по параметрам модели напомним основные положения и формулы метода наименьших квадратов, необходимые для дальнейшего понимания. Уравнение линейной регрессии для n переменных-факторов имеет вид: y = a 0 + a1 x1 + ... + a n x n + ε ,

x0 = 1 ,

(1.1.2)

где, как и в (1.1), y- выход или моделируемый показатель, x1 ,.., x n - входы или факторы, ε- случайная погрешность с математическим ожиданием 0 и постоянной дисперсией σ2. Задача состоит в том, чтобы по данных N точек

( y j , x1 j ,..., x nj ),

имеющемуся набору экспериментальных j = 1, N , n – число факторов, оценить

коэффициенты a 0 , a1 ,..., a n . В методе наименьших квадратов в качестве оценок aˆ i коэффициентов ai выбирают те, которые минимизируют квадратичную форму Q=

j=N

∑

j =1

(y j −

i =n

∑ ai xij ) 2 .

(1.1.3)

i =0

r Введем матричные обозначения для векторов-столбцов выхода Y - результатов r r наблюдений, погрешностей ε , коэффициентов a и матрицы эксперимента или

матрицы плана наблюдений X

ε1

y1 r Y=

. . yN

r ,ε =

a0

x 01 ,..., xi1 ,..., x n1

x ,..., xi 2 ,..., x n 2 . r a , a = 1 , X = 02 , x 0 j = 1 , j = 1, N (1.1.4) . ...................... . εN an x 0 N ,..., xiN ,..., x nN

Тогда оценки метода наименьших квадратов находятся как решения так называемой системы нормальных уравнений относительно вектора коэффициентов r a: 6

r r X T Xa = X T Y .

(1.1.5)

Если матрица X T X невырожденная, существует обратная ей матрица, и решение системы нормальных уравнений (1.5) описывается уравнением r r aˆ = ( X T X ) −1 X T Y

(1.1.6)

r r Оценки МНК несмещенные, то есть ( M aˆ = a ), и эффективные в классе линейных оценок.

Одной из характеристик оценок является их ковариационная

матрица, имеющая вид:

r r r r r D(aˆ ) = M (aˆ − a )(aˆ − a ) T = σ2 ( X T X ) −1

(1.1.7)

r На главной диагонали матрицы D(aˆ ) стоят дисперсии соответствующих r оценок параметров a , а вне неё значения ковариаций между соответствующими компонентами вектора оценок. Подстановка полученных оценок в уравнение регрессии дает аналитический вид модели зависимости переменной y от независимых переменных:

yˆ = aˆ + aˆ x + ... + aˆn xn 0 11

(1.1.8)

Неотъемлемой частью процедуры построения уравнения регрессии является проверка статистической значимости модели регрессии и проверка её адекватности имеющимся экспериментальным данным. В тех случаях, когда дисперсия ошибки измерения σ2 неизвестна априори или нет предварительных оценок для её величины sˆ2 ( y) , в регрессионном анализе при проверках гипотез по статистическим критериям применяется так называемая остаточная сумма квадратов, которая является суммой квадратов отклонений вычисленных

значений по уравнению регрессии (1.1.8)

выходной переменной yˆ от их наблюденных значений yэ:

2 = N Sост ∑ ( yˆ j − y jэ )2 j =1 7

(1.1.9)

2 Если выбранная модель верна, остаточная сумма Sост может являться оценкой. sˆ2 ( y) для дисперсии ошибки измерения σ2 зависимой переменной: 2 /(N-(n+1)) sˆ2 ( y ) = Sост

,

(1.10)

Здесь, как и ранее, N - общее число наблюдений, n+1 – число оцениваемых параметров модели, а число степеней свободы равно f= N-(n+1). Следует отметить, что эта форма оценивания ошибки эксперимента используется в расчетах уравнения регрессии по программам пакетов Statistica и Excel.. Подобное оценивание не всегда приводит к хорошим результатам. Более 2 в процедурах проверки гипотез будет рассмотрено ниже. подробно применение Sост В практических задачах построения зависимостей, как правило, значение дисперсии σ2 ошибки измерения неизвестно. Для оценивания дисперсии обычно прибегают к проведению экспериментальных замеров выходной величины y при фиксированном (неизменном)

наборе значений факторов, а затем

полученные

результаты используют для оценки ошибки, распространяя полученное значение на всю область независимых переменных или факторов. (Этот прием вследствие

исходных

предположений

регрессионного

оправдан

анализа

о

независимости и равной точности измерений). Пусть

при

фиксированном

наборе

факторов

выполнено

m

параллельных измерений величины yk , k=1÷m, тогда оценкой дисперсии ошибки измерения является эмпирическая дисперсия sˆ2 ( y ) : m

∑(y

sˆ2 ( y ) = k =1

k

− y) 2 (m − 1)

, y = (1 / m)∑ y k

(1.1.11)

k

а оценкой среднеквадратического отклонения величина m

sˆ ( y ) =

∑(y k =1

k

− y) 2 ( m − 1)

8

,

(1.1.12)

при этом число степеней свободы равно f=m-1. В процедурах проверки статистических гипотез о значимости параметров модели и проверки адекватности должен учитывать способ оценивания ошибки эксперимента, так как число степеней свободы f в описанных способах разное. 1.3 Пассивный и активный эксперименты. Планирование эксперимента

Традиционный

путь

экспериментального

изучения

объектов

и

процессов состоит в проведении экспериментов (наблюдений) в точках факторного пространства, выбранных исследователем либо из неформальных соображений, либо произвольно, а затем — в обработке полученных экспериментальных данных с применением методов математической статистики. Сюда же относятся сбор и обработка данных о процессе в режиме нормальной эксплуатации, когда исследователь не вмешивается в процесс, а лишь регистрирует интересующие его параметры. Такой метод сбора данных называют пассивным экспер иментом . Построение моделей по данным пассивного эксперимента имеет следующие недостатки [8,11]: 1. При пассивном многофакторном

эксперименте трудно оценить

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

(факторы)

{ x1 ,.., x n },

вообще

говоря,

оказываются

(статистически связанными), поэтому соответствующие им

влияния на выход невозможно разделить. 3. Если изучаемый производственный процесс налажен (регулируем), то, как правило, факторы варьируются в очень узких пределах. В этом случае исследователь оказывается перед неразрешимой задачей: нужно описать выход процесса какой-то функцией, отражающей поведение процесса в разных

9

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

алгебраических

равнений — нормальных уравнений. Коэффициенты системы выражаются через значения

факторов;

если

факторы

зависимы

детерминированная или статистическая), матрица плохо обусловленной

(будь системы

то

зависимость

может оказаться

и решение будет некорректным. Для получения

корректных результатов необходимо провести дополнительный

анализ

и

обработку экспериментальных данных [5]. Перечисленных недостатков лишены модели, построенные по данным активного эксперимента. В активном эксперименте исследователь имеет возможность указать в каждом опыте значения (уровни) факторов и назначить последовательность проведения опытов. При этом исследователь всегда заинтересован

в получении

нужной

информации при

наименьших

затратах:

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

поиске

законов

или

зависимостей

между

факторами

использовали

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

Укажем два требования, выполнение которых предполагает теория активного эксперимента. 1. Воспроизводимость. факторов и повторим неравные

интервалы

Зафиксируем

уровни (значения)

несколько раз эксперимент времени).

Разброс

(вообще говоря,

значений отклика

воспроизводимость результатов: если разброс мал

всех через

характеризует

(не превышает заданного

заранее значения), имеется воспроизводимость результатов; если он велик, воспроизводимости нет. 2. Управляемость.

Предполагается, что

в

каждом опыте возможен

выбор тех уровней факторов или, иначе говоря, значений управляющих параметров, которые представляют интерес. Соответствующий объект называют управляемым, и именно при таком условии возможен активный эксперимент. В дальнейшем будем рассматривать лишь управляемые объекты, для которых выполнено требование воспроизводимости. Теория планирования эксперимента опирается на методы математической статистики и подробно отражена в литературе [11,13]. Ниже приведены лишь необходимые для дальнейшего изложения сведения. Пусть Ω x — допустимая область, в которой могут быть осуществлены r эксперименты. В каждой точке x из Ω x можно, вообще говоря, произвести не одно, а несколько наблюдений. Планом

эксперимента

называют

r множество точек x j из Ω x

( j = 1, N ), в которых проводят наблюдения, с указанием их числа r j в каждой r точке x j = ( x1 j ,..., x nj ), общее число наблюдений N= ∑ r j . j

В предположении нормальности измерений точность оценок параметров математической модели полностью описывается их ковариационной матрицей. Поэтому критерии оптимальности плана должны отражать ряд желательных свойств ковариационной матрицы. План обычно характеризуется некоторым функционалом

от

матрицы.

Наиболее 11

естественным

функционалом

от

ковариационной

матрицы

служит

ее

детерминант

r r det D(aˆ ) =⎢ D(aˆ ) ⎢,

называемый обобщенной дисперсией . Для линейных по параметрам моделей ковариационная матрица зависит только от выбора координат точек плана и не зависит от истинных значений r вектора a = ( a 0 , a1 ,..., a n ) оцениваемых параметров. Все статистические критерии качества плана эксперимента можно разбить на две группы. К первой относят критерии, связанные с оценкой параметров, ко второй — критерии и свойства планов, связанные с погрешностью в оценке модели. Среди критериев оценки параметров наиболее часто упоминается так называемая D-о п т и м а л ь н о с т ь, состоящая в том, что выбирают план, которому соответствует наименьший на множестве планов детерминант ковариационной

матрицы.

Эллипсоид

рассеяния

оценок

параметров

(коэффициентов регрессии) в этом случае имеет наименьший объем [13]. Невырожденные планы для оценки параметров модели (1.2) i =n

y = a 0 + a1 x1 + ... + a n x n + ε = ∑ a i xi ; i =0

x0 = 1

(1.3.1)

называют п л а н а м и первого порядка , или линейными. Область

планирования

выбирают

на

основании

неформальных

соображений, относящихся к существу изучаемой реальной задачи. Сначала нужно оценить границы областей определения факторов, а затем, пользуясь априорной информацией, выбрать локальную область для планирования эксперимента. Чаще всего область планирования Ω x задают одним из трех следующих типов ограничений. 1.

Пусть

zi —

натуральные

значения

независимой

переменной

(измеряемые в единицах времени, температуры, объема и т. д.). Область Z задается системой неравенств ci min ≤ zi ≤ ci max .

12

Для моделирования удобно перейти к безразмерным переменным, приняв за единицу половину длины интервала варьирования переменной ( ci max -

ci min ):

xi = [2 z i − (ci max + ci min )] /(ci max − ci min ) ,

(1.3.2)

так что − 1 ≤ xi ≤ +1,

i =1÷ n

(1.3.3)

Область, задаваемая неравенством (3.3), - это n-мерный куб; все 2 n его вершин находятся в точках, координаты которых равны либо +1, либо -1. В теории

планирования

эксперимента

переменные

xi

называют

кодированными. 2. Область Z в кодированных переменных представляет собой n-мерный n

{∑ xi2 ≤ 1}

шар:

i =1

n

3. Область {x : ∑ xi = 1, i =1

xi ≥ 0, i = 1 ÷ n} - это (n-1)-мерный правильный

симплекс; при n = 3 правильный симплекс — равносторонний треугольник, при n= 4 — равносторонний тетраэдр. При планировании экспериментов

прежде всего выбирается точка, в

окрестности которой будет изучаться поверхность отклика, и определяются интервалы варьирования переменных. Координаты в натуральных переменных преобразуют с помощью соотношений (1.3.2) в кодированные переменные так, что новые координаты выбранной точки становятся равными нулю. Эти координаты называют основным или нулевым уровнем рассматриваемых факторов. Выбор основного уровня опирается на неформальную априорную информацию об изучаемом процессе или объекте. В частности, если эта точка расположена на границе области Ω x или вблизи нее, то в качестве основного уровня

приходится

выбирать

другую

первоначальной.

13

точку,

сдвинутую

относительно

Если нельзя указать заранее, где условия наилучшие, основной уровень выбирают либо в центре области Ω x , либо случайно. 1.4 Полный факторный эксперимент и дробные реплики

Для оценки параметров линейной математической модели достаточно варьировать факторы на двух уровнях - верхнем и нижнем. Эксперимент, в котором реализуются все возможные сочетания факторов на двух уровнях, называют полным факторным экспериментом (ПФЭ) . Очевидно, число всех сочетаний на двух уровнях для n факторов равно N = 2n . Поэтому часто полный факторный эксперимент на двух уровнях обозначают как эксперимент типа 2n. Ниже записаны план и матрица планирования ПФЭ типа 2 2 для двух факторов x1 и x 2 (в строке указаны значения соответствующих факторов и отклика y ) для линейной модели. Область планирования при геометрической интерпретации — это квадрат с вершинами в точках x1 = ±1, x 2 =±1, а нулевой уровень расположен в начале координат. В вершинах квадрата производятся опыты, их результаты обозначены через y . Нетрудно независимых

построить переменных

аналогичные

планы

в

модели

линейной

для

любого

(3.1)

ПФЭ

числа

n

типа 2n

соответствует гиперкуб в n-мерном пространстве, вершины которого имеют координаты xi = ±1, i=1,n. Есть несколько простых способов построения полных факторных планов на двух уровнях. Наиболее удобной представляется процедура, при которой производится пропорциональное присвоение уровней +1 и -1 с периодичностью, соответствующей номеру переменной – для первого факторы – через один ( +1, -1, +1, -1 и т. д.), для второй переменной, через 2 шага (+1,+1,-1,-1, +1,+1, -1.-1. и т. д.), для третьей переменной – через 4 шага, затем через восемь, шестнадцать

и так далее. При такой процедуре общее

число экспериментов для n факторов составляет N=2n. 14

Таблица 1 − Полный факторный эксперимент для двух независимых переменных, варьируемых на двух уровнях (планирование типа 22) Переменные № опыта

x1

x2

Матрица планирования Х

y

x0

x1

x2

+1

-1

-1

1

-1

-1

y1

2

+1

-1

y2

+1

+1

-1

3

-1

+1

y3

+1

-1

+1

4

+1

+1

y4

+1

+1

+1

В таблице 2 приведена матрица планирования полного факторного эксперимента на двух уровнях для трех факторов. Пусть эксперимент проводится для оценки коэффициентов регрессии линейного уравнения

y = a 0 + a1 x1 + a 2 x 2 .

(1.4.1)

Тогда ПФЭ типа 22 содержит четыре опыта, что позволяет не только оценить коэффициенты a 0 , a1 , a 2 , но и проверить их значимость. Проверка адекватности модели может быть осуществлена, если известна дисперсия воспроизводимости опыта. Для оценки коэффициентов линейной модели общего вида (1.3.1) при N опытах укажем основные свойства матрицы ПФЭ: 1) симметричность относительно центра эксперимента: сумма элементов столбца равна нулю: N

∑ xij

= 0; i = 0 ÷ n

(1.4.2)

j =1

2) нормированность: сумма квадратов опытов: 15

элементов столбца равна числу

N

∑ ( xij ) 2 = N , i = 0 ÷ n ;

(1.4.3)

j =1

3) Ортогональность - скалярное произведение

любых

двух векторов-

столбцов матрицы равно нулю: N

∑ xij * x kj

= 0, i ≠ k , i = 0 ÷ n,

k = 0 ÷ n.

(1.4.4)

j =1

Регрессионный анализ в условиях ПФЭ ввиду указанных свойств не представляет

труда,

ибо

система

нормальных

уравнений

МНК

(1.5)

распадается на независимые уравнения, и оценки коэффициентов задаются простыми формулами: aˆ i =

1 N ∑ xij y j , i = 0 ÷ n N j =1

(1.4.5)

Таблица 2. − Полный факторный эксперимент для трех переменных, варьируемых на двух уровнях (планирование типа 23, линейная модель) Матрица планирования Х x0

x1

x2

x3

+1

+1

+1

+1

+1

-1

+1

+1

+1

+1

-1

+1

+1

-1

-1

+1

+1

+1

+1

-1

+1

-1

+1

-1

+1

+1

-1

-1

+1

-1

-1

-1

Вследствие исходных предположений регрессионного анализа о независимости измерений и их равноточности получаем, что для дисперсий 16

оценок параметров модели выполняется и равенство их дисперсий и, следовательно, среднеквадратических ошибок

σ 2 (aˆ i ) = σ 2 ( y ) / N ,

σ (aˆ i ) = σ ( y ) / N

(1.4.6)

Ввиду диагональности ковариационной матрицы оценки коэффициентов

aˆ i регрессии оказываются некоррелированными:

cov(aˆ i , aˆ k ) = 0, i ≠ k ;

i, k = 0 ÷ n

(1.4.7)

Это означает, что коэффициенты регрессии определяются независимо друг от друга (в обычном предположении нормальности ошибок), и, следовательно, их значения соответствуют относительному вкладу в выходное значение y каждого из слагаемых модели (1.3.1). В результате расчетов оценок линейная по параметрам модель выходной переменной y получает вид:

yˆ = aˆ 0 + aˆ1 x1 + ... + aˆ n x n ,

(1.4.8)

где знак ^ означает оценку величины. Для любого набора независимых факторов оценка значения выходной переменной в соответствие с моделью (1.4.8) находится путем подстановки их величин факторов в построенное уравнение регрессии. 1.4.1. Критерий Стьюдента для проверки значимости параметров модели Оценку значимости (проверку нуль − гипотезы о том, что aˆ i = 0) выполняют с помощью t-критерия Стьюдента, который с учетом свойств ортогональности плана, независимости оценок параметров (1.4.5) и (1.4.6) и их равной точности (1.4.7) принимает вид: | ti | =

aˆ i (σ ( y ) / N )

17

, i=0÷n

(1.4.9)

Если для выбранного уровня значимости α

(или доверительного уровня

(1-α)) вычисленная по формуле (1.4.9) величина | t i | больше соответствующего табличного значения | tα , f |- критического значения t-статистики,

нуль-

гипотезу отвергают. Считается, что параметр aˆ i значимо отличен от нуля с

уровнем доверия не менее (1-α). Переменная xi рассматривается как значимо влияющая на величину отклика – выходную переменную y . В модели (1.4.8) правомочно присутствие члена aˆ i xi . В противном случае, фактор xi признаётся как незначимо влияющий на выходную переменную и из дальнейшего анализа и модели он может быть исключён. Исключение фактора из модели в силу ортогональности матрицы эксперимента не влечет изменения в величинах оценок оставшихся в модели коэффициентов. Это свойство в момент разработки теории планирования экспериментов было весьма существенным, так как позволяло экономить вычислительные ресурсы. Одно из главных достоинств многофакторного эксперимента состоит в том, что в опытах одновременно варьируются все переменные. Поэтому каждый коэффициент регрессии определяют по результатам всех N экспериментов, и дисперсия в оценке коэффициентов регрессии в N раз меньше дисперсии ошибки единичного опыта.

Так, что

в полном

факторном эксперименте дисперсия оценки коэффициента модели в N раз меньше, чем в традиционном однофакторном. При использовании критерия Стьюдента следует различать случаи, каким образом получена оценка

ошибки измерения σ ( y ) , которая, как правило,

априорно неизвестна. В регрессионном анализе статистических пакетов Statistica и Excel

в

качестве оценки σ 2 ( y ) используется остаточная сумма квадратов (1.1.9) и формула (1.1.10): 2 S ост sˆ { y} = N − ( n + 1) 2

18

с числом степеней свободы для статистики Стьюдента f= N − (n + 1) , где N число

экспериментальных точек,

n +1 – общее число параметров регрессионной

модели. В случае наличия оценки ошибки измерения σ ( y ) по параллельным наблюдениям (1.1.12) m

sˆ ( y ) =

∑(y k =1

k

− y) 2 ( m − 1)

,

число степеней свободы для статистики Стьюдента равно f=m-1. Эти различия в числе степеней свободы необходимо учитывать при определении критических значения t-статистики. Доверительные интервалы для параметров регрессии помощью критического значения t- статистики

строятся

с

для уровня значимости α,

соответствующему (1-α) доверительному уровню. Так, используя (1.4.9), (1-α) доверительный интервал для параметра aˆ i записывается в виде: ( a) i − ∆ a (α , f ), a) i + ∆ a (α , f ) )=( a) i − tα , f s(a) i ), a) i + tα , f s(a) i ) )= i

i

=( a) i − tα , f s( y ) / N , a) i + tα , f s( y ) / N ) ,

(1.4.10)

а полуширина доверительного интервала определяется формулой: ∆ ai (α , f ) = tα , f s( y ) / N

(1.4.11)

Второй способ проверки значимости параметра aˆ i по t-критерию состоит

в сравнении его величины с полушириной доверительного интервала, этот способ является модификацией критерия Стьюдента в проверке нулевой гипотезы aˆ i =0.. Если | aˆ i |≤ ∆ a (α , f ) = tα , f s( y ) / N i

,

(1.4.12)

то принимается нулевая гипотеза, коэффициент aˆ i признается незначимым с

уровнем доверия (1-α) и из модели соответствующий член исключается. 1.4.2. Проверка адекватности модели по критерию Фишера

19

До сих пор полагали, что зависимость выхода y описывается линейной моделью, то есть уравнением регрессии вида (1.4.8). Справедливость этого предположения помощью

критерия

Фишера,

если

можно известна

проверить оценка

с

дисперсии

s 2 ( y ) величины y , полученная в серии из m измерений y в одинаковых условиях, и число степеней свободы для оценки дисперсии ошибки опыта равно m-1 (1.1.11).

Критерий Фишера производит сравнение отклонений модельных значений yˆ j , вычисленных с помощью уравнения (1.4.8), от реально наблюденных значений y j в j-той экспериментальной точке. Как уже упоминалось выше, 2 : строится так называемая остаточная сумма квадратов S ост

N 2 S ост = ∑ ( yˆ j − y jэ ) 2 , j =1 которая представляет

(1.4.13)

собой сумму квадратов уклонений экспериментальных

точек от теоретической кривой регрессии. Если модель верна, то величина 2 S ост s { y} = N − (n + 1)

ˆ2

(1.4.14)

является оценкой для дисперсии величины y . Таким образом, имеются две оценки одной и той же величины, полученные разными путями.

В случае

2 не адекватности модели оценка sˆ 2 { y} по остаточной сумме квадратов S ост

может слишком сильно отличаться от оценки s 2 ( y ) (1.1.11)., непосредственно полученной в эксперименте. Сравнивая остаточную дисперсию sˆ 2 { y} с оценкой дисперсии s 2 ( y ) , можно проверить гипотезу об адекватности модели по F-критерию Фишера: Fˆ = sˆ 2 { y} / s 2 ( y ) 20

(1.4.15)

После вычисления (1.4.15), проводится сравнение полученной величины Fˆ с табличными значениями, устанавливаемыми в зависимости от выбранного уровня значимости α и от степеней свободы числителя

m1= N − (n + 1) и

знаменателя m2= m-1. Можно также воспользоваться возможностями Excel (функция FРАСПОБР) При уровне значимости α и

указанных выше степенях

свободы m1 и m2 находится Fα ,m1,m2 - предельное максимальное значение для величины Fˆ . Если Fˆ < Fα ,m1,m2 ,

принимается гипотеза об адекватном

представлении экспериментальных данных линейной моделью. В противном случае, если

Fˆ ≥ Fα ,m1,m2 , гипотеза адекватности должна быть отвергнута,

поскольку это означает, что sˆ 2 { y} значительно превышает s 2 ( y ) .

Допустим, что модель нелинейна. Можно ли оценить количественно нелинейность с помощью ПФЭ? Когда факторы функционально связаны (эффект одного фактора зависит от уровня, на котором находится другой), говорят о наличии эффекта взаимодействия. Простейший его вид — произведение эффектов xi x k ,

i ≠ k.

Если рассматривать модель в виде неполного полинома второго порядка

y = ∑ ai xi + ∑ aik xi x k ,

i≠k

(1.4.16)

(без квадратичных слагаемых), то ПФЭ позволяет оценить эффекты взаимодействия, описываемые членом модели aik xi x k . Для этого нужно расширить матрицу планирования, добавив столбцы, соответствующие всем учитываемым

произведениям

xi x k .

Эти

столбцы

получают

путем

поэлементного перемножения столбцов, соответствующих xi и x k .

После этого с новым вектор - столбцом можно обращаться как с вектор столбцом нового фактора. Полный факторный эксперимент позволяет провести оценивание

параметров моделей, содержащих и тройные взаимодействия

факторов, что иллюстрирует таблица 3. При этом для всех оценок параметров 21

моделей с парными и тройными взаимодействиями сохраняются свойства независимости оценок и ортогональности столбцов матрицы плана. Таблица 3. − Полный факторный эксперимент для трех переменных, варьируемых на двух уровнях (планирование типа 23, модель с парными и тройным взаимодействиями) Матрица планирования Х

x0

x1

x2

x3

x1 x 2

x1 x3

x 2 x3

x1 x 2 x 3

+1

+1

+1

+1

+1

+1

+1

+1

+1

-1

+1

+1

-1

-1

+1

-1

+1

+1

-1

+1

-1

+1

-1

-1

+1

-1

-1

+1

+1

-1

-1

+1

+1

+1

+1

-1

+1

-1

-1

-1

+1

-1

+1

-1

-1

+1

-1

+1

+1

+1

-1

-1

-1

-1

+1

+1

+1

-1

-1

-1

+1

+1

+1

-1

Пользуясь таким планированием, можно вычислить коэффициенты регрессии неполного квадратного уравнения (1.4.16), используя описанные выше формулы. 1.5 Дробные реплики полного факторного эксперимента Если в линейной модели n+1 коэффициентов, а ПФЭ содержит 2n опытов, то с ростом n число опытов все более превосходит число оцениваемых коэффициентов. Например, при n=7 число опытов равно 27= 128,

то есть из 128 опытов лишь 8 нужны для однозначной оценки

коэффициентов существенно

линейной сократить

модели. число

За

опытов,

счет

этой

причем

избыточности так,

чтобы

можно матрица

планирования не утратила оптимальных свойств. План эксперимента, сохраняющий основные свойства ПФЭ, но включающий значительно меньше 22

опытов н а з ы в а ю т д р о б н ы м ф а к т о р н ы м э к с п е р и м е н т о м [11,13,14]. За

сокращение

числа

опытов

приходится

«расплачиваться»:

оценки

коэффициентов некоторых линейных и нелинейных членов в уравнении регрессии оказываются смешанными (подробнее об этом см. в литературе [11,13,14]). Поясним идею дробных реплик на примере. Допустим, что нам нужно получить

линейное

приближение

некоторого

небольшого

участка

поверхности отклика при трех независимых переменных. Для решения этой задачи можно ограничиться четырьмя опытами, если в планировании для полного факторного эксперимента типа 22 произведение x1 x 2 приравнять третьему

фактору

x3 .

Будет

получена

матрица

планирования,

представленная в таблице 4 и называемая дробной репликой ПФЭ. Элементы этой матрицы в точности равны элементам матрицы, представленной в таблице 1. Но опыты здесь будут уже ставиться с включением третьего независимого переменного x3 .

Подобная замена

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

взаимодействия

отсутствуют,

и

соответствующие

коэффициенты

регрессии aik равны нулю.

Таблица 4. − Дробный факторный эксперимент для трех переменных, варьируемых на двух уровнях (планирование типа 23-1, линейная модель) Матрица планирования Х x0

x1

x2

x3 = x1 x 2

+1

+1

+1

+1

+1

-1

+1

-1

+1

+1

-1

-1

+1

-1

-1

+1

С новым вектор - столбцом обращаются как с вектор - столбцом нового фактора. Пользуясь таким планированием, можно оценить свободный член a 0 23

и три коэффициента регрессии при линейных членах. Если коэффициенты регрессии aik при парных произведениях не строго равны нулю, то найденные коэффициенты регрессии будут оценками для совместных эффектов

a1' → a1 + a 23 a 2' → a 2 + a13 a 3' → a3 + a12 Эти эффекты не могут быть раздельно оценены в планировании, состоящем всего из четырёх опытов, так как здесь неразличимы столбцы для линейных членов и парных произведений. Для решения вопросов о разрешающей способности дробных реплик, то есть возможности раздельного оценивания параметров модели, вводятся понятия генерирующих соотношений, при помощи которых задаются дробные реплики. Например, планирование типа 23-1

может быть представлено двумя разными

полурепликами, каждая из которых задается одним из следующих генерирующих соотношений:

x3 = x1 x 2

и

x3 = − x1 x 2

(1.5.1)

Будем обозначать элементы первого столбца матрицы планирования символов I (все они всегда равны единице). Найдём соотношения, определяющие элементы первого столбца для каждой из полуреплик. Умножая правые и левые части написанных выше соотношений (5.1) на x3 , получим: x32 = I = x1 x 2 x3

и

x32 = I = − x1 x 2 x3 ,

(1.5.2)

так как всегда xi2 = I . Будем называть определяющим контрастом соотношения, задающие элементы первого столбца.

Две рассмотренные полуреплики (1.5.1) имеют

определяющие контрасты

I = x1 x 2 x3

и

I = − x1 x 2 x3

24

,

(1.5.3)

Зная определяющие контрасты, легко найти соотношения, задающие все совместные оценки. Для этого надо последовательно помножить независимые переменные на определяющий контраст. В нашем случае совместные оценки будут задаваться соотношениями x1 = x 2 x3

и

x1 = − x 2 x3

x 2 = x1 x3

и

x 2 = − x1 x3

x3 = x1 x 2

и

x3 = − x1 x 2

(1.5.4)

Это значит, что коэффициенты регрессии двух полуреплик будут оценками a1' → a1 + a 23

a1' → a1 − a 23

a 2' → a 2 + a13

a 2' → a 2 − a13

a 3' → a3 + a12

a3' → a3 − a12

Две полуреплики в планировании типа 24-1

могут быть заданы

генерирующими соотношениями

x 4 = x1 x 2 x3

x 4 = − x1 x 2 x3

и

,

(1.5.5)

Определяющими контрастами для этих реплик будут соотношения I = x1 x 2 x3 x 4

и

I = − x1 x 2 x3 x 4 ,

(1.5.6)

Совместные оценки будут определяться соотношениями x1 = x 2 x3 x 4

и

x1 = − x 2 x3 x 4

x 2 = x1 x3 x 4

и

x 2 = − x1 x3 x 4

x3 = x1 x 2 x 4

и

x3 = − x1 x 2 x 4

x 4 = x1 x 2 x3

и

x 4 = − x1 x 2 x3

x1 x 2 = x3 x 4

и

x1 x 2 = − x3 x 4

x1 x3 = x 2 x 4

и

x1 x3 = − x 2 x 4

x1 x 4 = x 2 x3

и

x1 x 4 = − x 2 x3

Таблица 5. − Матрица плана дробной реплики 24-1 факторного эксперимента для четырех факторов с определяющим контрастом 25

(1.5.7)

I = x1 x 2 x3 x 4 Матрица плана X x0

x1

x2

x3

+1

+1

+1

+1

+1

+1

-1

+1

+1

-1

+1

+1

-1

+1

-1

+1

-1

-1

+1

+1

+1

+1

+1

-1

-1

+1

-1

+1

-1

+1

+1

+1

-1

-1

+1

+1

-1

-1

-1

-1

x 4 = x1 x 2 x 3

Если априори можно принять, что все тройные взаимодействия равны нулю, то, пользуясь любой из этих полуреплик, можно получить раздельные оценки для всех четырёх линейных эффектов и три совместные оценки для парных произведений. В отличие от планирования 23-1 здесь линейные эффекты определяются отдельно от парных взаимодействий. В

случае,

когда

линейная

аппроксимация

поверхности

отклика

оказывается неадекватна, переходят к описанию поверхности с помощью моделей второго или более высокого порядка. Для полинома второго порядка модель имеет вид: y = a0 +

n

∑ a i xi

i =1

+

n

∑ aik xi x k

+ ε,

(1.5.8)

l ,k =1

Она содержит 1+n+n(n+1)/2 параметров a i и a ik

и, следовательно,

для их оценки требуется большее, чем для линейного приближения, число опытов. В отличие от ПФЭ, различные планы второго порядка обладают неодинаковыми оптимальными свойствами, так что планы, оптимальные по одним критериям, оказываются неудовлетворительными по другим. 26

1.6 Примеры применения факторных планов и дробных реплик

Рассмотрим примеры применения методов планирования экспериментов в реальных практических задачах. 1.6.1 Планирование эксперимента и математическое моделирование при поиске оптимального состава консервационного масла. Полный факторный план

Для защиты металлических изделий от атмосферной коррозии при транспортировке и хранении применяются

консервационные масла. Обычно

основу консервационных масел составляют минеральные масла, полученные в результате переработки нефти. Эти масла, особенно высоковязкие, уже сами по себе

предотвращают коррозионное повреждение

металлической поверхности.

Однако длительность защитного механизма не велика и для усиления защитного эффекта в масляную основу вводят так называемые пакеты противокоррозионных присадок. В зависимости от условий выбираются

как

основы

масляных

транспортировки и хранения изделий фракций,

так

и

составы

пакетов

противокоррозионных присадок. Для создания консервационного масла, защищающего изделия трубного завода в процессе транспортировки на открытом воздухе, была выбрана масляная фракция с величиной кинематической вязкости 18 сСт. Испытания показали, что защитная

способность

только

основы

не

удовлетворяет

предъявляемым

требованиям. Согласно техническим условиям на эксплуатацию консервационное масло должно выдерживать испытания в камере соляного тумана не менее

2

суток, в то время как консервация только основой не выдерживает и одних суток. (При оценке защитной эффективности в лабораторных условиях используются соответствующие методики ГОСТ 9.054-75). Другая характеристика работоспособности масла – его стабильность, то есть стабильное удерживание в объеме масла введённых присадок. Согласно

27

техническим требованиям осадок в консервационном масле не должен превышать 3 % от объёма через месяц после выпуска продукции. Пакет присадок составляли две противокоррозионные присадки - А и В, хорошо зарекомендовавшие себя в других масляных композициях Таким образом рецептура разрабатываемого консервационного масла состояла из масляной основы и двух компонентов А и В, выполняющих функции защитных противокоррозионных присадок. Основная

цель

разработки

формулировалась

следующим

образом:

определить концентрации компонентов, при которых свойства или показатели качества консервационного масла удовлетворяют предъявляемым требованиям к основным рабочим характеристикам – защитным свойствам и стабильности масла. Роль

независимых переменных или факторов

концентрации двух присадок, а зависимыми,

x1 и

x2

здесь играют

выходными переменными или

откликами выступают технические характеристики масла, такие как количество циклов до появления коррозии при испытаниях в камере солевого тумана y1, количество осадка, выпадающее после 30 суток с момента производства y2, кинематическая вязкость масла y3 . По аналогии с похожими композициями присадок для диапазонов концентраций выбраны интервалы значений: (% масс.): 5–15 для присадки

А, и

1–2 для присадки В. Относительная узость диапазонов позволила предположить, что

поведение показателей качества масла

можно будет описать линейной

моделью их зависимости от концентраций присадок. Для линейной модели достаточно применить полный факторный план на двух уровнях для двух факторов с числом экспериментов (образцов) N=22=4, представленный в таблице 6.1. Варьирование концентрациями

присадок

А и В проводилось в

соответствии с планом таблицы 6.1, где план записан в кодированных переменных x. Таблица 6.1 − Полный факторный эксперимент

28

для двух независимых переменных на двух уровнях (планирование 22) План эксперимента x0

x1

x2

1

+1

-1

-1

2

+1

+1

-1

3

+1

-1

+1

4

+1

+1

+1

№ опыта

Диапазоны варьирования переменных, как указано выше, представляли собой интервалы (% масс.): для А С1,мин=5≤C1≤15=С1,макс, для В С2,мин=1≤C2≤2=С2,макс, а центральная точка имела координаты С1,0=(С1,мин+С1,макс)/2 =10; С2,0=(С2,мин+С2,макс)/2 =1,5; Cвязь кодированных переменных с натуральными переменными С описывается равенствами: х1=2(С1-С1,0)/(С1,макс-С1,мин), х2=2(С2-С2,0)/(С2,макс-С2,мин) или х1= (С1-10)/5,

х2=(С2-1.5)/ 0.5 .

По плану таблицы 6.1 составлялись композиции присадок, причём реализовывались сочетания концентраций со значениями концов интервалов, которым соответствуют значения кодированных переменных +1 и -1. Составы, а также результаты измерений характеристик качества образцов консервационного масла, составленных по описанному выше плану, представлены в Таблице 6.2. Данные Таблицы 6.2 позволяют

получить в аналитическом виде

зависимости различных свойств от концентраций компонентов и определить на основе этих зависимостей оптимальные сочетания концентраций. 29

− Показатели качества образцов консервационного масла с

Таблица 6.2 присадками N экспери-

Концентрации присадок C,

мента

% объема Присадка А,

Присадка В,

С1

С2

Технические характеристики Защита в соляном тумане, цикл y1 ,

Стабильность, % осадка, y2

y3

Вязкость при +40, сСт,

1.

5

1

2

10

20.39

2.

15

1

3

0

23.79

3.

5

2

1

25

20.39

4.

15

2

5

3

24.84

При моделировании зависимостей на основе полных факторных планов или дробных реплик

вычисления оценок параметров модели, проведение анализа

значимости коэффициентов и

проверка адекватности модели проводится в

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

запись x

модели

в

кодированных

переменных

использует

для независимых переменных, а для натуральных значений

переменных - обозначения C. Линейная регрессионная модель зависимости характеристики y

от

концентраций в кодированных переменных имеет вид: y = a0 + a1 x1+ a2 x2, где

(1.6.1)

y -значение моделируемой характеристики, x1, x2,,– концентрации

компонентов, а0, a1, a2–неизвестные параметры, численные оценки которых находятся по экспериментальным данным методом наименьших квадратов. Оценки для ортогональных планов находятся по формулам (1.4.5) 1 N aˆ i = ∑ xij y j , i = 0 ÷ n N j =1 30

,

j-номер эксперимента (здесь меняется от 1 до 4). Статистические

параметров модели а0, a1, a2 связаны с

свойства

дисперсией ошибки измерения выходной характеристики σ 2 ( y ) (см. 1.4.6):

σ 2 (aˆi ) = σ 2 ( y) / N

,

σ (aˆi ) = σ ( y) / N

В тех случаях, когда дисперсия σ 2 ( y ) ошибки наблюдения неизвестна,

2 . для ее оценки используется остаточная сумма квадратов Sост N 2 Остаточная сумма (1.4.13) S ост = ∑ ( yˆ j − y jэ )2 j =1

представляет

собой сумму квадратов уклонений экспериментальных точек от уравнения регрессии, она связана с дисперсией величины y соотношением (1.4.14):. 2 Sост s {y} = ≈ σ 2 ( y) N − (n + 1)

ˆ2

Оценка значимости (проверка нуль − гипотезы о том, что параметр незначим или aˆ i = 0) выполняется с помощью t-критерия Стьюдента, который с учетом свойств ортогональности плана, независимости оценок параметров и их равной точности (1.4.6) принимает вид (1.4.9) | ti | =

aˆ i (σ ( y ) / N )

≈

aˆi ( sˆ( y) / N )

i=0÷n

Ниже приведены модели показателей качества масла с присадками, а также статистические характеристики их параметров. Количество циклов в соляном тумане без коррозионного поражения:

y1 = 2.75 + 1.25 x1+ 0.25 x2 при этом sˆ( y) =1.5,

t (aˆ ) = 1.66 , 1

,

(1.6.2)

t (aˆ ) =0.33 и число степеней свободы для 2

критерия Стьюдента в этом случае равно 1 (N-(n+1)=4-3=1). Согласно полученным критическим точкам (функция СТЬЮДРАСПОБР в Excel) оба параметра должны быть признаны незначимыми на уровне значимости 31

α=0.05,

если

за оценку

дисперсии

остаточную

ошибку. В

действительности 5-ю параллельными опытами установлено, что

оценка для

дисперсии составляет s 2 ( y ) =0.16

выбирать

или s(y)=0.4. В этом случае t (aˆ ) = 6.25 , 1

t (aˆ ) =1.25 , так что незначимым для числа степеней свободы 4 оказывается 2

лишь параметр, стоящий при x2. Исключение соответствующего члена из модели приводит к зависимости: y1 = 2.75 + 1.25 x1

(1.6.3)

Величина осадка (%) после 30 суток описывается моделью: y2= 9.5 - 8 x1+4.5 x2, при этом sˆ( y) =6.0,

t (aˆ ) = -2.75 , 1

(1.6.4)

t (aˆ ) =1.5 , число степеней свободы для 2

критерия Стьюдента равно 1. В соответствии со статистическими таблицами (или воспользовавшись как и в предыдущем случае возможностями Excel) оба параметра должны быть признаны незначимыми на уровне значимости α=0.05, если за оценку дисперсии выбирать остаточную ошибку sˆ( y) . В действительности 4-мя параллельными опытами установлено, что

оценка для дисперсии ошибки измерения осадка

составляет s 2 ( y ) =3.6 или s(y)=1.9. В этом случае t (aˆ ) = 9.41 , 1

t (aˆ ) =5.29 и оба 2

коэффициента оказываются значимыми. Проверка адекватности модели по критерию Фишера следующему заключению:

(1.4.15) приводит к

для числа степеней свободы числителя m1=

N − (n + 1) =4-3=1 и знаменателя m2=4-1=3 и при вычисленном статистики

Фишера

Fˆ = sˆ 2 { y} / s 2 ( y ) =36/3.6=10

имеем

⎯

значении для

уровня

значимости α=0.05 критическая величина статистики Фишера равна 10.1

(функция FРАСПОБР с параметрами 0.05, 1,3) , то есть больше, чем полученное по экспериментальным данным значение Fˆ .

32

Следовательно,

можно считать, что линейная модель

(1.6.4) адекватно

описывает зависимость величины осадка от концентраций компонентов. Величина кинематической вязкости описывается уравнением y3= 22.35+1.96 x1+0.26 x2,, при этом sˆ( y) =0.52, t (aˆ ) = 7.47 , 1

(1.6.5)

t (aˆ ) =1. По критерию Стьюдента значимым 2

на уровне α=0.1 оказывается лишь коэффициент при

первой переменной x1,

поэтому модель для кинематической вязкости описывается уравнением y3= 22.35+1.96 x1

(1.6.6)

Подстановка в полученные модели натуральных значений концентраций переменных приводит к зависимостям: Количество циклов в соляном тумане без коррозионного поражения

y1 = 2.75 + 1.25 x1

=2.75+1.25 (С1-10)/5=0.25 +0.25 С1 ,

(1.6.7)

множественный коэффициент корреляции 0.845. Величина осадка (%) после 30 суток y2= 9.5 - 8 x1+4.5 x2,=9.5-8(С1-10)/5 +4.5 (С2-1.5)/ 0.5=12 -1.6 С1 +9.0 С2, (1.6.8) множественный коэффициент корреляции 0.95. Величина кинематической вязкости, сСт, y3= 22.35+1.96 x1 =22.35+1.96 (С1-10)/5=18.35+0.4 С1

(1.6.9)

множественный коэффициент корреляции 0.98. Из представленных моделей можно сделать следующие заключения. Максимальный

уровень

защитных

свойств,

а

также

стабильности

(наименьшие значения осадка (%)) достигаются при максимальном содержании в композиции присадки А, то есть когда концентрация А достигает 15%. В то же время

с

увеличением

концентрации

присадки

В

начинает

ухудшаться

стабильность масла, так при росте содержания в композиции компонента В увеличивается процент осадка. Поэтому наилучшей концентрацией для присадки В является её минимальное содержание в композиции, то есть 1 %. Таким образом, наилучшей композицией присадок для консервационного масла является композиция под номером 2 таблицы 6.2: присадка А -15 %, присадка В -1 %. 33

1.6.2 Разработка композиции присадок к консервационному маслу. Дробный факторный эксперимент

Теперь рассмотрим процесс применения дробного факторного эксперимента в той же задаче разработки консервационного масла, но с иными тремя противокоррозионными присадками D, C, F. Рецептура разрабатываемого консервационного масла состоит из основы – минерального масла с вязкостью 30 сСт- и компонентов D, C и F. Цель разработки - определить оптимальные концентрации компонентов, при которых свойства или показатели качества консервационного масла удовлетворяют

предъявляемым

требованиям

к

основным

рабочим

характеристикам. Предварительными

экспериментами

установлено,

что

диапазоны

концентраций компонентов представляются следующими интервалами (% масс.): 1-4 для D, 0,8-3 для C , 3-5 для F. Для получения в аналитическом виде зависимости

различных свойств

консервационного масла от концентраций компонентов и определения на основе этих зависимостей оптимальных композиций применяется метод планирования эксперимента и регрессионного анализа. Поскольку диапазон варьирования концентраций компонентов достаточно узкий, справедливо предположить, что в этом диапазоне зависимости могут быть описаны линейными регрессионными моделями, в которых независимыми переменными выступают концентрации компонентов: С1- концентрация D (% масс.), С2- концентрация C (% масс.), С3- концентрация (% масс.). Линейная

регрессионная

модель

изучаемой

характеристики

консервационного масла y от концентраций компонентов имеет вид: y = a0 + a1 C1+ a2 C2+ a3 C3 , где

(1.6.10)

y значение моделируемой характеристики, C1, C2, C3 – концентрации

компонентов – присадок соответственно D, C, F; а0, a1, a2 , a3 –неизвестные

34

параметры, численные оценки которых находятся из экспериментальных данных методом наименьших квадратов. В данной задаче мы имеем дело с 3-мя независимыми переменными – факторами- концентрациями присадок, а число неизвестных параметров модели (1.6.10) равно четырем. Поэтому число экспериментов, необходимое для нахождения оценок параметров, не может быть менее 4-х. линейной формы зависимости выходной

В условиях предположения

характеристики для оценки 4-х неизвестных параметров модели (1.6.10) можно применить наиболее экономически выгодный план – дробную реплику полного факторного эксперимента на 2-х уровнях для 3-х факторов с числом экспериментов N=23-1=4, определяющим контрастом I=-х1х2х3 и генерирующим соотношением х3=-х1х2 [11,13,14]. Для того, чтобы можно было осуществить проверку значимости параметров по критерию Стьюдента композиция

с

к дробной реплике добавлена центральная точка

координатами,

соответствующими

серединам

варьирования переменных. В данной задаче добавление

-

диапазона

центральной точки

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

интервалы

(%

масс.):

С1,мин=1≤C1≤4=С1,макс

для

присадки

D,

С2,мин=0,8≤C2≤3=С2,макс, для присадки C, С3,мин=3≤C3≤5=С3,макс для присадки F, а центральная точка имеет координаты С1,0=(С1,мин+С1,макс)/2 =2,5; С2,0=(С2,мин+С2,макс)/2 =1,9; С3,0=(С3,мин+С3,макс)/2 =4. План эксперимента 23-1

с генерирующим соотношением х3=-х1х2

и

центральной точкой представлен в кодированных переменных в Таблице 6.3. Связь кодированных и натуральных переменных описывается равенствами: х1=2(С1-С1,0)/(С1,макс-С1,мин),

х2=2(С2-С2,0)/(С2,макс-С2,мин), 35

х3=2(С3-С3,0)/(С3,макс-С3,мин). Таблица 6.3− План эксперимента- дробная реплика. 23-1 с определяющим контрастом I=-х1х2х3 номер эксперимента 1 2 3 4 5

х1 -1 1 -1 1 0

Кодированные переменные х2 х3 -1 -1 -1 1 1 1 1 -1 0 0

Условия эксперимента в натуральных переменных и результаты измерений характеристик качества масла представлены в Таблице 6.4 В силу особенностей выбранной дробной реплики для оценки остаточной дисперсии здесь

остаётся лишь одна степень свободы, и при проверке

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

показатели качества. Поэтому в данной задаче

целесообразно воспользоваться достаточно точными оценками ошибок измерения характеристик, полученными в результате 5-ти параллельных наблюдений в одних и тех же условиях

s 2 ( y ) =0,16 ( s ( y ) =0.4) для соляного тумана,

s 2 ( y ) =0,09 ( s ( y ) =0.3) для осадка, s 2 ( y ) = 0,09 ( s ( y ) =0.3) для вязкости. Как уже упоминалось в примере 6.1, моделирование зависимостей на основе полных факторных планов и дробных реплик переменных в силу

простоты расчетов и

производится

в кодированных

хороших статистических свойств

оценок. Запись модели в кодированных переменных использует обозначения x для независимых переменных, а для натуральных значений переменных будет использоваться обозначение C.

36

Таблица 6.4 − Составы образцов и показатели качества образцов консервационного масла с присадками N опыта

Технические характеристики масла

Составы образцов Концентрации присадок, % объема C1

C2

C3

1.

1

0,8

3

2

0

32,98

2.

4.

0,8

5

3

1

34,58

3.

1.

3

5

4

3

33,43

4.

4.

3

3

3

0,1

34,07

5.

2,5

1,9

4

2

0,3

34,04

Защита в соляном Стабильность, Вязкость при 40,сСт тумане, цикл % осадка, y3 y1 , y2

Линейная регрессионная модель характеристики консервационного масла y от концентраций 3-х компонентов в кодированных переменных имеет вид: y = a0 + a1 x1+ a2 x2,+ a3 x3

(1.6.11)

где y -значение моделируемой характеристики, x1, x2,, x3 – концентрации компонентов, а0, a1, a2 , a3 – неизвестные параметры, численные оценки которых находятся по экспериментальным данным таблицы 6.4. Оценки находятся по формулам (1.4.5): aˆ i =

1 N ∑ xij y j , i = 0 ÷ 4 , j-номер эксперимента N , изменяется здесь от 1 N j =1

до 5. Расчеты приводят к следующим моделям: Количество циклов в соляном тумане без коррозионного поражения:

y1 = 2.8 + 0.5 x2+ 0.5 x3 , при этом sˆ( y) =0.89,

(1.6.12)

а соответствующие остаточной ошибке статистики

Стьюдента равны t (aˆ2 ) = 1.12 , t (aˆ ) =1.12 при числе степеней свободы 1. Ясно, что 3 эта оценка грубая и надо воспользоваться более точной оценкой ошибки измерения s ( y ) =0.4.

37

aˆ

i В этом случае величины статистик Стьюдента | t i | = ( s( y ) / N )

с числом

опытов в плане N=5 и ошибкой 0.4 равны t (aˆ2 ) = 2.79 , t (aˆ ) =2.79 , значения 3

статистик

превосходят

критическую величину

2.77 для уровня значимости

α=0.05, что означает значимость коэффициентов регрессии. Величина осадка (%) после 30 суток описывается моделью: y2= 0.88 – 0.475 x1+0.525 x2 + 0.975 x3 , при этом sˆ( y) =0.648, t (aˆ ) = -1.465 , 1

(1.6.13)

t (aˆ ) =1.62 , t (aˆ ) =3.0 , а число степеней 3 2

свободы для критерия Стьюдента равно 1. В соответствии со статистическими таблицами все параметры должны быть признаны незначимыми на уровне α=0.05, если

за оценку дисперсии

выбирать остаточную ошибку sˆ( y) . В действительности использование более точных оценок дисперсий по серии наблюдений даёт значения статистик t (aˆ ) = 1 3.42 ,

t (aˆ ) =3.78 , t (aˆ ) =7.26, превосходящие по модулю критическое 3 2

значение для α=0.05, то есть все параметры значимы. Таким образом,

зависимость величины осадка от концентраций

компонентов описывается уравнением (1.6.13). Величина кинематической вязкости (сСт) описывается уравнением y3= 33.82+0.56 x1+0.24 x3,, при этом sˆ( y) =0.24, t (aˆ ) = 4.55 , 1

(1.6. 14 )

t (aˆ ) =1.95. 3

Оказывается, что при полученной остаточной ошибке и числе степеней свободы 1 значимым по критерию Стьюдента на уровне α=0.1 оказывается лишь коэффициент при первой переменной x1. Использование более точной оценки для ошибки измерения вязкости приводит к тому, что оба коэффициента модели оказываются статистически значимыми и в модели остаются. Модели (1.6.12) и (1.6.13) используются далее для выбора оптимальной композиции присадок.

38

Технологические требования, предъявляемые к разрабатываемому маслу, состоят в следующем: консервационное масло должно обеспечивать защиту в камере соляного тумана не менее 3-х циклов, масло должно быть стабильно, то есть осадок за 30 суток после изготовления не должен превышать 1% от объёма масла. Эти технологические требования можно записать в виде ограничений на зависимые переменные y: y1 ≥ 3;

y2 ≤ 1

(1.6.15)

Исходя из требований (1.6.15), необходимо найти оптимальные составы (рецептуры) для консервационного масла, то есть x1, x2 , x3

удовлетворяющие

требованиям: y1 = 2.8 + 0.5 x2+ 0.5 x3 ≥ 3;

(1.6.16)

y2= 0.88 – 0.475 x1+0.525 x2 + 0.975 x3

≤1

(1.6.17)

Анализ моделей показывает, что уменьшение осадка достигается только при увеличении содержания компонента x1, при этом согласно соотношению (6.16) изменение концентрации x1 в области экспериментирования не сказывается на потере защитной способности пакета. Поэтому концентрацию для x1 нужно выбрать максимально возможной, то есть x1=1. Подстановка

этого

значения

в

неравенство

(1.6.17)

приводит

к

соотношению: 0.4+0.525 x2 + 0.975 x3

≤1

(1.6.18)

Наибольший рост осадка вызывает увеличение концентрации x3 присадки F, поэтому разумно выбрать для неё уровень концентрации x3=0, так как при выборе для неё значения x3= -1 происходит потеря защитной способности ниже требуемых

3-х циклов. В предлагаемом варианте, если x3=0,. необходимый

уровень достигается при x2=1. Таким образом, оптимальной композицией оказывается сочетание x1=1, x2=1, x3=0 или в натуральном выражение 4 % присадки D, 3 % присадки C, 4% присадки F. Выводы по результатам приведённых примеров. 39

Как показывают результаты параметров

линейных моделей

статистического анализа значимости

для функций отклика, при использовании в

качестве оценок дисперсий опыта лишь остаточной ошибки аппроксимации, как правило, приводит к принятию нулевой гипотезы о не значимости параметров моделей. Этот эффект активного эксперимента связан с тем, что экономия на количестве необходимых экспериментов N, приводит к малому числу степеней свободы для остаточной ошибки, что в свою очередь приводит к росту критического порога для статистических критериев. В связи с этим рекомендуется предварительное

определение

ошибки

исследованиями.

40

эксперимента

дополнительными

1.7 Планирование эксперимента на диаграммах состав - свойство 1.7.1 Симплексная система координат

В настоящем разделе рассмотрено планирование эксперимента для систем, являющихся смесями q различных компонентов. Переменные xi (i=1,2, . . ., q) таких систем являются пропорциями, относительными содержаниями iтыx компонентов смеси и удовлетворяют условию: ∑ xi = 1, 1≤i ≤q

(xi ≥0)

(1.7.1)

Геометрическое место точек, удовлетворяющих условию нормирования суммы переменных (1.7.1), симплекс.

(Правильный

представляет собой симплекс

с

q

(q—1)-мерный

вершинами

в

правильный

(q—1)-мерном

пространстве, треугольник для q = 3, тетраэдр для q = 4 и т. д.). Каждой точке такого симплекса соответствует смесь определенного состава, и, наоборот, любой комбинации относительных содержаний q компонентов соответствует определенная точка симплекса. Так как в дальнейшем при планировании эксперимента и построении диаграмм состав—свойство [11,13,15]

придется оперировать факторным

пространством в виде симплексов, целесообразно определить координаты компонентов не в обычной системе координат, а в специальной — симплексной, в которой относительные содержания каждого компонента откладываются вдоль соответствующих сторон (граней) симплекса. Связь между двумя координатными системами — обычной декартовой и симплексной — для трехкомпонентных смесей показана на рисунке 1.7.1, а. Здесь изменению относительного содержания первого компонента x1 от 0 до 1 вдоль оси x1 (в долях от длины отрезка, равного единице) соответствует пропорциональное изменение координаты

x ′ ( x ′ = √2 x1 ) вдоль стороны ab 1 1

(длины √2) от точки а, где компонент x1 присутствует в пропорции 0, до точки 41

b, где содержание первого компонента равно 1, т. е. смесь состоит лишь из первого компонента. Легко заметить, что x1 долей (частей) отрезка оb, длиной 1, равно x ′ ( x ′ 1 1 = √2

x1) долям отрезка ab, длиной √2, т. е. относительное содержание

(пропорции) x1= x ′ . 1

Рисунок 1.7.1 − Переход к симплексной системе координат. Поэтому в дальнейшем, при рассмотрении симплексов штрих будем опускать, обозначая относительное содержание компонентов на его сторонах просто через xi. Аналогично, перемещению точки вдоль оси х2 от центра координат к точке с координатой x2=1 будет соответствовать пропорциональное перемещение точки вдоль стороны bc от точки b, где второй компонент отсутствует, к точке с, где имеется лишь второй компонент. Третий компонент на треугольной диаграмме (рисунок 1.7.1, b) откладывается вдоль стороны сa, начиная от точки с с нулевым содержанием данного компонента до точки а, где х 3 =1. Если в декартовой системе координат для определения уровня первого компонента, соответствующего какой-нибудь точке М смеси, достаточно взять отсекаемый на оси х 3 отрезок от, системе координат

то в симплексной (треугольной)

это аналогично проведению через точку М прямой, 42

параллельной стороне х3 (са), и взятию отрезка х1= am (рисунок 1.7.1, в и 1.7.1, г). Для определения второй координаты точки М проведем через эту точку прямую, параллельную стороне x1 и возьмем отсекаемый на стороне bc отрезок х2. Аналогично, пропорция третьего компонента определяется путем проведения через точку М прямой, параллельной стороне х2. (рисунок 1.7.1, г). В четырехкомпонентном случае для определения координаты х 1 какойнибудь точки трехмерного симплекса — правильного тетраэдра (рисунок 1.7.1, д) — необходимо провести через нее плоскость, параллельную двухмерной грани тетраэдра с ребром пропорций третьего компонента х3, и взять отсекаемый этой плоскостью на оси х 1 отрезок x( M ) . 1 Ввиду того, что переменные смеси, согласно (1.7.1), не являются независимыми, оценка коэффициентов модели (1.2) невозможна, так как матрица ( X T X ) оказывается вырожденной и не существует для неё обратной матрицы для уравнения (1.6). Исключение одной из смесевых переменных (обычно переменной, соответствующей основе смеси) и применение относительно q — 1 оставшихся независимых компонентов смеси планов, описанных в разделах 1.4 и 1.5, является одним из возможных вариантов решения задачи. Однако в тех случаях, когда выходной параметр системы зависит от всех q исследуемых смесевых переменных, следует пользоваться иными моделями, построенными с учетом (1.7.1), и строить для оценки коэффициентов этих полиномов оптимальные планы на области определения всех переменных. Впервые задача построения математической модели состав — свойство, включающей все компоненты системы, была решена Шеффе [15], который ввел каноническую (приведенную) форму полинома. Для простоты ограничимся рассмотрением обычных полиномов не выше второго порядка

y = a + ∑ ai xi + ∑ a xi x , 0 ik k 43

(1.7.2)

Полиномы Шеффе второго порядка получаются из полиномов обычного вида (1.7.2)

с учетом соотношения (1.7.1), например, для q=3 он имеет

следующее аналитическое выражение:

y =β x +β x +β x +β x x +β x x +β x x 11 2 2 3 3 12 1 2 13 1 3 23 2 3

(1.7.3)

Рассмотрим основные, наиболее широко применяемые планы на симплексе – симплекс - решетчатые планы. 1.7.2 Симплекс - решетчатые планы

Принципы построения

симплекс – решётчатых планов рассмотрим на

примерах полиномов первого и второго порядка. При построении этих планов используются представления полиномов в так называемой приведённой или канонической форме для областей переменных, заданных симплексом. Эти полиномы получаются из полиномов обычного вида подстановкой в них уравнения связи (1.7.1) Рассмотрим обычный полином первого порядка y = a + ∑ ai xi 0

(1.7.4)

При переходе к приведённой форме исключается свободный член a 0 умножением уравнения (1.7.1) на a0, подстановкой

его в (1.7.4)

и

последующей группировкой подобных членов: ∑ a0 xi =a0 , 1≤i ≤q

(1.7.5)

После обозначения как β коэффициентов при независимых переменных получаем: y=

∑

1≤ i ≤ q

(a + a x ) = ∑ βi xi i i 0 i

44

(1.7.6)

Для случая трёхкомпонентной системы линейное уравнение регрессии в приведённой форме принимает вид:

y =β x +β x +β x 11 2 2 3 3 Аналогичная подстановка (1.7.5) в

(1.7.7)

полином второй степени,

в

общем случае имеющий вид (1.7.2) при i=3, а также группировка членов и исключение квадратичных членов при помощи выражений: x2 = x − x x − x x , x2 = x − x x − x x , x2 = x − x x − x x , 2 12 2 3 3 1 1 12 13 2 3 13 2 3

приводит

(1.7.8)

к представлению уравнения регрессии на симплексе для трех

переменных в форме (1.7.3) после обозначения коэффициентов (параметров при переменных) через β:

y =β x +β x +β x +β x x +β x x +β x x 11 2 2 3 3 12 1 2 13 1 3 23 2 3

(1.7.9)

Для q переменных приведённый полином второго порядка имеет вид [11,13,14]: y= ∑ βi xi + β ij xi x j ∑ 1≤i≤q 1≤ i < j ≤ q

(1.7.10)

В случае трёхкомпонентной системы часто используется так называемая неполная кубическая модель

y = β x + β x + β x + β x x + β x x + β x x + β x x x (1.7.11) 11 2 2 3 3 12 1 2 13 1 3 23 2 3 123 1 2 3 Для оценки коэффициентов приведенного полинома Шеффе были предложены планы, обеспечивающие равномерный разброс экспериментальных точек по (q — 1)-мерному симплексу. Точками таких планов являются узлы {q,n}-симплексных (компонента)

решеток.

используется

В п+1

{q,n}-

решетке

для

каждого

равнорасположенных

фактора

уровней

в

интервале от 0 до 1 (x i = 0, 1/n, 2/n, . . ., 1) и берутся все возможные их

45

комбинации. Так как число таких комбинаций Cqn+ n −1 равно числу оцениваемых коэффициентов в приведенном полиноме степени п. Для n=1 набор точек представляет собой вершины симплекса, для n=2 набор точек состоит из вершин симплекса и середин граней как показано на рисунке 1.7.2

Рисунок 1.7.2 — Некоторые из {q ,n}-решеток: q= 3а — линейная {q= 3, n=1}; б — квадратичная {q= 3, n=2}; в — неполнокубическая; г — кубическая {q= 3, n=3); д — четвертой степени {q= 3, n=4}, с — восьмой степени {q= 3, n=8} q= 4ж — квадратичная {q=4, n=2}; 3 — кубическая {q=4, n=3} Примеры {q,n}- решетки приведены на рисунке 1.7.2. Ряд симплексрешетчатых планов низшего порядка (см. рисунок 1.7.2) входит составной частью в планы более высокого порядка. Так, например, квадратичная решетка может быть получена из линейных решеток добавлением срединных точек сторон, неполно-кубическая решетка — добавлением к квадратичной решетке всего одной точки в центре тяжести симплекса. Это свойство композиционности планов необходимо учитывать при изучении многокомпонентных систем, так как у исследователей на практике не всегда имеются достаточные априорные сведения о виде поверхности отклика каждой конкретной системы. Поэтому возможна и недооценка сложности поверхности отклика. Это может привести к тому, что найденный полином степени n, будет неадекватным и потребуется проведение дополнительных экспериментов для повышения степени полинома. 46

Рисунок 1.7.3 — Обозначение откликов в точках симплексных решеток Для оценки коэффициентов аппроксимирующего полинома степени п

(1.7.9) во всех точках плана, соответствующих узлам {q,n}-решетки, реализуются опыты и определяются отклики системы у Под

у

могут

подразумеваться

как

результаты

единичного

определения, так и средние значения нескольких определений. Удобно ввести специальные обозначения для этих откликов (рисунок 1.7.3). Отклик для смесей, содержащих только один ненулевой компонент (вершины симплекса, т. е. точки с координатами ( 01,..., 230 ; 1; 0, . . ., 0)), обозначим i −1

через y i , отклик для 1:1, бинарной смеси компонентов i и j — через уij (i