17428

ТЕОРИЯ ОПТИМИЗАЦИИ. Конспект лекций

Конспект

Информатика, кибернетика и программирование

ТЕОРИЯ ОПТИМИЗАЦИИ Конспект лекций Приведены программы методические указания по изучению курса контрольные вопросы характеристика лабораторных работ и задание на контрольную работу. Методические указания предназначены для студентов заочного отделения со с...

Русский

2013-07-01

1.22 MB

77 чел.

ТЕОРИЯ ОПТИМИЗАЦИИ

Конспект лекций

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

Методические указания предназначены для студентов заочного отделения со сроком обучения 6 лет специальности «190300» а также для студентов дневной формы обучения по специальнстям «190300» и «131000».

Подготовлены к публикации кафедрой______________ по рекомендации методической комиссии факультета _____________________ Санкт- Петербургского Государственного университета аэрокосмического приборостроения.

  1.  
    ЦЕЛИ И ОБЩАЯ ХАРАКТЕРИСТИКА УЧЕБНОЙ ДИСЦИПЛИНЫ

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

        Данная учебная дисциплина ставит своей целью:
              ознакомить пользователя с теоретическими аспектами методов оптимизации и их компьютерной реализацией на алгоритмических языках Паскаль, Си или Фортран;
              показать возможности разработанных методов, программ, программных пакетов и комплексов на примерах решения прикладных задач;
              научить пользоваться разработанными программными продуктами при решении проблем связанных с оптимизацией.


Целью преподавания дисциплины является изучение методов решения задач оптимизации.

В результате изучения дисциплины студент должен знать:

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

Студент должен уметь:

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

Тема №1
Общая характеристика оптимизационных задач

ПОСТАНОВКА ЗАДАЧИ ОПТИМИЗАЦИИ

[ Математическая постановка задач оптимизации. ] [ Критерии оптимальности: частный, аддитивный, мультипликативный, максиминный. ] [ Виды ограничений. ] [ Классификация задач. ] 

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

Оптимизация - целенаправленная деятельность, заключающаяся в получении наилучших результатов при соответствующих условиях.

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

При постановке задачи оптимизации необходимо:

    1. Наличие объекта оптимизации и цели оптимизации. При этом формулировка каждой задачи оптимизации должна требовать экстремального значения лишь одной величины, т.е. одновременно системе не должно приписываться два и более критериев оптимизации, т.к. практически всегда экстремум одного критерия не соответствует экстремуму другого.
    Типичный пример неправильной постановки задачи оптимизации:
      "Получить максимальную производительность при минимальной себестоимости". Ошибка заключается в том, что ставится задача поиска оптимума 2-х величин, противоречащих друг другу по своей сути.
    Правильная постановка задачи могла быть следующая:
       а) получить максимальную производительность при заданной себестоимости;
       б) получить минимальную себестоимость при заданной производительности;
    В первом случае критерий оптимизации - производительность, а во втором - себестоимость.

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

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

    4. Учет ограничений.


Математическая постановка задач оптимизации
Виды ограничений
 

Несмотря на то что прикладные задачи относятся к совершенно разным областям, они имеют общую форму. Все эти задачи можно классифицировать как задачи минимизации вещественнозначной функции f(x) N-мерного векторного аргумента x=(x1, x2,..., xn), компоненты которого удовлетворяют системе уравнений hk(x)=0, набору неравенств gi(x) 0, а также ограничены сверху и снизу, т.е. xi(u) xi xi(l). В последующем изложении функцию f(x) будем называть целевой функцией, уравнения hk(x)=0 - ограничениями в виде равенств, а неравенства gi(x) 0 - ограничениями в виде неравенств. При этом предполагается, что все фигурирующие в задаче функции являются вещественнозначными, а число ограничений конечно.

Задача общего вида:

минимизировать f(x)

при ограничениях

hk(x)=0, k=1, ..., K,

gj(x) 0 j=1, ..., J,

xi(u) xi xi(l), i=1, ..., N

называется задачей оптимизации с ограничениями или задачей условной оптимизации. Задача, в которой нет ограничений, т.е.

J=K=0;

xi(u)= - xi(l) = , i=1, ..., N,

называется оптимизационной задачей без ограничений или задачей безусловной оптимизации.


Критерии оптимальности

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

Критерием оптимальности называется количественная оценка оптимизируемого качества объекта.

          На основании выбранного критерия оптимальности составляется целевая функция, представляющая собой зависимость критерия оптимальности от параметров, влияющих на ее значение. Вид критерия оптимальности или целевой функции определяется конкретной задачей оптимизации.
          Таким образом, задача оптимизации сводится к нахождению экстремума целевой функции.
          Наиболее общей постановкой оптимальной задачи является выражение критерия оптимальности в виде экономической оценки (производительность, себестоимость продукции, прибыль, рентабельность).Однако в частных задачах оптимизации, когда объект является частью технологического процесса, не всегда удается или не всегда целесообразно выделять прямой экономический показатель, который бы полностью характеризовал эффективность работы рассматриваемого объекта. В таких случаях критерием оптимальности может служить технологическая характеристика, косвенно оценивающая экономичность работы агрегата (время контакта, выход продукта, степень превращения, температура). Например устанавливается оптимальный температурный профиль, длительность цикла - "реакция - регенерация". Но в любом случае любой критерий оптимальности имеет экономическую природу.
          Рассмотрим более подробно требования, которые должны предъявляться к критерию оптимальности.
               1. Критерий оптимальности должен выражаться количественно.
               2. Критерий оптимальности должен быть единственным.
               3. Критерий оптимальности должен отражать наиболее существенные стороны процесса.
               4. Желательно чтобы критерий оптимальности имел ясный физический смысл и легко рассчитывался.
          Любой оптимизируемый объект схематично можно представить следующим образом (рис. 5.1) .




Рис. 5.1

На рис. 5.1 обозначено:
Y - выходы объекта;
X - контролируемые входные параметры;
U - регулируемые входные параметры управляющие параметры;
Z - неконтролируемые воздействия;

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

R=R(X1, X2,...,XN, Y1,Y2,...,YN, U1,U2,..., UN);      (5.1)

Так как Y=f (U), то при фиксированных Х можно записать:
     (5.2)
При этом всякое изменение значений управляющих параметров двояко сказывается на величине R:
     - прямо, т.к управляющие параметры непосредственно входят в выражение критерия оптимизации;
     - косвенно - через изменение выходных параметров процесса, которые зависят от управляющих.
Если же случайные возмущения достаточно велики и их необходимо учитывать,то следует применять экспериментально - статистические методы, которые позволят получить модель объекта в виде функции
Y=j(Xi,Ui)      (5.3)
которая справедлива только для изученной локальной области. Тогда критерий оптимальности примет следующий вид:
     (5.4)
          В принципе, для оптимизации вместо математической модели можно использовать и сам объект, однако оптимизация опытным путем имеет ряд существенных недостатков:
               а) необходим реальный объект;
              б) необходимо изменять технологический режим в значительных пределах, что не всегда возможно;
               в) длительность испытаний и сложность обработки данных. Наличие математической модели (при условии, что она достаточно надежно описывает процесс) позволяет значительно проще решить задачу оптимизации аналитическим либо численным методами.

          В задачах оптимизации различают простые и сложные критерии оптимизации.
          Критерий оптимальности называется
простым, если требуется определить экстремум целевой функции без задания условий на какие-либо другие величины. Такие критерии обычно используются при решении частных задач оптимизации (например,определение максимальной концентрации целевого продукта,оптимального времени пребывания реакционной смеси в аппарате и т.п.).
          Критерий оптимальности называется
сложным,если необходимо установить экстремум целевой функции при некоторых условиях,которые накладываются на ряд других величин (например, определение максимальной производительности при заданной себестоимости,определение оптимальной температуры при ограничениях по термостойкости катализатора и др.).
          Процедура решения задачи оптимизации обязательно включает, помимо выбора управляющих параметров,еще и установление ограничений на эти параметры (термостойкость, взрывобезопасность, мощность перекачивающих устройств). Ограничения могут накладываться как по технологическим, так и по экономическим соображениям.
          Итак, для решения задачи оптимизации необходимо:
               а) составить математическую модель объекта оптимизации:
     Y=f(X,U)     (5.5)
               б) выбрать критерий оптимальности и составить целевую функцию:
              
               в) установить возможные ограничения, которые должны накладываться на переменные:
                 
                 
               г) выбрать метод оптимизации, который позволит найти экстремальные значения искомых величин.
          Принято различать задачи статической оптимизации для процессов, протекающих в установившихся режимах, и задачи динамической оптимизации.
          В первом случае решаются вопросы создания и реализации оптимальной модели процесса, во втором - задачи создания и реализации системы оптимального управления процессом при неустановившихся режимах эксп луатации.

Классификация задач 

Задачи оптимизации можно классифицировать в соответствии с видом функций f, hk, gi и размерностью вектора x. Задачи без ограничений, в которых x представляет собой одномерный вектор, называются задачами с одной переменной и составляют простейший, но вместе с тем весьма важный подкласс оптимизационных задач.Задачи условной оптимизации, в которых функции hk, и gi являются линейными, носят название задач с линейными ограничениями.

В таких задачах целевые функции могут быть либо линейными, либо нелинейными. Задачи, которые содержат только линейные функции вектора непрерывных переменных x, называются задачами линейного программирования; в задачах целочисленного прграммирования компоненты вектора x должны принимать только целые значения.

Задачи с нелинейной целевой функцией и линейными ограничениями иногда называют задачами нелинейного программирования с линейными ограничениями. Оптимизационные задачи такого рода можно классифицировать на основе структурных особенностей нелинейных целевых функций. Если f(x) - квадратичная функция, то мы имеем дело с задачей квадратичного программирования; если f(x) есть отношение линейных функций, то соответствующая задача носит название задачи дробно - линейного программирования, и т.д. Деление оптимизационных задач на эти классы представляет значительный интерес, поскольку специфические особенности тех или иных задач играют важную роль при разрабоке методов их решения.

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

1.Безусловная минимизация функций:

у(х)

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

2.Минимизация функций при ограничениях в виде равенств:

у(х)                                              (7)

W : fi(х)=0 (i=1, :. , m; m<n)             (8)

Область задания Wфункции (14)-бесконечное множество-решений системы уравнений (15), представляющее собой p-мерную поверхность в n-мерном пространстве. Поскольку любая точка хW лежит на поверхности fi (х) == 0 (i = 1, ... , m) и, следовательно, является граничной, то и минимум (в случае его существования) также лежит на границе области W.

3.Минимизация функций, заданных на гиперпараллелепипеде:

у(х)

W : fi(х) 0 (i=m+1, :. , m+k)

fi(х) = 0 (i=1, :.., m).

Здесь минимум функции может находиться внутри области и на ее границах

4.Минимизация функций при ограничениях в виде неравенств. 

у(х)

W : fi(х) 0 (i=1, :. , k)

х 0

Минимум может находиться внутри области W и на ее границе.

Задачи 1, 2 обычно относят к классическим задачам оптимизации; 3, 4 называют задачами математического программирования. Легко видеть, что задачи 1-3-частный вырожденный случай задачи 4.

Класификация задач математического программирования

УПРАЖНЕНИЯ И ВОПРОСЫ К ТЕМЕ №1 

1.Что подразумевается под термином "оптимизация"?

2.Что необходимо при постановке задачи оптимизации?

3.Запишите условие задачи оптимизации общего вида. Чем оно отличается от записи в другом виде (например)?

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

R=R(X1, X2,...,XN, Y1,Y2,...,YN, U1,U2,..., UN)

5.Составьте свою задачу оптимизации. К какой классификационной группе она относится?


Тема №2

Методы безусловной оптимизации (для функции одной переменной)

МЕТОДЫ БЕЗУСЛОВНОЙ ОПТИМИЗАЦИИ

[ Одномерная оптимизация.] [ Необходимые и достаточные условия оптимальности.]
[ Методы половинного деления], [ "золотого" сечения], [ Фибоначчи.]
[ Методы с использованием производных ]. [ Mетоды полиномиальной аппроксимации ] 

Методы безусловной оптимизации делятся на методы одномерной и многомерной оптимизации.

МЕТОДЫ ОДНОМЕРНОЙ ОПТИМИЗАЦИИ

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

Своеобразным индикатором важности методов оптимизации функции одной переменной является огромное множество реализованных алгоритмов, которые условно можно сгруппировать следующим образом:

  •  методы исключения интервалов: 
    •  метод половинного деления; 
      •  метод "золотого" сечения; 
      •  метод Фибоначчи; 
    •  методы полиномиальной аппроксимации; 
    •  методы с использованием производных. 

Необходимые и достаточные условия оптимальности.

МЕТОДЫ ИСКЛЮЧЕНИЯ ИНТЕРВАЛОВ 

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

Все методы одномерной оптимизации основаны на предположении, что исследуемая целевая функция в допустимой области по крайней мере обладает свойством унимодальности, так как для унимодальной функции W(x) сравнение значений W(t) в двух различных точках интервала поиска позволяет определить, в каком из заданных двумя указанными точками подынтервалов точки оптимума отсутствуют.

Правило исключения интервалов. Пусть W(x) унимодальна на отрезке [a,b], а ее минимум достигается в точке x*. Рассмотрим x1 и x2, расположенные a<x1<x2<b.

  •  Если W(x1)>W(x2), то точка минимума W(x) не лежит в интервале (a,x1), т.е. x*Є (x1,b).
    •  Если W(x1)<W(x2), то точка минимума W(x) не лежит в интервале (x2,b), т.е. x*Є (a,x2).

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

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

Процесс применения методов поиска на основе исключения интервалов включает два этапа:

  •  этап установления границ интервала;
    •  этап уменьшения интервала.

Этап установления границ интервала

Выбирается исходная точка, а затем на основе правила исключения строится относительно широкий интервал, содержащий точку оптимума. Обычно используется эвристический метод, например, Свенна, в котором (k+1) пробная точка определяется по рекуррентной формуле

xk+1 = xk + 2kD , k=0,1,2... (3.1)

где

xo - произвольно выбранная начальная точка;

D - подбираемая величина шага.

Знак D определяется путем сравнения значений W(x), W(xo + |D | ), W(xo -|D | ):

  •  если W(xo -|D| ) W(x) W(xo + |D| ), то D имеет положительное значение;
  •  если W(xo -|D| ) W(x) W(xo + |D| ), то D имеет отрицательное значение;
  •  если W(xo -|D| ) W(x) W(xo + |D| ), то точка минимума лежит между xo - |D| и xo + |D| и поиск граничных точек завершен;
  •  если W(xo -|D| ) W(x) W(xo + |D| ), то имеем противоречие предположению об унимодальности.

Пример 3.

W(x)=(100-x)2, xo=30, |D| =5.

Определим знак D :

W(30)=4900;

W(30+5)=4225;

W(30-5)=5625.

Выполняется условие W(xo -|D| ) W(x) W(xo + |D| ), следовательно, D имеет положительное значение; x*=30.

x1=xo+20D = 35;

x2=x1+21D = 45, W(45)=3025 < W(x1) ® x*>35;

x3=x2+22D = 65, W(65)=1225 < W(x2) ® x*>45;

x4=x3+23D = 105, W(105)=25 < W(x3) ® x*>65;

x5=x4+24D = 185, W(185)=7225 > W(x4) ® x*<185.

Искомый интервал 65<x*<185.

МЕТОД ДЕЛЕНИЯ ПОПОЛАМ

Найти W(x) на отрезке [a,b].

Шаг 1. xm=(a+b)/2; L=b-a; вычислить W(xm).

Шаг 2. x1=a+L/4; x2=b-L/4; вычислить W(x1) и W(x2).

Шаг 3.

  •  Если W(x1)<W(xm), то исключить (xm,b], т.е. b=xm, xm=x1. Перейти к шагу 5.
  •  Если W(x1) W(xm), то перейти к шагу 4.

Шаг 4.

  •  Если W(x2)<W(xm), то исключить [a,xm), т.е. a=xm, xm=x2. Перейти к шагу 5.
  •  Если W(x2) W(xm), то исключить [a,x1) и (x2,b], т.е. a=x1, b=x2. Перейти к шагу 5.

Шаг 5. L=b-a. Если | L| <e , то закончить поиск. В противном случае вернуться к шагу 2.

Продолжение примера 2.1. Оптимальный раскрой лесоматериалов (MathCAD) (рис.2.1)

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

Программа:

Рисунок 2.1 - Решение примера 1.

Метод золотого сечения 

Сущность метода. Интервал поиска делится на две части так, чтобы отношение длины большого отрезка к длине всего интервала было равно отношению

=. Учитывая, что z1+z2=z, имеем

W(x)

z

z12=z z2 = (z1+z2)z2 = z1z2 + z22;

 

z1

z2

 

z1z2 + z22 - z12 = 0,

 

 

 

 

откуда =.

а

x

b

x

Рисунок 2.2. Определение коэффициента дробления

Найти W(x) на отрезке [a, b].

Шаг 1. Вычисляем коэффициент дробления отрезка [a, b]    k=(- 1)/2.

Шаг 2. x1=a+(1-k)(b-a), вычислить W(x1).

Шаг 3. x2=a+k(b-a), вычислить W(x2).

Шаг 4.

  •  Если | x2-x1| < e, где e - заданная погрешность, то xm = (x1+x2)/2, вычислить W(xm) и закончить поиск.
  •  Если | x2-x1| >e, то перейти к шагу 5.

Шаг 5.

  •  Если W(x1)>W(x2), то исключить a=x1, x1=x2 и W(x1)=W(x2). Перейти к шагу 3, затем к шагу 4.
  •  Если W(x1) <W(x2), то b=x2, x2=x1 и W(x1)=W(x2). Перейти к шагу 2 и 4.

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

Метод золотого сечения 

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

Необходимо найти минимум функции f = x4-14x3+60x2-70x. Проведем вычисления методом золотого сечения на интервале (1;2) с погрешностью 1e-5.
Программа будет иметь вид:

Примечание

   program prim;
   uses crt, ustp, metodo, idvvod;

   {*------Подключение модулей ------*
   UstP       - Установка параметров и типов, открытие файла для работы
   MetodO - Методы нахождения оптимума функции 1- переменной без ограничений
   IDVvod - Процедуры ввода исходных данных для методов оптимизации  
   }

   var
      a,b,e:real;
      nom:integer;
      fil:text;
     {$F+}
{- обязательно подключать директиву!!! }

  
{--Описываем функцию, оптимум которой необходимо найти--}
   function q(x:real):real;
       begin
       
q:=x*x*x*x-14*x*x*x+60*x*x-70*x;
       end;

  
{---Основная программа ---}
    Begin
    Clrscr;
     
{--Выводим строку для удобства пользователя, т.е. уточняем метод--}
      writeln('Нахождение оптимума по методу "золотого" сечения');
     
{--Подключаем процедуру получения имени файла для ввода результатов --}
      GetInputName(fil);
     
{--Подключаем процедуру ввода исходных данных --}
      VvodIsx_GS(a,b,e,fil);
     
{--Подключаем процедуру вычисления оптимума --}
      GoldSech(q,a,b,e,fil);
     
{--Закрываем файл результатов--}
      close(fil);
   end.

