Решение. Сделаем копию данной последовательности и отсортируем ее по возрастанию

Первый способ

Сделаем копию данной последовательности и отсортируем ее по возрастанию. Тогда НВП будет также НВП для отсортированной последовательности. Поэтому достаточно найти наибольшую общую подпоследовательность исходной и отсортированной последовательности. Обратите внимание: на самом деле таким способом мы найдем не возрастающую, а неубывающую подпоследовательность. Чтобы найти возрастающую подпоследовательность, нужно сделать еще один шаг. Сложность этого алгоритма – O(n²).

Второй способ

Попробуем воспользоваться методом динамического программирования. Пусть мы знаем НВП для последовательности a₁, …, a_k. Можем ли мы найти НВП для последовательности a₁, …, a_k+1? Кажется, что это легко сделать: если a_k+1 больше, чем последний член ранее найденной НВП, то a_k+1 нужно добавить в НВП, в противном случае нужно оставить НВП неизменной. К сожалению, такое решение оказывается неверным. Рассмотрим такой пример: 1 2 3 6 4 5. Для первых пяти членов этой последовательности длина НВП равна 4, но есть две различных НВП: 1 2 3 6 и 1 2 3 4. Пусть мы выбрали первую из них. Тогда на следующем шаге мы не сможем дополнить ее числом 5 и получим более короткую НВП (1 2 3 6) для всей последовательности, чем должны был получить (1 2 3 4 5). Модифицируем предложенный алгоритм. Будем теперь запоминать не произвольную НВП, а ту, которая заканчивается на наименьшую возможную цифру. К сожалению, и этот способ не приведет нас к успеху. Рассмотрим такой пример: 1 4 2 3. Для последовательности из первых двух символов НВП будет иметь длину 2 (1 4). Но при добавлении символов длина НВП не увеличится, поскольку все оставшиеся числа меньше четырех. Для решения задачи нам потребуется более сложная система подзадач. Пусть для последовательности a₁, a₂, …, a_k мы знаем возрастающие последовательности длины 1, 2, 3,... Причем среди всех возрастающих последовательностей данной длины мы будем хранить ту, которая оканчивается на меньшее число. Тогда для последовательности a₁, a₂, …, a_k+1 мы сможем построить аналогичный набор последовательностей.

Третий способ

Будем постепенно строить возрастающие цепочки и в каждый момент времени для каждого элемента a_i хранить длину самой длинной построенной цепочки, заканчивающейся этим элементом. Приведенный алгоритм имеет сложность O(n²).

Четвертый способ

Будем поочередно обрабатывать члены исходной последовательности. Пусть после обработки k-го элемента нам известно, что возрастающая подпоследовательность длины 1 может оканчиваться числом e[1], ВП длины 2 может оканчиваться числом e[2] и т.д., причем указанные числа e[i] – наименьшие возможные. Заметим, что е[1] < e[2] < e[3] < … Рассмотрим теперь число a_k+1. Найдем такое i, что e[i-1] ≤ a_k+1 < e[i] (если изначально положить e[0] = -∞, а все остальные члены массива e равными +∞, то такое всегда найдется). Положим e[i]: = a_k+1, а остальные элементы массива e[i] оставим без изменений. Заметим, что сложность поиска в упорядоченном массиве есть O(log m), где m – длина НВП, поэтому сложность приведенного алгоритма – О(n log m).

Примеры задач, решаемых при помощи динамического программирования. Наикратчайшие пути через сети

Пусть имеется N узлов с номерами 1, 2, …, N. Пусть t_ij> 0 – время прохождения от узла i к узлу j. Пусть N – целевой узел, а u_i – минимальное время перехода от узла i к узлу N. Очевидно, что u_N=0. По принципу оптимальности получим:

u_i = min_{j¹ i}(t_ij + u_j), i = 1, 2, …, N-1.

Рассмотрим некоторые свойства полученных уравнений. Пусть (u₁, u₂, …, u_N) и (U₁, U₂, …, U_N) – два решения и m – индекс, для которого максимальна разность

U_j – u_j, j = 1, …, N.

Пусть

U_m = t_mn + U_n £ t_mr + U_r

u_m = t_mr + u_r Þ U_m - u_m £ U_r - u_r

Так как m – индекс для которого U – u = max, то получаем знаки равенства. Точно также можно найти s, s¹ m, s¹ r и U_m – u_m = U_r – u_r = U_s – u_s Þ U_N – u_N = 0.

Метод итераций. Пусть

u_i⁽⁰⁾ = t_iN, t_NN = 0,

если нет пути от i к N, то положим t_iN = ¥. Рассмотрим итерационный процесс:

u_i⁽¹⁾ = min_{i¹ j}(t_ij + u_j⁽⁰⁾), u_N⁽¹⁾ = 0, …

Очевидно, что u_i^(k+1) £ u_i^(k)). Так как оптимальный путь содержит не более N-2 узла (поскольку нет петель), то u_i^(N-2) есть решение Þ число итераций не более N-2.

Оптимизация метода динамического программирования (на примере задачи нахождения числа сочетаний )

