Журнал вычислительной математики и математической физики, 2019, T. 59, № 6, стр. 937-950

ВВЕДЕНИЕ

В работах [1]–[4] развита эффективная в вычислительном плане процедура обобщенного инвариантно-несмещенного оценивания (ОИНО) значений непрерывных линейных функционалов (НЛФ) от информационного сообщения (в дальнейшем эти значения будем называть числовыми характеристиками), реализация которого наблюдается в аддитивной смеси с флуктуационным шумом и сингулярной помехой. При этом сообщение и помеха рассматриваются в виде линейных комбинаций из соответствующих конечномерных функциональных пространств с заданными базисами, а спектральные коэффициенты при базисных функциях полагаются неизвестными. Помеха на интервале наблюдения может многократно превосходить (по абсолютной величине) сумму реализаций информационного сообщения и шума, т.е. речь идет о так называемых некорректных наблюдениях, с которыми мы зачастую сталкиваемся на практике [1]–[10]. В качестве значений НЛФ могут рассматриваться оценки сообщения и его производных различных порядков в фиксированных точках интервала наблюдения, определенный интеграл от данного сообщения на этом интервале, спектральные коэффициенты разложения по заданной системе функций и др. [1]–[4].

Для решения задачи оценивания значений НЛФ от информационного сообщения в условиях некорректных наблюдений можно использовать известные статистические методы оптимального оценивания [1]–[10]: максимума апостериорной плотности вероятности, максимального правдоподобия, наименьших квадратов и их многочисленные разновидности (см., например, [7]). В данных методах для учета фактора неопределенности, связанного с возможностью появления сингулярной помехи, предусматривается использование процедуры расширения пространства состояний, когда все неизвестные спектральные коэффициенты данной помехи включаются в число оцениваемых параметров. Однако при таком подходе могут существенно вырасти размерность решаемой задачи, объем вычислительных затрат и может наблюдаться известный эффект “размазывания точности” [8], [10]. Описываемые в работах [1]–[10] подходы к решению указанной задачи ориентированы на различные уровни априорной неопределенности, но при этом все факторы такой неопределенности учитываются заданием компактных множеств допустимости (например, в виде эллипсоидов) в рамках вероятностных и невероятностных критериев. Гораздо более сложную задачу представляет случай, когда наблюдения, помимо традиционной флуктуационной погрешности, содержат и сингулярную помеху, структура которой характеризуется заданным с точностью до параметров семейством возможных функциональных базисов, может принимать сколь угодно большие значения, а также иметь разрывы I рода с неизвестными моментами времени их появления. В этих условиях, когда отсутствует эффективный инструмент описания множеств допустимости, может быть применена вышеуказанная процедура ОИНО.

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

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

Вопросы, связанные с оценкой значений некоторых НЛФ для ДС, зачастую являются весьма актуальными. Например, это важно при анализе различных переходных процессов, решении задач оптимального управления, идентификации, селекции, распознавания и др. (см. [5]–[10]). Как правило, эти задачи должны решаться в реальном времени, поэтому развиваемый нами численно-аналитический метод исследования ДС, не требующий расширения пространства состояний, может быть эффективным инструментом достижения оперативности расчетов. Так, в известной задаче селекции истинных и ложных пеленгов в двухпозиционной триангуляционной системе необходимо в реальном времени проводить оценку производных (до третьего порядка включительно) по некорректным наблюдениям пеленгов излучающей цели [8]. В других задачах, связанных с синтезом управлений и оцениванием состояний ДС с учетом результатов наблюдений [5]–[7], необходимо иметь информацию не только о числовых характеристиках, но и прогнозировать траекторию движения ДС с учетом данных измерительных датчиков.

Таким образом, актуальна задача разработки заявленного выше метода исследования ДС без традиционного расширения пространства состояний, который бы удовлетворял условиям оптимальности оценивания, а именно обеспечивал бы несмещенность формируемых оценок, их инвариантность к сингулярной помехе и минимизацию дисперсии ошибок оценивания (сглаживающий эффект по отношению к флуктуационному шуму). Требуется получить удобные для вычислений явные векторно-матричные выражения, позволяющие исследовать поведение ДС с учетом возможности оценивания значений интересующих нас НЛФ, а также определить границы применимости развиваемого метода и дать полезные практические рекомендации. Необходимо также сформировать аналогичные выражения для анализа качества формируемых оценок и расчета методических ошибок, обусловленных применением конечномерных функциональных пространств, подтвердить возможность эффективного использования метода в ходе численного эксперимента на конкретном иллюстративном примере.

