Журнал вычислительной математики и математической физики, 2020, T. 60, № 7, стр. 1151-1169

ВВЕДЕНИЕ

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

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

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

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

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

Пусть ДС описывается векторным уравнением в нормальной форме

где есть решение уравнения с начальным условием ${{x}_{0}} = x\left( {{{t}_{0}},{{x}_{0}}} \right) = {{[{{x}_{{i0}}},i = \overline {1,I} ]}^{{\text{T}}}} \in {{G}_{{x0}}}$, ${{G}_{{x0}}} = \{ {{x}_{0}} \in {{\mathbb{R}}^{I}}{\kern 1pt} {\kern 1pt} :\;{{с}_{i}} \leqslant {{x}_{{i0}}} \leqslant {{d}_{i}},\;i = \overline {1,I} \} $ – гиперпараллелепипед возможных начальных условий, $f(t,x)$ – функция, имеющая ту гладкость (по аргументам $t$ и $x$), которая соответствует рассматриваемой ДС и требуется по смыслу проводимых в дальнейшем выкладок.

При гладкой правой части $f(t,x)$ (с учетом существования всевозможных производных определенного порядка по $t$ и $x$) решение $x\left( {t,{{x}_{0}}} \right)$ уравнения (1.1) и его соответствующие частные производные существуют и непрерывны в области ${{G}_{{xt}}} = \left\{ {{{G}_{{x0}}},{{G}_{t}}} \right\}$ (см. [8], с. 124, 125). Будем полагать, что функции ${{x}_{i}}\left( {t,\;{{x}_{0}}} \right)$ принадлежат пространству функций ${{C}^{r}}\left[ {{{G}_{{xt}}}} \right]$ с равномерной сходимостью (здесь число $r$ ограничивает длину соответствующего мультиндекса для всевозможных частных производных), при этом норме этого пространства соответствует обозначение ${{\left\| {\, \cdot \,} \right\|}_{{C\left[ {{{G}_{{xt}}}} \right]}}}$ (см. [9], с. 51, 52).

Уравнению (1.1) можно поставить в соответствие приближенное аналитическое решение $\tilde {x} = \tilde {x}\left( {t,{{x}_{0}}} \right)$, где ${{\tilde {x}}_{0}} = \tilde {x}\left( {{{t}_{0}},{{x}_{0}}} \right)$, при этом в общем случае $\mathop {\max }\limits_i \left| {{{{\tilde {x}}}_{{i0}}} - {{x}_{{i0}}}} \right| \ne 0,$ $i = \overline {1,I} $. Считаем, что

является ${{{\varepsilon }}_{0}}$ – приближенным по невязке решением [10], т.е. и при подстановке функций ${{\tilde {x}}_{i}}\left( {t,{{x}_{0}}} \right)$ в (1.1) получаем где невязки ${{\vartheta }_{i}}\left( t \right)$ удовлетворяют неравенству

В дальнейшем близость точной ${{x}_{i}} = {{x}_{i}}\left( {t,{{x}_{0}}} \right)$ и приближенной ${{\tilde {x}}_{i}} = {{\tilde {x}}_{i}}\left( {t,{{x}_{0}}} \right)$ фазовых координат ДС для фиксированного ${{x}_{0}}$ будем характеризовать невязкой

соответственно для оценки близости их производных первого порядка (по аргументу $t$) вводится невязка Для оценки близости векторных решений $x\left( {t,{{x}_{0}}} \right)$, $\tilde {x}\left( {t,{{x}_{0}}} \right)$ и их производных ${{x}^{{\left( 1 \right)}}}\left( {t,{{x}_{0}}} \right)$, ${{\tilde {x}}^{{\left( 1 \right)}}}\left( {t,{{x}_{0}}} \right)$ используются невязки ${{{\varepsilon }}_{x}}\left( {{{x}_{0}}} \right){\text{ = }}\mathop {{\text{max}}}\limits_i {{{\varepsilon }}_{{xi}}}\left( {{{x}_{0}}} \right)$ и ${\varepsilon }_{x}^{{\left( 1 \right)}}\left( {{{x}_{0}}} \right) = \mathop {\max }\limits_i {\varepsilon }_{{xi}}^{{\left( 1 \right)}}\left( {{{x}_{0}}} \right)$ соответственно.

Зададимся также векторным уравнением измеряемого выхода ДС [11], [12]:

при любом $t \in {{G}_{t}}$. Здесь $y\left( t \right) = {{[{{y}_{j}}\left( {t,{{x}_{0}}} \right),j = \overline {1,J} ]}^{{\text{T}}}}$– вектор измеряемых функций с координатами ${{y}_{j}}\left( {t,{{x}_{0}}} \right)$, которые также принадлежат ${{C}^{r}}\left[ {{{G}_{{xt}}}} \right]$, ${\varphi }(t,x)$ – гладкая функция своих аргументов (степень гладкости также определяется типом измеряемого выхода рассматриваемой ДС).

На траекториях ${{y}_{j}}\left( {t,{{x}_{0}}} \right) \in {{C}^{r}}\left[ {{{G}_{{xt}}}} \right]$ вводится семейство $\Omega $ ЛНФ ${\omega }$, которые необходимы при решении тех или иных целевых задач, связанных с наблюдением за конкретной ДС, при этом будем использовать представление

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

