CC BY
50
УЧЕНЫЕ ЗАПИСКИ Ц АГ И Том XXI 1990
№ 3
УДК 629.7.015.017.2
629.7.015.4:533.6.013.43
СРЕДНЕКВАДРАТИЧНАЯ АППРОКСИМАЦИЯ ЧЕТНЫМИ НЕОТРИЦАТЕЛЬНЫМИ ДРОБНО-РАЦИОНАЛЬНЫМИ ФУНКЦИЯМИ
С. Н. Супруненко
Рассматривается задача среднеквадратичной аппроксимации дробно-рациональными функциями, числитель и знаменатель которых являются четными неотрицательными вещественными полиномами. Такие задачи возникают, в частности, при построении моделей спектральных плотностей случайных процессов и амплитудно-частотных характеристик динамических систем. Предлагается метод решения таких задач, основанный на сведении их к определению минимально-фазового спектрального фактора аппроксимирующей функции.
1. Введение. Четные неотрицательные дробно-рациональные функции наиболее широко используются в теории случайных процессов и в теории автоматического управления, где они отождествляются соответственно с рациональной спектральной плотностью стационарного случайного процесса и с квадратом амплитудно-частотной характеристики линейной стационарной динамической системы. При исследовании систем управления летательных аппаратов может возникать потребность построения таких зависимостей по имеющимся теоретическим или экспериментальным данным. В связи с этим практический интерес представляет решение следующей аппроксимационной задачи. Для вещественной функции (ш), заданной на интервале вещественного аргумента шє[0,1], требуется определить дробно-рациональную функцию наилучшего приближения
5 («а, Ьх, 6Г) = П ^°*2’ = ежі + дх2 «>* + - + вхт о>2т~2
П(ш2, ву) ву і -|- буа ш*+ 0ул <о2”-2
исходя из минимизации среднеквадратичной ошибки аппроксимации
1
» (6х, Оу) = | [/(<») - 5(ш», вх, ву)рау (со) (1)
о
по коэффициентам полиномов числителя и знаменателя при условии неотрицательности этих полиномов:
П((й2, 0*)^О и П(ш», 0у)^;О ДЛЯ всех 10. (2)
В выражении (1) функция ф(со)—ограниченная и неубывающая на рассматриваемом интервале. Она задает распределение весов аппроксимации. Поскольку коэффициенты дробно-рациональной функции могут задаваться с точностью до общего множителя,
в постановке задачи используется также дополнительное условие однозначности решения, формулируемое в виде
ftr П (">*» 9лг) + М-у П (со2 , 0у) = 1, (3)
где ч>х, «у, (Ау^О — заданные параметры.
Отметим, что принятое определение среднеквадратичной ошибки аппроксимации через интеграл Стильтьеса дает возможность единообразно учитывать случаи формирования этой ошибки как обычным интегрированием (ср(со) •—непрерывна и дифференцируема), так и конечным суммированием (ф(со)—«ступенчатая» функция). Единичный интервал аппроксимации принят с целью упрощения записи и обеспечения благоприятного режима вычислений при решении задачи на ЭВМ. Получить этот интервал можно соответствующим масштабированием исходного аргумента.
Вследствии нелинейности аппроксимирующей функции по коэффициентам знаменателя и необходимости соблюдения условий (2) — (3) решение рассматриваемой задачи требует привлечения довольно сложных методов нелинейного программирования
[1]. Для условной минимизации ошибки аппроксимации принципиально возможно использование универсальных алгоритмов и соответствующего программного обеспечения. Однако из-за сложности преобразования условия неотрицательности (2) в эквивалентные ограничения на искомые коэффициенты (простые соотношения существуют только для т, я<3, а для т, п>5 явные формульные зависимости вообще отсутствуют [2]), этот подход оказывается неэффективным или даже проблематичным. Очевидно, в рассматриваемом случае необходим специализированный метод решения, учитывающий особую структуру аппроксимирующей функции. Разработке такого метода и посвящена данная статья.
Развиваемый в статье подход основывается на двух идеях, предложенных соответственно в [3, 4] и в [5]. Первая из них предусматривает замену сложной процедуры непосредственной минимизации по коэффициентам знаменателя дробно-рациональной аппроксимирующей функции на более простую итерационную процедуру уточнения искомых коэффициентов методом последовательных приближений. Однако, использование только этой идеи еще не гарантирует решения поставленной задачи. Численные расчеты для ряда примеров из практики показывают, что нарушение условия (2) в процессе итераций — явление довольно частое, и оно ведет к прекращению сходимости итерационного процесса. В связи с этим уместно указать на ненадежность предложенного в [4] метода аппроксимации амплитудно-частотных характеристик, в котором условие неотрицательности аппроксимирующей функции не учитывается.
Чтобы учесть условие (2), предлагается использовать идею [5] перехода от условной минимизации к безусловной путем нелинейного преобразования переменных минимизируемой функции на основе соотношений спектральной факторизации полиномов [6]:
П(и2, 0*) = Я(г<о, a)P(-ta), а), П (ш*, 0у) = Р(ш, b)P(—iu,b).
В этих разложениях полиномы
P(s, а) = а1 + я2 s + ат sm~l , P(s, b) = bx + b2 s + ... + bnsn~l
имеют вещественные коэффициенты и не имеют корней в правой полуплоскости комплексного аргумента s.
Для минимизации ошибки аппроксимации по коэффициентам полиномов P(s, а) и P(s, Ь)—■ могут быть построены эффективные итерационные схемы решения с применением быстродействующих ньютоновских алгоритмов безусловной минимизации. Предлагаемый метод решения ведет к получению минимально-фазового фактора U7(s, a, b)=P(s, a)/P(s, Ь) в спектральном разложении аппроксимирующей функции S(ю2, 0х, 0y) = W'(i(B, a, b) W(—/о), а, Ь), что весьма важно для многих практических приложений.
2. Минимизация ошибки аппроксимации. По аналогии с [3, 4] итерации уточнения коэффициентов аппроксимирующей функции методом последовательных приближений будем строить путем минимизации ошибки аппроксимации очередного приближения, формируемой в виде
1
5?=6(6?> ey) = j(n(“2- е?)/Н-П(“2. еЯ)2^9<“)- (4)
о
где d.i/q (ш) = d<f (<*)/n (u>2, 0^-1)2, q — номер итерации. Зависимость 8 0у)—
квадратичная и в стационарной точке итераций (т. е. при 0® = 0^-1) совпадает
9—Ученые записки № 3
115
с исходной функцией ошибки аппроксимации (1). Если на каждой итерации ц обеспечивается абсолютный минимум ошибки 69, и если последовательность этих минимумов оказывается монотонно убывающей, то пределом получаемой последовательности должно быть стационарное решение, совпадающее с локальным минимумом исходной функции ошибки аппроксимации е (0^., 0у).
Следует иметь в виду, что возможность существования множества локальных минимумов в рассматриваемой задаче является естественным следствием нелинейности аппроксимирующей функции по коэффициентам своего знаменателя. Сходимость итераций к конкретной точке локального минимума, равно как и сходимость вообще, зависит от выбора стартовой точки по этим коэффициентам. В частности, во многих практических случаях для задания нулевого приближения оказывается достаточным принять П(шз, 0®^ = 1. при таком задании стартовой точки на первой итерации будет получено решение, совпадающее с точным решением исходной аппроксимацион-ной задачи, но только с дополнительной весовой функцией, равной квадрату знаменателя искомой аппроксимирующей дробно-рациональной функции. Как правило, эта первая итерация дает основной вклад в аппроксимацию, а последующие уточняют решение. Разумеется, указанный способ задания стартовой точки не универсален и вопрос о ее выборе должен решаться с учетом физического содержания аппроксима-ционной задачи.
Необходимым условием сходимости метода является обеспечение абсолютного минимума ошибки аппроксимации очередного приближения. Так как минимизация 5, должна выполняться с соблюдением условий (2) — (3), требуется прежде всего выяснить существование и единственность нужного решения. Предварительно сделаем
некоторые преобразования. Используя соотношения П (<о2, = а? д (о>2, 0у) =
= Оу 0у, В которых Их = [1, (Й2.......ш2т—2]Т И 2у = [1, 0)2, ... , й)2п_2]т — векторы
размерности тип соответственно, преобразуем (4) к следующему виду (для упрощения дальнейших записей индекс <7 опущен):
5 (0„ 0у) = 0* Нхх 6Х - 20у Нух дх + в; Яуу 0у, (4а)
1 1 1 где Нхх = 2Х ЯТХ СО/ (Ш), Нух = J йу а;/(О,) аГф («), Нуу = § ^y Щ / (ш)1 Щ (<о) 0 0 о
Матрицы Нхх, Нух, Нуу (размерностью тХт, пХт и пХп соответственно) имеют ганкелеву структуру. Элементы матриц образованы последовательностями степенных четных моментов. Например:
1
луж(/, Л =- = J Ш2 (‘+/-2)/(“)^(“), (I = 1, , п; у = 1, ... , т).
О
[Ьх1
Объединение 0^ и 0^ в единый вектор вге= позволяет на основании (4а)
записать У
8(0г)=0^Ягг0г> (46)
[г_г _//т 1
*х У* —матрица блочно-ганкелевой структуры.
“уЛГ Нуу J
Из анализа знаковой определенности функции б (02) следует, что если на рассматриваемом интервале изменения аргумента со функция ф(со) имеет бесконечное число точек роста, то матрица Нгг положительно определена. В противном случае, гарантируется только неотрицательная определенность этой матрицы. Конечное число точек роста функции 'ф(ю) соответствует формированию ошибки аппроксимации в виде конечной суммы. Во избежание возможных затруднений для таких случаев будем полагать, что число слагаемых в сумме достаточно велико и, таким образом, строгая положительная определенность матрицы Нгг обеспечивается.
С использованием объединенного вектора перепишем в более компактной форме также и условие (3): /^0г = 1. В этой записи т+п—верный вектор
.Р* определяется в соответствии с формулами:
ч:;] ’ рхг=^х’ лг®*"-2]1* /гу=к> !*у®2г-. ^
В m-мерном евклидовом пространстве Ет с ортонормированным базисом, в котором задается вектор 0*, выделим область Qx, во внутренних точках которой полином П (w2, 0*) положителен, а в граничных точках может принимать и нулевые значения. Нетрудно показать, что эта область является выпуклым замкнутым конусом. Аналогичную область QvczEn выделим и для полинома П(ш2, 0У). В объединенном евклидовом пространстве Ет+п, образованном суммой ортогональных подпространств Ет и Еп, области Qx и Qy являются ортогональными проекциями некоторого выпуклого замкнутого конуса, Qz, во внутренних точках которого вещественная четная дробнорациональная функция положительна, а в граничных точках может также иметь вещественные нули и полюсы. Геометрическим образом условия (3) в Ет+п является некоторая гиперплоскость Rz. Сечение конуса Qz гиперплоскостью Rz дает выпуклое замкнутое множество Tz = Qzf\Rz. Так как функция б(02) в силу своей положительной определенности выпукла всюду в Ет+п, то на множестве Тг она имеет единственный локальный минимум, который является также и минимумом абсолютным [1].
3. Переход к задаче безусловной минимизации. При минимизации функции 8(0Z) на Tz наиболее просто учесть условие (3). Воспользуемся методом множителей Лагранжа и с учетом (4 б) составим функцию
Ф (вг) = ^Я22вг+ А/=-гт0г> где X — множитель Лагранжа.
Как и б(02), эта функция выпукла всюду в Ет+п и поэтому на выпуклом множестве Qz имеет единственный минимум, зависящий от Я. Очевидно, выбором X можно обеспечить принадлежность точки минимума также и множеству Rz.
Если игнорировать условие (2), то минимизируя Ф(0г) на открытом множестве и затем определяя К таким образом, чтобы выполнялось равенство F\ 0z=1, нетрудно получить следующее решение:
6* = arg min 5 (0г) = 8* HJ*FZ, z
8, = 8(0г*) = (^/у-1^)-1.
В геометрической интерпретации это решение определяется как точка касания гиперплоскости Rz с эллипсоидом минимального уровня 0* Hzz Qz = й*. Если эта точка принадлежит конусу Qz, то найденное решение будет искомым. Однако в общем случае 0г £ Qz и поэтому минимизировать функцию Ф(02) приходится с учетом
(2). Рассмотрим способ учета этого ограничения путем нелинейного преобразования аргумента Qz.
В соответствии с [5] введем для полиномов P(s, а) и P(s, b) вспомогательные полиномы P(s, х) и P(s, у), коэффициенты которых связаны соотношениями
xk = (— l)[kl2] ak, k=\,...,m\
У* = (-1 )[kmbk, * = 1......я.
(Квадратные скобки в операциях над индексами означают целую часть). Коэффициенты числителя и знаменателя функции 5(ш2, 0*, 0„) определяются через коэффициенты вспомогательных полиномов на основании соотношений
0,= Г(*)тлг, 0у=Г(у)ту, (5)
в которых Г (х) и Г (у) — так называемые расширенные матрицы Гурвица (РМГ) векторов хну. Элементы матрицы Г (х) определяются по правилу
ПРИ 1:S2/— гг£/и 1 7(*. -/) = ] i \ I, J= 1. ... . т.
{ 0 в остальных случаях J
Аналогичную структуру имеет матрица Г (у). Если ввести объединенный т+п-мерный вектор г = [ у ] ’ т0 два соотношеиия (5) объединяются в одно:
0г = Г BD(z)T г. (5а)
ГГ(лт) О
Здесь матрица ГBD (г) = ^
О Г(>).
— блочно-диагональная.
Нелинейное преобразование (5а) отображает всю открытую область определения вещественного вектора г в замкнутую область <2*. Но это преобразование свойством взаимно-однозначного соответствия не обладает, поскольку разным г^Ят+п может соответствовать одна и та же точка ЭгеЕЕОг. Используя (5а) как замену переменных В Ф(0г), получим
ф (о, (г)) = гт Тво {г) Нгг Гво (г)тг + Щ Гво (г)т г.
Теперь можно минимизировать функцию Ф(0г(г:)) на открытом множестве. Необходимо только учесть, что из-за нелинейности преобразования переменных эта функция может иметь множество локальных минимумов, не все из которых отображаются в точку абсолютного минимума исходной функции Ф(02) на множестве <32. Выясним, каким образом при безусловной минимизации функции Ф(02(г)) гарантировать получение абсолютного минимума.
В евклидовом пространстве Ет, в ортонормированием базисе которого задается вектор х, существует связная замкнутая область ГРХ, во внутренних точках которой корни полинома Р(э, а(х)) принадлежат левой полуплоскости комплексного аргумента 5, а в граничных точках среди этих корней имеются и чисто мнимые [7]. Эта область коническая, причем, если х £ Л® то и — х £ О®. Для определенности за О® примем только одну сторону конуса, соответствующую значениям а, >0. Аналогичную область й°;с£" выделим и для полинома р(я, Ъ(у\). Между точками областей О® и <2х, равно как и между точками областей О® и (2У, существует взаимно-однозначное соответствие [5]. В объединенном евклидовом пространстве Ет+п области /Э^. и Оу следует рассматривать как ортогональные проекции замкнутой связной конической области в точках которой спектральный фактор №(«, а(х) Ь(у)) является минимально фазовым. Поскольку между парами проекций (х £ у £ и
йу^Яу) имеется взаимно-однозначное соответствие, то такое же соответствие должно существовать и непосредственно между точками ге и 02е<2г. Но в таком случае в области Э® присутствует только одна точка г*, являющаяся прообразом точки локального абсолютного минимума функции Ф(0г) на множестве <32- Учитывая теперь непрерывность и дифференцируемость функции Ф(0*(2)). можно также заключить, что эта точка г*е /Э® должна быть и единственной точкой локального абсолютного минимума функции ф(0*(,г)) на множестве И®.
Таким образом, можно представить следующий способ решения. Необходимо искать локальный минимум функции Ф(0*(г)) каким-либо методом безусловной минимизации и, если найденная точка принадлежит Д® , принять ее за искомую. В противном случае перенести решение в О® и продолжить поиск локального минимума. Перенос решения в область выполняется путем спектральной факторизации аппроксимирующей функции с выделением ее минимально-фазового фактора. Поскольку при каждой такой операции переноса найденная точка 0г и соответствующее значение Ф(0г) не изменяются, то в конечном итоге должен быть достигнут локальный минимум функции Ф(0г(г)), обязательно принадлежащий области О0г.
Обратимся теперь непосредственно к минимизации Ф(02, г)). При выводе необходимых соотношений понадобится обозначение следующей блочно-диагональной матрицы:
Ово(г) о 1
вок [ 0 й (у)
По принятой в [5] терминологии й(х) и в (у)—четно-разреженные ганкелевы матрицы (ЧРГМ) векторов х и у. Элементы матрицы и (х) определяются правилом:
. . [хал-1 при I + /=четном )
«(*./ = (г+')/2 ... /, /= 1,.,.,/п.
I. 0 при I 4-у = нечетном )
Аналогичную структуру имеет и матрица О (у).
Введенные блочно-диагональные матрицы обладают теми же свойствами, что и матрицы РМГ-ЧРГМ {5]. Используя эти свойства, нетрудно разложить зависимость Ф(0г (г)) в ряд Тейлора в окрестности точки г. В результате оказывается возмож-
ным получить явные выражения для градиента и матрицы вторых производных минимизируемой функции:
У2 Ф (9* (*)) = 4Гво (г) £ (г), |
у2г Ф (бг (г)) = 40вд (,§■ (г)) + 8Гао (г) Нгг ГВ£) (г)т, |
где g(z) = НггГво(г)т г + I Рг = Ф (0г).
Точка локального минимума функции Ф(02(г)) очевидно должна удовлетворять уравнению
Гво(2)г(*)“°. (7)
которое получается приравниванием градиента нулю. Структура этого уравнения предполагает возможность существования решений двух типов. Решения первого типа удовлетворяют соотношению £(2)=0. Поскольку в данном случае у0 Ф(02)=О,
то получение решения эквивалентно минимизации функции Ф(0*) на открытом множестве. Соответствующие формулы расчета 0* уже были приведены в начале раздела. Может оказаться, что 6* £ фг и значит в %™.+п отсутствует точка г*,
обращающая вектор £(г) в ноль. Тогда необходимо искать второй тип решения, отвечающий условию g(z)фQ. В соответствии с требованием разрешимости уравнения (7) относительно нулевого вектора §(г) этот второй тип решения должен удовлетворять соотношению сЫ Гвп(г) =с1е1 Г (л:) с1е1 Г (г/) г 0. Отсюда можно заключить, что решение принадлежит границе области £)° (в граничных точках областей и
матрицы Г (я) и Г (у)—вырождаются [5]). В качестве иллюстрации такого решения на рис. 1 отображена ситуация, когда пограничное решение возникает по коэффициентам числителя аппроксимирующей функции (<1е(: Г(;с) =0).
Минимизация функции Ф(0г(г)) приводит к решению, зависящему от неизвестного множителя Лагранжа. Учитывая, что в точке минимума одновременно с выполнением (7) должно также выполняться равенство Гво (г*)т г* = 1, находим:
X = - 2г*т Гво (г*) Нгг Гво (2*)т г* = - 28*.
Полученное соотношение связывает X со значением ошибки аппроксимации б в точке минимума. Поскольку эта точка заранее неизвестна, то определить X до начала процесса минимизации не удается. Чтобы обойти данное затруднение, исключим X из процесса минимизации. С этой целью введем нормированные переменные 2 = 8^“1/2 г, 0г = 8Г1 в2 и вместо функции Ф(0г(г)) будем минимизировать функцию
Ф (г) = 6Г2 Ф (в, (ъУ2 г)) = !тГво (2) Нгг Гво (г)1 г — Тво (г)тг. (8)
Функция Ф(г) не зависит от Я и поэтому минимизация ее не требует знания этого параметра. Необходимые для определения точки локального минимума ге£т+" формулы расчета градиента и матрицы вторых производных функций Ф(г) точно
Рис. 1
такие же, как и у функции Ф(02(г)). При использовании соответствующих формул (8) — (9) следует только переобозначить г на г, а вместо вектора g(z) принять вектор g (г) = Нгг Гво (г)т г—Рг. При поиске абсолютного минимума функции Ф (г) следует иметь в виду, что масштабирование г не меняет статуса принадлежности г области £>“ и поэтому условия : ( и г ( /^ — эквивалентны.
После того, как найдена точка минимума г* функции Ф (г), значение 5* вычисляется по формуле
5* = (г’т Г во (*•) Нгг Гво (^)т ,
а точка минимума г* функции Ф(вг(г)) определяется как г* = 8У2г*.
4. Алгоритм решения. Рассмотрим вначале задачу минимизации ошибки аппроксимации очередного приближения. Наличие явных выражений для градиента и матрицы вторых производных минимизируемой функции дает возможность применить метод Ньютона, обладающий свойством квадратичной сходимости [1].
Пусть гк — известное приближение к точке локального минимума функции Ф (г) на шаге итерации Л. Решение на шаге к + 1 задается в виде гк + 1— = г*-|-аАДг*, где вектор Дг* и параметр ак задают соответственно направление и длину шага итерации. С учетом (6) для определения Дг* можно получить следующее уравнение:
J (Р) аР => - Гво (Р) $ (Р) , (9)
где У (Р) = 2Гво (Р) Игг Гво (Р)т + аво (г*)) ,
£ £к) « Нгг Гво (Ру Р-Р2.
Чтобы вектор Агк был направлением убывания функции Ф(г), требуется положительная определенность матрицы J(zh). В процессе итераций это требование может нарушаться, даже если итерации начинались из области положительной определенности. В связи с этим при решении (9) следует предусмотреть возможность модификации матрицы ./(г*), позволяющей восстановить свойство положительной определенности матрицы при его утрате. Из существующих различных способов для этой цели можно рекомендовать модификацию по методу Гилла и Мюррея [8,] предусматривающую коррекцию диагональных элементов системной матрицы при разложении ее на треугольные множители (модифицированная ЬЬТ или 1.йи факторизация) в процессе
решения (9). С целью обеспечения гарантированной сходимости итераций к точке локального минимума, особенно в условиях возможных модификаций матрицы J(zk). необходимо оптимизировать длину шага. Подставляя выражение для гк+1 в формулу (8) и учитывая (9), получим приращение функции Ф(г) на шаге к в зависимости от параметра длины шага ак (индекс к для упрощения записи формул далее опускается):
ДФ (а) = — 4а + 2А2 а2 4<43 а3 -г А4 а4,
где
Аг = Дгт I (г) Аг, А2 «= Дгт J (г) Дг ,
43 = ? Гво (Аг) Нгг Тво (Дг)т Дг, А4 = дР Гво (Д5) Н2г Гво (ДТу Д5 .
Через У(г) в соотношении для коэффициента А1 обозначена фактически используемая модификация матрицы J(z). Так как J{z) должна быть положительно определенной, то имеет место Л4>0. Кроме того, при Дг^=0 выполняется А',>0. Отсюда нетрудно заключить, что функция ДФ(“) имеет один или два минимума в области положительных а. Оптимальные значения а находятся из условия выполнения дДФ(“)/<5а=0, т. е. как корни уравнения
Л2 а + З<43 а3 А4а3 — Ах.
При наличии двух положительных корней следует взять наименьший.
Ускорить сходимость итераций можно также оптимальным масштабированием вектора г, а еще лучше — раздельным масштабированием векторов х и у, являющихся ортогональными составляющими вектора г. При раздельном масштабировании статус принадлежности г решения области не нарушается. Масштабные
множители введем преобразованием х : = ъЦ2 х, у : — т|/2 у. Используя формулу
(8), нетрудно получить приращение минимизируемой функции в зависимости от масштабов т* и т„:
дф ту) — —2В1тх — 2В2ту + В3т2 — 2В^хху + Вьх2у , где Вг = Р].ёх> °> в2 = Р],Ъу> О,
В3 — 0* Нхх Ьх, В4 *= 0у НуХ 6Х, Вс, = 0у Нуу 0у.
Абсолютный минимум функции ДФ(Тх, ху) достигается при
■гх = (б! В, + Я2 ВЖВг Въ - В\) > О,
1у = (В% В3 + В1 В4)/В3 В& — В2) 0.
С этими масштабами значение ошибки аппроксимации очередного приближения равно:
8 = 82 т* + В2 т*) .
ПР“ машинной реализации рассмотренной выше схемы минимизации функции Ф(г) необходимо учесть, что в случае малых значений 5* возможно переполнение арифметического устройства ЭВМ. по мере приближения итераций к точке решения (если ||г*||~1, то ||г*||~ 8~1/2). В связи с этим целесообразно ввести дополнительное масштабирование переменных, позволяющее оперировать с приемлемыми для ЭВМ числами. Можно предложить следующий способ такого масштабирования. На каждом шаге итераций после вычисления оптимальных масштабов т* и ту выполняется
переопределение гк : = р^/2 гк, 0* : = р* 6*, Рг : = р* Рг. Параметр дополнительного масштабирования р* вычисляется по формуле р*=8*/8А_в которой Ьк_1 и 8Й - оценки ошибки аппроксимации 6 предыдущего и текущего шагов итераций. Для вычисления 8* рекомендуется формула: 8* = Ч- 1/'К р1 вх + х
в*)- Выражение в скобках в этой формуле вычисляется до выполнения
дополнительного масштабирования. Если перед началом итераций принять 50 = 1, то по достижению стационарной точки решения окажется (если это точка минимума) 5* : = 8*, а место векторов г* и 0* займут вектора г* и 0*. Исходный вектор Р2 будет модифицирован значением 8*. Отметим, что рассмотренное дополнительное масштабирование не влияет на сходимость итераций.
На основании изложенного можно рекомендовать следующий алгоритм минимизации ошибки аппроксимации очередного приближения.
1. Подготовка. Задать стартовую точку 2° £ (если она не задана), принять г*: =2°, 8А_1: = 1, /: = 0.
2. Циклы переноса решений в /3° . Проверить условие />/доп и в случае его выполнения уйти на 3.
2.1. Ньютоновские итерации по минимизации 6. Установить к: =1.
2.1.1. Определить оптимальные масштабы гх, ъу и затем вычислить 8к, рь. Промасштабировать векторы 2й, 6*, Рг и принять 6ь_1: =6*.
2.1.2. При й>2 проверить условия окончания итераций: если р&>1 или 6»6ДОп, то уйти на 2.2.
2.1.3. Определить элементы матрицы J (г*) (нужен только нижний
или верхний треугольник матрицы) и вектора Гво (2*) (2й). Выполнить
разложение (по необходимости — модифицированное) матрицы на треугольные множители и определить Д2* из (9).
2.1.4. Вычислить оптимальную длину шага аи и определить новое решение гк : = гк + а^ Агк и соответствующий ему вектор 0* . Принять к : =&+1 и уйти на 2.1.1.
2.2. Проверка условия гк £ /?® . Используя критерий устойчивости Рауса, проверить выполнение хк £ Ойх и ук £ Оу. При одновременном выполнении — уйти на 3.
2.3. Спектральная факторизация. Перенести решение в С*® путем спектральной факторизации полиномов П 0*) и/или П (<о2, 0*^
(в соответствии с нарушением условия хк £ £>° и/или ук £ £)“) . Если для факторизации используется итерационный алгоритм, то в качестве стартовой точки итераций взять го и найденное решение опять запомнить в 2°. Принять I : =/+1 и уйти на 2.
3. Окончание решения. Выполнить преобразование ак = а (хк),
Ьк = Ь{ук) и выдать ай, 6*, 0*, 0у, 8* :=В* в качестве решения задачи.
Отметим, что при отсутствии априорных данных о стартовой точке г° £ й®’ для ее задания могут быть использованы биномиальные коэффициенты (я® = ,
и Ь°1 = . *=1..я), являющиеся коэффициентами устойчивых поли-
номов (1+5)ш_1 и (1+5)п“*. Ограничения /Доп и йдоп устанавливаются во избежание зацикливания алгоритма, которое может возникнуть, если требуемая точность решения не достигается из-за ошибок округления ЭВМ. Для спектральной факторизации
полиномов можно воспользоваться итерационным методом, изложенным в [5]. Как альтернативный вариант — возможно формирование спектрального фактора по корням факторизуемого полинома.
Представленный алгоритм минимизации б является ядром следующего алгоритма аппроксимации.
1. Подготовка. Вычислить элементы вектора Рг и задублировать /•'г
в . Задать нулевое приближение для коэффициентов знаменателя
апроксимирующей функции через вектор 0у—*. При отсутствии априорных данных по этим коэффициентам принять 0у—1= [1, 0,..., 0]т. Установить счетчик циклов последовательных приближений <7=1.
2. Минимизация функции е(0х, 0^).
2.1. Вычислить элементы матриц Нхх, НУх, Пуу, используя дополнительную весовую функцию п 0у9—1^-2.
2.2. Обратиться к алгоритму минимизации функции В? = В ( 0®, 0у).
При д=1 в качестве начальных условий ньютоновских итераций принять биномиальные коэффициенты (если нет априорных данных), при ?>1 для этой цели принять решение с предыдущего шага.
2.3. Проверить выполнение условий окончания итераций последовательных приближений: по малости величины б, (возможна также проверка условия 69»6*_1) и по числу итераций </><7доП (ограничение Яяоп устанавливается во избежание возможного зацикливания алгоритма). При выполнении условий'—уйти на 3.
2.4. Принять/7,* : = 0^~1 : = 0у, <7 : =<7+1 и уйти на 2.1.
3. Окончание решения. Выдать в качестве готового результата
а* : = а9, Ь*: = Ь«, 0* : = : = в*, Ч- = 5? •
Представленные алгоритмы реализованы в виде подпрограмм на языке программирования ФОРТРАН в арифметике двойной точности. Ниже рассматриваются примеры практических расчетов, выполненных на ЭВМ БЭСМ-6.
5. Примеры. Дробно-рациональная аппроксимация спектра Кармана.
Для статистического описания турбулентного движения воздушных масс хорошо подходит модель Кармана [9] со спектральной плотностью продольной и поперечной составляющих турбулентных пульсаций в виде:
(2) = а< V [1 + (1,339/;, о,у]516 ’
„ ,п 2 ь 1 + 8/3 (1,339/.п ^)2 *н(и)г [! (1>339/,л 2)2]п/6
где □ — пространственная частота, <71, ап — среднеквадратичное значение пульсаций скорости, /.(, /,„ — масштабы турбулентности. Эти спектральные плотности являются
иррациональными функциями, поэтому генерирование случайного процесса турбулентности по Карману, как правило, затруднено. Если для 5<(£2) и 5П(£3) построить аппроксимации в виде четных неотрицательных дробно-рациональных функций, то открывается возможность приближенного моделирования турбулентности Кармана с помощью легко генерируемого случайного процесса с рациональной спектральной плотностью. Чтобы не иметь дела с размерными физическими величинами, вместо (О) и (й) будем аппроксимировать функции
[1 +(1,339.)Ч5/6
(ч) = 2~т ^п (у1^п) ~
1 + 8/3 (1,339 м)*
[1 + (1,339 у)з]и/6 ’
где v=L<i „й — безразмерная частота.
В представленных ниже результатах вычислений коэффициенты аппроксимирующих дробно-рациональных функций 5(^2, 0Х, 0у) и 5„(г2, 0*,0И) определялись путем минимизации сумм квадратов невязок вида
е=1®2 0*))2’ к
взятых для 100 распределённых по логарифмической шкале в интервале у = 0,1-ь25 значений V*. Рассмотрены порядки аппроксимаций: т/п—1/2, 2/3, 3/4, 4/5. При минимизации е ставилось условие &Хт='1, т. е. принималось:
Рх = [0...-, 0,1 ]т, Ру — [0.0,0]т.
Число итераций последовательных приближений ограничивалось величиной <7Доп = П. Полученные в расчетах значения коэффициентов аппроксимирующих дробно-рациональных функций и коэффициентов минимально-фазовых факторов для этих аппроксимаций приведены в табл. 1 (продольные пульсации) и 2 (поперечные пульсации). Отметим, что минимально-фазовый фактор является передаточной функцией фильтра, с помощью которого белый шум преобразуется в случайный процесс с моделируемым спектром. Наглядное представление о качестве аппроксимации дают рис. 2, 3, где
м
>*»
Т а б л и ц а 1
Коэффициенты аппроксимации спектральной плотности атмосферной турбулентности (для продольной составляющей) модели Кармана
т/п ®X 1 ®Х 2 3 4 0у 1 2 з 0у 4 5 £
«1 ч Я3 «4 Ь\ Н К ь&
1/2 1,0 1,0 0,5049 0.7106 0,4929 0,7020 3,45-10-4
2/3 6,4590 2,5415 1,0 1,0 3,2324 1,7979 4,0408 2,4597 0,3122 0,5588 4,11-10-6
3/4 121,3773 11,0171 4,5954 8,2455 1.0 1,0 60,6963 7,7908 90,4719 13,3825 19,5134 5,6874 0,2299 0,4795 8,4Ы0-8
4/5 3379,576 58,1341 2123,162 64,2430 168,593 17,2360 1,0 1,0 1689,819 41,1074 2945,952 85,1989 1065,489 52,4589 56,2945 10,1081 0,1912 0,4373 1,53.10-э
10 Ученые записки № 3
Таблица 2
Коэффициенты аппроксимации спектральной плотности атмосферной турбулентности (для поперечной составляющей) модели Кармана
т}п 1 а1 9,2 -я2 3 а3 4 «4 0у 1 Ьг ®у 2 »2 ®у 3 Ьз 9у 4 ®у 5 Ь г
1/2 1,0 1.0 0,9018 0,9496 0.2198 0,4689 2,62-10-з
2/3 0,1591 0,3988 1,0 1,0 0,1607 0,4009 0.7522 1,0901 0,2959 0,5440 6,84-10"5
3/4 5,0372 2,2444 19,9085 4,9394 1,0 1,0 5,0419 2,2454 14,0801 5,4700 7,6276 3,5275 0,1937 0.4402 7,95-10-7
4/5 184,7791 702,7376 108,1172 1.0 184,8016 495,4352 325,3157 29,5863 0,1521 1,24-10-8
13.5933 32,5666 13,1625 1,0 13,5942 35,7649 28,8246 7,2158 . 0,3900
Ю
сл
показаны графики исходных нормированных функций спектральной плотности и их аппроксимации. Там же приведены графики относительной погрешности
8?М-(5(Л 0л еу)-5(м))/5(м).
Из характера поведения <55( (V) и 65п (\!) следует, что наибольшие относительные ошибки моделирования возникают на правом конце интервала аппроксимации. Они уменьшаются с ростом соотношения т/п и в случае т/га=4/5 не превышают 5-*-7% (меньшие ошибки относятся к продольной составляющей). В целях сопоставления полученных аппроксимаций с широко используемой моделью турбулентности Драйдена
[9] на рис. 2 и 3 приведены дополнительно графики нормированных функций спектральной плотности этой модели, рассчитанные по формулам:
_ 2 - 1 + 3^
А (V) =--------, 5- Ы =-------1---- .
’ 1 + V* " К (14- ^2)2
Видно, что к модели Драйдена наиболее близки аппроксимации с т/'/г=1/2 в случае продольной и с т/п=2/3 в случае поперечной составляющих турбулентности.
Упрощение передаточной функции. При исследовании управляемого движения летательного аппарата часто возникает потребность упрощения сложных передаточных функций, описывающих его динамику. Например, в [10] в связи с решением задачи проектирования робастной системы управления ставится вопрос об упрощении передаточной функции упругих колебаний корпуса ЛА, имеющей следующий вид:
0.686 (52-2809) (52—152.25+ 14 500) (52 + 153,8 5+ 14500)
(5^ = —■ .■■■■—■ . - I 11.
’ (в2-1-5 + 605) (52 + 4,55-5 + 2660) (52 + 2.515 + 3900) (в2 + 3.995 + 22 980) '
А
Попытаемся заменить ^(в) передаточной функцией более низкого порядка №(«), воспользовавшись методом среднеквадратичной аппроксимации функции [(ш) = | № (ш) |2. В качестве примера рассмотрим решение задачи для диапазона частот ш= 1-*-100 рад/с. Из анализа поведения /(со) следует, что в этом диапазоне присутствуют две резонирующие моды. Чтобы их учесть, знаменатель упрощенной передаточной функции должен иметь порядок не ниже четвертого. Поэтому численные расчеты проводились для значений я = 5 и т= 1н-5. В качестве критерия близости аппроксимации рассматривалась такая же сумма квадратов невязок, как и в предыдущем примере. Учитывая возможность существования локальных минимумов, расчеты для каждого значения т выполнялись четырехкратно с последовательным назначением единицы коэффициентам 6x1, 0хт, 0В 1, вУп и выбором наилучшего по точности аппроксимации варианта. По-
л
лученные, таким образом, для разных т передаточные функции Щя) сравнивались с № (в) по амплитудно-частотным и фазо-частотным характеристикам (АЧХ — ФЧХ). Выяснилось, что хотя увеличение т приводит к снижению ошибки модели по АЧХ, од-
нако может вести к повышению ошибки по ФЧХ._ Наилучшее соответствие по ФЧХ достигнуто при т=2 (при т = 3 получены практически те же результаты, так как •оказалось а3«0). Полученная для т = 2 передаточная функция имеет вид:
&(,)-________________________260 (26.1 + х)
2,364-10» + 5,662-103 s + 4,514-Юз s2 + 3,855 «з +
Графики АЧХ—ФЧХ этой функции и исходной функции Щя) приведены на рис. 4. Видно, что в рассматриваемом диапазоне частот обеспечивается удовлетворительное соответствие по АЧХ. Однако точность модели по ФЧХ недостаточна, поскольку с ростом частоты модель дает значительно меньший фазовый сдвиг, чем ее оригинал. С целью повышения точности моделирования ФЧХ к полученной передаточной функции л
Щз) было добавлено звено чистого запаздывания ехр(—тя). Оно не меняет АЧХ передаточной функции и поэтому оптимальную величину т легко определить из условия согласования только лишь ФЧХ. Путем минимизации суммы квадратов невязок по ФЧХ получена следующая формула для расчета т:
т = ^2 “* (аг8 ^ ~ аг2 Ш ш*))^2 “I •
В случае т=2 и использовании тех же узловых точек <оа, что и при формировании •ошибки аппроксимации функции /(со) получено т=0,047 с. Из рис. 4 видно, что ФЧХ
Л
передаточной функции И7(в) с добавлением запаздывания теперь уже вполне удовлетворительно согласуется с ФЧХ исходной функции №(х) и в такой комбинации позволяет сохранить основные динамические свойства упругого ЛА в заданном диапазоне частот.
ЛИТЕРАТУРА
1. Поляк Б. Т. Введение в оптимизацию. — М.: Наука, 1983.
2. Джу р и Э. Инноры и устойчивость динамических систем. — М.: Наука, 1979.
3. Sanathanan С. К-, KoernerJ. Transfer function synthesis as a ratio of two complex polynomials. — IEEE Transactions on Automatic Control, vol. AC-8, 1963, N 1.
4. J о n g М. Т., Shanmungam K- S. Determination of a Transfer Function from Amplitude Frequency Response Data. — Int. J. Control, 1977, vol. 25, N 6.
5. СупруненкоС. H. Среднеквадратичная аппроксимация четным неотрицательным полиномом. — Ученые записки ЦАГИ, 1988, т. 19, № 2.
6. К а л м а н Р. Е. Когда линейная система управления является оптимальной. — Теоретические основы инженерных расчетов. — Труды амер. общества инженеров механиков, 1964, № 2.
7. Чеботарев Н. Г., Мейман Н. И. Проблема Рауса—Гурви-ца для полиномов и целых функций. — Труды математического ин-та им. Стеклова, 1949, т. 26.
8. Гилл Ф., Мюррей У., Райт М. Практическая оптимизация.— М.: Мир, 1985.
9. Доброленский Ю. П. Динамика полета в неспокойной атмосфере.— М.: Машиностроение, 1969.
10. Lin J. М., Н а п К. W. Model reduction for control systems with parameter variations. — J. of the Franklin Institute, vol. 317, 1984, N 2.
Рукопись поступила 1/XIl 1988 г.
