close

Вход

Забыли?

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

?

AfanasenkoObuhova

код для вставкиСкачать
МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ
Федеральное государственное автономное образовательное учреждение
высшего профессионального образования
САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ
АЭРОКОСМИЧЕСКОГО ПРИБОРОСТРОЕНИЯ
А. С. Афанасенко, Н. А. Обухова, Б. С. Тимофеев
ОСНОВЫ ТЕОРИИ ОПТИМИЗАЦИИ
Учебное пособие
Санкт-Петербург
2013
УДК 004.94(075)
ББК 32.81я73
А94
Рецензенты:
доктор технических наук, профессор А. А. Гоголь;
доктор технических наук, профессор Л. Л. Полосин
Утверждено
редакционно-издательским советом университета
в качестве учебного пособия
А94
Афанасенко, А. С.
Основы теории оптимизации: учеб. пособие/ А. С. Афанасенко,
Н. А. Обухова, Б. С. Тимофеев. – СПб.: ГУАП, 2013. – 72 с.: ил.
ISBN 978-5-8088-0873-7
Даны теоретические материалы по основным разделам теории
оптимизации. Рассмотрены целевые функции различного вида. Изложены методы поиска минимума функции одной и многих переменных. Обсуждено решение задачи векторной оптимизации на
основе экспертных оценок. Приведены основы регрессионного анализа. Теоретические материалы сопровождаются заданиями на моделирование. Предложены примеры решения оптимизационных задач на основе системы инженерных расчетов Matlab.
Предназначено для студентов радиотехнического факультета,
изучающих курс «Основы оптимизации». Может быть полезно студентам других специальностей и специалистам, применяющим методы оптимизации для решения научно-технических задач.
УДК 004.94(075)
ББК 32.81я73
ISBN 978-5-8088-0873-7
© Санкт-Петербургский государственный
университет аэрокосмического
приборостроения (ГУАП), 2013
© А. С. Афанасенко, Н. А. Обухова,
Б. С. Тимофеев, 2013
1. ЭЛЕМЕНТЫ ТЕОРИИ ОПТИМИЗАЦИИ
1.1. Исследование целевых функций
1.1.1. Основные задачи оптимизации
Функциональные зависимости показателей качества kj от регулируемых параметров i называют целевыми функциями. Задачей
оптимизации является отыскание экстремума (минимума или максимума) целевой функции. В технических задачах целевая функция имеет физический смысл меры нежелательных свойств, которую нужно минимизировать.
В зависимости от числа регулируемых параметров различают
одномерную и многомерную оптимизацию, в зависимости от числа
показателей качества: скалярную (один показатель) и векторную
оптимизацию (несколько показателей качества в каждой точке
пространства параметров). При этом оптимизация может быть выполнена с соблюдением ряда ограничений (условная оптимизация)
или без них (безусловная оптимизация).
Будем рассматривать задачу оптимизации в конечномерном евклидовом пространстве, где целевая функция является скалярной
функцией конечного числа n переменных, составляющих вектор
T
параметров  = [1, 2 ,..., n ] ;  Î R n .
Функция имеет локальный минимум в точке *, если существует такая окрестность этой точки, что для всех значений  в этой
окрестности f(*)< f(). Если это неравенство справедливо для всей
области определения , то имеет место глобальный минимум.
а)
2
1.8
1.6
1.4
1.2
y 1
0.8
0.6
0.4
0.2
0
-1.5
б)
1
0.5
0
-0.5
-1
y
-1.5
-2
D0
D2
D3
D1
-1
-0.5
0
x
-2.5
-3
0.5
1
1.5
-3.5
-2.5
-2
-1.5
-1
x
-0.5
0
0.5
Рис. 1.1.1. Рельеф унимодальной (а) и мультимодальной (б) функций
3
Целевая функция может быть:
унимодальной – функция имеет единственный минимум (рис.
1.1.1, a);
мультимодальной – функция имеет много минимумов, из которых один является глобальным, либо все минимумы равнозначны
(рис. 1.1.1, б).
1.1.2. Численное нахождение минимума
целевой функции
Классический подход к задаче поиска минимума заключается в
поиске уравнения, которому должна удовлетворять точка *. Если
функция f() и ее производные непрерывны, то в точке минимума
производная функции равна нулю f(*)= 0 Это условие необходимое, но недостаточное для минимума. Его выполнение гарантирует
только то, что в этой точке существует перегиб функции (экстремум). Анализ знака второй производной позволяет выяснить, что
это – минимум или максимум (плюс – минимум, минус – максимум). Если вторая производная равна нулю, то ситуация остается
неопределенной.
В общем случае решение нелинейного уравнения f() = 0 не
всегда простое, поэтому используют численный метод Ньютона.
Рассмотрим метод Ньютона для одномерной функции. Мы находимся в точке p с координатой i (рис. 1.1.2). Нужно найти точку
*, которая является решением уравнения. Построим касательную
к кривой f() в точке p. Tочка i+1 будет лучшей аппроксимацией решения, чем i (см. рис. 1.1.2), так как находится ближе к *.
Отношение f()/(I–i+1)=tg=f() определяет вторую производную целевой функции в точке i. Тогда в конце первой итерации
g(D)
P
D
D1
T
M
D0
D
A
Рис. 1.1.2. Решение нелинейного уравнения
4
f ¢( i )
. Это формула Ньютона, которая позволяет итераf ¢¢( i )
ционно получить решение нелинейного в общем случае уравнения.
За исключением простых целевых функций, минимум которых
может быть найден аналитически, оптимизация осуществляется
численными методами в виде последовательности шагов, дающих
приближение к точке минимума:
 i+1 =  i -