ПРИМЕР ФАЙЛА РЕЗУЛЬТАТОВ

   Результаты нахождение оптимума
   Исходные данные для метода золотого сечения
   Интервал (a,b):  a = 0.0000000000,    b = 2.0000000000
   EPSILON (точность) E=  9.99999999999612E-0006
   Нахождение оптимума по методу  "золотого" сечения
   x =  0.763932022501   a =  0.000000000000   b =  1.236067977499
   x =  0.944271909999   a =  0.472135955000   b =  1.236067977499
   x =  0.763932022501   a =  0.472135955000   b =  0.944271909999
   x =  0.832815729997   a =  0.652475842498   b =  0.944271909999
   x =  0.763932022501   a =  0.652475842498   b =  0.832815729997
   x =  0.790243257493   a =  0.721359549996   b =  0.832815729997
   x =  0.806504495004   a =  0.763932022501   b =  0.832815729997
   x =  0.790243257493   a =  0.763932022501   b =  0.806504495004
   x =  0.780193260012   a =  0.763932022501   b =  0.790243257493
   x =  0.784032017463   a =  0.773982019981   b =  0.790243257493
   x =  0.780193260012   a =  0.773982019981   b =  0.784032017463
   x =  0.781659534883   a =  0.777820777433   b =  0.784032017463
   x =  0.780193260012   a =  0.777820777433   b =  0.781659534883
   x =  0.780753327176   a =  0.779287052304   b =  0.781659534883
   x =  0.781099467719   a =  0.780193260012   b =  0.781659534883
   x =  0.780753327176   a =  0.780193260012   b =  0.781099467719
   x =  0.780885541099   a =  0.780539400555   b =  0.781099467719
   x =  0.780967253797   a =  0.780753327176   b =  0.781099467719
   x =  0.780885541099   a =  0.780753327176   b =  0.780967253797
   x =  0.780916752572   a =  0.780835039875   b =  0.780967253797
   x =  0.780885541099   a =  0.780835039875   b =  0.780916752572
   x =  0.780897462821   a =  0.780866251348   b =  0.780916752572
   x =  0.780885541099   a =  0.780866251348   b =  0.780897462821
   x =  0.780890094791   a =  0.780878173070   b =  0.780897462821
   x =  0.780885541099   a =  0.780878173070   b =  0.780890094791
   x =  0.780882726763   a =  0.780878173070   b =  0.780885541099
   Минимум найден по методу "золотого" сечения
   При Xmin = 0.780881857085 значение функции  Fmin =-24.369601567218
   Количество вычислений k =26
   

Поиск Методом Фибоначчи 

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

Предположим, что имеется интервал неопределенности (x1; x3) и известно значение функции f(x2 ) внутри этого интервала (см. рис. 2.3.). Если можно вычислить функцию всего один раз в точке x4, то где следует поместить точку x4, для того чтобы получить наименьший возможный интервал неопределенности?

Рисунок 2.3 - Графическое представление функции

Положим x2-x1 = L и x3 -x2 = R, причем L > R, как показано на рис.2.3, и эти значения будут фиксированы, если известны x1, x2 и x3. Если x4 находится в интервале (x1; x2), то:

1.Если f (x4) < f (x2 ), то новым интервалом неопределенности будет (x1;x2) длиной x2- x1=L;

2.Если f (x4) >f (x2 ), то новым интервалом неопределенности будет (x4; x3) длиной x3-x4 ;

Поскольку не известно, какая из этих ситуаций будет иметь место, выберем x4 таким образом, чтобы минимизировать наибольшую из длин x3-x4 и x2-x1. Достигнуть этого можно, сделав длины x3-x4 и x2-x1 равными, т.е. поместив x4 внутри интервала симметрично относительно точки x2, уже лежащей внутри интервала. Любое другое положение точки x4 может привести к тому, что полученный интервал будет больше L. Помещая x4 симметрично относительно x2, мы ничем не рискуем в любом случае.

Если окажется, что можно выполнить еще одно вычисление функции, то следует применить описанную процедуру к интервалу (x1;x2), в котором уже есть значение функции, вычисленное в точке x4, или к интервалу (x4;x3), в котором уже есть значение функции, вычисленное в точке x2. Следовательно, стратегия ясна с самого начала. Нужно поместить следующую точку внутри интервала неопределенности симметрично относительно уже находящейся там точке. Парадоксально. Но, чтобы понять, как следует начинать вычисления, необходимо разобраться в том, как его следует кончать.

На n-м вычислении n-ю точку следует поместить симметрично по отношению к (n-1)-й точке. Положение этой последней точки в принципе зависит от нас. Для того, чтобы получить наибольшее уменьшение интервала на данном этапе, следует разделить пополам предыдущий интервал. Тогда точка xn будет совпадать с точкой xn-1. Однако при этом мы не получаем никакой новой информации. Обычно точки xn-1 и xn отстоят друг от друга на достаточном расстоянии, чтобы определить, в какой половине, левой или правой, находится интервал неопределенности. Они помещаются на расстоянии e/2 по обе стороны от середины отрезка Ln-1; можно самим задать величину e или выбрать эту величину равной минимально возможному расстоянию между двумя точками. (Предположим, что в нашем примере инженер может регулировать температуру интервалом в 1°С, поэтому e =1.)

Интервал неопределенности будет иметь длину Ln, следовательно,

Ln-1=2Ln-e (рис. 2.4, нижняя часть).

На предыдущем этапе точки xn-1 и xn-2 должны быть помещены симметрично внутри интервала Ln-2 на расстоянии Ln-1 от концов этого интервала. Следовательно,     Ln-2=Ln-1+ Ln (рис. 2.4, средняя часть).

Рисунок 2.4

Замечание. Из рисунка ясно, что на предпоследнем этапе xn-2 остается в качестве внутренней точки.

Аналогично

Ln-3=Ln-2+ Ln-1 (рисунок 2.4, верхняя часть).

В общем случае

Lj-1=Lj+ Lj+1 при 1 < j <n

Таким образом,

Ln-1=2Ln-e ,

Ln-2=Ln-1+Ln =3Ln -e

Ln-3=Ln-2+Ln-1=5Ln-2 -e

Ln-4=Ln-3 +Ln-2 =8Ln-3 -e и т.д.

Если определить последовательность чисел Фибоначчи следующим образом:

F0=1, F1=1 и Fk =Fk-1 +F k-2 для k=2,3,...,
то L
n-j =Fj +1Ln-Fj-1 e , j=1,2,...,n-1.

Если начальный интервал (a;b) имеет длину Lj=(b-a), то

Lj=FnLn-Fn-2e,

т.е.

Ln =L1/Fn +e Fn-2/Fn .

Следовательно, произведя n вычислений функции, мы уменьшим начальный интервал неопределенности в 1/Fn раз по сравнению с его начальной длиной (пренебрегая e ), и это - наилучший результат.

Если поиск начат, то его несложно продолжить, используя описанное выше правило симметрии. Следовательно, необходимо найти положение первой точки, которая помещается на расстоянии L2 от одного из концов начального интервала, причем не важно, от какого конца, поскольку вторая точка помещается согласно правилу симметрии на расстоянии L2 от второго конца интервала:

L2=Fn-1Ln- Fn-3e = Fn-1 L1 /Fn +e (Fn-1 Fn-2 - Fn Fn-3 )/Fn = Fn-1 L1 /Fn +e (-1)n /Fn.

После того как найдено положение первой точки, числа Фибоначчи больше не нужны. Используемое значение e может определяться из практических соображений. Оно должно быть меньше L1 /Fn+1, в противном случае мы будем напрасно тратить время на вычисление функции.

Таким образом, поиск методом Фибоначчи, названный так в виду появления при поиске чисел Фибоначчи, является итерационной процедурой. В процессе поиска интервала (x1; x2 ) с точкой x2, уже лежащей в этом интервале, следующая точка x4 всегда выбирается такой,


что x
3 - x4=x2-x1 или

x4 -x1 =x3 -x2, т.е.     x4 =x1 -x2 +x3.

Если f ( x2 ) =f 2 и f (x4 ) =f 4, то можно рассмотреть четыре случая:

  •  х4 < x2, f4 < f2 Новый интервал (х1 ; x2), содержащий точку х4 
    •  х4 > x2, f4 < f2 Новый интервал (х2 ; x3 ), содержащий точку х4 
    •  х4 < x2, f4 > f2 Новый интервал (х4 ; x3 ), содержащий точку х4 
    •  х4 > x2, f4 > f2 Новый интервал (х1 ; x4 ), содержащий точку х4 

В приведенной программе реализованы указанные выше случаи. В том виде, как она приведена здесь, эта программа позволяет производить до 40 вычислений функции. В программе исследуется функция f ( x) =x4 -14x3+60x2-70x. Методы полиномиальной аппроксимации 

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

Рассмотрим следующие вопросы:

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

Квадратичная аппроксимация 

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

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

Если целевая функция W(x) в точках x1, x2, x3 принимает соответствующие значения W1, W2, W3, то можно определить коэффициенты aо, a1, a2 таким образом, что значения квадратичной функции

q(x) = ao + a1(x - x1) + a2(x - x1)(x - x2)

совпадут со значением W(x) в трех указанных точках. Вычислим q(x) в трех указанных точках (см. Табл.2.1).

Таблица 2.1 - Вычисление значений

W1 = W(x1) = q(x1) = ao

ao = W1

W2 = W(x2) = q(x2) = W1 + a1(x2 - x1)

a1 =(W2 - W1)/(x2 - x1)

W3 =q(x3) = W1 + [(W2 - W1) (x3 - x1)]/ /(x2 - x1) + a2(x3 - x1) (x3 - 2)

a2 =

=0

= -

Метод Пауэлла 

Если извесны значения функции f(x) в трех разных точках a, b, g, равные соответвенно fa, fb, fg, то функция f(x) может быть аппроксимирована квадратичной функцией

j(x)=Ax2+Bx+C,

где А, В и С определяются из уравнений

Aa2+Ba+C=fa Ab2+Bb+C=fb Ag2+Bg+C=fg 

После преобразований этих уравнений получаем

A=[(g-b)fa+(a-g)fb+(b-a)fg] / D, B=[(b2-g2)fa+(g2-a2)fb+(a2-b2)fg] / D, C=[bg(g-b)fa+ga(a-g)fb+ab(b-g)fg] / D, 

где D=(a-b)(b-g)(g-a).Ясно, что j(x) будет иметь минимум в точке x=-B/2A, если А>0. Cледовательно можно аппроксимировать точку минимума функции f(x) значением

Этод метод может непосредственно применяться к функциям одной переменной. Он может быть очень полезен для выполнения линейного поиска в процедурах, описанных в теме №3. В этих процедурах требуется найти минимум функции f(x) в точках прямой x0+ld, где x0- заданная точка, а d определяет заданное направление. Значение функции f(x0+ld) на этой прямой является значениями функции одной переменной l:

j(l) = f(x0+ld).

Идеи и результаты, изложенные выше, преобразуются в вычеслительные процедуры, описанные далее. Предположим, что заданы унимодальная функция одной переменной f(x), начальная апроксимация положения минимума и длинна шага D, является величиной того же порядка, что и расстояние от точки А до точки истенного минимума x*(условие, которое не всегда просто удовлетворить). Вычислительная процедура имеет следующие шаги:

Шаг 1. x2 = x1 + D x.

Шаг 2. Вычислить W(x1) и W(x2).

Шаг 3.

  •  Если W(x1) > W(x2), то x3 = x1 + 2 D x.
    •  Если W(x1)< W(x2), то x3 = x1 - D x.

W(x1) > W(x2),

Шаг 4. Вычислить W(x3) и найти

Wmin = min{ W(x1),W(x2), W(x3)},

Xmin = xi, соответствующая Wmin.

Шаг 5. По x1, x2, x3 вычислить x*, используя формулу для оценивания с помощью квадратической аппроксимации.

Шаг 6. Проверка окончания

  •  Если | Wmin - W(x*)| < W, то закончить поиск. В противном случае к шагу 7.
    •  Если | Xmin - x*| < x, то закончить поиск. В противном случае к шагу 7.

Шаг 7. Выбрать Xmin или x* и две точки по обе стороны от нее. Обозначить в естественном порядке и перейти к шагу 4.

Данный метод часто называют методом Пауэлла, он реализован в программе

Заметим, что если точка Е задана слишком малой, то a, b, g, а также fa, fb, fg будут очень близко друг к другу и значение d (см. Уравнение (1) ) может стать вообще недостижимыми. Чтобы преодолеть эту трудность, перепишем уравнение(1) для второй и последней интерполяции:

Кубическая интерполяция

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

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

Рассмотрим задачу минимизации функции f(x) на прямой x0 + hd , то есть минимизацию функции

Предположим, что известные следующие значения:

Эту информацию можно использовать для построения в кубического полинома

который будет аппроксимировать функцию j(h) Если p=0 , то уравнения, которые определяют a, b, c, d имеют вид :

Точки перегиба кубического полинома являются решением уравнения

Следовательно, если r является точкой минимума кубического полинома,

Одно из значений этого выражения является минимумом. Друга производная равна 2c +6dh 

Если мы выберем положительный знак, то при

друга производная будет


(2.1)

Выбор точки q зависит от нас. Если Gp >0 , то надо выбрать значение q положительным, то есть сделать шаг в направлении уменьшения функции j (h) , в другом случае значения q надо выбирать отрицательным. Значение должно быть таким, чтобы интервал (0, ) включал в себя минимум. Это будет верным, если j q > j p или Gp >0.

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

Давидон, Флетчер и Пауэлл предложили выбирать q следующим чином:

где j m - оценка самого малого значения истинного минимумf j (h) , а h - константа, значение которой обычно берут 2 или 1.

итерационная процедура:
1. Найти  и 
2. Проверить, если не выполняется условие
Gp < 0 , то выполнять поиск в направлении -d . Выбрать q из выражения. При этом надо "угадать" j m .
3. Вычислить  и 
4. Если
Gq < 0 или j q > j p , то интервал, который содержит минимум, найден. В другом случае изменить q на 2q и возвратиться к шагу 3.
5. Использовать уравнение (2.1) для апроксимации точки минимума на интервале
(0, q) значением r.
6. Если  , где
e - заданная точность, то остановиться.
7. Возвратиться на шаг 5, используя интервал
(0,r) , если Gr > 0 , или используется интервал (r,q) , если Gr0.

Квадратичные функции 

Квадратичная функция

(1)

где a - стала; b - постоянный вектор и G - положительно определенная симметричная матрица - имеет минимум в точке x* , где x* определяется следующим путем:

(2)

оттуда x*= -G-1 b. Любую функцию можно аппроксимировать в окресности точки x0 функцией

(3)

где G(x0) - матрица Гессе, вычисленная в точке x0.

Аппроксимацией минимума функции f(x) может быть минимум функции j(x). Если последний находится в точке xm, то

оттуда

или

(4)

Таким образом точкой xи следующей аппроксимации минимума будет:

(5)

или,

(6)

где длина шага lи определяется одномерным поиском в направлении G-1(xи)g(xи). 

Метод Ньютона-Рафсона основан на последнем уравнении.

Уравнения (5) и (6) в том виде, как они записаны, требуют вычисления и транспонирования матрицы Гессе на каждом шаге, который чаще является главной частью вычислений. Если точка xи близко расположенная к точке x*, то сходимость будет быстрой, потому, что в общем случае функция j(x) будет хорошо аппроксимировать функцию f(x) в этой ячейке. Как норму градиента |g(xи+1)|, так и расстояние между точками | xи+1 - xи | надо проверить на выполнение критерий завершения. Интересно отметить, что сравнивая с простым методом самого быстрого спуска направлением спуска в данном случаю будет не -g(xи), а -G-1(xи)g(xи), где учитываются и вторые сменные. Методом Давидона - Флетчера - Пауэлла можно получить самый лучший результат, делая поиск на и-му этапе в направлении -Hиg(xи), где Hи- положительно определенная симметричная матрица, которая в конечном счете будет приравнивать -G-1(x*). Таким образом, этот метод обходит как вычисления, так и транспонирование матрицы G(xи) при каждой итерации.

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

Говорят, что два направления p и q сопряженные относительно симметричной положительно определенной матрицы G, если

(7)

Можно показать, что если p0, p1,:,pn-1 это n взаимно сопряженных направлений в n-мерном пространстве, то они ленейно независимые.

Перепишем функцию (1) в виде

причем ее минимум находится в точке

и

(8)

Допустим, что для поиска минимума функции (8) используется итерационная процедура.

Начнем из точки x0 и проведем поиск в направлении p0 с целью нахождения минимума в точке

(9)

где l0 - некоторая скалярная величина.

Отметим, что в точке x1 направление   ортогональный направлению p0 и

(10)

В общему случаю на шаге и происходит поиск из точки xи в направлении pи с целью нахождения минимума в точке

(11)

где для F(x) справедливые соотношения

(12) (13)

Повторным применением уравнения (11) после шагов получим

(14)

для всех j на промежутке 0j<n-1.

Таким образом, согласно уравнению (13)

(15)

Следовательно,

(16)

и из уравнения (16) получим:

Сейчас, если все векторы p0,p1,:,pn-1 взаимно сопряженные так, что   то из соотношения (12) получим

Но потому, что в этом случае векторы p0,p1,:,pn-1 ленейно независимые и, таким образом, образовывают базис, то

g(xn)=0 

оттуда

G(xn - x*)=0 

и

xn = x*.

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

Методы с использованием производных

Целесообразно предположить, что эффективность поисковых процедур существенно повысится, если в дополнение к условию непрерывности ввести требование дифференцируемости целевой функции. Необходимое условие существования оптимума целевой функции в точке x*

W'(x*) = dW/dx f x=x* = 0.

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

  •  метод хорд; 
    •  метод касательных; 
    •  метод средней точки. 

Метод хорд 

Сущность метода. Ориентирован на нахождение корня уравнения W'(x) в интервале [a,b], в котором имеются две точки N и P, в которых знаки производных различны. Алгоритм метода хорд позволяет аппроксимировать функцию W'(x) "хордой" и найти точку, в которой секущая графика W'(x) пересекает ось абсцисс.

Рисунок 2.5 - Схема метода хорд

Шаг 1. Следующее приближение к стационарной точке x* определяется по

формуле

R = P - .

Шаг 2. Вычислить W'(R).

Шаг 3. Если | W'(R)| < e , то закончить поиск. В противном случае необходимо выбрать одну из точек P или N, чтобы знаки производных в этой точке и точке R были различны. Вернуться к шагу 1.

Как видно из алгоритма, метод хорд реализован на исследовании как знака производной, так и ее значении. Поэтому он более эффективен, чем метод средней точки.

Пример.

Минимизировать W(x)=2x2+(16/x) в интервале 1 < x < 5.

W'(x) = dW(x)/dx = 4x-16/x2.

Итерация 1.

Шаг 1. N=1, P=5, W'(5)=19.36, W'(1) = -12.

Шаг 2. R=5-{19.36/[(19.36+12)/4]}=2.53.

Шаг 3. W'(2.53)=7.62>0; положить R=2.53.

Итерация 2.

Шаг 2. R=2.53-{7.62/[(7.62+12)/1.53]}=1.94.

Шаг 3. W'(1.94)=3.51>0; положить R=1.94.

Итерации продолжаются до тех пор, пока не будет выполняться неравенство | W'(R)| < e.

Метод касательных 

Сущность метода. Ориентирован на нахождение корня уравнения W'(x) в интервале [a,b], в котором имеются две точки N и P, в которых знаки производных различны. Работа алгоритма начинается из точки xo, которая представляет начальное приближение корня уравнения W'(x)=0. Далее строится линейная аппроксимация функции W'(x) в точке x1, и точка, в которой аппроксимирующая линейная функция обращается в нуль, принимается в качестве следующего приближения. Если точка xk принята в качестве текущего приближения к оптимальной точке, то линейная функция, аппроксимирующая функцию W'(x) в точке xk, записывается в виде

W'(x,xk) = W'(xk) + W''(xk)(x-xk).

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

Рис.2. Схема метода касательных

Шаг 1. Следующее приближение к стационарной точке x* определяется по формуле