Одна из проблем рекурсивных алгоритмов – длительная работа рекурсии за счет повторяющихся вызовов рекурсивной функции для одного и того же набора параметров. Чтобы не тратить машинное время на вычисление значений рекурсивной функции, которые мы уже когда-то вычислили, можно сохранить эти значения в массиве и вместо рекурсивного вызова функции мы будем брать это значение из массива. В задаче вычисления числа сочетаний создадим двумерный массив B и будем в элементе массива B[n][k] хранить величину C _n^k. В результате получим следующий массив (по строчкам – значения n, начиная с 0, по столбцам -- значения k, начиная с 0, каждый элемент массива B[n][k] равен сумме двух элементов: стоящего непосредственно над ним B[n-1][k] и стоящего над ним слева B[n-1][k-1]).

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

Приведенный алгоритм для вычисления C _n^k требует объем памяти O (n ²) и такая же его временная сложность (для заполнения одного элемента массива требуется O (1) операций, всего же элементов в массиве O (n ²)).

Подобный метод решения, когда одна задача сводится к меньшим подзадачам, вычисленные же ответы для меньших подзадач сохраняются в массиве или иной структуре данных, называется динамическим программированием. В рассмотренной задаче вычисления числа сочетаний использование метода динамического программирования вместо рекурсии позволяет существенно уменьшить время выполнения программы за счет некоторого увеличения объема используемой памяти. Тем не менее, приведенный алгоритм тоже имеет небольшие недостатки. Мы вычисляем значения всех элементов последней строки, хотя нам необходим только один из них – B[n][k].

Аналогично, в предпоследней строке нас интересуют только два элемента – B[n-1][k-1] и B[n-1][k], в строке, стоящей над ней – три элемента, и т. д., в то время, как мы вычисляем все элементы во всех строках. Кроме того, мы не используем почти половину созданного массива B[n+1][n+1], а именно все элементы, стоящие выше главной диагонали. От этих недостатков можно избавиться, более того, можно уменьшить объем используемой памяти до O (n).

Если же в программе часто возникает необходимость вычисления числа сочетаний C _n^k для каких-то ограниченных значений n, то лучше всего создать глобальный массив для хранения всевозможных значений C _n^k для нужных нам значений n, заполнить его в начале программы, а затем сразу же брать значение числа сочетаний из этого массива, не вычисляя его каждый раз заново.

Тема 4. Системы и процессы.

План лекции:

Системы.

Многошаговые процессы.

Редукция информации.

Реккурентные соотношения.

Бесконечные процессы.

Процессы с заданными правилами остановки.

Нестационарные процессы.

Краткое содержание лекции

Системы

Сначала нас будет интересовать изучение систем (термин, который здесь еще не определен, но интуитивно вполне понятен) и то, как они ведут себя с течением времени. Эта точка зрения позволяет ввести понятие «процесса» простым и естественным способом. Аналитически представим себе систему как вектор состояния х(t) и правило для определения его значения в любой момент времени t. Для простоты будем рассматривать только конечномерные векторы

Каждая компонента х(t) определяет различные свойства системы. Число N называется размерностью системы. Классически х(t) может быть задан следующими способами:

= g(x(t)) – дифференциальные уравнения,

x(t)=g(x(t-t)) – разностное уравнение,

= h(x(t), x(t-t)) – дифференциально-разностное уравнение,

x(t)=u(t) + ò g(x(s), t)ds – интегральное уравнение,

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

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

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

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

Многошаговые процессы

Введя чисто математическое определение понятия «система», уточним теперь другую интуитивную концепцию – понятие «процесс». Ради общности заменим символ х(t) символом р, считая р точкой множества или пространства R. Во всем последующем изложении R будет множеством точек М -мерного пространства (возможно с некоторой метрикой) или некоторого подмножества.

Рассмотрим теперь функцию Т(р) – преобразование, обладающее свойством, что преобразованная точка р₁ = Т(р) принадлежит R для всех точек р из R. Можно считать, что р изображает начальное состояние системы, р₁ = Т(р) – состояние на одну единицу времени позднее и вообще система векторов [ р, р₁, р₂, …, р_n, … ], где р₀ = р, р_{n+ 1}= Т(р_n), n =0, 1, 2, …, изображает временную историю системы, наблюдаемую в дискретные моменты времени n =0, 1, 2, …, т.е. последовательные состояния системы. Можно также написать р_n=Тⁿ(р), что обозначает n -кратное применение преобразования Т, и иногда будем применять это обозначение.

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

Удобно ввести обозначение [ р, Т(р) ] или просто [ р, Т ]. более точно, это есть многошаговый процесс дискретного детерминированного типа. Так как это единственный многошаговый процесс, который мы указали, то ограничимся пока более простым термином, а при рассмотрении других типов многошаговых процессов отметим их дополнительные особенности.

Редукция информации

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