1, 2 ,.... k ; k Î [1,¥[; f ( k )  min.
Процесс минимизации сходится, если lim  k = * , и расходитk¥
ся в остальных случаях. Сходимость зависит от применяемого метода и вида целевой функции.
Унимодальная целевая функция может быть удобной или неудобной для реализации процесса минимизации. При формулировке задачи оптимизации важно найти правильное представление
целевой функции и масштабы входящих в нее переменных. Поиск
минимума во многих случаях осуществляют путем обследования
функции в направлениях, определяемых ее частными производными. Разумеется, целевая функция в этих случаях должна иметь эти
производные.
1.1.3. Разложение функции в ряд Тейлора
Формула Тейлора позволяет представить аналитическую функцию f в окрестности некоторой точки в виде бесконечной суммы степенных функций. Коэффициенты ряда Тейлора пропорциональны
значениям производных функции в данной точке. Пусть функция
f() имеет непрерывные первую и вторую производные в окрестности некоторой точки 0. Тогда в этой окрестности она может быть
аппроксимирована первыми тремя членами ряда Тейлора:
f () = f (10 , 20 ,..., n0 ) +
n
¶f,
( - 0i ) +
0 i
¶

i=1
+å
+
1 n n ¶2 f
0
0
0 ( i -  i )( j -  j ) + R3 ,
åå
2 i=1 j=1 ¶ i ¶ j 
где R3 – остаточный член, который стремится к нулю при   0.
5
Введем вектор-градиент
T
é ¶f ¶f
¶f ùú
f () = ê
,
,...,
,
ê ¶ ¶
¶n úû
2
ë 1
составленный из первых частных производных целевой функции
по всем координатным осям. Из анализа известно, что векторный
дифференциал
d = (d1, d2 ,..., dn )T
на поверхности f() = const удовлетворяет соотношению (приращения функции на этой поверхности нет):
n
df
d i =áf (), dñ= 0.
d
i
i=1
df () = å
Вектор-градиент нормален в точке 0 к поверхности f() = const,
так как скалярное произведение равно нулю.
Вещественная матрица Гессе (гессиан)
é d2 f
ù
ú
H() = êê
0 ú , 1 £ i, j £ n.

d

d

êë i j
úû
Значение второй производной не зависит от порядка дифференцирования, поэтому матрица Гессе – вещественная симметричная
матрица.
Используя введенные величины, уравнение разложения в ряд
Тейлора можно переписать в виде:
( ) ( - 0 ) + 12 ( - 0 ) H(0 )( - 0 ) + R3.
f ( ) = f ( ) + f 0
T
T
Представление произвольной нелинейной функции первыми
тремя слагаемыми ряда Тейлора можно рассматривать как ее аппроксимацию квадратичной функцией в окрестности заданной
точки.
1.1.4. Квадратичная целевая функция
и ее особенности
Наиболее распространенным видом аппроксимации целевой
функции является квадратичная форма (рис. 1.1.3), обоснованием
6
f(x,y)
8
7
6
5
4
3
2
1
1
1
0.5
0.5
0
y
–0.5
0
–0.5 x
–1 –1
Рис. 1.1.3. Поверхность квадратичной целевой функции
которой является приведенное выше разложение в ряд Тейлора.
В n-мерном пространстве такая функция описывается формулой
f ( ) = T  + T a + c,
(1.1.1)
T
где  = [1, 2 ,..., n ] – вектор координат точки;  – симметричная, положительно определенная матрица размера nn; a = [a1, a2,
…, an]T – вектор параметров, определяющий положение минимума; с – свободный параметр. В двумерном случае формулу (1.1.1)
можно переписать в виде
é 11
f (x, y) = éë x yùû ê
ê 21
ë
éa ù
12 ù é xù
ú ê ú + é x yù ê 1 ú + c,
ë
û
ê a2 ú
22 úû êë y úû
ë û
где [ x, y ] = [1, 2 ] – вектор координат точки на плоскости.
Градиент квадратичной функции (вектор-столбец, составленный из первых частных производных):
é ¶f ù
ê ú
é211x + 212 y + a1 ù
ê ¶x ú
ú.
f (x, y) = ê ú = 2x + a = ê
ê221x + 222 y + a2 ú
ê ¶f ú
ë
û
ê ú
êë ¶y úû
7
Гессиан (матрица вторых частных производных) не зависит от
значений координат и равен:
é ¶2 f
¶2f ùú
ê
ê
2
¶x¶y úú
ê ¶x
Df (x, y) = ê 2
ú = 2.
ê ¶ f
¶2 f ú
ê
ú
ê ¶y¶x ¶y2 ú
ë
û
Все производные более высоких порядков равны нулю. Поэтому
разложение квадратичной функции в ряд Тейлора содержит всего
три слагаемых. Нахождение минимума квадратичной функции
принято называть линейной оптимизацией или линейным программированием, так как первые частные производные в этом случае линейно зависят от координат.
При известных значениях матрицы  и вектора a минимум квадратичной функции может быть найден аналитически. Для этого
достаточно приравнять к нулю ее градиент:
f (* ) = 2* + a = 0,
откуда
* = -(2-1 a.
(1.1.2)
Рельеф целевой функции определяется матрицей .. Если матрица  диагональная и коэффициенты на диагонали равны между
собой, то линии равных значений (геометрическое место точек, соответствующее f() = const) имеют вид окружностей.
Если матрица  диагональная и коэффициенты на диагонали не
равны друг другу, то линии равных значений имеют вид эллипсов,
вытянутых вдоль одной из координатных осей. Чем больше отношение максимального диагонального коэффициента к минимальному, тем сильнее искажена исходная окружность.
Если матрица  заполнена, то эллипсы будут повернуты относительно координатных осей (рис. 1.1.4).
Степень отличия линий равных значений от окружностей характеризует коэффициент обусловленности матрицы :
cond() =  max /  min ,
где  max и  min – собственные значения матрицы .
Элементы диагональной матрицы  совпадают с ее собственными значениями . В случае, если матрица  заполнена, применяют
более сложную методику определения собственных чисел.
8
а)
1 0¯
°
4 ¡¡
°
¢0 1±
б)
3 0¯
°
4 ¡¡
°
¢ 0 2±
в)
3 1¯
°
4 ¡¡
°
¢1 4±
Рис. 1.1.4. Рельеф квадратичной функции
при различных значениях матрицы 
Пусть А – квадратная матрица n-го порядка и выполняется равенство Av = v, при v  0. Тогда число  называют собственным
значением матрицы А, а ненулевой вектор v – соответствующим
ему собственным вектором. Перепишем задачу в виде
(A - 1)v = 0, v ¹ 0,
(1.1.3)
где 1 – единичная матрица. Для существования нетривиального решения задачи (при v  0) должно выполняться условие
det(A - 1) = 0.
Этот определитель является многочленом n-й степени от ; его
называют характеристическим многочленом. Значит, существует
n собственных значений – корней этого многочлена, среди которых
могут быть одинаковые – кратные.
Если найдено некоторое собственное значение, то, подставляя
его в однородную систему уравнений (1.1.3), определяют соответствующий собственный вектор. Будем нормировать собственные
векторы, т. е. делить на их длину (норму). Собственные векторы
определяют направления главных осей эллипсоида, заданного матрицей .. Собственные векторы линейно независимы и, следовательно, образуют базис. В большинстве случаев собственные векторы ортогональны.
Задачу на собственные значения легко решить для некоторых
видов матриц: диагональной, трехдиагональной, треугольной.
Для диагональной матрицы собственные числа лежат на диагонали. Поэтому для заполненных матриц решение ищут путем
9
приведения исходной матрицы к диагональному виду с помощью
ортогональной матрицы F. Диагональную матрицу G = F–1F называют подобной матрице Q (преобразование подобия). Пусть , v
есть собственное значение и собственный вектор матрицы G, тогда
v = Gv = F–1Fv. Умножим обе части этого уравнения слева на
матрицу F:(Fv) = (Fv). Отсюда видно, что  есть собственное значение матрицы , а Fv – собственный вектор. Преобразование подобия, или приведение матрицы к диагональному виду, не меняет
собственных значений матрицы и по определенному закону преобразует ее собственные вектора. Это свойство позволяет использовать преобразование подобия для расчета искомых величин.
Для двухмерного пространства параметров собственные числа
удовлетворяют квадратному уравнению
11 - 
12
= 0,
21
22 - 
или в развернутом виде
2 - (11 + 22 ) + (1122 - 1221 ) = 0,
причем это уравнение всегда имеет действительные корни.
1.1.5. Функция Розенброка
Функция Розенброка, также называемая «банановой» функцией, используется в качестве тестовой целевой функции для проверки сходимости различных методов минимизации. В двумерном
случае функция имеет вид
f (x, y) = 100(x - y2 )2 + (1 - x)2 .
(1.1.4)
Единственный минимум функции находится в точке с координатами {1,1}. Минимум расположен на дне изогнутого «оврага»
(рис. 1.1.5), что замедляет сходимость большинства методов минимизации при работе с этой функцией.
Задание на моделирование
1. Построить поверхность и линии равных значений квадратичной целевой функции.
2. Рассчитать значения коэффициента обусловленности cond ()
и собственные векторы матрицы .
10
11
Рис. 1.1.5. Поверхность функции Розенброка (а) и ее рельеф вблизи минимума (б)
3. Исследовать зависимость вида рельефа целевой функции от
численных значений параметров: коэффициентов матрицы  (диагональной и заполненной); выяснить влияние вектора a.
4. Получить поверхности и линии равных значений тестовой целевой функции Розенброка.
Контрольные вопросы
1. Что такое целевая функция?
2. Что такое вектор-градиент, каковы его свойства, позволяющие использовать его для поиска минимума целевой функции?
3. Что такое матрица Гессе, каковы ее свойства, позволяющие
использовать ее для поиска минимума целевой функции?
4. Квадратичная форма, роль матрицы  и вектора a.
5. Что такое коэффициент обусловленности cond , как его определить; его влияние на вид линий равного уровня целевой функции.
6. Каковы особенности функции Розенброка, позволяющие использовать ее в качестве тестовой?
1.2. Методы поиска минимума функций
одной переменной
1.2.1. Одномерное представление функции
многих переменных в выбранном направлении
Процедуры систематического получения последовательности
векторов a0, a1, a2,, ak в пространстве параметров Rn, таких, что
выполняются неравенства f(a0 f(a1f(a2 f(ak называют
методами спуска или методами поиска минимума.
В большинстве методов поиска минимума выбор следующей
точки приближения предполагает:
определение направления поиска;
нахождение локального минимума целевой функции по заданному направлению.
При этом на каждой итерации формируется одномерная целевая функция f(x), являющаяся проекцией n-мерной функции F(),
 = [1, 2, …, n]T на найденное направление поиска:
f (x) = F (0 + xd ),
(1.2.1)
где F(…) – исходная функция векторного аргумента; x – скалярная
переменная; 0 – начальная точка поиска; d = [d1, d2, …, dn]T – единичное направление поиска.
12
Если у функции F() существуют первые частные производные
по всем компонентам вектора , ее градиент также может быть использован при одномерном поиске. Тогда производная одномерной
функции f(x) будет проекцией градиента на направление поиска:
df
=< d,F (0 + xd ) >,
dx
(1.2.2)
где угловые скобки обозначают скалярное произведение векторов.
Следует различать случаи, когда градиент целевой функции действительно известен, и те, когда он аппроксимирован численными
методами. В простейшем случае
F ( ) » ( F ( + he) - F ( )) h ,
где h – малая положительная величина; e – единичный вектор. При
замене производных конечными разностями методы поиска минимума могут утрачивать сходимость из-за неустойчивости численного дифференцирования.
1.2.2. Общая структура методов одномерного поиска
Все рассматриваемые в данной работе методы одномерного поиска имеют схожую структуру (рис. 1.2.1).
1. Путем предварительного анализа значений функции локализуют отрезок, внутри которого находится минимум.
2. Выполняют итерационное уточнение положения минимума.
3. Поиск завершают при выполнении одного или нескольких
требований, характеризующих точность определения координаты
минимума, значения функции в точке минимума, количество выполненных итераций и пр.
Методы одномерной минимизации могу быть разделены на две
группы:
методы поиска, основанные на разбиении интервала неопределенности;
методы поиска, основанные на интерполяции целевой функции.
Методы первой группы делят интервал неопределенности по
некоторому правилу и отбрасывают ту часть интервала, в которой
функция принимает наибольшие значения. Методы второй группы
выполняют интерполяцию целевой функции на интервале неопределенности, используя небольшое число измерений, и ищут минимум интерполирующей функции (наиболее широко распространена интерполяция с помощью полинома 2-й или 3-й степени).
13
Исходные данные: целевая
функция, начальная точка
и направление поиска
Одномерный поиск
минимума
Определение начального
интервала неопределенности
На каждой итерации
последовательно уточняется
положение точки минимума
Вычисление положения точки
минимума
Уточнение границ интервала
неопределенности
Да
Первая итерация?
Нет
Проверка критериев
завершения поиска
Достигнута ли заданная
точность?
Нет
Да
Превышено ли допустимое
число итераций?
Нет
Да
Конец
Рис. 1.2.1. Обобщенная структура методов одномерного поиска
1.2.3. Выбор начального интервала неопределенности
Начальный интервал неопределенности, от которого зависит
успешность всей процедуры, выбирают путем грубого поиска. На
первом этапе поиска нужно задать интервал поиска a < x < b, внутри которого находится минимум целевой функции. После окончания предыдущей итерации получена точка xk., с которой начинают
процедуру грубого поиска (рис. 1.2.2). Из точки xk делают шаг h в
направлении d (или –d, как заранее определено). Если шаг неудачный, направление поиска изменяют на обратное. Если шаг удачный, делают удвоенный шаг 2h в том же направлении. Когда целевая функция стала возрастать, целесообразно произвести реверс и
14
f (x)
a
-h +h
0
b
+2h
+4h
-2h
x
Рис. 1.2.2. Определение начального интервала поиска
вернуться обратно на половину длины последнего шага. Минимизация произойдет успешно, если минимум будет заключен между
крайними точками.
1.2.4. Методы поиска, основанные
на разбиении интервала неопределенности
Бисекция
Метод бисекции известен как простейший численный метод решения нелинейных уравнений вида f(x) = 0. Он может быть применен для поиска минимума в том случае, когда целевая функция на
интервале неопределенности непрерывно дифференцируема (рис.
1.2.3), хотя метод может сходиться и при наличии у функции разрывов первого рода.
В начале поиска задается интервал [a, b], на границах которого
производная принимает значения противоположных знаков. Выполнение этого условия для функции с монотонной производной означает существование единственного экстремума в пределах интервала неопределенности. На каждой итерации измеряется значение
производной в середине интервала f((a+b)/2). Если знак f((a+b)/2)
не равен знаку f(a), новым интервалом является [a, (a+b)/2]. В противном случае новым интервалом является [(a+b)/2, b]. Поиск завершается тогда, когда ширина интервала становится меньше удвоенной величины погрешности x, либо абсолютное значение производной становится меньше заданной величины f(x). За N итераций
интервал неопределенности сужается в 2N раз, при этом требуется
N+2 раза вычислить значение производной целевой функции.
15
f’(x)
f ’((a+b )/ )2
f ’ (b )
a
(a+b )/2
b
x
f ’(a)
a*
b*
Рис. 1.2.3. Поиск минимума методом бисекции
Дихотомия
В методе дихотомии так же, как в методе бисекции, интервал
неопределенности сужается в 2 раза за одну итерацию, но не требуется вычисления производной.
Границы начального интервала неопределенности задаются
точками a, b (рис. 1.2.4). Значения целевой функции измеряются
в точках x1…x5, где
x1 = a,
1
x2 = a + (b - a),
4
(a + b )
x3 =
,
2
3
x4 = a + (b - a),
4
x5 = b.
Среди измеренных в указанных точках значений f(x1) … f(x5)
выбирается минимальное. Обозначив за i индекс точки, в которой
определено минимальное значение f(x), оставим только ту часть
исходного интервала, которая находится в окрестности найденной
точки минимума xi. Новый интервал неопределенности [xi–1, xi+1],
имеет длину в 2 раза меньше начальной длины. На приведенном
рисунке новым интервалом является [x2, x4]. С учетом этого на следующей итерации происходит переобозначение точек:
16
x1 = xi-1 = a* ,
x3 = xi ,
x5 = xi+1 = b* ,
1 *
b - a* ,
4
3
x4 = a* + b* - a* .
4
Значение целевой функции требуется измерять только для точек x2 и x4, так как остальные значения известны с предыдущей
итерации. Процесс сужения интервала неопределенности продолжается до тех пор, пока его ширина не станет меньше 2x.
Метод дихотомии имеет достаточно высокую скорость сходимости (за N итераций интервал неопределенности сокращается в 2N
раз), не использует производной, но при этом требует 2N+3 измерений целевой функции, что для многих задач делает его менее выгодным, чем метод «золотого сечения»
x2 = a* +
(
(
)
)
Метод «золотого сечения»
Правилом «золотого сечения» в искусстве принято называть
разделение некоторой величины на две части, при котором отношение меньшей части к большей и отношение большей части ко всей
величине равны между собой:
M
m
=
= , m < M,
M +m M
где «золотое число»   0,618.
e (w3 ) e (w2 ) e (w4 ) e (w5 ) e (w1)
f (x)
e (w1)
новый интервал
e (w4 )
e (w2 )
a
e (w5 )
e (w3 )
( a+ b)/2
b
x
Рис. 1.2.4. Поиск минимума методом дихотомии
17
При поиске минимума методом золотого сечения (рис. 1.2.5) на
интервале неопределенности выбираются четыре точки x1…x4, где
x1 = a,
x2 = b - 0,618(b - a),
x3 = a + 0,618(b - a),
x4 = b.
Среди измеренных в указанных точках значений f(x1) … f(x4)
выбирается минимальное. Очевидно, что минимум функции располагается на одном из двух прилегающих к этой точке отрезков,
поэтому интервалом неопределенности становится [a*, b*]. Оставшийся отрезок необходимо отбросить. Обозначив за i индекс точки,
в которой определено минимальное значение f(x), определим новые
точки x*1…x*4:
ìïx* = x = a* ,
ïï 1
1
ïï *
*
ïïx4 = x3 = b ,
åñëè i = 1 èëè i = 2,
í *
ïïx3 = x2 ,
ïï
ïïx* = b* - 0,618 b* - a* ,
ïîï 2
(
)
ìïx* = x = a* ,
ïï 1
2
ïï *
*
ïïx4 = x4 = b ,
åñëè i = 3 èëè i = 4.
í *
ïïx2 = x3 ,
ïï
ïïx* = a* + 0,618 b* - a* ,
ïïî 3
(
)
новый интервал
f (x)
e (w4 )
e (w1)
e (w3 )
e (w2 )
0,618 (b-a)
0,618 (b-a)
b
w1* w1 þ
w2*
w3*
w4*
x
Рис. 1.2.5. Поиск минимума методом золотого сечения
18
На следующей итерации появляется только одна новая точка, в
которой необходимо измерить значение целевой функции. Эта недостающая точка располагается относительно середины интервала
симметрично имеющейся внутренней точки. Новые точки x*1…x*4
расположены также в пропорции золотого сечения.
Несмотря на то, что интервал неопределенности сокращается
только в 1,618 раза за одну итерацию, поиск требует лишь N+3 измерения целевой функции, что в большинстве случаев делает его
более быстродействующим, чем метод дихотомии.
Метод Фибоначчи
Метод Фибоначчи является приближенным аналогом метода
«золотого сечения», но в нем исключены арифметические операции
с плавающей запятой, что может быть выигрышным для реализации этого метода на некоторых вычислительных архитектурах.
Особенностью метода Фибоначчи является необходимость предварительного определения числа точек, в которых вычисляются значения целевой функции. Каждая пара точек, в которых вычисляют
целевую функцию, располагается симметрично относительно центра интервала неопределенности, причем координаты этих точек
определяются из выражений:
Ô
Ôk
x2 = k-1 (x4 - x1 ) + x1, x3 =
(x4 - x1 ) + x1,
Ôk+1
Ôk+1
где Фk – числа Фибоначчи: Ф0 = Ф1 = 1; Фk+1 = Фk+Фk–1.
В методе Фибоначчи на каждой итерации вычисляют только
одно значение целевой функции, тогда как в методе дихотомии –
два. Интервал неопределенности с каждым шагом сокращается в
пределах от 1,6 до 2 раз.
1.2.5. Методы поиска,
интерполирующие целевую функцию
Квадратичная интерполяция
Метод применим в тех случаях, когда целевую функцию на интервале поиска можно аппроксимировать полиномом 2-й степени:
(x) = ax2 + bx + c.
Чтобы определить координату минимума аппроксимирующей
функции (x), нужно продифференцировать это выражение по x и
приравнять производную нулю:
19
d
b


= 2ax + b = 0, îòêóäà x = - .
dx
2a
Неизвестными в этом уравнении являются коэффициенты аппроксимирующего полинома a и b. Для того чтобы их найти, необходимо задать три точки на оси x (рис. 1.2.6) и решить систему
уравнений:
ì
ï
(x1 ) = ax12 + bx1 + c;
ï
ï
ï
ï(x ) = ax2 + bx + c;
í
2
2
2
ï
ï
2
ï
(x3 ) = ax3 + bx3 + c.
ï
ï
î
(1.2.3)
Результат измерения целевой функции в точках x1, x2 и x3 подставляется в левую часть (аппроксимирующий полином в этих точках точно совпадает с целевой функцией) и система принимает вид
é 2
é f (x1 ) ù ê x1
ê
ú ê 2
ê
ú ê
ê f (x2 )ú = ê x2
ê f (x )ú ê 2
êë 3 úû ê x3
ë
x1 1ùú é aù
ú ê ú
x2 1ú * êê b úú , èëè f = Ra.
ú
ê ú
x3 1úú êë c úû
û
(1.2.4)
Наиболее известным способом решения системы линейных
уравнений, подобной (1.2.3), является умножение на обратную
e (w), M(w)
f (x)
M(w)
w1
w2
ŵ
(w3 – w1)
j
w3
x
новый интервал
Рис. 1.2.6. Поиск минимума методом квадратичной интерполяции
20
матрицу: a = R–1f. Найденные коэффициенты полинома a = [a, b,
c]T подставляют в (1.2.3) и получают координату экстремума аппроксимирующей функции. Эта точка является минимумом лишь
в том случае, когда a > 0 (то есть вторая производная (x) положительна). В противном случае аппроксимация является неудачной
и следует выбрать другие точки x1, x2 и x3. Так как размерность
системы (1.2.4) невысока, решение можно записать в явном виде,
без обращения матрицы:
xˆ =
(x22 - x32 )f (x1 ) + (x32 - x12 )f (x2 ) + (x12 - x22 )f (x3 ) .
2 éë(x2 - x3 )f (x1 ) + (x3 - x1 )f (x2 ) + (x1 - x2 )f (x3 )ùû
(1.2.5)
Еще проще выглядит решение в случае равноотстоящих друг от
друга точек x1 < x2 < x3. Путем замены переменной в (5) получаем:
xˆ = x2 +
f (x1 ) - f (x3 )
h
*
, h = x2 - x1 = x3 - x2 . (1.2.6)
2 f (x1 ) - 2f (x2 ) + f (x3 )
Для того чтобы квадратичная аппроксимация позволяла найти
минимум, необходимо, чтобы аппроксимирующая функция (x)
была строго выпуклой вниз. Это требование эквивалентно требованию a > 0 в (1.2.3) и требованию положительного знака знаменателя в (1.2.6). Найденная точка минимума x̂ отличается от истинного минимума целевой функции. На следующей итерации исходный
интервал неопределенности [x1, x3] сокращается в k раз и аппроксимация повторяется с использованием нового набора точек. Обычно выбирают k в диапазоне 4…9.
Кубическая интерполяция без использования производных
Метод поиска минимума, когда целевую функцию аппроксимируют полиномом 3-й степени, несколько более сложен в реализации по сравнению с квадратичной интерполяцией, но сходится для
более широкого класса целевых функций. Аппроксимирующая
функция в этом случае имеет 4 коэффициента:
(x) = ax3 + bx2 + cx + d.
Для нахождения минимума (x), аналогично (3), приравнивают
к нулю первую производную по x:
d
= 3axˆ 2 + 2bxˆ + c = 0.
dx
21
Если данное квадратное уравнение имеет действительные корни, они являются точками, подозрительными на экстремум:
xˆ1,2 =
-2b  (2b)2 -12ac
.
6a
(1.2.7)
Чтобы выяснить, какое из решений является точкой минимума,
проверяют знак второй производной:
d2 
dx2
= 6ax + 2b,
(1.2.8)
которая в точке минимума должна быть положительной. Для нахождения четырех неизвестных коэффициентов полинома a, b, c, d
необходимо решить систему из четырех уравнений:
ì
ï
(x1 ) = ax13 + bx12 + cx1 + d;
ï
ï
ï
ï
(x2 ) = ax23 + bx22 + cx2 + d;
ï
ï
í
ï
(x3 ) = ax33 + bx32 + cx3 + d;
ï
ï
ï
ï
3
2
ï
ï
î(x4 ) = ax4 + bx4 + cx4 + d.
(1.2.9)
В выбранных точках x1...x4 значения аппроксимирующего полинома в точности совпадают со значениями целевой функции:
(xi) = f(xi). Положение точек x2 и x3 в пределах интервала неопределенности может быть произвольным, но удобнее распределить
точки равномерно (рис. 1.2.7). После подстановки значений f(xi) и
перевода в матричную форму система (1.2.9) принимает вид
é 3
é f (x1 ) ù ê x1
ê
ú ê 3
ê f (x2 )ú ê x2
ê
ú ê
ê f (x ) ú = ê 3
ê 3 ú ê x3
ê f (x )ú ê
ëê 4 ûú ê x43
ë
x12
x22
x32
x42
x1 1ùú é a ù
ú ê ú
x2 1ú ê b ú
ú * ê ú , èëè f = Ra.
ê ú
x3 1úú ê c ú
ú êd ú
x4 1úû ëê ûú
(1.2.10)
Одним из способов решения системы (1.2.10) является умножение на обратную матрицу: a = R–1f. Но обращение матрицы – достаточно медленная и неустойчивая операция, поэтому на практике системы линейных уравнений решают, используя QR- или SVDразложение матрицы R. Найденные коэффициенты полинома a =
[a, b, c, d]T подставляют в уравнение (1.2.7) для нахождения экстремумов, а затем согласно (1.2.8) выбирают минимум (при усло22
e (w), M(w)
f (x)
M(w)
w1
ŵ w2
w3
(w4 – w1)
j
w
новый интервал
Рис. 1.2.7. Поиск минимума методом кубической интерполяции
вии его существования). Полученное значение x̂ является текущей
оценкой минимума функции и переходит на следующую итерацию
в качестве центра интервала неопределенности. Ширина интервала
при этом сокращается в k раз (обычно k выбирают в диапазоне 4…9).
На рис. 1.2.7 показана целевая функция f(x), измеренная в
точках x1…x4. Аппроксимирующий полином (x) имеет минимум
в точке xˆ, который отличается от истинного минимума целевой
функции. На следующей итерации исходный интервал неопределенности [x1, x4] сокращается в k раз и аппроксимация повторяется с использованием нового набора точек.
Кубическая интерполяция с использованием производных
Если для измерения доступна не только целевая функция, но
и ее производная, то аппроксимирующий полином может быть
построен c использованием значений целевой функции на краях
интервала [x1, x2] и производных, измеренных в этих же точках.
Удобно в интерполирующей функции осуществить замену переменной x на , изменяющуюся от 0 до 1 на заданном интервале:
() = a3 + b2 + c + d,  = (x - x1 ) (x2 - x1 ).
На краях интервала значения полинома () и его производной
определяются целевой функцией:
() = f ((x - x1 ) (x2 - x1 )),
¢() = f ¢(x) (x2 - x1 ).
23
Для нахождения коэффициентов интерполирующего полинома
необходимо решить систему уравнений:
ìï(0) = f (x1 ) = d;
ïï
ïï¢(0) = f ¢(x ) (x - x ) = c;
1
2
1
ï
í
ïï(1) = f (x2 ) = a + b + c + d;
ïï
ïï¢(1) = f ¢(x2 ) (x2 - x1 ) = 3a + 2b + c.
î
Полученная система легко решается методом подстановки и коэффициенты полинома записываются в явном виде:
ì
ï
a = (f ¢(x1 ) + f ¢(x2 ))(x2 - x1 ) - 2(f (x2 ) - f (x1 ));
ï
ï
ï
ï
b = 3(f (x2 ) - f (x1 )) - (2f ¢(x1 ) + f ¢(x2 ))(x2 - x1 );
ï
ï
í
ï
c = f ¢(x1 )(x2 - x1 );
ï
ï
ï
ï
ï
ï
îd = f (x1 ).
После того, как вычислены коэффициенты полинома, точку минимума opt находят таким же образом, как и в случае кубической
интерполяции по четырем точкам. С учетом обратной замены переменной:
xopt = x1 + (x2 - x1 )opt .
Достаточным условием сходимости описанного метода является
строгая выпуклость функции f(x) вниз на заданном интервале, то
есть выполнение условий f(a)<0, f(b)>0 (рис. 1.2.8, а). Однако существование минимума возможно и в других случаях (рис. 1.2.8, б).
а)
e a þ 0
б)
e a a
0
a
b
e a þ 0
в)
e a a
0
a
b
e a þ
0
г)
e a a
0
a
e a þ 0
e a a
0
b
b
a
Рис. 1.2.8. Примеры существования (а, б) и отсутствия (в, г)
минимума кубической параболы в зависимости от знаков производных
на концах интервала
24
1.2.6. Критерии остановки счета
Если целевая функция не квадратичная, то любым из предложенных методов минимум за конечное число шагов может быть
найден лишь приближенно. Чтобы выполнить поиск минимума за
конечное число шагов, счет останавливают согласно заранее установленному критерию. Различают критериальные и аргументные
задачи оптимизации. В первом случае важно лишь обеспечить минимальное значение целевой функции. Во втором случае надо точно найти координаты минимума.
Критерий остановки счета в первом случае: f (xk+1 ) ³ f (xk ), т. е.
счет прекращают, если целевая функция перестает уменьшаться. Во втором случае критерий (аргументный) остановки счета
xk+1 - xk £ , где  – малая заданная величина.
Для методов, использующих значение первой производной, счет
может быть прекращен, если абсолютное значение производной
приблизилось к нулю: f ¢(xk ) < f ¢(x).
Дополнительными критериями остановки счета могут быть:
превышение допустимого числа итераций: i > imax;
достижение минимальным значение целевой функции заданного предела: fk < fmin .
Задание на моделирование
Разработать и написать программу, реализующую:
нахождение минимума целевой функции методом одномерной
минимизации (вид функции и метод указаны в приложении);
вывод на экран графика целевой функции, а также результатов
поиска минимума на каждой итерации.
Контрольные вопросы
1. Поясните, каким образом сокращают интервал неопределенности на каждой итерации в методах линейного поиска минимума
целевой функции?
2. Какой из методов имеет большую скорость сходимости – дихотомии или золотого сечения?
2. Какие методы одномерного поиска используют аппроксимацию целевой функции?
3. Приведите алгоритмы методов квадратичной и кубической
интерполяции.
25
1.3. Методы поиска минимума функций
многих переменных
1.3.1. Классификация методов поиска минимума функций
многих переменных
Методы минимизации функции n переменных делят на три
группы. В зависимости от того, какая информация о целевой функции используется, различают:
методы нулевого порядка (методы прямого поиска), использующие только значения целевой функции;
методы первого порядка, использующие значения целевой
функции и ее первых производных;
методы второго порядка, использующие значения целевой
функции, ее первых и вторых частных производных. Использование вторых частных производных (матрицы Гессе) означает, что
методы данной группы пытаются локально аппроксимировать целевую функцию в квадратичной форме. На практике из методов
этой группы применяют в основном квазиньютоновские методы,
заменяющие матрицу Гессе ее аппроксимацией, и, строго говоря,
относящиеся к методам первого порядка.
Разработка значительного числа разнообразных методов поиска
минимума обусловлено тем, что наилучшего метода, пригодного во
всех случаях практического применения, не существует.
К факторам, которые следует учитывать при выборе метода поиска минимума целевой функции можно отнести: число итераций;
число арифметических операций (умножений или делений), выполняемых на каждой итерации; точность при заданном критерии
остановки счета; размерность пространства параметров. Так как в
подавляющем большинстве практических случаев целевая функция задана численно, быстродействие метода, в основном, определяется числом обращений к датчику целевой функции, а также
числом выполненных арифметических операций.
1.3.2. Методы нулевого порядка
Симплекс-метод Нелдера–Мида относится к методам прямого
поиска, использующим только значения целевой функции. Этот
метод, хотя и не обеспечивает высокой скорости сходимости, широко применяется на практике в тех случаях, когда нет априорных
сведений о рельефе целевой функции, либо ее производные недо26
ступны, либо требуется ценой небольших вычислительных затрат
получить решение, хоть сколько-нибудь лучшее, чем начальное
приближение.
Симплекс в n-мерном пространстве задается набором из n+1
вершины {1, 2 ,..., n+1 }. На каждой итерации симплекс модифицируется таким образом, чтобы приблизиться к точке локального
минимума. Один из алгоритмов, реализующих симплекс-метод,
описывается следующей последовательностью шагов.
1. Вершины симплекса сортируют в порядке возрастания значений целевой функции:
f (1 ) £ f (2 ) £  £ f (n+1 ).
2. Операции отражения, растяжения и сжатия производят путем переноса «наихудшей» вершины n+1 в сторону центра тяжести симплекса. Центр тяжести рассчитывают по всем вершинам
симплекса, кроме наихудшей:
0 =
1 n
å i .
n 1
3. Вычисляют точку отражения по формуле
r = 0 + r (0 -  n+1 ),
где стандартным является значение r = 1. Если в найденной точке
значение функции меньше, чем в следующей наихудшей точке n ,
но не меньше, чем в наилучшей точке 1:
f (1 ) £ f (r ) < f (n ),
то вершину n+1 в симплексе заменяют на r и выполняют переход к шагу 1.
4. Если значение целевой функции f (r ) > f (1 ), то переходят к
шагу 5. В противном случае вычисляют точку растяжения по формуле
 e = 0 +  (0 - n+1 ),
