close

Вход

Забыли?

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

?

FreeBASIC09

код для вставкиСкачать
Пособие для начинающих программировать на языке высокого уровня FreeBASIC. Работа состоит из нескольких глав. В главе 9 "Методы оптимизации" – продолжение пособия, которое будет интересно для учащихся школ, студентов институтов, а также преподавател
 1
Евгений Рыжов, инженер
Программирование на языке FreeBASIC
пособие для начинающих
Пособие для начинающих программировать на языке высокого уровня FreeBASIC.
Работа состоит из нескольких глав. В главе 9 "Методы оптимизации" –
продолжение пособия, котор
ое будет интересно для учащихся школ,
студентов институтов, а также преподавателей.
Фрагмент 9. Методы оптимизации.
Февраль. Достать чернил и плакать!
Писать о феврале навзрыд,
Пока грохочущая слякоть
Весною черною горит.
Достать пролетку. За шест
ь гривен,
Чрез благовест, чрез клик колес,
Перенестись туда, где ливень
Еще шумней чернил и слез.
2
Где, как обугленные груши,
С деревьев тысячи грачей
Сорвутся в лужи и обрушат
Сухую грусть на дно очей.
Под ней проталины чернеют,
И ветер криками изрыт,
И чем случайней, тем вернее
Слагаются стихи навзрыд.
Борис Пастернак, 1912
Еще даже не весна, а картинка просится осеняя. Называется картина "17 октября 1905 года", автор ее Илья Ефимович Репин
(1844
—
1930) —
русский художник, живописец, мастер портретов, исторических и бытовых сцен, мемуарист, автор ряда очерков, составивших книгу воспоминаний "Далёкое близкое".
Картина "17 октября 1905 года" (1907, 1911)
Вот что об этой картине писал русский религиозный философ, литературный критик и публицист Васили
й Васильевич Розанов (1856 —
1919):
Сколько понимания, сколько верности! Конечно, все жившие в 1905
-
1906 годах в Петербурге скажут о картине: "Это -
так! Это -
верно!" Несут на плечах маньяка, с сумасшедшим выражением лица и потерявшего шапку. "До шапки ли
тут, когда конституция". Лицо его не ясно в мысли, как именно у сумасшедшего, и видны только "глаза в одну точку" и расклокоченная борода. Это "назарей" революции, к шевелюре которого вообще никогда не притрагивались ножницы, бритва, гребенка и щетка. Умс
твенная роль его небольшая: самого его несут на плечах, и он в свою очередь высоко держит над толпою "венок победы". Таким образом маньяк, как ему и следовало, вышел в простую деревянную подставку для плаката. Впереди всей процессии два гимназиста, и один не старше IV или V класса, но и другой, старший, 3
ближайший к зрителю, тоже не VIII класса, а класса VI или VII. Кто видал массы гимназистов, не ошибется, взглянув на лицо, к которому классу относится "питомец школы". Два эти гимназиста и стоящий позади шес
тиклассника студент в фуражке, положивший ему руки на плечи, -
"инструктор" пения и идей -
какая это опера!! Боже, до чего все это -
так!!
Прошли годы, но все те
же лица, те
же проблемы... Вот и вспомнился русский философ, писатель и публицист Иван Алекса
ндрович Ильин (1883 -
1954). Приведу только две цитаты
из его сочинений:
И.А.Ильин "Зависть как источник бедствий" (1 июля 1952)
...всплывают новые силы -
новые диктаторы, новые классы, новые нации. Эти диктаторы принадлежат к полуинтеллигенции, думают уп
рощающе, не ведают ни правосознания, ни чувства ответственности, но одержимы волею к необузданной власти. Эти новые классы не имеют ни малейшего представления о религии, о душе и о культуре...
И.А.Ильин "Что сулит миру расчленение России" (15 июля 1950)
.
..когда после падения большевиков мировая пропаганда бросит во всероссийский хаос лозунг "Народы бывшей России, расчленяйтесь!", то откроются две возможности: или внутри России встанет русская национальная диктатура, которая возьмет в свои руки крепкие "бр
азды правления", погасит этот пробельный лозунг и поведет Россию к единству, пресекая все и всякие сепаратистские движения в стране; или же такая диктатура не сложится, и в стране начнется непредставимый хаос передвижений, возвращений, отмщений, погромов, развала транспорта, безработицы, голода, холода и безвластия. Тогда Россия будет охвачена анархией и выдаст себя с головой своим национальным, военным, политическим и вероисповедным врагам...
Думаете, это не имеет никакого отношения к обсуждаемой проблеме
"Методы оптимизации"? Ошибаетесь –
имеет! За неимением времени об этом будет вскользь упомянуто ниже...
1.
Почтовый ящик.
В этот раз "почтовый ящик" оказался пустым, хотя "темп скачивания
"
брошюр не стал меньше...
2.
Вместо введения.
Ключевым слов
ом в средствах массовой информации все чаще становятся ободряющие слова "инновация", "оптимизация". В математике и информатике это слово "оптимизация" достаточно конкретно –
оно означает решение задачи нахождения экстремума (минимума или максимума) некотор
ой целевой функции в известной области конечномерного векторного пространства, ограниченной набором линейных и/или нелинейных равенств и/или неравенств. Теорию и методы решения задачи оптимизации изучает научное направление "Математическое программирование
". При этом оптимизационные методы делятся на три большие группы
:
4
аналитические методы;
численные методы;
графические методы.
Как вы уже, наверное, догадались, в настоящей брошюре речь пойдет о численных методах (иначе, зачем бы нужны были эти монстры -
компьютеры) и, значит, прежде всего, о математическом программировании. Это огромная область математики, разрабатывающая общую теорию и численные методы (алгоритмы и программы) решения многомерных задач с ограничениями. В отличие от классической математики
, математическое программирование занимается математическими методами решения задач нахождения наилучших вариантов из всех возможных
. В процессе проектирования систем перед инженерами обычно ставится задача поиска наилучших, в некотором смысле, структуры о
бъекта и/или значений его параметров. Если оптимизация связана с задачей расчета оптимальных значений параметров при заданной структуре объекта, то говорят о параметрической оптимизации. Как правило, именно эту задачу подразумевают СМИ (нельзя трогать Конс
титуционный строй). Задача выбора оптимальной структуры объекта является структурной оптимизацией. Здесь следует заметить, что существует большой класс структурно неустойчивых систем, например, большая экономическая система, сложившаяся в стране. Понятно, что без оптимизации структуры подбор параметров практически лишен смысла...
На протяжении всей своей истории люди при необходимости принимать решения прибегали к сложным ритуалам. Они устраивали торжественные церемонии, приносили в жертву животных, гадали
по звездам и следили за полетом птиц. Они полагались на народные приметы и старались следовать примитивным правилам, облегчающим им трудную задачу принятия решений. В настоящее время для принятия решения используют новый и, по
-
видимому, более научный "рит
уал", основанный на применении компьютеров. Даже не буду пытаться привести строгую классификацию этих научных "ритуалов"
. Сделаю только несколько замечаний для дальнейших изысканий.
Для обозначения научного "ритуала" часто используется термин "математичес
кое программирование", который предложил Роберт Дорфманом приблизительно в 1950 году.
Роберт Дорфман
(Robert Dorfman; 1916 —
2002) —
американский экономист.
Теперь термин объединяет линейное программирование, целочисленное программирование, выпуклое прогр
аммирование, нелинейное программирование и программирование при наличии неопределенности. Нелинейное программирование имеет дело с оптимизацией нелинейных функций при линейных и
\
или нелинейных ограничениях... Однако, многие классические алгоритмы в список не попали. Попытка добавить в эту кашу немного масла не удалась. Но, если очень коротко, то должны быть рассмотрены следующие оптимизационные методы:
Одномерные:
Метод золотого сечения (д
ихотомия
)
, Метод парабол, Перебор по сетке, Метод Фибоначчи, Троичный
поиск.
5
Прямые методы:
Метод Гаусса, Метод Нелдера
—
Мида, Метод Хука
-
Дживса, Метод конфигураций, Метод Розенброка (метод вращающихся координат).
Первого порядка:
Градиентный спуск, Метод Зойтендейка, Покоординатный спуск, Метод сопряжённых градиентов, Квази
ньютоновские методы, Алгоритм Левенберга
—
Марквардта.
Второго порядка:
Метод Ньютона, Метод Ньютона
—
Рафсона.
Стохастические:
Метод Монте
-
Карло, Имитация отжига, Эволюционные алгоритмы, Дифференциальная эволюция, Муравьиный алгоритм, Метод роя частиц.
Методы
линейного программирования:
Симплекс
-
метод, Алгоритм Гомори, Метод эллипсоидов, Метод потенциалов.
Методы нелинейного программирования:
Последовательное квадратичное программирование.
Должны
-
то, должны, но за неимением времени рассмотрим только некоторые
из них. Так что –
простите! И не пугайтесь слова "Экстремум"
(лат. extremum —
крайний) в математике —
максимальное или минимальное значение функции на заданном множестве. Точка, в которой достигается экстремум, называется точкой экстремума. Соответственно
, если достигается минимум —
точка экстремума называется точкой минимума, а если максимум —
точкой максимума. В математическом анализе выделяют также понятие локальный экстремум (соответственно минимум или максимум). Наверное, уже заметили, что все методы можно классифицировать по следующим признакам:
1.1. методы решения линейных задач
1.2. методы решения нелинейных задач
2.1. методы одномерной оптимизации
2.2. методы многомерная оптимизация
3.1. методы оптимизация без ограничений
3.2. методы оптимизация с ограничениями
Научным методом "случайного тыка"
было получено, что первыми будут рассмотрены методы многомерной оптимизации линейных задач с ограничениями. Ну, не спорить же с наукой!
3.
Симплекс
-
метод решения задачи линейного программирования
.
Симпле
кс
-
метод представляет собой алгоритм решения оптимизационной задачи линейного программирования. Он применяется непосредственно к общей задаче линейного программирования сформулированной в стандартной форме:
найти набор (наборы) действительных чисел X = (X(
1), Х(2),... , X(N),
доставляющий экстремум (максимум или чаще минимум) линейной целевой функции
L(X) = C(1)*X(1) + C(2)*X(2) +...+ C(N)*X(N)
при условии неотрицательности переменных X(1), X(2),..., X(N) >= 0
,
и ограничений типа линейных равенств:
A(1,1)*X
(1) + A(1,2)*X(2) + A(1,N)*X(N) = B(1),
A(2,1)*X(1) + A(2,2)*X(2) + A(2,N)*X(N) = B(2),
...
A(M,1)*X(1) + A(M,2)*X(2) + A(M,N)*X(N) = B(M).
6
Метод был разработан советским математиком Канторовичем Л.В. еще в 1937 году. В 1939 вышла брошюра Канторовича "Ма
тематические методы организации и планирования производства"
, в которой зафиксировано открытие линейного программирования -
направления, оказавшего большое влияние на развитие экономической науки. В этой работе впервые давалась математическая постановка пр
оизводственных задач оптимального планирования и предлагались эффективные методы их решения. Во введении Канторович писал: "Существуют два пути повышения эффективности работы цеха, предприятия и целой отрасли промышленности. Один путь —
это различные улучш
ения в технике, т.е. новые приспособления в отдельном станке, изменение технологического процесса, нахождение новых, лучших видов сырья. Другой путь, пока гораздо меньше используемый,
—
это улучшения в организации производства и планировании"
... Тем самы
м идея оптимальности в экономике была поставлена на прочный научный фундамент... Но линейное программирование не используется в управлении экономикой до сих пор и это представляется огромным недостатком. Было бы логичнее, если его применение началось где
-
н
ибудь около 1758 года, когда экономисты впервые начали математически описывать экономические системы. Простейшие примеры моделей линейного программирования можно найти уже в Экономической таблице французского экономиста Франсуа Кенэ
(1694 –
1774), который пытался определить взаимосвязи между землевладельцами, крестьянами и ремесленниками... Позже можно обнаружить, что Леон Вальрас
(1834 –
1910; ему принадлежит формулировка критерия рыночного равновесия: спрос равен предложению) разработал в 1874 сложную мат
ематическую модель, частью которой являлись постоянные технологические коэффициенты. Однако до 1930 года все эти работы носили разрозненный характер, и в области использования моделей линейного типа было сделано очень мало.
Канторович считал необходимым п
родолжать исследования в следующих направлениях:
1.
развитие алгоритмов линейного программирования и их конкретизация для отдельных классов задач;
2.
обобщение предложенных методов с целью изучения более широких классов экстремальных задач с ограничениями, включ
ая нелинейные задачи и задачи в функциональных пространствах; приложение таких методов к экстремальным задачам математики, механики и техники;
3.
распространение новых методов экономического анализа отдельных производственных задач на общие экономические сист
емы; приложение этих методов к задачам планирования и анализа структуры экономических показателей на уровне отрасли, региона и народного хозяйства в целом.
Многие экономисты не понимали идей математика Канторовича, а ортодоксы заявляли, что они несовмести
мы с марксизмом и вообще его опровергают. Алексей Николаевич Косыгин
(председатель Совмина СССР, 1964 —
1980) предлагал провести, как сейчас принято говорить, модернизацию СССР, "научная поддержка" 7
была поручена Леониду Витальевичу. Концепция реформ была и
зложена в докладе на Пленуме ЦК КПСС в 1965 году под названием "об улучшении управления промышленностью, совершенствовании планирования". Но тогда "власть" не нуждалась в укреплении социализма и реформа не состоялась. Видимо, страна ждала еще не оперившихс
я реформаторов
-
революционеров..
. и вот настало время и они пришли...
Человек разумный всегда был, есть и будет человеком хозяйствующим.
Практическая экономика для каждого из нас и наших предков —
это арена здравого смысла. Здравый смысл представляет собой
особую способность человека к мгновенным оценочным суждениям. Понимание выше здравого смысла и проявляется как осознанная адаптивность поведения. Понимание не наследуется и, стало быть, не принадлежит к числу врожденных свойств. Уникальной особенностью че
ловека является способность пониманием делиться, превращая оценки в материальные и идеальные артефакты.
Для детального ознакомления с предметом могу посоветовать сборник:
http://www.ma
th.nsc.ru/LBRT/g2/english/ssk/selecta.pdf
Канторович Л.В. Математико
-
экономические работы. Новосибирск, 2011, 760 с.
Достаточно легко читается книга:
Дж. Данциг. Линейное программирование, его применения и обобщения. М.: 1966.
Джордж Бернард Данциг (Ge
orge Bernard Dantzig; 1914 —
2005) —
американский математик, известен как разработчик симплексного алгоритма, применяемого в решениях задач симплекс
-
методом. Наряду с советским математиком Леонидом Канторовичем, считается основоположником линейного програм
мирования.
Монография Данцига удачно сочетает в себе предельно элементарное изложение основных, исходных вопросов линейного программирования, которое будет доступно даже совсем неискушенному в математике читателю, с главами, посвященными таким глубоким и м
атематически тонким теориям, как принцип разложения или дискретное (целочисленное) программирование. Эту книгу можно, с известным основанием, считать своего рода энциклопедией
линейного программирования, в которой содержится в той или иной форме описание б
ольшинства основных вопросов, относящихся к этой дисциплине.
Для неподготовленного читателя можно порекомендовать книгу:
Гасс С. Путешествие в Страну Линейного Программирования. М., 1973. 176 стр. ил.
Почему самые разные специалисты вынуждены прибегать к математическим методам оптимального управления и, в частности, к линейному программированию? Как от сугубо практической задачи перейти к ее математической модели? Как соотносится эта модель с реальной действительностью? Каковы возникающие при этой трудност
и? На все эти вопросы в доступной н занимательной форме отвечает в настоящей книге крупный американский ученый.
8
Возможны варианты постановки приведенной выше стандартной формы задачи линейного программирования. Например, вместо равенств в ограничениях могу
т быть записаны неравенства.
Это легко исправить, внося в задачу дополнительные (вспомогательные) переменные. Они вносятся также и в целевую функцию с коэффициентами, равными нулю.
Будем решать пример из книги (кстати, содержит хорошие примеры
):
Лунгу К. Н. Линейное программирование. Руководство к решению задач. —
М.: 2005.
Стр. 45: § 4. Симплексный метод решения задач линейного программирования
Пример 4.7.
Целевая функция:
L(X) = 1*x(1) -
1*x(2) –
3*x(3) -
> min
Смешанные ограничения:
2*x(1) -
1*x(2) + 1*x
(3) <= 3
4*x(1) -
2*x(2) + 1*x(3) >= -
6
3*x(1) + 1*x(3) <= 15
В стандартной форме с неотрицательными дополнительными переменными ограничениями задача принимают вид:
+ 1*x(1) -
1*x(2) –
3*x(3) + 0*x(4) + 0*x(5) + 0*x(6) -
> min
+ 2*x(1) -
1*x(2) + 1*x(3) + 1*x(4) + 0*x(5) + 0*x(6) = 3
+ 4*x(1) -
2*x(2) + 1*x(3) + 0*x(4) -
1*x(5) + 0*x(6) = -
6
+ 3*x(1) -
0*x(2) + 1*x(3) + 0*x(4) + 0*x(5) + 1*x(6) = 15
Текст программы на языке FreeBASIC приведен ниже:
' P R O G R A M "Extremum_ESML"
' 08.01.88
'
-------------------------------------------------------------
' Решение задачи линейного программирования
'
-------------------------------------------------------------
'
' Для решения задачи и
спользуется симплексный метод (метод
' последовательного улучшения плана). В качестве примера
' рассматривается задача:
' Целевая функция:
' Y = L(X) = 1*x(1) -
1*x(2) –
3*x(3) -
> min
' Смешанные ограничения:
' 2*x(1) -
1*x(2) + 1*x(3) <= 3
' 4*x(1) -
2*x(2) + 1*x(3) >= -
6
' 3*x(1) + 1*x(3) <= 15
' кроме того: X1 >= 0, X2 >= 0, X3 >= 0.
'
' Ввод исходных данных из блока данных.
' Вывод результатов на экран или в файл.
' Lmin = L(1,11,12) = -
46
'
#Lang "FB" ' режим совм
естимости с FreeBASIC
' Описание используемых подпрограмм (процедур)
Declare Sub ESML(N As Integer, M As Integer, V As Single,_
C() As Single, E() As Single, X() As Single,_
ByRef Y As Single, ER As Integer)
9
Data 3, 3, -
1 ' П
араметры задачи
Data 1, -
1, -
3
Data 2, -
1, 1, 1, 3
Data 4, -
2, 1, -
1, -
6
Data 3, 0, 1, 1, 15
'
Dim As Single C(8), E(10, 10), X(18), Y, V
Dim As Integer N, M, I, J, ER
Read N, M, V ' Ввод параметров задачи
For J = 1 To N: Read C(J)
: Next J
For I = 1 To M
For J = 1 To N+2: Read E(I, J): Next J
Next I
'
ESML(N, M, V, C(), E(), X(), Y, ER)
'
If ER = 0 Then
' Open "ESML.res" For Output As #2
' Open "CON" For Output As #2
Print " ERR = "; ER: Print
Print " Ex
tremum = "; Y: Print
Print " Actual X(I):": Print
For J = 1 To N
Print Using " X(#) = ###.#####"; J; X(J)
Next J
' Close #2
Else
If ER = -
1 Then Print " Solution does NOT exist"
If ER = -
2 Then Print " Solu
tion is NOT limited"
End If
'
Sleep
'
Sub ESML(N As Integer, M As Integer, V As Single,_
C() As Single, E() As Single, X() As Single,_
ByRef Y As Single, ER As Integer)
'
-------------------------------------------------------------
' Симплексный метод решения задачи линейного программирования
'
-------------------------------------------------------------
'
' N -
количество неизвестных переменных;
' M -
количество линейных ограничений;
' V -
вариант поиска: mi
n (
-
1) / max (+1);
' C() -
коэффициенты уравнения функции;
' E() -
коэффициенты уравнений ограничений;
' X() -
значения переменных в точке экстремума;
' Y -
значение функции в точке экстремума;
' ER -
код завершения по
дпрограммы (0 -
норм.);
' A() -
рабочий массив.
'
Dim As Single A(M+2,N+M+1)
Dim As Integer I, J, K, R, S, T
10
ER = 0 ' Пока ошибок нет
K = N + M + 1
For J = 1 To N ' Заполнение матрицы
A(1,J) = -
V*C(J): X(J) = 1
Next J
Fo
r I = 2 To M+1
For J = 1 To N+2: A(I,J) = E(I
-
1,J): Next J
A(I,K) = A(I,N+2): A(I,N+2) = 0
If I <> 2 Then A(I,N+I
-
1) = A(I,N+1): A(I,N+1) = 0
Next I
'
R = 1
For I = 2 TO M+1
If A(I, N+I
-
1) <> -
1 Then
X(N+I
-
1) = 0
Else
X(N+I
-
1) = 1
R = M+2
For J = 1 To N+M
A(R,J) = A(R,J) -
A(I,J)
Next J
EndIf
Next I
'
Do
S = 1: T = 1
For J = 2 To N+M
If A(R,J) >= A(R,S) Then S = J
If A(R,J) < A(R,T) Then T = J
Next J
Print R, S, T
'
If A(R,T) < 0 Then
S = 1
For I = 2 To M+1
If A(I,T) > 0 Then
Y = A(I,K)/A(I,T)
If S = 1 Then
S = I
Else
If Y < A(S,K)/A(S,T) Then S = I
EndIf
EndIf
Next I
If S = 1 Then ER = -
2: Exit Do
For J = 1 To N+M
If X(J
) <> 1 Then If A(S,J) = 1 Then Exit For Next J
X(J) = 1: X(T) = 0: Y = A(S,T)
For J = 1 To K: A(S,J) = A(S,J)/Y: Next J
For I = 1 To M+2
If I <> S Then
11
Y = A(I,T)
For J = 1 To K
A(I,J) = A(I,J) -
Y*A(S,J)
Next J
EndIf
Next I
Else
If R = 1 Then
For J = 1 To N
If X(J) <> 0 Then
X(J) = 0
Else
For I = 2 To M+1
If A(I,J) = 1 Then Exit For
Next I
X(J) = A(I,K)
EndIf
Next J
Y = V*A(1,K)
Exit Do Else
If A(R,S) > 0.0001 Then ER = -
1: Exit Do
R = 1
EndIf
EndIf
Loop
'
End Sub
4.
Методы одномерной оптимизации
.
В этом разделе рассмотрим несколько простых численных процедур, непосредственно локализирующих минимум функции F(X). С помощью численных методов непосредственно ищут минимум функции F(X) в некотором интервале A < X < B, в котором, как предполагается, лежит минимум, выборочно вычисляя значения функции в точках данного интервала. Иногда это единственно возможная стратегия поиска.
4.1.
Метод Ньютона (поиск точек перегиба функции)
Для функций одной переменной классический
подход при поиске значений X в точках перегиба (экстремума) функции F(X) состоит в решении уравнения:
F'(X) = 0.
Решить такое уравнение не всегда просто. Поэтому часто используют численный метод его решения (конечно, после построения графика рассматриваем
ой функции). Приблизительный эскиз кривой Y = F'(X) позволит получить приближенное решение. Если можно найти два значения аргумента A и B таких, что F'(A) и F'(B) имеют противоположные знаки, то тогда, в силу очевидных предположений о непрерывности функции
, будет существовать корень уравнения F'(X) в указанном интервале.
12
Ниже приводится текст программы на языке FreeBASIC. Не думаю, что она будет слишком популярной, но для разминки –
вполне пригодна:
' P R O G R A M "Peregib_Fun"
' 08.01.88
'
-------------------------------------------------------------
' Программа поиска точек перегиба функции F(x)
'
-------------------------------------------------------------
'
' Ординарная переработка подпрогр
аммы на языке QBasic
' Рассматривается функция: F(X) = X^2/2
-
Sin(X)
' E -
требуемая точность расчета = 0.00001
' Z -
начальное значение аргумента = 0.5
' Результаты вычисления (последовательные аппроксимации): ' 0.5000000 0.7552225
' 0.7552225 0.7391416
' 0.7391416 0.7390851
' 0.7390851 0.7390851
' Minimum is -
0.4004886 at the point 0.7390851
#Lang "FB" ' режим совместимости с FreeBASIC
Function Fun(X As Single) As Single
' Вычисление функции F(X) арг
умента X
Fun = X*X/2
-
Sin(X)
End Function
'
Function F1p(X As Single) As Single
' Вычисление первой призводной функции F(X) аргумента X
F1p = X
-
Cos(X)
End Function
'
Function F2p(X As Single) As Single
' Вычисление второй производной функции F(X) аргумент
а X
F2p = 1+Sin(X)
End Function
'
Dim As Single E, Z, X
Print "Search for points of inflection functions F(x)"
Data 0.00001, 0.5
Read E, Z ' Ввод исходных данных Print " ": Print " Successive Approximation:"
Do
X = Z
Z = X
-
F1p(X)/F2p(X)
Print X
, Z
Loop While Abs(Z
-
X)> E
'
Print " "
If F2p(X)>0 Then Print "Minimum is ";Fun(X);" at the point"; X
If F2p(X)<0 Then Print "Maximum is ";Fun(X);" at the point"; X
If F2p(X)= 0 Then Print "It rarely happens" ' редкий
случай
!
Sleep
13
4.2.
Поиск минимума ме
тодом Фибоначчи
Предположим, что нужно определить минимум как можно точнее, т.е. с наименьшим возможным интервалом неопределенности, но при этом можно выполнить только N вычислений функции. Как следует выбрать N точек, в которых вычисляется функция? С пер
вого взгляда кажется ясным, что не следует искать решение для всех точек, получаемых в результате эксперимента. Напротив, надо попытаться сделать так, чтобы значения функции, полученные в предыдущих экспериментах, определяли положение последующих точек. Де
йствительно, вычисленные значения функции уже подсказывают информацию о самой функции, которую можно использовать в дальнейшем поиске.
Предположим, что имеется интервал неопределенности (X1,X3) и известно значение функции F(X2) внутри этого интервала. Есл
и можно вычислить функцию всего один раз в точке X4, то где следует поместить эту точку, для того чтобы получить наименьший возможный интервал неопределенности? Можно, например, выбрать точку X4 таким образом, чтобы минимизировать наибольшую из длин (X3 –
X4) и (X2 –
X1). Достигнуть этого можно, сделав длины (X3 -
X4) и (X2 –
X1) равными т.е. поместив X4 внутри интервала симметрично относительно точки X2, уже лежащей внутри интервала.
Если окажется, что можно выполнить еще одно вычисление функции, то следу
ет применить описанную процедуру к интервалу (X1, X2), в котором уже есть значение функции, вычисленное в точке X4, или к интервалу (X4, X3), в котором уже есть значение функции, вычисленное в точке X2. Следовательно, стратегия ясна с самого начала. Нужно поместить следующую точку внутри интервала неопределенности симметрично относительно уже находящейся там точке... Алгоритм выбора элементов последовательности {XK} называют стратегией поиска. При заданном количестве вычислений функции F(X) оптимальной явля
ется стратегия, которая приводит к наименьшему интервалу неопределенности. Установлено, что стратегия поиска минимума будет оптимальной, если для построения последовательности {XK} использовать числа Фибоначчи. Об этом ученом уже было упомянуто в брошюре "
Фрагмент 3. Некое техническое устройство". Итерационная процедура, реализующая указанную стратегию, называется методом Фибоначчи. Метод "Золотого сечения" незначительно уступает оптимальному методу Фибоначчи по времени поиска минимума. Установлено, что п
ри одинаковом числе итераций метод "Золотого сечения" приводит к отрезку неопределенности, лишь на 17% большему, чем метод Фибоначчи. Метод золотого сечения иногда сочетают с одним из методов аппроксимации минимизируемой функции, которую в области минимума
обычно приближают параболой, минимум которой определяется по аналитической формуле. Программа на языке Фортран, реализующая очень похожий алгоритм, приведена ниже, в разделе 4.3.
Текст программы на языке FreeBASIC для двух вариантов расчета приведен ниже
:
' P R O G R A M "Extremum_Fibonacci"
' 10.10.1996
'
-------------------------------------------------------------
' Поиск экстремума функции методом Фибоначчи
'
------------------------------------
-------------------------
14
' Ординарная переработка подпрограммы на языке QBasic
' В программе производится поиск отрезка (А, В), содержащего
' точку минимума унимодальной функции F(X), за N вычислений.
' Функция F(X) вычисляется в подпрограмме
-
функ
ции.
' Необходимые числа Фибоначчи вычисляются при входе.
' Для выбора варианта примера необходимо ввести NV = 1/2
' Примеры для программы:
' Пример 1:
' F(X) = 2*X^2 -
Exp(X) ' исследуемая функция ' Data 10, 0 ' количество вычислений, эпсило
н
' Data 0, 1 ' границы интервала
' Результаты расчета:
' конечный интервал: 0.3595500 0.3707850 ' значение функции: -
1.1741321
' Пример 2:
' F(X) = X*X*X*X
-
14*X*X*X+60*X*X
-
70*X ' исследуемая функция
' Data 20, 0 ' количество вычислений, эпсилон
' Data 0, 2 ' границы интервала
' Результаты расчета:
' конечный интервал: 0.7808039 0.7809753 ' значение функции: -
24.3696022 #Lang "FB" ' режим совместимости с FreeBASIC
Function Fun(X As Single, NV As Integer) As Single
' Вычисле
ние функции F(X) аргумента X
If NV = 1 Then Fun = 2*X^2 -
Exp(X)
If NV = 2 Then Fun = X*X*X*X
-
14*X*X*X+60*X*X
-
70*X End Function
'
Dim As Integer NV, N, K, I
Dim As Single E, A, B, X, X1, X2, X3, X4, F2, F4
Print
Print " === Extremum Fibonacci ==="
Print "
Enter number of the option ";: Input NV
If NV = 1 Then
' данные для варианта 1
N = 10: E = 0 ' количество вычислений, эпсилон
A = 0: B = 1 ' границы интервала
EndIf
If NV = 2 Then
' данные для варианта 2
N = 20: E = 0 ' количество выч
ислений, эпсилон
A = 0: B = 2 ' границы интервала
EndIf
Dim As Single F(N)
F(0) = 1: F(1) = 1
For I = 2 To N
F(I) = F(I
-
1)+F(I
-
2)
Next I
X1 = A: X2 = A+((B
-
A)*F(N
-
1)+E*(
-
1)^N)/F(N): X3 = B
X = X2: F2 = Fun(X, NV)
Print " Actual Interval (A,B)"
K = 1: Print Using " ###.###### ###.######"; X1; X3
15
L2: X4 = X1
-
X2+X3
X = X4: F4 = Fun(X, NV)
If F4 > F2 Then GoTo L4
If X2 < X4 Then GoTo L3
X3 = X2: X2 = X4: F2 = F4
Print Using " ###.###### ###.######"; X1; X3
GoTo L6
L3: X1 = X2
: X2 = X4: F2 = F4
Print Using " ###.###### ###.######"; X1; X3
GoTo L6
L4: If X2 < X4 Then GoTo L5
X1 = X4: Print Using " ###.###### ###.######"; X1; X3
GoTo L6
L5: X3 = X4: Print Using " ###.###### ###.######"; X1; X3
L6: K = K + 1
If K <= N Then GoTo L2
Print
Print Using " Interval: A = ###.####### B = ###.#######";_
X1; X3 Print Using " Min F(X) = ###.#######"; F2
Sleep
4.3.
Поиск минимума методом Брента
Вот и добрались до моей любимой книжки:
Фо
рсайт Дж., Малькольм М., Моулер К., Машинные методы математических вычислений.
М., Мир, 1980, 280 с., ил.
На странице 197 находим подходящий материал:
8.1. Одномерная оптимизация
Задача состоит в том, чтобы найти множество абсцисс X1, X2,..., Xk, в которы
х вычисляется функция, такое, что оптимальное значение F() лежит при некотором i в интервале X(i
-
1) <= Xo <= X(i+1). Такой интервал называется интервалом неопределенности. Алгоритм выбора абсцисс X(i) где i = 1,..., k называется планом поиска. Если извест
но только то, что F() в интервале унимодальна, то какова оптимальная стратегия для нахождения X? При заданном количестве вычислений функции оптимальным планом поиска будет тот, который приводит к наименьшему интервалу неопределенности... Заманчиво использо
вать метод "Золотого сечения" (он аналогичен методу бисекции для нахождения действительноrо нуля функции) -
он с гарантией работает даже в самом худшем случае и что платой за эту гарантию является медленная, всего лишь линейная сходимость. Метод "Золотого сечения" не извлекает никаких выгод из возможной гладкости функции F{).
Если известно, что функция F() имеет непрерывные производные и можно начать достаточно близко к Х, то можно сконструировать следующий итерационный процесс: начинаем с трех произвольн
ых действительных чисел V1, V2, V3. В общем случае это могут быть числа V(k
-
2), V(k
-
1) и V(k). Пусть V(i+1) -
абсцисса вершины параболы (с вертикальной осью), проходящей через точки (Vi, F(Vi)) где i = k
-
2, k
-
1, k. Продолжаем итерации с числами V(k
-
1), V(k
) и V(k+1). Этот процесс называется последовательной параболической интерполяцией. Можно доказать, что для достаточно хороших 16
начальных приближений итерации сходятся со скоростью порядка ~1.324..., при условии, что F"(X) > 0 в точке минимума (т.е. F"(X) им
еет простой нуль).
В книге Брент
(Richard Brent, 1973, Algorithms for Minimization) предложен алгоритм нахождения минимума, в основе которого комбинация метода "Золотого сечения" и последовательной параболической интерполяции. Алгоритм этот совершенно ана
логичен методу нахождения нуля в программе "Zeroin_Fors.bas" (см. "Фрагмент 8."), который использует комбинацию бисекции и обратной квадратичной интерполяции. Для приводимой ниже программы "FMin_Fors.bas" выбрана тестовая функция F(X) = X^3 –
2*X –
5
, та ж
е, что и для указанной выше программы.
' P R O G R A M "FMin_Fors"
' 28.01.2013
'
-------------------------------------------------------------
' Нахождение минимума функции методом Брента
'
--------
-----------------------------------------------------
' Ординарная переработка программ из книги:
' Форсайт Дж., Малькольм М., Моулер К., Машинные методы
' математических вычислений. М., Мир, 1980, 280 с., ил.
#Lang "FB" ' режим FreeBASIC
-
совместимо
сти
'
' Иллюстрирующая подпрограмма
-
функция
Function F(X As Single) As Single
' X -
аргумент функции;
' F -
значение функции.
F = X*(X*X -
2) -
5
End Function
'
Function FMIN(AX As Single, BX As Single, TOL As single) As Single
' Подпрограмма вычисляет приближение Х к точке, где функция
' F(X) достигает минимума на интервале (АХ, ВХ)
'
' Входные данные:
' АХ –
левая граница исходного интервала
' ВХ –
правая граница исходного интервала
' F –
подпрограмма
-
функци
я, которая вычисляет F(X)
' для любого Х в интервале (АX, ВХ)
' TOL –
желаемая длина интервала неопределенности
' конечного результата ( >= 0)
'
' Выходные данные:
' FMIN -
абсцисса, аппроксимирующая точку, где функция F
' достигает минимума
'
' Метод использует комбинацию поиска золотого сечения и
' последовательной параболической интерполяции. Сходимость
' никогда не бывает значительно хуже, чем для поиска методом
' Фибоначчи. Если функция F() имеет непр
ерывную вторую
' производную, положительную в точке минимума (не
' совпадающей ни с точкой AX, ни с точкой BX), то сходимость
' сверхлинейная и обычно имеет порядок примерно 1.324...
17
' Функция F() никогда не вычисляется в двух точках,
' отстоя
щих друг от друга менее чем EPS*ABS(X)+(TOL/3), где
' EPS примерно равно квадратному корню из относительной
' машинной
точности
.
'
Dim As Single A, B, C, D, E, EPS, XM, P, Q, R, TOL1, TOL2, U, V, W
Dim As Single FU, FV, FW, FX, X, Sign
'
' Константа C есть возведенная в квадрат величина,
' обратная "Золотому сечению"
C = 0.5*(3.0
-
Sqr(5.0))
'
' Вычисление значения EPS
EPS = 1.0
L10: EPS = EPS/2.0
TOL1 = 1.0 + EPS
If TOL1 > 1.0 Then GoTo L10
EPS = Sqr(EPS)
'
' Сохранение (
приравнивание) начальных значений интервала
A = AX: B = BX
V = A+C*(B
-
A): W = V
X = V
E = 0.0
FX = F(X)
FV = FX: FW = FX
'
' Начало основного цикла вычислений
L20: XM = 0.5*(A+B)
TOL1 = EPS*ABS(X)+TOL/3.0
TOL2 = 2.0*TOL1
'
' Проверка критерия окончания
If ABS(X
-
XM) <= (TOL2
-
0.5*(B
-
A)) Then GoTo L90
'
' Проверка необходимости "Золотого сечения"
If ABS(E) <= TOL1 Then GoTo L40
'
' Построение аппроксимирующей параболы
R=(X
-
W)*(FX
-
FV)
Q=(X
-
V)*(FX
-
FW)
P=(X
-
V)*Q
-
(X
-
W)*R
Q=2.0*(Q
-
R)
If Q > 0 Then P = -
P
Q = Abs(Q)
R=E
E=D
'
' Проверка приемлемости параболы
L30: If ABS(P) >= Abs(0.5*Q*R) Then GoTo L40
If P <= Q*(A
-
X) Then GoTo L40
18
If P >= Q*(B
-
X) Then GoTo L40
'
' Шаг параболической интерполяции
D = P/Q: U = X+D
'
' Проверка близости вычисления F() точке AX или BX
If (U
-
A) < TOL2 Then
If (XM
-
X) >= 0 Then Sign = Abs(TOL1)
If (XM
-
X) < 0 Then Sign = -
Abs(TOL1)
D = Sign
EndIf
If (B
-
U) < TOL2 Then
If (XM
-
X) >= 0 Then Sign = Abs(TOL1)
If (XM
-
X) < 0 Then Sign = -
Abs(TOL1)
D = Sign
EndIf
GoTo L50
'
' Очередной шаг "Золотого сечения"
L40: If X >= XM Then E = A
-
X
I
f X < XM Then E = B
-
X
D = C*E
'
' Проверка близости вычисления F() к точке X
L50: If Abs(D) >= TOL1 Then
U = X+D
Else
If D >= 0 Then Sign = Abs(TOL1)
If D < 0 Then Sign = -
Abs(TOL1)
U = U + Sign
EndIf
FU=F
(U)
'
' Присвоение новых значений параметрам
If FU > FX Then GoTo L60
If U >= X Then A = X
If U < X Then B = X
V = W: FV = FW
W = X: FW = FX
X = U: FX = FU
GoTo L20
'
L60: If U < X Then A = U
If U >= X Th
en B = U
If FU <= FW Then GoTo L70
If W = X Then GoTo L70
If FU <= FV Then GoTo L80
If V = X Then GoTo L80
If V = W Then GoTo L80
GoTo L20
'
L70: V = W: FV = FW
19
W = U: FW = FU
GoTo L20
'
L80: V = U
FV
= FU
GoTo L20
'
' Завершение подпрограммы
L90: FMIN=X
End Function
' 28.01.2013
'
-------------------------------------------------------------
' Иллюстрирующая программа для подпрограммы FMIN(
)
'
-------------------------------------------------------------
'
' Используется внешняя функция F(X)
Dim As Single A, B, Z, TOL, X
A = 0.0 ' левая граница интервала
B = 1.0 ' правая граница интервала
TOL = 1.0E
-
5 ' интервал нео
пределенности
Z = FMIN(A, B, TOL)
Open "CON" For Output As #2
'Open "FMin_Fors.res" For Output As #2
Print #2, Using " A = ###.####### B = ###.#######"; A; B
Print #2, Using " Z = ###.#######"; Z
Close #2
Sleep
' Результат вычислений:
' A = 0.0000000
B = 1.0000000
' Z = 0.8165590
5.
Многомерная оптимизация
.
Локальная минимизация функции N переменных настолько важная задача, что алгоритмы для ее решения были изобретены еще более 130 лет назад. В 1845 году Коши предложил метод, называемый тепер
ь методом скорейшего спуска... В дальнейшем теоретический анализ предсказал, что метод Коши в некоторых ситуациях должен сходиться крайне медленно, и опыт подтвердил это во многих практических случаях даже для скромных значений N, равных 2, 3 или 4. Причин
у можно понять, рассматривая при N = 2 функцию F(), для которой линиями уровня являются, например, банальные эллипсы для которых локальные градиенты даже приблизительно не указывает направление к центру эллипсов. В связи с этим предлагаю посмотреть замечат
ельную книгу:
Химмельблау Д. Прикладное нелинейное программирование. М., 1975, 536 с.
Книга посвящена методам оптимального управления системами с нелинейными целевыми функциями. Описаны методы нелинейного программирования как при отсутствии ограничений на
управляющие переменные, так и при наличии ограничений. Рассматриваются такие вопросы, как возможность получения решения, время оптимизации, точность решения и т.д. Конечно, есть и другие достойные книги, они 20
легко доступны в сети, а это хороший повод суще
ственно ограничить "повествовательную часть" настоящей брошюры..
.
5.1.
Численные методы оптимизации нулевого порядка
Далее рассмотрим два метода поиска локального экстремума функций многих переменных и относящиеся к прямым методам (методам безусловной оптимизации нулевого порядка), т.е. опирающихся непосредственно на значения функции:
метод прямого поиска Хука
—
Дживса (Hooke
—
Jeeves),
метод деформируемого многогранника Нелдера
-
Мида (Nelder
-
Mead).
5.1.1. Метод прямого поиска (метод Хука
-
Дживса)
Этот мет
од был разработан в 1961 году, но до сих пор является весьма эффективным и оригинальным. Суть этого метода состоит в следующем. Задаются некоторой начальной точкой X[0]. Изменяя компоненты вектора F[0], обследуют окрестность данной точки, в результате чего
находят направление, в котором происходит уменьшение минимизируемой функции Fun(X). В выбранном направлении осуществляют спуск до тех пор, пока значение функции уменьшается. После того как в данном направлении не удается найти точку с меньшим значением фу
нкции, уменьшают величину шага спуска. Если последовательные дробления шага не приводят к уменьшению функции, от выбранного направления спуска отказываются и осуществляют новое обследование окрестности и т.д.
Достоинством метода прямого поиска является пр
остота его программирования на компьютере. Он не требует знания целевой функции в явном виде, а также легко учитывает ограничения на отдельные переменные, и даже "сложные ограничения" на область поиска. Недостаток метода прямого поиска состоит в том, что в
случае сильно вытянутых, изогнутых или обладающих острыми углами линий уровня целевой функции он может оказаться неспособным обеспечить продвижение к точке минимума. Действительно, в этих случаях, каким бы малым ни был шаг в направлении X1, X2,... из точк
и X’ нельзя получить уменьшения значения целевой функции.
' P R O G R A M "Extremum_Hooke_Jeeves"
' 10.10.1996
'
-------------------------------------------------------------
' Поиск экстремума фун
кции методом Хука и Дживса
'
-------------------------------------------------------------
' Ординарная переработка подпрограммы на языке QBasic
#Lang "FB" ' режим совместимости с FreeBASIC
'
Function Fun(X() As Single, ByRef FE As Integer) As Sing
le
' Подпрограмма вычисления значений функции Fun(X1,X2,X3)
Fun = (X(1)
-
2)^2 + (X(2)
-
5)^2 + (X(3)+2)^4
' Увеличить счетчик количества вычислений функции
FE = FE + 1
End Function
'
21
Dim As Integer I, J, N, FE, PS, BS
Dim As Single H, FI, FB, Z
' Блок и
сходных данных
Data 3 ' Количество переменных
Data 4, -
2, 3 ' Координаты начальной точки
Data 1 ' Исходная величина шага
' Функция вычисляется в подпрограмме
-
функции Fun()
Read N ' Ввод количества переменных
Dim As Single
X(N), B(N), Y(N), P(N), K
For I = 1 To N: Read X(I): Next I ' Ввод
исходных
координат
Read H ' Ввод величины шага
K = H: FE = 0
For I = 1 To N: Y(I) = X(I): P(I) = X(I): B(I) = X(I): Next I
Z = Fun(X(), FE): FI = Z
Print " Sourc
e value F0 = "; Z
For I = 1 To N: Print Using " ###.##### "; X(I);: Next I
Print
PS = 0: BS = 1
' Исследование вокруг базисной точки
J = 1: FB = FI
L2:X(J) = Y(J)+K
Z = Fun(X(), FE)
If Z < FI Then GoTo L3
X(J) = Y(J)
-
K
Z = Fun(X(), FE)
If Z < FI Then GoTo L3
X(J) = Y(J)
GoTo L4
L3:Y(J) = X(J)
L4:Z = Fun(X(), FE)
FI = Z
Print " Exploring Search"; Z
For I = 1 To N: Print Using " ###.##### "; X(I);: Next I
Print
If J = N Then GoTo L5
J = J+1
GoTo L2
L5:If FI < FB
-
0.0000001 Then GoTo L8
' Если функция уменьшилась, то произвести поиск по образцу
If PS = 1 And BS = 0 Then GoTo L6
' ...но если исследование производилось вокруг точки шаблона
' РТ и уменьшение функции не было достигнуто, то изменить
' базис
ную точку в операторе 420...
' в противном случае уменьшить длину шага
GoTo L7
L6:For I = 1 To N: P(I)=B(I): Y(I)=B(I): X(I)=B(I): Next I
Z = Fun(X(), FE): BS = 1: PS = 0
FI = Z: FB = Z
Print " Replacing reference point "; Z
For I = 1 To N
: Print X(I);: Next I: Print
' ...и провести исследование вокруг новой базисной точки
J = 1: GoTo L2
22
L7:K = K/10
Print " Reduce length of step "
If K < 0.00000001 Then GoTo L9
' Если поиск не закончен, то произвести новое
' исследование вокр
уг новой базисной точки
J = 1: GoTo L2
' Произвести поиск по образцу
L8:For I = 1 To N
P(I)=2*Y(I)
-
B(I)
B(I)=Y(I): X(I)=P(I): Y(I)=X(I)
Next I
Z = Fun(X(), FE): FB=FI: PS=1: BS=0: FI=Z
Print " Sample search "; Z
For I = 1 To N: Print
Using " ###.##### "; X(I);: Next I
Print
' ...после этого произвести исследование
' вблизи последней точки образца
J = 1: GoTo L2
L9:Print " Found minimum!"
For I = 1 To N: Print Using " ###.##### "; P(I);: Next I
Print
Print Using " Minimum is equal F(X) = ###.#####"; FB
Print Using " Number of calculations = ###"; FE
Print
Sleep
Приведенная выше программа реализует процедуру поиска минимума. Одной или двух точек бывает недостаточно для определения начальной точки. Первая точк
а всегда должна выбираться осмотрительно. Подпрограмма, вычисляет значение минимизируемой функции:
Fun(X1, X2, X3) = (X1 -
2)^2 + (X2 -
5)^2 + (X3 + 2)^4.
Минимум, очевидно, находится в точке (2, 5, —
2).
Для начальной точки (4, —
2, 3) и начального шага длино
й 1 на экран монитора выводятся промежуточные результаты. По ним можно легко проследить за изменениями в ходе исследований и поиска по образцу. Количество выполненных вычислений функции запоминается в счетчике. Это часто используется как средство сравнения
эффективности различных методов поиска. Чем лучше метод, тем меньше в общем случае требуется вычислений значений функции. Нетрудно дописать подпрограмму
-
функцию операторами проверки условий и присвоения функции искусственного очень большого значения с тем
, чтобы решить задачу оптимизации с ограничениями.
5.1.2. Метод деформируемого многогранника (метод Нелдера
—
Мида)
Данный метод является развитием симплексного метода Спендли, Хекста и Химсворта. Он состоит в том, что для минимизации функции N переменных
Fun(X) в N
-
мерном пространстве строится многогранник, содержащий (N+1) вершину (в двумерном пространстве симплексом является равносторонний треугольник). Очевидно, что каждая вершина соответствует некоторому вектору X. Вычисляются значения целевой функции
Fun(X) в каждой из вершин многогранника, определяются максимальное из этих значений и соответствующая ему вершина Xh. Через эту 23
вершину и центр тяжести остальных вершин проводится проецирующая прямая, на которой находится точка Xq с меньшим значением целе
вой функции, чем в вершине Xh. Затем исключается вершина Xh. Из оставшихся вершин и точки Xq строится новый многогранник, с которым повторяется описанная процедура. В процессе выполнения данных операций многогранник изменяет свои размеры, что и обусловило название метода.
С помощью операции растяжения и сжатия размеры и форма деформируемого многогранника адаптируются к топографии целевой функции. В результате многогранник вытягивается вдоль длинных наклонных поверхностей, изменяет направление в изогнутых в
падинах, сжимается в окрестности минимума, что определяет эффективность рассмотренного метода. В результате получился очень надежный и быстрый метод прямого поиска (не градиентный), оказавшимся одним из самых эффективных при проверке до значений N <=6.
Ни
же приведена программа на языке FreeBASIC, использующая метод Нелдера
-
Мида для минимизации функции Розенброка (о ней будем говорить еще, но в следующих выпусках):
F(X1,X2) = 100*(X(2) –
X(1)^2)^2 + (1 –
X(1))^2
в ней шаг SL = 0,5, а в качестве начального п
риближения выбрана точка (1,5; 2), для коэффициентов деформации используются значения ALFA = 1.0, BETA = 0.5, GAMMA = 2,0.
Фрагмент функции Розенброка
' P R O G R A M "Extremum_Nelder_Mead"
' 10.02.
2013
'
-------------------------------------------------------------
' Поиск экстремума функции методом Нелдера и Мида
'
-------------------------------------------------------------
' Ординарная переработка подпрограммы на языке QBasic
#Lang "FB" '
режим совместимости с FreeBASIC
24
'
Function Fun(X() As Single, ByRef FEX As Integer) As Single
' Подпрограмма вычисления значений функции Розенброка
Fun = 100*(X(2)
-
X(1)^2)^2+(1
-
X(1))^2 ' функции Розенброка
' Нарастить счетчик количества вычислений фун
кции
FEX = FEX + 1 ' нарастить счетчик числа вычислений функции
' Functions EXecution
End Function
' Текст
основной
программы
Dim As Single SL, AL, BE, GA, SIG, S1, S2
Dim As Single Z, FO, FH, FG, FL, FR, FE, FC
Dim As Integer FEX, N, I, J,
H, G, L
Data 2 ' число переменных
Data 1.5, 2 ' начальное приближение
Data 0.5 ' величина шага
Data 1, 0.5, 2 ' рекомендация Нелдера и Мида
Print
Print " === Extremum Nelder Mead ==="
Print
FEX = 0 ' Функция еще не вычис
лялась
Read N ' Ввести число переменных
Dim As Single S(N+1,N) ' матрица
вершин
симплекса
' I -
номер вершины, J -
номер координаты
Dim As Single X(N), XH(N), XG(N), XL(N), XO(N)
Dim As Single XR(N), XC(N), XE(N), F(N+1)
For J =
1 To N
Read S(1,J) ' начальное
приближение
Next J
Read SL ' Step Length -
длина
шага
' Построение первого симплекса вокруг исходной точки
For I = 2 To N+1
For J = 1 To N
If J = I
-
1 Then
S(I,J)= S(1,J)+SL
Else
S(I,J)=S(1,J)
En
d If
Next J
Next I
' ввод
параметров
ALFA, BETA, GAMMA
Read AL, BE, GA
' Вычисление значения функции For I = 1 TO N+1
For J = 1 TO N
X(J) = S(I,J)
Next J
Z = Fun(X(), FEX): F(I) = Z
Next I
' Наибольшее и наименьшее значения функции
'
и соответствующие им координаты точек
M1: FH = -
1.0E+20: FL = 1.0E+20
For I = 1 To N+1
25
If F(I) > FH Then FH = F(I): H = I
If F(I) < FL Then FL = F(I): L = I
Next I
' Следующее за наибольшим значение функции
' и соответствующая ему точка
FG = -
1.0E+20
For I = 1 To N+1
If I = H Then GoTo M2
If F(I) > FG Then FG = F(I): G = I
M2: Next I
For J = 1 To N
XO(J) = 0
For I = 1 To N+1
If I <> H Then XO(J) = XO(J)+S(I,J)
Next I
' цент
р тяжести XO, наибольшая XH, большая XG, наименьшая XL
XO(J) = XO(J)/N
XH(J) = S(H,J): XG(J) = S(G,J): XL(J) = S(L,J)
Next J
For J = 1 To N
X(J) = XO(J)
Next J
Z = Fun(X(), FEX): FO = Z
Print " Calculation of the centr
e of gravity"
' Выполнение отражения
For J = 1 To N
XR(J) = XO(J)+AL*(XO(J)
-
XH(J)): X(J) = XR(J)
Next J Z = Fun(X(), FEX): FR=Z
Print " Performing reflection"; Z
' Если
FR<FL, то
выполнить
растяжение
If FR < FL Then GoTo M3
' Если FR>FL и FR>FG, то проверить значения FR и FH
' в противном случае заменить точку ХН на XR
If FR > FG Then GoTo M6
GoTo M5
' Выполнение растяжения
M3: For J = 1 To N
XE(J) = GA*XR(J)+(1
-
GA)*XO(J)
X(J) = XE(J)
Next J
Z = Fun
(X(), FEX): FE=Z
If FE < FL Then GoTo M4
GoTo M5
M4: For J = 1 To N
S(H,J) = XE(J)
Next J
F(H) = FE
Print " Implementation of the tensile"; Z
' Переход к проверке сходимости процесса
GoTo M8
26
M5: For J = 1 TO N
S(H,J) = XR
(J)
Next J
F(H) = FR
Print " Performing reflection"
GoTo M8
M6: If FR <= FH Then
For J = 1 To N
XH(J) = XR(J)
Next J
F(H) = FR
EndIf
' Далее выполнить сжатие
For J = 1 To N
XC(J) = BE*XH(J)+(1
-
BE)*XO(J)
X(J
) = XC(J)
Next J
Z = Fun(X(), FEX): FC = Z
If FC > FH Then GoTo M7
For J = 1 To N
S(H,J) = XC(J)
Next J F(H) = FC
Print " Implementation of compression "; Z
GoTo M8
' Начало редукции (преобразования) симплекса
M7: For
I = 1 To N+1
For J = 1 To N
S(I,J) = (S(I,J)+XL(J))/2
X(J) = S(I,J)
Next J
Z = Fun(X(), FEX): F(I) = Z
Next I
Print " Implementation of the reduction"
' Проверка сходимости процесса
M8: S1 = 0: S2 = 0
For I = 1 To N+1
S1 = S1+F(I) ' сумма
значений
функции
S2 = S2+F(I)*F(I) ' сумма квадратов значений
Next I SIG = S2
-
S1*S1/(N+1): SIG = SIG/(N+1)
If SIG >= 1E
-
10 Then GoTo M1
' Вывод результатов поиска минимума функции
Print: Print " Minimum at the point:"
For J = 1 To N
Print " X"; J; " ="; XL(J)
Next J: Print " "
Print " Value of minimum ="; FL
Print " Number of calculations ="; FEX
Sleep
27
При выполнении программы минимум найден в точке:
X1 = 1.000634
X2 = 1.00135
7
Значение функции минимума = 1.193644e
-
006
Количество вычислений функции = 108
5.2.
Градиентные численные методы оптимизации
Ниже рассмотрим два метода поиска, в которых наряду со значениями функции используется и ее градиент:
Метод наискорейшего спус
ка
Метод Давидона
-
Флетчера
-
Пауэлла
5.2.1.
Метод наискорейшего спуска
Наиболее известным среди градиентных методов является метод наискорейшего спуска. Идея этого метода проста: поскольку вектор градиента указывает направление наискорейшего возрастания функции, то минимум следует искать в обратном направлении.
Используя приведенную ниже программу, найдем минимум функции:
f(xl,x2,x3) = (X(1)
-
1)^2 + (X(2)
-
3)^2 + 4*(X(3)+5)^2
Минимум, очевидно, равен нулю при X(1) = 1, X(2)= 3, X(3) = -
5.
Поиск проводится из начальной точки (4; -
1; 2) с начальным шагом длиной 4.
' P R O G R A M "Extremum_FastDescent"
' 10.02.2013
'
-------------------------------------------------------------
' Поиск экстремума функц
ии методом наискорейшего спуска
'
-------------------------------------------------------------
' Ординарная переработка подпрограммы на языке QBasic
' Одномерный поиск методом квадратичной интерполяции.
#Lang "FB" ' режим совместимости с FreeBASIC
'
Function Fun(X() As Single, ByRef FEX As Integer) As Single
' Подпрограмма вычисления исследуемой функции
Fun = (X(1)
-
1)^2 + (X(2)
-
3)^2 + 4*(X(3)+5)^2
FEX = FEX + 1 ' нарастить счетчик
End Function
'
Sub Gfn(X() As Single, N As Integer, G() As S
ingle, ByRef G0 As Single)
' Подпрограмма вычисления частных производных
Dim As Integer I
G0 = 0
G(1) = 2*(X(1)
-
1)
G(2) = 2*(X(2)
-
3)
G(3) = 8*(X(3)+5)
For I = 1 To N: G0 = G0+G(I)*G(I): Next I
28
G0 = Sqr(G0)
End Sub
' Основная про
грамма поиска минимума
' Блок исходных данных
Data 3 ' число переменных
Data 4, -
1, 2 ' начальное приближение
Data 4 ' начальный шаг
'
Dim As Integer N, I, J, K, FEX
Dim As Single LL, Z, ZZ, G0, G2, LO, FO, S1, S2, S3, DN, F
Read N ' Ввод количества переменных
Dim As Single X(N), Y(N), G(N), D(N), L(4), FF(4)
' Массивы L(4) и FF(4) нужны для квадратичной интерполяции
' Ввод координат начальной точки
For I = 1 To N: Read X(I): Next I
' Ввод начального шага
Read LL ' Ввод начального шага
Print " *** Current Value ***"
' Заголовок промежуточного вывода
For I = 1 To N: Y(I)=X(I): Next I
Z = Fun(X(), FEX): Gfn(X(), N, G(), G0)
If G0 < 0.000001 Then GoTo M6
M1: For I = 1 To N: D(I) = -
G(I)/G0: Next I
' Задать направление наискорейшего спуска из точки Y(I)
' с координатами, равными координатам начальной точки
L(1) = 0: FF(1) = Z: ZZ = Z
M2: L(3) = LL
For I = 1 To N: X(I)=Y(I)+L(3)*D(I): Next I
Z = Fun(X(), FEX): FF(3)=Z
Gfn(X(), N, G(), G0)
G2 = 0
For I = 1 To N: G2 = G2+G(I)*D(I): Next I
If FF(3) >= FF(1) Or G2 >= 0 Then GoTo M3 ' Удвоить длину шага LL, чтобы "накрыть" минимум
LL = 2*LL: GoTo M2
'
M3: L(2) = LL/2
For I = 1 To N: X(I) = Y(I)+L(2)*D(I): Next I
Z = Fun(X(), FEX): FF(2) = Z
' Выполнить первую квадратичную интерполяцию
L(4)=LL*(FF(2)
-
0.75*FF(1)
-
0.25*FF(3))/(2*FF(2)
-
FF(1)
-
FF(3))
If L(4) < 0 Then Print " Attention!"
For I = 1 TO N: X(I) = Y(I)+L(4)*D(I): Next I
Z = Fun(X(), FEX): FF(4
) = Z
' Упорядочить 4 значения LAMBDA и 4 значения функции ' в порядке убывания значений функции
M4: For J = 1 To 3
For K = J+1 To 4
If FF(J) > FF(K) Then
' Поменять местами FF(J) и FF(K), если они не упорядочены
LO = L(J): L(J) = L(K): L(K) = LO ' LO и FO -
буферы
29
FO = FF(J): FF(J) = FF(K): FF(K) = FO
EndIf
Next K
Next J
' Если заданная точность достигнута, то закончить поиск
' в текущем направлении
If Abs(L(1)
-
L(2)) < 0.00005 Then GoTo M5
' Запомнить три лу
чшие точки
S1 = Sgn(L(2)
-
L(1)): S2 = Sgn(L(3)
-
L(1))
S3 = Sgn(L(4)
-
L(1))
If S1 = S2 And S1 = -
S3 Then L(3)=L(4): FF(3)=FF(4)
' Вторая и последняя интерполяция DN = (L(2)
-
L(3))*FF(1)+(L(3)
-
L(1))*FF(2)+(L(1)
-
L(2))*FF(3)
F = (FF(1)
-
FF(2)
)/(2*DN)
F = F*(L(2)
-
L(3))*(L(3)
-
L(1))
L(4) = (L(1)+L(2))/2+F
For I = 1 To N: X(I) = Y(I)+L(4)*D(I): Next I
Z = Fun(X(), FEX): FF(4)=Z
' Повторить вторую интерполяцию GoTo M4
M5: For I = 1 To N: X(I) = Y(I)+L(1)*D(I)
Y(I) = X(I):
Print " X";I; "="; X(I)
Next I
Z = Fun(X(), FEX): Gfn(X(), N, G(), G0)
Print " F ="; Z: Print
LL = LL/2
If G0 > 0.00001 Then GoTo M1
' Переход к началу итерации поиска из текущей точки
M6: Print: Print
For I = 1 To N: Print " X"
;I;"=";X(I): Next I
Print " Minimum of Function F(X1,X2,...,XN) = "; Z
Print Using " Number of calculations = ###"; FEX
Print
Sleep
Результаты вычисления:
X1 = 1.000003
X2 = 2.999996
X3 =
-
5.000000
Минимум функции
= 1.966782e
-
011
Ко
личество вычислений функции
= 56
5.2.2.
Метод Давидона
-
Флетчера
-
Пауэлла
Первоначально метод был предложен Дэвидоном (Davidon, 1959), а затем развит Флетчером и Пауэллом (Fletcher, Powell, 1963). Метод Дэвидона -
Флетчера -
Пауэлла называют также и мето
дом переменной метрики. Он попадает в общий класс квазиньютоновских процедур... Поиск из начальной точки X0 начинают, взяв в качестве начальной матрицу Н0 (обычно единичную матрицу, хотя в этом случае может подойти любая симметрическая положительно определ
енная матрица). Итерационная процедура заключается в выполнении одномерного поиска в 30
направлении градиента с целью определения "наилучшей" точки. После оценки "качества" точки определяется сопряженное направление и повторяется одномерный поиск. При выполне
нии описываемого алгоритма поиск после первой попытки ведется в тех направлениях, в которых целевая функция в ближайшей окрестности имеет значения, приближающиеся к оптимальному.
В качестве "тестового" используется пример предыдущего раздела –
найти миним
ум функции трех переменных:
f(xl,x2,x3) = (X(1)
-
1)^2 + (X(2)
-
3)^2 + 4*(X(3)+5)^2
из начальной точки (4; -
1; 2) с начальной единичной матрицей. Текст программы на языке FreeBASIC приведен ниже:
' P R O G R A M "Extremum_DavFlePow"
' 10.02.2013
'
-------------------------------------------------------------
' Поиск экстремума функции методом ДФП
'
-------------------------------------------------------------
' Ординарная переработка подпрограммы на яз
ыке QBasic
#Lang "FB" ' режим совместимости с FreeBASIC
'
Function Fun(X() As Single, ByRef FEX As Integer) As Single
' Подпрограмма вычисления исследуемой функции
Fun = (X(1)
-
1)^2 + (X(2)
-
3)^2 + 4*(X(3)+5)^2
FEX = FEX + 1 ' нарастить счетчик
End Function
'
Sub Gfn(X() As Single, N As Integer, G() As Single, ByRef G0 As Single)
' Подпрограмма вычисления частных производных
Dim As Integer I
G0 = 0
G(1) = 2*(X(1)
-
1)
G(2) = 2*(X(2)
-
3)
G(3) = 8*(X(3)+5)
For I = 1 To N: G0 = G0
+G(I)*G(I): Next I
G0 = Sqr(G0)
End Sub
' Основная программа поиска минимума
' Блок исходных данных
Data 3 ' число переменных
Data 4, -
1, 2 ' начальное приближение
' Ввод количества переменных
Dim As Integer N, I, J, K, FEX
Dim As Single
CC, Z, ZZ, FP, FQ, FR, FO, G0, G1, G2, G3, GP, GQ, GR
Dim As Single HH, QX, BB, W, WW, DD, KK, WK, DK
Read N ' Ввод числа переменных
Dim As Single X(N), P(N), Q(N), R(N), D(N), G(N), U(N), V(N), Y(N), M(N)
Dim As Single H(N,N)
' Построение начал
ьной (единичной) матрицы H()
For I = 1 To N
31
For J = 1 To N: H(I,J) = 0: Next J
H(I,I) = 1
Next I
CC = 0: FEX = 0 ' Сброс
всех
счетчиков
' Ввод координат начальной точки
For I = 1 To N: Read X(I): Next I
Print " *** Current Val
ue ***"
' Просто заголовок промежуточного вывода
M1: For I = 1 To N: P(I)=X(I): Y(I)=X(I)
Print " X";I;" ";X(I): Next I
Z = Fun(X(), FEX)
Print " Iteration"; CC; " value "; Z
Print
FP = Z: Gfn(X(), N, G(), G0): G1 = G0
' Градиент за
помнить в U и выбрать начальное направление D
For I = 1 To N
U(I) = G(I): D(I) = 0
For J = 1 To N
D(I) = D(I)
-
H(I,J)*G(J)
Next J Next I M2: GP = 0
For I = 1 To N: GP = GP+G(I)*D(I): Next I
If GP < 0 Then GoTo M3
' Найти начальный шаг и, если необходимо, ' изменить направление спуска на противоположное QX = Abs(2*FP/GP): If QX > 1 Then QX = 1
For I = 1 To N
X(I) = P(I)
-
QX*D(I): P(I) = X(I): Next I
Z = Fun(X(), FEX): FP = Z: Print " Instability?"
Gfn
(X(), N, G(), G0): G1 = G0: GoTo M2 M3: QX = Abs(2*FP/GP): If QX > 1 Then QX = 1
HH = QX
' Найти следующую точку Q
M4: BB = HH
For I=1 To N Q(I)=P(I)+BB*D(I): X(I)=Q(I)
Next I
Z = Fun(X(), FEX): FQ = Z
Gfn(X(), N, G(), G0): G2 =
G0
GQ = 0
For I=1 To N
GQ = GQ+G(I)*D(I)
Next I
If GQ > 0 Or FO > FP Then GoTo M5
' Выполнить кубическую интерполяцию ' ил удвоить шаг, чтобы "накрыть" минимум HH = 2*HH: GoTo M4
M5: ZZ = 3*(FP
-
FQ)/HH: ZZ = ZZ+GP+GQ
WW = ZZ*ZZ
-
GP*GQ: If WW < 0 Then WW = 0
W = Sqr(WW)
DD = HH*(1
-
(GQ+W
-
ZZ)/(GQ
-
GP+2*W))
32
For I=1 To N: X(I)=P(I)+DD*D(I): Next I
Z = Fun(X(), FEX): FR = Z
Gfn(X(), N, G(), G0): G3 = G0
' Поиск градиента в новой точке GR = 0
For I = 1
To N: GR = GR+G(I)*D(I): Next I If Z <= FP And Z <= FQ Then GoTo M7 If GR > 0 Then GoTo M6
HH = HH
-
DD
For I = 1 To N: P(I) = X(I): Next I
FP = Z: GP = GR: G1 = G0: GoTo M5
M6: HH = DD For I = 1 To N: Q(I) = X(I): Next I ' Обн
овление матрицы Н() M7: KK = 0: WK = 0: DK = 0
For I = 1 To N
U(I)=G(I)
-
U(I): V(I)=X(I)
-
Y(I)
Next I
For I = 1 To N: M(I) = 0
For J = 1 To N
M(I) = M(I)+H(I,J)*U(J)
Next J
KK = KK+M(I)*U(I): WK = WK+V(I)*U(I)
DK = DK+V(I
)*V(I)
Next I
If KK = 0 Or WK = 0 Then GoTo M8
For I = 1 To N
For J = 1 To N
H(I,J) = H(I,J)
-
M(I)*M(J)/KK+V(I)*V(J)/WK
Next J
Next I
M8: CC = CC+1
' Проверка критерия завершения вычислений
If SQR(DK) < 0.00005 Or G3 < 0.0
0001 Then GoTo M9
' Переход к началу итерации поиска
GoTo M1
' Переход сюда, если минимум функции найден
M9: Print
Print " Minimum at the point:"
For I = 1 To N
Print " X"; I; " "; X(I)
Next I
Print Using " Minimum is equal F
(X) = ###.#####"; Z
Print Using " Number of calculations = ### ###"; FEX; CC
Print
Sleep
Результаты поиска минимума функции:
X1 = 0.999999
X2 = 3.000000
X3 =
-
5.000000
Минимум функции = 1.966782e
-
011
Количество вычислений функции FEX = 7 CC = 2
33
Честно признаюсь –
не ожидал такого успешного испытания!
Конечно, вместо всего этого великолепия человеческой мысли можно искать экстремум простым перебором по достаточно "густой" координатной сетке. Но вступает в силу "проклятье размерности", де
лающую задачу трудноисполнимой даже на современных вычислительных машинах за приемлемое время. Однако существуют методы случайного поиска использующие различные стратегии генерации случайных точек. Есть даже удачные попытки глобальной оптимизации на основе
гибридизации методов роя частиц, эволюции разума и клональной селекции... Это –
интереснейшая область и надеюсь в ближайшее время уделить ей внимание.
В заключение хочу сказать следующее:
"Буду благодарен за всяческие улучшени
я
и упрощения представленн
ых в этой брошюре программ, а также за внятные описания (для школьников) используемых в них алгоритмов".
На этом позвольте закруглиться.
До новых встреч!
Пишите: eugene_r@mail.ru
Документ
Категория
Информатика
Просмотров
206
Размер файла
599 Кб
Теги
программы, freebasic, алгоритмы
1/--страниц
Пожаловаться на содержимое документа