где В этом уравнении: $x_{0}^{'} \in {{G}_{{x0}}}$ – неизвестное начальное условие для ДС, ${{s}_{j}}\left( t \right)$ – сингулярная помеха, ${{{\xi }}_{j}}\left( t \right)$ – флуктуационная погрешность. Введение индекса $j$ во временную сетку $\left\{ {{{t}_{{jn}}}} \right\}_{{n = 0}}^{{{{N}_{j}}}}$ позволяет рассмотреть случай асинхронных измерений различных функций выхода исследуемой ДС.

Помеху опишем в виде

где ${{C}_{j}} = {{[{{c}_{{jp}}},p = \overline {1,{{M}_{{sj}}}} ]}^{{\text{T}}}}$ – вектор неизвестных спектральных коэффициентов, ${{\Theta }_{j}}\left( t \right) = {{[{{{\theta }}_{{jp}}}\left( t \right),p = \overline {1,{{M}_{{sj}}}} ]}^{{\text{T}}}}$ ${{\Theta }_{j}}\left( t \right) = {{[{{{\theta }}_{{jp}}}\left( t \right),p = \overline {1,{{M}_{{sj}}}} ]}^{{\text{T}}}}$ – вектор заданных линейно независимых функций, ${{M}_{{sj}}}$ – заданное число степеней свободы (на практике зачастую прибегают к упрощению ${{M}_{{sj}}} = {{M}_{s}}$ и ${{\Theta }_{j}}\left( t \right) = \Theta \left( t \right)$).

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

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

Введем векторные обозначения:

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

Для фиксированного $j \in \overline {1,J} $ из множества $\Omega $ выберем произвольный функционал ${\omega }$: ${{y}_{j}}\left( {t,{{x}_{0}}} \right) \to {{\mathbb{R}}^{1}}$, в другой записи ${\omega }\left\{ {{{y}_{j}}\left( {t,{{x}_{0}}} \right)} \right\} = {{{\omega }}_{j}} \in {{\mathbb{R}}^{1}}$. На основе наблюдения ${{H}_{j}}$ мы должны дать оценку значения ${{{\omega }}_{j}}$, т.е. вычислить числовую характеристику измеряемой функции ${{y}_{j}}\left( {t,{{x}_{0}}} \right)$, соответствующую функционалу ${\omega } \in \Omega $. Данную оценку будем искать в классе линейных оценок (здесь и далее индекс $ * $ соответствует некоторой оценке) [13]:

где ${{P}_{j}} = {{[{{p}_{{jn}}},\;n = \overline {0,\;N} ]}^{{\text{T}}}}$ – вектор искомых весовых коэффициентов.

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

Выбор оптимального значения $P_{j}^{*}$ для вектора ${{P}_{j}}$ соответствует критерию

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

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

– сформировать алгоритм построения численно-аналитического решения $\tilde {x}\left( {t,{{x}_{0}}} \right)$ для многомерной ДС (1.1) применительно к области ${{G}_{{xt}}}$,

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

– сформировать алгоритм построения численно-аналитического выражения для измеряемого многомерного выхода (1.5) применительно к области ${{G}_{{xt}}}$,

– сформировать алгоритм нахождения оптимального вектора $P_{j}^{*}$ и оценки ${\omega }_{j}^{*}$, обеспечивающих выполнение критерия (1.12) при условиях (1.13) и (1.14), рассмотреть возможность оценивания различных НЛФ на наблюдениях (1.9),

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

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

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

2. ЧИСЛЕННО-АНАЛИТИЧЕСКОЕ РЕШЕНИЕ ДИФФЕРЕНЦИАЛЬНОГО УРАВНЕНИЯ ДИНАМИЧЕСКОЙ СИСТЕМЫ (ЭТАП 1)

Для формирования численно-аналитического решения уравнения (1.1) в области ${{G}_{{xt}}}$ можно использовать широко применяющуюся в общей теории приближенных методов [14] математическую конструкцию

где $\sum\nolimits_k^{} {\kern 1pt} $ – одномерная сумма по индексу $k$ с заданным числом неизвестных коэффициентов $\beta _{{ik}}^{x}$ и известных линейно независимых функций $\varsigma _{{ik}}^{x}\left( {t,{{x}_{0}}} \right)$. Очевидно, что с учетом гладкости решения ${{x}_{i}}\left( {t,\;{{x}_{0}}} \right)$ число слагаемых в сумме можно подобрать так, чтобы для любого сколь угодно малого ${{\delta }_{{xi}}} \geqslant 0$ обеспечивалось выполнение условия ${{{\varepsilon }}_{{xi}}}\left( {{{x}_{0}}} \right) < {{\delta }_{{xi}}}$ не зависимо от конкретного значения ${{x}_{0}} \in {{G}_{{x0}}}$.

В частном случае можно перейти от (2.1) к приближению на основе многомерного многочлена (см. [15], с. 152)

${\text{z}} = \left( {t,{{x}_{{10}}},...,{{x}_{{I0}}}} \right)$, $\sum\nolimits_{\text{m}}^{} {\kern 1pt} $ – многомерная сумма по мультиндексу ${\text{m}}$ с заданным числом неизвестных коэффициентов $\alpha _{{i{\text{m}}}}^{x}$. Представление (2.2) используется в теореме Вейерштрасса (для многомерного случая, см. [15], теорема 5, с. 152), подтверждающей возможность наилучшего приближения к ${{x}_{i}}\left( {t,{{x}_{0}}} \right)$ в равномерной метрике для области ${{G}_{{xt}}}$.

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

где $\Psi _{i}^{x}\left( t \right)$ и $\Lambda _{i}^{x}\left( {{{x}_{0}}} \right)$ – векторы заданных гладких линейно независимых функций соответствующей размерности, $B_{i}^{x}$ – матрица неизвестных коэффициентов. Если в качестве координат векторов $\Psi _{i}^{x}\left( t \right)$ и $\Lambda _{i}^{x}\left( {{{x}_{0}}} \right)$ использовать многомерные многочлены с соответствующими мультиндексами (по аналогии с [15], с. 152), то представление (2.3) несложно получить из (2.1). То есть, в данном частном случае, теорема Вейерштрасса подтверждает наличие наилучшего приближения к ${{x}_{i}}\left( {t,{{x}_{0}}} \right)$ в рамках конструкции (2.3). В общем случае численно-аналитическое решение (2.3) можно обосновать так. Пусть имеется частное решение ${{x}_{i}}(t,x_{0}^{'})$, соответствующее начальному условию $x_{0}^{'} \in {{G}_{{x0}}}$. Согласно [14], [15] для фиксированного $x_{0}^{'} \in {{G}_{{x0}}}$ и любого сколь угодно малого ${{\delta }_{{xi}}} \geqslant 0$ можно подобрать приближение такое, что будет выполнено условие ${{{\varepsilon }}_{{xi}}}(x_{0}^{'}) < {{\delta }_{{xi}}}$, где $b_{{ik}}^{x}(x_{0}^{'})$ – спектральные коэффициенты, значения которых зависят от выбора $x_{0}^{'} \in {{G}_{{x0}}}$. Так как точное решение ${{x}_{i}}\left( {t,{{x}_{0}}} \right)$ является гладким, а число слагаемых под знаком суммы (для приближенного решения) не зависит от выбора начального условия из области ${{G}_{{x0}}}$ (согласно предположению), то можно подобрать такое ${{x''}_{0}} \in {{G}_{{x0}}}$, что для приближения будет обеспечиваться выполнение условия ${{{\varepsilon }}_{{xi}}}(x_{0}^{{''}}) < {{\delta }_{{xi}}}$. При этом, если невязка ${\text{|}}x_{0}^{'} - x_{0}^{{''}}{\text{|}}$ мала, то и приращения ${\text{|}}b_{{ik}}^{'} - b_{{ik}}^{{''}}{\text{|}}$ будут также малыми. Следовательно, можно подобрать такие гладкие функции $b_{{ik}}^{x}\left( {{{x}_{0}}} \right)$, что для всех ${{x}_{0}} \in {{G}_{{x0}}}$ и $t \in {{G}_{t}}$ данному условию будет удовлетворять и численно-аналитическое решение При этом, если коэффициенты $b_{{ik}}^{x}\left( {{{x}_{0}}} \right)$ представить соответствующей линейной комбинацией гладких функций ${\lambda }_{{ir}}^{x}\left( {{{x}_{0}}} \right)$ с необходимым числом степеней свободы , то можно добиться выполнения условия ${{{\varepsilon }}_{{xi}}}\left( {{{x}_{0}}} \right) < {{\delta }_{{xi}}}$ и для приближения (2.3).

В формуле (2.4) использованы обозначения: $B_{i}^{x}\left( {x{}_{0}} \right)$ – вектор-столбец спектральных коэффициентов $b_{{ik}}^{x}\left( {{{x}_{0}}} \right)$, зависящих от ${{x}_{0}} \in {{G}_{{x0}}}$, $\Psi _{i}^{k}(t)$ – вектор-столбец линейно независимых функций ${\psi }_{{ik}}^{x}\left( t \right)$.

На практике конструкция (2.3) широко распространена для описания большого класса ДС (см., например, [5], [6]).

Для расчета коэффициентов, фигурирующих в математических конструкциях (2.1)–(2.4), воспользуемся идеей, основанной на использовании экспериментального семейства ОИК [1], [2], [5]. Как показывает опыт, применение этой идеи оказалось продуктивным в различных областях, в том числе при анализе и синтезе управляемых и неуправляемых ДС [16], при этом большие затраты на построение данного семейства компенсируются возможностью формирования качественных численно-аналитических решений для сравнительно больших областей ${{G}_{{xt}}}$. Пример неудачного построения решения на базе известного метода Тейлора продемонстрирован в работе [2], в то время как метод ОИК дал хороший результат. По существу, данный метод позволяет за счет учета гладкой зависимости решения от времени, начальных условий, а также других параметров (которые могут присутствовать в правой части дифференциального уравнения) рационально “сжать” всю информацию, содержащуюся в семействе ОИК, не прибегая к сложным аналитическим операциям (как, например, в методе Тейлора и его различных модификациях) над правой частью уравнения (1.1).

Обобщим метод ОИК [2] на векторное уравнение (1.1). Для этого в области ${{G}_{{x0}}}$ рассмотрим сетку, состоящую из узлов ${{x}_{{0(r)}}} \in {{G}_{{x0}}}$, $r \in \overline {1,{{M}_{x}}} $. Этим узлам поставим в соответствие семейство ОИК

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

Теперь в области ${{G}_{t}}$ задается временная сетка $\left\{ {{{t}_{{\left( {ik} \right)}}}} \right\}$, $k = \overline {1,\;{{M}_{{ti}}}} $, объема которой достаточно для представления фазовых координат ${{x}_{i}}\left( {t,\;{{x}_{{0\left( r \right)}}}} \right)$ в области ${{G}_{{xt}}}$ с требуемой точностью. Далее с использованием семейства ОИК (2.5) формируется массив чисел

Для нахождения неизвестных коэффициентов в формулах (2.1)–(2.4) можно применять известные операторы ${\rm E}_{1}^{x}$, ${\rm E}_{2}^{x}$ и ${\rm E}_{3}^{x}$ (интерполяции или аппроксимации [15], [17], [18]), при этом с учетом (2.6), используя эти операторы, вычисляются указанные коэффициенты При построении численно-аналитического решения $\tilde {x}\left( {t,{{x}_{0}}} \right)$ уравнения (1.1) на базе (2.7) в общем случае применяется неравномерная сетка $\left\{ {{{t}_{{i\left( k \right)}}}} \right\}_{{k = 1}}^{{{{M}_{{ti}}}}}$ по аргументу $t$. Кроме того, для многих ДС зависимость фазовой координаты ${{x}_{i}}\left( {t,{{x}_{0}}} \right)$ от ${{x}_{0}}$ является слабо выраженной, что существенно снижает объем ${{M}_{x}}$ вычислительной сетки по начальному условию ${{x}_{0}}$ и упрощает рассмотренную процедуру (интерполяционную или аппроксимационную) построения численно-аналитического решения $\tilde {x}\left( {t,{{x}_{0}}} \right)$ для уравнения (1.1).

Таким образом, априорно на первом этапе (до получения измерительных данных) строится приближенное аналитическое решение $\tilde {x}\left( {t,{{x}_{0}}} \right)$ уравнения (1.1), обеспечивающее в области ${{G}_{{xt}}}$ требуемую для практики точность анализа ДС. В число параметров, позволяющих достичь такой точности, можно отнести объемы используемых сеток (${{M}_{{ti}}}$ и ${{M}_{x}}$), кроме того важен выбор систем линейно независимых функций в приближениях (2.1)–(2.4) с учетом вида правой части уравнения (1.1).

Конкретизируем описанную процедуру на случай, когда при построении решения на базе (2.3) используются фундаментальные многочлены интерполяции (как тензорное произведение одномерных фундаментальных многочленов [15]). С этой целью в области ${{G}_{{x0}}}$ зададим многомерную вычислительную сетку

где ${\text{m}} = \left( {{{m}_{1}},...,{{m}_{I}}} \right)$ – мультиндекс (многомерный узел), (M_x – объем сетки по ${{x}_{0}}$) и для всех ее узлов построим семейство ОИК

На отрезке ${{c}_{i}} \leqslant {{x}_{{i0}}} \leqslant {{d}_{i}}$ (см. пояснения к формуле (1.1)) задаются узлы ${{x}_{{i0\left( 1 \right)}}},...,{{x}_{{i0\left( {{{M}_{{xi}}}} \right)}}}$ и принимается

В этом случае решение (2.3) можно представить в виде

где $\sum\nolimits_{\text{m}}^{} {\kern 1pt} $ – многомерная сумма,

Кроме того, иногда целесообразно использовать неравномерные сетки, например такие, в которых узлы интерполяции ${{x}_{{i0\left( {{{m}_{i}}} \right)}}}$ и ${{t}_{{\left( {ik} \right)}}}$ совпадают с корнями многочленов Чебышёва:

Если выполнены условия (1.2)–(1.4), то для невязки ${{{\varepsilon }}_{{xi}}}$ между фазовыми координатами ${{x}_{i}}\left( {t,\;{{x}_{0}}} \right)$ и ${{\tilde {x}}_{i}}\left( {t,\;{{x}_{0}}} \right)$ справедлива оценка (по аналогии с [10])

где $t \in {{G}_{t}}$, ${{x}_{0}} \in {{G}_{{x0}}}$, ${{L}_{0}}$ – константа Липшица для правой части уравнения (1.1).

Соответственно для оценки близости двух решений $x\left( {t,\;{{x}_{0}}} \right)$ и $\tilde {x}\left( {t,\;{{x}_{0}}} \right)$ можно воспользоваться оценкой

где $t \in {{G}_{t}}$, ${{x}_{0}} \in {{G}_{{x0}}}$. При $I = 1$ из (2.14) следует (2.13).

Опираясь на работы [15], [17], [18], можно также оценить вклад наследственной погрешности и погрешности округлений в результирующую ошибку построения приближенного решения уравнения (1.1) для области ${{G}_{{xt}}}$.

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

Соотношения (2.1)–(2.14) составляют суть численно-аналитического описания решения дифференциального уравнения ДС для заданной области допустимости ${{G}_{{xt}}}$.

3. ЧИСЛЕННО-АНАЛИТИЧЕСКОЕ ОПИСАНИЕ ИЗМЕРЯЕМОГО ВЫХОДА ДИНАМИЧЕСКОЙ СИСТЕМЫ (ЭТАП 1)

Точный измеряемый выход ДС с учетом (1.5) можно представить в виде

при этом ${\varphi }{\kern 1pt} :\;{{G}_{{xt}}} \to {{G}_{{yt}}}$ – гладкая функция по $t$ и ${{x}_{0}}$.

Соответственно для приближенного выхода ДС имеем

По аналогии с разд. 2, для численно-аналитического описания измеряемого выхода ДС можно использовать семейство ОИК, при этом также применяем принцип гладкой зависимости $y\left( t \right) = {\gamma }\left( {t,{{x}_{0}}} \right)$ от своих аргументов. Для этой цели семейству ОИК ${{x}_{{\left( r \right)}}} = x(t,{{x}_{{0\left( r \right)}}}),$ $r = \overline {1,{{M}_{x}}} $, поставим в соответствие семейство выходных траекторий

которые представим массивом чисел

Теперь по аналогии с (2.7) можно получить оценки неизвестных коэффициентов в приближенных выражениях для измеряемого выхода ДС, если в формулах (2.1)–(2.4) заменить индекс $x$ на $y$, а индекс $i$ на $j$:

По аналогии с (2.3) и (2.8)–(2.11) выход ДС можно представить в виде

Соотношения (3.1)–(3.7) составляют суть численно-аналитического описания измеряемого выхода ДС, при этом формулы (3.6) и (3.7) являются наиболее удобными в численно-аналитических расчетах.

4. ОЦЕНИВАНИЕ ЧИСЛОВЫХ ХАРАКТЕРИСТИК ИЗМЕРЯЕМЫХ ПАРАМЕТРОВ ДИНАМИЧЕСКОЙ СИСТЕМЫ (ЭТАП 2)

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

Распространим известную процедуру обобщенного инвариантно-несмещенного оценивания (ОИНО) [3] на задачу вычисления числовых характеристик измеряемого выхода ДС. С учетом (3.2) и по аналогии с (2.4) выходную траекторию ДС с заданной точностью можно представить в виде

где ${{\bar {Y}}_{j}} = {{[{{y}_{j}}\left( {{{t}_{{\left( k \right)}}}} \right),k = \overline {1,{{M}_{t}}} ]}^{{\text{T}}}} = B_{i}^{y}(x_{0}^{'})$ – вектор временных отсчетов траектории ${{\tilde {y}}_{j}}(t) = {{{\tilde {\gamma }}}_{j}}(t,x_{0}^{'})$ на вычислительной сетке $\left\{ {{{t}_{{\left( k \right)}}}} \right\}_{{k = 1}}^{{{{M}_{t}}}}$, $x_{0}^{'} \in {{G}_{{x0}}}$– неизвестное начальное условие, фигурирующее в уравнении наблюдения (1.8), $\Psi _{j}^{y}(t) = {{[{\psi }_{{jk}}^{y}(t),k = \overline {1,{{M}_{t}}} ]}^{{\text{T}}}}$ – вектор линейно независимых функций.

Полагаем, что

при этом вектор временных отсчетов координаты ${{y}_{{jn}}} = {{{\gamma }}_{j}}\left( {{{t}_{{jn}}},{{x}_{0}}} \right)$ на измерительной сетке $\left\{ {{{t}_{{jn}}}} \right\}_{{n = 0}}^{{{{N}_{j}}}}$ (где ${{N}_{j}} \gg {{M}_{t}}$) приближенно может быть задан в виде где $\Psi _{i}^{y} = [{\psi }_{{jkn}}^{y},k = \overline {1,\;{{M}_{t}}} ,n = \overline {0,\;{{N}_{j}}} ]$ – матрица отсчетов ${\psi }_{{jkn}}^{y} = {\psi }_{{jk}}^{y}\left( {{{t}_{{jn}}}} \right)$ функций ${\psi }_{{jk}}^{y}(t)$.

В силу линейности функционала ${\omega } \in \Omega $, условие несмещенности (1.13) принимает вид

где ${\omega }\{ \Psi _{j}^{y}(t)\} = {{[{\omega }\{ {\psi }_{{jk}}^{y}\left( t \right)\} ,k = \overline {1,\;{{M}_{t}}} ]}^{{\text{T}}}}$ – вектор значений функционала ${\omega }$ на линейно независимых функциях ${\psi }_{{jk}}^{y}(t)$.

Непосредственно из (4.4) вытекает окончательное условие несмещенности

Поскольку необходимо, чтобы

то условие инвариантности (1.14) можно преобразовать к виду где ${{\left[ 0 \right]}_{{{{M}_{{sj}}} \times 1}}}$ – нулевой вектор-столбец размерности ${{M}_{{sj}}} \times 1$, ${{\Theta }_{j}} = [{{{\theta }}_{{jsn}}},s = \overline {1,{{M}_{{sj}}}} ,n = \overline {0,{{N}_{j}}} ]$ – матрица отсчетов ${{{\theta }}_{{jsn}}} = {{{\theta }}_{{js}}}\left( {{{t}_{{jn}}}} \right)$ линейно независимых функций сингулярной помехи.

Задача нахождения оптимальной оценки $P_{j}^{*}$ вектора ${{P}_{j}}$ решается методом условной оптимизации Лагранжа на базе критерия (1.12) с учетом ограничений (1.13) и (1.14). В итоге получаем следующую процедуру автокомпенсационного оценивания НЛФ применительно к измеряемому выходу ДС (по аналогии с [3]):

где ${{Z}_{{j1}}} = {{Z}_{{j2}}}K_{{j\Xi }}^{{ - 1}},$ ${{Z}_{{j2}}} = {{E}_{j}} - K_{{j\Xi }}^{{ - 1}}\Theta _{j}^{{\text{T}}}{{({{\Theta }_{j}}K_{{j\Xi }}^{{ - 1}}\Theta _{j}^{{\text{T}}})}^{{ - 1}}}{{\Theta }_{j}},$ ${{E}_{j}}$ – единичная матрица.

В выражениях (4.8), (4.9) имеем следующие размерности для матриц:

Из (4.8) непосредственно вытекает, что матрицы ${{K}_{{j\Xi }}}$, ${{\Theta }_{j}}K_{{j\Xi }}^{{ - 1}}\Theta _{j}^{{\text{T}}}$ и $\Psi _{j}^{y}{{Z}_{{j1}}}{{(\Psi _{j}^{y})}^{{\text{T}}}}$должны быть невырожденными и хорошо обусловленными, что соответствует получению единственного и устойчивого решения $P_{j}^{*}$. Данный вопрос относится к сфере оптимального планирования измерительного эксперимента и в данной статье не рассматривается. При рациональном выборе основных параметров процедуры (4.8), (4.9) можно избежать необходимости решения некорректной (с вычислительной точки зрения) задачи оценивания.

Для дисперсии оптимальной оценки ${\omega }_{j}^{*}$ числовой характеристики ${{{\omega }}_{j}}$ имеем следующее выражение (по аналогии с [3]):

где

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

Дадим оценку усредненной методической погрешности ${{{\delta }}_{j}}$ оценивания характеристики ${\omega }$ применительно к представлению (3.6). Поскольку в (3.6) не учитывается “хвост”

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

Так как ${\omega }\{ {{[\Psi _{j}^{y}(t)]}^{{\text{T}}}}B_{j}^{y}\Lambda _{j}^{y}\left( {{{x}_{0}}} \right)\} = {\omega }_{{yj}}^{*}$ (условие несмещенности), ${\omega }_{{Sj}}^{*} = P_{j}^{{\text{T}}}{{S}_{j}} = 0$ (условие инвариантности) и ${\text{M}}\left[ {{{\Xi }_{j}}} \right] = 0$ (шум является центрированным), то непосредственно из (4.11) вытекает равенство

Таким образом, методическая погрешность оценивания в среднем определяется значением НЛФ на “хвосте” бесконечного ряда и линейной оценкой этого значения.

Необходимые и достаточные условия существования и единственности решения (4.8), (4.9) задачи оценивания числовой характеристики измеряемого выхода ДС требуют невырожденности и соблюдения некоторых ограничений на ранги ряда матриц. Выполнение этих условий на практике обеспечивается рациональным выбором векторов $\Psi _{j}^{y}(t)$ и ${{\Theta }_{j}}(t)$ линейно-независимых функций, числа измерительных узлов $\left\{ {{{t}_{{jn}}}} \right\}_{{n = 0}}^{{{{N}_{j}}}}$ и их расположением, числа степеней свободы в моделях приближенного решения дифференциального уравнения и измеряемого выхода ДС, а также сингулярной помехи. Все эти вопросы относятся к планированию измерительного эксперимента. Используя результаты работ [4]–[6], [11], можно получить оптимальные оценки параметров предлагаемого метода с учетом статистических характеристик флуктуационных погрешностей измерений. Мы не рассматриваем эти вопросы, чтобы не перегружать объем статьи хотя и важными для практики, но второстепенными для теории результатами.

Теперь рассмотрим вопрос, связанный с нахождением оптимальной траектории $j$-го измеряемого выхода ДС (где $j \in \overline {1,J} $). Под ней понимается траектория (4.1), вектор спектральных коэффициентов которой ${{\bar {Y}}_{j}}$ находится оптимальным образом (с учетом (1.12)–(1.14)) в соответствии с алгоритмом ОИНО (4.8), (4.9). Согласно данному алгоритму, достаточно найти оптимальные оценки $y_{{j\left( k \right)}}^{*}$ отсчетов траектории ${{y}_{j}}(t) = {{{\gamma }}_{j}}(t,x_{0}^{'})$ в узлах ${{t}_{{\left( k \right)}}},\;k = \overline {1,{{M}_{t}}} $, с учетом некорректных наблюдений (1.10). Если истинные значения ${{y}_{{j\left( k \right)}}}$ рассматривать как НЛФ ${\omega }$ на траектории ${{y}_{j}}(t) = {{{\gamma }}_{j}}(t,x_{0}^{'})$, т.е. ${{y}_{{j\left( k \right)}}} = {\omega \{ }{{{\gamma }}_{j}}(t,x_{0}^{'}){\text{\} }}$, то, принимая во внимание (4.6), (4.7), можно воспользоваться оценками

положив здесь ${\omega }\left\{ {{{\Psi }_{{yj}}}\left( t \right)} \right\} = {{\left[ {0,0,...,1,0,...,0} \right]}^{{\text{T}}}}$ – нулевой столбец с единицей в $k$-й позиции.

Зная оценку $\bar {Y}_{j}^{*} = {{[y_{{j\left( k \right)}}^{*},k = \overline {1,{{M}_{t}}} ]}^{{\text{T}}}}$ вектора ${{\bar {Y}}_{j}} = {{[{{y}_{j}}({{t}_{{\left( k \right)}}}),k = \overline {1,{{M}_{t}}} ]}^{{\text{T}}}}$ с учетом (4.1)–(4.3) находим оптимальные скалярные траектории

для всех измеряемых функций ДС.

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

Формулы (4.1)–(4.14) могут быть использованы для решения многих прикладных задач (указанных во Введении), например, для прогнозирования эволюции измеряемых функций ДС во времени.

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

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

где ${\eta } = {{[{{{\eta }}_{i}},\;i = \overline {1,\;{{I}_{2}}} ]}^{{\text{T}}}}$ – вектор произвольных постоянных параметров.

Здесь полагаются выполненными ограничения:

где $\left[ {{{c}_{i}},{{d}_{i}}} \right] \subset {{\mathbb{R}}^{1}},$ $i = \overline {1,{{I}_{1}}} $, $\left[ {{{r}_{i}},{{d}_{i}}} \right] \subset {{\mathbb{R}}^{1}},$ $i = \overline {1,{{I}_{2}}} $.

Вводя расширенный вектор $x = {{[{{z}^{{\text{T}}}},{{{\eta }}^{{\text{T}}}}]}^{{\text{T}}}} \in {{\mathbb{R}}^{{{{I}_{1}} + {{I}_{2}}}}} = {{\mathbb{R}}^{I}}$ и учитывая, что $d{\eta /}dt \equiv 0$, приходим к уравнению (1.1), в котором

где ${{\bar {\bar {f}}}_{i}}\left( {t,\;z,\;{\eta }} \right) \equiv 0,\;\;i = \overline {1,{{I}_{2}}} $, при этом ${{G}_{{xt}}} = \left\{ {{{G}_{z}},{{G}_{{\eta }}},{{G}_{t}}} \right\}$.

Таким образом, класс ДС вида (5.1), которые широко используются на практике (например, в баллистике [19], [20], и ряде задач идентификации и оптимального управления [16]), также попадает в сферу применения развитого метода.

Задачу ОИНО, рассмотренную в разд. 4, также можно обобщить на многомерный случай, который связан с оценкой векторной числовой характеристики ${\omega } = {{[{{{\omega }}_{d}},\;d = \overline {1,\;D} ]}^{{\text{T}}}}$ на выходных траекториях ДС. В этом случае линейная оценка векторного НЛФ ${\omega }\left\{ {{{y}_{j}}} \right\}$ для $j$-го измеряемого параметра находится в векторно-матричном виде ${\omega }_{j}^{*} = P_{j}^{*}{{H}_{j}}$, где $P_{j}^{*} = [p_{{jdn}}^{*},d = \overline {1,D} ,n = \overline {0,N} ]$ – оценка матрицы ${{P}_{j}}$ весовых коэффициентов. Данная матрица находится из условия минимизации следа корреляционной матрицы ${{K}_{j}} = {{P}_{j}}{{K}_{{j\Xi }}}P_{j}^{{\text{T}}}$ этой оценки, кроме того, также должны быть выполнены условия несмещенности и инвариантности, рассмотренные в разд. 4 по отношению ко всем координатам вектора ${\omega }_{j}^{*}$. Несложный анализ показывает, что данная задача развивается на $D$ подзадач, связанных с ОИНО по каждой координате вектора ${\omega } = {{[{{{\omega }}_{d}},\;d = \overline {1,\;D} ]}^{{\text{T}}}}$, т.е. достигается максимально возможная декомпозиция вычислительной процедуры.

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

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

В качестве НЛФ могут рассматриваться выборочное значение $j$-й измеряемой функции в произвольной точке $t' \in {{G}_{t}}$ (в этом случае ${\omega }\{ {{\Psi }_{y}}(t)\} = {{[{{{\psi }}_{{yk}}}(t{\text{'}}),k = \overline {1,{{M}_{t}}} ]}^{{\text{T}}}}$), значение производной $r$-го порядка в точке $t{\text{'}}$ (в этом случае ${\omega }\{ {{\Psi }_{y}}(t)\} = {{[{\psi }_{{yk}}^{{\left( r \right)}}(t{\text{'}}),k = \overline {1,\;{{M}_{t}}} ]}^{{\text{T}}}}$), значение определенного интеграла на отрезке $\left[ {{{t}_{0}},{{t}_{0}} + T} \right]$ (в этом случае ${\omega }\left\{ {{{\Psi }_{y}}(t)} \right\} = {{\left[ {\int\limits_{{{t}_{0}}}^{{{t}_{0}} + T} {{{{\psi }}_{{yk}}}()d\tau } ,k = \overline {1,{{M}_{t}}} } \right]}^{{\text{T}}}}$) и др. Несмещенные оценки этих и других НЛФ, инвариантные к сингулярной помехе, можно получить по результатам некорректных наблюдений (1.10), если в алгоритме (4.8), (4.9) раскрыть вектор ${\omega }\left\{ {{{\Psi }_{y}}\left( t \right)} \right\}$ с учетом специфики используемого НЛФ (по аналогии с приведенными примерами).

