Крутое восхождение по поверхности отклика

Новый подход к решению задачи поиска оптимума был предложен Боксом и Уильсоном в 1951 г. [2]. Они предложили использовать последовательный метод изучения поверхности отклика, напоминающий итерационный метод решения задач вычислительной математики.

Согласно концепции Бокса–Уилсона задача оптимизации условно разбиваете на два этапа. Первый этап – крутое восхождение с целью скорейшего достижения области оптимума. При этом используется линейное планирование. Линейный план может использоваться один или несколько раз в зависимости от интенсивности продвижения. Второй этап – описание области оптимума методами нелинейного планирования [2]. Если область оптимума не достигнута, в этом случае строится линейный план следующего цикла и исследование продолжается.

Алгоритм поиска оптимальных условий протекания процесса изображен на рис. 1.

Опишем процедуру поиска экстремума функции отклика. Ставится ДФЭ с шагом варьирования Δz_i , i=1,2,..к . Получают линейное описание поверхности отклика:

(1)

где коэффициенты регрессии представляют собой частные производные:

(2)

Начало

1.Описание небольших участков поверхности отклика линейным уравнением

2. Нахождение частных производных

3. Движение в направлении крутого восхождения

4. Попали в почти стационарную область?

5. Описание почти стационарной области

6. Нахождение оптимальных значений факторов

Конец

Рисунок 1 – Алгоритм поиска оптимальных условий протекания

процесса

Градиент функции отклика y задается выражением:

(3)

где i, j, ….k – единичные векторы (орты) в направлении координатных осей.

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

i =1..k, (4)

где с – константа, шаг изменения факторов при крутом восхождении.

При постановке эксперимента всегда приходится переходить к натуральным переменным z_i. В натуральных переменных величина шага принимается равной:

c_i^’ = сΔz_i. (5)

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

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

Пример. Исследуется система состав стали и семи компонентов [2]. Необходимо найти такое соотношение факторов, обеспечивающее наибольшее значение показателя «у» - сопротивление при испытаниях на разрыв.

Реализован ДФЭ 2^7-4, в первых строках табл. 1 записаны интервалы варьирования и указаны кодовые обозначения переменных. В строках 7-14 дана матрица планирования. В строке 15 указаны численные значения коэффициентов регрессии.

Линейное приближение функции отклика имеет вид:

b₀=3,893.

Поскольку линейные эффекты оценивались независимо от парных взаимодействий, крутое восхождение проводилось по градиенту линейного приближения. Опыты крутого восхождении реализуются с использованием фактических значений факторных переменных. В строке 16 приведены величины, полученные после умножения коэффициентов регрессии на интервал варьирования. Пользуясь этими величинами, выполнена серия мысленных опытов (строки 19-22) для крутого восхождения по поверхности отклика. Результаты мысленных опытов вычислялись по уравнению регрессии, выраженному через фактические переменные.

Опыт 9, соответствующий строке 23, был реализован и дал хороший результат. Далее были реализованы 25-28 . Лучший результат получен в опыте 11 (строка 26). Дальнейшее продвижение по градиенту приводит к уменьшению сопротивления при испытаниях на разрыв.