1. ПОСТАНОВКА ЗАДАЧИ

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

Скалярная ДС задается дифференциальным уравнением

где $x({{t}_{0}},{{x}_{0}}) = {{x}_{0}} \in \left[ {c,d} \right]$, $f(t,x)$ – функция, удовлетворяющая для всех $x = x(t,{{x}_{0}})$ из прямоугольника $G = \left\{ {{{t}_{0}} \leqslant t \leqslant {{t}_{0}} + T,\;c \leqslant {{x}_{0}} \leqslant d} \right\}$ известным условиям существования и единственности решения уравнения (1.1) с постоянной Липшица ${{L}_{0}}$ (см. [11]).

В прямоугольнике $G$ для заданных натуральных чисел ${{M}_{t}}$ и ${{M}_{x}}$ (смысл которых устанавливается в дальнейшем) справедливы ограничения:

где Данные ограничения показывают, что в дальнейшем мы оперируем с уравнениями (1.1), степень гладкости решений которых зависит от выбора значений $j$ и $k$. На практике значения констант ${{\eta }_{{0j}}}$, ${{\eta }_{{1k}}}$ и ${{\eta }_{{2jk}}}$ находятся либо в ходе численно-аналитических расчетов с использованием уравнения (1.1), либо при наблюдении за конкретной ДС в ходе различных испытаний, либо при ее эксплуатации.

Поскольку точное решение $x = x(t,{{x}_{0}})$, как правило, неизвестно, то в дальнейшем будем оперировать с некоторым приближенным решением $\tilde {x} = \tilde {x}(t,{{x}_{0}})$, которое с учетом (1.2) может быть представлено в виде

где есть вектор функций, характеризующих зависимость решения $\tilde {x}(t,{{x}_{0}})$ от начального условия ${{x}_{0}} = \tilde {x}({{t}_{0}},{{x}_{0}})$, суть векторы коэффициентов, подлежащих оцениванию, есть вектор заданных базисных функций от аргумента ${{x}_{0}}$, есть вектор заданных базисных функций от аргумента $t$.

Выбор вектора $\Psi (t)$ зачастую представляет собой достаточно непростую и трудоемкую задачу [7], при этом учитываются возможность и адекватность описания интегральных кривых уравнения (1.1), соответствие условиям наблюдения ДС и различные вычислительные аспекты [7].

Далее, для начального условия ${{x}_{0}} \in \left[ {c,d} \right]$ будет использоваться вычислительная сетка $\left\{ {{{x}_{{0\left( m \right)}}}} \right\}_{{m = 1}}^{{{{M}_{x}}}}$, на которой базисные функции ${{\lambda }_{1}}({{x}_{0}})$, ${{\lambda }_{2}}({{x}_{0}})$, …, ${{\lambda }_{{{{M}_{x}}}}}({{x}_{0}})$ мы полагаем линейно независимыми. Кроме того, аналогичное предположение относится также к временной координате $t \in \left[ {{{t}_{0}},{{t}_{0}} + T} \right]$ и базисным функциям ${{\psi }_{1}}(t)$, ${{\psi }_{2}}(t)$, …, ${{\psi }_{{{{M}_{t}}}}}(t)$, для которых в дальнейшем применяется вычислительная сетка $\left\{ {{{t}_{{\left( k \right)}}}} \right\}_{{k = 1}}^{{{{M}_{t}}}}$, т.е. при разработке метода мы будем оперировать с двумя конечномерными базисами для ${{x}_{0}}$ и $t$ соответственно. Наличие таких базисов гарантирует невырожденность и выполнение ограничений на ранги некоторых матриц, которые встречаются при изложении основных положений метода.

В области $G$ для заданных значений ${{\varepsilon }_{0}} > 0$, ${{M}_{t}}$ и ${{M}_{x}}$ на точном и приближенном решениях должно выполняться условие