где стандартным является значение  = 2. Если f ( l ) < f (r ), то вершину n+1 в симплексе заменяют на  l , в противном случае заменяют n+1 на r . Переходят к шагу 1.
5. На данном шаге f (r ) ³ f (n ). Вычисляют точку сжатия по
формуле
27
 c = 0 + (0 - n+1 ),
где стандартным является значение  = –0,5. Если новая точка приводит к уменьшению целевой функции: f ( c ) < f (n+1 ), то вершину
n+1 в симплексе заменяют на  c , и переходят к шагу 1.
6. Переход на данный шаг, как правило, означает, что точка минимума оказалась внутри фигуры, образованной вершинами симплекса. Производят уменьшение размеров симплекса по формуле
 i =  i + ( i - 1 ) äëÿ âñåõ i Î [2,n + 1],
где стандартным является значение  = 0,5. Вычисляют новые значения целевой функции в точках 2...n+1 и переходят к шагу 1.
Поиск завершается при выполнении одного из условий:
максимальное и минимальное значения функции в вершинах
симплекса отличаются не более, чем на заданную величину:
f (n+1 ) – f (1 ) < f ;
размеры симплекса уменьшились до заданного предела:
max  i - 1 <  x , i Î [2,n + 1];
превышено максимально допустимое число итераций.
Критическое влияние на сходимость метода оказывают начальное положение и начальный размер симплекса. Как правило, для
их задания используют некоторые априорные данные о решаемой
задаче.
Графическое пояснение операций отражения, растяжения, сжатия и уменьшения симплекса для двумерного случая приведено
на рис. 1.3.1. Полная схема описанного алгоритма приведена на
рис. 1.3.2.
1.3.3. Методы первого порядка
Градиентные методы
Градиентным называют метод спуска, в котором следующее
значение вектора параметров k+1 выбирают по отношению к предыдущему k в направлении, противоположном вектору-градиенту
f ( k ) = (df / d1,df / d2 ,...,df / dn )Ò , составленному из первых
частных производных целевой функции. Вектор-градиент f ( k )
показывает направление наибыстрейшего роста функции f() из
точки k. Если поделить вектор-градиент на его длину (норму),
28
а)
x0
xc
x3
б)
x1
x1
xe
xr
x0
x2
в)
xr
x3
x2
г)
x1
x1
xe
x0
x3
x3
x0
xc
x2
x2
д)
x1
w1 V(w3 ”w1)
w1 V(w2 – w1)
x3
x2
Рис. 1.3.1. Операции над двумерным симплексом:
положение промежуточных точек (а), отражение (б),
растяжение (в), сжатие (г) и уменьшение (д)
то получим единичный вектор-градиент, указывающий лишь направление
rk = f ( k ) / f ( k ) .
При минимизации целевой функции градиентными методами
следующая k+1 точка при спуске  k+1 =  k -  k rk , где k – величина шага в направлении – rk, противоположном направлению вектора-градиента. Варианты градиентного метода различаются выбором метрики в пространстве Rn и величины шага .
В методе Гаусса–Зейделя движение к точке минимума происходит по координатным осям (прямоугольной сетке). Направление
поиска совпадает с той координатной осью, вдоль которой абсолютная величина частной производной df / d j максимальна. Следующая точка при спуске
 k+1 =  k - k sign[df ( k ) / d j ]; k > 0;