Для начала будем рассматривать только часть многошагового процесса, N -шаговый процесс, последовательность векторов [ p, p₁, p₂, …, p_N ], где, как и прежде, р_{k +1}= T(p_k), k =0, 1, 2, …, N -1. Рассмотрим свойства различных скалярных функций, зависящих от этого процесса, обозначим их { g(p, p₁, …, p_N) }.

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

, …

Независимость от прошлого

Мы предполагаем, что в начальный момент времени система определяется вектором состояния р, а процесс описывается последовательностью:

р ₁= Т(р), р ₂= Т(р ₁ ), …

Как следствие этого предположения последующее для k -й точки состояние зависит только от p_k, и не зависит от предыдущих состояний. Другой путь введения этих понятий состоит в том, что p_N можно рассматривать как N -е состояние на N -м шаге процесса, начавшегося из начального положения p, или (N-k) -e состояние процесса на (N-k) -м шаге процесса, начавшегося из состояния p_k в момент времени k. Символически это можно записать в виде уравнения T^N = T^N-kT^k. Это соотношение может рассматриваться как аналитическая формулировка причинной связи, или, другими словами, того факта, что будущее единственным образом определяется настоящим (теорема единственности).

Реккурентные соотношения

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

,

где h(p) – заданная функция состояния, и заметим, что рассматриваемая функция полностью определена начальным состоянием р и числом шагов N, если задано преобразование Т. Запишем последовательность функций { f_N(p) }, определяемых соотношением

f_N(p) = = h(p) + h(T(p)) + h(T ² (p)) + … + h(T^N(p)),

где N =0, 1, 2, … и р принадлежит множеству R. Рассмотрим для N ³ 1 частичную сумму

h(T(p)) + h(T ² (p)) + … + (T^N(p)).

Эта сумма эквивалентна

h(p ₁ ) + h(T(p ₁ )) + … + h(T^{N- 1} (p ₁ )),

а следовательно, и функции f_N-1(p₁), определяемой последовательностью функций {f_N(p)}. Из этого следует, что f_N(p) удовлетворяет функциональному (реккурентному) уравнению

Аналогично, если обратится к функции , то рассмотрим

f_N(p) = ,

и получим функциональное уравнение

f_N(p) = h(p)f_N-1(T(p)).

Предоставим читателю вывод соотношения f_N(p) = max [ h(p), f_N-1(T(p)) ], N ³ 1,

для функции f_N(p) = , определенной выше под пунктом (с).

Таким же образом соотношение

f_N(p) = h(p, T(p)) + f_N-1(T(p))

может быть установлено для функции под пунктом (d).

Бесконечные процессы

Соотношения, полученные в предыдущем пункте, значительно упрощаются для бесконечных процессов, т.е. процессов, которыми становятся ранее рассмотренные многошаговые процессы при N =¥. Поступим формально и, не рассматривая пока вопросов сходимости, напишем f(p) = . Заметим, что бесконечность числа шагов исключает зависимость от N и таким образом устраняет одну переменную, что весьма существенно. При этом аналогично предыдущему получим

f(p) = h(p) + , f(p) = h(p) + f(T(p)).

Подобным же образом, если f(p) = , то, учитывая очевидное соотношение

max [ x₁, x₂, …]= max [ x₁, max [ x₂, …]],

получаем

f(p) = max [ h(p), f(T(p)) ].

Во многих случаях функция f_N(p) не имеет предела при N ® ¥. В этом случае иногда имеет смысл величина вида с асимптотическим соотношением вида ~ g₁(p) при N ® ¥ или соотношением вида f_N(p) ~ g₂(p)^Nh(p).

Процессы с заданными правилами остановки

Во многих важных ситуациях в рассматриваемых конечных процессах число шагов заранее не фиксировано, а зависит от начального состояния, т.е. N = N(p). При этом приходим к соответствующим функциям вида f(p) = .

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

£ a.

Другой немаловажный пример дает нам исследование траекторных процессов. Объект продолжает движение до тех пор, пока он не окажется на некотором заранее обусловленном расстоянии от другого объекта. Аналитически, вводя метрику в R, обозначая через d(p, q) расстояние между двумя элементами p и q в R, мы ставим условие, что процесс продолжается до тех пор, пока d(T^Np, q) £ e, где q и e заданы заранее. Рассмотрим траекторный процесс, который оканчивается, как только d(p_N, q) £ e. Если ввести функцию f(p)= , где N = N(p) определяется приведенным выше правилом остановки, то получаем

Нестационарные процессы

Нестационарные процессы имеют вид:

(p_m, р_m+1, р_m+2, …, р_n, …),

где

р_{n+ 1}= Т_n(р_n), n = m, m+1, m+2, …

Выполнив редукцию для величины получаем:

f_{N, m}(p) = h_m(p) + f_{N, m+1}(T_m(p))

Если N = const, то j_m(p) = f_{N, m}(p) и получаем:

j_m(p) = h_m(p) + j_m+1(T_m(p)), 0£ m£ N

j_N(p) = h_N(p)

Тема 5. Дискретные и непрерывные задачи.

План лекции:

Непрерывные процессы.

Траекторный процесс.

Неоднородная атмосфера.