xk+1 = xk - [W'(xk)/W''(xk)].

Шаг 2. Вычислить W'(xk+1), W''(xk+1)

Шаг 3. Если |W'(xk+1)| < e , то закончить поиск. В противном случае необходимо вернуться к шагу 1.

Как явствует из алгоритма, целевая функция W(x) должна быть дважды дифференцируема.

Пример 3.

Минимизировать W(x)=2x2+(16/x), положив x1=1.

W'(x) = dW(x)/dx = 4x-16/x2,

W"(x)= 4+(32/x3).

Итерация 1.

x1=1, W'(1)=-12, W"(1)=36, x2=1-(-12/36)=1.33

Итерация 2.

x2=1.33, W'(1.33)=-3.73, W"(1.33)=17.6, x3=1.33-(-3.76/17.6)=1.54

Итерации продолжаются до тех пор, пока не будет выполняться неравенство |W'(xk)| <e.

Метод средней точки 

Сущность метода. Основан на алгоритме исключения интервалов, на каждой итерации которого рассматривается одна пробная точка R. Если в точке R выполняется неравенство W'(R) < 0, то в следствии унимодальности функции точка оптимума не может лежать левее точки R. Аналогично, если W'(R) > 0, то интервал x>R можно исключить.

Пусть в интервале [a,b] имеются две точки N и P, в которых производные W'(N)<0 и W'(P)>0. Оптимальная точка x* расположена между N и P.

Шаг 1. Положить P=b, N=a, причем W'(a)<0 и W'(b)>0.

Шаг 2. Вычислить R=(P+N)/2 и W'(R).

Шаг 3.

  •  Если | W'(R)| < e , то закончить поиск. В противном случае, если W'(R)<0, положить N=R, и перейти к шагу 2.
  •  Если | W'(R)| > e , положить P=R и перейти к шагу 2.

Как следует из логической структуры, процедура поиска по методу средней точки основана на исследовании только знака производной.

Пример.

Минимизировать W(x)=2x2+(16/x) в интервале 1 < x < 5.

W'(x) = dW(x)/dx = 4x-16/x2.

Итерация 1.

Шаг 1. N=1, P=5, W'(5)=19.36, W'(1) = -12.

Шаг 2. R=(5+1)/2=3.

Шаг 3. W'(3)=10.22>0; положить P=3.

Итерация 2.

Шаг 2. R=(3+1)/2=2.

Шаг 3. W'(2)=4>0; положить P=2.

Итерации продолжаются до тех пор, пока не будет выполняться неравенство | W'(R)| <e.

Функции одной переменной. Необходимые и достаточные условия.

Функция f(x) имеет локальный минимум в точке x0, если существует некоторая положительная величина d, такая, что если |x - x0| < d, то f(x)>=f(x0), если существует окрестность точки x0, такая, что для всех значений x в этой окрестности f(x) > f(x0). Функция f(x) имеет глобальный минимум в точке x*, если для всех x справедливо неравенство f(x)>=f(x*).

На рис.2.8 дано графическое представление функции f(x), которая имеет локальный минимум в точке x0 и глобальный минимум в точке x*.

Рисунок 2.8 - Графическое представление функции

Классический подход к задаче нахождения значений x0 и x* состоит в поиске уравнений, которым они должны удовлетворять. Представленная на рис.2.8 функция и ее производные непрерывны, и видно, что в точках x0 и x* будут решениями уравнения

f'(x)=0     (2.1)

Точка xm, в которой достигается локальный минимум, и точка xc, в которой имеется точка горизонтального перегиба функции, также удовлетворяет этому уравнению. Следовательно, уравнение (2.1) является только необходимым условием, но не является достаточным условием.

Заметим, однако, что в точках x0 и x* производная f'(x) меняет знак с отрицательного на положительный. В точке xm знак меняется с положительного на отрицательный, в то время как в точке xc он не меняется. Следовательно, производная в минимуме является возрастающей функцией, а поскольку степень возрастания f'(x) измеряется второй производной, можно ожидать, что f''(x0)>0,
f'' (x*)>0, тогда как f'' (x) <0.

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

Полученные выше результаты могут найти надежное обоснование, если рассмотреть разложение функции f(x) в ряд Тейлора в окрестности точки x0 (или x*, или xm), что, конечно, требует непрерывности функции f(x) и ее производных:

f(x0 + h) - f(x0) = hf'(x)+(h2/2!) f''(x0)+...   (2.2)

Если в точке x0 достигается минимум, то левая часть (2.2) будет неотрицательной для любого достаточно малого h (|h|<d). Следовательно, первая производная f'(x0) должна быть равна нулю, и это и есть достаточное условие (см. уравнение 2.1). Если бы она была положительной, то достаточно малое отрицательное значение h делало бы правую часть (2.2) отрицательной, а если бы она была отрицательной, то достаточно малое положительное значение h делало бы ее отрицательной.

Так как в следующем члене (2.2) всегда h2 > 0, то, если

f'' (x*)>0,       (2.3)

в точке x0 достигается минимум. Если f'(xm) = 0 и f''(xm) < 0, то из аналогичных соображений в точке xm достигается максимум. Для определения различия между локальным и глобальным минимумами необходимо сравнить значения функций f(x0) и f(x*).

УПРАЖНЕНИЯ И ВОПРОСЫ К ТЕМЕ №2

1. Найти значения максимума и минимума функции f(x) = x(x-1)2.

2. Найти значения максимума и минимума функции f(x) = x/(x2+1).

3.Показать, что минимальным значением функции acosq + bsinq является -

4. Равнобедренный треугольник с вертикальным углом 2q вписан в окружность радиуса r. Найти выражение для площади треугольника как функции от q, показать, что она максимальна, когда треугольник равносторонний.

5. Исследовать функцию f(x) = x2/3 - 1. Нарисовать ее график. Показать, что f(x) имеет минимум при x = 0. Чему равно значение f'(x) при x = 0? Меняет ли знак f'(x), если x возрастает при прохождении через 0?

6. Исследовать функцию f(x) = |x|. Найти ее минимум. Что можно сказать относительно поведения f'(x) в точке минимума?

7. Найти минимум функции f(x) = -e-xsh(x/2).

8. В процессе производства определенного количества некоторого товара его цена устанавливается равной LK. Товар хранится на складе до тех пор, пока не будет использован, и стоимость хранения единицы товара составляет LS в единицу времени. Показать, что если товар производится регулярно в количестве x в течение времени x/R, то стоимость функционирования такой системы в единицу времени равна

С = KR/x + Sx/2

Покажите, что стоимость С достигает минимума при х =

9. Исследовать точки перегиба функции f(x) = x4 - 14x3 + 60x2 - 70x. [Этот кажущийся простым пример иллюстрирует одну из классических задач. Требуется решить уравнение f'(x) = 0. В данном случае оно является кубическим уравнением, которое не так просто раскладывается на множители. Необычным является способ решения, при котором используется один из численных методов, описанных в следующей главе: он предназначен для минимизации функции j (x) = [f'(x)]2. Минимум функции j(x) равен нулю, а это означает, что получено решение уравнения f'(x) =0]

10. Если F0 = 1 и F1 = 1, а Fn = Fn-1 + Fn-2 для n>=2, то показать, что

Fn = (tn+1 - (- t)-(n+1))/,

где t = (1 + )/2 = 1,618033989.

Используя это рекуррентное соотношение, покажите, что F2= 2; F3=3, F4=5, F5=8, F6=13, F10=89, F19=6765, F20=10946.

11. Показать, что

1) Fn-1Fn-2 - Fn Fn - 3 = (-1) n.

2) Fn = t;n+1/ для больших n.

3) Fn-1/Fn = 1/t для больших n.

12. Если расстояние между точками равно по крайней мере Eps, то 2Eps будет наименьшим интервалом неопределенности; следовательно, неравенство Ln >= 2Eps устанавливает границу для n - числа полезных экспериментов. Показать, что неравенство приводит к соотношению Eps < L1 /Fn+1.

13. Показать, что для сокращения интервала неопределенности на 1% от начальной величины необходимо выполнить 11 вычислений функции при использовании метода Фибоначчи. Если размещение точек задать вначале, то каково будет минимальное число необходимых точек? (Не рассчитывайте на получение удачных результатов).

14. Применить метод Фибоначчи при наличии 10 вычислений функции для определения минимума функции 2x2 + 3e-x на интервале (0; 1).

15. Воспользоваться методом Фибоначчи для нахождения минимума функции x4 - 14x3 + 60x2 - 70x. на интервале (5; 7) с точностью 0,01. Сколько раз необходимо вычислить функцию?

16. Использовать метод "золотого сечения" для определения минимума функции 2x2 + 3e-x с точностью до двух десятичных знаков. В качестве начального интервала неопределенности использовать интервал (0; 1). Сколько раз необходимо вычислить функцию? Сравните с упр.14.

17. Применить квадратичную интерполяцию для определения точки минимума функции -e-x ln(x) на интервале (1, 3) с точностью до 0,001. Проверьте, не используете ли вы отрицательных значений x, ибо в противном случае получите ошибку при работе программы.

18. Применить квадратичную интерполяцию для поиска минимума функции f(x1, x2) = x1 2 + 2x2 2 + 2x1 x2 на прямой a + ld, где a =1 и d=2/3 c точностью до двух десятичных знаков.

19. Решить уравнение exsinx = 1. Попытайтесь минимизировать функцию f(x) = (1 - exsin x)2.

Тема №3
Методы безусловной оптимизации (для функции n переменных)

МЕТОДЫ БЕЗУСЛОВНОЙ ОПТИМИЗАЦИИ

[ Многомерная оптимизация.] [ Необходимые и достаточные условия оптимальности.]
[ Методы нулевого порядка.] [ Методы первого порядка.] 

Методы безусловной оптимизации делятся на методы одномерной и многомерной оптимизации.

К методам многомерной оптимизации относятся:

  1.  Методы нулевого порядка:
    •  покоординатного спуска, 
    •  Хука-Дживса, 
    •  симплексный метод Нелдера-Мида. 
  2.  Методы первого порядка:
  •  градиентный; 
  •  наискорейшего спуска; 
  •  сопряженных градиентов:
    •  метод Давидона-Флетчера-Пауэлла; 
    •  метод Флетчера-Ривса. 

Необходимые и достаточные условия оптимальности.

Метод покоординатного спуска

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

Рисунок 3.1 - Метод прямого поиска

Рассмотрим функцию двух переменных. Ее линии уровня представлены на рис.3.1, а минимум лежит в точке (x1*,x2*). Простейшим методом поиска является метод покоординатного спуска. Из точки А произведем поиск минимума вдоль направления оси х1 и, таким образом, находим точку В, в которой касательная к линии постоянного уровня параллельна оси х1. Затем, производя поиск из точки В в направлении оси х2, получаем точку С, производя поиск параллельно оси х2, получаем точку D, и т.д. Таким образом, мы приходим к оптимальной точке. Очевидным образом эту идею можно применить для функции n переменных.

Данный метод проиллюстрирован программой.

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

Описание метода Хука-Дживса

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

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

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

Xp(k+1)=X(k)+(X(k)-X(k-1))

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

X(k) - текущая базовая точка,

X(k-1)- предыдущая базовая точка,

Xp(k+1)- точка, построенная при движении по образцу,

X(k+1)- следующая (новая) базовая точка.

Приведенная последовательность характеризует логическую структуру поиска по методу Хука-Дживса.

Описание алгоритма метода 

А. Выбрать начальную базисную точку b1 и шаг длиной hj для каждой переменной xj, j=1,2,...,n.

Б. Вычислить f(x) в базисной точке b1 с целью получения сведений о локальном поведении функции f(x). Функция f(x) в базисной точке b1 находится следующим образом:

Вычисляется значение функции f(b) в базисной точке b1.

Каждая переменная по очереди изменяется прибавлением длины шага. Таким образом, мы вычисляем значение функции f(b1+h1*e1), где e1-единичный вектор в направлении оси х1.

  •  Если это приводит к уменьшению значения функции, то d1 заменяется на b1+h1*e1. В противном случае вычисляется значение функции f(b1-h1*e1), и если ее значение уменьшилось, то b1 заменяем на b1-h1*e1.
  •  Если ни один из проделанных шагов не приводит к уменьшению значения функции, то точка b1 остается неизменной и рассматриваются изменения в направлении оси х2, т.е. находится значение функции f(b1+ h2*e2) и т.д.   Когда будут рассмотрены все n-переменныx,   мы будем иметь новую базисную точку b2.
  •  Если b2=b1, т.е. уменьшение функции не было достигнуто, то исследование повторяется вокруг той же базисной точки b1, но с уменьшенной длиной шага.
  •  Если b2 b1, то производится поиск по образцу.

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

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

P1=b1+2*(b2-b1).

В общем случае

Pi=bi+2*(b(i+1)-bi).

2.Затем исследование следует продолжать вокруг точки P1(Pi).

Если наименьшее значение на шаге B,2 меньше значения в базисной точке b2(в общем случае b(i+1)), то получают новую базисную точку b3 (b(i+2)), после чего следует повторить шаг B,1 . В противном случае не производить поиск по образцу из точки b2(b(i+1)), а продолжить исследования в точке b2(b(i+1)).

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

Данный метод реализован в программе

Метод Нелдера-Мида 

Метод Нелдера-Мида (поиск по деформируемому многограннику) является развитием симплексного метода Спендли, Хекста и Химсворта. Множество (n+1)-й равноудаленной точки в n-мерном пространстве называется регулярным симплексом. Эта конфигурация рассматривается в методе Спендли, Хекста и Химсворта. Следовательно, в двумерном пространстве симплексом является равносторонний треугольник, а в трехмерном пространстве правильный тетраэдр. Идея метода состоит в сравнении значений функции в (n+1)вершинах симплекса и перемещении в направлении оптимальной точки с помощью итерационной процедуры. В симплексном методе, предложенном первоначально, регулярный симплекс использовался на каждом этапе. Нелдер и Мид предложили несколько модификаций этого метода, допускающих, чтобы симплексы были неправильными. В результате получился очень надежный метод прямого поиска, являющийся одним из самых эффективных, если n<=6. 

В методе Спендли, Хекста и Химсворта симплекс перемещается с помощью трех основных операций: отражения, растяжения и сжатия. Смысл этих операций станет понятным при рассмотрении шагов процедуры.

А. Найдем значения функции

f1=f( x1),     f2=f(x2) ... fn+1=f(xn+1)

в вершинах симплекса.

Б. Найдем наибольшее значение функции fh, следующее за наибольшим значением функции fg , наименьшее значение функции fl и соответствующие им точки xh, xg и xl.

В. Найдем центр тяжести всех точек, за исключением точки xh. Пусть центром тяжести будет

x0 = xi (3.1)

и вычислим f(x0)=f0.

Г. Удобнее всего начать перемещение от точки xh. Отразив точку xh относительно точки x0, получим точку xr и найдем f(xr) = fr.

Операция отражения иллюстрируется рис.3.2. Если a>0 - коэффициент отражения, то положение точки xr определяется следующим образом:

xr-x0=a (x0-xh),

т.е.

xr=(1+a)x0 -axh.     (3.2)

Замечание. a= |xr-x0| / |x0 -xh|.

Д. Сравним значения функций fr и fl.

1.Если fr<fl, то мы получили наименьшее значение функции. Направление из точки x0 в точку xr наиболее удобно для перемещения. Таким образом, мы производим растяжение в этом направлении и находим точку xe и значение функции fe=f(xe). Рисунок 3.3. иллюстрирует операцию растяжения симплекса. Коэффициент растяжения g >1 можно найти из следующих соотношений:

xe-x0=g (xr-x0),

т.е.

xe=gxr+ (1-g)x0. (3.3)

Рисунок 3.2 - Операция отражения

Рисунок 3.3 - Операция растяжения

Замечание. g = | xe-x0 | / |xr-x0 |.

а) Если fe>=fl, то заменяем точку xh на точку xe и проверяем (n+1)-ую точку симплекса на сходимость к минимуму ( см. шаг З). Если сходимость достигнута, то процесс останавливается; в противном случае возвращаемся на шаг Б.

б) Если fe=fl , то отбрасываем точку xe. Очевидно, мы переместились слишком далеко от точки x0 к точке xr. Поэтому следует заменить точку xh на точку xr, в которой было получено улучшение (шаг Д, 1) проверить сходимость и, если она достигнута, вернуться на шаг В.

2.Если fr>fl, но fr <=fgто xr является лучшей точкой по сравнению с другими двумя точками симплекса и мы заменяем точку xh на точку xr и, если сходимость не достигнута, возвращаемся на шаг Б, т.е. выполняем пункт 1,б, описанный выше.

3.Если fr>fl и fr>fgто перейдем на шаг Е.

Е. Сравним значения функций fr и fh.

1. Если fr<fh, то заменяем точку xh на точку xr и значение функции fh на значение функции fr. Запоминаем значение fr>fg из шага Д.2, приведенного выше. Затем переходим на шаг Е.2

2. Если fr>fh, ясно, что мы переместились слишком далеко от точки xh к точке x0. Пытаемся исправить это, найдя точку xc (а затем и f c) с помощью шага сжатия, показанного на рис. 3.4.

Если fr>f h, то сразу переходим к шагу сжатия и находим точку xc из соотношения

xc-x0=b(x0-xh),

где b(0<b<1)- коэффициент сжатия. Тогда

xc=(1+b)x0-bxh.    (3.4)

Если fr<f h, то сначала заменим точку xh на точку xr, а затем произведем сжатие. Тогда точку xc найдем из соотношения

xc-x0=b(xr-x0),

т.е.

xc=bxr+(1-b)x0. (3.5)

(рис. 3.5).

Ж. Сравниваем значения функций fc и fh.

1.Если fc<f h, то заменяем точку xh на точку xc, и если сходимость не достигнута ,то возвращаемся на шаг Б.

2.Если fc>f h, то очевидно, что все наши попытки найти значение меньшее fh закончились неудачей, поэтому мы переходим на шаг З.

З. На этом шаге мы уменьшаем размерность симплекса делением пополам расстояния от каждой точки симплекса до xl-точки, определяющей наименьшее значение функции.

Рисунок 3.4. - Шаг сжатия для fr>fh

Рисунок 3.5 - Шаг сжатия для fr<fh

Таким образом, точка xi заменяется на точку xi-( xi-xl), т.е. заменяем точку xi точкой

(xi-xl) (3.6)

Затем вычисляем fi для i=1,2,...,(n+1), проверяем сходимость и, если она достигнута, возвращаемся на шаг В.