Предложенный алгоритм ОИНО можно выполнять на подвижной сетке (“скользящем окне”) $\left\{ {{{t}_{{j,r + l - 1}}}} \right\}_{{l = - {{{\bar {n}}}_{j}}/2}}^{{{{{\bar {n}}}_{j}}/2}} \in {{G}_{t}}$, где $\left( {{{{\bar {n}}}_{{j + 1}}} + 1} \right)$ – объем “скользящего окна”, индекс $r$ характеризует положение его центрального узла в момент времени ${{t}_{{jr}}} \in {{G}_{t}}$, а индекс $j$ позволяет занумеровать узлы “скользящего окна” при фиксированном $r$. Очевидно, что для формирования полномерного “скользящего окна” объемом ${{\bar {n}}_{{j + 1}}} + 1$ должны выполняться ограничения $1 + {{\bar {n}}_{j}}{\text{/}}2 \leqslant r \leqslant 1 + {{N}_{j}} - {{\bar {n}}_{j}}{\text{/}}2$ и $ - {{\bar {n}}_{j}}{\text{/}}2 \leqslant l \leqslant {{\bar {n}}_{j}}{\text{/}}2$, в противном случае алгоритм ОИНО реализуется на неполномерном “скользящем окне” нарастающего объема (это относится к краям временного отрезка $\left[ {{{t}_{0}},T} \right]$). Но и в этом случае имеется ограничение – этого объема должно быть достаточно для обеспечения невырожденности задачи оценивания на базе алгоритма (4.8), (4.9).

