Методы оптимизации

МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ ПОЛИТЕХНИЧЕСКИЙ ИНСТИТУТ (ФИЛИАЛ) ФЕДЕРАЛЬНОГО ГОСУДАРСТВЕННОГО БЮДЖЕТНОГО ОБРАЗОВАТЕЛЬНОГО УЧРЕЖДЕНИЯ ВЫСШЕГО ОБРАЗОВАНИЯ «ДОНСКОЙ ГОСУДАРСТВЕННЫЙ ТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ» В Г. ТАГАНРОГЕ РОСТОВСКОЙ ОБЛАСТИ ПИ (филиал) ДГТУ в г. Таганроге Автомобилестроение и сервис транспортных средств МЕТОДЫ ОПТИМИЗАЦИИ Конспект лекций для студентов – заочников направления 09.03.02 «Информационные технологии» по курсу «Методы оптимизации» (бакалавриат) 2016 г. 1 СОДЕРЖАНИЕ ВВЕДЕНИЕ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1. ПРИМЕРЫ ЗАДАЧ ОПТИМИЗАЦИИ . . . . . . . . . . . . . . . . . . . . . . . 1.1. Задача о распределении ресурсов . . . . . . . . . . . . . . . . . . . . . . . 1.2. Задача о перевозках . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3. Задача о загрузке . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4. Задача о раскрое материала . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5. Задача об использовании мощностей. (Загрузка предприятия) 1.6. Выбор плана лечения больного пациента . . . . . . . . . . . . . . . . 1.7. Метод наименьших квадратов для обработки результатов экспериментов . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2. ЛИНЕЙНЫЕ ЗАДАЧИ ОПТИМИЗАЦИИ . . . . . . . . . . . . . . . . . . . . 2.1. Постановка и анализ задач ЛО . . . . . . . . . . . . . . . . . . . . . . . . . 2.2. Графическое решение задач ЛО. . . . . . . . . . . . . . . . . . . . . . . . . 3. ВЫПУКЛЫЕ МНОГОГРАННИКИ И ЛО. . . . . . . . . . . . . . . . . . . . . . 3.1. Выпуклые множества и свойства допустимой области. . . . . . 3.2. Свойства допустимых решений ОЗЛО. . . . . . . . . . . . . . . . . . . . 4. СИМПЛЕКС–МЕТОД РЕШЕНИЯ ЗАДАЧ ЛО. . . . . . . . . . . . . . . . . 4.1. Оптимальный допустимый базис и оптимальное решение. . . . 4.2. Табличный симплекс–метод. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3. Симплекс–алгоритм. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4. Начальный допустимый базис. . . . . . . . . . . . . . . . . . . . . . . . . . . . 5. ДВОЙСТВЕННОСТЬ В ЛИНЕЙНОЙ ОПТИМИЗАЦИИ . . . . . . . . 5.1. Двойственные задачи . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6. ТРАНСПОРТНАЯ ЗАДАЧА . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1. Постановка задачи, свойства и нахождение опорного решения 6.2. Распределительный метод решения ТЗ . . . . . . . . . . . . . . . . . . . . 6.3. Метод потенциалов решения транспортной задачи . . . . . . . . . . 6.4. Двойственная транспортная задача . . . . . . . . . . . . . . . . . . . . . . . 6.5. Особые случаи решения ТЗ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7. НЕЛИНЕЙНАЯ ОПТИМИЗАЦИЯ . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.1. Безусловная оптимизация . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2. Выпуклые функции и их свойства . . . . . . . . . . . . . . . . . . . . . . . . 7.3. Условная оптимизация при ограничениях–равенствах . . . . . . . Стр. 3 3 3 4 6 6 7 8 9 10 10 13 18 18 21 27 27 33 36 38 44 44 50 50 56 63 70 72 76 76 80 85 2 7.4. Условная оптимизация при ограничениях–неравенствах . . . . . 8. ЧИСЛЕННЫЕ МЕТОДЫ БЕЗУСЛОВНОЙ ОПТИМИЗАЦИИ . . . . 8.1. Общая характеристика численных методов оптимизации . . . . 8.2. Методы нулевого порядка одномерной оптимизации . . . . . . . . 8.2.1. Метод деления интервала пополам . . . . . . . . . . . . . . . . . . . . . . 8.2.2. Метод дихотомии . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.2.3. Метод золотого сечения . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.3. Методы нулевого порядка многомерной оптимизации . . . . . . 8.3.1. Метод конфигураций . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.3.2. Метод деформируемого многогранника . . . . . . . . . . . . . . . . . 8.3.3. Метод сопряженных направлений . . . . . . . . . . . . . . . . . . . . . . 8.3.4. Адаптивный метод случайного поиска . . . . . . . . . . . . . . . . . . 8.4. Методы первого порядка . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.4.1. Метод градиентного спуска с постоянным шагом . . . . . . . . . 8.4.2. Метод наискорейшего градиентного спуска . . . . . . . . . . . . . . 8.4.3. Метод Гаусса–Зейделя . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.4.4. Метод Дэвидона–Флетчера–Пауэлла . . . . . . . . . . . . . . . . . . . . 8.5. Методы второго порядка . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.5.1. Метод Ньютона . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.5.2. Метод Ньютона–Рафсона . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ЛИТЕРАТУРА . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 93 93 98 100 102 103 105 105 107 110 111 113 113 115 117 120 122 122 124 127 3 ВВЕДЕНИЕ Содержание курса ―Методы оптимизации‖ составляют методы нахождения экстремумов функций нескольких переменных на множествах, которые задаются линейными или нелинейными равенствами и неравенствами. Другое распространѐнное название курса ―Математическое программирование‖, исторически более раннее, объясняется тем, что первые математические исследования и их приложения развивались совместно с задачами экономики и планирования, т.е. созданием в некотором, смысле оптимальной программы–плана. Сначала для такого рода задач с линейными функциями Дж. Данцигом (1949 г.) и был предложен термин ―линейное программирование‖. Затем для другого типа задач с нелинейными функциями Х. Кун и А. Таккер (1951) предложили использовать название ―нелинейное программирование‖. Впоследствие, в зависимости от используемых в задачах классов функций, появились термины ―целочисленное программирование‖, ―выпуклое программирование‖, ―динамическое программирование‖ и в 1959 году появился термин ―математическое программирование‖. В дальнейшем, в 70-х годах ХХ столетия, чтобы отличать эту дисциплину от программирования для вычислительных машин, чаще стал использоваться термин ―методы оптимизации‖. Методы оптимизации являются одним из наиболее развивающихся разделов прикладной математики в силу большого многообразия и важности их приложений. В качестве примеров можно указать следующие области применения методов оптимизации: технико-экономические системы: планирование, управление, транспортные задачи, распределение кадров, ресурсов; численный анализ: аппроксимация, регрессия, решение систем уравнений, численные методы решения задач математической физики; технические системы: распознавание образов, оптимальное управление, роботы, управление производством, проектирование систем связи; математическая экономика: анализ больших макроэкономических моделей, микроэкономических моделей или моделей предпринимательства, теория принятия решений, теория игр. 1. ПРИМЕРЫ ЗАДАЧ ОПТИМИЗАЦИИ 1.1. Задача о распределении ресурсов Многие практические задачи начинаются оптимизацией распределения ресурсов. Пусть необходимо производить товары T1 ,...,Tn T j , j 1, n . Для этого имеются ресурсы R1 ,...,Rm (эл. энергия, оборудование, рабочая сила и т.д.) в количествах b1,...,bm bi , i 1, m единиц. Стоимость единицы ресурса 4 bi равна d i д.е. Для производства одной единицы товара T j надо aij единиц ресурса Ri . Каждая единица товара T j может быть реализована по цене c j . Известно, что рынок не может поглотить более, чем k j единиц товара T j . Необходимо составить план производства товара T j , таким образом, чтобы, хватило ресурсов, а полученная прибыль была максимальной. Составим математическую модель. Обозначим x j , j 1, n количество товара j-го типа, которое надо произвести. Тогда, по условию xj 0, xj kj , причѐм ресурсов должно хватать, т.е. (1.1) n aij x j bi , (i 1, 2,..., m). (1.2) j 1 Себестоимость единицы товара T j равна sj a1 j d1 a2 j d 2 ... amj d m m aij d i . i 1 Чистая прибыль от одной единицы товара T j есть q j прибыль будет равна: cj s j , а общая n z q1x1 ...qn xn q j x j. (1.3) j 1 Задача состоит в следующем: необходимо выбрать такие значения переменных, x1, x2 ,...,xn которые удовлетворяют неравенствам (1.1), (1.2) и обращают в максимум функцию n переменных z из (1.3). 1.2. Задача о перевозках Имеется m складов C1, C2 , C3 ,...,Cm и n пунктов потребления П1 , П 2 , ..., П n . На складах имеются запасы товаров одного типа в количествах a1, a2 ,...,am единиц. Пункты потребления подали заявки соответственно на b1 , b2 ,...,bn , единиц товара. Заявки выполнимы, т.е. сумма всех заявок не превосходит суммы всех имеющихся запасов n bj j 1 m i 1 ai . Стоимость перевозки одной единицы товара со склада Ci в пункт П j равна cij , i 1, m; j 1, n . Требуется составить план перевозок таким образом, чтобы все заявки были выполнены, а общие расходы на перевозки были бы минимальными. 5 Составим математическую модель. Обозначим через xij количество товара, отправленного со склада Ci в пункт П j . Решение (план перевозок) будет состоять из m n чисел, т.е. является матрицей: x11 x12  x1n x 21 x 22  x 2 n .     x m1 x m 2  x mn Необходимо выбрать такие неотрицательные значения переменных, xij чтобы были выполнены следующие условия: 1. Ёмкость склада не должна быть превышена X n xij a i , ( i 1, m ). (1.4) j 1 2. Поданные заявки должны быть удовлетворены, т.е. m i 1 xij b j , ( j 1, n ). (1.5) Общая стоимость перевозок будет равна: z c11 x11 c12 x12 ... c1n x1n c 21 x 21 c 22 x 22 ... c 2n x 2n ............................................... c m1 x m1 c m 2 x m 2 ... c mn x mn или m n z cij xij . (1.6) i 1j 1 Необходимо составить план перевозок, т.е. найти матрицу xij такую, чтобы переменные xij были неотрицательными и удовлетворяли условиям (1.4) и (1.5), а при этом функция z принимала минимальное значение. Если сумма всех заявок равна сумме запасов n bj j 1 n i 1 ai , то из каждого склада будут вывезены все товары и неравенства (1.4) обращаются в равенства. Задача (1.4)-(1.6) называется транспортной задачей: 6 m n z cij xij min , i 1j 1 xij m i 1 n 0, xij bj, xij ai . j 1 1.3. Задача о загрузке (задача о ранце, рюкзаке …) Пусть некоторое транспортное средство загружается предметами n типов. Каждый предмет j-го типа j 1, 2,..., n имеет вес (массу) m j кг. Максимальная грузоподъѐмность транспорта M кг. Каждый предмет j-го типа имеет определѐнную ценность c j (стоимость единицы, калорийность единицы и т.д.). Сколько предметов каждого типа следует погрузить, чтобы суммарная ценность z была максимальной, а вес не превышал допустимую грузоподъѐмность М. Составим математическую модель. Обозначим через x j количество предметов j-го типа, которое погружается. Суммарная ценность груза будет равна z c1 x1 c2 x2 ... cn xn n cjxj . (1.7) j 1 Ограничение по массе задается неравенством n mjxj j 1 m1 x1 m2 x 2 ...mn x n M. (1.8) Все переменные должны быть целыми неотрицательными x j 0, x j N , j 1, n . Тогда задача состоит в следующем. Необходимо определить целые неотрицательные переменные x1, x2 ,..., xn которые, доставляют максимум функции (1.7) при ограничении-неравенстве (1.8). 1.4. Задача о раскрое материала На раскрой (распил, обработку) имеется материал одного образца в количестве a единиц. Требуется изготовить из него m разных комплектующих изделий в количествах пропорциональных числам b1, b2 ,...,bm (условие комплектности). Каждая единица материала может быть раскроена n способами, 7 причѐм использование j-го способа j 1, n дает aij единиц i-го изделия i 1, m . Необходимо найти план раскроя, обеспечивающий максимальное число комплектов. Составим математическую модель. Обозначим через x j число единиц материала, раскраиваемых j-м способом, а x - число изготавливаемых комплектов изделий. Так как общее количество материала равно сумме его единиц, раскраиваемых различными способами, то n xj a. (1.9) j 1 Требование комплектности выразится равенством n aij x j bi x, i 1, m . (1.10) j 1 Очевидно, что x j 0, i 1, n . Требуется найти неотрицательные значения переменных x1, x2 ,...,xn , удовлетворяющие системе уравнений (1.9) при которых функция z x из (1.10) принимала максимальное значение. 1.5. Задача об использовании мощностей. (Загрузка предприятия) Предприятию необходимо за время T выпустить k1, k 2 ,...,k n единиц продукции P1, P2 ,...,Pn . Продукция производится на станках S1 , S 2 ,...,S m . Для каждого станка известна производительность aij (т.е. число единиц продукции Pj , которое можно произвести на станке S i ) и затраты bij на изготовление продукции Pj на станке S i в единицу времени. Необходимо составить такой план работы станков (т.е. надо распределить выпуск продукции между станками), чтобы затраты на производство всей продукции были минимальными. Составим математическую модель. Обозначим через xij время, в течении которого станок Si , i 1, m будет занят изготовлением продукции Pj , j 1, n . Так как время работы каждого станка ограниченно и не превышает T , то справедливы неравенства: x11 x12 ... x1n T , x21 x22 ... x2 n T , .............................., xm1 xm 2 ... xmn T . (1.11) 8 Для выполнения плана выпуска по номенклатуре необходимо, чтобы выполнялись неравенства: a11x11 a21x21 ... xm1 k1, a12 x12 a22 x22 ... xm 2 k2 , (1.12) .........................................., a1n x1n a2 n x2 n ... xmn kn . Затраты на производство всей продукции выразится функцией (1.13) z b11x11 b12 x12 ... bmn xmn . Необходимо найти неотрицательные значения переменных xij , i 1, m; j 1, n , удовлетворяющие системам (1.11) и (1.12) при которых функция (1.13) принимала бы минимальное значение. 1.6. Выбор плана лечения больного пациента Пусть комплекс заболеваний больного пациента однозначно определѐн. Известны методы лечения M1, M 2 ,...,M n , которые представляют собой различные лекарственные препараты и их дозы, параметры схем применения лекарственных средств, физиотерапевтические процедуры и их продолжительность, облучение и т.п. Состояние больного описывается набором показателей f1, f 2 ,..., f m (температура тела, давление крови, частота дыхания, параметры различных анализов), которые зависят от параметров метода лечения M j , j 1, n . Пусть допустимый диапазон f i изменяется для здорового человека ai ; ai . Аналогичным способом должны быть установлены зависимости токсичности, показателей побочных эффектов и других вредных воздействий q1, q2 ,...,q p на больного от параметров метода лечения M j , j 1, n . Предположим, что допустимое значение нежелательного фактора k-го типа k 1, p ограничено сверху величиной bk . Кроме того, из практического опыта известен диапазон изменения параметров метода лечения M j . С каждым набором лекарств и процедур связывается некоторый критерий качества лечения z , который зависит от параметров выбранного метода лечения, характеризующий стоимость и продолжительность лечения. Это может быть, например, совместная вероятность того, что различного рода затраты и время лечения не превысят заданных допустимых значений. Составим математическую модель такой детерминированной задачи выбора метода лечения при фиксированном диагнозе. Обозначим через x1, x2 ,...,xn параметры методов лечения M1, M 2 ,...,M n . Тогда объективные 9 показатели fi , i 1, m состояния больного будут функциями x j , j 1, n , т.е. fi fi x1, x2 ,..., xn , i 1, m и удовлетворяют неравенствам ai fi x1, x2 ,...xn ai , i 1, m . (1.14) Токсичность и побочные эффекты qk , k 1, p также будут функциями x j и будут удовлетворять неравенствам: qk x1, x2 ,..., xn bk , k 1, p . (1.15) Кроме того, переменные x j должны лежать в известном заранее диапазоне изменения: (1.16) xj xj xj , а критерий качества (1.17) z z x1, x2 ,..., xn , должен быть максимальным. Таким образом, необходимо выбрать такие параметры x1, x2 ,...,xn методов лечения M1, M 2 ,...,M n чтобы выполнялись системы неравенств (1.14)(1.16), а функция качества z принимала максимальное значение. 1.7. Метод наименьших квадратов для обработки результатов экспериментов Пусть проведено n опытов, при которых изменялась величина y в зависимости от изменений x . Получилась таблица наблюдений: x x1 x2 … xn y где ~ yj ~ y1 ~ y2 … ~ yn значение измереной величины. Необходимо построить многочлен степени m y f ( x) m a k x k , ко- k 0 торый бы аппроксимировал результаты наблюдений, так чтобы значения y j при y f ( x) в некотором смысле были близки к результатам наблюдений ~ значениях xj. Согласно методу наименьших квадратов надо найти минимум следующей функции n переменных a0 , a1, a2 ,..., am : z a0 , a1,...,an n j 1 ~ yj m ak x kj k 0 2 . 10 Все приведѐнные примеры позволяют сформулировать задачу в общем виде: найти совокупность переменных x1, x2 ,...,xn , которые удовлетворяют некоторой системе неравенств (равенств) и доставляют min или max функции n переменных z f x1, x2 ,...,xn . z f ( x1, x2 ,..., xn ) min(max), аi xj bj , j 1, n , gi ( x1, x2 ,..., xn ) 0, i 1, m. Такая задача называется задачей оптимизации. Функция z f ( x1, x2 ,...,xn ) называется целевой функцией, а все неравенства и равенства, которым должны удовлетворять неизвестные называются ограничениями, накладываемыми на переменные. В зависимости от класса функций f ( x1, x2 ,..., xn ) и gi ( x1, x2 ,..., xn ) для решения подобных задач используются различные методы и по-разному называются соответствующие разделы курса. 2. ЛИНЕЙНЫЕ ЗАДАЧИ ОПТИМИЗАЦИИ Изучение методов оптимизации начнѐм с линейных задач. Часто методы линейной оптимизации (ЛО) служат основой для решения других, более сложных задач оптимизации. 2.1. Постановка и анализ задач ЛО В общем виде задача ЛО имеет следующую формулировку: найти значения n переменных x1, x2 ,...,xn , которые обращают в min (max) линейную функцию z c1x1 c2 x2 ... cn xn n cjxj , j 1 сj R . (2.1) и удовлетворяют системе линейных неравенств a11x1 a12 x2 ... a1n xn b1, ................................................, am1x1 am 2 x2 ... amn xn bm , (2.2) где ai , b j R , причѐм x j 0 , ( j 1, n ). (2.3) В дальнейшем, для определѐнности, всегда будем решать задачу минимизации z , т.к. если по содержанию задачи надо найти max z , то такая задача сводится к нахождению min ( z ). 11 Определение 1. Допустимым решением задачи ЛО (2.1)-(2.3) называется любая совокупность переменных x j удовлетворяющих условиям (2.2), (2.3). Определение 2. Оптимальным решением задачи ЛО называется такое решение из множества допустимых решений, при котором линейная функция z обращается в минимум. Определение 3. Множество точек пространства Rn, координаты которых удовлетворяют ограничениям (2.2), (2.3), называется областью допустимых решений (ОДР). Система ограничений (2.2) не обязательно полностью должна состоять из неравенств, т.е. она может быть смешанной и состоять как из неравенств, так и равенств. Если все ограничения задачи ЛО заданы в виде равенств, то такая форма ограничений называется канонической, а сама задача называется основной задачей линейной оптимизации (ОЗЛО). Запишем ОЗЛО в матричной (векторной) форме. Обозначим a11 a12 .. a1n c1 b1 x1 a 21 a 22 .. a 2n c2 b2 x2 ,c ,b , x . .. .. .. .. .. .. .. a m1 a m 2 .. a mn cn bm xn Если все ограничения являются равенствами, то можно записать следующую задачу ОЗЛО А z cT x min, ОЗЛО : Ax b, x 0. (2.4) От общей задачи ЛО (2.1)-(2.3) всегда можно перейти к ОЗЛО. Пусть все ограничения (2.2) заданы в виде неравенств n aij x j j 1 bi , i 1, m . Введѐм новые переменные xn 1, xn 2 ,..., xn xn 1 b1 (a11x1 ... a1n xn ), xn 2 b2 (a21x1 ... a2n xn ), ..............................................., xn m bm (am1x1 ... amn xn ). m следующим образом (2.5) Новые, ‖добавочные‖, переменные xn i , i 1, m , согласно ограничениям должны быть неотрицательными. Эти переменные также называют ―слабыми‖ переменными. 12 Теперь получили ОЗЛО: найти такие неотрицательные переменные x1, x2 ,...,xn xn 1, xn 2 ,...,xn m , которые бы обращали в min линейную функn цию z c j x j и удовлетворяют системе уравнений (2.5). j 1 Пример 2.1. Привести следующую задачу ЛО к ОЗЛО. z x1 2 x2 3x3 min 2 x1 x2 3 x3 6, x3 3 x2 1, x5 2 x4 x1 x5 x1 0, 1, x1, x2 , x3 , x4 , x5 0. Решение. Приведѐм неравенства к виду (2.2) введя добавочные переменные x6 , x7 , x8 , x9 2 x1 x2 3x3 6, x6 x3 3x2 1, x7 x5 1, x8 1 x1 2 x4 x5 6 2 x1 x2 3x3 , 1 3 x2 x3 , x1 2 x4 x5 , (2.6) x1 0, x9 0 x1 x5 . Таким образом, необходимо найти неотрицательные значения переменных x1 , x2 , x3 , x4 , x5 , x6 , x7 , x8 , x9 , которые удовлетворяют системе уравнений (2.6) и обращают в min функцию z . Итак, в дальнейшем будем анализировать и искать решения только ОЗЛО с ограничениями в виде равенств, т.к. общую задачу с ограниченияминеравенствами всегда можно свести к ОЗЛО. Рассмотрим каноническую систему ограничений - уравнений a11x1 a12 x2 ... a1n xn b1, n aij x j b j , (i 1, m). ..............................................., (2.7) j 1 am1x1 am 2 x2 ... amn xn bn , Для нахождения решения ОЗЛО необходимо, чтобы она имела решение, т.е. совместна. Из линейной алгебры известно (теорема КронекераКапелли), что для совместности системы линейных уравнений необходимо и достаточно, чтобы ранг матрицы системы A был равен рангу расширенной матрицы A r r , где a11 .. a1n b1 A' .. .. .. .. . am1 .. amn bm 13 Ранг матрицы есть максимальное число линейно независимых строк или столбцов матрицы, поэтому ясно, что ранг системы (2.7) не больше, чем число уравнений или переменных, т.е. r m или r n . Возможны два случая: m n и m n . 1. Пусть m n и r n . Отбросим в системе ―лишние‖ уравнения ( m n ), которые являются линейной комбинацией остальных и система примет вид: n aij x j j 1 bi , i 1, n . Поскольку r r n , то система имеет единственное решение x ( x1, x2 ,..., xn ) . Если в этом решении хотя бы одно x j 0 , то ОЗЛО решений не имеет. Если j : x j 0 , то найденное решение является допустимым. Оно же является и оптимальным. Отметим, что если здесь r n , то этот случай сводится к б). 2. Если m n , тогда очевидно, что r n , т.е. число ограничений меньше чем число переменных. В этом случае, если система совместна, то у неѐ существует бесчисленное множество решений. Если ( n m ) переменным придавать произвольные значения (свободные переменные), то остальные переменных будут через них выражаться (базисные переменные). Пример 2.2. Выделить свободные и базисные переменные. 2 x1 x2 x3 x4 1, 2 1 1 1 , r 2. 1 1 2 1 x1 x2 2 x3 x4 2, T Решение. Пусть х1, х2 – базисные; х3, х4 - свободные переменные, тогда: 2 x1 x2 1 x3 x4 , x1 3 x3 , x1 x2 2 2 x3 x4 , x2 5 3x3 x4 . В дальнейшем, записывая ограничения в канонической форме, будем считать их линейно-независимыми т.е., r m n , n m 0 . 2.2. Графическое решение задач ЛО Предполагая, что n m , положим в начале n m 1 , r m . В этом случае свободных переменных одна, поэтому целевую функцию z и остальные ( n 1 ) переменные можно выразить через неѐ. Пусть эта переменная будет x1 . Тогда z f ( x1 ) 1 x1 0 min(max) , x2 2 x1 2 0, x3 3 x1 3 0, ............................, xn n x1 n 0. 14 Если система ограничений не противоречива, то она определяет отрезок на оси Ox1, а z min или z max достигаются на концах отрезка. Теперь пусть n m 2 . В этом случае две переменные, например x1 и x2 можно взять в качестве свободных переменных и выразить через них остальные x3 , x4 ,...,xn (это можно сделать методом элементарных преобразований), которые будут базисными. x3 31x1 32 x2 3 0, x4 41x1 42 x2 4 0, (2.8) ........................................., xn n1x1 n 2 x2 n 0. На плоскости Ox1 x2 покажем допустимую область изменения переменных x1 и x2 . Это первая четверть, т.к. x1 0 и x2 0 . Остальные переменные в силу (2.8) тоже неотрицательные. Пусть 31 x1 32 x 2 3 0. Это прямая плоскость Ox1 x2 делит на две полуплоскости, причѐм в одной полуплоскости x3 0 , а в другой x3 0 . Отметим штрихами область где x1, x2 , x3 0 . Аналогичным образом построим и все остальные прямые x4 0, x5 0, x6 0,...,xn 0 : 0; 41 х1 42 х2 4 ……………………. n1 х1 n 2 х2 n 0. Получим n прямых: две оси координат x1 0 и x2 0 , а также n-2 прямые x3 0 ,…, xn 0 . Каждое неравенство х j 0 определяет полуплоскость. Часть плоскости Ox1 x2 принадлежащая одновременно всем полуплоскостям х j 0 образует область допустимых решений ОЗЛО. Возможны три случая: ОДР не существует. Это будет тогда, когда система уравнений-ограничений несовместна (противоречива), например, как показано на рис.2.1. В этом случае нет общей области пересечения прямых и ОЗЛО не имеет допустимого решения. Оптимального решения тем более нет. ОДР ограничена и существует. Пересечение полуплоскостей x j 0, j 1, n представляет собой замкнутый многоугольник, например АBCDE (рис. 2.2.). Это и будет ОДР ОЗЛО, которая есть ограниченное замкнутое множество в R2. ОДР существует, но неограниченна. Общей областью пересечения полуплоскостей х j 0 является неограниченная область на Ox1 x2 . Здесь ОДР существует, но неограниченна (рис. 2.3). 15 x2 x2 В С А ОДР x1 E Рис. 2.1 D x1 Рис. 2.2 x2 ОДР x1 Рис. 2.3 Ясно, что ОЗЛО может иметь решение только в случаях 2 и 3. Поскольку целевая функция z n c j x j есть функция n переменных, то по j 1 второй теореме Вейерштрасса функция n переменных f x1, x2 ,...,xn , непрерывная в ограниченной замкнутой области достигает в ней минимального и максимального значения, которые могут быть достигнуты либо в одной из стационарных точек, либо на границе области. Стационарные точки находятz z c j . Поскольку, c j од0 , j 1, n . В нашем случае ся из условия xj xj новременно не равны 0, то стационарных точек нет. Следовательно, минимум находится на границе ОДР. Если х1, х2 – свободные переменные, то используя (2.8) выразим z через x1 и x2 : z 1 x1 2 х2 0. 16 Обозначим z ' 1 x1 2 х2 ; z ' z 0 , где 0 не зависит от x1 и x2 . Очевидно, что z и z достигают минимального и максимального значений при одних и тех же значениях x1 и x2 . Запишем z ' 1x1 2 х2 С и получим уравнение прямой на Ox1 x2 с нормальным вектором n ( 1, 2 ) grad z z. Различным значениям z соответствуют различные прямые на плоскости Ox1 x2 , и все они параллельны между собой. Поэтому достаточно изобразить на плоскости одну прямую (основную) z 0 , проходящую через начало координат. Перемещая еѐ параллельно самой себе, получаем различные значения z . Если перемещать в направлении нормали n z , то значения z , а следовательно и z возрастают, а в противоположном направлении убывают. Двигаясь по ОДР можно определить, имеет ли ОЗЛО единственное решение, или его нет, или решений бесконечно много. Например, из рис. 2.4. следует, что z min z B , zmax z A . X2 В ОДР А n z n z z C2 X1 C1 Рис. 2.4 Пример 2.3. Найти оптимальное решение задачи. z 16 х1 х 2 х3 5х 4 5х5 min , 2 x1 x2 x3 10, 2 x1 3x2 2 x1 4 x2 x4 6, x5 8, x1, x2 , x3 , x4 , x5 0. Решение. Здесь m 3 , n 5 , следовательно, три переменные базисные, две переменные свободные. Пусть x1 и x2 – свободные переменные. Выразим переменные x3 , x4 , x5 и целевую функцию z через переменные x1 , x2 , причѐм x3 0, x4 0, x5 0 . 17 х3 10 2 x1 x2 0, 2 x1 x4 0, 2 x1 3 x2 6 2 x1 3 x2 x2 10, 6, х5 8 2 x1 4 x2 0, 2 x1 4 x2 8. Построим на плоскости Ox1 x2 прямые x3 0, x4 0, x5 0 : 2 х1 3х2 6 – II 2 х1 х 2 10 – I 2 х1 4 х 2 8 – III и находим ОРД из условий x j 0, j 1,2,3,4,5 . x2 I II A ОДР Z Z III x1 n ( 2, 3) Рис. 2.5 Подставим x3 , x4 , x5 в z z 16 x1 x2 (10 2 x1 x2 ) 5 (2 x1 3x2 z 6) 5 ( 8 2 x1 4 x2 ), 2 x1 3x2 , z z 2; 3; n 2; 3 grad z. x1 x2 Строим прямую z 0 и вектор n на плоскости Ox1 x2 . Перемещая еѐ в направлении вектора n получаем, что min значение z находится в точке A (рис. 2.5.) пересечение прямых, I и II. Решая систему их уравнений, получаем x1 3, x2 4 . Подставляя в x 3 , x4 , x5 находим оптимальное решение 18 . x 3,4,0,0,14 , z Из геометрических соображений в случае n m 2 можно сделать следующие предварительные выводы: Оптимальное решение ОЗЛО, если оно существует, то лежит обязательно на границах ОДР. Оптимальное решение ЛО может быть не единственным. Это будет, когда прямая z const , z C параллельна одной из сторон ОДР. 18 Оптимального решения может не существовать, даже когда ОДР существует. Это будет, когда вектор n z направлен в сторону неограниченности ОДР. Оптимальное решение достигается всегда в одной из вершин ОДР (если оно достигается на целой стороне, то оно же достигается и в каждой из вершин, через которые проходит сторона). Решение лежащее на одной из вершин ОДР называется опорным (базисным) решением, а сама вершина опорной (базисной) точкой. Для того, чтобы найти оптимальное решение ОЗЛО достаточно перебрать все вершины ОДР и выбрать из них ту, где достигается min (max) z . Оптимальное решение достигается в одной точке, где две из переменных x j обращаются в 0. В заключение, отметим, что, если n m 3 , то свободных переменных три и геометрическая интерпретация ОЗЛО должна строиться в трѐхмерном пространстве. Смысл нахождения ОДР как части трѐхмерного пространства, образованного пересечением полупространств сохраняется, но в силу невозможности построения теряется наглядность. Поэтому, графически решают ОЗЛО только тогда, когда n m 2 . 3. ВЫПУКЛЫЕ МНОГОГРАННИКИ И ЛИНЕЙНАЯ ОПТИМИЗАЦИЯ В задачах оптимизации, как с ограничениями, так и без ограничений, важную роль играет понятие выпуклости. Для большинства методов решения задач оптимизации существование решения может доказать только лишь при условиях выпуклости. 3.1. Выпуклые множества и свойства допустимой области х1( 2) х1(1) Пусть х (1) х2(1) ... , х ( 2) х2( 2) ... х1( r ) , …, х х2( r ) (r ) некоторые точки ... хn(1) хn( 2) хn( r ) n -мерного пространства Rn. Определение 1. Выпуклой линейной комбинацией точек x 1 ,...,x вается точка х 1х (1) 2х ( 2) ... rх (r ) r i 1 ix (i ) , r назы- 19 где i 0и r i 1 1. i Определение 2. Выпуклая линейная комбинации двух точек Rn называется отрезком, соединяющим эти точки. x (1) х ,х ( 2) R 1x n 1 2 1, 2 (1) 2x (2) отрезок, 1, 0. х (1) (1 ) х ( 2) , где Часто отрезок записывают в виде х 0,1 . n Определение 3. Множество S R называется выпуклым, если оно содержит вместе с любыми двумя своими точками и соединяющий их отрезок, x x 2 x 1 x 1 x x 2 x 2 1 1 x x x 2 1 2 Рис. 3.1 0.1 : x (1) (1 ) x ( 2) S . х (1) , х (2) S Выпуклые множества в Rn: куб, шар, пирамида и т.п. Невыпуклые множества: сфера, окружность, кольцо, тор. Примеры выпуклых и невыпуклых множеств на плоскости приведены на рис. 3.1. Теорема 1. Множество S Rn выпуклое тогда и только тогда, когда любая выпуклая линейная комбинация точек из S принадлежит S .  Необходимость. Докажем по индукции. Пусть S – выпуклое множество и x1 , x 2 ...x r S . Если r 1, то х (1) S ; Если r 2 , то т.е. х (1) k i 1 (1 (i ) ix ) х ( 2) S . Пусть утверждение справедливо для k S, ( i 0: i i 1 1) . r k , т.е. 20 Докажем r для k 1, т.е. докажем, k 1 что i 1 i k 1 i 1 k k i 1 k 1 ix Т.к. S ix (i ) k 1x k ( k 1) ( i 1 i x k 1 1, i 1 k (i ) i , i 1 k 1 i 1 k (i ) 1. i Обозначим x ix i 1 ) x( k (1 1) . i 1 k 1 i ) x(i ) i x (i ) S , то тогда x k ( i ) x (i ) 1 i 1 i 1 S по определению выпуклого множества. Достаточность доказать самостоятельно.  Определение 4. Полупространством в Rn называется множество точек xT x1 ,...,xn координаты которых удовлетворяют неравенству i n ajxj b ( j 1 n ajxj j 1 b) , a j , b R. Определение 5. Гиперплоскостью в Rn называется множество точек xT x1 ,...,xn координаты, которых удовлетворяют линейному уравнению. n ajxj b. j 1 Ясно, что любая гиперплоскость может быть представлена пересечением двух полупространств. Определение 6. Многогранным множеством называется пересечение конечного числа полупространств. Ограниченное многогранное множество называется многогранником. Из постановки задачи ЛО, ясно, что ОДР, если она существует, является многогранным множеством. Теорема 2. Полупространство является выпуклым множеством.  Пусть x 1 и x 2 две произвольные точки полупространства, т.е. для координат выполняются неравенства: n a j x (j1) j 1 bи n a j x (j2) j 1 b. 21 x (1) Покажем, что отрезок x странству. xj Для ) x 2 следующее: n n ajxj n ) a j [ x (j1) (1 j 1 a j x (j2) b 1 xT точек x j (1) (1 j 1 (1 координат ) x ( 2) принадлежит полупро- (1 ) x (j2) ] x1 ,...,xn n имеем в силу a j x (j1) j 1 b b. j 1 Следовательно, x x (1) (1 ) x (2) — точка полупространства.  Теорема 3. Пересечение конечного числа выпуклых множеств есть выпуклое множество.  Докажем для двух множеств. Пусть S1 и S 2 - выпуклые множества, и x 1 ,x x (1) 2 (1 S1  S 2 . Т.к. x (1) , x ( 2) ) x ( 2) S1 и x (1) S1 и x 1 , x (1 ) x ( 2) 2 S2 , S 2 , а S1 и S 2 выпуклые, то 0,1 . Следовательно, x (1) (1 ) x ( 2) S1  S 2 . Для любого конечного числа множеств доказывается по индукции.  Следствие 1. Гиперплоскость является выпуклым множеством. Следствие 2. Многогранное множество, если оно не пусто, является выпуклым множеством. Определение 7. Точка x выпуклого замкнутого множества называется крайней, если она не может быть представлена в виде выпуклой линейной комбинации каких–либо двух различных, отличных от x , точек этого множества. На рис. 3.2 S - выпуклое многогранное пространстX во в R2. Точка x не может быть представлена как выпуклая линейная комбинация других точек S ,т.е. лежать на S отрезке, целиком лежащим в S . Из предыдущего следует, что крайние точки ОДР ОЗЛО могут находиться на пересечении гиперплоскостей. Рис.3.2 Ясно также, что ОДР может иметь лишь конечное число крайних точек. Крайние точки многогранника называются вершинами. 3.2. Свойства допустимых решений ОЗЛО Рассмотрим следующую задачу (ОЗЛО): 22 n z z cT x min, или Ax b, x 0, cjxj n aij x j cТ bТ (c1 , c2 ...cn ) , bi (i 1, m), (3.1) j 2 xj где min, j 1 0. (b1 , b2 ...bm ) , xТ ( x1 , x2 ...xn ) , A (aij ), i 1, m; j 1, n . T ( x1* , x2* ,...,xn* ) . Пусть Решение задачи (3.1) будем обозначать x* система Ax b совместна и r A r m . Определение 8. Базисом называется всякая не особая квадратная подматрица порядка m матрицы A . Ясно, что такая подматрица существует, т.к. r m . Пусть B — некоторый базис, тогда перестановкой столбцов матрицы A еѐ всегда можно привести к виду: a11 .. a1m a1m 1 .. a1n А (В N ) ; B .. .. .. ; N .. .. .. ; a m1 .. a mm a mm 1 .. a mn где N – подматрица, состоящая из столбцов матрицы A , не принадлежащих базису B . Размерность N : m n m . Точно также можно произвести разбиение векторов x и c : T T xB х xB x N x1, x2 ,..., xm xm 1,..., xn ; xN T T cB cB c N c1,..., cm cm 1,..., cn . cN Переменные x1,...,xm соответствуют базису, и векторы представлены как клеточные. Тогда систему уравнений ограничений можно записать в виде xB BN b; (3.2.) xN c BxB Nx N b, x1 ,..., xm , x TN x m 1 ,...,x n . где x TB Переменные входящие в x B называются базисными, а в x N – свободными. Система уравнений (3.2) имеет бесчисленное множество решений. 23 Определение 9. Базисным решением называется частное решение системы уравнений (3.2), полученное при условии, что x N 0 . Тогда при x N 0 из (3.2) получаем базисное решение (3.3) x B B 1b . Ясно, что базисное решение будет допустимым, т.е. удовлетворять системе ограничений x j 0 , если х В 0 , т.е. х1 0,..., х m 0 . В этом случае его также называют опорным. Соответствующий базис называется допустимым базисом. Допустимое базисное решение называется невырожденным, если оно имеет ровно m положительных компонент, т.е. xB 0 (в x B нет нулевых компонент). Умножим (3.2) слева на B 1 , получим Х В В 1Nx N B 1b , хВ B 1b B 1Nx N B 1 b Nx N . (3.4) Таким образом, выразили базисные переменные через свободные. Обычно на практике это делается методом Гаусса. Теорема 4. Допустимое решение x является крайней точкой многогранного множества допустимых решений тогда и только тогда, когда оно базисное (множества крайних точек и опорных решений совпадают).  Достаточность. Пусть x – опорное решение (допустимое базисное B xB B 1 , b 0 , xT ( x1B , x2B ,..., xm ,0,...,0), х N 0, решение). Тогда B x1B , x2B ,..., xm 0 , m n . Докажем, что x – крайняя точка ОДР. От противного. Предположим, что x не является крайней точкой, тогда, еѐ можно представить х х1 1 х 2 , где х 1 , х 2 ОДР и х1 x х 2 T 2 . Но 0и1 0 , тогда, так как x(1) T (2) x1(2) , x2(2) ,..., xm ,0...,0 , то точки x 1 и x ниям ограничениям: х (1) (2) В 1 1 b, x(1) x(2) 2 (1) x1(1) , x2(1) ,..., xm ,0,...,0 , удовлетворяют уравнеx(1) x(2) , т.е. они сов- x B b, падают. Противоречие. Необходимость. Пусть x — крайняя точка, тогда часть переменных обращается в ноль, хT x1 , x2 ,..., xk ,0...0 — крайняя точка ОДР, причѐм яс- но, что x j 0 , т.к. это допустимое решение. Покажем, что x — базисное решение, то есть надо доказать, что k m. От противного, пусть k m . Т.к. x крайняя точка, то она удовлетворяет системе уравнений ограничений (3.2) которую можно записать в виде 24 (3.3), причѐм, учитывая, что хk 1 xk 2 .. xn получим: 1 х1 0 х2 0 х3 ... 0 хk 1, 0 x1 1 x2 0 x3 .. 0 xk 2, 0 x1 0 x2 1 x3 ... 0 xk 3, ........, 0 x1 0 x2 0 x3  1 xk k. 1 , , Ak  1 2 где B 1b  A1 , k 0 , то в развѐрнутом виде ,  1 или x1 A1 x 2 A2  x k Ak B 1 b . Т.к. r m и k m , то система столбцов A1, A2 ,, Ak симая, а следовательно существуют такие 1, 2 ,, k , что k j 1 j Aj k 0 , где 2 j j 1 линейно зави- 0. (3.5) С другой стороны k B 1b . х j Aj (3.6) j 1 0 и сложим, а потом вычтем из (3.6). Умножим равенство (3.5) на k xj aj Aj B 1 b, j 1 (3.7) k xj aj Aj B 1 b. j 1 Рассмотрим две точки x 1 и x x1 x1 1  x 1 xk с компонентами x j j и xj j 1  k  2 , x 2 xk k  . (3.8) 25 Первые k компонент за счѐт выбора можно сделать положительными, следовательно, х 1 , х 2 ОДР, т.к. х 1 B 1b , x 2 B 1b . Тогда 1 1 1 2 х х (точка решения системы ограничений) – являиз (3.7) точка х 2 2 ется выпуклой линейной комбинацией точек x 1 и x 2 , которая не является крайней точкой. Противоречие. Следовательно, k m и x допустимое базисное (опорное) решение. Следствие. ОДР ОЗЛО, т.е. множество вида Х х : Ax b, x 0 имеет конечное число крайних точек (вершин), которое не превосходит C nm . Действительно, максимальное число базисов, извлечѐнных из A есть число возможностей выбора m столбцов из n , т.е. C nm , но не все они могут быть допустимы. Теорема 5. Любая точка ОДР ОЗЛО есть выпуклая линейная комбинация еѐ крайних точек.  Пусть Х х : Ax b, x 0 есть ОДР ОЗЛО. Возьмѐм х Х . Тогда, х1, х2 ,, хk ,0,,0 , причѐм x1, x2 ,, xk 0, и k m , то x базисесли хТ ное решение и по теореме 4 это крайняя точка, а, следовательно, теорема доказана. Если x1, x2 ,, xk k и m, то существует y (1) точка T (1) y1 1 , y2 1 ,..., ym ,0,...,0 у которой все еѐ компоненты y1(1) , y2(1) ,, ym со1 держатся среди ненулевых компонент точки x и выполняется условие Ay 1 b, y 1 0, y 1 y (1) j x, x1, x2 ,, xk , то есть y 1 является крайней точкой. Обозначим x 1y 1 ;x 1y 1 xj min 1 0, x 1y j 1, m . 0, y (1) j j 1 2 такую, что y 2 0, Ay 2 2 y1 0 . Отсюда x 2y через конечное число шагов p 1 точек y , y x 1y 1 2 ,, y 2y 2 p ... b, 1y p что 1 y 2 2y 1y y 1 и число 1 , найдѐм точку 2, такое, что 2 0 . Продолжая этот процесс, m (так как r A m ) получим p крайних и чисел py ясно, Rn – новая точка. Действуя точно также относительно точки x y Тогда 1, 2 ,, p , удовлетворяющих равенству: Следовательно, x 1 1y 2y 2  26 p . Т.е. x – линейная комбинация этих крайних точек. Покажем, что это выпуклая линейная комбинация. Во-первых, 1, 2 ,, p 0 по построению, p тогда из Ax b следует A( 1 y (1) 2 y (2)  1 Ay (1) 1b 2 Ay 2b (2) py p Ay  pb  ( p) ) b; ( p) b; p b, 1. j j 1 Следовательно, x есть выпуклая линейная комбинация p крайних точек.  Теорема 6. (Необходимое условие оптимального решения). Решение ОЗЛО достигается в вершине многогранника ОДР или в крайней точке ОДР. Если оптимум достигается не в одной вершине, а в нескольких, то он достигается на всѐм многограннике, порождѐнном этими вершинами. 1 2  Пусть х , х ОЗЛО, т.е. z* z x* p ,, х - вершины ОДР, а x оптимальное решение x ОДР . Если x – вершина ОДР, то теорема z x доказана. Пусть x не является вершиной ОДР (или крайней точкой ОДР), p тогда по теореме 5 имеем х* p i ix , i 1 i 1, 0 . Т.к. z x линейная i i 1 функция, то zx * p z ix i 1 pz m Пусть z( x m zxm x n i p cj j 1 i 1 ) min z x 1 ... p 1 p i ixj c j x ji i i 1 j 1 ,, z x zxm , p zx m i 1 z x* . Тогда zx p n iz 1z xi . x m  . Но по определению, для любого x принадлежащего ОДР, z* z x откуда следует, что z zx m , т.е. оптимальное значение z достигается в вершине x m . Допустим, что z принимает наименьшее значение более, чем в одной точке, например, в x 1 , x 2 ,...,x s , т.е. z(x(1))=…=z(x(s))=z*. Рассмотрим произвольную выпуклую линейную комбинацию этих точек. Это будет некоторая точка x многогранного множества, порождѐнная этими точками: x 1х (1) 2х ( 2)  sx (s) s ; i 1 i 1 i 0. 27 Тогда zx z S i 1 ix S (i ) i 1 iz x (i ) S i 1 iz * z S * i 1 i z* . Следовательно, z достигает минимума в произвольной точке многогранного множества с вершинами x 1 , x 2 ,...,x s .  Итак, если ОДР – многогранник, то для решения ОЗЛО надо найти значения z во всех вершинах ОДР и выбрать наименьшее. Следовательно, теорема 6 даѐт необходимое условие оптимального решения. Теперь задача состоит в нахождении такого допустимого базиса B , при котором базисное решение x B являлось бы оптимальным. 4. СИМПЛЕКС – МЕТОД РЕШЕНИЯ ЗАДАЧ ЛИНЕЙНОЙ ОПТИМИЗАЦИИ 4.1. Оптимальный допустимый базис и оптимальное решение Пусть имеется ОЗЛО: z cT x min, Ax b, x 0. Обозначим через B базис матрицы А размерностью m m и базисные переменные x1,...,xm , а свободные – xm 1 , xm 2 , ..., xn . По (3.4) имеем хВ B 1b B 1NxN , n k m, где B c 1 (4.1) – обратная матрица для B , b b1  bm , xN xm xm  xm два 1 2 k x1 x2 , xB .  xm вектора, т.е. Если вектор с также разбит на cB , cT cTB cTN , то целевую функцию можно представить в виде: cN z cT x cTB cTN , z cTB B 1b Обозначим xB xN cTB xB B 1Nx N cTN xN . Подставим сюда хВ из (4.1), получим cTN x N cTB B 1b cTN cTB B 1 N x N (4.2) 28 cTB B 1 ; В 1b ; cTN N cTN cTB B 1N ; 0 b. Тогда из (4.2) (4.3) z b cTN N xN . После разделения на свободные и базисные переменные получили следующую задачу ЛО: z xN min, (4.4) xB B 1Nx N 0. Вектор cTB B 1 1, 2 ,, m называется вектором симплексных cTN N множителей. Компоненты вектора 1, 2 ,, n m называются приведѐнными значениями или приведѐнными ценами свободных переменных. Обсудим идею метода решения задач ЛО. Для этого вначале запишем задачу (4.4) в развѐрнутом координатном виде: (4.5) z 0 1xm 1 2 xm 2  S xs m  k xn m k , x1 1 x2 2  xl 11xm 1 21xm 1 12 xm 2 22 xm 2   1s xm s 2 s xm s   1k xn , 2 k xn ,             , l l1xm 1 l 2 xm 2  ls xm s  lk xn , (4.6)              , xm m m1xm 1 m 2 xm 2  ms xm s  mk xn . Матрица системы (4.6) размерности ( m k ) имеет вид 11 12  1k 21 22  2k B 1N .     m1 m2  mk xm 1 xm 2  xn 0 и вектор столбец Пусть x N 0 (4.7) 0 , т.е. xB является опорxN ным (допустимым базисным) решением, но неизвестно будет ли оно оптимальным. Рассмотрим различные ситуации. Если все компоненты вектора γ в 0 , то всякое изменение опорного решения, а имен(4.5) положительны но, увеличение каких либо свободных переменных сверх нуля приводит к 1 0,, m 0 . Это означает, что вектор-столбец x 29 тому, что z может только возрастать. Следовательно, исходное опорное реxB шение x является оптимальным. xN Пусть некоторые коэффициенты в (4.5) отрицательные, допустим ys 0 . Это означает, что, если в (4.6) переменную xm s увеличивать сверх нуля, то целевая функция z будет уменьшаться, и исходное опорное решение xB x не является оптимальным, а следовательно базис В не являxN ется оптимальным. Поэтому необходимо изменить базис и найти новое опорное решение. Причѐм ясно, что надо увеличить xm s сверх нуля, т.е. перевести еѐ из свободных переменных в базисные, а некоторую другую переменную из базисных в свободные. Если среди уравнений системы ограничений (4.6) нет уравнения с положительным коэффициентом при xm s т.е. столбец матрицы B 1N , соответствующий xm s, отрицательный ( B 1N ) s 0 , то, увеличивая xm s сколько угодно, все переменные x1, x2 ,..., xm останутся положительными и z неограниченно уменьшается. Следовательно, оставаясь в ОДР, z min не существует, и оптимального решения нет. Если среди уравнений системы (4.6) есть уравнение с положительными коэффициентами при xm s , например, в уравнении xl l l1 xm 1 ... ls xs .... lk xk коэффициент ls 0 . Тогда, увеличивая xm s до значения l lS (иначе xm s станет отрицательным), а остальные свободные пе- ременные оставляя равными нулю, получим, что xl 0 , а хm s l . lS Таким образом, переменные x1, x2 ,..., xi 1, xi 1, xm s 0 станут базисными, а переменные xm 1 ,..., xn , xl – свободными. Теперь надо найти новый базис В и выразить целевую функцию z через новые свободные переменные. Если все новые коэффициенты в z положительны, то получаем оптимальное решение. Если есть отрицательные, то полученное решение является опорным, но не оптимальным, и процедура улучшения решения продолжается пока не будет найден оптимальный базис или будет показано, что ОЗЛО решения не имеет. В заключение заметим, что если в столбце B 1N s имеется несколько положительных элементов, а не только x1, x2 ,..., xm для которой отношение ls , то выбирают ту переменную из i is , i 1, m будет минимальным, т.е. 30 i min находят i: is 0 , и ту переменную из базисных, которая быстрее об- is ратится в ноль. Приведѐнное рассуждение поиска оптимального решения обосновывается следующими двумя теоремами. Теорема 1. (об оптимальном базисе). Для того, чтобы матрица В была оптимальным допустимым базисом необходимо и достаточно чтобы (4.8) 0 cTN N 0. xB любое решеxN ние, удовлетворяющее ограничениям ОЗЛО (не обязательно базисное) и z(х) соответствующее значение целевой функции z. Если В некоторый базис, тоB 1Nx N ; гда после разделения переменных получаем (4.2)-(4.3), т.е. xB  Достаточность. Пусть выполняется (4.8), x z ( х) xN , где cTN b, Т. к. выполняется (4.8), т.е. Значение cTB N; 0 и xN cTB B 1 . 0 , то z ( x) b cTB B 1b . получается при допустимом базисном решении (опорном решении), когда xN 0, xB B 1b 0 в точке x * . При дру- гих х z (x) 0 , следовательно z ( x ) z оптимальное решение, а B B допустимый оптимальный базис. Необходимость. Пусть В* — оптимально допустимый базис, которому соответствует оптимальное решение z* и которое достигается в точке x * х *В х *N . Покажем, что если существует свободная переменная с индексом m s , удовлетворяющая условию S 0 , то можно выбрать решение с меньшей ценой, чем 0 , а следовательно z* не является оптимальным решением. Сделаем замену свободных и базисных переменных, т.е. по-другому выберем x B и x N , так, чтобы только переменная xm s изменила своѐ значение и вместо 0 она станет равной некоторой θ>0. xm 0, в решении х*; 0, в решении х; s Тогда новая точка x xN x*N eS xB xN удовлетворяет равенству eS ; eTS 0 0  1 0 0 ; 31 где es – вектор-столбец той же размерности, что и x N , т.е. k=n–m, в котором все компоненты равны нулю, за исключением компоненты s , которая равна xB 1. Число надо подобрать таким, чтобы x оставалась решением, т.е. xN x B B 1b B 1 Nx N B 1b B 1 NeS 0 . Произведение NeS есть столбец АS матрицы N, который соответствует переменной xm s Обозначим As B 1 As , тогда m 1  mm 1 Если   1n 1m s , xB .    1 NeS AS As 0  mn mm s  всегда можно подобрать так чтобы x B 0 , тогда для 0 , то eS z * x 0 , по z x будем иметь z ( х) 0 х N z ( x* ) S . Т.к. S 0 , а отсюда следует, что z(x)>0. z cjxj My1 j 1 My2 My3 min . ―Запретительная цена‖ (оценка) делает невозможным присутствия хотя бы одной из искусственных переменных в оптимальном решении. 2. Пусть ограничения заданы в виде неравенств x1 x2 x3 5, 3 x1 x2 2 x3 4, 2 x1 x3 1. Введѐм в неравенства дополнительные переменные y1, y2, y3 чтобы получить ограничения типа равенств. x1 x2 x3 y1 5, 3x1 x2 2 x3 y2 4, 2 x1 x3 y3 1. На главной диагонали единичной подматрицы будут стоять элементы (–1.) Умножать уравнения на (–1), нельзя, т.к. в правой части системы уравнений получаются отрицательные свободные члены, а, следовательно, базисное решение с y1, y2, y3 в качестве базисных переменных не будет допустимым. Поэтому вводим три искусственных переменных y4, y5, y6 с коэффициентами (+1), которые образуют искусственный единичный базис, а именно x1 x2 x3 y1 0 y2 0 y3 1y4 0 y5 0 y6 5, 3 x1 x2 2 x3 0 y1 2 x1 0 x2 xi вид: 0, y2 0 y3 0 y4 1 y5 x3 0 y1 0 y2 yj y3 0 y4 0 y6 4, 0 y5 1 y6 1, 0. Здесь y4, y5, y6 –базисные переменные, а целевая функция будет иметь z 3 cjxj j 1 0 y1 0 y2 0 y3 My4 My5 My6 min . 3. Пусть ограничение смешанные. В этом случае базис образуется так же одновременно введением дополнительных и искусственных переменных. 41 2 x1 7 x2 x3 37, 5 x1 9 x2 x3 24, x1 x2 x3 10. В первое ограничение введѐм дополнительную переменную y1, в третью – дополнительную переменную y3 со знаком минус. Получаем ограничения-равенства. Для получения базиса во второе и третье уравнение вводим искусственные переменные y2 и y4 с коэффициентами (+1). Получим задачу: 2 x1 7 x2 x3 0 y3 y1 0 y2 0 y4 37, 5 x1 9 x2 x1 x2 xi 3 z x3 0 y3 0 y1 1y2 0 y4 24, x3 1y3 0 y1 0 y2 1y4 10, 0, y j cjxj 0. 0 y1 j 1 0 y3 My2 My4 min . Метод элементарных преобразований Пусть имеется ОЗЛО: z c T x min; Ax b ; x 0 . 1. Преобразование уравнений-ограничений к стандартному виду. Если нет естественного допустимого базиса, полученного с помощью дополнительных переменных, то методом элементарных преобразований матрицы преобразуем расширенную матрицу системы ограничений ( A b) B N b к виду E N ; N B 1N . Без ограничения общности можно считать, что Е — единичная подматрица соответствующая первым m столбцам матрицы А. Таким образом, ограничения будут представлены в виде xi N n ij x j j m 1 i , i 1, m, j m 1, n, . Следовательно, х1, х2,…,хm — базисные переменные, хm+1, хm+2 ,…,хn – свободные. Целевая функция выражается через свободные переменные одновременно при выполнении процедуры Гаусса: ij z n jxj j m 1 и записывается в виде 2. Поиск опорного решения. 1 xm 1  n m xn z 42 2а. Если i 1, m : i 0, то решение x 1, 2 ,, m ,0,0,,0 является опорным, Е — допустимый базис, и переход к поиску оптимального решения (п. 4.3). 2б. Если i 1, m : i 0 и j m 1, n : ij 0 (т.е. в i-ой строке матрицы N все элементы ij 0 ) то, система ограничений x j 0, i 1, m, x j 0 , j m 1, n несовместна, а значит, что допустимое решение отсутствует, оптимального решения нет. 2в. Если i 1, m : i 0 и j m 1, n : ij 0 т.е. в i-ой строке и в jом столбце матрицы N есть отрицательный элемент, то необходимо продолжение поиска опорного решения (замена базиса). 3. Замена базиса. В i-ой строке, где i 0 и ij 0 , выбирается отрицательный iq max ij ij элемент iq с наибольшей абсолютной величиной и столбец q, в котором он находится, будет разрешающим. В этом столбце ищется разрешающий элемент. Для этого рассматриваются все элементы данного столбца, имеющие одинаковый знак со свободным членом. Из них выбирается тот элемент pq , для которого отношение к нему свободного члена минимально. Выбором элемента строка p : q pq min i: ij q 0 q pq определяется разрешающая . iq Производится переход к новому базису B заменой x p возxq вращение в 2а. Выход из программы в 2а. или 2б. Замечание 1. Замену свободных и базисных переменных при табличном алгоритме можно делать методом элементарных преобразований (Гаусса). В программном варианте лучше использовать матричное произведение В В F (замечание в 4.3). Замечание 2. Если в 2б есть несколько строк, где i 0 (т.е. отрицательные свободные члены в ограничении), то выбирается строка, в которой i max . Пример. 4.2 43 z x1 2 x2 3 x3 min x1 2 x2 2 x1 x2 x1 x2 x3 1, x3 5, 2, x2 x3 1, x1, x2 , x3 0 . Решение. Преобразуем систему уравнений-неравенств в систему уравнений-равенств y1 x1 2 x2 x3 1, y2 y3 2 x1 x2 x3 5, x1 x2 0 x3 2, y4 0 x1 x2 x3 1. Чтобы получить единичную подматрицу Е умножим второе уравнение на (–1) и заполним симплекс-таблицу. 44 N Б x 1 z y 10 x 2 1 1 y 2 - y 1 01 z y 3 3 2 - 2 2 x y 1 3 1 - 1 1 1 - 1 2 3 y 2 3 1 1 Таблица 4.3. y y z 4 1 5 1 i 1 - 5/-2=2,5 2 2/ 1 1=2 y4 11 x y 1 2 1 01 1 y - 1 - 1 - 3 1 1 - 1 1 2 1 1 1 2 1 1 3 5 3 1 4 z y 1 x 32 x 1 1 13 3/ 1=3 2 1/-1=1 1 1/ 1=1 - 5 2 1 2 1 1 1 1 3 2 1 1 2 1 2 1 2 y 4 Получим y1=2, x3=1, x1=2, y4=0, x2=0, y2=0, y3=0. Это допустимое решение. Оно является и оптимальным, поскольку в строке z все коэффициенты неотрицательные и z * 5 . В заключении этой главы отметим особые случаи в симплекс – методе, т.е. сделаем анализ. 1.Отсутствие допустимых решений. В этом случае ограничения противоречивы, и это означает, что модель построена некорректно. Необходимо построить модель, имеющую иную структуру. 2. Неограниченность. В этом случае возможно неограниченное увеличение переменных без нарушения наложенных ограничений, т.е. пространство решений, по крайней мере, в одном направлении неограниченно и целевую функцию z можно сделать как угодно малой. Разработанная модель не- 45 достаточно точна. Ясно, например, что бессмысленно использовать модели, прогнозирующие ―бесконечную прибыль‖. Обычно в этом случае: а) не учтены одно или несколько ограничений (не хватает ограничений); б) неточно оценены параметры, входящие в ограничения. Необходимо корректировать модель. 3. Вырожденность. Так называется ситуация, когда один или больше свободных членов в уравнениях-ограничениях равны нулю, т.е. некоторые базисные переменные в опорном решении равны нулю. В этом случае может оказаться, что замена х j xi приводит только к перестановке переменных без уменьшения функции z и через конечное число замен получается тоже, что уже было. Это явление называется ―зацикливанием‖. Наличие вырожденности не свидетельствует, о какой-то ―опасности‖ для исследователя. С практической точки зрения это означает, что в модели есть, по крайней мере, одно ―избыточное‖ ограничение. Чтобы избежать зацикливания используют следующие приѐмы: а) надо взять на каком-либо этапе другой разрешающий столбец; б) несколько изменить коэффициенты ограничений; в) применить лексиграфическую процедуру, т.е. расположить базисы в лексиграфическом порядке. Каждый раз, когда возникает вырожденность, приходят к другому базису в соответствии с лексиграфически возрастающим последовательностям базисов. Таким образом, никогда не встретится базис уже встречавшийся ранее. Это обеспечивает, что через конечное число итераций z снова будет убывать. 4. Альтернативность. Альтернативные решения возникают тогда, когда оптимальное решение достигается не в одной крайней точке, а на некотором многограннике, порожденном крайними точками. Здесь z принимает одно и тоже значение во всех точках этого многогранника. Существование альтернативного решения на практике может быть удобно, т.к. x* можно выбрать используя ещѐ какие-либо дополнительные соображения. 5. ДВОЙСТВЕННОСТЬ В ЛИНЕЙНОЙ ОПТИМИЗАЦИИ 5.1. Двойственные задачи Двойственность является важным понятием в ЛО, приводящим к результатам теоретического и практического характера: теоремам двойственности. Для каждой задачи ЛО можно поставить другую задачу ЛО в совершенно другом пространстве переменных. Если известна взаимосвязь между этими задачами, то можно получить решение одной из них, когда известно решение другой. Пусть необходимо найти такой вектор - столбец x, что z cT x min; Ax b; x 0 . (5.1) 46 b1 c1 x1 А (aij ); i 1, m; j 1, n; b  ; c  ; x  . bm cn xn Эту задачу будем называть прямой или исходной. Двойственная задача к (5.1) формулируется следующим образом: найти такой вектор-столбец переменных u , что (5.2) w bT u max; AT u c, u o, где u T (u1 , u 2 ,...,u m ) . Переменные u1,u2,…,um называются двойственными оценками или теневыми ценами. Такую форму постановки прямой и двойственной задач называют стандартной. В двойственной задаче вектору правых частей ограничений в (5.1) ставится вектор коэффициентов целевой функции в (5.2), столбцы матрицы А ограничений (5.1) становятся строками матрицы ограничений в (5.2). Каждая двойственная переменная u1, u2, …, um связана с соответствующим ограничением общего вида прямой задачи, т.е. каждому i-му ограничению ставится в соответствие новая переменная ui Пример. Привести к стандартной форме задачу ЛО и сформулировать двойственную задачу. z x1 4 x2 3x3 max 3 x1 4 x2 x3 x1 2 x2 x3 x3 7; 6; 4; x1 0, x2 0, x3 0. Решение. Так как x1 2 x2 то можно записать: z x1 4 x2 3 x3 3 x1 x1 4 x2 2 x2 x3 x3 x3 6 x1 2 x2 x3 6, x1 2 x2 x3 6, min 7; 6; x1 2 x2 x3 6; 0 x1 0 x2 x3 4; x1 0 x2 0 x3 0. Это стандартная форма прямой задачи. Двойственная задача состоит в следующем: найти такие u1 0, u2 0, u2 0, u3 0 , чтобы выполнялись неравенства: 47 3u1 u2 u2 4u1 2 u2 0u3 2u2 1, 0u3 4, u1 u2 u2 u3 3, и функция w 7u1 6u2 6u2 4u3 принимала max значение. Видно, что u 2 u 2 u 2 – может считаться единственной переменной из-за вида коэффициентов, и задача может быть сформулирована следующим образом u1 0, u2 0, u3 0 такие, что w 7u1 6u 2 4u3 max и выполняются неравенства: 3u1 u2 1, 4u1 2u2 4, u1 u2 u3 3. Двойственную задачу можно определить для любой линейной задачи. Для этого можно пользоваться следующей таблицей 5.1 соответствий. Таблица 5.1. Прямая (исходная) задача Двойственная задача Коэффициенты целевой функции Правая часть (min) Коэффициенты целевой Правая часть функции (max) Т A — матрица ограничений А — матрица ограничений i-e ограничение: Переменная ui 0 i-e ограничение:= Переменная x j 0 Переменная Переменная x j j-e ограничение: = ui j-e ограничение: Имеют место следующие теоремы двойственности. Теорема 1. Двойственная задача к двойственной есть прямая задача.  Запишем двойственную задачу (5.2) в стандартном виде как прямую задачу (5.1) и построим для неѐ соответствующую двойственную задачу. Тогда w AT u u 0,  bT u min, z c, c AT x 0, T T x x max, bT T cT x max, Ax b, x 0, z , z cT x min, Ax b, x 0. 48 Теорема 2. Если х и u допустимые решения соответствующих исходной и двойственной задач, то z c T x w b T u , т.е. значение функции z соответствующее любому допустимому решению прямой задачи не меньше значения w, соответствующего допустимому решению двойственной задачи.  Пусть имеются прямая (5.1) и двойственная (5.2) задачи. Тогда, из Ax b, u 0 u T Ax u T b b T u w , а из условий AT u c, x 0 xT AT u xT c cT x z . Т.к. u T Ax x T AT u — число, то z xT AT u u T Ax w z w . Неравенство (5.3) cT x bT u называется основным неравенством двойственности.  Из теоремы следует, что если имеется допустимое значение функции z, равное допустимому значению функции w, то это должно быть минимальное значение функции w. Однако, ясно, что пока ниоткуда не следует, что такие значения существуют. Следствие. Если х* и u* – допустимые решения прямой и двойственной задач, и выполняется равенство с T x * b T u * , то х* есть оптимальное решение задачи (5.1), а u* – оптимальное решение задачи (5.2).  Согласно основному неравенству двойственности х, и : z( x) w(u) . Пусть х — произвольное допустимое решение (5.1), а u* — оптимальное решение (5.2) z ( x) w(u * ) z ( x * ) х* — оптимальное решение (5.1).  Наоборот, пусть u — произвольное решение (5.2); х* — оптимальное решение (5.1) w(u ) z ( x * ) w(u * ) u* — оптимальное решение (5.2). Это следствие называется критерием оптимальности Канторовича. Теорема 3. (первая теорема двойственности). Пусть прямая задача (5.1) имеет конечное оптимальное решение и * — строка симплексмножителей, соответствующая этому оптимальному решению. Тогда оптимальное решение двойственной задачи существует и * оно равно * ,т.е. u * .  По теореме о необходимом и достаточном условии оптимума имеем zx b cTN N x ; cTB B 1; Если * – симплекс множители, соответствующие оптимальному ре* * cTN N 0 и тогда j ( )Tj шению х* задачи (5.1), то Aj 0 , T * (j=m+1,m+2,…,n), для свободных переменных и Aj 0 , j (cB ) j j 1,2...m для базисных переменных, где Аj – j-й столбец матрицы A. Следо- вательно * Аj (cT ) j , j 1,2n. Тогда, допустимое решение задачи (5.2). * А cТ АТ ( * T ) c * - 49 Докажем, что * – оптимальное решение (5.2). Пусть В — оптимальный базис задачи (5.1), т.е. базис на котором реализуется оптимальное решение х*. Тогда, т.к. * c TB B 1 , то *b c TB B 1b c TB x B c T x * z * В силу следствия теоремы 2 это означает, что * – оптимальное решение (5.2).  Теорема 3 соответствует случаю, когда обе задачи (5.1) и (5.2) имеют решение. Следующая теорема охватывает все случаи, которые возможны. Эта теорема даѐт возможность находить решение двойственной задачи не * решая еѐ: u * cTB B 1 . Теорема 4. (вторая теорема двойственности). Пусть заданы линейная задача (5.1) в общем случае с ограничением типа неравенств и двойственная ей задача (5.2). Тогда а) Если задача (5.1) или (5.2) имеет решение, то имеет решение и другая из них, причѐм z * zmin w* wmax . б) Если одна из этих задач имеет неограниченное решение, то другая задача не имеет решения (противоречива система ограничений).  а) Если существует решение задачи (5.1), то справедливость следует из теоремы 3, т.к. при ограничениях типа неравенств можно перейти к стандартной форме введением дополнительных переменных. Если существует решение (5.2), то теорема справедлива в силу теоремы 2 и предыдущего рассуждения. б) Пусть решение задачи (5.2) неограниченно. Это означает, что для любого M 0 существует двойственное решение u : w bT u M . Тогда, если (5.1) имела решение x , то по теореме 2 выполняется неравенство: w bT u cT x при любом двойственном решении а это противоречие, т.к. с одной стороны w>M, а с другой стороны w cT x .  Теперь пусть ограничение в прямой задаче (5.1) все являются равенствами. Имеет место следующая теорема. Теорема 5. (о дополняющей нежѐсткости). Для того, чтобы х и u — решения прямой (ОЗЛО) и двойственной задачи были оптимальными, необходимо и достаточно, чтобы ATj u c j x j 0, j 1, n , (5.4) где Aj — j-й столбец матрицы А.  Необходимость. Пусть х и u — оптимальные решения. Тогда Ax b u T Ax u T b u T Ax c T x u T b c T x . Т.к. х и u оптимальные решения, то по теореме 5.2 cT x bT u u T b а, AT u c T x 0 . следовательно, u T Ax c T x 0, u T A c T x 0 50 Поскольку x 0 и AT u cT 0 (из ограничения (5.2)), то скалярное произведение будет равно нулю, если каждое слагаемое в нѐм равно нулю, т.е. ( ATj u ci ) x j 0 . Достаточность. Пусть выполняется (5.4). Тогда AT u cT x 0 , т.к. все слагаемые равны нулю, то скалярное произведение равно нулю. Откуда u T A c T x 0 u T Ax c T x 0 u T b c T x . Следовательно, u и х — оптимальные решения двойственных задач.  Ясно, что если все ограничения в двойственной задаче равенства, то можно доказать аналогичную теорему с необходимым и достаточным условиями оптимальности в виде: (5.5) ( Ai x bi )ui 0 i 1, m , Ai строка А. Условия (5.4) и (5.5) называются условиями дополняющей нежѐсткости. Эти условия означают: а) если какое-либо ограничение одной из задач с еѐ оптимальным решением обращается в строгое неравенство, то соответствующая компонента оптимального решения двойственной задачи должна равняться нулю; б) если же какая-либо компонента оптимального решения одной из задач положительна, то соответствующее ограничение в двойственной задаче с еѐ оптимальным решением должно обращаться в строгое равенство. ui* 0 , если, ATj u* c j Таким образом, если Ai x * bi x*j 0 . Ai x * bi . В случае б), если x*j 0 ATj u* c j , а если ui* 0 Экономическая интерпретация условий дополняющей нежесткости: если по некоторому оптимальному плану х* производства расход i-го ресурса строго меньше его запаса bi, то в оптимальном плане соответствующая двойственная оценка равна нулю. Если в некотором оптимальном плане оценок iя компонента строго больше нуля, то в оптимальном плане производства расход соответствующего ресурса равен его запасу. Таким образом, двойственные оценки могут быть мерой дефицитности ресурсов. Дефицитный ресурс (т.е. полностью используемый по оптимальному плану производства) имеет положительную оценку, а ресурс избыточный (не полностью используемый) имеет нулевую оценку. Теперь выясним содержательный смысл двойственных оценок (двойственных переменных ui). Решение задачи (5.1) можно рассматривать как функцию правых частей системы ограничений, т.е. z z b1, b2 ,...,bm . Исследуем эту функцию. Пусть для некоторой фиксированной правой части В — * оптимальный базис задачи (5.1) и пусть u * c TB B 1 есть оптимальные двойственные переменные. Для j-го индекса свободной переменной имеем 51 (сTN ) j uA j b b b b1 ( j m 1, m 2,, n) . b1 ;...;bm bm T Зададим b приращение b и рассмотрим задачу z cT x min, Ax b , x 0. Заметим, что u * c TB B 1 не зависят от b и всегда остаются двойственb 0 , то базис B остаѐтся одными решениями. Если новый вектор b b новременно прямым и двойственным допустимым базисом, а следовательно и оптимальным. Значит B будет оптимальным для любого b , и, следоваb) 0 , если решение не вырожденное, т.е. B 1b 0 Тогда тельно, B 1 (b для малых b имеем z (b ) = u* bl T (b z (b b) w(u * ) (b )T (u * ) (u * )T b b) (u * )T b (u * )T b , z (b b) z (b) (u * )T b . Пусть все bi 0 , кроме bl 0 . Тогда переходя к пределу при z (b) u l* . Если записать в конечных приращениях 0 получим bl z (b bl ) z (b) ul* ; bl z (b bl ) z (b) z (bl ) ul* bl . Таким образом, величина двойственной оценки (двойственной переменной) численно равна изменению целевой функции при изменении соответствующего свободного члена ограничения на единицу. Рассмотрим экономическую интерпретацию двойственной задачи на примере. Пусть для производства трѐх видов изделий A, B, и C используются три различных вида сырья. Каждый из видов сырья может быть использован в количествах 180, 210, 244 кг. Нормы затрат каждого из видов сырья на единицу продукции данного вида и цена единицы продукции каждого вида заданы в таблице. Определить план выпуска продукции, при котором обеспечивается еѐ максимальная стоимость при условиях, что используемое сырьѐ каждого вида не превосходит заданный объѐм. Таблица 5.2 Вид Нормы затрат сырья (кг) на единицу проОбъсырья дукции ѐм сырья А В С I. 4 2 1 180 52 II. III. Цена ед. продукции 3 1 1 2 3 5 10 14 12 210 244 Обозначим: x1 – произведено изделий A, х2 – произведено изделий B, х3 – произведено изделий C. Исходная (прямая задача) имеет вид: z 10 x1 14 x2 12 x3 max 4 x1 2 x2 3 x1 x3 180, x2 3 x3 210, x1 2 x2 5 x3 244, x1, x2 , x3 0. Формально построим двойственную задачу по описанному выше правилу. Обозначим двойственные переменные u1, u2, u3. w 180u1 210u2 224u3 min 4u1 3u2 u3 10, 2u1 u2 2u3 14, u1 3u2 5u3 12, u1, u2 , u3 0. Формулировка двойственной задачи: найти такую совокупность u1, u2, u3 цен каждого вида сырья (оценки сырья каждого вида), при которых общая стоимость (общая оценка) всего сырья, используемого для производства продукции, была минимальной при условиях, что суммарная цена сырья всех трѐх видов, используемого на производство единицы продукции данного типа (A, B, C) была не меньше цены единицы продукции данного типа. Можно дать экономическую интерпретацию любой двойственной задачи линейной оптимизации. 6. ТРАНСПОРТНАЯ ЗАДАЧА 6.1 Постановка задачи, свойства и нахождение опорного решения Транспортная задача (ТЗ) формулируется следующим образом. В m пунктах отправления A1, A2,…, Am имеется однородных груз в количествах соответственно a1 , a 2 ,...,a m единиц. Этот груз необходимо доставить в n пунктов назначения B1 , B2 ,, Bn , которые подали заявки в количествах 53 b1 , b2 ,...,bn соответственно. Известна стоимость cij перевозки единицы груза из i-го i 1, m пункта отправления в j-й j 1, n пункт назначения. Требуется составить план перевозок при котором заявки были бы удовлетворены и при этом общая стоимость всех перевозок была бы минимальной. Такую задачу называют транспортной задачей по критерию стоимости. Обозначим через x ij i 1, m ; j 1, n количество единиц груза которое необходимо доставить из i-го пункта отправления в j-й пунк назначения. x11 x12  x1n x21 x22  x2 n .     xm1 xm 2  xmn называется матрицей перевозок или план перевозок. Предполагается, что i, j : xij 0 . Удельные транспортные издержки (расходы) запишем в виде матрицы тарифов c11 c12  c1n c 21 c 22  c 2 n cij .     c m1 c m 2  c mn Математическая модель ТЗ должна отражать все условия и цели задачи в математической форме. Переменные xij должны удовлетворять ограничениям по запасам и по потребностям. Это можно записать следующим образом: x11 x12  x1n a1, xij x21 x22  x2n a2 , , xm1 xm 2  xmn am, x11 x21  xm1 b1, x12 x22  xm 2 b2 , n xij ai , (i 1, m) . (6.1) j 1 m xij b j , ( j 1, n) . (6.2) , i 1 xnm x2n  xmn bn , Суммарная стоимость z всех перевозок должна быть минимальной 54 a11x11 a12 x12  a1n x1n z a21x21 a22 x22  a2n x2n m n  z aij xij min, (6.3) i 1j 1 am1xm1 am 2 xm 2  amn xmn . План перевозок xij называется допустимым, если он удовлетворяет ограничениям (6.1), (6.2). Допустимый план перевозок называется оптимальным, если он доставляет минимум функции z. Теорема 1. Для того, чтобы ТЗ имела допустимый план, необходимо и достаточно, чтобы выполнялось равенство m i 1 ai n bj . (6.4) j 1 □ Необходимость. Пусть ТЗ имеет допустимый план xij , т. е. выполняются равенства (6.1) и (6.2). Просуммируем все переменные xij по строкам и по столбцам, т. е. по i и по j. Получим m n i 1j 1 xij n bj ; j 1 m n j 1i 1 xij m ai i 1 m i 1 ai n bj . j 1 Очевидно, что суммируются все элементы xij как по строкам, так по столбцам, различие лишь в перестановке этих элементов. Но при перестановке сумма не изменяется, т. е. равенство (6.4) является необходимым условием решения ТЗ. Достаточность. Пусть выполняется условие (6.4). Обозначим m n ai b j ai b j s ; xij , ( i 1, m, j 1, n ). (6.5) s i 1 j 1 Покажем, что это допустимое решение. Поскольку ai 0, bj 0 , то и xij 0 . Покажем, что переменные (6.5) составляют допустимый план. Этот набор неотрицательных чисел будет составлять допустимый план, тогда, когда он удовлетворяет системе ограничений (6.1), (6.2). Просуммируем равенства (6.5) по индексу i. Получим m m ab m bj m i j xij xij ai b j 0 . s s i 1 i 1 i 1 i 1 Аналогично, просуммируем (6.5) по индексу j, получим s>0 55 n xij j 1 n ai b j s j 1 ai s n bj ai j 1 Следовательно, числа xij 0, i 1, n . ai b j удовлетворяют системе (6.1), (6.2) огs раничений задачи, и следовательно являются допустимым планом.  Транспортная задача. называется закрытой, если суммарный объем груза в ПО равен суммарным заявкам ПН, т. е. если выполняется равенство (6.4). Если для ТЗ выполняется одно из условий m ai i 1 n bj j 1 или m ai i 1 n bj , j 1 то транспортная задача называется открытой. Вначале методы решения ТЗ будут рассмотрены для закрытой модели, а затем методы решения открытых задач. Транспортную задачу как линейную задачу можно решать симплексметодом. Но, видно, что особенностью ТЗ является то, что все коэффициенты при переменных (6.1) и (6.2) равны единице, а также, то что система уравнений ограничений из (m+n) уравнений является линейно зависимой. Теорема 2. Ранг матрицы A системы уравнений (6.1)-(6.2) равен r=m+n-1. □ Неизвестные в (6.1)-(6.2) x11x12 ,, x1n x21x22 ,, x2n x31x32 ,, x3n,  , xm1x m2 ,, xmn . Тогда матрица системы ограничений (7.1), (7.2) имеет вид: 1 1  1 0 0  0  0 0  0 0 0  0 1 1  1  0 0  0              0 0  0 0 0  0  1 1  1 A 1 0  0 1 0  0  1 0  0 0 1  0 0 1  0  0 1  0              0 0  1 0 0  1  0 0  1 В каждом столбце матрицы A содержится только два элемента, равных единице, остальные элементы равны 0. При этом, если сложить первые m строк матрицы, получим строку, элементы которой будут 1. Этот же результат получим если сложить последние n строк. Обозначим i-ю строку через Ai , получаем A1 A2  Am Am 1 Am 2  Am n Отсюда следует, что любая строка есть линейная комбинация остальных строк, например 56 A1 Am 1 Am 2  Am n ( A2 Am 2 ... Am ) . Значит не меняя ранга матрицы можно вычеркнуть, например, последнюю строку, r=m+n-1. Таким образом, в системе уравнений ограничений m+n-1 базисных переменных, а остальные свободные. Их число равно k=mn-(m+n-1)=(m n-m)-(n-1)=(m-1)(n-1) Следовательно, для оптимального плана перевозок по крайней мере (m-1)(n-1) значений xij должны бать равны нулю. Допустимый план перевозок будет опорным, если в нем отличны от нуля не более r=m+n-1 базисных перевозок (переменных). Условия транспортной задачи записываются в виде транспортной таблицы (ТТ). Таблица 6.1 ПН ПО А1 A2 … Am заявки bj В1 В2 С 11 C 12 C 21 C 22 … cm1 b1 … C m2 b2 ... … … … … b3 Вn c1n C 2n запасы ai a1 а2 … C mn bn … аm ai bj Решение транспортной задачи ТЗ всегда существует и, поскольку множество значений z , ограничено снизу z 0 , то среди допустимых планов имеется оптимальный. Решение ТЗ начинается с построения базисного плана. Для этого существуют разные способы, из которых простейшим является диагональный метод или метод ―северо-западного угла‖. В этом методе перевозки начинают расставлять с клетки (1, 1) и заканчивают в клетке (m, n), при этом перевозки должны удовлетворять условиям (6.1) и (6.2). Пример 1. Условие ТЗ задано транспортной таблицей 6.2. Полученное решение является не только допустимым, но и опорным. Клетки таблицы, в которых стоят ненулевые перевозки являются базисными, их число Остальные клетки свободные, их r m n 1 8. (m 1)(n 1) 3 4 12 . Таблица 6.2 57 ПН B1 ПО B2 10 B3 8 A1 27 B4 5 запасы B5 6 аi 9 3 48 18 6 A2 7 8 A3 7 7 A4 заявки 18 bj 5 27 8 30 10 9 4 42 6 5 8 7 12 6 6 8 12 30 27 20 20 26 125=125 z 10 18 8 27 5 3 8 30 10 9 8 12 7 6 8 20 1039 Это есть стоимость плана перевозок. Пример 2. Дана транспортная таблица 6.3. Составить опорный план перевозок. Таблица 6.3 ПН запасы ПВ В2 В3 В4 В5 1 ai ПО A1 A2 A3 A4 заявки bj 10 10 20- 10 25 20- 10 10 20 35 20 20+ 30 25 2095=95 Применив метод северо-западного угла, получим опорный план. Но в нем шесть отличных от нуля переменных, две базисных переменных оказались равными нулю. Это, так называемый, вырожденный случай. Такие случаи вырождения возможны не только при составлении данного плана, но и при его оптимизации. В дальнейшем удобно всегда иметь в транспортной таблице m n 1 базисных ненулевых клеток. Для этого достаточно в нужных местах изменить запасы или заявки на величину , а после нахождения оптимального решения положить =0. Например, такое изменение сделано в таблице 6.3. Другим методом построения опорного решения является метод минимального элемента. При построении опорного плана, очевидно, выбор 58 пунктов назначения и отправления целесообразно производить, ориентируясь на тарифы перевозок, а именно, на каждом шаге следует выбирать какуюлибо клетку, отвечающую минимальному тарифу (если таких клеток несколько, то следует выбирать любую), и рассмотреть пункты назначения и отправления, соответствующие выбранной клетке. Сущность метода минимального элемента и состоит в выборе клетки с минимальным тарифом. Этот метод, как правило, позволяет найти опорный план транспортной задачи, при котором общая стоимость перевозок груза меньше, чем общая стоимость перевозок при плане, полученном методом северо-западного угла. Таблица 6.4 ПН П В1 В2 В3 ЗАПА СЫ В4 О ai A1 7 A2 4 8 1 2 9 8 160 5 120 20 9 A3 ЗАЯВ КИ 120 2 3 6 50 30 90 50 190 110 160 140 170 470 bj Стоимость перевозок по методу минимальных перевозок будет равна z 160 120 4 20 8 50 2 30 3 90 6 1530 6.2 Распределительный метод решения ТЗ Рассмотрим транспортную таблицу (ТТ), имеющую m строк и n столбцов. Определение 1. Циклом в ТТ называется ломаная с вершинами в клетках и звеньями, лежащими вдоль строк и вдоль столбцов таблицы, удовлетворяющая условиями: 1) ломаная должна быть связанной, т. е. из любой ее вершин можно попасть в любую другую вершину по звеньям ломаной; 59 2) в каждой вершине ломаной встречаются ровно два звена, причем одно из них располагается по строке, а другое по столбцу. Отметим, что циклом может быть и самопересекающаяся ломаная, но точки ее самопересечений не могут в силу условия второго оказаться ее вершинами. В табл. 6.5. построены следующие циклы: Ц1 : цикл с вершинами в клетках (1,1), (1,2), (3,2), (3,1). Ц 2 : цикл с вершинами в клетках (2,3), (2,5), (4,5), (4,4), (3,4), (3,6), (5,6), (5,3). Таблица 6.5 ПН B1 ПО A1 + A2 A3 A4 A5 заявки bj B2 B3 B4 B5 B6 C 12 C 13 C 14 C 15 C 16 C 21 C 22 C 23 C 24 C 25 C 26 C 31 + C =+ 32 C 33 C 34 C 35 C 36 C 42 + C 43=+ C 44 С 11 C 41 C 51 b1 - + =+ C 52 b2 C 53 =+ b3 =+ C 54 b4 =+ C 45 + =+ =+ C 55 b5 C 46 C 56 + =+b запасы ai a1 a2 a3 a4 a5 6 Лемма 1. Пусть в ТТ из m строк и n столбцов отмечено m n клеток mn m n . Тогда существует цикл, вершины которого лежат в отмеченных клетках (может быть и не во всех). □ Если m n 4 , то m n m n , что не удовлетворяет условию. Если m n 4 и m n m n , то m 2, n 2 . В этом случае все клетки должны быть отмечены и существование цикла очевидно. 60 По методу полной математической индукции, предположим, что лемма справедлива для ТТ с суммой числа строк и столбцов m n 1, и, следовательно, с m n 1 отмеченными клетками. Докажем, что она будет справедливой и для ТТ с суммой строк и столбцов m n . Рассмотрим отдельно два случая. 1. Допустим, что имеется отмеченная клетка, для которой в строке (или столбце), в которой она расположена, нет других отмеченных клеток. Удалим эту строку. В оставшейся ТТ сумма строк и столбцов m n 1. По предположению, в такой ТТ имеется цикл с вершинами в отмеченных клетках, но тогда цикл имеется и в исходной ТТ. 2. Какую бы отмеченную из числа m n клеток не взять, и в строке, и в столбце, где она расположена, имеются другие отмеченные клетки. Начав с какой-нибудь из отмеченных клеток и двигаясь по столбцу, можно перейти от нее к другой отмеченной клетке, от этой второй клетки, двигаясь по строке, можно перейти к третьей отмеченной клетке и т.д. Так как число строк и столбцов конечно, то через конечное число шагов выйдем на строку, где уже были. Пусть это i -я строка, и здесь мы вышли на отмеченную клетку. Соединив стрелкой эти клетки, мы замкнем цикл. Лемма доказана. Лемма 2. Число вершин в каждом цикле четно. Доказательство очевидно из определения. Зададим в транспортной таблице произвольный цикл. Некоторую вершину примем за начальную и припишем ей знак +. Затем, обходя цикл в каком-либо направлении, приписываем вершинам поочередно знаки + и - (см. табл. 6.5). Цикл с указанной расстановкой знаков называется означенным. Перенести (―перебросить‖) какое-то количество k единиц груза по означенному циклу – это значит увеличить перевозки (переменные) стоящие в положительных вершинах цикла на k единиц, а перевозки стоящие в отрицательных вершинах уменьшить на то же количество. Очевидно, что при переносе любого числа единиц по означенному циклу равновесие между запасами и заявками не меняется, следовательно при любом циклическом переносе, оставляющем перевозки неотрицательными, допустимый план остается допустимым. Определение 2. Ценой означенного цикла называется увеличение стоимости перевозок при перемещении по этому циклу одной единицы груза. Цена цикла 1 c11 c12 c32 c31 (см. табл. 6.5) 2 c34 c36 c56 c53 c23 c25 c45 c44 . 61 Если цена цикла и перебрасывается по циклу k единиц груза, то k . Очевидно, что для опстоимость перевозок z изменяется на величину тимизации плана имеет смысл перемещать перевозки только по циклам с отрицательной ценой. Метод последовательного улучшения плана перевозок и состоит в том, что в таблице означаются только циклы с отрицательной ценой, по ним перемещаются перевозки и план улучшается до тех пор пока циклов с отрицательной ценой не останется. При улучшении плана циклическими переносами пользуются идеями, заимствованиями из симплекс-метода: при каждом шаге (цикле) заменяют одну свободную переменную на базисную, то есть заполняют одну свободную клетку и взамен освобождают одну из базисных. Лемма 3. В транспортной таблице ТТ не существует цикла, все вершины которого находятся в базисных клетках. □ Дадим свободным переменным какие-либо определенные числовые значения. Тогда и базисные переменные примут определенные значения. Пусть лемма неверна, т.е. существует цикл, все вершины которого находятся в базисных клетках. Перенесем по этому циклу какое-либо количество единиц груза. При этом снова получим решение системы (6.1) (6.2). Поскольку вершины цикла находятся в базисных клетках, при переносе изменились только базисные переменные, а ни одна из свободных не изменилась. Но это невозможно. Противоречие. ■ Определение 3. Циклом пересчета данной свободной клетки называется цикл, одна вершина которого находится в этой свободной клетке, а все остальные вершины – в базисных. Свободной клетке всегда будем приписывать знак +. Теорема 3. В транспортной таблице для каждой свободной клетки существует, и притом единственный, цикл пересчета. □ Выделим произвольную свободную клетку ТТ. Вместе с ней рассмотрим все базисные клетки, их m+n-1. Следовательно, вместе со свободной клеткой выделим m+n клеток. Согласно лемме 1 найдется цикл с вершинами в выделенных клетках, причем в силу леммы 3 среди этих вершин содержится рассматриваемая свободная клетка. Для доказательства единственности предположим противное. Пусть для данной свободной клетки существуют два различных цикла пересчета. Простым перебором возникающих при этом возможностей 62 свободная клетка можно показать, что в этом случае всегда можно выделить цикл с вершинами только в базисных клетках, что противоречит лемме 3.  Таким образом, можно сделать вывод: если цена цикла пересчета с плюсом в свободной клетке отрицательна, то план можно улучшить перемещением перевозок по этому циклу. Количество единиц груза k , которое можно переместить, определяется минимальным значением перевозок, стоящих в отрицательных вершинах цикла (если перемещать большее число единиц груза, возникнут отрицательные перевозки). Пример 3. Найти оптимальный план ТЗ приведенной в таблице 6.6. Таблица 6.6 63 ПН Ы П В1 В2 В3 В4 запасы ai О А 10 7 22 A1 А A2 А A3 - 8 9 5 + 25 8 + 6 31 6 5 4 23 7 48 6 - 7 38 18 20 41 20 З ЗАЯВ КИ 22 34 117=117 bj Составляем опорный план методом северо-западного угла. План невырожденный. Число положительных базисных переменных r=m+n-1=3+4-1=6. Стоимость этого плана z=10 22+7 9+6 25+5 23+6 18+7 20=796. Попробуем улучшить план, заняв свободную клетку (2, 4) с минимальной стоимостью c24 4 . Цена цикла с вершиной в точке (2, 4) равна 2 . Перемещаем по этому циклу 20 единиц груза. Полу24 4 7 6 5 чаем новый план (табл. 6.7). Стоимость этого плана z=756. Таблица 6.7 64 ПН З П В1 В2 В3 ЗАПАС Ы В4 О ai А A1 - 10 22 А A2 7 + 6 8 9 5 31 6 25 5 3 4 20 48 - + А 8 7 6 7 38 A3 38 З ЗАЯВК И 22 34 41 20 117=117 bj Теперь возьмем свободную клетку с наименьшей стоимостью c21 = 5. Цикл пересчета этой клетки имеет цену 21 5 10 7 6 4 . Переместим 22 единицы груза. Получим таблицу 6.8. Таблица 6.8 65 ПН П В1 В2 В3 ЗЗАПА СЫ В4 О ai А 10 7 6 8 31 A1 А 5 22 A2 А 31 6 3 8 5 3 7 4 20 6 7 38 A3 48 38 З ЗАЯВ КИ 22 34 41 20 117=117 bj Дальнейший перебор свободных клеток показывает, что цены их циклов пересчета либо положительные, либо нулевые. Таким образом, план Т3 является оптимальным. Решение будет следующим: 0 31 0 0 22 3 3 20 0 0 38 0 zm / n 668 . Пример 4. Найти оптимальный план перевозок транспортной задачи приведенной в транспортной таблице 6.9. Таблица 6.9 ПН В1 В2 В3 ai ПО – 5 10 4 + АA1 40+ 20 20 + 6 4 5 2 АA2 23 3 7 3 6 2 АA3 – 200+ bj 20 20 43 83=83 66 Начальный план получился вырожденный, т.к. число ненулевых базисных переменных равно 4, а должно быть m+n-1=5=r. Чтобы получить невырожденный план изменим запасы в первой и третьей строках. Перемещаем единиц груза по циклу пересчета для клетки (3; 2), в которой 32 = 320 5+4-6= - 4. Получим табл. 6.10. Таблица 6.10 ПН П В1 В2 В3 ai О А- 20 A1 А A2 5 4 7 5 3 20 20 23 6 20- 20- A3 40+ 23 - + 4 20 6 А bj 10 + 43 83=83 Свободная клетка с наименьшей стоимостью (2, 2). Цена цикла пересчета 22 4 5 4 5 2 . Можно переместить единиц груза (так как 23 нельзя, будут отрицательными перевозки). Результат занесем в таблицу 6.11. Таблица 6.11 67 ПН П В1 В2 В3 ai О АA1 4 + 40+ 20+ 6 4 5 - + А 7 20 3 20 23 236 20- 20- A3 bj 5 20 А A2 10 43 83=83 Отрицательную цену имеет цикл пересчета свободной клетки (2, 1) 5 . Переносим по нему 20 единиц груза. 21 6 10 4 5 Таблица 6.12 68 ПН П В1 В2 В3 ai О А 10 5 4 40+ 40+ A1 А 6 4 20 A2 А 5 7 3 6 20- 20- A3 bj 23 3- 20 20 43 83=83 В этой таблице все циклы пересчета имеют неотрицательную стоимость. Поэтому план является оптимальным. Полагая в нем 0 , получим окончательный оптимальный план, причем zmin 4 40 6 20 5 3 3 20 355 . Рассмотренный метод решения транспортной задачи называется распределительным. 6.3 Метод потенциалов решения транспортной задачи Пусть имеется ТЗ с балансовыми условиями (6.1), (6.2): n xij ai , i 1, m, xij bj, j 1, n, j 1 m i 1 причем выполняется условие (6.4) m n ai i 1 bj . j 1 Стоимость перевозки единицы груза из Ai в B j равна cij . Требуется найти план перевозок ( xij ), который удовлетворял бы балансовым условиям 69 (6.1), (6.2) и при этом стоимость всех перевозок z m n cij xij была бы ми- i 1j 1 нимальной. Поставим в соответствие каждому из пунктов отправления Ai некоторое число i , i 1, m и каждому пункту назначения B j искомое число j , j 1, n . Обозначим c~ . Будем называть величину c~ ―псевдостоиi j ij ij мостью‖ перевозок. Числа i и j будем называть платежами. Пусть в базисных клетках (6.5) i j сij . т. е. в базисных клетках псевдостоимость равна стоимости. Теорема 4. Если система (6.1), (6.2) невырождена, то ранг r системы (6.5) равен рангу системы (6.1) (6.2), т. е. r=m+n-1. (Без доказательства). Уравнений в (6.5) всего m+n-1, а число неизвестных i и j равно m+n. Следовательно одну из этих неизвестных можно задать произвольно. После этого из m+n-1 уравнений можно найти все остальные платежи i , j , а по ним вычислить псевдостоимости для каждой свободной клетки. Пример 5. Найти псевдостоимости ТЗ заданой табл. 6.13. Таблица 6.13 ПН П В1 В2 В3 В4 заявки платежи ai i В5 О A1 A2 A3 А4 заявки bj 10 10 8 17 5 8 9 6 4 13 9 7 4 3 7 5 14 12 17 4 25 32 -2 40 -1 20 20 4 5 4 12 10 21 24 117=117 3 2 2 19 22 14 5 5 8 8 6 9 6 5 4 3 14 10 9 41 3 4 12 8 14 8 70 платежи 10 8 6 5 4 j В правом верхнем углу будем писать стоимости перевозок cij , а в левом верхнем – псевдостоимости c~ij . По способу северо-западного угла строим опорный план. Положим 1 1 10 1 2 8 2 2 6 2 3 4 3 3 5 3 4 3 5 1=0, и применяя далее условие (6.5), получим: 1 10 1 2 8 3 6 4 4 5 3 5 4 2 2 1 3 8 4 4. Далее находим псевдостоимости c~ij i j для свободных клеток, например, 1 3 0 6 и заполняем этими значениями левые углы клеток. Теорема 5. Для данной совокупности платежей ( i , j ) суммарная стои4 5 z мость перевозок с псевдостоимостями ~ m n c~ij xij при любом i 1j 1 допустимом плане перевозок ( xij ) сохраняет одно и тоже значение ~ z c const . z □ ~ m n ( i 1j 1 m n i 1 i xij j 1 m n j ) xij i m n i 1 j i 1j 1 m xij j 1 i 1 i xij i ai m n i 1j 1 m i 1 jb j j xij . (6.6) Правая часть не зависит от плана перевозок ( xij ), а зависит только от заказов a i , заявок b j и платежей ( i , j ). Теорема 6. (Достаточное условие оптимальности плана перевозок). Если для всех базисных клеток плана ( xij 0 ) выполняется условие c~ c , i j ij ij а для всех свободных клеток ( xij c~ c , i j ij ij то план является оптимальным. 0 ) выполняется условие 71 □ Пусть ( xij ) – план перевозок, ( i , j ) система платежей, удовлетворяющих условиям теоремы. Стоимость плана будет m n z cij xij . i 1j 1 В этой сумме xij отличны от нуля только слагаемые, соответствующие базисным клеткам, в них c~ij cij . Поэтому m n ~ z z c~ij xij . i 1j 1 Пробуем изменить план ( xij ), заменив его каким-то другим планом ( xij ). Стоимость нового плана будет m n z cij xij . i 1j 1 Значения плана ( xij ) отличны от нуля, вообще говоря, в других клетках, чем плана ( xij ). Некоторые из этих клеток совпадают с прежними – базисными для плана ( xij ), а другие, со свободными для плана ( xij ). В первых по-прежнему а во вторых Поэтому c~ c , c c~ . z m n cij xij i 1j 1 ij m n ij c~ij xij ij ~ z , то есть z ij ~ z . Но по теореме 5 при таком же i 1j 1 z ~ z z , т. е. z z. Следовательно, никаким изнаборе платежей ( i , j ) ~ менением плана ( xij ) его стоимость не может быть уменьшена. Следовательно, план ( xij ), удовлетворяющий условию теоремы 6, является оптимальным. Теорема верна для вырожденного случая.  Таким образом, признаком оптимального плана ( xij ) является выполнение двух условий cij cij для всех базисных клеток , cij cij для всех свободных клеток. План, обладающий таким свойством называется потенциальным, а соответствующие ему платежи ( i , j ) потенциалами пунктов Ai и B j . Теорему 6 можно переформулировать теперь так: всякий потенциальный план является оптимальным. Теорема 7. Какова бы ни была система платежей ( i , j ), удовлетворяющая условию теоремы 6, цена цикла пересчета каждой свободной клетки равна разности между стоимостью cij и псевдостоимостью c~ в данной клетке ij 72 cij c~ij . □ Начнем обходить цикл пересчета выйдя из свободной клетки (i, j) например по столбцу. Сначала мы попадем в базисную клетку (k, j), расположенную в k-й строке и j-м столбце. Двигаясь по k-й строке попадем в клетку (k, l), расположенную в l-м столбце и той же строке k и т. д. Наконец, выйдя из некоторой базисной клетки (i, p), лежащей в p-м столбце мы вернемся и притом по i-й строке в свободную клетку (i, j). Таким образом, вершинами цикла пересчета служат следующие клетки (i, j) (k, j) (k, l) (s, l) (s, t) … (r, p) (i, p) Подсчитаем цену цикла пересчета (6.9) ij cij ckj ckl csl cst  crp cip . Стоимость cip войдет в сумму со знаком ―-‖, т. к. число вершин в цикле четно. Все клетки, кроме (i, j) в цикле пересчета являются базисными, поэтому в сумме (6.9) все стоимости, кроме cij равны соответствующим псевдостоимостям. В силу этого, ij cij ( k j) ( k l) ( s l) ( s t)  ( r p) ( i p) Раскрывая скобки, получим ~ ij cij ( i j ) cij cij .  Теперь ясно, что при использовании метода потенциалов для решения ТЗ отпадает наиболее трудоемкий элемент распределительного метода: поиск циклов с отрицательной ценой. ij Алгоритм решения ТЗ методом потенциалов: 1. Составить опорный план перевозок, в котором m+n-1 базисная клетка. 2. Определить для этого плана платежи, исходя из условия, что в каждой базисной клетке псевдостоимости равны стоимостям, причем один из платежей можно выбрать произвольно, например, положив равными нулю. 3. Подсчитать псевдостоимости всех свободных клеток ( c~ij i j ). Если окажется, что все они не превышают стоимостей, то имеющийся план оптимален. 4. Если хотя бы в одной свободной клетке псевдостоимость c~ij больше стоимости cij , следует приступить к улучшению плана путем переноса перевозок по циклу пересчета для любой свободной клетки с отрицательной ценой (для которой c~ij cij ). 5. После этого заново пересчитываются платежи и псевдостоимости. Если план все еще не оптимальный, процедура улучшения продолжается до тех пор, пока не будет найден оптимальный план. Пример 6. Решить ТЗ, данную в табл. 7.13 (Пример 5). 73 Полученный с помощью метода северо-западного угла план не оптимален. Попробуем улучшить его, переводя в базисную одну из свободных клеток, для которых c~ij cij , например клетку (2, 1). Цена цикла пересчета для этой клетки с с~ 5 8 3 . Перенесем по циклу 13 единиц груза 21 21 21 уменьшив тем самым стоимость плана на 39 и получим след. таблицу. Таблица 6.14 ПН П В1 В2 В3 В4 заявки платежи ai i В5 О 10 А A1 10 8 - 4 5 A2 А 8 9 21 5 6 А 8 5 7 + 3 4 13 + 9 9 8 4 3 3 8 6 7 A3 5 4 22 + 9 8 14А 11 10 4 3 14 - 10 8 9 8 8 20 заявки bj платежи 17 21 41 14 24 10 8 9 8 7 j Вычислим платежи. Положив снова 1 32 -5 40 -4 20 1 3 4 A4 2 19 4 5 25 117 74 1 1 10 1 2 8 2 1 5 2 3 4 3 3 5 3 4 3 4 1 2 8 3 9 4 4 8 5 3 5 7 5 8 1 2 3 4 5 4 1. Располагаем в таблице 3 псевдостоимости. Еще есть свободная клетка, в которой c~ij cij , например (1, 4). Перенос по циклу пересчета для этой клетки приводит к табл. 6.15. Таблица 6.15 ПН П В1 В2 В3 В4 заявки платежи ai i В5 О A1 A2 10 8 8 8 9 21 5 17 9A -5 5 6 4 6 7 6 5 15 заявки bj платежи 14 5 25 32 –3 40 –2 20 20 3 117 +4 3 3 8 2 5 4 4 3 3 26 + 11 10 6 5 4 A3 А4 7 6 10 10 8 11 8 4 9 8 8 17 21 41 14 24 8 8 7 6 5 j Снова вычислим платежи и псевдостоимости. Перенеся по циклу, соответствующему свободной клетке (4, 3) 20 единиц груза, получает новый план с новыми платежами и псевдостоимостями. Таблица 6.16 75 ПН П В1 В2 В3 В4 заявки платежи ai i В5 О 10А 8 8 8 7 9 6 21 A1 A2 5А 5 6 17 5 4 9А 6 5 4 3 3 8 5 4 6 9 10 9 8 4 3 10 8 8 платежи 32 –3 40 –2 20 1 3 24 7 8 6 20 заявки bj 25 2 15 6 7 14 5 4 A3 А4 6 5 17 21 41 14 24 8 8 7 6 5 117=117 j В этой таблице псевдостоимости не превосходят соответствующих стоимостей, следовательно, этот план оптимален. zmin = 8 21+6 4+5 17+4 15+5 6+4 10+3 24+8 20 = 639. 6. 4. Двойственная транспортная задача Транспортной задаче, как всякой задаче линейной оптимизации можно поставить двойственную задачу (ДЗ) по общему правилу построения ДЗ. Запишем подробнее основную ТЗ m n z aij xij min, i 1j 1 cij c11 c 21  c m1 c12 c 22  cm2  c1n  c 2n .    c mn 76 x11 x12  x1n x21 x22  x2 n a1, a2 , , xm1 xm 2  xmn am , x11 x21  xm1 b1, x12 x22  xm 2 b2 , , xnm x2n  xmn bn , xij 0 . Каждому из (m+n) ограничений поставим в соответствие двойственную переменную и получим вектор u T (u1 um , um 1 um n ) . Согласно правилу построения получим следующую двойственную задачу: w a1u1  amum b1um 1  bnum n w m i 1 ai ui n b j um j 1 Ограничения: u1 um 1 c11, u1 um 2 c12 ,  u1 um n c1n ,  um um 1 cm1, um um 2 j max ui um j cij , i 1, m, j 1, n, cm 2 ,  um um 1 cmn , u j 0 или u j 0 , um j 0 или um j 0 . Из постановки двойственной транспортной задачи ясно, что двойственные переменные являются платежами i , j . То есть ui i , um j j , u ( 1 m , 1 n ) . Двойственную задачу можно записать в виде: 77 m n w ai bj i i 1 i max, j j 1 сij , j i 1, m, j 1, n. Для двойственной транспортной задачи справедливы все теоремы теории двойственности. В частности, если основная ТЗ или двойственная ТЗ имеют решения w или: xij X и u , то z m n cij xij i 1j 1 m i 1 m i ai i 1 jb j , где i , j - потенциалы. Смысл этого равенства состоит в том, что суммарные транспортные расходы при оптимальном плане перевозок равны суммарной оценке стоимости продукции при полном удовлетворении спроса. Смысл достаточных условий оптимальности (Теорема 6) состоит в следующем. В свободных клетках i xij 0 , т. е. Перевозки из Ai в j cij B j не рентабельные. В базисных клетках i j cij xij 0 , и, следовательно, перевозки рентабельные. Величину i + j можно интерпретировать как приращение ценности единицы продукции при перевозке из Ai в B j . i – ценность единицы продукции, содержащейся в Ai ; j – ценность единицы продукции, которая окажется после перевозки в B j . 6.5. Особые случаи решения ТЗ 1. Открытая ТЗ. Модель ТЗ называется открытой, если выполняется одно из условий m ai i 1 n bj или j 1 m ai i 1 n bj . j 1 Для решения открытой ТЗ необходимо преобразовать ее в закрытую, то есть, чтобы сумма запасов и заявок была равна. При выполнении первого условия, когда запасов больше чем заявок надо ввести (n+1)-й пункт назначения Bn 1 , то есть в ТТ введется дополнительный столбец. Спрос фиктивного потребителя полагается равным небалансу, то есть bn m 1 i 1 ai n bj , j 1 78 а все тарифы одинаковыми, чаще всего равными нулю, то есть ci,n 1 0 i 1, m . При выполнении второго условия, когда заявки больше, чем запасы, аналогично, вводится фиктивный поставщик Am 1 , запас груза у которого равен ai n 1 bj j 1 m i 1 ai , а тарифы дополнительной строки ТТ полагают равными нулю, то есть cm 1, j 0 j 1, n . Пример. В трех хранилищах A1 , A2 , A3 имеется соответственно 70 т, 90 т и 50 т топлива. Требуется спланировать перевозку топлива четырем потребителям спрос которых соответственно равен 50, 70, 40 и 40 т., так, чтобы затраты на транспортировку были минимальными. Стоимость одной тонны перевозки из Ai в B j задается матрицей стои- 5 2 3 6 мостей. (cij ) 4 3 5 7 . 2 4 1 5 Решение. ai 210 и b j 200 . Составим ТТ и, поскольку, запасы топлива в хранилищах превышают спрос потребителей, то вводится фиктивный потребитель bs ai b j 210 200 10 Все затраты для фиктивного потребителя положим равной нулю: cis 0 . Построим ТТ Найдем опорный план методом минимального элемента. 79 Таблица 6.17 ПО П В1 В2 В3 В4 В5 ai i Н 3А 5 2 2 2 3 6 6 0 60 A1 4 A2 4 3 40 A3 bj - 10 3 3 5 7 4 1 1 5- 5 -1 40 3 2 2 + 90 1 50 –1 40 70 70 40 50 j 7 1 10 2+А 2 1 10 40 10 6 210 Наименьшими являются нулевые тарифы для клеток (1, 5), (2, 5); (3, 5), то загрузим сначала клетку (1, 5) 10; Второй загрузим клетку (3, 3) 40, далее загружаем клетки (1, 2), (3, 1), (2, 2), (2, 1), (2, 4) 60, 10, 10, 40, 10. Получили невырожденный план: m+n-1=3+5-1=7. Далее решаем методом потенциалов. Из условия i j сij получим систему и найдем платежи. 1 2 2 2 0 1 0 1 5 5 2 1 4 1 4 1 3 2 2 3 2 4 7 3 2 2 3 3 1 2 3 cij 3 2 1 2 3 4 7 1 6 3 1 ( 1) 2 cij i j. 1 Найдем оценки свободных клеток: sij сij cij c~ij . Среди i j 1. свободных клеток одна отрицательная имеет отрицательную оценку s25 Следовательно, план перевозок можно улучшить за счет загрузки клетки (2, 80 5). Наибольшее количество топлива в отрицательных клетках (вершинах) цикла равно 10. Перемещаем. Получим табл. 6.18. Таблица 6.18 ПН П В1 В2 В3 В4 В5 ai i О 3 A1 5 2 2 2 3 6 6 0 704 A2 4 3 3 3 5 7 7 0 40 2 A3 40 2 1 4 1 10 bj j 5 -2 70 3 40 3 40 7 90 50 –2 40 50 4 –1 10 1 5 70- 10 210 Полученный план оказался вырожденным. m+n-1=6. Изменим запасы a1 на , a2 на . 1 1 2 2 1 1 4 2 1 4 2 0 2 2 3 2 3 2 4 7 4 7 2 5 5 3 1 2 3 2 1 3 3. Найдем оценки свободных клеток: ij cij c~ij . Они все неотрицательны ij 0 . Поскольку среди оценок имеются равные нулю 14 0 ; 34 0 , то можно получить новые планы, но значение целевой функции не изменится. Это случай бесчисленного множества оптимальных планов. Полагая 0, получаем 0 70 0 0 3 xij 3 40 10 0 40 40 0 81 z 2 70 4 40 7 40 2 10 1 40 640 . 10 т. груза на топлива в хранилище A2 осталось невостребованным. 2. При некоторых реальных условиях перевозки из некоторого пункта Ai в B j не могут быть осуществлены. Для нахождения оптимальных планов в этих случаях предполагают, что тариф перевозки единицы груза cij является сколь угодно большим, то есть cij M 1. При этих условиях находят решения известными методами новой ТЗ. Такой подход решения ТЗ называется запрещением перевозок или блокированием соответствующей клетки ТЗ. 3. В отдельных ТЗ дополнительным условием является обеспечение перевозки по заданным маршрутам определенного количества груза. Пусть, например, из Ai в B j требуется обязательно перевезти d ij единиц груза. Тогда в клетку ТЗ, стоящей на пересечении i-й строки и j-го столбца записывают это число d ij , и в дальнейшем эту клетку считают свободной со сколь угодно большим тарифом перевозок M . Находят оптимальный план новой задачи, который является оптимальным планом исходной задачи. 4. В некоторых случаях необходимо найти решение ТЗ, при котором из пункта Ai в B j должно быть завезено не менее заданного количества груза d ij . Для определения оптимального плана такой ТЗ считают, что запасы пункта Ai и потребности пункта B j меньше фактических на d ij единиц. Определяют оптимальный план новой ТЗ, из которого исходит оптимальный план исходной задачи. 5. В некоторых ТЗ требуется найти такое решение, при котором из пункта Ai в пункт назначения B j должно быть завезено не более заданного количества груза d ij , т. е. xij dij (*). Такую задачу можно решить следующим образом. В исходной ТЗ для каждого j-го ограничения предусматривают дополнительный столбец, т. е вводят дополнительный пункт назначения. В этом столбце записывают те же тарифы, что и в столбце B j , за исключения тарифа, находящегося в i-й строке. В дополнительном столбце в этой строке тариф считается сколь угодно большим M. При этом потребности пункта B j считают равными d ij , а потребности вновь введенного пункта назначения полагают равными b j d ij . Задачу решают методом потенциалов и находят оптимальный план или устанавливается ее неразрешимость. Заметим, что исходная ТЗ разрешима лишь тогда, когда для нее существует хотя бы один опорный план. Замечание. Методами ТЗ могут быть решены и другие экономические задачи, не связанные с транспортной. 82 7. НЕЛИНЕЙНАЯ ОПТИМИЗАЦИЯ 7.1. Безусловная оптимизация. Условия оптимальности Пусть f ( x ) f ( x1 , x 2 ,...,x n ) – функция, заданная в R n . Требуется решить задачу: найти точку x , где f (x) достигает min (max). Дальше всюду будем искать min f ( x) , т.к. max f ( x) min( f ( x)) . По определению точки минимума надо найти такую точку ( x* )T ( x1* , x2* ,...,xn* ) R n , что x R n : f ( x* ) f ( x) (8.1), т.е. найти точку глобального минимума функции f (x ) в R n . Если выполняется строгое неравенство f ( x * ) f ( x ) x x* , то x * – единственная точка глобального минимума. Существуют случаи и неединственного глобального минимума (рис.7.2). f ( x) f ( x) x Рис. 7.1 x x x x x Рис. 7.2 Задача нахождения min f ( x) x R n называется задачей нелинейной оптимизации без ограничений. Для f ( x) многих оптимизационных задач без ограничений основные известные методы не дают возможности определить глобальный минимум и находят локальный минимум, т.е. точки, удовлетворяющие условию (7.1) только в некоторой окрестности точки x * (рис.7.3). x x x x На рис. 7.3. обозначены: Рис. 7.3 x * – точка глобального минимума; x , x – точки локального минимума. 83 Пусть функция f (x ) дифференцируема в некоторой точке x Rn. Вектор f f f T f ( x) ; ;...; f x1 , f x2 ,..., f xn , x1 x 2 xn называется градиентом функции f (x ) в точке x . Если f ( x ) f ( x1 , x 2 ,...,x n ) дважды дифференцируема в точке n x R , то матрица вторых производных f (x) называется матрицей Гессе, а соответствующий ей определитель – гессианом. 2 2 2 f f f ... x1 x 2 x1 x n x12 2 . f ( x) ... ... ... ... 2 2 2 f f f ... 2 x n x1 xn x2 xn Матрица Гессе является матрицей второго дифференциала, n n 2 f ( x) 2 d f ( x) dxi dx j как квадратичной формы относительно dxi . x x i j i 1j 1 Действительно, если обозначить 2 f ( x) aij ; dxi li ; dx j l j , xi x j то d 2 f ( x) n aij li l j , i, j 1 или d2 f lT 2 f ( x) l . lT (l1, l2...ln ) Квадратичная форма называется положительно определѐнной, если l R n : l T 2 f l 0 . Теорема 1. (Необходимое условие оптимальности). Пусть функция f (x ) непрерывна и имеет непрерывные частные производные первого и второго порядка в окрестности точки x * . Для того, чтобы точка x * была точкой минимума (глобального или локального) функции f (x ) , необходимо, чтобы выполнялись условия: а) f ( x * ) 0 (условие стационарности); 84 б) матрица Гессе 2 f ( x * ) является положительно полуопре- делѐнной матрицей, т.е. xT 2 f ( x* ) x 0, x Rn . □ Пусть x * – точка локального минимума функции f (x ) . Поскольку f (x ) дважды непрерывно дифференцируема, то еѐ можно разложить по формуле Тейлора в окрестности точки x * : n f (x* ) 1 n n 2 f (x* ) * * f (x x) f ( x ) xi xi x j xi 2! i 1 j 1 x i x j i 1 где, xT x1, x2... xn x x1 , x2 ...xn x12 x22 ... x x* x, 2 Или в векторном виде: 1 T f ( x) f ( x* ) T f ( x* ) x x 2 где O( 2 ) 0 при x ), xn2 . 2 f ( x* ) x O( 2 ), (7.2) 0. x* гда, если выбрать точку x f (x* ) 2 , f ( x* ) 0 . То- Предположим, что не выполняется условие а), т.е. T o( T f (x* ) f (x* ) , x 2 0 , f ( x) f ( x * ) где f ( x* ) T f ( x* ) 0 , то 2 0( ) . Следовательно, из (7.2) следует, что f ( x ) f ( x * ) , т.е. x * не является точкой минимума (противоречие). Теперь пусть не выполняется б), т.е. матрица Гессе 2 f ( x * ) не является положительно полуопределѐнной. В силу того, что условие а) является необходимым, (7.2) можно записать 1 T 2 f ( x) f ( x * ) x f ( x * ) x o( 2 ) . 2 По предположению, существует вектор d R n : d T 2 f ( x* )d 0 . Если выбрать x x * d для достаточно малого 0 получим, 1 2 T 2 d f ( x * )d o( 2 ) . что f ( x ) f ( x * ) 2 * x x x d . Т.к. d T 2 f ( x * )d 0 , то f ( x ) f ( x * ) , т.е. x * не является точкой минимума (противоречие). ■ Заметим, что а) и б) являются необходимыми условиями минимума. Тогда точка x * , удовлетворяющая условию а), называется стационарной точкой функции f (x ) . Стационарность не является достаточным условием локальной оптимальности. Например, для функции одной пе- 85 ременной (рис.7.4) касательная в точке x * параллельна оси Ох, т.е. df ( x * ) 0 . Однако точки минимума нет. dx f x Рис. 7.4 Теорема 2. (Достаточное условие оптимальности). Пусть функция f (x ) дважды непрерывно дифференцируема в окрестности точки x * , а x * - стационарная точка. Тогда, если матрица Гессе 2 f ( x * ) в точке x * является положительно определѐнной матрицей, т.е. выполняется условие x Rn , x 0 : xT 2 f ( x*) x 0 , то x * является точкой ло- кального минимума функции f (x ) . □ Рассмотрим точку x * , которая удовлетворяет условиям теоремы 1, x * –- стационарная точка. Т. к. f (x ) дважды непрерывно дифференцируема, запишем формулу Тейлора в окрестности точки x 1 T 2 f ( x) f ( x * ) x f ( x * ) x o( 2 ) . 2 2 2 x1 x 2 ... x n2 , x O (x* ) , Пусть d T xT x1 x1* , x 2 x 2* ,...,x n (d1, d 2 ,..., d n ) – произвольное единичное приращение (перемещение) из точки x * , причѐм d x x* x n* . d12 d; d 22 1. ... d n2 d x ; x 2 Тогда f ( x * d) * f (x ) dT 2 o(1) – бесконечно малое, т.е. o(1) 2 f (x* ) d 2 0 при x 0. o(1) ; x* d. 86 По условию матрица Гессе 2 f ( x * ) положительно определѐнная, то d T f (x* ) d 0 . Тогда, если – достаточно малое число, то, отбрасывая второе и третье слагаемое в правой части, получим f ( x * d ) f (x* ), f ( x) f (x* ) x O (x* ) , т.е. x * – точка локального минимума. ■ Ясно, что на практике теоремой 2 можно пользоваться, если f (x ) дважды непрерывно дифференцируема. Это очень большое требование к f (x) . Оно не всегда выполняется. Если f (x) непрерывна и имеет непрерывные частные производные первого порядка, то очевидный путь нахождения точки минимума функции f (x ) состоит в нахождении стационарных точек, т.е. проверки условия f (x) 0 , которое равносильно f ( x) 0 , i 1, n . Решение этой системы, как прависистеме уравнений xi ло, можно найти только численно. Однако существует класс достаточно распространѐнных функций, для которых необходимое условие теоремы 1 f (x) 0 является и достаточным, а также все локальные минимумы такой функции f (x ) являются глобальными. Это выпуклые функции. 7.2. Выпуклые функции и их свойства Пусть X Rn – выпуклое множество (и дальше везде). Определение 1. Функция f (x ) , определѐнная на выпуклом множестве X Rn , называется выпуклой, если x , x ется неравенство f x 1 x f (x ) 1 f (x ) . X 0,1 выполня- (8. 3) Если выполняется строгое неравенство, то f (x ) называется строго выпуклой функцией. В одномерном случае ясна геометрическая интерпретация этого понятия (рис.7.5). 87 f ( x) z f (x ) 1 f (x ) f (x ) f ( x) f (x ) z x x f ( x, y ) x (1 x x x x2 y2 a2 b2 )x Рис. 7.5 Рис. 7.6 Выпуклая функция на отрезке x , x не может принимать значение больше, чем линейная функция, интерполирующая значения f (x ) и f (x ) . Из определения 1 следует, что свойство выпуклости функции обязательно предполагает выпуклость множества, на котором задана функция, т.к. x (1 )x X. Свойства выпуклых функций Теорема 3. Если функции f 1 ( x ), f 2 ( x ),..., f k ( x ) выпуклые на некотором выпуклом множестве X , то выпуклой будет на этом множестве и их неотрицательная линейная комбинация, т.е. функция k f ( x) i fi ( x), i 0, i 1, k . i 1 В частности будет выпуклой сумма двух функций. □ x ,x X и x 0,1 имеем: k f( x (1 )x ) i i 1 x0 k fi ( x (1 )x ) i f i ( x ) (1 ) f i (x ) i 1 f ( x ) (1 ) f ( x ) . Функция f (x ) по определению является выпуклой. ■ Теорема 4. Если g (x ) – выпуклая функция x 0 , то множество решений системы g ( x) b, (7.4) x 0, тоже есть выпуклое множество. □ Проверим, что вместе с решениями x и x системы (7.4) точка x 1 x , 0,1 также принадлежит множеству решений системы 88 0 . Остаѐтся показать, что g ( x 0 ) b . Т. к. g ( x ) g( x ) b . g ( x ) b , то 0,1 : 1 g( x ) 1 b По свойству выпуклости функции g (x ) имеем: (4). Ясно, что x 0 b и g ( x 0 ) g ( x (1 ) x ) g ( x ) (1 ) g ( x ) b (1 )b b g (x0 ) b ■ Справедлива теорема 7.5, которую оставим без доказательства. Теорема 5. Всякая выпуклая функция f (x ) , определѐнная на выпуклом множестве X , непрерывна в любой внутренней точке этого множества. Теорема 6. (Необходимое и достаточное условие выпуклости). Для того чтобы функция f (x ) , дифференцируемая во всех внутренних точках множества X , была выпуклой, необходимо и достаточно, чтобы x , x X выполнялось условие f (x ) T f (x ) f (x ) x x , (7.5) f (x ) f (x ) f (x ) , ,..., . x1 x2 xn □ Необходимость. Пусть f (x ) – дифференцируемая на множестве X и выпуклая функция. Тогда справедливо неравенство (7.3): f (1 )x x 1 f (x ) f (x ) f (x ) f (x ) f (x ) , откуда f x (x x) f (x ) f (x ) f (x ) . где T f (x ) (7.6) Запишем формулу Тейлора с точностью o( ) , т.е. при n 1 : f ( x0 f ( x0 ) x) Полагая x 0 x, T x f ( x0 x x) x, 0 1. x , будем иметь: T f x (x x ) f (x ) f x Подставляя в (7.6), получим: T f x (x x ) x x (x f (x ) x) x x . f (x ) . Т. к. это верно 0 , в силу непрерывно0,1 , переходя к пределу при сти f (x ) получим неравенство (7.5). Необходимость доказана. Достаточность. Пусть выполняется условие (7.6). Возьмѐм x x 1 x , 0,1 . x x 1 x x , x x x x . 89 Тогда в силу (7.6) можно записать: T f ( x ) f ( x) f ( x) x x , T f ( x ) f ( x) f ( x) x x . Умножим первое равенство на , а второе на 1 f ( x ) f ( x) 1 f ( x ) f ( x) f ( x) T T f ( x) x T f ( x) x T f ( x) x и сложим: T f ( x) x f ( x) x T f ( x) x T f ( x) x x x x 0, (7.3). ■ f (x ) 1 f ( x ) f ( x) Теорема 7. Если функция f (x ) дважды непрерывно дифференцируема на выпуклом множестве X , то утверждения а), б), в) эквивалентны: а) f (x ) – выпуклая на множестве X функция; б) в) x ,x x X : f (x ) T f (x ) X матрица Гессе 2 f (x ) x x ; f ( x ) в точке x – положительно по- 2 луопределѐнная матрица, т.е. x X : xT f ( x) x 0 (второй дифференциал является положительно полуопределѐнной квадратичной формой). □ Утверждения а) и б) эквивалентны по теореме 7.6. Для доказательства теоремы надо доказать, что утверждение в) эквивалентно либо а), либо б). Докажем, что в) эквивалентно б). x . ТоПусть x , x X - произвольные точки. Обозначим x x x x . Запишем формулу Тейлора в точке x с точностью до гда x малых второго порядка, т.е. o( 2 n 2 ), x i2 : i 1 f (x ) f (x ) 1 df ( x ) 1! 1 2 d f (x ) 2! или f (x ) f (x ) T f (x ) (x x) 1 (x 2 x) 2 f (x ) (x x ). (7.7) Если выполняется б), то xT 2 форма положительная полуопределенная. Верно и наоборот. ■ f ( x) x 0 , т.е. квадратичная 90 Следствие. Квадратичная форма f ( x) n n aij xi x j xT Ax как i 1j 1 функция n переменных является выпуклой функцией, если она положительно полуопределѐнная. Следует из (в), т.к. 2 f ( x ) 2 a ij 2 A . Теорема 8. Выпуклая дифференцируемая функция f (x ) , определѐнная на выпуклом множестве X , достигает своего глобального минимума в каждой точке x , где f (x) 0 . Т.е. если f (x) 0 , то x – точка глобального минимума функции f (x) (глобальный оптимум). □ Пусть в точке x 0 X: f ( x 0 ) 0 . Покажем, что x 0 - точка гло- бального минимума. Покажем, что Из теоремы 6 следует (7.5): x X : f ( x) f ( x) f ( x0 ) f ( x0 ) . T f ( x0 ) x x 0 . Т.к. f ( x 0 ) 0 , то f ( x) f ( x 0 ) 0 f ( x) f ( x 0 ) . ■ Теорема 9. Всякий локальный минимум выпуклой функции f (x ) , определѐнной на выпуклом множестве X , является еѐ глобальным минимумом на этом множестве. □ Пусть x * – точка локального минимума функции f (x ) , т.е. O (x* ) x O (x* ) : f (x* ) f ( x) . (7.8) Докажем, что (7.8) верно x X : f (x* ) f (x ) , x X . От противного. Пусть (7. 9) т.е. x * не является точкой глобального минимума. Рассмотрим отрезок, соединяющий точки x * и x , т.е. множество точек x : x x (1 ) x * , 0,1 . Т. к. X – выпуклое множество, то точки отрезка принадлежат множеству X , т.е. x X . Возьмѐм на нѐм точку x O ( x * ) , причѐм x x * . Поскольку точка x лежит на отрезке, то ей соответствует некоторое значение и выполняется неравенство x x (1 )x * , 0,1 . В силу выпуклости функции f (x ) для точек x и x * имеет место неравенство: 91 f (x ) или f (x ) f x )x * (1 f (x* ) f ( x ) (1 ) f (x* ) f (x ) f (x* ) (7.1 0) Но по предположению выполняется (7.9), т.е. f ( x ) f ( x * ) 0 . Тогда, усиливая неравенство (7.10), отбросим второе слагаемое: f ( x ) f ( x * ), f ( x * ) f ( x ) . А это противоречит (7.8), следовательно, x * - точка глобального минимума ( x принадлежит исходной окрестности O ( x * ) ). ■ Замечание. Пусть f (x ) – выпуклая дифференцируемая функция, заданная на выпуклом множестве X . Запишем условие (7.5) из теоремы 6: T x, x * X : f ( x ) f ( x * ) f (x* ) x x* . (7.1 1) Если x * – точка локального минимума функции f (x ) на множестf (x* ) ве X , то f ( x ) f (x* ) f ( x) Тогда x T из X. f (x* ) (x (7.11) x* ) , m следует, f ( x) что f (x* ) 0 . T f (x* ) (x x* ) m 0. (7.1 2) Это условие в силу теоремы 6 является необходимым и достаточным условием локального минимума функции f (x ) . В силу теоремы 9 для выпуклой функции это условие является и условием глобального минимума. Геометрически условие (7.12) означает, что, градиент f ( x * ) , если он отличен от нуля, составляет не тупой угол 1 с вектором направленным из точки x * в+ любую точку x X (рис. 7.7). Точка ~ x не удовле- ряет условию (7.12): из Xнеѐ можно сдвинуться, уменьшая ваясь в X . Т.е. для точки ~ x существует точка x , где f (x) , но . Если 2 x * X – внутренняя точка, тоf x f ( x * ) 0 . (Знаками «+» и «-» указано направление возрастания и убывания значений линий уровня) 2 + - x f (x ) 1 x x Рис. 7.7 x 2 92 Итак, если f (x ) - дифференцируемая и выпуклая функция, то для нахождения глобального оптимума (минимума f (x ) ), т.е. точки минимума, надо найти f ( x) f x1 , f x2 ,..., f xn , приравнять его к нулю и получить систему n нелинейных уравнений, решение которой по теореме 9 является точкой глобального минимума. Для того, чтобы определить, является ли дважды непрерывно дифференцируемая функция f (x ) выпуклой, надо согласно теореме 7 проверить, является ли еѐ матрица Гессе ределѐнной матрицей. 2 f ( x ) положительно полуоп- 7.3. Условная оптимизация при ограничениях-равенствах Задача ставится следующим образом: найти точку (точки) x Rn , в которой функция f ( x ) f ( x1 , x 2 ,...,x n ) принимает минимальное (максимальное) значение и при этом выполняются следующие неравенства (равенства): g k ( x) g k ( x1 , x 2 ,...,x n ) 0, k 1, m , которые называются ограничениями. f ( x) f ( x1, x2 ,..., xn ) min, g k ( x) g k ( x1, x2 ,..., xn ) 0, k 1, m . (7.1 3) Множество точек x , удовлетворяющих системе неравенств g i (x ) 0 , называется областью допустимых решений (ОДР), а точка x принадлежащая ОДР и доставляющая min (max) функции f (x ) , называется точкой условного экстремума. Пример 1. Минимизировать функцию f ( x ) x12 x 22 2x 2 при ограничениях g 1 ( x ) x12 x 22 1 0 g 2 ( x ) 2 x1 4 x 2 g 3 ( x ) x1 0 g 4 ( x) x 2 0 1 0 ОДР f ( x) x12 (x2 1) 2 1 Линии уровня x12 ( x 2 1) 2 1 c const (рис.7.8). Оптимальное решение находится в углах многоугольника, ограничивающего ОДР (на рис. 7.8. заштрихована). Оптимальное решение: x1 1, x 2 0 : x* (1/ 2, 0) , f ( x* ) 1/ 4. 93 x2 g1 ( x) 0 ОДР 1 x1 x -1 g 2 ( x) 0 Рис. 7.8 Будем рассматривать сначала задачи нелинейного программирования, когда все ограничения являются равенствами f ( x) f ( x1, x2 ,..., xn ) min, g k ( x) 0, k 1, m (m n). (7.14) Очевидный путь решения такой задачи следующий: из m ограничений-равенств выразить через ( n m) переменных остальные m и подставить их в f (x ) . Получим задачу безусловной оптимизации в пространстве меньшего измерения Rn m . И еѐ надо решать. Практически эта задача никогда неразрешима. Оригинальный метод решения задачи (7.3) предложил Лагранж. Сначала рассмотрим геометрическую интерпретацию. Рассмотрим когда n 2, m 1. z f ( x1, x2 ) min, g ( x1, x2 ) 0, z f ( x1 , x2 ) – некоторая поверхность. Уравнение g ( x1 , x 2 ) 0 определяет кривую на плоскости ( x1 , x 2 ) . Решением задачи минимизации f (x ) при g ( x1 , x 2 ) 0 будет точка x * ( x1* , x 2* ) , лежащая на кривой L : g ( x1 , x 2 ) 0 , которой соответствует точка B , лежащая на цилиндрической поверхности, у которой образующей является кривая g ( x1 , x 2 ) 0 , а направляющая параллельна оси Oz . Отметим, что условный минимум в точке B не совпадает с глобальным минимумом в точке A (рис.7.9). 94 z z f ( x1, x2 ) C B A x2 x g ( x1, x2 ) 0 L x1 Рис. 7.9 В математическом анализе доказывается теорема: для того чтобы точка x * была точкой условного экстремума функции f (x) , т.е. решением задачи (7.12), необходимо, чтобы выполнялись условия: f ( x* ) m k g k ( x* ) 0, i 1, n, k 1 g k ( x* ) 0, k 1, m, где k – некоторые постоянные, одновременно не равные нулю. Эти условия удобнее записать через функцию Лагранжа: m F ( x, ) f ( x) k gk ( x) , k 1 ( 1, 2 ,..., m ) ( 1, 2 ,..., m ) . Метод Лагранжа состоит в следующем: 1. Составляют функцию Лагранжа F (x, ) ; F F , ( j 1, n, k 1, m) ; 2. Находят частные производные xj k 95 3. Решают F xj систему F 0, и отыскивают точку g k ( x), k ( x10 , x20 ,..., xn0 ) , удовлетворяющую этой системе. Найденные точки исследуют далее на минимум (или на максимум). Метод множителей Лагранжа является классическим методом решения задач нелинейной оптимизации. К сожалению, при практическом применении метода могут быть большие трудности, связанные с решением системы уравнений, которые сужают область его применения. Однако он является аппаратом, который используется для обоснования различных современных численных методов решения практических задач оптимизации. Пример 2. Найти оптимальное распределение ограниченного ресурса в a единиц между n потребителями, если прибыль, получаемая при выделении j - му потребителю x j единиц ресурса, вычисляется по форx муле c j x j . Математическая модель задачи: n max z n cj xj , j 1 xj j 1 n Обозначим g ( x ) a. x j – ограничение. a j 1 Составим функцию Лагранжа n F ( x1 , x 2 ,..., x n , ) n cj xj j 1 a xj . j 1 Находим частные производные и приравниваем к нулю. cj F 0, j 1, n, xj 2 xj n F a xj 0. j 1 Решаем систему: cj a j c 2j 14 2 1 4 n c 2j 2 j 1 * c 2j xj 4x j 2 xj n c 2j 2 1 n 2 * cj ; xj 4a j 1 2 4 ac 2j n c 2j j 1 . , 96 Т.о., если j - му потребителю будет выделено ac 2j n единиц ресурса, c 2j j 1 то суммарная прибыль достигает максимальной величины и будет равна z * n cj j 1 ac 2j n c 2j j 1 n c 2j a n . j 1 c2 j j 1 В заключение дадим множителям Лагранжа экономическую интерпретацию. Рассмотрим простую задачу оптимизации z f ( x1 , x 2 ) min (max) . (7.15) g ( x1 , x 2 ) b . (7.1 6) Предположим, оптимум достигается в точке x * ( x1* , x 2* ) и оптимальное значение z есть z * f ( x1* , x 2* ) . Предположим, что правая часть в ограничении (7.16) меняется. Тогда координаты точки оптимума x * будут функциями b , т.е. x1* x1* (b), x 2* x 2* (b), z * f x1* (b), x 2* (b) . Найдѐм производную по b от z * f dx1* f dx 2* dz * . db x1 db x 2 db (7.17) С другой стороны, g x1* (b), x 2* (b) b . После дифференцирования этого тождества по b получим: g dx1* g dx 2* 1. x1 db x 2 db (7.18) Кроме того, в точке x * выполняется необходимое условие экстремума функции Лагранжа: f ( x1 , x 2 ) венств при n 2 и m 1 получаем: g ( x1 , x 2 ), f x1 g ; x1 xj f x2 0 . Из этих раg . x2 97 Подставив последние равенства в (7.17) и учитывая (7.18), получим: g dx1* x1 db dz * db dz * db g dx 2* x 2 db g dx1* x1 db g dx 2* x 2 db , . dz * i , i 1, m . dbi Если z интерпретировать как доход или стоимость, а bi - как объѐм некоторых ресурсов, то множители Лагранжа i показывают, как изменится максимальный доход (или минимальная стоимость), если количество ресурса i -го вида увеличится на единицу (если перейти к конечным b, z ). Аналогично для общей задачи получили бы 7.4. Условная оптимизация при ограничениях-неравенствах f ( x) Рассмотрим задачу min, gi ( x) bi , i 1, m. Ограничения в виде неравенств могут быть преобразованы в ограничения-равенства добавлением к каждому неравенству неотрицательной переменной, которая называется ослабленной (слабой). gi ( x) bi gi ( x) yi2 bi или (7.19) g i ( x) y i2 bi 0 . Следовательно, задача сводится к минимизации функции f (x ) при наличии ограничений (8.19) в виде равенств. По правилу, которое сформулировано в 7.3, составим функцию Лагранжа m F ( x, , y ) f ( x) i gi ( x ) yi2 bi . i 1 (7.20) Необходимые условия должны выполняться в стационарной точке, где частные производные по всем переменным равны нулю. m gi F f (7.21) 0; j 1, n , i xj xj i 1 xj F i gi ( x ) yi2 bi 0; i 1, m , (7.22) 98 F yi 2 i yi Умножив (7.23) на bi 0; i 1, m . yi , получим (7.23) 2 i yi 0 . Тогда из (7.22) y i2 следует: gi n i bi gi ( x ) 0, i 1, m . (7.24) i 1 Уравнения (7.21), (7.22) и (7.24) являются необходимым условием минимума в точке x * при наличии ограничений. Условия (7.24) называются условиями дополняющей ………… и означают, что либо i 0 , либо bi gi ( x* ) 0 . Если i 0 , то g i ( x * ) bi , и ограничения являются активными типа равенств. С другой стороны, если ограничение является ограничением в виде строгого неравенства g i ( x * ) bi , то соответствующий множитель Лагранжа i 0 . В самом деле, если g i ( x * ) bi , то рассматривается минимум, удовлетворяющий ограничению, которое является неактивным (пассивным) и которым можно пренебречь, а соответствующий множитель Лагранжа i 0 . Но предварительно не известно, какими ограничениями можно пренебречь. Есть также дополнительные условия, которые должны быть выполнены в точке минимума при наличии ограничений, а именно i 0 . Покажем это. Предположим, что уравнения (7.21), (7.22), (7.23) справедливы в точке ( x* , * , y* ) . Фактический минимум функции при наличии ограничений есть z f ( x * ) и z можно рассматривать как функцию от bi , т.к. изменения bi будут изменять ограничения, и, таким образом, изменять z * саму функцию z . Покажем, что i . По формуле дифференцироваbi ния сложной функции имеем: n z f xj , i 1, m , bi j 1 x j bi (7.25) где частные производные вычисляются в точке x * . Поскольку gk ( x) bk yk2 , то g k n g k x j 0, если i k . bi j 1 x j bi 1, если i k (7.26) Составим линейную комбинацию: 99 m z z * gk * k i. bi k 1 bi bi (7.27) С другой стороны из (7.25) и (7.26): m n m xj z z * gk * k bi bi k 1 n n z xj j 1 j 1 j 1 * k gk xj xj bi xj bi n k k 1 j 1 gk xj xj bi 0. (7.28) z z * * 0 или i i . bi bi Из неравенства g ( x) bi следует, что с возрастанием bi область ограничений расширяется, а это не может привести к увеличению значения zmin функции f (x ) , находящегося внутри области ограничений, а z может лишь уменьшить его. Поэтому, 0 , а следовательно, i 0 . b Условия оптимальности (7.21), (7.23), (7.24) в стационарной точке x можно записать в векторном виде: Из (7.27) и (7.28) следует, что m f (x ) xi gi ( x ) 0 , (7.29) i 1 m i (bi gi ( x )) 0 , (7.30) i 1, m . (7.31) i 1 i При рассмотрении задачи минимизации f ( x) при условиях gi ( x) bi может случиться так, что не существует таких i , i 1, m , при которых без дополнительных предположениях о природе функций gi ( x) были бы справедливы условия (7.23) – (7.31), где x – оптимальное решение. Эти дополнительные предположения называют условиями регулярности ограничений. В частности при ограничениях равенства такими условиями является линейная независимость векторов–градиентов ограничений. В теории нелинейной оптимизации центральное место занимает теорема Куна–Таккера, обобщающая классический метод множителей Лагранжа на случай, когда в нелинейной задаче помимо ограничений – равенств имеются также неравенства. 100 Теорема 10. Пусть функция f ( x) gi ( x) i 1, m имеют непрерывные частные производные на некотором открытом множестве Rn , содержащим x . Если x является точкой минимума функции f ( x) при ограничениях gi ( x) bi , удовлетворяющих условию регулярности в виде линейной независимости градиентов, то существуют такие 1, 2 ,..., m , что справедливы условия (7.29)–(7.31). Заметим, что множители Лагранжа в задаче нелинейной оптимизации с ограничениями равенствами являются знаконеопределенными, тогда как в теореме Куна–Таккера они должны быть неположительными. Условие регулярности было впервые предложено Куном и Таккером в 1951 году и имеет различные формы. В частном случае, когда f ( x) и все gi ( x) являются выпуклыми функциями, то условие регулярности состоит в существовании такой точки x , где gi ( x) 0 для всех i 1, m . Это условие называется условием регулярности Слейтера. Рассмотрим функцию F ( x, ) от векторов x и . Определение 2. Пара называется седловой точкой функции x , F ( x, ) , если при всех 0 и x Rn выполняется условие F ( x , ) F ( x , ) F ( x, ) . (7.32) Это неравенство называется неравенством для седловой точки и означает, что n -мерная точка x является точкой минимума для функции F ( x, ) , а m –мерная точка – точкой максимума для функции F ( x , ) , т.е. F (x , ) max min F ( x, ) Между x Rn понятием min max F ( x, ) . x Rn седловой точки функции Лагранжа m F ( x, ) f ( x) i gi ( x ) и решением задачи нелинейной оптимизации i 1 существует связь, которую устанавливает следующая теорема. Теорема 11. Пусть f ( x) и все функции gi ( x) выпуклые, а также gi ( x) удовлетворяют условию регулярности Слейтера. Тогда точка x является решением задачи нелинейной оптимизации 101 f ( x) min, gi ( x) bi , (i 1, m), (7.33) тогда и только тогда, когда существует такой вектор 0 , что выполняется неравенство для седловой точки (7.32) и равенство (7.30). (Без доказательства). Таким образом, при выполнении условий теоремы II задача нелинейной оптимизации оказывается эквивалентной задаче об отыскании седловой точки функции Лагранжа. Пример 3. Написать условие Куна–Таккера для минимума функции f ( x) 3x 2 4 x1x2 5 x 2 при ограничениях x1 0, x2 0, x1 x 2 4 . Постановка задачи имеет вид: f ( x) 3x 2 4 x1x2 5 x2 min, x1 0, x2 0, x1 x2 4. Составим функцию Лагранжа F ( x, , y) 3x12 4 x1x2 5x22 1( y12 x1) 2 ( y2 2 x2 ) 2 3 ( y3 x1 x2 . Необходимыми условиями оптимальности будут: 6 x1 4 x2 1 3 0, 4 x1 10 x2 x1 0, x2 2 3 0, 0, x1 x2 0, 1x1 0 , 2 x2 0, 3 4 x1 x2 1, 2 , 3 x1 Эти 3, x2 1, 0, 0. 1 0, 2 условия 0, 3 22, fmin выполняются 44 . при 8. ЧИСЛЕННЫЕ МЕТОДЫ БЕЗУСЛОВНОЙ ОПТИМИЗАЦИИ 8.1. Общая характеристика численных методов оптимизации Решать задачу 4) 102 (8.1) f (x) min , x R n ( X ) , в практическом большинстве случаев приходится численно. При этом с помощью необходимых условий оптимальности можно свести эту задачу к другой и пытаться применять известные численные методы для еѐ решения. Так, например, если функция f (x ) дифференцируема, то воспользоваться каким-либо методом для решения системы уравнений f x 0 i 1, n . Однако наиболее эффективными являются методы, i разработанные специально для оптимизационных задач, т.к. они наиболее полно учитывают их специфику. Любой численный алгоритм (метод) решения задач оптимизации основан на точном или приближѐнном вычислении еѐ характеристик (значений f (x ) , f (x ) ,…, значений функций, задающих дополнительные области и т.д.). На основании этого строится приближение к решению задачи, т.е. искомой точки минимума x * , или, если такая точка не единственная, то к множеству точек минимума. Для каждой конкретной задачи вопрос о том, какие характеристики надо выбрать для вычисления, зависит от свойств функции и имеющихся возможностей по хранению информации. Так, если функция f (x ) недифференцируема, то нельзя пользоваться алгоритмами с производными f (x ) . Также, если доступен небольшой объѐм памяти ЭВМ, то для решения задач высокой размерности нельзя пользоваться алгоритмами, которые требуют на каждом шаге хранения всех матриц, используемых в методе. Алгоритмы, использующие лишь информацию о значениях минимизируемой функции, называются алгоритмами нулевого порядка; алгоритмы, использующие информацию о первых производных, называются алгоритмами первого порядка; алгоритмы, использующие информацию о вторых производных, называются алгоритмами второго порядка. Основное отличие этих алгоритмов друг от друга состоит в следующем: методы нулевого порядка не требуют знания производных исследуемой функции, следовательно их можно применять для наиболее широкого класса функций, в том числе и для тех, у которых невозможно вычислить производные в конкретных точках; методы первого порядка требуют знания первых производных, гарантированно сходятся при исследовании выпуклых функций и обеспечивают сходимость со скоростью геометрической прогрессии; 103 методы второго порядка требуют знание вторых производных, гарантированно сходятся при исследовании выпуклых функций и обеспечивают сходимость с квадратичной скоростью. На первый взгляд может показаться, что наиболее предпочтительными оказываются методы нулевого порядка как самые универсальные, однако это не так. За их универсальность приходится расплачиваться уменьшение точности расчета по сравнению с методами первого, а, тем более, второго порядка. Если же исследовать в программе ЭВМ задает высокую точность расчета, то он проигрывает в скорости вычисления по сравнению с методами первого и второго порядка, использующими ту же самую точность при намного меньшем числе итераций. С другой стороны, возможно, при решении конкретной задачи исследователю не требуется высокая точность и увеличение времени расчета на несколько секунд малозначимо, а процесс нахождения производных очень затруднителен, то имеет смысл обратиться к более универсальным методам. Иначе говоря, нельзя рекомендовать, не только тот или иной численный алгоритм, но даже целый класс методов. Выбирать тот или иной алгоритм для поиска безусловного экстремума следует пользователю, которых владеет информацией об исследуемой функции, ее производных, требуемой точности расчета и так далее. Работа любого алгоритма состоит, как правило, из двух этапов. На первом этапе вычисляются все необходимые характеристики задачи, а на втором – строится приближение к решению. Если все точки, в которых вычисляется функция f (x ) на втором этапе, выбираются одновременно до начала вычислений, то алгоритм минимизации называется пассивным. Такие алгоритмы используются в связи с применением многопроцессорных ЭВМ, в связи с условиями постановки и использования физических экспериментов, результатом которых является значение минимума функции и т.д. Однако в большинстве задач точки вычисления выбираются поочерѐдно, т.е. точка x i 1 выбирается, когда уже выбраны предыдущие точки x 1 , x 2 ,…, x i и в каждой произведены необходимые расчѐты. Эти 104 алгоритмы называются последовательными и при записи таких методов оптимизации обычно пользуются соотношениями вида xk 1 xk k 0,1,2,... . k h, k R, (8. 2) При этом конкретный алгоритм определяется заданием точки x0 и правилами выбора векторов h k и чисел k . Вектор h k определяет направление (k 1) шага, а коэффициент k – длину этого шага. Ясно, что при h k 1 длина отрезка, соединяющего точки x k и x k 1 не равна k . Обычно название метода оптимизации определяется выбором h k , а его различные варианты связываются с выбором k . Наряду с термином «шаг метода» используется термин «итерация метода». Важной характеристикой итерационных методов является сходимость. Определение 1. Метод (8.2) называется сходящимся, если x k x * при , где x * – решение задачи (8.1). k Если f ( x k ) f ( x * ) , то иноf ( x) гда также говорят, что метод (8.2) сходится (по функции), а последовательность x k называется минимизирующей. Минимизирующая последовательность может и не сходиться к точке минимума. Так для функции на рис.8.1 минимизирующая последовательx2 2 ность x k k не сходится к точке x x 0 x1 1 x 4 4минимума x * 0 . Если x * не единственная, то Рис. 8.1 под сходимостью понимают сходимость последовательности x k к множеству X * точек минимума функции f (x ) Эффективность сходящегося метода характеризуется скоростью сходимости. Различают следующие скорости сходимости: 1) линейная сходимость: xk 1 x* q xk x* при k k0 , q const . 2) квадратичная сходимость: x k 1 x* q xk 2 x* , k k0 . (8.3) 105 3) сверхлинейная сходимость: xk 1 x* qk 1 xk x* , где 0 при k . Большинство теорем о сходимости методов оптимизации доказано в предположении выпуклости целевой функции. Для невыпуклых задач численные методы оптимизации обычно позволяют находить лишь локальные решения, т.е., вообще говоря, стационарные точки. Задача отыскания глобального оптимума чрезвычайно сложна даже для функций одной переменной. Получение достаточно точного решения многомерных задач глобальной оптимизации с помощью существующих в настоящее время методов оказывается вообще невозможным. При выборе подходящего метода решения реальных задач надо руководствоваться здравым смыслом, интуицией, опытом, а также результатами численных экспериментов. Для создания конкретного вычислительного алгоритма бесконечно шаговый метод надо дополнять условием остановки. Оно часто определяется имеющимися вычислительными ресурсами. Теоретически остановку можно производить, используя (8.3). На практике часто используются следующие условия остановки: qk 1 xk 1 xk 1, f ( x k 1) f ( xk ) f ( xk ) 3, 2, (8.4) причѐм, как правило, необходимо выполнение сразу двух или даже всех трѐх условий. Третье условие относится лишь к безусловной оптимизации дифференцируемых функций (равносильно условию стационарности). Вместо критериев (8.4), основанных на понятии абсолютной погрешности, можно использовать аналогичные критерии, основанные на понятии относительной погрешности. xk f (x k 1 ) 1 xk 1 f (x k ) f (x k 1 ) 2 3 xk 1 1 1 1 1 , f (x k 1 ) , f (x k 1 ) . 106 Отметим, что выполнение указанных критериев и других подобных им эвристических условий остановки в общем случае не гарантирует доставление необходимой точности решения задачи. Многие методы оптимизации относятся к числу методов спуска, в которых движение к минимуму на каждом шаге выбирается из числа направлений убывания функции f (x ) . Говорят, что вектор h задаѐт направление убывания функции f (x ) в точке x , если f ( x h) f ( x) для любых достаточно малых 0 . Сам вектор h называется направлением убывания. Векторов h множество в каждой точке. Обозначим это множество U ( x, f ) . Т.о., если любой достаточно мелкий сдвиг из точки x в направлении h приводит к уменьшению f , то h U ( x, f ) . Очевидно, что из формулы Тейлора следует необходимое и доста1 T 2 T f h h f h o( 2 ) . точное условие убывания: f 2 Если функция f (x ) дифференцируема в точке x R n , то для убывания f (x ) в направлении вектора h необходимо, чтобы h U ( x, f ) . T f (x ) 0 Если h U ( x, h) , то T f ( x )h 0 . Проверить самостоятельно. Геометрически это означает, что вектор h составляет тупой угол с градиентом f (x ) . Метод (9.2), т.е. x k 1 xk kh k называется методом спуска, если вектор h k задаѐт направление убывания функции f (x ) в точке x k , т.е. hk U ( x k , f ), k 0,1,2,..., а число k 0 и f ( x k 1 ) f ( x k ) . Примерами метода спуска являются метод покоординатного спусf (x k ) и hk ка и градиентный метод, в которых h k f ( x k ) соотxi ветственно. Важным в методах спуска является выбор параметра спуска k . Имеются следующие способы выбора k . 1. Минимизация функции вдоль заданного направления. Коэффициk енты k можно определить из условия f ( x k min f ( x k hk ) , kh ) (8.5) где минимум берѐтся по 0 . Этот метод является в некотором смысле наилучшим, т.к. он обеспечивает достижение наименьшего значения f вдоль заданного направления. Однако он требует решения на каждом шаге одномерной задачи минимизации, которую часто прихо- 107 дится решать численно, а, следовательно, приводит к увеличению объема вычислений. В простейших случаях можно получить явное решение для , например, для квадратичной функции 1 T f ( x) x Ax b x c . 2 2. Априорный выбор k , определяющих длину шага. Иногда величины k выбираются до начала вычислений. Например, требуют выполнения условий 0 ( k 0,1,2,...), k k 0 k , k 0 Можно взять 2 k . c . k 1 Интуитивно это можно объяснить следующим образом. Условие k сходимости ряда k 0 2 k накладывается, чтобы добиться достаточно бы- строй сходимости последовательности k к нулю с целью обеспечения сходимости метода в окрестности точки экстремума x * . Условие расходимости ряда k 0 k призвано обеспечить достижение точки x * даже при неудачном начальном приближении x0 , т.е. при большом расстоянии от x0 до x * . 3. Дробление шага. Если h k – направление убывания, то дробление 0 и шага состоит в следующем. Выбираются некоторые константы 1 ). Для коэффициента проверяется условие 1 (часто 2 f (xk hk ) f (xk ) . (8.6) Если оно выполняется, то полагают k , если нет, то производится дробление шага, т.е. полагают и вновь проверяется выполнение условия (8.6). Процесс дробления, т.е. умножения текущего значения на продолжается до тех пор, пока условие (8.6) не окажется выполненным. Этот процесс не может быть бесконечным, потому что h k – направление убывания. Первое , при котором условие выполнено, принимается за k . Эти три способа выбора длины шага в методах спуска не являются единственными, на практике есть и другие. 108 8.2. Методы нулевого порядка одномерной оптимизации Сначала рассмотрим задачу одномерной оптимизации, которая состоит в нахождении безусловного минимума функции одной переменной f ( x) , т.е. точка x R , такой, что f ( x ) min f ( x) . R Для методов одномерной оптимизации обычно необходимо задавать начальную (априорную) информацию о положении точки минимума с помощью начального интервала неопределенности L0 a0 , b0 , которому принадлежит искомая точка x , но ее точное значение неизвестно. Большинство известных методов одномерной оптимизации применяется для унимодальных на интервале функций, т.е. функций, которые достигают глобального минимума на L a, b в единственной точке x , причем слева от x f ( x) строго убывает, а справа от x f ( x) строго возрастает. Ясно, что непрерывная строго выпуклая функция на a, b является унимодальной на этом отрезке. Однако, конечно унимодальными функциями могут быть и разрывные, и невыпуклые на a, b функции. Стратегия поиска точки минимума включает в себя три этапа: 1. Выбор начального интервала L0 a0 , b0 , границы которого должны быть такими, чтобы функция f ( x) была на нем унимодальной. 2. Уменьшение интервала L0 по некоторому правилу. 3. Проверка условия окончания. Поиск заканчивается, если длина текущего интервала будет меньше заданной величины. Ответом является множество точек последнего интервала неопределенности, среди которых выбирается решение задачи x . В некоторых методах заранее задается и определяется количество N вычислений функции и, тогда продолжительность поиска ограничена. Для выбора начального интервала неопределенности существуют различные алгоритмы, например, алгоритм Сенна. 1. Произвольно из каких либо соображений задаются следующие параметры: x0 – некоторая точка, t 0 – величина шага. Положим k 0 . 2. Вычисляются значения f ( x) в трех точках x0 t , x0 , x0 t . 3. Проверка условия окончания. 3а. Если f ( x0 t ) f ( x0 ) f ( x0 t ) , то начальный интервал неопределенности найден a0 , b0 x0 t ; x0 t . 109 3б. Если f ( x0 t ) f ( x0 ) f ( x0 t ) , то функция не является унимодальной и интервал неопределенности не может быть найден. Необходимо задать другую начальную точку x0 . 4. Определяется величина . 4а. Если f ( x0 t ) f ( x0 ) f ( x0 t ) , то t , a0 x0 , x1 x0 t , k 1 . 4б. x1 Если f ( x0 t ) f ( x0 ) f ( x0 t ) , то t , a0 x0 t , k 1 . 5. Находится следующая точка x k 1 x k 2k . 6. Проверка условия убывания функции. 6а. Если f ( x k 1) f ( x k ) и t , то a0 x k , если f ( x k 1) t , то b0 x k , в обоих случаях положить k 6б. Если f ( x k 1) b0 x0 , k 1 и перейти к п.5. f ( x k ) , то процесс завершается. При x k 1 , а при t положить a0 x k 1 . В результате имеем искомый интервал L0 8.2.1 Метод деления интервала пополам f ( xk ) и t положить a0 , b0 . Метод относится к последовательным алгоритмам и позволяет исключить из дальнейшего рассмотрения на каждой итерации в точности половину текущего интервала неопределенности. Задается начальный интервал неопределенности, а алгоритм уменьшения интервала, основан на анализе величин функции в трех точках, равномерно распределенных на текущем интервале (делящих его на четыре равные части). Условия окончания процесса поиска стандартные: поиск заканчивается, когда длина текущего интервала неопределенности оказывается меньше установленной величины. Алгоритм 1. Задать начальный интервал неопределенности L0 требуемую мощность. 2. Положить k 0 . ak bk , L2k bk 3. Вычислить среднюю точку: xk c 2 a0 , b0 и l ak , f x c k . 0 – 110 4. Вычислить точки: yk L2k L2k и f ( yk ) , f ( zk ) . Заме4 4 тим, что точки yk , xkc , zk делят интервал L0 [ak , bk ] на четыре равные части. 5. Сравнение значений f ( yk ) и f ( xkc ) . 1 , zk bk f ( xkc ) , исключить интервал xkc , bk , положив bk 1 xkc , a k . Средней точкой нового интервала становится точка yk : xkc 1 yk 5а. Если f ( yk ) ak ak (рис. 8.2, а). Перейти к п. 7. 5б. Если f ( yk ) f ( xkc ) , перейти к п. 6. 6. Сравните f ( zk ) с f ( xkc ) . 6а. ak xkc , bk 1 f ( zkc ) , Если f ( zk ) 1 исключить интервал a , xkc , k положив bk . Средней точкой нового интервала становится точка zk : xck 1 zk (рис. 8.2, б). Перейти к п. 7. 6б. Если f ( zk ) f ( xkc ) , исключить интервалы ( a k , y k ] , ( z k , bk ] , положив c zk . Средней точкой нового интервала остается xk : ak 1 yk , bk xkc 1 xkc (рис. 8.2, в). 1 7. Вычислить L2( k 1) 7а. k 1 x* L2( k Если 1) L2 bk ak 1 l, 1 и проверить условия окончания. процесс поиска завершается и [ak 1, bk 1] . В качестве приближенного решения можно взять * середину последнего интервала: x 7б. Если L2( k 1) xkc l , то положить k 1. k 1 и перейти к п. 111 f f f ( zk ) f ( yk ) f f ( yk ) ( xkc ) f ( xkc ) f ( zk ) ak yk xkc zk x y k xkc ak bk Новый интервал Текущий интервал а zk x bk Новый интервал Текущий интервал б f f ( yk ) f ( xkc ) f ( zk ) ak yk xkc zk bk x Новый интервал Текущий интервал в Рис. 8.2 – Геометрическая интерпретация метода деления интервала пополам Для метода деления интервала пополам характеристика относительного уменьшения начального интервала неопределенности находится по формуле LN 1 , где N – количество вычислений функции. R( N ) N L0 22 Средняя точка последовательно получаемых интервалов всегда совпадает с одной из трех пробных точек, найденных на предыдущей итерации. Следовательно, на каждой итерации требуется два новых вычисления функции. 112 Если задана величина R (N ) , то требуемое для достижения желаемой точности количество вычислений функции находится как наименьшее целое, 2 ln R( N ) удовлетворяющее условию N . ln 0.5 Текущие интервалы имеют четные номера L0 , L2 , L4 , …, где индекс указывает на сделанное количество вычислений функций. 8.2.2. Метод дихотомии Метод относится к последовательным стратегия. Задается начальный интервал неопределенности и требуемая точность. Алгоритм опирается на анализ значений функций в двух точках. Для их нахождения текущий интервал неопределенности делится пополам и в обе стороны от середины откладывается по , где – малое положительное число. Условия окончания 2 процесса поиска стандартные: поиск заканчивается, когда длина текущего интервала неопределенности ок5азвывается меньше установленной величины. Алгоритм 1. Задать начальный интервал неопределенности L0 a0 , b0 и 0 – малое число, l 0 – точность. 2. Положить k 0 . ak bk ak bk , f ( yk ), zk , f ( zk ) . 3. Вычислить yk 2 2 4. Сравнить f ( yk ) с f ( zk ) : 4а. Если f ( yk ) f ( zk ) , положить ak 1 yk , bk 1 bk (рис. 8.3, а) и перейти к п. 5. 4б. Если f ( yk ) f ( zk ) , положить ak 1 ak , bk 1 zk (рис. 8.3, б). 5. Вычислить L2(k 1) 5а. Если L2( k 1) bk 1, 1 ak 1 и проверить условие окончания: процесс поиска завершается и . В качестве приближенного решения можно взять ak 1 bk 1 середину последнего интервала x . 2 5б. Если L2(k 1) 1, положить k k 1 и перейти к п. 3. x L2(k 1) ak 1, b k 1 113 f f f ( zk ) f ( yk ) f ( zk ) f ( yk ) 2 2 yk ak zk bk x yk ak 2 Новый интервал 2 zk bk x Новый интервал Текущий интервал а Текущий интервал б Ри с. 8.3. Для метода дихотомии характеристика относительного уменьшения начального интервала неопределенности находится по формуле R( N ) 1 , N 22 где N – количество вычислений функции. Текущие интервалы неопределенности L0 , L2 , L4 ,... имеют четные номера, указывающие на количество сделанных вычислений функции, как и в методе деления интервала пополам. Эффективность методов дихотомии и деления интервала пополам при малых можно считать одинаковой. 8.2.3. Метод золотого сечения Для построения конкретного метода одномерной минимизации, работающего по принципу последовательного сокращения интервала неопределенности, следует задать правило выбора на каждом шаге двух внутренних точек. Конечно, желательно, чтобы одна из них всегда использовалась в качестве внутренней и для следующего интервала. Тогда число вычислений функции сократится вдвое и одна итерация потребует расчета только одного нового значения функции. В методе золотого сечения в качестве двух внутренних точек выбираются точки золотого сечения. Говорят, что точка производит «золотое сечение» отрезка, если отношение длины всего отрезка к большей части равно отношению большей части к меньшей. На отрезке a0 , b0 имеются две симметричные относительно его концов точки y0 и z0 : b0 b0 a0 y0 b0 y0 y0 a0 b0 a0 z0 a0 z0 a0 b0 z0 1 5 2 1,618 114 Кроме того, точка y0 производит золотое сечение отрезка a0 , z0 , а точка z0 – отрезка y0 , b0 . Метод относится к последовательным стратегиям. Задается начальный интервал неопределенности и требуемая точность. Алгоритм уменьшения интервала опирается на анализ значений функции в двух точках. В качестве точек вычисления функции выбираются точки золотого сечения. Тогда с учетом свойств золотого сечения на каждой итерации, кроме первой, требуется только одно новое вычисление функции. Условия окончания процесса поиска стандартные: поиск заканчивается, когда длина текущего интервала неопределенности оказывается меньше установленной величины. Алгоритм 1. Задать начальный интервал неопределенности L0 a0 , b0 , точность l 0 yk yk 1 1 2. Положить к=0. 3. Вычислить 3 5 3 5 y0 a0 (b0 a0 ); z0 a0 b0 y0 , 0,38196 . 2 2 4. Вычислить f ( yk ), f ( zk ) . 5. Сравнить f ( yk ), f ( zk ) : f ( yk ) f ( zk ) , то положить ak 1 ak , bk 5а. Если ak 1 bk 1 yk , zk 1 yk (рис. 8.4, а). Перейти к п. 6. 5б. Если положить ak 1 yk , bk f ( yk ) f ( zk ) , zk , zk 1 ak 1 bk 1 zk (рис. 8.4, б). f f 1 zk и 1 bk и f ( zk ) f ( yk ) f ( yk ) f ( zk ) yk ak yk 1 zk zk 1 zk bk Новый интервал Текущий интервал а 6. Вычислить x ak yk zk yk 1 1 x bk Новый интервал Текущий интервал б Ри с. 8.4. ak 1 bk 1 и проверить условие окончания: 115 ak 1, bk 1 . В каче6а. Если l , процесс поиска завершается и x стве приближенного решения можно взять середину последнего интервала: ak 1 bk 1 x . 2 6б. Если l , положить к=к+1 и перейти к п. 4. Для метода золотого сечения характеристика относительно уменьшения начального интервала неопределенности находится по формуле R( N ) (0,618) N 1 , где N – количество вычислений функции. Текущие интервалы неопределенности имеют следующий вид: L0 , L2 , L3 , L4 ,... Они отражают тот факт, что на первой итерации производится два вычисления функции, на последующих – по одному. Сокращение длины L L3 L2 интервала неопределенности постоянно: 0 ... 1,618 . L2 L3 L4 Если задана величина R( N ) , то требуемое для достижения желаемой точности количество вычислений функции находится как наименьшее целое ln R( N ) число, удовлетворяющее условию N 1 . ln 0,618 8.3. Методы нулевого порядка многомерной оптимизации 8.3.1. Метод конфигураций Задача многомерной оптимизации состоит в следующем: требуется найти безусловный минимум функции f ( x) f x1, x2 ,..., xn многих переменных, то есть найти такую точку x* R n , что f ( x*) min f ( x) . x Rn Метод конфигураций (метод Хука-Дживса) представляет собой комбинацию исследующего поиска с циклическим изменением переменных и ускоряющего поиска по образцу. Исследующий поиск ориентирован на выявление локального поведения целевой функции и определение направления ее убывания вдоль «оврагов». Полученная информация используется при поиске по образцу при движении вдоль «оврагов» . Исследующий поиск начинается в некоторой начальной точке x 0 , называемой старым базисом. В качестве множества направлений поиска выбирается множество координатных направлений. Задается величина шага, которая может быть различной для разных координатных направлений и переменной в процессе поиска. Фиксируется первое координатное направление и делается шаг в сторону увеличения соответствующей переменной. Если значение функции в пробной точке меньше значения функции в исходной точке, шаг считается удачным. В противном случае необходимо вернутся в исход- 116 ную точку и сделать шаг в противоположном направлении с последующей проверкой поведения функции. После перебора всех координат исследующий поиск завершается. Полученная точка называется новым базисом (на рис. 8.2 в точке x 0 произведен исследующий способ и получена точка x 1 – новый базис). Если исследующий поиск с данной величиной шага неудачен, то она уменьшается и процедура продолжается. Поиск заканчивается, когда текущая величина шага станет меньше некоторой величины. x2 3 f ( x) ( x1 1) 2 x22 4 3 2 8 1 x x 14 6 2 5 x 1 2 7 10 x1 -1 15 x1 9 f ( x) 1 17 1 2 4 1 x3 2 -1 13 11 12 16 -2 Рис. 8.5 – Геометрическая интерпретация метода конфигураций Поиск по образцу заключается в движении по направлению от старого базиса к новому (от точки x 0 через точку x 1 , из точки x 1 через точку x 2 , из x 2 через x3 на рисунке 8.5). Величина ускоряющего шага задается ускоряющим множителем . Успех поиска по образцу определяется с помощью исследующего поиска из полученной точки (например, из точек 6, 11, 15 на рисунке 8.5). Если при этом значение в наилучших размерах меньше, чем в точке предыдущего базиса, то поиск по образцу удачен (точки 6, 11 – результат удачного поиска по образцу, а точка 15 – неудачного). Если поиск по образцу неудачен, происходит возврат в новый базис, где продолжается исследующий поиск с уменьшенным шагом. На рис. 8.5 удачный поиск отображается сплошными линиями, а неудачный – пунктирными, числа соответствуют порождаемым алгоритмом точки. Обозначим через d 1 ,...,d n – координатные направления. 1 d1 , d2 ... 1 ,..., d n ... . ... 1 117 При поиске по направлению di меняется только переменная x i , а остальные переменные остаются зафиксированными. Алгоритм 1. Задать начальную точку x 0 , число 0 для остановки алгоритма, начальные величины шагов по координатным направлениям 1 ,..., n , ускоряющим множитель 0 , коэффициент уменьшения шага 1 . Положить y1 x0 , i 1, k 0 . 2. Осуществить исследующий поиск по выбранному координатному направлению i i 2а. Если f ( yi f ( y1i ,..., yii i di ) f ( y ) , то есть i ,..., yn ) , f ( y1i ,..., y ii ,..., yni ) yi шаг считается удачным. В этом случае положить 1 yi i di и перейти к п. 3. 2б. Если в пункте 2а шаг неудачен, то делается шаг в противоположном i f ( y i ) , то есть f ( y1i ,..., yii направлении. Если f ( y i idi ) i ,..., yn ) , f ( y1i ,..., yii ,..., yni ) шаг считается удачным. В этом случае следует положить yi 1 yi i di и перейти к п. 3. 2с. Если в пунктах 2а и 2б шаги неудачны, положить y i 1 y i . 3. Проверить условия: 3а. Если i n , то положить i i 1 и перейти к п. 2 (продолжить исследующий поиск по оставшимся направлениям). 3б. Если i n , проверить успешность исследующего поиска: если f ( y n 1) f ( x k ) , перейти в п. 4; если f ( y n 1) 4. y1 Провести f ( x k ) , перейти в п. 5. поиск по образцу. Положить xk 1 yn 1 , i . xk 1 ( xk 1 xk ), i 1, k k 1 и перейти к п. 2. 5. Проверить условие окончания. 5а. Если все i , то поиск закончить: x* x k . 5б. Для тех i , для которых i , уменьшить величину шага: i Положить y1 x k , x k 1 x k , k k 1, i 1 и перейти к п. 2. В алгоритме можно использовать одинаковую величину шага по координатным направлениям, то есть вместо 1 , 2 ,..., n применять . 118 Существует модификация метода, где при исследующем поиске и поиске по образцу используется одномерная минимизация. Тогда если функция f ( x ) дифференцируема, метод сходится к стационарной точке. 8.3.2. Метод деформируемого многогранника В основу метода деформируемого многогранника (метода НелдераМида) положено построение последовательности систем n 1 точек xi (k ), i 1,..., n 1 , которые являются вершинами выпуклого многогранника. Точки системы xi (k 1), i 1,..., n 1 , на k 1 итерации совпадают с точками системы xi (k ), i 1,..., n 1 , кроме i h , где точка x h (k ) - наихудшая в системе xi (k ), i 1,..., n 1 , то есть f ( x h (k )) max f ( xi (k )) . Точка 1 i n 1 x h (k ) заменяется на другую точку по специальным правилам, описанным ниже. В результате многогранники деформируются в зависимости от структуры линий уровня целевой функции, вытягиваясь вдоль длинных наклонных плоскостей, изменяя направление в изогнутых впадинах и сжимаясь в окрестности минимума. Построение последовательности многогранников заканчивается, когда значения функции в вершинах текущего многогранника отличаются от значения функции в центре тяжести системы xi (k ), i 1,..., n 1 ; 0. i h не более чем на Алгоритм 1. Задать координаты вершин многогранника x1,..., x n 1 ; параметры отражения , сжатия , растяжения ; число 0 для остановки алгоритма. Положить k 0 (последующие пункты 2-6 соответствуют текущему номеру k системы точек). 2. Среди вершин ―наилучшую‖ x l и ―наихудшую‖ x h , соответствующие минимальному и максимальному значениям функции: l j h j f (x ) f (x ) max f ( x ) , min f ( x ) ; j 1,...,n 1 j 1,...,n 1 а также точку x S , в которой достигается второе по величине после максимального значения функции f ( x S ) . 3. Найти ―центр тяжести‖ всех вершин многогранника за исключением наихудшей x h : xn 2 1 n n 1 j 1 xj xh 1n 1 j x . nj 1 j h 4. Проверка условия окончания. 119 4а. Если 1 n 1 n 1j f (x j ) 1 2 2 f ( xn 2 ) , процесс поиска можно 1 завершить и в качестве приближенного решения взять наилучшую точку те* l x . кущего многогранника: x 4б. Если , продолжать процесс. 5. Выполнить операцию отражения наихудшей вершины через центр тяжести x n 2 (рис. 8.6, а): xn 3 xn 2 ( xn 2 xh ) . 6. Проверить выполнение условий: 6а. Если f ( x n 3 ) f ( xl ) , выполнить операцию растяжения (рис.8.6, б): xn 4 xn 2 ( xn 3 xn 2 ) . Найти вершины нового многогранника: если f ( x n 4 ) f ( xl ) , то вершина x h заменяется на x n 4 (при n 2 многогранник будет содержать вершины x1, x3 , x6 ). Затем следует положить k k 1 и перейти в п. 2; f ( x n 4 ) f ( xl ) , то вершина x h заменяется на x n 3 (при n 2 многогранник будет содержать вершины x1, x3 , x6 ). Далее следует положить k k 1 и перейти к шагу 2; 6б. Если f ( x S ) f ( x n 3 ) f ( x h ) , то выполнить операцию сжатия (рис.8.6, в): xn 5 xn 2 ( xh xn 2 ) . Следует заменить вершину x h на x n к п. 2 (при n 5 , положить k k 1 и перейти 2 многогранник будет содержать вершины x1, x3 , x7 ); 6в. Если f ( xl ) f ( x n 3 ) f ( x S ) , то вершина x h заменяется на x n 3 . Затем следует положить k k 1 и перейти к п. 2; f ( x n 3 ) f ( x h ) , выполняем операцию редукции 6г. Если (рис 8.6, г). Формируется новый многогранник с уменьшенными вдвое сторонами и вершиной x l : x j xl 0.5( x j xl ), j 1,..., n 1 . При этом следует положить k k 1 и перейти к п. 2. 120 x2 x x1 2 x xn x1 x1 xn 2 x4 x3 xn 1 x n 3 x n 3 x1 xn x f ( xn 4 ) x4 x3 xn 3 xn x5 1 xs 1 xs f ( x1 ) f ( x1 ) x6 x2 x1 2 4 x7 x1 xn xn f ( xn 4 ) 5 5 5 x4 x3 xn xh 2 x1 xs x x2 xh h xh x1 xn 2 x4 x3 xn 3 xn 1 xs x5 Рис. 8.6 – Геометрическая интерпретация метода деформируемого многогранника В литературе авторы в результате вычислительных экспериментов рекомендуют различные значения параметров. Так Нелдер и Мид рекомендуют использовать параметры Павиани – 1; 0.5; 2; Паркинсон и Хатчинсон – 1; 0.4 0.6; 2.8 3; 2; 0.25; 2.5 . В последнем случае в рамках операции фактически выполняется растяжение. 8.3.3. Метод сопряженных направлений В методе сопряженных направлений (метод Пауэлла) используется факт, что минимум квадратичной функции может быть найден не более чем за n шагов при условии, что поиск ведется вдоль сопряженных относительно матрицы Гессе направлений. Так как достаточно большой класс целевых функций может быть представлен в окрестности точки минимума своей квадратичной аппроксимацией, описанная идея применяется и для неквадратичных функций. Задается начальная точка и направления d1, d2 ,..., dn , совпадающие с координатами. Находится минимум f ( x ) при последовательном движении по (n 1) направлениям с помощью одного из методов одномерной минимизации. При этом полученная ранее точка минимума берется в качестве исходной для поиска по следующему направлению, а направление dn используется как при первом (d0 dn ) , так и при последнем поиске. Находится новое направление поиска, сопряженное с dn . Оно проходит через точки, полученные при первом и последнем поиске. Заменяется d1 на d 2 , d 2 на d 3 и т.д. Направление dn заменяется сопряженным направлением, после 121 чего повторяется поиск по (n 1) направлениям, уже не содержащим старого направления d 1 . Для квадратичных функций последовательность n 2 одномерных поисков приводит к точке минимума (если все операции выполнены точно). Построение сопряженного направления для квадратичной функции при n 2 изображено на рис. 8.7. Оно проходит через точки 1 и 3. x2 x0 2 1 3 x 4 d2 1 1 x1 d1 Рис. 8.7 – метода сопряженных направлений Геометрическая интерпретация Алгоритм 1. Задать начальную точку x 0 , число начальные направления поиска 1 d1 dn , i 0, y 0 Положим d0 2. Найти yi ! yi , d2 ... x0 , k 1 , d3 ... 0 для окончания алгоритма, . ... 1 0. ti di , где шаг ti находится в результате поиска ми- нимума функции f ( yi ti di ) по t i одним из методов одномерной минимизации. 3. Проверить выполнение условий: 3а. Если i n 1 , положить i i 1 и перейти к п. 2. 3б. Если i n 1 , проверить успешность поиска по n первым направx* y n , иначе положить лениям. Если y n y 0 , то поиск завершить: i i 1 и перейти к п. 2. 3с. Если i n , проверить успешность поиска по n последним направлениям. Если y n 1 y1 , поиск завершить: x* y n 1 , иначе перейти к п. 4 для построения сопряженного направления. 122 4. Положить x k 1 4а. Если xk 1 xk yn 1 и проверить условие окончания: то поиск заверить: x* xn 1 . положить d0 dn y n 1 y1 (новое направление); di di 1, i 1,..., n 1 (исключается старое направление). Если при этом rang (d1,..., dn ) n , то новая система направлений линейно независима. В этом случае положить k 1 и перейти к п. 2. di di , i 0,1,..., n; k k 1, i 0, y x 4б. Иначе Если rang (d1,..., dn ) n , то новая система направлений линейно зависима. Тогда следует продолжить поиск в старых направлениях. Для этого положить: di di , i 0,1,..., n; k k 1, i 0, y 0 x k 1 и перейти к п. 2. 8.3.4. Адаптивный метод случайного поиска Задается начальная точка x 0 . каждая последующая точка находится по формуле x k 1 x k tk k где tk 0 – величина шага; k – случайный вектор единичной длины, определяющий направление поиска; k – номер итерации. На текущей итерации при помощи генерирования случайных векторов k получаются точки, лежащие на гиперсфере радиуса tk с центром в точке k x (рис. 8.8). Если значение функции в полученной точке не меньше, чем в центре, шаг считается неудачным (точки y1, y 2 при поиске из x0 ; y1, y 3 при поиске из x 1 ). Если число неудачных шагов из текущей точки достигает некоторого числа М, дальнейший поиск продолжается из той же точки, но с меньшим шагом до тех пор, пока он не станет меньше заранее заданной величины R . Если же значение функции в полученной точке меньше, чем в центре, шаг считается удачным и в найденном направлении делается увеличенный шаг, играющий роль ускоряющего шага (как при поиске по образцу в методе конфигураций). Если при этом значение функции снова меньше, чем в центре, направление считается удачным и дальнейший поиск продолжается из этой точки (точки z 3 x1 при поиске из x 0 , z 4 x 2 при поиске из x 1 ). Если же значение функции стало не меньше, чем в центре, направление считается неудачным и поиск продолжается из старого центра (в точке y 2 при поиске из x1 функция меньше, чем в x 1 , а в точке z 2 уже не меньше, поэтому направление ( z 2 x1) неудачное). 123 x z3 x1 y1 z4 y3 y4 y2 y3 x2 z2 y1 x0 y2 Рис. 8.8 – Геометрическая адаптивного метода случайного поиска Алгоритм интерпретация 1. Задать начальную точку x 0 , коэффициенты расширения 1 и сжа1 , М – максимальное число неудачно выполненных испытаний на тия 0 текущей итерации, t0 1 – начальную величину шага, R – минимальную величину шага, N – максимальное число итераций. Положить k 0, j 1 . 2. Получить случайный вектор j ( 1j ,..., nj ) , где i j – случайная величина, равномерно распределенная на интервале [ 1, 1] . 3. Вычислить y j x k j tk j . 4. Проверить выполнение условий. 4а. Если f ( y j ) f ( x k ) , шаг удачный. Положить z j Определить, является ли текущее направление y j если xk 1 z j , tk 1 f (z j ) tk , k xk (y j xk ) . x k удачным: f ( x k ) , направление поиска удачное. Положить k 1 и проверить условие окончания. Если k N, положить j 1 и перейти к п. 2. Если k N , поиск завершить: x* x k ; если f ( z j ) f ( x k ) , направление поиска неудачное, перейти к п. 5; 4б. Если f ( y j ) f ( x k ) , шаг неудачный и перейти к п. 5. 5. Оценить число неудачных шагов из текущей точки: 5а. Если j M , следует положить j j 1 и перейти к п. 2; 5б. Если j M , проверить условие окончания: 124 R , процесс закончить: x* x k , f ( x* ) f ( xk ) ; R , положить tk tk , j 1 и перейти к п. 2. если tk если tk Величина i j , равномерно распределенная на интервале [ 1, 1] , генерируется обычно с помощью датчиков псевдослучайных чисел на ЭВМ. Вырабатывается случайная величина i j , равномерно распределенна на [0, 1] , а затем используется линейное преобразование: i j 2 ij 1 . В литературе рекомендуют следующие параметры алгоритма: =1,618; = 0,618; M 3n . При =1 точка z j на шаге 4 совпадает с y j , то есть аналог поиска по образцу не производится. Начальный шаг t0 R можно задать произвольно. Если выполнено условие окончания tk R , то в качестве ответа можно использовать любую точку внутри шара с радиусом tk и центром в точке x k . 8.4 Методы первого порядка Постановка задачи. Пусть функция f ( x ) , ограниченна снизу на множестве Rn и имеет непрерывные частные производные во всех его точках. Требуется найти локальный минимум функции f ( x ) на множестве допусти- Rn , мых решений X f ( x* ) min f ( x) то есть найти такую точку x* Rn , что x Rn 8.4.1 Метод градиентного спуска с постоянным шагом Стратегия решения задачи состоит в построении последовательности точек {x k }, k 0,1,..., таких, что f ( x k 1) f ( x k ), k 0,1,... . Точки последовательности {x k } вычисляются по правилу xk 1 xk tk f ( x k ), k 0,1,..., где точка x 0 задается пользователем; f ( x k ) - градиент функции f ( x ) , вычисленный в точке x k ; величина шага tk задается пользователем и остается постоянной до тех пор, пока функция убывает в точках последовательности, что контролируется путем проверки выполнения условия f ( x k 1) f ( x k ) 0 или f ( x k 1) последовательности f ( xk ) 1 , где 1 2 f ( xk ) , 0 f ( x) {x k } заканчивается в точке 1 . Построение xk , для которой - заданное малое положительное число, или k M , где M - предельное число итераций, или при двукратном одновременном выпол- 125 нении двух неравенств x k 1 x k 2 , f ( x k 1) f (xk ) 2 где 2 - малое положительное число. Алгоритм 1. Задать x 0 , 0 0 , M - предельное число итераций. f ( x) f ( x) T ,..., ) . Найти градиент функции в произвольной точке f ( x ) ( x1 xn 2. Положить k 0 . 3. Вычислить f ( x k ) . 4. Проверить выполнение критерия окончания f ( x k ) 1. 1, 1 0, 2 4а. Если критерий выполнен, расчет закончен, x* x k . 4б. Если критерий не выполнен, то перейти к п. 5. 5. Проверить выполнение неравенства k M . 5а. Если неравенство выполнено, то расчет окончен, x* 5б. Если нет, то перейти к п. 6. 6. Задать величину шага tk . xk . 7. Вычислить x k 1 x k tk f ( x k ) . 8. Проверить выполнение условия f ( x k 1) f ( x k ) 0 (или f ( x k 1) f ( xk ) f ( xk ) 2 ). 8а. Если условие выполнено, то перейти к п. 9. 8б. Если условие не выполнено, положить tk 9. Проверить выполнение условий xk 1 xk 2 , f ( x k 1) f ( x k 1) 2 tk и перейти к п. 7. 2 : 9а. Если оба условия выполнены при текущем значении k и k k 1, то расчет окончен, x* x k 1 ; 9б. Если хотя бы одно из условий не выполнено, положить k k 1 и перейти к п 3. Геометрическая интерпретация метода для n 2 приведена на рис. 8.7. 126 C1 C2 C3 f ( x) C1 x0 x f ( x ) C3 f ( x0 ) f ( x ) C2 Рис. 8.9 – Геометрическая интерпретация метода градиентного спуска с постоянным шагом Метод градиентного спуска гарантирует сходимость последовательности к точке минимума для сильно выпуклых {x k } функций. Итерационный процесс подбора удачной величины tk отражается в индексации п. 7 и 8. Первый индекс совпадает с номером k, а второй с числом делений текущей величины tk пополам. Оценки скорости сходимости получены только для сильно выпуклых функций, когда последовательность {x k } сходится к точке минимума f ( x ) со скоростью геометрической прогрессии: f ( xk ) f ( x* ) q k ( f ( x0 ) f ( x* )), xk x* C ( q )k где q ( 0,1), C 0 - константы. 8.4.2 Метод наискорейшего градиентного спуска Стратегия поиска x состоит в построении последовательности точек xk , k 0,1,..., таких, что f ( x k 1) f ( x k ), k 0,1,... . Точки последователь- ности {x k } вычисляются по правилу x k 1 x k tk f ( x k ), где точка x 0 задается пользователем; величина шага tk определяется для каждого значения k из условия (tk ) f ( x k tk f ( x k )) min . tk Решение задачи может осуществляться с использованием необходимоd 0 с последующей проверкой достаточного услого условия минимума dt k d2 0 . Такой путь может быть использован либо при доставия минимума 2 dt k точно простой минимизируемой функции (tk ) , либо при предварительной аппроксимации достаточно сложной функции (tk ) f ( xk tk f ( x k )) поли- 127 номом P(tk ) (как правило, второй или третьей степени), и тогда условие d dР d2 0 замещается условием 0 , а условие 0 - условием dt k dt k dt 2 k d 2Р 0. dt 2 k Другой путь решения задачи связан с использованием численных методов, когда ищется min (tk ) min f ( x k tk f ( x k )) . Границы интервала t a;b t a;b [a , b] задаются пользователем. При этом степень близости найденного значеd 0, ния tk к оптимальному значению tk* , удовлетворяющему условиям dt k d2 0 , зависит от задания интервала [a , b] и точности методов одномерной dt 2 k минимизации. Построение последовательности {x k }, k 0,1,..., заканчивается в точке x k , для которой f ( xk ) 1 , где 1 – заданное число, или, если k M, M – предельное число итераций, или при двукратном одновременном выполнении неравенств xk 1 xk 2 , f ( x k 1) f (xk ) 2 , где 2 – малое поло- жительное число. Алгоритм 1. Задать x 0 , 1 , 2 0 , предельное число итераций М. найти градиент функции в произвольной точке f ( x) df ( x) df ( x) ,..., dx1 dxn T . 2. Положить k 0 . 3. Вычислить f ( x k ) . 4. Проверить выполнение критерия окончания f ( xk ) 4а. Если критерий выполнен, то x* x k ; 4б. Если критерий не выполнен, то перейти к п. 5. 5. Проверить выполнение неравенства k M . 5а. Если неравенство выполнено, то x* x k ; 5б. Если нет, то перейти к п. 6. 6. Вычислить величину шага tk* из условия (tk ) f ( xk tk f ( x k )) 7. Вычислить x k 1 xk min . tk tk f ( x ) . 1 . 128 8. Проверить выполнение условий xk 1 xk 2 , f ( x k 1) f (xk ) 2 . 8а. Если оба условия выполнены при текущем значении k и k k 1, то расчет окончен, x* x k 1 ; 8б. Если хотя бы одно из условий не выполнено, то положить k k 1 и перейти к п. 3. Геометрическая интерпретация метода для n 2 приведена на рис. 8.10. x2 C1 f ( x) C1 C2 x0 f ( x0 ) x1 f ( x ) C2 x1 Рис. 8.10 – Геометрическая интерпретация метода наискорейшего градиентного спуска Метод наискорейшего спуска гарантирует сходимость последовательности {x k } к точке минимума для сильно выпуклых функций. Оценки скорости сходимости получены только для сильно выпуклых функций, когда последовательность {x k } сходится к точке минимума функции f ( x ) со скоростью геометрической прогрессии (линейная сходимость): k 1 M m k x xk x x* , где M и m – оценки наибольшего и наименьшеM m го собственных значений матрицы Гессе функции f ( x ) . 8.4.3 Метод Гаусса-Зейделя Стратегия метода Гаусса-Зейделя состоит в построении последовательности точек xk , k 0,1,..., таких, что f ( x k 1) следовательности {x k } вычисляются по правилу x jk 1 x jk tk f ( x) xk 1 ek x x jk 1 , f ( x k ), k 0,1,... Точки по- 129 где j – номер цикла вычислений, j 0,1,2,... ; k – номер итерации внутри цикла, k 0,1,..., n 1 ; ek 1 – единичный вектор, (k 1) – я проекция которого равна 1; точка x00 задается пользователем, величина шага tk выбирается из условия f x jk (tk ) tk f ( x) xk 1 * ek min . 1 tk x x jk Эта задача является задачей одномерной минимизации функции (tk ) f x jk f ( x) xk 1 tk ванием условий d dtk * ek x x 2 1 и может быть решена либо с использо- jk d 0 , либо численно с использованием методов dt 2k одномерной минимизации, как задача (tk ) min . 0, tk a;b d 0 имеет высокую степень и корни его трудно dtk определить, можно аппроксимировать функцию (tk ) полиномом P(tk ) Если уравнение dР d 2Р 0, 2 0. второй или третей степени и определить tk из условий dtk dt k При численном решении задача определения величины шага степень близости найденного значения tk к оптимальному значению tk *, удовлетвоd d2 0, 2 ряющему условиям 0 , зависит от задания интервала [a , b] и dtk dt k точности методов одномерной минимизации. Легко видеть, что при фиксированном j за одну итерацию с номером * k изменяется только одна проекция точки x jk , имеющая номер k 1 , а в течение всего цикла с номером j , то есть начиная с k 0 и кончая k n 1 , изменяются все n проекций точки x j 0 . После этого точке x jn присваивается номер x j 1,0 ( j 1) м цикле. и она берется за начальную точку для вычислений в Расчет заканчивается в точке x jk при выполнении по крайней мере одного из трех критериев окончания счета: f ( xk ) кратного выполнения неравенств x jk 1 x jk Здесь 1 , 2 1 2 , или k , f ( x jk 1) M , или двуf ( x jk ) 2 . – малые положительные числа, M – предельное число циклов 130 итераций. Вопрос о том, является ли точка x jk найденным приближением искомой точки, решается путем проведения дополнительного исследования. Алгоритм 1. Задать x00 , 1 0, 2 0 ; предельное число M циклов счета, кратное n , где n – размерность вектора x . Найти градиент f ( x ) . 2. Задать номер цикла j 0 . 3. Проверить условие j M : 3а. Если j M , то расчет окончен и x * x jk ; 3б. Если j M , то перейти к п. 4. 4. Задать k 0 . 5. Проверить условие k n 1: 5а. Если k n 1, то перейти к п. 6; 5б. Если k n , то положить j j 1 и перейти к п. 3. f ( x jk ) . 6. Вычислить f ( x jk ) 7. Проверить выполнение условия 1 : 7а. Если условие выполнено, то расчет окончен и x* 7б. Если нет, то перейти к п. 8. 8. Вычислить tk* из условия (tk ) f x jk tk 9. Вычислить x jk 1 f ( x) xk 1 x x x jk tk* * ek min . 1 tk jk f ( x) xk 1 x jk . ek x x 1 . jk 10. Проверить выполнение условий xk 1 xk 2 , f ( x k 1) f (xk ) 2 : 10а. Если оба условия выполнены в двух последовательных циклах с номерами j и j 1 , то расчет окончен, найдена точка x* x jk 1 ; 10б. Если не выполняется хоты бы одно условие, положить k k 1 и перейти к п. 5. Геометрическая интерпретация метода для n 2 приведена на рис. 8.11. 131 x2 x 00 x 01 x 02 x1 Рис. 8.11 – Геометрическая интерпретация метода Гаусса-Зейделя 8.4.4 Метод Дэвидона-Флетчера-Пауэлла Стратегия метода Дэвидона-Флетчера-Пауэлла состоит в построении последовательности xk , k точек 0,1,..., таких, что f ( x k 1) f ( x k ), k 0,1,... . Точки последовательности {x k } вычисляются по правилу x k 1 x k tk Ak f ( x k ) , k 0,1,..., где Ak – матрица размера n n , которая вычисляется по правилу Ak k A 1 Ak Ack , A0 E, T Ak g k k x ( xk ) c x где x k k T xk g 1 k g k T xk , g k gk T Ak g k , k A f xk 1 f xk . Точка x 0 задается пользователем, величина шага tk определяется из условия kt f xk tk Ak f x k min . tk Решение задачи может осуществляться как из условий d2 задачи kt min t k a;b 0, d 2Р 0 , где P(tk ) – полином, аппроксиdt 2 k (tk ) , так и численно, то есть путем поиска решения 0 или из условий dt 2k мирующий функцию d dtk dР dtk 0, методами одномерной минимизации. 132 Построение последовательности {x k } заканчивается в точке x k , для которой f ( xk ) 1 , где – заданное число, или при k 1 M ( M – пре- дельное число итераций), или при двукратном одновременном выполнении двух неравенств x k 1 x k , f ( x k 1) 2 f (xk ) 2 , где 2 – малое по- ложительное число. Алгоритм 1. Задать x0 , диент f ( x) . 0, 1 0, M - предельное число итераций. Найти гра- 2 2. Положить k 0, A0 3. Вычислить f ( xk ) . E. f ( xk ) 4. Проверить критерий окончания 1 : 4а. Если критерий выполнен, то расчет окончен и x* x k . 4б. Если нет, то перейти к п. 5. 5. Проверить условие k M : 5а. Если неравенство выполнено, то расчет окончен и x* x k . 5б. Если нет, перейти при k 0 к шагу 10, а при k 1 к п. 6. 6. Вычислить gk 7. Вычислить xk f xk xk f xk . 1 xk . k T x ( xk ) 8. Вычислить Ak c x 9. Вычислить Ak 1 1 10. Определить d k k T Ak g k g k T gk T Ak g k . k A Ack . Ak f ( x k ) . 11. Вычислить t *k из условия 12. Вычислить x k Ak g k 1 xk kt f xk tk Ak f x k min . tk tk Ak f ( x k ) . 13. Проверить условия x k 1 x k 2 , f ( x k 1) f (xk ) 2 : 13а. Если оба условия выполнены в двух последовательных циклах с номерами k и k 1 , то расчет окончен, найдена точка x* x k 1 ; 13б. Если не выполняется хоты бы одно условие, положить k k 1 и перейти к п. 3. 133 Алгоритм метода Дэвидона-Флетчера-Пауэлла в применении к квадра1 2 ( fx, x) (b, x) с положительно определенной маттичной функции f ( x) 2 рицей Гессе 2 f обеспечивает отыскание минимума x* ( 2 ) 1b не более чем за n шагов. Для неквадратичных функций алгоритм не является конечным и его скорость сходимости зависит от точности решения задачи. 8.5 Методы второго порядка Постановка задачи. Пусть дана функция f ( x ) , ограниченная снизу на множестве x R n и имеющая непрерывные частные производные во всех его точках. Требуется найти локальный минимум функции f ( x ) на множестве допустимых решений X f ( x* ) Rn , т. о. есть найти такую точку x* R n , что min f ( x), f ( x) C 2 . x Rn 8.5.1 Метод Ньютона Стратегия метода Ньютона состоит в построении последовательности точек {x k }, k 0,1,..., таких, что f ( x k 1) f ( x k ), k 0,1,... . Точки последовательности вычисляются по правилу x k 1 x k d k , k 0,1,..., где x 0 - задается пользователем, а направление спуска d k определяется для каждого значения k по формуле dk 2 f ( xk ) 1 f ( xk ) . Выбор d k гарантирует выполнение требования f ( x k 1) f ( x k ) при условии, что 2 f ( x k ) 0 . Формула для d k получена из следующих соображений: 1. Функция f ( x ) аппроксимируется в каждой точке последовательности x k квадратичной функцией Fk f ( x k ) ( f ( x k ), d k ) 0.5(d k , 2 f ( xk ) d k ) . 2. Направление d k определяется из необходимого условия экстремума dFk первого порядка: 0 . Таким образом, при выполнении требования dd k 2 f ( x k ) 0 последовательность является последовательностью точек минимумов квадратичных функций Fk , k 0,1,... (рис. 8.12). Чтобы обеспечить выполнение требования f ( x k 1) f ( x k ), k 0,1,... , даже в тех случаях, когда 134 для каких-либо значений матрица Гессе 2 f ( x k ) не окажется положительно определенной, рекомендуется для соответствующих значений k вычислить точку x k 1 по методу градиентного спуска x величины шага tk из условия f ( x k tk f ( x k )) k 1 xk tk f ( xk ) с выбором f ( xk ) . Построение последовательности {x k } заканчивается в точке x k , для которой f ( xk ) 1 , где - заданное малое положительное число, или 1 при k M ( M - предельное число итераций), или при двукратном одновреk 1 менном выполнении двух неравенств x k 1 x k ) f (xk ) 2 , f (x 2 , где 2 - малое положительное число. y F0 F1 f(x) x2 x x0 x1 x Рис. 8.12 - Геометрическая интерпретация метода Ньютона Алгоритм 1. Задать x 0 , 1 , 2 , M – предельное число итераций. Найти градиент f ( x ) и матрицу Гессе 2 f ( x) . 2. Положить k 0 . 3. Вычислить f ( x k ) . 4. Проверить выполнение критерия окончания f ( xk ) 1 . Если не- равенство выполнено, то расчет окончен и x* x k , в противном случае перейти к п. 5. 5. Проверить выполнение неравенства k M . Если неравенство выполнено, расчет окончен и x* x k , если нет, перейти к п. 6. 6. Вычислить матрицу 2 f ( xk ) . 135 2 7. Вычислить матрицу 8. 2 Проверить 1 f ( xk ) выполнение 2 условия f 1 ( xk ) 1 0. Если f ( xk ) . 9. Определить d k 10. Найти 2 dk . 0 , то перейти к п. 9. В противном случае перейти к п. 10, по- ложив d k dk 1 f ( xk ) f ( xk ) 1 2 f ( xk ) xk точку 1 f ( xk ) . 1 xk tk d k , положив f ( xk ) или выбрав tk из условия f ( x k 1) tk 1, f ( x k ) , если f ( xk ) . 11. Проверить выполнение условий xk 1 xk 2 f ( x k 1) , f ( xk ) 2 : 11а. Если оба условия выполнены при текущем значении k и k k 1, то расчет окончен, x* x k 1 ; 11б. В противном случае положить k k 1 и перейти к п. 3. Сходимость метода Ньютона доказана лишь для сильно выпуклых функций и для достаточно хорошего начального приближения. При практическом использовании метода Ньютона следует: 2 f ( x k ) на выполнение условия а) анализировать матрицу 2 f ( xk ) 0 k на формулу x k 0,1,... и заменять формулу x k 1 1 xk xk 2 f ( xk ) 1 f ( xk ) tk f ( x k ) в случае его невыполнения; б) производить анализ точки x k с целью выяснения, является ли она найденным приближением искомой точки x* . 8.5.2 Метод Ньютона-Рафсона Стратегия метода Ныотона-Рафсона состоит в построении последовательности точек {x k }, k 0,1,..., таких, что f ( x k 1) f ( x k ), k 0,1,... . Точки последовательности вычисляются по правилу xk 1 xk tk 2 f ( xk ) 1 f ( x k ), k 0,1,... , где x 0 задается пользователем, а величина шага tk определяется из условия 136 (tk ) f ( xk tk 2 1 f ( xk ) f ( x k )) min . tk Эта задача может решаться либо аналитически с использованием необd 0 с последующей проверкой достаточного ходимого условия минимума dtk d2 условия min , где интервал 0 , либо численно как задача (tk ) 2 dtk t k [ a ,b ] [a , b] задается пользователем. Если функция (tk ) достаточно сложна, то возможна ее замена полиномом P(tk ) второй или третьей степени и тогда шаг tk может быть опредеdP d 2P 0 при выполнении условия лен из условия 0. 2 dtk dtk При численном решении задачи определения величины шага степень близости найденного значения tk к оптимальному значению tk * , удовлетвоd d2 0, ряющему УСЛОВИЯМ 0 зависит от задания интервала [a , b] и 2 dtk dtk точности методов одномерной минимизации. Построение последовательности {x k } заканчивается в точке x k , для которой f ( xk ) 1 , где – заданное число, или при k 1 M ( M - пре- дельное число итераций), или при двукратном одновременном выполнении двух неравенств x k 1 x k 2 , f ( x k 1) f (xk ) 2 , где 2 - малое поло- жительное число. Алгоритм 1. Задать x 0 , 1 , 2 , M - предельное число итераций. Найти градиент f ( x ) и матрицу Гессе 2 f ( x) . 2. Положить k 0 . 3. Вычислить f ( x k ) . 4. Проверить выполнение критерия окончания f ( xk ) 1 . Если не- равенство выполнено, то расчет окончен и x* x k , в противном случае перейти к п. 5. 5. Проверить выполнение неравенства k M . Если неравенство выполнено, расчет окончен и x* x k , если нет, перейти к п. 6. 6. Вычислить матрицу 2 f ( xk ) . 137 2 7. Вычислить матрицу f ( xk ) 1 . 2 8. Проверить выполнение условия полняется, dk то найти d k f ( xk ) H 1( xk ) f ( xk ) , 1 если 0 . Если условие вынет, то положить f ( xk ) . 9. Определить x k 1 xk tk d k . 10. Найти шаг t *k из условия k 1 k (tk ) f ( xk tk d k ) min . tk * k x tk d . 11. Вычислить x 12. Проверить выполнение неравенств xk 1 xk 2 , f ( x k 1) f ( xk ) 2 . Если оба условия выполнены при текущем значении k и k k 1, то расчет окончен, x* x k 1 , в противном случае положить k k 1 и перейти к п. 3. Сходимость к точке минимума метода Ньютона-Рафсона гарантируется независимо от выбора начального приближения лишь для сильно выпуклых функций. Поэтому при практическом использовании метода НьютонаРафсона следует: а) анализировать матрицу Гессе 2 f ( x k ) на выполнение условия 2 f ( xk ) 0 k 0,1,... и заменять формулу x k 1 xk 2 f ( xk ) 1 f ( xk ) на формулу метода градиентного спуска x k 1 x k tk f ( x k ) в случае его невыполнения; б) производить анализ точки x k с целью выяснения, является ли она найденным приближением искомой точки x* . 138 ЛИТЕРАТУРА 1.Акулич И.Л. математическое программирование в примерах и задачах.–М.: Высшая школа, 1986. 2. Базара М., Шетти К. Нелинейное программирование. Теория и алгоритмы. – М.: Мир, 1982. 3. Банди Б. Методы оптимизации. Вводный курс. – М., Радио и связь, 1989. 4. Банди Б. Основы линейного программирования. – М., Радио и связь, 1989. 5. Вентцель Е.С. Исследование операций. – М.: Советское радио, 1972. 6. Грешилов А.А. Прикладные задачи математического программирования. – М.: Изд-во МГТУ, 1990. 7. Мину М. Математическое программирование. Теория и алгоритмы. – М.: Наука, 1990. 8. Муханова Э.А., Рубенштейн Г.Ш. Математическое программирование. – Новосибирск., Наука СО, 1987. 9. Понтелеев А.В, Летова Т.А. Методы оптимизации в примерах и задачах. – М.: Высшая школа, 2002 . 10. Поляк Б.Г. Введение в оптимизацию. – М.: Наука, 1983. 11. Сухарев А.Г., Тимохов А.В., Федоров В.В. Курс методов оптимизации. – М.: Наука, 1986.