close

Вход

Забыли?

вход по аккаунту

код для вставкиСкачать
Министерство образования и науки Российской Федерации
Новосибирский государственный технический университет
Кафедра программных систем и баз данных
КУРСОВОЙ ПРОЕКТ
по дисциплине “Планирование и анализ эксперимента”
Факультет:
Группа:
Студент:
Преподаватель:
Вариант:
ПМИ
ПМ-83
Пустосмехов В.А.
Гультяева Т.А.
23
Новосибирск
2012
1. Цель работы
Изучение методики оптимизации функции отклика по методу Бокса-Уилсона
2. Постановка задачи
y( x1, x2 )  x2 sin x12  e
 
Начальная точка: x(0)   0.4, 0.1
Ограничения на факторы: 2  x1  2 , 2  x2  2
Условие: y( x1, x2 )  max
1. Реализовать программу, которая позволяет выбирать значение параметров и просматривать
протокол решения.
2. Используя алгоритм движения к почти стационарной области на основе использования метода
крутого восхождения найти такую область.
3. Описать найденную область уравнением регрессии 2-ого порядка
4. Привести уравнение регрессии к канонической форме.
5. Определить тип поверхности.
3. Ход работы
Процедура оптимизации по методу Бокса-Уилсона делится на 2 этапа:
1. движение к «почти стационарной области», близкой к экстремальному
2. уточнение положения экстремума в ней.
Алгоритм движения к «почти стационарной области» следующий:
1. Выбор одного из факторов xl за базовый и выбор для него шага hl движения.
Выбор интервалов варьирования  i , i  1, m .
Выбор числа k параллельных измерений в точках ПФЭ. В программе проводилось по 3
параллельных измерения.
Номер итерации j = 0.
( j)
Выбор начального приближения x .
( j)
( j)
2. Проведение k измерений функции отклика в точке x , рассчитать y  y , sb2  sb2( j )
, где
2
sb2
s
s2y
y
(или sb2 
– если есть параллельные опыты),

N
N
s2y
2
s 
,
y k 1
1 N
s2y   s2j , ,
N j 1


s2j 
2
1 k
y ji  y j ,

k  1 i 1
yj 
1 k
y ,
k i 1 ji
j  1, N ,
j  1, N ,
N - число серий опытов (то есть в случае измерения дисперсии в центре N = 1), k – число опытов
в каждой серии.
2
( j)
3. Перенос начала координат факторного пространства (т.е. кодированных переменных) в точку x
– это точка будет центром плана. Переход к кодированным координатам заключается во введении
новых кодированных переменных:
x x
X i  i 0i , i  1, n
i
,
где  i – масштаб по оси xi .
В точках ПФЭ провести по
k измерений функции отклика.
Проведение ПФЭ в пространстве кодированных переменных. Расчет ОМП оценок bˆi , проверка
значимости bˆi , i  1, m .
̂ = ( )−1  
Значимость каждого коэффициента уравнения регрессии устанавливают с помощью критерия
Стьюдента, вычисляя его расчетное значение:
tp 
b
sb2
.
где b - коэффициент уравнения регрессии, для которого устанавливается значимость. Каждое
рассчитанное значение сравнивают с табличным значением критерия Стьюдента, которое выбирают для
заданного уровня значимости p при числе степеней свободы f  N  k 1 . Если выполняется условие
t p  t , то коэффициент считается значимым.
4. Если все bˆi незначимы или y  bˆ0 , то необходимо перейти к шагу 10.
(Если все коэффициенты, соответствующие влиянию факторов (но не взаимодействию),
незначимы, то градиент равен нулю и никуда не направлен. Следовательно, поиск в направлении
градиента невозможен и следует перейти ко второму этапу метода Бокса-Уилсона)
Иначе рассчитать шаги hi , i  l, i  1, m
При выборе шага движения один из факторов принимают за базовый и для него выбирают шаг hl
движения. Для всех остальных факторов шаг движения рассчитывают по формуле:
b
hi  hl i i
b l l ,
где hi - шаг движения для i -ого фактора; hl - шаг движения для базового фактора l; b i и b l –
коэффициенты регрессионного уравнения,  i , l – интервал варьирования.
5. Проверить адекватность модели, в зависимости от которой будут приниматься решения о
проведении о реализации опытов.
Проверка его адекватности осуществляется с помощью критерия Фишера, который представляет
собой отношение:
2
sад
Fp  2
s ,
y
2 – оценка дисперсии адекватности, которая вычисляется как:
где sад
2
1 N
2
sàä 
y  yˆ j
N  m j 1 j


3
где m – число коэффициентов регрессии искомого уравнения, включая и свободный член; y j , yˆ j
– среднее экспериментальное и расчетное значение функции отклика в j-м опыте; N – число опытов
ПФЭ. С оценкой дисперсии адекватности связано число степеней свободы fад  N  m . Уравнение
регрессии считается адекватным, если выполняется условие:
Fp  F , где F – значение критерия Фишера .
6.
l  0,
0
j
z   x  .
i
i
l 1
l 
Сделать мысленный шаг zi
 zi  hi , i  1, m. Принять решение о его реализации.
