CC BY
68
Пушков С.Г., Горелик А.А.
Оренбургский государственный университет E-mail: psg@mail.osu.ru; anna_gmn3@rambler.ru
ИСПОЛЬЗОВАНИЕ МЕТОДОВ ИНТЕРВАЛЬНОГО АНАЛИЗА ДЛЯ ВЫЧИСЛЕНИЯ РАЗМЕРНОСТИ КОНЕЧНОМЕРНОЙ РЕАЛИЗАЦИИ ЛИНЕЙНОЙ ДИНАМИЧЕСКОЙ СИСТЕМЫ
В статье рассматривается задача построения модели в пространстве состояний на основе данных о поведении вход-выход линейной стационарной динамической системы с дискретным временем. Исследуется проблема получения численных алгоритмов вычисления конечномерной реализации в пространстве состояний отображения вход-выход, представленного импульсной последовательностью матриц. Предложен метод вычисления размерности такой реализации, основанный на использовании интервального анализа.
Ключевые слова: динамические системы, пространство состояний, проблема реализации, численные методы, интервальный анализ.
Введение
Решение многих прикладных задач сводится к построению математической модели некоторого процесса или явления. Зачастую моделирование основывается на результатах наблюдения за поведением некоторой динамической системы и сводится к представлению экспериментальных данных в форме, удобной для дальнейшего исследования. Целью такого моделирования может быть прогнозирование реакции объекта на различные входные возмущения и управление объектом. Наблюдения за поведением системы позволяют задать ее с помощью отображения «вход-выход». В этом случае для определения реакции системы необходимо знать и использовать всю предысторию входов и выходов, что не всегда удобно осуществить на практике. Поэтому часто ставится задача задания динамической системы, используя ее свойства (пространство состояний) таким образом, чтобы реакция системы определялась через ее состояние и входной сигнал в данный момент времени. В теории систем эта проблема известна как задача реализации.
Для класса линейных стационарных динамических систем с дискретным временем задача построения конечномерной реализации успешно решается в рамках алгебраического подхода к проблеме [1]. Одним из первых алгоритмов вычисления конечномерных реализаций для систем над полями является алгоритм Б.Л. Хо, который основан на приведении ганкелевой матрицы отображения вход-выход-системы к диагональному виду [1]. В [2, 3] представлены варианты численной реализации этого алгоритма.
В настоящей работе освещены проблемы численной реализации алгоритма Б.Л. Хо, связанные с неточностью представления данных в ЭВМ. В частности, рассматривается вопрос о вычислении размерности пространства состояний линейной динамической системы с помощью методов интервальной арифметики.
1. Задача реализации импульсной
последовательности
Поведение вход-выход линейной стационарной многомерной управляемой системы с дискретным временем может быть охарактеризовано импульсной последовательностью матриц
Кл2,...|, (1.1)
являющейся представлением импульсной характеристики системы. В этом случае для заданной последовательности векторов управлений (входной последовательности) и( 0), и(1),... выходная последовательность векторов у( 1 ),у( 2),... определяется соотношением:
г
уО) = '^А1и(1-г), г = 1,2,к .
г=1
Соответствующая модель в пространстве состояний определяется по следующим формулам:
х(г +1) = Fx(г) + Ои(г), у(г) = Нх(г), (1.2)
где г е 2 - момент времени; х(г),х(г+1) - векторы состояния системы в моменты времени г и г+1; и(г) - вектор входного сигнала и у(г) - соответствующий ему вектор выходного сигнала; F,G,H - линейные отображения.
Проблема построения модели с пространством состояний на основе данных вход-выход в теории систем известна как задачареализации [1].
Математическая формулировка этой задачи для класса линейных динамических систем с дискретным временем состоит в следующем:
- для заданной последовательности матриц (1.1) размера pхш над полем Kпостроить тройку матриц (Д G, Я) над тем же полем К, таких,что
д. = HFi~lG,(i = 1,2,...;, (1.3)
где для некоторого п (которое также нужно определить) Я - матрица размера р х п , F - матрица размера п х п , G - матрица размера п х ш.
Критерий реализуемости импульсной последовательности матриц дает следующая теорема, доказанная в [1]:
Теорема 1.1. Импульсная последовательность матриц { ,а2 ,...} реализуема тогда и только тогда, когда она рекуррентна, то есть когда существует такое г и такие коэффициенты Р1 ,Р2,...рг е я, что
Г
аг+/+1 = -£р.д+. (] = 0,1,...). (1.4)
.=1
Данное условие называют условием рекуррентности.
Последовательности (1.1), а следовательно, и отображению вход-выход / можно поставить в соответствие бесконечную в двух направлениях матрицу:
Ч А2
в(1) = '
которая является блочно симметричной и называется ганкелевой матрицей исходного отображения. Введем также в рассмотрение ганкеле-вы блоки размера д'хд:
' Ai А2 " Aq )
Bq'q(f) = А2 А3 " + q
, А + • Aq+q'-1
На B(f) и Bq,q(f) можно ввести оператор сдвига о следующим образом:
о kB(f) = B( okf) =
A2+k A3+k
Введение ганкелевой матрицы отображения вход-выход в случае систем над полями позволяет сформулировать следующий критерий реализуемости [1]:
- отображение f реализуемо тогда и только тогда, когда матрица B(f) имеет конечный ранг. Этот критерий также можно сформулировать в другой форме:
- отображение f реализуемо тогда и только тогда, когда существует целое п, такое, что rank Brr (f) < n для любых положительных r.
2. Алгоритм Б.Л. Хо построения конечномерной реализации
Для формулировки алгоритма построения реализации введем в рассмотрение матрицу Em, которая содержит единицы на главной диагонали и нули на всех внедиагональных позициях.
Одним из известных алгоритмов вычисления конечномерных реализаций для систем над полями является алгоритм Б.Л. Хо [1]. Этот алгоритм основан на следующей теореме:
Теорема 2.1. Пусть для последовательности матриц {a1 ,a2,...} размера p хm над полем K существует такое r, что выполняется условие рекуррентности (1.2). Тогда, если матрицы P размера pr х pr и M размера mr х mr над Kудов-летворяют условию
P[Brr(f)]M =
ln On
n mr-n
OPr-n oPr-
^ УІ ^ in V—
= E^Elr (2.1)
для n = rank Brr(f), то каноническую реализацию отображения f можно вычислить по формулам:
F = EnprP[( oBUf)]MEnmr,
G = EnprP[Brr(f)]Emr,
H = EpPr[Brr(f)]ME'mr. (2.2)
Таким образом, алгоритм Б.Л. Хо вычисления конечномерной реализации состоит из следующих этапов [1]:
1) выбор числа r, удовлетворяющего условию рекуррентности (1.4) импульсной последовательности матриц;
2) построение невырожденных матриц P и M, удовлетворяющих условию (2.1);
3) вычисление матриц F, G, H канонической реализации по указанным в теореме формулам (2.2).
При численной реализации этого алгоритма на ЭВМ возникают трудности, связанные с бессмысленностью понятия бесконечной матрицы для машинного представления, а также с ограниченной точностью представления чисел, которые могут сильно исказить решение. В ра-
боте [3] предложен следующий численный алгоритм вычисления конечномерной реализации:
Алгоритм 2.1
Шаг 1. Выбор числа г, удовлетворяющего условию (1.2), эквивалентен выбору пары чисел (д ' ,д), таких, что выполняется условие:
гапкВч ,д (/ )= гапкВд,+1? (/ )= гапкВд ,д+1 (/ ).
Шаг 2. С помощью элементарных преобразований над строками и столбцами (метода Гаусса) ганкелева матрица Вд,9 (/) приводится к виду в правой части соотношения (2.1). При этом преобразованиям над строками матрицы Вд 'д (/) будет соответствовать матрица Р, а преобразованиям над столбцами - матрица М.
Шаг 3. Вычисляются матрицы г, G,H по приведенным в теореме формулам (2.2). Умножение матриц Р, В и М происходит по правилам матричного умножения, а умножению на матрицы епрг,е:г, Етг,Еррг эквивалентно выделение подматриц соответствующего размера.
Для реализации первого шага описанного выше алгоритма из ганкелевой матрицы выделяется блок Вд 'д (/), содержащий максимально возможное количество строк и столбцов. Можно показать, что если количество исходных мат-
до р "К '
-----, а д = д0 - д. За-
риц А. равно ч0, то ч = тем вычисляются ранги матриц вч ,д (/), вч ,+ц (/) и вд 'ч _ц(/) с помощью метода Гаусса приведения матрицы к треугольному виду, и если ранги всех трех матриц совпадут, то это означает, что существует реализация первых ч + ч' членов импульсной последовательности. В противном случае вычисляется одна из частичных реализаций размерности шт(р,ч 'ш).
3. Особенности машинной арифметики
Большинство проблем, связанных с машин-
ной реализацией алгоритма Б.Л. Хо, возника-
ют на первом и втором этапах, поскольку в про-
цессе диагонализации матрицы с помощью ал-
горитма Гаусса и его модификаций требуется
производить большое количество делений, ухуд-
шающих точность вычислений.
Общеизвестно [4], что в вычислительной
машине может быть представлено лишь конеч-
ное множество чисел, каждое из которых имеет
конечное число разрядов. Чаще всего они за-
писываются в полулогарифмической форме, а точнее - в форме с плавающей точкой: х = ш ■ Ье, где т - мантисса, Ь - основание степени, е -
порядок. Как правило, для внутримашинного представления выбирается основание b, равное 2, а мантисса нормализуется, то есть ее абсолютное значение помещается в интервал
2,2 . Согласно стандарту IEEE (Institute of lectrical and Electronics Engineers) FloatingPoint Standard числа с плавающей точкой представляются в памяти 32 двоичными разрядами, из них 24 бита отводится для мантиссы и 8 битов для показателя. Это означает, что наибольшее число с плавающей точкой равно 1.11к12 х10127 = 1038, наименьшее число с плавающей точкой приблизительно равно 10-38, а мантисса из 24 битов позволяет хранить примерно 7 значащих десятичных разрядов. Кроме того, даже точные исходные данные, заданные в десятичной системе счисления, могут искажаться при переводе в двоичную систему счисления для хранения в компьютере и выполнения различных операций. Например, числа 1/3, 0,1 не имеют точного представления в двоичной системе счисления. При выполнении многочисленных арифметических операций малые погрешности исходных данных могут катастрофически увеличиваться.
Также особенностью выполнения машинных арифметических операций является то, что при сложении машинных чисел различной величины результат может оказаться точно равен одному из слагаемых. Например, если выполнить над двумя числами, имеющими 7 значащих разрядов, действие 1.234 х1015 + 2.345 Х10-10, то получится 1.234х1015. Наименьшее число с плавающей точкой, которое при сложении с числом 1,0 дает результат, больший чем 1,0, называется машинным эпсилоном [4] и обозначается е маш. Машинный эпсилон определяет относительную погрешность арифметики компьютера: если x и у - два положительных числа с плавающей точкой и x > у , то их сумму можно
і+У
. Очевидно, что
записать в виде х+у = х
у ' и
при — < є маш сумма с плавающей точкой чисел хи у совпадет с х. Для 32-битовой арифметики с плавающей точкой, удовлетворяющей стандарту ІЕЕЕ, ємаш = 2-22 = 2.38х 10-7 при использовании округлений. Поэтому бессмысленно рассчитывать более чем на семь верных десятичных знаков в любом результате таких вычислений с плавающей точкой.
4. Использование интервального анализа
Поскольку при численной реализации алгоритма Б.Л. Хо происходит вычисление ранга матрицы, для которого требуется сравнивать результаты вычислений с нулем, неправильное округление результата может привести не только к искажению матриц реализации F, о, н, но и к неправильному определению размерности пространства состояний. Причиной неправильного определения ранга ганкелевой матрицы является то, что элементы, которые должны быть равны нулю при точных вычислениях, при машинных вычислениях ограниченной точности могут оказаться малыми, но не нулевыми элементами. Для решения этой проблемы можно ввести достаточно малую константу и все элементы матрицы, меньшие данной константы, считать нулями. Но чем больше вычислений будет производиться, тем больше такие элементы будут отличаться от нуля. Поэтому одно и то же малое число может на начальных этапах процедуры диагонализации соответствовать ненулевому элементу, а на заключительных этапах нулю.
Для того чтобы динамически следить за увеличением погрешности результата вычислений и оценивать ее накопление, можно использовать методы интервальной арифметики [5]. Для этого все исходные данные (элементы ган-келевой матрицы) заключаются в интервалы, имеющие своими границами машинные числа.
Обозначим множество машинных чисел через ям . Тогда любому действительному чис-
лу х из интервала
жіпу,таху У€&м У€&м
для хранения в
ЭВМ ставится в соответствие некоторое машинное число ~ е ям с помощью некоторого отображения ~ = /1(х). Это отображение называется округлением, и для него выполняется свойство монотонности: х < у ^ /1(х) < /I (у).
Для замены действительных чисел интервалами используется направленное округление [5]. Для округления вниз (4 х ) справедлива импликация: х е я ^4 х < х, а округление вверх ( Т х ) определяется следующим образом: Т х = -(4- (-х)) для любого действительного числа х, причем 4 х е ям и Т х е ям. Округление /I связано с направленным округлением следующим двойным неравенством: 4 х < /1(х) <Т х. Тогда переход от
интервала X = [х1 ,х2 ] к его машинному представлению X = [Зс1 ,~2 ] происходит с помощью направленного округления интервала ( Ь X ) по правилу: X =Ь X =Ь [х1, х2 ] = [4 х1, Т х2 ] .Соответственно при реализации алгоритма Б.Л. Хо каждый элемент х исходной последовательности импульсных матриц можно заменить на интервал [ х,Т х].
Для нахождения ранга ганкелевой матрицы отображения вход-выход, используя алгоритм Гаусса, необходимо выполнять с элементами матрицы все четыре арифметические операции. Если над двумя машинными числами х и у из ям производится машинная операция *, где * е {+,- ■,} ,то ее результатом оказывается число, которое может не принадлежать ям, но оно будет записано в память ЭВМ как число 7 е ям, полученное округлением: г = /1(х * у), при этом справедливо [5], что г = /I(х * у)е г =Ь ([х,х]*[у,у]). Аналогично определяется результат машинной операции над интервалами.
В интервальном анализе доказано [5], что если А и В - машинные интервалы, то результат операции а * в с с =Ь (а * в), и если а е А,Ь е в , то а * Ь е с =Ь (а * в). Таким образом, если все вычисления, производимые над элементами ганкеле-вой матрицы, производить над соответствующими интервалами, то точные результаты вычислений всегда будут входить в полученные интервалы. Поэтому элементам, которые при точных вычислениях должны были обнулить-ся, будут соответствовать интервалы, содержащие нуль, то есть интервалы, нижняя и верхняя границы которых отличаются знаком. Следовательно, замена нульсодержащих интервалов нулями на каждой итерации приведения матрицы к диагональному виду позволяет уменьшить погрешность вычисления ранга ганкеле-вой матрицы.
Данный метод дает необходимое условие для идентификации элемента как нуля, то есть если элемент будет нулевым, то ему обязательно будет соответствовать интервал, содержащий нуль. Однако достаточным это условие не является, так как длины интервалов постоянно растут и интервал может содержать нуль за счет своего расширения, когда соответствующее ему число не нуль. Причем ширина интервала возрастает в большей степени, чем это обусловлено исходными данными. Таким образом, в ре-
зультирующие интервалы будут входить числа, которые не могут быть получены из чисел, принадлежащих исходным интервалам. Причина этого эффекта в том, что определение интервальных операций предполагает предысторию получения операндов неизвестной. Например, по определению интервального частного имеем [2,4]: [2,4]= [0.5,2]. Хотя если эти интервалы соответствуют одному и тому же значению х = [2,4], то частное х: х равно 1 и результирующий интервал будет [1,1]. В связи с тем, что ганкелева матрица является блочно симметричной и элементы в ней повторяются, применение стандартного метода диагонализации Гаусса неизбежно приводит к расширению интервалов «лишними» значениями. Рассмотрим два примера.
Пример 1. Рассмотрим отображение вход-выход скалярной системы (с одним входом и одним выходом), заданное последовательностью чисел (матриц размера 1x1) { ,^2,а^,^4,^5 }. Соответствующая ганкелева матрица будет иметь вид:
г \
а #2
В3 3 = #2 #4
0-3 @4 ас
345
После первой итерации метода Гаусса получим матрицу:
0
0 a4 - a5
a> a2 a3 a>
a> a;a3 a>
Предположим, что с элементом а2 связан интервал [-1,2] (с шириной, равной 3), тогда операция интервального умножения а2 ■ а2 даст результат [-1,2] ■ [-1,2] = [-2,4] (с шириной, равной 6), хотя на самом деле множеством значений данного произведения является интервал [-1,2] = [0,4] (с шириной, равной 4). Таким образом, мы получили, что вместо естественного расширения интервала с ширины 3 до ширины 4 интервал пополняется «лишними» значениями от -2 до 0 и его ширина увеличивается до 6. Аналогично расширится интервал для произведения а3 ■ а3, расположенного в последнем элементе матрицы. На следующем шаге вторую строку необходимо умножить на значение
a4
a2 aз a>
aз '
2
a2
a>
которое можно преобразовать к следующему
виду:
a>a4 — a2 a;
2"~
a>a;— a2
Пусть с элементами а1 ,а2,а3,а4 связаны соответственно интервалы [2,3] [1,2] [3,4] [4,5], тогда, выполняя последовательно операции первого из представленных выражений, получим интервал [0,4] , а вычисляя значения интервальных операций по второму выражению, получим интервал [0,6], то есть второй интервал содержит значения, которые не могут быть получены при подстановке в него чисел, принадлежащих заданным интервалам.
Пример 2. Рассмотрим импульсную последовательность системы с одним входом и одним
выходом из 22 матриц Л1 =1:
1,
111111111 1 1
2 3 4 5 6 7 8 9 10 11 12
.1 2. — — — — .! — — .1
13 ’ 14’ 15’ 16’ 17’ 18’ 19’ 20’ 21’ 22.
Соответствующая ганкелева матрица отображения вход-выход будет иметь вид:
B =
1
1
1
111
12
Из
114
^Xl Х2 Хз /14 ••• Yll
V у
Матрица B является матрицей Гильберта и поэтому имеет полный ранг, то есть ранг матрицы B должен быть равен 11.
Для вычисления ранга матрицы B на ЭВМ все ее элементы представим действительными числами с 14 знаками после запятой (используя тип double двойной точности), а машинную погрешность (начальную ширину интервалов) будем считать равной 10-14. Применяя к данной матрице В и к соответствующей ей матрице интервалов алгоритм Гаусса, на последней итерации должны получить элемент
1
B
11,11
716830370256
(1.39503017937 10-12 X
который не равен нулю. Но полученный соответствующий ему интервал
[-1.47622562214 10-12 ! 8.89842221169 10-12 ] кроме истинного значения содержит нуль, поэтому элемент %п будет ошибочно считаться нулем. Ошибка происходит из-за увеличения ширины интервала с 10-14 до 10-12. Данная ошибка приводит к тому, что мы находим конечномерную реализацию, размерность которой меньше минимально возможной, то есть такой реализации для данной системы не существует.
Заключение
В работе рассмотрен численный алгоритм построения реализации линейной динамической
системы с дискретным временем над полем действительных чисел, основанный на алгоритме Б.Л. Хо. Предложен вариант решения проблемы его машинной разработки, связанной с неточностью представления данных в ЭВМ. Показано, что использование интервальной арифметики позволяет точнее выявлять нулевые элементы матрицы для правильного определения ранга ганкелевой матрицы отображения вход-выход и, соответственно, размерности пространства состояний. Для большей эффективности применения интервальной арифметики требуется модификация стандартных методов определения ранга матрицы с учетом свойств ганке-левой матрицы отображения вход-выход.
4.05.2010
Список использованной литературы:
1. Калман, Р. Очерки по математической теории систем / Р. Калман, П. Фалб, М. Арбиб. - М.: Едиториал УРСС, 2004. - 400 с.
2. Пушков, С.Г. О вычислении конечномерной реализации / С.Г. Пушков // Кибернетика и системный анализ. - 1991, №6. - С. 107-112.
3. Пушков, С.Г. Представление динамических систем в пространстве состояний: точная и приближенная реализация / С.Г. Пушков. - Барнаул: Изд-во АлтГТУ, 2003. - 270 с.
4. Каханер, Д. Численные методы и математическое обеспечение / Д. Каханер, К. Моулер, С. Нэш. -М.: Мир, 1998. - 575 с., ил.
5. Алефельд, Г. Введение в интервальные вычисления / Г. Алефельд, Ю. Херцбергер. - М.: Мир, 1987. - 360 с.
Cведения об авторах:
Пушков Cергей Григорьевич, заведующий кафедрой математической кибернетики Оренбургского государственного университета, доктор технических наук, профессор 460018, г. Оренбург, пр-т Победы, 13, тел. (3532)372535, e-mail: psg@mail.osu.ru
Горелик Анна Александровна, преподаватель кафедры администрирования информационных систем Оренбургского государственного университета, Оренбургского государственного университета, доктор технических наук, профессор 460018, г. Оренбург, пр-т Победы, 13, ауд. 1502, тел. (3532)372539, e-mail: anna_gmn3@rambler.ru
Pushkov S.G., Gorelik A.A.
The use of methods of interval analysis for enumerating the dimensionality of the finite-dimensional realization of linear dynamic system
The article examined the task of constructing the model in the state space on the basis of the data about the behavior the input-output of linear stationary dynamic system with the discrete time. The authors investigated the problem of obtaining the numerical algorithms of the calculation of finite-dimensional realization in the state space of mapping input-output, represented by the pulse sequence of matrices. They proposed the method of enumerating the dimensionality of this realization, based on the use of interval analysis.
The key words: dynamic systems, state space, the problem of realization, numerical methods, interval analysis.
Bibliography:
1. Kalman, R. Essays on the mathematical theory of systems / R. Kalman, P Falb, MA Arbib. - M.: Editorial URSS, 2004. -400 pp.
2. Pushkov, SG On the computation of finite realization / SG GUN / / Cybernetics and system analysis. - 1991 № 6. - S. 107-112.
3. Pushkov, SG Representation of dynamic systems in state space: an exact and approximate realization / SG Pushkov. -Barnaul: Altai State Technical University Publishing House, 2003. - 270 pp.
4. Kahaner, D. Numerical methods and software / D. Kahaner, C. Mowler, C. Nash. -M.: Mir, 1998. - 575 pp., Il.
5. Alefeld, G. Introduction to interval computation / G. Alefeld, J. Hertsberger. - M.: Mir, 1987. - 360 pp.