df ( k ) / d j > df ( k ) / d i , j ¹ i,
где символ sign означает «знак».
29
30
Изменение симплекса
Нет
Размеры симплекса
меньше допустимых?
Нет
Превышено
допустимое
число итераций?
Нет
fn+1 – f1 < tol?
f1 <= f2 <= ... <= fn+1
Сортировка вершин
по возрастанию ц. ф.
Проверка критериев
завершения
Задание начальных
координат x1...xn+1
Симплекс-метод
Проверка критериев
завершения
Растяжение
Да
fe < fr?
fr <= f1?
Да
Вычисление fe
Нет
Нет
fr < fn?
Да
Нет
Отражение
Сжатие
Нет
fc < fn+1?
Вычисление fc
Уменьшение
Вычисление f2…fn+1
Да
Конец
Вычисление центра
тяжести
Вычисление fr
Завершение поиска
Изменение симплекса
Рис. 1.3.2. Схема симплекс-метода в нотации «Дракон»
Завершение поиска
Да
Да
Да
Симплекс состоит из n+1
вершины (n - размерность
пространства параметров)
В методе наискорейшего спуска величину шага k определяют после решения задачи одномерной оптимизации по направлению – rk:
f ( k - k rk ) = min;
 k > 0.
В случае квадратичной функции это приводит к тому, что поиск
в n-мерном пространстве осуществляется по n взаимно-перпендикулярным направлениям (рис. 1.3.3).
Устойчивость и сходимость градиентных методов при квадратичной целевой функции f () =á, ñ+áà, ñ+ c зависят от матрицы . Если матрица положительно определена (все собственные
значения больше нуля), градиентные методы сходятся к стационарной точке, где вектор-градиент равен нулю. Скорость сходимости градиентных методов зависит от коэффициента обусловленности cond
 Главным недостатком градиентных методов считают
их медленную сходимость для плохо обусловленных задач, когда
cond >> 1.
Метод сопряженных градиентов
Медленную сходимость градиентных методов преодолевают за
счет использования для выбора следующего направления поиска информации, полученной на предыдущей итерации. В методе
сопряженных градиентов Флетчера–Ривса направление следующего поиска определено вектором dk,  -сопряженным по отношению к предыдущему направлению поиска dk–1. Два вектора dk
и dk–1 называют  -сопряженными, если скалярное произведение
ád k , d k-1 ñ= 0. Новое направление поиска является линейной ком-
y
1
0
-1
-3
x0
xk
-2
x
-1
0
Рис. 1.3.3. Траектория поиска минимума квадратичной функции
методом наискорейшего спуска
31
бинацией старого направления dk–1 и вектора-градиента в новой
точке rk
d k = -rk + k-1d k-1,
(1.3.1)
где скаляр k–1 выбирают таким, чтобы ád k , d k-1 ñ= 0.
Для расчета k-1 образуем скалярное произведение
0 = dk , dk-1 = rk , dk-1 + k-1 dk-1, dk-1 .
Тогда
k-1 =
rk , d k-1
d k-1, d k-1
.
(1.3.2)
В подавляющем большинстве случаев целевая функция задана
численно и матрица  неизвестна. Чтобы исключить эту матрицу
воспользуемся следующим равенством:
d k-1 = ( k –  k-1 ) = (rk – rk-1 ).
Тогда выражение (1.3.2) можно записать в следующем виде:
k-1 =
rk , rk – rk-1
dk-1, rk – rk-1
=
rk , rk – rk , rk-1
dk-1, rk – dk-1, rk-1
.
Но векторы-градиенты на соседних шагах всегда ортогональны
árk , rk-1 ñ= 0. Член в знаменателе ád k-1, rk ñ также равен нулю, так
как новый вектор-градиент всегда ортогонален предыдущему направлению поиска. Теперь нужно исключить d k-1 из знаменателя.
Умножим слева на rk использованное ранее равенство
d k = -rk + k-1dr -1.
Тогда
árk ,d k ñ= -árk , rk ñ+ k árk ,d k-1 ñ.
Последний член этой формулы равен нулю, так как новый вектор градиент rk всегда ортогонален предыдущему направлению поиска d k-1. Следовательно
árk ,d k ñ= -árk , rk ñ.
Равенство справедливое на k-м шаге, справедливо и для k-1 шага,
значит
árk-1,d k-1 ñ= -árk-1, rk-1 ñ.
32
С учетом этих равенств
k-1 =
árk , rk ñ
.
árk-1, rk-1 ñ
(1.3.3)
Оказывается, чтобы найти новое направление поиска нужно
лишь вычислить новый вектор градиент, умножить его скалярно
самого на себя и подставить в числитель формулы (1.3.3). Заметим,
что скалярное произведение, стоящее в знаменателе, уже найдено
на предыдущем шаге.
Для расчета k-1 может быть использована формула Полака–
Рибьера:
T
æ
ö
ç < r ,(rk - rk-1 ) >÷÷
k-1 = maxçç0, k
÷÷.
çç
< rkT-1, rk-1 > ø÷
è
(1.3.4)
Алгоритм нахождения минимума методом сопряженных градиентов состоит в следующей последовательности шагов.
1. Задают начальную точку в n-мерном пространстве.
2. Определяют направление поиска. На первой итерации и далее
через каждые n итераций поиск ведется в направлении, обратном
градиенту. На остальных итерациях направление вычисляется по
формуле (1.3.1), причем значение k-1 может быть найдено по
(1.3.3) или (1.3.4).
3. Выполняют одномерный поиск по выбранному направлению.
В результате получают значение оптимального шага, дающее новую оценку точки минимума:  k+1 =  k + k d k .
4. В найденной точке вычисляют значение целевой функции и ее
градиента. Проверяют критерии завершения поиска. Если один из
критериев выполняется, завершают поиск в текущей точке. В противном случае возвращаются к шагу 2.
Метод сопряженных градиентов позволяет найти минимум квадратичной целевой функции за n итераций, где n – размерность
пространства параметров.
1.3.4. Методы второго порядка
Метод Ньютона
В методе Ньютона последовательность точек спуска
1
 k+1 =  k - k Hk ( k )f ( k ),
33
é d2 f ù
ú , 1 £ i, j £ n – матрица вторых частных прогде H( k ) = êê
ú


d
d
ëê ki kj ûú
изводных целевой функции (матрица Гессе), взятых в точке  k .
Применение метода Ньютона для поиска минимума квадратичной целевой функции позволяет ограничиться одним шагом. Однако этому методу присущи серьезные недостатки: при неквадратичной целевой функции найти минимум за один шаг не удается, а чтобы выработать новое направление поиска надо заново заполнить и
обратить матрицу Гессе. При значительной размерности пространства параметров на эти операции затрачивается много времени.
Квазиньютоновские методы
Метод Давидона-Пауэлла
Для преодоления этого и других недостатков метода Ньютона его
модифицируют. Матрицу вторых производных H-1
k ( k ) заменяют
на ее аппроксимацию Hk , которую затем уточняют рекуррентным
способом на основе информации, полученной на каждом последу1
ющем шаге итерационного процесса, так что Hk - Hk ()  0 при
k  ¥. Этим приемом устраняют операции обращения матрицы
Hk ( k ), а также существенные затраты времени на определение
численными методами n2 / 2 вторых частных производных целевой функции.
Направление поиска на текущей итерации определяют решением матричного уравнения по формуле, аналогичной формуле Ньютона:
Hk d k = -f ( k ),
(1.3.5)
где dk – направление поиска на k-й итерации; Hk – текущая аппроксимация матрицы Гессе.
На первом шаге матрицу H0 принимают равной единичной; первое направление поиска из начальной точки 0 – противоположным вектору-градиенту – r0 . Далее на каждом шаге вычисляют
вектор-градиент и обновляют матрицу Hk .
Следующее приближение матрицы
Hk+1 = Hk + k
T
d k dT
k - Hk p k p k Hk ,
ád k , pk ñ ápk , Hk pk ñ
где вектор pk = f ( k+1 ) -f ( k ).
34
На практике желательно избегать обращения матрицы при решении (1.3.5): эта операция не только сложна в вычислительном
отношении (объем вычислений растет пропорционально третьей
степени размерности пространства параметров), но и может дать
неточный результат для плохо обусловленной матрицы. Поэтому
на шаге 4 сразу вычисляют матрицу, обратную Hk (формула Шермана–Моррисона):
-1
1
Hk+1 = Hk +
(sTk pk + pTk H-k 1pk )(sksTk ) - H-k 1pksTk + skpTk H-k 1 , (1.3.6)
2
sT
k pk
(sTk pk )
где sk = k d k .
Несмотря на кажущуюся сложность, последняя формула предпочтительнее, чем обращение матрицы Hk ( k ), по соображениям
числовой устойчивости, хотя при небольшой размерности пространства параметров вычисление (1. 3. 6) более затратное.
Методы, основанные на замене матрицы H-1
k ( k ) на ее аппроксимацию Hk , получили название квазиньютоновских и должны
быть отнесены к градиентным методам, так как не требуют вычисления вторых частных производных целевой функции. Это
весьма мощная процедура, очень эффективная при оптимизации
большинства функций независимо от того, квадратичные они или
нет. Скорость сходимости этих методов несколько выше, чем у метода сопряженных градиентов Флетчера–Ривса. Главный их недостаток заключается в необходимости хранить и пересчитывать
матрицу Hk размерностью n ´n, что при больших n требует значительного объема памяти и выполнения многих арифметических
операций.
Задание на моделирование
1. Реализовать методы многомерного поиска:
симплекс-метод;
метод наискорейшего спуска;
метод Гаусса–Зейделя;
метод сопряженных градиентов;
квазиньютоновский метод.
2. Протестировать указанные методы на наборе функции, определенном вариантом задания (таблица).
3. Сведите результаты тестирования в таблицу.
35
Вид
Метод Начальная Число
Точка
Значение
Число
целевой поиска
точка
итераций минимума функции обращений
функции
в точке к датчику
минимума
Контрольные вопросы
1. Дайте классификацию методов поиска минимума для функций многих переменных.
2. Как выбирается направление поиска минимума в методах Гаусса–Зейделя, наискорейшего спуска, Флетчера–Ривса ?
3. За какое число итераций сходятся исследуемые в лабораторной работе методы в случае квадратичной целевой функции?
4. Как меняется сходимость методов при неквадратичной целевой функции и от чего она зависит?
1.4. Метод экспертных оценок
1.4.1. Получение экспертных оценок
Экспертные оценки получают в результате статистической обработки результатов опроса группы или нескольких групп экспертов.
К опросу экспертов прибегают в случаях, когда математическими методами задачу не решить. Методом экспертных оценок:
определяют совокупности исходных данных – показателей качества, параметров, условий;
выявляют структуру предпочтений – коэффициенты замещения
для разных точек пространства показателей качества;
выбирают параметры критерия предпочтения – структуру обобщенной (интегральной) целевой функции; назначают весовые коэффициенты и т. д.
выбирают наилучшую систему, характеризуемую вектором показателей качества.
Существуют различные методы опроса экспертов применительно к различным решаемым задачам. В большинстве случаев они
сводятся к следующему.
1. Составляют таблицу экспертного опроса (анкету), которую
эксперты должны заполнить, а также пояснения к ней, указывающие как следует заполнять таблицу и какие факторы учитывать.
2. Выбирают состав экспертов из числа специалистов, имеющих
опыт разработки и эксплуатации аналогичных или близких систем. Важно обеспечить репрезентативность выборки и независимость высказываний экспертов.
36
3. Опрашиваемые эксперты заполняют таблицу опроса: проставляют оценки в баллах по 100-5-балльной шкалам, а также ранжируют (присваивают ранги или места от 1-го и далее) всем исследуемым системам или предлагаемым решениям.
4. Производят статистическую обработку полученных данных
для определения средних оценок (например, средний балл или сумма мест).
5. Определяют надежность полученных оценок, а заодно и добросовестность, и квалификацию экспертов.
1.4.2. Обработка экспертных оценок
Пусть N экспертов должны дать r-балльные оценки решениям
Р1, Р2, …, Рm. При обработке результатов r-балльнные оценки часто
переводят в нормализованные оценки wji, занимающие интервал от
нуля до единицы. Их связь с r-балльными оценками линейна:
M
wji = w1ji / å w1ji ,
i=1
здесь w1ji – ненормированная оценка i-го показателя для эксперта
с номером j, причем усреднение проводят по сумме оценок, данных
этим экспертом.
Разброс оценок характеризуют величиной вариации
vi = 2 / i ,
где  и 2 – среднее значение и дисперсия оценки, данных решению
Рi, причем усреднение производят по оценкам N экспертов (по параметру j).
При математической обработке результатов оценок обычно принимают нормальный (гауссов) закон распределения оценок. Параметром, определяющим разброс оценок в ансамбле N, является
среднеквадратическое отклонение .
О надежности оценок судят, например, по разбросу оценок отдельных экспертов и групп экспертов и (или) по степени согласованности
оценок. Степень согласованности оценок, даваемых по совокупности
всех оцениваемых решений, характеризуют так называемым коэффициентом конкордации (согласия) G и соответствующим уровнем
значимости . Здесь коэффициент конкордации 0  G  1, а  есть
вероятность того, что значение согласованности G есть случайное
совпадение. Согласованность оценок считают удовлетворительной,
если G  0,5 при a  0,01, и хорошей, если G  0,7 при a  0,001.
37
Для оценки систем или решений назначают две группы экспертов. Например, в первую группу включают разработчиков систем, а
во вторую – специалистов по эксплуатации. Степень согласованности решений этих двух групп экспертов характеризуют коэффициентом ранговой корреляции , 0    1, причем, чем ближе  к единице, тем лучше согласованы решения. Обычно согласованность
считают удовлетворительной при   0,85 и хорошей при   0,95.
Для вычисления этих величин предварительно находят ранги R1,
R2, …, RM соответствующих решений Р1, Р2, …, РM. В математической статистике ранжированием называют присвоение элементам
порядковых номеров (мест) 1, 2, …, M в зависимости от убывания
какого-либо количественного признака, например качества. Присваиваемый элементу порядковый номер называют рангом этого
элемента. Например, в строю ранжируют по росту, гимнастов – по
степени мастерства. Если несколько сравниваемых элементов имеют одинаковый количественный признак, им присваивают один и
тот же ранг, равный среднему арифметическому. (Например, если
солдаты имеют рост 173, 177, 177, 177 и 180 см, то им присваивают
ранги 5, 3, 3, 3, и 1).
Коэффициент конкордации
é
ù
M
N
ê
ú
G = 12å Si2 / ê N 2 (M 3 - M) - N å Tj ú,
ê
i=1
j=1 úû
ë
где
N
Si = Si - Si , Si = å Rji ,
j=1
здесь R1i – ранг, присвоенный решению Рi j-м экспертом; Si – сумма
рангов, присвоенная решению Рi всеми экспертами:
Si = (S1 + S2 + ... + SM ) / M
– среднее значение суммы рангов, где усреднение производят по
всем M решениям;
Lj
3
Tj = å (tjl
- tjl ),
i=1
где Lj – число групп решений с совпавшими рангами в совокупности
решений j-го эксперта; l – номер группы (с совпавшими рангами);
tjl – число совпавших решений j-го эксперта в группе с номером l.
38
Коэффициент ранговой корреляции
é M
ù
 = 1 - êê6å (Ri1 - Ri2 )2 úú / M (M2 -1),