Особого внимания заслуживает порядок реализации мысленных опытов (опытов, которые связаны с
движением по градиенту)
При использовании адекватной линейной модели целесообразно реализовать только те мысленные
опыты, которые заходят за область ПФЭ хотя бы по одному из учитываемых факторов. Если крутое
восхождение ведется с помощью линейной модели, которая не является адекватной, то часть (1–3)
мысленных опытов осуществляют в области ПФЭ, а остальные, если это разумно, за пределами этой
области
Значения критерия оптимизации в области эксперимента находят простым расчетом yˆ , а за ее
пределами реализуют часть мысленных опытов.
l  l 1
Повторять шаг до тех пор, пока не примется решение о реализации мысленного опыта.
 j 1 , которая соответствует zl 1 , k
i
7. Реализовать мысленный опыт, проведя в точке xi
( j 1) 2( j 1)
измерений функции отклика, рассчитать y
, sb
( j 1)
( j)
8. Если шаг удачен, т.е. y
 y , j  j  1 , то перейти на шаг 6.
Иначе: если число удачных шагов меньше 2, то перейти к шагу 10. Если число удачных шагов больше
2, то перейти на шаг 9.
( j)
9. Подготовка к новой итерации ( j – номер последнего удачного шага): y  y , sb2  sb2( j ) .
Перейти на шаг 3.
10. Второй этап метода Бокса-Уилсона – описание области оптимума методами нелинейного
( j)
планирования с переносом начала координат факторного пространства в точку x .
4
Текст программы:
#!/usr/bin/env python
# coding=utf8
from itertools import chain
from math import sin, sqrt
from numpy import array, average, matrix, random
from rpy2 import robjects
sign = lambda x: 0 if x == 0 else 1 if x > 0 else -1
σ = 0.01
f = lambda x1, x2: x2 * sin(x1 ** 2)
x0 = (0.4, -0.1)
x1 = [-2, 2]
x2 = [-2, 2]
# Матрица ПФЭ
X = matrix([[
[
[
[
[
[
[
[
[
[
[
[
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.1, 0.1) # Интервал варьирования
x = array(x0)
# Начальное приближение
points = [x]
# Точки плана
while True:
def measure(x):
return f(*x) + random.normal(0, σ)
y = [measure(x), measure(x), measure(x)]
# Измерения в точке x
sb2 = sum([(yi - average(y)) ** 2 for yi in y]) / 2 # Вычисление sb2
FFE_points = [x + ε, x - ε, x + (-ε[0], ε[1]), x + (ε[0], -ε[1])]
# Точки ПФЭ
y = list(chain(*[[measure(x), measure(x), measure(x)] for x in FFE_points])) # Делаем в них по 3
# измерения
points = points + FFE_points
# Добавляем их в план
b = (X.T * X).I * X.T * matrix(y).T
b = b.T.tolist()[0]
# Вычисляем b
# Проверям их значимость
significant = lambda b: abs(b) / sqrt(sb2) > robjects.r("qt(0.95, 2)")[0] if sb2 != 0 else True
if not significant(b[1]) and not significant(b[2]):
# Конец первого этапа
break
#
h
#
h
Выбираем шаг движения. Базовый фактор второй,
= (0.2 * sign(b[2]) * b[1] / b[2], 0.2 * sign(b[2]))
Если фактор незначим, то по нему не надо двигаться
= (0 if not significant(b[1]) else h[0], 0 if not significant(b[2]) else h[1])
average3 = lambda v: [average(x) for x in list(zip(*[v[i::3] for i in range(3)]))] # Усреднить по 3
ẙ = average3(y)
# Среднее экспериментальное значение функции отклика
ỹ = average3(X * matrix(b).T)
# Расчётное значение функции отклика
sad2 = sum([(ỹ[i] - ẙ[i]) ** 2 for i in range(0, 4)])
# Оценка дисперсии адекватности
adequate = sad2 / sb2 <= robjects.r("qf(0.95, 1, 2)")[0] if sb2 != 0 else
print("Адекватность: ", adequate, "b: ", b, "h: ", h)
def stay_in_area(x):
x[0] = max(x[0],
x[0] = min(x[0],
x[1] = max(x[1],
x[1] = min(x[1],
return x
x1[0])
x1[1])
x2[0])
x2[1])
x_next = stay_in_area(x + h)
if adequate:
# Выходим за область ПФЭ
5
while max(x[0] - ε[0], x1[0]) < x_next[0] and x_next[0] < min(x[0] + ε[1], x1[1]) and max(x[1]
- ε[0], x2[0]) < x_next[1] and x_next[1] < min(x[1] + ε[1], x2[1]):
x_next = stay_in_area(x_next + h)
# Производим мысленные опыты
step = 0
break_main_loop = False
while step < 100:
points = points + [x_next]
y_next = [measure(x_next), measure(x_next), measure(x_next)]
# Добавляем точку в план
# Проводим мысленный опыт
if average(y_next) > average(y):
# Шаг удачен
x_next = stay_in_area(x_next + h)
y = y_next
else:
if step == 0:
break_main_loop = True
else:
# Подготовка новой итерации
x = stay_in_area(x_next - h)
break
step = step + 1
if break_main_loop:
break
print(x)
print(points)
Незашумлённая функция имеет максимум в точке (±1.2533, 2), равный 2.
В качестве базового фактора был выбран второй.
1 = 0.1, 2 = 0.1
Протокол решения программы:
Итерация Адекватность b0
b1
b2
1
-0.0248
-0.0148
0.0154
2
+
1.4173
0.2158
0.0862
Найдена точка (-1.3367, 1.9), значение функции 1.8558
h1
-0.1929
0.5006
h2
0.2
0.2
Точки, в которых проводились измерения функции отклика:
2.5
2
1.5
1
0.5
0
-2
-1.5
-1
-0.5
0
-0.5
6
0.5
1
Описание найденной области уравнением регрессии второго порядка:
Для адекватного математического описания здесь требуется многочлен более высокой степени,
например, отрезок ряда Тейлора, содержащий члены с квадратами переменных:
y  b0  b1X1  b2 X 2  ...  bn X n  b12 X1X 2  ...  b(n1)n X n1X n 
b11X12  b22 X 22  ...  bnn X n2.
При ротатабельном ЦКП для вычисления коэффициентов регрессии и соответствующих оценок
дисперсий находят следующие константы:
где
n–
число факторов;
N – общее число опытов ротатабельного ЦКП;
N0 – число опытов в центре плана.
Далее в формулах вместо s 2y используется s 2 . В каждой точке проводится k  3 параллельных
y
опытов.
На основании результатов эксперимента вычисляют следующие суммы:
Формулы для расчета коэффициентов регрессии имеют вид:
Оценки дисперсий в определении коэффициентов регрессии вычисляют по следующим формулам:
7
Воспользовавшись данными формулами, получим уравнение регрессии следующего вида:
 = −0.0564 + 0.01231 + 0.10942 − 0.02411 2 − 0.034512 − 0.037622
Приведем его к каноническому виду. Сначала запишем систему:

= 0.0123 + 0.06901 + 0.02412 = 0
1
,

= 0.1094 + 0.02411 + 0.07532 = 0;
{2
Решим её:
1 = −1,3708; 2 = 1.9640; (1 , 2 ) = 1,8727
Составим характеристическое уравнение вида:
 −  0.512
| 11
| = 0, 11 = 0.0345, 12 = 21 = −0.0241, 22 = 0.0376
0.521 22 − 
Раскрыв определитель, получим:
−0.1684 − 0.1879 + 2 = 0
11 = −0.0254, 12 = 0.0975
Выполнено условие
11 + 22 = 0.0721
11 + 22 = 0.0721
значит, коэффициенты канонической формы найдены верно.
Уравнение регрессии в канонической форме имеет вид:
 = 1,8727 − 0.025412 + 0.066222
Коэффициенты канонической формы имеют разный знак, поверхность отклика имеет вид «седловины».
Чтобы увеличивать значение Y, следует двигаться от центра поверхности по направлению 22 .
Найдем соотношение между переменными 1 , 2 и 1 , 2 . Составим систему уравнений:
(b11  B11 )m11  0.5b12  m12  0,

0.5b12  m11  (b22  B11 )m12  0.
Подставив известные значения коэффициентов, получим
11 = 1, 12 = −1.1377
Теперь найдем величины:
11
11 =
= 0.6601,
2
2
√11
+ 12
12
12 =
= −0.7511
2
2
√11
+ 12
Составим вторую систему уравнений:
(b11  B22 )m21  0.5b12  m22  0,

0.5b12  m21  (b22  B22 )m22  0.
Подставив известные значения коэффициентов, получим
21 = 1, 22 = −0.8789
Теперь найдем величины:
21
21 =
= 0.7511,
2
2
√21
+ 22
22
22 =
= 0.6601
2
2
√21
+ 22
Теперь можно представить связь между координатами 1 , 2 и 1 , 2 .
Z1  M11 ( X1  X1s )  M12 ( X 2  X 2 s ),
Z2  M 21( X1  X1s )  M 22 ( X 2  X 2 s ).
8
Подставив коэффициенты, получим:
1 = 0.6601(1 + 1.3780) − 0.7511(2 − 1.9640),
2 = 0.7511(1 + 1.3780) − 0.6601(2 − 1.9640).
Для нахождения угла поворота новых координатных осей относительно старых вычислим величину
12
(2) =  −
= 0.09, откуда  = 5.14. Угол положителен, координаты повернуты против часовой
11
22
стрелки.
Вывод по проделанной работе:
Процедура оптимизации зашумлённой функции недетерминирована: кажущееся увеличение или
уменьшение её значения может быть вызвано шумом, что существенно осложняет поиск экстремума.
Поэтому, к подбору таких параметров алгоритма, как интервалы варьирования факторов и величина
шага следует относиться с особой тщательностью: для исследованной функции при слишком малой
величине шага план содержал сотни точек, при слишком большой останавливался на первой же
итерации ввиду невозможности совершить удачный шаг.
9
1/--страниц
Пожаловаться на содержимое документа