Известно, что точность оценивания сглаженного значения измеряемой функции наибольшая в середине “скользящего окна”, что также в полной мере относится к нахождению различных НЛФ (например, производных), которые соответствуют некоторой произвольной точке $t' \in {{G}_{t}}$, являющейся одним из узлов подвижной сетки.

Алгоритм ОИНО несложно реализовать и на адаптивном “скользящем окне”, объем которого меняется оптимальным образом с учетом условий наблюдения за ДС.

6. ИЛЛЮСТРАТИВНЫЙ ПРИМЕР

Опираясь на модель (5.1) общего вида, оценим эффективность предлагаемого метода применительно к ракете, движущейся в однородном поле силы тяжести, при линейном законе расхода массы и отсутствии сопротивления среды. Как показано в [20, с. 31–34], рассматриваемый случай является хорошим иллюстративным примером для тестирования различных численно-аналитических методов.

Криволинейное движение ракеты описывается системой дифференциальных уравнений (здесь время измеряется в секундах)

где ${{z}_{1}} = {{x}_{1}}$ – фазовая координата (безразмерная), соответствующая углу ${\beta }$ (рад) наклона вектора скорости к горизонту, ${{z}_{2}} = {{x}_{2}}$ – величина скорости ракеты (м/с), $g{}_{0}$ – ускорение силы тяжести (м/с²), ${\eta } = {{x}_{3}}$ – удельный расход массы ракеты (с^–1), ${{c}_{v}}$ – относительная скорость отбрасываемых частиц (м/с).

