Пособие ТОКН

~ ~ ~	~ ~ ~ ~	~ ~ ~
1 = x1x2 x4	= x1x2 x3x5	= x3x4 x5 .

Система смешиваний при этом имеет довольно сложный вид и определяется умножением обобщающего определяющего контраста на ~x1, ~x2 , ~x3 и т.д.:

~	~ ~	~ ~ ~		~ ~ ~ ~		,
x1	= x2 x4	= x2 x3x5		= x1x3x4 x5
~	~ ~ ~ ~		~ ~ ~		~ ~	,
x3	= x1x2 x3x4		= x2 x1x5		= x4 x5
~	~ ~ ~ ~		~ ~ ~		~ ~	,
x5	= x1x2 x4 x5		= x1x2 x3		= x3x4

~x1~x5 = ~x2 ~x5~x4

~		~ ~	~ ~ ~		~ ~ ~ ~		,
x2	= x1x4		= x1x3x5 = x2 x3x4 x5
~		~ ~	~ ~ ~ ~ ~			~ ~	, (1.17)
x4	= x2 x1		= x1x2 x3x4 x5			= x3x5
~ ~		~ ~ ~		~ ~		~ ~ ~	,
x1x3		= x2 x3x4		= x2 x5	= x1x4 x5
~ ~		~ ~ ~
= x2 x3		= x1x4 x3.

В рассмотренной полуреплике все линейные эффекты смешаны с парными, и если коэффициенты β12 и β23 окажутся значимыми, то возникает сомне-

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

перимент –	четвертьреплика с		~ ~ ~	и
		определяющими контрастами 1 = −x1x2 x4
~ ~ ~ ~	. Обобщающий определяющий контраст тогда будет таким:
1 = x1x2 x3x5
	~ ~ ~	~ ~ ~ ~	~ ~ ~
	1 = −x1x2 x4 = x1x2 x3x5		= −x3x4 x5 .

Он дает следующую систему смешиваний для линейных членов:

b1 →β1	−β24 + β235 −β1345 ,		b2 →β2 −β14 +β135 −β2345 ,
b3 →β3	−β1234 +β125 −β45 ,		b4 →β4 −β12 + β12345		−β35 , (1.18)
		b5 →β5 −β1245 + β123 −β34 .
Усредняя (1.18) и систему смешиваний, полученную из (1.17), получим
b1 →β1 + β235 ,		b2 →β2 + β135 ,		b3 →β3	+β125 ,
	b4 →β4 + β12345 ,		b5 →β5 + β123 .

Система смешиваний имеет более простой вид, а главное, линейные эффекты уже не смешаны с парными.

31

На основе плана 23 можно построить еще два плана. Один из них – план 26−3 , другой – план с максимальным числом факторов – 7, при этом получается 1/16-реплика, или план 27−4 . Эти планы, однако, будут иметь очень сложную систему смешиваний, кроме того, план 27−4 будет насыщенным.

1.5.Учет влияния неконтролируемых внешних факторов

1.5.1.Полная рандомизация

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

температура окружающей среды, влажность, давление, радиация и прочее; их можно отнести к естественным (природным) факторам;

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

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

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

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

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

1.5.2. Блочная рандомизация В некоторых случаях оказывается выгодным проводить не полную рандоми-

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

32

Рассмотрим это на конкретном примере. Пусть при проведении эксперимента по

плану 22 в течение одного дня можно ставить только два опыта (или имеются две партии сырья, каждой из которых хватит не более чем на два опыта). В этом случае возникает опасность, что условия второго дня (или состав второй партии) несколько отличаются от первого и это скажется на результатах опыта. Желательно спланировать эксперимент так, чтобы возможные изменения в отклике не отразились на линейных оценках. Для этого проведем 1-й и 4-й опыты в первый день (с первой партией сырья), а 2-й и 3-й – во второй (со второй партией). Тем самым матрица разбивается на два блока, и рандомизацию теперь следует проводить внутри этих блоков. В этом и состоит ограничение на рандомизацию.

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

b0 = 14 (y1 + y2 + ε + y3 + ε + y4 )= 14 (y1 + y2 + y3 + y4 )+ 2ε , b1 = 14 (y1 + y2 + ε − y3 − ε − y4 )= 14 (y1 + y2 − y3 − y4 ), b2 = 14 (y1 − y2 − ε + y3 + ε − y4 )= 14 (y1 − y2 + y3 − y4 ), b12 = 14 (y1 − y2 − ε − y3 − ε + y4 )= 14 (y1 − y2 − y3 + y4 )− 2ε .