И. Проверка сходимости основана на том, чтобы стандартное отклонение (n+1) -го значения функции было меньше некоторого заданного малого значения ( В этом случае вычисляется

s2=(fi-`f)2/(n+1) (3.7)

где `f= fi / n+1.

Если s< e, то все значения функции очень близки друг к другу, и поэтому они, возможно, лежат вблизи точки минимума функции xl. Исходя из этого, такой критерий сходимости является разумным, хотя Бокс, Дэвис и Свенн [1] предлагают то, что они считают более "безопасной" проверкой.

Коэффициенты a, b, g в вышеприведенной процедуре являются соответственно коэффициентами отражения, сжатия и растяжения. Нелдер и Мид рекомендуют брать a=1, b=0.5 и g=2. Рекомендация основана на результатах экспериментов с различными комбинациями значений. Эти значения параметров позволяют методу быть эффективным, но работать в различных сложных ситуациях.

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

x2=x1+ke1

x3=x1+ke2 (3.8)

xn+1=x1+ken

где k - произвольная длина шага,
   а e
j - единичный вектор.


Данный метод реализован в
программе 

Градиентный метод 

В этом методе наряду со значением функции используется и ее градиент.

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

Это свойство может быть обосновано следующим образом. Предположим, что осуществляется перемещение из точки x в следующую точку x+h*d, где d - некоторое направление, а h - шаг некоторой длины. Следовательно, перемещение производится из точки (x1,x2,..,xn) в точку (x1+dx1, x2+dx2, ..., xn+dxn) , где

dxi =hdi,           (3.9)

а di - косинусы направления , такие, что

di2 = 1.           (3.10)

Изменение значений функции определяется соотношениями

df =f(x1+dx1, x2+dx2, ..., xn+dxn) - f(x1,x2,..,xn) = dx1+dx2+...+dxn (3.11)

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

f(x0+h) - f(x0) = hi(x1, x2, ... , xn) + hihj(x1,x2,..,xn)+...

=hT+1/2*hTG(x0)h+...    (3.11*)

Как следует выбрать направления di, удовлетворяющие уравнению (3.11), чтобы получить наибольшее значение изменения функции df?

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

j(d1,d2,..,dn)= df+l(di2 -1).

Величина df, удовлетворяющая ограничению (3.9), достигает максимума, когда функция

j(d1,d2,..,dn)=h(d1+d2 +...+dn) + l( d12 +d22 +...+dn2 - 1)

достигает максимума. Ее производная

=h+2ldj, при j=1, 2,...,n.      (3.12)

Если   =0 , то dj= - .      (3.13)

Следовательно,

       (3.14)

Тогда di ~ и направление d параллельно направлению в точке х.

Таким образом, наибольшее локальное возрастание функции для заданного малого шага h имеет место, когда d есть направление или g(x). Поэтому направлением наискорейшего спуска является направление

- или - g(x).         (3.15)

В более простом виде уравнение (3.11) можно записать так:

df= |||dx| cosq ,

где q- угол между векторами и dx. Для заданной величины мы минимизируем df, выбирая q=180О, чтобы направление dx совпадало с направлением -.

Замечание. Направление градиента перпендикулярно в любой точке линии постоянного уровня, поскольку вдоль этой линии функция постоянна. Таким образом, если (d1,d2,..,dn) - малый шаг вдоль линии уровня, то

f(x1+d1, x2+d2, ..., xn+dn) = f(x1,x2,..,xn)

и, следовательно,

df= dj= []Td=0.       (3.16)

(рис.3.6)

Рисунок 3.6 - Метод градиентного спуска

Данный метод представлен программой и может быть применен для решения прикладных задач

Описание метода наискорейшего спуска 

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

xi+1=xi - li,(3.17)

где li - значение l, минимизирующее функцию

j (li)=f[xi- l].(3.18)

Значение li может быть найдено с помощью одного из методов одномерного поиска (например, методом квадратичной интерполяции).

В приложении приведена программа, позволяющая реализовать метод наискорейшего спуска. В ней множитель Лагранжа обозначен через h. Вектор di является единичным.

Для поиска минимума функции

j( l)=f(xi- ldi) (3.19)

в направлении di из точки xi используется метод квадратичной интерполяции.

В точке xil= 0, и мы выбираем длину шага l такой, чтобы шаг "перекрыл " минимум функции j(l) . Производная

= g(xi-ldi)Tdi(3.20)

Данный оператор for(i=0;i<n;i++) g2+=g[i]*d[i]; - вычисляет выражение

f(x0-hdi)Td=g(x0-hdi)Td (3.21)

Оператор if (ff[2]>=ff[0] || g2>=0) проверяет условие "перекрытия" минимума, которое выполняется при выполнении либо одного, либо другого условия. Если минимум не попал в отрезок (0,l), то l удваивается, и это повторяется столько раз, сколько необходимо для выполнения условия "перекрытия".

Удостоверившись, что отрезок (0,l) содержит минимум, в качестве третьей точки возьмем точку l/ 2. Минимальную точку сглаживающего квадратичного полинома находим в соответствии с соотношением

, (3.22)

что отражено следующими операторами

l[3]=h*(ff[1]-.75*ff[0]-.25*ff[2]);

l[3]/=2*ff[1]-ff[0]-ff[2];

Оператор for(i=0;i<n;i++)

{ x[i]=y[i]+l[0]*d[i]; y[i]=x[i]; }

производит присваивание xi+1=xi, и если |g(xi+1)| достаточно мало, то процесс заканчивается. В процессе поиска предполагается сходимость к экстремуму, поэтому для эффективности процедуры разумно уменьшить длину шага. При этом деление шага пополам выбрано произвольно.

Данный метод реализован в программе и может быть применен для решения прикладных задач.

Метод Давидона-Флетчера-Пауэлла 

Метод Давидона-Флетчера-Пауэлла (ДФП) основан на использовании соотношений (3.23) и (3.24).

Из уравнения

f(x0+h)-f(x0) =

    =hT*Ñ f(x0)+* hTG(x0)* h+...

следует, что при выполнении условий непрерывности любую функцию можно аппроксимировать в окрестности точки x0 функцией

j(x)=f(x)+(x-x0)T* Ñ f(x0)+* (x-x0)T*G(x0)* (x-x0),   (3.23)

где G(x0)-матрица Гессе, вычисленная в точке x0.

Разумной аппроксимацией минимума функции f(x) может быть минимум функции j(x). Если последний находится в точке xm, то

Ñ f(x0)+G(x0)(xm-x0)=0,

откуда

xm=x0-G-1(x0)* Ñ f(x0),

или

xm=x0-G-1* g(x0)    (3.24)

В методе Давидона-Флетчера-Пауэлла не требуется на каждом шаге вычислять обратный гессиан G-1(xi) так как направление поиска на шаге i является направлением -Hig(xi), где Hi - положительно определенная симметричная матрица, которая обновляется на каждом шаге, как это будет описано ниже. В пределе матрица H становится равной обратному гессиану.

Начнем поиск из начальной точки x0, взяв в качестве начальной матрицу H0 (обычно единичную матрицу, хотя в этом случае может подойти любая симметричная положительно определенная матрица). Итерационная процедура может быть представлена следующим образом (вместо g(xi) удобнее писать gi):

  1.  На шаге i имеется точка xi и положительно определенная симметрическая матрица Hi.
  2.  В качестве направления поиска взять направление di=-Hi* gi 
  3.  Чтобы найти функцию l i, минимизирующую функцию f(xi+l idi), произвести одномерный поиск вдоль прямой xi+l* di.
  4.  Положить

vi=li* di.

5. Положить

xi+1=xi+vi.

6.Найти f(xi+1) и gi+1.Завершить процедуру, если величины |gi+1| или |vi| достаточно малы. В противном случае продолжить.

Замечание. Из соотношения ( g(xi+1)Tpi=0, где pi-направление с целью нахождения минимума в точке) следует, что

gi+1T* vi=0.

7.Положить

ui=gi+1 - gi.

8.Обновить матрицу H следующим образом:

Hi+1=Hi+Ai+Bi,

где

Ai=vi* viT/(viT* ui);

Bi=-Hi*ui uiT*Hi/ (uiT*Hi ui).

9.Увеличить i на единицу и вернуться на шаг 2.

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

Данный метод представлен программой и может быть применен для решения прикладных задач

Метод Флетчера-Ривса 

Метод Флетчера-Ривса основан на том, что для квадратичной функции n переменных n одномерных поисков вдоль взаимно сопряженных направлений позволяют найти минимум.

Рассмотрим функцию

f ( x) = a+ bTx+ xTGx.

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

В качестве первого направления поиска из первой точки x1 возьмем направление наискорейшего спуска

d1= - g1   (1)

и найдем значение l 1, минимизирующее функцию

f ( x1+ l d1).

Положим

x2= x1+ l 1d1    (2)

и произведем поиск в направлении d2, сопряженном направлению d1 (выберем вектор d2 как линейную комбинацию векторов d1 и - g2), и найдем

x3= x2+ l 2d2     (3)

минимизацией функции f ( x2+ l d2). Направление поиска d2 из точки x3 выбирается сопряженным направлениям d1 и d2. На (k + 1 ) - м шаге выбираем dk + 1 в виде линейной комбинации - gk + 1 , d1, d2,...,dk , сопряженной всем направлениям d1, d2,...,dk .

Таким образом, dk + 1 = - gk + 1 + kr=1 ar dr, k = 1, 2 . . . Оказывается, все a r равны нулю, за исключением a k, так что

dk+1=-gk+1+a kdk    (4)

и

a k=g2k+1/g2k.    (5)

Прежде чем перейти к индуктивным рассуждениям, докажем справедливость соотношений (4) и (5) при k=1.Поскольку f ( x2) = f (x1+ l 1d1) является минимумом функции f ( x1+ l d1) на прямой, то

gT2d1=-gT2g1=0.    (6)

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

g2=b+Gx2, g1=b+Gx1.

Тогда, если d1 и d2=-g2+a 1d1 сопряжены, то

dT2Gd1=0,

т.е.

-gT2Gd1+a 1dT1Gd1=0,

следовательно,

(-gT2-a 1gT1)G(x2-x1)/l 1=0,

откуда

(-gT2-a 1gT1)(g2-g1)/l 1=0.

Таким образом,

-g22+a 1g21=0.

Остальные члены исчезают из соотношения (6), и, следовательно

a 1=g22/g21,

что и требовалось доказать. Это как раз и есть соотношение (5) при k=1.

Теперь перейдем к доказательству соотношений (4) и (5) по индукции, полагая, что векторы d1,d2,..,dk получены описанным выше способом и являются взаимно сопряженными.

Точка xk+1=xk+l kdk является минимумом функции f ( xk+ l dk) на прямой xk+ l kdk.

Тогда

gTk+1dk=0.    (7)

Имеем

xk+1=xk+l kdk= xk-1+l k-1dk-1+l kdk и т.д.

Таким образом,

xk+1=xj+1ki=j+1l idi; при 1< = j< = k-1,    (8)

следовательно,

Gxk+1=Gxj+1ki=j+1l iGdi,

тогда

gTk+1=gTj+1ki=j+1l idTiG при 1< = j< = k-1,

откуда

gTk+1dj=gTj+1djki=j+1l idTiGdj.

В результате преобразований имеем gTj+1dj=0 (в соответствии с соотношениями (6) и (7)) и из-за взаимной сопряженности dTiGdj=0 при j<i. Таким образом, каждое слагаемое в правой части равно нулю.

Следовательно

gTk+1dj=0 при j=1,2,...,k-1   (9)

и из соотношения (7) окончательно имеем

gTk+1dj=0 при j=1,2,...,k.    (10)

Таким образом, было доказано, что вектор gk+1 ортогонален каждому из векторов d1,d2,...,dk.

Можно также показать, что вектор gk+1 ортогонален векторам g1,g2,...,gk.

Из соотношения (10) имеем

gTk+1dj=0 при j=1,2,...,k.

Так как из предположения в начале доказательства по индукции

dj=-gj+a j-1dj-1,

то приведенное выше соотношение принимает вид

-gTk+1gj+aj-1gTk+1dj-1=0,

следовательно, -gk+1gj=0, поскольку gTk+1dj-1=0 из соотношения (10).

Таким образом

gTk+1gj=0 при j=1,2,...,k.    (11)

Доказательство по индукции будет закончено, если показать, что вектор dk+1, определенный в соотношении (4), сопряжен с векторами d1,d2,...,dk.

Для j=1,2,...,k-1 имеем

dTk+1Gdj=-gTk+1Gdj+a kdTk Gdj=

=-gTk+1Gdj

в силу взаимной сопряженности.

Тогда

-gTk+1Gdj=-gTk+1G (xj+1-xj)/l j= gTk+1G (gj+1-gj)/l j=0

с учетом соотношения (11).

Таким образом, dTk Gdj=0 при j=1,2,...,k-1, и это справедливо для любого a k. Для завершения доказательства необходимо определить ak так, чтобы выполнялось равенство

dTk Gdk=0:

Следовательно,

dTk Gdk=(-g2k+1+a kg2k)/ l k,

поскольку все другие члены из правой части исчезают в силу соотношений (10) и (11).

Следовательно, направление dk+1 будет сопряжено с направлением dk, если a k=g2k+1/g2k, что и требовалось доказать.

Таким образом, направления поиска в методе Флетчера-Ривса являются взаимно сопряженными и в данном методе минимум квадратичной функции n переменных можно найти не более чем за n шагов. Это означает, что одномерный поиск производится с нужной точностью и устраняются любые ошибки округления, которые могут возникнуть.

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

В общем случае данный метод не столь эффективен, как метод ДФП, но, тем не менее, бывает весьма полезен.

Данный метод представлен программой и может быть применен для решения прикладных задач

Функции n-переменных. Необходимые и достаточные условия.

Рассмотрим функцию n действительных переменных

f(x1,x2,...,xn) = f(x)

Точка в n-мерном евклидовом пространстве с координатами x1,x2,x3,...,xn обозначается вектором-столбцом x. Градиент функции, т.е. вектор с компонентами

f/x1, f/x2,...,f/xn,

обозначается Ñf (x) или, иногда, g(x). Матрица Гессе (гессиан) функции f(x) обозначается как G(x) и является симметричной матрицей n x n элементов вида

Gij = 2 f/xixj

Функция f(x) имеет локальный минимум в точке x0, если существует окрестность точки x0, такая, что f(x) > f(x0) во всех точках этой окрестности, т.е. существует положительная величина d, такая, что для |x - x0| < d справедливо неравенство f(x) >= f(x0).

В случае глобального минимума в точке x* для всех x справедливо неравенство f(x) >= f(x*).

При таких определениях и очевидных предположениях относительно дифференцируемости можно обобщить уравнение (2.2) из темы №2 и получить:

f(x0+h)-f(x0)= hi (x1,x2,..,xn)+ hihj (x1,x2,..,xn)+... =h + hTG(x0)h+...    (4)

Тогда, если x0 является точкой минимума функции f(x), то каждая первая частная производная f/xi (i = 1,...,n) должна обращаться в нуль в точке x0. Если это не так, то соответствующим выбором hi можно добиться того, что разность f(x0 + h) - f(x0 ) будет отрицательна.

Следовательно, необходимым условием минимума в точке x0 является уравнение

Ñf(x0)= 0;(5)

т.е.

f(x0) / xi = 0, (i=1,...,n), (6)

Тогда знак разности f(x0 + h) - f(x0) определяется членом:

1/2 hTG(x0)h    (7)

Если матрица G(x0) положительно определена, то этот член положителен для всех h. Следовательно, необходимыми и достаточными условиями минимума являются:

Ñf(x0) = 0, G(x0) положительно определена.  (8)

Необходимыми и достаточными условиями максимума являются:

Ñf(xm) = 0, G(xm) отрицательно определена.  (9)

УПРАЖНЕНИЯ И ВОПРОСЫ К ТЕМЕ №3

1. Исследовать точки перегиба функции f(x) = x12 + 4x1x2 + 5x22.

2. Исследовать точки перегиба функции

f(x) = -x12 - 6x22 -23 x32 - 4x1x2 + 6x1x3 + 20x2x3.

3. Пусть f(x) есть квадратичная функция вида f(x) = a + bTx + 1/2xTGx, где a - константа, b - вектор, не зависящий от x, а G - положительно определенная симметричная матрица, не зависящая от x. Показать, что x* = -G-1b.

4. Показать, что функция f(x) = (x1 - a)2 + (x2 - b)2 + (x2 - c)2 имеет минимум в точке (a; b; c).

5. Фирма выпускает два аналогичных товара 1 и 2. Прибыль от реализации товара составляет ciqi (i=1, 2), где ci - константа, а qi - объем реализации товара. Последний зависит от цен (p1 и p2) двух товаров. Анализ последних данных продажи дает следующие эмпирические зависимости:

где a1, a2, b1, b2 - некоторые положительные константы. Требуется определить цены p1 и p2, которые максимизируют общую прибыль.

Найти уравнения для p1 и p2, которые максимизируют прибыль, и решите их. Показать, что если решения этих уравнений положительны и 4b1 b2 > (a1 + a2)2, то эти уравнения дают оптимальные цены.

6. Найти минимум функции e-x - cos x.

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

8. Используйте программу, реализующую метод покоординатного спуска для поиска минимума функции (x1-1)2+ (x2-2)2 + (x3-3)2 и функции Розенброка 100(x2 - x12)2 + (1 - x1)2.

9. Повторите упражнения 7 и 8 с разными начальными точками.

10. Используя программу Нелдера-Мида, минимизируйте функцию Розенброка, функцию Пауэлла, двумерную экспоненциальную функцию. Используйте различные начальные значения и сравните количество вычислений функции. Особенно неудобными для минимизации являются точки (-1,2;1) для функции Розенброка и (-3;-1;0;1) для функции Пауэлла.

11. Модифицируйте метод Нелдера-Мида, заменив коэффициент, уменьшающий длину шага, т.е. вместо 10 используйте (а) 2, (б) 4, (в) 8, (г) 100. Выполните повторно упражнения 7 и 8. Сравните количество вычислений функции, необходимых для достижения конечного результата.

12. Ознакомьтесь с работой Бокса, Дэвиса и Свена, используйте предложенную ими проверку сходимости в программе Нелдера-Мида и примените результаты к функциям, упомянутым в упражнении 10.

13. Модифицируйте программу Нелдера-Мида, используя различные значения для a, b, g, проведите оценку получающихся результатов.

14. Минимизируйте функцию:

f(x1, x2) = x14 + x24 + 2x12x22 - 4x1 + 3.

15. Минимизируйте функцию:

f(x1, x2) = (x12 + x2 - 11)2 + (x1 + x22 - 7)2.

16. Решите уравнения.

x12 + x2 = 11

x1 + x22 = 7

Как можно использовать упр.15, чтобы помочь в решении этих уравнений?

17. Решите систему уравнений:

x + y + z = 6,

x2 + y2 + z2 = 14,

x3 + y3 + z3 = 36.

18. Переменные F и C связаны соотношением F = a + bC, но при измерении величины F появляется ошибка.

Определите значения a и b по следующим данным:

F

51

68

84

103

121

141

C

10

20

30

40

50

60

Используйте принцип наименьших квадратов и найдите a и b по результатам минимизации функции

S = (Fi - a - bCi)2, где i = 1:6.

19. Известно, что переменные Q и h связаны соотношением (обусловленным ошибкой в измерении Q) Q = ahn, где a и n - постоянные.

Используя приведенные ниже данные, покажите, что соотношение имеет смысл, и получите значения a и n:

h

4

6

8

10

12

Q

650

1740

3640

6360

9790

Замечание. Если Q = ahn, то ln(Q) = ln(a) + n*ln(h). Таким образом, настоящую задачу можно свести к задаче линейной регрессии с преобразованными переменными ln(h) и ln(Q). Решите ее с помощью преобразования и без него, основываясь на принципе наименьших квадратов, т.е. минимизируйте функции

1) S1 = [ln(Qi) - ln(a) - n*ln(hi)]2,

2) S2 = [Qi - ahin]2,

изменяя значения a и n. Получились ли значения одинаковыми?

20. Измените программу метода наискорейшего спуска, так, чтобы она могла выполнять линейный поиск методом Фибоначчи, и методом "золотого сечения". Примените полученную программу в решении задач, приведенных в примерах. Обратите внимание на любые возникающие трудности.

21. Найти минимум функции f(x1, x2, x3) = 3(x1 -1)2+2(x2-2)2 + (x3-3)2. В качестве начальной точки выбрать точку (9; -7; 11).

22. Попробуйте решить контрольный пример, используя различные начальные точки. Появились ли какие-нибудь трудности?

23. Пусть P0, P1 и P2 - следующие друг за другом точки, полученные методом наискорейшего спуска. Предполагается, что расчеты будут производиться быстрее, если следующий поиск произвести в направлении P0P2. Нарисуйте линии постоянного уровня вблизи минимума для иллюстрации приведенных выше соображений. Отразите эту идею в программе наискорейшего спуска.

