CC BY
44
УЧЕНЫЕ ЗАПИСКИ ЦА Г И Том XIX 198 8
№ 2
УДК 629.7.051.062.2—52
СРЕДНЕКВАДРАТИЧНАЯ АППРОКСИМАЦИЯ ЧЕТНЫМ НЕОТРИЦАТЕЛЬНЫМ ПОЛИНОМОМ
С. Н. Супруненко
Предлагается метод решения задачи среднеквадратичной аппроксимации вещественной функции вещественным четным неотрицательным полиномом, основанный на минимизации ошибки аппроксимации по коэффициентам спектрального фактора аппроксимирующего полинома. Показано, что решение может быть найдено с использованием алгоритмов безусловной минимизации. Дается алгоритм решения с применением известных итерационных алгоритмов ньютоновского типа. Рассматриваемая полиномиальная аппроксимация может найти применение при исследовании систем управления ЛА в частотной области.
1. В некоторых задачах, связанных с анализом и синтезом систем управления летательных аппаратов в частотной области требуется полиномиальная аппроксимация следующего специального вида. Для вещественной функции /(со), заданной на отрезке изменения аргумента частоты сое{0, 1], необходимо определить вещественный четный полином наилучшего приближения
П (ш, 0) = 01 -{- 02 О)2 ... + Вп СВ2"-2
исходя из минимизации взвешенной среднеквадратичной ошибки аппроксимации
1
£ (6) = ]“(/(«.) _пк, е))*<*р(а)) (1)
и
по коэффициентам полинома при дополнительном условии
П(ш2, 0)>О для всех м. (2)
В (1) функция ф(со) —ограниченная и неубывающая на отрезке
со е [0, 1] — задает распределение весов аппроксимации. Относительно
аппроксимируемой функции /(<д) предполагается, что существуют степенные четные моменты
1
Рк == |/(ш) О)»-2 <*<р (со), £ = 1, ... , П. (3)
0
Отметим, что приведенная здесь запись выражения для среднеквадратичной ошибки аппроксимации с использованием интеграла Стильтьеса обеспечивает единообразный учет возможных на практике случаев формирования е(0) как путем обычного интегрирования (если ф(со) дифференцируемая функция), так и путем конечного суммирования (если ф(со) —ступенчатая разрывная функция). Что касается аргумента со, то для удобства записи в качестве него принята частота безразмерная, получающаяся нормировкой истинной частоты значением правого конца диапазона аппроксимации.
Особенностью рассматриваемой аппроксимационной задачи является необходимость соблюдения условия (2), в силу чего минимизация функции е(0) в общем случае требует привлечения алгоритмов условной минимизации [1]. Однако использование последних в данном случае встречает затруднения из-за сложности преобразования условия (2) к соответствующим нелинейным ограничениям на коэффициенты полинома П(со2, ф) (см. приведенные в [2] громоздкие соотношения для п = 4 и 5). Другой альтернативой могла бы быть минимизация е(0) в классе решения задач квадратичного программирования с континуумом ограничений (2). Применение этого подхода так же проблематично из-за необходимости поиска активных ограничений среди бесконечного числа заданных.
Предлагаемый в статье метод решения позволяет обойти данные затруднения путем преобразования исходной задачи условной минимизации к минимизации безусловной на основе использования соотношения спектральной факторизации. Известно [3], что вещественный четный полином, если он неотрицателен, может быть представлен в виде произведения двух комплексно сопряженных полиномов, имеющих вещественные коэффициенты, т. е. П(со2, 0) =Р(т, а) Р(—/со, а), где Р(з, а) =а1 + а2з + а382+...+ansn~i — так называемый спектральный фактор полинома П(со2, 0). При подстановке этого соотношения в (1) функция Ошибки аппроксимации преобразуется в зависимость от коэффициентов спектрального фактора и ее минимизация по этим коэффициентам может выполняться на открытом множестве. Вычислительная эффективность предлагаемого метода достигнута главным образом за счет того, что оказалось возможным получить в явной компактной зависимости от степенных моментов аппроксимируемой функции и от коэффициентов спектрального фактора все необходимые для построения практических алгоритмов минимизации соотношения. В статье определены условия единственности решения преобразованной задачи и показана возможность отыскания его с применением известных итерационных алгоритмов ньютоновского типа.
Отметим ряд конкретных примеров практического использования аппроксимации рассматриваемого вида. Так, при синтезе оптимальных линейных динамических систем на основне решения интегрального уравнения Винера—Хопфа, требуется спектральная факторизация числителя и знаменателя дробно-рациональной функции спектральной плотности [4]. Спектральная факторизация в рамках решения аппроксимационной задачи (см. соответствующий пример в данной статье) позволяет сделать точность решения частотно зависимой, а так же дает возможность расширить класс факторизуемых функций, включив в него в том числе и функции знаконеопределенные. Указанные обстоятельства отделяют соответствующий алгоритм решения от уже известных (см. например [5]). Аппроксимация четным неотрицательным полиномом может быть применена при синтезе линейных фильтров с желаемыми амплитудно-частотными характеристиками при условии, что заданы
нули или полюса искомой передаточной функции фильтра. Наконец, рассматриваемая аппроксимация может быть применена при упрощении передаточных функций систем управления с целью, например, облегчения последующего анализа системы, либо с целью придания ей определенных свойств робастности,[6].
2. Минимизация ошибки аппроксимации. Рассмотрим вначале исходную задачу минимизации функции е(0). Используя соотношение ПК, 6) = 2Т0, в котором й = [1, со2, ..., ш2п~2]т — я-мерный вектор, преобразуем (1) к виду
в(9) = 80 — 2Ртъ + втне,
где
1 1 «0 = ]7и2<*9>(®), Р= <*?(“),
о о
Элементы вектора Р определяются по формуле (3). Матрица Н — ган-келева, элементы ее образованы последовательностью степенных четных моментов:
1
К = °;+у-1 = [ “2{Ш~2) о??(ш). у'=1, •, п..
6
В силу принятых для функции ф(со) ограничений эта последовательность позитивна и соответственно матрица Н — положительно определена [7].
Выделим в евклидовом пространстве Яп, прямоугольными координатами которого являются элементы вектора 0, область С}, во внутренних точках которой полином П((о2, 0) положителен, а в граничных точках может принимать и нулевые значения. Вследствие линейности (по 0) ограничения (2) область С} представляет собой выпуклый конус [8]. Для решения рассматриваемой задачи существенно следующее утверждение: функция е(0) имеет в области С} единственный локальный минимум, являющийся также и абсолютным. Если бы не было ограничения (2), этот минимум легко можно было определить, приравнивая нулю градиент функции е(0). Соответствующее решение, полученное с учетом (4), имеет вид:
0*=//~1/7 е(9*) = е0 — Р* Н~1 Р.
В общем случае 0* (£ С}, поэтому для решения задачи необходимо применение более сложных методов поиска минимума, учитывающих (2).
Возможным способом получения решения с нужными свойствами является сведение исходной задачи условной минимизации к минимизации безусловной путем нелинейного преобразования аргумента минимизируемой функции е(0). Предлагаемое преобразование базируется на использовании вспомогательного полинома и двух специальных матриц. Поэтому дадим вначале соответствующие определения.
Вспомогательным к исходному полиному Р(з, а) назовем полином Р(5, х), коэффициенты которого определяются по правилу:
** = (—ак, 1................
(Квадратные скобки здесь означают символ целой части числа.)
(4)
1
н= | аат с(<р (со).
Расширенной матрицей Гурвица (РМГ) вектора х = [хи ..., х„]т назовем пХп матрицу Г (х), элементы которой определены правилом:
Г Ху-1 при 1<2/ —/<я,
~ 1 0 в остальных случаях.
Обычная матрица Гурвица, рассматриваемая в задаче устойчивости полинома [9], отличается от РМГ отсутствием первых строки и столбца.
Четно-разреженной ганкелевой матрицей (ЧРГМ) вектора х назовем пХп ганкелеву матрицу (/(х), элементы которой определены правилом:
Г Х(1+п12 при /+7 —четном,
^1’ 1 0 при Ь -{- 7 — нечетном.
В матрице ЧРГМ каждая четная косая диагональ нулевая.
Матрицы РМГ и ЧРГМ обладают двумя важными свойствами, которые будут использоваться при выводе ряда формул алгоритма минимизации:
Г(у)тх = Г(л;)т_у, Г (y)x = G(x)y.
Здесь х и у— произвольные векторы одинаковой размер.ности.
С учетом введенных понятий можно получить следующее представление для четного неотрицательного полинома:
П (ш2, 6)= ЙТГ (х)т х.
Здесь х — вектор коэффициентов вспомогательного к спектральному фактору P(s, а) полинома.
Из этого представления нетрудно усмотреть формулу взаимосвязи для коэффициентов полиномов П(ю2, 0) и P(s, х):
6 = Г {хУ х. (5)
Подставляя (5) в (4), преобразуем тем самым функцию ошибки аппроксимации к явной зависимости от х:
г (х) = е0 — 2FTГ (х)тх + х'Г(х)НГ(х)тх. (6)
Нелинейное преобразование (5), назовем его здесь отображением Ф, отображает всю открытую область Rn определения вещественного вектора х в замкнутую область Q, в которой коэффициенты полинома П(со2, 0) удовлетворяют условию его неотрицательности, т. е. имеет место #: Rn-+Q. Поэтому поиск минимума функции е(х) может проводиться на открытом множестве. Так как между коэффициентами полиномов P(s, а) и P(s, х) имеется простое взаимно однозначное соответствие, то рассматриваемое преобразование по своей сути означает переход к минимизации ошибки аппроксимации непосредственно по коэффициентам спектрального фактора.
Следует, однако, отметить, что отображение 0 свойством гомеоморфизма, т. е. взаимно-однозначного соответствия не обладает. Это слепляющее отображение, поскольку разным x^Rn может соответствовать
одна и та же точка 0е<2. Как следствие, несмотря на то, что функция б(0) имеет в области С? единственный локальный абсолютный минимум, функция г(х) в открытой области Яп может иметь множество не обязательно равнозначных локальных минимумов. Необходимо выяснить, какие из локальных минимумов функции г{х) в Ип гарантированно отображаются в абсолютный минимум функции е(0) в <3.
Выделим в евклидовом пространстве Яп, прямоугольными координатами которого являются элементы вектора х, п связных областей Од, (? = 0, 1 — 1, во внутренних точках которых корни полинома
Р(5, а(х)) принадлежат правой полуплоскости комплексного аргумента 5, а в граничных точках среди этих корней появляются чисто мнимые [10]. Области Од являются коническими, причем, если хе/),, то и —х^Од. Так как смена знака х мало что добавляет к процессу определения минимума функции е(х), то в дальнейшем за область Од будем принимать только часть конуса, соответствующую положительным значениям XI.
У каждой области существует образ (2д = 0(/)д), который при не обязательно совпадает со всей областью ф, а возможно составляет только часть ее. Например, при п = 3 имеет место
62;> — 2/Мз, в80). В то же время область =
{л:3 ^0} преобразованием & отображается в (31 = {б1;>0, 62;>
>2^01 03 , 93 >0} с:<2. Заметим, что область к тому же не выпукла.
Нарушение условия $(Од)=С} может привести к тому, что среди точек Од не найдется прообраза абсолютного минимума &(0) на С}. При отсутствии выпуклости <Зд возможно появление не абсолютных, а только локальных минимумов &(0) на <Зд а <3, для которых имеются прообразы в областях Од. Эти ситуации наглядно показаны на рис. 1.
Данные обстоятельства приводят к необходимости поиска минимума функции е(х) только в тех областях Од, которые гарантированно отображаются на всю область (?. Известно [3], что таким свойством обладают области А) и /)«-1- Так как в практических приложениях основной интерес представляет устойчивый спектральный фактор [3, 4], то в дальнейшем будем рассматривать только одну область £>0. В связи с этим введем ограниченное отображение которое как и Ф задается формулой (5), но действие которого ограничено пределами области О,. Это отображение обладает тем важным свойством, что оно обеспечи-
т>(Л0) = !} ц!-
вает взаимно однозначное соответствие между точками областей и С?. Так, если 0* — точка абсолютного минимума е(0) на <3, то найдется прообраз этой точки х* = '&о(0*)-1, причем единственный, принадлежащий О0. В силу непрерывности и дифференцируемости функции е(х) этот прообраз должен быть точкой локального минимума функции.
Обратимся к задаче отыскания локальных минимумов функции е(х). Раскладывая г(х) в ряд Тейлора в окрестности произвольной точки х и выполняя преобразования, учитывающие свойства матриц РМГ-ЧРГМ, нетрудно получить следующие выражения для вектора градиента и матрицы вторых производных минимизируемой функции:
V*6 (*)£(*); (7)
у*, г (х) = АО (ё (х)) + 8Г (х) ЯГ (х)т, (8)
где
g(x) = HГ(x')x-F^±Ъe(e).
В экстремальных точках градиент функции е(х) равен нулю, т. е.:
Г(х)£(х) = 0. (9)
Структура уравнения (9) предполагает возможность существования решений, удовлетворяющих условию £(х) = 0, т. е. уравнению
НГ(х)тх — Р = 0. (10)
Нетрудно заметить, что решениям уравнения (10) соответствует образ 0* = Н~1 Р, являющийся глобальным минимизатором функции е(0) на открытом множестве. Поэтому все вещественные решения уравнения (10), если они имеются, обеспечивают абсолютный минимум ошибки аппроксимации с соблюдением условия (2). Среди таких решений найдется и точка, принадлежащая области £)0.
Рассмотрим теперь решения (9), которые уравнению (10) не удовлетворяют. Из условия разрешимости (9) относительно £(-*0=0 следует, что такие экстремальные точки должны принадлежать гиперповерхности в Яп, задаваемой соотношением с!е1:Г(л;) = 0. Можно показать, что если эти точки принадлежат также£)0, то спектральный фактор находится на границе устойчивости (т. е. граница области £>0 является частью указанной гиперповерхности). В зависимости от знаковой
определенности матрицы экстремальная точка может быть сед-
лом, минимумом или максимумом. Из анализа выражения (8) следует, что локальный максимум возможен только в точке х = 0. Так как в пределах области /)0 имеется только одна точка локального минимума, то экстремальными точками, принадлежащими границе области Б0, являются седла и, если не существует вещественных решений (10), одна точка минимума.
На рис. 2, а, б показаны изолинии уровня функции е(х) при п = 2 с вариацией параметров функции таким образом, чтобы минимум принадлежал либо внутренности £>0 (рис. 2, а), либо границе Д, (рис. 2, б). В данном примере п = 2 и поэтому областью Д, является правый ниж-
Рис. 2
ний квадрант в координатах х1х2. Точечным пунктиром на рисунках отмечены линии уровня, проходящие через седловые точки.
3. Итерационные алгоритмы поиска минимума. Для отыскания точек локального минимума функции г(х) можно применить большинство известных методов безусловной минимизации. Если удается обеспечить принадлежность решения области Д, то найденный таким образом локальный минимум будет абсолютным. Наличие явных выражений для вектора градиента (7) и матрицы вторых производных (8) минимизируемой функции делает привлекательным применение ньютоновских 'итерационных алгоритмов минимизации, обладающих, как известно [1]. свойством квадратичной сходимости. Рассмотрим эти алгоритмы.
Пусть хк — известное приближение к точке локального минимума на шаге итерации к. Решение на шаге й+1 представляется в виде хк+1 = хк -\-аькхк, где вектор Ахк и параметр задают соответственно направление и длину шага итерации.
Алгоритм Ньютона (АН). В исходном варианте этого алгоритма, именуемом иногда алгоритмом Ньютона — Рафсона, принимается аь=1 и Ахк = — (у^8(л*)ГЧ^(л*).
С учетом соотношений (7) и (8) в данном случае имеем:
Дх* = — У(х*)-1 Г (л:*) ^(л:4), (И)
где
J(xk) = 2Г(xk)HГ(xky + G[g(xk)\, g{xk) = H(Y){xkyxk -Р.
Для того чтобы вектор Ахк был направлением убывания е (л:), необходима положительная определенность матрицы J{xk) [1]. Анализ струк-
туры этой матрицы показывает, что областью положительной определенности является все открытое пространство исключая центральную его часть, имеющую форму звезды. Благодаря этому всегда можно обеспечить положительную определенность матрицы J (х) путем необходимого удлинения вектора х. Другим способом обеспечения положительной определенности матрицы У (х) является модификация этой матрицы в случае нарушения данного свойства. Вопросы восстановления положительной определенности матрицы путем некоторой модификации ее элементов достаточно широко освящен в соответствующей литературе. Отметим здесь только наиболее удобный в вычислительном отношении метод Гилла и Мюррея [11], который предусматривает автоматическую коррекцию диагональных элементов матрицы в процессе решения системы уравнений типа (11) с разложением матрицы системы на треугольные множители (т. е. с Ьи или 1^Е)и факторизацией матрицы). Ввиду универсальности и простоты метода, применение его в алгоритме АН наиболее предпочтительно.
Для обеспечения гарантированной сходимости итераций к локальному минимуму, особенно в случае возможных модификаций матрицы J (х), необходимо так же оптимизировать длину шага итераций. Практическая проверка алгоритма АН с модификацией матрицы J (х) и оптимальным выбором параметра а показала, что обеспечивается сходимость итераций к локальному минимуму функции г{х) независимо от стартовой точки. Однако найденное решение не обязательно принадлежит Дь даже если итерации начинались из этой области.
Алгоритм Гаусса—Ньютона (АГН). От рассмотренного выше этот алгоритм отличается только тем, что вместо полной матрицы вторых производных используется ее «укороченный» аналог. В данном случае это матрица J(х) = 2Г (х) НГ (х)с. Из (11) с использованием такой замены можно получить:
2Г (л:*)т Ахк — 6* — Ьк, (12)
где 0й = Г (хкУ хк, 0* = Н~х Р.
Специальная структура матрицы Г(х) дает возможность, используя метод исключения Гаусса [9], построить быстродействующий алгоритм решения (12). Подобный алгоритм рассмотрен в [5]. Однако эффективность данного подхода теряется, если точка 0* близка к границе области (2, поскольку в точке решения матрица Г(х) близка к вырождению. Получение решения непосредственно из уравнения типа (11) более надежно, так как имеется возможность обеспечения строгой положительной определенности матрицы /(*) путем ее модификации как и в алгоритме АН.
Алгоритм АГН обладает весьма важным свойством, выделяющим его среди всех прочих алгоритмов поиска минимума е(х). А именно если 0*еС, то итерации АГН с параметром длины шага а<1, начавшиеся в области А>, эту область не покидают.
Данное свойство позволяет использовать алгоритм АГН для целей спектральной факторизации четного неотрицательного полинома с получением устойчивого спектрального фактора. Отметим, что условие 0* е <2 является весьма существенным для работоспособности алгоритма. Если оно нарушается (т. е. существуют значения со, при которых полином П(со2, 0*) отрицателен), то при подходе итераций АГН к точке минимума возникают осцилляции итераций с пересечением границы области А». Подобный пример показан на рис. 3.
Рассмотрим возможность повышения эффективности ньютоновских итераций путем оптимизации длины шага и масштабирования х.
Оптимизация длины шага. Используя (6) и (11), получим величину приращения е(х) вдоль направления Ах в зависимости от а:
Ае (а) = — 4А1 а + 2Л2 а2 + 4А3 а® -}- А4 а*,
где
А1 = ^x!J (х) Ах, А2 = Ахт/(х)Ах,
Л3 = хтГ (Ах)//Г (Ах)тАх, Л4 = ДхТ (Ах) ЯГ (Ах)1 Ах.
А
Через 3(х) обозначена фактически используемая модификация матрицы вторых производных У (х). В области сходимости итераций матрица
Л
У (х) должна быть положительно определена, поэтому Л^О (при АхФО). Очевидно, так же и Л4>0, поэтому функция Де(а) обязательно имеет один или два минимума при положительных а. Эти минимумы удовлетворяют условию дЛе(а)/<За = 0 и поэтому могут быть определены как корни уравнения Л2а+ЗЛза2-!-Л4Сс3=Л1. Рекомендуется использовать наименьший положительный корень.
Оптимальное масштабирование. Путем выбора масштабного множителя в преобразовании х-^-хх можно обеспечить дополнительную минимизацию ошибки аппроксимации без существенных вычислительных затрат. Используя (4) получим
г (х) = е0 — 2т2 Г 6 + х4 6Т Я0.
/ /?т 0 \1/2
Минимум функции г (г) достигается при тор4 = I "еТ/лГ) и Равен
(рч 0)2
г(?6Р1) = г0—-рЖ' Масштабирование х удобно тем, что оно не ме-
няет статуса принадлежности х области О0.
Комбинированный алгоритм. Рассмотренные алгоритмы АН и АГН вместе с соответствующими модификациями 'можно объединить в
рамках единого алгоритма поиска минимума е(х) с гарантированным решением в области Д>. В предлагаемом ниже комбинированном алгоритме с помощью итераций АН отыскивается какой-либо локальный минимум е(х) в открытой области Яп, а итерациями АГН делается перенос решения в область Д, если оно этой области не принадлежало. Перенос решения в область Д по существу означает спектральную факторизацию четного неотрицательного полинома с получением устойчивого (возможно и на границе устойчивости) спектрального фактора. В укрупненном виде комбинированный алгоритм состоит из следующих шагов.
1. Задать стартовую точку х°еД и установить счетчик циклов переноса решений <7=1.
2. Проверить выполнение условия окончания счета (например, по малости нормы градиента IIУхе(х°) ||). В случае выполнения условия окончания уйти на 6.
3. Найти локальный минимум е(х), используя итерации АН с оптимальным выбором т, а и по необходимости с модификацией матрицы У(х). Итерации начинать с точки х°. При <7=1 выполнять проверку принадлежности итераций области А> с запоминанием в х° последней итерации, еще принадлежащей Д (учитывая возможность выхода итераций АН из области О0). При <7>1 такой проверке подвергается только найденное решение х*. Если найденное решение х* (<7^1) принадлежит В0, то принять х°=х* и уйти на 6.
4. Перенести решение в Д> путем спектральной факторизации. Для этой цели применить алгоритм АГН с ао<1 (возможно с минимизацией по т и а, а так же с модификацией матрицы У (х) как и в алгоритме АН), используя в качестве стартовой точку х° и минимизируя функцию е(х), в которой принято Р—НГ(х*)тх*. Запомнить найденное решение в х°. Если <7 = <7дсга (ограничение <7д0п устанавливается во избежание возможного зацикливания алгоритма), уйти на 6.
5. Принять <7 = <7+1 и уйти на 2.
6. Выполнить преобразование а = а(х°) и выдать а в качестве готового решения.
Так как на шаге 3 рассматриваемого алгоритма обеспечивается убывание ошибки аппроксимации е(х), а на шаге 4 при переносе решения в область Д> величина этой ошибки не возрастает (при условии, что спектральная факторизация выполняется достаточно точно), то на каждом цикле итераций ц происходит непрерывное убывание е(х) с получением решения в области Д. В целом, сходимость итераций гарантируется сходимостью используемой модификации АН.
Отметим, что в качестве стартовой точки х° на шаге 1 может быть использована точка х(а), полученная с использованием спектрального фактора нулевого приближения следующего вида: Р(в, а) = (1 + 5)п-1. Коэффициенты этого полинома являются биномиальными коэффициентами, т. е.: = &=1,...,я. Требуемая на шаге 3 проверка при-
надлежности итераций области А, т. е. проверка устойчивости спектрального фактора Р{в, а (х)) может производиться непосредственно по коэффициентам вспомогательного полинома Р{з, х) (см. приложение) .
4. Примеры. Предложенный выше комбинированный алгоритм реализован в виде программы на алгоритмическом языке ФОРТРАН. Рассмотрим конкретные примеры, иллюстрирующие, работоспособность этого алгоритма.
Спектральная факторизация. Для четного полинома П (<о2, 6) = = 2 — 6ш2 —[— ЗОси4—104ш6 100и>8, который на интервале ве[0, 67,
О, 80] отрицателен (на рис. 4 соответствующая зависимость показана сплошной линией) и поэтому не имеет спектрального фактора с действительными коэффициентами, выполнена приближенная спектральная факторизация путем минимизации среднеквадратичной ошибки аппроксимации следующего вида:
1
е = ^[П(о>2, 6) — Р(ш, а)Р(—ш, а)]2о)й?ш.
о
В данном случае использовалась функция распределения весов аппроксимации ср (св) = -1— со2, обеспечивающая повышение точности аппроксимации с ростом аргумента частоты со. При поиске минимума е(х) в качестве стартовой точки итераций АН принято х°=[1; —4; —6; 4; 1]т. Останов итераций производился по достижению заданной величины относительного приращения ошибки аппроксимации (принята величина 10-11). В процессе решения потребовалось выполнить 12 итераций АН, 18 итераций АГН (перенос решения в £>0), затем две итерации АН
(стартовой точкой этих итераций было решение предыдущих итераций АГН) и одна итерация АГН (перенос решения в Б0 с использованием той же стартовой точки, что и в предыдущих итерациях АН). Получено:
{0Й} = {2,273; —12,587; 62,190; —152,666; 122,844},
[ак] — {1,508; 3,784; 8,923; 6,718; 11,083).
Поведение аппроксимирующего четного полинома показано пунктиром на рис. 4. Спектральный фактор находится на границе устойчивости.
Упрощение системы управления ЛА. При разработке систем управления летательных аппаратов в настоящее время широко применяются методы современной теории управления, базирующиеся на идеях синтеза оптимальных динамических систем в пространстве состояний. Однако в ряде случаев эти методы приводят к усложнению системы управления, к повышению чувствительности ее к возможным изменениям
характеристик объекта управления. В [12] для модели динамики ЛА по крену с передаточной функцией
Гл* (5) - 4500 [і (1 + 5,2) (і + тет- + -да) (> + 2ШГ + 2да) х:
приведена следующая передаточная функция оптимального регулятора:
Для работы этого регулятора необходимо обеспечить измерение восьми переменных состояния, что может оказаться не выполнимым при практической его, реализации. Попытаемся уменьшить число измеряемых переменных, заменив исходный полином №су($) полиномом
л
№ су (5) более низкого порядка, достаточно хорошо аппроксимирующим №су(5) в основном диапазоне рабочих частот регулятора. Воспользуемся для этой цели методом среднеквадратичной аппроксимации функции /(со) = | Шсу(ш) |2. В результате аппроксимации должен быть получен
четный неотрицательный полином | №су(«й) |2 и его устойчивый (воз-
л
можно на границе устойчивости) спектральный фактор И7су(«). В качестве примера рассмотрим результаты расчетов, полученные с использованием среднеквадратичной ошибки аппроксимации в виде конечной суммы:
Суммирование в этой формуле выполнялось для последовательности значений частоты т, к= 1,..., 250, распределенных по логарифмической шкале в заданных диапазонах аппроксимации. В расчетах рассматри-
вались полиномы №су($) 1, 2, 3-го порядков, для которых установлены разные диапазоны аппроксимации: от 1 рад/с до соответственно 50, 100, 150 рад/с. Запись ошибки аппроксимации в виде конечной суммы означает, что функция распределения весов аппроксимации <р((о) является разрывной ступенчатой со скачками в точках со^. Поэтому при вычислении степенных моментов, из которых формируются вектор Р и матрица Н, интегралы в соответствующих формулах расчета степенных моментов заменялись на конечные суммы с суммированием по указанным выше значениям частоты &>&. В конечном итоге получены следующие аппроксимации исходного регулятора:
Гоф) = 0,05 (1 -М/12.1)(і + + -§г)(і + + -^і) X
Л
£ = Х(|№су (*»,) і 2 - і і^су (*Ч) і *)*•
ь
Л
НРсу («) = 0,05 (1 + 8,2-10-*2 5),
ІГсу (5) = 0,05 (1,152 + 8,1 • 10-г х + 8,2-10~4 я2),
ГСУ (5) = 0,05 (1,074 + 9,1 • 10~2 Я + 1,18 • Ю-Зв2 + 6,5 • 10"5 53).
6—«Ученые записки» № 2
81
На рис. 5 показаны амплитудно-фазово-частотные характеристики разомкнутого контура ЛА+СУ с исходным регулятором 7-го порядка и с регуляторами пониженного порядка, полученными путем аппроксимации. Видно, что даже с регулятором 1-го порядка могут быть обеспечены вполне удовлетворительные динамические характеристики ЛА с системой управления (в данном случае запасы устойчивости по амплитуде и фазе равны соответственно 5,03 дБ и 33,5°).
ПРИЛОЖЕНИЕ
В данном приложении приводятся доказательства наиболее существенных для рассматриваемой задачи утверждений, упоминаемых в статье.
1. Единственность локально-абсолютного минимума е(в) на ()■
Множество ф, во внутренних точках которого полином П(<в2, 9) положителен, образовано пересечением замкнутых полупространств {6 £ А?” : Ут 6 ;> 0] и поэтому является выпуклым замкнутым конусом [8]. В силу положительной определенности матрицы Н функция е(0) является выпуклой всюду в Яп. Выпуклая функция имеет на выпуклом замкнутом множестве единственную точку локального минимума, являющуюся так же и минимумом абсолютным [1].
2. Векторно-матричное представление для четного неотрицательного полинома П(со2, 0).
Находим:
п п
п («3, в) = Р(т, а)Р(— Ы, а) = ^ (т)к+1-2 ак аг ^ “6+,~2-*й хг .
к, г-1 к, 1 = 1
к+1=2т к+1=2т
Полученное выражение есть квадратичная форма хТ0(&)х, в чем легко убедиться, если учесть, что элементы матрицы 0(2) равны g/г^ = юк+1~2 при к+1 четном и gkl — 0 при к-\-1 нечетном. Учитывая, что О (2) х = Г (х) 2, получим П (о>2, 6) = = 2Т Г (*)т х.
3. Основное свойство итераций АГН при минимизации е(х).
Справедливо следующее утверждение. Если Н~1 У7 = 6* £ и хк <= £>0 (здесь
символ М означает внутреннюю часть области без ее граничных точек), то для всех значений параметра длины шага итерации а/* е [0, 1) выполняется хк+1 £ шЮ0> 3 для 01/1=1 выполняется Хк + 1£Оо.
а| Г (Д**)т bxk = Г {xk+xy xk+l - 2Г (xkf xk+1 + Г (xky xk
и с учетом (12) получим:
а\ Г (Дл:*)т Д** =0*+I-(l — ak) 6* — я* 6*.
Так как для всех Axk£Rn справедливо 2Т Г (Длс*)т Дат*^0, то должно выполняться
QT6*+l^aA QT0* + (1 _ а/г) QT0* > (! _aft)2T0ft.
В рассматриваемом случае xk£inl D0, поэтому и 0* ^ int Q. Но тогда для всех 1) должно выполняться 2Т0Й+1>О, т. е. 6ft+1£intQ. Так как в данном случае граница Q не достигается, то в силу непрерывности зависимости л'*+1 от я* и в силу свойства гомеоморфизма отображения *>0 не достигается и граница D0, а значит xk+1£int D0. Если же а* = 1, то х*+* является предельной точкой последовательности {xft+x (aft) ^ int D0} при ak 1. Очевидно, что самое многое эта точка может принадлежать границе D0. Поэтому в общем случае xk+lcDo.
4. Проверка устойчивости спектрального фактора по коэффициентам вспомогательного полинома.
Докажем вначале следующее утверждение.
Главные диагональные определители k-ro порядка матриц РМГ векторов исходного и вспомогательного полиномов отличаются только знаками, причем имеет место
det (Г (a))ft = (-l)^det(r (x))k,
где
р(Л) = [(А + 2)/4]-[(*+1)/4].
Доказательство. Рассмотрим знаковую структуру матрицы, полученной из матрицы Г (а) заменой элементов вектора а элементами вектора х. Она имеет вид:
“ + — + — + — ■ • •
0 — 4- - + - • • •
0 + — + — + ...
0 0 - + — + ...
0 0 + — + — ...
Если на —1 умножить все элементы четных столбцов этой матрицы, а затем все элементы строк 3, 4, 7, 8,..., то знаковая структура матрицы станет полностью положительной. Элементарный подсчет показывает, что общее по числу строк и столбцов количество таких умножений нечетно только для порядков матрицы, равных 2, 6, 10, 14,.... Этому циклическому чередованию соответствует приведенная выше функция четности ц(6). Применение теперь известного правила умножения строк и столбцов определителя на число [9] и завершает доказательство.
Используя данное утверждение, преобразуем детерминантный критерий Гурвица [9] проверки устойчивости полинома P(s, а) к системе неравенств, включающих только коэффициенты вспомогательного полинома P(s. х):
(— 1 )•*■<*>det (Г (х))й > 0, * = п.
На основании этих неравенств нетрудно получить и модификацию критерия Рауса [9]:
С — 1 >^*/21 £>0, k=\.п. Здесь Я, j, — k-ft элемент первого столбца таблицы
Рауса, составленной по коэффициентам вспомогательного полинома P(s, х). Отметим, что если исходный полином P(s, а) находится на границе устойчивости, т. е. среди его корней имеются чисто мнимые, то должно выполняться det (Г (л:))п = 0.
1. Поляк Б. Т. Введение в оптимизацию. — М.: Наука, 1983.
2. Д ж у р и Э. Инноры и устойчивость динамических систем. — М.: Наука, 1979.
3. К а л м а н Р. Е. Когда линейная система управления является оптимальной. — Теоретические основы инженерных расчетов. Труды амер. общества инженеров-механиков, 1964, № 1.
4. Чанг Ш. А. С. Синтез оптимальных систем автоматического управления.— М.: Машиностроение, 1964.
5. V о s t г у Z. ‘New algorithm for polynomial spectral factorization with quadratic convergence. — Kybernetika (CSSR), vol. 12, 1976.
6. Б e с e к e p с к и й В. А., Н е б а л о в А. В. Робастные системы автоматического управления. — М.: Наука, 1983.
7. Ахиезер Н. И. Классическая проблема моментов.— М.: Физ-матгиз, 1966.
8. Никайдо X. Выпуклые структуры и математическая экономика.— М.: Мир, 1972.
9. Гантмахер Ф. Р. Теория матриц. — М.: Наука, 1967.
10. Чеботарев Н. Г., М е й м а н Н. Н. Проблема Рауса—Гурвица для полиномов и целых функций. — Труды математич. ин-та им. Стеклова, 1949, т. 26.
И. Гилл Ф., Мюррей У., Райт М. Практическая оптимизация.— М.: Мир, 1985.
12. Несли Ф. В., Зархан П. Причины возникновения неустойчивости при практическом применении систем, разработанных с помощью современной теории управления. — Аэрокосмическая техника, 1985, т. 3. № 3.
Рукопись поступила 25/XII 1986