Близость точного $x = x\left( {t,{{x}_{0}}} \right)$ и приближенного $\tilde {x} = \tilde {x}\left( {t,{{x}_{0}}} \right)$ решений в области $G$ будем характеризовать невязкой $\varepsilon = \left\| {x - \tilde {x}} \right\|$, соответственно для оценки близости производных вводится невязка ${{\varepsilon }^{{\left( 1 \right)}}} = \left\| {dx{\text{/}}dt - d\tilde {x}{\text{/}}dt} \right\|$ между точной ${{x}^{{\left( 1 \right)}}} = dx{\text{/}}dt$ и приближенной ${{\tilde {x}}^{{\left( 1 \right)}}} = d\tilde {x}{\text{/}}dt$ производными.

На множестве интегральных кривых $W[{{t}_{0}},{{t}_{0}} + T]$ уравнения (1.1) вводится семейство $\Omega $ функционалов $\omega $, которые характеризуют рассматриваемую ДС.

Воспользуемся дискретным уравнением наблюдения

где В этом уравнении: ${{x}_{0}}$ – неизвестное начальное условие, $h(t)$ – сингулярная помеха,$\xi (t)$ – флуктуационный шум.

Для описания помехи воспользуемся представлением

где есть вектор неизвестных коэффициентов, есть вектор заданных базисных функций, ${{M}_{h}}$ – заданное число степеней свободы.

По аналогии с [3], [4] вектор $\Theta (t)$ может формироваться из семейства возможных базисов, которые в наибольшей степени соответствуют всем возможным факторам неопределенности, соответствующим различным вероятным условиям наблюдения ДС. Такое семейство должно предусматривать возможность появления сингулярных помех самого разного типа, с которыми можно столкнуться на практике при исследовании конкретной ДС в тех или иных условиях наблюдения.

Шум $\xi (t)$ на сетке $\left\{ {{{t}_{n}}} \right\}_{{n = 0}}^{N}$ характеризуется нулевым математическим ожиданием и соответствующей корреляционной матрицей ${{K}_{\Xi }}$.

Введем векторы

Это позволяет перейти от (1.5) к векторному уравнению наблюдения

Из множества $\Omega $ выберем произвольный функционал

в другой записи: На основе наблюдения $Y$ мы должны дать оценку значения ${{\omega }_{0}}$, т.е. вычислить числовую характеристику ДС, соответствующую функционалу $\omega \in \Omega $. Данную оценку будем искать в классе линейных оценок: где есть вектор искомых весовых коэффициентов.

В силу линейности процедуры оценивания (1.8), при фиксированном $P$ дисперсия ошибки вычисления характеристики ${{\omega }_{0}}$ находится как

Выбор оптимального значения $P{\kern 1pt} *$ соответствует критерию

при этом должны выполняться условия несмещенности оценки и ее инвариантности к помехе

Требуется с учетом (1.1)–(1.12):

– сформировать алгоритм построения приближенного решения $\tilde {x}(t,{{x}_{0}})$,

– рассмотреть вопрос оптимизации выбора параметров данного алгоритма,

– исследовать методическую погрешность алгоритма, обусловленную применением усеченных базисов в модели (1.3),

– сформировать алгоритм нахождения оптимального вектора $P{\kern 1pt} *$ и оценки $\omega _{Y}^{*}$, обеспечивающих выполнение критерия (1.10) при условиях (1.11) и (1.12),

– исследовать методическую погрешность данного алгоритма,

– определить условия и границы применимости развиваемого метода,

– дать комментарии к достигаемому вычислительному эффекту.

2. ЧИСЛЕННО-АНАЛИТИЧЕСКОЕ ОПИСАНИЕ ПОВЕДЕНИЯ ДИНАМИЧЕСКОЙ СИСТЕМЫ

Для описания поведения ДС в области $G$ (что эквивалентно построению приближенного решения (1.3)) воспользуемся принципом непрерывной зависимости решения от начального условия и идеей использования совокупности опорных интегральных кривых (ОИК) для численно-аналитических построений [12]. С этой целью на отрезке $\left[ {c,d} \right]$ возможных начальных условий зададим сетку $\left\{ {{{x}_{{0\left( m \right)}}}} \right\}_{{m = 1}}^{{{{M}_{x}}}}$ и для каждого фиксированного ${{x}_{{0\left( m \right)}}}$ построим (любым из известных высокоточных численных методов [13–15] на первичной временной сетке большого объема) ОИК

Мы используем здесь непрерывное (по $t$) представление, хотя на практике эта кривая задается таблично.