24. Если f(x) - положительно определенная квадратичная функция двух переменных, минимум которой находится в точке (0; 0), а x0, x1 и x2 - три точки, полученные последовательно методом наискорейшего спуска, то покажите, что прямая, проходящая через точки x0 и x2, проходит через точку минимума функции (см. метод ускорения расчетов, описанный в упр.23).

25. Минимизируйте функцию

f(x) = 3(x1 - 4)2 + 5(x2 + 3)2 + 7(2x3 + 1)2,

используя: метод ДФП, метод Флетчера - Ривса. Покажите, что независимо от начальной точки для нахождения минимума в первом случае потребовалось три итерации, во втором - три поиска.

26. Минимизируйте функцию

f(x1, x2) = 1 - 2x1 - 2x2 - 4x1x2 + 10x12 + 2x22.

27. Минимизируйте функцию

x14 + x24 + 2x12x22 - 4x +3.

28. Минимизируйте функцию

(x12 + x2 - 11)2 + (x1 + x22 - 7)2.

29. Минимизируйте функцию

f(x1, x2) = x13 + x22 - 3x1 - 2x2 +2.

ПРИКЛАДНЫЕ ЗАДАЧИ К ТЕМЕ

Пример 3.1

Пример 3.2

Пример 3.3

Пример 3.1. Определение оптимальных размеров контейнера.

Пусть требуется переправить 400 м3 сыпучего материала через большую реку. Для перевозки груза необходимо построить контейнер. Известны следующие данные: стоимость каждого рейса на противоположный берег реки и обратно равна 4.2 грн.; стоимость материалов для изготовления дна контейнера составляет 20.00 грн./м2; боковых стенок контейнера - 5.00 грн./м2; крышки контейнера - 20.00 грн./м2.

Сконструируйте контейнер таким образом, чтобы минимизировать полные затраты на перевозку груза.

Для составления математической модели используем такие обозначения:

x - длина ребра дна контейнера;

y - ширина дна контейнера;

h - высота контейнера;

Sд = xy - площадь дна;

Sк = xy - площадь крышки;

Sб = 2xh+2yh - площадь боковых стенок;

Сод = 20*Sд - общая стоимость материалов для изготовления дна контейнера;

Сок = 20*Sк - общая стоимость материалов для изготовления крышки контейнера;

Соб = 5*Sб = 10h(x+y) - общая стоимость материалов для изготовления боковых стенок контейнера;

Ср - стоимость рейса;

Со = Сод + Сок + Соб + Ср (3.1.1) - общая стоимость материалов для изготовления контейнера и его перевозки, т.е. полные затраты на перевозку груза.

Таким образом, конечная функция, которую необходимо минимизировать имеет вид

Со = 20хy + 20хy + 10*h(x+y) + 4.2 min      (3.1.2)

Выразив h через ограничение xyh=400 по объему перевозимого груза получим функцию для оптимизации

Со = 40хy + 4000*(x+y)/ (xy) + 4.2 min      (3.1.3)

При аналитическом решении (первые частные производные равны нулю y=100/x2; x=100/y2) была получена точка (4,6415; 4,6415; 18,5667).

Полученные значения после вычислений по
    · наискорейшего спуска - (4,6416; 4,6416; 18,5667), Со=2589,52;
    · покоординатного спуска - (5,00; 4,50; 18,1818), Со=2593,09;
    ·модифицированному методу Хука-Дживса - (4,6415; 4,6416; 18,5667), Со=2589,52;

 

Пример 3.2. 

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

Рисунок 3.2 - Система подачи газа

C (D, P, L, r ) = 7.84 D2P1 + 450000 + 36900 + 657*106/L + 657*106/L * (r0,219 - 1) (долл./год),    (3.2.1)

где

D - внутренний диаметр труб, дюйм;

P - давление на выходе компрессора, фунт.дюйм2;

L - расстояние между компрессорными станциями, миля;

r =P1/P2 - отношение давлений на выходе и входе компрессора.

Преположим, что расход газа в единицу времени можно описать функцией

(фут3/ч),     (3.2.2)

где     f = 0.008D-1/3 - коэффициент трения.

Пусть расход газа составляет 100*106 фут3/день. Воспользуйтесь формулой (3.2.2) для исключения переменной Р1 из (3.2.1). Затем с помощью метода поиска по симплексу Нелдера-Мида и метода Пауэла найдите такие значения параметров системы, которым соответствует минимум суммарных эксплуатационных затрат в единицу времени.

Тема №4
Методы условной оптимизации. Нелинейное программирование

НЕЛИНЕЙНОЕ ПРОГРАММИРОВАНИЕ

[ Необходимые и достаточные условия оптимальности.] [ Методы решения задач нелинейного программирования.] 

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

Методы решения задач нелинейного программирования (условной многопараметрической оптимизации) подразделяют следующим образом:

  •  методы прямого поиска: 
    •  модифицированный метод Хука-Дживса; 
      •  метод комплексов; 
      •  метод случайного поиска. 
    •  методы штрафных функций; 
    •  методы линеаризации. 

Методы прямого поиска

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

 

Перед непосредственным применением методов прямого поиска необходимо провести ряд мероприятий по подготовке задачи к решению, а именно

  •  исключить ограничения в виде равенств;
    •  определить начальную допустимую точку.

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

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

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

После проведения процедуры подготовки задачи к решению следует приметить один из методов условной оптимизации. Рассмотрим методы прямого поиска:

  •  модифицированный метод Хука-Дживса; 
    •  метод комплексов; 
    •  метод случайного поиска. 

Модифицированный поиск Хука-Дживса

Метод Хука-Дживса нетрудно модифицировать и для учета ограничений.

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

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

В тексте программы модифицированного метода прямого поиска Хука-Дживса сделана попытка реализовать такую процедуру. Рассматриваемая задача формулируется следующим образом:

минимизировать f(x1, x2) = 3x12+4x1x2+5x22

при ограничениях x1³ 0, x2³ 0, x1+x2³ 4.

Программа

Минимум, равный 44, достигается в точке (3; 1) при ограничении х12=4

Для начальной точки (4;3) и при длине шага, равной единице, программой успешно решена задача минимизации.

Для начальной точки (3;4) и длины шага =1 программой также успешно решена задача минимизации.

Для начальной точки (5;6) и длины шага =1 задача не решена, так как программа остановилась в точке (1;3), т.е. на активном ограничении, и выдала неверный результат.

Аналогичные неутешительные результаты были получены для начальной точки (5;6) и длины шага =0,5. Неверное решение было найдено в точке (1,5;2,5). Для начальной точки (4;3) и длины шага 0,5 программа работала нормально, но было получено неверное решение в точке (2,5;1,5).

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

Метод комплексов

Алгоритм:

Заданы границы значений всех переменных xiL, xiU, i=1,2,..., N (размерность задачи), допустимая точка xo, параметр отображения a (рекомендуется a =1,3) и параметры окончания вычислений e и d .

Шаг 1. Построение начального комплекса, состоящего из P допустимых точек (рекомендуется P=2N). Для каждой точки p = 1, 2,...,P-1

  •  случайным образом определить координаты xp;
    •  если xp - недопустимая точка, то найти центр тяжести xцт уже найденных точек и положить xp = xp + (xцт - xp); повторять процедуру до тех пор, пока xp не станет допустимой;
    •  если xp - допустимая точка, повторять до тех пор, пока p=P;
    •  вычислить W(xp) для p=0,1,...,P-1.

Шаг 2. Отражение комплекса:

  •  выбрать точку xR, для которой W(xR) = max W(xp) = Wmax (решается задача минимизации);
    •  найти центр тяжести xцт и новую точку xm = xцт + a (xцт - xR);
    •  если xm - допустимая точка и W(xm)< Wmax, уменьшить в два раза расстояние между xm и центром тяжести xцт, продолжать поиск, пока W(xm)<Wmax;
    •  если xm - допустимая точка и W(xm)<Wmax, то перейти к шагу 4;
    •  если xm - недопустимая точка, то перейти к шагу 3.

Шаг 3. Корректировка для обеспечения допустимости:

  •  если xim<xiL(нижняя граница допускаемой области), то положить xim = xiL;
  •  если xim>xiU(верхняя граница допускаемой области), то положить xim = xiU;
  •  если xm - недопустимая точка, то уменьшить в два раза расстояние до центра тяжести; продолжать до тех пор, пока xm не станет допустимой точкой.

Шаг 4. Проверка условий окончания вычислений:

положить и   ;

если и , то вычисления прекратить; в противном случае перейти к шагу 2a.

Программа

Методы случайного поиска

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

xi = xiL +ri (xiU - xiL), i=1,2,...,N,

где

ri - случайная величина, равномерно распределенная на интервале [0,1].

После проверки каждой точки на допустимость вычисляются значения целевой функции, которые сравниваются с наилучшим значением, полученным к данному моменту. Если текущая точка дает лучшее значение, то она запоминается, в противном - отбрасывается. Процесс прекращается после заданного числа итераций N или по исчерпанию заданного машинного времени. Для получения 90% доверительного интервала величиной e i (xiU - xiL), где 0<e <1, для переменной xi совместный случайный поиск требует

испытаний. При N=5, e i=0,01 число испытаний оценивается в 2,3 1010, что исключает возможность непосредственного использования метода.

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

Заданы границы значений всех переменных xiL, xiU, i=1,2,..., N (размерность задачи), начальные допустимая точка xo и интервал поиска D xo = xiU - xiL, количество серий Q, количество точек в серии P и параметр окончания вычислений e . Для каждой из серий, начиная с q = 1, необходимо выполнить следующие действия.

Шаг 1. Для i = 1,2,...,N вычислить

xip = xiq-1 +ri D xq-1,

где

ri - случайная величина, равномерно распределенная на интервале [-0,5;0,5].

Шаг 2.

  •  Если xp - недопустимая точка и p < P, то повторить шаг 1.
    •  Если xp - допустимая точка, то запомнить xp и W(xp) и
      •  если p < P, то повторить шаг 1.
      •  если p = P, то найти среди всех допустимых точек xp точку с наименьшим значением W(xp) и положить ее равной xq.

Шаг 3. Уменьшить интервал поиска, полагая D xiq = e i D xiq-1.

Шаг 4.

  •  Если q > Q, то закончить вычисления.
  •  В противном случае увеличить q = q + 1 и продолжить вычисления, начиная с шага 1.

Описание штрафных функций.

Основная идея метода штрафной функции состоит в преобразовании задачи минимизации функции Z=f(x) c соответствующими ограничениями, наложенными на X, в задачу поиска минимума без ограничений функции

Z=f(x)+p(x).

Функция P(x) является штрафной. Необходимо, чтобы при нарушении ограничений она "штрафовала" функцию Z, т.е. увеличивала ее значение. В этом случае минимум Z будет находиться внутри области ограничений. Функция P(x), удовлетворяющая этому условию, может быть не единственной.

Задачу минимизации можно сформулировать следующим образом:

минимизировать функцию

z=f(x)

при ограничениях cj(x)>0, j=1,2,..,m.

Замечание.Ограничение вида "меньше или равно", h(x)0, всегда может быть записано как -h(x)0, поэтому в приведенной выше формулировке нет потери общности.

Функцию P(x) удобно записать следующим образом:

P(x)=r.

где r- положительная величина. Тогда функцию Z= j(x, r) принимает вид

Z=(x, r)=f(x)+r. (1)

Если x принимает допустимые значения, т.е. значения, для которых cj(x)0, то Z принимает значения, которые больше соответствующих значений f(x) (истинной целевой функции данной задачи), и разность можно уменьшить за счет того, что r может быть очень малой величиной. Но если x принимает значения, которые хотя и являются допустимыми, но близки к границе области ограничений, и по крайней мере одна из функций cj(x) близка к нулю, тогда значения функции P(x) и, следовательно, значения функции Z станут очень велики. Таким образом, влияние функции P(x)состоит в создании "гребня с крутыми краями" вдоль каждой границы области ограничений. Следовательно, если поиск начинается из допустимой точки и осуществляется поиск минимума функции (x, r) без ограничений, то минимум, конечно, будет достигаться внутри допустимой области для задачи с ограничениями. Полагая r достаточно малой величиной, для того чтобы влияние P(x) было малым в точке минимума, мы можем сделать точку минимума функции (x, r) без ограничений совпадающей с точкой минимума функции f(x) с ограничениями.

В общем случае невозможно аналитически определить положение минимума функции x, (x, r), рассматривая ее как обычную функцию от r. Для его определения необходимо обратиться к численным методам.

Следует отметить, что если целевая функция f(x) выпукла, а функция cj(x) вогнута, то функция (x, r), заданная уравнением (1), также является выпуклой функцией в области ограничений, которая сама является выпуклой. Следовательно, (x, r) имеет для данного значения r единственный минимум.

Если x1 и x2-точки, принадлежащие допустимой области, т. е. Cj(x1)0 и cj(x2)0 для j=1,2,...,m, то при 0<<1 справедливо неравенство

cj( x2+(1-) x1)cj(x2)+(1-)cj(x1) 0,

так как функция cj(x) выпуклая. Следовательно, допустимая область является выпуклой.

Таким образом, точка x2+(1-)x1 при 0<<1 также является допустимой.

Кроме того, функция 1/cj(x) является выпуклой для всех x, которые удовлетворяют неравенству cj(x) 0. Если h(x)=1/cj(x), то

h(x)=.

Следовательно, гессиан функции h(x) имеет вид

H(x)=-+,

где c(x)i k= есть гессиан функции cj(x). Тогда, если p(x)- произвольный вектор, то справедливо равенство

pth(x)p=- ,

где всегда pth(x)p>0, так как c(x)- отрицательно определенная матрица ввиду того, что cj(x)- выпуклая функция и cj(x) 0. Тогда матрица H(x) положительно определена и 1/cj(x) выпукла во всей области.

Предположим, что x1*,x2*,...,xk*- минимальные точки функции (x, r) для убывающей последовательности значений r1, r2,...,rk,..., стремящейся к нулю. Тогда последовательность точек x1*,x2*,...,xk*,... сходится к оптимальному решению задачи с ограничениями при rk0.

Следовательно, lim xk*=x* и lim при rk0,

где x* - минимальная точка функции f(x) при наличии ограничений.

Полученный результат можно доказать следующим образом. Поскольку f(x) непрерывная функция и f(x*)f(x) для всех допустимых точек, то, задав произвольное достаточно малое значение eps, можно найти допустимую точку x', такую, что

f(x')<f(x*)+eps/2.  (2)

Поскольку rk- убывающая последовательность, стремящаяся к нулю, можно найти такое значение K, что при kK справедливо неравенство

rk (3)

Поскольку P(x)>0 из определения функции (x,r), имеем

f(x*)min(x, rk)=(xk*,rk),

где xk*- минимальная точка функции (x, rk) для задачи без ограничений.

Кроме того, если k>K, то rk<rK и справедливо неравенство

(xk*,rk)(xk*,rk).

Это следует из того, что, поскольку xk* минимизирует функцию (x, rk) и в любой другой точке области x, в частности в точке xK*, функция будет принимать значение, больше чем (xk*,rk).

Поэтому

(xk*,rk)=f(xk*)+rk>f(xk*)+rk,

поскольку rk<rK.

Следовательно,

(xk*,rk)>(xk*,rk).

Тогда

f(x*)(xk*,rk)(xK*,rk)<(xK*,rK). (4)

Но так как значение xK* минимизирует функцию (x, rK), то

(xK*,rK)(x',rK)=f(x')+rk. (5)

Следовательно, из уравнений (4) и (5) получим

f(x*)(xk*,rk)f(x')+rk.

Из уравнения (3) следует, что

f(x')+rkf(x')+eps/2.

Тогда из уравнения (2) следует, что

f(x*)(xk*,rk)<f(x*)+eps/2+eps/2.

Поскольку eps может быть выбрано произвольно малым, всегда можно найти такое значение k, при котором

f(x*)<(xk*,rk)<f(x*)+eps.

Таким образом, при k (rk0)

lim (xk*,rk)=f(x*) при rk0.

Из приведенного выше доказательства следует, что при rk0

f(xk*)f(x*) и rk0. (6)

f(x1*),f(x2*),...,f(xk*)образуют убывающую последовательность, такую, что f(x*k+1)<f(xk*).

Очевидно, что если функция f(x) выпукла, а функция cj(x) при j=1,2,...,n вогнута, то функция f(x) при наличии ограничений имеет единственный минимум.

Штрафы. Виды штрафов.

С помощью штрафных функций

P(x,R) = W(x) + W (R,g(x),h(x)),

где
      R - набор штрафных параметров;
     
W - штраф,

исходная задача условной оптимизации преобразуется в последовательность задач безусловной оптимизации. Штраф W определяется так, чтобы допустимые точки задачи имели преимущество перед недопустимыми в отношении безусловной оптимизации штрафной функции. Здесь штраф как бы создает вдоль границы допустимой области барьер из бесконечно больших значений функции P.

К штрафу выдвигаются следующие требования:

  •  решение подзадач должны стремиться к решению исходной задачи нелинейного программирования ;
  •  сложность оптимизации P(x, R) должна быть такого же порядка, что и W(x).

Методы штрафных функций классифицируются в соответствии со способами учета ограничений - неравенств g(x), так как ограничения-равенства h(x) учитываются во всех методах одинаково с помощью квадратичного штрафа

W = R{h(x)}2.

При рассмотрении любой штрафной функции требуется выбрать начальное значение R и изменять его после решения каждой подзадачи безусловной оптимизации с тем, чтобы обеспечить сходимость. Для квадратичного штрафа, учитывающего ограничения - равенства, представляется целесообразным начинать с R=0, а затем последовательно увеличивать R на некоторое D R или использовать возрастающие степени какого-либо числа, например 10. В результате получаемые точки будут все точнее и точнее удовлетворять ограничениям.

Для учета ограничений - неравенств используют следующие штрафы:

"Бесконечный" штраф

W = 1020,

где

- множество индексов нарушенных ограничений gj(x)<0 при j Є J.

Логарифмический штраф

W = -R ln[g(x)].

Отрицательный штраф исключают положив W = 0 для таких x, что g(x)>1.

Логарифмический штраф - барьерная функция, не определенная в недопустимых точках. Итерационный процесс следует начинать из допустимой начальной точки при положительном начальном значении R (R=10 или R=100). После решения каждой подзадачи условной оптимизации параметр R следует уменьшать и в пределе устремить к нулю.

Штраф обратной функции

W = R [1/g(x)].

Итерации следует начинать с начальной допустимой точки при положительном R, значения которого в пределе должно стремиться к нулю.

Штраф квадрата срезки

W = R [g(x)]2 ,

где

В данном методе недопустимые точки не создают проблем (в отличие от предыдущих), поэтому он весьма удобен. Кроме того функция P(x,R) определена и непрерывна всюду. Вычисления следует проводить с положительными Ri; после решения очередной подзадачи безусловной оптимизации R необходимо увеличивать.

Алгоритм методов штрафных функций:

Шаг 1. Задать значения N, J, K, e1, e2, e3, xo, Ro,

где

e1, e2, e3 - соответственно, параметры окончания процедур одномерного и многомерного поиска безусловной оптимизации а также работы алгоритма штрафных функций;

xo - начальное приближение для x*;

Ro - начальный выбор штрафных параметров.

Шаг 2. Построить P(x, R) = W(x) + W (R, g(x), h(x)).

Шаг 3. Найти xt+1 минимизирующее значение P(xt+1,Rt) при фиксированном Rt. В качестве начальной точки использовать xt, а в качестве параметра окончания шага - константу e2 (возможно и e1).

Шаг 4. Проверить, выполняется ли условие

| P(xt+1, Rt) - P(xt, Rt-1) | < e3.

если "да" - положить xt+1=xT и закончить процесс решения;

если "нет" - перейти к следующему шагу.

Шаг 5. Положить Rt+1=Rt+D Rt в соответствии с используемым правилом пересчета, после чего вернуться к шагу 2.

Методы штрафных функций.

Идея метода штрафных функций состоит в следующем. Исходную задачу математического программирования

                                    min j (x) = m                                  (1)
x
є X

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

M(x,b ) = j (x) + (1/b )y (x),

либо задачей минимизации функции M(x,b) на некотором множестве Г, более простой структуры с точки зрения численных методов минимизации, чем исходное множество Х.
При этом функцию y (x) выбирают таким образом, чтобы решение задачи

min M(x, b ) = m(b ) (2)
x є Г

сходилось при b(r) 0 к решению исходной задачи, или если это обеспечить не удается, то по крайней мере чтобы

  m(b) m при b(r) 0.

Итак, метод штрафных функций состоит в построении функции М(х, b), в выборе последовательности {bi} и в решении соответствующих этой последовательности задач вида (2).
Имеются два главных, в основе своей различных, подхода к построению построению штрафных функций. Первый подход, известный как метод внешних штрафных функций, был предложен Курантом в 1943 г., когда почти все хорошо известные сейчас алгоритмы нелинейного программирования еще только разрабатывались. В методе внешних штрафных функций штрафы выбираются так, что точки вне С (замкнутое множество в пространстве R
n) последовательно по мере возрастания и становятся все более "невыгодными" и либо mi(b) = 0 для всех bєC, либо mi(b) 0 при i (r) для всех bє C.

Второй подход называется методом внутренних штрафных функций. В нем штрафные функции выбираются так, чтобы оптимальные для задач (2) точки bi принадлежали внутренности С0 множества С. Благодаря этому свойству штрафных функций методы внутренних штрафных функций известны также как барьерные методы, поскольку в процессе решения задачи они "отталкивают" от границы множества С точки b . Поэтому начальная точка b0, используемая при решении задачи (2), должна выбираться из С0.
Основные свойства методов штрафных функций: при достаточно слабых предположениях они строят последовательности, сходящиеся к точкам, оптимальным для задачи (2) определяемой (1). В то же время все результаты опираются на значительную идеализацию ситуации, поскольку мы предполагаем, что можем решать задачи минимизации без ограничений (2). Возникает две проблемы. Первая из них состоит в решении вопроса о "начале и конце", или, точнее, о прерывании процесса построения минимизирующей последовательности {b
i}. Вторая проблема касается достоинств различных методов штрафных функций.
Пусть задано число e > 0 и мы хотим найти точку . Все методы штрафных функций строят точки b
i, являющиеся оптимальными решениями для задачи (2).
Таким образом, в принципе, существует такое целое число j, что точка b
j удовлетворяет нашим требованиям. Предположим, что это целое число j, которое на практике оказывается очень большим, нам известно, и  допустим, что мы хотим найти минимум функции
f
0(b ) + mj(b ), исходя из некоторого начального приближения b0, с помощью одного из алгоритмов, требующего вычисления градиента. Весьма вероятно, что f0(b0) будет слишком мало по сравнению с mj(b ), если mj - внешняя штрафная функция, или слишком велико, если mj(b ) - внутренняя штрафная функция.
Обычная практика состоит в том, чтобы начать со штрафа умеренной величины и затем увеличивать его в процессе вычислений.
Существует много ситуаций, когда имеет смысл комбинировать методы внешних и внутренних штрафных функций и создавать смешанные методы.

Прикладной пример 1     

Постановка задачи

В качестве примера применения метода штрафных функций рассмотрим задачу оптимального проектирования щита люка. Отверстие прямоугольной формы шириной l0=600 см должно быть покрыто полыми брусками из алюминия длиной l0 и шириной в 6 см. Покрытие должно выдерживать максимальную удельную нагрузку до 1000 кг/м3 .

a)

б)

