Известия РАН. Теория и системы управления, 2021, № 2, стр. 25-34

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

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

Для класса линейных стационарных динамических систем с дискретным временем методы оценивания параметров импульсной характеристики могут оказаться также полезными при построении моделей с пространством состояний на основе данных об измерениях входных и выходных сигналов системы. В частности, эти методы оказываются необходимой частью при осуществлении двухэтапной процедуры построения модели с пространством состояний [5]. Здесь задачей первого этапа является получение отображения вход-выход, хорошим представлением которого является импульсная характеристика, для систем с дискретным временем это есть импульсная последовательность матриц. На втором этапе из решения задачи реализации вычисляются матрицы эквивалентной системы в виде модели с пространством состояний. Некоторые методы решения задачи второго этапа для различных подходов к представлению параметрической неопределенности можно найти в следующих работах: [6] – при статистическом подходе, [7] – для интервальных систем, [8] – для систем над нечеткими числами.

Настоящая статья посвящена задаче оценивания параметров импульсной характеристики на основе данных об измерении входных и выходных сигналов системы. С этой целью после краткого анализа проблемы оценивания импульсной характеристики для класса линейных стационарных динамических систем с дискретным временем рассматриваются постановки задач и методы их решения для различных подходов к описанию неопределенностей в данных. Наряду с методом, основанным на минимизации среднеквадратичного отклонения (СКО), представлены варианты методов интервального и нечеткого оценивания матриц импульсной последовательности. Работа завершается сравнительным анализом и численной иллюстрацией методов.

1. Импульсная характеристика динамической системы и ее оценивание. Рассмотрим линейную стационарную динамическую систему $\Sigma $ с дискретным временем с $m$ входами $n$ состояниями и p выходами над полем действительных чисел $\mathbb{R}$, которую можно отождествить с четверкой линейных отображений $\Sigma = (F,G,H,J)$:

где $U = {{\mathbb{R}}^{m}}$, $X = {{\mathbb{R}}^{n}}$ и $Y = {{\mathbb{R}}^{p}}$ – векторные пространства над полем $\mathbb{R}$ размерностей $m$, $n$ и $p$ (пространства входных сигналов, состояний и выходных сигналов соответственно). Динамическое поведение системы $\Sigma $ описывается разностными уравнениями: где $x(t) \in X$, $u(t) \in U$, $y(t) \in Y$ для всех дискретных моментов времени t. Отображения $F$, $G$, $H$ и J в уравнениях (1.1) могут быть представлены матрицами с вещественными коэффициентами размеров $n \times n$, $n \times m$, $p \times n$ и $p \times m$ соответственно. Поэтому далее мы будем обозначать эти отображения и соответствующие им матрицы одинаковыми символами.

Пусть на вход системы $\Sigma = (F,G,H,J)$ с динамическим поведением (1.1), находящейся в начальный момент времени $t = 0$ в состоянии x(0), поступает последовательность входных сигналов $u(0),u(1),u(2),...$ Тогда по индукции легко установить, что значения состояния и выходного сигнала в любой последующий момент времени $t$ можно вычислить по формулам

Если $x(0) = 0$, то $y(0) = 0$ (в правых частях последних двух уравнений $0$ следует понимать в векторном смысле) и, согласно (1.2),

Вводя обозначения

получаем уравнения

Таким образом поведение вход-выход линейной стационарной динамической системы с дискретным временем, находящейся в начальный момент времени t = 0 в нулевом состоянии, $x(0) = 0$, наряду с отображением вход-выход ${{f}_{\Sigma }}:\left( {u(0),u(1),...} \right) \mapsto \left( {y(0),y(1),...} \right)$ [1] может быть охарактеризовано импульсной последовательностью матриц (называемой также последовательностью весовых матриц или последовательностью марковских параметров)

размера p × m, являющейся представлением импульсной характеристики системы. Это название связано с тем, что данную последовательность можно получить в результате специального эксперимента с системой, подав на ее вход единичный импульс и наблюдая значения выходного сигнала. Уравнения (1.3) иногда называют уравнениями свертки, а поскольку матрицы (1.4) определяют веса, с которыми входные сигналы в предшествующие моменты времени представлены в выходном сигнале, становится понятным другое название этой последовательности.

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

Линейная стационарная динамическая система с дискретным временем $\Sigma = (F,G,H,J)$ называется асимптотически устойчивой, если