Таблица 1.10

Разбиение матрицы на 2 блока

	~		~		~		~ ~		y
	x0		x1		x2		x1x2
1	+		+		+		+		y1
2	+		+		–		–		y2 + ε
3	+		–		+		–		y3 + ε
4	+		–		–		+		y4

Можно заметить, что искажены оценки b0 и b12 . Это объясняется выбран-

ным способом разбиения матрицы на блоки. В принципе, можно рассматривать день (или сырье) как новый фактор x3 , и тогда рассмотренный план будет

представлять собой не что иное, как полуреплику 23−1 с определяющим кон-

= ~ ~ ~

трастом 1 x1x2 x3. Говорят также, что в этом случае матрица разбита на орто-

гональные блоки.

33

Если известно, что коэффициент β12 незначим, то можно оценить величи-

ну эффекта, вызванного неоднородностью (разными днями или разными партиями сырья), – так называемый межблоковый эффект. Как это сделать, ясно из выражения для оценки b12 .

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

план можно рассматривать как четвертьреплику 25−2 . Определяющие контрасты возьмем в виде 1 = ~x1~x2~x4 и 1 = ~x1~x2~x3~x5 , где ~x4 и ~x5 – факторы, «отвечающие» за день и сырье соответственно.

Предположим, что постановка опытов во второй день дает дополнительный вклад, равный γ, а вторая партия сырья – вклад δ. План и результаты эксперимента приведены в табл. 1.11. Согласно этому плану в первый день проводятся опыты 1, 2, 7, 8, во второй – 3, 4, 5, 6; первая партия сырья используется в опытах 1, 4, 6, 7, вторая – в опытах 2, 3, 5, 8.

Разбиение матрицы на 4 блока

Таблица 1.11

~

~

~

~

~ ~

~ ~ ~

y

x0

x1

x2

x3

x1x2

x1x2 x3

1

+

+

+

+

+

+

y1

2

+

+

+

–

+

–

y2+δ

3

+

+

–

+

–

–

y3+γ+δ

4

+

+

–

–

–

+

y4+γ

5

+

–

+

+

–

–

y5+γ+δ

6

+

–

+

–

–

+

y6+γ

7

+

–

–

+

+

+

y7

8

+

–

–

–

+

–

y8+δ

Для проверки подсчитаем линейные оценки:

b =	1	(y + y	2	+δ+ y					3	+ γ +δ+ y						4	+ γ − y					5		− γ −δ− y							6	− γ − y				7	− y	8	−δ)=
1	8	1
					1
				=			(y + y				2	+ y		3	+ y			4	− y		5		− y			6	− y		7	− y			8	),
				8				1
			b =			1		(y + y				2	− y		3	− y			4	+ y			5		+ y		6	− y		7		− y			),
				2		8			1																									8
			b =			1		(y − y				2	+ y		3	− y			4	+ y			5		− y		6	+ y		7		− y			).
				3		8			1																									8

Видно, что эффекты дня и сырья не влияют на линейные оценки, что и следовало ожидать.

34

1.6.Планирование эксперимента при поиске оптимальных условий

1.6.1.Общие понятия

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

Первый – метод Гаусса-Зейделя. Он заключается в поочередном варьировании всех переменных до достижения условного экстремума по каждой переменной; такая методика характерна для однофакторного эксперимента. Метод Гаусса-Зейделя не всегда решает поставленную задачу, а в случае успеха зачастую приходится делать много шагов.

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

Боксом и Уилсоном было предложено соединить градиентный метод поиска экстремума с ортогональным планированием. Суть этого соединения состоит

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

Метод Бокса-Уилсона напоминает лазание по горе с завязанными глазами

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

Координатные составляющие градиента можно определить, разделив оценки коэффициентов модели на интервалы варьирования соответствующих факторов. Покажем, что это действительно так для n = 2.

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

~	~	,
y = b0 + b1x1	+ b2 x2

к модели в обычных факторах, используя формулы перехода (1.4):

35

	x	− 0,5(xв + x		н)		x	2	− 0,5(xв		+ xн)
y =b + b	1	1	1		+ b				2	2	.
								0,5(xв	− xн)
0 1		0,5(xв − xн)			2
		1	1					2		2

Найдем составляющие градиента:

	∂y ∂xi =	bi	,	i = 1, 2.	(1.19)
		0,5(xв
			− xн)
		i	i
	Простая проверка на размерность подтверждает правильный характер этих
	формул, так как нетрудно видеть, что размерность коэффициентов b0 ,				b1 и b2

