§ 4. Оптимизация планов. Поиск экстремума функции отклика

1. Построение линейных оптимальных планов. Пусть функция отклика является полиномом первой степени от k переменных:

. (1)

Определение 1. План называется линейным, если ему соответствует функция отклика (1) и он позволяет получить несмещенные МНК-оценки параметров .

Если матрица линейного плана есть

, (2)

то матрица планирования для этого плана имеет вид

.

Определение 2. Линейный план называется линейным ортогональным планом, если столбцы матрицы X попарно ортогональны, т. е.

. (3)

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

Пусть матрица планирования X имеет ранг k + 1 и вектор наблюдений удовлетворяет условиям ; , тогда для ковариационной матрицы МНК-оценок вектора β имеем равенство

.

Будем предполагать, что

, , (4)

где величины c_i заданы; положим также .

Требуется выбрать матрицу планирования X или матрицу линейного плана , так чтобы дисперсии МНК-оценок параметров в классе линейных планов с ограничениями (4) были минимальными.

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

Теорема Бокса. Пусть функция отклика имеет вид (1), столбцы матрицы линейного плана удовлетворяют условиям (4) и ранг матрицы X равен k + 1. Тогда для МНК-оценок параметров выполняется неравенство

( ), (5)

причем минимум дисперсий в классе линейных планов с ограничениями (4) достигается тогда и только тогда, когда столбцы матрицы X попарно ортогональны.

В частности, поскольку первый столбец матрицы X состоит из единиц (x_0u  1), из условия ортогональности (3) при j  0 следует . Таким образом, получаем

Следствие. Необходимым условием минимума дисперсий МНК-оценок является выполнение при всех i  1, 2, ..., k равенства (условие симметрии плана).

2. Градиентный метод. Будем предполагать, что функция отклика непрерывна и имеет непрерывные частные производные 1-го порядка в ограниченной замкнутой области G  ^k.

Определение 1. Функция f называется унимодальной в области G, если она в этой области имеет единственный экстремум.

Пусть изменение функции многих переменных вдоль некоторой траектории, проведенной из точки в точку и определенной в области G, задается уравнением , где  – параметр, пробегающий числовой отрезок [0; 1].

Траекторию назовем строго возрастающей (строго убывающей), если функция () строго возрастает (строго убывает) на отрезке [0; 1].

Определение 2. Пусть – функция отклика, определенная в области G и имеющая во внутренней точке максимум или минимум. Функция называется строго унимодальной, если отрезок, проведенный из любой точки в точку , является строго возрастающей траекторией в случае максимума или строго убывающей траекторией в случае минимума.

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

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

, (6)

где – градиент функции f в точке ;  – некоторое положительное число. В этом случае функция f предполагается строго унимодальной.

Различные варианты градиентного метода отличаются друг от друга способом выбора . Одним из вариантов градиентного метода является метод наискорейшего спуска (или подъема). При оптимизации технологических процессов в теории планирования эксперимента используется его статистический аналог – метод Бокса и Уилсона. В этом методе используется не сам градиент, а его оценка.

3. Оценивание градиента. Пусть функция отклика

(7)

определена в области G  ^k.

Выберем произвольно точку . Используя эту точку, построим полный факторный эксперимент. Выберем для каждого i  1, 2, ..., k нижний уровень и верхний уровень так, чтобы было . Положим , введем кодированные переменные и выразим функцию отклика (7) через кодированные переменные:

. (8)

Под задачей оценивания градиента будем понимать определение оценки градиента функции отклика (8) в точке x⁰  (0, 0, ..., 0). Предположим, что в окрестности точки x⁰ функция (8) допускает разложение по формуле Маклорена:

(9)

где . Введем обозначения: , , , . Тогда равенство (9) перепишется в виде

.

Так как , задача оценивания градиента сводится к нахождению МНК-оценок параметров .

Пусть матрица полного факторного плана с центром в точке задана равенством (2).

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

. (10)