Далее на отрезке $\left[ {{{t}_{0}},{{t}_{0}} + T} \right]$ задается вторичная (интерполяционная) сетка $\left\{ {{{t}_{{\left( k \right)}}}} \right\}_{{k = 1}}^{{{{M}_{t}}}}$ малого объема, на которой формируется массив чисел

Проводя интерполяцию (по $t$) данного массива (для фиксированного $m$), получаем приближенное представление m-й ОИК в виде где есть вектор коэффициентов, которые находятся путем решения следующей системы линейных алгебраических уравнений (характеризуют условия интерполяции):

Далее, учитывая непрерывную зависимость решения уравнения (1.1) от ${{x}_{0}}$, проводим для фиксированного момента времени ${{t}_{{\left( k \right)}}}$, где $k \in \overline {1,{{M}_{t}}} $, интерполяцию (по ${{x}_{0}}$) коэффициентов $\left\{ {{{a}_{k}}({{x}_{{0\left( m \right)}}})} \right\}_{{m = 1}}^{{{{M}_{x}}}}$:

где координаты вектора ${{B}_{k}} = {{[{{b}_{{km}}},m = \overline {1,{{M}_{x}}} ]}^{{\text{т }}}}$ находятся согласно условиям интерполяции Решая систему линейных алгебраических уравнений (2.4), с учетом (1.3), (2.1)–(2.4), получаем приближенное решение уравнения (1.1): где есть матрица коэффициентов двумерной интерполяции.

Следует помнить, что на этапе применения высокоточного численного метода, как правило, используется сетка большого объема, который многократно превышает объем $\left( {{{M}_{t}}} \right)$ сетки (в случае выполнения условия гладкости (1.2) решения $x\left( {t,\;{{x}_{0}}} \right)$), используемой на этапе временной интерполяции (2.1), (2.2). Кроме того, для многих ДС зависимость решения $x\left( {t,\;{{x}_{0}}} \right)$ от ${{x}_{0}}$ является слабо выраженной, что существенно снижает объем $\left( {{{M}_{x}}} \right)$ сетки интерполяции по начальному условию ${{x}_{0}}$ и упрощает процедуру (2.3), (2.4).

Таким образом, априорно на первом этапе (до получения измерительных данных) строится приближенное решение (2.5) уравнения (1.1) для области G, которое на втором этапе (с учетом наблюдения (1.5)) используется для вычисления оптимальных (с минимальной дисперсией, несмещенных и инвариантных к сингулярной помехе) оценок $\omega _{Y}^{*}$ числовых характеристик и оценки $x_{0}^{*}$ начального условия ${{x}_{0}}$, а также построения оптимальной ОИК

Узлы интерполяции ${{x}_{{0\left( m \right)}}}$ и ${{t}_{{\left( k \right)}}}$ можно выбирать таким образом, чтобы они совпадали с корнями многочленов Чебышева [12]:

Если при построении решения использовать интерполяцию на основе многочленов Лагранжа, то процедура (2.1)–(2.5) приводит к представлению (в данном случае ${{b}_{{km}}} = {{x}_{{\left( {km} \right)}}}$):

где

Пусть выполняется условие (1.4), и, кроме того, полагаем, что для любого${{x}_{0}}$, не совпадающего ни с одним из узлов интерполяции $\left\{ {{{x}_{{0\left( m \right)}}}} \right\}_{{m = 1}}^{{{{M}_{x}}}}$, точное и приближенное решения имеют при ${{t}_{0}}$ невязку ${{\Delta }_{0}}$.Тогда для невязки $\varepsilon $ между этими решениями справедлива оценка (по аналогии с [16])

где

Соотношения (2.6)–(2.8) позволяют оптимизировать выбор параметров метода построения приближенного решения $\tilde {x}(t,{{x}_{0}})$, описывающего поведение ДС, на базе ОИК.

В данном разделе процедура построения численно-аналитического решения уравнения (1.1) реализована на интерполяционных алгоритмах, предполагающих выполнение соответствующих условий интерполяции. Однако данную процедуру можно реализовать и на аппроксимационных алгоритмах (например, наименьших квадратов), что потребует решения задач минимизации соответствующих невязок для времени $t$ и начального условия ${{x}_{0}}$. При учете погрешностей аппроксимации достаточно воспользоваться оценками, которые приводятся в работах [13]–[15].