Рисунок 1 - Условные обозначения размеров щита.

В качестве переменных задачи (конструктивных параметров)- приняты следующие: x1=tf - толщина стенок и x2=h - высота бруска (рис. 1а,б).

Приняты также следующие предположения:

  1.  Коэффициент Пуассона v=0,3 ;
  2.  Материал обладает линейной упругостью с модулем Юнга E=7*105 кг/см2 
  3.  Остаточная деформация должна отсутствовать.

При построении модели формируются следующие ограничения:

1) На максимальное касательное напряжение tmax. Полагая, что максимум поперечной силы Q = 1800 кг, для максимального касательного напряжения имеем:

= 450 (1)

2) На максимальное изгибающее напряжение sbmax . На основе предварительного анализа конструкции получен максимальный изгибающий момент M = 2,7*104кг/см. Для изгибающего напряжения, действующего на стенки, имеем

(2)

3) На максимальный прогиб dmaxpl0 , где p=0,025. Величина прогиба определяется по формуле (3)

где q1=q0b, и для величины I принята аппроксимация (4)

4) На изгибающее напряжение в стенках. Изгибающее напряжение в стенках выражается формулой

(5)

В результате задача оптимального проектирования формируется с учетом следующих ограничений:


;
;
; (7)
;
h >= 0; t
f >= 0;

Критерием оптимальности является масса щита на единицу длины F(tf, h) = h +120 tf (8)

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

Алгоритм нахождения оптимального решения методом случайного поиска.

Введем обозначения:

n- количество испытаний;

i- текущий номер испытания;

temp-текущее наилучшее значение функции цели;

f-функция цели;

h,tf- конструктивные параметры;

h0,tf0- начальные значения, удовлетворяющие условию (7).

  1.  Вычисляем значение функции f = h0+120*tf0; присваиваем temp=f;
  2.  Увеличиваем i и находим новые случайные значения для h,tf.
  3.  Если h,tf удовлетворяют условию (7) то находим значение функции f. Если f < temp, то присваиваем temp=f;
  4.  Если i <= n переходим к пункту 2.

Алгоритм нахождения оптимального решения методом штрафных функций.

Рассмотрим задачу поиска локального минимума критерия оптимальности F(tf, h) в области, ограниченной системой неравенств (7). Введение обобщенного критерия оптимальности по методу штрафных функций производится с помощью непрерывной функции

(9)

Обобщенным критерием оптимальности согласно методу штрафных функций является выражение TT(tf, h, t) = h +120tf + tQ(x), где t - некоторое положительное число, называемое коэффициентом штрафа.

Рассматривается некоторая неограниченная, монотонно возрастающая последовательность {tk}, k=1,2,... положительных чисел. Для первого элемента этой последовательности с помощью метода покоординатного спуска отыскивается безусловный локальный минимум функции TT(tf,h,t). Пусть этот минимум достигается при значениях (tf1,h1).

Вектор (tf1,h1) используется как начальное приближение для решения задачи поиска минимума функции TT(tf,h,t2) где t2>t1 и т.д. Таким образом, решается последовательность задач минимизации функций TT(tf,h,th), k=1,2,..., причем результат предыдущей оптимизации используется в качестве начального приближения для поиска последующей.

Блок-схема метода случайного поиска

Текст программы

unit cover;
interface
uses

Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs,
StdCtrls, Menus, ComCtrls;

type

TForm1 = class(TForm)
Button1: TButton;
Label1: TLabel;
Label2: TLabel;
Label3: TLabel;
Label4: TLabel;
Label5: TLabel;
Label6: TLabel;
MainMenu1: TMainMenu;
exit: TMenuItem;
Label7: TLabel;
Label8: TLabel;
ProgressBar1: TProgressBar;
ProgressBar2: TProgressBar;
procedure Button1Click(Sender: TObject);
procedure exitClick(Sender: TObject);

private
{ Private declarations }
public
{ Public declarations }
end;

var

Form1: TForm1;

implementation

{$R *.DFM}
const

l0=600;
v=0.3;
e=7e5;
q=1800;
ts=0.5;
m=2.7e4;
b=6;
sbmax=700;
p=0.025;
q0=0.1;
tmax=450;

var h,tf:real;
t:integer;

function f1:real;

begin
f1:=1-tmax*2*ts*h/q;
end;


function f2:real;

begin
f2:=1-sbmax*(b*tf*h)/m;
end;


function f3:real;

begin
f3:=1-p*l0*e*0.5*b*tf*sqr(h)/(0.0158*q0*b*sqr(sqr(l0))) ;
end;


function f4:real;

begin
f4:=1-e*sqr(tf)*b*tf*h/(1000*m);
end;


function f:real;

begin
f:=h+120*tf;
end;


function max(x1,x2:real):real;

begin
if x1>x2 then max:=x1
else max:=x2;
end;


function min(x1,x2:real):real;

begin
if x1
else min:=x2;
end;



function TT:REAL;

begin
tt:=f+t*(sqr((f1+abs(f1))/2)+sqr((f2+abs(f2))/2)+sqr((f3+abs(f3))/2)+sqr((f4+abs(f4))/2));
end;


Procedure lmin;
label l1,l2;
const eps=1e-9;
var temp:real;
step:real;
count:longint;

begin
count:=0;
for t:=1 to 2000 do

begin
inc(count);
repeat

step:=1;
repeat

repeat

temp:=tt;
h:=h+step;
if h<0 then begin
h:=0;
goto l1;
end;

until tt>temp;
step:=-step/2

until abs(temp-tt)
l1:

step:=0.01;
repeat

repeat

temp:=tt;
tf:=tf+step;
if tf<0 then begin
tf:=0;
goto l2;
end;

until tt>temp;
step:=-step/2;

until abs(temp-tt)
l2:

until abs(temp-tt)
Form1.ProgressBar2.Position:=round(count/20);

end;
form1.Label4.caption:='H= '+floattostr(h);
form1.Label5.caption:='Tf= '+floattostr(tf);
form1.Label6.caption:='F= '+floattostr(f);

end;


procedure rlmin;
var i,c,count:longint;
temp,temp1,ht,tft:real;
hmin,hmax,tfmin,tfmax,minold:real;

begin
randomize;
hmin:=0;
hmax:=100;
tfmin:=0;
tfmax:=1;
count:=0;

for c:=1 to 10 do

begin

temp:=f;
ht:=h;
tft:=tf;

for i:=1 to 10000 do
begin
inc(count);
h:=hmin+random*(hmax-hmin);
tf:=tfmin+random*(tfmax-tfmin);
if (f1<=0 )and (f2<=0) and (f3<=0) and (f4<=0) then
begin
temp1:=f;
if (f
begin
temp:=f;
ht:=h;
tft:=tf;
end

end;


end;


h:=ht;
tf:=tft;

minold:=hmin;
hmin:=h-(hmax-minold)/20;
hmax:=h+(hmax-minold)/20;

minold:=tfmin;
tfmin:=tf-(tfmax-minold)/20;
tfmax:=tf+(tfmax-minold)/20;

Form1.ProgressBar1.Position:=round(count/1000);
end;


form1.Label1.caption:='H= '+floattostr(h);
form1.Label2.caption:='Tf= '+floattostr(tf);
form1.Label3.caption:='F= '+floattostr(f);
end;


procedure TForm1.Button1Click(Sender: TObject);

begin
h:=40;tf:=0.7;
rlmin;
h:=40;tf:=1;
lmin;

end;


procedure TForm1.exitClick(Sender: TObject);

begin
Close;
end;


end.

Инструкция для пользователя

Программа написана на языке Delphi 3.0. Для запуска программы необходимо запустить на исполнение файл project1.exe. В откомпилированном виде программа занимает 195 Кб. Для нормального функционирования программе необходима операционная система Windows'95. Поэтому минимальные требования к аппаратному обеспечению следующие:

  •  · Процессор 80386/80387 и выше.
  •  · Объем ОЗУ более 8МБ.

Программа не требует для своей работы ввода данных. После запуска программы нужно щелкнуть мышью по кнопке "PROCESS" и ждать до появления результата. Выход из программы осуществляется с помощью пункта меню "Exit".

Контрольный пример. 

В качестве контрольного примера возьмем следующие значения: l0=600; v=0.3; E=7*105; Q=1800; ts=0.5; M=2.7*104; b=6; sbmax=700; p=0.025; q0=0.1; tmax=450.

Вычисление методом штрафных функций дает следующие результаты: h=25.29   tf=0.632 (величины даны в сантиметрах) F=101.14

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

Таблица 5.1 

Количество итераций

Полученные значения

1000

h=24.81 tf =0.638 F=101.32

10 000

h=25.49 tf = 0.632 F=101.44

100 000

h=25.44 tf =0.628 F=101.33

1 000 000

h=25.45 tf = 0.632 F= 101.30

10 000 000

h=25.39 tf = 0.632 F= 101.30

Заключение

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

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

  Прикладной пример 2 

Оптимальное проектирование диска турбины

Необходимо спроектировать диск турбины со следующими заданными свойствами:
· Масса диска W должна быть минимальна
· Должно удовлетворяться требование по прочности

где s0 - предельно допустимое давление, st и sq- радиальное и тангенциальное давление соответственно.
· Заданы ограничения на пределы изменения толщины диска. Для ее решения использовать метод случайного поиска и метод штрафных функций.

Постановка задачи оптимального проектирования


Диск турбины можно рассматривать как вращающийся круглый диск переменной толщины (рис.1). Масса такого диска выражается с помощью формулы:

Где a1 и am - внутренний и внешний радиусы соответственно;
р- плотность; h(r) -толщина диска на расстоянии r от оси вращения.

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

Уравнение равновесия сил для вращающегося диска имеет вид:

где st и sq - радиальное и тангенциальное давление соответственно, а w - частота вращения диска. Это уравнение справедливо в предположении радиальной симметрии сил в плоскости, ортогональной к оси вращения.

Cилы могут быть выражены через величину радиального смещения u(r) c помощью следующих соотношений:

er=du/dr; eq=u/r;       (3.4)


где e
r и eq - радиальная и тангенциальная деформации; E- модуль Юнга и n- отношение Пуассона.

Подставляя уравнения (3) и (4) в (2), получаем:

Отсюда следует, что для определения u(r) необходимо знать геометрию диска, т.е. h=h(r), a1<=r<=am. В этом случае u(r) определится как решение краевой задачи для уравнения (3.5) с граничными условиями:

В результате решения краевой задачи и последующего применения формул (3.3) и (3.4) можно определить давления st и sq.

Кроме условия равновесия сил при вращении диска, должны удовлетворяться требования по прочности. Соответствующее ограничение имеет вид:

   (3.7)

где s0 - предельно допустимое давление.

Наконец, ограничения на пределы изменения толщины диска заданы следующим образом:


где
e - нижний предел толщины, заданный из технологических соображений; am-am-1=const.

Полученная в результате задача минимизации выражения (3.1) при ограничениях (3.5)-(3.8) может быть решена методами оптимального проектирования.

Описание алгоритма

Значения функций ограничений и критерия оптимальности полностью определены, если задана зависимость h(r). Чтобы представить задачу оптимального проектирования как задачу варьирования многомерного вектора конструктивных параметров, производится аппроксимация зависимости h(r) некоторой кусочно-линейной функцией. С этой целью производится разбиение радиуса диска на m частей:

a1 < a2 <...< am-1 < am      (3.9)

На каждом отрезке строится линейная функция, аппроксимирующая соответствующий участок изменения h(r) (см. рис. 1)

Если толщина диска задана в точках разбиения величинами

h(aj)=bj, j=1,2,...,m,      (3.10)


то кусочно-линейная аппроксимация полностью определена формулой

Подставляя формулу (3.11) в критерий оптимальности (3.1), получаем следующее выражение:

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

Задание разбиения диска фиксированными точками a1, a3, a4,..., am и переменной a2 дает возможность перейти к следующей постановке задачи оптимального проектирования диска.

Предположим, что значения b1 = b2 и bm = bm-1 фиксированы. Тогда остальные проектируемые параметры составляют вектор x = (b3, b4,..., bm-2, a2) и размерности m-3.

Ограничения на геометрию диска являются следствием (3.8) и записываются в виде

bj e1, j = 3...m-2; a1+ e3 a2 a3 - e2     (3.13)

где e1, e2, e3 - заданные величины. В векторной форме эти ограничения могут быть представлены неравенством

l x u           (3.14)



Чтобы представить ограничение (3.5) в соответствии с принятым разбиением, запишем его отдельно для каждого участка. При этом будем предполагать разбиение таким, что изменением h(r) внутри каждого участка можно пренебречь.

Уравнение (3.5) для каждого участка принимает более простую форму:

что позволяет получить общее решение в виде

где с1 и с2 - постоянные, получаемые следующим методом:

В начальной точке


Пусть rисходн = a[1]; E= 2*1010; n=0,3
Тогда приравняв e
qисходн и er исходн и подставив числовые значения, получим:

следовательно с1=1,4*10-3 и c2=0,084*10-3

Отсюда с помощью (3.3) и (3.4) получаем выражение для сил вращения:

Для проектируемой конструкции определяются также силы sr, sq в точках r1,r2,...,rn. Ограничение (3.7) записывается отдельно для каждого интервала [rj-1,rj]. Имеем:

откуда ввиду (3.7) получаем:

Вектор функций ограничений может быть записан в виде

Пределы его изменения в соответствии с условием (3.16) записываются с помощью неравенства

где L=(0, 0,..., 0) и U=(t0, t0, t0).

Таким образом, задача оптимального проектирования диска паровой турбины сводится к поиску минимума критерия оптимальности (3.12) при ограничениях (3.16) и (3.17) путем варьирования вектора x=(b3,..., bn-2, a2).

Алгоритм метода случайного поиска

Введем обозначения:

numsteps - количество испытаний;
i - текущий номер испытания;
temp -текущее наилучшее значение функции цели; (вначале присвоим temp=1E+38)
W -функция цели;
- параметры проектирования;
- текущий наилучший набор параметров проектирования.


1. Случайным образом генерируем набор элементов массива b.
2. Проверяем, удовлетворяет ли данный набор ограничениям.
3. Если "да", то находим значение целевой функции, соответствующей данному набору;     переходим к пункту 4. Иначе если "нет" и i < numsteps, увеличиваем i, переходим к пункту 1.
4. Если W < temp то temp=W и запоминаем текущее лучшее значение =
5. Если i < numsteps, увеличиваем i, переходим к пункту 1.
6. Вывод результатов. Останов.
Алгоритм метода штрафных функций

Рассмотрим задачу поиска локального минимума критерия оптимальности W в области, ограниченной системой неравенств (3.16)-(3.17). Введение обобщенного критерия оптимальности по методу штрафных функций производится с помощью непрерывной функции

Обобщенным критерием оптимальности согласно методу штрафных функций является выражение T = W + RQ(x), где R - некоторое положительное число, называемое коэффициентом штрафа.

Рассматривается некоторая неограниченная, монотонно возрастающая последовательность {Rk}, k=1,2,... положительных чисел. Для первого элемента этой последовательности с помощью метода покоординатного спуска отыскивается локальный минимум функции T. Пусть этот минимум достигается при значениях (b*,R1).

Вектор (b*,R1) используется как начальное приближение для решения задачи поиска минимума функции T где R2>R1 и т.д. Таким образом, решается последовательность задач минимизации функций T(b*,Rk), k=1,2 ..., причем результат предыдущей оптимизации используется в качестве начального приближения для поиска последующей

Блок-схема метода случайного поиска

Методы внешних штрафных функций.

Причина, по которой используется множественное число (методы, а не метод), состоит в том, что выборы различных семейств штрафных функций приводят к различным алгоритмам. Таким образом, методы внешних штрафных функций образуют целый класс.
Определение. Пусть С-замкнутое множество в Rn. Последовательность непрерывных функций mi є Rn(r) R1, i = 0, 1, 2,..., называется последовательностью внешних штрафных функций для множества С, если:
    1) m
i'(b ) = 0 для всех b є C, i = 0, 1, 2,...,
    2) m
i'(b ) > 0 для всех b C, i = 0, 1, 2,...,
    3) m
i+1'(b ) > miє(b ) для всех b C, i = 0, 1, 2,...,
    4) m
