Использование методов интервального анализа для учета погрешности вычисления конечномерных реализаций линейных динамических систем

Горелик А.А., Пушков С.Г.

Решение многих прикладных задач математики сводится к построению математической модели некоторого процесса или явления. Зачастую моделирование основывается на результатах наблюдения за поведением некоторой динамической системы, и сводится к представлению экспериментальных данных в форме, удобной для дальнейшего исследования. Наблюдения за поведением системы позволяют задать ее с помощью отображения «вход»-«выход». В этом случае для определения реакции системы необходимо знать и использовать всю предысторию входов и выходов, что не всегда удобно осуществить на практике. Поэтому часто ставится задача задания динамической системы, используя ее свойства (пространство состояний), таким образом, чтобы реакция системы определялась через ее состояние и входной сигнал в данный момент времени. Для линейных динамических стационарных систем с дискретным временем эта реакция определяется по следующим формулам:

где – момент времени; – векторы состояния системы в моменты времени и ; – вектор входного сигнала и – соответствующий ему вектор выходного сигнала; – линейные отображения.

Проблема построения модели с пространством состояний на основе данных «вход»-«выход» в теории систем известна как задача реализации [1].

Математическая формулировка этой задачи для класса линейных динамических систем с дискретным временем состоит в следующем:

для заданной последовательности матриц размера pxm над полем K

построить тройку матриц (F, G, H) над тем же полем K таких, что

где для некоторого n (которое также нужно определить) H – матрица размера pxn, F – матрица размера nxn, G – матрица размера nxm.

Одним из известных алгоритмов вычисления конечномерных реализаций для систем над полями является алгоритм Б.Л.Хо [1]. Численная реализация данного алгоритма [2,3] основана на приведении ганкелевой матрицы отображения вход-выход системы к диагональному виду, которое можно осуществить с помощью метода Гаусса и его модификаций. Однако при реализации на ЭВМ алгоритма диагонализации матрицы возникают проблемы, связанные с погрешностью машинных вычислений. Поскольку при вычислении ранга матрицы требуется сравнивать результаты вычислений с нулем, то неправильное округление результата может привести не только к искажению матриц реализации F, G, H, но и к неправильному определению размерности пространства состояний.

Для решения этой проблемы можно завести достаточно малую константу, и все элементы матрицы, меньшие нее, считать нулями. Но чем больше вычислений будет производиться, тем больше такие элементы будут отличаться от нуля. Поэтому одно и то же малое число может на начальных этапах метода Гаусса соответствовать ненулевому элементу, а на заключительных этапах нулю.

Для того, чтобы динамически следить за увеличением погрешности результата вычислений, можно использовать методы интервальной арифметики [3]. Для этого все исходные данные (элементы ганкелевой матрицы) представляются интервалами с радиусом, равным начальной машинной погрешности. Далее все вычисления, производимые над элементами ганкелевой матрицы, производятся над соответствующими интервалами. В результате вычислений радиусы интервалов динамически растут. Элемент ганкелевой матрицы считается равным нулю, если соответствующий ему интервал содержит нуль.

Данный метод дает необходимое условие для идентификации элемента как нуля, то есть если элемент будет нулевым, то ему обязательно будет соответствовать интервал, содержащий нуль. Однако достаточным это условие не является, так как длины интервалов постоянно растут, и интервал может содержать нуль за счет расширения интервала, когда соответствующее ему число не нуль.

Для примера рассмотрим импульсную последовательность системы с одним входом и одним выходом из 22 матриц :

Соответствующая ганкелева матрица отображения вход-выход будет иметь вид:

Матрица B является матрицей Гильберта и поэтому имеет полный ранг, то есть ранг матрицы B должен быть равен 11.

Для вычисления ранга матрицы B на ЭВМ все ее элементы представим действительными числами с 14 знаками после запятой, а машинную погрешность (начальную ширину интервалов) будем считать равной 10^-16. Применяя к данной матрице B и к соответствующей ей матрице интервалов алгоритм Гаусса, на последней итерации должен получиться элемент (), который не равен нулю. Но полученный соответствующий ему интервал [-;] кроме истинного значения содержит нуль, поэтому элемент будет ошибочно считаться нулем. Ошибка происходит из-за увеличения ширины интервала с 10^-16 до 10^-12. Данная ошибка приводит к тому, что мы находим конечномерную реализацию, которой в данной системе быть не может.

Использование интервальной арифметики при построении на ЭВМ реализации линейной динамической системы с дискретным временем над полем действительных чисел позволяет точнее выявлять нулевые элементы матрицы за счет повышения трудоемкости вычислений.

1. 
   Калман, Р. Очерки по математической теории систем / Р. Калман, П. Фалб, М.Арбиб. – М.: Едиториал УРСС, 2004. – 400 с.
   
2. 
   Пушков, С.Г. О вычислении конечномерной реализации / С.Г. Пушков // Кибернетика и системный анализ. – 1991, №6. – С.107-112.
   
3. 
   Пушков, С.Г. Представление динамических систем в пространстве состояний: точная и приближенная реализация / С.Г. Пушков. – Барнаул: Изд-во АлтГТУ, 2003. – 270 с.
   
4. 
   Алефельд, Г. Введение в интервальны вычисления / Г. Алефельд, Ю. Херцбергер. – М.: Мир, 1987. – 360 с.