Таблица 1 – Планирование эксперимента при выборе оптимального состава легированной стали

	Наименование факторов	Cr	Ni	Mo	V	Nb	Mn	C	у3*103
	Основной уровень Zi0 %			0,1	0,02	0,1	0,4	0,4	6,809
	Интервал варьирования Δzi			0,1	0,02	0,1	0,1	0,1
	Верхний уровень Zimax			0,2	0,04	0,2	0,5	0,5
	Нижний уровеньZimin						0,3	0,3
	Кодовые обозначения переменных	х1	х2	х3	х4	х5	х6	х7
	Опыт 1	-	-	-	-	-	-	-	1,5
	Опыт 2	+	+	-	-	+	+	-	3,5
	Опыт 3	+	-	+	-	+	-	+	6,2
	Опыт 4	-	+	+	-	-	+	+	3,2
	Опыт 5	+	-	-	+	-	+	+	5,3
	Опыт 6	-	+	-	+	+	-	+	5,1
	Опыт 7	-	-	+	+	+	+	-	5,3
	Опыт 8	+	+	+	+	-	-	-	5,8
	bi	0,71	-0,09	0,64	0,89	0,54	-0,16	0,46
	bi* Δхi	0,71	-0,09	0,064	0,018	0,054	-0,016	0,046
	C=1
	Шаг (с округлением)	0,8	-0,1	0,07	0,02	0,06	-0,02	0,05
	Мысленный шаг 1	4,8	1,9	0,17	0,04	0,16	-0,38	0,45	7,507
	Шаг 2	5,6	1,8	0,24	0,06	0,22	0,36	0,5
	Шаг 3	6,4	1,7	0,31	0,08	0,28	0,34	0,55
	Шаг 4	7,2	1,6	0,38	0,1	0,34	0,32	0,6	9,62
	Реализованный опыт 9	8,0	1,5	0,45	0,12	0,4	0,3	0,65	10,3
	Мысленный опыт	8,8	1,4	0,52	0,14	0,46	0,28	0,7	11.0
	Реализованный опыт 10	9,6	1,3	0,59	0,16	0,52	0,26	0,75	11,0
	Реализованный опыт 11	10,4	1,2	0,66	0,18	0,58	0,24	0,8	11,5
	Реализованный опыт 12	11,2	1,1	0,73	0,2	0,64	0,22	0,85	11,2
	Реализованный опыт 13		1,0	0,8	0,22	0,7	0,2	0,9	10,1

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

Контрольные вопросы

1. Как определяется градиент функции отклика?

2. От чего зависит величина составляющих градиента?

3. Как выбирается шаг движения по градиенту?

4. Что представляет мысленный опыт и алгоритмы его реализации?

5. Из какой точки плана начинают движение по градиенту?

6. В каком случае движение по градиенту считается эффективным?

7. Какие возможны исходы при крутом восхождении?

8. Как учитывается возможный временной дрейф при крутом восхождении?

Лекция4 Ротатабельное планирование второго порядка

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

Ротатабельным является такое планирование, у которого корреляционая матрица (Х^тХ)^-1 инвариантна к ортогональному вращению координат [2]. Это условие удовлетворяется для плана второго порядка, если все нечеткие моменты до четвертого порядка равны нулю, а для четных моментов имеет место соотношение:

(1)

где λ₂, λ₄ – некоторые константы, удовлетворяющие неравенству ;

N=2^k+2k+n₀ – общее количество опытов в плане; k – количество факторов.

Свойство ротатабельности не инвариантно к изменению масштабности независимых переменных. Вид информационного профиля зависит от λ₄ .

Иногда интерес представляет информация о функции отклика в некоторой окрестности центра плана. В этом случае следует добиться одинаковой погрешности модели внутри гиперсферы единичного радиуса. План, обеспечивающий такое свойство функции отклика, называется униформ-ротатабельным (λ₄<=1). При этом теряется ортогональность:

Для его формирования достаточно обеспечить равенство дисперсии в центре плана (ρ = 0) и на поверхности гиперсферы радиуса ρ = 1. Этого добиваются подбором числа наблюдений n₀ в центре плана, а именно параметр λ₄ следует взять равным положительному корню квадратного уравнения [11]:

2λ₄ (λ₄ – 1)(k + 2) + λ₄ (k + 1) – (k – 1) = 0.

Ротатабельные планы с равномерным расположением точек на сфере приводят к выражденным матрицам (Х^тХ). Для устранения этого принимают комбинации ротатабельных планов с различным радиусом сферы.

Точки ротатабельного центрального композиционного плана (РЦКП) Бокса второго порядка располагают на концентрических гиперсферах, количество которых не менее двух. Первая гиперсфера может быть вырожденной, т.е. представлять собой центральную точку плана, ее радиус ρ₁ = 0. Именно такая сфера часто используется на практике.

Вторая гиперсфера соответствует вписанному в нее кубу, выбранному в качестве ядра плана. Для ядра х_i = 1, следовательно, радиус этой гиперсферы равен:

ρ₂ = (х₁² + х₂² + … + х_k²)^1/2 = (k)^1/2.