Далее полагаем, что ${{g}_{0}} = 9.80665$ (взято из справочника), ${{c}_{v}} = 2290$, ${{t}_{0}} = 0$, $T = 70$, ${{x}_{{i0}}} = {{x}_{i}}\left( 0 \right) \in {{G}_{{xi}}} = \left[ {{{c}_{i}},\;{{d}_{i}}} \right]$, $i = \overline {1,3} $, ${{c}_{1}} = 1.0107$, ${{d}_{1}} = 2.4362$, ${{c}_{2}} = 100$, ${{d}_{2}} = 200$, ${{c}_{3}} = r = 0.0060$, ${{d}_{3}} = p = 0.0100$, в качестве измеряемой функции $y$ используем угол ${\beta }$, т.е. $J = 1$, $y = {\beta }$ (с учетом скалярного выхода ДС в используемых далее формулах индекс $j$ опускается). Для построения численно-аналитического решения в формуле (5.1) принято ${{I}_{1}} = 2$ и ${{I}_{2}} = 1$, следовательно, $I = {{I}_{1}} + {{I}_{2}} = 3$, при этом: ${{f}_{1}}(t,x) = - {{g}_{0}}x_{2}^{{ - 1}}$, ${{f}_{2}}(t,x) = {{c}_{v}}{{x}_{3}}{{\left( {1 - {{x}_{3}}t} \right)}^{{ - 1}}} - {{g}_{0}}\operatorname{th} \left( {{{x}_{1}}} \right)$, ${{f}_{3}}(t,x) \equiv 0$, ${\text{m}} = \left( {{{m}_{1}},{{m}_{2}},{{m}_{3}}} \right)$, ${{M}_{{x1}}} = 4$, ${{M}_{{x2}}} = 6$, ${{M}_{{x3}}} = 5$, ${{M}_{x}} = {{M}_{{x1}}}{{M}_{{x2}}}{{M}_{{x3}}} = 120$, ${{M}_{t}} = 11$, ${{M}_{s}} = 2$, ${{{\psi }}_{{xk}}}(t) = {{L}_{{\left( k \right)}}}(t)$, ${{L}_{{\left( {\text{m}} \right)}}}\left( {{{x}_{0}}} \right) = {{L}_{{\left( {{{m}_{1}},{{m}_{2}},{{m}_{3}}} \right)}}}\left( {{{x}_{{01}}},{{x}_{{02}}},{{x}_{{03}}}} \right)$ = ${{L}_{{\left( {{{m}_{1}}} \right)}}}\left( {{{x}_{{10}}}} \right){{L}_{{\left( {{{m}_{2}}} \right)}}}\left( {{{x}_{{20}}}} \right){{L}_{{\left( {{{m}_{3}}} \right)}}}\left( {{{x}_{{30}}}} \right)$, ${{t}_{{\left( k \right)}}} = 7\left( {k - 1} \right)$, $k = \overline {1,11} $.