совпадает с размерностью отклика. Кроме того, увеличение интервала варьирования приводит к приблизительно пропорциональному увеличению соответствующих коэффициентов; составляющие градиента при этом практически не изменяются, что и следует из соотношения (1.19).

В ряде руководств (например, [2, 3, 8]) составляющие градиента определяются умножением оценок коэффициентов на интервалы варьирования, что представляется неоправданным в силу приведенных выше аргументов.

Формулы (1.19) используются в случае, когда адекватна линейная модель; когда же она неадекватна, поступают обычно так же, но при этом нужно иметь в виду, что дальнейшее крутое восхождение может оказаться не столь эффективным.

Шаг крутого восхождения выбирается пропорционально градиенту и таким образом, чтобы одна из его координат (в принципе – любая) была примерно равна интервалу варьирования. Остальные коэффициенты шага обычно округляются до ближайших удобных для проведения эксперимента значений. Крутое восхождение обычно начинают из той точки, которая хотя бы по одной координате выходит за пределы предыдущей области экспериментирования. Крутое восхождение особенно эффективно в том случае, когда коэффициенты модели приблизительно одинаковы. Добиться этого можно надлежащим выбором интервалов варьирования каждого фактора. Напомним, что эффект фактора (т.е. вклад фактора в отклик) равен его удвоенному коэффициенту, а величина эффекта приблизительно пропорциональна интервалу варьирования.

1.6.2. Пример крутого восхождения Для иллюстрации изложенного рассмотрим искусственно синтезирован-

ный пример (табл. 1.12). Здесь Sад2 =1,0 , S y2 =1,0 (последняя подсчитана на ос-

новании опытов в центре плана), F = Sад2 / S y2 =1,0 , т.е. линейная модель адек-

ватна. Обе координаты градиента слишком велики – они равны –4,0 и –4,5 соответственно, поэтому уменьшим их в 20 раз и округлим результат. Крутое

36

восхождение начинается из точки x1 = 0,8 , x2 = 5,7 , координаты 5-8-го опытов

записываются в реальном масштабе. Более подробно крутое восхождение и его особенности рассмотрены в работах [2, 3, 8].

			Таблица 1.12
Расчет градиента и план крутого восхождения (план 22)
		x1	x2	y
Основной уровень		1,5	7
Интервал варьирования		0,5	1
Верхний уровень		2	8
Нижний уровень		1	6
	~	~	~
	x0	x1	x2
1	+	+	+	82
2	+	+	–	90
3	+	–	+	85
4	+	–	–	95
Оценка bi	88	–2	–4,5
bi /интервал варьирования		–4	–4,5
Шаг при изменении фактора x1 на 0,2		–0,2	–0,25
Округление		–0,2	–0,3
		x1	x2
5		0,8	5,7	–
6		0,6	5,4	–
7		0,4	5,1	–
8		0,2	4,8	–

1.7.Планы второго порядка

1.7.1.Общие вопросы

Рассмотренные выше планы для построения линейных моделей называются также планами первого порядка. Как было отмечено, эти планы не позволяют учесть квадратичные эффекты, которые могут потребоваться для адекватного описания поверхности отклика в сильно искривленной области, например в области оптимума. Для описания подобных областей необходимы более сложные планы, в которых помимо точек с координатами ±1 должны присутствовать другие точки. Обычно в качестве таких точек выбираются центральная точка с координатами (0,0,...,0) и так называемые звездные точки на осях фак-

торного пространства с координатами (±α,0,...,0) , (0,±α,0,...,0),..., (0,0,...,± α). В

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

37

1.7.2. Ортогональные планы второго порядка Рассмотрим подробно случай двухфакторного эксперимента (рис. 1.6 и

табл. 1.13), который описывается моделью

~

~

~

~2

~2

y =β′0 + β1x1

+ β2 x2

+ β12 x1x2

+ β11x1 + β22 x2 .

~

x2

(+1,–1)

× (0,α)

(+1,+1)

*

*

(–α,0)

(0,0)

(α,0)

~

×

×

x1

(–1,–1) *

× (0, –α)

* (+1,–1)

Рис. 1.6. Размещение точек в факторном пространстве

Двухфакторный план

Таблица 1.13

~

~

~

y

x0

x1

x2

План 22

1

+

+

+

у1

2

+

+

–

у2

3

+