Причинность.

Стохастичность и корреляция.

Вариационное исчисление.

Геометрическая интерпретация.

Уравнение эйконала (в R²).

Краткое содержание лекции

Непрерывные процессы

Рассмотрим дискретный процесс, заданный соотношениями:

p_{i +1} = T (p_i)

Пусть этим значениям i соответствует значения t = 0, D, 2D, …

Тогда в силу непрерывности T(p)» p, а точнее

T(p) = p + S(p)D + o(D)

Пусть соответствующая целевая функция имеет вид суммы

f_t(p) = g(p)D + g(T(p))D + g(T²(p))D + … + g(Tⁿ(p))D

И предположим, промежуток D достаточно мал. Если обозначим t = nD, то получим

f_t+D(p) = g(p)D + f_t(p + S(p)D + o(D))

Рассмотрим теперь M-мерный случай, получим

Тогда уравнение для целевой функции примет вид

Если перейдем к пределу при D®0, то получим линейное уравнение в частных производных

Формально

p¢ (t) = S(p),

p(0) = p₀

Траекторный процесс

Рассмотрим задачу: в пустоте с малой скоростью v брошено вертикально вверх тело. Найти наибольшую высоту подъёма тела. Классическое решение задачи имеет вид:

Способ динамического программирования. Пусть f(v) – максимальная высота

Последовательно получаем:

D ® 0

T(v) = v - gD

f(v) = vD + f(v - gD) + o(D)

D ® 0 Þ 0 = v - f¢ (v)g, f(0) = 0 Þ

Задача (объект А догоняет объект Б). Пусть f(r, d) – время до встречи, g(r, d) – расстояние до точки встречи.

Получаем уравнение

Если D ® 0, то разложив f в правой части по формуле тейлора, получим уравнение в частных производных для функции f и граничное условие для f:

Найти f(r, d); составить уравнение для g(r, d); найти g(r, d) (методом характеристик).

Неоднородная атмосфера

Пусть теперь задача о теле, брошенном вверх, имеет вид (движение с сопротивлением):

x² = h(x, x¢),

x(0) = 0,

x¢ (0) = v

Следовательно нужна дополнительная переменная состояния – высота в начале процесса (b). Пусть f(b, v) – максимальная высота.

Последовательно получаем

b₁ = T₁(b, v) = b + vD + o(D)

v₁ = T₂(b, v) = v + h(b, v)D + o(D)

f(b, v) – max высота Þ

f(b, v) = vD + f(b+vD, v+h(b, v)D) + o(D)

D ® 0 Þ

,

f = -b + j(const),

где const – первый интеграл.

Самостоятельно рассмотреть конкретный пример h(x, x¢) = -g - kx¢.

Причинность

Рассмотрим последовательно идущие стадии процесса

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

f(b, t+t) = f(f(b, t), t)

это и есть принцип детерминизма. Например для матрицы A получаем

e^A(t+t) = e^Ate^At

Если e^At – есть решение уравнения x¢ = Ax, x(0) = E.

Самостоятельно рассмотреть случаи sin(t+t) и cos(t+t)

Стохастичность

Пусть T не полностью известно, тогда (p₀, р₁, р₂, …) – случайный процесс, следовательно получаем систему с описанием

p_k = T(p_k-1, r_k), k³ 1,

p₀ = p,

где r_k - независимые случайные величины с распределением dG(r).

Требуется найти математическое ожидание.

M(g(p)) = M(g(p₁) + g(p₂) + …+g(p_N))

Обозначим

f_N(p) = M_r(p_m, р_m+1, р_m+2, …, р_n, …), N³ 1,

где математическое ожидание берётся относительно независимых случайных величин r₁, r₂, …, r_N. Тогда получим