Тогда матрице плана и функции отклика (10) соответствует матрица планирования ( j  0, 1, 2, ..., k; u  1, 2, ..., N; x_0u  1) и для МНК-оценки параметра имеем равенство

( j  0, 1, 2, ..., k),

где – наблюдения в точках плана.

Поскольку являются оценками компонент градиента , т. е. частных производных , то МНК-оценка градиента функции отклика в точке определяется равенством

.

Пример. Пусть   f (x₁, x₂) – функция отклика (x₁, x₂ –кодированные переменные), матрица плана и результаты наблюдений . Используя аппроксимацию вида   ₀ + ₁x₁ + ₂x₂ + ₃x₁x₂ ( ), найдем оценку градиента в центре плана (0, 0). Так как и , то планирование ортогонально,  (–y₁ + y₂ – y₃ + y₄) : 4  –3,  (–y₁ – y₂ + y₃ + y₄) : 4  –10 и, следовательно, получаем .

4. Метод Бокса и Уилсона. Будем предполагать, что функция отклика (7) в области G строго унимодальна. Поиск ее максимума может быть осуществлен методом, предложенным Боксом и Уилсоном.

Пусть – начальная точка поиска максимума. При переходе к кодированным переменным точке соответствует точка , а функция отклика (7) принимает вид (8).

Пусть – оценка градиента в точке , полученная с использованием факторного эксперимента, матрица плана которого задана равенством (2). Для поиска максимума сделаем некоторый шаг из точки в направлении оценки градиента . Положим

, (11)

где параметр , – МНК-оценка градиента, , (i  1, 2, ..., k).

Точке в системе координат кодированных переменных соответствует точка в системе координат исходных переменных, при этом имеем .

Выполним в точке наблюдения и найдем в ней оценку функции отклика

. (12)

Так как измерение функции отклика происходит с ошибкой, то в точке можно найти только ее оценку, а не точное значение. Предположим, что оценка значимо больше оценки функции отклика в точке , где согласно (10) .

В этом случае делаем второй шаг в направлении оценки градиента и т. д. В общем виде, учитывая (11), для l-го шага имеем:

, (i  1, 2, ..., k), (13)

где . Оценка функции отклика в точке аналогична (12):

,

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

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

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

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

. (14)

Задача проверки гипотезы адекватности этой модели для случая, когда матрица плана (i  1, 2, ..., k; u  1, 2, ..., N₀) является матрицей плана с кратными повторными наблюдениями (l  1, 2, ..., n; s  1, 2, ..., m; всего наблюдений) и когда наблюдения в центре плана отсутствуют, рассматривалась в п. 3 § 3.

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

Оценка дисперсии , связанная с неадекватностью модели, есть , где r – ранг матрицы (см. § 3, п. 6). Учитывая условия (5) – (7) § 3, ортогональность планирования и то, что ранг r матрицы равен , получаем ,

где

, ( j  0, 1, ..., p₀).

Гипотеза H₀ адекватности модели (14) принимается, если , и отклоняется, если .

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

Определение. План, которому соответствует функция отклика вида

,

являющаяся полиномом степени d от переменных x₁, x₂, ..., x_k, называется планом порядка d, если он позволяет получить несмещенные МНК-оценки неизвестных параметров .

Для аппроксимации функции отклика сначала используется полином 2-й степени:

. (15)

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

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

Но полный факторный эксперимент типа использовать для построения планов второго порядка нецелесообразно, так как избыточность опытов   N – ( p + 1), где – число неизвестных параметров в (15), очень велика.

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

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

.

Для их оценивания построим композиционный план. Полный факторный эксперимент типа 2³ образует ядро композиционного плана. В качестве дополнительных точек для наблюдений возьмем еще шесть так называемых «звездных» точек с координатами (–₁, 0, 0), (₁, 0, 0), (0, –₂, 0), (0, ₂, 0), (0, 0, –₃), (0, 0, ₃). Кроме этих точек при построении композиционного плана используют повторных опытов в его центре. Наблюдения в нем или уже имеются до построения композиционного плана, или их выполняют. Они необходимы для проверки гипотезы адекватности модели, а также для получения информации о центре плана.