i'(b ) при i для всех b C.

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

Алгоритм:

Шаг 0. Выбрать e> 0, a О (0,1/2), g > 1, r > 0 и b0 є Rn.
Шаг 1. Положить e
0 = e , положить j = 0 и положить i = 0.
Шаг 2. Вычислить

Шаг 3. Если ||h(bj; ei)||> ei то перейти к шагу 4; иначе положить
e
i+1 + 1=ei /b, положить bi = bj и перейти к шагу 2.
Шаг 4. Используя один из изученных ранее методов (например, один из градиентных методов) , вычислить такое число
lj, что
-
l (1- a )|| h(bj; ei)||2 <f0(bj +ljh(bj; ei)) + (1/ei)m''(bj +lj h(bj; ei))- f0(bj) - -(1/ei )m (bj) < -lja|| h(bj; ei)||2.
Шаг 5. Положить b
j+1 = bj + ljh(bj; ei), положить j = j + 1 и перейти к шагу 2.    

Методы внутренних штрафных функций.

Изложим теперь методы внутренних штрафных функций для решения задачи (2), определяемой (1). Как и в случае методов внешних штрафных функций, представленные результаты тривиально обобщаются на бесконечномерную задачу (2), если заменить Rn на L во всех приведенных ниже утверждениях.

Определение. Пусть С-множество в Rn. Последовательность непрерывных функций m ''i :C0 R1, i = 0, 1, 2, ..., называется последовательностью внутренних, шрафных функций для множества С, если

0 < m ''i+1(b ) < m ''i(b ) для всех b є C0,i = 0, 1, 2,...,
m ''
i(b ) 0 при i для всех b є C0,
m(b
j) при j для любой последовательности {bj} є C0, для которой bj(r)b* є dN (граница множества С).

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

Алгоритм:

Шаг 0. Выбрать e> 0, a є (0,1/2), g > 1, r > 0 и b0 є C0.
Шаг 1. Положить e
0 = e , положить j = 0 и положить i = 0.
Шаг 2. Вычислить

Шаг 3. Если ||h(bj; ei)||> ei то перейти к шагу 4; иначе положить ei+1 = ei /b, положить bi = bj и перейти к шагу 2.
Шаг 4. Используя один из изученных ранее методов (например,один из градиентных методов) , вычислить такое число
lj, что bj + lj h(bj; ei) є C и
-
l (1- a )|| h(bj; ei)||2< f0(bj +ljh(bj; ei)) + eim''(bj + lj h(bj; ei))- f0(bj) - ei m'' (bj)< -lja|| h(bj; ei)||2.

Шаг 5. Положить bj+1 = bj + ljh(bj; ei), положить j = j + 1 и перейти к шагу 2.

Комбинированные методы штрафных функций.

При соответствующих предположениях методы внешних и внутренних штрафных функций можно комбинировать, в результате чего получится смешанный метод решения задачи (2), определяемой (1).
Наилучшие результаты дает комбинация методов внутренних и внешних штрафных функций. Задавая начальное приближение b
0 для задачи (2) и полагая С = { b| f j (b ) <? 0, j = 1,2,...,m, rj(b ) = 0, j = 1,2,...,l}, можно назначить внешние штрафы за нарушение ограничений rj(b ) = 0, j = 1,2,...,l и за нарушение ограничений f j (b ) < 0, когда f j (b ) > 0. Можно, кроме того, назначить внутренние штрафы для тех ограничений f j (b )< 0, для которых

f j (b) < 0.

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

Алгоритм:

Примечание. Этот алгоритм решает задачу

min {f° (b ) | fj (b ) < 0, i = 1, 2, ... m, rj(b } = 0,j = 0, 1, 2,...,l},

где функции f и r непрерывно дифференцируемы. Для него требуется предположение о том, что множество С = {b є Rn | fi (b )< 0,
r =1, 2,..., m, r
j(b ) = 0, j=1, 2,...,l} удовлетворяет условиям регулярности Куна-Таккера.
Шаг 0. Выбрать b
0 є Rn , выбрать a, a',a'' є (0; 0,5), выбрать b є (0,5; 0,8).
Шаг 1. Положить b = b
0 и положить i=0.
Шаг 2. Определить множества индексов

I' = {j є {l, 2,..., m} | fj(b ) > 0},
I'' = { j є {l, 2,..., m} | f
j(b ) < 0}.

Шаг 3. Определить внешние и внутренние штрафные функции:


Шаг 4. Если i=0, то перейти к шагу 5; иначе перейти к шагу 8.
Шаг 5. Вычислитьf
0(b ), p' (b ), p'' (b ).
Шаг 6. Положить

Шаг 7. Выбрать e0 є (0,l; 1).
Шаг 8. Вычислить:

Шаг 9. Если ||h(b )|| > ei, то перейти к шагу 10; иначе - положить
e
i+1=aei,
e' =
a'e' , e''=a''e'' , bi = b , i = i + 1 и перейти к шагу 2.
Примечание. Вычисление величины шага.
Шаг 10. Положить
l= 1.
Шаг 11. Вычислить
= f
0(b + l h(b )) + (1/e')p' (b +l h(b )) + e'' p'' (b + l h(b )) - f0(b ) -
-(1/e' )p' (b ) - e'' p'' (b )+ + 0.5|| h(b )||
2.
Шаг 12. Если < 0, то положить b= b +
lhb и перейти к шагу 8; иначе положить l=bl и перейти к шагу 11.

Весьма общий подход к ускорению сходимости штрафных функций состоит в использовании при решении последовательности задач сверхлинейно сходящегося алгоритма. Следующий алгоритм основан на методе сопряженных градиентов Полака-Рибьера с восстановлением, поскольку при исключении ограничений типа равенств с помощью штрафных функций могут не выполняться предположения о выпуклости. В случае, когда метод сопряженных градиентов восстанавливается через каждые n + 1 итерации, его сходимость к стационарной точке можно установить, не предполагая выпуклости целевой функции.

Алгоритм:

Примечание. Этот алгоритм решает задачу

min {f° (b ) | fj (b ) 0, i = 1, 2, ... m, rj(b } = 0,j = 0, 1, 2,...,l},

где функции f и r непрерывно дифференцируемы. Для него требуется предположение о том, что условия регулярности Куна-Таккера выполнены.

Шаг 0. Выбрать b0 є Rn , выбрать a, a' , a'' є (0; 0,5), выбрать b є (0,5; 0,8).
Шаг 1. Положить b = b
0 и положить i=0.
Шаг 2. Определить множества индексов

I' = {j є {l, 2,..., m} | fj(b ) > 0},
I'' = { j є {l, 2,..., m} | f
j(b ) < 0}.

и определить внешние и внутренние штрафные функции:


Шаг 3. Если i=0, то перейти к шагу 4; иначе перейти к шагу 7.
Шаг 4 Вычислитьf
0(b), p' (b ),p'' (b ).
Шаг 5 Положить

Шаг 6 Выбрать e0 є (0,l; 1).
Шаг 7. Определить(b ) =
f0(b ) + (1/e' )p' (b ) - e' p' (b ).
Шаг 8. Вычислитьf
0'(b ).
Шаг 9. Положить g = h = -f
0'(b ); положить j = 0.
Шаг 10. Положить q = 0.
Шаг 11. Определить q: R
1 (r) Rn равенством

q (x) = f0'(b + xh) - f0'(b ).

Шаг 12. Положить x = 0.
Шаг 13. Вычислить

q' (x) = <f0'(b + xh),b >.

Шаг 14. Если q'(x) = 0, то положить m = x и перейти к шагу 19; иначе перейти к шагу 15.
Шаг 15. Положить
l= 1.
Шаг 16. Вычислить  = q (x -
lq' (x)) - q (x) + 0.5lq' (x)2.
Шаг 17. Если < 0 то положить q = q + 1 и перейти к шагу 18; иначе положить
l=bl , и перейти к шагу 16.
Шаг18. Если q < 2, то положить x = x -
lq' (x) и перейти к шагу 13; иначе положить m= x - lq' (x) и перейти к шагу 19.
Шаг 19. Положитьb = b + m h.
Шаг 20. Вычислитьf
0(b ).
Шаг 21. Если ||f
0(b )II > ei, то перейти к шагу 22; иначе положить ei+1=aei, e' = a'e' , e''=a'e' , bi = b , i = i + 1 и перейти к шагу 2.
Шаг 22. Положить j = j + 1.
Шаг 23. Если j < n, то перейти к шагу 24; иначе перейти к шагу 9.
Шаг 24. Положить
g = - f
0(b ), 
Шаг 25. Положить g = g' , положить h = h' и перейти к шагу 10.
 

К комбинированным методам штрафных функций относится метод SUMT Фиакко и Маккормика

Методы линеаризации

Идея методов заключается в сведении задачи нелинейного программирования к задаче линейного программирования. С этой целью нелинейные функции целевой функции W(x) и ограничений g(x), h(x) раскладывают в ряд Тейлора до членов первого порядка в окрестности точки линеаризации xt, что позволяет W(x), g(x), h(x) аппроксимировать линейными функциями и свести общую задачу нелинейного программирования

к следующей задаче линейного программирования

Решая ее при помощи методов линейного программирования, находим новое приближение xt+1. В случае нелинейных функций точка xt+1 обычно недопустимая точка. Однако для сходимости к оптимуму достаточно, чтобы последовательность точек {xt}, полученных в результате решения последовательности подзадач линейного программирования, выполнялось следующее условие: значение целевой функции W и невязки по ограничениям в xt+1 меньше их значений в точке xt.

Описание метода

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

Рассмотрим задачу математического программирования вида:

f(x)  min,

gi(x)  0,i = 1,...,m, (1)

где функции f, g1,...,gm дифференцируемы, а функции g1,...,gm, кроме того, выпуклы на Rn.

Пусть Х - допустимое множество задачи (1). Будем считать, что Х - непусто.
Произвольной точке хR
n поставим в соответствие следующую задачу квадратичного программирования относительно h:

< f' (x),h > + 0.5||h||2  min,

< gi(x),h > + gi(x)  0, i = 1,...,m. (2)

Эта задача, как мы видим, сформирована на основе линейных частей разложений

f(x + h) = f(x) + < f' (x),h > + 0(||h||),

gi(x + h) = gi(x) + < gi'(x),h > + 0(||h||), i = 1,...,m,

причем в целевую функцию добавлен квадратичный член 0.5||h||2, а константа f(x) отброшена.

Заметим, что, в силу выпуклости функций g1,...,gm и условия Х , допустимое множество задачи (2) всегда не пусто. В самом деле, для любой точки Х, имеем:

< gi(x) , - x > + gi(x)  gi( ) 0, i = 1,...,m. (3)

Стало быть, вектор h =-x удовлетворяет ограничениям этой задачи.

Далее, благодаря квадратичному члену целевая функция задачи (2) сильно выпукла. Поэтому данная задача имеет решение, причем единственное.
Обозначим его через h(x).
В методе линеаризации последовательность приближений {x
k} генерируется по правилу

xk+1 = xk + akhk, k = 0,1,..., (4)

где hk - решение задачи (2) при x = xk,т.е. hk = h(xk), а ak выбирается тем или иным способом. При этом (в отличие от метода условного градиента) последовательность {xk} не обязана принадлежать множеству Х, а последовательность {f(xk)} - убывать.

Случай линейных ограничений

Задача НЛП с линейными ограничениями имеет следующий вид: 

f(x) min

при ограничениях Ах  b, при x 0.

Она отличается от задачи ЛП в стандартной форме нелинейностью целевой функции. Также как для задачи ЛП, допустимая область является многогранником. Поскольку функция f(x) нелинейная, оптимальное решение может не совпадать с вершиной или угловой точкой допустимой области. Кроме того, если f(x) не выпуклая, задача НЛП с линейными ограничениями может иметь несколько локальных минимумов. Рассмотрим с учетом этих особенностей последовательность действий при формулировке и решении линеаризованной в некоторой допустимой точке x0 задачи НЛП с линейными ограничениями.

Ясно, что задача вида

(x; x0) min

при ограничениях Ax b, x 0

представляет собой задачу ЛП и поэтому имеет оптимальное решение в допустимой угловой точке. Важен вопрос о соотношении между решением линеаризованной задачи * и решением исходной задачи x* .При решении линеаризованной задачи получается единственное решение, тогда как исходная задача может иметь несколько минимумов. Выбор точки локального минимума вблизи * зависит от выбора начальной точки линеаризации. Следовательно, при невыпуклой целевой функции нельзя быть уверенным в том, что найдено приближение к глобальному минимуму. Однако даже в случае выпуклой функции f(x) нет гарантии, что точка * хорошо "приближает" x* . Истинное решение x* может лежать внутри допустимой области, а точка * должна быть угловой. Из этого следует, что и в случае выпуклой функции f(x) для получения хорошего приближения x*необходима дополнительная корректировка * .

В силу допустимости точек х0,  * имеет место неравенство

(x0;x0) >(*; x0) .

Следовательно, если воспользоваться формулой для  , то получится следующее неравенство:

f(x0) > f(x0) + f(x0)(  * - x0),

или f(x0)(  * - x0) < 0.

Очевидно, что вектор (* - x0) задает направление спуска. Одномерный поиск из точки х0 в направлении ( * - x0) приводит к точке*. Поскольку * является угловой точкой допустимой области и x0 допустима, все точки на отрезке прямой между ними допустимы (так как допустимая область выпуклая). Более того, поскольку * представляет собой угловую точку, точки прямой, лежащие вне этого отрезка, недопустимы. Таким образом, одномерный поиск должен производиться на следующем отрезке прямой:

x = x0 + a (  * - x0), 0 a1.

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

В результате решения задачи

минимизировать f(x0 + a (  * - x0)), 0 a 1,

находится допустимая точка x1, такая, что f(x1) < f(x0). 

Поскольку величина f(x1) вообще не равна нулю, полученная точка х1 может служить точкой линеаризации для построения следующей аппроксимации. Решение последовательности задач ЛП и одномерный поиск продолжаются до тех пор, пока последовательные оптимумы хt, получаемые при помощи одномерного поиска, не станут достаточно близкими (т. е. расстояние между ними не станет меньше заранее выбранной малой величины). Ниже описан соответствующий алгоритм Франка - Вульфа .

Алгоритм Франка - Вульфа

Заданы х° и два числа e >0 и d>0.

Шаг 1. Вычислить f(xt). Если ||f(xt)|| e, то прекратить вычисления. В противном случае перейти к шагу 2.

Шаг 2. Решить следующую задачу ЛП:

минимизировать f(xt)y

при ограничениях Ау  b, у 0. 

Пусть уt - оптимальное решение этой задачи.

Шаг 3. Найти at, представляющее собой решение задачи

минимизировать f(х(t) +a(t) - x(t)))

при 0 al.

Шаг 4. Вычислить:

x(t+1) = х(t) +a(t)(t) - x(t))

Шаг 5. Проверить близость к решению. Если

|| x(t+1) - х(t)|| < d || x(t+1)||

и || f(x(t+1) ) - f(х(t))||  e ||f(x(t+1) )||,

то прекратить вычисления. В противном случае перейти к шагу 1.

В формулировке подзадачи ЛП опущены постоянные члены целевой функции (x, x(t)).

В случае оптимума, лежащего внутри допустимой области, алгоритм Франка - Вульфа весьма похож на градиентный метод оптимизации при отсутствии ограничений.
Для граничных точек работа алгоритма Франка-Вульфа сводится к проектированию градиента на грани допустимой области. При выполнении некоторых условий: непрерывной дифференцируемости f(x)
, ограниченности допустимой области, возможности точного решения подзадач ЛП и одномерного поиска легко может быть доказана сходимость алгоритма к точке Куна - Таккера при выборе любой допустимой точки в качестве начальной. Если функция f(x) выпуклая, решение представляет собой глобальный минимум. Поскольку определение направления спуска зависит от угловых точек, трудно дать оценку скорости сходимости алгоритма к оптимуму. Однако ясно, что скорость сходимости не может быть больше, чем у обычного градиентного метода.
Необходимо заметить, что при условии выпуклости целевой функции f(x) для алгоритмов Франка - Вульфа (а также любых алгоритмов, в которых используются подзадачи с линеаризованной целевой функцией) после решения каждой подзадачи можно получить очень удобную оценку разности между достигнутым и истинным значениями целевой функции. Для выпуклой функции f(x), линеаризованной в точке x
(t), выполняется неравенство

f(х)  f(х(t)) +  f(xt)(х-х(t)) =(x; х(t))
для всех х. Отсюда следует, что для всех допустимых точек

min f(х) min f (х; х(t)) = (y(t); х(t)).

Поэтому (y(t); х(t)) представляет собой нижнюю оценку значения f(x) в точке минимума. С другой стороны, f(y(t)), а лучше f(х(t+1))служит верхней оценкой минимального значения f(x). Таким образом, на каждой итерации разность

f(x(t+1)) -  (y(t); х(t))

может служить в качестве оценки улучшения целевой функции. Ее удобно использовать в качестве критерия останова на шаге 5 алгоритма Франка - Вульфа.

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

Общая задача нелинейного программирования

Рассмотрим общую задачу НЛП 

минимизировать f (х) 

при ограничениях gi(х) 0, j = 1, ..., J,

hk(x) = 0, k = 1, ..., K,

xi(U) xi xi(L), i = 1, ..., N.

При некоторой заданной оценке решения x(t) использование непосредственной линеаризации приводит к следующей задаче:

минимизировать f(х(t)) +  f(xt)(х-х(t))

при ограничениях gi(t)) +  gi(xt)(х-х(t)) 0, j = 1, ...,J,

hk(t)) + hk(xt)(х-х(t)), k = 1, ... K,

xi(U) xi xi(L), i = 1, ..., N.

Эта задача представляет собой задачу ЛП. Решая ее при помощи методов ЛП, находим новое приближение x(t+1). Заметим, что если xt - допустимая точка исходной задачи, то нет уверенности, что новая точка xt+1 будет допустимой. В действительности в случае нелинейных функций, входящих в ограничения, точка xt+1 обычно недопустима. Если xt+1 - недопустимая точка, то оптимальное значение линеаризованной целевой функции, удовлетворяющее неравенству

(x(t+1); х(t)) <  (x(t); х(t)),

может не быть точной оценкой истинного оптимума. Для сходимости к оптимуму достаточно, чтобы для последовательности точек {х(t)}, полученных в результате решения последовательности подзадач ЛП, выполнялось следующее условие: значения целевой функции и невязки по ограничениям в х(t+1) меньше их значений в точке х(t).

Правило изменения шага

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

-di  xi - x(t)i di, i = 1,..., N,

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

УПРАЖНЕНИЯ И ВОПРОСЫ К ТЕМЕ №4

1. Покажите, что функция f(x1, x2) = x12 + x22, где x1, x2=4 имеет минимум, равный 8 при (x1, x2) = ( +2, +2).

2. Покажите, что функция f(x, y) = x2 + y2 при ограничении x - y = 5 достигает минимума при x = 2,5, y = -2,5.

3. Если константы a, b, c, k - положительны, то найдите положительные x, y, z, такие, что функции k = x+y+z    и    w=ax2+by2+cz2 минимальны.

4. Если a, b, c - отрицательны и f(x1, x2, x3) = ax2x3 + bx3x1 + cx1x2, где x1 + x2 + x3 = 1, то покажите, что f(x1, x2, x3) имеет минимальное значение, равное abc/[2(ab + bc + ac) - (a2 + b2 + c2)] при условии положительности знаменателя.

5. Найдите стационарное значение (значения) функции f(x, y, z) = xy2z3 при ограничении x + y + z = 6.

