5.4. Методы многопараметрической оптимизации
Под большими системами будем понимать системы, описываемые моделями с большим числом управляемых параметров. Если степень овражности соответствующих критериев оптимальности достаточно высока, то стандартные вычислительные средства оказываются неэффективными в силу изложенных в 3.5.4 причин. Методы ОПС, а также ЭР-методы неприменимы, так как их вычислительные схемы содержат заполненные матрицы размерности п х и, что при больших (порядка 1000) п определяет чрезмерные требования к объему необходимой памяти компьютера.
Наиболее часто в указанной ситуации рекомендуется применять различные нематричные формы метода сопряженных градиентов (СГ). Однако далее будет показано, что в классе матричных градиентных схем (5.2) существуют более эффективные для рассматриваемых задач алгоритмы, чем методы СГ.
Пусть оптимизируемая система может быть представлена как совокупность взаимосвязанных подсистем меньшей размерности, в этом случае требования к выходным параметрам системы могут быть сформулированы в виде следующих неравенств:
где xj есть ^-мерный частный вектор управляемых параметров; хд — п^-мерный вектор управляемых параметров, влияющий на все q выходных параметров и
осуществляющий связь отдельных подсистем оптимизируемой системы. Размерность полного вектора управляемых параметров х = [х1, х2, хч] равна
Используя технику оптимизации, представленную в 2.8.2, можно привести задачу решения системы неравенств (5.54) к виду
где критерий (5.56) является сглаженным вариантом критерия минимального запаса работоспособности.
Функционалы (5.56) возникают и при других постановках задач оптимального параметрического синтеза, не основанных непосредственно на критериях минимального запаса работоспособности. Поэтому задача (5.56) имеет достаточно общий характер.
Далее будут рассмотрены методы решения задачи (5.56) при следующих предположениях.
Критерий (5.56) обладает относительно высокой степенью овражности, а его выпуклость гарантируется только в окрестности точки минимума.
Размерность (5.55) полного вектора управляемых параметров х велика, что, с одной стороны, затрудняет применение стандартных методов оптимизации из-за нехватки доступной памяти компьютера, а с другой — не позволяет реализовать предельно возможные характеристики сходимости алгоритмов.
Решение задачи анализа оптимизируемой системы требует значительных вычислительных затрат. Поэтому в процессе оптимизации требуется минимизировать количество обращений к вычислению значений/ (х).
Коэффициент заполнения у матрицы G (х) = J"(x) достаточно мал. Обычно можно полагать у ~ \/q.
Структура матрицы G(x) не зависит от точки х.
Подматрица Gy имеет размеры щ х rip а общее число ненулевых элементов равно
Таким образом, учитывая симметричность матрицы G (х), в памяти компьютера необходимо хранить
ненулевых элементов. Необходимые сведения о схемах хранения разреженных матриц содержатся, например, в [12], [23].
5.4.1. Методы с чебышевскими функциями релаксации
Обратимся снова к классу матричных градиентных схем.
Пусть Х{ (Gk) є [-rn, At], М» 777 > 0. В силу приведенных ранее предположений и сформулированных в 5.1 требований к функциям релаксации, наиболее рациональный метод должен иметь функцию релаксации, значения которой резко снижаются от R = 1 при X = 0, оставаясь малыми во всем диапазоне [0, М]. И напротив, при X < 0 функция R (X) должна интенсивно возрастать. Кроме того, отвечающая R (X) матричная функция Н должна строиться без матричных умножений для сохранения свойства разреженности матрицы Gk - J"(xk).
Покажем, что в качестве такой R (X) с точностью до множителя могут быть использованы смещенные полиномы Чебышева второго рода Ps (X), удовлетворяющие следующим соотношениям
Данное утверждение вытекает из известного факта равномерней сходимости по-Г Р (X)]
следовательности < ——-1 к нулю при s -» оо на открытом промежутке (0, 1). Да-
лее будем предполагать, что собственные числа матрицы Gk нормированы к промежутку (0, 1). Для этого достаточно вместо матрицы G рассматривать матрицу
G g -—-, а вместо, вектора g — вектор -р--.
Отвечающая принятой R (X) зависимость Н (X) имеет вид
Н(Х) = [1 - R(X)] / X = [1 - PL(X) / L] / X. (5.59)
Построение методов (5.2) непосредственно с функцией (5.59) возможно, но приводит к необходимости решения на каждом шаге по k больших линейных систем уравнений с разреженной матрицей. Далее показано, что существуют более эффективные приемы реализации.
Из (5.59) следует, что Н (X) является полиномом степени L - 2, в то время как R (X) имеет степень 1-1. Поэтому для реализации матричного градиентного метода с указанной функцией Н (X), вообще говоря, нет необходимости решать линейные системы. Метод будет выглядеть следующим образом:
хы =xk -(axE + a2Gk + ...+ aL_iG^2)gk =xk -H(Gk)gk. (5.60)
Реализация метода (5.60) может быть основана на методах вычисления коэффициентов а, для различных степеней L. При этом число L должно выбираться из условия наиболее быстрого убывания J(x). Далее обсуждается альтернативный подход, основанный на других соображениях.
Таким образом, при фиксированной квадратичной аппроксимации f(x) функционала / (л:) в окрестности х = хк мы имеем возможность переходить от Ps kPs+x за счет одного умножения матрицы Е - 2Gk на вектор 0S, в полной мере используя свойство разреженности матрицы Gk и не прибегая к дополнительным вычислениям градиента. Эффективность алгоритма (5.62) при больших значениях г) определяется множителями релаксации для малых собственных значений матрицы Gk. Рассмотрим положительную часть спектра (X > 0), что особенно важно в окрестности оптимума, где матрица G (х) положительно определена. Основное
достоинство метода с RS(X) = — состоит в том, что уже при малых s происхо-
поэтому значения производных R's (0) в последней строке таблицы характеризуют множители релаксации для отрицательных слагаемых в (5.58). Вычисление производных R's (0) может быть выполнено, исходя из следующих рекуррентных соотношений:
р;=о,р~ =-4;p;+l =2p;-4S-p;_,;i?;(o)=^-.
Значения а5, (35 для s > 8 (при X > 0) могут быть вычислены по асимптотической формуле
а5 = *^, р5 =1-а5; (5.63)
при этом Rs < 0,22.
Соотношение (5.63) получается из следующего представления полиномов Чебы-шева:
Л(Я) = т^, X = sm^,X, Се[0,1]. LsmC, 2
Действительно, при достаточно малых С, имеем [34]
дит заметное подавление слагаемых из (5.58) в широком диапазоне значений X. Далее представлены значения Rs для внутреннего максимума Rs (X) и границы диапазонов а5 < X < (3S, где |/?S(A,)| < Rs:
График функции Ф (^) представлен на рис. 5.7. Полагая х = ^ получим
|
Последняя функция и ее
числовые характеристики показаны на рис. 5.8.
Ф<£) й Ф($кр) при 4 > $кр,
где ^ = xl = 6523; Ф(^ ) = <р(хкр ) = 022. Таким образом, полагая
получим следующее утверждение: если для наименьшего (положительного) собственного числа т выполняется неравенство
Укрупненная схема алгоритма, построенного на основе соотношения (5.62), может быть реализована с помощью следующей последовательности шагов. При этом предполагается, что все переменные задачи надлежащим образом нормализованы (см. 4.3.1). Предполагается также, что переменные пронумерованы некоторым оптимальным способом, гарантирующим эффективное хранение разреженной матрицы Е - Gk в памяти компьютера.
Алгоритм RELCH.
Шаг 1. Задать начальную точку х; вычислить J := J(x); задать I, определяющее количество пересчетов по формуле (5.62) (об априорном выборе L см. далее).
Шаг 2. Вычислить g :=/(*), G :=f(x)\ принять g := g/ \\G\\, G := G/ \\G\\\ a : = 1.
Шаг 3. По формуле (5.62) построить 0L; принять хс: = х + QL.
Шаг 4. Вычислить Jt :=J(xc). Если Jt > J, перейти к шагу 5, иначе — к шагу 6.
Шаг 5. Принять a := xc: = x + a0L и перейти к шагу 4. Шаг 6. Принять х := x\J :=Jt и перейти к шагу 2.
Критерий окончания процесса здесь
не указан. Как правило, вычисления закан-
чиваются по исчерпании заданного количества вычислений функционала либо
при явной остановке алгоритма. Число пересчетов I по формуле (5.62) является
параметром, задаваемым пользователем. Согласно (5.63) первоначально целесо-
образно полагать L = / =
1,3 Jr\, где ц — оценка степени овражности миними-
зируемого функционала. При таком выборе L множители релаксации в положи-
тельной части спектра будут гарантированно меньше 0,23. При конструировании
алгоритмических способов задания L необходимо учитывать, что последователь-
ность {Js}, где Js = J(xk + Qs) не будет при 5 -> оо убывать монотонно. На
шаге 5
алгоритма применена регулировка нормы вектора продвижения с целью предотвращения выхода из области справедливости локальной квадратичной модели функционала.
5.4.2. Характеристики сходимости и сравнение с методами сопряженных градиентов
Дадим оценку эффективности метода (5.62) по сравнению с методами сопряженных градиентов (СГ-методами). Для задач большой размерности (когда число итераций меньше размерности) можно гарантировать сходимость СГ-мето-дов только со скоростью геометрической прогрессии даже для сильно выпуклых квадратичных функционалов [48].
Действительно, рассмотрим случай
и оценим скорость сходимости метода СГ к экстремальной точке х = 0. Итерация я* полученная методом СГ, может быть представлена в виде [65]
где Pk(G) — матричный полином k-й степени. При этом из свойств метода СГ следует, что коэффициенты ct, ck полинома Pk(G) на каждой итерации принимают такие значения, чтобы минимизировать величину /(xk), только множителем отличающуюся от функции ошибки. Иначе говоря, k-e приближение минимизирует f(xk) среди векторов х° + V, где вектор V является элементом подпространства, натянутого на векторы Gx°, G2x°, Gkx°. Полагая
Выберем в качестве полинома Pk (X) близкий к оптимальному полином, наименее уклоняющийся от нуля на промежутке [ти, М], содержащем все собственные
значения положительно определенной матрицы G и нормированный так, что Р4(0) = 1.
Линейной заменой переменных
задача сводится к построению полинома наименее отклоняющегося от нуля на
г л
41
4. М + т ,
промежутке t є [-1, 1] и принимающего в точке
t0 ,
(соответствующей
М-т
X = 0) значение 1. Решение последней задачи дается полиномом [7]
f (0_ ад _ rt(t) ^
cos (k arccos t0) Tk (t0)' где Tk (t) = cos (k arccos t) — полином Чебышева. При этом
тахІ7;(ґ)|=1—-—rmax|71(*)l
Очевидно,
Так как справедливо представление
где q = 1 —р. ТакимЧ)бразом, сходимость метода СГ со скоростью геометрической прогрессии доказана. Точное значение Ц, справедливое для любых к, будет при этом равно
Из (5.68) следует, что при г] » 1 сходимость может быть очень медленной.
«Конечность» метода СГ, то есть точное решение задачи минимизации квадратичной функции за п шагов, где п — размерность пространства поиска, проявляется только при достаточно большом количестве итераций. При этом степень полинома Pk (X) в (5.65) будет равна п, и оптимальный выбор этого полинома сводится к локализации его п корней в точках Хи Х2, Хп, что приведет к точному решению задачи (f(xn) = 0).
Легко видеть, что оценка (5.68) для квадратичной функции общего вида
преобразуется к виду
где х* — оптимальная точка, не совпадающая в общем случае с началом координат.
Важная особенность алгоритмов типа RELCH заключается в том, что соответствующие множители релаксации будут определяться только числом итераций L и степенью г] обусловленности задачи независимо от размерности п. В то же время в схемах методов СГ для завершения каждого цикла спуска требуется порядка п итераций; в противном случае, согласно последней формуле, скорость сходимости может быть очень малой. Кроме того, каждая итерация метода С Г даже для квадратичного случая требует нового вычисления градиента, то есть дополнительных вычислительных затрат по анализам функционирования оптимизируемой системы.
Будем далее полагать, что алгоритм RELCH реализован с постоянным L = 1,3^п, имея в области X > 0 множители релаксации, не превышающие значения 0,23.
Рассмотрим задачу минимизации квадратичного функционала f (х) = - (Gx, х) с
положительно определенной матрицей G. Оценим количество вычислений / (х), требуемое для достижения контрольного вектора х! с нормой \ х! || < 0,23 методом СГ и алгоритмом RELCH из начальной точки хР с || хР || = 1. По достижении точки xf вся ситуация повторяется, поэтому полученные далее сравнительные оценки эффективности имеют достаточно общий характер.
Будем предполагать, что для вычисления производных применяются двусторонние конечноразностные соотношения, что в нижеследующем анализе дает дополнительные преимущества методу СГ.
Для достижения вектора xf алгоритму RELCH требуется вычислить в точке хР слабо заполненную матрицу Гессе и вектор градиента /' (хР). При коэффициенте заполнения у это потребует около 2уп2 вычислений /. Далее выполняется L = ХЗ-у/ц итераций по формуле (5.62), не требующих дополнительных вычислений целевого функционала /.
Чтобы получить вектор xf, методу СГ потребуется N итераций, где число N определяется из условия:
то есть N= . Для выполнения каждой итерации необходимо обновление век-
Inq
тора градиента, что связано с 2п вычислениями f(x). Общее число вычислений /
равно ^п Относительный выигрыш в количестве вычислений / методом Inq
RELCH по сравнению с методом СГ задается функцией ¥ (г|) = -2,2/ (уп In q). Очевидно, при г| -> оо имеем q (г|) -> 1 и ¥ (г|) -> оо. Характерные значения ¥ для у = 0,01 и п = 1000 даны ниже:
Таким образом, для получения сравнимых результатов при г] = 104 алгоритму RELCH потребуется приблизительно в 11 раз меньше вычислений /, чем методу СГ. Следует, однако, учитывать, что при увеличении г| возрастает количество L пересчетов по формуле (5.62). Это может приводить к возрастанию влияния вычислительных погрешностей при вычислении 9S с большими номерами S.
Пример. Рассмотрим модельную задачу минимизации квадратичного функционала f(x) с п = 200, г| = 1500, у = 0,025. Для определенности примем, что время однократного вычисления / (х) эквивалентно выполнению 102гг операций умножения с плавающей точкой. Время выполнения одной операции умножения для определенности и чисто условно положим равным ty = 3-Ю"5 с. Вычисление значения/^) занимает при этом t/= 0,6 с процессорного времени. Для вычисления /' и /" с помощью общих конечноразностных формул потребуется, соответствен
но, t' = 2nt/=A мин, if' = 2yn2t = 20 мин. Число пересчетов по формуле (5.62) равно L = i,3yfr\ = 50. При каждом пересчете производится умножение слабо заполненной матрицы Е - 2Gk на вектор 0S, что требует yn2ty = 3-Ю"2 с машинного времени. Время построения вектора 0^ без учета вычисления /', /" составит около 50-3-10"2 = 1,5 с и может в расчет не приниматься.
В результате получается, что для построения контрольного вектора xf с || xf || < 0,23 методом RELCH потребуется около £" = 20 мин машинного времени. Метод СГ затратит, соответственно, Ч* (1500) • 20 = 1,3 ч.
При повторном применении алгоритма RELCH .к построенному вектору xf мы получим вектор xf' с || xf' Н < 0,231| xf\, и т. д. Следовательно, если обозначить соответствующую последовательность векторов через {хт}, то норма вектора х будет убывать по закону геометрической прогрессии || хт \\ < dm \ х° ||, где d < 0,23 независимо от величины Г] и п.
Важным дополнительным преимуществом алгоритма RELCH по сравнению с методом СГ является его достаточно высокая эффективность в невыпуклом случае, так как функция релаксации метода в левой полуплоскости целиком расположена в разрешенной области и множители релаксации для X < 0 быстро растут по абсолютной величине при переходе от 0S к 0S+1. Характеристики роста были приведены ранее.
л
Рис. 5.9. Области работоспособности: 1—е^ = £^\2 — ем = £^>е^
Так же как и в случае ЭР-методов, можно показать, что эффективность рассматриваемого подхода сохраняется при степенях рвражности, удовлетворяющих I
неравенству г| < . Области работоспособности алгоритмов RELEX, RELCH
в плоскости (я, г|) представлены на рис. 5.9. Ясно, однако, что при умеренных размерностях п более эффективными, вообще говоря, оказываются алгоритмы
|
типа RELEX. Они позволяют
за меньшее число Ny операций умножения матрицы на вектор получить заданные
значения множителей релаксации. При больших г| это приводит к существенному
уменьшению накопленной вычислительной погрешности. Для подтверждения данного
замечания достаточно проанализировать характер изменения множителей релаксации
при применении формул пересчета (5.39) и (5.62). Характерные зависимости для
рассмотренных случаев (для фиксированного Х> < 0) представлены на рис.
5.10.
Видно, что если область локальной квадратичности функционала J(x) невелика (C,k мало), то необходимые значения = 1 и более эффективными могут оказаться методы типа RELCH.