В этом случае матрица плана при имеет вид: , где обозначено

, .

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

Аналогично строятся центральные композиционные пла­ны 2-го порядка для произвольного числа факторов k, при этом каждый из факторов варьируется на пяти уровнях , – 1, 0, 1, (i  1, 2, ..., k). Число наблюдений выражается, очевидно, равенством

N  N₀ + 2k + n₀, (16)

где N₀ – число наблюдений в точках ядра плана; 2k – число «звездных» точек; n₀ – число наблюдений в центре плана.

В дальнейшем рассматриваются только такие планы, для которых все _i (i  1, 2, ..., k) равны одному и тому же числу . Число  в этом случае называют плечом плана.

Среди центральных композиционных планов наиболее широко применяются планы Бокса.

Определение. Центральный композиционный план 2-го порядка называется планом Бокса, если его ядром является полный факторный эксперимент типа 2^k.

8. Ортогональные планы 2-го порядка. В общем случае центральный композиционный план второго порядка не является ортогональным. Однако если он представляет собой план Бокса, то изменением звездного плеча  и преобразованием функции отклика (15) его можно сделать ортогональным.

Приведем пример таблицы, содержащей матрицу планирования плана Бокса для k = 3 и п₀  3, не все столбцы которой ортогональны друг другу.

План	1 x1 x2 x3				x1x2 x1x3 x2x3
Полный факторный эксперимент типа 23	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 1 1 1 1 1 1 1	1 1 1 1 1 1 1 1
	1	–	0	0	0	0	0	2	0	0
Звездный план	1 1 1 1 1	 0 0 0 0	 –  0 0	0 0  – 	   0 0	   0 0	   0 0	2 0 0 0 0	0 2 2 0 0	0 0 0 2 2
Наблюдения в центре плана	1 1 1	0 0 0	0 0 0	0 0 0	0 0 0	0 0 0	0 0 0	0 0 0	0 0 0	0 0 0

Здесь имеем при i  1, 2, 3; при i, j  1, 2, 3, i  j .

Пусть (i  1, 2, ..., k; и  1, 2, ..., N) – матрица плана Бокса, где N определяется формулой (16), а X – матрица планирования, соответствующая матрице плана и функции отклика (15).

Для того чтобы сделать план Бокса ортогональным, введем переменные

(i  1, 2, ..., k).

В результате от модели (15) перейдем к модели

(17)

где

. (18)

У моделей (15) и (17) все неизвестные коэффициенты, кроме первого, совпадают. Очевидно, что

(i  1, 2, ..., k). (19)

Обозначим

. (20)

Матрица планирования X для модели (15) будет отличаться от матрицы планирования для модели (17). Таблица, приведенная на стр. 140, содержит матрицу для случая k = 3 и п₀ = 3.

Для модели (17) столбец матрицы , соответствующий переменной x₀, будет ортогонален каждому столбцу, соответствующему переменной . Действительно,

(i  1, 2, ..., k).

Неортогональными друг другу будут лишь последние k столбцов этой матрицы: , i, j  1, 2, ..., k; i  j.

Требуется подобрать величину  так, чтобы эти столбцы стали попарно ортогональными. Для этого нужно ре-

План	x0 x1 x2 x3				x1x2 x1x3 x2x3
Полный факторный эксперимент типа 23	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– c 1– c 1– c 1– c 1– c 1– c 1– c 1– c	1– c 1– c 1– c 1– c 1– c 1– c 1– c 1– c	1– c 1– c 1– c 1– c 1– c 1– c 1– c 1– c
Звездный план	1 1 1 1 1 1	–  0 0 0 0	0  –  0 0	0 0 0  – 	0    0 0	0    0 0	0    0 0	2– c 2– c – c – c – c – c	– c – c 2 – c 2 – c – c – c	– c – c – c – c 2 – c 2 – c
Наблюдения в центре плана	1 1 1	0 0 0	0 0 0	0 0 0	0 0 0	0 0 0	0 0 0	– c – c – c	– c – c – c	– c – c – c