–

+

у3

4

+

–

–

у4

Звёздные точки

5

+

+α

0

у5

6

+

–α

0

у6

7

+

0

+α

у7

8

+

0

–α

у8

Центральные точки

9

+

0

0

у9

10

+

0

0

у10

..

..

..

..

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

~ ~2		1	N ~2
xi′ = xi	−		∑xij	,	i = 1,2,…,n,
		N j=1

38

где N = 2n + 2n + n	– общее число опытов. Тогда, очевидно, получим
0
N ~	~	N ~2	N ~2	= 0 ,	i = 1,2,…,n.
∑x0 j xij′		= ∑xij	− ∑xij
j=1		j=1	j=1

Новая модель будет выглядеть следующим образом:

′	~		~	~		~′	~′	(1.20)
y =β0	+β1x1 +β2 x2			+ β12 x1x2 + β11x1			+ β22 x2 .
Понятно, что коэффициенты β′0 и β0 связаны зависимостью
	β′0 =β0	−	β11	N ~2	β22	N ~2
			N	∑x1 j −	N	∑x2 j .
				j=1		j=1

Столбцы линейных членов и столбец взаимодействия ортогональны столбцу факторов xi′. Для того чтобы были ортогональны столбцы факторов xi′ и x′j (i≠j),

необходимо соответствующим образом подобрать величину звездного плеча α и число опытов в центре плана n0 . Приведем соответствующие выкладки:

N ~ ~

N

~2

1 N

~2

~2

1

N ~2

∑xik′ x′jk = ∑ (xik

−

∑xil

)(x jk

−

∑x jl ) =

k =1

k =1

N l =1

N l =1

N ~2 ~2

1

N ~2

N ~

2

n

(2n + 2α2 )2

= ∑x x

−

∑x

∑x

= 2

−

= 0,

k =1 ik

jk

N l =1

il

l =1

jl

2n + 2n + n0

откуда α = 0,5( 2n N − 2n ) . Выбирая различные значения n0 , получим соответствующие α. Если нет особых требований, можно положить n0 =1. Значения α в этом случае приведены в табл. 1.14.

Приведем матрицу планирования для случая n = 2 (табл. 1.15). Можно проверить, что её столбцы симметричны и попарно ортогональны. Эти свойства обусловливают независимость оценок. Свойство нормированности, однако, не выполняется, а именно

9

~ 2

= 9 ,

9

~ 2

~ 2

= 6 ,

9

~

2

~

2

9

~ 2 ~ 2

) = 4 .

∑ x0 j

∑ x1 j

= ∑x2 j

∑(x'1 j )

= ∑(x'2 j )

= 2 , ∑

(x1 j x2 j

j =1

j =1

j =1

j =1

39

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

1

~

1

(y1

− y3 − y4 + y5 − y6 ).

b1

=

∑y j x1 j =

+ y2

~2

6

∑x1 j

Значения α и N

Таблица 1.14

Число факторов n

2

3

4

5 (план 24–1)

Величина плеча α

1

1,215

1,414

1,517

Число опытов N

9

15

25

27

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

Таблица 1.15

~

~

~

~ ~

~

~

y

x0

x1

x2

x1x2

x1′

x2′

1

+

+

+

+

1/3

1/3

у1

2

+

+

–

–

1/3

1/3

у2

3

+

–

+

–

1/3

1/3

у3

4

+

–

–

+

1/3

1/3

у4

5

+

+1

0

0

1/3

–2/3

у5

6

+

–1

0

0

1/3

–2/3

у6

7

+

0

+1

0

–2/3

1/3

у7

8

+

0

–1

0

–2/3

1/3

у8

9

+

0

0

0

–2/3

–2/3

у9

NZ

9

6

6

4

2

2

Дисперсия воспроизводимости S y2 рассчитывается так же, как в случае линей-

ного планирования. Методика проверки значимости коэффициентов остается прежней, несколько меняется формула (1.12) для t-отношения, в которой число опытов N заменяется соответствующим нормировочным коэффициентом. Это объясняетсятем, чтодисперсииразныхоценокотличаютсядруготдруга, например:

1

~

D y

D y

D y

D[bi ]= D

∑yi xij

=

=

,

D[b0 ]=

.

~2

~2

n

2

rN

r(2

+ 2α

)

∑xij

r∑xij

Точно так же обстоит дело с проверкой адекватности, общая схема которой сохраняется, но конечная формула для дисперсии адекватности выглядит несколько по-другому:

40