Линейное, нелинейное и квадратичное программирование

С. М. БОРОДАЧЁВ ТЕОРИЯ СИСТЕМНОГО АНАЛИЗА И ПРИНЯТИЯ РЕШЕНИЙ Лекции s-m-borodachev.blogspot.com 1 Оглавление ВВЕДЕНИЕ .................................................................................................................. 4 1. ЭЛЕМЕНТЫ ЛИНЕЙНОГО ПРОГРАММИРОВАНИЯ ..................................... 5 1.1. Канонический вид задачи линейного программирования............................ 6 1.2. Типичные применения линейного программирования................................. 8 Оптимальное использование ресурсов .................................................................. 8 Планирование инвестиций ...................................................................................... 8 Транспортная задача ................................................................................................ 9 Задача о назначениях ............................................................................................. 10 Задача коммивояжёра ............................................................................................ 11 2. НЕЛИНЕЙНОЕ И КВАДРАТИЧНОЕ ПРОГРАММИРОВАНИЕ ................... 14 2.1. Выбор инвестиционного портфеля (задача Марковица) ............................ 16 3. ДИНАМИЧЕСКОЕ ПРОГРАММИРОВАНИЕ .................................................. 18 Уравнение Р. Беллмана.......................................................................................... 20 3.1. Стохастические модели динамического программирования ..................... 21 Задача о наилучшем выборе ................................................................................. 24 4. ПРИНЯТИЕ РЕШЕНИЙ В УСЛОВИЯХ НЕОПРЕДЕЛЁННОСТИ И РИСКА ..................................................................................................................................... 28 4.1. Условия неопределённости............................................................................ 28 Некоторые нестандартные критерии ................................................................... 28 4.2. Условия риска (критерий Байеса – Лапласа) ............................................... 29 Ожидаемая ценность точной информации (EVPI) ............................................. 30 4.3. Статистические решающие правила ............................................................. 31 Минимаксное решающее правило ....................................................................... 32 Байесовское решающее правило .......................................................................... 33 Свойства минимаксных и байесовских стратегий. ............................................ 35 Выборочный контроль качества продукции ....................................................... 38 5. ТЕОРИЯ СИСТЕМ ................................................................................................ 41 5.1. Система как математический объект ............................................................ 41 2 5.2. Математическая классификация систем ...................................................... 43 5.3. Синтез конечных автоматов .......................................................................... 44 5.4. Основные проблемы теории систем ............................................................. 50 5.5. Передаточная матрица системы .................................................................... 51 5.6. Макроэкономические модели ........................................................................ 54 Модель Леонтьева (~1930 г) ................................................................................. 54 Двойственная модель Леонтьева .......................................................................... 56 Дискретно-временная динамическая модель Леонтьева ................................... 57 БИБЛИОГРАФИЧЕСКИЙ СПИСОК...................................................................... 57 3 Введение Любая сфера человеческой деятельности связана с принятием решений. В зависимости от предметной области дисциплины, посвящённые этой проблематике, назывались исследование операций, экономико-математические методы, математические методы в управлении. Для постановки задачи принятия управленческого решения необходимо выполнить два условия: 1. Должна быть возможность выбора, т. е. по крайней мере, 2 варианта действий. 2. Вариант выбирается по определённому принципу. Известны два принципа выбора решения: a. Волевой – применяют при отсутствии формализованных моделей. b. Критериальный – заключается в принятии некоторого(-ых) критерия(-ев) (мерило) и сравнении возможных вариантов по этому критерию. Принятый критерий называют целевой функцией. Вариант, для которого целевая функция принимает наилучшее значение, называют оптимальным и на нём останавливают выбор. Классификация задач принятия решений может быть проведена по нескольким признакам: 1. По степени полноты информации об условиях, в которых принимается решение. a. Детерминированные задачи – относительно каждого действия известно, что оно неизменно приводит к некоторому конкретному исходу, внешние факторы также известны. b. Вероятностные задачи (или задачи с риском) – включают в своей постановке параметры, являющиеся случайными величинами, для которых известны или хотя бы экспертно могут быть оценены распределения вероятностей. c. Задачи в условиях неопределённости – результаты ваших действий и/или внешние факторы могут быть различными, но их вероятности совершенно неизвестны или не имеют смысла. 2. По используемым математическим методам для их решения (хотя, как правило, степень полноты информации и определяет используемый метод). 4 1. Элементы линейного программирования Здесь слово программирование используется в смысле определения программы действий, программы выпуска изделий, т.е. как планирование. Пример 1.1[9] Фирма производит скороварки (А) и кофеварки (Б). Основные операции, и производственные мощности (если производить что-то одно, в шт. за неделю): Операция Штамповка Отделка Сборка А Сборка Б Прибыль от шт. Вид продукции А 25000 33333 22500 – 15 ₽ Б 35000 16667 – 15000 12.5 ₽ Какова должна быть программа выпуска на следующую неделю x1 – скороварок, шт., x2 – кофеварок, шт., чтобы прибыль была максимальной? Запишем функцию прибыли и ограничения по мощностям. Пусть полная производственная мощность цеха штамповки = 1 и т. д. f ( x1 , x2 ) = f ( x) = 15 x1 + 12.5 x2 → max x2  x1  25000 + 35000 ≤ 1   x1 + x2 ≤ 1   33333 16667  x1 x2 ≤ 1, ≤1  15000  22500  x1 ≥ 0, x2 ≥ 0 Изобразим ограничения и целевую функцию на рисунке. 5 С учётом знаков неравенств, получаем, что каждая точка допустимого плана выпуска  x1  x =   должна находиться внутри заштрихованной области или на её границе.  x2  Рассмотрим семейство прямых 15 x1 + 12.5 x2 = f (при различных значениях f). Нормальный  15  вектор n =   у них одинаков, значит они параллельны. Перемещаем прямую вдоль n для  12.5  получения максимального значения прибыли f (оставаясь в допустимой области!). Находим,  x1* = 20370  что максимум достигается в угловой точке x* =   . При этом фирма получит  x2 * = 6481  максимальную прибыль f max = 15 x1 * +12.5 x2 * = 386562.5 руб. Кстати, использование мощностей по сборке кофеварок будет значительно ниже 100%. Поэтому разумно часть работников отсюда перевести на другие участки. 1.1. Канонический вид задачи линейного программирования В общем случае число переменных может быть произвольным, x1 , x2 ,..., xn и подчинены они m ограничениям вида 6 a11 x1 + a12 x2 + ... + a1n xn ≤ b1  ... a x + a x + ... + a x ≤ b mn n m  m1 1 m 2 2 (1) Введя вектор переменных (плана)  x1    x x = 2,  ..     xn  вектор ограничений (2)  b1    . b =   , матрицу коэффициентов  .     bm   a11 ... a1n    A =  ... ...  , a   m1 ... amn  систему ограничений (1) можно записать в матрично-векторной форме Ax ≤ b . Целевую функцию можно записать f ( x ) = c1 x1 + ... + cn xn = c T x , тогда задача линейного программирования примет вид x * = arg max c T x ; x∈E (3) E = {x : Ax ≤ b , x ≥ 0}, где E – допустимая область, любая x ∈ E – допустимая точка или план. Пример 1.2 Понятие argmax изображено на рисунке. Ограничения в виде неравенств можно свести к ограничениям-равенствам вводя дополнительные вспомогательные переменные, поиск максимума целевой функции можно заменить поиском минимума функции, умноженной на -1. Тогда приходим к определению. 7 Канонический вид задачи линейного программирования – найти вектор плана x * , для которого целевая функция f ( x ) = c T x достигает минимума на множестве точек E0, удовлетворяет системе ограничений Ax = b , где b ≥ 0 . x * = arg min c T x ; x∈E0 (4) E0 = {x : Ax = b , x ≥ 0}. Вспомним метод Гаусса решение систем линейных уравнений. В нём расширенная матрица системы приводится к ступенчатой форме, определяются базисные и свободные переменные. Если все свободные переменные положить равными нулю, то полученное частное решение системы называется базисным решением. Теорема Каждое допустимое базисное решение системы Ax = b есть угловая точка области E0. Обратно: каждая угловая точка E0 – допустимое базисное решение. 1.2. Типичные применения линейного программирования Оптимальное использование ресурсов Предприятие выпускает n видов изделий. Для их производства используются m видов ресурсов (разное сырье, людские ресурсы, финансовые ресурсы и т. п.). Эти ресурсы ограничены, их запасы составляют в планируемый период b1, b2, …, bm условных единиц. Известны также технологические коэффициенты: aij – сколько единиц i-го ресурса требуется в производстве единицы j-го вида продукции. cj – прибыль от реализации единицы j-го вида продукции. Требуется составить такой план выпуска продукции ( x1* – единицы изделий первого вида, x2* – единицы изделий второго вида и т. д.), при котором прибыль предприятия была бы наибольшей. Пусть будем выпускать план (2), тогда a12 x2 – единиц 1-го ресурса пойдёт на всю программу выпуска 2-го вида продукции. Аналогично с другими видами продукции, поэтому общий расход a11 x1 + a12 x2 + ... + a1n xn 1-го ресурса не должен превышать его запас. Записав это для всех ресурсов, имеем (1). Суммарная прибыль – c T x , в результате получаем задачу линейного программирования (3). Планирование инвестиций Пример 1.3 8 Плановый горизонт – 3 года. Моменты времени их начала и окончания: k = 0, k = 1, k = 2, k = 3. Существуют 5 инвестиционных проектов с характеристиками: Момент 1 2 3 Приход, руб., в момент k на вложенный рубль A B C D -1 -1 -1 0.3 -1 1.1 1 0.3 1 1.75 E -1 1.4 Отрицательный приход означает инвестирование. В любой момент деньги, вернувшиеся из какого-либо проекта, можно вкладывать в другие, открывающиеся в тот же момент. Всегда деньги можно положить в банк под 6% годовых. Максимальное вложение в проект A – 500000 рублей. В момент k = 0 имеется 1000000 рублей для инвестирования. Цель: составить план инвестиций, максимизирующий сумму денег в k = 3. Обозначим ak, bk, ck, dk, ek, sk – инвестиции в проекты и в банк в момент k. Ограничения: a0 ≤ 500000 a + c + d + s = 1000000 (k = 0)  0 0  b + s = 0.3 a + 1.1c0 + 1.06 s0 (k = 1) 1 1 e2 + s2 = a0 + 0.3b1 + 1.06 s1 (k = 2) Целевая функция (приход в момент k = 3) f = b1 + 1.75d0 + 1.4e2 + 1.06s2 → max . Имеем ЗЛП с 8 переменными, 4 ограничениями. Найдём решение fmax = =1797600 руб. Т. е. прирост 79.76%. Инвестиции a0 = 500000, d0 = 500000, e2 = =659000, s1 = 150000. Остальные = 0. Транспортная задача Имеется m пунктов отправления однородного, бесконечно делимого груза с запасом в каждом из них ai и n пунктов потребления этого груза с потребностями bj в каждом. Затраты на перевозку единицы груза из i-го пункта отправления в j пункт потребления известны и равны cij. Задача: как спланировать перевозки, т.е. какие количества xij груза перевести от i-го поставщика к j-му потребителю, чтобы все потребности были удовлетворены, а суммарные затраты на перевозки были минимальны. Поскольку нельзя вывезти больше запаса, постольку 9  x11 + x12 + ... + x1n ≤ a1  .................................  x + x + ... + x ≤ a mn m  m1 m 2 . (5) Все потребности должны быть удовлетворены:  x11 + x21 + ... + xm1 = b1  ................................. .  x + x + ... + x = b mn n  1n 2 n (6) Суммарные затраты на перевозки m n f = c11 x11 + c12 x12 + ... + cmn xmn =  cij xij → min . (7) i =1 j =1 Уравнения (5)–(7) – открытая транспортная задача, не весь груз должен быть вывезен. Замкнутая транспортная задача – сумма запасов равна сумме потребностей. Задача о назначениях Пример 1.4 Пусть имеется 5 исполнителей и 5 работ Исполнители 1 2 3 4 5 Работы 1 2 6 4 4 5 2 9 8 6 2 3 3 2 7 5 7 9 4 7 6 3 3 5 5 1 1 1 1 1 В таблице указано, какое время cij , ч. требуется каждому исполнителю для выполнения каждой работы. Задача: назначить исполнителей на работы так, чтобы суммарное время выполнения всех работ было минимально. Обозначим xij – назначенность i-го исполнителя на j-ю работу. 1, i − й назначен на j − ю работу xij =  . 0, i − й не назначен на j − ю работу Каждый работник должен быть назначен один раз:  x11 + x12 + ... + x1n = 1  ................................. .  x + x + ... + x = 1 nn  n1 n 2 (8) 10 Каждая работа должна быть выполнена один раз:  x11 + x21 + ... + xn1 = 1  ............................... .  x + x + ... + x = 1 2n nn  1n (9) n n f = c11 x11 + c12 x12 + ... + cnn xnn =  cij xij → min . i =1 j =1 Видно, что задача о назначениях есть частный случай замкнутой транспортной задачи. Задача коммивояжёра Постановка задачи: есть N+1 городов, стоимости переезда из любого i-го города в любой jй известны (cij). Коммивояжёр должен побывать в каждом городе один раз и вернуться в исходный город маршрута (0-й), затратив минимум средств на переезды. Сколько всего возможно маршрутов? Ясно, что каждый из них можно записать в виде последовательности проезжаемых городов, где первый и последний члены «0». Все они отличаются перестановкой N элементов (1, 2, …, N) между этими «0». Известно, что число перестановок из N элементов равно N! Например, при 10 городах 9! = 362880 – число возможных маршрутов. Выбрать самый дешёвый путём перебора сложно. При N > 70 полный перебор вариантов невозможен, ни на каких мыслимых компьютерах за миллиарды лет! К задаче коммивояжёра сводятся разнообразные задачи производственного менеджмента. Пример 1.5 Выбор порядка обработки деталей на автоматической линии, при котором затраты времени на наладку были бы минимальны. Поскольку время наладки очередной детали зависит от того, какая деталь обрабатывалась до неё, порядок обработки деталей влияет на суммарные потери времени на наладку станков линии. Время обработки деталей неизменно и не влияет на решение задачи. Т. о. матрица (cij) – матрица времён переналадки. Пример 1.6 Строительство n объектов. Строительство каждого состоит из m последовательных стадий (земляные работы, фундамент, кладка стен и т. п.). Продолжительность каждой стадии на каждом объекте известна. Мощность строительной организации не позволяет вести один и тот же вид работ одновременно на нескольких объектах. Считая, что работа на каждом объекте должна продолжаться непрерывно до окончания его строительства, требуется определить сроки начала строительства каждого объекта (а значит и последовательность объектов) так, чтобы суммарный срок строительства всех объектов был минимальным. Для математической постановки задачи коммивояжёра введём переменные xij : 1, маршрут включает переезд из i в j xij =  . 0, маршрут не включает переезд из i в j 11 N Целевая функция N  c x i =0 j =0 ij ij → min . В силу замкнутости маршрута, коммивояжёр из каж- дого города обязательно в какой-то уедет (8), и в каждый из какого-то приедет (9). Чтобы исключить отсутствие переезда ( xii = 1 ), примем все cii = +∞ . Т. о. задача коммивояжёра имеет тот же вид, что и задача о назначениях. Но этих условий здесь недостаточно, так как может получиться решение, состоящее из нескольких несвязанных замкнутых маршрутов. Как изображено на рисунке. В общем случае задачи коммивояжёра не существует метода её решения лучшего, чем полный перебор вариантов. Разработаны разнообразные методы поиска решения (возможно, не самого лучшего) с сокращением числа перебираемых вариантов. Рассмотрим один из них на примере. Пример 1.7 Пусть N = 4 (всего пять городов). Матрица расстояний между ними (в км): 300 250 200 400 300 500 350 600 250 500 250 200 200 350 250 250 400 600 200 250 Будем поэтапно рассматривать под маршруты, начинающиеся в исходном городе и включающие всё большее число городов. Начнём с под маршрутов из 3 участков. В пункт 1 при этом можно попасть через всевозможные последовательности 2-х промежуточных городов 6ю способами: Вариант движения 0231 Расстояние, км 250+250+350=850 12 Перспективно или нет да 0321 0241 0421 0341 0431 … 0124 0214 0134 0314 0234 0324 200+250+500=950 1050 1100 1050 1000 … 1000 1350 900 1150 750 650 нет да нет нет да … да нет да нет нет да Аналогично во все другие конечные пункты. Рассмотрим первые две строки. Оба варианта оканчиваются в одном и том же пункте (1) и имеют уже пройдёнными одно и то же множество (2, 3). В дальнейшем оба могут развиваться одинаково, но ясно, что проигрыш в 100 км делает второй бесперспективным, по сравнению с первым при любом его дальнейшем развитии. Поэтому отбросим его, а значит и все его возможные продолжения, число которых может быть огромным. Теперь строим под маршруты из 4 участков. Нужно попасть в каждую точку через всевозможные последовательности 3-х промежуточных городов кратчайшим путем при каждом составе промежуточных городов. Поэтому, во всех перспективных вариантах 1-го этапа добавляем в конец отсутствующую в нем точку. Из всех оканчивающихся в одной точке выбираем кратчайший. Например, см. выбор из подчёркнутых для точки 4. И будем знать кратчайшие маршруты попадания во все точки, при всевозможных наборах промежуточных 3-х! Вариант движения 02314 Расстояние, км 1450 13 Перспективно или нет нет 02413 04312 01324 01423 03421 01234 01432 02431 01243 01342 03241 нет нет да нет нет нет нет да да да нет 1400 1500 1100 1350 1150 1350 1400 1050 1250 1100 1250 Эти варианты дополняем пунктом возвращения (0). Варианты движения 013240 024310 012430 013420 Расстояние, км 1500 1350 1450 1350 Видно, что имеется два оптимальных маршрута следования 024310 и 013420, имеющие минимальную длину, равную 1350 км. 2. Нелинейное и квадратичное программирование Задача нелинейного программирования: x * = arg min F ( x ); x∈M M = {x : f j ( x ) ≤ 0, j = 1,..., m; x ≥ 0}, где F( x ), fj( x ) – произвольные нелинейные функции. Пример 2.1 При каких размерах канал с сечением большим или равным S в форме равнобочной трапеции будет иметь наименьшую поверхность смачивания (наименьшие потери)? 14 Поверхность смачивания будет минимальна b ⋅ h + d ⋅ h ≥ S u = b + 2l = b + 2 h2 + d 2 → min при сохранении сечения  . b, h, d ≥ 0 Эта задача нелинейного программирования имеет решение h* = 4 S S 2 S , , d* = , b* = 4 3 4 3 3 3 3 α* = 60°. Задача квадратичного программирования (частный случай задачи нелинейного программирования, когда целевая функция квадратична, а ограничения линейны): x * = arg min(c T x + x T Kx ) x∈M M = { x : Ax ≤ b , x ≥ 0} . K – симметричная матрица размера [n*n]. Для решения задач нелинейного и квадратичного программирования широко применяются численные методы нахождения решения (приближенного решения). Например, в MathCAD есть функция Minimize( f, x, y, …), (Maximize( f, x, y, …)) дающая набор значений неизвестных x, y, …, минимизирующих функцию f , при ограничениях (уравнениях и неравенствах), указанных в блоке Given. Т. е. это фактически argmin f(x). Пример 2.2 Решить задачу F ( x ) = −2 x1 − 4 x2 + x12 + 2 x22 → min;  x1 + 2 x2 ≤ 8  2 x1 − x2 ≤ 12. x , x ≥ 0  1 2 Решение в MathCAD ORIGIN:=1 F ( x ) := −2 ⋅ x1 − 4 ⋅ x2 + ( x1 ) 2 + 2 ⋅ ( x2 ) 2 x1 := 0 x2 := 0 Given x1 + 2 ⋅ x2 ≤ 8 2 ⋅ x1 − x2 ≤ 12 x1 ≥ 0 x2 ≥ 0 15 xopt := Minimize( F , x)  1 xopt =   F ( xopt ) = −3  1 В PTC Mathcad Prime начальные приближения для решения, ограничения и функция поиска минимума записываются в соответствующие поля Блока решения в закладке Математика: 2.1. Выбор инвестиционного портфеля (задача Марковица) Рентабельность некоторого i-го актива (акция, облигация и т.п.) за определённый период Ri = Di + Pi1 − pi0 , pi0 где Di – дивиденды или текущий доход от актива (случайная величина), pi0 – цена актива в начале периода (известна), Pi1 – цена актива в конце периода (случайная величина). Т. о. Ri случайная величина и её смысл – прибыль на рубль затрат. Часто рентабельность выражают в процентах, умножая это выражение на 100%. Пусть имеются w – средства для инвестирования, на рынке имеются n видов активов. wi – средства, инвестированные в i-й актив. Тогда прибыль от портфеля: R1w1 + R2 w2 + ... + Rn wn , 16 n рентабельность портфеля Rπ = Rw i =1 i i w n =  Ri xi = x T R , xi = i =1 wi – доля средств, инвестированw  x1  x    ных в i-й актив, вектор портфеля x =  2  .  ...   xn  Г. Марковиц в 1952 г. поставил и решил задачу: найти портфель x * , обеспечивающий ожидаемую рентабельность портфеля MRπ ≥ r0 при наименьшем риске (дисперсии рентабельности портфеля DRπ → min). n n i =1 i =1 MRπ = M ( Ri xi ) =  xi MRi , ai = MRi – ожидаемая рентабельность i–го актива (можно оценить по прошлым годам, считаем известными). DRπ = D ( x T R ) = M ( x T R − Mx T R ) 2 = M ( x T ( R − M R )) ⋅ ( x T ( R − M R ))T = = M ( x T ( R − M R ) ⋅ ( R − M R )T x ) = x T M [( R − M R ) ⋅ ( R − M R )T ] x = x T K R x , где K R – ковариационная матрица случайного вектора рентабельностей активов R : DR1 cov( R1 , R2 )   cov( R2 , R1 ) DR2 K R = M [( R − M R)( R − M R )T ] =   ... ...   cov( Rn , R1 ) cov( Rn , R2 ) ⋯ cov( R1 , Rn )   ... cov( R2 , Rn )   ... ⋱  ... DRn  cov( Ri , R j ) = M ( Ri − MRi )( R j − MR j ) , (можно оценить по прошлым годам, считаем извест- ными). Т. о. получилась задача квадратичного программирования:  n  xi ai ≥ r0  i =1  n F ( x ) = x T K R x → min,  xi ≤ 1 . i = 1   x ≥ 0, i = 1, n  i  Пример 2.3 Пусть на рынке есть 3 актива с ожидаемыми рентабельностями a1 = 11%, a2 = 15%, a3 = 8%. 0.5 −0.7   1.5  Оценка ковариационной матрицы K =  0.5 2.5 −0.3  [%]2. Как правило, чем доход −0.7 −0.3 1   нее актив, тем больше дисперсия его доходности. Найти оптимальный портфель, имеющий ожидаемую рентабельность не менее 11%, если доля вложения в третий актив должна быть не более 0.25. 17 F ( x ) = 1.5 x12 + x1 x2 − 1.4 x1 x3 + 2.5 x22 − 0.6 x2 x3 + x32 → min; 11 − 11x1 − 15 x2 − 8 x3 ≤ 0 x + x + x ≤ 1  1 2 3 .   x3 − 0.25 ≤ 0  x1 , x2 , x3 ≥ 0  0.436    F ( xopt ) = 0.472[%]2 Решение этой задачи в MathCAD даст xopt =  0.28  , . Риск портфеля  0.25  σ = 0.687%   меньше, чем риск самого мало рискового актива в нём. Это эффект диверсификации – комбинирование ценных бумаг в нужных пропорциях. 3. Динамическое программирование Пример 3.1 Пусть туристам нужно добраться из пункта A в пункт D (к железной дороге) за K (= 2) дня (перехода, этапа) кратчайшим путём. 18 Обозначим x k – состояние в момент времени k (k = 0, 1, 2 = K) т. е. где они. u k – управление (принятое решение по какой дороге пойти), применённое в момент k, lk ( x k , u k ) – путь, прой- дённый при этом (на k+1-м этапе). Например, l1 ( B, ε ) = 2 . Тогда суммарный путь 1  l (x k =0 k k , т.е. надо решить, какие лучше применить управления? , u k ) → min 0 1 {u , u } Следует ли сразу взять u 0 = β (с кратчайшим путём)? Нет, это было бы недальновидно! Попробуем анализировать с конца. Так как в итоговую точку можно прийти на последнем этапе путём разных управлений, то последний участок нужно выбрать покороче, но сам по себе, а вместе с тем путём, который по минимуму придётся пройти, чтобы попасть в стартовую точку последнего этапа. Обозначим zk ( x k ) – минимальное значение целевой функции (пути) для достижения состояния x k . Тогда zk ( x k ) = min {lk −1 ( x k −1 , u k −1 ) + zk −1 ( x k −1 )} , при k = K это даст минимум (по всем u k −1 , ведущим в x k . x k −1 соответствующие) итогового пути. Но в каждую точку x k −1 можно попасть тоже путём нескольких управлений, lk − 2 ( x k −2 , u k −2 ) + zk −2 ( x k −2 )} и т. д. до z1 ( x1 ) = min{ l0 ( x 0 , u 0 )} , т.к. поэтому zk −1 ( x k −1 ) = min{ k −2 u u z0 ( x ) = 0 . Т. о.  z0 ( x 0 ) = 0,  k lk −1 ( x k −1 , u k −1 ) + zk −1 ( x k −1 )}, k = 1, 2,..., K  zk ( x ) = min{ u k −1 (10) Пример 3.1(продолжение)  x1 = C , 0min {3, 2} = 2 u ∈{α , β }  k = 1. z1 ( x1 ) = min l0 ( x 0 , u 0 ) =  1 . u x = B , min{5} = 5  u 0 ∈{γ }  Или эти же вычисления в виде таблицы x1 C u0 α β γ A A A l0 ( x , u ) 3 2 5 z1 ( x1 ) 2 x B 5 Обозначим u k −1 ( x k ) – выявленное условно-оптимальное управление в k-1-й момент (в зависимости от будущего состояния x k ). В нашем примере это подчёркнутые управления во второй строке. Итоги таблицы 1-го этапа заносим в основную таблицу: 19 k=1 xk k=2 1 C B D z1 ( x ) u (x ) z2 ( x 2 ) u1 ( x 2 ) 2 5 – β γ – – – 7 – – ε 1 Теперь второй (последний) этап: l1 ( x1 , u1 ) + z1 ( x1 )} k = 2. z2 ( x 2 ) = min{ 1 u x 2 D u1 δ ε μ x1 C B B l1 ( x1 , u1 ) 6 2 3 z1 ( x1 ) 2 5 5 Σ 8 7 7 8 2 z2 ( x ) Итог в основной таблице. Видим, что минимально пункт x 2 = D можно достичь через 7 км. Осталось обратным ходом проследить, как это удалось: u*1 = u1 ( x*2 = D ) = ε  значит x*1 = B  u*0 = u 0 ( x*1 = B ) = γ  x*0 = A . Ответ: x*k A 1 B u*k γ ε zk ( x*k ) 5 k 2 D 7 Уравнение Р. Беллмана Есть некоторая система (например, предприятие) и дискретные моменты времени (начала очередных финансовых годов k = 0, 1, …, K). 20 Состояние системы в каждый из этих моментов характеризуется некоторой величиной x k (или набором величин – вектором состояния x k ). В моменты k = 0, 1, 2, …, K-1 вы должны применять управления u k (или векторные u k ), при этом несёте затраты (потери) lk ( x k , u k ) зависящие от наличного состояния и выбранного управления, а система переходит в новое состояние x k +1 = f k ( x k , u k ) . Задача динамического программирования: найти оптимальный набор управлений {u*0 , u*1 ,..., u*K −1} , обеспечивающий нужное конечное состояние x K с минимумом суммарных по- терь z K ( x K ) = K −1 min 2 K −1 {u 0 ,u ,...,u }  l (x k =0 k k , u k ) . Решение уже найдено, это рекуррентное уравнение Р. Бел- лмана (10) (в прямой форме). По нему находится последовательность функций {u k −1 ( x k )} , дающих условно-оптимальное управление для каждого состояния x k (см. основная таблица). Затем, обратный ход: • k = K , оптимальное состояние известно: x*k = x K (заданное) ; • находим оптимальное управление u*k −1 , ведущее в известное оптимальное состояние x*k : u*k −1 = u k −1 ( x*k ) ; • находим оптимальное состояние предыдущего момента: x*k −1 = f k−−11 ( x*k , u*k −1 ) ; • k – 1 = 0? если да – Stop, известны все оптимальные управления u*K −1 , u*K − 2 ,..., u*0 , оптималь- ная траектории x*K , x*K −1 ,..., x*0 , и динамика критерия качества (целевой функции) z K ( x*K ), z K −1 ( x*K −1 ),..., z0 ( x*0 ) ; • если нет, k = k – 1 и возврат к п. 2. 3.1. Стохастические модели динамического программирования Стохастический – от греческого слова «благоразумный», т.е. модели применяются к случайным, вероятностным явлениям в системах. Здесь после вашего управления «природа» добавляет своё случайное управление и новое состояние системы становится случайным. От этого и выигрыш на k – м этапе становится, как правило, случайным Wk −1 ( x k −1 , u k −1 ) , и оптимизировать можно только математическое ожидание суммарного выигрыша за все этапы. 21 Т. е. в этом примере Wk −1 ( x k −1 , u k −1 = 1) может принять значение w1 с вероятностью p1, а система оказаться в состоянии x1k ; w2 с вероятностью p2 , а система оказаться в состоянии x2k ; w3 с вероятностью p3 , а система оказаться в состоянии x3k . Дерево решений – графическое изображение последовательности возможных состояний системы, переходов между ними под влиянием решений человека и “природы”, с указанием соответствующих вероятностей и выигрышей. Пример 3.2 Фирме предстоит принять решение: внедрять или нет в массовое производство новый продукт. В начальном состоянии у фирмы есть 3 возможности: отказаться от проекта вообще (D), приступить к массовому производству (L), произвести пробный сбыт небольшой партии на экспериментальном рынке (T). При отказе – дальнейший выигрыш = 0, при L с вероятностью 0.3 может быть «успех» (прибыль 3 млн. долл.) и с вероятностью 0.7 «неудача» (убыток 0.25 млн. долл. (прибыль 0.25)). Пробный сбыт требует затрат 22500 долл. и с вероятностями 0.5, 0.25, 0.25 приводит к состояниям: a) продукт продегустировали < 10% потребителей; b) продукт продегустировали > 10%, купили повторно < 50% из них; c) продукт продегустировали > 10%, купили повторно > 50% из них. После каждого из них можно отказаться от проекта или запустить массовое производство, вероятности «успеха» при этом 0.06, 0.28, 0.8 соответственно. (Эти вероятности берутся из накопленной статистики по аналогичным примерам или из субъективных оценок экспертов). Что делать в начальном состоянии? Как поступать в дальнейшем, чтобы ожидаемый суммарный выигрыш за все этапы был максимальным? 22 Ясно, что в таких задачах мы не станем заранее связывать себя обязательствами по нашим управлениям на все шаги вперёд, а наоборот, подождём, в каком состоянии окажемся в момент принятия решения и, исходя из этого, примем решение, т.е. желательно узнать условно оптимальные управления u k −1 ( x k −1 ) - что делать в момент (k-1) если в этот момент мы оказались в состоянии x k −1 . Такие условно оптимальные управления даёт уравнение Беллмана в обратной форме. Обозначим f k −1 ( x k −1 ) – максимальное математическое ожидание прибыли от момента (k – 1) до K, в зависимости от имеющегося состояния x k −1 . Ясно, что оно складывается из математического ожидания прибыли на ближайшем этапе и ожидания того, что будет после него до конца:  f k −1 ( x k −1 ) = max{ M [Wk −1 ( x k −1 , u k −1 ) + f k ( X k )]}, k = K , K − 1,...,1 u k −1   f K ( x K ) = 0 (11) Пример 3.2 (продолжение) MW1 ( x1 , u1 )} k = 2 = K. f1 ( x1 ) = max{ 1 u 1 x O S F a b u1 – – – D L x2 – – – O S MW1 ( x1 , u1 ) f1 ( x1 ) 0.06*3+0.94*(- 0 0.66 0.25) = - 0.055 0.66 F c D L O S F D L O S 2.35 2.35 M [W0 ( x 0 , u 0 ) + f1 ( X 1 )]} k = 1. f 0 ( x 0 ) = max{ u x x0 u0 D L x1 O S T F a b 23 c F MW0 ( x 0 , u 0 ) 0.725 – 0.0225 Mf1 ( X 1 ) 0.5f1(a)+0.25f1(b)+0.25f1(c) = 0.7525 Σ 0.725 0.73 f0 ( x ) 0.73 Т. о. u*0 ( x 0 ) = T (пробный сбыт), а мы и есть в x 0 , значит это безусловно оптимальное управление. Ожидаемая прибыль за все этапы f 0* = f 0 ( x 0 ) = 0.73 млн. долл. Из приведённых таблиц следует, если x*1 = a → u*1 = u1 (a) = D; x*1 = b → u*1 = u1 (b) = L; x*1 = c → u*1 = u1 (c) = L. k x0 1 a x*k b c u*k T D L L f k ( x*k ) 0.73 0.66 2.35 Ожидаемое значение прибыли с момента 1 может быть больше, чем с момента 0 (2.35> 0.73). Задача о наилучшем выборе Пример 3.3 Невеста выбирает себе жениха из N претендентов, среди которых, ясно, один наилучший. Претенденты появляются по одному и отвергнутого вернуть уже нельзя. Процесс выбора заканчивается, как только невеста останавливается на каком-либо очередном претенденте. Как ей поступить? Ограничимся процедурами отбора следующего типа: первый жених либо принимается, либо отвергается. Если он отвергается, то автоматически отвергаются и все последующие женихи, которые хуже, чем первый, и следующее решение принимается лишь при появлении жениха лучшего, чем первый. Он (второй лидер) либо принимается, либо отвергается. Если он отвергается, то автоматически отвергаются и все последующие женихи, которые хуже, чем второй лидер, и следующее решение принимается лишь при появлении жениха лучшего, чем второй лидер. И так далее. Ясно, что если взять раннего лидера, то вероятно, что самый лучший ещё впереди, если многих лидеров пропустить, то вероятно, что остальные все хуже и выбор вообще не состоится. Обозначим k = 0, 1, 2, …– моменты появления очередных лидеров (очевидно, всё лучших и лучших). Состояние системы в каждый из этих моментов описывается x k – номер жениха 24 (число уже появившихся), ставшего лидером в момент k. Т. о. x 0 = 1 и т.д. В каждый момент k мы можем применить управление u k (= 1 – принять соответствующего лидера за жениха, = 2 – продолжить осмотр следующих). Цель – максимизировать мат. ожидание вероятности того, что выбранный в момент k* жених является наилучшим из всех. Если выбор останавливается на x k -м по счёту женихе, то вероятность, что он и есть наилучший: w( x k , u k = 1) = xk . Т.к это вероятность, что один конкретный (наилучший воN обще) из N попадёт в первые x k человек. k −1 Каковы вероятности переходов Px(ku−1 , x=k 2) ? k −1 Px(ku−1 , x=k 2) = P( x k | x k −1 ) = P(следующий лидер будет x k - м|предыдущий был x k −1 - м) = P[( x k −1 был лидером, в диапазоне x k −1 + 1, x k −1 + 2,..., x k − 1 лидера не было) ⋅ ⋅( x k - м будет лучший из x k )] = x k −1 1 ⋅ . xk − 1 xk Первый множитель: лидер из x k − 1 людей вошёл в x k −1 число. Т. к. они зависят от момента (k – 1), цепь не однородна. Пусть f k −1 ( x k −1 ) – максимальное математическое ожидание вероятности выбора наилучшего на шагах с k – 1 до последнего в зависимости от состояния x k −1 . Уравнение Беллмана f k −1 ( x k −1 ) = max{ w( x k −1 , u k −1 = 1), k −1 u N  x k = x k −1 +1 k −1 f k ( x k ) Px(ku−1 , x=k 2) } , второе выражение соответствует управлению «продолжить». Слагаемое с x k = ∅ равно 0 т.к. f k (∅) = 0 . 25 Если x k −1 = N − 1 , то сумма состоит из одного слагаемого, в нём, ясно f k ( x k = N ) = 1 , т.к. это означает, что последний лидер – абсолютный и вероятность успеха = 1! Но вероятность перехода в такое «счастье» очень мала. Сравним: N  1⋅ Pxk−1xk = x =N k x k −1 x k −1 < w( x k −1 , u k −1 = 1) = , ( N − 1) ⋅ N N значит u k −1 ( x k −1 = N − 1) = 1 – «брать», и x k −1 . N Так будет, пока второе выражение не превысит первое, тогда впервые (с конца) условно оптимальным управлением будет «продолжить», т. е. k-1=k* f k −1 ( x k −1 ) = N N xk* x k *+1 xk* x k* 1 x k *+1 ≤  Px k* , xk*+1 ≤  или , N x k*+1 = x k* +1 N N x k*+1 = xk* +1 x k *+1 − 1 x k *+1 N N 1  l − 1 ≥ 1 . Или l = m +1 1 1 1 + + ... + ≥ 1 . Отсюда находим максимальное m, при котором m m +1 N −1 это выполняется: mmax = x k * – номер человека-лидера, при котором u k * = «продолжить» (при меньших x k u k также будет «продолжить»). Т.о. нужно последним пропустить лидера с номером mmax (если он, очередной лидер, появится с этим номером). Следующего лидера уже надо брать. Это эквивалентно такой стратегии: пропустить первые mmax женихов, затем взять первого, который лучше всех предшествующих. f * = N=2 N=3 1 N N mmax l = mmax +1 l − 1  1 mmax = 1 ≥1 1 1 1 mmax = 1 + ≥1 1 2 N=4 1 1 1 + + ≥1 1 2 3 N=5 1 1 1 + + ≥ 1 mmax = 2 2 3 4 … … N = 10 mmax = 1 mmax = 3 f* = 0.5 1 1 1 ( + ) = 0.5 3 1 2 f* = 11/24=0.458 f* = … … … … N = 20 … mmax = 7 … … … N=∞ … N/e … 1/e 26 27 4. Принятие решений в условиях неопределённости и риска Пример 4.1 Вы собираетесь пойти на работу и думаете: надеть плащ или нет. Погода – неконтролируемый, внешний фактор может быть: θ = 1 – сухо, θ = 2 – дождь, θ = 3 – ливень. Фактор, которым вы распоряжаетесь (контролируемый), может быть: a = 1 – выйти без плаща, a = 2 – надеть плащ. M0 = {1, 2} – множество (чистых) стратегий. Всевозможные сочетания этих факторов ведут к следующим вашим потерям или неудобствам: L(a, θ) – функция или матрица потерь (в условных денежных единицах). 0 5 7 L ( a, θ ) =   , т. е. наибольшие потери 7 будут, если вы без плаща попали в ливень. 3 2 1 Иногда удобнее говорить не о потерях, а о выигрышах  0 −5 −7  W ( a, θ ) = − L ( a, θ ) =   , и тогда оперирующая сторона стремится максимизиро −3 −2 −1  вать выигрыш. 4.1. Условия неопределённости Никакой дополнительной информации о неконтролируемом факторе нет. Назовём оценкой эффективности стратегии a величину W (a ) = min W (a,θ ) – то есть θ наименьший выигрыш, который может быть при данной стратегии a (гарантированный результат). W (1) = min W (1, θ ) = min(0, −5, −7) = −7; θ θ W (2) = min(−3, −2, −1) = −3. θ Стратегия a* называется оптимальной, если она приносит наилучший гарантированный результат: a* = arg max W (a) . W ( M 0 ) = W (a*) – наилучший гарантированный результат во a∈M 0 множестве стратегий M0. W (1)   −7  a* = arg max  = arg max   = 2 , т.е. нужно надеть плащ. W ( M 0 ) = W (2) = −3 наилуч a a W (2)   −3  шее, что можно гарантировать правильным выбором стратегии. Это стандартный (весьма пессимистичный) взгляд (образ действий) называют максиминным подходом [поскольку a* = arg max min W (a,θ ) ] или критерием Вальда (A. Wald). a θ Некоторые нестандартные критерии 28 А. Критерий оптимизма (макси-макса). В этом случае каждая стратегия характеризуется W (a) = max W (a,θ ) наибольшим возможным выигрышем θ и оптимальной считается aO = arg max W (a) . a Если применить этот критерий в примере 4.1: W (1) = max W (1,θ ) = max(0, −5, −7) = 0; θ θ W (2)= max(−3, −2, −1) = −1; θ 0 aO = arg max   = 1 – без плаща. a  −1 Б. Критерий пессимизма-оптимизма Гурвица: { } aH = arg max α min W (a,θ ) + (1 − α ) max W (a, θ ) , где α – доля пессимизма, 0 ≤ α ≤ 1 . a θ θ Продолжение примера: пусть α = 0.5   −7   0   −3.5  aH = arg max 0.5   + 0.5    = arg max   = 2 – надеть плащ. a a  −1   −2    −3  4.2. Условия риска (критерий Байеса – Лапласа) Критерий Вальда (с оглядкой на худшее состояние природы) выглядит уж очень осторожным. Если у вас есть представление о том, как редко это состояние природы бывает, и что чаще бывают не столь губительные состояния, то может быть стоит оценивать стратегии по среднему выигрышу, который она дает при многократном применении. Другими словами, внешний фактор θ считается случайной величиной Θ, имеющей известное распределение pΘ (θ ) – плотность распределения или закон, и оценкой эффективности стратегии a называют среднее по состояниям природы θ значение выигрыша при данной стратегии (математическое ожидание)   W ( a, θ ) pθ  , EMV ( a ) = MW (a, Θ) =  θ   W ( a, θ ) pΘ (θ ) dθ  в дискретном и непрерывном случае соответственно. EMV – expected monetary value. И тогда, естественно, aBL = arg max EMV (a) . (11a) a EMV (aBL ) - наилучший ожидаемый денежный результат. Пример 4.1 (продолжение) Считаем известным распределение погоды θ pθ*100% 1 60 2 30 29 3 10 EMV (1) = W (1,1) p1 + W (1, 2) p2 + W (1,3) p3 = 0 ⋅ 0.6 + (−5) ⋅ 0.3 + (−7) ⋅ 0.1 = −2.2; EMV (2) = −3 ⋅ 0.6 + (−2) ⋅ 0.3 + (−1) ⋅ 0.1 = −2.5;  −2.2  aBL = arg max   = 1 – без плаща, наилучший ожидаемый денежный результат EMV (aBL ) = a  −2.5  -2.2. Сравним это с максиминным решением a* = 2 – надеть плащ. При том решении средний выигрыш меньше (-2.5), но хуже, чем наилучший гарантированный результат (-3) не будет. При a* = 1 средний выигрыш лучше, но можно и сильно промокнуть (-7). Ожидаемая ценность точной информации (EVPI) Пример 4.1 (продолжение) Пусть предлагается услуга – утром перед выходом из дома вам сообщат, какая на улице погода (точная информация). Сколько (максимально) следует заплатить за услугу накануне вечером? Если проинформируют, что сухо, то ясно – пойдёте без плаща и w = 0, если дождь, то в плаще и w = -2, если ливень, то в плаще и w = -1. А с какими вероятностями так сообщат – с теми, как оно бывает. То есть ожидаемый оптимальный выигрыш при покупке информации 0 ⋅ 0.6 + ( −2) ⋅ 0.3 + ( −1) ⋅ 0.1 = −0.7 . Разность: ожидаемый оптимальный выигрыш (при покупке информации) – оптимальный ожидаемый выигрыш (в условиях риска) и составляет ожидаемую ценность информации (EVPI – expected value of perfect information) EVPI = −0.7 − ( −2.2) = 1.5 . Вплоть до такой цены можно заплатить за информацию. Очевидно, EVPI = M max W (a, Θ) − max MW (a, Θ) . a a Пример 4.2 Фирма продаёт скоропортящийся продукт. Нереализованный в данном периоде продукт не может быть продан в следующем периоде. k1 – себестоимость и затраты на хранение пропавшего товара – издержки излишка (у.е. на 1 тонну). k2 – потеря прибыли из-за неудовлетворённого спроса – издержки недостатка (у.е. на 1 тонну). Если a – запас и θ - предстоящий спрос то, очевидно, функция потерь k1 (a − θ ), при θ ≤ a, L ( a, θ ) =  k2 (θ − a), при θ > a. Пусть предстоящий спрос Θ случаен, с распределением p(θ ) . Найти оптимальный запас a (в тоннах) по критерию Байеса – Лапласа. Подставим функцию потерь в выражение для EMV(a). 30 a ∞ a EMV (a ) = k1  (a − θ ) p (θ )dθ + k2  (θ − a ) p (θ )dθ = a a = k1[aF (a ) −  θ p (θ )dθ ] + k2 {M(Θ) −  θ p (θ )dθ − − a[1 − F (a )]} = a = (k1 + k2 )aF (a ) + k2 M(Θ) − k2 a − (k1 + k2 )  θ p (θ )dθ . Здесь F (θ ) = P(Θ ≤ θ ) - функция распределения вероятностей спроса. Минимизируем по a: ∂EMV (a ) k2 . Т. о. опти= (k1 + k2 ) F (a ) − k2 = 0 , откуда F ( aB ) = ∂a k1 + k2 мальный запас aBL есть квантиль порядка k2 распределения спроса k1 + k2 aB = q ( k2 ). k1 + k2 Напомню, квантиль порядка p распределения есть то минимальное значение аргумента функции распределения, при котором она ≥ p. Учитывая монотонно возрастающий характер квантильной функции, найденный результат конкретизирует очевидные качественные соображения относительно влияния соотношения издержек k1 и k2 на выбор запаса. 4.3. Статистические решающие правила Резюмируем суть «игр против природы»: оперирующая сторона или лицо, принимающее решение (ЛПР), не обладает полной информацией о состоянии природы (конъюнктура рынка, замыслы конкурентов, погодные условия, состояние оборудования и т. п.) – о стратегии природы θ. Однако ЛПР должен предпринять какое-то действие (стратегию-константу) a. Известна функция потерь L ( a , θ ) , которые несёт оперирующая сторона, предприняв действие a в условиях природы θ. Новым в «статистических играх» будет возможность оперирующей стороне перед принятием решения провести эксперимент, в котором она получает значение x случайной величины X, которая связана с θ, причём известно её условное распределение p X ( x, θ ) ≜ p ( x | θ ) . Функция a = d ( x ) , определяющая действия ЛПР при всех возможных исходах эксперимента x, называется решающей функцией (решающим правилом, стратегией-функцией). Принятие решения по правилу d, даже при не случайном θ, в силу случайности X ведёт к случайным потерям L ( d ( X ), θ ) . Ожидаемое их значение называется риском. Т. о. риск — это ожидаемые потери, которые вы понесёте, приняв решение (а его принимать всё равно придётся) в соответствии с таким-то решающим правилом (d) в таких-то внешних условиях (θ). R(d , θ ) = M {L[d ( X ), θ ] | θ } =  L[d ( x), θ ]p( x | θ )dx . 31 (12) Для всех сочетаний правил и условий - функция риска. Т. о. при использовании стратегий-функций акцент переводится с функции потерь на функцию риска, и поиск оптимальной функции d будет связан с минимизацией риска. Минимаксное решающее правило Применим стандартный подход – критерий Вальда, который состоит в следующем: т. к. θ неизвестен, то любое правило d ( x ) лучше охарактеризовать максимально возможным при нём риском – rm (d ) = max R(d ,θ ) (максимальный риск правила d). Следует остановиться на праθ виле, где этот риск минимален. Это определяет минимаксное решающее правило: d mm = arg min max R(d , θ ) . (13) θ d Часто его можно найти перебором всех возможных d по минимуму rm ( d ) . rm * = rm (d mm ) – минимаксный риск. Пример 4.3 Обычная цена месторождения 2 у. е. На торги выставлено 3 месторождения (пакет) за 4 у. е. Эта дешевизна неспроста: из них одно (состояние «природы» θ1 ) или два (состояние «природы» θ 2 ) исчерпаны. У вас есть 2 варианта действий: a1 – покупать, a2 – не покупать. Функθ1 θ2 a 0 2 ция потерь L ( a, θ ) = 1  . Ваши потери в 1 у. е. будут, если неплохие (при θ1 ) месторожa2  1 0  дения достанутся конкурентам. Продавец согласен допустить вас проверить одно любое месторождение. Результатом мо- жет быть: x1 – месторождение годное или x2 – месторождение исчерпанное. Распределение результата X в зависимости от состояния природы заранее известно (и очевидно таково): p( x | θ ) = θ1 θ2 x1   x2  13 1 3 2 3 2 3  .  Возможны всего 4 различных решающих правил: a1 , x = x1 , a1 , x = x1 , a2 , x = x1 , a2 , x = x1 , d1 ( x) =  d 2 ( x) =  d3 ( x ) =  d 4 ( x) =  a1 , x = x2 ; a2 , x = x2 ; a1 , x = x2 ; a2 , x = x2 . Можно обозначить их в виде цепочек d1 = 11 ; d 2 = 12 ; d3 = 21 ; d4 = 22 . На первом месте указано, что делать при первом возможном исходе эксперимента, на втором – что при втором исходе и т. д. В общем случае если число возможных действий – Na и число возможных исходов эксперимента – Nx, то имеется ( Na ) Nx различных решающих правил-цепочек длины Nx (для одного исхода – Na, для двух – Na⋅Na вариантов и т. д.). Найдём значения функции риска правила d 2 при всех θ. По формуле (12) 32 R(d 2 , θ ) =  L(d 2 ( x), θ ) p( x | θ ) = L(d 2 ( x1 ), θ ) p( x1 | θ ) + L(d 2 ( x2 ), θ ) p ( x2 | θ ) = x = L(a1 , θ ) p( x1 | θ ) + L(a2 , θ ) p( x2 | θ ) = =| сразу для двух θ |= (0, 2) ⋅ ( 23 , 13 ) + (1, 0) ⋅ ( 13 , 23 ) = (0, 23 ) + ( 13 , 0) = ( 13 , 32 ). Поэтому максимальный риск этого правила rm (d 2 ) = max( 13 , 32 ) = 2 / 3 . Точно так же нахоθ дим rm (d1 ) = 2 , rm (d3 ) = 4 / 3 , rm (d 4 ) = 1 . По (13) минимаксное решающее правило d mm 2    2 / 3  = arg min = d2 . d  4 / 3   1  Минимаксный риск rm * = rm (d mm ) = rm (d 2 ) = 2 / 3 . В отличие от стратегических игр, минимаксный риск 2/3 в нашей статистической игре не является гарантированным результатом при применении минимаксной стратегии d 2 . Например, получив в эксперименте x1, будем покупать пакет месторождений, но при этом может быть θ = θ 2 и потери будут 2 у. е. Байесовское решающее правило Будем предполагать, что состояние природы Θ – случайная величина с известным априорным распределением p (θ ) , тогда любое правило d ( x ) можно охарактеризовать средним по состояниям природы риском (байесовский риск правила d), аналог EMV(a): rB (d ) = MR (d , Θ) =  R(d , θ ) p(θ )dθ (14) (ранее каждое d характеризовали худшим (максимальным) по состояниям природы риском). Следует остановиться на правиле, где этот риск минимален, – байесовское решающее пра- вило: d B = arg min rB (d ) . d rB * = rB (d B ) – минимальный байесовский риск. Подчеркнём, что байесовское решающее правило зависит от предположенного априорного распределения p (θ ) . Пример 4.4 (продолжение примера 4.3) Будем считать, что p(θ1 ) = α , p(θ 2 ) = 1 − α составляет априорное распределение. Функцию риска уже рассчитывали. Найдём байесовские риски всех решающих правил. rB (d1 ) =  R (d1 , θ ) p (θ ) = R (d1 , θ1 )α + R (d1 , θ 2 )(1 − α ) = 0α + 2 ⋅ (1 − α ) = 2 ⋅ (1 − α ); θ rB (d 2 ) = ... = 13 α + 32 ⋅ (1 − α ) = 13 ⋅ (2 − α ); rB (d3 ) = ... = 23 α + 43 ⋅ (1 − α ) = 23 ⋅ (2 − α ); rB (d 4 ) = ... = 1α + 0 ⋅ (1 − α ) = α . 33 Пусть α = 0.4 (т.е. вероятность «хорошего» пакета меньше половины),  1.2   1.2      0.533  0.533    rB (d ) = , d B = arg min rB (d ) = arg min = d 4 , минимальный байесовский риск d d  1.06   1.06       0.4   0.4  rB * = rB (d B ) = 0.4 у. е. Применённый подход легко реализуем, когда стратегий-функций не много и они явно известны. Однако хотелось бы уметь автоматически строить d B как функцию от x, т. е. для каждого возможного значения x заготовить a = d B ( x) . Теорема (байесовский принцип) Байесовская решающая функция имеет вид: a = d B ( x) = arg min EMV (a | x) , a (15) где EMV ( a | x ) – ожидаемая денежная ценность (потери) действия a при известном x, (т. е. следует делать то, что несёт наименьшие ожидаемые при данном x потери). EMV (a | x) = M [ L(a, Θ) | x] =  L(a,θ ) p(θ | x)dθ , (16) Здесь p (θ | x ) – апостериорное распределение состояния природы, p(θ | x) = p(θ ) p( x | θ )  p(θ ) p( x | θ )dθ – формула Байеса. Доказательство. Подставляем в байесовский риск (14) функцию риска (12), меняем местами интегрирование, и, т. к. по формуле Байеса p (θ ) p ( x | θ ) = p ( x ) p (θ | x ) , имеем, rB (d ) =  { L[d ( x), θ ]p (θ | x)dθ } p( x)dx . (17) Минимизация rB ( d ) достигается минимизацией внутреннего интеграла при каждом x, значит нужно брать действие d ( x ) = a дающее интеграл поменьше, что и означает (15). Пример 4.5 (продолжение примера 4.4) Построим байесовскую решающую функцию по формуле (15). Вначале рассчитаем ожидаемую денежную ценность каждого действия при каждом x: p (θ ) p( x | θ ) EMV (a | x) =  L(a,θ ) . Пусть x = x2. p ( x) θ EMV (a | x2 ) = = 1 [ L(a, θ1 )α p ( x2 | θ1 ) + L(a, θ 2 )(1 − α ) p( x2 | θ 2 ) = p( x2 ) 1 0... + 2 ⋅ (1 − α ) ⋅ 23  1  43 ⋅ (1 − α )  =   1 . α p ( x2 ) 1α ⋅ 13 + 0... p ( x ) ⋅ 2  3  34 Пусть α = 0.4, EMV (a | x2 ) = d B ( x2 ) = arg min EMV (a | x2 ) = arg min a a 1 0.4  = a2 . p ( x1 ) 0.266  d B ( x1 ) = arg min a 1 0.8  , байесовское действие при этом x: p( x2 ) 0.133 1 0.8  = a2 . Аналогично p ( x2 ) 0.133 Т. е. байесовская решающая функция a2 , x = x1 d B ( x) =  = d 4 . Минимальный байесовский риск по формуле (17) a2 , x = x2 rB * = rB (d B ) =  EMV (d B ( x) | x) p ( x) =  (min EMV (a | x)) p ( x) =  min( EMV (a | x) p ( x)) (18) x x a x a = 0.133 + 0.266 = 0.4 у.е. Результаты совпадают с прежним решением (путём перебора байесовских рисков всех решающих правил). Свойства минимаксных и байесовских стратегий. Байесовское решение, построенное для любого априорного распределения p (θ ) , является допустимым, т. е. нет другого, лучшего правила с риском меньшим при всех θ. Если существует априорное распределение p (θ ) , для которого байесовское правило dB имеет постоянную функцию риска R(d B ,θ ) = const, то такое dB есть минимаксное правило. Доказательства см. [14]. Резюмируем ситуации информированности о внешнем факторе: 1. Если знаем θ, выбираем решение a = arg min L(a,θ ) . a 2. Если не знаем θ. a) Нет эксперимента. Критерии Вальда, Байеса – Лапласа. b) Есть эксперимент. По косвенной информации о θ, содержащейся в результатах эксперимента x, выбираем решение a = d ( x) . Но по какой d? Даже зная свойства эксперимента (связь x с θ), единой, лучшей всегда, решающей функции нет, т. к. риск решающей функции зависит от значения θ. Т. е. опять всё упирается в θ! Но при одинаковом незнании θ, по сравнению со стратегическими играми (Вальд, Байес – Лаплас) проведение эксперимента (и затем d mm ( x) , d B ( x) ) должно давать преимущество. Если оптимальная решающая функция уже найдена, то можно рассчитать ожидаемую ценность информации от эксперимента EVEI. EVEI = мера (зависящая от информированности о θ) потерь оптимального решения без эксперимента – та же мера риска оптимальной решающей функции. Больше этого тратить на эксперимент нерационально. Пример 4.6 (продолжение примера 4.3) • Нет информации о распределении состояний природы. 35  2 Без эксперимента: a* = arg min L (a) = arg min   = 2 т.е. оптимальное решение - не покупать. a∈M 0 a∈M 0 1   Мера потерь здесь L ( M 0 ) = L ( a∗ ) = 1 . С экспериментом: минимаксный риск оптимальной решающей функции ( d 2 = 12 ) rm * = rm (d mm ) = rm (d 2 ) = 2 / 3 . Т.о. EVEI = 1 - 2/3 = 1/3 ≈ 0.333 у.е. • Известно априорное распределение состояний природы: p(θ1 ) = α , p(θ 2 ) = 1 − α , α = 0.4. Без эксперимента:  0 ∗ 0.4 + 2 ∗ 0.6   1.2  aBL = arg min EMV (a) = arg min  = arg min   = 2 т.е. оптимальное решение  a∈M 0 a∈M 0 1 ∗ 0.4 + 0 ∗ 0.6 a∈M 0 0.4     не покупать. Мера потерь здесь EMV (aBL ) = 0.4 . С экспериментом: минимальный байесовский риск оптимальной решающей функции ( d B = d 4 = 22 ) rB * = rB (d B ) = 0.4 . Т.о. EVEI = 0.4 – 0.4 = 0 у.е.! Понятно, т.к. этот эксперимент всегда посоветует то, что и без него бы сделали (не покупать). Ясно, что крайний случай эксперимента – точная информация и EVPI = EMV (aBL ) - [ 0  2 0.4 ∗ min   + 0.6 ∗ min   ]= 0.4 – 0 = 0.4 у.е. 1 0 Пример 4.7 На технологическую линию может поступать сырье с малым θ1 и большим θ2 количеством примесей. Известно, что в среднем поступает 60 % сырья первого вида и 40 % сырья второго вида. Для использования различных видов сырья предусмотрено три режима работы технологической линии: a1, a2, a3. Потери, связанные с качеством выпускаемой продукции и расходом сырья, в зависимости от качества сырья и режима работы технологической линии составляют: a1 a2 a3 θ1 1 3 θ2 5 3 2 Допускается произвести грубый предварительный анализ содержания примесей (точный лабораторный анализ невозможен, т. к. требует затраты значительного времени, а значит, простоя оборудования). Результаты эксперимента могут быть следующими: x1 – примесей не обнаружено; x2 – примеси в небольшом количестве; x3 – примесей много. Условное распределение результатов известно: x1 x2 θ1 0.6 0.25 θ2 0.2 0.3 36 x3 0.15 0.5 Найдём байесовское решающее правило. По формуле (15) при заданном априорном распределении (0.6, 0.4) получим байесовскую стратегию d B = 123 . Минимальный байесовский риск равен 1.58 у. е. (Самостоятельный разбор решения). Если бы не знали априорное распределение ((0.6, 0.4)), то искали бы минимаксное решающее правило. Всего имеется 27 решающих правил: d1 = 111; d 2 = 112; ...; d 27 = 333 (в каждом указаны действия при исходах эксперимента x1, x2, x3 соответственно). Максимальные риски всех правил: 5 4 3.5 4.4 3.4 2.9 4.1 3.1 2.6 4.6 3.6 3.1 4 3 2.5 3.7 2.7 2.2 4.4 3.4 2.9 3.8 2.8 2.5 3.5 2.7 3 Поэтому минимаксным оказывается решающее правило d mm = d18 = 233. Минимаксный риск равен 2.2 у.е. Пример 4.8 (Продолжение примера 4.2) Пусть теперь о предстоящем случайном спросе Θ косвенно говорит x - инфляция в предыдущем периоде (в %). Найти оптимальный запас a (в тоннах) по байесовскому решающему правилу. По байесовскому принципу, оптимальный запас по байесовскому решающему правилу при наличии результата эксперимента x, находится по (15), т.е. также, как стратегия-константа по критерию Байеса – Лапласа (11a), только с использованием апостериорного распределения. Поэтому, преобразования примера 4.2 справедливы с заменой априорного распределения на k2 апостериорное. Т. о. оптимальный запас aB есть квантиль порядка апостериорного k1 + k2 распределения спроса aB = q ( k2 | x) . k1 + k2 Пусть известно, что x = 1.38. По прошлым данным, при таком же уровне инфляции в предыдущем периоде, наблюдались следующие спросы с частотами: Спрос (в тоннах) Частота <= 0 0–4 0.029 4–8 0.095 8 – 12 0.057 12 – 16 0.15 16 – 20 0.205 20 – 24 0.069 37 24 – 28 0.138 28 – 32 0.207 32 – 36 0.05 > 36 По этим группированным данным строим эмпирическую функцию апостериорного распределения вероятностей спроса Fˆ (θ | x) = Pˆ (Θ ≤ θ | x ) . Считаем распределение спроса равномерным в пределах интервалов группировки, т.е. эмпирическая функция распределения линейно нарастает в интервалах группировки от накопленной частоты быть <= начального значения интервала до накопленной частоты быть <= конечного значения интервала. См. рисунок. Например, в интервале от 20 до 24 функция изменяется от 0+0.029+0.095+0.057+0.15+0.205 = 0.536 до 0.536+0.069 = 0.605. Пусть k1 = 0.81, k2 = 1.15, k2 = 0.587. По графику видно, что квантиль порядка 0.587 k1 + k2 примерно равна 23 тоннам. Чтобы найти её точное значение, составим уравнение нужного куска прямой Fˆ = 0.0173θ + 0.191 , обратная функция θ = 57.97 Fˆ − 11.072 по ней aB = 57.97 *0.587 − 11.072 = 22.94 тонн. Выборочный контроль качества продукции Поставщик должен поставить потребителю партию изделий объёмом M = 100 шт. По договору поставщик обязуется заменять дефектные изделия, выявленные в процессе эксплуатации партии за свой счёт. Издержки замены одного дефектного изделия c2 = 2000 р. Перед 38 отправкой поставщик проводит выборочную проверку партии. Объем выборки N = 30 шт. Издержки проверки одного изделия c1 = 180 р. Если число обнаруженных дефектных изделий в выборке X ≤ m , где m – приёмное число, то партия принимается (решение a = 1). Выявленные в выборке дефектные изделия заменяются на годные, это не даёт дополнительных издержек. Если X > m, партия отвергается (решение a = 2), все изделия в ней проверяются, годные остаются в резерве (например, для замены дефектных в выборке в принятых партиях). Требуется найти оптимальное приёмное число m, минимизирующее риск (средние потери (издержки) в случае применения этого решающего правила). Обозначим θ – доля дефектных изделий в партии (характеризует состояние «природы»). 0 ≤ θ ≤ 1 . Функция потерь c1 N + c2 ( M − N )θ , при a = 1, L( a, θ ) =  c1 N + c1 ( M − N ), при a = 2. Мыслимые решающие правила: 1, x ≤ m, a = d m ( x) =  m = 0,1,..., N . Каждое из них нумеруется применяемым в нём приём2, x > m, ным числом, всего (N + 1) правил. N Функция риска правила m: R (m, θ ) =  L[d m ( x), θ ] p ( x | θ ) ; условное распределение числа x=0 дефектных изделий в выборке p ( x | θ ) для большой партии можно принять биномиальным: p ( x | θ ) = C Nx θ x (1 − θ ) N − x . Применим байесовский подход для нахождения оптимального решающего правила. Считаем известным априорное распределение θ в партиях p (θ ) . При отлаженном производстве доля дефектных изделий обычно мала, поэтому можно для него принять бета-распределение с параметрами формы 1 и 5: 39 1 Перебором: байесовские риски (14) всех правил: rB ( m ) =  R ( m, θ ) p (θ ) dθ , и байесовское ре0 шающее правило (оптимальное приёмное число mB): mB = arg min rB (m) .Расчёт в MathCAD по m этим формулам дал оптимальное приёмное число mB = 2. Средние потери (издержки) в случае применения решающего правила с приёмным числом 2 составят 16040 р. на партию (минимальный байесовский риск). Те же результаты можно найти построением байесовского решающего правила по (15), (18). 40 5. Теория систем Система – совокупность элементов, объединённых некоторой формой регулярного взаимодействия или взаимозависимости для выполнения заданной функции. Отдельные элементы системы или её подсистемы называются компонентами. Системный анализ [systems analysis] — 1. Научная дисциплина, разрабатывающая общие принципы исследования сложных объектов с учётом их системного характера. 2. Методология исследования объектов посредством представления их в качестве систем и анализа этих систем. Пример 5.1 Система управления запасами. Предприятие покупает, хранит и отпускает потребителю некоторый товар. Впереди T периодов времени, их начала в моменты t = 0, 1, …, T-1. Обозначим x t – запас в момент t, u t – объём закупки в момент t (или сразу после), d t – известный (или случайный) объём отпуска товара в момент t (или сразу после). Ясно, x t +1 = x t + u t − d t . Хранение запаса x t +1 стоит y t +1 = g t +1 ( xt +1 ) – известные функции. Т. о. процесс u t закупок определяет процесс y t издержек на хранение. Сами закупки стоят p t (u t ) . Желательно оптимизировать функционирование системы, т.е. выбрать процесс u t так, чтобы обеспечить необходимые объёмы отпуска с наименьшими суммарными затратами T −1 (y t +1 + p t (u t )) . t =0 5.1. Система как математический объект Рассмотрим подробнее основные понятия, входящие в определение системы. Прежде всего, это понятие взаимодействия, которое означает зависимость значений, принимаемых одним процессом (процесс-следствие) в текущий момент времени t, от того, какие значения принимал другой процесс (причина). Множество моментов времени будем обозначать T. В теории систем причинный процесс называют входом (input), а процесс-следствие называют выходом (output). Множество значений входа будем обозначать U, а значений выхода – через Y, элементы этих множеств – через u и y, соответственно. Зависимость выхода от значений входа в предшествующие моменты не может быть прямой (причина, погасшая до возникновения следствия, не может его вызвать…), это значит, что имеются объекты, связывающее всю предысторию входов-причин до момента t и выход в этот момент. В теории систем эти объекты называются состояниями. Множество состояний будем обозначать через X, а его элементы как x. Итак, в каждый момент t система характеризуется 41 некоторым состоянием. Это состояние однозначно определяет значение выхода в этот момент t. В состоянии накапливаются все причины, реализовавшиеся в прошлом, которые определяют событие в настоящем. Это аналогично памяти. В понятиях математической модели это означает, что y(t) является определённой функцией от x(t) y (t ) = g ( x (t )) . В общем случае характер этой зависимости может меняться во времени, т.е. g также является функцией от t, и тогда y (t ) = g (t , x (t )) (1) g называют оператором выхода. Наконец, для определения понятия системы требуется понятие, отражающее зависимость состояния от входа и предыдущих состояний. Будем считать, что состояние в будущий момент t однозначно определяется значениями, принимаемыми входом, начиная с начального момента t 0 и до момента t, при условии, что в момент t 0 было определённое состояние x(t0 ) . Пусть u(.) является отображением из T в U. Тогда x(t ) = σ (t; t0 , x(t0 ), u (.)), t > t0 . (2) σ называют оператором динамики состояния. Система определена, если заданы множества T, U, Y, X, оператор выхода g, оператор динамики σ и уравнения (1) и (2). Они представляют в математической форме закон развития (поведения, движения) системы. При известном начальном состоянии x(t0 ) и входе u (.) эти уравнения позволяют найти выход y(t) для любого t > t0 по формуле y (t ) = g (t , σ (t; t0 , x(t0 ), u (.))) . В теории систем помимо математического также используется схематическое представление системы, как показано на рисунке. 42 В этом представлении стрелки указывают входы и выходы, а прямоугольник – их причинно-следственную связь, определяемую законом преобразования входа в выход. Можно рассмотреть две системы, которые взаимодействуют между собой. Это означает, что какие-то входы одной из систем одновременно являются выходами другой. Соединение или композиция этих систем образуют новую систему. Такое соединение может быть последовательным или соединением с обратной связью, когда каждая система имеет в качестве входа выход другой системы. Беря другие системы Σ1 ,..., Σ m и отождествляя входы одних с выходами других, будем получать всё более разнообразные и сложные системы. В частности, так получают системы с иерархическим строением причинно-следственных связей. Такая композиция опять будет системой в соответствии с введённым определением: у неё будут свои входы, свои выходы и свои операторы выхода и динамики состояний. 5.2. Математическая классификация систем Математическое строение входящих в определение системы множеств и отображений лежит в основе математической классификации систем. По отношению к структуре множества T системы разбиваются на два класса. Если T совпадает с множеством действительных чисел R, то уравнения (1) и (2) определяют класс непрерывных по времени систем. Если же T = {t : t = −∞ ,..., −2, −1, 0,1, 2,... + ∞} , то система принадлежит классу дискретных по времени систем. Для дискретных систем (используя обозначения xt ≡ x(t ) и т. д.) если в (2) за t 0 и t взять соседние моменты дискретного времени t и t+1 между которыми только одно управление u t , то (1) и (2) можно переписать xt +1 = f t ( xt , u t ) , y t = g t ( xt ) . (3) Системы можно разбить ещё на два класса – детерминированные и стохастические. Модель детерминированных систем не содержит случайных величин. Если в уравнения (3) входят величины, имеющие вероятностный характер, то такие системы называются стохастиче- скими. Рассмотрим детерминированные системы и проведём их классификацию. Система (3) называется конечным автоматом, если она дискретна по времени и множества U, Y и X имеют конечное число элементов. В теории конечных автоматов множества U и Y называются входным и выходным алфавитами. 43 Система (3) называется конечномерной, если X, Y и U являются конечномерными линейными пространствами. Система называется стационарной, если функции f и g не зависят от t. В этом случае модель системы задаётся уравнениями xt +1 = f ( x t , u t ) , y t = g ( xt ) (4) Система называется линейной, если множества X, Y и U являются линейными пространствами, отображение f линейно по x и u, отображение g линейно по x. В этом случае математическая модель имеет следующий вид (если в пространствах X, Y и U выбраны базисы и их элементы представлены арифметическими векторами, а линейные операторы – матрицами) x t +1 = F t x t + G t u t + wt yt = H t x t + v t , где x t ∈ R N , y t ∈ R L , u t ∈ R M . Вход u t может формироваться человеком, воздействующим на систему (оператор, пилот, директор фирмы и т.п.). Случайный вектор w t является воздействием окружающей среды на эту систему и управлять его значениями невозможно (шум системы), v t - случайный вектор (шум наблюдения). Тогда структура линейной системы может быть изображена следующим образом 5.3. Синтез конечных автоматов 44 Конечное число состояний (а также значений входа и выхода) конечного автомата удобно обозначать (кодировать) с помощью набора двоичных (1 или 0) (булевых, логических) переменных. Например, если 4 возможных состояния, то достаточно двухмерных векторов: 0 0 1   1 S1 =   , S 2 =   , S 3 =   , S 4 =   . 0 1  0  1 Это удобно и с точки зрения его физической (электронной цифровой) реализации, когда используют компоненты, производящие операции над булевыми сигналами (напряжение «есть» = 1 или «нет» = 0). Поэтому алгоритм работы автомата целесообразно представлять в виде булевой функции. Рассмотрим частный случай стационарного конечного автомата, когда состояние x t +1 однозначно определяется входом u t . Тогда f ( xt , u t ) не зависит от x t и y t +1 = g ( f (u t )) = h(u t ) . Если в реальном автомате время срабатывания пренебрежимо мало по сравнению с циклом изменения времени, то можно считать xt + 0 = f (u t ) и он может быть описан моделью y t + 0 = h(u t ) . Так мы приходим к определению. Комбинационным автоматом (K-автоматом) или автоматом без памяти называется объект или модель, который имеет m входов и l двоичных выходов, причём каждой комбинации значений входов в любой момент времени соответствует одна определённая комбинация значений выходов. Вспомним операции над булевыми переменными. Дизъюнкция (логическое сложение «или – или») ∨ (+), конъюнкция (логическое умножение «и – и») ∧ ( ⋅ или знак вообще опускается) и дополнение (логическое отрицание «не») ¬ (––). Правила: 0+1=1, 1+1=1, 0+0=0; 11=1, 01=0, 00=0, 0 = 1, 1 = 0 . Легко проверить тождества (u + v) w = uw + vw , u + v = uv , u + vu = u + v . Докажем второе. Пусть л.ч. = 1, → u + v = 0 → u = 0, v = 0 → u = 1, v = 1 → п.ч. = 1. Аналогично доказывается, если п.ч. = 1 то л.ч. = 1. Если любая из них = 0, то и другая = 0, т.к. если бы была 1, то и первая должна была быть 1. Пример 5.2 Построить K -автомат, реагирующий (дающий на выходе 1) не менее, чем на два сигнала  u1    (= 1) из трёх, поступающих в качестве входа u =  u2  , y – его выход. Ясно, что u   3 y = h(u ) = u1u2u3 + u1u2u3 + u1u2u3 + u1u2u3 . 45 Пример 5.3 С помощью блоков, изображающих компоненты производящие операции сложения и умножения, построить физический K-автомат, реализующий функцию y = h(u ) = u1 (u2u3 + u4 ) . Или другая реализация: считая управления (= 1) - замыканием соответствующих выключателей, создать цепь, управляющую y (=1) – ток идёт. 1   Пусть u =   . Каков y?.... Ответ : 0. 1   0 Пример 5.4 Булевские уравнения конечного автомата x1t +1 = x1t x2t + u t , x2t +1 = x2t u t , y t = x2t . Будет ли это стационарная система? Линейная? K-автомат? Построить его физическую реализацию (используя в том числе блоки сдвига назад по времени (z-1) и отрицания ( ¬ )). 46 0 1 Пусть x 0 =   , u 0 = 1. Каков y 0 ? y1 ? x1 ? Пусть u1 = 1. Каков y 2 ? Ответы: 1, 0,   , 0. 1 0 Однако в практических ситуациях вначале даётся словесное описание закона поведения автомата, затем строится геометрическое представление поведения в виде сети из кругов и соединяющих их дуг, которая называется графом. Круги называются вершинами графа и изображают возможные состояния автомата, а дуги указывают переходы из состояния в состояние. Когда задан граф автомата, желательно построить соответствующие логические функции. Эта задача называется синтезом автомата. Пример 5.5 В результате наблюдений за «чёрным ящиком» с одним входом и одним выходом было установлено, что у него возможны 4 состояния, между которыми происходят переходы под влиянием управления u = (0, 1). А именно: S1→0→S2, S1→1→S4; S2→0→S1, S2→1→S3; S3→0→S2, S3→1→S4; S4→0→S4, S4→1→S4. Получаем граф 47 Теперь можно решать задачи управления системой. Например, автомат находится в состоянии S2, как его перевести в состояние S4 и поддерживать в нём? Видно, что можно применить управление 0, в следующем такте 1 и в следующем такте система будет в требуемом состоянии S4. Затем нужно постоянно подавать управление 0 или 1. Для решения задачи синтеза автомата вначале необходимо закодировать его состояния векторами, каждая координата которых является логической переменной. Так как состояний 4 то  x1t  0 0 1   1 достаточно двухмерных векторов состояний: x =  t  , S1 =   , S 2 =   , S 3 =   , S 4 =   . x  0 1  0  1  2 t Далее строятся функции f1 ( x t , u t ) и f 2 ( x t , u t ) . Первая координата. При каких условиях x1t +1 примет значение равное единице? Т.е. после управления перейдём в состояние, в котором первая координата = 1, а это S3 или S4. По графу видно, что это будет когда предыдущее состояние было S2 и поступило управление 1 (логическая функция x1 x2u ), либо предыдущее состояние было S3 и поступило управление 1, либо предыдущее состояние было S1 и поступило управление 1, либо предыдущее состояние было S4 и поступило управление 0 или 1. Т.о. 48 f1 ( x , u ) = x1 x2u + x1 x2u + x1 x2u + x1 x2u + x1 x2u = x1 x2u + ( x1 + x1 ) x2u + x1 x2u + x1 x2u = x2u + x2u + x1 x2u = u + x1 x2u = u + x1 x2 . Аналогично вторая координата. Т.о. результат будет x1t +1 = x1t x2t + u t , x2t +1 = x1t + x2t . 0 Проверим, x t = S1 =   , u t = 1 . Куда перейдёт система? 0  00 + 1 1 x t +1 =   =   = S 4 т . д.  0 + 1   1 Пример 5.6 На конвейере детали A и B образуют беспорядочную последовательность. Для последующего использования детали следует рассортировать так, чтобы они на конвейере образовывали тройки ABA ABA ABA … Сконструировать автомат, который делает это, сбрасывая ненужную деталь на вспомогательный конвейер, который возвратит детали на главный. Воздействия – приходящие детали (A, B), состояния S1 – исходное (готов пропустить A), S2 – только что пропустил первую A тройки, S3 – только что пропустил B тройки. Со штрихом – те же самые, но люк сброса открыт (деталь сбросилась). 49 Пример 5.7 (Росс Эшби У. Введение в кибернетику. М., Изд-во иностр. лит., 1959, с. 92— 93) Дорогой друг! Некоторое время назад я купил старый дом, но обнаружил, что он посещается двумя призрачными звуками: Пением и Смехом. В результате он мало подходит для жилья. Однако я не отчаиваюсь, ибо я установил путём практической проверки, что их поведение подчиняется определёнными законам, непонятным, но непререкаемым, и что я могу воздействовать на них, играя на Органе или сжигая Ладан. В течение каждой минуты каждый из этих звуков либо звучит, либо молчит; никаких переходов они не обнаруживают. Поведение же их в последующую минуту зависит только от событий предыдущей минуты, и эта зависимость такова: Пение в последующую минуту ведёт себя так же, как и в предыдущую (звучит или молчит), если только в эту предыдущую минуту не было игры на Органе при молчащем Смехе. В последнем случае оно меняет своё поведение на противоположное (звучание на молчание и наоборот). Что касается Смеха, то если в предыдущую минуту горел Ладан, Смех будет звучать или молчать в зависимости от того, звучало или молчало Пение (так что Смех копирует Пение минутой позже). Если, однако, Ладан не горел, Смех будет делать противоположное тому, что делало Пение. В ту минуту, когда я пишу Вам это, Смех и Пение оба звучат. Прошу Вас сообщить мне, какие действия с Ладаном и Органом должен я совершить, чтобы установить и поддерживать тишину в доме? 5.4. Основные проблемы теории систем В теории систем изучаются три группы проблем: (1) это выявление и представление в тех или иных формах законов динамики систем, (2) выявление или оценивание текущих состояний, в частности, для решения проблем прогнозирования эволюции систем и (3) формирование входов, обеспечивающих требуемое поведение систем. Целью при этом является как выяснение условий, при которых эти проблемы разрешимы, так и получение методов из решения. Первую группу проблем называют проблемами представления или идентификации си- стем. Вторую группу называют проблемами наблюдаемости и третью – управляемости. Проблема идентификации, по существу, связана с описанием множеств Y и U, используя имеющуюся информацию о реальных процессах и их взаимодействиях. Далее встаёт задача построения множества состояний X и отображений f и g. 50 Если X, f и g заданы или идентифицированы, то можно приступить к изучению одной из основных проблем – проблеме прогнозирования. Обычно представляет интерес прогноз наблюдаемого поведения, т.е. выхода y(t). Однако, учитывая связь (1) эту задачу можно свести к прогнозу процесса в пространстве состояний. По (3) имеем, что прогнозирование состояния в будущий момент окажется возможным, если в состояние известно в предыдущий момент. Поэтому при известной модели проблема прогнозирования сводится к задаче нахождения (или оценки) состояния системы в требуемые моменты времени. Подчеркнём, что реальная информация о системе связана лишь с измерениями процессов входа и выхода. Пример: пусть для нахождения состояния линейной системы в момент t1 используется информация, даваемая выходом только в тот же момент. Тогда задача сводится к разрешимости относительно x t1 уравнения y t1 = H t1 x t1 + v t1 . Если размерности X и Y одинаковы, и матрица H не вырождена, то из него можно определить единственное решение x t1 . Однако обычно размерность X больше размерности Y. Тогда решение этого уравнения не единственно, т.е. наблюдаемому значению y t1 будут соответствовать различные состояния, при которых оно может реализоваться. В таком случае для различения состояний необходимо расширить обследуемую информацию, наблюдая выход более, чем в один момент, например, на интервале [t0, t1]. Система, для которой при этом возможно восстановить процесс x[t0, t1] называется наблюдаемой. Следующая основная проблема теории систем – формирование специального поведения систем, что вызывается необходимостью удовлетворения заданных требований, накладываемых на протекающие в ней процессы. Эти требования называют целью, которая ставится перед системой, или её назначением. Требования на поведение системы, которые обычно задаются в терминах процесса-выхода, целесообразно формулировать в виде требований на процесс в пространстве состояний. Одной из важнейших задач теории управления является так называемая двухточечная задача. Пусть цель состоит в том, чтобы в момент t1 > t0 система находилась в состоянии x1, если она была в состоянии x0. Если её можно привести в любое состояние – она называется дости- жимой. Если её можно привести из любого состояния – она называется управляемой. Наконец, проблема устойчивости системы. Будет ли выход системы, подвергшейся возмущению сходиться к невозмущённому выходу при t → ∞ ? 5.5. Передаточная матрица системы 51 Рассмотрим линейную стационарную систему x t +1 = Fx t + Gu t , (5) y t = Hx t . (6) Хотелось бы найти замкнутое (аналитическое) выражение того, как выходной процесс y t зависит от входного (управления) u t . Для этого нужно решить разностное уравнение (5) и подставить решение в (6). Разностные уравнения решаются с помощью z-преобразования ∞ X Z ( z ) =  xt z −t , (7) t =0 xt = s + iπ 1 X Z (e p )e pt dp .  2π i s −iπ (8) Применим (7) к обеим частям (5) и (6) (стрелки у векторов будем опускать). ∞ ∞ ∞ t =0 t =0 t =0  xt +1 z −t = F  xt z −t + G  u t z −t , ∞ ∞ t =0 t ′= 0 YZ ( z ) = HX Z ( z ) z ( x t +1 z − (t +1) + x 0 − x 0 ) = z ( xt ′ z − t ′ − x 0 ) , т .о . , л.ч. первого выражения zX Z ( z ) − zx 0 = FX Z ( z ) + GU Z ( z ) , равна отсюда X Z ( z ) = ( zI − F ) −1 GU Z ( z ) + z ( zI − F ) −1 x 0 и YZ ( z ) = H ( zI − F ) −1 GU Z ( z ) + Hz ( zI − F ) −1 x 0 . (9) Для x 0 = 0 получаем YZ ( z ) = W ( z )U Z ( z ) , (10) W ( z ) = H ( zI − F )−1 G (11) где называется передаточной матрицей системы (5) и (6). Если только u 0 = 1 а все остальные = 0 (единичный импульс), то по (7) U Z ( z ) ≡ 1 и YZ ( z ) = W ( z ) , т.е. W ( z ) - образ отклика системы на единичный импульс (характеристическая функция). Её обратное z-преобразование (8) (отклик) называется импульсной характеристикой и если содержит конечное число ненулевых отсчётов, то говорят о системе с конечной импульсной характеристикой (КИХ), в противном случае с бесконечной – БИХ. Комплексная частотная характеристика КЧХ ( ω ) = W (eiω ) показывает в, соответствии с (10), как составляющие различных частот управления передаются системой на выход. Её модуль называется амплитудно-частотной характеристикой - АЧХ(ω). Пример 5.8 52 Вход-выход модель порядка N = 2. a0 y t + a1 y t +1 + y t + 2 = c0u t + c1u t +1 (12) Её можно переписать в пространстве состояний (5), (6)  0 −a0  F =   1 −a1  c  c0 + c1 z G =  0  , H = ( 0 1) . По (11) W ( z ) = . Пусть a0 = 0.7, a1 = −0.5 a0 + a1 z + z 2  c1  c0 = −0.5, c1 = 1 . АЧХ - см. на рисунке. На следующих рисунках изображён отклик этой системы на различные воздействия. По входной последовательности по формуле (7) вычислялось её z-преобразование, находился об-  0 раз выхода по (9), x 0 =   . Выходная последовательность вычислялась по формуле  0 yt = 1 2π π π Y Z (e s +iω )e ( s +iω ) t d ω численно при s = 0.1. Расчёты проведены в PTC Mathcad Prime 4.0. − При подаче единичного импульса в момент 0, имеем. Т.е. это БИХ - система. Если 2 одиночных импульса. 53 Максимум АЧХ при частоте ω0 = arccos( −a1 ) = 1.267 (период = 4.958) соответствует резо2 a0 нансной частоте осциллятора в л.ч. (12). Поэтому управления с близким периодом (5, частота ω = 2 π/5 = 1.257) дают большой выход. В PTC Mathcad Prime 4.0 есть встроенная функция response (u, C1, τ), вычисляющая τ значений выходной величины y t по входным отсчётам u и коэффициентам C1 вход-выход модели (12) (или, что тоже, дробно-рациональной характеристической функции W ( z ) ). Так в этом 0  примере C1 =  c1 c  0 1  a1  , и получим те же результаты что и на последних 3-х рисунках. a0  5.6. Макроэкономические модели Модель Леонтьева (~1930 г) 54 Пусть есть n отраслей экономики, каждая производит по одному продукту. yi – валовый выпуск i-го продукта (в натуральных единицах тонн, шт. и т.п.), vij – производственные заn траты i-го продукта в j-й отрасли. vi =  vij – затраты i-го продукта во всех отраслях (т.е. на j =1 производство и самого i-го продукта и всех других) - прямые затраты. Ясно yi = vi + xi , xi – конечный продукт i-й отрасли. Будем считать vij , естественно, пропорциональными объёмам y j . Т.о. vij = aij y j или aij = vij / y j – затраты единиц i-го продукта на производство единицы jго продукта. A = (aij ) – матрица технологических коэффициентов. Подставляя, имеем n  y =  i  aij y j + xi , или j =1  a ≥ 0, y ≥ 0  ij  y = Ay + x   A ≥ 0, y ≥ 0 (13) – статическая модель Леонтьева. Модель Леонтьева называется продуктивной, если ∃x > 0 и ∃y ≥ 0 удовлетворяющие (13). Продуктивность означает, что есть, что потреблять по всем видам продукции при некотором валовом выпуске. Матрица продуктивной модели называется продуктивной матрицей. Пример 1.5 A = 0  x = y – продуктивна. Пример 2.4 1 ... 1   A = 1 ... 1 ? (Указание: от противного, предположив ∃x > 0 , расписать (13) покомпо1 ... 1   нентно и сложить равенства). Теорема Для продуктивности модели Леонтьева необходимо и достаточно ∃( I − A) −1 ≥ 0 . Доказательство достаточности. Из (13) следует ( I − A) y = x , ∃( I − A) −1 ≥ 0  y = ( I − A)−1 x  ∀x > 0∃! y = ( I − A)−1 x ≥ 0. Т.о. для продуктивной модели можно задавшись вектором конечных продуктов найти необходимый валовый выпуск y = Nx , 55 (14) где N = ( I − A)−1 -- матрица полных затрат. Если хотим увеличить конечный продукт на ∆x , то каково соответствующее увеличение вала? y + ∆y = A( y + ∆y ) + x + ∆x  ∆y = A∆y + ∆x  ∆y = ( I − A) −1 ∆x Пример 1.6 Две отрасли 1-я сельское хозяйство (лён), 2-я промышленность (ткань). Матрица межотраслевого баланса: прямые затраты конечное потребле- валовая сельское промыш- ние (домохозяйства) хозяйство ленность сельское хозяйство 25 20 ? 100 кг промышленность 14 6 ? 50 м продукция Каково конечное потребление? Как должен измениться валовый выпуск, если потребление ткани должно вырасти на 20%? Двойственная модель Леонтьева Пусть по-прежнему есть n отраслей экономики, каждая производит по одному продукту. Обозначим pi – цена за единицу i-го продукта. Каковы должны быть цены в равновесной экономике? Ясно, цена = производственные затраты на сырьё, комплектующие, энергию и т.п. (продукты других отраслей) + добавленная стоимость (на зарплату, налоги, …). pi = a1i p1 + a2i p2 + ... + ani pn + wi , wi – добавленная стоимость в единице i-го продукта. pi =  a ji p j + wi =  ( AT )ij p j + wi , j j T  p = A p + w   A ≥ 0, p ≥ 0 (15) – двойственная модель Леонтьева. Двойственная модель Леонтьева называется продуктивной, если ∃w > 0 и ∃p ≥ 0 , удовлетворяющие (15). Теорема 56 Прямая и двойственная модели Леонтьева продуктивны одновременно (или не продуктивны). ( I − AT )−1 = [( I − A)T ]−1 = [( I − A)−1 ]T ≥ 0 ⇔ ( I − A)−1 ≥ 0 . Это соответствие между натуральной и стоимостной сторонами экономики. Дискретно-временная динамическая модель Леонтьева В ней предполагается, что конечный продукт в году t x t используется для инвестиций i t и конечного потребления c t xt = i t + c t . (16) Инвестиции дают прирост валового выпуска ∆y t = ( y t +1 − y t ) к следующему году i t = B ∆y t . (17) B – матрица капиталоёмкости прироста валовых выпусков. ii = bi1∆y1 + bi 2 ∆y2 + ... + bin ∆yn , i = 1,..., n , т.е её элементы bij – инвестиции единиц i –го продукта, необходимые для прироста единицы j –го продукта за единицу времени (год). Подставляя (16) и (17) в (13) получим y t = Ay t + B( y t +1 − y t ) + c t . Считая, что матрицы A и B могут меняться со временем (технический прогресс), и выражая y t через x t по (14) имеем x t +1 = F t x t + G t u t yt = H t x t , – вид уравнений управляемой системы. Управление u t = c t – конечное потребление, наблюдаемая выходная величина – валовый выпуск, состояние – конечный продукт, используемый для инвестиций и конечного потребления (ВВП). F t = ( Bt N t +1 ) −1 + ( N t +1 )−1 N t , G t = −( Bt N t +1 )−1 , Ht = Nt . Библиографический список 1. Акоф Р. Основы исследования операций / Р. Акоф, М. Сасиени. М.: Мир, 1971. 2. Бородачёв С. М. Имитационное моделирование в экономике : учебное пособие / С. М. Бородачёв. Екатеринбург : УрФУ, 2010. 3. Вагнер Г. Основы исследования операций. Т. 1-3. / Г. Вагнер. М.: Мир, 1973. 4. Глухов В. В. Математические методы и модели для менеджмента / В. В. Глухов, М. Д. Медников, С. Б. Коробко. С. Пб.: Лань, 2000. 57 5. Дегтярев Ю. И. Исследование операций / Ю. И. Дегтярев. М.: Высшая школа, 1986. 6. Дубров А. М. Моделирование рисковых ситуаций в экономике и бизнесе / А. М. Дубров, Б. А. Лагоша, Е. Ю. Хрусталев. М.: Финансы и статистика, 2000. 7. Исследование операций в экономике / Под ред. Н. Ш. Кремера / М.: ЮНИТИ, 1997. 8. Коршунов Ю. М. Математические основы кибернетики / Ю. М. Коршунов. М.: Энергоатомиздат, 1987. 9. Кофман А. Займёмся исследованием операций / А. Кофман, Р. Фор. М.: Мир, 1966. 10. Морозов В. В. Исследование операций в задачах и упражнениях / В. В. Морозов, А. Г. Сухарев, В. В. Фёдоров. М.: Высшая школа, 1986. 11. Александр Мороз. Математические основы менеджмента /А. И. Мороз. М.: Academia, 1997. 12. Росс Эшби У. Введение в кибернетику /М.: Изд-во иностр. лит., 1959. 13. Бородачёв С. М. Методы математической статистики : учебное пособие / С. М. Бородачёв. Екатеринбург: УрФУ, 2012. 14. Ивченко Г. И., Медведев Ю. И. Математическая статистика. – М.: Высшая школа, 1984. 58