ëê i=1
ûú
где Ri1 и Ri2 ранги, присвоенные решению Рi, соответственно, 1- и
2-й группами экспертов.
Величину уровня значимости определяют по таблицам распределения случайной величины 2 с (N – 1) степенями свободы, где
é
ù
M
N
ê
ú
2 = 12å Si2 / ê NM (M + 1) - (1 / (M -1)) å Tj ú.
ê
i=1
j=1 úû
ë
Задание на моделирование
Провести экспертную оценку предлагаемых для продажи автомобилей и решить задачу выбора автомобиля. Основные задачи:
1) определить совокупность показателей для оценки имеющихся в продаже автомобилей (6 или более показателей);
2) определить по 100-балльной шкале весовые коэффициенты
для этих показателей;
3) провести оценку автомобилей по 100-балльной шкале в соответствии с выбранными показателями;
4) ввести обобщенный показатель в виде суммы взвешенных
оценок показателей и рассчитать его значение для каждого из автомобилей;
5) ранжировать автомобили в соответствии с полученными значениями обобщенного показателя;
6) оценить степень достоверности проведенных оценок.
Контрольные вопросы
1. В каких случаях применяют экспертные оценки?
2. По каким принципам формируют группы экспертов?
3. Какие системы оценок приняты для проведения экспертиз?
4. Как выявить некомпетентного (недобросовестного) эксперта?
5. Как оценить степень надежности решений экспертов?
6. Какие составляющие входят в формулу для расчета коэффициента конкордации?
7. С какой целью рассчитывают коэффициент ранговой корреляции?
39
8. Какие составляющие входят в формулу для расчета коэффициента ранговой корреляции?
9. Как оценивают уровень значимости полученных оценок?
10. Какими приемами улучшают достоверность оценок?
1.5. Основы регрессионного анализа
1.5.1. Регрессионные модели
В инженерной и научной деятельности важное значение имеют организация (планирование) и обработка результатов эксперимента. Предположим, необходимо выявить закономерность между
некоторой величиной x (предсказывающей переменной – предиктором или независимой переменной – фактором) и откликом y.
Вполне возможно, что предикторов будет много, и мы будем иметь
дело с вектором x, а отклик y есть скалярная величина. Сам эксперимент заключается в следующем.
Составляют план эксперимента, под которым понимают конкретные значения x1, причем задают как число экспериментов для
получения N откликов, так и конкретные значения x1.
Назначают x1 и проводят измерения откликов y1.
На основании эмпирических соображений и инженерного чутья
предлагают математическую модель зависимости y = f(x).
Идентифицируют параметры модели.
Проверяют адекватность модели.
Такой подход к эксперименту называют регрессией, так как анализ результатов делают после того, как получены опытные данные. В задачах такого рода обычно предполагают, что значения независимых переменных x заданы точно, а отклики y подвержены
помехам.
Зависимость математического ожидания M какой-либо случайной величины y от параметра(ов) x называют уравнением регрессии,
а методы экспериментального нахождения уравнения регрессии –
регрессионным анализом.
Точное уравнение регрессии можно записать, зная средние M
для всех допустимых значений x. В практических наблюдениях
такая ситуация невозможна. Более того, даже отдельные средние
значения M не могут быть найдены точно, а допускают лишь приближенные оценки. Поэтому можно искать лишь уравнения приближенной регрессии, оценивая тем или иным способом величину
и вероятность этой приближенности.
40
В общем случае задача ставится так: по данной выборке {(x1,y1),
(x2,y2 ),..., (x,y )} найти уравнение приближенной регрессии (рис.
1.5.1) и оценить допущенную при этом ошибку. В дальнейшем
уравнение приближенной регрессии мы будем записывать в виде
y = f(x), понимая под этим, что на самом деле M = f(x).
Уравнение приближенной регрессии существенно зависит от
выбранного критерия приближения. В методе наименьших квадратов (МНК) за такой критерий принимают дисперсию и оценивают
степень приближения в норме l2. Подбирать это уравнение нужно
по значениям {xi, yi}. Здесь могут возникнуть:
«чистая ошибка», связанная с тем, что измеренные значения
y подвержены случайным воздействиям и не совпадают с соответствующими средними M;
«ошибка интерполяции», ибо число уравнений, по которым
определяют значения коэффициентов регрессии, обычно больше
числа этих коэффициентов, что позволяет сгладить ошибки эксперимента, но не гарантирует прохождения регрессионной кривой
(поверхности) через экспериментальные точки (система переопределенная: число уравнений больше, чем неизвестных).
Из сказанного видно, что бессмысленно гоняться за такими уравнениями, которые на всех парах чисел (x, y) давали бы точное равенство. Последним, к сожалению, грешат очень многие экспериментаторы: построив на координатной плоскости несколько найденных
из опыта точек, они в качестве графика зависимости берут плавную
кривую, непременно проходящую через все построенные точки.
y
x a0 a1w
x
Рис. 1.5.1. Линия регрессии
41
Итак, метод интерполяции, весьма полезный при обработке неслучайных данных, в теории случайных величин неприменим. Необходимо искать принцип приближенности, который бы использовал существенные статистические свойства самого уравнения
регрессии.
Принцип МНК вытекает из общего положения математической
статистики, согласно которому в качестве меры рассеяния всегда
берут дисперсию (среднее из суммы квадратов отклонений). Этот
принцип позволяет сравнить, какое из двух произвольных уравнений дает лучшее приближение регрессии. Обычно сравнивают
функции фиксированного типа, но с произвольными коэффициентами: например, совокупность всех многочленов степени r. Предположим, линия регрессии имеет вид прямой – мы предложили линейную модель: y = 0 + 1x + , что означает – для данного x соответствующее значение y состоит из 0 + 1x плюс некоторая ошибка . Порядок модели соответствует наивысшей степени x, так что
это модель первого порядка. Мы могли бы построить модель второго порядка: y = 0 + 1x + 2 x2 + , которая нелинейна по x, но линейна по параметрам . Мы не умеем находить точно параметры ,
но можем получить их оценки b0 и b1 по результатам наблюдений –
это и называют идентификацией параметров модели. Сумма квадратов отклонений от найденных из опытов значений откликов yi
N
N
i=1
i=1
S = å 2i = å (yi - 0 - 1xi )2 .
Нужно выбрать значения оценок b0 и b1 такими, чтобы их подстановка вместо коэффициентов  в уравнение давала наименьшее
возможное (минимальное) значение S. Для этого продифференцируем выражения по 0 и по 1:
N
dS
= -2å (yi - 0 - 1xi );
d0
i=1
N
dS
= -2å xi (yi - 0 - 1xi ),
d1
i=1
так что для оценок b0 и b1 имеем
N
N
i=1
i=1
å (yi - b0 - b1xi ) = 0; å xi (yi - b0 - b1xi ) = 0.
42
(1.5.1)
Это система так называемых нормальных уравнений, которую
можно преобразовать:
N
N
N
N
N
i=1
i=1
i=1
i=1
i=1
b0 å xi + b1 å xi2 = å xi yi ; b0 N + b1 å xi = å yi .
(1.5.2)
Не будем получать конкретные формулы для расчета коэффициентов b0 и b1, а найдем решение в матричной форме, чтобы обобщить формулу на большее число параметров и модель более высокой степени. Для этого снова запишем систему уравнений:
b0 + b1x1 = y1;
b0 + b1x2 = y2 ;
.
.
b0 + b1xN = yN
или в матричной форме
é1 x1 ù
é y1 ù
ê
ú
ê ú
ê1 x2 ú é b0 ù ê y2 ú
ê
ú * ê ú = ê ú , Xb = y.
ê.
. úú êë b1 úû êê . úú
ê
ê1 x ú
êy ú
N úû
ëê
ëê N ûú
Чтобы найти решение, нужно умножить правую и левую части
уравнения слева на обратную матрицу X–1. Однако эта матрица не
обращается, так как не является квадратной. Поэтому проведем
следующие действия. Умножим слева матрицу X на XТ:
N
ù
é1 x1 ù éê
xi úú
ê
ú ê N
å
é 1 1 ... 1 ù ê1 x2 ú ê
i=1 ú
ú=ê
ê
ú*ê
ú.
ú
ê x1 x2 ... xN ú ê .
N
N
. ú ê
ú
ë
û ê
2
ê1 x ú êê å xi å xi úú
N ûú ê i=1
ëê
i=1
ë
ûú
Затем умножим матрицу XТ на вектор y:
ù
é y1 ù éê N
ê ú ê å yi úú
é 1 1 ... 1 ù ê y2 ú ê i=1 ú
ê
ú*ê ú=
ú.
ê x1 x2 ... xN ú ê . ú êê N
ú
ë
û ê ú ê
ê y ú ê å xi yi úú
êë N úû ê i=1
ë
ûú
43
Теперь нормальные уравнения (1.5.2) можно записать в следующем векторно-матричном виде:
XÒ Xb = XÒ y,
для решения которых нужно умножить правую и левую части слева на обратную матрицу [XТX]–1, которая является квадратной.
Тогда окончательно
b = [XÒ X]-1XÒ y,
что является оценкой параметров модели yˆ = b0 + b1x + b2 x2 + ... +
+ br xr по методу наименьших квадратов, где ŷ – предсказанное по
модели значение y; r – порядок модели.
Итак, МНК основан на двух предположениях:
сумма остатков i = 0;
сумма квадратов остатков 2i = min.
1.5.2. Статистический анализ уравнения регрессии
После вычисления оценок коэффициентов регрессии нужно произвести статистический анализ найденного уравнения регрессии.
Целью статистического анализа является проверка адекватности
уравнения регрессии и значимости его отдельных коэффициентов.
Для проверки гипотезы адекватности необходимо сравнить точность предсказания, достигнутую с помощью модели, с точностью
измерений откликов yi во время эксперимента.
Точность модели можно охарактеризовать величиной дисперсии
адекватности, показывающей степень разброса значений отклика относительно найденного уравнения регрессии
DA =
1
A
N
å (yi - yˆ)2 ,
i=1
где yi – результат i-го измерения; ŷ – значение отклика, предсказанного найденным уравнением регрессии; А = N – (r + 1) – число
степеней свободы, равное числу слагаемых в сумме, уменьшенному
на число линейных связей между ними.
Число связей r + 1 равно числу коэффициентов в уравнении регрессии. Число степенейB свободы равно разности между числом
опытов и параметров (коэффициентов) модели, которые определяют на основании этих опытов.
Точность измерений характеризуется дисперсией 2 случайной помехи . Так как  обычно неизвестна, вместо нее используют
44
«чистую ошибку» – дисперсию воспроизводимости DВ, представляющую выборочную дисперсию отклика, рассчитанную на основании повторных опытов, поставленных в одной, нескольких или
даже всех N точках плана. Если априори известно, что дисперсия
2 измеренных значений y во всех точках плана одинакова, тогда
безразлично в какой из них ставить повторные опыты. Дисперсия
воспроизводимости
1 n
DB =
å (yij -yiñð )2 ,
B j=1
где yij – значение отклика в точке xi при j-м измерении, число степеней свободы B = n – 1.
Если это предположение неоправдано, следует проводить повторные опыты во многих точках плана (но необязательно во всех,
так как мы все время имеем дело с приближенными величинами).
В этом случае формула для дисперсии воспроизводимости усложняется:
DB =
n
N
1 N i
(yij - yiñð ) 2, ãäå B = å ni - N
åå
B i=1 j=1
i=1
и предполагается, что в каждой из точек плана делают ni измерений (число ni  1, но никак не менее единицы, иначе получается,
что мы вовсе и не использовали эту точку плана).Если найденное
уравнение регрессии правильно (адекватно), то DА будет обусловлена только случайным разбросом результатов измерений y относительно их математического ожидания – кривой регрессии y(x).
Она будет несущественно отличаться от DВ. Если же уравнение
регрессии найдено неверно (неадекватно), к случайному разбросу
точек относительно математического ожидания добавится систематическая составляющая, обусловленная разностью между истинным уравнением регрессии и его оценкой (DА будет значимо отличаться от DВ). Это положение иллюстрирует рис. 1.5.2, где показан случай, когда отклик y зависит от одной скалярной величины
x. Для сравнения DА и DВ используют их отношение F = DÀ / DÂ
которое представляет собой так называемый F-критерий проверки гипотезы об адекватности модели. Представленные в числителе и знаменателе дроби гауссовы величины независимы и имеют
2-распределение соответственно с А и В степенями свободы.
Следовательно статистика F подчиняется F-распределению Фишера–Снедекера с А и В степенями свободы. Для проверки гипоте45
x
a0
a
1w
a2 w
2
y
x
a0
w
a1
x
Рис. 1.5.2. Адекватность модели
зы адекватности задают уровень значимости  и находят квантиль
F(1–)(А, В) F-распределения уровня (1–). Таблицы квантилей
F(1–)(А, В) для  = 0,01 и 0,05 приведены в справочниках. Если
F < F(1–)(А, В), можно утверждать, что экспериментальные данные не противоречат проверяемой гипотезе; если же F > F(1–)(А,
В), гипотеза адекватности должна быть отвергнута.
Чтобы выбрать лучшую модель, можно увеличить ее порядок и
добавить квадратичный член b2 x2. Затем снова производят дисперсионный анализ и проверяют адекватность модели по F-критерию:
это так называемый последовательный F-критерий. Процесс прекращают, после того, как гипотеза об адекватности модели может
быть принята.
Вполне возможно, что для описания результатов эксперимента
будет предложена модель нелинейная по параметрам. Например,
это экспоненциальная функция y = 0e–x, или участок синусоиды.
Тогда МНК применять нельзя, а коэффициенты модели следует
искать методом нелинейного программирования. С этой целью составляют целевую функцию и находят ее минимум.
Задание на моделирование
1. Сформировать вектор откликов, соответствующий виду функции. Вид функции указан в приложении.
2. Ввести шумы с равномерным законом распределения в сформированный вектор откликов.
46
3. Найти уравнение регрессии по полученным результатам эксперимента.
4. Проверить адекватность найденного уравнения регрессии.
В случае если модель неадекватна, увеличить порядок модели; повторить пп.1–4.
Контрольные вопросы
1. Что такое регрессионный анализ, уравнение регрессии, факторы, факторное пространство?
2. В чем особенности метода наименьших квадратов?
3. Статистические оценки и их свойства.
4. Что понимается под адекватностью уравнения регрессии?
5. F-критерий и методика его применения.
47
2. МЕТОДИЧЕСКИЕ УКАЗАНИЯ
К ВЫПОЛНЕНИЮ ЗАДАНИЙ НА МОДЕЛИРОВАНИЕ
2.1. Изучение особенностей рельефа
квадратичных целевых функций и функции Розенброка
Работа выполняется в среде Matlab версии 7.1 или выше. Программа для построения целевых функций разбита на несколько
функциональных модулей, которые будут повторно использованы
при выполнении следующих работ.
2.1.1. Датчик квадратичной целевой функции
Создайте подпрограмму – датчик целевой функции. Искомая
функция в общем случае зависит от нескольких переменных.
Вектор-столбец, составленный из значений этих переменных, будет определять точку, в которой необходимо вычислить значение
функции. Если несколько таких столбцов объединить в двумерный массив и передать их датчику целевой функции, за один вызов
можно рассчитать значения искомой функции на множестве точек.
Для этого создайте файл quadratic.m со следующим содержанием.
function f = quadratic(x, theta, a, c)
% f = quadratic(x,theta,a,c) – многомерная квадратичная функция
% входные аргументы:
% x - массив RxN, где R - размерность пространства,
% N - число точек (одна точка занимает один столбец)
% theta – матрица параметров размером RxR
% a – вектор параметров размером Rx1
% с - скаляр
% выходные значения:
% f - вектор-строка 1xN, содержащий значения функции
% Число столбцов равно числу точек
num_points = size(x,2);
% Предварительное выделение памяти для записи значений функции
f = zeros(1,num_points);
% Гарантия того, что вектор параметров a является столбцом
a = a(:);
% Вычисление значений функции в каждой точке
for i=1:num_points
% Извлечение столбца
pp = x(:,i);
f(i) = pp’*theta*pp + pp’*a + c;
end
48
2.1.2. Построение графиков квадратичной функции
Разработайте подпрограмму для отображения целевой функции
в различных видах. Для этого создайте файл lab_one_quadratic.m и
занесите в него следующий код.
function lab_one_quadratic
% Исследование квадратичной функции
%% ===========================================================
% Пределы и шаг изменения координат
x_limits = -1:0.05:1;
y_limits = -1:0.05:1;
% Создание координатной сетки
[cx,cy] = meshgrid(x_limits, y_limits);
% Объединение cx и cy в единый массив,
% содержащий координаты каждой точки в виде столбца
xy = [cx(:) cy(:)]’;
% Задание параметров квадратичной функции
theta = [3 1; 1 2];
a = [0; 0];
c = 1;
% Вычисление значений целевой функции
f = quadratic(xy, theta, a, c);
% Превращение f из вектора-строки в двумерный массив
f = reshape(f, size(cx));
%% ===========================================================
% Построение графика функции в виде поверхности
subplot(2,2,1)
surf(cx, cy, f)
% Построение графика функции в виде сетки
subplot(2,2,2)
mesh(cx, cy, f)
% Построение изолиний
subplot(2,2,3)
contour(cx, cy, f)
% Установка равных масштабов по осям
% (чтобы круги не выглядели, как эллипсы)
axis equal
% Отсечение неиспользуемой части координатного пространства
axis tight
На рис. 2.1.1. показаны примеры графиков, получаемых в результате выполнения приведенного кода.
49
1.5
-0.8
-1
-1
4.5
1.5
3.5
-1 -1
4
-0.5
-0.6
3
2.5
0
y
0 x
-0.5
4.5
0.5
2
1
0.5
55
5. 6
1
1
2.5
1.5
1.5
3 .5
3
2
2
3.5
-0.4
6
4
3
-0.2
3
4.5
5.55
4
2.5
0.2
y
0
3.5
2
5
3
1.5
4
6
4.5
0.4
3.5
0.6
7
2.5
2
3
8
3
2.5
1
0.8
3.5
б)
2.5
а)
2
4
-0.5
2
2.5
3
2.5
3
0
x
0.5
1
Рис. 2.1.1. Поверхность квадратичной функции, построенная функцией
mesh (а), и линии равного уровня, построенные функцией contour (б)
Дополните контурный график квадратичной функции собственными векторами матрицы , как это показано на рис. 2.1.1, б. Собственные векторы постройте в виде единичных отрезков, выходящих из точки минимума, рассчитанной по формуле (1.1.2). Для
этого дополните функцию lab_one_quadratic следующим фрагментом кода.
%% ===========================================================
% Собственные значения и собственные векторы заданной матрицы
% вычисляет функция eig. собственные значения находятся
% на главной диагонали матрицы L,
% собственные векторы расположены в столбцах матрицы V.
[V,L] = eig(theta);
% Извлечение главной диагонали:
L = diag(L);
% Вывод собственных значений и коэффициента обусловленности
% в командное окно:
disp(‘собственные значения:’);
disp(L);
disp(‘коэффициент обусловленности:’);
disp(max(L)/min(L));
% Расчет точки минимума
x0 = -inv(2*theta) * a;
% Построение собственных векторов в виде отрезков,
% выходящих из точки x0:
50
for i = 1:numel(L)
vect = V(:,i);
line([x0(1) x0(1)+vect(1)], [x0(2) x0(2)+vect(2)]);
end
2.1.3. Исследование влияния параметров
квадратичной функции на ее рельеф
Аналогично п. 2.1.2 проведите расчеты при различных значениях параметров, взятых из таблицы. Значение параметра c – произвольное.
Исследуемые комбинации параметров квадратичной функции
№
1
2
3
4
5