При расчетах погрешность операций над числами составляла $2.2 \times {{10}^{{ - 16}}}$, а результаты вычислений представлены с точностью до четвертого знака после запятой.

Этап 1. Для построения семейства ОИК использовался метод Рунге–Кутты 4 -го порядка с контролируемой точностью ${{10}^{{ - 5}}}$. С целью наглядности результатов моделирования, для оценки точности использовалась частная невязка ${{{\bar {\varepsilon }}}_{{xi}}} = \left| {{{x}_{i}}\left( {t,{{x}_{0}}} \right) - {{{\tilde {x}}}_{i}}\left( {t,{{x}_{0}}} \right)} \right|$, $i \in \left\{ {1,2} \right\}$, как функция от времени $t$ и начального условия ${{x}_{0}}$. Так, на фиг. 1 представлена зависимость частной невязки ${{{\bar {\varepsilon }}}_{{x1}}}$ от $t$ и ${{x}_{{10}}}$ для фиксированных значений ${{x}_{{20}}} = 104$ и ${{x}_{{30}}} = 0.0063$. При этом на исходной области $\left\{ {{{G}_{{x0}}},{{G}_{t}}} \right\}$ имеем общую невязку

что составляет $0.3237\% $. Данный максимум соответствует узлу

На усеченной области $\left\{ {{{{\bar {G}}}_{{x0}}},{{{\bar {G}}}_{t}}} \right\}$ (где ${{\bar {G}}_{{x0}}} = \left\{ {\left[ {1.3170,\;1.7454} \right],\;\left[ {120,\;180} \right],\;\left[ {0.0070,\;0.0090} \right]} \right\}$ и ${{G}_{t}} = \left[ {7,63} \right]$) получили общую невязку ${{{\varepsilon }}_{{x1}}} = 0.0004$, что составляет $0.1853\% $. Данный максимум соответствует узлу