f_N(p) = g(p) + ò f_N-1(T(p, r₁)dG(r₁)

f₀(p) = g(p)

(т.к. f_N(p) = M_r(g(p) + f_N-1(T(p, r₁)))

Корреляция

Случайные величины r_k могут быть зависимы, простейший случай, когда r_k зависит только от r_k-1. В этом случае известно распределение dG(r_k, r_k-1). Тогда

f_N(p, r₀) = g(p) + ò f_N-1(T(p, r₁))dG(r₁, r₀)

Вариационное исчисление

Простейшая классическая задача вариационного исчисления имеет вид: найти max J(u), где

J(u) = ò ₀^t0g(u, u¢)dt,

u(t) Î C_{[0, t0]},

u(0) = c.

Обозначим формально

f(c, t₀) = max_uJ(u)

и воспользуемся свойством интеграла ò ₀^t0 = ò ₀^D + ò _D^t0.

Учтя, что

ò ₀^Dg(u, u¢)dt = g(c, v)D + o(D), (где v = u¢ (0)) получим:

f(c, t₀) = max_v(g(c, v)D + f(c+vD, t₀-D)) + o(D)

И далее учитывая, что D ® 0 получаем:

f¢ _t0 = max_v(g(c, v) + vf¢ _c), f(c, 0) = 0

Замечание. v = v(c, t) называется функцией стратегии

Замечание. В классическом подходе решается уравнение:

При этом если, например, дополнительно задано условие |u¢ (t)|£ k, 0£ t£ t₀, то в классическом подходе возникают серьёзные трудности, а метод динамического программирования легко даёт уравнение:

f¢ _t0 = max_{|v|£ k}(g(c, v) + vf¢ _c),

f(c, 0) = 0

Геометрическая интерпретация

Вариационное исчисление: ищем кривую u(t), tÎ [0, T], такую, что J(u) ® max, а u – точка в пространстве функций. Динамическое программирование: для каждой точки ищем направление, которое является оптимальным.

Иначе: кривая – геометрическое место точек (В. И.)

кривая – огибающая семейства касательных (Д. П.)

Уравнение эйконала (в R²)

Пусть дано распределение коэффициента преломления среды в виде n(x, y), тогда скорость луча v(x, y) = c/n(x, y).

По принципу Ферма путь луча такой, что время ® min.

Пусть это время от точки (x₀, y₀) до точки (x, y) есть t(x, y).

Пусть угол наклона q, тогда заменим

(x, y) ®(x+Dcosq, y+Dsinq) и Dt = D/v(x, y).

Последовательно получаем:

t(x, y) = min_q(D/v(x, y) + t(x+Dcosq, y+Dsinq) + o(D)) Þ

0 = min_q(D/v(x, y) + D t¢ _x cosq + D t¢ _y sinq + o(D)) Þ

-v^-1(x, y) = min_q(t¢ _x cosq + t¢ _y sinq)

q_min: -t¢ _x sinq + t¢ _y cosq = 0 Þ tgq = t¢ _y/t¢ _x Þ

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

Выводы

Процедура решения непрерывных задач методом динамического программирования для задач, описываемых в терминах траекторных процессов приводит к уравнению в частных производных. Сравним классический подход и метод ДП по сложности численного решения в одномерном случае (N – число узлов)

Таким образом при решении серии задач метод ДП конкурирует с классическими.

Тема 6. Динамическое программирование и линейное программирование

План лекции:

Задачи оптимизации производства услуг.

Задача оптимального распределения ресурсов.

Краткое содержание лекции

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

В распределительных задачах с большим числом различных результатов производственной деятельности (i=1, …, n) и видов ресурсов (j=1, …, m) общее решение задачи оптимизации может быть при определенных условиях заменено совокупностью последовательно решаемых менее сложных частных задач оптимизации, например, по каждому из отдельных видов производственной деятельности. При этом важными понятиями ДП являются: последовательность шагов оптимизации, состояние системы распределения ресурсов и варианты решения (области изменения оптимизируемых переменных)

Рассмотрим приемы и методы динамического программирования и введенных выше понятий на примере общей задачи линейного программирования:

z =c₁x₁+c₂x₂+…+c_nx_n Þ max,

a₁₁x₁+a₁₂x₂+…+a_1nx_n £ b₁,

a₂₁x₁+a₂₂x₂+…+a_2nx_n £ b₂,

…

a_m1x₁+a_m2x₂+…+a_mnx_n £ b_m,

x₁, …, x_n³ 0,

где z - целевая функция, подлежащая максимизации;

x_i - оптимизируемые переменные;

i=1, n-номер оптимизируемой переменной;

с_i - доход от реализации i-го вида производственной деятельности;

j=1, m-номер ограничений на значения переменных;

а_ij - коэффициенты уравнений-ограничений;

b_j - величина j-го ресурса (правая часть ограничений).

Здесь каждый вид производственной деятельности i может рассматриваться как отдельный шаг оптимизации; множество возможных значений переменных (допустимая область решений) x_i как варианты решений, а количество каждого j-го вида ресурса (B_i1, …, B_ij, … B_im), 0£ b_j£ B_ij, доступного для распределения на текущем и предыдущих (либо текущем и последующих) шагах (по i-м типам деятельности) как состояние модели.

Тогда оптимальное значение целевой функции z для шагов i, i+1, …, n при заданных состояниях {B_ij} может быть записано в виде следующей рекуррентной функции Беллмана (алгоритма прямой прогонки):

f_i(B_i1, …, B_im) = max {c_i x_i + f_i-1(B_i1-a_i1x_i, …, B_im-a_imx_i)}, i =1, n; j =1, m,

где 0£ a_ijx_i£ B_ij

с начальными условиями f₀(B₀₁, …, B_0m)=0.

Оптимальное значение целевой функции в обратном времени для шагов n, …, i, i-1, …, 1 при заданных состояниях {B_ij} может быть записано в виде следующего алгоритма обратной прогонки:

f_n(B_n1, …, B_nm)=max{c_nx_n}, 0£ a_njx_n£ B_nj

f_i(B_i1, …, B_im)=max {c_ix_i+f_i+1(B_i1-a_i1x_i, …, B_im-a_imx_i)}, i=1, n; j=1, m,

где 0£ a_ijx_i£ B_ij 0£ B_ij£ b_j.

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

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

Процесс решения задачи методом динамического программирования включает два этапа. На первом этапе пошаговые задачи оптимизации приводят к условно-оптимальным по ресурсу (состояниям) решениям и одному (конечному) безусловно-оптимальному решению. На втором этапе формируется окончательная безусловно-оптимальная стратегия (x₁^опт, …, x_i^опт, …, x_n^опт) путем учета полученного на первом этапе конечного решения и затрат ресурсов на его реализацию, а также обратного по шагам анализа множества условно-оптимальных решений и выделения из него оптимальной стратегии.

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

max z = 3x_a + 2x_б, при ограничениях

x_а + 2x_б £ 6,

2x_а + x_б £ 8,

-x_а + x_б £ 1,

x_б £ 2,

[x_а], [x_б] ³ 0.

Здесь шаги оптимизации определяются порядком оптимизации различных видов услуг связи: на 1^ом шаге i=1 оптимизируется число услуг типа а, на втором шаге i=2 оптимизируется число услуг типа б.

Вектор состояния (B_i1, B_i2) определяется двумя видами ресурсов: числом телеграфных и телефонных каналов, подлежащих распределению на i-м шаге, т.е.0£ B_i1£ 6 и 0£ B_i2£ 8. Варианты решений определяются допустимой областью определения переменных x_i, подлежащих оптимизации. Далее рассмотрим реализацию алгоритма ДП для случая обратного алгоритма.

Критерий оптимальности для первого шага (оптимизация производства услуги x_б) принимает вид:

f₂(B₂₁, B₂₂)= max{2x_б},

0 £ 2x_б £ B₂₁

0 £ x_б £ B₂₂

Далее определяем допустимые границы изменения состояния на первом шаге:

Процесс поиска условно оптимальных решений по величине производства услуги x_б отобразим в форме таблицы.

Таблица

		Оптимальное решение
Xб=0	Xб=1	Xб=2	Xб=3	Xб=4	f*2(B2)	x*б
(2, 0)		-	-	-	-
(3, 2)			-	-	-
(4, 4)				-	-
(5, 6)				-	-
(6, 8)					-

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

Целевая функция на втором шаге имеет следующий вид:

f₁(B₁) = max{R₁(x_а) + f₂(B₁ – C₂(x^*_б))}, при 0£ x_а£ B₁₁; 0£ 2x_а£ B₁₂

Решение задачи второго этапа представлено таблицей.

Таблица

	R1(xа) + f2(B1 – C2(x*б))	Оптимальное решение
Xa=0	Xa=1	Xa=2	Xa=3	Xa=4	f*1(B1)	x*а
(6, 8)	0+6=6	3+4=7	6+4=10	9+2=11	12+0

Выполнение второго этапа предполагает определение на основе полученного безусловного решения x^*_а = 4 состояния ресурса, оставшегося только на первый шаг , и отыскание из таблицы 2 второго безусловного решения, удовлетворяющего этому ограничению x^*_б = 0.

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

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

f₂(B₂₁, B₂₂)= max{2x_б}

0£ 2x_б£ B₂₁

0£ x_б£ B₂₂

Так как из ограничений следует, что x_б £ min{B₂₁/2, B₂₂}, а f₂(x_б|B₂₁, B₂₂) = 2x_б, то, подставляя первое во второе и переходя к максимуму целевой функции, получим

f₂(B₂₁, B₂₂) = max{f₂(x_б|B₂₁, B₂₂) = 2x_б = 2min{B₂₁/2, B₂₂},

откуда x^*_б = min{B₂₁/2, B₂₂}.

Т.е. оптимальное значение числа услуг типа б, получаемое на шаге 1, равно минимальному из двух видов ресурсов (числу телеграфных B₂₁/2 или телефонных каналов B₂₂), распределяемых на первом шаге.

Далее для шага 2 имеем:

f₁(B₁₁, B₁₂) = max {3x_a+f₂(B₁₁-x_a, B₁₂-2x_a)}= max {3x_a+2min{(B₁₁-x_a)/2, B₁₂-2x_a}}, 0£ x_a£ B₁₁; 0£ 2x_a£ B₁₂, где B₁₁=6 и B₁₂=8 для первого шага оптимизации.

Подставляя значение ресурсов ТГ и ТФ каналов в ограничения, получим обобщенное ограничение в виде x_a £ min(B₁₁, B₁₂/2)=4. Учитывая пропорциональную зависимость значения целевой функции от значения x_a последнее неравенство превращается в равенство, а соответствующее решению x_a^*= 4 оптимальное (максимальное) ее значение f₁^*(4) = 12.

Второй этап. Безусловное решение на первом шаге оптимизации числа предоставляемых услуг x_б необходимо проводить для следующего состояния по ресурсам:

B₂₁ = B₁₁- x_a^*=6-4=2; B₂₂ = B₁₂- 2x_a^* =8-8=0.

Откуда из условного решения следует безусловное: x_б^*= min { B₂₁/2, B₂₂ }= min { 1, 0 }= 0. Таким образом, пошаговое решение задачи линейного программирования методом динамического программирования обеспечило в целом оптимальную стратегию производства услуг в одну единицу времени x_a= 4, x_б = 0 при максимальном для случая целочисленного решения значении дохода сети спутниковой связи z = 12 у.е./ед.вр.

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

,

где P_i –вероятность безотказной работы i-го интервала.

Для повышения надежности данной последовательной системы используется резервирование станций на каждом из отдельных интервалов РРЛ, вероятность безотказной работы которых определяется выражением:

,

где q_i и p_i –вероятность отказа и вероятность безотказной работы элемента на i-м интервале соответственно; x_i- число резервных станций на i-м интервале; 1+x_i –общее число (одна рабочая и x_i резервная) станций в i-м интервале. Пусть также введены ограничения на число резервных станций x_i £ 2; i=1, n. При этом суммарная стоимость резервных элементов с учетом известных ограничений на число станций в подразделении, развертывающем РРЛ не может превысить величины С_S = 1200у.ед. Остальные исходные данные, содержащие сведения о надежности p_i(x_i) и стоимости резервных средств c_i(x_i) i-го интервала, даны в таблице.

В указанной таблице значения вероятностей безотказной работы i-х интервалов при использовании в них x_i резервных станций определяются из выражения p_i(x_i) = 1-(1-p_i)^1+x_i, где p_i-вероятность безотказной работы интервала при отсутствии резерва. Необходимо определить оптимальное количество резервных станций на каждом интервале x_i^опт, обеспечивающих максимальную надежность РРЛ, т.е.

,

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

Для решения этой нелинейной задачи применим метод динамического программирования и, в частности, алгоритм обратной прогонки. При этом номер шага соответствует номеру интервала, под состоянием s_i (0£ s_i£ C_S) понимается суммарная стоимость основного и резервного оборудования, задействованного на i-м и последующих интервалах, а под вариантами решения x_i понимаем допустимое число резервных элементов в i-м интервале. Рекуррентное соотношение для функции Беллмана в данном случае может быть записано в виде:

F_n(s_n) = ,

F_i (s_i) = ,

s_i = C_S - , k = 1, i = 1, 2, …, n-1.

Из таблицы найдем границы изменения состояния s_i на каждом шаге:

s_3min=c₃(x₃=0)=300у.е.,

s_3max=С_S - = C_S - [c₁(x₁ = 0) + c₂(x₂ = 0)] = 1200 – 300 = 900у.е.,

s_2min=c₂(x₂ = 0) + c₃(x₃ = 0) = 200 + 300 = 500у.е.,

s_2max=C_S - c₁(x₁ = 0)=1200 – 100 = 1100у.е., s₁=C_S = 1200у.е.

На первом этапе обратного алгоритма рассматривают шаг 3. Найденные на этом и последующих этапах значения критерия и условно-оптимальные решения для всех допустимых значений s_i и оптимизируемых переменных x_i=0, 1, 2 представлены в нижеследующих таблицах.

Шаг 3

Шаг 2

Здесь значения функции Беллмана на втором шаге определяются с учетом ее значения на предыдущем шаге согласно выражению F_i(c_i)= p_i(x_i)F_i+1(C_S - c_i).

Шаг 3

Безусловное оптимальное решение на первом этапе получено лишь для шага 3 (числа резервных станций на первом интервале x₁^опт=1), поэтому окончательная стратегия относительно необходимых резервных станций в каждом из трех интервалов может быть получена лишь на втором этапе – этапе анализа результатов пошаговой оптимизации.

Так как x₁^опт=1, то величина ресурса, подлежащая распределению на первом и втором шагах будет равно S₂ = C_S - c(x₁) = 1200 - 200 = 1000у.е. Тогда из предпоследней таблицы для S₂=1000 получим оптимальное число резервных станций во втором интервале: x₂^опт=1. Следовательно, состояние S₃=C_S - ë c(x₁^опт) + c(x₂^опт)û = 1200 – (200 + 400) = 600у.е. Наконец, из последней таблицы для S₃=600 имеем оптимальное число резервных станций на третьем интервале: x₃^опт=1.

Окончательно, оптимальная стратегия распределения резерва имеет вид: , т.е. на каждом радиорелейном интервале должно находиться по одной резервной радиорелейной станции. Вероятность безотказной работы РРЛ при этом составляет P_РРЛ= 0, 5733; суммарные затраты и затраты по интервалам составляют C_S=1200; С₁=200; С₂=400; С₃=600у.е. Интервальные вероятности безотказной работы составляют p₁(x₁=1)=0, 91; p₂(x₂=1)=0, 84; p₃(x₃=1)=0, 75.

Тема 7. Задачи поиска оптимальных путей на графах с неотрицательными весами ребер.

План лекции:

Алгоритм Флойда – Уоршелла.

Алгоритм Дейкстры.

Краткое содержание лекции

Алгоритм Флойда - Уоршелла

Алгоритм Флойда - Уоршелла – динамический алгоритм для нахождения кратчайших расстояний между всеми вершинами взвешенного ориентированного графа. Разработан в 1962 году Робертом Флойдом и Стивеном Уоршеллом. Более строгая формулировка этой задачи следующая: есть ориентированный граф G = (V, Е) каждой дуге v –> w этого графа сопоставлена неотрицательная стоимость C[v, w]. Общая задача нахождения кратчайших путей заключается в нахождении для каждой упорядоченной пары вершин (v, w) любого пути от вершины v к вершине w, длина которого минимальна среди всех возможных путей от v к w.

Формальное описание. Для определенности положим, что вершины графа последовательно пронумерованы от 1 до n. Алгоритм Флойда использует матрицу А размера n´ n, в которой вычисляются длины кратчайших путей. Вначале A[i, j] = C[i, j] для всех i ≠ j. Если дуга i –> j отсутствует, то C[i, j] = infinity. Каждый диагональный элемент матрицы А равен 0. Над матрицей А выполняется n итераций. После k-й итерации A[i, j] содержит значение наименьшей длины путей из вершины i в вершину у, которые не проходят через вершины с номером, большим k.

Другими словами, между концевыми вершинами пути i и j могут находиться только вершины, номера которых меньше или равны k. На k-й итерации для вычисления матрицы А применяется следующая формула:

A_k[i, j] = min(A_k-1[i, j], A_k-1[i, k] + A_k-1[k, j])

Нижний индекс k обозначает значение матрицы А после k-й итерации, но это не означает, что существует n различных матриц, этот индекс используется для сокращения записи. Для вычисления A_k[i, j] проводится сравнение величины A_k-i[i, j] (т.е. стоимость пути от вершины i к вершине j без участия вершины k или другой вершины с более высоким номером) с величиной A_k-1[i, k] + A_k-1[k, j] (стоимость пути от вершины i до вершины k плюс стоимость пути от вершины k до вершины j). Если путь через вершину k дешевле, чем A_k-1[i, j], то величина A_k[i, j] изменяется.

Оценка сложности. Время выполнения этой программы, очевидно, имеет порядок O(n³), поскольку в ней практически нет ничего, кроме вложенных друг в друга трех циклов. Доказательство правильности работы этого алгоритма также очевидно и выполняется с помощью математической индукции по k, показывая, что на k-й итерации вершина k включается в путь только тогда, когда новый путь короче старого.

Пример графа и таблица расстояний для него.

Алгоритм Дейкстры.

Алгоритма Дейкстры находит оптимальные маршруты и их длину между одной конкретной вершиной (источником) и всеми остальными вершинами графа. Недостаток данного алгоритма в том, что он будет некорректно работать если граф имеет дуги отрицательного веса.

Формальное определение. Дан взвешенный ориентированный граф G (V, E) без петель и дуг отрицательного веса. Найти кратчайшие пути от некоторой вершины a графа G до всех остальных вершин этого графа.

Неформальное описание. Каждой вершине из V сопоставим метку – минимальное известное расстояние от этой вершины до a. Алгоритм работает пошагово – на каждом шаге он «посещает» одну вершину и пытается уменьшать метки. Работа алгоритма завершается, когда все вершины посещены.

Инициализация. Метка самой вершины a полагается равной 0, метки остальных вершин – бесконечности. Это отражает то, что расстояния от a до других вершин пока неизвестны. Все вершины графа помечаются как непосещённые.

Шаг алгоритма. Если все вершины посещены, алгоритм завершается. В противном случае, из ещё не посещённых вершин выбирается вершина u, имеющая минимальную метку. Мы рассматриваем всевозможные маршруты, в которых u является предпоследним пунктом. Вершины, в которые ведут рёбра из u, назовем соседями этой вершины. Для каждого соседа вершины u, кроме отмеченных как посещённые, рассмотрим новую длину пути, равную сумме значений текущей метки u и длины ребра, соединяющего u с этим соседом. Если полученное значение длины меньше значения метки соседа, заменим значение метки полученным значением длины. Рассмотрев всех соседей, пометим вершину u как посещенную и повторим шаг алгоритма.

Рассмотрим алгоритм на примере, возьмем ориентированный граф G:

Этот граф можно представить в виде матрицы С:

Возьмем в качестве источника вершину 1. Это значит что мы будем искать кратчайшие маршруты из вершины 1 в вершины 2, 3, 4 и 5. Данный алгоритм пошагово перебирает все вершины графа и назначает им метки, которые являются известным минимальным расстоянием от вершины источника до конкретной вершины. Рассмотрим этот алгоритм на примере. Присвоим 1-й вершине метку равную 0, потому как эта вершина – источник. Остальным вершинам присвоим метки равные бесконечности.

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

После того как мы рассмотрели все вершины, в которые есть прямой путь из W, вершину W отмечаем как посещённую, и выбираем из ещё не посещенных такую, которая имеет минимальное значение метки, она и будет следующей вершиной W. В данном случае это вершина 2 или 5. Если есть несколько вершин с одинаковыми метками, то не имеет значения какую из них мы выберем как W.

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

Исходя из картинки мы можем увидеть, чт