6. Открытая коробка, изготовляемая из тонкого листа железа, имеет высоту z и прямоугольное основание с размерами x и y. Основание и стороны имеют толщину, равную d (малая величина), а стороны длиной y имеют толщину 2d. Если количество материала фиксировано, то покажите, что объем коробки максимален при x=2y=4z.

7. Решите задачу: минимизировать функцию f(x, y) = x2 + y2 при ограничениях x>=0, y>=0, x+y>=5.

8. Решите задачу: минимизировать функцию f(x, y) = x2 + 6xy - 4x - 2y при ограничениях x2-2y<=1, 2x-2y<=1.

9. Найдите минимальное и максимальное значения функции f(x,y)=x2+y2 при ограничении 3x2 + 4xy +6y2 = 140. Дайте геометрическую интерпретацию результата.

10. Посылка, которую должны отправить по почте, имеет форму прямоугольника с размерами x1, x2, x3. При отправлении посылки накладываются следующие ограничения: x1<=120, x2<=11, x3 <= 42   и   x1+x2+2x3<=72. Найдите размеры, при которых объем будет максимальным. (Это модифицированная Розенброком задача о почтовой рассылке).

11. Стандартная задача о почтовой рассылке аналогична приведенной выше задаче, за исключением того, что ограничения имеют вид xi<=42 при i=1, 2, 3, x1+2x2+2x3<=72.

12. Используйте комплексный метод для минимизации функции

f(х, у) = х22 

при ограничениях: x0, у0, х+у5.

13. Минимизируйте функцию

f(х, у) = х2+ 6ху - 4х - 2у

при ограничениях х2+2у1 2x-2у 1.

14. Минимизируйте функцию

f(x,у) = 3x12 +4х1х2 + 5х22 

при ограничениях х1>0, x2 0 и х1 + х,24

15. Проэкспериментируйте с программой, реализующей комплексный метод, заменив k - число точек в комплексе ; а - коэффициент отражения. Решите задачи, приведенные ниже и рассмотренные в тексте данной темы, с новыми значениями.

16. Решите задачу минимизации функции

f(x1, x2) = +(x1 + x2)

при ограничениях х10,    х20,    х12S. Рассмотрите два случая S=6 и S=4.

17. Минимизируйте функцию

f(х) =-[9-(х1- 3)223/27

при ограничениях x10,    0 x2 x1/,    0x1+x26.

18. Минимизируйте функцию

f(х) = х14 + х22 

при ограничениях х1 0,    х2 0,    х1*x28.

19. Минимизируйте функцию

f(х) = x12 + x2232,

при ограничениях х123 3,    х123 3.

(В качестве начальной точки возьмите точку (1; 2; 3))

20. Минимизируйте функцию

f(х1, х2,) = (х1-3)2+(x2-4)

при ограничениях 2*х12 + х22 34, 2*х1+З*х2 18, х1 0, x2 0.

21. Минимизируйте функцию

f(х1, х2, x3) = -х12*x3,

при ограничениях х1 0, х2 0, x3 0, 2*х1222+3*х32 51.

22. Используя метод SUMT, минимизируйте функцию f(x1, x2) = 3x12 + 4x1x2 + 5x22 при ограничениях x1 >=0, x2 >= 0, x1 + x2 >= 4.

23. Минимизируйте функцию f(x1, x2) = -x12 - x22 при ограничениях x1 >= 0, x2 >= 0 и x1 + 2x2 <= 3.

24. Используя метод SUMT, решите задачу о почтовых посылках: минимизируйте функцию V = -x1x2x3 при ограничениях 0 <= xi <= 42, где i=1, 2, 3 и x1 + 2x2 + 2x3 <= 72.

25. Минимизируйте функцию f(x) = x12 + x22 + x32 при ограничениях x1 + x2 + x3 >= 3, x1x2x3 >= 3, x1, x2, x3, >= 0.

26. Поэкспериментируйте с программой SUMT, изменив начальный выбор r , определение l, уменьшение r .

27. Поэкспериментируйте с программой метода SUMT, изменив критерий сходимости В последнем случае попробуйте условие if(R<IE-12).

28. Минимизируйте функцию (x1 - 1)4 + (x2 - 3)2, если x1, x2 >= 0 и 3x12 + 2x22 <= 21, 4x1 + 5x2 <= 20.

29. Минимизируйте функцию f(x1, x2) = x12 + x22 при ограничениях x1>= 2, x12 - x22 <= 1.

ПРИКЛАДНЫЕ ЗАДАЧИ К ТЕМЕ

Пример 4.1

Пример 4.2

Пример 4.3

Пример 4.4

Рассмотрим конкретный пример задачи условного нелинейного программирования.

Пример 4.1. План выпуска мебели.

Предприятие может выпускать два вида корпусной мебели. На их изготовление идет древесина трех видов. Запасы древесины на предприятии, нормы их расхода aij (i=1,2,3; j=1,2), себестоимость сj и оптовые цены указаны в табл.№1. Из-за брака в процессе производства расход древесины зависит от объема xj производства изделий и в первом приближении выражается линейной функцией aij + xj , а себестоимость продукции - функцией cj + 0,1xj . Изделия могут выпускаться в любых соотношениях, так как сбыт обеспечен. Предприятие обязано с целью изучения спроса населения выпустить не менее двух комплектов каждого вида мебели. Составить план выпуска изделий, обеспечивающий получение максимальной прибыли.


Таблица 1.
 

Порода

Запас сырья, м3

Нормы расхода на изделие вида

 

1

2

Сосна

100

10

20

Береза

120

20

20

Дуб

150

20

20

Себестоимость, млн. грн.

5

10

Цена, млн. грн.

7

13

Постановка задачи.

1. В качестве показателя эффективности целесообразно взять прибыль предприятия за операцию (в млн. руб.).

2. В качестве управляемых переменных задачи следует взять:

  •  x1 - количество комплектов корпусной мебели 1;
  •  x2 - количество комплектов корпусной мебели 2.

3. Целевая функция:

W = [7 - (5+0,1x1)]x1 + [13 - (10+0,1x2)] x2max,

или

W = 2x1 - 0,1x12 + 3x2 - 0,1x22max.

4. Ограничения:

4.1. По использованию сосны, м3:

(10+x1)x1 + (20+x2)x2100.

4.2. По использованию березы, м3:

(20+x1)x1 + (10+x2)x2 120.

4.3. По использованию дуба, м3:

(20+x1)x1 + (10+x2)x2 150.

4.4. Обязательства по контракту, шт.:

x1 2.

4.5. Областные ограничения:

x2 2.

Задача сводится к нахождению неотрицательных x1* и x2*, удовлетворяющих нелинейным ограничениям и доставляющих максимум нелинейной целевой функции.

Пример 4.2. Распределение топлива в силовых установках.

Каждый из двух электрических генераторов может работать на жидком топливе и газе или в любом их сочетании. Суммарная мощность, развиваемая обоими генераторами, должна составлять 50 МВт. Потребление газа ограничено 10 ед./ч. Желательно выбрать схему работы каждого генератора так, чтобы минимизировать потребление жидкого топлива.

Согласно кривой, отражающей эксплуатационные характеристики генератора 1, затраты топлива на обеспечение выходной мощности х1 МВт составляют

w1 = 1.4609 + 0.15186x1 + 0.00145x12,     (4.2.1)

w2 = 1.5742 + 0.1631x1 + 0.001358x12,    (4.2.2)

где w1 и w2 - затраты жидкого топлива и газа в час соответственно.

Аналогично для генератора 2 затраты топлива на обеспечение выходной мощности х2 МВт составляют

у1 = 0.8008 + 0.2031х2 + 0.000916х22,     (4.2.3)

у2 = 0.7266 + 0.2256х2 + 0.000778х22,    (4.2.4)

где у1 и у2 - затраты жидкого топлива и газа в час соответственно.

По эксплуатационным характеристикам генераторов х1 может изменяться в диапазоне

18 х1 30,    (4.2.5)

а х2 - в диапазоне

14 х2 25,    (4.2.6)

Предположим, что разные виды топлива могут комбинироваться аддитивным образом, т.е. для любого значения выходной мощности х1* любая линейная комбинация топливных компонент l1w1(x1*)+(1-l1)w2(x1*),

где 0l11, также дает мощность, равную х1*. Подобное предположение верно и для генератора 2.

Задача теперь формулируется следующим образом: требуется определить мощность каждого генератора х1 и х2 и доли использования разных видов топлива l1 и l2 так, чтобы минимизировать общее потребление жидкого топлива:

f = l1w1(x1)+ l2y1(x1)

при ограничении на использование газа

(1-l1)w2(x1)+(1-l2)y2(x2)10,    (4.2.7)

ограничении на суммарную мощность

х12=50     (4.2.8)

и ограничениях на переменные (4.2.5) и (4.2.6), а также

0l11 и 0l21.    (4.2.9)

Используя равенство (4.2.8), можно исключить переменную х2 и свести исходную задачу к задаче, содержащей три переменные, одно нелинейное неравенство (4.2.7), неравенства (4.2.9) и

25 х1 30;

последнее получается из (4.2.5), (4.2.6) и (4.2.8)

В работе [1] показано, что при начальном решении

z1 =20, l1=1/2, l2=1/2

и оценках границ начального интервала

z1=20, z2=1, z3=1

решение с точностью 0.1% достигается при реализации 55 серий по 100 пробных точек в каждой. Для получения точности в 0.01% необходимо 86 серий или 86000 пробных точек. Оптимальное решение имеет вид

f=3.05 ед.кол.нефти/час,

х1=30, х2=20,

l1=0 и l2=0.58.

Очевидно, что количество пробных точек 5500 итераций для решения задачи с целевой функцией и единственным нелинейным ограничением третьей степени слишком велико. Кроме того, требование дополнительных 3100 пробных точек для повышения точности на один знак неприемлемо.

Пример 4.3. 

Задача о нахождении размеров балки

Описание задачи.

Задана строительная конструкция, состоящая из балки А и жесткой опоры В. Балка А крепится на жесткой опоре В с помощью сварного соединения. Балка изготавливается из стали и должна выдержать нагрузку F = 2721,5 кг. Размеры балки, толщина, ширина сварного шва, а также марка стали выбираются таким образом, чтобы полные затраты были минимальными.

Рисунок 1 - Нагруженная балка.

Система состоит из балки А и сварного шва, необходимого для крепления балки к опоре В. Независимыми, или управляемыми, переменными служат размеры:

x=(x1, x2, x3, x4 )= (h, l, t, b )
х
1= h - ширина сварного шва;
x
2= l - глубина заделки балки;
x
3= t - высота балки;
x
4= b - ширина балки;
k - марка стали.

Математическая задача сводится к поиску минимума функции цели F(x)(1), представляющую полные затраты на реализацию строительной конструкции.

F(x)=(1+C3)X12X2+C4[K]X3X4(L+X2) min     (1)

Независимы, или управляемыми, переменными служат размеры

h- ширина сварного шва;
l- глубина заделки балки;
t- высота балки;
b- ширина балки.

Длина балки предполагается равной 35,56 см. Для удобства записи представим введенные переменные как компоненты неизвестного вектора Х, где X1=h;   X2=l;  X3=t;   X4=b;

X = [ X1, X2, X3, X4 ]= [ h, l, t, b ]

Характеристическим показателем качества проекта служат затраты на построение сварной группы. Основными стоимостными характеристиками такой группы являются С0 затраты на подготовительные работы, С1 затраты на сварочные работы и С2 стоимость материалов, т. е.

F(x)=C0+C1+C2,

где F(x)-функция затрат

Затраты на подготовительные работы C0.

Решение о проектировании рассматриваемой строительной конструкции в виде сварного узла было основано на наличии аппарата для сварки прямолинейным швом у заказчика проекта. Далее предположим, что все устройства, необходимые для установки и поддержания балки A в положении, удобном для проведения сварочных работ, также имеются в наличии. В случае затраты C0 можно не учитывать при построении модели.

C0 - затраты на подготовительные работы; C0=0;

Затраты на сварочные работы C1.

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

C1=(10 долл/4)*(1/60 * 4/мин)*(6 мин/дюйм)Vw = 1(долл/дюйм)*Vw,

где Vw - объем сварного шва в кубических дюймах.

Стоимость материалов С2.

С23Vw+C4[K]VB,

где С3=стоимость/объем сварного шва = (0,37*0,283)*163.87(долл./мм3);

С4=стоимость/объем балки = (0,17*0,283)*163,87(долл./мм3);

Vw= объем балки А

Vw=2(1/2h2l)=h2l, VB=tb(L+l) .


Имеем

С23*h2*l+C4[K]*t*b*(L+l).

Следовательно, функция затрат принимает следующий вид:

F(x)=h2*l+C3*h2*l+C4[K]*t*b*(L+l), (1.1)

или, (если выразить ее через переменные x)

x=[x1, x2, x3,x4]=[h, l, t, b], то:

F(x)=(1+C3)x12x2+C4[K]x3x4(L+x2) min    (1)

Заметим, что все комбинации значений x1,x2,x3,x4 могут оказаться допустимыми, если балка выдерживает заданную нагрузку. Минимум функции цели необходимо искать при условии, что выполняется ряд требований (ограничений). Указанные ограничения аналитически представлены следующими ограничениями:


g
1(x) = td[K] - t(x) >= 0;   (1.2)
g
2(x) = sd[K] - s(x) >= 0;   (1.3)
g
3(x) = x4 - x1 >= 0;   (1.4)
g
4(x) = x2 >= 0;   (1.5)
g
5(x) = x3 >= 0;   (1.6)
g
6(x) = Pc(x) - F >= 0;   (1.7)
g
7(x) = x1 - 0.125 >= 0;   (1.8)
g
8(x) = 0.25 - d(x) >= 0.   (1.9)

где

td- расчетное напряжение в сварном шве при сдвиге;
t(x)- максимальное напряжение в сварном шве при сдвиге, функция х;
sd- расчетное нормальное напряжение для материала балки;
s(x)-максимальное нормальное напряжение в балке, функция х ; Рс(х)- критическая нагрузка на балку, функция х;
d(х)- величина прогиба конца балки, функция х.
K=1,2,3...10 - номер марки стали.

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

Напряжение в сварном шве t(х).

Напряжение в сварном шве можно разложить на две составляющие t' и t",

где t'- первичное напряжение в плоскости поперечного сечения сварного шва,
а
t"- вторичное напряжение при кручении:

t'= и t"= .

При этом

M=F[L+(x2/2)],

R={

J=2{0.707x1x2[ ,

где
   M - момент силы F относительно центра тяжести сварной группы,
   J - полярный момент инерции сварной группы.


Напряжение
t в сварном шве вычисляется в соответствии с формулой:

t(x)=[(t')2+2t't" cos q+(t")2]1/2,

где cosq=x2/2R..

Напряжение при изгибе балки s(х).

Можно показать, что максимальное напряжение при изгибе равно

s(x)=6FL/x4x32 .

Критическая нагрузка на балку Pc(x).

С ростом отношения t/b=x3/x4 наблюдается тенденция к потере устойчивости балки. Те комбинации значений x3 и х4 при которых возможна потеря устойчивости, необходимо исключить из числа допустимых. Критическая нагрузка, для балок близких к прямоугольным, приближенно описывается следующим выражением:

Pc(x)= где

E- модуль Юнга = 30*106*64.52 *4,54(г/мм2),

I=(1/12)x3x43,

a=(1/3)Gx3x42,

Gk- модуль сдвига=12*106*64,52*4,54(г/мм2)

Прогиб балки d(х).

Чтобы вычислить величину прогиба, предположим, что балка представляет собой консоль с вылетом L. Тогда

d(x)=(4FL3)/(Ex33x4).

Остальные неравенства интерпретируются следующим образом:

g3 >=0 - устанавливает практическую невозможность получения сварного шва, ширина которого превышает ширину балки;

g4 >=0 и g5 >=0 отражает требования неотрицательности x2 и x3. Заметим, что неотрицательность x1 и x4 следует из неравенств g3 >=0 и g7 >=0.

Ограничение g6 >=0 гарантирует, что критическая нагрузка на балку не будет превышена.

Неравенство g7 >=0 отражает тот факт, что физически невозможно сделать сварной шов, ширина которого меньше некоторого порогового значения.

Наконец параметры td и sd, фигурирующие в ограничениях g1 и g2, зависят от материала конструкции.

Таким образом, надо найти xe1, xe2, xe3, xe4 и ke, при котором F(x) (функция цели) будет минимальна.

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

Оптимизационная задача проектирования включает функции затрат (1) и сложную систему неравенств, которая получается путем подстановки приведенных выше формул (1.2)-(1.9).

При этом все функции оказываются выраженными через четыре независимые переменные x1,x2,x3,x4, которые носят непрерывный характер, а k- переменная дискретная, то из этого следует, что данная задача нелинейного целочисленного программирования достаточно сложна и, как легко видеть, не может быть решена графическим способом. Данная задача решена методом простого случайного поиска.

Решение задачи

Выбираем случайным образом марку стали, а также неизвестные параметры x1,x2,x3,x4. Проверяем, является ли полученное решение допустимым, т.е. удовлетворяет ли оно перечисленным ограничениям. Если нет - генерируем новое случайное решение. Если да - находим значение функции цели F(x) и сравниваем его с предыдущим экстремальным Fe. Для начала его значение взято самым невыгодным (1038 долларов). Если найденное решение лучше предыдущего, мы экстремальному значению функции цели присваиваем только что найденное, соответственно запоминаем параметры, при которых оно получено xei=x; ke=k, где k- марка стали.

Так поступаем N раз, где N- число испытаний, после чего находим оптимальное решение xei, где i=1,2,3,4; ke и ему соответствующее значение функции цели Fei. В данной задаче функцией цели является функция полных затрат (1).

Логическое описание программы

Данная программа состоит из:

1. Ввода исходных данных и N - количества испытаний.

2. Создания 3 файлов для занесения в них исходных данных таких, как :

а) расчетное напряжение в сварном шве при сдвиге;
б) расчетное нормальное напряжение для материала балки ;
в) стоимости материала.

3.Цикла со счетчиком j=1,N для определения fei=fi т.е. определения минимальной функции затрат.

4. Оператора цикла for со счетчиком i=0,3 для генерирования случайного числа sl, определения первичного и вторичного напряжения при кручении t, R радиуса при кручении, I полярного момента инерции, M момента силы F относительно центра тяжести, cosq и

5. Проверки ограничений.

6. Функции random поиска псевдослучайного числа sl .

7. Оператора if, для поиска экстремального значения функции цели (fei) - минимальной функции затрат. Если fi > fei, продолжаем испытания , а если fi < fei, тогда fei = fi и xei=xi.

8. Создания файла Margo.dan и вывода результатов fei и xei, ke, где i=1, 2, 3, 4 (параметры балки), а ke - марка стали, в файл.

9. По данным результатов строим график зависимости функции цели от количества испытаний.

ПРОГРАММНАЯ РЕАЛИЗАЦИЯ

Данные алгоритма представлены в виде программы на языке СИ. Программа состоит из исходного файла и файла запуска программы.

Исходный файл составляет 6305 Кбт,
Файл запуска программы 32351 Кбт.

Программа работает на машинах IBM. Для работы подходит любая машина, даже 286, что для практического применения очень выгодно.

Предварительно, для работы программы необходимо создать 3 файла, куда заносятся данные:

1) расчетное напряжение в сварном шве при сдвиге;

2) расчетное нормальное напряжение для материала балки;

3) стоимость материала.

Далее на языке СИ программа решает задачу оптимизации методом простого случайного поиска, в которой определяет функцию цели - минимальную функцию затрат. Результаты заносим в таблицу "Соотношение функции полных min затрат от количества испытаний". Затем строим график зависимости функции цели от числа испытаний.

Таблица - Соотношение функции полных min затрат от количества испытаний

Кол-во испытаний N

Функция цели fi

100

49,537

200

35,035