для некоторой матричной нормы ${\text{||}}\, \cdot \,{\text{||}}$. Легко установить, что если система $\Sigma = (F,G,H,J)$ устойчива, то для ее матриц импульсной последовательности будет выполняться условие и тогда для некоторого ограниченного τ, называемого временем затухания переходных процессов, вместо (1.3) мы можем записать приближенное равенство

Более того, поскольку в этом случае $\mathop {\lim }\limits_{j \to \infty } {\text{||}}H{{F}^{j}}x(0){\text{||}} = 0$ для любого начального состояния $x(0)$, то мы можем использовать уравнения (1.5) и для систем с ненулевым начальным состоянием.

С учетом вышеизложенного задача оценивания импульсной характеристики асимптотически устойчивой линейной стационарной динамической систем с дискретным временем приобретает следующий вид. Пусть N – число экспериментов – пар $\left( {u(i),y(i)} \right)$, $i = \overline {1,N} $, $N > \tau $. Будем определять A_j для $j = \overline {0,\tau } $ из решения задачи минимизации некоторого функционала:

который характеризует меру несогласованности (ошибки) между моделью (1.5) и наблюдаемыми экспериментальными данными вход-выход. Конкретный вид функционала $\Phi $ в значительной степени определяет метод решения этой задачи.

2. Статистический подход. Наиболее изучена методология оценивания параметров на основе минимизации СКО с использованием статистического подхода к проблеме [2, 9]. Необходимые оценки и свойства этих оценок можно получить с помощью статистических понятий и методов.

Будем определять A_j для $j = \overline {0,\tau } $ из решения задачи минимизации квадратичного функционала:

где ${\text{||}}\, \cdot \,{\text{||}}$ – матричная евклидова норма.

Введем в рассмотрение блочные матрицы

где S_jt – матрица размера $m \times m$, ${{R}_{t}}$ – матрица размера $p \times m$,

В этих обозначениях исходная задача сводится к задаче вычисления параметров множественной линейной регрессии.

Предложение. Для асимптотически устойчивых линейных стационарных динамических систем с дискретным временем решение задачи (2.1) удовлетворяет матричному уравнению

Вариант доказательства этого утверждения без обращения к статистическим понятиям можно найти в работе [5].

Таким образом, оценки матриц импульсной последовательности асимптотически устойчивой системы мы можем получить, решив матричное уравнение (2.5). В результате решения этой задачи оценки будут точечными. Если мы обладаем достаточной информацией о законе распределения ошибок в исходных данных, то на основе вероятностно-статистических методов можно получить информацию о свойствах этой оценки, в том числе найти оценки доверительных интервалов.

3. Интервальное оценивание. При данном подходе мы изначально ищем не точечные оценки, а интервальные. Причем эти интервальные оценки следует интерпретировать в ином по сравнению со статистическим подходом смысле, не как доверительные интервалы, а как множество равноправных точечных моделей, которые объясняют наблюдаемые данные.

Итак, пусть по-прежнему мы имеем в распоряжении данные измерений входного и выходного сигнала некоторой системы:

а также модель, описывающую поведение вход-выход:

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

То, что матрицы A_j являются интервальными, в частности, означает, что все элементы матриц A_j – интервалы. Эти матрицы мы будем представлять в виде

где ${{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{A} }_{j}}$ и ${{\bar {A}}_{j}}$ – точечные матрицы, элементы которых состоят из нижних и верхних границ элементов матрицы A_j  соответственно. Уравнения (3.1) с интервально заданными матрицами A_j, $j = \overline {0,\tau } $, следует понимать как семейство уравнений с точечными матрицами ${{B}_{j}} \in \left[ {{{{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{A} }}_{j}},{{{\bar {A}}}_{j}}} \right]$, $j = \overline {0,\tau } $.

Для решения задачи оценивания матриц импульсной характеристики в виде (3.2) образуем допустимое множество, состоящее из всех наборов точечных матриц $\left\{ {{{B}_{0}},...,{{B}_{\tau }}} \right\}$, удовлетворяющих уравнениям (3.3):

Таким образом требование соответствия модели (1.5) наблюдаемым данным вход-выход, которое ранее выражалось в задаче минимизации (2.1) функционала Φ, теперь состоит в принадлежности допустимому множеству ${\mathbf{B}}$. Сами интервальные оценки будем находить из решения $2(\tau + 1)pm$ задач линейного программирования на нахождение нижних и верхних границ всех элементов матриц A_j:

где ${{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{a} }_{{kql}}}$, ${{\bar {a}}_{{kql}}}$, ${{b}_{{kql}}}$ – элементы, расположенные в q-й строке и l-м столбце, матриц ${{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{A} }_{k}}$, ${{\bar {A}}_{k}}$, ${{B}_{k}}$ соответственно. В каждой из этих $2(\tau + 1)pm$ задач на допустимом множестве B, определенном, согласно (3.4), в качестве функционала, используется один из $(\tau + 1)pm$ элементов матриц ${{B}_{j}}$ (остальные свободны и могут принимать любые значения в рамках ограничивающего множества B), для которых находится минимум и максимум. Они же являются и искомыми параметрами.

Задачи (3.5) на допустимом множестве ${\mathbf{B}}$, определенном согласно (3.4), в большинстве случаев не имеют решения, поскольку несовместной оказывается система ограничений. Эти задачи гарантировано имеют решения только для точных данных и точного указания порядка модели. Поэтому при решении задач (3.5) вместо допустимого множества, определенного согласно (3.4), целесообразно рассматривать более широкое множество, учитывающее погрешность определения выхода $\Delta y(t)$. При этом ограничения принимают форму неравенств:

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

Здесь $\mathbb{I}\mathbb{R}$ – множество всех вещественных интервалов, соответственно $\mathbb{I}{{\mathbb{R}}^{p}}$ – множество всех интервальных векторов размерности p. Через $\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{y} (t)$ и $\bar {y}(t)$ будем обозначать векторы нижних и верхних границ вектора выхода соответственно. В итоге получается задача (точнее, семейство задач) (3.5) на допустимом множестве

Для случая

с неотрицательными интервальными входами $u(t) = \left[ {\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{u} (t),\bar {u}(t)} \right]$ можно рассматривать задачи (3.5) на допустимом множестве

Наконец, возможно применение интервального варианта метода, описанного в разд. 2. Тогда потребуется решение матричного уравнения

с интервальной матрицей $R$ и, в общем случае, даже с интервальной матрицей $S$. Здесь при формировании интервальных матриц $R$ и $S$ в (2.2) вместо уравнений (2.3) и (2.4) нужно использовать их интервальные аналоги:

Последний случай с интервальной матрицей $S$ может привести к весьма сложной и не всегда полезной задаче (из-за очень широких оценочных интервалов), поскольку потребует обращения интервальной матрицы. В то время как первый случай (только с интервальной матрицей $R$ и точечной $S$) всегда имеет решение (когда $S$ является невырожденной) с адекватной шириной оценочных интервалов. Причем когда $\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{y} (t) = y(t) - \Delta y(t)$ и $\bar {y}(t) = y(t) + \Delta y(t)$, точечное решение, полученное на основе метода, описанного в разд. 2, оказывается в точности центром решения интервальной задачи.

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

В отличие от подхода разд. 3, где интервал представлялся нижней и верхней границами, $a = \left[ {\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{a} ,\bar {a}} \right]$, здесь будем описывать интервалы с использованием его центра c и ширины $\delta $, применяя запись $a = (c,\delta )$, где $\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{a} = c - \delta $, $\bar {a} = c + \delta $. Интервал $a = (c,\delta )$ будем считать нечетким числом с треугольной функцией принадлежности (треугольное нечеткое число), в то время как интервалы $a = \left[ {\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{a} ,\bar {a}} \right]$ целесообразно рассматривать как нечеткие числа с постоянной на $\left[ {\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{a} ,\bar {a}} \right]$ функцией принадлежности.

Итак, как и прежде, мы имеем в распоряжении данные измерений:

и импульсная характеристика описывается моделью но теперь A_j есть матрицы с нечеткими треугольными числами, которые мы будем представлять в виде где C_j – матрицы центров и ${{\Delta }_{j}}$ – матрицы ширин матриц A_j. Таким образом формула (4.1) для вычисления нечеткого (интервального) выхода принимает вид где $\left| {u(t - j)} \right|$ понимается как вектор, все координаты которого состоят из абсолютных значений соответствующих координат вектора $u(t - j)$:

В соответствии с формулами (4.2) суммарная ширина оценочного интервала для p-мерного вектора выхода $y(t)$ определяется как