Ядро представляет собой ПФЭ вида 2^k или ДФЭ вида 2^{k – p} , причем должно соблюдаться условие (k – p)/4 > 3/4. Следовательно, с учетом ограничений на ЦКП Бокса, если k ≥ 5, то в качестве ядра можно использовать полуреплику, если k ≥ 8, ядром может служить четверть реплики.

Третья гиперсфера имеет радиус ρ₃ = 2 ^{k / 4} для ядра в виде ПФЭ и радиус ρ₃=2^(k-p)/4 для ядра в виде ДФЭ.

Таким образом, каждый фактор в РЦКП Бокса варьируется на пяти уровнях. В некоторых случаях радиусы второй и третьей гиперсферы совпадают:

k = 2:

ρ₂ = 2 ^1/2, ρ₃ = 2 ^2/4 = 2^1/2;

k = 8 и p = 2:

ρ₂ = 8 ^1/2 = 2 ^3/2, ρ₃ = 2 ^{(8 – 2)/4} = 2^3/2.

Для каждого из факторов кроме нижнего, верхнего и основного уровней устанавливаются два дополнительных уровня “±a”, где a вычисляется по формуле:

,

где k – количество факторов.

В таблице 1 приведены необходимые сведения для составления ротатабельных центральных композиционных планов (РЦКП).

Таблица 1 - Сведения для составления ротатабельных центральных композиционных планов

Количество факторов	Число точек ПФЭ	Число звездных точек	Число точек в центре эксперимента	Значение α
				1,414
				1,682
				2,000
				2,378
5, полуреплика				2,000
				2,828
6, полуреплика				2,378
				3,333
7, полуреплика				2,828

Рассмотрим для примера случай объекта с двумя входными переменными (k = 2) [10]. Точки исходного симплексного плана находятся в вершинах равностороннего треугольника (рис.1).

x1

x2

Рисунок 1 – Построение РЦКП с двумя факторными переменными

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

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

Пример. В качестве примера рассмотрим матрицу рототабельного униформ-планирования второго порядка для k=3 [2]. При использовании полинома в качестве математической модели процесса факторы кодируют по формуле

где x_i – кодовое значение i-го фактора, z_i – натуральное текущее значение i-го фактора, z_i0 – начальный уровень фактора, ∆z_i – интервал варьирования i-го фактора.

Таблица 2 - Матрица планирования РЦКП (униформ-планирование)

второго порядка для k=3

x0	x1	x2	x3	x12	x22	x32	x1x2	x1x3	x2x3
	-1 -1 -1 -1 -1,682 1,682	-1 -1 -1 -1 -1,682 1,682	-1 -1 -1 -1 -1,682 1,682	2,828 2,828	2,828 2,828	2,828 2,828	-1 -1 -1 -1	-1 -1 -1 -1	-1 -1 -1 -1

Коэффициенты модели и их дисперсии рассчитываются по формулам [11]:

,

;

;

;

;

;

;

.

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

cov(β ₀, β _ii) = – 2σ²(ỹ) λ₄ A/N ;

cov(β _ii, β _ij) = σ² (ỹ) (1–λ₄ )A/N.

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

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

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

отличается от нуля. Здесь – величина, предсказанная уравнением модели, – найденная экспериментально.

Величина

σ_R² =S²_R / [N–(k+1)(k+2)/2]

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

На основании рассчитанных величин можно провести все необходимые проверки коэффициентов и модели в целом. Адекватность модели проверяется по критерию Фишера:

.

Значимость коэффициентов модели проверяется по критерию Стьюдента аналогично определению значимости при ортогональном планировании эксперимента:

.

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

Контрольные вопросы

1. Что такое «композиционный план второго порядка»?

2. Какое планирование обеспечивает равномерность распределения информации в уравнении функции отклика по всем направлениям?

3. Свойства униформ-ротатабельных планов?

4. Расположение точек ротатабельного центрального композиционного плана (РЦКП)?

5. На скольких уровнях варьируется каждый фактор в РЦКП Бокса?

6. Как выбирается величина звездного плеча α и чему равно их количество в РЦКП?

7. Из каких блоков состоят композиционные планы, предложенные Боксом и Уильсоном?