Под большими системами будем понимать системы, описываемые моделями с большим числом управляемых параметров. Если степень овражности соответст­вующих критериев оптимальности достаточно высока, то стандартные вычис­лительные средства оказываются неэффективными в силу изложенных в 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.