Найдем оценки матриц A_j, которые дают минимальную суммарную величину ширины оценочных интервалов всех векторов выхода, которая в данном случае вычисляется по формуле

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

В отличие от традиционной регрессионной оценки для данного метода характерно то, что увеличение числа данных расширяет оценочный интервал.

5. Сравнительный анализ и численная иллюстрация методов. Проведем сравнительный анализ представленных выше методов оценивания матриц импульсной характеристики и проиллюстрируем их на численных примерах. Таким образом в сравнении будут участвовать следующие четыре метода оценивания импульсной характеристики:

1) метод, который будем обозначать через StaImp, основанный на минимизации квадратичной ошибки путем решения задачи (2.1), которое получается из решения матричного уравнения (2.5) (его условно можем считать статистическим);

2) интервальный метод IntImp1, базирующийся на решении $2(\tau + 1)pm$ задач линейного программирования (3.5) на допустимом множестве B, определенном согласно (3.6);

3) интервальный метод IntImp2 из решения интервального матричного уравнения (3.7) с интервальной матрицей R;

4) нечеткий (интервальный) метод FuzImp на основе решения задачи линейного программирования (4.3).

Метод StaImp является точечным в том смысле, что мы получаем точечные (неинтервальные) оценки параметров импульсной характеристики системы. Поскольку этому методу можно дать статистическую интерпретацию, то, используя те или иные предположения о законах распределения искажений выходного сигнала и соответствующие им статистические методы, удается получить интервальные оценки, например 95%- или 99%-ные доверительные интервалы. Но проблема заключается в том, что эти доверительные интервалы в большинстве реальных ситуаций оказываются настолько большими, что теряют всякий практический смысл. Достоинством данного метода является относительно небольшие вычислительные затраты. Точность оценок улучшается при увеличении длины временного ряда.

Интервальный метод IntImp1 требует значительно больших (по сравнению с методом StaImp) вычислительных затрат. Достоинство метода состоит в том, что получаемые интервальные оценки содержат все точечные решения, которые описывают наблюдаемый временной ряд с учетом заранее заданных погрешностей измерений или наблюдений. Также плюсом этого метода является то, что в систему ограничений задачи линейного программирования можно добавить ограничения, касающиеся априорных знаний об оцениваемых параметрах. В отличие от метода StaImp, ширина полученных интервальных оценок растет с увеличением длины временного ряда. Недостаток данного метода состоит также в том, что при достаточно большой длине временного ряда задача линейного программирования, лежащая в основе метода, не всегда имеет решение (из-за несовместимости системы ограничений). Это затруднение может быть преодолено путем взятия больших величин $\Delta y(t)$ в задачах (3.6), что в свою очередь приведет к увеличению ширины оценочных интервалов.

Интервальный метод IntImp2 дает в итоге более широкие интервалы в качестве оценок параметров, и эти оценки являются практически гарантированными. Решение задачи существует всегда, когда невырожденной оказывается матрица $S$ в уравнении (3.7). Вычислительные затраты на реализацию этого метода будут существенно меньшими, чем при реализации интервального метода IntImp1. Когда $\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{y} (t) = y(t) - \Delta y(t)$ и $\bar {y}(t) = y(t) + \Delta y(t)$, центры интервальных оценок, полученных методом IntImp2, в точности совпадают с точечными оценками, полученными методом StaImp. Недостаток метода – часто избыточно большая ширина полученных интервалов.

Нечеткий (интервальный) метод FuzImp требует меньших по сравнению с методом IntImp1 вычислительных затрат, поскольку решается только одна задача линейного программирования. Полученные оценочные интервалы являются несколько более широкими, чем при использовании метода IntImp1. В отличие от метода IntImp1, задача линейного программирования метода FuzImp практически всегда имеет решение. Увеличение числа данных (длины временного ряда наблюдения) расширяет оценочные интервалы. Данную ситуацию, как и в случае с методом IntImp1, можно интерпретировать как приобретение новой информации, что приводит к расширению возможной оценки. Метод FuzImp имеет тот же плюс, что и метод IntImp1, а именно возможность дополнять систему ограничений априорными знаниями специалистов об оцениваемых параметрах.

Пример 1. Рассмотрим линейную стационарную систему $\Sigma = (F,G,H,J)$ с одним входом, одним выходом и тремя состояниями ($m = p = 1$, n = 3):