3. ОБОБЩЕННОЕ ИНВАРИАНТНО-НЕСМЕЩЕННОЕ ОЦЕНИВАНИЕ

Используя найденное решение (2.5), распространим известную процедуру ОИНО на задачу оценивания НЛФ, начального условия ${{x}_{0}}$ и интегральной кривой $x(t,{{x}_{0}})$, соответствующей данному условию применительно к ДС (1.1). Так как вектор отсчетов этой интегральной кривой может быть представлен в виде

то в силу линейности функционала $\omega \in \Omega $ имеем где есть вектор значений функционала $\omega $ на базисных функциях интегральной кривой, а есть матрица их отсчетов ${{\psi }_{{kn}}} = {{\psi }_{k}}\left( {{{t}_{n}}} \right)$.

Непосредственно из (3.1) следует, что для выполнения условия несмещенности (1.11), необходимо

Поскольку применительно к сингулярной помехе

то с учетом этого для соблюдения условия инвариантности (1.12) достаточно выполнения равенства где ${{\left[ 0 \right]}_{{{{M}_{h}} \times 1}}}$ – нулевой вектор столбец размерности ${{M}_{h}} \times 1$, есть матрица отсчетов ${{\theta }_{{hn}}} = {{\theta }_{h}}({{t}_{n}})$ базисных функций сингулярной помехи.

Таким образом решается экстремальная задача нахождения оптимальной оценки $P{\kern 1pt} *$ вектора $P$:

для которой использован метод условной оптимизации Лагранжа. В итоге получили следующую процедуру ОИНО применительно к ДС (по аналогии с [2, см. формулу (2.5)]: где здесь $E = {{E}_{{\left( {N + 1} \right) \times \left( {N + 1} \right)}}}$ – единичная матрица размерности $\left( {N + 1} \right) \times \left( {N + 1} \right)$.

Из (3.4) непосредственно вытекает, что матрицы ${{K}_{\Xi }}$, $\Theta K_{\Xi }^{{ - 1}}{{\Theta }^{{\text{т }}}}$ и $\Psi {{Z}_{1}}{{\Psi }^{{\text{т }}}}$ должны быть невырожденными, что соответствует получению единственного решения ${{P}^{*}}$.

Для дисперсии оптимальной оценки $\omega _{Y}^{*}$ числовой характеристики ${{\omega }_{0}}$ имеем следующее выражение (см. [2, см. формулу (2.20)])

где

Выражения (3.4)–(3.6) имеют явный векторно-матричный вид, что существенно упрощает их применение в любой вычислительной среде. Оценка НЛФ (3.5) формируется сразу без промежуточного вычисления спектральных коэффициентов интегральной кривой и сингулярной помехи, т.е. без расширения пространства состояний. Кроме того, эта оценка является несмещенной, поскольку достигается инвариантность результата оценивания к сингулярной помехе. В силу этого алгоритм нахождения значения НЛФ на базе выражений (3.4) и (3.5) можно отнести к классу линейных автокомпенсационных алгоритмов.

Усредненная методическая погрешность ${{\Delta }_{\omega }}$ оценивания характеристики ${{\omega }_{0}}$, обусловленная ошибкой

аппроксимации точного решения его моделью (2.5), равна (используя символ математического ожидания${\text{M[}}\, \cdot \,{\text{]}}$): где

Поскольку

(условие несмещенности), $\omega _{H}^{*} = {{P}^{{\text{т }}}}H = 0$ (условие инвариантности) и ${\text{M}}[\Xi ] = 0$ (шум является центрированным), непосредственно из (3.7) вытекает равенство ([2, см. формулу (3.23)]):

Таким образом, методическая погрешность оценивания в среднем определяется точным значением функционала на ошибке $\Delta (t,{{x}_{0}})$ и его линейной оценкой на векторе отсчетов $\Delta \left( {{{x}_{0}}} \right)$ этой ошибки.

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

числа измерительных узлов $\left\{ {{{t}_{n}}} \right\}_{{n = 0}}^{N}$ и их расположением, числа степеней свободы в моделях приближенного решения дифференциального уравнения ДС и сингулярной помехи.

Для нахождения оптимальной интегральной кривой, характеризующей траекторию ДС, достаточно найти оптимальную (сглаженную) оценку $x_{0}^{*}$ начального условия ${{x}_{0}}$ для некорректных наблюдений, которую надо подставить в решение (2.5), т.е. $x{\kern 1pt} {\text{*}} = \tilde {x}(t,x_{0}^{*})$. Если ${{x}_{0}}$ рассматривать как значение некоторого функционала $\omega $ на решении $x(t,{{x}_{0}})$, т.е. ${{x}_{0}} = \omega \left\{ {x(t,{{x}_{0}})} \right\}$, то можно воспользоваться оценкой (3.4), (3.5), положив здесь $\omega \left\{ {\Psi (t)} \right\} = {{\left[ {1,0,0, \ldots ,0} \right]}^{{\text{т }}}}.$ В итоге получим

и несмотря на то что рассматриваемая нами ДС является нелинейной, оценивание как ее начального условия, так и самой интегральной кривой осуществляются в рамках линейной процедуры без привлечения традиционно используемых процедур линеаризации нелинейных функций [6]–[9].

4. НЕКОТОРЫЕ ОБОБЩЕНИЯ, ВЫЧИСЛИТЕЛЬНЫЕ АСПЕКТЫ

Предложенный в разд. 2 подход к построению приближенного решения уравнения (1.1) для скалярной ДС, несложно обобщить и на многомерную ДС, описываемую системой уравнений вида (в них правая часть удовлетворяет условию Липшица по $x$)

где есть вектор фазовых координат ДС, есть вектор произвольных постоянных параметров. Полагаются выполненными ограничения:

По аналогии с (2.1) для многомерной ДС в области ${{G}_{{x\eta }}} = \left\{ {{{G}_{x}},{{G}_{\eta }}} \right\}$ вводится семейство узлов, для каждого из которых строится своя опорная интегральная кривая. Далее такая кривая представляется соответствующим интерполяционным полиномом (на сетке $\left\{ {{{t}_{{\left( k \right)}}}} \right\}_{{k = 1}}^{{{{M}_{t}}}}$), а для всех опорных интегральных кривых, с учетом непрерывной зависимости решения от параметров ${{x}_{0}}$ и $\eta $, строится (с использованием процедуры многомерной интерполяции [13]–[15]) искомое численно-аналитическое решение:

где $\alpha ({{x}_{0}},\eta )$ – векторная функция, удовлетворяющая соответствующим характеристическим условиям интерполяции по параметрам ${{x}_{0}}$ и $\eta $. Используя [13]–[15], по аналогии с разд. 2, несложно оценить невязку точного и приближенного решений для области $G = \left\{ {{{G}_{t}},{{G}_{{x\eta }}}} \right\}$.

Задачу ОИНО, рассмотренную в разд. 3, также можно обобщить на многомерный случай, который связан с оценкой векторной числовой характеристики $\omega = {{[{{\omega }_{i}},i = \overline {1,I} ]}^{{\text{т }}}}$ наблюдаемой ДС. В этом случае линейная оценка этой характеристики находится в виде $\omega _{Y}^{*} = {{P}^{*}}Y$, где $P{\kern 1pt} {\text{*}} = [p_{{in}}^{*},i = \overline {1,I} ,n = \overline {0,N} ]$ – оценка матрицы $P$ весовых коэффициентов. Данная матрица находится из условия минимизации следа корреляционной матрицы ${{K}_{Y}} = P{{K}_{\Xi }}{{P}^{{\text{т }}}}$ этой оценки, кроме того, также должны быть выполнены условия несмещенности и инвариантности, рассмотренные в разд. 2 по отношению ко всем координатам вектора $\omega _{Y}^{*}$.

Несложный анализ показывает, что данная задача развивается на $I$ подзадач, связанных с ОИНО по каждой координате вектора, т.е. достигается максимально возможная декомпозиция вычислительной процедуры.

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

5. ИЛЛЮСТРАТИВНЫЙ ПРИМЕР ПОСТРОЕНИЯ ПРИБЛИЖЕННОГО РЕШЕНИЯ

Для оценки эффективности предлагаемого метода ограничимся ДС, описываемой уравнением Риккати специально подобранного вида. Известно, что уравнения Риккати широко используются на практике при тестировании различных методов, а также решении задач аналитического конструирования регуляторов, оценки параметров и состояния различных систем, математической физики и др. [5]–[10], [17].

Пример выбран достаточно простым, чтобы иметь точное решение для последующего сравнительного анализа. Пусть в уравнении (1.1) правая часть задается выражением

Точное решение уравнения Риккати с такой правой частью имеет вид

где ${{\mu }_{0}} = {{x}_{0}}\cos {{t}_{0}} - 1$.

Вычислительный эксперимент реализован при следующих параметрах:

Все используемые вычислительные сетки полагались равномерными (сетка по ${{x}_{0}}$ для построения пяти ОИК с шагом ${{\Delta }_{x}} = 1{\text{/}}{{M}_{x}} = 0.25$; сетка по $t$ для построения самих ОИК на базе численного метода (гарантирующего погрешность, не превышающую ${{10}^{{ - 6}}}$) с шагом ${{\Delta }_{t}} = {{10}^{{ - 2}}}$; интерполяционная сетка по $t$ для построения приближенного решения $\tilde {x}(t,{{x}_{0}})$ с шагом ${{\tilde {\Delta }}_{t}} = 1{\text{/}}{{M}_{t}} = 0.16(6)$).

На фиг. 1 представлены зависимости, характеризующие ДС, а именно, семейства ОИК и их производных.

Для рассматриваемого примера была построена матрица ${\mathbf{B}}$, фигурирующая в (2.5):

Невязка точного $x\left( {t,\;{{x}_{0}}} \right)$ и приближенного $\tilde {x}\left( {t,\;{{x}_{0}}} \right)$ (построенного с учетом (2.5) и данной матрицы ${\mathbf{B}}$) решений в рабочей области $G = \left\{ {\left[ {0,1} \right],\left[ {0,1} \right]} \right\}$ составляет $\delta = 3.8 \times {{10}^{{ - 3}}}$ (или $\delta = 0.3885\% $). Соответственно в усеченной области интерполяции имеем $\delta = 5.439 \times {{10}^{{ - 4}}}$ (или $\delta = 0.2112\% $). Для невязки производных от приближенного решения в области $G$ получаем значение ${{\delta }^{{(1)}}} = 1.646 \times {{10}^{{ - 1}}}$ (относительную погрешность вычислить нельзя, поскольку точное значение производной в соответствующей точке $\left( {{{t}_{0}} = 0,\;{{x}_{0}} = 1} \right)$ равно нулю). Переходя к усеченной области $\bar {G}$, имеем ${{\delta }^{{(1)}}} = 1.117 \times {{10}^{{ - 2}}}$ (или ${{\delta }^{{(1)}}} = 1.705\% $). Эта невязка соответствует точке $\left( {t = 0.21,{{x}_{0}} = 2{\text{/}}3} \right)$. Действительно, усечение области $G$ существенно повышает точность расчетов.

Зависимости невязки $\varepsilon $ от ${{x}_{0}}$ (для ряда фиксированных значений $t$) в области $G$ представлена на фиг. 2 (для детерминированного случая).

Наблюдаем некоторый традиционный рост погрешности на краях сетки по ${{x}_{0}}$. Это необходимо учитывать при практической реализации предлагаемого метода.

6. ИЛЛЮСТРАТИВНЫЙ ПРИМЕР ОЦЕНИВАНИЯ ИНТЕГРАЛЬНОЙ КРИВОЙ ДИНАМИЧЕСКОЙ СИСТЕМЫ И ПРОИЗВОДНОЙ ЕЕ ДВИЖЕНИЯ

Ниже используются исходные данные из разд. 5 и дополнительно следующие:

Для построения оптимальной интегральной кривой $x{\kern 1pt} * = \tilde {x}(t,x_{0}^{*})$ было конкретизировано уравнение наблюдений (1.7), в котором флуктуационный шум формировался датчиком случайных чисел (при ${{\sigma }^{2}} = {{10}^{{ - 2}}}$), а сингулярная помеха, представляемая двумя базисными функциями: задавалась в виде кривой, изображенной на фиг. 3.

В результате реализации алгоритма (3.4), (3.5) и (3.9) было показано, что невязка межу точным $x(t,{{x}_{0}})$ и оптимальным $x{\kern 1pt} {\text{*}}(t,{{x}_{0}})$ решениями в области $G$ составляет $\varepsilon = 7.073 \times {{10}^{{ - 4}}}$ (или $\varepsilon = 0.212\% $). Для усеченной области $\bar {G}$ имеем меньшую невязку $\varepsilon = 5.439 \times {{10}^{{ - 4}}}$ (или $\varepsilon = 0.211\% $).

На фиг. 4 представлены соответственно зависимости невязок $\varepsilon $ и ${{\varepsilon }^{{(1)}}}$ от ${{x}_{0}}$ (при использовании экспериментальных данных), которые характеризуют влияние флуктуационных помех на точность идентификации. Кроме того, было показано, что увеличение числа узлов интерполяции по ${{x}_{0}}$ нецелесообразно. Оптимальным является число ${{M}_{x}} \in \left[ {5,6} \right]$. Соответственно для $t$ такое число ${{M}_{t}} \in \left[ {6,\;7} \right]$.

Были также сформированы наблюдения (с учетом сингулярных и флуктуационных помех, а также принятых исходных данных), относящиеся к интегральной кривой с начальным условием ${{x}_{0}} = 3.8$. Ставилась задача оценить числовую характеристику, а именно: производную от этой кривой (по результатам наблюдений) в точке $t = 0.5$. На базе алгоритма

была получена оценка производной $ - 0.5594$, что соответствует невязке ${{\varepsilon }^{{\left( 1 \right)}}} = 0.0360$ (или ${{\varepsilon }^{{\left( 1 \right)}}} = 0.6416\% $). При расчете невязки помимо данной оценки учитывалось точное значение производной ($ - 0.5630$). Было проверено также, что невязка, возникающая при оценивании производной для всей области $\bar {G} \subset G$, равна ${{\varepsilon }^{{\left( 1 \right)}}} = 0.0117$ (или ${{\varepsilon }^{{\left( 1 \right)}}} = 1.7047\% $).

Данный пример подтверждает возможность применения развитого численно-аналитического метода при исследовании ДС в условиях некорректных наблюдений даже при относительно небольших значениях параметров ${{M}_{x}}$ и ${{M}_{t}}$.

ЗАКЛЮЧЕНИЕ

Развитый численно-аналитический метод исследования ДС ориентирован, в первую очередь, на случаи, связанные с решением широкого круга прикладных задач [5]–[10], требующих вычислений в реальном времени. Данный метод предполагает вынесение основного объема вычислений на первый (предварительный) этап, не связанный непосредственно с наблюдением ДС. К концу этого этапа на базе семейства ОИК мы должны сформировать приближенное аналитическое решение (зависящее от всех характерных параметров, значения которых могут принадлежать заданному множеству допустимости), а также сформировать вектор-столбец или матрицу весовых коэффициентов ОИНО, учитывающих спектральный состав интегральной кривой и сингулярной помехи, а также статистические характеристики флуктуационного шума. На втором этапе для построения оптимальной интегральной кривой ДС (по результатам некорректных наблюдений), достаточно подставить в сформированное решение оптимальную оценку начального условия, найденную для частного случая общей процедуры ОИНО. Любые другие числовые характеристики поведения ДС находятся в виде скалярного произведения вектора наблюдений и вектора весовых коэффициентов общего вида, не прибегая к непосредственному вычислению соответствующего НЛФ для найденной оптимальной траектории ДС. Предполагаемый метод не требует выполнения процедур линеаризации (не принимая во внимание этап, связанный с построением ОИК, даже если ДС является нелинейной), не требует расширения пространства состояний и не связан с необходимостью самостоятельного расчета спектральных коэффициентов в линейных комбинациях интегральной кривой и сингулярной помехи. Метод допускает возможность комбинирования с традиционными подходами к решению прикладных задач, связанных с оптимальной обработкой измерений [6], [7].

Дальнейшее совершенствование метода может идти по следующим направлениям:

– рассмотреть многомерную ДС при наличии в правой части соответствующего дифференциального уравнения вектора коэффициентов с априорно неизвестными значениями, но заданными множествами допустимости;

– ввести в уравнение наблюдения векторный измеряемый параметр, который нелинейно связан с вектором состояния ДС;

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

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

– рассмотреть возможность применения развитого метода для решения задач пассивной и активной идентификации ДС на основе заранее спланированного эксперимента.

Для решения задач идентификации многомерных ДС можно воспользоваться идеей построения численно-аналитического решения дифференциального уравнения на базе экспериментального семейства опорных интегральных кривых [18].