шить уравнение (см. приведенную таблицу)

,

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

. (21)

Учитывая (16), (18), (19) и (20), получаем с²  N₀ / N. Из уравнения (21) следует, что .

Зависимость  и c от числа факторов при n₀  3 показана в следующей таблице.

Число факторов	2	3	4
Ядро
c	0,603	0,686	0,770
	1,147	1,353	1,547

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

, j  1, 2, ..., p.

Дисперсия оценки определяется по теореме Бокса:

, j  1, 2, ..., p. (22)

Оценка параметра определяется формулой

,

ее дисперсия

. (23)

Согласно (18) и (19) и, следовательно, .

Дисперсии оценок коэффициентов регрессии для планов второго порядка в отличие от линейных планов являются различными (см. (22) и (23)).

Оценкой функции отклика в точке является

,

а оценка ее дисперсии есть

.

Пример. Пусть k  3, N₀  2³ и n₀  3. Построим центральный композиционный план 2-го порядка. Так как в этом случае   1,35, c  0,69, то матрица планирования для ортогонального плана Бокса и функции отклика (17) имеет вид

Таблица 1. Значения функции

x	0	1	2	3	4	5	6	7	8	9
0,0	0,3989	3989	3989	3988	3986	3984	3982	3980	3977	3973
0,1	3970	3965	3961	3956	3951	3945	3939	3932	3925	3918
0,2	3910	3902	3894	3885	3876	3867	3857	3847	3836	3925
0,3	3814	3802	3790	3778	3765	3752	3739	3726	3712	3697
0,4	3683	3668	3653	3637	3621	3605	3589	3572	3555	3538
0,5	3521	3503	3485	3467	3448	3429	3410	3391	3372	3352
0,6	3332	3312	3292	3271	3251	3230	3209	3187	3166	3144
0,7	3123	3101	3079	3056	3034	3011	2989	2966	2943	2920
0,8	2897	2874	2850	2827	2803	2780	2756	2732	2709	2685
0,9	2661	2637	2613	2589	2565	2541	2516	2492	2468	2444
1,0	0,2420	2396	2371	2347	2323	2299	2275	2251	2227	2203
1,1	2179	2155	2131	2107	2083	2059	2036	2012	1989	1965
1,2	1942	1919	1895	1872	1849	1826	1804	1781	1758	1736
1,3	1714	1691	1669	1647	1626	1604	1582	1561	1539	1518
1,4	1497	1476	1456	1435	1415	1394	1374	1354	1334	1315
1,5	1295	1276	1257	1238	1219	1200	1182	1163	1145	1127
1,6	1109	1092	1074	1057	1040	1023	1006	0989	0973	0957
1,7	0940	0925	0909	0893	0878	0863	0848	0833	0818	0804
1,8	0790	0775	0761	0748	0734	0721	0707	0694	0681	0669
1,9	0656	0644	0632	0620	0608	0596	0584	0573	0562	0551
2,0	0,0540	0529	0519	0508	0498	0488	0478	0468	0459	0449
2,1	0440	0431	0422	0413	0404	0396	0387	0379	0371	0363
2,2	0355	0347	0339	0332	0325	0317	0310	0303	0297	0290
2,3	0283	0277	0270	0264	0258	0252	0246	0241	0235	0229
2,4	0224	0219	0213	0208	0203	0198	0194	0189	0184	0180
2,5	0175	0171	0167	0163	0158	0154	0151	0147	0143	0139
2,6	0136	0132	0129	0126	0122	0119	0116	0113	0110	0107
2,7	0104	0101	0099	0096	0093	0091	0088	0086	0084	0081
2,8	0079	0077	0075	0073	0071	0069	0067	0065	0063	0061
2,9	0060	0058	0056	0055	0053	0051	0050	0048	0047	0046
3,0	0,0044	0043	0042	0040	0039	0038	0037	0036	0035	0034
3,1	0033	0032	0031	0030	0029	0028	0027	0026	0025	0025
3,2	0024	0023	0022	0022	0021	0020	0020	0019	0018	0018
3,3	0017	0017	0016	0016	0015	0015	0014	0014	0013	0013
3,4	0012	0012	0012	0011	0011	0010	0010	0010	0009	0009
3,5	0009	0008	0008	0008	0008	0007	0007	0007	0007	0006
3,6	0006	0006	0006	0005	0005	0005	0005	0005	0005	0004
3,7	0004	0004	0004	0004	0004	0004	0003	0003	0003	0003
3,8	0003	0003	0003	0003	0003	0002	0002	0002	0002	0002
3,9	0002	0002	0002	0002	0002	0002	0002	0002	0001	0001