Импульсная последовательность для этой системы имеет вид

Оценивание импульсной последовательности будем проводить на основе временного ряда

порождаемого уравнениями где в качестве входного сигнала $u(t)$ используется равномерно распределенная случайная величина из интервала [–1, 1], а выходной сигнал $y(t)$ искажен равномерно распределенным на интервале $[ - 0.1,\;0.1]$ или нормально распределенным случайным шумом $\eta (t)$ с нулевым средним и дисперсией, подобранной так, чтобы соотношение сигнал/шум было близким к 10.

Для достаточно короткого временного ряда длины N = 30, искажения выходного сигнала равномерным на интервале $[ - 0.1,\;0.1]$ шумом и временем затухания переходных процессов τ = 3 результаты оценивания импульсной характеристики по описанным выше четырем методам приведены в табл. 1. Здесь для метода StaImp при вычислении средней ширины оценочных интервалов использовались 95%-ные доверительные интервалы, вычисленные на основе соответствующих процентилей распределения Стьюдента. СКО в интервальных методах вычислялось для центров полученных интервальных оценок суммарно по всем четырем параметрам. Для более длинного временного ряда длины $N = 100$ при этих же условиях результаты оценивания импульсной характеристики дают более высокую точность оценивания по всем методам (если сравнивать СКО центров оценок), при этом ширина оценочных интервалов, полученных методом FuzImp, существенно возрастет.

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

Пример 2. Теперь рассмотрим линейную стационарную систему с двумя выходами Σ = = $(F,G,H,J)$ ($m = 1$, $p = 2$, $n = 4$):

В отличие от системы примера 1 данная система хотя и является асимптотически устойчивой, матрицы ее импульсной последовательности не превращаются в нулевые, начиная с четвертого члена.

Оценивание матриц импульсной последовательности также будем проводить на основе временного ряда (5.1), порождаемого уравнениями (5.2), где выходной сигнал $y(t)$ искажен равномерно и нормально распределенным векторным случайным шумом $\eta (t)$ с теми же параметрами распределений, что и в примере 1.

Результаты оценивания импульсной характеристики на основе временного ряда длины N = 30 с равномерными искажениями (при соотношении сигнал/шум около 10) приведены в табл. 2 (τ = 3).

Вообще при небольшом уровне зашумления выходного сигнала метод IntImp1 довольно часто может не давать результата либо полученные оценочные интервалы окажутся настолько узкими, что могут не содержать точные значения. Но по мере уменьшения соотношения сигнал/шум картина существенно меняется (особенно для равномерно распределенного шума). При высоком уровне зашумления выходного сигнала этот метод может стать наиболее предпочтительным по сравнению с другими.

Заключение. Рассмотрены методы оценивания матриц импульсной характеристики линейных стационарных динамических систем с дискретным временем, основанные на различных подходах к представлению параметрической неопределенности. Наряду со стандартным подходом, допускающим вероятностно-статистическую интерпретацию и основанным на минимизации СКО, рассматриваются подходы в рамках описания неопределенностей на языках интервального анализа и теории нечетких множеств.

В работе представлены четыре метода оценивания матриц импульсной характеристики: метод StaImp, основанный на минимизации квадратичной ошибки; интервальный метод IntImp1, базирующийся на вычислении минимальных интервалов, которые содержат наблюдаемые данные; интервальный метод IntImp2 можно считать обобщением метода StaImp на интервальный случай; нечеткий (интервальный) метод FuzImp на основе решения задачи минимизации области неопределенности (метод нечеткой регрессии).

Представленные методы дают разные по качеству оценки и имеют различную вычислительную сложность. Метод StaImp имеет преимущества по сравнению с другими методами, когда в нашем распоряжении имеются достаточно длинные временные ряды наблюдений, а характер шума, искажающего выходные сигналы, близок к нормальному распределению. Для искажений случайным шумом без тенденции к центрированию (в особенности, равномерном распределении) наиболее предпочтительным становятся интервальные методы IntImp1 и IntImp2. Метод IntImp1 имеет преимущества для коротких временных рядов и высоких уровней “зашумления” выходного сигнала, а метод IntImp2 дает лучшие результаты для более длинных временных рядов и низких уровней зашумления. Метод FuzImp может с успехом применяться для случаев искажения шумами с тенденцией к центрированию, а также при нестационарном характере шума.