é 1 0ù
ê
ú
ê0 1ú
ë
û
é 3 0ù
ê
ú
ê0 1ú
ë
û
é4 1ù
ê
ú
ê1 2ú
ë
û
é2 1 ù
ê
ú
ê1 3ú
ë
û
é15 3ù
ê
ú
ê 3 1ú
ë
û
a
é 0ù
ê ú
ê 0ú
ë û
é2ù
ê ú
ê3ú
ë û
é 0ù
ê ú
ê 0ú
ë û
é3ù
ê ú
ê–1ú
ë û
é–1ù
ê ú
ê2ú
ë û
При построении графиков следите за тем, чтобы в область построения попадали точка минимума и несколько эллипсов. Этого можно добиться подбором значений переменных x_limits и y_
limits, но лучше задавать пределы области построения относительно точки минимума.
2.1.4. Исследование функции Розенброка
Создайте новый файл banana.m, в котором разместите датчик
значений функции Розенброка по формуле (1.1.4).
function f = banana(x)
% f = banana(x) – двумерная тестовая функция Розенброка («банан»)
% входные аргументы:
% x - массив 2xN, где N - число точек
% (одна точка содержит 2 координаты и занимает один столбец)
% выходные значения:
% f - вектор-строка 1xN, содержащий значения функции
if size(x,1) ~= 2
error(‹Входной аргумент должен иметь размерность 2xN›)
end
f = 100 * (x(2,:) - x(1,:).^2).^2 + (1-x(1,:)).^2;
51
Для отображения поверхности и контурного графика этой функции создайте файл lab_one_banana.m, ориентируясь на ранее написанную функцию lab_one_quadratic.m. Расчет точки минимума,
собственных значений и собственных векторов в данном случае не
требуется. Добейтесь того, чтобы точку с координатами (1, 1) можно было визуально определить на графике как минимум функции
(см. рис. 1.1.5, б). Для этого ознакомьтесь со справкой по функции
contour.
Содержание отчета
1. Исходный код программы.
2. Графики квадратичной функции при различных значениях
матрицы и вектора a.
3. Результаты расчета коэффициента обусловленности, собственных чисел и собственных векторов для каждой конфигурации
квадратичной функции. Собственные векторы нанести на соответствующие контурные графики в виде отрезков, выходящих из точки минимума функции.
4. График поверхности функции Розенброка. С помощью контурного графика показать истинное положение минимума. Для
этого необходимо подобрать пределы изменения координат такими, чтобы область минимума хорошо просматривалась, а также
принудительно задать значения функции, для которых будут построены линии равного уровня, в функции contour.
5. Выводы об изученных свойствах целевых функций.
2.2. Исследование методов поиска минимума функций
одной переменной
2.2.1. Варианты заданий
Необходимо разработать в среде Matlab программу, реализующую заданный преподавателем метод одномерного поиска. Получите у преподавателя вариант задания, содержащий метод поиска,
целевую функцию, начальную точку и направление поиска, а также дополнительные параметры.
В таблице использованы следующие обозначения: a0 – начальная точка, d – направление поиска, h – начальный шаг, x – точность по аргументу, f – точность по значению функции, g – точность по значению производной.
52
Варианты заданий
Начальная точка,
Параметры
направление
поиска
поиска
№
Метод
Целевая
функция
1
Бисекция
Функция Розенброка
a0 = [0; 0]T
d = [1; 1]T
h = 0,01
x = 10–4
Дихотомия
Квадратичная,
é5 1ù
é-1ù
ú, a = ê ú, c = 1
=ê
ê1 2ú
ê2ú
ë
û
ë û
a0 = [–2; 1]T
d = [3; 2]T
h = 0,1
x = 10–5
Функция Била
a0 = [2; 0,7]T
d = [1; –0,2]T
h = 0,1
x = 10–5
Функция Химельблау
a0 = [0; 3]T
d = [1; –1]T
h = 0,3
x = 10–8
Функция Химельблау
a0 = [–2; –3]T
d = [–1; –0,5]T
h = 0,3
x = 10–8
Функция Била
a0 = [1; 0]T
d = [3; 0,3]T
h = 0,2
x = 10–8
Функция Розенброка
a0 = [0; 1]T
d = [–1; –1]T
h = 0,5
x = 10–8
Квадратичная,
é 6 -2ù
é2ù
ú, a = ê ú, c = 3
=ê
ê-2 3 ú
ê5ú
ë
û
ë û
a0 = [0; 0]T
d = [–0,5;1]T
h = 0,5
x = 10–3
2
3 Золотое сечение
Квадратичная
интерполяция
без производной
Квадратичная
5 интерполяция
с производной
Кубическая
6 интерполяция
без производных
Кубическая
7 интерполяция
с производными
4
8
Фибоначчи
Функция Била от двух переменных:
2
(
(
f (x, y) = (1,5 - x (1 - y)) + 2,25 - x 1 - y2
2
2
)) + (2,625 - x(1- y3 )) .
Функция Химмельблау от двух переменных:
(
2
) (
2
)
f (x, y) = x2 + y -11 + x + y2 - 7 .
2.2.2. Структура программы
Процедуру одномерного поиска необходимо проектировать таким образом, чтобы было возможно использовать ее для нахождения минимума различных целевых функций с различными на53
lab_two.m
датчик ц. ф.
датчик градиента
параметры поиска
Точка входа в программу
координата
минимума
minimize_linear.m
function x_min =
minimize_linear(…
funobj, gradobj, options)
funobj.m
function
f=funobj(x)
gradobj.m
function
g=gradobj(x)
Рис. 2.2.1. Структура программы одномерного поиска
чальными условиями. Для этого подходит иерархическая структура, показанная на рис. 2.2.1.
Метод поиска минимума выделен в m-функцию minimize_linear,
которая обращается к датчикам значений целевой функции и ее
градиента (при необходимости). Точкой входа в программу служит
функция верхнего уровня lab_two, которая выполняет следующие
действия:
определяет целевую функцию и ее градиент;
задает условия и параметры поиска;
вызывает функцию поиска минимума;
обрабатывает результаты поиска.
2.2.3. Датчики целевой функции и ее градиента
Используйте датчики значений целевой функции из разд. 2.1,
для других вариантов создайте аналогичные. Для методов, использующих производные, создайте в отдельном файле датчик
градиента, вычисляющий первые частные производные заданной
функции. Например, для квадратичной функции создайте файл
quadratic_grad.m со следующим содержанием.
function g = quadratic_grad(x, theta, a)
% g = quadratic(x,theta,a,c) – градиент квадратичной функции
% входные аргументы:
% x - массив RxN, где R - размерность пространства,
% N - число точек (одна точка занимает один столбец)
% theta – матрица параметров размером RxR
54
%
%
%
%
%
a – вектор параметров размером Rx1
с - скаляр
выходные значения:
g - матрица RxN, содержащая первые частные производные
квадратичной функции в строках
% Размерность пространства параметров
num_dims = size(x,1);
% Число столбцов равно числу точек
num_points = size(x,2);
% Предварительное выделение памяти
g = zeros(num_dims,num_points);
% Гарантия того, что вектор параметров a является столбцом
a = a(:);
% Вычисление значений функции в каждой точке
for i = 1:num_points
% Извлечение столбца
pp = x(:,i);
g(:,i) = 2*theta*pp + a;
end
Датчик значений целевой функции можно объединить в одном
файле с датчиком градиента, используя m-функцию с переменным
числом выходных аргументов.
2.2.4. Работа с функциональным типом данных
языка Matlab
Функциональный тип данных в данной работе помогает унифицировать работу с различными целевыми функциями и различными начальными условиями. Пусть датчиком является m-функция
с именем “fun”, принимающая аргумент x в виде вектора-столбца (см. рис. 2.2.1). Для превращения многомерной функции в
одномерную создается вспомогательная переменная funobj типа
function_handle, осуществляющая проекцию скалярного аргумента t на направление поиска в многомерном пространстве:
funobj = @(t) fun(a0(:) + t*d(:));
где a0 – начальная точка, d – направление поиска, а оператор (:)
превращает матрицу произвольной размерности в вектор-столбец.
Данный прием, известный также как «связывание аргументов»,
описан в справке Matlab в разделе «anonymous functions». К созданному объекту funobj можно обращаться так же, как к любой
другой функции среды Matlab. Например, график изменения це55
левой функции вдоль заданного направления строится следующим
фрагментом кода:
t = 0:h/100:h;
f = zeros(size(t));
for i = 1:numel(t)
f(i) = funobj(t(i));
end
plot(t,f)
Если требуется получить функцию одной переменной, а исходная функция имеет более одного аргумента, к остальным аргументам привязываются фиксированные значения. Например, при работе с квадратичной функцией:
funobj = @(t) quadratic(a0 + t*d, [4 1;1 2], [2;-1], 0);
Производная по направлению поиска формируется аналогичным образом:
gradobj = @(t) d’ * quadratic_grad(a0+t*d, [4 1;1 2], [2;-1]);
2.2.5. Реализация заданного метода
Ознакомьтесь с примерным содержанием функции верхнего
уровня lab_two.
function lab_two
% lab_two – одномерный поиск минимума
% начальная точка
a0 = [-2; 1];
% направление поиска
d = [3; 2];
% нормировка
d = d./ norm(d);
% датчики
funobj = @(t) quadratic(a0 + t*d, [5 1; 1 2], [-1;2], 0);
gradobj = @(t) d› * quadratic_grad(a0 + t*d, [5 1; 1 2], [-1;2]);
% параметры поиска
options.max_iter = 100; % макс. число итераций
options.h = 0.1; % начальный шаг
options.eps = 1e-8; % точность
% запуск одномерного поиска
t_min = minimize_linear(funobj, gradobj, options);
% вычисление точки минимума и значения функции в ней
a_min = a0 + t_min*d;
f_min = funobj(t_min);
56
% вывод результатов в командное окно
fprintf(‘\nТочка минимума: (%f, %f)\nЗначение функции: %f\n’,...
a_min(1), a_min(2), f_min);
Входными данными для минимизации выбранным методом являются датчик значений целевой функции, датчик градиента целевой функции, начальная точка и направление поиска в многомерном пространстве, а также параметры точности, необходимые
для критериев завершения поиска. Алгоритм работы m-функции
minimize_linear должен повторять рис. 1.2.1 с учетом специфики
конкретного метода. Для облегчения реализации рассмотрим процедуру поиска начального интервала неопределенности.
function tmin = minimize_linear(fun, grad, options)
% начальный шаг
step = options.h;
% границы интервала
a = 0;
b = step;
% значения функции
fa = fun(a);
fb = fun(b);
% цикл не более options.max_iter итераций
for it = 1:options.max_iter
if fb > = fa
% функция возрастает
% возврат на половину шага
c = b - step/2;
fc = fun(c);
if fc < = fa
break;
else
b = c;
a = a - step/2;
end
break
else
% функция убывает
step = step*2;
a = b;
b = a + step;
end
end
fa = fb;
fb = fun(b);
57
% далее производится минимизация на интервале a < = t < = b,
% заканчивающаяся присвоением значения переменной tmin
Приведенный пример, в отличие от алгоритма из п. 1.2.3, использует начальную точку t = 0, соответствующую точке a0 в многомерном
пространстве, и всегда выполняет первый шаг в сторону возрастания
переменной t (так как направление поиска уже задано вектором d).
2.2.6. Графическое представление результатов
Для облегчения отладки необходимо ввести в программу графическое отображение результатов поиска. Для этого график одномерной целевой функции на интервале неопределенности строится
с помощью приведенного выше фрагмента кода, а затем, по мере
выполнения итераций, на этот график наносятся в виде маркеров
оценки точки минимума. Очевидно, что при отсутствии ошибок в
программе оценка точки минимума будет постепенно приближаться к истинному минимуму функции.
2.2.7. Тестирование разработанной программы
Проведите тестирование разработанной программы, варьируя
начальную точку и направление поиска. Зафиксируйте результаты 4–5 запусков, отмечая успешность или не успешность поиска, а
также количество выполненных итераций.
Содержание отчета
1. Функциональная схема заданного алгоритма поиска минимума.
2. Исходный код программы.
3. Результаты работы программы в текстовом и графическом виде.
4. Зависимость числа итераций от задаваемой точности поиска.
5. Выводы о сходимости реализованного метода в различных условиях.
2.3. Исследование методов поиска минимума функций
многих переменных
2.3.1. Варианты заданий
Реализовать в среде Matlab методы многомерного поиска: cимплекс-метод; метод наискорейшего спуска; метод Гаусса-Зейделя;
метод сопряженных градиентов; квазиньютоновский метод.
58
За исключением симплекс-метода, во всех вариантах задания
потребуется подпрограмма одномерного поиска minimize_linear,
которая является результатом выполнения лабораторной работы
из разд. 2.2.
Получите у преподавателя вариант задания, содержащий целевые функции, на которых метод должен быть протестирован.
№
1
2
3
4
5
Квадратичная
да
да
да
нет
да
Функции для тестирования
Розенброка
Экспоненциальная
нет
нет
да
да
да
да
да
да
да
да
Била
нет
нет
нет
да
нет
Экспоненциальная функция имеет следующую запись
f (x, y) = ye-x
2
-y2
.
Данная функция является примером неквадратичной функции,
имеет один минимум в точке [0; –0,707] и один максимум в точке
[0; 0,707]. Формулы для остальных функций приведены в разд. 2.2.
2.3.2. Реализация заданного метода
Симплекс-метод реализуется в соответствии с алгоритмом, приведенным на рис. 1.3.2. В качестве примера градиентного метода
рассмотрим реализацию метода покоординатного спуска
function gauss_seidel_demo
% минимизация методом покоординатного спуска
% датчики многомерной функции и градиента
fun_multi = @(x) quadratic(x, [5 1; 1 2], [-1;2], 0);
grad_multi = @(x) quadratic_grad(x, [5 1; 1 2], [-1;2]);
% параметры многомерного поиска
a0 = [1;-1.2]; % начальная точка
tol = 1e-8; % точность
max_iter = 50; % максимальное число итераций
% параметры одномерного поиска
options.max_iter = 100; % макс. число итераций
options.h = 0.01; % начальный шаг
options.eps = 1e-8; % точность
% выделение памяти для построения траектории
track = zeros(max_iter,3);
59
% основной цикл
for iter = 1:max_iter
f0 = fun_multi(a0);
g0 = grad_multi(a0);
track(iter,:) = [a0(1) a0(2) f0];
if iter > 1
fprintf(‘Итерация %d, точка (%.3f, %.3f), f = %.3f\n’,...
iter-1, a0(1), a0(2), f0);
end
% проверка критериев завершения поиска
if iter > 1 && track(iter-1,3) - track(iter,3) < tol
break;
end
% формирование направления поиска
num_dims = numel(g0);
% направлением поиска является частная производная
% по одной из координат, номер которой циклично меняется
% в зависимости от номера итерации
k = 1 + mod(iter-1, num_dims);
d = zeros(size(g0));
d(k) = -g0(k);
d = d./ norm(d);
funobj = @(t) fun_multi(a0 + t*d);
gradobj = @(t) d’ * grad_multi(a0 + t*d);
t_min = minimize_linear(funobj, gradobj, options);
a0 = a0 + t_min*d;
end
% построение траектории
px = track(1:iter,1);
py = track(1:iter,2);
line(px, py, ‘Color’, ‘red’, ‘Marker’,’o’);
Дополните программу, построив линии равного уровня исследуемой целевой функции вместе с траекторией поиска. Это позволит более наглядно представить траекторию поиска относительно
истинного положения точки минимума. Введите дополнительные
критерии завершения поиска по малому изменению координат минимума и по уменьшению нормы градиента.
2.3.3. Тестирование программы
Протестируйте программу, запуская ее на заданных целевых
функциях. Пробуя различные начальные точки, не слишком далеко удаленные от минимума, выявите характерные черты поведе60
ния метода на квадратичных и неквадратичных функциях. Зафиксируйте свои наблюдения.
Содержание отчета
1. Исходный код программы.
2. Траектории поиска минимума реализованным методом для
5–6 характерных случаев.
3. Число итераций и координаты минимума для каждого эксперимента.
4. Выводы о сходимости данного метода. Анализ случаев отсутствия сходимости, если таковые имели место.
2.4. Организация экспертного совещания
и обработка результатов экспертных оценок
Работа выполняется в пакете Microsoft Exell.
1. Определим список показателей для оценки имеющихся в продаже автомобилей.
Список показателей:
k1 – скоростные качества (приемистость, максимальная скорость); k2 – потребление горючего; k3 – мощность двигателя; k4 –
состояние (год выпуска, необходимость ремонта и пр.); k5 – внешний вид и уровень комфорта; k6 – стоимость.
Предположим, в результате работы экспертов получены следующие весовые коэффициенты для показателей (табл. 1):
Таблица 1
Показатель
k1
k2
k3
k4
k5
k6
Вес wj
0,111
0,166
0,083
0,278
0,139
0,22
2. Организуем две группы экспертов и проведем оценку предлагаемых к продаже автомобилей каждым экспертом.
Результатом работы одного из экспертов являются данные, приведенные в табл. 2 и 3.
3. Приведем оценки к стандартному виду (от 0 до 1) путем нормирования по следующей формуле:
M
kji = Kji / å Kji ,
i=1
где Kji – ненормированная оценка i-го показателя для эксперта с
номером j, причем усреднение проводят по сумме оценок, данных
этим экспертом.
61
Таблица 2
Автомобили
Волга
Москвич 2141
ВАЗ 2107
ВАЗ 2110
Опель-Вектра
Мерседес 190
Сумма
Первичные оценки экспертов
K2
K3
K4
K5
K1
35
20
50
40
70
70
285
50
90
85
30
60
60
375
70
50
70
50
70
50
360
60
80
60
100
90
45
435
K6
30
35
80
37
80
50
312
70
90
80
60
75
30
405
Таблица 3
Автомобили
Волга
Москвич 2141
ВАЗ 2107
ВАЗ 2110
Опель-Вектра
Мерседес 190
Нормированные оценки
k1
0,123
0,070
0,175
0,140
0,246
0,246
k2
0,133
0,240
0,227
0,080
0,160
0,160
k3
0,194
0,139
0,194
0,139
0,194
0,139
k4
0,138
0,184
0,138
0,230
0,207
0,103
k5
0,096
0,112
0,256
0,119
0,256
0,160
k6
0,173
0,222
0,198
0,148
0,185
0,074
Таблица 4
Автомобили
Волга
Москвич 2141
ВАЗ 2107
ВАЗ 2110
Опель-Вектра
Мерседес 190
Взвешенные оценки
Cумма Ранги
w1*k1 w2*k2 w3*k3 w4*k4 w5*k5 w6*k6
0,01
0,01
0,02
0,02
0,03
0,03
0,02
0,04
0,04
0,01
0,03
0,03
0,02
0,01
0,02
0,01
0,02
0,01
0,04
0,05
0,04
0,06
0,06
0,03
0,01
0,02
0,04
0,02
0,04
0,02
0,04
0,05
0,04
0,03
0,04
0,02
0,14
0,18
0,19
0,15
0,20
0,13
5
3
2
4
1
6
Примечание: Все импортные машины имеют возраст 5–10 лет.
Такие таблицы необходимо получить для каждого эксперта,
участвующего в решении задачи о приобретении оптимальной автомашины.
4. Cоставим сводные таблицы и определим коэффициент конкордации.
Для примера приведем гипотетические результаты экспертизы
(табл. 5). О надежности этих оценок можно судить по величине вариаций. Видно, что она колеблется от 0,177 до 0,386. Обычно считают удовлетворительным, если вариация меньше 0,3 и хорошим,
62
если меньше 0,2. В данном примере результаты оценок можно считать удовлетворительными только с большой натяжкой.
Таблица 5
№
№ эксгруппы перта Волга
1
1
2
3
4
5
6
7
2
8
9
10
Средние
значения
Дисперсии
Вариации
Оценки экспертов
Сумма
для нормированных весов показателей
абсолютных
МоскВАЗ
ВАЗ
Опель- Мерсеотклонений
вич
2107
2110
Вектра дес 190
0,14
0,256
0,156
0,111
0,263
0,18
0,128
0,156
0,083
0,053
0,19
0,128
0,156
0,139
0,184
0,15
0,103
0,222
0,278
0,211
0,2
0,256
0,222
0,222
0,158
0,13
0,128
0,089
0,166
0,132
0,182
0,245
0,173
0,304
0,230
0,25
0,2
0,103
0,196
0,286
0,156
0,171
0,128
0,087
0,086
0,188
0,086
0,231
0,217
0,143
0,125
0,286
0,256
0,13
0,114
0,219
0,143
0,205
0,196
0,2
0,063
0,114
0,077
0,174
0,171
0,250
0,300
0,282
0,200
0,272
0,196 0,123
0,166
0,188
0,202
0,124
0,004 0,002
0,339 0,353
0,002
0,265
0,005
0,379
0,001
0,161
0,002
0,315
Для вычисления коэффициента конкордации составим табл. 6
рангов для оценок, полученных из таблиц, предоставленных экспертами (см. табл. 4 результатов эксперта под номером 1).
Таблица 6
№ группы № эксперта
1
2
3
1
4
5
6
7
8
2
9
10
Сумма мест
Отклонения
от среднего DS
Ранги, присвоенные экспертами
5
1,5
4
5
1
1
2
5
2,5
1
25
3
4
4
6
6
4
3
4
6
6
28
2
4
4
4
3
3
6
2
1
4
46
4
6
1,5
1
2
5
1
1
5
5
33
1
1,5
1,5
2
4
2
4
3
2,5
2
31,5
6
4
6
3
5
6
5
6
4
3
23,5
–10
–7
11
–2
–3,5
–11,5
Среднее
значение
DS
35
487,5
63
Данные, приведенные в табл. 6, позволяют рассчитать коэффициент конкордации G и оценить его значимость по таблице распределения случайной величины 2 с (N–1) степенями свободы. Первый эксперт все ранги сделал различными. Поэтому L1 = 0, T1 = 0.
У второго эксперта имеются две группы (L2 = 2) совпадающих оценок: первая группа содержит два совпадающих ранга R1 и R5; вторая – три совпадающих ранга R2, R3, R6. Поэтому t21 = 2, t22 = 3 и
T2 = (2^3–2)+(3^3–3) = 30; T3 = 30; T4–T8 = 0; T9 = 6; T10 = 0; Сумма
T = 66. Подставляя полученные значения в формулу для коэффициента конкордации, получим
G = 0,287611; 2 = 14,38.
По таблицам распределения 2 находим, что величина 14,38
соответствует величине вероятности  < 0,01. Значит с вероятностью 1 –  = 0,99 можно считать, что согласованность мнений экспертов, характеризуемая величиной коэффициента конкордации
G = 0,344, неслучайна. Однако величина коэффициента конкордации недостаточно велика для того, чтобы считать совпадение мнений экспертов удовлетворительным, а оценки – надежными.
Таблица 7
№
группы
№
эксперта
1
2
1
3
4
5
Сумма мест
Ранги по группе
6
7
2
8
9
10
Сумма мест
Ранги по группе
Ранги, присвоенные экспертами
5
1,5
4
5
1
16,5
3
1
2
5
2,5
1
11,5
1
3
4
4
6
6
23
5
4
3
4
6
6
23
1
2
4
4
4
3
17
4
3
6
2
1
4
16
5
4
6
1,5
1
2
14,5
2
5
1
1
5
5
17
3
1
1,5
1,5
2
4
10
1
2
4
3
2,5
2
13,5
4
6
4
6
3
5
24
6
6
5
6
4
3
24
2
Причиной ненадежности оценок может быть некомпетентность
отдельных экспертов. Первым признаком ненадежного эксперта
64
является сумма абсолютных значений отклонений оценок данного эксперта от усредненных значений, данная в правом крайнем
столбце табл. 5. Из табл. 5 видно, что мнения экспертов 4, 7 и 8-го
существенно отличаются от усредненных оценок. Затем нужно исключать экспертов по одному и вычислять для оставшихся экспертов коэффициент конкордации. Если после исключения какого-то
эксперта окажется, что коэффициент конкордации существенно
возрос, то его нужно вывести из состава группы или выяснить причины его поведения. Можно рассмотреть вариант с исключением
еще одного эксперта, если коэффициент конкордации увеличится
недостаточно.
Найдем теперь, насколько совпадают мнения экспертов двух
групп. По данным табл. 7 рассчитаем коэффициент ранговой корреляции  = 0,71. Значение r показывает, что мнения групп экспертов в достаточной мере совпадают.
Реализуем выбор автомобиля на основе полученных данных.
Подготовительную работу можно считать оконченной и приступить
непосредственно к выбору автомобиля для покупки. С этой целью
проведем оценку автомобилей по разработанным показателям. Затем сформируем обобщенный показатель, например, в виде суммы
взвешенных значений отдельных показателей
N
kîáù = å
M
å wi kji ,
j=1 i=1
причем на первое место поставим автомобиль, для которого сумма
M взвешенных показателей по совокупности N экспертов приобретает максимум.
Содержание отчета
1. Таблицы экспертных оценок, выполненные в программном
пакете Microsoft Excel.
2. Обоснованное решение о покупке конкретного автомобиля.
2.5. Обработка результатов эксперимента
с использованием регрессионного анализа
Работа выполняется в среде Matlab.
Рассмотрим алгоритм моделирования на примере двумерной
гауссоиды (рис. 2.5.1).
65
Рис. 2.5.1. График исходной функции
Введем датчик исследуемой функции:
function f = gaussoid(x, x0, sigma)
% Многомерная гауссоида
% Координаты каждой точки должны лежать в матрице x в виде столбцов
% число столбцов - это число точек
num_points = size(x,2);
% Размерность
n = size(x,1);
% Массив значений функции
f = zeros(1,num_points);
% Нормировочный множитель
c = 1 / ((2*pi)^(n/2) * sqrt(det(sigma)));
sigma = inv(sigma);
% Вычисление значений функции в каждой точке
for i=1:num_points
pp = x(:,i) - x0(:);
f(i) = c * exp(-0.5 * pp’ * sigma * pp);
end
В двумерном случае приведенный код соответствует вычислениям
по формуле
f (x, y) =
66
1
2 
1
- xT -1x
e 2
, x = [ x - x0 y - y0 ],
где x0, y0 – координаты центра гауссоиды,  – ковариационная
матрица. Для построения поверхности исследуемой функции
необходимо задать параметры гауссоиды и координаты точек.
% параметры гауссоиды
x0 = [0; -1];
sigma = diag([5 8]);
% центральная точка
% ковариационная матрица
% набор точек на плоскости
R = 5; % диапазон изменения координат
small_step = 0.5; % шаг изменения координат
[x,y] = meshgrid(-R:small_step:R, -R:small_step:R);
% вычисление функции в заданном наборе точек
p = [x(:) y(:)]’;
f = gaussoid(p, x0, sigma);
f = reshape(f, size(x));
% построение поверхности
surf(x, y, f, ‘EdgeColor’, ‘none’);
Введем нормированные переменные,
пределах 1
которые изменяются в
nx = x(:) / R;
ny = y(:) / R;
Нормированные переменные нужны, чтобы план эксперимента
был представлен в стандартном виде. В исходную функцию нужно
подставлять истинные (ненормированные) переменные.
При использовании линейной модели исходная функция аппроксимируется плоскостью:
fˆ(x, y) = ax + by + c.
Матрица плана эксперимента для линейной модели имеет вид:
é x1
ê
ê x2
A = êê
ê ...
êx
ëê N
y1
y2
...
yN
1ù
ú
1ú
ú,
1úú
1úûú
что записывается в Matlab как
Plan = [nx(:) ny(:) ones(numel(nx),1)];
Смоделируем погрешность измерений путем добавления шума к
значениям исходной функции:
67
noise_std = 0.001; % СКО шума
Response = Response(:) + noise_std * randn(numel(Response),1);
Вектор параметров модели вычисляется по формуле МНК:
Cf = inv(Plan’ * Plan) * Plan’ * Response;
что дает результат
>> Cf =
>>
0.0000
>> -0.0004
>>
0.0064
(значения могут меняться под действием шума). В среде Matlab для
решения системы линейных уравнений предпочтительно использовать оператор «\», который, в зависимости от свойств матрицы,
применяет один из способов факторизации (например, разложение
по сингулярным числам):
Cf = Plan’ * Plan \ Plan’ * Response;
В результате использования линейной модели получен
совершенно неадекватный результат: плоскость, имеющая мало
общего с исходной функцией (рис. 2.5.2). Для проверки адекват-
0,16
0,14
f(x, y)
0,12
0,01
0,008
0,006
0,004
0,002
0
5
0
y
–5
–5
3 4
2
0 1
3 4
2
1
x
–5
Рис. 2.5.2. График исходной функции и ее линейная модель
68
ности модели нужно рассчитать дисперсии адекватности DА и
воспроизводимости DВ. Число степеней свободы  A = N - (r + 1)
равно шести, так имеется три параметра модели r+1=3 и проведено
9 опытов N=9.
Рассчитаем значения функции в точках плана эксперимента по
полученной модели:
Model = Plan * Cf;
Дисперсия адекватности рассчитывается как средний квадрат
отклонения экспериментальных данных от откликов, полученных
по модели:
DA = mean((Model – Response).^2);
Для расчета дисперсии чистой ошибки в реальных условиях
проведения регрессионного эксперимента нужны повторные опыты в точках плана. Поскольку при моделировании ошибки введены искусственно, дисперсия воспроизводимости равна дисперсии
шума, подмешанного к экспериментальным данным:
DB = noise_std^2;
Значение F-критерия для проведенного моделирования
вычисляется как отношение дисперсии адекватности к дисперсии
воспроизводимости:
Fval = DA/DB
>> Fval =
>> 40.6315
Число степеней свободы В равно числу точек в плане
эксперимента минус 1. Число степеней свободыА равно числу точек
в плане эксперимента, уменьшенному на число параметров модели
(для плоскости – 3). Квантиль F-распределения, соответствующая
заданному уровню доверия, рассчитывается с помощью датчика
обратной функции распределения icdf, имеющегося в среде
Matlab.
fi_a = size(Plan,1) - size(Plan,2);
fi_b = size(Plan,1) - 1;
qF = icdf(‘f’,0.99, fi_a, fi_b) % уровень доверия 0,99
>> qF =
>>
2.6953
Поскольку Fval > qF, считаем, что гипотеза о равенстве дисперсий неверна и нужно повысить степень модели – например,
применить квадратичную функцию.
69
Таблица 8
Вариант
1
2
Вид функции
Пространственная Гауссоида:
æ (r 2 + r 2 ) ö÷
ç
1
x
y ÷
exp ççR (r) =
÷÷÷
2
ç
 2