На фиг. 2 представлена зависимость частной невязки ${{{\bar {\varepsilon }}}_{{x2}}}$ от $t$ и ${{x}_{{30}}}$ для фиксированных значений ${{x}_{{10}}} = 2.1466$ и ${{x}_{{20}}} = 124$. При этом на исходной области $\left\{ {{{G}_{{x0}}},{{G}_{t}}} \right\}$ имеем общую невязку ${{{\varepsilon }}_{{x2}}} = \mathop {\max }\limits_{{{x}_{0}}} {{{\varepsilon }}_{{x2}}}\left( {{{x}_{0}}} \right) = 4.0542$ (где ${{{\varepsilon }}_{{x2}}}\left( {{{x}_{0}}} \right) = \mathop {\max }\limits_t \left| {{{{\tilde {x}}}_{2}}\left( {t,{{x}_{0}}} \right) - {{x}_{2}}\left( {t,{{x}_{0}}} \right)} \right|$), что составляет $0.1744\% $. Данный максимум соответствует узлу

На усеченной области $\left\{ {{{{\bar {G}}}_{{x0}}},{{{\bar {G}}}_{t}}} \right\}$ получили общую невязку ${{{\varepsilon }}_{{x2}}} = 0.6355$, что составляет $0.0367\% $. Данный максимум соответствует узлу

Данные результаты наглядно демонстрируют возможность существенного повышения точности численно-аналитического решения дифференциального уравнения при усечении исходной области $\left\{ {{{G}_{{x0}}},{{G}_{t}}} \right\}$, что соответствует уменьшению влияния краевых эффектов.

Учитывая, что измеряемая функция ДС (6.1) может быть представлена в виде

на фиг. 3 представлен трехмерный график частной невязки ${{{\bar {\varepsilon }}}_{{\beta }}} = \left| {{\gamma }\left( {t,{{x}_{0}}} \right) - {\tilde {\gamma }*}\left( {t,{{x}_{0}}} \right)} \right|$ как функции от $t$ и ${{x}_{{10}}}$ для фиксированных значений ${{x}_{{20}}} = 108$ и ${{x}_{{30}}} = 0.0063$.

Соответственно на фиг. 4 рассмотрен случай зависимости ${{{\bar {\varepsilon }}}_{{\beta }}}$ от $t$ и ${{x}_{{30}}}$ для фиксированных значений ${{x}_{{10}}} = 2.1466$ и ${{x}_{{20}}} = 108$. При этом на исходной области $\left\{ {{{G}_{{x0}}},{{G}_{t}}} \right\}$ имеем общую невязку ${{{\varepsilon }}_{{\beta }}} = \mathop {\max }\limits_{{{x}_{0}}} {{{\varepsilon }}_{{\beta }}}({{x}_{0}}) = 0.0988$ (где ${{{\varepsilon }}_{{\beta }}}({{x}_{0}}) = \mathop {\max }\limits_t \left| {\tilde {\beta }\left( {t,{{x}_{0}}} \right) - \beta \left( {t,{{x}_{0}}} \right)} \right|$), что составляет $0.3303\% $. Данный максимум соответствует узлу