Таблица 2. Значения функции

x	0	1	2	3	4	5	6	7	8	9
0,0	0,00000	00399	00798	01197	01595	01994	02392	02790	03188	03586
0,1	03983	04380	04776	05172	05567	05962	06356	06749	07142	07535
0,2	07926	08317	08706	09095	09483	09871	10257	10642	11026	11409
0,3	11791	12172	12552	12930	13307	13683	14058	14431	14803	15173
0,4	15542	15910	16276	16640	17003	17364	17724	18082	18439	18793
0,5	19146	19497	19847	20194	20540	20884	21226	21566	21904	22240
0,6	22575	22907	23237	23565	23891	24215	24537	24857	25175	25490
0,7	25804	26115	26424	26730	27035	27337	27637	27935	28230	28524
0,8	28814	29103	29389	29673	29955	30234	30511	30785	31057	31327
0,9	31594	31859	32121	32381	32639	32894	33147	33398	33646	33891
1,0	34134	34375	34614	34850	35083	35314	35543	35769	35993	36214
1,1	36433	36650	36864	37076	37286	37493	37698	37900	38100	38298
1,2	38493	38686	38877	39065	39251	39435	39617	39796	39973	40147
1,3	40320	40490	40658	40824	40988	41149	41309	41466	41621	41774
1,4	41924	42073	42220	42364	42507	42647	42786	42922	43056	43189
1,5	43319	43448	43574	43699	43822	43943	44062	44179	44295	44408
1,6	44520	44630	44738	44845	44950	45053	45154	45254	45352	45449
1,7	45543	45637	45728	45818	45907	45994	46080	46164	46246	46327
1,8	46407	46485	46562	46638	46712	46784	46856	46926	46995	47062
1,9	47128	47193	47257	47320	47381	47441	47500	47558	47615	47670
2,0	47725	47778	47831	47882	47932	47982	48030	48077	48124	48169
2,1	48214	48257	48300	48341	48382	48422	48461	48500	48537	48574
2,2	48610	48645	48679	48713	48745	48778	48809	48840	48870	48899
2,3	48928	48956	48983	49010	49036	49061	49086	49111	49134	49158
2,4	49180	49202	49224	49245	49266	49286	49305	49324	49343	49361
2,5	49379	49396	49413	49430	49446	49461	49477	49492	49506	49520
2,6	49534	49547	49560	49573	49585	49598	49609	49621	49632	49643
2,7	49653	49664	49674	49683	49693	49702	49711	49720	49728	49736
2,8	49744	49752	49760	49767	49774	49781	49788	49795	49801	49807
2,9	49813	49819	49825	49831	49836	49841	49846	49851	49856	49861
3,0 3,5 4,0 4,5 5,0	0,49865 0,49977 0,499968 0,499997 0,49999997		3,1 3,6	49903 49984	3,2 3,7	49931 49989	3,3 3,8	49952 49993	3,4 3,9	49966 49995

Таблица 3. Корни уравнений P(² < x)  q и P(² > y) 