2
÷ø
çè
rx = x–x1, ry = y–y1, (x1,y1) – центр функции
Кривая Релея:
æ (r 2 + r 2 ) ö÷
rx2 + ry2
ç
x
y ÷
÷÷
exp ççR (r) =
2
çç
÷
2
2

è
ø÷
3
Пространственная синусоида:
E(x, y) = Em sin (2x + 2y + 0 )
4
Пространственная косинусоида:
E(x, y) = Em cos (2x + 2y + 0 )
Содержание отчета
1. Задание.
2. Исходные данные для эксперимента (вектора откликов и матрица плана эксперимента).
3. Полученные уравнения регрессии и результаты их статистического исследования.
4. Выводы.
Библиографический список
1. Тимофеев Б. С. Автоматическая настройка ТВ-систем: учеб.
пособие. Л.: ЛИАП. 1984.
2. Банди Б. Методы оптимизации. Вводный курс. М.: Радио и
связь, 1988.
3. Аоки М. Введение в методы оптимизации. М.: Наука. 1977.
4. Тимохов А. В. Теория оптимизации в задачах и упражнениях:
учеб. пособие. СПб.: Лань, 2012.
5. Дрейпер Н., Смит Г. Прикладной регрессионный анализ. М.:
Финансы и статистика, 1973.
6. Браунли К. А. Статистическая теория и методология в науке и
технике. М.: Наука, 1977.
7. Румшинский Л. З. Математическая обработка результатов
эксперимента. М.: Наука, 1981.
70
8. Гроп Д. Методы идентификации систем. М.: Мир, 1979.
9. Формалев В. Ф., Ремизников Д. Л. Численные методы. М.:
Физматлит, 2004.
10. Измаилов А. Ф., Солодов М. В. Численные методы оптимизации: учеб. пособие. М.: Физматлит, 2005.
11. Тимофеев Б. С. Автоматическая настройка телевизионных
систем с помощью микроЭВМ. М.: Радио и связь, 1988.
СОДЕРЖАНИЕ
1. Элементы теории оптимизации ........................................
1.1. Исследование целевых функций ...............................
1.2. Методы поиска минимума функций одной
переменной .................................................................
1.3. Методы поиска минимума функций многих
переменных .................................................................
1.4. Метод экспертных оценок ........................................
1.5. Основы регрессионного анализа ................................
2. Методические указания к выполнению заданий
на моделирование .............................................................
2.1. Изучение особенностей рельефа квадратичных
целевых функций и функции Розенброка .........................
2.2. Исследование методов поиска минимума функций
одной переменной ........................................................
2.3. Исследование методов поиска минимума функций
многих переменных ......................................................
2.4. Организация экспертного совещания и обработка
результатов экспертных оценок ......................................
2.5. Обработка результатов эксперимента с использованием
регрессионного анализа .................................................
Библиографический список ................................................
3
3
12
26
36
40
48
48
52
58
61
65
70
71
Учебное издание
Афанасенко Арсений Сергеевич
Обухова Наталья Александровна
Тимофеев Борис Семенович
ОСНОВЫ ТЕОРИИ ОПТИМИЗАЦИИ
Учебное пособие
Редактор А. В. Подчепаева
Верстальщик С. Б. Мацапура
Сдано в набор 9.09.13. Подписано к печати 19.12.13.
Формат 6084 1/16. Бумага офсетная. Усл. печ. л. 4,2.
Уч.-изд. л. 4,5. Тираж 100 экз. Заказ № 648.
Редакционно-издательский центр ГУАП
190000, Санкт-Петербург, Б. Морская ул., 67
Документ
Категория
Без категории
Просмотров
0
Размер файла
858 Кб
Теги
afanasenkoobuhova
1/--страниц
Пожаловаться на содержимое документа