На усеченной области $\left\{ {{{{\bar {G}}}_{{x0}}},{{{\bar {G}}}_{t}}} \right\}$ получили общую невязку ${{{\varepsilon }}_{{\beta }}} = 0.0199$, что составляет $0.2197\% $. Данный максимум соответствует узлу

Соответственно для производной ${{y}^{{\left( 1 \right)}}} = {{{\beta }}^{{\left( 1 \right)}}}$ на исходной области $\left\{ {{{G}_{{x0}}},{{G}_{t}}} \right\}$ имеем общую невязку ${\varepsilon }_{{\beta }}^{{\left( 1 \right)}} = \mathop {\max }\limits_{{{x}_{0}}} {\varepsilon }_{{\beta }}^{{\left( 1 \right)}}({{x}_{0}}) = 5.4919$ (где ${\varepsilon }_{{\beta }}^{{\left( 1 \right)}}({{x}_{0}}) = \mathop {\max }\limits_t \left| {{{{\tilde {\beta }}}^{{\left( 1 \right)}}}\left( {t,{{x}_{0}}} \right) - {{\beta }^{{\left( 1 \right)}}}\left( {t,{{x}_{0}}} \right)} \right|$), что составляет $6.9520\% $. Данный максимум соответствует узлу

На усеченной области $\left\{ {{{{\bar {G}}}_{{x0}}},{{{\bar {G}}}_{t}}} \right\}$ получили общую невязку ${\varepsilon }_{{\beta }}^{{\left( 1 \right)}} = 3.1612$, что составляет $5.5860\% $. Данный максимум соответствует узлу

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

Этап 2. При реализации данного этапа использовались также следующие данные: $x_{{10}}^{'} = 1.7354$, $x_{{20}}^{'} = 112$, $x_{{30}}^{'} = 0.0800$, $s(t) = {{c}_{1}}{{{\theta }}_{1}}(t) + {{c}_{2}}{{{\theta }}_{2}}(t),$ ${{c}_{1}} = 5$, ${{{\theta }}_{1}}(t) = \sin (0.64t)$, ${{c}_{2}} = 15$, ${{{\theta }}_{2}}(t) = \exp ( - 0.07t)$, ${{K}_{\Xi }} = {\text{diag}}\left( {{\sigma },...,{\sigma }} \right),$ ${\sigma } = {\pi /}360$, $\left\{ {{{t}_{n}}} \right\}_{{n = 0}}^{N}$, $N = 70$, ${{t}_{n}} - {{t}_{{n - 1}}} = 1$, $\bar {n} = 20$, $11 \leqslant r \leqslant 61$. Алгоритм ОИНО реализовывался на полномерном “скользящем окне” с 21 узлом, при этом первая оценка измеряемого параметра и его производной относятся к текущему моменту времени ${{t}_{{20}}} = 20$, а последняя – к моменту ${{t}_{{70}}} = 70$. Именно при ${{t}_{{20}}} = 20$ было сформировано первое полномерное “скользящее окно”, а при ${{t}_{{70}}} = 70$ – последнее. Поскольку все оценки привязываются к середине “скользящего окна”, то первая полученная оценка соответствовала моменту ${{t}_{{10}}} = 10$, а последняя – моменту ${{t}_{{60}}} = 60$. Следовательно, приводимые ниже характеристики точности относятся к усеченному отрезку времени ${{\bar {G}}_{t}} = \left[ {10,60} \right]$.

На фиг. 5 приведены графики зависимости, характеризующие сингулярную помеху $s\left( t \right)$.

С использованием датчика случайных чисел и принятых исходных данных было сформировано уравнение наблюдения, при этом на фиг. 6 представлены графики истинного значения измеряемой функции $y = {\beta }$ и наблюдения $h = y + s + {\xi }$.

Начиная с момента времени ${{t}_{n}} = 20$ для каждого положения “скользящего окна” (всего 51 положение, т.е. $11 \leqslant r \leqslant 61$) формировались текущие оценки ${\beta }_{n}^{*}$ и ${\beta }_{n}^{{\left( 1 \right)}}{\text{*}}$, соответствующие центральному узлу. Для описания выходной координаты $y(t) = {\beta }(t)$ в пределах “скользящего окна” использовалась линейная комбинация из трех первых полиномов Лежандра, т.е. в формуле (4.1) принималось:

Для перехода от параметра ${\tau } \in \left[ { - 1,1} \right]$ к временной координате $t \in [{{t}_{{j,r - {{{\bar {n}}}_{j}}/2 - 1}}},{{t}_{{j,r + {{{\bar {n}}}_{j}}/2 - 1}}}] \in {{G}_{t}}$ применялось линейное преобразование: $t = {{2}^{{ - 1}}}[{{t}_{{j,r - {{{\bar {n}}}_{j}}/2 - 1}}} + {{t}_{{j,r + {{{\bar {n}}}_{j}}/2 - 1}}} - {\tau }({{t}_{{j,r - {{{\bar {n}}}_{j}}/2 - 1}}} - {{t}_{{j,r + {{{\bar {n}}}_{j}}/2 - 1}}})]$.

Ниже представлены графики зависимости частной невязки ${\bar {\varepsilon }}_{{\beta }}^{ * }$ оценивания выходной координаты (фиг. 7) и ее производной ${\bar {\varepsilon }}_{{\beta }}^{{\left( 1 \right)}}{\text{*}}$ (фиг. 8) от времени при фиксированном начальном условии $x_{0}^{'} = {{\left[ {1.7354,112,0.0800} \right]}^{{\text{T}}}}$, при этом общая невязка для выходной координаты составляет ${\varepsilon }_{{\beta }}^{*} = \mathop {\max }\limits_t {\bar {\varepsilon }}_{{\beta }}^{*} = 0.8162$ (соответственно $1.5717\% $), а для производной ${\varepsilon }_{{\beta }}^{{\left( 1 \right)}}{\text{*}} = \mathop {\max }\limits_t {\bar {\varepsilon }}_{{\beta }}^{{\left( 1 \right)}}{\text{*}} = 0.0518$ (соответственно $5.0601\% $).

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

ЗАКЛЮЧЕНИЕ

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

Развитый метод не предполагает выполнения процедур линеаризации (даже если ДС и ее выход являются нелинейными), не требует расширения пространства состояний и не связан с необходимостью самостоятельного расчета спектральных коэффициентов в линейных комбинациях выходных траекторий и сингулярной помехи. Метод допускает возможность комбинирования с традиционными подходами к решению прикладных задач, связанных с оптимальной обработкой наблюдений, идентификацией и планированием измерительного эксперимента [4], [9], [11], [15], [18], [19].

Выражаю свою благодарность моим аспирантам П.Ю. Раду и А.Г. Кондрашову за помощь в проведении вычислительного эксперимента.